CoupledSolver.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDSOLVER_HPP_
#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDSOLVER_HPP_
 
#include "physicsSolvers/PhysicsSolverBase.hpp"
#include "physicsSolvers/LogLevelsInfo.hpp"
 
#include <tuple>
 
namespace geos
{
 
template< typename ... SOLVERS >
class CoupledSolver : public PhysicsSolverBase
{
 
public:
 
  CoupledSolver( const string & name,
                 Group * const parent )
    : PhysicsSolverBase( name, parent )
  {
    forEachArgInTuple( m_solvers, [&]( auto solver, auto idx )
    {
      using SolverType = TYPEOFPTR( solver );
      string const key = SolverType::coupledSolverAttributePrefix() + "SolverName";
      registerWrapper( key, &m_names[idx()] ).
        setRTTypeName( rtTypes::CustomTypes::groupNameRef ).
        setInputFlag( dataRepository::InputFlags::REQUIRED ).
        setDescription( "Name of the " + SolverType::coupledSolverAttributePrefix() + " solver used by the coupled solver" );
    } );
 
    this->getWrapper< string >( PhysicsSolverBase::viewKeyStruct::discretizationString() ).
      setInputFlag( dataRepository::InputFlags::FALSE );
 
    addLogLevel< logInfo::Coupling >();
  }
 
  CoupledSolver( CoupledSolver const & ) = delete;
 
  CoupledSolver( CoupledSolver && ) = default;
 
  CoupledSolver & operator=( CoupledSolver const & ) = delete;
 
