Computes the hysteresis order for use by the hysteretic capillary pressure and relative-permeability objects. More...

#include <PorousFlowHysteresisOrder.h>

Inheritance diagram for PorousFlowHysteresisOrder:

Public Types
typedef DerivativeMaterialPropertyNameInterface::SymbolName	SymbolName

Public Member Functions
	PorousFlowHysteresisOrder (const InputParameters &parameters)

virtual void	initQpStatefulProperties () override

virtual void	computeQpProperties () override

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

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

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="")

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 Attributes
const unsigned	_liquid_ph_num
	Liquid phase number.

const std::string	_liquid_phase
	Stringified liquid phase number.

const unsigned	_initial_order
	Initial order.

const std::vector< Real >	_previous_turning_points
	Previous turning points that were encountered prior to the simulation.

const MaterialProperty< std::vector< Real > > &	_sat_old
	Old value of saturation.

const MaterialProperty< std::vector< Real > > &	_sat_older
	Older value of saturation.

MaterialProperty< unsigned > &	_hys_order
	Computed hysteresis order at the nodes or quadpoints.

const MaterialProperty< unsigned > &	_hys_order_old
	Old value of hysteresis order at the nodes or quadpoints.

MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > &	_hys_sat_tps
	Recorded saturation values at the turning points.

const MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > &	_hys_sat_tps_old
	Old value of recorded saturation values at the turning points.

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

Computes the hysteresis order for use by the hysteretic capillary pressure and relative-permeability objects.

Definition at line 20 of file PorousFlowHysteresisOrder.h.

Constructor & Destructor Documentation

◆ PorousFlowHysteresisOrder()

PorousFlowHysteresisOrder::PorousFlowHysteresisOrder ( const InputParameters & parameters )

Definition at line 41 of file PorousFlowHysteresisOrder.C.

  : DerivativeMaterialInterface<PorousFlowMaterial>(parameters),
    _liquid_ph_num(getParam<unsigned>("liquid_phase")),
    _liquid_phase(Moose::stringify(_liquid_ph_num)),
    _initial_order(getParam<unsigned>("initial_order")),
    _previous_turning_points(getParam<std::vector<Real>>("previous_turning_points")),
    _sat_old(_nodal_material
                 ? getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")
                 : getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_qp")),
    _sat_older(_nodal_material
                   ? getMaterialPropertyOlder<std::vector<Real>>("PorousFlow_saturation_nodal")
                   : getMaterialPropertyOlder<std::vector<Real>>("PorousFlow_saturation_qp")),
    _hys_order(_nodal_material ? declareProperty<unsigned>("PorousFlow_hysteresis_order_nodal")
                               : declareProperty<unsigned>("PorousFlow_hysteresis_order_qp")),
    _hys_order_old(_nodal_material
                       ? getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_nodal")
                       : getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_qp")),
    _hys_sat_tps(_nodal_material
                     ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
                           "PorousFlow_hysteresis_saturation_tps_nodal")
                     : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
                           "PorousFlow_hysteresis_saturation_tps_qp")),
    _hys_sat_tps_old(
        _nodal_material
            ? getMaterialPropertyOld<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
                  "PorousFlow_hysteresis_saturation_tps_nodal")
            : getMaterialPropertyOld<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
                  "PorousFlow_hysteresis_saturation_tps_qp"))
{
  if (_liquid_ph_num >= _dictator.numPhases())
    paramError("liquid_phase",
               "The Dictator proclaims that the number of fluid phases is ",
               _dictator.numPhases(),
               " while you have foolishly entered ",
               _liquid_ph_num,
               " for your liquid phase number.  Remember that indexing starts at 0.  Be aware that "
               "the Dictator does not tolerate mistakes.");
  if (_initial_order > PorousFlowConstants::MAX_HYSTERESIS_ORDER)
    paramError("initial_order",
               "The initial order must be less than the max order, which is hard-coded to ",
               PorousFlowConstants::MAX_HYSTERESIS_ORDER);
  if (_previous_turning_points.size() != _initial_order)
    paramError("previous_turning_points",
               "initial_order is ",
               _initial_order,
               " so there must be this many previous_turning_points");
  for (const auto & tp : _previous_turning_points)
    if (tp < 0.0 || tp > 1.0)
      paramError("previous_turning_points",
                 "Each previous_turning_point must lie in the range [0.0, 1.0]");
  // check that the previous_turning_points are ordered correctly
  Real previous_tp = 1.0;
  Real previous_jump = std::numeric_limits<Real>::max();
  bool drying = true;
  for (unsigned i = 0; i < _initial_order; ++i)
  {
    const Real this_jump = _previous_turning_points[i] - previous_tp;
    if (drying)
    {
      if (this_jump >= 0)
        paramError("previous_turning_points",
                   "The previous_turning_points vector, {a, b, c, ...}, must be obey the following "
                   "order: a < 1; b > a; a < c < b, etc");
    }
    else
    {
      if (this_jump <= 0)
        paramError("previous_turning_points",
                   "The previous_turning_points vector, {a, b, c, ...}, must be obey the following "
                   "order: a < 1; b > a; a < c < b, etc");
    }
    if (std::abs(this_jump) >= previous_jump)
      paramError("previous_turning_points",
                 "The previous_turning_points vector, {a, b, c, ...}, must be obey the following "
                 "order: a < 1; b > a; a < c < b, etc");
    drying = !drying;
    previous_tp = _previous_turning_points[i];
    previous_jump = std::abs(this_jump);
  }
}

