CompositionalMultiphaseBase.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_FLUIDFLOW_COMPOSITIONALMULTIPHASEBASE_HPP_
#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEBASE_HPP_
 
#include "physicsSolvers/fluidFlow/FlowSolverBase.hpp"
#include "fieldSpecification/FieldSpecificationManager.hpp"
#include "constitutive/fluid/multifluid/MultiFluidBase.hpp"
#include "constitutive/solid/CoupledSolidBase.hpp"
#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp"
#include "physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp"
 
namespace geos
{
 
enum class CompositionalMultiphaseFormulationType : integer
{
  ComponentDensities, 
  OverallComposition  
};
 
ENUM_STRINGS( CompositionalMultiphaseFormulationType,
              "ComponentDensities",
              "OverallComposition" );
 
//START_SPHINX_INCLUDE_00
class CompositionalMultiphaseBase : public FlowSolverBase
{
public:
 
  CompositionalMultiphaseBase( const string & name,
                               Group * const parent );
 
  CompositionalMultiphaseBase() = delete;
 
  CompositionalMultiphaseBase( CompositionalMultiphaseBase const & ) = delete;
 
  CompositionalMultiphaseBase( CompositionalMultiphaseBase && ) = default;
 
  CompositionalMultiphaseBase & operator=( CompositionalMultiphaseBase const & ) = delete;
 
  CompositionalMultiphaseBase & operator=( CompositionalMultiphaseBase && ) = delete;
 
  virtual ~CompositionalMultiphaseBase() override = default;
 
//START_SPHINX_INCLUDE_01
 
  virtual void registerDataOnMesh( Group & meshBodies ) override;
 
  virtual void registerDataForCFL( Group & meshBodies ) { GEOS_UNUSED_VAR( meshBodies ); }
 
  virtual void
  implicitStepSetup( real64 const & time_n,
                     real64 const & dt,
                     DomainPartition & domain ) override;
 
  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;
 
  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;
 
  template< typename SUBREGION_TYPE >
  void accumulationAssemblyLaunch( DofManager const & dofManager,
                                   SUBREGION_TYPE const & subRegion,
                                   CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                   arrayView1d< real64 > const & localRhs );
 
  virtual void
  resetStateToBeginningOfStep( DomainPartition & domain ) override;
 
  virtual void
  implicitStepComplete( real64 const & time,
                        real64 const & dt,
                        DomainPartition & domain ) override;
 
  void updateGlobalComponentFraction( ObjectManagerBase & dataGroup ) const;
 
  real64 updatePhaseVolumeFraction( ObjectManagerBase & dataGroup ) const;
 
  void updateFluidModel( ObjectManagerBase & dataGroup ) const;
 
  void updateRelPermModel( ObjectManagerBase & dataGroup ) const;
 
  void updateCapPressureModel( ObjectManagerBase & dataGroup ) const;
 
  void updateCompAmount( ElementSubRegionBase & subRegion ) const;
 
  void updateEnergy( ElementSubRegionBase & subRegion ) const;
 
  void updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const;
 
  virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const = 0;
 
  real64 updateFluidState( ElementSubRegionBase & subRegion ) const;
 
  virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override final;
 
  virtual void saveSequentialIterationState( DomainPartition & domain ) override final;
 
  virtual void updateState( DomainPartition & domain ) override final;
 
  integer numFluidComponents() const { return m_numComponents; }
 
  integer numFluidPhases() const { return m_numPhases; }
 
  string referenceFluidModelName() const { return m_referenceFluidModelName; }
 
  virtual units::Unit getMassUnit() const override
  { return m_useMass ? units::Unit::Mass : units::Unit::Mole; }
 
  void assembleLocalTerms( DomainPartition & domain,
                           DofManager const & dofManager,
                           CRSMatrixView< real64, globalIndex const > const & localMatrix,
                           arrayView1d< real64 > const & localRhs ) const;
 
  virtual void
  assembleFluxTerms( real64 const dt,
                     DomainPartition const & domain,
                     DofManager const & dofManager,
                     CRSMatrixView< real64, globalIndex const > const & localMatrix,
                     arrayView1d< real64 > const & localRhs ) const = 0;
 
  virtual void
  assembleStabilizedFluxTerms( real64 const dt,
                               DomainPartition const & domain,
                               DofManager const & dofManager,
                               CRSMatrixView< real64, globalIndex const > const & localMatrix,
                               arrayView1d< real64 > const & localRhs ) const = 0;
 
  struct viewKeyStruct : FlowSolverBase::viewKeyStruct
  {
    static constexpr char const * elemDofFieldString() { return "compositionalVariables"; }
 
    // inputs
 
    static constexpr char const * useMassFlagString() { return "useMass"; }
    static constexpr char const * formulationTypeString() { return "formulationType"; }
 
    // time stepping controls
 
    static constexpr char const * solutionChangeScalingFactorString() { return "solutionChangeScalingFactor"; }
    static constexpr char const * targetRelativePresChangeString() { return "targetRelativePressureChangeInTimeStep"; }
    static constexpr char const * targetRelativeTempChangeString() { return "targetRelativeTemperatureChangeInTimeStep"; }
    static constexpr char const * targetPhaseVolFracChangeString() { return "targetPhaseVolFractionChangeInTimeStep"; }
    static constexpr char const * targetRelativeCompDensChangeString() { return "targetRelativeCompDensChangeInTimeStep"; }
    static constexpr char const * targetCompFracChangeString() { return "targetCompFracChangeInTimeStep"; }
    static constexpr char const * targetFlowCFLString() { return "targetFlowCFL"; }
 
 
    // nonlinear solver parameters
 
    static constexpr char const * maxCompFracChangeString() { return "maxCompFractionChange"; }
    static constexpr char const * maxRelativePresChangeString() { return "maxRelativePressureChange"; }
    static constexpr char const * maxRelativeTempChangeString() { return "maxRelativeTemperatureChange"; }
    static constexpr char const * maxRelativeCompDensChangeString() { return "maxRelativeCompDensChange"; }
    static constexpr char const * allowLocalCompDensChoppingString() { return "allowLocalCompDensityChopping"; }
    static constexpr char const * useTotalMassEquationString() { return "useTotalMassEquation"; }
    static constexpr char const * useSimpleAccumulationString() { return "useSimpleAccumulation"; }
    static constexpr char const * minCompDensString() { return "minCompDens"; }
    static constexpr char const * minCompFracString() { return "minCompFrac"; }
    static constexpr char const * maxSequentialCompDensChangeString() { return "maxSequentialCompDensChange"; }
    static constexpr char const * minScalingFactorString() { return "minScalingFactor"; }
 
    static constexpr char const * relPermNamesString() { return "relPermNames"; }
    static constexpr char const * capPressureNamesString() { return "capPressureNames"; }
    static constexpr char const * diffusionNamesString() { return "diffusionNames"; }
    static constexpr char const * dispersionNamesString() { return "dispersionNames"; }
  };
 
  virtual void initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) override;
 
  virtual void initializeThermalState( MeshLevel & mesh, string_array const & regionNames ) override;
 
  virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override;
 
  void applyDirichletBC( real64 const time,
                         real64 const dt,
                         DofManager const & dofManager,
                         DomainPartition & domain,
                         CRSMatrixView< real64, globalIndex const > const & localMatrix,
                         arrayView1d< real64 > const & localRhs ) const;
 
  void applySourceFluxBC( real64 const time,
                          real64 const dt,
                          DofManager const & dofManager,
                          DomainPartition & domain,
                          CRSMatrixView< real64, globalIndex const > const & localMatrix,
                          arrayView1d< real64 > const & localRhs ) const;
 
  virtual void applyAquiferBC( real64 const time,
                               real64 const dt,
                               DofManager const & dofManager,
                               DomainPartition & domain,
                               CRSMatrixView< real64, globalIndex const > const & localMatrix,
                               arrayView1d< real64 > const & localRhs ) const = 0;
 
  void keepVariablesConstantDuringInitStep( real64 const time,
                                            real64 const dt,
                                            DofManager const & dofManager,
                                            DomainPartition & domain,
                                            CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                            arrayView1d< real64 > const & localRhs ) const;
 
 
  void chopNegativeDensities( DomainPartition & domain );
  void chopNegativeDensities( ElementSubRegionBase & subRegion );
 
  void chopNegativeCompFractions( DomainPartition & domain );
 
  virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt,
                                              DomainPartition & domain ) override;
 
 
  virtual void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL )
  {
    GEOS_UNUSED_VAR( domain, dt, maxPhaseCFL, maxCompCFL );
    GEOS_ERROR( GEOS_FMT( "{}: computeCFLNumbers is not implemented for {}", getDataContext(), getCatalogName() ) );
  }
 
  virtual void initializePostInitialConditionsPreSubGroups() override;
 
  integer useSimpleAccumulation() const { return m_useSimpleAccumulation; }
 
  integer useTotalMassEquation() const { return m_useTotalMassEquation; }
 
  virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const override;
 
