diff --git a/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ALM_multiphasePoromechanics_curvedFrac_smoke.xml b/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ALM_multiphasePoromechanics_curvedFrac_smoke.xml index b74322d75c0..bd3162df17b 100644 --- a/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ALM_multiphasePoromechanics_curvedFrac_smoke.xml +++ b/inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ALM_multiphasePoromechanics_curvedFrac_smoke.xml @@ -63,6 +63,7 @@ discretization="fluidTPFA" targetRegions="{ Region, Reservoir, Fracture}" temperature="368.15" > + + targetRegions="{ Region, Reservoir, Fracture }" > + + targetRegions="{ Region, Fault}" > + 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(), diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index 8d05dfd6958..cf51f7872a1 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/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,33 +922,21 @@ 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() ); } } @@ -953,6 +944,15 @@ real64 PhysicsSolverBase::nonlinearImplicitStep( real64 const & time_n, 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(); diff --git a/src/coreComponents/physicsSolvers/SolverStatistics.hpp b/src/coreComponents/physicsSolvers/SolverStatistics.hpp index 334f1672c98..c3110f4612d 100644 --- a/src/coreComponents/physicsSolvers/SolverStatistics.hpp +++ b/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 */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp index 356661e434e..95e40473898 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp @@ -383,6 +383,8 @@ TestSet getTestSet() discretization="singlePhaseTPFA" targetRegions="{ reservoir }" > + @@ -619,11 +621,11 @@ TestSet getTestSet() discretization="fluidTPFA" targetRegions="{ reservoir }" temperature="366.483" - useMass="1" - logLevel="1" > + useMass="1" + logLevel="1" > + allowNonConverged="1" /> @@ -891,11 +893,11 @@ TestSet getTestSet() discretization="fluidTPFA" targetRegions="{ reservoir }" temperature="366.483" - useMass="0" - logLevel="1" > + useMass="0" + logLevel="1" > + allowNonConverged="1" /> diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp index c5cddf6acd1..996e9ca150b 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp +++ b/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() ); } } diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index deee39412d2..bb6ccc35163 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -2452,7 +2452,8 @@ Information output from lower logLevels is added with the desired log level - +