Base class Material designed to provide the porosity. More...

#include <PorousFlowPorosityExponentialBase.h>

Inheritance diagram for PorousFlowPorosityExponentialBaseTempl< is_ad >:

Public Types
typedef DerivativeMaterialPropertyNameInterface::SymbolName	SymbolName

Public Member Functions
	PorousFlowPorosityExponentialBaseTempl (const InputParameters &parameters)

const GenericMaterialProperty< U, is_ad > &	getDefaultMaterialProperty (const std::string &name)

const GenericMaterialProperty< U, is_ad > &	getDefaultMaterialPropertyByName (const std::string &name)

Moose::Kokkos::MaterialProperty< U, dimension >	declareKokkosPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

void	validateDerivativeMaterialPropertyBase (const std::string &base)

const MaterialPropertyName	derivativePropertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const

const MaterialPropertyName	derivativePropertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const

const MaterialPropertyName	derivativePropertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const

const MaterialPropertyName	derivativePropertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, unsigned int v2, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, unsigned int v1, unsigned int v2=libMesh::invalid_uint, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, unsigned int v2, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, unsigned int v1, unsigned int v2=libMesh::invalid_uint, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

Moose::Kokkos::MaterialProperty< U, dimension >	getKokkosMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

Moose::Kokkos::MaterialProperty< U, dimension >	getKokkosMaterialPropertyDerivative (const std::string &base, const SymbolName &c1)

Moose::Kokkos::MaterialProperty< U, dimension >	getKokkosMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

Moose::Kokkos::MaterialProperty< U, dimension >	getKokkosMaterialPropertyDerivative (const std::string &base, const SymbolName &c1)

Moose::Kokkos::MaterialProperty< U, dimension >	getKokkosMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)

Moose::Kokkos::MaterialProperty< U, dimension >	getKokkosMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)

void	validateCoupling (const MaterialPropertyName &base, const std::vector< VariableName > &c, bool validate_aux=true)

void	validateCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

void	validateCoupling (const MaterialPropertyName &base, const std::vector< VariableName > &c, bool validate_aux=true)

void	validateCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

void	validateNonlinearCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

void	validateNonlinearCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

const MaterialPropertyName	propertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const

const MaterialPropertyName	propertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const

const MaterialPropertyName	propertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const

const MaterialPropertyName	propertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const

const MaterialPropertyName	propertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const

const MaterialPropertyName	propertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const

const MaterialPropertyName	propertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const

const MaterialPropertyName	propertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const

Static Public Member Functions
static InputParameters	validParams ()

Protected Member Functions
virtual void	initQpStatefulProperties () override

virtual void	computeQpProperties () override

virtual GenericReal< is_ad >	atNegInfinityQp () const =0
	Returns "a" at the quadpoint (porosity = a + (b - a) * exp(decay))

virtual Real	datNegInfinityQp (unsigned pvar) const =0
	d(a)/d(PorousFlow variable pvar)

virtual GenericReal< is_ad >	atZeroQp () const =0
	Returns "b" at the quadpoint (porosity = a + (b - a) * exp(decay))

virtual Real	datZeroQp (unsigned pvar) const =0
	d(b)/d(PorousFlow variable pvar)

virtual GenericReal< is_ad >	decayQp () const =0
	Returns "decay" at the quadpoint (porosity = a + (b - a) * exp(decay))

virtual Real	ddecayQp_dvar (unsigned pvar) const =0
	d(decay)/d(PorousFlow variable pvar)

virtual RealGradient	ddecayQp_dgradvar (unsigned pvar) const =0
	d(decay)/d(grad(PorousFlow variable pvar))

void	applyPorosityMin ()
	If _porosity[_qp] < porosity_min, set _porosity[_qp] = porosity_min and multiply the derivatives of porosity by _zero_modifier (the derivative material properties for non-AD objects, and the AD derivatives for AD objects).

Protected Attributes
const bool	_strain_at_nearest_qp
	When calculating nodal porosity, use the strain at the nearest quadpoint to the node.

const bool	_ensure_positive
	for decayQp() > 0, porosity can be negative when using porosity = a + (b - a) * exp(decay).

	usingPorousFlowPorosityBaseMembers

