3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
parkervangenuchten.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef PARKER_VANGENUCHTEN_3P_HH
25#define PARKER_VANGENUCHTEN_3P_HH
26
27#include <algorithm>
33
34namespace Dumux::FluidMatrix {
35
37{
44 template<class Scalar>
45 struct Params
46 {
47 Params(const Scalar swr = 0.0, const Scalar snr = 0.0, const Scalar sgr = 0.0)
48 : swr_(swr), snr_(snr), sgr_(sgr)
49 {}
50
54 Scalar swr() const
55 { return swr_; }
56
60 void setSwr(Scalar v)
61 { swr_ = v; }
62
66 Scalar snr() const
67 { return snr_; }
68
72 void setSnr(Scalar v)
73 { snr_ = v; }
78 Scalar sgr() const
79 { return sgr_; }
80
84 void setSgr(Scalar v)
85 { sgr_ = v; }
86
87 bool operator== (const Params& p) const
88 {
89 return Dune::FloatCmp::eq(swr(), p.swr(), 1e-6)
90 && Dune::FloatCmp::eq(snr(), p.snr(), 1e-6)
91 && Dune::FloatCmp::eq(sgr(), p.sgr(), 1e-6);
92 }
93 private:
94 Scalar swr_;
95 Scalar snr_;
96 Scalar sgr_;
97 };
98
103 template<class Scalar>
104 static Params<Scalar> makeParams(const std::string& paramGroup)
105 {
106 Params<Scalar> params;
107 params.setSwr(getParamFromGroup<Scalar>(paramGroup, "Swr", 0.0));
108 params.setSnr(getParamFromGroup<Scalar>(paramGroup, "Snr", 0.0));
109 params.setSgr(getParamFromGroup<Scalar>(paramGroup, "Sgr", 0.0));
110 return params;
111 }
112
122 template<class Scalar>
123 static Scalar swToSwe(const Scalar sw, const Params<Scalar>& params)
124 {
125 return (sw - params.swr())/(1.0 - params.swr()); // TODO other residual saturations?
126 }
127
137 template<class Scalar>
138 static Scalar sweToSw(const Scalar swe, const Params<Scalar>& params)
139 {
140 return swe*(1.0 - params.swr()) + params.swr(); // TODO other residual saturations?
141 }
142
151 template<class Scalar>
152 static Scalar dswe_dsw(const Params<Scalar>& params)
153 {
154 return 1.0/(1.0 - params.swr()); // TODO other residual saturations?
155 }
156
165 template<class Scalar>
166 static Scalar dsw_dswe(const Params<Scalar>& params)
167 {
168 return 1.0 - params.swr(); // TODO other residual saturations?
169 }
170
180 template<class Scalar>
181 static Scalar snToSne(const Scalar sn, const Params<Scalar>& params)
182 {
183 return sn; // sne equals sn // TODO other residual saturations?
184 }
185
195 template<class Scalar>
196 static Scalar stToSte(const Scalar st, const Params<Scalar>& params)
197 {
198 return (st-params.swr()) / (1-params.swr()); // TODO other residual saturations?
199 }
200
210 template<class Scalar>
211 static Scalar steToSt(const Scalar ste, const Params<Scalar>& params)
212 {
213 return ste*(1.0 - params.swr()) + params.swr(); // TODO other residual saturations?
214 }
215
224 template<class Scalar>
225 static Scalar dste_dst(const Params<Scalar>& params)
226 {
227 return 1.0/(1.0 - params.swr() /*- params.snr() - params.sgr()*/); // TODO other residual saturations?
228 }
229
238 template<class Scalar>
239 static Scalar dst_dste(const Params<Scalar>& params)
240 {
241 return 1.0 - params.swr(); // TODO other residual saturations?
242 }
243};
244
253{
254
255public:
270 template<class Scalar>
271 struct Params
272 {
273 Params(Scalar alpha, Scalar n, Scalar swr = 0.0, Scalar snr = 0.0,
274 Scalar betaNw = 1.0, Scalar betaGn = 1.0, Scalar betaGw = 1.0, bool regardSnr = false)
275 : alpha_(alpha), n_(n), m_(1.0 - 1.0/n), swr_(swr), snr_(snr)
276 , betaNw_(betaNw), betaGn_(betaGn), betaGw_(betaGw), regardSnr_(regardSnr)
277 {}
278
279 Scalar alpha() const { return alpha_; }
280 void setAlpha(Scalar alpha) { alpha_ = alpha; }
281
282 Scalar m() const { return m_; }
283 void setM(Scalar m) { m_ = m; n_ = 1.0/(1.0 - m); }
284
285 Scalar n() const{ return n_; }
286 void setN(Scalar n){ n_ = n; m_ = 1.0 - 1.0/n; }
287
288 Scalar swr() const { return swr_; }
289 void setSwr(Scalar swr) { swr_ = swr; }
290
291 Scalar snr() const { return snr_; }
292 void setSnr(Scalar swr) { snr_ = swr; }
293
294 Scalar betaNw() const { return betaNw_; }
295 void setBetaNw(Scalar betaNw) { betaNw_ = betaNw; }
296
297 Scalar betaGn() const { return betaGn_; }
298 void setBetaGn(Scalar betaGn) { betaGn_ = betaGn; }
299
300 Scalar betaGw() const { return betaGw_; }
301 void setBetaGw(Scalar betaGw) { betaGw_ = betaGw; }
302
303 bool regardSnrForKrn() const { return regardSnr_; }
304 void setRegardSnrForKrn(bool v) {regardSnr_ = v; }
305
306 bool operator== (const Params& p) const
307 {
308 return Dune::FloatCmp::eq(alpha_, p.alpha_, 1e-6)
309 && Dune::FloatCmp::eq(n_, p.n_, 1e-6)
310 && Dune::FloatCmp::eq(m_, p.m_, 1e-6)
311 && Dune::FloatCmp::eq(swr_, p.swr_, 1e-6)
312 && Dune::FloatCmp::eq(snr_, p.snr_, 1e-6)
313 && Dune::FloatCmp::eq(betaNw_, p.betaNw_, 1e-6)
314 && Dune::FloatCmp::eq(betaGn_, p.betaGn_, 1e-6)
315 && Dune::FloatCmp::eq(betaGw_, p.betaGw_, 1e-6)
316 && regardSnr_ == p.regardSnr_;
317 }
318
319 private:
320 Scalar alpha_, n_, m_, swr_, snr_, betaNw_, betaGn_, betaGw_;
321 bool regardSnr_;
322 };
323
328 template<class Scalar = double>
329 static Params<Scalar> makeParams(const std::string& paramGroup)
330 {
331 const auto n = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenN");
332 const auto alpha = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenAlpha");
333 const auto swr = getParamFromGroup<Scalar>(paramGroup, "Swr", 0.0);
334 const auto snr = getParamFromGroup<Scalar>(paramGroup, "Snr", 0.0);
335 const auto betaNw = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenBetaNw", 1.0);
336 const auto betaGn = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenBetaGn", 1.0);
337 const auto betaGw = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenBetaGw", 1.0);
338 const auto regardSnr = getParamFromGroup<bool>(paramGroup, "ParkerVanGenuchtenRegardSnrForKrn", false);
339 return Params<Scalar>(alpha, n, swr, snr,
340 betaNw, betaGn, betaGw, regardSnr );
341 }
342
348 template<class Scalar>
349 static Scalar pcgw(Scalar swe, const Params<Scalar>& params)
350 {
351 assert(0 <= swe && swe <= 1);
352 return pc_(swe, params);
353 }
354
360 template<class Scalar>
361 static Scalar pcnw(Scalar swe, const Params<Scalar>& params)
362 {
363 assert(0 <= swe && swe <= 1);
364 return pc_(swe, params)/params.betaNw();
365 }
366
372 template<class Scalar>
373 static Scalar pcgn(const Scalar ste, const Params<Scalar>& params)
374 {
375 assert(0 <= ste && ste <= 1);
376 return pc_(ste, params)/params.betaGn();
377 }
378
384 template<class Scalar>
385 static Scalar pcAlpha(Scalar sne, const Params<Scalar>& params)
386 {
387 /* regularization */
388 if (sne <= 0.001)
389 sne = 0.0;
390 if (sne >= 1.0)
391 sne = 1.0;
392
393 if (sne > params.snr())
394 return 1.0;
395 else
396 {
397 if (params.snr() >= 0.001)
398 return sne/params.snr();
399 else
400 return 0.0;
401 };
402 }
403
410 template<class Scalar>
411 static Scalar dpcgw_dswe(const Scalar swe, const Params<Scalar>& params)
412 {
413 using std::pow;
414 const Scalar powSeRegu = pow(swe, -1/params.m());
415 return - 1.0/params.alpha() * pow(powSeRegu - 1, 1.0/params.n() - 1)/params.n()
416 * powSeRegu/swe/params.m()/params.betaGw();
417 }
418
425 template<class Scalar>
426 static Scalar dpcnw_dswe(const Scalar swe, const Params<Scalar>& params)
427 {
428 using std::pow;
429 const Scalar powSeRegu = pow(swe, -1/params.m());
430 return - 1.0/params.alpha() * pow(powSeRegu - 1, 1.0/params.n() - 1)/params.n()
431 * powSeRegu/swe/params.m()/params.betaNw();
432 }
433
440 template<class Scalar>
441 static Scalar dpcgn_dste(const Scalar ste, const Params<Scalar>& params)
442 {
443 using std::pow;
444 const Scalar powSeRegu = pow(ste, -1/params.m());
445 return - 1.0/params.alpha() * pow(powSeRegu - 1, 1.0/params.n() - 1)/params.n()
446 * powSeRegu/ste/params.m()/params.betaGn();
447 }
448
461 template<class Scalar>
462 static Scalar krw(const Scalar swe, const Params<Scalar>& params)
463 {
464 using std::pow;
465 using std::sqrt;
466 const Scalar r = 1.0 - pow(1 - pow(swe, 1/params.m()), params.m());
467 return sqrt(swe)*r*r;
468 }
469
487 template<class Scalar>
488 static Scalar krn(const Scalar swe, const Scalar sn, const Scalar ste, const Params<Scalar>& params)
489 {
490 Scalar krn;
491 using std::pow;
492 krn = pow(1 - pow(swe, 1/params.m()), params.m());
493 krn -= pow(1 - pow(ste, 1/params.m()), params.m());
494 krn *= krn;
495
496 using std::clamp;
497 using std::sqrt;
498 if (params.regardSnrForKrn())
499 {
500 // regard Snr in the permeability of the n-phase, see Helmig1997
501 const Scalar resIncluded = clamp(sn - params.snr()/ (1-params.swr()), 0.0, 1.0);
502 krn *= sqrt(resIncluded);
503 }
504 else
505 krn *= sqrt(sn / (1 - params.swr())); // Hint: (ste - swe) = sn / (1-Swr)
506
507 return krn;
508 }
509
522 template<class Scalar>
523 static Scalar krg(const Scalar ste, const Params<Scalar>& params)
524 {
525 assert(0 <= ste && ste <= 1);
526 using std::cbrt;
527 using std::pow;
528 return cbrt(1 - ste) * pow(1 - pow(ste, 1/params.m()), 2*params.m());
529 }
530
540 template<class Scalar>
541 static Scalar dkrg_dste(const Scalar ste, const Params<Scalar>& params)
542 {
543 assert(0 < ste && ste <= 1);
544
545 using std::pow;
546 const Scalar x = pow(ste, 1.0/params.m());
547 return
548 -pow(1.0 - x, 2*params.m())
549 *pow(1.0 - ste, -2.0/3)
550 *(1.0/3 + 2*x/ste);
551 }
552
560 template<class Scalar>
561 static Scalar kr(const int phaseIdx, const Scalar swe, const Scalar sne, const Params<Scalar>& params)
562 {
563 switch (phaseIdx)
564 {
565 case 0:
566 return krw(params, swe, sne);
567 case 1:
568 return krn(params, swe, sne);
569 case 2:
570 return krg(params, swe, sne);
571 }
572 DUNE_THROW(Dune::InvalidStateException,
573 "Invalid phase index ");
574 }
575
576private:
577
584 template<class Scalar>
585 const static Scalar pc_(const Scalar se, const Params<Scalar>& params)
586 {
587 using std::pow;
588 return pow(pow(se, -1/params.m()) - 1, 1/params.n())/params.alpha();
589 }
590
591};
592
600template <class Scalar>
602{
603 using BaseLawParams = typename ParkerVanGenuchten3P::Params<Scalar>;
604
605public:
607 template<class S>
608 struct Params
609 {
617 { pcLowSwe_ = pcLowSwe; }
618
622 Scalar pcLowSwe() const
623 { return pcLowSwe_; }
624
629 { pcHighSwe_ = pcHighSwe; }
630
637 Scalar pcHighSwe() const
638 { return pcHighSwe_; }
639
645 { krnLowSwe_ = krnLowSwe; }
646
651 Scalar krnLowSwe() const
652 { return krnLowSwe_; }
653
659 { krgLowSte_ = krgLowSte; }
660
665 Scalar krgLowSte() const
666 { return krgLowSte_; }
667
673 { krwHighSwe_ = krwHighSwe; }
674
679 Scalar krwHighSwe() const
680 { return krwHighSwe_; }
681
682
688 void setConstRegularization(const bool input)
689 { constRegularization_ = input; }
690
696 { return constRegularization_; }
697
698private:
699 S pcLowSwe_ = 0.01;
700 S pcHighSwe_ = 0.99;
701 S krnLowSwe_ = 0.1;
702 S krwHighSwe_ = 0.9;
703 S krgLowSte_ = 1e-3;
704 bool constRegularization_ = false;
705 };
706
708 template<class MaterialLaw>
709 void init(const MaterialLaw* m, const std::string& paramGroup)
710 {
711 pcLowSwe_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenPcLowSweThreshold", 0.01);
712 pcHighSwe_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenPcHighSweThreshold", 0.99);
713 krwHighSwe_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenKrwHighSweThreshold", 0.9);
714 krnLowSwe_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenKrnLowSweThreshold", 0.1);
715 krgLowSte_ = getParamFromGroup<Scalar>(paramGroup, "ParkerVanGenuchtenKrgLowSteThreshold", 1e-3);
716 constRegularization_ = getParamFromGroup<bool>(paramGroup, "VanGenuchtenConstantRegularization", false);
717
718 initPcParameters_(m, pcLowSwe_, pcHighSwe_);
719 initKrParameters_(m, krnLowSwe_, krwHighSwe_);
720 }
721
722 template<class MaterialLaw, class BaseParams, class EffToAbsParams>
723 void init(const MaterialLaw* m, const BaseParams& bp, const EffToAbsParams& etap, const Params<Scalar>& p)
724 {
725 pcLowSwe_ = p.pcLowSwe();
726 pcHighSwe_ = p.pcHighSwe();
727 krwHighSwe_ = p.krwHighSwe();
728 krnLowSwe_ = p.krnLowSwe();
729 krgLowSte_ = p.krgLowSte();
730 constRegularization_ = p.constRegularization();
731
732 initPcParameters_(m, pcLowSwe_, pcHighSwe_);
733 initKrParameters_(m, krnLowSwe_, krwHighSwe_);
734 }
735
740 {
741 return Dune::FloatCmp::eq(pcLowSwe_, o.pcLowSwe_, 1e-6)
742 && Dune::FloatCmp::eq(pcHighSwe_, o.pcHighSwe_, 1e-6)
743 && Dune::FloatCmp::eq(krwHighSwe_, o.krwHighSwe_, 1e-6)
744 && Dune::FloatCmp::eq(krnLowSwe_, o.krnLowSwe_, 1e-6)
745 && constRegularization_ == o.constRegularization_;
746 }
747
756 OptionalScalar<Scalar> pcgw(Scalar swe) const
757 {
758 // if specified, a constant value is used for regularization
759 using std::clamp;
760 if (constRegularization_)
761 swe = clamp(swe, 0.0, 1.0);
762
763 // make sure that the capillary pressure observes a derivative
764 // != 0 for 'illegal' saturations. This is favourable for the
765 // newton solver (if the derivative is calculated numerically)
766 // in order to get the saturation moving to the right
767 // direction if it temporarily is in an 'illegal' range.
768 if (swe < pcLowSwe_)
769 return pcgwLowSwePcgwValue_ + pcgwDerivativeLowSw_*(swe - pcLowSwe_);
770
771 else if (swe > 1.0)
772 return pcgwDerivativeHighSweEnd_*(swe - 1.0);
773
774 else if (swe > pcHighSwe_)
775 return pcgwSpline_.eval(swe);
776
777 else
778 return {}; // no regularization
779 }
780
789 OptionalScalar<Scalar> pcnw(Scalar swe) const
790 {
791 // if specified, a constant value is used for regularization
792 using std::clamp;
793 if (constRegularization_)
794 swe = clamp(swe, 0.0, 1.0);
795
796 // make sure that the capillary pressure observes a derivative
797 // != 0 for 'illegal' saturations. This is favourable for the
798 // newton solver (if the derivative is calculated numerically)
799 // in order to get the saturation moving to the right
800 // direction if it temporarily is in an 'illegal' range.
801 if (swe < pcLowSwe_)
802 return pcnwLowSwePcnwValue_ + pcnwDerivativeLowSw_*(swe - pcLowSwe_);
803
804 else if (swe > 1.0)
805 return pcnwDerivativeHighSweEnd_*(swe - 1.0);
806
807 else if (swe > pcHighSwe_)
808 return pcnwSpline_.eval(swe);
809
810 else
811 return {}; // no regularization
812 }
813
822 OptionalScalar<Scalar> pcgn(Scalar ste) const
823 {
824 // if specified, a constant value is used for regularization
825 using std::clamp;
826 if (constRegularization_)
827 ste = clamp(ste, 0.0, 1.0);
828
829
830 // make sure that the capillary pressure observes a derivative
831 // != 0 for 'illegal' saturations. This is favourable for the
832 // newton solver (if the derivative is calculated numerically)
833 // in order to get the saturation moving to the right
834 // direction if it temporarily is in an 'illegal' range.
835 const Scalar pcLowSte = pcLowSwe_;
836 const Scalar pcHighSte = pcHighSwe_;
837 if (ste < pcLowSte)
838 return pcgnLowStePcgnValue_ + pcgnDerivativeLowSt_*(ste - pcLowSte);
839
840 else if (ste > 1.0)
841 return pcgnDerivativeHighSteEnd_*(ste - 1.0);
842
843 else if (ste > pcHighSte)
844 return pcgnSpline_.eval(ste);
845
846 else
847 return {}; // no regularization
848 }
849
855 {
856 // no regularization
857 return {};
858 }
859
864 OptionalScalar<Scalar> krw(const Scalar swe) const
865 {
866 if (swe < 0.0)
867 return 0.0;
868 else if (swe > 1.0 - std::numeric_limits<Scalar>::epsilon())
869 return 1.0;
870 else
871 return {}; // no regularization
872 }
873
874
881 OptionalScalar<Scalar> krn(Scalar swe, const Scalar sn, Scalar ste) const
882 {
883 using std::clamp;
884 swe = clamp(swe, 0.0, 1.0);
885 ste = clamp(ste, 0.0, 1.0);
886
887 if (ste - swe <= 0.0)
888 return 0.0;
889 else
890 return ParkerVanGenuchten3P::krn(swe, sn, ste, *baseLawParamsPtr_);
891 }
892
897 OptionalScalar<Scalar> krg(const Scalar ste) const
898 {
899 //return 0 if there is no gas
900 if (ste > 1.0 - std::numeric_limits<Scalar>::epsilon())
901 return 0.0;
902
903 // use linear regularization for very high gas saturations
904 // to avoid a kink in the curve and to maintain a slope for
905 // the Newton solver
906 if (ste <= krgLowSte_)
907 return krgLowStkrgValue_ + krgDerivativeLowSt_*(ste - krgLowSte_);
908 else
909 {
910 // For very low gas saturations:
911 // We use a scaling factor that decreases the gas phase permeability quite fast a very low gas phase
912 // saturations, thus making that phase virtually immobile.
913 // This prevents numerical issues related to the degeneration of the gas phase mass balance for the 3p3c model
914 // at very low gas phase saturations.
915
916 // get the absolute gas phase saturation
917 const Scalar st = ste*(1 - swr_) + swr_;
918 const Scalar sg = 1.0 - st;
919
920 // do not regularize
921 if (sg > 0.1)
922 return {};
923
924 // return original curve scaled by factor
925 using std::max;
926 const Scalar scalFact = max(0.0, (sg - sgr_)/(0.1 - sgr_));
927
928 return ParkerVanGenuchten3P::krg(ste, *baseLawParamsPtr_) * scalFact;
929 }
930 }
931
938 OptionalScalar<Scalar> kr(const int phaseIdx, const Scalar swe, const Scalar sne) const
939 {
940 switch (phaseIdx)
941 {
942 case 0:
943 return krw(swe, sne);
944 case 1:
945 return krn(swe, sne);
946 case 2:
947 return krg(swe, sne);
948 }
949 DUNE_THROW(Dune::InvalidStateException,
950 "Invalid phase index ");
951 }
952
953private:
954 template<class MaterialLaw>
955 void initPcParameters_(const MaterialLaw* m, const Scalar lowSwe, const Scalar highSwe)
956 {
957 const auto lowSw = MaterialLaw::EffToAbs::sweToSw(lowSwe, m->effToAbsParams());
958 const auto highSw = MaterialLaw::EffToAbs::sweToSw(highSwe, m->effToAbsParams());
959 const auto dsw_dswe = MaterialLaw::EffToAbs::dsw_dswe(m->effToAbsParams());
960 const auto dst_dste = MaterialLaw::EffToAbs::dst_dste(m->effToAbsParams());
961
962 baseLawParamsPtr_ = &m->basicParams();
963
964 // pcgw
965 pcgwLowSwePcgwValue_ = m->template pcgw<false>(lowSw, 0.0);
966 pcgwDerivativeLowSw_ = m->template dpcgw_dsw<false>(lowSw, 0.0)*dsw_dswe;
967 pcgwHighSwePcgwValue_ = m->template pcgw<false>(highSw, 0.0);
968 pcgwDerivativeHighSweThreshold_ = m->template dpcgw_dsw<false>(highSw, 0.0)*dsw_dswe;
969 pcgwDerivativeHighSweEnd_ = 2.0*(0.0 - m->template pcgw<false>(highSw, 0.0))/(1.0 - highSwe);
970 pcgwSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1
971 pcgwHighSwePcgwValue_, 0, // y0, y1
972 pcgwDerivativeHighSweThreshold_, pcgwDerivativeHighSweEnd_); // m0, m1
973
974 // pcnw
975 pcnwLowSwePcnwValue_ = m->template pcnw<false>(lowSw, 0.0);
976 pcnwDerivativeLowSw_ = m->template dpcnw_dsw<false>(lowSw, 0.0)*dsw_dswe;
977 pcnwHighSwePcnwValue_ = m->template pcnw<false>(highSw, 0.0);
978 pcnwDerivativeHighSweThreshold_ = m->template dpcnw_dsw<false>(highSw, 0.0);
979 pcnwDerivativeHighSweEnd_ = 2.0*(0.0 - m->template pcnw<false>(highSw, 0.0))/(1.0 - highSwe);
980 pcnwSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1
981 pcnwHighSwePcnwValue_, 0, // y0, y1
982 pcnwDerivativeHighSweThreshold_, pcnwDerivativeHighSweEnd_); // m0, m1
983
984 // pcgn
985 pcgnLowStePcgnValue_ = m->template pcgn<false>(lowSw, 0.0);
986 pcgnDerivativeLowSt_ = m->template dpcgn_dst<false>(lowSw, 0.0)*dst_dste;
987 pcgnHighSwePcgnValue_ = m->template pcgn<false>(highSw, 0.0);
988 pcgnDerivativeHighSteThreshold_ = m->template dpcgn_dst<false>(highSw, 0.0);
989 pcgnDerivativeHighSteEnd_ = 2.0*(0.0 - m->template pcgn<false>(highSw, 0.0))/(1.0 - highSwe);
990 pcgnSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1
991 pcgnHighSwePcgnValue_, 0, // y0, y1
992 pcgnDerivativeHighSteThreshold_, pcgnDerivativeHighSteEnd_); // m0, m1
993
994 }
995
996 template<class MaterialLaw>
997 void initKrParameters_(const MaterialLaw* m, const Scalar lowSwe, const Scalar highSwe)
998 {
999 krgLowStkrgValue_ = ParkerVanGenuchten3P::krg(krgLowSte_, *baseLawParamsPtr_);
1000 krgDerivativeLowSt_ = ParkerVanGenuchten3P::dkrg_dste(krgLowSte_, *baseLawParamsPtr_);
1001
1002 swr_ = m->effToAbsParams().swr();
1003 sgr_ = m->effToAbsParams().sgr();
1004 }
1005
1006 Scalar krgLowStkrgValue_;
1007 Scalar krgDerivativeLowSt_;
1008
1009 Scalar pcLowSwe_, pcHighSwe_;
1010 Scalar krwHighSwe_, krnLowSwe_, krgLowSte_;
1011
1012 // pcgw
1013 Scalar pcgwLowSwePcgwValue_;
1014 Scalar pcgwHighSwePcgwValue_;
1015 Scalar pcgwDerivativeLowSw_;
1016 Scalar pcgwDerivativeHighSweThreshold_;
1017 Scalar pcgwDerivativeHighSweEnd_;
1018
1019 // pcgn
1020 Scalar pcgnLowStePcgnValue_;
1021 Scalar pcgnHighSwePcgnValue_;
1022 Scalar pcgnDerivativeLowSt_;
1023 Scalar pcgnDerivativeHighSteThreshold_;
1024 Scalar pcgnDerivativeHighSteEnd_;
1025
1026 // pcnw
1027 Scalar pcnwLowSwePcnwValue_;
1028 Scalar pcnwHighSwePcnwValue_;
1029 Scalar pcnwDerivativeLowSw_;
1030 Scalar pcnwDerivativeHighSweThreshold_;
1031 Scalar pcnwDerivativeHighSweEnd_;
1032
1033 Spline<Scalar> pcgwSpline_;
1034 Spline<Scalar> pcnwSpline_;
1035 Spline<Scalar> pcgnSpline_;
1036 Spline<Scalar> krwSpline_;
1037 Spline<Scalar> krnSpline_;
1038
1039 Scalar swr_, sgr_;
1040
1041 bool constRegularization_;
1042
1043 const BaseLawParams* baseLawParamsPtr_;
1044};
1045
1046template<class ScalarType,
1047 class BaseLaw,
1048 class Regularization = NoRegularization,
1049 class EffToAbsPolicy = ParkerVanGenuchten3PEffToAbsPolicy>
1050class ParkerVanGenuchtenMaterialLaw : public Adapter<ParkerVanGenuchtenMaterialLaw<ScalarType, BaseLaw, Regularization, EffToAbsPolicy>, ThreePhasePcKrSw>
1051{
1052public:
1053
1054 using Scalar = ScalarType;
1055
1056 using BasicParams = typename BaseLaw::template Params<Scalar>;
1057 using EffToAbsParams = typename EffToAbsPolicy::template Params<Scalar>;
1058 using RegularizationParams = typename Regularization::template Params<Scalar>;
1059
1060 using EffToAbs = EffToAbsPolicy;
1061
1065 static constexpr bool isRegularized()
1066 { return !std::is_same<Regularization, NoRegularization>::value; }
1067
1073
1078 explicit ParkerVanGenuchtenMaterialLaw(const std::string& paramGroup)
1079 : basicParams_(makeBasicParams(paramGroup))
1080 , effToAbsParams_(makeEffToAbsParams(paramGroup))
1081 {
1082 if constexpr (isRegularized())
1083 regularization_.init(this, paramGroup);
1084 }
1085
1091 const EffToAbsParams& effToAbsParams = {},
1092 const RegularizationParams& regParams = {})
1093 : basicParams_(baseParams)
1094 , effToAbsParams_(effToAbsParams)
1095 {
1096 if constexpr (isRegularized())
1097 regularization_.init(this, baseParams, effToAbsParams, regParams);
1098 }
1099
1103 template<bool enableRegularization = isRegularized()>
1104 Scalar pcgw(const Scalar sw, const Scalar /*dummySn*/) const
1105 {
1106 const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_);
1107 if constexpr (enableRegularization)
1108 {
1109 const auto regularized = regularization_.pcgw(swe);
1110 if (regularized)
1111 return regularized.value();
1112 }
1113
1114 return BaseLaw::pcgw(swe, basicParams_);
1115 }
1116
1120 template<bool enableRegularization = isRegularized()>
1121 Scalar pcnw(const Scalar sw, const Scalar /*dummySn*/) const
1122 {
1123 const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_);
1124 if constexpr (enableRegularization)
1125 {
1126 const auto regularized = regularization_.pcnw(swe);
1127 if (regularized)
1128 return regularized.value();
1129 }
1130
1131 return BaseLaw::pcnw(swe, basicParams_);
1132 }
1133
1139 template<bool enableRegularization = isRegularized()>
1140 Scalar pcgn(const Scalar sw, const Scalar sn) const
1141 {
1142 const auto swe = EffToAbs::swToSwe(sw + sn, effToAbsParams_);
1143 if constexpr (enableRegularization)
1144 {
1145 const auto regularized = regularization_.pcgn(swe);
1146 if (regularized)
1147 return regularized.value();
1148 }
1149
1150 return BaseLaw::pcgn(swe, basicParams_);
1151 }
1152
1156 template<bool enableRegularization = isRegularized()>
1157 Scalar pcAlpha(const Scalar /*dummySw*/, const Scalar sn) const
1158 {
1159 const auto sne = EffToAbs::snToSne(sn, effToAbsParams_);
1160 if constexpr (enableRegularization)
1161 {
1162 const auto regularized = regularization_.pcAlpha(sne);
1163 if (regularized)
1164 return regularized.value();
1165 }
1166 return BaseLaw::pcAlpha(sne, basicParams_);
1167 }
1168
1172 template<bool enableRegularization = isRegularized()>
1173 Scalar dpcgw_dsw(const Scalar sw, const Scalar /*dummyS*/) const
1174 {
1175 const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_);
1176 if constexpr (enableRegularization)
1177 {
1178 const auto regularized = regularization_.dpcgw_dswe(swe);
1179 if (regularized)
1180 return regularized.value()*EffToAbs::dswe_dsw(effToAbsParams_);
1181 }
1182
1183 return BaseLaw::dpcgw_dswe(swe, basicParams_)*EffToAbs::dswe_dsw(effToAbsParams_);
1184 }
1185
1189 template<bool enableRegularization = isRegularized()>
1190 Scalar dpcnw_dsw(const Scalar sw, const Scalar /*dummySw*/) const
1191 {
1192 const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_);
1193 if constexpr (enableRegularization)
1194 {
1195 const auto regularized = regularization_.dpcnw_dswe(swe);
1196 if (regularized)
1197 return regularized.value()*EffToAbs::dswe_dsw(effToAbsParams_);
1198 }
1199
1200 return BaseLaw::dpcnw_dswe(swe, basicParams_)*EffToAbs::dswe_dsw(effToAbsParams_);
1201 }
1202
1206 template<bool enableRegularization = isRegularized()>
1207 Scalar dpcgn_dst(const Scalar st, const Scalar /*dummySw*/) const
1208 {
1209 const auto ste = EffToAbs::stToSte(st, effToAbsParams_);
1210 if constexpr (enableRegularization)
1211 {
1212 const auto regularized = regularization_.dpcgn_dste(ste);
1213 if (regularized)
1214 return regularized.value()*EffToAbs::dswte_dst(effToAbsParams_);
1215 }
1216
1217 return BaseLaw::dpcgn_dste(ste, basicParams_)*EffToAbs::dste_dst(effToAbsParams_);
1218 }
1219
1225 template<bool enableRegularization = isRegularized()>
1226 Scalar krw(const Scalar sw, const Scalar sn) const
1227 {
1228 const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_);
1229 if constexpr (enableRegularization)
1230 {
1231 const auto regularized = regularization_.krw(swe);
1232 if (regularized)
1233 return regularized.value();
1234 }
1235
1236 return BaseLaw::krw(swe, basicParams_);
1237 }
1238
1244 template<bool enableRegularization = isRegularized()>
1245 Scalar krn(const Scalar sw, const Scalar sn) const
1246 {
1247 const auto swe = EffToAbs::swToSwe(sw, effToAbsParams_);
1248 const auto ste = EffToAbs::stToSte(sw + sn, effToAbsParams_);
1249 if constexpr (enableRegularization)
1250 {
1251 const auto regularized = regularization_.krn(swe, sn, ste);
1252 if (regularized)
1253 return regularized.value();
1254 }
1255
1256 return BaseLaw::krn(swe, sn, ste, basicParams_);
1257 }
1258
1264 template<bool enableRegularization = isRegularized()>
1265 Scalar krg(const Scalar sw, const Scalar sn) const
1266 {
1267 const auto ste = EffToAbs::stToSte(sw + sn, effToAbsParams_);
1268 if constexpr (enableRegularization)
1269 {
1270 const auto regularized = regularization_.krg(ste);
1271 if (regularized)
1272 return regularized.value();
1273 }
1274
1275 return BaseLaw::krg(ste, basicParams_);
1276 }
1277
1284 template<bool enableRegularization = isRegularized()>
1285 Scalar kr(const int phaseIdx, const Scalar sw, const Scalar sn) const
1286 {
1287 switch (phaseIdx)
1288 {
1289 case 0:
1290 return krw(sw, sn);
1291 case 1:
1292 return krn(sw, sn);
1293 case 2:
1294 return krg(sw, sn);
1295 }
1296 DUNE_THROW(Dune::InvalidStateException,
1297 "Invalid phase index ");
1298 }
1299
1304 template<bool enableRegularization = isRegularized()>
1305 Scalar dkrg_dst(const Scalar st) const
1306 {
1307 const auto ste = EffToAbs::stToSte(st, effToAbsParams_);
1308 if constexpr (enableRegularization)
1309 {
1310 const auto regularized = regularization_.dkrg_dste(ste);
1311 if (regularized)
1312 return regularized.value()*EffToAbs::dste_dst(effToAbsParams_);
1313 }
1314
1315 return BaseLaw::dkrg_dste(ste, basicParams_)*EffToAbs::dste_dst(effToAbsParams_);
1316 }
1317
1322 {
1323 return basicParams_ == o.basicParams_
1324 && effToAbsParams_ == o.effToAbsParams_
1325 && regularization_ == o.regularization_;
1326 }
1327
1332 static BasicParams makeBasicParams(const std::string& paramGroup)
1333 {
1334 return BaseLaw::template makeParams<Scalar>(paramGroup);
1335 }
1336
1341 { return basicParams_; }
1342
1347 static EffToAbsParams makeEffToAbsParams(const std::string& paramGroup)
1348 {
1349 return EffToAbs::template makeParams<Scalar>(paramGroup);
1350 }
1351
1356 { return effToAbsParams_; }
1357
1358private:
1359 BasicParams basicParams_;
1360 EffToAbsParams effToAbsParams_;
1361 Regularization regularization_;
1362};
1363
1368template<class Scalar>
1370
1375template<class Scalar>
1377
1378} // end namespace Dumux::FluidMatrix
1379
1380#endif
Provides 3rd order polynomial splines.
A wrapper that can either contain a valid Scalar or NaN.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
A tag to turn off regularization and it's overhead.
Wrapper type to combine an arbitrary number of different laws for fluid-matrix interaction (e....
Definition: brookscorey.hh:35
A 3rd order polynomial spline.
Definition: spline.hh:55
Definition: parkervangenuchten.hh:37
static Scalar snToSne(const Scalar sn, const Params< Scalar > &params)
Convert an absolute nonwetting saturation to an effective one.
Definition: parkervangenuchten.hh:181
static Scalar dswe_dsw(const Params< Scalar > &params)
Derivative of the effective saturation w.r.t. the absolute saturation.
Definition: parkervangenuchten.hh:152
static Params< Scalar > makeParams(const std::string &paramGroup)
Construct from a subgroup from the global parameter tree.
Definition: parkervangenuchten.hh:104
static Scalar steToSt(const Scalar ste, const Params< Scalar > &params)
Convert an effective wetting saturation to an absolute one.
Definition: parkervangenuchten.hh:211
static Scalar dste_dst(const Params< Scalar > &params)
Derivative of the effective saturation w.r.t. the absolute saturation.
Definition: parkervangenuchten.hh:225
static Scalar stToSte(const Scalar st, const Params< Scalar > &params)
Convert an absolute total liquid saturation to an effective one.
Definition: parkervangenuchten.hh:196
static Scalar swToSwe(const Scalar sw, const Params< Scalar > &params)
Convert an absolute wetting saturation to an effective one.
Definition: parkervangenuchten.hh:123
static Scalar dsw_dswe(const Params< Scalar > &params)
Derivative of the absolute saturation w.r.t. the effective saturation.
Definition: parkervangenuchten.hh:166
static Scalar sweToSw(const Scalar swe, const Params< Scalar > &params)
Convert an effective wetting saturation to an absolute one.
Definition: parkervangenuchten.hh:138
static Scalar dst_dste(const Params< Scalar > &params)
Derivative of the absolute saturation w.r.t. the effective saturation.
Definition: parkervangenuchten.hh:239
The parameter type.
Definition: parkervangenuchten.hh:46
Scalar swr() const
Return the residual wetting saturation.
Definition: parkervangenuchten.hh:54
Params(const Scalar swr=0.0, const Scalar snr=0.0, const Scalar sgr=0.0)
Definition: parkervangenuchten.hh:47
void setSwr(Scalar v)
Set the residual wetting saturation.
Definition: parkervangenuchten.hh:60
Scalar snr() const
Return the residual nonwetting saturation.
Definition: parkervangenuchten.hh:66
void setSgr(Scalar v)
Set the residual gas phase saturation.
Definition: parkervangenuchten.hh:84
bool operator==(const Params &p) const
Definition: parkervangenuchten.hh:87
void setSnr(Scalar v)
Set the residual nonwetting saturation.
Definition: parkervangenuchten.hh:72
Scalar sgr() const
Return the residual gas phase saturation.
Definition: parkervangenuchten.hh:78
Implementation of Parker/vanGenuchten's capillary pressure <-> saturation relation for three phases....
Definition: parkervangenuchten.hh:253
static Params< Scalar > makeParams(const std::string &paramGroup)
Construct from a subgroup from the global parameter tree.
Definition: parkervangenuchten.hh:329
static Scalar krn(const Scalar swe, const Scalar sn, const Scalar ste, const Params< Scalar > &params)
The relative permeability for the nonwetting phase after the Model of Parker et al....
Definition: parkervangenuchten.hh:488
static Scalar pcnw(Scalar swe, const Params< Scalar > &params)
The capillary pressure-saturation curve for the non-wettigng and wetting phase.
Definition: parkervangenuchten.hh:361
static Scalar krg(const Scalar ste, const Params< Scalar > &params)
The relative permeability for the nonwetting phase of the medium implied by van Genuchten's parameter...
Definition: parkervangenuchten.hh:523
static Scalar dpcnw_dswe(const Scalar swe, const Params< Scalar > &params)
Returns the partial derivative of the capillary pressure to the effective saturation.
Definition: parkervangenuchten.hh:426
static Scalar kr(const int phaseIdx, const Scalar swe, const Scalar sne, const Params< Scalar > &params)
The relative permeability for a phase.
Definition: parkervangenuchten.hh:561
static Scalar dkrg_dste(const Scalar ste, const Params< Scalar > &params)
The derivative of the relative permeability for the gas phase in regard to the total liquid saturatio...
Definition: parkervangenuchten.hh:541
static Scalar pcgw(Scalar swe, const Params< Scalar > &params)
The capillary pressure-saturation curve for the gas and wetting phase.
Definition: parkervangenuchten.hh:349
static Scalar dpcgn_dste(const Scalar ste, const Params< Scalar > &params)
Returns the partial derivative of the capillary pressure to the effective saturation.
Definition: parkervangenuchten.hh:441
static Scalar krw(const Scalar swe, const Params< Scalar > &params)
The relative permeability for the wetting phase of the medium implied by van Genuchten's parameteriza...
Definition: parkervangenuchten.hh:462
static Scalar pcAlpha(Scalar sne, const Params< Scalar > &params)
This function ensures a continuous transition from 2 to 3 phases and vice versa.
Definition: parkervangenuchten.hh:385
static Scalar dpcgw_dswe(const Scalar swe, const Params< Scalar > &params)
Returns the partial derivative of the capillary pressure to the effective saturation.
Definition: parkervangenuchten.hh:411
static Scalar pcgn(const Scalar ste, const Params< Scalar > &params)
The capillary pressure-saturation curve for the gas and nonwetting phase.
Definition: parkervangenuchten.hh:373
The parameter type.
Definition: parkervangenuchten.hh:272
Scalar betaGn() const
Definition: parkervangenuchten.hh:297
Scalar snr() const
Definition: parkervangenuchten.hh:291
Params(Scalar alpha, Scalar n, Scalar swr=0.0, Scalar snr=0.0, Scalar betaNw=1.0, Scalar betaGn=1.0, Scalar betaGw=1.0, bool regardSnr=false)
Definition: parkervangenuchten.hh:273
Scalar swr() const
Definition: parkervangenuchten.hh:288
void setAlpha(Scalar alpha)
Definition: parkervangenuchten.hh:280
void setSwr(Scalar swr)
Definition: parkervangenuchten.hh:289
void setM(Scalar m)
Definition: parkervangenuchten.hh:283
void setBetaGw(Scalar betaGw)
Definition: parkervangenuchten.hh:301
void setBetaNw(Scalar betaNw)
Definition: parkervangenuchten.hh:295
void setBetaGn(Scalar betaGn)
Definition: parkervangenuchten.hh:298
Scalar alpha() const
Definition: parkervangenuchten.hh:279
void setN(Scalar n)
Definition: parkervangenuchten.hh:286
bool regardSnrForKrn() const
Definition: parkervangenuchten.hh:303
Scalar betaNw() const
Definition: parkervangenuchten.hh:294
void setRegardSnrForKrn(bool v)
Definition: parkervangenuchten.hh:304
Scalar m() const
Definition: parkervangenuchten.hh:282
Scalar n() const
Definition: parkervangenuchten.hh:285
void setSnr(Scalar swr)
Definition: parkervangenuchten.hh:292
Scalar betaGw() const
Definition: parkervangenuchten.hh:300
bool operator==(const Params &p) const
Definition: parkervangenuchten.hh:306
A regularization for the ParkerVanGenuchten3PRegularization material law.
Definition: parkervangenuchten.hh:602
void init(const MaterialLaw *m, const BaseParams &bp, const EffToAbsParams &etap, const Params< Scalar > &p)
Definition: parkervangenuchten.hh:723
OptionalScalar< Scalar > pcgn(Scalar ste) const
The regularized capillary pressure-saturation curve for the gas and nonwetting phase regularized part...
Definition: parkervangenuchten.hh:822
OptionalScalar< Scalar > pcnw(Scalar swe) const
The regularized capillary pressure-saturation curve for the nonwetting and wetting phase regularized ...
Definition: parkervangenuchten.hh:789
OptionalScalar< Scalar > pcgw(Scalar swe) const
The regularized capillary pressure-saturation curve for the gas and wetting phase regularized part:
Definition: parkervangenuchten.hh:756
bool operator==(const ParkerVanGenuchten3PRegularization &o) const
Equality comparison with another instance.
Definition: parkervangenuchten.hh:739
OptionalScalar< Scalar > kr(const int phaseIdx, const Scalar swe, const Scalar sne) const
The relative permeability for a phase.
Definition: parkervangenuchten.hh:938
void init(const MaterialLaw *m, const std::string &paramGroup)
Initialize the spline.
Definition: parkervangenuchten.hh:709
OptionalScalar< Scalar > krn(Scalar swe, const Scalar sn, Scalar ste) const
The regularized relative permeability for the nonwetting phase.
Definition: parkervangenuchten.hh:881
OptionalScalar< Scalar > krg(const Scalar ste) const
The regularized relative permeability for the gas phase.
Definition: parkervangenuchten.hh:897
OptionalScalar< Scalar > pcAlpha(Scalar sne) const
This function ensures a continuous transition from 2 to 3 phases and vice versa.
Definition: parkervangenuchten.hh:854
OptionalScalar< Scalar > krw(const Scalar swe) const
The regularized relative permeability for the wetting phase.
Definition: parkervangenuchten.hh:864
Regularization parameters.
Definition: parkervangenuchten.hh:609
Scalar krnLowSwe() const
Threshold saturation below which the relative permeability of the nonwetting phase gets regularized.
Definition: parkervangenuchten.hh:651
void setKrnLowSwe(Scalar krnLowSwe)
Set the threshold saturation below which the relative permeability of the nonwetting phase gets regul...
Definition: parkervangenuchten.hh:644
Scalar pcLowSwe() const
Threshold saturation below which the capillary pressure is regularized.
Definition: parkervangenuchten.hh:622
void setPcLowSwe(Scalar pcLowSwe)
Set the threshold saturation below which the capillary pressure is regularized.
Definition: parkervangenuchten.hh:616
Scalar pcHighSwe() const
Threshold saturation above which the capillary pressure is regularized.
Definition: parkervangenuchten.hh:637
void setConstRegularization(const bool input)
Choose whether to use a constant value for regularization of the pc-S curves or not.
Definition: parkervangenuchten.hh:688
void setKrwHighSwe(Scalar krwHighSwe)
Set the threshold saturation above which the relative permeability of the wetting phase gets regulari...
Definition: parkervangenuchten.hh:672
bool constRegularization() const
Returns whether to use a constant value for regularization of the pc-S curves or not.
Definition: parkervangenuchten.hh:695
Scalar krwHighSwe() const
Threshold saturation above which the relative permeability of the wetting phase gets regularized.
Definition: parkervangenuchten.hh:679
void setPcHighSwe(Scalar pcHighSwe)
Set the threshold saturation above which the capillary pressure is regularized.
Definition: parkervangenuchten.hh:628
void setKrgLowSte(Scalar krgLowSte)
Set the threshold saturation below which the relative permeability of the nonwetting phase gets regul...
Definition: parkervangenuchten.hh:658
Scalar krgLowSte() const
Threshold saturation below which the relative permeability of the nonwetting phase gets regularized.
Definition: parkervangenuchten.hh:665
Definition: parkervangenuchten.hh:1051
typename Regularization::template Params< Scalar > RegularizationParams
Definition: parkervangenuchten.hh:1058
Scalar kr(const int phaseIdx, const Scalar sw, const Scalar sn) const
The relative permeability for the nonwetting phase.
Definition: parkervangenuchten.hh:1285
ParkerVanGenuchtenMaterialLaw(const std::string &paramGroup)
Construct from a subgroup from the global parameter tree.
Definition: parkervangenuchten.hh:1078
Scalar pcnw(const Scalar sw, const Scalar) const
The capillary pressure-saturation curve for the nonwetting and wetting phase.
Definition: parkervangenuchten.hh:1121
typename EffToAbsPolicy::template Params< Scalar > EffToAbsParams
Definition: parkervangenuchten.hh:1057
const EffToAbsParams & effToAbsParams() const
Return the parameters of the EffToAbs policy.
Definition: parkervangenuchten.hh:1355
ParkerVanGenuchtenMaterialLaw()=delete
Deleted default constructor (so we are never in an undefined state)
Scalar dpcnw_dsw(const Scalar sw, const Scalar) const
The partial derivative of the capillary pressure w.r.t. the saturation.
Definition: parkervangenuchten.hh:1190
ScalarType Scalar
Definition: parkervangenuchten.hh:1054
Scalar krn(const Scalar sw, const Scalar sn) const
The relative permeability for the nonwetting phase.
Definition: parkervangenuchten.hh:1245
static BasicParams makeBasicParams(const std::string &paramGroup)
Create the base law's parameters using input file parameters.
Definition: parkervangenuchten.hh:1332
Scalar pcgn(const Scalar sw, const Scalar sn) const
The capillary pressure-saturation curve for the gas and nonwetting phase.
Definition: parkervangenuchten.hh:1140
Scalar krw(const Scalar sw, const Scalar sn) const
The relative permeability for the wetting phase.
Definition: parkervangenuchten.hh:1226
Scalar dpcgw_dsw(const Scalar sw, const Scalar) const
The partial derivative of the capillary pressure w.r.t. the saturation.
Definition: parkervangenuchten.hh:1173
Scalar krg(const Scalar sw, const Scalar sn) const
The relative permeability for the nonwetting phase.
Definition: parkervangenuchten.hh:1265
Scalar pcgw(const Scalar sw, const Scalar) const
The capillary pressure-saturation curve for the gas and wetting phase.
Definition: parkervangenuchten.hh:1104
static EffToAbsParams makeEffToAbsParams(const std::string &paramGroup)
Create the parameters of the EffToAbs policy using input file parameters.
Definition: parkervangenuchten.hh:1347
const BasicParams & basicParams() const
Return the base law's parameters.
Definition: parkervangenuchten.hh:1340
Scalar dkrg_dst(const Scalar st) const
The derivative of the relative permeability for the nonwetting phase w.r.t. saturation.
Definition: parkervangenuchten.hh:1305
ParkerVanGenuchtenMaterialLaw(const BasicParams &baseParams, const EffToAbsParams &effToAbsParams={}, const RegularizationParams &regParams={})
Construct from parameter structs.
Definition: parkervangenuchten.hh:1090
static constexpr bool isRegularized()
Return whether this law is regularized.
Definition: parkervangenuchten.hh:1065
bool operator==(const ParkerVanGenuchtenMaterialLaw &o) const
Equality comparison with another instance.
Definition: parkervangenuchten.hh:1321
typename BaseLaw::template Params< Scalar > BasicParams
Definition: parkervangenuchten.hh:1056
Scalar pcAlpha(const Scalar, const Scalar sn) const
This function ensures a continuous transition from 2 to 3 phases and vice versa.
Definition: parkervangenuchten.hh:1157
Scalar dpcgn_dst(const Scalar st, const Scalar) const
The partial derivative of the capillary pressure w.r.t. the saturation.
Definition: parkervangenuchten.hh:1207
EffToAbsPolicy EffToAbs
Definition: parkervangenuchten.hh:1060
Adapter to inherit from, allowing the inheriting class to be wrapped by the makeFluidMatrixInteractio...
Definition: fluidmatrixinteraction.hh:67