Member Function Documentation

◆ computeQpProperties()

void PorousFlowHysteresisOrder::computeQpProperties ( )

overridevirtual

Definition at line 131 of file PorousFlowHysteresisOrder.C.

{
  // NOTE to readers: the computed properties depend on old and older quantities only.  Hence, it
  // would be nice if we could compute them only at the start of each timestep instead of every
  // nonlinear iteration
 
  // For finite elements involved in boundary-condition calculations, number of nodes > number of
  // quadpoints.  Material values on these elements are not computed during initial setup, which
  // means that for timestep = 1, _sat_older.size() = number of quadpoints.  Since
  // computeQpProperties() is computed with 0 <= _qp < number of nodes (for nodal materials), this
  // results in an out-of-bounds error on _sat_older.  PorousFlow's implementation of hysteresis
  // lags one timestep behind, and assumes that _sat_older = _sat_old at timestep = 1 anyway, so
  // just use it:
  const MaterialProperty<std::vector<Real>> & sat_older =
      (_nodal_material && _sat_old.size() != _sat_older.size()) ? _sat_old : _sat_older;
 
  // whether saturation has been reducing ("drying" or "draining"):
  const bool drying = (_hys_order_old[_qp] % 2 == 0);
  // whether have been drying but are now increasing saturation ("wetting" or "imbibing"):
  const bool dry2wet = drying && (_sat_old[_qp][_liquid_ph_num] > sat_older[_qp][_liquid_ph_num]);
  // whether have been wetting but are now decreasing saturation:
  const bool wet2dry = !drying && (_sat_old[_qp][_liquid_ph_num] < sat_older[_qp][_liquid_ph_num]);
  if ((dry2wet || wet2dry) && _hys_order_old[_qp] < PorousFlowConstants::MAX_HYSTERESIS_ORDER)
  {
    _hys_order[_qp] = _hys_order_old[_qp] + 1;
    _hys_sat_tps[_qp] = _hys_sat_tps_old[_qp];
    _hys_sat_tps[_qp].at(_hys_order_old[_qp]) = sat_older[_qp][_liquid_ph_num];
  }
  else
  {
    // no change in wetting/drying direction
    _hys_order[_qp] = _hys_order_old[_qp];
    _hys_sat_tps[_qp] = _hys_sat_tps_old[_qp];
  }
  if (_sat_old[_qp][_liquid_ph_num] >= 1.0)
    _hys_order[_qp] = 0; // assume the system gets back to the original drying curve
 
  // reduce order by 1 if _hys_order = max and saturation exceeds bounds of last turning-point
  if (_hys_order[_qp] == PorousFlowConstants::MAX_HYSTERESIS_ORDER)
  {
    // the maxiumum-order curve is a drying curve, and the saturation exceeds the previous turning
    // point:
    const bool drying_and_high_sat =
        (_hys_order[_qp] % 2 == 0) &&
        (_sat_old[_qp][_liquid_ph_num] >= _hys_sat_tps[_qp].at(_hys_order[_qp] - 1));
    // the maxiumum-order curve is a wetting curve, and the saturation is less than the previous
    // turning point:
    const bool wetting_and_low_sat =
        (_hys_order[_qp] % 2 == 1) &&
        (_sat_old[_qp][_liquid_ph_num] <= _hys_sat_tps[_qp].at(_hys_order[_qp] - 1));
    if (drying_and_high_sat || wetting_and_low_sat)
      _hys_order[_qp] -= 1;
  }
 
  // reduce order by 2 if saturation exceeds the bounds of the second-to-last turning-point
  // encountered
  bool can_reduce_order = (_hys_order[_qp] > 1);
  while (can_reduce_order)
  {
    can_reduce_order = false;
    // are drying and the saturation is <= a previously-encountered turning point:
    const bool drying_and_low_sat =
        (_hys_order[_qp] % 2 == 0) &&
        (_sat_old[_qp][_liquid_ph_num] <= _hys_sat_tps[_qp].at(_hys_order[_qp] - 2));
    // are wetting and the saturation is >= a previously-encountered turning point:
    const bool wetting_and_high_sat =
        (_hys_order[_qp] % 2 == 1) &&
        (_sat_old[_qp][_liquid_ph_num] >= _hys_sat_tps[_qp].at(_hys_order[_qp] - 2));
    if (drying_and_low_sat || wetting_and_high_sat)
    {
      _hys_order[_qp] -= 2;
      can_reduce_order = (_hys_order[_qp] > 1);
    }
  }
}

