PhysicsSolverBase.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_PHYSICSSOLVERBASE_HPP_
#define GEOS_PHYSICSSOLVERS_PHYSICSSOLVERBASE_HPP_
 
#include "codingUtilities/traits.hpp"
#include "common/DataTypes.hpp"
#include "common/format/LogPart.hpp"
#include "dataRepository/ExecutableGroup.hpp"
#include "dataRepository/RestartFlags.hpp"
#include "linearAlgebra/interfaces/InterfaceTypes.hpp"
#include "linearAlgebra/utilities/LinearSolverResult.hpp"
#include "linearAlgebra/DofManager.hpp"
#include "mesh/MeshBody.hpp"
#include "physicsSolvers/NonlinearSolverParameters.hpp"
#include "physicsSolvers/LinearSolverParameters.hpp"
#include "physicsSolvers/SolverStatistics.hpp"
 
#include <limits>
 
namespace geos
{
 
class DomainPartition;
 
class PhysicsSolverBase : public ExecutableGroup
{
public:
 
  enum class StatsOutputType : integer
  {
    none, iteration, convergence, all
  };
 
  explicit PhysicsSolverBase( string const & name,
                              Group * const parent );
 
  PhysicsSolverBase( PhysicsSolverBase && ) = default;
 
  virtual ~PhysicsSolverBase() override;
 
  PhysicsSolverBase() = delete;
 
  PhysicsSolverBase( PhysicsSolverBase const & ) = delete;
 
  PhysicsSolverBase & operator=( PhysicsSolverBase const & ) = delete;
 
  PhysicsSolverBase & operator=( PhysicsSolverBase && ) = delete;
 
  virtual string getCatalogName() const = 0;
 
 
  virtual void registerDataOnMesh( Group & MeshBodies ) override;
 
  virtual void initialize_postMeshGeneration() override;
 
  void generateMeshTargetsFromTargetRegions( Group const & meshBodies );
 
  virtual void cleanup( real64 const time_n,
                        integer const cycleNumber,
                        integer const eventCounter,
                        real64 const eventProgress,
                        DomainPartition & domain ) override;
 
  virtual bool execute( real64 const time_n,
                        real64 const dt,
                        integer const cycleNumber,
                        integer const eventCounter,
                        real64 const eventProgress,
                        DomainPartition & domain ) override;
 
  ParallelMatrix & getSystemMatrix() { return m_matrix; }
 
  ParallelMatrix const & getSystemMatrix() const { return m_matrix; }
 
  ParallelVector & getSystemRhs() { return m_rhs; }
 
  ParallelVector const & getSystemRhs() const { return m_rhs; }
 
  ParallelVector & getSystemSolution() { return m_solution; }
 
  ParallelVector const & getSystemSolution() const { return m_solution; }
 
  DofManager & getDofManager() { return m_dofManager; }
 
  DofManager const & getDofManager() const { return m_dofManager; }
 
  CRSMatrix< real64, globalIndex > & getLocalMatrix() { return m_localMatrix; }
 
  CRSMatrixView< real64 const, globalIndex const > getLocalMatrix() const { return m_localMatrix.toViewConst(); }
 
  virtual real64 solverStep( real64 const & time_n,
                             real64 const & dt,
                             integer const cycleNumber,
                             DomainPartition & domain );
 
  virtual real64 setNextDt( real64 const & currentTime,
                            real64 const & currentDt,
                            DomainPartition & domain );
 
  virtual real64 setNextDtBasedOnIterNumber( real64 const & currentDt );
 
  virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt,
                                              DomainPartition & domain );
 
  virtual real64 explicitStep( real64 const & time_n,
                               real64 const & dt,
                               integer const cycleNumber,
                               DomainPartition & domain );
 
  virtual real64 nonlinearImplicitStep( real64 const & time_n,
                                        real64 const & dt,
                                        integer const cycleNumber,
                                        DomainPartition & domain );
 
