bugfix: allowNonConverged warning instead of aborting

inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ALM_multiphasePoromechanics_curvedFrac_smoke.xml

@@ -63,6 +63,7 @@
       discretization="fluidTPFA"
       targetRegions="{ Region, Reservoir, Fracture}" 
       temperature="368.15" >
+      <!-- 'allowNonConverged' is a workaround, DO NOT USE IN PRODUCTION WITHOUT CONSIDERATION. -->
       <NonlinearSolverParameters
         newtonTol="1.0e-2"
         newtonMaxIter="8"

inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ALM_singlephasePoromechanics_curvedFrac_smoke.xml

@@ -63,7 +63,8 @@
       name="flowSolver"
       logLevel="1"
       discretization="fluidTPFA"
-      targetRegions="{ Region, Reservoir, Fracture }" > 
+      targetRegions="{ Region, Reservoir, Fracture }" >
+      <!-- 'allowNonConverged' is a workaround, DO NOT USE IN PRODUCTION WITHOUT CONSIDERATION. -->
       <NonlinearSolverParameters
         newtonTol="1.0e-2"
         newtonMaxIter="8"

inputFiles/poromechanicsFractures/Contact/LagrangianMultipliers/PoroElastic_dfm_PEBICrack_smoke.xml

@@ -37,7 +37,8 @@
       name="flowSolver"
       logLevel="1"
       discretization="singlePhaseTPFA"
-      targetRegions="{ Region, Fault}" > 
+      targetRegions="{ Region, Fault}" >
+      <!-- 'allowNonConverged' is a workaround, DO NOT USE IN PRODUCTION WITHOUT CONSIDERATION. -->
       <NonlinearSolverParameters
         newtonTol="1.0e-2"
         newtonMaxIter="8"

src/coreComponents/physicsSolvers/NonlinearSolverParameters.cpp

@@ -104,8 +104,9 @@ NonlinearSolverParameters::NonlinearSolverParameters( string const & name,
   registerWrapper( viewKeysStruct::allowNonConvergedString(), &m_allowNonConverged ).
     setApplyDefaultValue( 0 ).
     setInputFlag( InputFlags::OPTIONAL ).
-    setDescription( "Allow non-converged solution to be accepted. "
-                    "(i.e. exit from the Newton loop without achieving the desired tolerance)" );
+    setDescription( "When enabled, the simulation can continue even with a non-converged solution (residuals exceeding the tolerance during the "
+                    "newton loops) at the deepest level of sub-timesteps (maxSubSteps).\n"
+                    "When disabled, the simulation will stop with an error when no converged solution has been found at the last sub-step level." );
 
   registerWrapper( viewKeysStruct::timeStepDecreaseIterLimString(), &m_timeStepDecreaseIterLimit ).
     setApplyDefaultValue( 0.7 ).
@@ -214,6 +215,14 @@ NonlinearSolverParameters::NonlinearSolverParameters( string const & name,
 
 void NonlinearSolverParameters::postInputInitialization()
 {
+  if( m_allowNonConverged > 0 )
+  {
+    GEOS_WARNING( GEOS_FMT( "Solver '{}' has '{}' enabled.\nThe simulation will not stop on non-converged solutions, "
+                            "result may contain inaccurate solutions that VIOLATE CONVERGENCE TOLERANCES.",
+                            getParent().getName(), viewKeysStruct::allowNonConvergedString() ),
+                  getWrapperDataContext( viewKeysStruct::allowNonConvergedString() ) );
+  }
+
   GEOS_ERROR_IF_LE_MSG( m_timeStepDecreaseIterLimit, m_timeStepIncreaseIterLimit,
                         GEOS_FMT( "{} Value should be smaller than {}",
                                   viewKeysStruct::timeStepIncreaseIterLimString(),

src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp

@@ -217,17 +217,17 @@ void PhysicsSolverBase::registerDataOnMesh( Group & meshBodies )
   forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &,
                                                     MeshLevel & mesh,
                                                     string_array const & regionNames )
-  {
-    ElementRegionManager & elemManager = mesh.getElemManager();
-    elemManager.forElementSubRegions< ElementSubRegionBase >( regionNames,
-                                                              [&]( localIndex const,
-                                                                   ElementSubRegionBase & subRegion )
     {
-      setConstitutiveNamesCallSuper( subRegion );
-      setConstitutiveNames( subRegion );
-    } );
+      ElementRegionManager & elemManager = mesh.getElemManager();
+      elemManager.forElementSubRegions< ElementSubRegionBase >( regionNames,
+                                                                [&]( localIndex const,
+                                                                     ElementSubRegionBase & subRegion )
+      {
+        setConstitutiveNamesCallSuper( subRegion );
+        setConstitutiveNames( subRegion );
+      } );
 
-  } );
+    } );
 
 }
 