  CoupledSolver & operator=( CoupledSolver && ) = delete;
 
 
  void
  setSubSolvers()
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto idx )
    {
      using SolverPtr = TYPEOFREF( solver );
      using SolverType = TYPEOFPTR( SolverPtr {} );
      auto const & solverName = m_names[idx()];
      auto const & solverType = LvArray::system::demangleType< SolverType >();
      solver = this->getParent().template getGroupPointer< SolverType >( solverName );
      GEOS_THROW_IF( solver == nullptr,
                     GEOS_FMT( "{}: Could not find solver '{}' of type {}",
                               getDataContext(),
                               solverName, solverType ),
                     InputError );
      GEOS_LOG_LEVEL_RANK_0( logInfo::Coupling,
                             GEOS_FMT( "{}: found {} solver named {}",
                                       getName(), solver->getCatalogName(), solverName ) );
    } );
  }
 
 
  virtual void
  setupCoupling( DomainPartition const & domain,
                 DofManager & dofManager ) const
  { GEOS_UNUSED_VAR( domain, dofManager ); }
 
  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 )
  { GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); }
 
  void
  setupDofs( DomainPartition const & domain,
             DofManager & dofManager ) const override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->setupDofs( domain, dofManager );
    } );
 
    setupCoupling( domain, dofManager );
  }
 
  virtual void
  implicitStepSetup( real64 const & time_n,
                     real64 const & dt,
                     DomainPartition & domain ) override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->implicitStepSetup( time_n, dt, domain );
    } );
  }
 
  virtual void
  implicitStepComplete( real64 const & time_n,
                        real64 const & dt,
                        DomainPartition & domain ) override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->implicitStepComplete( time_n, dt, domain );
    } );
  }
 
  // general version of assembleSystem function, keep in mind many solvers will override it
  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
  {
 
    // 1. Assemble matrix blocks of each individual solver
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->assembleSystem( time_n, dt, domain, dofManager, localMatrix, localRhs );
    } );
 
    // 2. Assemble coupling blocks
    assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs );
  }
 
  virtual void
  applySystemSolution( DofManager const & dofManager,
                       arrayView1d< real64 const > const & localSolution,
                       real64 const scalingFactor,
                       real64 const dt,
                       DomainPartition & domain ) override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->applySystemSolution( dofManager, localSolution, scalingFactor, dt, domain );
    } );
  }
 
  virtual void
  updateState( DomainPartition & domain ) override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->updateState( domain );
    } );
  }
 
  virtual void
  resetStateToBeginningOfStep( DomainPartition & domain ) override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->resetStateToBeginningOfStep( domain );
    } );
  }
 
  real64
  solverStep( real64 const & time_n,
              real64 const & dt,
              int const cycleNumber,
              DomainPartition & domain ) override final
  {
    GEOS_MARK_FUNCTION;
 
    if( getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::FullyImplicit )
    {
      return fullyCoupledSolverStep( time_n, dt, cycleNumber, domain );
    }
    else if( getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::Sequential )
    {
      return sequentiallyCoupledSolverStep( time_n, dt, cycleNumber, domain );
    }
    else
    {
      GEOS_ERROR( getDataContext() << ": Invalid coupling type option." );
      return 0;
    }
 
  }
 
 
  virtual void
  updateAndWriteConvergenceStep( real64 const & time_n, real64 const & dt,
                                 integer const cycleNumber, integer const iteration ) override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      if( m_writeStatisticsCSV >= 2 )
      {
        solver->updateAndWriteConvergenceStep( time_n, dt, cycleNumber, iteration );
      }
    } );
  }
 
  virtual real64
  calculateResidualNorm( real64 const & time_n,
                         real64 const & dt,
                         DomainPartition const & domain,
                         DofManager const & dofManager,
                         arrayView1d< real64 const > const & localRhs ) override
  {
    real64 norm = 0.0;
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      real64 const singlePhysicsNorm = solver->calculateResidualNorm( time_n, dt, domain, dofManager, localRhs );
      norm += singlePhysicsNorm * singlePhysicsNorm;
    } );
 
    return sqrt( norm );
  }
 
  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
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->applyBoundaryConditions( time_n, dt, domain, dofManager, localMatrix, localRhs );
    } );
  }
 
  virtual bool
  checkSystemSolution( DomainPartition & domain,
                       DofManager const & dofManager,
                       arrayView1d< real64 const > const & localSolution,
                       real64 const scalingFactor ) override
  {
    bool validSolution = true;
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      bool const validSinglePhysicsSolution = solver->checkSystemSolution( domain, dofManager, localSolution, scalingFactor );
      if( !validSinglePhysicsSolution )
      {
        GEOS_LOG_RANK_0( GEOS_FMT( "      {}/{}: Solution check failed. Newton loop terminated.", getName(), solver->getName()) );
      }
      validSolution = validSolution && validSinglePhysicsSolution;
    } );
    return validSolution;
  }
 
  virtual real64
  scalingForSystemSolution( DomainPartition & domain,
                            DofManager const & dofManager,
                            arrayView1d< real64 const > const & localSolution ) override
  {
    real64 scalingFactor = PhysicsSolverBase::scalingForSystemSolution( domain, dofManager, localSolution );
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      real64 const singlePhysicsScalingFactor = solver->scalingForSystemSolution( domain, dofManager, localSolution );
      scalingFactor = LvArray::math::min( scalingFactor, singlePhysicsScalingFactor );
    } );
    return scalingFactor;
  }
 
  virtual real64
  setNextDt( real64 const & currentTime,
             real64 const & currentDt,
             DomainPartition & domain ) override
  {
    real64 nextDt = PhysicsSolverBase::setNextDt( currentTime, currentDt, domain );
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      real64 const singlePhysicsNextDt =
        solver->setNextDt( currentTime, currentDt, domain );
      nextDt = LvArray::math::min( singlePhysicsNextDt, nextDt );
    } );
    return nextDt;
  }
 
  virtual void cleanup( real64 const time_n,
                        integer const cycleNumber,
                        integer const eventCounter,
                        real64 const eventProgress,
                        DomainPartition & domain ) override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain );
    } );
    PhysicsSolverBase::cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain );
  }
 
  virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const override
  {
    bool isConverged = true;
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      isConverged &= solver->checkSequentialSolutionIncrements( domain );
    } );
    return isConverged;
  }
 
  virtual bool updateConfiguration( DomainPartition & domain,
                                    integer const configurationLoopIter ) override
  {
    bool result = true;
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      result &= solver->updateConfiguration( domain, configurationLoopIter );
    } );
    return result;
  }
 
  virtual void outputConfigurationStatistics( DomainPartition const & domain ) const override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->outputConfigurationStatistics( domain );
    } );
  }
 
  virtual void resetConfigurationToBeginningOfStep( DomainPartition & domain ) override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->resetConfigurationToBeginningOfStep( domain );
    } );
  }
 
  virtual bool resetConfigurationToDefault( DomainPartition & domain ) const override
  {
    bool result = true;
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      result &= solver->resetConfigurationToDefault( domain );
    } );
    return result;
  }
 