  virtual bool
  lineSearch( real64 const & time_n,
              real64 const & dt,
              integer const cycleNumber,
              integer const newtonIter,
              DomainPartition & domain,
              DofManager const & dofManager,
              CRSMatrixView< real64, globalIndex const > const & localMatrix,
              ParallelVector & rhs,
              ParallelVector & solution,
              real64 const scaleFactor,
              real64 & lastResidual );
 
  virtual bool
  lineSearchWithParabolicInterpolation ( real64 const & time_n,
                                         real64 const & dt,
                                         integer const cycleNumber,
                                         integer const newtonIter,
                                         DomainPartition & domain,
                                         DofManager const & dofManager,
                                         CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                         ParallelVector & rhs,
                                         ParallelVector & solution,
                                         real64 const scaleFactor,
                                         real64 & lastResidual,
                                         real64 & residualNormT );
 
  virtual real64 linearImplicitStep( real64 const & time_n,
                                     real64 const & dt,
                                     integer const cycleNumber,
                                     DomainPartition & domain );
 
  virtual void
  implicitStepSetup( real64 const & time_n,
                     real64 const & dt,
                     DomainPartition & domain );
 
  virtual void
  setupDofs( DomainPartition const & domain,
             DofManager & dofManager ) const;
 
  virtual void
  setupSystem( DomainPartition & domain,
               DofManager & dofManager,
               CRSMatrix< real64, globalIndex > & localMatrix,
               ParallelVector & rhs,
               ParallelVector & solution,
               bool const setSparsity = true );
 
  virtual void
  setSparsityPattern( DomainPartition & domain,
                      DofManager & dofManager,
                      CRSMatrix< real64, globalIndex > & localMatrix,
                      SparsityPattern< globalIndex > & pattern );
 
  virtual std::unique_ptr< PreconditionerBase< LAInterface > >
  createPreconditioner( DomainPartition & domain ) const;
 
  virtual void
  assembleSystem( real64 const time,
                  real64 const dt,
                  DomainPartition & domain,
                  DofManager const & dofManager,
                  CRSMatrixView< real64, globalIndex const > const & localMatrix,
                  arrayView1d< real64 > const & localRhs );
 
  virtual void
  applyBoundaryConditions( real64 const time,
                           real64 const dt,
                           DomainPartition & domain,
                           DofManager const & dofManager,
                           CRSMatrixView< real64, globalIndex const > const & localMatrix,
                           arrayView1d< real64 > const & localRhs );
 
  void
  debugOutputSystem( real64 const & time,
                     integer const cycleNumber,
                     integer const nonlinearIteration,
                     ParallelMatrix const & matrix,
                     ParallelVector const & rhs ) const;
 
  void
  debugOutputSolution( real64 const & time,
                       integer const cycleNumber,
                       integer const nonlinearIteration,
                       ParallelVector const & solution ) const;
 
  virtual void
  updateAndWriteConvergenceStep( real64 const & time_n,
                                 real64 const & dt,
                                 integer const cycleNumber,
                                 integer const iteration );
 
  virtual real64
  calculateResidualNorm( real64 const & time,
                         real64 const & dt,
                         DomainPartition const & domain,
                         DofManager const & dofManager,
                         arrayView1d< real64 const > const & localRhs );
 
  virtual void
  solveLinearSystem( DofManager const & dofManager,
                     ParallelMatrix & matrix,
                     ParallelVector & rhs,
                     ParallelVector & solution );
 
  virtual bool
  checkSystemSolution( DomainPartition & domain,
                       DofManager const & dofManager,
                       arrayView1d< real64 const > const & localSolution,
                       real64 const scalingFactor );
 
  virtual real64
  scalingForSystemSolution( DomainPartition & domain,
                            DofManager const & dofManager,
                            arrayView1d< real64 const > const & localSolution );
 
  virtual void
  applySystemSolution( DofManager const & dofManager,
                       arrayView1d< real64 const > const & localSolution,
                       real64 const scalingFactor,
                       real64 const dt,
                       DomainPartition & domain );
 
  virtual bool updateConfiguration( DomainPartition & domain,
                                    integer configurationLoopIter );
 
  virtual void outputConfigurationStatistics( DomainPartition const & domain ) const;
 
  virtual void resetConfigurationToBeginningOfStep( DomainPartition & domain );
 
