https://mooseframework.inl.gov
Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
PorousFlow2PhaseHysPP Class Reference

Material designed to calculate the 2 porepressures and 2 saturations, as well as derivatives of them, as well as capillary pressure, in two-phase situations with hysteretic capillary pressure, assuming the phase porepressures as the nonlinear variables. More...

#include <PorousFlow2PhaseHysPP.h>

Inheritance diagram for PorousFlow2PhaseHysPP:
[legend]

Public Types

typedef DerivativeMaterialPropertyNameInterface::SymbolName SymbolName
 

Public Member Functions

 PorousFlow2PhaseHysPP (const InputParameters &parameters)
 
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 Member Functions

virtual void initQpStatefulProperties () override
 
virtual void computeQpProperties () override
 
void buildQpPPSS ()
 Assemble std::vectors of porepressure and saturation at the nodes and quadpoints. More...
 
Real landSat (Real slDel) const
 
void computeTurningPointInfo (unsigned tp_num, Real tp_sat, Real tp_pc)
 Compute all relevant quantities at the given turning point. More...
 
Real capillaryPressureQp (Real sat) const
 
Real dcapillaryPressureQp (Real sat) const
 
Real d2capillaryPressureQp (Real sat) const
 
Real liquidSaturationQp (Real pc) const
 
Real dliquidSaturationQp (Real pc) const
 
Real d2liquidSaturationQp (Real pc) const
 

Protected Attributes

MaterialProperty< Real > & _pc
 Computed nodal or quadpoint values of capillary pressure. More...
 
const VariableValue_phase0_porepressure
 Nodal or quadpoint value of porepressure of the zero phase (eg, the water phase) More...
 
const VariableGradient_phase0_gradp_qp
 Gradient(phase0_porepressure) at the qps. More...
 
const unsigned int _phase0_porepressure_varnum
 Moose variable number of the phase0 porepressure. More...
 
const unsigned int _p0var
 PorousFlow variable number of the phase0 porepressure. More...
 
const VariableValue_phase1_porepressure
 Nodal or quadpoint value of porepressure of the one phase (eg, the gas phase) More...
 
const VariableGradient_phase1_gradp_qp
 Gradient(phase1_porepressure) at the qps. More...
 
const unsigned int _phase1_porepressure_varnum
 Moose variable number of the phase1 porepressure. More...
 
const unsigned int _p1var
 PorousFlow variable number of the phase1 porepressure. More...
 
const Real _alpha_d
 van Genuchten alpha parameter for the primary drying curve More...
 
const Real _alpha_w
 van Genuchten alpha parameter for the primary wetting curve More...
 
const Real _n_d
 van Genuchten n parameter for the primary drying curve More...
 
const Real _n_w
 van Genuchten n parameter for the primary wetting curve More...
 
const Real _s_l_min
 Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity) More...
 
const Real _s_lr
 Liquid saturation below which the liquid relative permeability is zero. More...
 
const Real _s_gr_max
 Residual gas saturation: 1 - _s_gr_max is the maximum saturation for which the van Genuchten expression is valid for the wetting curve. More...
 
const Real _pc_max
 Maximum capillary pressure: for Pc above this value, a "lower" extension will be used. More...
 
const Real _high_ratio
 The high-saturation extension to the wetting will commence at _high_ratio * (1 - _s_gr_del) More...
 
const PorousFlowVanGenuchten::LowCapillaryPressureExtension::ExtensionStrategy _low_ext_type
 Type of low-saturation extension. More...
 
const Real _s_low_d
 Saturation on the primary drying curve where low-saturation extension commences. More...
 
const Real _dpc_low_d
 d(Pc)/dS on the primary drying curve at S = _s_low_d More...
 
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
 Parameters involved in the low-saturation extension of the primary drying curve. More...
 
const Real _s_low_w
 Saturation on the primary wetting curve where low-saturation extension commences. More...
 
const Real _dpc_low_w
 d(Pc)/dS on the primary wetting curve at S = _s_low_w More...
 
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_w
 Parameters involved in the low-saturation extension of the primary wetting curve. More...
 
const PorousFlowVanGenuchten::HighCapillaryPressureExtension::ExtensionStrategy _high_ext_type
 Type of high-saturation extension of the wetting curves. More...
 
const Real _s_high
 Saturation at the point of high-saturation extension. More...
 
const Real _pc_high
 Pc at the point of high-saturation extension. More...
 
const Real _dpc_high
 d(Pc)/d(S) at the point of high-saturation extension More...
 
const PorousFlowVanGenuchten::HighCapillaryPressureExtension _high_ext
 Parameters involved in the high-saturation extension of the primary wetting curve. More...
 
const MaterialProperty< unsigned > & _hys_order
 Hysteresis order, as computed by PorousFlowHysteresisOrder. More...
 
const MaterialProperty< unsigned > & _hys_order_old
 Old value of hysteresis order, as computed by PorousFlowHysteresisOrder. More...
 
const MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _hys_sat_tps
 Saturation values at the turning points, as computed by PorousFlowHysteresisOrder. More...
 
const MaterialProperty< Real > & _pc_older
 Older value of capillary pressure. More...
 
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _pc_tps
 Nodal or quadpoint values of Pc at the turning points. More...
 
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _s_d_tps
 Computed nodal or quadpoint values of saturation on the drying curve at _pc_tps. More...
 
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _s_gr_tps
 Computed nodal or quadpoint values of S_gr_Del, ie, the Land expression, at the turning points. More...
 
MaterialProperty< std::array< PorousFlowVanGenuchten::LowCapillaryPressureExtension, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _w_low_ext_tps
 Nodal or quadpoint values of the low extension of the wetting curve defined by _s_gr_tps. More...
 
MaterialProperty< std::array< PorousFlowVanGenuchten::HighCapillaryPressureExtension, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _w_high_ext_tps
 Nodal or quadpoint values of the high extension of the wetting curve defined by _s_gr_tps. More...
 
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _s_w_tps
 Computed nodal or quadpoint values of liquid saturation on the wetting curve defined by _s_gr_del, at pc = _pc_d_tps. More...
 
const unsigned int _num_phases
 Number of phases. More...
 
const unsigned int _num_components
 Number of components. More...
 
const unsigned int _num_pf_vars
 Number of PorousFlow variables. More...
 
GenericMaterialProperty< std::vector< Real >, is_ad > & _porepressure
 Computed nodal or quadpoint values of porepressure of the phases. More...
 
MaterialProperty< std::vector< std::vector< Real > > > *const _dporepressure_dvar
 d(porepressure)/d(PorousFlow variable) More...
 
GenericMaterialProperty< std::vector< RealGradient >, is_ad > *const _gradp_qp
 Grad(p) at the quadpoints. More...
 
MaterialProperty< std::vector< std::vector< Real > > > *const _dgradp_qp_dgradv
 d(grad porepressure)/d(grad PorousFlow variable) at the quadpoints More...
 
MaterialProperty< std::vector< std::vector< RealGradient > > > *const _dgradp_qp_dv
 d(grad porepressure)/d(PorousFlow variable) at the quadpoints More...
 
