FlowSolverBase.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_FINITEVOLUME_FLOWSOLVERBASE_HPP_
#define GEOS_PHYSICSSOLVERS_FINITEVOLUME_FLOWSOLVERBASE_HPP_
 
#include "physicsSolvers/PhysicsSolverBase.hpp"
#include "common/Units.hpp"
#include "finiteVolume/BoundaryStencil.hpp"
#include "fieldSpecification/AquiferBoundaryCondition.hpp"
 
namespace geos
{
 
class FlowSolverBase : public PhysicsSolverBase
{
  template< typename VIEWTYPE >
  using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
 
public:
 
  static string coupledSolverAttributePrefix() { return "flow"; }
 
  FlowSolverBase( const string & name,
                  Group * const parent );
 
 
  FlowSolverBase() = delete;
 
  FlowSolverBase( FlowSolverBase const & ) = delete;
 
  FlowSolverBase( FlowSolverBase && ) = default;
 
  FlowSolverBase & operator=( FlowSolverBase const & ) = delete;
 
  FlowSolverBase & operator=( FlowSolverBase && ) = delete;
 
  virtual void registerDataOnMesh( Group & MeshBodies ) override;
 
  struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct
  {
    // misc inputs
    static constexpr char const * isThermalString() { return "isThermal"; }
    static constexpr char const * inputTemperatureString() { return "temperature"; }
    static constexpr char const * allowNegativePressureString() { return "allowNegativePressure"; }
    static constexpr char const * maxAbsolutePresChangeString() { return "maxAbsolutePressureChange"; }
    static constexpr char const * maxSequentialPresChangeString() { return "maxSequentialPressureChange"; }
    static constexpr char const * maxSequentialTempChangeString() { return "maxSequentialTemperatureChange"; }
 
    static constexpr char const * fluidNamesString() { return "fluidNames"; }
    static constexpr char const * solidNamesString() { return "solidNames"; }
    static constexpr char const * permeabilityNamesString() { return "permeabilityNames"; }
    static constexpr char const * solidInternalEnergyNamesString() { return "solidInternalEnergyNames"; }
    static constexpr char const * thermalConductivityNamesString() { return "thermalConductivityNames"; }
  };
 
  void prepareStencilWeights( DomainPartition & domain ) const;
 
  void updateStencilWeights( DomainPartition & domain ) const;
 
  void enableFixedStressPoromechanicsUpdate() { m_isFixedStressPoromechanicsUpdate = true; }
 
  void enableJumpStabilization() { m_isJumpStabilized = true; }
 
  void updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const;
 
  virtual void updatePorosityAndPermeability( SurfaceElementSubRegion & subRegion ) const;
 
  virtual void saveSequentialIterationState( DomainPartition & domain ) override;
 
  integer & isThermal() { return m_isThermal; }
 
  virtual units::Unit getMassUnit() const { return units::Unit::Mass; }
 
  void allowNegativePressure() { m_allowNegativePressure = 1; }
 
  void setKeepVariablesConstantDuringInitStep( bool const keepVariablesConstantDuringInitStep )
  { m_keepVariablesConstantDuringInitStep = keepVariablesConstantDuringInitStep; }
 
  virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const override;
 
  void enableLaggingFractureStencilWeightsUpdate(){ m_isLaggingFractureStencilWeightsUpdate = 1; };
 
  real64 sumAquiferFluxes( BoundaryStencil const & stencil,
                           AquiferBoundaryCondition::KernelWrapper const & aquiferBCWrapper,
                           ElementViewConst< arrayView1d< real64 const > > const & pres,
                           ElementViewConst< arrayView1d< real64 const > > const & presOld,
                           ElementViewConst< arrayView1d< real64 const > > const & gravCoef,
                           real64 const & timeAtBeginningOfStep,
                           real64 const & dt );
 