  virtual bool resetConfigurationToDefault( DomainPartition & domain ) const;
 
 
  virtual void updateState( DomainPartition & domain );
 
  virtual void
  resetStateToBeginningOfStep( DomainPartition & domain );
 
  virtual void
  implicitStepComplete( real64 const & time,
                        real64 const & dt,
                        DomainPartition & domain );
 
 
  virtual real64 getTimestepRequest( real64 const ) override
  {return m_nextDt;};
  real64 getTimestepRequest()
  {return m_nextDt;};
 
  virtual Group * createChild( string const & childKey, string const & childName ) override;
 
  using CatalogInterface = dataRepository::CatalogInterface< PhysicsSolverBase, string const &, Group * const >;
 
  static CatalogInterface::CatalogType & getCatalog();
 
  struct viewKeyStruct
  {
    static constexpr char const * cflFactorString() { return "cflFactor"; }
 
    static constexpr char const * initialDtString() { return "initialDt"; }
 
    static constexpr char const * minDtIncreaseIntervalString() { return "minDtIncreaseInterval"; }
 
    static constexpr char const * discretizationString() { return "discretization"; }
 
    static constexpr char const * targetRegionsString() { return "targetRegions"; }
 
    static constexpr char const * writeLinearSystemString() { return "writeLinearSystem"; }
 
    static constexpr char const * usePhysicsScalingString() { return "usePhysicsScaling"; }
 
    static constexpr char const * allowNonConvergedLinearSolverSolutionString() { return "allowNonConvergedLinearSolverSolution"; }
 
    static constexpr char const * writeStatisticsCSVString() { return "writeStatistics"; }
 
    static constexpr char const * numTimestepsSinceLastDtCutString() { return "numTimestepsSinceLastDtCut"; }
  };
 
  struct groupKeyStruct
  {
    static constexpr char const * linearSolverParametersString() { return "LinearSolverParameters"; }
 
    static constexpr char const * nonlinearSolverParametersString() { return "NonlinearSolverParameters"; }
 
    static constexpr char const * solverStatisticsString() { return "SolverStatistics"; }
  };
 
  Timestamp getSystemSetupTimestamp() const { return m_systemSetupTimestamp; }
 
  Timestamp getMeshModificationTimestamp( DomainPartition & domain ) const;
 
  void setSystemSetupTimestamp( Timestamp timestamp );
 
  R1Tensor const gravityVector() const;
 
  virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const;
 
  virtual void saveSequentialIterationState( DomainPartition & domain );
 
  LinearSolverParameters & getLinearSolverParameters()
  {
    return m_linearSolverParameters.get();
  }
 
  LinearSolverParameters const & getLinearSolverParameters() const
  {
    return m_linearSolverParameters.get();
  }
 
  NonlinearSolverParameters & getNonlinearSolverParameters()
  {
    return m_nonlinearSolverParameters;
  }
 
  NonlinearSolverParameters const & getNonlinearSolverParameters() const
  {
    return m_nonlinearSolverParameters;
  }
 
  virtual void
  synchronizeNonlinearSolverParameters()
  { /* empty here, overriden in CoupledSolver */ }
 
  localIndex targetRegionIndex( string const & regionName ) const;
 
  string_array const & getTargetRegionNames() const {return m_targetRegionNames;}
 
 
  template< typename LAMBDA >
  void forDiscretizationOnMeshTargets( Group const & meshBodies, LAMBDA && lambda ) const
  {
    for( auto const & target: m_meshTargets )
    {
      string const meshBodyName = target.first.first;
      string const meshLevelName = target.first.second;
      string_array const & regionNames = target.second;
      MeshBody const & meshBody = meshBodies.getGroup< MeshBody >( meshBodyName );
 
      MeshLevel const * meshLevelPtr = meshBody.getMeshLevels().getGroupPointer< MeshLevel >( meshLevelName );
      if( meshLevelPtr==nullptr )
      {
        meshLevelPtr = meshBody.getMeshLevels().getGroupPointer< MeshLevel >( MeshBody::groupStructKeys::baseDiscretizationString() );
      }
      lambda( meshBodyName, *meshLevelPtr, regionNames );
    }
  }
 
