Time step cut support for sequential solver (#2861)

Currently the simulation would stop if outer loop didn't converge. This commit adds time step cut logic.
Bonus:
- print next dt in sub-timestepping process
- add `setNextDtBasedOnNewtonIter` for coupled solver to avoid subsolvers doing a lot of iters
- when pressure change is computed - switch from relative to absolute for pres < 1

src/coreComponents/physicsSolvers/SolverBase.cpp

@@ -273,7 +273,7 @@ bool SolverBase::execute( real64 const time_n,
 
     if( getLogLevel() >= 1 && dtRemaining > 0.0 )
     {
-      GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "{}: sub-step = {}, accepted dt = {}, remaining dt = {}", getName(), subStep, dtAccepted, dtRemaining ) );
+      GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "{}: sub-step = {}, accepted dt = {}, next dt = {}, remaining dt = {}", getName(), subStep, dtAccepted, nextDt, dtRemaining ) );
     }
   }
 
@@ -820,7 +820,7 @@ real64 SolverBase::nonlinearImplicitStep( real64 const & time_n,
     }
     else
     {
-      GEOS_ERROR( getDataContext() << ": Nonconverged solutions not allowed. Terminating..." );
+      GEOS_ERROR( "Nonconverged solutions not allowed. Terminating..." );
     }
   }
 

src/coreComponents/physicsSolvers/SolverBase.hpp

@@ -153,7 +153,7 @@ class SolverBase : public ExecutableGroup
    * @param[in] currentDt the current time step size
    * @return the prescribed time step size
    */
