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:
CoupledSolidBase const & porousSolid =
getConstitutiveModel< 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.