www.mooseframework.org
RichardsSeff1BWsmall.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
18 
21 {
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  params.addRangeCheckedParam<Real>(
31  "Ss",
32  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.");
38  "C", "C > 1", "BW's C parameter. Must be > 1. Typical value would be 1.05.");
39  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  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  return params;
50 }
51 
53  : RichardsSeff(parameters),
54  _sn(getParam<Real>("Sn")),
55  _ss(getParam<Real>("Ss")),
56  _c(getParam<Real>("C")),
57  _las(getParam<Real>("las"))
58 {
59  if (_ss <= _sn)
60  mooseError("In BW effective saturation Sn set to ",
61  _sn,
62  " and Ss set to ",
63  _ss,
64  " but these must obey Ss > Sn");
65 }
66 
67 Real
68 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  if (z < 1.0)
122  {
123  /* series near 0 */
124  p = std::sqrt(2.0 * (2.7182818284590452353602874713526625 * z + 1.0));
125  w = -1.0 + p * (1.0 + p * (-0.333333333333333333333 + p * 0.152777777777777777777777));
126  }
127  else
128  w = std::log(z); /* asymptotic */
129  if (z > 3.0)
130  w -= std::log(w); /* useful? */
131  for (i = 0; i < 10; i++)
132  {
133  /* Halley iteration */
134  e = std::exp(w);
135  t = w * e - z;
136  p = w + 1.0;
137  t /= e * p - 0.5 * (p + 1.0) * t / p;
138  w -= t;
139  if (std::abs(t) < eps * (1.0 + std::abs(w)))
140  return w; /* rel-abs error */
141  }
142  /* should never get here */
143  mooseError("LambertW: No convergence at z= ", z, "\n");
144 }
145 
146 Real
147 RichardsSeff1BWsmall::seff(std::vector<const VariableValue *> p, unsigned int qp) const
148 {
149  Real pp = (*p[0])[qp];
150  if (pp >= 0)
151  return 1.0;
152 
153  Real x = (_c - 1.0) * std::exp(_c - 1 - _c * pp / _las);
154  Real th = _c / (1.0 + LambertW(x)); // use branch 0 for positive x
155  return _sn + (_ss - _sn) * th;
156 }
157 
158 void
159 RichardsSeff1BWsmall::dseff(std::vector<const VariableValue *> p,
160  unsigned int qp,
161  std::vector<Real> & result) const
162 {
163  result[0] = 0.0;
164 
165  Real pp = (*p[0])[qp];
166  if (pp >= 0)
167  return;
168 
169  Real x = (_c - 1) * std::exp(_c - 1.0 - _c * pp / _las);
170  Real lamw = LambertW(x);
171  result[0] = Utility::pow<2>(_c) / _las * lamw / Utility::pow<3>(1 + lamw);
172 }
173 
174 void
175 RichardsSeff1BWsmall::d2seff(std::vector<const VariableValue *> p,
176  unsigned int qp,
177  std::vector<std::vector<Real>> & result) const
178 {
179  result[0][0] = 0.0;
180 
181  Real pp = (*p[0])[qp];
182  if (pp >= 0)
183  return;
184 
185  Real x = (_c - 1) * std::exp(_c - 1 - _c * pp / _las);
186  Real lamw = LambertW(x);
187  result[0][0] = -Utility::pow<3>(_c) / Utility::pow<2>(_las) * lamw * (1.0 - 2.0 * lamw) /
188  Utility::pow<5>(1 + lamw);
189 }
RichardsSeff1BWsmall(const InputParameters &parameters)
void dseff(std::vector< const VariableValue *> p, unsigned int qp, std::vector< Real > &result) const
derivative of effective saturation as a function of porepressure
void addRequiredRangeCheckedParam(const std::string &name, const std::string &parsed_function, const std::string &doc_string)
const Real eps
Base class for effective saturation as a function of porepressure(s) The functions seff...
Definition: RichardsSeff.h:18
static InputParameters validParams()
Definition: RichardsSeff.C:15
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
"Broadbridge-White" form of effective saturation for Kn small as a function of porepressure (not capi...
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
const std::vector< double > x
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 _las
BW&#39;s lambda_s parameter multiplied by (fluiddensity*gravity)
static InputParameters validParams()
Real _sn
BW&#39;s initial effective saturation.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real _c
BW&#39;s C parameter.
Real LambertW(const double z) const
LambertW function, returned value satisfies W(z)*exp(W(z))=z.
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
Real seff(std::vector< const VariableValue *> p, unsigned int qp) const
effective saturation as a function of porepressure
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
Real _ss
Effective saturation where porepressure = 0.
registerMooseObject("RichardsApp", RichardsSeff1BWsmall)