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