GenericMaterialProperty< Real, is_ad > &	_porosity
	Computed porosity at the nodes or quadpoints.

MaterialProperty< std::vector< Real > > *const	_dporosity_dvar
	d(porosity)/d(PorousFlow variable)

MaterialProperty< std::vector< RealGradient > > *const	_dporosity_dgradvar
	d(porosity)/d(grad PorousFlow variable)

const Real	_porosity_min
	Minimum allowed porosity (default zero, so porosity is never negative)

const Real	_zero_modifier
	When the porosity_min floor is active, the porosity derivatives are multiplied by this.

const unsigned int	_num_phases
	Number of phases.

const unsigned int	_num_components
	Number of fluid components.

const unsigned int	_num_var
	Number of PorousFlow variables.

Private Member Functions
bool	haveMaterialProperty (const std::string &prop_name)

std::vector< VariableName >	buildVariableVector (const VariableName &c1, const VariableName &c2, const VariableName &c3)

void	validateCouplingHelper (const MaterialPropertyName &base, const std::vector< VariableName > &c, const System &system, std::vector< VariableName > &missing)

bool	isNotObjectVariable (const VariableName &name)

Private Attributes
FEProblemBase &	_dmi_fe_problem

Detailed Description

template<bool is_ad>
class PorousFlowPorosityExponentialBaseTempl< is_ad >

Base class Material designed to provide the porosity.

In this class porosity = a + (b - a) * exp(decay) where a = atNegInfinityQp() b = atZeroQp() decay = decayQp() Since this expression can become negative for decay > 0, if ensure_positive = true then for decay > 0 the following expression is used instead: porosity = a + (b - a) * exp(c * (1 - Exp(- decay / c))) where c = log(a/(b-a)) This latter expression is C1 continuous at decay=0 with the former expression. It is monotonically decreasing with "decay" and is positive.

Definition at line 32 of file PorousFlowPorosityExponentialBase.h.

Constructor & Destructor Documentation

◆ PorousFlowPorosityExponentialBaseTempl()

template<bool is_ad>

PorousFlowPorosityExponentialBaseTempl< is_ad >::PorousFlowPorosityExponentialBaseTempl ( const InputParameters & parameters )

Definition at line 45 of file PorousFlowPorosityExponentialBase.C.

  : PorousFlowPorosityBaseTempl<is_ad>(parameters),
    _strain_at_nearest_qp(this->template getParam<bool>("strain_at_nearest_qp")),
    _ensure_positive(this->template getParam<bool>("ensure_positive"))
{
}

Member Function Documentation

◆ applyPorosityMin()

template<bool is_ad>

void PorousFlowPorosityBaseTempl< is_ad >::applyPorosityMin ( )

protectedinherited

If _porosity[_qp] < porosity_min, set _porosity[_qp] = porosity_min and multiply the derivatives of porosity by _zero_modifier (the derivative material properties for non-AD objects, and the AD derivatives for AD objects).

Otherwise do nothing. Derived classes call this at the end of computeQpProperties(), after the unfloored porosity and its derivatives have been computed.

Definition at line 58 of file PorousFlowPorosityBase.C.

{
  if (MetaPhysicL::raw_value(_porosity[_qp]) >= _porosity_min)
    return;
 
  if constexpr (!is_ad)
  {
    _porosity[_qp] = _porosity_min;
    for (unsigned int v = 0; v < _num_var; ++v)
    {
      (*_dporosity_dvar)[_qp][v] *= _zero_modifier;
      (*_dporosity_dgradvar)[_qp][v] *= _zero_modifier;
    }
  }
  else
  {
    // The AD path carries its derivatives inside the value and has no derivative material
    // properties, so the value and the derivatives are softened separately.  Assigning
    // _porosity_min to the ADReal would instead discard the derivatives entirely.
    _porosity[_qp].value() = _porosity_min;
    _porosity[_qp].derivatives() *= _zero_modifier;
  }
}

Referenced by PorousFlowPorosityHMBiotModulus::computeQpProperties(), and PorousFlowPorosityLinear::computeQpProperties().

