Constitutive models

In GEOS, all constitutive models defining fluid and rock properties are implemented in the namespace constitutive and derived from a common base class, ConstitutiveBase. All objects are owned and handled by the ConstitutiveManager.

Standalone models

Standalone constitutive models implement constitutive laws such as:

  • mechanical material models (linear elasticity, plasticity, etc.),

  • PVT fluid behaviors,

  • relative permeability relationships,

  • porosity and permeability dependencies on state variables,

  • contact laws.

Storage, allocation, and update of properties

Each constitutive model owns, as member variables, LvArray::Array containers that hold the properties (or fields) and their derivatives with respect to the other fields needed to update each property. Each property is stored as an array with the first dimension representing the elementIndex and the second dimension storing the index of the integration point. These dimensions are determined by the number of elements of the subregion on which each constitutive model is registered, and by the chosen discretization method. Vector and tensor fields have an additional dimension to identify their components. Similarly, an additional dimension is necessary for multiphase fluid models with properties defined for each component in each phase. For example, a single-phase fluid model where density and viscosity are functions of the fluid pressure has the following members:

  array2d< real64 > m_density;
  array2d< real64 > m_dDensity_dPressure;
  array2d< real64 > m_dDensity_dTemperature;

  array2d< real64 > m_density_n;

  array2d< real64 > m_viscosity;
  array2d< real64 > m_dViscosity_dPressure;
  array2d< real64 > m_dViscosity_dTemperature;

  array2d< real64 > m_internalEnergy;
  array2d< real64 > m_internalEnergy_n;
  array2d< real64 > m_dInternalEnergy_dPressure;
  array2d< real64 > m_dInternalEnergy_dTemperature;

  array2d< real64 > m_enthalpy;
  array2d< real64 > m_dEnthalpy_dPressure;
  array2d< real64 > m_dEnthalpy_dTemperature;

Resizing all fields of the constitutive models happens during the initialization phase by the ConstitutiveManger through a call to ConstitutiveManger::hangConstitutiveRelation, which sets the appropriate subRegion as the parent Group of each constitutive model object. This function also resizes all fields based on the size of the subregion and the number of quadrature points on it, by calling CONSTITUTIVE_MODEL::allocateConstitutiveData. For the single phase fluid example used before, this call is:

void SingleFluidBase::allocateConstitutiveData( Group & parent,
                                                localIndex const numConstitutivePointsPerParentIndex )
{
  ConstitutiveBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex );

  resize( parent.size() );

  m_density.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_dDensity_dPressure.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_dDensity_dTemperature.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_density_n.resize( parent.size(), numConstitutivePointsPerParentIndex );

  m_viscosity.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_dViscosity_dPressure.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_dViscosity_dTemperature.resize( parent.size(), numConstitutivePointsPerParentIndex );

  m_internalEnergy.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_internalEnergy_n.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_dInternalEnergy_dPressure.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_dInternalEnergy_dTemperature.resize( parent.size(), numConstitutivePointsPerParentIndex );

  m_enthalpy.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_dEnthalpy_dPressure.resize( parent.size(), numConstitutivePointsPerParentIndex );
  m_dEnthalpy_dTemperature.resize( parent.size(), numConstitutivePointsPerParentIndex );
}

Any property or field stored on a constitutive model must be updated within a computational kernel to ensure that host and device memory in GPUs are properly synced, and that any updates are performed on device. Some properties are updated within finite element kernels of specific physics (such as stress in a mechanics kernel). Consequently, for each constitutive model class, a corresponding nameOfTheModelUpdates, which only contains LvArray::arrayView containers to the data, can be captured by value inside computational kernels. For example, for the single phase fluid model Updates are:

/**
 * @brief Base class for single-phase fluid model kernel wrappers.
 */
class SingleFluidBaseUpdate
{
public:

  /**
   * @brief Get number of elements in this wrapper.
   * @return number of elements
   */
  GEOS_HOST_DEVICE
  localIndex numElems() const { return m_density.size( 0 ); }

  /**
   * @brief Get number of gauss points per element.
   * @return number of gauss points per element
   */
  GEOS_HOST_DEVICE
  localIndex numGauss() const { return m_density.size( 1 ); };

protected:

  /**
   * @brief Constructor.
   * @param density     fluid density
   * @param dDens_dPres derivative of density w.r.t. pressure
   * @param viscosity   fluid viscosity
   * @param dVisc_dPres derivative of viscosity w.r.t. pressure
   */
  SingleFluidBaseUpdate( arrayView2d< real64 > const & density,
                         arrayView2d< real64 > const & dDens_dPres,
                         arrayView2d< real64 > const & viscosity,
                         arrayView2d< real64 > const & dVisc_dPres )
    : m_density( density ),
    m_dDens_dPres( dDens_dPres ),
    m_viscosity( viscosity ),
    m_dVisc_dPres( dVisc_dPres )
  {}

  /**
   * @brief Copy constructor.
   */
  SingleFluidBaseUpdate( SingleFluidBaseUpdate const & ) = default;

  /**
   * @brief Move constructor.
   */
  SingleFluidBaseUpdate( SingleFluidBaseUpdate && ) = default;

  /**
   * @brief Deleted copy assignment operator
   * @return reference to this object
   */
  SingleFluidBaseUpdate & operator=( SingleFluidBaseUpdate const & ) = delete;

