https://mooseframework.inl.gov
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
PorousFlow1PhaseHysP Class Reference

Material designed to calculate fluid phase porepressure and saturation for the single-phase partially-saturation situation with hysteretic capillary pressure and assuming porepressure is a nonlinear variable. More...

#include <PorousFlow1PhaseHysP.h>

Inheritance diagram for PorousFlow1PhaseHysP:
[legend]

Public Types

typedef DerivativeMaterialPropertyNameInterface::SymbolName SymbolName
 

Public Member Functions

 PorousFlow1PhaseHysP (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.
 
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.
 
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.
 
const VariableValue_porepressure_var
 Nodal or quadpoint value of porepressure of the fluid phase.
 
const VariableGradient_gradp_qp_var
 Gradient(_porepressure at quadpoints)
 
const unsigned int _porepressure_varnum
 Moose variable number of the porepressure.
 
const unsigned int _p_var_num
 The PorousFlow variable number of the porepressure.
 
const Real _alpha_d
 van Genuchten alpha parameter for the primary drying curve
 
const Real _alpha_w
 van Genuchten alpha parameter for the primary wetting curve
 
const Real _n_d
 van Genuchten n parameter for the primary drying curve
 
const Real _n_w
 van Genuchten n parameter for the primary wetting curve
 
const Real _s_l_min
 Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity)
 
const Real _s_lr
 Liquid saturation below which the liquid relative permeability is zero.
 
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.
 
const Real _pc_max
 Maximum capillary pressure: for Pc above this value, a "lower" extension will be used.
 
const Real _high_ratio
 The high-saturation extension to the wetting will commence at _high_ratio * (1 - _s_gr_del)
 
const PorousFlowVanGenuchten::LowCapillaryPressureExtension::ExtensionStrategy _low_ext_type
 Type of low-saturation extension.
 
const Real _s_low_d
 Saturation on the primary drying curve where low-saturation extension commences.
 
const Real _dpc_low_d
 d(Pc)/dS on the primary drying curve at S = _s_low_d
 
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
 Parameters involved in the low-saturation extension of the primary drying curve.
 
const Real _s_low_w
 Saturation on the primary wetting curve where low-saturation extension commences.
 
const Real _dpc_low_w
 d(Pc)/dS on the primary wetting curve at S = _s_low_w
 
const PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_w
 Parameters involved in the low-saturation extension of the primary wetting curve.
 
const PorousFlowVanGenuchten::HighCapillaryPressureExtension::ExtensionStrategy _high_ext_type
 Type of high-saturation extension of the wetting curves.
 
const Real _s_high
 Saturation at the point of high-saturation extension.
 
const Real _pc_high
 Pc at the point of high-saturation extension.
 
const Real _dpc_high
 d(Pc)/d(S) at the point of high-saturation extension
 
const PorousFlowVanGenuchten::HighCapillaryPressureExtension _high_ext
 Parameters involved in the high-saturation extension of the primary wetting curve.
 
const MaterialProperty< unsigned > & _hys_order
 Hysteresis order, as computed by PorousFlowHysteresisOrder.
 
const MaterialProperty< unsigned > & _hys_order_old
 Old value of hysteresis order, as computed by PorousFlowHysteresisOrder.
 
const MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _hys_sat_tps
 Saturation values at the turning points, as computed by PorousFlowHysteresisOrder.
 
const MaterialProperty< Real > & _pc_older
 Older value of capillary pressure.
 
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _pc_tps
 Nodal or quadpoint values of Pc at the turning points.
 
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.
 
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< 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.
 
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.
 
const unsigned int _num_phases
 Number of phases.
 
const unsigned int _num_components
 Number of components.
 
const unsigned int _num_pf_vars
 Number of PorousFlow variables.
 
GenericMaterialProperty< std::vector< Real >, is_ad > & _porepressure
 Computed nodal or quadpoint values of porepressure of the phases.
 
MaterialProperty< std::vector< std::vector< Real > > > *const _dporepressure_dvar
 d(porepressure)/d(PorousFlow variable)
 