◆ initQpStatefulProperties()

void PorousFlowHysteresisOrder::initQpStatefulProperties ( )

overridevirtual

Definition at line 123 of file PorousFlowHysteresisOrder.C.

{
  _hys_order[_qp] = _initial_order;
  for (unsigned i = 0; i < _initial_order; ++i)
    _hys_sat_tps[_qp].at(i) = _previous_turning_points[i];
}

◆ validParams()

InputParameters PorousFlowHysteresisOrder::validParams ( )

static

Definition at line 16 of file PorousFlowHysteresisOrder.C.

{
  InputParameters params = PorousFlowMaterial::validParams();
  params.addPrivateParam<std::string>("pf_material_type", "hysteresis_order");
  params.addParam<unsigned>("liquid_phase", 0, "The phase number of the liquid phase");
  params.addParam<unsigned>(
      "initial_order",
      0,
      "The initial order.  0 means the simulation will begin on the primary drying curve.  1 means "
      "the simulation will begin on the first-order wetting curve.  2 means the simulation will "
      "begin on the second-order drying curve, etc.");
  params.addParam<std::vector<Real>>(
      "previous_turning_points",
      std::vector<Real>(),
      "The turning points (liquid saturation values) at occured prior to the simulation.  There "
      "must be exactly initial_order of these values.  The first value is the liquid saturation at "
      "the transition from primary-drying to first-order-wetting; the second value is the liquid "
      "saturation at the transition from first-order-wetting to second-order-drying; and so on.  "
      "You must ensure that your initial saturation field is commensurate with the "
      "previous_turning_points.");
  params.addClassDescription("Computes hysteresis order for use in hysteretic capillary pressures "
                             "and relative permeabilities");
  return params;
}