GenericMaterialProperty< std::vector< Real >, is_ad > & _saturation
 Computed nodal or qp saturation of the phases. More...
 
MaterialProperty< std::vector< std::vector< Real > > > *const _dsaturation_dvar
 d(saturation)/d(PorousFlow variable) More...
 
GenericMaterialProperty< std::vector< RealGradient >, is_ad > *const _grads_qp
 Grad(s) at the quadpoints. More...
 
MaterialProperty< std::vector< std::vector< Real > > > *const _dgrads_qp_dgradv
 d(grad saturation)/d(grad PorousFlow variable) at the quadpoints More...
 
MaterialProperty< std::vector< std::vector< RealGradient > > > *const _dgrads_qp_dv
 d(grad saturation)/d(PorousFlow variable) at the quadpoints More...
 

Detailed Description

Material designed to calculate the 2 porepressures and 2 saturations, as well as derivatives of them, as well as capillary pressure, in two-phase situations with hysteretic capillary pressure, assuming the phase porepressures as the nonlinear variables.

Definition at line 19 of file PorousFlow2PhaseHysPP.h.

Constructor & Destructor Documentation

◆ PorousFlow2PhaseHysPP()

PorousFlow2PhaseHysPP::PorousFlow2PhaseHysPP ( const InputParameters parameters)

Definition at line 32 of file PorousFlow2PhaseHysPP.C.