GenericMaterialProperty< std::vector< RealGradient >, is_ad > *const _gradp_qp
 Grad(p) at the quadpoints.
 
MaterialProperty< std::vector< std::vector< Real > > > *const _dgradp_qp_dgradv
 d(grad porepressure)/d(grad PorousFlow variable) at the quadpoints
 
MaterialProperty< std::vector< std::vector< RealGradient > > > *const _dgradp_qp_dv
 d(grad porepressure)/d(PorousFlow variable) at the quadpoints
 
GenericMaterialProperty< std::vector< Real >, is_ad > & _saturation
 Computed nodal or qp saturation of the phases.
 
MaterialProperty< std::vector< std::vector< Real > > > *const _dsaturation_dvar
 d(saturation)/d(PorousFlow variable)
 
GenericMaterialProperty< std::vector< RealGradient >, is_ad > *const _grads_qp
 Grad(s) at the quadpoints.
 
MaterialProperty< std::vector< std::vector< Real > > > *const _dgrads_qp_dgradv
 d(grad saturation)/d(grad PorousFlow variable) at the quadpoints
 
MaterialProperty< std::vector< std::vector< RealGradient > > > *const _dgrads_qp_dv
 d(grad saturation)/d(PorousFlow variable) at the quadpoints
 

Private Member Functions

Real firstOrderWettingPc (Real sat) const
 
Real dfirstOrderWettingPc (Real sat) const
 
Real d2firstOrderWettingPc (Real sat) const
 
Real secondOrderDryingPc (Real sat) const
 
Real dsecondOrderDryingPc (Real sat) const
 
Real d2secondOrderDryingPc (Real sat) const
 
Real firstOrderWettingSat (Real pc) const
 
Real dfirstOrderWettingSat (Real pc) const
 
Real d2firstOrderWettingSat (Real pc) const
 
Real secondOrderDryingSat (Real pc) const
 
Real dsecondOrderDryingSat (Real pc) const
 
Real d2secondOrderDryingSat (Real pc) const
 
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

Material designed to calculate fluid phase porepressure and saturation for the single-phase partially-saturation situation with hysteretic capillary pressure and assuming porepressure is a nonlinear variable.

Definition at line 19 of file PorousFlow1PhaseHysP.h.

Constructor & Destructor Documentation

◆ PorousFlow1PhaseHysP()

PorousFlow1PhaseHysP::PorousFlow1PhaseHysP ( const InputParameters parameters)

Definition at line 27 of file PorousFlow1PhaseHysP.C.

