37 Real n = 1.0 / (1.0 - m);
40 Real dseff_dp = -m *
std::pow(inner, -m - 1) * dinner_dp;
52 Real n = 1.0 / (1.0 - m);
57 m *
std::pow(inner, -m - 1.0) * d2inner_dp2;
79 if (seff <= 0.0 || seff >= 1.0)
95 if (seff <= 0.0 || seff >= 1.0)
107 d2pc *= (1.0 - m) / m /
alpha;
117 if (seff <= 0.0 || seff >= 1.0)
121 const Real da = -1.0 / m *
std::pow(seff, 1.0 / m - 1.0);
132 if (seff <= 0.0 || seff >= 1.0)
136 const Real da = -1.0 / m *
std::pow(seff, 1.0 / m - 1.0);
137 const Real d2a = -(1.0 / m) * (1.0 / m - 1.0) *
std::pow(seff, 1.0 / m - 2.0);
142 return -0.25 *
std::pow(seff, -1.5) * Utility::pow<2>(
b) + 2.0 *
std::pow(seff, -0.5) *
b *
db +
150 if (seff <= 0.0 || seff >= 1.0)
154 const Real da = -1.0 / m *
a / (1.0 - seff);
156 const Real db = -2.0 * m *
b / (1.0 -
a) * da;
165 if (seff <= 0.0 || seff >= 1.0)
169 const Real da = -1.0 / m *
a / (1.0 - seff);
170 const Real d2a = 1.0 / m * (1.0 / m - 1) *
std::pow(1.0 - seff, 1.0 / m - 2.0);
172 const Real db = -2.0 * m *
b / (1.0 -
a) * da;
174 -2.0 * m * (
db / (1.0 -
a) * da +
b * Utility::pow<2>(da / (1.0 -
a)) +
b / (1.0 -
a) * d2a);
194 pc = low_ext.
Pc + low_ext.
dPc * 0.5 * (sl * sl - low_ext.
S * low_ext.
S) / low_ext.
S;
197 pc = low_ext.
Pc * std::exp(low_ext.
dPc * (sl - low_ext.
S) / low_ext.
Pc);
214 const Real expon = -high_ext.
dPc / high_ext.
Pc * (1.0 - high_ext.
S);
215 pc = high_ext.
Pc *
std::pow((1.0 - sl) / (1.0 - high_ext.
S), expon);
224 const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
227 else if (seff <= 0.0)
252 dpc = low_ext.
dPc * sl / low_ext.
S;
255 dpc = low_ext.
dPc * std::exp(low_ext.
dPc * (sl - low_ext.
S) / low_ext.
Pc);
272 const Real expon = -high_ext.
dPc / high_ext.
Pc * (1.0 - high_ext.
S);
273 dpc = high_ext.
dPc *
std::pow((1.0 - sl) / (1.0 - high_ext.
S), expon - 1.0);
282 const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
285 else if (seff <= 0.0)
290 const Real dseff = 1.0 / (1.0 - sgrdel - slmin);
292 std::pow(seff, n / (1.0 - n) - 1.0);
293 dpc = dpc_dseff * dseff;
313 d2pc = low_ext.
dPc / low_ext.
S;
317 std::exp(low_ext.
dPc * (sl - low_ext.
S) / low_ext.
Pc);
334 const Real expon = -high_ext.
dPc / high_ext.
Pc * (1.0 - high_ext.
S);
335 d2pc = high_ext.
dPc * (1.0 - expon) / (1.0 - high_ext.
S) *
336 std::pow((1.0 - sl) / (1.0 - high_ext.
S), expon - 2.0);
345 const Real seff = (sl - slmin) / (1.0 - sgrdel - slmin);
348 else if (seff <= 0.0)
353 const Real dseff = 1.0 / (1.0 - sgrdel - slmin);
354 const Real d2pc_dseff =
355 (1.0 /
alpha / (1.0 - n)) *
357 (n / (1.0 - n) - 1.0) *
std::pow(
a, 1.0 / n - 1.0) *
std::pow(seff, n / (1.0 - n) - 2.0));
358 d2pc = d2pc_dseff * dseff * dseff;
381 const Real s2 = low_ext.
S * low_ext.
S + 2.0 * (pc - low_ext.
Pc) * low_ext.
S / low_ext.
dPc;
392 const Real ss = low_ext.
S + std::log(pc / low_ext.
Pc) * low_ext.
Pc / low_ext.
dPc;
406 if (pc < high_ext.
Pc)
412 const Real expon = -high_ext.
dPc / high_ext.
Pc * (1.0 - high_ext.
S);
413 s = 1.0 -
std::pow(pc / high_ext.
Pc, 1.0 / expon) * (1.0 - high_ext.
S);
421 if (pc == std::numeric_limits<Real>::max())
427 s = (1.0 - sgrdel - slmin) * seff + slmin;
450 const Real s2 = low_ext.
S * low_ext.
S + 2.0 * (pc - low_ext.
Pc) * low_ext.
S / low_ext.
dPc;
457 const Real ds2 = 2.0 * low_ext.
S / low_ext.
dPc;
464 const Real s = low_ext.
S + std::log(pc / low_ext.
Pc) * low_ext.
Pc / low_ext.
dPc;
470 ds = low_ext.
Pc / pc / low_ext.
dPc;
478 if (pc < high_ext.
Pc)
484 const Real expon = -high_ext.
dPc / high_ext.
Pc * (1.0 - high_ext.
S);
485 ds = -(1.0 - high_ext.
S) / pc / expon *
std::pow(pc / high_ext.
Pc, 1.0 / expon);
493 if (pc == std::numeric_limits<Real>::max())
498 const Real dseffpow = n * (seffpow - 1.0) / pc;
500 const Real dseff = (1.0 - n) / n * seff / seffpow * dseffpow;
501 ds = (1.0 - sgrdel - slmin) * dseff;
524 const Real s2 = low_ext.
S * low_ext.
S + 2.0 * (pc - low_ext.
Pc) * low_ext.
S / low_ext.
dPc;
531 const Real ds2 = 2.0 * low_ext.
S / low_ext.
dPc;
532 d2s = -0.25 * ds2 * ds2 /
std::pow(s2, 1.5);
538 const Real s = low_ext.
S + std::log(pc / low_ext.
Pc) * low_ext.
Pc / low_ext.
dPc;
552 if (pc < high_ext.
Pc)
558 const Real expon = -high_ext.
dPc / high_ext.
Pc * (1.0 - high_ext.
S);
559 d2s = -(1.0 - high_ext.
S) * (1.0 / expon) * (1.0 / expon - 1.0) /
568 if (pc == std::numeric_limits<Real>::max())
573 const Real dseffpow = n * (seffpow - 1.0) / pc;
574 const Real d2seffpow = (n - 1.0) * dseffpow / pc;
576 const Real dseff = (1.0 - n) / n * seff / seffpow * dseffpow;
579 (dseff * dseffpow - seff * dseffpow * dseffpow / seffpow + seff * d2seffpow) / seffpow;
580 d2s = (1.0 - sgrdel - slmin) * d2seff;
592 Real upper_liquid_param,
600 const Real sl_bar = (sl - slr) / (1.0 - slr);
605 if (sgrdel == 0.0 || sl <= sldel)
612 const Real my_sldel = (sldel < slr) ? slr : sldel;
613 const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
614 if (sl >= 1.0 - 0.5 * my_sgrdel)
620 else if (sl > upper_liquid_param * (1.0 - my_sgrdel))
625 sl, upper_liquid_param * (1.0 - my_sgrdel), y0, y0p, 1.0 - 0.5 * my_sgrdel, y1, y1p);
630 const Real sl_bar_del = (my_sldel - slr) / (1.0 - slr);
631 const Real s_gt_bar =
632 my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
633 a = (1 - s_gt_bar / (1.0 - sl_bar_del)) *
635 b = s_gt_bar / (1.0 - sl_bar_del) *
std::pow(1.0 -
std::pow(sl_bar_del, 1.0 / m), m);
638 return std::sqrt(sl_bar) * Utility::pow<2>(1.0 -
a -
b);
648 Real upper_liquid_param,
657 return std::numeric_limits<Real>::max();
658 const Real sl_bar = (sl - slr) / (1.0 - slr);
659 const Real sl_bar_prime = 1.0 / (1.0 - slr);
666 if (sgrdel == 0.0 || sl <= sldel)
669 const Real dc_dsbar =
c / m / sl_bar;
671 const Real da_dsbar = -m *
a / (1.0 -
c) * dc_dsbar;
672 a_prime = da_dsbar * sl_bar_prime;
679 const Real my_sldel = (sldel < slr) ? slr : sldel;
680 const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
681 if (sl >= 1.0 - 0.5 * my_sgrdel)
686 const Real dc_dsbar =
c / m / sl_bar;
688 const Real da_dsbar = -m *
a / (1.0 -
c) * dc_dsbar;
689 a_prime = da_dsbar * sl_bar_prime;
691 else if (sl > upper_liquid_param * (1.0 - my_sgrdel))
696 sl, upper_liquid_param * (1.0 - my_sgrdel), y0, y0p, 1.0 - 0.5 * my_sgrdel, y1, y1p);
701 const Real sl_bar_del = (my_sldel - slr) / (1.0 - slr);
702 const Real s_gt_bar =
703 my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
704 const Real s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
706 const Real c_prime =
c / m / (sl_bar + s_gt_bar) * (sl_bar_prime + s_gt_bar_prime);
707 a = (1 - s_gt_bar / (1.0 - sl_bar_del)) *
std::pow(1.0 -
c, m);
709 -s_gt_bar_prime / (1.0 - sl_bar_del) *
std::pow(1.0 -
c, m) - m *
a / (1.0 -
c) * c_prime;
710 b = s_gt_bar / (1.0 - sl_bar_del) *
std::pow(1.0 -
std::pow(sl_bar_del, 1.0 / m), m);
711 b_prime = s_gt_bar_prime *
b / s_gt_bar;
715 return 0.5 * kr / sl_bar * sl_bar_prime -
716 std::sqrt(sl_bar) * 2.0 * (1.0 -
a -
b) * (a_prime + b_prime);
737 if (sl > 1.0 - sgrdel)
739 const Real sl_bar = (sl - slr) / (1.0 - slr);
741 if (sgrdel != 0.0 && sl > sldel)
747 const Real my_sldel = (sldel < slr) ? slr : sldel;
748 const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
749 s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
752 if (sl_bar + s_gt_bar < 1.0)
758 kr = k_rg_max *
a *
b;
781 if (sl > 1.0 - sgrdel)
783 const Real sl_bar = (sl - slr) / (1.0 - slr);
784 const Real sl_bar_prime = 1.0 / (1.0 - slr);
786 Real s_gt_bar_prime = 0.0;
787 if (sgrdel != 0.0 && sl > sldel)
793 const Real my_sldel = (sldel < slr) ? slr : sldel;
794 const Real my_sgrdel = (sldel < slr) ? sgrmax : sgrdel;
795 s_gt_bar = my_sgrdel * (sl - my_sldel) / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
796 s_gt_bar_prime = my_sgrdel / (1.0 - slr) / (1.0 - my_sldel - my_sgrdel);
799 if (sl_bar + s_gt_bar < 1.0)
803 const Real a_prime = -gamma *
a / (1.0 - (sl_bar + s_gt_bar)) * (sl_bar_prime + s_gt_bar_prime);
806 (
c == 0 ? 0.0 :
c / m / (sl_bar + s_gt_bar) * (sl_bar_prime + s_gt_bar_prime));
808 const Real b_prime = -2.0 * m *
b / (1.0 -
c) * c_prime;
809 kr_prime = k_rg_max * (
a * b_prime + a_prime *
b);
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.
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.
Real drelativePermeabilityHys(Real sl, Real slr, Real sgrdel, Real sgrmax, Real sldel, Real m, Real upper_liquid_param, Real y0, Real y0p, Real y1, Real y1p)
Derivative of Hysteretic relative permeability for liquid, with respect to liquid saturation...
Real relativePermeabilityNWHys(Real sl, Real slr, Real sgrdel, Real sgrmax, Real sldel, Real m, Real gamma, Real k_rg_max, Real y0p)
Hysteretic relative permeability for gas.
Real cubic(Real x, Real x0, Real y0, Real y0p, Real x1, Real y1, Real y1p)
Cubic function f(x) that satisfies f(x0) = y0 f'(x0) = y0p f(x1) = y1 f'(x1) = y1p.
Parameters associated with the extension of the hysteretic wetting capillary pressure function to hig...
ExtensionStrategy strategy
ExtensionStrategy strategy
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
van Genuchten effective saturation, capillary pressure and relative permeability functions.
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...
Real drelativePermeabilityNWHys(Real sl, Real slr, Real sgrdel, Real sgrmax, Real sldel, Real m, Real gamma, Real k_rg_max, Real y0p)
Derivative of hysteretic relative permeability for gas with respect to the liquid saturation...
Parameters associated with the extension of the hysteretic capillary pressure function to low saturat...
const GeochemicalDatabaseReader db("database/moose_testdb.json", true, true, false)
Real d2RelativePermeabilityNW(Real seff, Real m)
Second derivative of relative permeability for a non-wetting phase with respect to effective saturati...
Real relativePermeabilityHys(Real sl, Real slr, Real sgrdel, Real sgrmax, Real sldel, Real m, Real upper_liquid_param, Real y0, Real y0p, Real y1, Real y1p)
Hysteretic relative permeability for liquid.
Real d2CapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
Second derivative of capillary pressure wrt effective saturation.
Real dRelativePermeability(Real seff, Real m)
Derivative of relative permeability with respect to effective saturation.
Real dcubic(Real x, Real x0, Real y0, Real y0p, Real x1, Real y1, Real y1p)
Derivative of cubic function, f(x), with respect to x.
Real effectiveSaturation(Real p, Real alpha, Real m)
Effective saturation as a function of porepressure.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string alpha
Real dCapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
Derivative of capillary pressure wrt effective saturation.
Real d2EffectiveSaturation(Real p, Real alpha, Real m)
Second derivative of effective saturation wrt porepressure.
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.
Real d2RelativePermeability(Real seff, Real m)
Second derivative of relative permeability with respect to effective saturation.
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.
Real dEffectiveSaturation(Real p, Real alpha, Real m)
Derivative of effective saturation wrt porepressure.
Real dRelativePermeabilityNW(Real seff, Real m)
Derivative of relative permeability for a non-wetting phase with respect to effective saturation...
MooseUnits pow(const MooseUnits &, int)
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 capillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
Capillary pressure as a function of effective saturation.