  template< typename LAMBDA >
  void forDiscretizationOnMeshTargets( Group & meshBodies, LAMBDA && lambda ) const
  {
    for( auto const & target: m_meshTargets )
    {
      string const meshBodyName = target.first.first;
      string const meshLevelName = target.first.second;
      string_array const & regionNames = target.second;
      MeshBody & meshBody = meshBodies.getGroup< MeshBody >( meshBodyName );
 
      MeshLevel * meshLevelPtr = meshBody.getMeshLevels().getGroupPointer< MeshLevel >( meshLevelName );
      if( meshLevelPtr==nullptr )
      {
        meshLevelPtr = meshBody.getMeshLevels().getGroupPointer< MeshLevel >( MeshBody::groupStructKeys::baseDiscretizationString() );
      }
      lambda( meshBodyName, *meshLevelPtr, regionNames );
    }
  }
 
  string getDiscretizationName() const {return m_discretizationName;}
 
  virtual bool registerCallback( void * func, const std::type_info & funcType ) final override;
 
  IterationsStatistics & getIterationStats()
  {
    return m_solverStatistics.m_iterationsStats;
  }
  IterationsStatistics const & getIterationStats() const
  {
    return m_solverStatistics.m_iterationsStats;
  }
  ConvergenceStatistics & getConvergenceStats()
  {
    return m_solverStatistics.m_convergenceStats;
  }
  ConvergenceStatistics const & getConvergenceStats() const
  {
    return m_solverStatistics.m_convergenceStats;
  }
 
  SolverStatistics & getSolverStatistics() { return m_solverStatistics; }
 
  SolverStatistics const & getSolverStatistics() const { return m_solverStatistics; }
 
#if defined(GEOS_USE_PYGEOSX)
  virtual PyTypeObject * getPythonType() const override;
#endif
 
  map< std::pair< string, string >, string_array > const & getMeshTargets() const
  {
    return m_meshTargets;
  }
 
  bool detectOscillations() const;
 
protected:
 
  virtual void postInputInitialization() override;
 
  real64 eisenstatWalker( real64 const newNewtonNorm,
                          real64 const oldNewtonNorm,
                          LinearSolverParameters::Krylov const & krylovParams );
 
  template< typename CONSTITUTIVE_BASE_TYPE >
  static string getConstitutiveName( ElementSubRegionBase const & subRegion );
 
  template< typename CONSTITUTIVE_BASE_TYPE >
  static string getConstitutiveName( ParticleSubRegionBase const & subRegion ); // particle overload
 
  template< typename CONSTITUTIVE >
  void setConstitutiveName( ElementSubRegionBase & subRegion, string const & wrapperName, string const & constitutiveType ) const;
 
  virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const { GEOS_UNUSED_VAR( subRegion ); }
 
 
  template< typename BASETYPE = constitutive::ConstitutiveBase, typename LOOKUP_TYPE >
  static BASETYPE const & getConstitutiveModel( dataRepository::Group const & dataGroup, LOOKUP_TYPE const & key )
  {
    dataRepository::Group const & constitutiveModels = dataGroup.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() );
    return constitutiveModels.getGroup< BASETYPE >( key );
  }
 
  template< typename BASETYPE = constitutive::ConstitutiveBase, typename LOOKUP_TYPE >
  static BASETYPE & getConstitutiveModel( dataRepository::Group & dataGroup, LOOKUP_TYPE const & key )
  {
    dataRepository::Group & constitutiveModels = dataGroup.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() );
    return constitutiveModels.getGroup< BASETYPE >( key );
  }
 
  template< typename CONSTITUTIVE_TYPE >
  static CONSTITUTIVE_TYPE & getConstitutiveModel( ElementSubRegionBase & subRegion )
  {
    return getConstitutiveModel< CONSTITUTIVE_TYPE >( subRegion, getConstitutiveName< CONSTITUTIVE_TYPE >( subRegion ) );
  }
 
  real64 m_cflFactor;
 
  real64 m_nextDt;
 
  integer m_allowNonConvergedLinearSolverSolution;
 
  integer m_numTimestepsSinceLastDtCut;
 
  string m_discretizationName;
 
  DofManager m_dofManager;
 
  ParallelMatrix m_matrix;
 
  ParallelVector m_rhs;
 
  ParallelVector m_solution;
 
  ParallelVector m_scaling;
 
  integer m_usePhysicsScaling;
 
  CRSMatrix< real64, globalIndex > m_localMatrix;
 
  std::unique_ptr< LinearSolverBase< LAInterface > > m_linearSolver;
 
  std::unique_ptr< PreconditionerBase< LAInterface > > m_precond;
 
  integer m_writeLinearSystem;
 
  StatsOutputType m_writeStatisticsCSV;
 
  LinearSolverParametersInput m_linearSolverParameters;
 
  LinearSolverResult m_linearSolverResult;
 
  NonlinearSolverParameters m_nonlinearSolverParameters;
 
  SolverStatistics m_solverStatistics;
 
  Timestamp m_systemSetupTimestamp;
 
  std::function< void( CRSMatrix< real64, globalIndex >, array1d< real64 > ) > m_assemblyCallback;
 
  std::map< std::string, std::chrono::system_clock::duration > m_timers;
 
  ArrayOfArrays< real64 > m_solutionHistory;
 
private:
  string_array m_targetRegionNames;
 
  map< std::pair< string, string >, string_array > m_meshTargets;
 
  virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const { GEOS_UNUSED_VAR( subRegion ); }
 
  bool solveNonlinearSystem( real64 const & time_n,
                             real64 const & dt,
                             integer const cycleNumber,
                             DomainPartition & domain );
 
  void logEndOfCycleInformation( integer const cycleNumber,
                                 integer const numOfSubSteps,
                                 stdVector< real64 > const & subStepDts ) const;
};
 
 
template< typename CONSTITUTIVE_BASE_TYPE >
string PhysicsSolverBase::getConstitutiveName( ElementSubRegionBase const & subRegion )
{
  string validName;
  dataRepository::Group const & constitutiveModels = subRegion.getConstitutiveModels();
 
  constitutiveModels.forSubGroups< CONSTITUTIVE_BASE_TYPE >( [&]( dataRepository::Group const & model )
  {
    GEOS_ERROR_IF( !validName.empty(), "A valid constitutive model was already found." );
    validName = model.getName();
  } );
 
  return validName;
}
 
template< typename CONSTITUTIVE_BASE_TYPE >
string PhysicsSolverBase::getConstitutiveName( ParticleSubRegionBase const & subRegion ) // particle overload
{
  string validName;
  dataRepository::Group const & constitutiveModels = subRegion.getConstitutiveModels();
 
  constitutiveModels.forSubGroups< CONSTITUTIVE_BASE_TYPE >( [&]( dataRepository::Group const & model )
  {
    GEOS_ERROR_IF( !validName.empty(), "A valid constitutive model was already found." );
    validName = model.getName();
  } );
  return validName;
}
 
template< typename CONSTITUTIVE >
void PhysicsSolverBase::setConstitutiveName( ElementSubRegionBase & subRegion, string const & wrapperName, string const & constitutiveType ) const
{
  subRegion.registerWrapper< string >( wrapperName ).
    setPlotLevel( dataRepository::PlotLevel::NOPLOT ).
    setRestartFlags( dataRepository::RestartFlags::NO_WRITE ).
    setSizedFromParent( 0 );
 
  string & constitutiveName = subRegion.getReference< string >( wrapperName );
  constitutiveName = getConstitutiveName< CONSTITUTIVE >( subRegion );
  GEOS_ERROR_IF( constitutiveName.empty(), GEOS_FMT( "{}: {} constitutive model not found on subregion {}",
                                                     getDataContext(), constitutiveType, subRegion.getName() ) );
}
 
ENUM_STRINGS( PhysicsSolverBase::StatsOutputType,
              "none",
              "iteration",
              "convergence",
              "all" );
 
} // namespace geos
 
 
#endif /* GEOS_PHYSICSSOLVERS_PHYSICSSOLVERBASE_HPP_ */