29 _pc(_nodal_material ? declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_nodal")
30 : declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_qp")),
31 _porepressure_var(_nodal_material ? coupledDofValues("porepressure")
32 : coupledValue("porepressure")),
33 _gradp_qp_var(coupledGradient("porepressure")),
34 _porepressure_varnum(coupled("porepressure")),
35 _p_var_num(_dictator.isPorousFlowVariable(_porepressure_varnum)
36 ? _dictator.porousFlowVariableNum(_porepressure_varnum)
37 : 0)
38{
39 if (_num_phases != 1)
40 mooseError("The Dictator proclaims that the number of phases is ",
41 _dictator.numPhases(),
42 " whereas PorousFlow1PhaseHysP can only be used for 1-phase simulations. Be aware "
43 "that the Dictator has noted your mistake.");
44}
void mooseError(Args &&... args)
const VariableGradient & _gradp_qp_var
Gradient(_porepressure at quadpoints)
const VariableValue & _porepressure_var
Nodal or quadpoint value of porepressure of the fluid phase.
const unsigned int _p_var_num
The PorousFlow variable number of the porepressure.
const unsigned int _porepressure_varnum
Moose variable number of the porepressure.
MaterialProperty< Real > & _pc
Computed nodal or quadpoint values of capillary pressure.
Base material designed to calculate and store quantities relevant for hysteretic capillary pressure c...
const unsigned int _num_phases
Number of phases.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

Member Function Documentation

◆ buildQpPPSS()

void PorousFlow1PhaseHysP::buildQpPPSS ( )
protected

Assemble std::vectors of porepressure and saturation.

Definition at line 85 of file PorousFlow1PhaseHysP.C.

86{
87 _porepressure[_qp][0] = _porepressure_var[_qp];
88 _pc[_qp] = -_porepressure_var[_qp];
89 _saturation[_qp][0] = liquidSaturationQp(_pc[_qp]);
90}
GenericMaterialProperty< std::vector< Real >, is_ad > & _porepressure
Computed nodal or quadpoint values of porepressure of the phases.
GenericMaterialProperty< std::vector< Real >, is_ad > & _saturation
Computed nodal or qp saturation of the phases.

Referenced by computeQpProperties(), and initQpStatefulProperties().

◆ 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.

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 PorousFlowVanGenuchten::LowCapillaryPressureExtension _low_ext_d
Parameters involved in the low-saturation extension of the primary drying curve.
const MaterialProperty< unsigned > & _hys_order
Hysteresis order, as computed by PorousFlowHysteresisOrder.
const Real _s_l_min
Minimum liquid saturation for which the van Genuchten expression is valid (Pc(_s_l_min) = infinity)
const MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _hys_sat_tps
Saturation values at the turning points, as computed by PorousFlowHysteresisOrder.
const Real _n_d
van Genuchten n parameter for the primary drying curve
const Real _alpha_d
van Genuchten alpha 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...

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

◆ computeQpProperties()

void PorousFlow1PhaseHysP::computeQpProperties ( )
overrideprotectedvirtual

Reimplemented from PorousFlowHystereticCapillaryPressure.

Definition at line 54 of file PorousFlow1PhaseHysP.C.

55{
56 // size stuff correctly and prepare the derivative matrices with zeroes
58
60 const Real pc = -_porepressure_var[_qp];
61 const Real ds = -dliquidSaturationQp(pc);
62
63 if (!_nodal_material)
64 {
65 (*_gradp_qp)[_qp][0] = _gradp_qp_var[_qp];
66 (*_grads_qp)[_qp][0] = ds * _gradp_qp_var[_qp];
67 }
68
69 // _porepressure is only dependent on _porepressure, and its derivative is 1
70 if (_dictator.isPorousFlowVariable(_porepressure_varnum))
71 {
72 // _porepressure is a PorousFlow variable
73 (*_dporepressure_dvar)[_qp][0][_p_var_num] = 1.0;
74 (*_dsaturation_dvar)[_qp][0][_p_var_num] = ds;
75 if (!_nodal_material)
76 {
77 (*_dgradp_qp_dgradv)[_qp][0][_p_var_num] = 1.0;
78 (*_dgrads_qp_dgradv)[_qp][0][_p_var_num] = ds;
79 (*_dgrads_qp_dv)[_qp][0][_p_var_num] = d2liquidSaturationQp(pc) * _gradp_qp_var[_qp];
80 }
81 }
82}
void buildQpPPSS()
Assemble std::vectors of porepressure and saturation.

◆ 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.

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
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}
const Real _high_ratio
The high-saturation extension to the wetting will commence at _high_ratio * (1 - _s_gr_del)
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.
const PorousFlowVanGenuchten::HighCapillaryPressureExtension::ExtensionStrategy _high_ext_type
Type of high-saturation extension of the wetting curves.
MaterialProperty< std::array< Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER > > & _pc_tps
Nodal or quadpoint values of Pc at the turning points.
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 Real _alpha_w
van Genuchten alpha parameter for the primary wetting 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.
const Real _pc_max
Maximum capillary pressure: for Pc above this value, a "lower" extension will be used.
const Real _n_w
van Genuchten n parameter for the primary wetting curve
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 PorousFlowVanGenuchten::LowCapillaryPressureExtension::ExtensionStrategy _low_ext_type
Type of low-saturation extension.
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.
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...
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.
Parameters associated with the extension of the hysteretic wetting capillary pressure function to hig...
Parameters associated with the extension of the hysteretic capillary pressure function to low saturat...

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

◆ 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.

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.

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

