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,
254 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
256 solver->updateAndWriteConvergenceStep( time_n, dt, cycleNumber, iteration );
268 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
270 real64 const singlePhysicsNorm = solver->calculateResidualNorm( time_n, dt, domain, dofManager, localRhs );
271 norm += singlePhysicsNorm * singlePhysicsNorm;
285 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
287 solver->applyBoundaryConditions( time_n, dt, domain, dofManager, localMatrix, localRhs );
295 real64 const scalingFactor )
override
297 bool validSolution =
true;
298 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
300 bool const validSinglePhysicsSolution = solver->checkSystemSolution( domain, dofManager, localSolution, scalingFactor );
301 if( !validSinglePhysicsSolution )
303 GEOS_LOG_RANK_0( GEOS_FMT(
" {}/{}: Solution check failed. Newton loop terminated.",
getName(), solver->getName()) );
305 validSolution = validSolution && validSinglePhysicsSolution;
307 return validSolution;
316 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
318 real64 const singlePhysicsScalingFactor = solver->scalingForSystemSolution( domain, dofManager, localSolution );
319 scalingFactor = LvArray::math::min( scalingFactor, singlePhysicsScalingFactor );
321 return scalingFactor;
330 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
332 real64 const singlePhysicsNextDt =
333 solver->setNextDt( currentTime, currentDt, domain );
334 nextDt = LvArray::math::min( singlePhysicsNextDt, nextDt );
342 real64 const eventProgress,
345 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
347 solver->cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain );
356 bool isConverged =
true;
357 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
359 isConverged &= solver->checkSequentialSolutionIncrements( domain );
365 integer const configurationLoopIter )
override
368 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
370 result &= solver->updateConfiguration( domain, configurationLoopIter );
377 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
379 solver->outputConfigurationStatistics( domain );
385 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
387 solver->resetConfigurationToBeginningOfStep( domain );
394 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
396 result &= solver->resetConfigurationToDefault( domain );
404 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
407 solver->synchronizeNonlinearSolverParameters();
424 int const cycleNumber,
449 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
451 if( meshModificationTimestamp > solver->getSystemSetupTimestamp() )
453 solver->setupSystem( domain,
454 solver->getDofManager(),
455 solver->getLocalMatrix(),
456 solver->getSystemRhs(),
457 solver->getSystemSolution() );
458 solver->setSystemSetupTimestamp( meshModificationTimestamp );
469 bool isConverged =
false;
476 for( dtAttempt = 0; dtAttempt < maxNumberDtCuts; ++dtAttempt )
481 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
483 solver->resetStateToBeginningOfStep( domain );
484 solver->getIterationStats().resetCurrentTimeStepStatistics();
496 startSequentialIteration( iter, domain );
499 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto idx )
502 GEOS_FMT(
" Iteration {:2}: {}", iter + 1, solver->getName() ) );
503 real64 solverDt = solver->nonlinearImplicitStep( time_n,
511 solver->saveSequentialIterationState( domain );
516 if( solverDt < stepDt )
524 isConverged = checkSequentialConvergence( cycleNumber,
539 finishSequentialIteration( iter, domain );
546 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
548 solver->getIterationStats().iterateTimeStepStatistics();
556 stepDt *= dtCutFactor;
562 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
564 solver->getIterationStats().updateTimeStepCut();
579 GEOS_ERROR(
"Nonconverged solutions not allowed. Terminating..." );
600 virtual bool checkSequentialConvergence(
integer const cycleNumber,
607 bool isConverged =
true;
615 GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, GEOS_FMT(
" Iteration {:2}: outer-loop convergence check", iter + 1 ) );
622 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
625 solver->getLocalMatrix().toViewConstSizes().zero();
626 solver->getSystemRhs().zero();
627 arrayView1d< real64 >
const localRhs = solver->getSystemRhs().open();
630 solver->assembleSystem( time_n,
633 solver->getDofManager(),
634 solver->getLocalMatrix().toViewConstSizes(),
636 solver->applyBoundaryConditions( time_n,
639 solver->getDofManager(),
640 solver->getLocalMatrix().toViewConstSizes(),
642 solver->getSystemRhs().close();
645 real64 const singlePhysicsNorm =
646 solver->calculateResidualNorm( time_n,
649 solver->getDofManager(),
650 solver->getSystemRhs().values() );
651 residualNorm += singlePhysicsNorm * singlePhysicsNorm;
655 residualNorm = sqrt( residualNorm );
657 GEOS_FMT(
" ( R ) = ( {:4.2e} )", residualNorm ) );
661 isConverged = ( residualNorm < params.
m_newtonTol );
667 forEachArgInTuple(
m_solvers, [&](
auto & solver,
auto )
669 NonlinearSolverParameters
const & singlePhysicsParams = solver->getNonlinearSolverParameters();
670 if( singlePhysicsParams.m_numNewtonIterations > singlePhysicsParams.m_minIterNewton )
688 GEOS_FMT(
"***** The iterative coupling has converged in {} iteration(s) *****", iter + 1 ) );
704 GEOS_FMT(
"{}: line search is not supported by the coupled solver when {} is set to `{}`. Please set {} to `{}` to remove this error",
706 NonlinearSolverParameters::viewKeysStruct::couplingTypeString(),
708 NonlinearSolverParameters::viewKeysStruct::lineSearchActionString(),
714 validateNonlinearAcceleration();
718 virtual void validateNonlinearAcceleration()
720 GEOS_THROW ( GEOS_FMT(
"{}: Nonlinear acceleration {} is not supported by {} solver '{}'",
721 getWrapperDataContext( NonlinearSolverParameters::viewKeysStruct::nonlinearAccelerationTypeString() ),
738 virtual void startSequentialIteration(
integer const & iter,
744 virtual void finishSequentialIteration(
integer const & iter,
745 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 void initializePreSubGroups() override
Called by Initialize() prior to initializing sub-Groups.
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()
NonlinearSolverParameters & getNonlinearSolverParameters()
accessor for the nonlinear solver parameters.
NonlinearSolverParameters m_nonlinearSolverParameters
Nonlinear solver parameters.
virtual void initializePreSubGroups()
Called by Initialize() prior to initializing sub-Groups.
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 updateAndWriteConvergenceStep(real64 const &time_n, real64 const &dt, integer const cycleNumber, integer const iteration)
Update the convergence information and write then into a CSV file.
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()