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,
257 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
259 real64 const singlePhysicsNorm = solver->calculateResidualNorm( time_n, dt, domain, dofManager, localRhs );
260 norm += singlePhysicsNorm * singlePhysicsNorm;
274 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
276 solver->applyBoundaryConditions( time_n, dt, domain, dofManager, localMatrix, localRhs );
284 real64 const scalingFactor )
override
286 bool validSolution =
true;
287 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
289 bool const validSinglePhysicsSolution = solver->checkSystemSolution( domain, dofManager, localSolution, scalingFactor );
290 if( !validSinglePhysicsSolution )
292 GEOS_LOG_RANK_0( GEOS_FMT(
" {}/{}: Solution check failed. Newton loop terminated.",
getName(), solver->getName()) );
294 validSolution = validSolution && validSinglePhysicsSolution;
296 return validSolution;
304 real64 scalingFactor = 1e9;
305 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
307 real64 const singlePhysicsScalingFactor = solver->scalingForSystemSolution( domain, dofManager, localSolution );
308 scalingFactor = LvArray::math::min( scalingFactor, singlePhysicsScalingFactor );
310 return scalingFactor;
319 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
321 real64 const singlePhysicsNextDt =
322 solver->setNextDt( currentTime, currentDt, domain );
323 nextDt = LvArray::math::min( singlePhysicsNextDt, nextDt );
331 real64 const eventProgress,
334 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
336 solver->cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain );
345 bool isConverged =
true;
346 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
348 isConverged &= solver->checkSequentialSolutionIncrements( domain );
356 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
358 result &= solver->updateConfiguration( domain );
365 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
367 solver->outputConfigurationStatistics( domain );
373 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
375 solver->resetConfigurationToBeginningOfStep( domain );
382 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
384 result &= solver->resetConfigurationToDefault( domain );
402 int const cycleNumber,
420 int const cycleNumber,
427 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
429 if( meshModificationTimestamp > solver->getSystemSetupTimestamp() )
431 solver->setupSystem( domain,
432 solver->getDofManager(),
433 solver->getLocalMatrix(),
434 solver->getSystemRhs(),
435 solver->getSystemSolution() );
436 solver->setSystemSetupTimestamp( meshModificationTimestamp );
447 bool isConverged =
false;
454 for( dtAttempt = 0; dtAttempt < maxNumberDtCuts; ++dtAttempt )
459 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
461 solver->resetStateToBeginningOfStep( domain );
462 solver->getSolverStatistics().initializeTimeStepStatistics();
475 startSequentialIteration( iter, domain );
478 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto idx )
481 GEOS_FMT(
" Iteration {:2}: {}", iter + 1, solver->getName() ) );
482 real64 solverDt = solver->nonlinearImplicitStep( time_n,
490 solver->saveSequentialIterationState( domain );
495 if( solverDt < stepDt )
504 isConverged = checkSequentialConvergence( iter,
518 finishSequentialIteration( iter, domain );
525 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
527 solver->getSolverStatistics().saveTimeStepStatistics();
535 stepDt *= dtCutFactor;
541 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
543 solver->getSolverStatistics().logTimeStepCut();
558 GEOS_ERROR(
"Nonconverged solutions not allowed. Terminating..." );
579 virtual bool checkSequentialConvergence(
int const & iter,
585 bool isConverged =
true;
593 GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, GEOS_FMT(
" Iteration {:2}: outer-loop convergence check", iter + 1 ) );
600 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
603 solver->getLocalMatrix().toViewConstSizes().zero();
604 solver->getSystemRhs().zero();
605 arrayView1d< real64 >
const localRhs = solver->getSystemRhs().open();
608 solver->assembleSystem( time_n,
611 solver->getDofManager(),
612 solver->getLocalMatrix().toViewConstSizes(),
614 solver->applyBoundaryConditions( time_n,
617 solver->getDofManager(),
618 solver->getLocalMatrix().toViewConstSizes(),
620 solver->getSystemRhs().close();
623 real64 const singlePhysicsNorm =
624 solver->calculateResidualNorm( time_n,
627 solver->getDofManager(),
628 solver->getSystemRhs().values() );
629 residualNorm += singlePhysicsNorm * singlePhysicsNorm;
633 residualNorm = sqrt( residualNorm );
635 GEOS_FMT(
" ( R ) = ( {:4.2e} )", residualNorm ) );
636 isConverged = ( residualNorm < params.
m_newtonTol );
642 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
644 NonlinearSolverParameters
const & singlePhysicsParams = solver->getNonlinearSolverParameters();
645 if( singlePhysicsParams.m_numNewtonIterations > singlePhysicsParams.m_minIterNewton )
663 GEOS_FMT(
"***** The iterative coupling has converged in {} iteration(s) *****", iter + 1 ) );
677 GEOS_FMT(
"{}: line search is not supported by the coupled solver when {} is set to `{}`. Please set {} to `{}` to remove this error",
679 NonlinearSolverParameters::viewKeysStruct::couplingTypeString(),
681 NonlinearSolverParameters::viewKeysStruct::lineSearchActionString(),
691 validateNonlinearAcceleration();
694 virtual void validateNonlinearAcceleration()
696 GEOS_THROW ( GEOS_FMT(
"{}: Nonlinear acceleration {} is not supported by {} solver '{}'",
697 getWrapperDataContext( NonlinearSolverParameters::viewKeysStruct::nonlinearAccelerationTypeString() ),
706 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
712 virtual void startSequentialIteration(
integer const & iter,
718 virtual void finishSequentialIteration(
integer const & iter,
719 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.
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 real64 sequentiallyCoupledSolverStep(real64 const &time_n, real64 const &dt, int const cycleNumber, DomainPartition &domain)
Sequentially coupled solver step. It solves a nonlinear system of equations using a sequential approa...
virtual bool checkSequentialSolutionIncrements(DomainPartition &domain) const override
Check if the solution increments are ok to use.
virtual void postInputInitialization() override
virtual bool updateConfiguration(DomainPartition &domain) override
updates the configuration (if needed) based on the state after a converged Newton loop.
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...
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
SolverStatistics m_solverStatistics
Solver statistics.
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
NonlinearSolverParameters & getNonlinearSolverParameters()
accessor for the nonlinear solver parameters.
NonlinearSolverParameters m_nonlinearSolverParameters
Nonlinear solver parameters.
void logTimeStepCut()
Tell the solverStatistics that there is a time step cut.
void logNonlinearIteration(integer const numLinearIterations)
Tell the solverStatistics that we are doing a nonlinear iteration.
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 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 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()