◆ d2firstOrderWettingPc()

Real PorousFlowHystereticCapillaryPressure::d2firstOrderWettingPc ( Real  sat) const
privateinherited
Returns
d^2(capillary pressure on the first-order wetting curve)/d(sat)^2
Parameters
satliquid saturation

Definition at line 467 of file PorousFlowHystereticCapillaryPressure.C.

468{
469 const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
471 ? 1.0
472 : 1.0 - _s_gr_tps[_qp].at(0);
473 const Real sat_to_use = _s_w_tps[_qp].at(0) + (max_s - _s_w_tps[_qp].at(0)) *
474 (sat - _hys_sat_tps[_qp].at(0)) /
475 (max_s - _hys_sat_tps[_qp].at(0));
476 const Real dsat_to_use = (max_s - _s_w_tps[_qp].at(0)) / (max_s - _hys_sat_tps[_qp].at(0));
478 _s_l_min,
479 _s_gr_tps[_qp].at(0),
480 _alpha_w,
481 _n_w,
482 _w_low_ext_tps[_qp].at(0),
483 _w_high_ext_tps[_qp].at(0)) *
484 dsat_to_use * dsat_to_use;
485}

Referenced by PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp().

◆ d2firstOrderWettingSat()

Real PorousFlowHystereticCapillaryPressure::d2firstOrderWettingSat ( Real  pc) const
privateinherited
Returns
d^2(saturation on the first-order wetting curve)/d(pc)^2
Parameters
pccapillary pressure

Definition at line 572 of file PorousFlowHystereticCapillaryPressure.C.

573{
574 const Real d2sat_to_use = PorousFlowVanGenuchten::d2saturationHys(pc,
575 _s_l_min,
576 _s_gr_tps[_qp].at(0),
577 _alpha_w,
578 _n_w,
579 _w_low_ext_tps[_qp].at(0),
580 _w_high_ext_tps[_qp].at(0));
581 const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
583 ? 1.0
584 : 1.0 - _s_gr_tps[_qp].at(0);
585 if (pc <= 0.0 && _w_high_ext_tps[_qp].at(0).strategy ==
587 return 0.0;
588 return d2sat_to_use * (max_s - _hys_sat_tps[_qp].at(0)) / (max_s - _s_w_tps[_qp].at(0));
589}
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.

Referenced by PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp().

◆ 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.

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}

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

◆ d2secondOrderDryingPc()

Real PorousFlowHystereticCapillaryPressure::d2secondOrderDryingPc ( Real  sat) const
privateinherited
Returns
d^2(capillary pressure on the second-order drying curve)/d(sat)^2
Parameters
satliquid saturation

Definition at line 517 of file PorousFlowHystereticCapillaryPressure.C.

518{
519 const Real tp0 = _hys_sat_tps[_qp].at(0);
520 const Real tp1 = _hys_sat_tps[_qp].at(1);
521 const Real s1 = _s_d_tps[_qp].at(1);
522 const Real sat_to_use = (sat >= tp0) ? tp0 + (sat - tp0) * (s1 - tp0) / (tp1 - tp0) : sat;
523 const Real dsat_to_use = (sat >= tp0) ? (s1 - tp0) / (tp1 - tp0) : 1.0;
525 sat_to_use, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d) *
526 dsat_to_use * dsat_to_use;
527}

Referenced by PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp().

◆ d2secondOrderDryingSat()

Real PorousFlowHystereticCapillaryPressure::d2secondOrderDryingSat ( Real  pc) const
privateinherited
Returns
d^2(saturation on the second-order drying curve)/d(pc)^2
Parameters
pccapillary pressure

Definition at line 617 of file PorousFlowHystereticCapillaryPressure.C.

618{
619 const Real sat_to_use =
621 const Real d2sat_to_use =
623 const Real tp0 = _hys_sat_tps[_qp].at(0);
624 const Real tp1 = _hys_sat_tps[_qp].at(1);
625 const Real s1 = _s_d_tps[_qp].at(1);
626 return (sat_to_use >= tp0) ? d2sat_to_use * (tp1 - tp0) / (s1 - tp0) : d2sat_to_use;
627}

Referenced by PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp().

◆ 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.

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}

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

◆ dfirstOrderWettingPc()

Real PorousFlowHystereticCapillaryPressure::dfirstOrderWettingPc ( Real  sat) const
privateinherited
Returns
d(capillary pressure on the first-order wetting curve)/d(sat)
Parameters
satliquid saturation

Definition at line 446 of file PorousFlowHystereticCapillaryPressure.C.

447{
448 const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
450 ? 1.0
451 : 1.0 - _s_gr_tps[_qp].at(0);
452 const Real sat_to_use = _s_w_tps[_qp].at(0) + (max_s - _s_w_tps[_qp].at(0)) *
453 (sat - _hys_sat_tps[_qp].at(0)) /
454 (max_s - _hys_sat_tps[_qp].at(0));
455 const Real dsat_to_use = (max_s - _s_w_tps[_qp].at(0)) / (max_s - _hys_sat_tps[_qp].at(0));
457 _s_l_min,
458 _s_gr_tps[_qp].at(0),
459 _alpha_w,
460 _n_w,
461 _w_low_ext_tps[_qp].at(0),
462 _w_high_ext_tps[_qp].at(0)) *
463 dsat_to_use;
464}

Referenced by PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp(), and PorousFlowHystereticCapillaryPressure::dcapillaryPressureQp().

◆ dfirstOrderWettingSat()

Real PorousFlowHystereticCapillaryPressure::dfirstOrderWettingSat ( Real  pc) const
privateinherited
Returns
d(saturation on the first-order wetting curve)/d(pc)
Parameters
pccapillary pressure

Definition at line 552 of file PorousFlowHystereticCapillaryPressure.C.

553{
554 const Real dsat_to_use = PorousFlowVanGenuchten::dsaturationHys(pc,
555 _s_l_min,
556 _s_gr_tps[_qp].at(0),
557 _alpha_w,
558 _n_w,
559 _w_low_ext_tps[_qp].at(0),
560 _w_high_ext_tps[_qp].at(0));
561 const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
563 ? 1.0
564 : 1.0 - _s_gr_tps[_qp].at(0);
565 if (pc <= 0.0 && _w_high_ext_tps[_qp].at(0).strategy ==
567 return 0.0;
568 return dsat_to_use * (max_s - _hys_sat_tps[_qp].at(0)) / (max_s - _s_w_tps[_qp].at(0));
569}
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.

Referenced by PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp(), and PorousFlowHystereticCapillaryPressure::dliquidSaturationQp().

◆ 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.

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}

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

◆ dsecondOrderDryingPc()

Real PorousFlowHystereticCapillaryPressure::dsecondOrderDryingPc ( Real  sat) const
privateinherited
Returns
d(capillary pressure on the second-order drying curve)/d(sat)
Parameters
satliquid saturation

Definition at line 504 of file PorousFlowHystereticCapillaryPressure.C.

505{
506 const Real tp0 = _hys_sat_tps[_qp].at(0);
507 const Real tp1 = _hys_sat_tps[_qp].at(1);
508 const Real s1 = _s_d_tps[_qp].at(1);
509 const Real sat_to_use = (sat >= tp0) ? tp0 + (sat - tp0) * (s1 - tp0) / (tp1 - tp0) : sat;
510 const Real dsat_to_use = (sat >= tp0) ? (s1 - tp0) / (tp1 - tp0) : 1.0;
512 sat_to_use, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d) *
513 dsat_to_use;
514}

Referenced by PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp(), and PorousFlowHystereticCapillaryPressure::dcapillaryPressureQp().

◆ dsecondOrderDryingSat()

Real PorousFlowHystereticCapillaryPressure::dsecondOrderDryingSat ( Real  pc) const
privateinherited
Returns
d(saturation on the second-order drying curve)/d(pc)
Parameters
pccapillary pressure

Definition at line 604 of file PorousFlowHystereticCapillaryPressure.C.

