21 #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDSOLVER_HPP_
22 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDSOLVER_HPP_
32 template<
typename ... SOLVERS >
44 Group *
const parent )
47 forEachArgInTuple(
m_solvers, [&](
auto solver,
auto idx )
49 using SolverType = TYPEOFPTR( solver );
50 string const key = SolverType::coupledSolverAttributePrefix() +
"SolverName";
52 setRTTypeName( rtTypes::CustomTypes::groupNameRef ).
54 setDescription(
"Name of the " + SolverType::coupledSolverAttributePrefix() +
" solver used by the coupled solver" );
60 addLogLevel< logInfo::Coupling >();
82 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto idx )
84 using SolverPtr = TYPEOFREF( solver );
85 using SolverType = TYPEOFPTR( SolverPtr {} );
86 auto const & solverName =
m_names[idx()];
87 auto const & solverType = LvArray::system::demangleType< SolverType >();
88 solver = this->
getParent().template getGroupPointer< SolverType >( solverName );
90 GEOS_FMT(
"{}: Could not find solver '{}' of type {}",
92 solverName, solverType ),
95 GEOS_FMT(
"{}: found {} solver named {}",
96 getName(), solver->getCatalogName(), solverName ) );
127 {
GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); }
140 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
142 solver->setupDofs( domain, dofManager );
153 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
155 solver->implicitStepSetup( time_n, dt, domain );
164 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
166 solver->implicitStepComplete( time_n, dt, domain );
182 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
184 solver->assembleSystem( time_n, dt, domain, dofManager, localMatrix, localRhs );
194 real64 const scalingFactor,
198 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
200 solver->applySystemSolution( dofManager, localSolution, scalingFactor, dt, domain );
207 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
209 solver->updateState( domain );
216 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
218 solver->resetStateToBeginningOfStep( domain );
227 int const cycleNumber,
253 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
257 solver->updateAndWriteConvergenceStep( time_n, dt, cycleNumber, iteration );
270 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
272 real64 const singlePhysicsNorm = solver->calculateResidualNorm( time_n, dt, domain, dofManager, localRhs );
273 norm += singlePhysicsNorm * singlePhysicsNorm;
287 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
289 solver->applyBoundaryConditions( time_n, dt, domain, dofManager, localMatrix, localRhs );
297 real64 const scalingFactor )
override
299 bool validSolution =
true;
300 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
302 bool const validSinglePhysicsSolution = solver->checkSystemSolution( domain, dofManager, localSolution, scalingFactor );
303 if( !validSinglePhysicsSolution )
305 GEOS_LOG_RANK_0( GEOS_FMT(
" {}/{}: Solution check failed. Newton loop terminated.",
getName(), solver->getName()) );
307 validSolution = validSolution && validSinglePhysicsSolution;
309 return validSolution;
318 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
320 real64 const singlePhysicsScalingFactor = solver->scalingForSystemSolution( domain, dofManager, localSolution );
321 scalingFactor = LvArray::math::min( scalingFactor, singlePhysicsScalingFactor );
323 return scalingFactor;
332 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
334 real64 const singlePhysicsNextDt =
335 solver->setNextDt( currentTime, currentDt, domain );
336 nextDt = LvArray::math::min( singlePhysicsNextDt, nextDt );
344 real64 const eventProgress,
347 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
349 solver->cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain );
358 bool isConverged =
true;
359 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
361 isConverged &= solver->checkSequentialSolutionIncrements( domain );
367 integer const configurationLoopIter )
override
370 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
372 result &= solver->updateConfiguration( domain, configurationLoopIter );
379 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
381 solver->outputConfigurationStatistics( domain );
387 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
389 solver->resetConfigurationToBeginningOfStep( domain );
396 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
398 result &= solver->resetConfigurationToDefault( domain );
416 int const cycleNumber,
441 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
443 if( meshModificationTimestamp > solver->getSystemSetupTimestamp() )
445 solver->setupSystem( domain,
446 solver->getDofManager(),
447 solver->getLocalMatrix(),
448 solver->getSystemRhs(),
449 solver->getSystemSolution() );
450 solver->setSystemSetupTimestamp( meshModificationTimestamp );
461 bool isConverged =
false;
468 for( dtAttempt = 0; dtAttempt < maxNumberDtCuts; ++dtAttempt )
473 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
475 solver->resetStateToBeginningOfStep( domain );
476 solver->getIterationStats().resetCurrentTimeStepStatistics();
488 startSequentialIteration( iter, domain );
491 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto idx )
494 GEOS_FMT(
" Iteration {:2}: {}", iter + 1, solver->getName() ) );
495 real64 solverDt = solver->nonlinearImplicitStep( time_n,
503 solver->saveSequentialIterationState( domain );
508 if( solverDt < stepDt )
516 isConverged = checkSequentialConvergence( cycleNumber,
531 finishSequentialIteration( iter, domain );
538 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
540 solver->getIterationStats().iterateTimeStepStatistics();
548 stepDt *= dtCutFactor;
554 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
556 solver->getIterationStats().updateTimeStepCut();
571 GEOS_ERROR(
"Nonconverged solutions not allowed. Terminating..." );
592 virtual bool checkSequentialConvergence(
integer const cycleNumber,
599 bool isConverged =
true;
607 GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, GEOS_FMT(
" Iteration {:2}: outer-loop convergence check", iter + 1 ) );
614 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
617 solver->getLocalMatrix().toViewConstSizes().zero();
618 solver->getSystemRhs().zero();
619 arrayView1d< real64 >
const localRhs = solver->getSystemRhs().open();
622 solver->assembleSystem( time_n,
625 solver->getDofManager(),
626 solver->getLocalMatrix().toViewConstSizes(),
628 solver->applyBoundaryConditions( time_n,
631 solver->getDofManager(),
632 solver->getLocalMatrix().toViewConstSizes(),
634 solver->getSystemRhs().close();
637 real64 const singlePhysicsNorm =
638 solver->calculateResidualNorm( time_n,
641 solver->getDofManager(),
642 solver->getSystemRhs().values() );
643 residualNorm += singlePhysicsNorm * singlePhysicsNorm;
647 residualNorm = sqrt( residualNorm );
649 GEOS_FMT(
" ( R ) = ( {:4.2e} )", residualNorm ) );
653 isConverged = ( residualNorm < params.
m_newtonTol );
659 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
661 NonlinearSolverParameters
const & singlePhysicsParams = solver->getNonlinearSolverParameters();
662 if( singlePhysicsParams.m_numNewtonIterations > singlePhysicsParams.m_minIterNewton )
680 GEOS_FMT(
"***** The iterative coupling has converged in {} iteration(s) *****", iter + 1 ) );
696 GEOS_FMT(
"{}: line search is not supported by the coupled solver when {} is set to `{}`. Please set {} to `{}` to remove this error",
698 NonlinearSolverParameters::viewKeysStruct::couplingTypeString(),
700 NonlinearSolverParameters::viewKeysStruct::lineSearchActionString(),
711 validateNonlinearAcceleration();
715 virtual void validateNonlinearAcceleration()
717 GEOS_THROW ( GEOS_FMT(
"{}: Nonlinear acceleration {} is not supported by {} solver '{}'",
718 getWrapperDataContext( NonlinearSolverParameters::viewKeysStruct::nonlinearAccelerationTypeString() ),
727 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
733 virtual void startSequentialIteration(
integer const & iter,
739 virtual void finishSequentialIteration(
integer const & iter,
740 DomainPartition & domain )
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
#define GEOS_THROW(msg, TYPE)
Throw an exception.
#define GEOS_ERROR(msg)
Raise a hard error and terminate the program.
#define GEOS_LOG_RANK_0(msg)
Log a message on screen on rank 0.
#define GEOS_THROW_IF(EXP, msg, TYPE)
Conditionally throw an exception.
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
virtual void assembleCouplingTerms(real64 const time_n, real64 const dt, DomainPartition const &domain, DofManager const &dofManager, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs)
Utility function to compute coupling terms.
CoupledSolver & operator=(CoupledSolver const &)=delete
deleted assignment operator
virtual real64 fullyCoupledSolverStep(real64 const &time_n, real64 const &dt, int const cycleNumber, DomainPartition &domain)
Fully coupled solution approach solution step.
virtual bool resetConfigurationToDefault(DomainPartition &domain) const override
resets the configuration to the default value.
CoupledSolver(const string &name, Group *const parent)
main constructor for CoupledSolver Objects
virtual void outputConfigurationStatistics(DomainPartition const &domain) const override
CoupledSolver(CoupledSolver &&)=default
default move constructor
virtual void resetConfigurationToBeginningOfStep(DomainPartition &domain) override
resets the configuration to the beginning of the time-step.
virtual real64 sequentiallyCoupledSolverStep(real64 const &time_n, real64 const &dt, integer const cycleNumber, DomainPartition &domain)
Sequentially coupled solver step. It solves a nonlinear system of equations using a sequential approa...
CoupledSolver(CoupledSolver const &)=delete
deleted copy constructor
virtual void synchronizeNonlinearSolverParameters() override
synchronize the nonlinear solver parameters.
virtual void setupCoupling(DomainPartition const &domain, DofManager &dofManager) const
Utility function to set the coupling between degrees of freedom.
virtual bool checkSequentialSolutionIncrements(DomainPartition &domain) const override
Check if the solution increments are ok to use.
virtual bool updateConfiguration(DomainPartition &domain, integer const configurationLoopIter) override
updates the configuration (if needed) based on the state after a converged Newton loop.
virtual void postInputInitialization() override
std::array< string, sizeof...(SOLVERS) > m_names
Names of the single-physics solvers.
std::tuple< SOLVERS *... > m_solvers
Pointers of the single-physics solvers.
CoupledSolver & operator=(CoupledSolver &&)=delete
deleted move operator
void setSubSolvers()
Utility function to set the subsolvers pointers using the names provided by the user.
virtual void mapSolutionBetweenSolvers(DomainPartition &domain, integer const solverType)
Maps the solution obtained from one solver to the fields used by the other solver(s)
The DoFManager is responsible for allocating global dofs, constructing sparsity patterns,...
Partition of the decomposed physical domain. It also manages the connexion information to its neighbo...
void updateNonlinearIteration(integer const numLinearIterations)
Tell the solverStatistics that we have done a newton iteration.
void updateTimeStepCut()
Tell the solverStatistics that we cut the time step and we increment the cumulative counters for disc...
SequentialConvergenceCriterion sequentialConvergenceCriterion() const
Getter for the sequential convergence criterion.
integer m_allowNonConverged
Flag to allow for a non-converged nonlinear solution and continue with the problem.
@ FullyImplicit
Fully-implicit coupling.
@ Sequential
Sequential coupling.
real64 m_newtonTol
The tolerance for the nonlinear convergence check.
NonlinearAccelerationType m_nonlinearAccelerationType
Type of nonlinear acceleration for sequential solver.
@ None
Do not use line search.
integer m_maxIterNewton
The maximum number of nonlinear iterations that are allowed.
real64 m_timeStepCutFactor
Factor by which the time step will be cut if a timestep cut is required.
integer m_numNewtonIterations
The number of nonlinear iterations that have been exectued.
integer m_numTimeStepAttempts
Number of times that the time-step had to be cut.
integer m_maxTimeStepCuts
Max number of time step cuts.
CouplingType couplingType() const
Getter for the coupling type.
LineSearchAction m_lineSearchAction
Flag to apply a line search.
@ ResidualNorm
convergence achieved when the residual drops below a given norm
@ NumberOfNonlinearIterations
convergence achieved when the subproblems convergence is achieved in less than minNewtonIteration
@ SolutionIncrements
convergence achieved when the solution increments are small enough
integer m_subcyclingOption
Flag to specify whether subcycling is allowed or not in sequential schemes.
Base class for all physics solvers.
virtual string getCatalogName() const =0
IterationsStatistics & getIterationStats()
integer m_numTimestepsSinceLastDtCut
Number of cycles since last timestep cut.
virtual void cleanup(real64 const time_n, integer const cycleNumber, integer const eventCounter, real64 const eventProgress, DomainPartition &domain) override
Called as the code exits the main run loop.
Timestamp getMeshModificationTimestamp(DomainPartition &domain) const
getter for the timestamp of the mesh modification on the mesh levels
virtual void postInputInitialization() override
ConvergenceStatistics & getConvergenceStats()
integer m_writeStatisticsCSV
NonlinearSolverParameters & getNonlinearSolverParameters()
accessor for the nonlinear solver parameters.
NonlinearSolverParameters m_nonlinearSolverParameters
Nonlinear solver parameters.
Wrapper< TBASE > & registerWrapper(string const &name, wrapperMap::KeyIndex::index_type *const rkey=nullptr)
Create and register a Wrapper around a new object.
DataContext const & getDataContext() const
string const & getName() const
Get group name.
Group & getParent()
Access the group's parent.
DataContext const & getWrapperDataContext(KEY key) const
#define GEOS_LOG_LEVEL_RANK_0(logInfoStruct, msg)
Output messages (only on rank 0) based on current Group's log level.
virtual void implicitStepSetup(real64 const &time_n, real64 const &dt, DomainPartition &domain) override
function to perform setup for implicit timestep
virtual void cleanup(real64 const time_n, integer const cycleNumber, integer const eventCounter, real64 const eventProgress, DomainPartition &domain) override
Called as the code exits the main run loop.
virtual real64 scalingForSystemSolution(DomainPartition &domain, DofManager const &dofManager, arrayView1d< real64 const > const &localSolution) override
Function to determine if the solution vector should be scaled back in order to maintain a known const...
virtual void updateAndWriteConvergenceStep(real64 const &time_n, real64 const &dt, integer const cycleNumber, integer const iteration) override
Update the convergence information and write then into a CSV file.
virtual void updateState(DomainPartition &domain) override
Recompute all dependent quantities from primary variables (including constitutive models)
virtual void implicitStepComplete(real64 const &time_n, real64 const &dt, DomainPartition &domain) override
perform cleanup for implicit timestep
virtual void applyBoundaryConditions(real64 const time_n, real64 const dt, DomainPartition &domain, DofManager const &dofManager, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs) override
apply boundary condition to system
virtual real64 scalingForSystemSolution(DomainPartition &domain, DofManager const &dofManager, arrayView1d< real64 const > const &localSolution)
Function to determine if the solution vector should be scaled back in order to maintain a known const...
virtual real64 calculateResidualNorm(real64 const &time_n, real64 const &dt, DomainPartition const &domain, DofManager const &dofManager, arrayView1d< real64 const > const &localRhs) override
calculate the norm of the global system residual
real64 solverStep(real64 const &time_n, real64 const &dt, int const cycleNumber, DomainPartition &domain) override final
virtual real64 setNextDt(real64 const ¤tTime, real64 const ¤tDt, DomainPartition &domain) override
function to set the next time step size
virtual real64 solverStep(real64 const &time_n, real64 const &dt, integer const cycleNumber, DomainPartition &domain)
entry function to perform a solver step
void setupDofs(DomainPartition const &domain, DofManager &dofManager) const override
Populate degree-of-freedom manager with fields relevant to this solver.
virtual void applySystemSolution(DofManager const &dofManager, arrayView1d< real64 const > const &localSolution, real64 const scalingFactor, real64 const dt, DomainPartition &domain) override
Function to apply the solution vector to the state.
virtual real64 setNextDt(real64 const ¤tTime, real64 const ¤tDt, DomainPartition &domain)
function to set the next time step size
virtual void assembleSystem(real64 const time_n, real64 const dt, DomainPartition &domain, DofManager const &dofManager, CRSMatrixView< real64, globalIndex const > const &localMatrix, arrayView1d< real64 > const &localRhs) override
function to assemble the linear system matrix and rhs
virtual bool checkSystemSolution(DomainPartition &domain, DofManager const &dofManager, arrayView1d< real64 const > const &localSolution, real64 const scalingFactor) override
Function to check system solution for physical consistency and constraint violation.
virtual void resetStateToBeginningOfStep(DomainPartition &domain) override
reset state of physics back to the beginning of the step.
@ FALSE
Not read from input.
@ REQUIRED
Required in input.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
unsigned long long int Timestamp
Timestamp type (used to perform actions such a sparsity pattern computation after mesh modifications)
std::string string
String type.
double real64
64-bit floating point type.
LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
int integer
Signed integer type.
Provides enum <-> string conversion facilities.
static constexpr char const * discretizationString()