protected:
 
  virtual void postInputInitialization() override;
 
  virtual void initializePreSubGroups() override;
 
 
  void validateConstitutiveModels( DomainPartition const & domain ) const;
 
  void initializeAquiferBC( constitutive::ConstitutiveManager const & cm ) const;
 
  template< typename OBJECT_TYPE >
  void applyFieldValue( real64 const & time_n,
                        real64 const & dt,
                        MeshLevel & mesh,
                        char const logMessage[],
                        string const fieldKey,
                        string const boundaryFieldKey ) const;
 
  integer m_numPhases;
 
  integer m_numComponents;
 
  integer m_useMass;
 
  CompositionalMultiphaseFormulationType m_formulationType;
 
  bool m_hasCapPressure;
 
  bool m_hasDiffusion;
 
  bool m_hasDispersion;
 
  real64 m_maxCompFracChange;
 
  real64 m_maxRelativePresChange;
 
  real64 m_maxRelativeTempChange;
 
  real64 m_maxRelativeCompDensChange;
 
  real64 m_solutionChangeScalingFactor;
 
  real64 m_targetRelativePresChange;
 
  real64 m_targetRelativeTempChange;
 
  real64 m_targetPhaseVolFracChange;
 
  real64 m_targetRelativeCompDensChange;
 
  real64 m_targetCompFracChange;
 
  real64 m_minScalingFactor;
 
  integer m_allowCompDensChopping;
 
  integer m_useTotalMassEquation;
 
  integer m_useSimpleAccumulation;
 
  real64 m_minCompDens;
 
  real64 m_minCompFrac;
 
  string m_referenceFluidModelName;
 
  real64 m_sequentialCompDensChange;
  real64 m_maxSequentialCompDensChange;
 
private:
 
  bool validateDirichletBC( DomainPartition & domain,
                            real64 const time ) const;
 
  virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override;
 
};
 
template< typename OBJECT_TYPE >
void CompositionalMultiphaseBase::applyFieldValue( real64 const & time_n,
                                                   real64 const & dt,
                                                   MeshLevel & mesh,
                                                   char const logMessage[],
                                                   string const fieldKey,
                                                   string const boundaryFieldKey ) const
{
  FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance();
 
  fsManager.apply< OBJECT_TYPE >( time_n + dt,
                                  mesh,
                                  fieldKey,
                                  [&]( FieldSpecificationBase const & fs,
                                       string const & setName,
                                       SortedArrayView< localIndex const > const & lset,
                                       OBJECT_TYPE & targetGroup,
                                       string const & )
  {
    if( fs.getLogLevel() >= 1 && m_nonlinearSolverParameters.m_numNewtonIterations == 0 )
    {
      globalIndex const numTargetElems = MpiWrapper::sum< globalIndex >( lset.size() );
      GEOS_LOG_RANK_0( GEOS_FMT( logMessage,
                                 getName(), time_n+dt, fs.getCatalogName(), fs.getName(),
                                 setName, targetGroup.getName(), fs.getScale(), numTargetElems ) );
    }
 
    // Specify the bc value of the field
    fs.applyFieldValue< FieldSpecificationEqual,
                        parallelDevicePolicy<> >( lset,
                                                  time_n + dt,
                                                  targetGroup,
                                                  boundaryFieldKey );
  } );
}
 
template< typename SUBREGION_TYPE >
void CompositionalMultiphaseBase::accumulationAssemblyLaunch( DofManager const & dofManager,
                                                              SUBREGION_TYPE const & subRegion,
                                                              CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                                              arrayView1d< real64 > const & localRhs )
{
  constitutive::MultiFluidBase const & fluid =
    getConstitutiveModel< constitutive::MultiFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) );
  constitutive::CoupledSolidBase const & solid =
    getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );
 
  string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
 
  using namespace isothermalCompositionalMultiphaseBaseKernels;
 
  BitFlags< KernelFlags > kernelFlags;
  if( m_useTotalMassEquation )
    kernelFlags.set( KernelFlags::TotalMassEquation );
  if( m_useSimpleAccumulation )
    kernelFlags.set( KernelFlags::SimpleAccumulation );
 
  if( m_isThermal )
  {
    thermalCompositionalMultiphaseBaseKernels::
      AccumulationKernelFactory::
      createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
                                                 m_numPhases,
                                                 dofManager.rankOffset(),
                                                 kernelFlags,
                                                 dofKey,
                                                 subRegion,
                                                 fluid,
                                                 solid,
                                                 localMatrix,
                                                 localRhs );
  }
  else
  {
    isothermalCompositionalMultiphaseBaseKernels::
      AccumulationKernelFactory::
      createAndLaunch< parallelDevicePolicy<> >( m_numComponents,
                                                 m_numPhases,
                                                 dofManager.rankOffset(),
                                                 kernelFlags,
                                                 dofKey,
                                                 subRegion,
                                                 fluid,
                                                 solid,
                                                 localMatrix,
                                                 localRhs );
  }
}
 
} // namespace geos
 
#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEBASE_HPP_