34  _pc(_nodal_material ? declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_nodal")
35  : declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_qp")),
36  _phase0_porepressure(_nodal_material ? coupledDofValues("phase0_porepressure")
37  : coupledValue("phase0_porepressure")),
38  _phase0_gradp_qp(coupledGradient("phase0_porepressure")),
39  _phase0_porepressure_varnum(coupled("phase0_porepressure")),
40  _p0var(_dictator.isPorousFlowVariable(_phase0_porepressure_varnum)
41  ? _dictator.porousFlowVariableNum(_phase0_porepressure_varnum)
42  : 0),
43 
44  _phase1_porepressure(_nodal_material ? coupledDofValues("phase1_porepressure")
45  : coupledValue("phase1_porepressure")),
46  _phase1_gradp_qp(coupledGradient("phase1_porepressure")),
47  _phase1_porepressure_varnum(coupled("phase1_porepressure")),
48  _p1var(_dictator.isPorousFlowVariable(_phase1_porepressure_varnum)
49  ? _dictator.porousFlowVariableNum(_phase1_porepressure_varnum)
50  : 0)
51 {
52  if (_num_phases != 2)
53  mooseError("The Dictator announces that the number of phases is ",
54  _dictator.numPhases(),
55  " whereas PorousFlow2PhaseHysPP can only be used for 2-phase simulation. When you "
56  "have efficient government, you have dictatorship.");
57 }
const VariableGradient & _phase0_gradp_qp
Gradient(phase0_porepressure) at the qps.
void mooseError(Args &&... args)
const unsigned int _phase0_porepressure_varnum
Moose variable number of the phase0 porepressure.
PorousFlowHystereticCapillaryPressure(const InputParameters &parameters)
const VariableValue & _phase1_porepressure
Nodal or quadpoint value of porepressure of the one phase (eg, the gas phase)
const VariableGradient & _phase1_gradp_qp
Gradient(phase1_porepressure) at the qps.
const unsigned int _num_phases
Number of phases.
const unsigned int _phase1_porepressure_varnum
Moose variable number of the phase1 porepressure.
MaterialProperty< Real > & _pc
Computed nodal or quadpoint values of capillary pressure.
const VariableValue & _phase0_porepressure
Nodal or quadpoint value of porepressure of the zero phase (eg, the water phase)
const unsigned int _p0var
PorousFlow variable number of the phase0 porepressure.
const unsigned int _p1var
PorousFlow variable number of the phase1 porepressure.

Member Function Documentation

◆ buildQpPPSS()

void PorousFlow2PhaseHysPP::buildQpPPSS ( )
protected

Assemble std::vectors of porepressure and saturation at the nodes and quadpoints.

Definition at line 132 of file PorousFlow2PhaseHysPP.C.

Referenced by computeQpProperties(), and initQpStatefulProperties().

133 {
134  _porepressure[_qp][0] = _phase0_porepressure[_qp];
135  _porepressure[_qp][1] = _phase1_porepressure[_qp];
136  _pc[_qp] = _phase1_porepressure[_qp] - _phase0_porepressure[_qp]; // this is >= 0
137  const Real sat = liquidSaturationQp(_pc[_qp]);
138  _saturation[_qp][0] = sat;
139  _saturation[_qp][1] = 1.0 - sat;
140 }
GenericMaterialProperty< std::vector< Real >, is_ad > & _porepressure
Computed nodal or quadpoint values of porepressure of the phases.
const VariableValue & _phase1_porepressure
Nodal or quadpoint value of porepressure of the one phase (eg, the gas phase)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MaterialProperty< Real > & _pc
Computed nodal or quadpoint values of capillary pressure.
const VariableValue & _phase0_porepressure
Nodal or quadpoint value of porepressure of the zero phase (eg, the water phase)
GenericMaterialProperty< std::vector< Real >, is_ad > & _saturation
Computed nodal or qp saturation of the phases.

◆ capillaryPressureQp()

Real PorousFlowHystereticCapillaryPressure::capillaryPressureQp ( Real  sat) const
protectedinherited
Returns
the value of capillary pressure, given the liquid saturation. This uses _hys_order[_qp]
Parameters
satliquid saturation

Definition at line 300 of file PorousFlowHystereticCapillaryPressure.C.

Referenced by PorousFlow2PhaseHysPS::buildQpPPSS(), PorousFlowHystereticInfo::computeQpInfo(), PorousFlowHystereticInfo::computeQpProperties(), and PorousFlowHystereticInfo::initQpStatefulProperties().

301 {
302  Real pc = 0.0;
303  if (_hys_order[_qp] == 0) // on primary drying curve
305  sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
306  else if (_hys_order[_qp] == 1) // first-order wetting
307  pc = firstOrderWettingPc(sat);
308  else if (_hys_order[_qp] == 2) // second-order drying
309  pc = secondOrderDryingPc(sat);
310  else // third order drying and wetting
311  {
312  const Real tp1 = _hys_sat_tps[_qp].at(1);
313  const Real tp2 = _hys_sat_tps[_qp].at(2);
314  const Real pc1 = firstOrderWettingPc(sat);
315  const Real pc2 = secondOrderDryingPc(sat);
316  // handle cases that occur just at the transition from 3rd to 2nd order, or 3rd to 1st order
317  if (sat < tp2)
318  pc = pc2;
319  else if (sat > tp1)
320  pc = pc1;
321  else if (pc1 >= pc2)
322  pc = pc1;
323  else if (pc1 <= 0.0 || pc2 <= 0.0)
324  pc = 0.0;
325  else
326  pc = std::exp(std::log(pc1) + (sat - tp1) * (std::log(pc2) - std::log(pc1)) / (tp2 - tp1));
327  }
328  return pc;
329 }
const MaterialProperty< unsigned > & _hys_order
Hysteresis order, as computed by PorousFlowHysteresisOrder.
const Real _alpha_d
van Genuchten alpha parameter for the primary drying curve
const Real _n_d
van Genuchten n parameter for the primary drying curve
Real capillaryPressureHys(Real sl, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
Hysteretic capillary pressure function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of Dou...
const MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _hys_sat_tps
Saturation values at the turning points, as computed by PorousFlowHysteresisOrder.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
Parameters involved in the low-saturation extension of the primary drying curve.
const Real _s_l_min
Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity) ...

◆ computeQpProperties()

void PorousFlow2PhaseHysPP::computeQpProperties ( )
overrideprotectedvirtual

Reimplemented from PorousFlowHystereticCapillaryPressure.

Definition at line 67 of file PorousFlow2PhaseHysPP.C.

68 {
69  // size stuff correctly and prepare the derivative matrices with zeroes
71 
72  buildQpPPSS();
73  const Real pc = _pc[_qp]; // >= 0
74  const Real ds = dliquidSaturationQp(pc); // dS/d(pc)
75 
76  if (!_nodal_material)
77  {
78  (*_gradp_qp)[_qp][0] = _phase0_gradp_qp[_qp];
79  (*_gradp_qp)[_qp][1] = _phase1_gradp_qp[_qp];
80  (*_grads_qp)[_qp][0] = ds * ((*_gradp_qp)[_qp][1] - (*_gradp_qp)[_qp][0]);
81  (*_grads_qp)[_qp][1] = -(*_grads_qp)[_qp][0];
82  }
83 
84  // the derivatives of porepressure with respect to porepressure
85  // remain fixed (at unity) throughout the simulation
86  if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
87  {
88  (*_dporepressure_dvar)[_qp][0][_p0var] = 1.0;
89  if (!_nodal_material)
90  (*_dgradp_qp_dgradv)[_qp][0][_p0var] = 1.0;
91  }
92  if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
93  {
94  (*_dporepressure_dvar)[_qp][1][_p1var] = 1.0;
95  if (!_nodal_material)
96  (*_dgradp_qp_dgradv)[_qp][1][_p1var] = 1.0;
97  }
98 
99  if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
100  {
101  (*_dsaturation_dvar)[_qp][0][_p0var] = -ds;
102  (*_dsaturation_dvar)[_qp][1][_p0var] = ds;
103  }
104  if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
105  {
106  (*_dsaturation_dvar)[_qp][0][_p1var] = ds;
107  (*_dsaturation_dvar)[_qp][1][_p1var] = -ds;
108  }
109 
110  _pc[_qp] = _phase1_porepressure[_qp] - _phase0_porepressure[_qp]; // this is >= 0
111  if (!_nodal_material)
112  {
113  const Real d2s_qp = d2liquidSaturationQp(pc); // d^2(S_qp)/d(pc_qp)^2
114  if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
115  {
116  (*_dgrads_qp_dgradv)[_qp][0][_p0var] = -ds;
117  (*_dgrads_qp_dv)[_qp][0][_p0var] = -d2s_qp * (_phase1_gradp_qp[_qp] - _phase0_gradp_qp[_qp]);
118  (*_dgrads_qp_dgradv)[_qp][1][_p0var] = ds;
119  (*_dgrads_qp_dv)[_qp][1][_p0var] = d2s_qp * (_phase1_gradp_qp[_qp] - _phase0_gradp_qp[_qp]);
120  }
121  if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
122  {
123  (*_dgrads_qp_dgradv)[_qp][0][_p1var] = ds;
124  (*_dgrads_qp_dv)[_qp][0][_p1var] = d2s_qp * (_phase1_gradp_qp[_qp] - _phase0_gradp_qp[_qp]);
125  (*_dgrads_qp_dgradv)[_qp][1][_p1var] = -ds;
126  (*_dgrads_qp_dv)[_qp][1][_p1var] = -d2s_qp * (_phase1_gradp_qp[_qp] - _phase0_gradp_qp[_qp]);
127  }
128  }
129 }
const VariableGradient & _phase0_gradp_qp
Gradient(phase0_porepressure) at the qps.
const unsigned int _phase0_porepressure_varnum
Moose variable number of the phase0 porepressure.
const VariableValue & _phase1_porepressure
Nodal or quadpoint value of porepressure of the one phase (eg, the gas phase)
const VariableGradient & _phase1_gradp_qp
Gradient(phase1_porepressure) at the qps.
void buildQpPPSS()
Assemble std::vectors of porepressure and saturation at the nodes and quadpoints. ...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const unsigned int _phase1_porepressure_varnum
Moose variable number of the phase1 porepressure.
MaterialProperty< Real > & _pc
Computed nodal or quadpoint values of capillary pressure.
const VariableValue & _phase0_porepressure
Nodal or quadpoint value of porepressure of the zero phase (eg, the water phase)
const unsigned int _p0var
PorousFlow variable number of the phase0 porepressure.
const unsigned int _p1var
PorousFlow variable number of the phase1 porepressure.

◆ computeTurningPointInfo()

void PorousFlowHystereticCapillaryPressure::computeTurningPointInfo ( unsigned  tp_num,
Real  tp_sat,
Real  tp_pc 
)
protectedinherited

Compute all relevant quantities at the given turning point.

Parameters
tp_numThe turning point number. Eg, tp_num=0 upon transition from zeroth-order drying to first-order wetting
tp_satLiquid saturation at the turning point
tp_pcCapillary pressure at the turning point

Definition at line 240 of file PorousFlowHystereticCapillaryPressure.C.

Referenced by PorousFlowHystereticCapillaryPressure::computeQpProperties(), and PorousFlowHystereticCapillaryPressure::initQpStatefulProperties().

243 {
244  _pc_tps[_qp].at(tp_num) = tp_pc;
245 
246  // Quantities on the drying curve:
247  // pc on the drying curve at the turning point saturation
249  tp_sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
250  // saturation on the drying curve, at tp_pc
251  _s_d_tps[_qp].at(tp_num) =
253 
254  // Quantities relevant to the wetting curve defined by the Land expression.
255  // s_gr_tps is the Land expression as a function of the turning point saturation
256  _s_gr_tps[_qp].at(tp_num) = landSat(tp_sat);
257  // the low extension of the wetting curve defined by _s_gr_tps
258  const Real s_w_low_ext = PorousFlowVanGenuchten::saturationHys(
259  _pc_max, _s_l_min, _s_gr_tps[_qp].at(tp_num), _alpha_w, _n_w);
261  s_w_low_ext, _s_l_min, _s_gr_tps[_qp].at(tp_num), _alpha_w, _n_w);
263  _low_ext_type, s_w_low_ext, _pc_max, dpc_w_low_ext);
264  // the high extension of the wetting curve defined by _s_gr_tps
265  const Real s_w_high_ext =
267  ? 1.0 - _s_gr_tps[_qp].at(tp_num)
268  : _high_ratio *
269  (1.0 -
270  _s_gr_tps[_qp].at(
271  tp_num)); // if NONE then use the vanGenuchten all the way to 1 - _s_gr_tps
272  const Real pc_w_high_ext =
274  _s_l_min,
275  _s_gr_tps[_qp].at(tp_num),
276  _alpha_w,
277  _n_w,
278  _w_low_ext_tps[_qp].at(tp_num));
279  const Real dpc_w_high_ext =
281  _s_l_min,
282  _s_gr_tps[_qp].at(tp_num),
283  _alpha_w,
284  _n_w,
285  _w_low_ext_tps[_qp].at(tp_num));
287  _high_ext_type, s_w_high_ext, pc_w_high_ext, dpc_w_high_ext);
288 
289  // saturation on the wetting curve defined by _s_gr_tps, at pc = pc_d_tps
290  _s_w_tps[_qp].at(tp_num) = PorousFlowVanGenuchten::saturationHys(pc_d_tps,
291  _s_l_min,
292  _s_gr_tps[_qp].at(tp_num),
293  _alpha_w,
294  _n_w,
295  _w_low_ext_tps[_qp].at(tp_num),
296  _w_high_ext_tps[_qp].at(tp_num));
297 }
Real dcapillaryPressureHys(Real sl, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
Derivative of capillaryPressureHys with respect to sl.
MaterialProperty< std::array< PorousFlowVanGenuchten::LowCapillaryPressureExtension, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _w_low_ext_tps
Nodal or quadpoint values of the low extension of the wetting curve defined by _s_gr_tps.
MaterialProperty< std::array< PorousFlowVanGenuchten::HighCapillaryPressureExtension, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _w_high_ext_tps
Nodal or quadpoint values of the high extension of the wetting curve defined by _s_gr_tps.
const PorousFlowVanGenuchten::LowCapillaryPressureExtension::ExtensionStrategy _low_ext_type
Type of low-saturation extension.
const Real _high_ratio
The high-saturation extension to the wetting will commence at _high_ratio * (1 - _s_gr_del) ...
const Real _pc_max
Maximum capillary pressure: for Pc above this value, a "lower" extension will be used.
Parameters associated with the extension of the hysteretic wetting capillary pressure function to hig...
const PorousFlowVanGenuchten::HighCapillaryPressureExtension::ExtensionStrategy _high_ext_type
Type of high-saturation extension of the wetting curves.
const Real _alpha_d
van Genuchten alpha parameter for the primary drying curve
const Real _n_d
van Genuchten n parameter for the primary drying curve
Real capillaryPressureHys(Real sl, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
Hysteretic capillary pressure function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of Dou...
Parameters associated with the extension of the hysteretic capillary pressure function to low saturat...
const Real _alpha_w
van Genuchten alpha parameter for the primary wetting curve
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _pc_tps
Nodal or quadpoint values of Pc at the turning points.
const Real _n_w
van Genuchten n parameter for the primary wetting curve
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
Parameters involved in the low-saturation extension of the primary drying curve.
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _s_gr_tps
Computed nodal or quadpoint values of S_gr_Del, ie, the Land expression, at the turning points...
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _s_w_tps
Computed nodal or quadpoint values of liquid saturation on the wetting curve defined by _s_gr_del...
const Real _s_l_min
Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity) ...
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _s_d_tps
Computed nodal or quadpoint values of saturation on the drying curve at _pc_tps.
Real saturationHys(Real pc, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
Hysteretic saturation function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of Doughty2008...

◆ d2capillaryPressureQp()

Real PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp ( Real  sat) const
protectedinherited
Returns
d^2(capillary pressure)/d(sat)^2. This uses _hys_order[_qp]
Parameters
satliquid saturation

Definition at line 371 of file PorousFlowHystereticCapillaryPressure.C.

Referenced by PorousFlowHystereticInfo::computeQpInfo(), and PorousFlow2PhaseHysPS::computeQpProperties().

372 {
373  Real d2pc = 0.0;
374  if (_hys_order[_qp] == 0) // on primary drying curve
376  sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
377  else if (_hys_order[_qp] == 1) // first-order wetting
378  d2pc = d2firstOrderWettingPc(sat);
379  else if (_hys_order[_qp] == 2) // second-order drying
380  d2pc = d2secondOrderDryingPc(sat);
381  else // third order drying and wetting
382  {
383  const Real tp1 = _hys_sat_tps[_qp].at(1);
384  const Real tp2 = _hys_sat_tps[_qp].at(2);
385  const Real pc1 = firstOrderWettingPc(sat);
386  const Real pc2 = secondOrderDryingPc(sat);
387  // handle cases that occur just at the transition from 3rd to 2nd order, or 3rd to 1st order
388  if (sat < tp2)
389  d2pc = d2secondOrderDryingPc(sat);
390  else if (sat > tp1)
391  d2pc = d2firstOrderWettingPc(sat);
392  else if (pc1 >= pc2)
393  d2pc = d2firstOrderWettingPc(sat);
394  else if (pc1 <= 0.0 || pc2 <= 0.0)
395  d2pc = 0.0;
396  else
397  {
398  const Real pc =
399  std::exp(std::log(pc1) + (sat - tp1) * (std::log(pc2) - std::log(pc1)) / (tp2 - tp1));
400  const Real dpc1 = dfirstOrderWettingPc(sat);
401  const Real dpc2 = dsecondOrderDryingPc(sat);
402  const Real dpc = pc * (dpc1 / pc1 + (std::log(pc2) - std::log(pc1)) / (tp2 - tp1) +
403  (sat - tp1) * (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1));
404  const Real d2pc1 = d2firstOrderWettingPc(sat);
405  const Real d2pc2 = d2secondOrderDryingPc(sat);
406  d2pc =
407  dpc * (dpc1 / pc1 + (std::log(pc2) - std::log(pc1)) / (tp2 - tp1) +
408  (sat - tp1) * (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1)) +
409  pc *
410  (d2pc1 / pc1 - dpc1 * dpc1 / pc1 / pc1 + (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1) +
411  (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1) +
412  (sat - tp1) *
413  (d2pc2 / pc2 - dpc2 * dpc2 / pc2 / pc2 - d2pc1 / pc1 + dpc1 * dpc1 / pc1 / pc1) /
414  (tp2 - tp1));
415  }
416  }
417  return d2pc;
418 }
Real d2capillaryPressureHys(Real sl, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
Second derivative of capillaryPressureHys with respect to sl.
const MaterialProperty< unsigned > & _hys_order
Hysteresis order, as computed by PorousFlowHysteresisOrder.
const Real _alpha_d
van Genuchten alpha parameter for the primary drying curve
const Real _n_d
van Genuchten n parameter for the primary drying curve
const MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _hys_sat_tps
Saturation values at the turning points, as computed by PorousFlowHysteresisOrder.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
Parameters involved in the low-saturation extension of the primary drying curve.
const Real _s_l_min
Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity) ...

◆ d2liquidSaturationQp()

Real PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp ( Real  pc) const
protectedinherited
Returns
d^2(liquid saturation)/d(pc)^2. This uses _hys_order[_qp]
Parameters
pccapillary pressure

Definition at line 698 of file PorousFlowHystereticCapillaryPressure.C.

Referenced by PorousFlowHystereticInfo::computeQpInfo(), computeQpProperties(), and PorousFlow1PhaseHysP::computeQpProperties().

699 {
700  Real d2sat = 0.0;
701  if (_hys_order[_qp] == 0) // on primary drying curve
703  else if (_hys_order[_qp] == 1) // first-order wetting
704  d2sat = d2firstOrderWettingSat(pc);
705  else if (_hys_order[_qp] == 2) // second-order drying
706  d2sat = d2secondOrderDryingSat(pc);
707  else // third order drying and wetting
708  {
709  const Real pc_tp1 = _pc_tps[_qp].at(1); // pc on first-order wetting at TP_1
710  const Real pc_tp2 = _pc_tps[_qp].at(2); // pc on second-order drying at TP_2
711  const Real sat1 = firstOrderWettingSat(pc);
712  const Real sat2 = secondOrderDryingSat(pc);
713  const Real dsat1 = dfirstOrderWettingSat(pc);
714  const Real dsat2 = dsecondOrderDryingSat(pc);
715  const Real d2sat1 = d2firstOrderWettingSat(pc);
716  const Real d2sat2 = d2secondOrderDryingSat(pc);
717  if (pc_tp1 <= 0.0 || pc <= 0.0 || pc_tp2 <= 0.0)
718  d2sat = d2sat2;
719  else if (pc > pc_tp2)
720  d2sat = d2sat2;
721  else if (pc < pc_tp1)
722  d2sat = d2sat1;
723  else
724  d2sat = d2sat1 - 1.0 / pc / pc * (sat2 - sat1) / std::log(pc_tp2 / pc_tp1) +
725  2.0 / pc * (dsat2 - dsat1) / std::log(pc_tp2 / pc_tp1) +
726  (std::log(pc) - std::log(pc_tp1)) * (d2sat2 - d2sat1) / std::log(pc_tp2 / pc_tp1);
727  }
728  return d2sat;
729 }
const MaterialProperty< unsigned > & _hys_order
Hysteresis order, as computed by PorousFlowHysteresisOrder.
const Real _alpha_d
van Genuchten alpha parameter for the primary drying curve
const Real _n_d
van Genuchten n parameter for the primary drying curve
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _pc_tps
Nodal or quadpoint values of Pc at the turning points.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
Parameters involved in the low-saturation extension of the primary drying curve.
Real d2saturationHys(Real pc, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
Second derivative of Hysteretic saturation function with respect to pc.
const Real _s_l_min
Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity) ...

◆ dcapillaryPressureQp()

Real PorousFlowHystereticCapillaryPressure::dcapillaryPressureQp ( Real  sat) const
protectedinherited
Returns
d(capillary pressure)/d(sat). This uses _hys_order[_qp]
Parameters
satliquid saturation

Definition at line 332 of file PorousFlowHystereticCapillaryPressure.C.

Referenced by PorousFlowHystereticInfo::computeQpInfo(), and PorousFlow2PhaseHysPS::computeQpProperties().

333 {
334  Real dpc = 0.0;
335  if (_hys_order[_qp] == 0) // on primary drying curve
337  sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
338  else if (_hys_order[_qp] == 1) // first-order wetting
339  dpc = dfirstOrderWettingPc(sat);
340  else if (_hys_order[_qp] == 2) // second-order drying
341  dpc = dsecondOrderDryingPc(sat);
342  else // third order drying and wetting
343  {
344  const Real tp1 = _hys_sat_tps[_qp].at(1);
345  const Real tp2 = _hys_sat_tps[_qp].at(2);
346  const Real pc1 = firstOrderWettingPc(sat);
347  const Real pc2 = secondOrderDryingPc(sat);
348  // handle cases that occur just at the transition from 3rd to 2nd order, or 3rd to 1st order
349  if (sat < tp2)
350  dpc = dsecondOrderDryingPc(sat);
351  else if (sat > tp1)
352  dpc = dfirstOrderWettingPc(sat);
353  else if (pc1 >= pc2)
354  dpc = dfirstOrderWettingPc(sat);
355  else if (pc1 <= 0.0 || pc2 <= 0.0)
356  dpc = 0.0;
357  else
358  {
359  const Real pc =
360  std::exp(std::log(pc1) + (sat - tp1) * (std::log(pc2) - std::log(pc1)) / (tp2 - tp1));
361  const Real dpc1 = dfirstOrderWettingPc(sat);
362  const Real dpc2 = dsecondOrderDryingPc(sat);
363  dpc = pc * (dpc1 / pc1 + (std::log(pc2) - std::log(pc1)) / (tp2 - tp1) +
364  (sat - tp1) * (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1));
365  }
366  }
367  return dpc;
368 }
Real dcapillaryPressureHys(Real sl, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
Derivative of capillaryPressureHys with respect to sl.
const MaterialProperty< unsigned > & _hys_order
Hysteresis order, as computed by PorousFlowHysteresisOrder.
const Real _alpha_d
van Genuchten alpha parameter for the primary drying curve
const Real _n_d
van Genuchten n parameter for the primary drying curve
const MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _hys_sat_tps
Saturation values at the turning points, as computed by PorousFlowHysteresisOrder.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
Parameters involved in the low-saturation extension of the primary drying curve.
const Real _s_l_min
Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity) ...

◆ dliquidSaturationQp()

Real PorousFlowHystereticCapillaryPressure::dliquidSaturationQp ( Real  pc) const
protectedinherited
Returns
d(liquid saturation)/d(pc). This uses _hys_order[_qp]
Parameters
pccapillary pressure

Definition at line 667 of file PorousFlowHystereticCapillaryPressure.C.

Referenced by PorousFlowHystereticInfo::computeQpInfo(), computeQpProperties(), and PorousFlow1PhaseHysP::computeQpProperties().

668 {
669  Real dsat = 0.0;
670  if (_hys_order[_qp] == 0) // on primary drying curve
672  else if (_hys_order[_qp] == 1) // first-order wetting
673  dsat = dfirstOrderWettingSat(pc);
674  else if (_hys_order[_qp] == 2) // second-order drying
675  dsat = dsecondOrderDryingSat(pc);
676  else // third order drying and wetting
677  {
678  const Real pc_tp1 = _pc_tps[_qp].at(1); // pc on first-order wetting at TP_1
679  const Real pc_tp2 = _pc_tps[_qp].at(2); // pc on second-order drying at TP_2
680  const Real sat1 = firstOrderWettingSat(pc);
681  const Real sat2 = secondOrderDryingSat(pc);
682  const Real dsat1 = dfirstOrderWettingSat(pc);
683  const Real dsat2 = dsecondOrderDryingSat(pc);
684  if (pc_tp1 <= 0.0 || pc <= 0.0 || pc_tp2 <= 0.0)
685  dsat = dsat2;
686  else if (pc > pc_tp2)
687  dsat = dsat2;
688  else if (pc < pc_tp1)
689  dsat = dsat1;
690  else
691  dsat = dsat1 + 1.0 / pc * (sat2 - sat1) / std::log(pc_tp2 / pc_tp1) +
692  (std::log(pc) - std::log(pc_tp1)) * (dsat2 - dsat1) / std::log(pc_tp2 / pc_tp1);
693  }
694  return dsat;
695 }
const MaterialProperty< unsigned > & _hys_order
Hysteresis order, as computed by PorousFlowHysteresisOrder.
const Real _alpha_d
van Genuchten alpha parameter for the primary drying curve
const Real _n_d
van Genuchten n parameter for the primary drying curve
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _pc_tps
Nodal or quadpoint values of Pc at the turning points.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
Parameters involved in the low-saturation extension of the primary drying curve.
Real dsaturationHys(Real pc, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
Derivative of Hysteretic saturation function with respect to pc.
const Real _s_l_min
Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity) ...

◆ initQpStatefulProperties()

void PorousFlow2PhaseHysPP::initQpStatefulProperties ( )
overrideprotectedvirtual

Reimplemented from PorousFlowHystereticCapillaryPressure.

Definition at line 60 of file PorousFlow2PhaseHysPP.C.

61 {
63  buildQpPPSS();
64 }
void buildQpPPSS()
Assemble std::vectors of porepressure and saturation at the nodes and quadpoints. ...

◆ landSat()

Real PorousFlowHystereticCapillaryPressure::landSat ( Real  slDel) const
protectedinherited
Returns
the value of gas saturation (called S_gr^Delta in the markdown documentation) using the Land expression. (This is a function of the liquid saturation at the turning point)
Parameters
slDelthe value of the liquid saturation at the turning point

Definition at line 233 of file PorousFlowHystereticCapillaryPressure.C.

Referenced by PorousFlowHystereticCapillaryPressure::computeTurningPointInfo().

234 {
235  const Real a = 1.0 / _s_gr_max - 1.0 / (1.0 - _s_lr);
236  return (1.0 - slDel) / (1.0 + a * (1.0 - slDel));
237 }
const Real _s_gr_max
Residual gas saturation: 1 - _s_gr_max is the maximum saturation for which the van Genuchten expressi...
const Real _s_lr
Liquid saturation below which the liquid relative permeability is zero.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ liquidSaturationQp()

Real PorousFlowHystereticCapillaryPressure::liquidSaturationQp ( Real  pc) const
protectedinherited
Returns
the value of liquid saturation, given the capillary pressure. This uses _hys_order[_qp]
Parameters
pccapillary pressure

Definition at line 630 of file PorousFlowHystereticCapillaryPressure.C.

Referenced by PorousFlow1PhaseHysP::buildQpPPSS(), buildQpPPSS(), and PorousFlowHystereticInfo::computeQpInfo().

631 {
632  Real sat = 0.0;
633  if (_hys_order[_qp] == 0) // on primary drying curve
635  else if (_hys_order[_qp] == 1) // first-order wetting
636  sat = firstOrderWettingSat(pc);
637  else if (_hys_order[_qp] == 2) // second-order drying
638  sat = secondOrderDryingSat(pc);
639  else // third order drying and wetting
640  {
641  // NOTE: this is not the exact inverse of the third-order capillary-pressure formula, but only
642  // the approximate inverse. In any simulation, only liquidSaturationQp or capillaryPressureQp
643  // are used (not both) so having a slightly different formulation for these two functions is OK
644  // Rationale: when pc is close to pc_tp1 then use the first-order wetting curve; when pc is
645  // close to pc_tp2 then use the second-order drying curve
646  const Real pc_tp1 = _pc_tps[_qp].at(1); // pc on first-order wetting at TP_1
647  const Real pc_tp2 = _pc_tps[_qp].at(2); // pc on second-order drying at TP_2
648  const Real sat1 = firstOrderWettingSat(pc);
649  const Real sat2 = secondOrderDryingSat(pc);
650  // handle cases that occur just at the transition from 3rd to 2nd order, or 3rd to 1st order
651  if (pc_tp1 <= 0.0 || pc <= 0.0 ||
652  pc_tp2 <= 0.0) // only the first condition is strictly necessary as cannot get pc==0 or
653  // pc_tp2=0 without pc_tp1=0 in reality. The other conditions are added in
654  // case of numerical strangenesses
655  sat = sat2;
656  else if (pc > pc_tp2)
657  sat = sat2;
658  else if (pc < pc_tp1)
659  sat = sat1;
660  else
661  sat = sat1 + (std::log(pc) - std::log(pc_tp1)) * (sat2 - sat1) / std::log(pc_tp2 / pc_tp1);
662  }
663  return sat;
664 }
const MaterialProperty< unsigned > & _hys_order
Hysteresis order, as computed by PorousFlowHysteresisOrder.
const Real _alpha_d
van Genuchten alpha parameter for the primary drying curve
const Real _n_d
van Genuchten n parameter for the primary drying curve
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _pc_tps
Nodal or quadpoint values of Pc at the turning points.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
Parameters involved in the low-saturation extension of the primary drying curve.
const Real _s_l_min
Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity) ...
Real saturationHys(Real pc, Real slmin, Real sgrdel, Real alpha, Real n, const LowCapillaryPressureExtension &low_ext=LowCapillaryPressureExtension(), const HighCapillaryPressureExtension &high_ext=HighCapillaryPressureExtension())
Hysteretic saturation function (Eqn(1) of Doughty2007) with extensions (page5 and Fig1 of Doughty2008...

◆ validParams()

InputParameters PorousFlow2PhaseHysPP::validParams ( )
static

Definition at line 15 of file PorousFlow2PhaseHysPP.C.

16 {
18  params.addRequiredCoupledVar("phase0_porepressure",
19  "Variable that is the porepressure of phase 0. This is assumed to "
20  "be the liquid phase. It will be <= phase1_porepressure.");
21  params.addRequiredCoupledVar("phase1_porepressure",
22  "Variable that is the porepressure of phase 1 (the gas phase)");
23  params.addParam<unsigned>(
24  "liquid_phase", 0, "Phase number of the liquid phase (more precisely, the phase of phase0)");
25  params.addClassDescription(
26  "This Material is used for 2-phase situations. It calculates the 2 saturations given the 2 "
27  "porepressures, assuming the capillary pressure is hysteretic. Derivatives of these "
28  "quantities are also computed");
29  return params;
30 }
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
void addClassDescription(const std::string &doc_string)

Member Data Documentation

◆ _alpha_d

const Real PorousFlowHystereticCapillaryPressure::_alpha_d
protectedinherited

◆ _alpha_w

const Real PorousFlowHystereticCapillaryPressure::_alpha_w
protectedinherited

◆ _dgradp_qp_dgradv

template<bool is_ad>
MaterialProperty<std::vector<std::vector<Real> > >* const PorousFlowVariableBaseTempl< is_ad >::_dgradp_qp_dgradv
protectedinherited

d(grad porepressure)/d(grad PorousFlow variable) at the quadpoints

Definition at line 51 of file PorousFlowVariableBase.h.

◆ _dgradp_qp_dv

template<bool is_ad>
MaterialProperty<std::vector<std::vector<RealGradient> > >* const PorousFlowVariableBaseTempl< is_ad >::_dgradp_qp_dv
protectedinherited

d(grad porepressure)/d(PorousFlow variable) at the quadpoints

Definition at line 54 of file PorousFlowVariableBase.h.

◆ _dgrads_qp_dgradv

template<bool is_ad>
MaterialProperty<std::vector<std::vector<Real> > >* const PorousFlowVariableBaseTempl< is_ad >::_dgrads_qp_dgradv
protectedinherited

d(grad saturation)/d(grad PorousFlow variable) at the quadpoints

Definition at line 66 of file PorousFlowVariableBase.h.

◆ _dgrads_qp_dv

template<bool is_ad>
MaterialProperty<std::vector<std::vector<RealGradient> > >* const PorousFlowVariableBaseTempl< is_ad >::_dgrads_qp_dv
protectedinherited

d(grad saturation)/d(PorousFlow variable) at the quadpoints

Definition at line 69 of file PorousFlowVariableBase.h.

◆ _dpc_high

const Real PorousFlowHystereticCapillaryPressure::_dpc_high
protectedinherited

d(Pc)/d(S) at the point of high-saturation extension

Definition at line 124 of file PorousFlowHystereticCapillaryPressure.h.

◆ _dpc_low_d

const Real PorousFlowHystereticCapillaryPressure::_dpc_low_d
protectedinherited

d(Pc)/dS on the primary drying curve at S = _s_low_d

Definition at line 108 of file PorousFlowHystereticCapillaryPressure.h.

◆ _dpc_low_w

const Real PorousFlowHystereticCapillaryPressure::_dpc_low_w
protectedinherited

d(Pc)/dS on the primary wetting curve at S = _s_low_w

Definition at line 114 of file PorousFlowHystereticCapillaryPressure.h.

◆ _dporepressure_dvar

template<bool is_ad>
MaterialProperty<std::vector<std::vector<Real> > >* const PorousFlowVariableBaseTempl< is_ad >::_dporepressure_dvar
protectedinherited

d(porepressure)/d(PorousFlow variable)

Definition at line 45 of file PorousFlowVariableBase.h.

◆ _dsaturation_dvar

template<bool is_ad>
MaterialProperty<std::vector<std::vector<Real> > >* const PorousFlowVariableBaseTempl< is_ad >::_dsaturation_dvar
protectedinherited

d(saturation)/d(PorousFlow variable)

Definition at line 60 of file PorousFlowVariableBase.h.

◆ _gradp_qp

template<bool is_ad>
GenericMaterialProperty<std::vector<RealGradient>, is_ad>* const PorousFlowVariableBaseTempl< is_ad >::_gradp_qp
protectedinherited

Grad(p) at the quadpoints.

Definition at line 48 of file PorousFlowVariableBase.h.

◆ _grads_qp

template<bool is_ad>
GenericMaterialProperty<std::vector<RealGradient>, is_ad>* const PorousFlowVariableBaseTempl< is_ad >::_grads_qp
protectedinherited

Grad(s) at the quadpoints.

Definition at line 63 of file PorousFlowVariableBase.h.

◆ _high_ext

const PorousFlowVanGenuchten::HighCapillaryPressureExtension PorousFlowHystereticCapillaryPressure::_high_ext
protectedinherited

Parameters involved in the high-saturation extension of the primary wetting curve.

Definition at line 126 of file PorousFlowHystereticCapillaryPressure.h.

◆ _high_ext_type

const PorousFlowVanGenuchten::HighCapillaryPressureExtension::ExtensionStrategy PorousFlowHystereticCapillaryPressure::_high_ext_type
protectedinherited

Type of high-saturation extension of the wetting curves.

Definition at line 118 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::computeTurningPointInfo().

◆ _high_ratio

const Real PorousFlowHystereticCapillaryPressure::_high_ratio
protectedinherited

The high-saturation extension to the wetting will commence at _high_ratio * (1 - _s_gr_del)

Definition at line 102 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::computeTurningPointInfo().

◆ _hys_order

const MaterialProperty<unsigned>& PorousFlowHystereticCapillaryPressure::_hys_order
protectedinherited

◆ _hys_order_old

const MaterialProperty<unsigned>& PorousFlowHystereticCapillaryPressure::_hys_order_old
protectedinherited

Old value of hysteresis order, as computed by PorousFlowHysteresisOrder.

Definition at line 132 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::computeQpProperties().

◆ _hys_sat_tps

const MaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER> >& PorousFlowHystereticCapillaryPressure::_hys_sat_tps
protectedinherited

◆ _low_ext_d

const PorousFlowVanGenuchten::LowCapillaryPressureExtension PorousFlowHystereticCapillaryPressure::_low_ext_d
protectedinherited

◆ _low_ext_type

const PorousFlowVanGenuchten::LowCapillaryPressureExtension::ExtensionStrategy PorousFlowHystereticCapillaryPressure::_low_ext_type
protectedinherited

Type of low-saturation extension.

Definition at line 104 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::computeTurningPointInfo().

◆ _low_ext_w

const PorousFlowVanGenuchten::LowCapillaryPressureExtension PorousFlowHystereticCapillaryPressure::_low_ext_w
protectedinherited

Parameters involved in the low-saturation extension of the primary wetting curve.

Definition at line 116 of file PorousFlowHystereticCapillaryPressure.h.

◆ _n_d

const Real PorousFlowHystereticCapillaryPressure::_n_d
protectedinherited

◆ _n_w

const Real PorousFlowHystereticCapillaryPressure::_n_w
protectedinherited

◆ _num_components

template<bool is_ad>
const unsigned int PorousFlowVariableBaseTempl< is_ad >::_num_components
protectedinherited

◆ _num_pf_vars

template<bool is_ad>
const unsigned int PorousFlowVariableBaseTempl< is_ad >::_num_pf_vars
protectedinherited

Number of PorousFlow variables.

Definition at line 39 of file PorousFlowVariableBase.h.

◆ _num_phases

template<bool is_ad>
const unsigned int PorousFlowVariableBaseTempl< is_ad >::_num_phases
protectedinherited

◆ _p0var

const unsigned int PorousFlow2PhaseHysPP::_p0var
protected

PorousFlow variable number of the phase0 porepressure.

Definition at line 46 of file PorousFlow2PhaseHysPP.h.

Referenced by computeQpProperties().

◆ _p1var

const unsigned int PorousFlow2PhaseHysPP::_p1var
protected

PorousFlow variable number of the phase1 porepressure.

Definition at line 54 of file PorousFlow2PhaseHysPP.h.

Referenced by computeQpProperties().

◆ _pc

MaterialProperty<Real>& PorousFlow2PhaseHysPP::_pc
protected

Computed nodal or quadpoint values of capillary pressure.

Definition at line 31 of file PorousFlow2PhaseHysPP.h.

Referenced by buildQpPPSS(), and computeQpProperties().

◆ _pc_high

const Real PorousFlowHystereticCapillaryPressure::_pc_high
protectedinherited

Pc at the point of high-saturation extension.

Definition at line 122 of file PorousFlowHystereticCapillaryPressure.h.

◆ _pc_max

const Real PorousFlowHystereticCapillaryPressure::_pc_max
protectedinherited

Maximum capillary pressure: for Pc above this value, a "lower" extension will be used.

Definition at line 100 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::computeTurningPointInfo().

◆ _pc_older

const MaterialProperty<Real>& PorousFlowHystereticCapillaryPressure::_pc_older
protectedinherited

Older value of capillary pressure.

Definition at line 139 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::computeQpProperties().

◆ _pc_tps

MaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER> >& PorousFlowHystereticCapillaryPressure::_pc_tps
protectedinherited

◆ _phase0_gradp_qp

const VariableGradient& PorousFlow2PhaseHysPP::_phase0_gradp_qp
protected

Gradient(phase0_porepressure) at the qps.

Definition at line 42 of file PorousFlow2PhaseHysPP.h.

Referenced by computeQpProperties().

◆ _phase0_porepressure

const VariableValue& PorousFlow2PhaseHysPP::_phase0_porepressure
protected

Nodal or quadpoint value of porepressure of the zero phase (eg, the water phase)

Definition at line 40 of file PorousFlow2PhaseHysPP.h.

Referenced by buildQpPPSS(), and computeQpProperties().

◆ _phase0_porepressure_varnum

const unsigned int PorousFlow2PhaseHysPP::_phase0_porepressure_varnum
protected

Moose variable number of the phase0 porepressure.

Definition at line 44 of file PorousFlow2PhaseHysPP.h.

Referenced by computeQpProperties().

◆ _phase1_gradp_qp

const VariableGradient& PorousFlow2PhaseHysPP::_phase1_gradp_qp
protected

Gradient(phase1_porepressure) at the qps.

Definition at line 50 of file PorousFlow2PhaseHysPP.h.

Referenced by computeQpProperties().

◆ _phase1_porepressure

const VariableValue& PorousFlow2PhaseHysPP::_phase1_porepressure
protected

Nodal or quadpoint value of porepressure of the one phase (eg, the gas phase)

Definition at line 48 of file PorousFlow2PhaseHysPP.h.

Referenced by buildQpPPSS(), and computeQpProperties().

◆ _phase1_porepressure_varnum

const unsigned int PorousFlow2PhaseHysPP::_phase1_porepressure_varnum
protected

Moose variable number of the phase1 porepressure.

Definition at line 52 of file PorousFlow2PhaseHysPP.h.

Referenced by computeQpProperties().

◆ _porepressure

template<bool is_ad>
GenericMaterialProperty<std::vector<Real>, is_ad>& PorousFlowVariableBaseTempl< is_ad >::_porepressure
protectedinherited

◆ _s_d_tps

MaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER> >& PorousFlowHystereticCapillaryPressure::_s_d_tps
protectedinherited

◆ _s_gr_max

const Real PorousFlowHystereticCapillaryPressure::_s_gr_max
protectedinherited

Residual gas saturation: 1 - _s_gr_max is the maximum saturation for which the van Genuchten expression is valid for the wetting curve.

Definition at line 98 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::landSat(), and PorousFlowHystereticCapillaryPressure::PorousFlowHystereticCapillaryPressure().

◆ _s_gr_tps

MaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER> >& PorousFlowHystereticCapillaryPressure::_s_gr_tps
protectedinherited

◆ _s_high

const Real PorousFlowHystereticCapillaryPressure::_s_high
protectedinherited

Saturation at the point of high-saturation extension.

Definition at line 120 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::PorousFlowHystereticCapillaryPressure().

◆ _s_l_min

const Real PorousFlowHystereticCapillaryPressure::_s_l_min
protectedinherited

Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity)

Definition at line 94 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::capillaryPressureQp(), PorousFlowHystereticCapillaryPressure::computeTurningPointInfo(), PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp(), PorousFlowHystereticCapillaryPressure::d2firstOrderWettingPc(), PorousFlowHystereticCapillaryPressure::d2firstOrderWettingSat(), PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp(), PorousFlowHystereticCapillaryPressure::d2secondOrderDryingPc(), PorousFlowHystereticCapillaryPressure::d2secondOrderDryingSat(), PorousFlowHystereticCapillaryPressure::dcapillaryPressureQp(), PorousFlowHystereticCapillaryPressure::dfirstOrderWettingPc(), PorousFlowHystereticCapillaryPressure::dfirstOrderWettingSat(), PorousFlowHystereticCapillaryPressure::dliquidSaturationQp(), PorousFlowHystereticCapillaryPressure::dsecondOrderDryingPc(), PorousFlowHystereticCapillaryPressure::dsecondOrderDryingSat(), PorousFlowHystereticCapillaryPressure::firstOrderWettingPc(), PorousFlowHystereticCapillaryPressure::firstOrderWettingSat(), PorousFlowHystereticCapillaryPressure::initQpStatefulProperties(), PorousFlowHystereticCapillaryPressure::liquidSaturationQp(), PorousFlowHystereticCapillaryPressure::PorousFlowHystereticCapillaryPressure(), PorousFlowHystereticCapillaryPressure::secondOrderDryingPc(), and PorousFlowHystereticCapillaryPressure::secondOrderDryingSat().

◆ _s_low_d

const Real PorousFlowHystereticCapillaryPressure::_s_low_d
protectedinherited

Saturation on the primary drying curve where low-saturation extension commences.

Definition at line 106 of file PorousFlowHystereticCapillaryPressure.h.

◆ _s_low_w

const Real PorousFlowHystereticCapillaryPressure::_s_low_w
protectedinherited

Saturation on the primary wetting curve where low-saturation extension commences.

Definition at line 112 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::PorousFlowHystereticCapillaryPressure().

◆ _s_lr

const Real PorousFlowHystereticCapillaryPressure::_s_lr
protectedinherited

Liquid saturation below which the liquid relative permeability is zero.

Definition at line 96 of file PorousFlowHystereticCapillaryPressure.h.

Referenced by PorousFlowHystereticCapillaryPressure::landSat(), and PorousFlowHystereticCapillaryPressure::PorousFlowHystereticCapillaryPressure().

◆ _s_w_tps

MaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER> >& PorousFlowHystereticCapillaryPressure::_s_w_tps
protectedinherited

◆ _saturation

template<bool is_ad>
GenericMaterialProperty<std::vector<Real>, is_ad>& PorousFlowVariableBaseTempl< is_ad >::_saturation
protectedinherited

◆ _w_high_ext_tps

MaterialProperty<std::array<PorousFlowVanGenuchten::HighCapillaryPressureExtension, PorousFlowConstants::MAX_HYSTERESIS_ORDER> >& PorousFlowHystereticCapillaryPressure::_w_high_ext_tps
protectedinherited

◆ _w_low_ext_tps

MaterialProperty<std::array<PorousFlowVanGenuchten::LowCapillaryPressureExtension, PorousFlowConstants::MAX_HYSTERESIS_ORDER> >& PorousFlowHystereticCapillaryPressure::_w_low_ext_tps
protectedinherited

The documentation for this class was generated from the following files: