WellNewtonSolver.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_WELLNEWTONSOLVER_HPP_
#define GEOS_PHYSICSSOLVERS_WELLNEWTONSOLVER_HPP_
 
#include "codingUtilities/traits.hpp"
#include "common/DataTypes.hpp"
#include "common/format/LogPart.hpp"
 
#include "dataRepository/RestartFlags.hpp"
#include "linearAlgebra/interfaces/InterfaceTypes.hpp"
#include "linearAlgebra/utilities/LinearSolverResult.hpp"
#include "linearAlgebra/DofManager.hpp"
#include "mesh/DomainPartition.hpp"
#include "mesh/MeshBody.hpp"
#include "physicsSolvers/NonlinearSolverParameters.hpp"
#include "physicsSolvers/LinearSolverParameters.hpp"
#include "physicsSolvers/SolverStatistics.hpp"
#include "physicsSolvers/LogLevelsInfo.hpp"
 
#include "physicsSolvers/fluidFlow/SolutionCheckHelpers.hpp"
#include "common/Timer.hpp"
#include <limits>
 
namespace geos
{
 
class DomainPartition;
 
class WellNewtonSolver : public dataRepository::Group
{
public:
 
  enum class StatsOutputType : integer
  {
    none, iteration, convergence, all
  };
 
  explicit WellNewtonSolver( string const & name,
                             Group * const parent );
 
  WellNewtonSolver( WellNewtonSolver && ) = default;
 
  virtual ~WellNewtonSolver() override;
 
  WellNewtonSolver() = delete;
 
  WellNewtonSolver( WellNewtonSolver const & ) = delete;
 
  WellNewtonSolver & operator=( WellNewtonSolver const & ) = delete;
 
  WellNewtonSolver & operator=( WellNewtonSolver && ) = delete;
 
  static string catalogName() { return "WellNewtonSolver"; }
 
 
 
  void generateMeshTargetsFromTargetRegions( Group const & meshBodies );
 
 
  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(); }
 
 
  template< typename T >
  void  setupSystem( T & well, DomainPartition & domain,
                     std::string const & meshBodyName,
                     MeshLevel const & meshLevel,
                     WellElementRegion & wellElementRegion,
                     bool const setSparsity =true );
 
  template< typename T >
  bool
  solveNonlinearSystem( T & well, real64 const & time_n,
                        real64 const & stepDt,
                        integer const cycleNumber,
                        DomainPartition & domain,
                        MeshLevel & mesh,
                        ElementRegionManager & elemManager,
                        WellElementSubRegion & subRegion );
 
 
  virtual void
  setupDofs( DomainPartition const & domain,
             DofManager & dofManager ) const;
 
  virtual void
  setSparsityPattern( DomainPartition & domain,
                      DofManager & dofManager,
                      CRSMatrix< real64, globalIndex > & localMatrix,
                      SparsityPattern< globalIndex > & pattern );
 
  virtual std::unique_ptr< PreconditionerBase< LAInterface > >
  createPreconditioner( DomainPartition & domain ) const;
 
 
  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 void
  solveLinearSystem( DofManager const & dofManager,
                     ParallelMatrix & matrix,
                     ParallelVector & rhs,
                     ParallelVector & solution );
 
 
 
  virtual Group * createChild( string const & childKey, string const & childName ) override;
 
  using CatalogInterface = dataRepository::CatalogInterface< WellNewtonSolver, 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"; }
 
    static constexpr char const * activeCoupledIterationsString() { return "activeCoupledIterations"; }
    static constexpr char const * estimateWellSolutionString() { return "estimateWellSolution"; }
    static constexpr char const * enableIsoThermalEstimatorString() { return "enableIsoThermalEstimator"; }
  };
 
  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; }
 
 
  void setSystemSetupTimestamp( Timestamp timestamp );
 
 
  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;}
 
 
 
  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; }
 
 
 
  bool detectOscillations() const;
 
 
  void enableThermalEffects ( bool enable ) { m_thermalEffectsEnabled = enable; };
 
  bool thermalEffectsEnabled() const { return m_thermalEffectsEnabled; }
 
  bool isoThermalEstimatorEnabled() const { return m_enableIsoThermalEstimator; }
 
  bool getNumActiveCoupledIterations() const { return m_activeCoupledIterations; }
 
protected:
 
  virtual void postInputInitialization() override;
 
  integer m_allowNonConvergedLinearSolverSolution;
 
 
 
  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;
 
  stdMap< 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;
 
  integer m_activeCoupledIterations;
 
  bool m_thermalEffectsEnabled;
  integer m_enableIsoThermalEstimator;
 
 
  void logEndOfCycleInformation( integer const cycleNumber,
                                 integer const numOfSubSteps,
                                 stdVector< real64 > const & subStepDts ) const;
};
 
template< typename T >
void WellNewtonSolver::setupSystem( T & well, DomainPartition & domain,
                                    std::string const & meshBodyName,
                                    MeshLevel const & meshLevel,
                                    WellElementRegion & wellElementRegion,
                                    bool const setSparsity )
{
  GEOS_MARK_FUNCTION;
 
  map< std::pair< string, string >, string_array > meshTargets;
  string_array regions;
 
  meshTargets.clear();
  regions.clear();
  regions.emplace_back( wellElementRegion.getName() );
  auto const key = std::make_pair( meshBodyName, meshLevel.getName() );
  meshTargets[key] = std::move( regions );
 
  m_dofManager.setDomain( domain );
  m_dofManager.addField( well.wellElementDofName(),
                         FieldLocation::Elem,
                         well.numDofPerWellElement(),
                         meshTargets );
 
  m_dofManager.addCoupling( well.wellElementDofName(),
                            well.wellElementDofName(),
                            DofManager::Connector::Node );
 
  m_dofManager.reorderByRank();
  if( setSparsity )
  {
    SparsityPattern< globalIndex > pattern;
    setSparsityPattern( domain, m_dofManager, m_localMatrix, pattern );
    m_localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) );
  }
  m_localMatrix.setName( this->getName() + "/matrix" );
 
  m_rhs.setName( this->getName() + "/rhs" );
  m_rhs.create( m_dofManager.numLocalDofs(), MPI_COMM_GEOS );
 
  m_solution.setName( this->getName() + "/solution" );
  m_solution.create( m_dofManager.numLocalDofs(), MPI_COMM_GEOS );
}
 
 
template< typename T >
bool WellNewtonSolver::solveNonlinearSystem( T & well, real64 const & time_n,
                                             real64 const & stepDt,
                                             integer const cycleNumber,
                                             DomainPartition & domain,
                                             MeshLevel & mesh,
                                             ElementRegionManager & elemManager,
                                             WellElementSubRegion & subRegion )
{
  integer const maxNewtonIter = m_nonlinearSolverParameters.m_maxIterNewton;
  integer const minNewtonIter = m_nonlinearSolverParameters.m_minIterNewton;
  real64 const newtonTol = m_nonlinearSolverParameters.m_newtonTol;
 
// keep residual from previous iteration in case we need to do a line search
 
  integer newtonIter = 0;
  real64 scaleFactor = 1.0;
 
  bool isNewtonConverged = false;
 
  for( newtonIter = 0; newtonIter < maxNewtonIter; ++newtonIter )
  {
    if( m_nonlinearSolverParameters.getLogLevel() > 4 )
      GEOS_LOG_LEVEL_RANK_0( logInfo::NonlinearSolver,
                             GEOS_FMT( " Well: {}   Est Attempt: NewtonIter: {:2}", subRegion.getName(), stepDt, newtonIter ));
 
    {
      Timer timer( m_timers.get_inserted( "assemble" ) );
 
// We sync the nonlinear convergence history. The coupled solver parameters are the one being
// used. We want to propagate the info to subsolvers. It can be important for solvers that
// have special treatment for specific iterations.
      synchronizeNonlinearSolverParameters();
 
// zero out matrix/rhs before assembly
      m_localMatrix.zero();
      m_rhs.zero();
 
      arrayView1d< real64 > const localRhs = m_rhs.open();
 
// call assemble to fill the matrix and the rhs
      well.assembleSystem( time_n,
                           stepDt,
                           cycleNumber,
                           elemManager,
                           subRegion,
                           m_dofManager,
                           m_localMatrix.toViewConstSizes(),
                           localRhs );
 
// apply boundary conditions to system
      well.applyWellBoundaryConditions( time_n,
                                        stepDt,
                                        elemManager,
                                        subRegion,
                                        m_dofManager,
                                        localRhs,
                                        m_localMatrix.toViewConstSizes() );
 
      m_rhs.close();
 
      if( m_assemblyCallback )
      {
// Make a copy of LA objects and ship off to the callback
        array1d< real64 > localRhsCopy( m_rhs.localSize() );
        localRhsCopy.setValues< parallelDevicePolicy<> >( m_rhs.values() );
        m_assemblyCallback( m_localMatrix, std::move( localRhsCopy ) );
      }
    }
 
    // well.outputSingleWellDebug( time_n, stepDt, 0, newtonIter, 0,
    //                             mesh, subRegion, dofManager, m_localMatrix.toViewConstSizes(), m_rhs.values()  );
    real64 residualNorm = 0;
    {
      Timer timer( m_timers.get_inserted( "convergence check" ) );
 
// get residual norm
      residualNorm = well.calculateWellResidualNorm( time_n, stepDt, m_nonlinearSolverParameters, subRegion, m_dofManager, m_rhs.values() );
      if( m_nonlinearSolverParameters.getLogLevel() > 4 )
        GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence,
                               GEOS_FMT( "        ( R ) = ( {:4.2e} )", residualNorm ) );
    }
    //auto iterInfo = currentIter( time_n, dt );
    //outputSingleWellDebug( time_n, stepDt, 0, newtonIter, 0,
    //                       mesh, subRegion, dofManager, m_localMatrix.toViewConstSizes(), m_rhs.values()  );
    // if the residual norm is less than the Newton tolerance we denote that we have
    // converged and break from the Newton loop immediately.
    if( residualNorm < newtonTol && newtonIter >= minNewtonIter )
    {
      isNewtonConverged = true;
      break;
    }
 
// if the residual norm is above the max allowed residual norm, we break from
// the Newton loop to avoid crashes due to Newton divergence
    if( residualNorm > m_nonlinearSolverParameters.m_maxAllowedResidualNorm )
    {
      string const maxAllowedResidualNormString = NonlinearSolverParameters::viewKeysStruct::maxAllowedResidualNormString();
      if( m_nonlinearSolverParameters.getLogLevel() > 4 )
        GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence,
                               GEOS_FMT( "    The residual norm is above the {} of {}. Newton loop terminated.",
                                         maxAllowedResidualNormString,
                                         m_nonlinearSolverParameters.m_maxAllowedResidualNorm )  );
      isNewtonConverged = false;
      break;
    }
 
    {
      Timer timer( m_timers.get_inserted( "linear solver total" ) );
 
// TODO: Trilinos currently requires this, re-evaluate after moving to Tpetra-based solvers
      if( m_precond )
      {
        m_precond->clear();
      }
 
      {
        Timer timer_setup( m_timers.get_inserted( "linear solver create" ) );
 
// Compose parallel LA matrix/rhs out of local LA matrix/rhs
//
        m_matrix.create( m_localMatrix.toViewConst(), m_dofManager.numLocalDofs(), MPI_COMM_GEOS );
      }
 
// Output the linear system matrix/rhs for debugging purposes
      //string tag = "_"+std::to_string( my_ctime );  tjb
      //debugOutputSystem( time_n, cycleNumber, newtonIter, m_matrix, m_rhs, tag );
 
      debugOutputSystem( time_n, cycleNumber, newtonIter, m_matrix, m_rhs );
// Solve the linear system
      solveLinearSystem( m_dofManager, m_matrix, m_rhs, m_solution );
 
// Increment the solver statistics for reporting purposes
      getIterationStats().updateNonlinearIteration( m_linearSolverResult.numIterations );
 
// Output the linear system solution for debugging purposes
      debugOutputSolution( time_n, cycleNumber, newtonIter, m_solution );
      //debugOutputSolution( time_n, cycleNumber, newtonIter, m_solution, tag );
    }
 
    {
      Timer timer( m_timers.get_inserted( "apply solution" ) );
 
// Compute the scaling factor for the Newton update
      scaleFactor = well.scalingForWellSystemSolution( subRegion, m_dofManager, m_solution.values() );
      if( m_nonlinearSolverParameters.getLogLevel() > 4 )
        GEOS_LOG_LEVEL_RANK_0( logInfo::Solution,
                               GEOS_FMT( "        {}: Global solution scaling factor = {}", getName(), scaleFactor ) );
 
      real64 minPressure = 0.0, minDensity = 0.0, minTotalDensity = 0.0;
      bool const solutionLogActive = isLogLevelActive< logInfo::Solution >( getLogLevel() );
      bool const solutionDetailsLogActive = isLogLevelActive< logInfo::SolutionDetails >( getLogLevel() );
      ElementsReporterBuffer rankNegPressureIds{ solutionLogActive, solutionDetailsLogActive ? 16 : 0 };
      ElementsReporterBuffer rankNegDensityIds{ solutionLogActive, solutionDetailsLogActive ? 16 : 0 };
      // output only total density sum, not cell details
      ElementsReporterBuffer rankTotalNegDensityIds{ solutionLogActive, 0 };
      if( !well.checkWellSystemSolution( subRegion, m_dofManager, m_solution.values(), scaleFactor, minPressure, minDensity, minTotalDensity, rankNegPressureIds, rankNegDensityIds,
                                         rankTotalNegDensityIds ) )
      {
// TODO try chopping (similar to line search)
        if( m_nonlinearSolverParameters.getLogLevel() > 4 )
          GEOS_LOG_RANK_0( GEOS_FMT( "    {}: Solution check failed. Newton loop terminated.", getName()) );
        break;
      }
 
// apply the system solution to the fields/variables
      well.applyWellSystemSolution( m_dofManager, m_solution.values(), scaleFactor, stepDt, domain, mesh, subRegion );
    }
 
    {
      Timer timer( m_timers.get_inserted( "update state" ) );
 
      // update derived variables (constitutive models)
      well.updateWellState( domain.getMeshBody( mesh.getParent().getParent().getName() ), elemManager, subRegion );
    }
 
  }
  return isNewtonConverged;
}
 
 
ENUM_STRINGS( WellNewtonSolver::StatsOutputType,
              "none",
              "iteration",
              "convergence",
              "all" );
 
} // namespace geos
 
 
#endif /* GEOS_PHYSICSSOLVERS_WELLNEWTONSOLVER_HPP_ */