◆ atNegInfinityQp()

template<bool is_ad>

virtual GenericReal< is_ad > PorousFlowPorosityExponentialBaseTempl< is_ad >::atNegInfinityQp ( ) const

protectedpure virtual

Returns "a" at the quadpoint (porosity = a + (b - a) * exp(decay))

Implemented in PorousFlowPorosityTempl< is_ad >.

◆ atZeroQp()

template<bool is_ad>

virtual GenericReal< is_ad > PorousFlowPorosityExponentialBaseTempl< is_ad >::atZeroQp ( ) const

protectedpure virtual

Returns "b" at the quadpoint (porosity = a + (b - a) * exp(decay))

Implemented in PorousFlowPorosityTempl< is_ad >.

◆ computeQpProperties()

template<bool is_ad>

void PorousFlowPorosityExponentialBaseTempl< is_ad >::computeQpProperties ( )

overrideprotectedvirtual

Reimplemented in PorousFlowPorosityHMBiotModulus.

Definition at line 81 of file PorousFlowPorosityExponentialBase.C.

{
  // Unqualified calls so ADL finds the MetaPhysicL overloads for the AD (DualNumber) instantiation
  using std::exp;
  using std::log;
  using std::pow;
 
  const GenericReal<is_ad> a = atNegInfinityQp();
  const GenericReal<is_ad> b = atZeroQp();
  const GenericReal<is_ad> decay = decayQp();
  GenericReal<is_ad> exp_term = 1.0; // set appropriately below
 
  if (decay <= 0.0 || !_ensure_positive)
  {
    exp_term = exp(decay);
    _porosity[_qp] = a + (b - a) * exp_term;
  }
  else
  {
    const GenericReal<is_ad> c = log(a / (a - b));
    // note that at decay = 0, we have expx = 1, so porosity = a + b - a = b
    // and at decay = infinity, expx = 0, so porosity = a + (b - a) * a / (a - b) = 0
    const GenericReal<is_ad> expx = exp(-decay / c);
    exp_term = exp(c * (1.0 - expx));
    _porosity[_qp] = a + (b - a) * exp_term;
  }
 
  if constexpr (!is_ad)
  {
    // The AD path carries the derivatives with respect to the PorousFlow variables inside the AD
    // porosity value, so only the non-AD path builds these derivative material properties.
    Real deriv = 0.0; // = d(porosity)/d(decay)
    if (decay <= 0.0 || !_ensure_positive)
      deriv = _porosity[_qp] - a;
    else
    {
      const Real c = log(a / (a - b));
      const Real expx = exp(-decay / c);
      deriv = (_porosity[_qp] - a) * expx;
    }
 
    (*_dporosity_dvar)[_qp].resize(_num_var);
    (*_dporosity_dgradvar)[_qp].resize(_num_var);
    for (unsigned int v = 0; v < _num_var; ++v)
    {
      (*_dporosity_dvar)[_qp][v] = ddecayQp_dvar(v) * deriv;
      (*_dporosity_dgradvar)[_qp][v] = ddecayQp_dgradvar(v) * deriv;
 
      const Real da = datNegInfinityQp(v);
      const Real db = datZeroQp(v);
      (*_dporosity_dvar)[_qp][v] += da * (1 - exp_term) + db * exp_term;
 
      if (!(decay <= 0.0 || !_ensure_positive))
      {
        const Real c = log(a / (a - b));
        const Real expx = exp(-decay / c);
        const Real dc = (a - b) * (da * b / a - db) / Utility::pow<2>(a);
        (*_dporosity_dvar)[_qp][v] += (b - a) * exp_term * dc * (1 - expx - expx / c);
      }
    }
  }
 
  // Apply the porosity floor last, after the unfloored derivatives above have been
  // formed (see PorousFlowPorosityBaseTempl::applyPorosityMin)
  applyPorosityMin();
}

◆ datNegInfinityQp()

template<bool is_ad>

virtual Real PorousFlowPorosityExponentialBaseTempl< is_ad >::datNegInfinityQp ( unsigned pvar ) const

protectedpure virtual