605{
606 const Real sat_to_use =
608 const Real dsat_to_use =
610 const Real tp0 = _hys_sat_tps[_qp].at(0);
611 const Real tp1 = _hys_sat_tps[_qp].at(1);
612 const Real s1 = _s_d_tps[_qp].at(1);
613 return (sat_to_use >= tp0) ? dsat_to_use * (tp1 - tp0) / (s1 - tp0) : dsat_to_use;
614}

Referenced by PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp(), and PorousFlowHystereticCapillaryPressure::dliquidSaturationQp().

◆ firstOrderWettingPc()

Real PorousFlowHystereticCapillaryPressure::firstOrderWettingPc ( Real  sat) const
privateinherited
Returns
capillary pressure on the first-order wetting curve
Parameters
satliquid saturation

Definition at line 421 of file PorousFlowHystereticCapillaryPressure.C.

422{
423 // Simplest version is to just use the wetting curve defined by _s_gr_tps[0], but want
424 // continuity at sat = _hys_sat_tps[0] (the turning point), so use the following process. The
425 // wetting curve is defined for S <= max_s, where
426 const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
428 ? 1.0
429 : 1.0 - _s_gr_tps[_qp].at(0);
430 // define an interpolation: s_to_use smoothly transitions from _s_w_tps[0] (the value of liquid
431 // saturation on the wetting curve defined by _s_gr_tps[0]) when sat = _hys_sat_tps[0], to max_s
432 // when sat = max_s
433 const Real sat_to_use = _s_w_tps[_qp].at(0) + (max_s - _s_w_tps[_qp].at(0)) *
434 (sat - _hys_sat_tps[_qp].at(0)) /
435 (max_s - _hys_sat_tps[_qp].at(0));
437 _s_l_min,
438 _s_gr_tps[_qp].at(0),
439 _alpha_w,
440 _n_w,
441 _w_low_ext_tps[_qp].at(0),
442 _w_high_ext_tps[_qp].at(0));
443}

Referenced by PorousFlowHystereticCapillaryPressure::capillaryPressureQp(), PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp(), PorousFlowHystereticCapillaryPressure::dcapillaryPressureQp(), and PorousFlowHystereticCapillaryPressure::initQpStatefulProperties().

◆ firstOrderWettingSat()

Real PorousFlowHystereticCapillaryPressure::firstOrderWettingSat ( Real  pc) const
privateinherited
Returns
saturation on the first-order wetting curve
Parameters
pccapillary pressure

Definition at line 530 of file PorousFlowHystereticCapillaryPressure.C.

531{
532 // this is inverse of firstOrderWettingPc: see that method for comments
533 const Real sat_to_use = PorousFlowVanGenuchten::saturationHys(pc,
534 _s_l_min,
535 _s_gr_tps[_qp].at(0),
536 _alpha_w,
537 _n_w,
538 _w_low_ext_tps[_qp].at(0),
539 _w_high_ext_tps[_qp].at(0));
540 const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
542 ? 1.0
543 : 1.0 - _s_gr_tps[_qp].at(0);
544 if (sat_to_use > max_s) // this occurs when using no high extension and pc = 0
545 return max_s;
546 return (sat_to_use - _s_w_tps[_qp].at(0)) * (max_s - _hys_sat_tps[_qp].at(0)) /
547 (max_s - _s_w_tps[_qp].at(0)) +
548 _hys_sat_tps[_qp].at(0);
549}

Referenced by PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp(), PorousFlowHystereticCapillaryPressure::dliquidSaturationQp(), and PorousFlowHystereticCapillaryPressure::liquidSaturationQp().

◆ initQpStatefulProperties()

void PorousFlow1PhaseHysP::initQpStatefulProperties ( )
overrideprotectedvirtual

◆ 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.

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.

Referenced by PorousFlowHystereticCapillaryPressure::computeTurningPointInfo().

◆ 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.

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}

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

◆ secondOrderDryingPc()

Real PorousFlowHystereticCapillaryPressure::secondOrderDryingPc ( Real  sat) const
privateinherited
Returns
capillary pressure on the second-order drying curve
Parameters
satliquid saturation