protected:
 
  virtual real64 fullyCoupledSolverStep( real64 const & time_n,
                                         real64 const & dt,
                                         int const cycleNumber,
                                         DomainPartition & domain )
  {
    return PhysicsSolverBase::solverStep( time_n, dt, cycleNumber, domain );
  }
 
  virtual real64 sequentiallyCoupledSolverStep( real64 const & time_n,
                                                real64 const & dt,
                                                integer const cycleNumber,
                                                DomainPartition & domain )
  {
    GEOS_MARK_FUNCTION;
 
    // Only build the sparsity pattern if the mesh has changed
    Timestamp const meshModificationTimestamp = getMeshModificationTimestamp( domain );
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      if( meshModificationTimestamp > solver->getSystemSetupTimestamp() )
      {
        solver->setupSystem( domain,
                             solver->getDofManager(),
                             solver->getLocalMatrix(),
                             solver->getSystemRhs(),
                             solver->getSystemSolution() );
        solver->setSystemSetupTimestamp( meshModificationTimestamp );
      }
    } );
 
    implicitStepSetup( time_n, dt, domain );
 
    NonlinearSolverParameters & solverParams = getNonlinearSolverParameters();
    integer const maxNumberDtCuts = solverParams.m_maxTimeStepCuts;
    real64 const dtCutFactor = solverParams.m_timeStepCutFactor;
    integer & dtAttempt = solverParams.m_numTimeStepAttempts;
 
    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;
 
    // outer loop attempts to apply full timestep, and managed the cutting of the timestep if
    // required.
    for( dtAttempt = 0; dtAttempt < maxNumberDtCuts; ++dtAttempt )
    {
      // TODO configuration loop
 
      // Reset the states of all solvers if any of them had to restart
      forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
      {
        solver->resetStateToBeginningOfStep( domain );
        solver->getIterationStats().resetCurrentTimeStepStatistics(); // initialize counters for subsolvers
      } );
      resetStateToBeginningOfStep( domain );
 
      integer & iter = solverParams.m_numNewtonIterations;
 
      for( iter = 0; iter < solverParams.m_maxIterNewton; iter++ )
      {
        // Increment the solver statistics for reporting purposes
        getIterationStats().updateNonlinearIteration( 0 );
 
        startSequentialIteration( iter, domain );
 
        // Solve the subproblems nonlinearly
        forEachArgInTuple( m_solvers, [&]( auto & solver, auto idx )
        {
          GEOS_LOG_LEVEL_RANK_0( logInfo::NonlinearSolver,
                                 GEOS_FMT( "  Iteration {:2}: {}", iter + 1, solver->getName() ) );
          real64 solverDt = solver->nonlinearImplicitStep( time_n,
                                                           stepDt,
                                                           cycleNumber,
                                                           domain );
 
          // save fields (e.g. pressure and temperature) after inner solve
          if( solver->getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::Sequential )
          {
            solver->saveSequentialIterationState( domain );
          }
 
          mapSolutionBetweenSolvers( domain, idx() );
 
          if( solverDt < stepDt ) // subsolver had to cut the time step
          {
            iter = 0; // restart outer loop
            stepDt = solverDt; // sync time step
            m_numTimestepsSinceLastDtCut = 0;
          }
        } );
        // Check convergence of the outer loop
        isConverged = checkSequentialConvergence( cycleNumber,
                                                  iter,
                                                  time_n,
                                                  stepDt,
                                                  domain );
 
        if( isConverged )
        {
          // we still want to count current iteration
          ++iter;
          // exit outer loop
          break;
        }
        else
        {
          finishSequentialIteration( iter, domain );
        }
      }
 
      if( isConverged )
      {
        // Save time step statistics for the subsolvers
        forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
        {
          solver->getIterationStats().iterateTimeStepStatistics();
        } );
        // 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( !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, stepDt, domain );
 
    return stepDt;
  }
 
  virtual void mapSolutionBetweenSolvers( DomainPartition & domain,
                                          integer const solverType )
  {
    GEOS_UNUSED_VAR( domain, solverType );
  }
 
  virtual bool checkSequentialConvergence( integer const cycleNumber,
                                           integer const iter,
                                           real64 const & time_n,
                                           real64 const & dt,
                                           DomainPartition & domain )
  {
    NonlinearSolverParameters const & params = getNonlinearSolverParameters();
    bool isConverged = true;
 
    if( params.m_subcyclingOption == 0 )
    {
      GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, "***** Single Pass solver, no subcycling *****" );
    }
    else
    {
      GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, GEOS_FMT( "  Iteration {:2}: outer-loop convergence check", iter + 1 ) );
 
      if( params.sequentialConvergenceCriterion() == NonlinearSolverParameters::SequentialConvergenceCriterion::ResidualNorm )
      {
        real64 residualNorm = 0;
 
        // loop over all the single-physics solvers
        forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
        {
 
          solver->getLocalMatrix().toViewConstSizes().zero();
          solver->getSystemRhs().zero();
          arrayView1d< real64 > const localRhs = solver->getSystemRhs().open();
 
          // for each solver, we have to recompute the residual (and Jacobian, although not necessary)
          solver->assembleSystem( time_n,
                                  dt,
                                  domain,
                                  solver->getDofManager(),
                                  solver->getLocalMatrix().toViewConstSizes(),
                                  localRhs );
          solver->applyBoundaryConditions( time_n,
                                           dt,
                                           domain,
                                           solver->getDofManager(),
                                           solver->getLocalMatrix().toViewConstSizes(),
                                           localRhs );
          solver->getSystemRhs().close();
 
          // once this is done, we recompute the single-physics residual
          real64 const singlePhysicsNorm =
            solver->calculateResidualNorm( time_n,
                                           dt,
                                           domain,
                                           solver->getDofManager(),
                                           solver->getSystemRhs().values() );
          residualNorm += singlePhysicsNorm * singlePhysicsNorm;
        } );
 
        // finally, we perform the convergence check on the multiphysics residual
        residualNorm = sqrt( residualNorm );
        GEOS_LOG_LEVEL_RANK_0( logInfo::ResidualNorm,
                               GEOS_FMT( "        ( R ) = ( {:4.2e} )", residualNorm ) );
        getConvergenceStats().setResidualValue( "R", residualNorm );
        updateAndWriteConvergenceStep( time_n, dt, cycleNumber, iter );
 
        isConverged = ( residualNorm < params.m_newtonTol );
 
      }
      else if( params.sequentialConvergenceCriterion() == NonlinearSolverParameters::SequentialConvergenceCriterion::NumberOfNonlinearIterations )
      {
        // TODO also make recursive?
        forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
        {
          NonlinearSolverParameters const & singlePhysicsParams = solver->getNonlinearSolverParameters();
          if( singlePhysicsParams.m_numNewtonIterations > singlePhysicsParams.m_minIterNewton )
          {
            isConverged = false;
          }
        } );
      }
      else if( params.sequentialConvergenceCriterion() == NonlinearSolverParameters::SequentialConvergenceCriterion::SolutionIncrements )
      {
        isConverged = checkSequentialSolutionIncrements( domain );
      }
      else
      {
        GEOS_ERROR( getDataContext() << ": Invalid sequential convergence criterion." );
      }
 
      if( isConverged )
      {
        GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence,
                               GEOS_FMT( "***** The iterative coupling has converged in {} iteration(s) *****", iter + 1 ) );
      }
    }
    return isConverged;
  }
 
  virtual void
  postInputInitialization() override
  {
    setSubSolvers();
 
    PhysicsSolverBase::postInputInitialization();
 
    bool const isSequential = getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::Sequential;
    bool const usesLineSearch = getNonlinearSolverParameters().m_lineSearchAction != NonlinearSolverParameters::LineSearchAction::None;
    GEOS_THROW_IF( isSequential && usesLineSearch,
                   GEOS_FMT( "{}: line search is not supported by the coupled solver when {} is set to `{}`. Please set {} to `{}` to remove this error",
                             getNonlinearSolverParameters().getWrapperDataContext( NonlinearSolverParameters::viewKeysStruct::couplingTypeString() ),
                             NonlinearSolverParameters::viewKeysStruct::couplingTypeString(),
                             EnumStrings< NonlinearSolverParameters::CouplingType >::toString( NonlinearSolverParameters::CouplingType::Sequential ),
                             NonlinearSolverParameters::viewKeysStruct::lineSearchActionString(),
                             EnumStrings< NonlinearSolverParameters::LineSearchAction >::toString( NonlinearSolverParameters::LineSearchAction::None ) ),
                   InputError );
 
    if( !isSequential )
    {
      synchronizeNonlinearSolverParameters();
    }
 
    if( m_nonlinearSolverParameters.m_nonlinearAccelerationType != NonlinearSolverParameters::NonlinearAccelerationType::None )
    {
      validateNonlinearAcceleration();
    }
  }
 
  virtual void validateNonlinearAcceleration()
  {
    GEOS_THROW ( GEOS_FMT( "{}: Nonlinear acceleration {} is not supported by {} solver '{}'",
                           getWrapperDataContext( NonlinearSolverParameters::viewKeysStruct::nonlinearAccelerationTypeString() ),
                           EnumStrings< NonlinearSolverParameters::NonlinearAccelerationType >::toString( m_nonlinearSolverParameters.m_nonlinearAccelerationType ),
                           getCatalogName(), getName()),
                 InputError );
  }
 
  virtual void
  synchronizeNonlinearSolverParameters() override
  {
    forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
    {
      solver->getNonlinearSolverParameters() = getNonlinearSolverParameters();
    } );
  }
 
  virtual void startSequentialIteration( integer const & iter,
                                         DomainPartition & domain )
  {
    GEOS_UNUSED_VAR( iter, domain );
  }
 
  virtual void finishSequentialIteration( integer const & iter,
                                          DomainPartition & domain )
  {
    GEOS_UNUSED_VAR( iter, domain );
  }
 
protected:
 
  std::tuple< SOLVERS *... > m_solvers;
 
  std::array< string, sizeof...( SOLVERS ) > m_names;
};
 
} /* namespace geos */
 
#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_COUPLEDSOLVER_HPP_ */