  /**
   * @brief Deleted move assignment operator
   * @return reference to this object
   */
  SingleFluidBaseUpdate & operator=( SingleFluidBaseUpdate && ) = delete;


  /// Fluid density
  arrayView2d< real64 > m_density;

  /// Derivative of density w.r.t. pressure
  arrayView2d< real64 > m_dDens_dPres;

  /// Fluid viscosity
  arrayView2d< real64 > m_viscosity;

  /// Derivative of viscosity w.r.t. pressure
  arrayView2d< real64 > m_dVisc_dPres;

Because Updates classes are responsible for updating the fields owned by the constitutive models, they also implement all functions needed to perform property updates, such as:

private:
  /**
   * @brief Compute fluid properties and derivatives at a single point.
   * @param[in]  pressure the target pressure value
   * @param[out] density fluid density
   * @param[out] dDensity_dPressure fluid density derivative w.r.t. pressure
   * @param[out] viscosity fluid viscosity
   * @param[out] dViscosity_dPressure fluid viscosity derivative w.r.t. pressure
   */
  GEOS_HOST_DEVICE
  virtual void compute( real64 const pressure,
                        real64 & density,
                        real64 & dDensity_dPressure,
                        real64 & viscosity,
                        real64 & dViscosity_dPressure ) const = 0;

  /**
   * @brief Compute fluid properties and derivatives at a single point.
   * @param[in]  pressure the target pressure value
   * @param[in]  temperature the target temperature value
   * @param[out] density fluid density
   * @param[out] dDensity_dPressure fluid density derivative w.r.t. pressure
   * @param[out] dDensity_dTemperature fluid density derivative w.r.t. temperature
   * @param[out] viscosity fluid viscosity
   * @param[out] dViscosity_dPressure fluid viscosity derivative w.r.t. pressure
   * @param[out] dViscosity_dTemperature fluid viscosity derivative w.r.t. temperature
   * @param[out] internalEnergy fluid internal energy
   * @param[out] dInternalEnergy_dPressure fluid internal energy derivative w.r.t. pressure
   * @param[out] dInternalEnergy_dTemperature fluid internal energy derivative w.r.t. temperature
   * @param[out] enthalpy fluid enthalpy
   * @param[out] dEnthalpy_dPressure fluid enthalpy derivative w.r.t. pressure
   * @param[out] dEnthalpy_dTemperature fluid enthalpy derivative w.r.t. temperature
   */
  GEOS_HOST_DEVICE
  virtual void compute( real64 const pressure,
                        real64 const temperature,
                        real64 & density,
                        real64 & dDensity_dPressure,
                        real64 & dDensity_dTemperature,
                        real64 & viscosity,
                        real64 & dViscosity_dPressure,
                        real64 & dViscosity_dTemperature,
                        real64 & internalEnergy,
                        real64 & dInternalEnergy_dPressure,
                        real64 & dInternalEnergy_dTemperature,
                        real64 & enthalpy,
                        real64 & dEnthalpy_dPressure,
                        real64 & dEnthalpy_dTemperature ) const = 0;

  /**
   * @brief Update fluid state at a single point.
   * @param[in] k        element index
   * @param[in] q        gauss point index
   * @param[in] pressure the target pressure value
   */
  GEOS_HOST_DEVICE
  virtual void update( localIndex const k,
                       localIndex const q,
                       real64 const pressure ) const = 0;

  /**
   * @brief Update fluid state at a single point.
   * @param[in] k           element index
   * @param[in] q           gauss point index
   * @param[in] pressure    the target pressure value
   * @param[in] temperature the target temperature value
   */
  GEOS_HOST_DEVICE
  virtual void update( localIndex const k,
                       localIndex const q,
                       real64 const pressure,
                       real64 const temperature ) const = 0;

};

Compound models

Compound constitutive models are employed to mimic the behavior of a material that requires a combination of constitutive models linked together. These compound models do not hold any data. They serve only as an interface with the individual models that they couple.

Coupled Solids

CoupledSolid models are employed to represent porous materials that require both a mechanical behavior and constitutive laws that describe the dependency of porosity and permeability on the primary unknowns.

The base class CoupledSolidBase implements some basic behaviors and is used to access a generic CoupledSolid in a physics solver:

  geos::constitutive::CoupledSolidBase const & solid =
    getConstitutiveModel< geos::constitutive::CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );

Additionally, a template class defines a base CoupledSolid model templated on the types of solid, porosity, and permeability models:

template< typename SOLID_TYPE,
          typename PORO_TYPE,
          typename PERM_TYPE >
class CoupledSolid : public CoupledSolidBase

While physics solvers that need a porous material only interface with a compound model, this one has access to the standalone models needed:

protected:
  SOLID_TYPE const & getSolidModel() const
  { return this->getParent().template getGroup< SOLID_TYPE >( m_solidModelName ); }

  PORO_TYPE const & getPorosityModel() const
  { return this->getParent().template getGroup< PORO_TYPE >( m_porosityModelName ); }

  PERM_TYPE const & getPermModel() const
  { return this->getParent().template getGroup< PERM_TYPE >( m_permeabilityModelName ); }

There are two specializations of a CoupledSolid:

  • CompressibleSolid: this model is used whenever there is no need to define a full mechanical model, but only simple correlations that compute material properties (like porosity or permeability). This model assumes that the solid model is of type NullModel and is only templated on the types of porosity and permeability models.

  • PorousSolid: this model is used to represent a full porous material where the porosity and permeability models need to be aware of the mechanical response of the material.