-  real64 setNextDtBasedOnNewtonIter( real64 const & currentDt );
+  virtual real64 setNextDtBasedOnNewtonIter( real64 const & currentDt );
 
   /**
    * @brief function to set the next dt based on state change

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp

@@ -2035,8 +2035,9 @@ real64 CompositionalMultiphaseBase::setNextDtBasedOnStateChange( real64 const &
       {
         if( ghostRank[ei] < 0 )
         {
-          subRegionMaxPresChange.max( LvArray::math::abs( pres[ei] - pres_n[ei] ) / LvArray::math::max( pres_n[ei], LvArray::NumericLimits< real64 >::epsilon ) );
-          subRegionMaxTempChange.max( LvArray::math::abs( temp[ei] - temp_n[ei] ) / LvArray::math::max( temp_n[ei], LvArray::NumericLimits< real64 >::epsilon ) );
+          // switch from relative to absolute when pressure less than 1
+          subRegionMaxPresChange.max( LvArray::math::abs( pres[ei] - pres_n[ei] ) / LvArray::math::max( LvArray::math::abs( pres_n[ei] ), 1.0 ) );
+          subRegionMaxTempChange.max( LvArray::math::abs( temp[ei] - temp_n[ei] ) / LvArray::math::max( LvArray::math::abs( temp_n[ei] ), 1.0 ) );
           for( integer ip = 0; ip < numPhase; ++ip )
           {
             subRegionMaxPhaseVolFracChange.max( LvArray::math::abs( phaseVolFrac[ei][ip] - phaseVolFrac_n[ei][ip] ) );

src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp

@@ -318,6 +318,18 @@ class CoupledSolver : public SolverBase
     return nextDt;
   }
 
+  virtual real64 setNextDtBasedOnNewtonIter( real64 const & currentDt ) override
+  {
+    real64 nextDt = SolverBase::setNextDtBasedOnNewtonIter( currentDt );
+    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
+    {
+      real64 const singlePhysicsNextDt =
+        solver->setNextDtBasedOnNewtonIter( currentDt );
+      nextDt = LvArray::math::min( singlePhysicsNextDt, nextDt );
+    } );
+    return nextDt;
+  }
+
   virtual void cleanup( real64 const time_n,
                         integer const cycleNumber,
                         integer const eventCounter,
@@ -369,10 +381,6 @@ class CoupledSolver : public SolverBase
   {
     GEOS_MARK_FUNCTION;
 
-    real64 dtReturn = dt;
-
-    real64 dtReturnTemporary;
-
     Timestamp const meshModificationTimestamp = getMeshModificationTimestamp( domain );
 
     // First call Coupled Solver setup  (important for poromechanics initialization for sequentially coupled)
@@ -397,73 +405,120 @@ class CoupledSolver : public SolverBase
     } );
 
     NonlinearSolverParameters & solverParams = getNonlinearSolverParameters();
-    integer & iter = solverParams.m_numNewtonIterations;
-    iter = 0;
-    bool isConverged = false;
+    integer const maxNumberDtCuts = solverParams.m_maxTimeStepCuts;
+    real64 const dtCutFactor = solverParams.m_timeStepCutFactor;
+    integer & dtAttempt = solverParams.m_numTimeStepAttempts;
 
-    // Reset the states of all solvers if any of them had to restart
-    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
-    {
-      solver->resetStateToBeginningOfStep( domain );
-      solver->getSolverStatistics().initializeTimeStepStatistics();     // initialize counters for subsolvers
-    } );
-    resetStateToBeginningOfStep( domain );
+    bool isConverged = false;
+    // dt may be cut during the course of this step, so we are keeping a local
+    // value to track the achieved dt for this step.
+    real64 stepDt = dt;
 
-    /// Sequential coupling loop
-    while( iter < solverParams.m_maxIterNewton )
+    // outer loop attempts to apply full timestep, and managed the cutting of the timestep if
+    // required.
+    for( dtAttempt = 0; dtAttempt < maxNumberDtCuts; ++dtAttempt )
     {
-      // Increment the solver statistics for reporting purposes
-      // Pass a "0" as argument (0 linear iteration) to skip the output of linear iteration stats at the end
-      m_solverStatistics.logNonlinearIteration( 0 );
+      // TODO configuration loop
 
-      startSequentialIteration( iter, domain );
+      // Reset the states of all solvers if any of them had to restart
+      forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
+      {
+        solver->resetStateToBeginningOfStep( domain );
+        solver->getSolverStatistics().initializeTimeStepStatistics();     // initialize counters for subsolvers
+      } );
+      resetStateToBeginningOfStep( domain );
 
-      // Solve the subproblems nonlinearly
-      forEachArgInTuple( m_solvers, [&]( auto & solver, auto idx )
+      integer & iter = solverParams.m_numNewtonIterations;
+      iter = 0;
+      /// Sequential coupling loop
+      while( iter < solverParams.m_maxIterNewton )
       {
-        GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "  Iteration {:2}: {}", iter+1, solver->getName() ) );
-        dtReturnTemporary = solver->nonlinearImplicitStep( time_n,
-                                                           dtReturn,
+        // Increment the solver statistics for reporting purposes
+        // Pass a "0" as argument (0 linear iteration) to skip the output of linear iteration stats at the end
+        m_solverStatistics.logNonlinearIteration( 0 );
+
+        startSequentialIteration( iter, domain );
+
+        // Solve the subproblems nonlinearly
+        forEachArgInTuple( m_solvers, [&]( auto & solver, auto idx )
+        {
+          GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "  Iteration {:2}: {}", iter + 1, solver->getName() ) );
+          real64 solverDt = solver->nonlinearImplicitStep( time_n,
+                                                           stepDt,
                                                            cycleNumber,
                                                            domain );
 
-        mapSolutionBetweenSolvers( domain, idx() );
+          mapSolutionBetweenSolvers( domain, idx() );
 
-        if( dtReturnTemporary < dtReturn )
+          if( solverDt < stepDt ) // subsolver had to cut the time step
+          {
+            iter = 0; // restart outer loop
+            stepDt = solverDt; // sync time step
+          }
+        } );
+
+        // Check convergence of the outer loop
+        isConverged = checkSequentialConvergence( iter,
+                                                  time_n,
+                                                  stepDt,
+                                                  domain );
+
+        if( isConverged )
         {
-          iter = 0;
-          dtReturn = dtReturnTemporary;
+          // Save Time step statistics for the subsolvers
+          forEachArgInTuple( m_solvers, [&]( auto & solver,
+                                             auto )
+          {
+            solver->getSolverStatistics().saveTimeStepStatistics();
+          } );
+          break;
+        }
+        else
+        {
+          finishSequentialIteration( iter, domain );
         }
-      } );
 
-      // Check convergence of the outer loop
-      isConverged = checkSequentialConvergence( iter,
-                                                time_n,
-                                                dtReturn,
-                                                domain );
+        ++iter;
+      }
 
       if( isConverged )
       {
-        // Save Time step statistics for the subsolvers
-        forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
-        {
-          solver->getSolverStatistics().saveTimeStepStatistics();
-        } );
+        // get out of time loop
         break;
       }
       else
       {
-        finishSequentialIteration( iter, domain );
+        // cut timestep, go back to beginning of step and restart the Newton loop
+        stepDt *= dtCutFactor;
+        GEOS_LOG_LEVEL_RANK_0 ( 1, GEOS_FMT( "New dt = {}", stepDt ) );
+
+        // notify the solver statistics counter that this is a time step cut
+        m_solverStatistics.logTimeStepCut();
+        forEachArgInTuple( m_solvers, [&]( auto & solver,
+                                           auto )
+        {
+          solver->getSolverStatistics().logTimeStepCut();
+        } );
       }
-      // Add convergence check:
-      ++iter;
     }
 
-    GEOS_ERROR_IF( !isConverged, getDataContext() << ": sequentiallyCoupledSolverStep did not converge!" );
+    if( !isConverged )
+    {
+      GEOS_LOG_RANK_0( "Convergence not achieved." );
+
+      if( m_nonlinearSolverParameters.m_allowNonConverged > 0 )
+      {
+        GEOS_LOG_RANK_0( "The accepted solution may be inaccurate." );
+      }
+      else
+      {
+        GEOS_ERROR( "Nonconverged solutions not allowed. Terminating..." );
+      }
+    }
 
-    implicitStepComplete( time_n, dt, domain );
+    implicitStepComplete( time_n, stepDt, domain );
 
-    return dtReturn;
+    return stepDt;
   }
 
   /**

src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp

@@ -335,7 +335,7 @@ void MultiphasePoromechanics< FLOW_SOLVER >::updateState( DomainPartition & doma
     } );
   } );
 
-  GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "        {}: Max phase volume fraction change: {}",
+  GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "        {}: Max phase volume fraction change = {}",
                                       this->getName(), GEOS_FMT( "{:.{}f}", maxDeltaPhaseVolFrac, 4 ) ) );
 }