@@ -358,10 +358,7 @@ bool PhysicsSolverBase::execute( real64 const time_n,
                                         GEOS_FMT( "{}: shortening time step to {} to match remaining time", getName(), nextDt ),
                                         m_nonlinearSolverParameters );
       }
-    }
 
-    if( dtRemaining > 0.0 )
-    {
       GEOS_LOG_LEVEL_RANK_0( logInfo::TimeStep,
                              GEOS_FMT( "{}: sub-step = {}, accepted dt = {}, next dt = {}, remaining dt = {}",
                                        getName(), subStep, dtAccepted, nextDt, dtRemaining ) );
@@ -827,10 +824,10 @@ real64 PhysicsSolverBase::eisenstatWalker( real64 const newNewtonNorm,
   return krylovTol;
 }
 
-real64 PhysicsSolverBase::nonlinearImplicitStep( real64 const & time_n,
-                                                 real64 const & dt,
-                                                 integer const cycleNumber,
-                                                 DomainPartition & domain )
+StepResult PhysicsSolverBase::nonlinearImplicitStep( real64 const & time_n,
+                                                     real64 const & dt,
+                                                     integer const cycleNumber,
+                                                     DomainPartition & domain )
 {
   GEOS_MARK_FUNCTION;
   // dt may be cut during the course of this step, so we are keeping a local
@@ -854,13 +851,19 @@ real64 PhysicsSolverBase::nonlinearImplicitStep( real64 const & time_n,
   // required.
   for( dtAttempt = 0; dtAttempt < maxNumberDtCuts; ++dtAttempt )
   {
-    // reset the solver state, since we are restarting the time step
     if( dtAttempt > 0 )
     {
+      // reset the solver state, since we are restarting the time step
       Timer timer( m_timers.get_inserted( "reset state" ) );
-
       resetStateToBeginningOfStep( domain );
       resetConfigurationToBeginningOfStep( domain );
+
+      // previous attempt failed, we try cuttting the timestep
+      stepDt *= dtCutFactor;
+      m_numTimestepsSinceLastDtCut = 0;
+      GEOS_LOG_LEVEL_RANK_0 ( logInfo::TimeStep, GEOS_FMT( "New dt = {}", stepDt ) );
+      getIterationStats().updateTimeStepCut();
+      getIterationStats().writeIterationStatsToTable();
     }
 
     // it's the simplest configuration that can be attempted whenever Newton's fails as a last resource.
@@ -919,40 +922,37 @@ real64 PhysicsSolverBase::nonlinearImplicitStep( real64 const & time_n,
 
     if( isConfigurationLoopConverged )
     {
-      // get out of outer loop
-      break;
-    }
-    else
-    {
-      // cut timestep, go back to beginning of step and restart the Newton loop
-      stepDt *= dtCutFactor;
-      m_numTimestepsSinceLastDtCut = 0;
-      GEOS_LOG_LEVEL_RANK_0 ( logInfo::TimeStep, GEOS_FMT( "New dt = {}", stepDt ) );
-
-      // notify the solver statistics counter that this is a time step cut
-      getIterationStats().updateTimeStepCut();
-      getIterationStats().writeIterationStatsToTable();
+      break; // we converged, get out of outer loop
     }
   } // end of outer loop (dt chopping strategy)
+  // if not converged & no timestep cut are permitted, we exit the loop
 
   if( !isConfigurationLoopConverged )
   {
-    GEOS_LOG_RANK_0( "Convergence not achieved." );
-
     if( allowNonConverged )
     {
-      GEOS_LOG_RANK_0( "The accepted solution may be inaccurate." );
+      if ( MpiWrapper::commRank() == 0 )
+        recordNonConvergedTimestep( time_n, cycleNumber );
     }
     else
     {
-      GEOS_ERROR( "Nonconverged solutions not allowed. Terminating...", getDataContext()  );
+      GEOS_ERROR( "Convergence not achieved. Terminating...", getDataContext() );
     }
   }
 
   // return the achieved timestep
   return stepDt;
 }
 
+void PhysicsSolverBase::recordNonConvergedTimestep( globalIndex const cycleNumber,
+                                                    real64 const time_n )
+{
+  GEOS_WARNING( dtRemaining > 0.0 && MpiWrapper::commRank() == 0,
+                "Maximum allowed number of sub-steps reached.\nNon-converged solutions are allowed, SIMULATION WILL CONTINUE WITH INACURATE RESULTS.",
+                getDataContext(), getWrapperDataContext( NonlinearSolverParameters::viewKeysStruct::allowNonConvergedString()) );
+  // TODO record non-converged timestep, cycle, residual...
+}
+
 bool PhysicsSolverBase::solveNonlinearSystem( real64 const & time_n,
                                               real64 const & stepDt,
                                               integer const cycleNumber,
@@ -1619,6 +1619,33 @@ void PhysicsSolverBase::cleanup( real64 const GEOS_UNUSED_PARAM( time_n ),
   getIterationStats().closeFile();
   getConvergenceStats().closeFile();
 
+  if( m_nonlinearSolverParameters.m_allowNonConverged && m_nonConvergedCycleNumbers.size() > 0 )
+  { // TODO : relocate where non-convergence data is stored
+    globalIndex const numNonConverged = m_nonConvergedCycleNumbers.size();
+
+    TableLayout nonConvLayout( GEOS_FMT( "{} NON-CONVERGED TIMESTEPS", getName()),
+                               { "Cycle", "Time (s)", "Residual" } );
+
+    TableTextFormatter const nonConvFormatter( nonConvLayout );
+
+    TableData nonConvData;
+    for( globalIndex i = 0; i < numNonConverged; ++i )
+    {
+      real64 const timeVal = i < m_nonConvergedTimes.size() ? m_nonConvergedTimes[i] : 0.0;
+      real64 const residVal = i < m_nonConvergedResiduals.size() ? m_nonConvergedResiduals[i] : 0.0;
+      nonConvData.addRow( m_nonConvergedCycleNumbers[i], timeVal, residVal );
+    }
+
+    std::string const tableStr = nonConvFormatter.toString( nonConvData );
+
+    GEOS_WARNING( GEOS_FMT( "WARNING: Simulation ended with non-converged timesteps:\n"
+                            "- Results contain NON-PHYSICAL values and are likely invalid for use.\n"
+                            "- The simulation has not been aborted because 'allowNonConverged' was active on this solver.\n"
+                            "{}\n",
+                            tableStr ),
+                           );
+  }
+
   for( auto & timer : m_timers )
   {
     real64 const time = std::chrono::duration< double >( timer.second ).count();

src/coreComponents/physicsSolvers/SolverStatistics.hpp

@@ -152,6 +152,13 @@ class IterationsStatistics : public dataRepository::Group
   void setFilename( string_view filename )
   { m_iterationsFilename = filename; }
 
+  /**
+   * @brief Set the table title in the log.
+   * @param name The table as a string_view.
+   */
+  void setTableName( string_view name )
+  { m_tableIterationName = name; }
+
 
   /**
    * @return  A const string iteration filename
@@ -166,13 +173,6 @@ class IterationsStatistics : public dataRepository::Group
   { return m_iterationsFilename;  }
 
 
-  /**
-   * @brief Set the filename output file.
-   * @param name The filename as a string_view.
-   */
-  void setTableName( string_view name )
-  { m_tableIterationName = name; }
-
   /**
    * @brief Output the statistics to the console in table format
    */

src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp

@@ -383,6 +383,8 @@ TestSet getTestSet()
                     discretization="singlePhaseTPFA"
                     targetRegions="{ reservoir }" >
 
+      <!-- 'allowNonConverged' is a workaround to speedup the test on some configurations,
+           DO NOT USE IN PRODUCTION WITHOUT CONSIDERATION. -->
       <NonlinearSolverParameters newtonMaxIter="40"
                                  allowNonConverged="1" />
       <LinearSolverParameters directParallel="0" />
@@ -619,11 +621,11 @@ TestSet getTestSet()
                                 discretization="fluidTPFA"
                                 targetRegions="{ reservoir }"
                                 temperature="366.483"
-                                useMass="1"
-                                logLevel="1" >
+                                 useMass="1"
+                                 logLevel="1" >
       <NonlinearSolverParameters newtonMaxIter="8"
                                  maxTimeStepCuts="8"
-                                 allowNonConverged="1" />
+                                 allowNonConverged="1" /> <!-- workaround for fracture contact convergence in smoke test -->
       <LinearSolverParameters directParallel="0" />
     </CompositionalMultiphaseFVM>
   </Solvers>
@@ -891,11 +893,11 @@ TestSet getTestSet()
                                 discretization="fluidTPFA"
                                 targetRegions="{ reservoir }"
                                 temperature="366.483"
-                                useMass="0"
-                                logLevel="1" >
+                                 useMass="0"
+                                 logLevel="1" >
       <NonlinearSolverParameters newtonMaxIter="8"
                                  maxTimeStepCuts="8"
-                                 allowNonConverged="1" />
+                                 allowNonConverged="1" /> <!-- workaround for fracture contact convergence in smoke test -->
       <LinearSolverParameters directParallel="0" />
     </CompositionalMultiphaseFVM>
   </Solvers>

src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp

@@ -505,6 +505,20 @@ class CoupledSolver : public PhysicsSolverBase
     // required.
     for( dtAttempt = 0; dtAttempt < maxNumberDtCuts; ++dtAttempt )
     {
+      if (dtAttempt > 0) {
+        // cut timestep, go back to beginning of step and restart the Newton loop
+        stepDt *= dtCutFactor;
+        m_numTimestepsSinceLastDtCut = 0;
+        GEOS_LOG_LEVEL_RANK_0( logInfo::TimeStep, GEOS_FMT( "New dt = {}", stepDt ) );
+
+        // notify the solver statistics counter that this is a time step cut
+        getIterationStats().updateTimeStepCut();
+        forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
+        {
+          solver->getIterationStats().updateTimeStepCut();
+        } );
+      }
+
       // TODO configuration loop
 
       // Reset the states of all solvers if any of them had to restart
@@ -580,33 +594,18 @@ class CoupledSolver : public PhysicsSolverBase
         // get out of the time loop
         break;
       }
-      else
-      {
-        // cut timestep, go back to beginning of step and restart the Newton loop
-        stepDt *= dtCutFactor;
-        m_numTimestepsSinceLastDtCut = 0;
-        GEOS_LOG_LEVEL_RANK_0( logInfo::TimeStep, GEOS_FMT( "New dt = {}", stepDt ) );
-
-        // notify the solver statistics counter that this is a time step cut
-        getIterationStats().updateTimeStepCut();
-        forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
-        {
-          solver->getIterationStats().updateTimeStepCut();
-        } );
-      }
     }
+    // if still not converged & no timestep cut are permitted, we exit the loop
 
-    if( !isConverged )
+    if( !isConfigurationLoopConverged )
     {
-      GEOS_LOG_RANK_0( "Convergence not achieved." );
-
-      if( m_nonlinearSolverParameters.m_allowNonConverged > 0 )
+      if( allowNonConverged )
       {
-        GEOS_LOG_RANK_0( "The accepted solution may be inaccurate." );
+        GEOS_LOG_RANK_0( "Convergence not achieved. The accepted solution residuals are out of tolerance bounds." );
       }
       else
       {
-        GEOS_ERROR( "Nonconverged solutions not allowed. Terminating...", getDataContext() );
+        GEOS_ERROR( "Convergence not achieved. Terminating...", getDataContext() );
       }
     }
 

src/coreComponents/schema/schema.xsd

@@ -2452,7 +2452,8 @@ Information output from lower logLevels is added with the desired log level
 		</xsd:restriction>
 	</xsd:simpleType>
 	<xsd:complexType name="NonlinearSolverParametersType">
-		<!--allowNonConverged => Allow non-converged solution to be accepted. (i.e. exit from the Newton loop without achieving the desired tolerance)-->
+		<!--allowNonConverged => When enabled, the simulation can continue even with a non-converged solution (residuals exceeding the tolerance during the newton loops) at the deepest level of sub-timesteps (maxSubSteps).
+When disabled, the simulation will stop with an error when no converged solution has been found at the last sub-step level.-->
 		<xsd:attribute name="allowNonConverged" type="integer" default="0" />
 		<!--configurationTolerance => Configuration tolerance-->
 		<xsd:attribute name="configurationTolerance" type="real64" default="0" />