d(a)/d(PorousFlow variable pvar)

Implemented in PorousFlowPorosityTempl< is_ad >.

◆ datZeroQp()

template<bool is_ad>

virtual Real PorousFlowPorosityExponentialBaseTempl< is_ad >::datZeroQp ( unsigned pvar ) const

protectedpure virtual

d(b)/d(PorousFlow variable pvar)

Implemented in PorousFlowPorosityTempl< is_ad >.

◆ ddecayQp_dgradvar()

template<bool is_ad>

virtual RealGradient PorousFlowPorosityExponentialBaseTempl< is_ad >::ddecayQp_dgradvar ( unsigned pvar ) const

protectedpure virtual

d(decay)/d(grad(PorousFlow variable pvar))

Implemented in PorousFlowPorosityTempl< is_ad >.

◆ ddecayQp_dvar()

template<bool is_ad>

virtual Real PorousFlowPorosityExponentialBaseTempl< is_ad >::ddecayQp_dvar ( unsigned pvar ) const

protectedpure virtual

d(decay)/d(PorousFlow variable pvar)

Implemented in PorousFlowPorosityTempl< is_ad >.

◆ decayQp()

template<bool is_ad>

virtual GenericReal< is_ad > PorousFlowPorosityExponentialBaseTempl< is_ad >::decayQp ( ) const

protectedpure virtual

Returns "decay" at the quadpoint (porosity = a + (b - a) * exp(decay))

Implemented in PorousFlowPorosityTempl< is_ad >.

◆ initQpStatefulProperties()

template<bool is_ad>

void PorousFlowPorosityExponentialBaseTempl< is_ad >::initQpStatefulProperties ( )

overrideprotectedvirtual

Definition at line 55 of file PorousFlowPorosityExponentialBase.C.

{
  // Unqualified calls so ADL finds the MetaPhysicL overloads for the AD (DualNumber) instantiation
  using std::exp;
  using std::log;
 
  const GenericReal<is_ad> a = atNegInfinityQp();
  const GenericReal<is_ad> b = atZeroQp();
  mooseAssert(a > b, "PorousFlowPorosityExponentialBase a must be larger than b");
  const GenericReal<is_ad> decay = decayQp();
 
  if (decay <= 0.0 || !_ensure_positive)
    _porosity[_qp] = a + (b - a) * exp(decay);
  else
  {
    const GenericReal<is_ad> c = log(a / (a - b));
    const GenericReal<is_ad> expx = exp(-decay / c);
    _porosity[_qp] = a + (b - a) * exp(c * (1.0 - expx));
  }
 
  if (_porosity[_qp] < _porosity_min)
    _porosity[_qp] = _porosity_min;
}

◆ validParams()

template<bool is_ad>

InputParameters PorousFlowPorosityExponentialBaseTempl< is_ad >::validParams ( )

static

Definition at line 16 of file PorousFlowPorosityExponentialBase.C.

{
  InputParameters params = PorousFlowPorosityBaseTempl<is_ad>::validParams();
  params.addParam<bool>("strain_at_nearest_qp",
                        false,
                        "When calculating nodal porosity that depends on strain, use the strain at "
                        "the nearest quadpoint.  This adds a small extra computational burden, and "
                        "is not necessary for simulations involving only linear lagrange elements. "
                        " If you set this to true, you will also want to set the same parameter to "
                        "true for related Kernels and Materials");
  params.addParam<bool>("ensure_positive",
                        true,
                        "Modify the usual exponential relationships that "
                        "governs porosity so that porosity is always "
                        "positive");
  params.setDocString(
      "porosity_min",
      "Minimum allowed value of the porosity: if the computed porosity is less than this value, "
      "porosity is set to this value instead.  The default of zero means porosity is never "
      "negative.  The ensure_positive transform only acts for decay > 0, so chemistry-driven "
      "(precipitation) porosity would otherwise become negative once pore space is filled by "
      "mineral; a small positive value may be preferable in that case.");
  params.addClassDescription("Base class Material for porosity that is computed via an exponential "
                             "relationship with coupled variables (strain, porepressure, "
                             "temperature, chemistry)");
  return params;
}