From 39e75de52fbb5b36d862ddb7adc95c16250771e1 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 21 Dec 2023 00:15:57 -0600 Subject: [PATCH] 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 --- .../physicsSolvers/SolverBase.cpp | 4 +- .../physicsSolvers/SolverBase.hpp | 2 +- .../fluidFlow/CompositionalMultiphaseBase.cpp | 5 +- .../multiphysics/CoupledSolver.hpp | 147 ++++++++++++------ .../multiphysics/MultiphasePoromechanics.cpp | 2 +- 5 files changed, 108 insertions(+), 52 deletions(-) diff --git a/src/coreComponents/physicsSolvers/SolverBase.cpp b/src/coreComponents/physicsSolvers/SolverBase.cpp index 76b9ae227b2..5ea5eea1754 100644 --- a/src/coreComponents/physicsSolvers/SolverBase.cpp +++ b/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..." ); } } diff --git a/src/coreComponents/physicsSolvers/SolverBase.hpp b/src/coreComponents/physicsSolvers/SolverBase.hpp index 8129e754402..f6174cd9b57 100644 --- a/src/coreComponents/physicsSolvers/SolverBase.hpp +++ b/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 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index f4897c4b9b6..c2edc85b3c8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/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] ) ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp index 0ae7961357f..ad49c26ea2c 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp +++ b/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; } /** diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp index 5e2f1690581..899254b1c64 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp +++ b/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 ) ) ); }