  virtual void assembleHydrofracFluxTerms( real64 const time_n,
                                           real64 const dt,
                                           DomainPartition const & domain,
                                           DofManager const & dofManager,
                                           CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                           arrayView1d< real64 > const & localRhs,
                                           CRSMatrixView< real64, localIndex const > const & dR_dAper )
  {
    GEOS_UNUSED_VAR ( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper );
    GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." );
  }
 
  void initializeState( DomainPartition & domain );
 
  virtual void initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); }
 
  virtual void initializeThermalState( MeshLevel & mesh, string_array const & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); }
 
  void findMinMaxElevationInEquilibriumTarget( DomainPartition & domain, // cannot be const...
                                               std::map< string, localIndex > const & equilNameToEquilId,
                                               arrayView1d< real64 > const & maxElevation,
                                               arrayView1d< real64 > const & minElevation ) const;
 
  void computeSourceFluxSizeScalingFactor( real64 const & time,
                                           real64 const & dt,
                                           DomainPartition & domain, // cannot be const...
                                           std::map< string, localIndex > const & bcNameToBcId,
                                           arrayView1d< globalIndex > const & bcAllSetsSize ) const;
 
  integer numberOfDofsPerCell() const { return m_numDofPerCell; }
 
protected:
 
  virtual void saveAquiferConvergedState( real64 const & time,
                                          real64 const & dt,
                                          DomainPartition & domain );
 
  virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const;
 
  virtual void validatePoreVolumes( DomainPartition const & domain ) const;
 
  virtual void precomputeData( MeshLevel & mesh,
                               string_array const & regionNames );
 
  virtual void initializePreSubGroups() override;
 
  void checkDiscretizationName() const;
 
  virtual void initializePostInitialConditionsPreSubGroups() override;
 
  virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) { GEOS_UNUSED_VAR( domain ); }
 
  void initializePorosityAndPermeability( MeshLevel & mesh, string_array const & regionNames );
 
  void initializeHydraulicAperture( MeshLevel & mesh, string_array const & regionNames );
 
  void saveInitialPressureAndTemperature( MeshLevel & mesh, string_array const & regionNames );
 
  virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override;
 
  integer m_numDofPerCell;
 
  integer m_isThermal;
 
  real64 m_inputTemperature;
 
  bool m_keepVariablesConstantDuringInitStep;
 
  bool m_isFixedStressPoromechanicsUpdate;
 
  bool m_isJumpStabilized;
 
  integer m_allowNegativePressure;
 
  real64 m_maxAbsolutePresChange;
 
  real64 m_sequentialPresChange;
  real64 m_maxSequentialPresChange;
 
  real64 m_sequentialTempChange;
  real64 m_maxSequentialTempChange;
 
  class BCMessage
  {
public:
    static string pressureConflict( string_view regionName, string_view subRegionName,
                                    string_view setName, string_view fieldName );
 
    static string temperatureConflict( string_view regionName, string_view subRegionName,
                                       string_view setName, string_view fieldName );
 
    static string missingPressure( string_view regionName, string_view subRegionName,
                                   string_view setName, string_view fieldName );
 
    static string missingTemperature( string_view regionName, string_view subRegionName,
                                      string_view setName, string_view fieldName );
 
    static string conflictingComposition( int comp, string_view componentName,
                                          string_view regionName, string_view subRegionName,
                                          string_view setName, string_view fieldName );
 
    static string invalidComponentIndex( int comp,
                                         string_view fsName, string_view fieldName );
 
    static string notAppliedOnRegion( int componentIndex, string_view componentName,
                                      string_view regionName, string_view subRegionName,
                                      string_view setName, string_view fieldName );
private:
    static string generateMessage( string_view baseMessage,
                                   string_view fieldName, string_view setName );
  };
 
private:
  virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override;
 
  // flag to determine whether or not to apply lagging update for the fracture stencil weights
  integer m_isLaggingFractureStencilWeightsUpdate;
 
};
 
 
}
 
#endif //GEOS_PHYSICSSOLVERS_FINITEVOLUME_FLOWSOLVERBASE_HPP_