Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : // "Broadbridge-White" form of effective saturation for Kn small (P Broadbridge and I White
11 : // ``Constant rate rainfall infiltration: A versatile nonlinear model 1. Analytic Solution'', Water
12 : // Resources Research 24 (1988) 145-154)
13 : //
14 : #include "RichardsSeff1BWsmall.h"
15 : #include "libmesh/utility.h"
16 :
17 : registerMooseObject("RichardsApp", RichardsSeff1BWsmall);
18 :
19 : InputParameters
20 25 : RichardsSeff1BWsmall::validParams()
21 : {
22 25 : InputParameters params = RichardsSeff::validParams();
23 50 : params.addRequiredRangeCheckedParam<Real>(
24 : "Sn",
25 : "Sn >= 0",
26 : "Low saturation. This must be < Ss, and non-negative. This is BW's "
27 : "initial effective saturation, below which effective saturation never goes "
28 : "in their simulations/models. If Kn=0 then Sn is the immobile saturation. "
29 : "This form of effective saturation is only correct for Kn small.");
30 75 : params.addRangeCheckedParam<Real>(
31 : "Ss",
32 50 : 1.0,
33 : "Ss <= 1",
34 : "High saturation. This must be > Sn and <= 1. Effective saturation "
35 : "where porepressure = 0. Effective saturation never exceeds this "
36 : "value in BW's simulations/models.");
37 50 : params.addRequiredRangeCheckedParam<Real>(
38 : "C", "C > 1", "BW's C parameter. Must be > 1. Typical value would be 1.05.");
39 50 : params.addRequiredRangeCheckedParam<Real>("las",
40 : "las > 0",
41 : "BW's lambda_s parameter multiplied "
42 : "by (fluiddensity*gravity). Must be "
43 : "> 0. Typical value would be 1E5");
44 25 : params.addClassDescription("Broadbridge-white form of effective saturation for negligable Kn. "
45 : "Then porepressure = -las*( (1-th)/th - (1/c)Ln((C-th)/((C-1)th))), "
46 : "for th = (Seff - Sn)/(Ss - Sn). A Lambert-W function must be "
47 : "evaluated to express Seff in terms of porepressure, which can be "
48 : "expensive");
49 25 : return params;
50 0 : }
51 :
52 13 : RichardsSeff1BWsmall::RichardsSeff1BWsmall(const InputParameters & parameters)
53 : : RichardsSeff(parameters),
54 13 : _sn(getParam<Real>("Sn")),
55 26 : _ss(getParam<Real>("Ss")),
56 26 : _c(getParam<Real>("C")),
57 39 : _las(getParam<Real>("las"))
58 : {
59 13 : if (_ss <= _sn)
60 1 : mooseError("In BW effective saturation Sn set to ",
61 1 : _sn,
62 : " and Ss set to ",
63 1 : _ss,
64 : " but these must obey Ss > Sn");
65 12 : }
66 :
67 : Real
68 5501270 : RichardsSeff1BWsmall::LambertW(const Real z) const
69 : {
70 : /* Lambert W function.
71 : Was ~/C/LambertW.c written K M Briggs Keith dot Briggs at bt dot com 97 May 21.
72 : Revised KMB 97 Nov 20; 98 Feb 11, Nov 24, Dec 28; 99 Jan 13; 00 Feb 23; 01 Apr 09
73 :
74 : Computes Lambert W function, principal branch.
75 : See LambertW1.c for -1 branch.
76 :
77 : Returned value W(z) satisfies W(z)*exp(W(z))=z
78 : test data...
79 : W(1)= 0.5671432904097838730
80 : W(2)= 0.8526055020137254914
81 : W(20)=2.2050032780240599705
82 : To solve (a+b*R)*exp(-c*R)-d=0 for R, use
83 : R=-(b*W(-exp(-a*c/b)/b*d*c)+a*c)/b/c
84 :
85 : Test:
86 : gcc -DTESTW LambertW.c -o LambertW -lm && LambertW
87 : Library:
88 : gcc -O3 -c LambertW.c
89 :
90 : Modified trially by Andy to use MOOSE things
91 : */
92 : mooseAssert(z > 0, "LambertW function in RichardsSeff1BWsmall called with negative argument");
93 :
94 : int i;
95 : const Real eps = 4.0e-16; //, em1=0.3678794411714423215955237701614608;
96 : Real p, e, t, w;
97 :
98 : /* Uncomment this stuff is you ever need to call with a negative argument
99 : if (z < -em1)
100 : mooseError("LambertW: bad argument ", z, "\n");
101 :
102 : if (0.0 == z)
103 : return 0.0;
104 : if (z < -em1+1e-4)
105 : {
106 : // series near -em1 in sqrt(q)
107 : Real q=z+em1,r=std::sqrt(q),q2=q*q,q3=q2*q;
108 : return
109 : -1.0
110 : +2.331643981597124203363536062168*r
111 : -1.812187885639363490240191647568*q
112 : +1.936631114492359755363277457668*r*q
113 : -2.353551201881614516821543561516*q2
114 : +3.066858901050631912893148922704*r*q2
115 : -4.175335600258177138854984177460*q3
116 : +5.858023729874774148815053846119*r*q3
117 : -8.401032217523977370984161688514*q3*q; // error approx 1e-16
118 : }
119 : */
120 : /* initial approx for iteration... */
121 5501270 : if (z < 1.0)
122 : {
123 : /* series near 0 */
124 226264 : p = std::sqrt(2.0 * (2.7182818284590452353602874713526625 * z + 1.0));
125 226264 : w = -1.0 + p * (1.0 + p * (-0.333333333333333333333 + p * 0.152777777777777777777777));
126 : }
127 : else
128 5275006 : w = std::log(z); /* asymptotic */
129 5501270 : if (z > 3.0)
130 5239606 : w -= std::log(w); /* useful? */
131 17389859 : for (i = 0; i < 10; i++)
132 : {
133 : /* Halley iteration */
134 17389859 : e = std::exp(w);
135 17389859 : t = w * e - z;
136 17389859 : p = w + 1.0;
137 17389859 : t /= e * p - 0.5 * (p + 1.0) * t / p;
138 17389859 : w -= t;
139 17389859 : if (std::abs(t) < eps * (1.0 + std::abs(w)))
140 5501270 : return w; /* rel-abs error */
141 : }
142 : /* should never get here */
143 0 : mooseError("LambertW: No convergence at z= ", z, "\n");
144 : }
145 :
146 : Real
147 3662140 : RichardsSeff1BWsmall::seff(std::vector<const VariableValue *> p, unsigned int qp) const
148 : {
149 3662140 : Real pp = (*p[0])[qp];
150 3662140 : if (pp >= 0)
151 : return 1.0;
152 :
153 3643142 : Real x = (_c - 1.0) * std::exp(_c - 1 - _c * pp / _las);
154 3643142 : Real th = _c / (1.0 + LambertW(x)); // use branch 0 for positive x
155 3643142 : return _sn + (_ss - _sn) * th;
156 : }
157 :
158 : void
159 1173182 : RichardsSeff1BWsmall::dseff(std::vector<const VariableValue *> p,
160 : unsigned int qp,
161 : std::vector<Real> & result) const
162 : {
163 1173182 : result[0] = 0.0;
164 :
165 1173182 : Real pp = (*p[0])[qp];
166 1173182 : if (pp >= 0)
167 : return;
168 :
169 1167464 : Real x = (_c - 1) * std::exp(_c - 1.0 - _c * pp / _las);
170 1167464 : Real lamw = LambertW(x);
171 1167464 : result[0] = Utility::pow<2>(_c) / _las * lamw / Utility::pow<3>(1 + lamw);
172 : }
173 :
174 : void
175 696382 : RichardsSeff1BWsmall::d2seff(std::vector<const VariableValue *> p,
176 : unsigned int qp,
177 : std::vector<std::vector<Real>> & result) const
178 : {
179 696382 : result[0][0] = 0.0;
180 :
181 696382 : Real pp = (*p[0])[qp];
182 696382 : if (pp >= 0)
183 : return;
184 :
185 690664 : Real x = (_c - 1) * std::exp(_c - 1 - _c * pp / _las);
186 690664 : Real lamw = LambertW(x);
187 690664 : result[0][0] = -Utility::pow<3>(_c) / Utility::pow<2>(_las) * lamw * (1.0 - 2.0 * lamw) /
188 690664 : Utility::pow<5>(1 + lamw);
189 : }
|