Definition at line 488 of file PorousFlowHystereticCapillaryPressure.C.

489{
490 // Simplest version is to just use the primary drying curve, but want
491 // continuity at sat = _hys_sat_tps[0] (the dry-to-wet turning point) and sat = _hys_sat_tps[1]
492 // (the wet-to-dry turning point), so use the following process.
493 const Real tp0 = _hys_sat_tps[_qp].at(0);
494 const Real tp1 = _hys_sat_tps[_qp].at(1);
495 const Real s1 = _s_d_tps[_qp].at(1);
496 const Real sat_to_use =
497 (sat >= tp0) ? tp0 + (sat - tp0) * (s1 - tp0) / (tp1 - tp0)
498 : sat; // final case can occur just at transition from 2nd to 0th order
500 sat_to_use, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
501}

Referenced by PorousFlowHystereticCapillaryPressure::capillaryPressureQp(), PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp(), PorousFlowHystereticCapillaryPressure::dcapillaryPressureQp(), and PorousFlowHystereticCapillaryPressure::initQpStatefulProperties().

◆ secondOrderDryingSat()

Real PorousFlowHystereticCapillaryPressure::secondOrderDryingSat ( Real  pc) const
privateinherited
Returns
saturation on the second-order drying curve
Parameters
pccapillary pressure

Definition at line 592 of file PorousFlowHystereticCapillaryPressure.C.

593{
594 // This is the inverse of secondOrderDryingPc: see that method for comments
595 const Real sat_to_use =
597 const Real tp0 = _hys_sat_tps[_qp].at(0);
598 const Real tp1 = _hys_sat_tps[_qp].at(1);
599 const Real s1 = _s_d_tps[_qp].at(1);
600 return (sat_to_use >= tp0) ? (sat_to_use - tp0) * (tp1 - tp0) / (s1 - tp0) + tp0 : sat_to_use;
601}

Referenced by PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp(), PorousFlowHystereticCapillaryPressure::dliquidSaturationQp(), and PorousFlowHystereticCapillaryPressure::liquidSaturationQp().

◆ validParams()

InputParameters PorousFlow1PhaseHysP::validParams ( )
static

Definition at line 15 of file PorousFlow1PhaseHysP.C.

16{
19 "porepressure", "Variable that represents the porepressure of the single liquid phase");
21 "This Material is used for unsaturated single-phase situations "
22 "where porepressure is the primary variable and the capillary pressure is hysteretic. The "
23 "hysteretic formulation assumes that the single phase is a liquid");
24 return params;
25}
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.

◆ _gradp_qp_var

const VariableGradient& PorousFlow1PhaseHysP::_gradp_qp_var
protected

Gradient(_porepressure at quadpoints)

Definition at line 41 of file PorousFlow1PhaseHysP.h.

Referenced by computeQpProperties().

◆ _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

◆ _p_var_num

const unsigned int PorousFlow1PhaseHysP::_p_var_num
protected

The PorousFlow variable number of the porepressure.

Definition at line 45 of file PorousFlow1PhaseHysP.h.

Referenced by computeQpProperties().

◆ _pc

MaterialProperty<Real>& PorousFlow1PhaseHysP::_pc
protected

Computed nodal or quadpoint values of capillary pressure.

Definition at line 31 of file PorousFlow1PhaseHysP.h.

Referenced by buildQpPPSS().

◆ _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

◆ _porepressure

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

◆ _porepressure_var

const VariableValue& PorousFlow1PhaseHysP::_porepressure_var
protected

Nodal or quadpoint value of porepressure of the fluid phase.

Definition at line 39 of file PorousFlow1PhaseHysP.h.

Referenced by buildQpPPSS(), and computeQpProperties().

◆ _porepressure_varnum

const unsigned int PorousFlow1PhaseHysP::_porepressure_varnum
protected

Moose variable number of the porepressure.

Definition at line 43 of file PorousFlow1PhaseHysP.h.

Referenced by computeQpProperties().

◆ _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: