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
9 #include "libmesh/utility.h"
10
12 {
13 Real
14 LambertW(Real z)
15 {
16  /* Lambert W function.
17  Was ~/C/LambertW.c written K M Briggs Keith dot Briggs at bt dot com 97 May 21.
18  Revised KMB 97 Nov 20; 98 Feb 11, Nov 24, Dec 28; 99 Jan 13; 00 Feb 23; 01 Apr 09
19
20  Computes Lambert W function, principal branch.
21  See LambertW1.c for -1 branch.
22
23  Returned value W(z) satisfies W(z)*exp(W(z))=z
24  test data...
25  W(1)= 0.5671432904097838730
26  W(2)= 0.8526055020137254914
27  W(20)=2.2050032780240599705
28  To solve (a+b*R)*exp(-c*R)-d=0 for R, use
29  R=-(b*W(-exp(-a*c/b)/b*d*c)+a*c)/b/c
30
31  Test:
32  gcc -DTESTW LambertW.c -o LambertW -lm && LambertW
33  Library:
34  gcc -O3 -c LambertW.c
35
36  Modified trivially by Andy to use MOOSE things
37  */
38  mooseAssert(z > 0, "LambertW function in RichardsSeff1BWsmall called with negative argument");
39
40  const Real eps = 4.0e-16; //, em1=0.3678794411714423215955237701614608;
41  Real p, e, t, w;
42
43  /* Uncomment this stuff is you ever need to call with a negative argument
44  if (z < -em1)
45  mooseError("LambertW: bad argument ", z, "\n");
46
47  if (0.0 == z)
48  return 0.0;
49  if (z < -em1+1e-4)
50  {
51  // series near -em1 in sqrt(q)
52  Real q=z+em1,r=std::sqrt(q),q2=q*q,q3=q2*q;
53  return
54  -1.0
55  +2.331643981597124203363536062168*r
56  -1.812187885639363490240191647568*q
57  +1.936631114492359755363277457668*r*q
58  -2.353551201881614516821543561516*q2
59  +3.066858901050631912893148922704*r*q2
60  -4.175335600258177138854984177460*q3
61  +5.858023729874774148815053846119*r*q3
62  -8.401032217523977370984161688514*q3*q; // error approx 1e-16
63  }
64  */
65  /* initial approx for iteration... */
66  if (z < 1.0)
67  {
68  /* series near 0 */
69  p = std::sqrt(2.0 * (2.7182818284590452353602874713526625 * z + 1.0));
70  w = -1.0 + p * (1.0 + p * (-0.333333333333333333333 + p * 0.152777777777777777777777));
71  }
72  else
73  w = std::log(z); /* asymptotic */
74  if (z > 3.0)
75  w -= std::log(w); /* useful? */
76  for (unsigned int i = 0; i < 10; ++i)
77  {
78  /* Halley iteration */
79  e = std::exp(w);
80  t = w * e - z;
81  p = w + 1.0;
82  t /= e * p - 0.5 * (p + 1.0) * t / p;
83  w -= t;
84  if (std::abs(t) < eps * (1.0 + std::abs(w)))
85  return w; /* rel-abs error */
86  }
87  /* should never get here */
88  mooseError("LambertW: No convergence at z= ", z, "\n");
89 }
90
91 Real
92 effectiveSaturation(Real pressure, Real c, Real sn, Real ss, Real las)
93 {
94  if (pressure >= 0.0)
95  return 1.0;
96  const Real x = (c - 1.0) * std::exp(c - 1.0 - c * pressure / las);
97  const Real th = c / (1.0 + LambertW(x)); // use branch 0 for positive x
98  return sn + (ss - sn) * th;
99 }
100
101 Real
102 dEffectiveSaturation(Real pressure, Real c, Real sn, Real ss, Real las)
103 {
104  if (pressure >= 0.0)
105  return 0.0;
106  const Real x = (c - 1.0) * std::exp(c - 1.0 - c * pressure / las);
107  const Real lamw = LambertW(x);
108  return (ss - sn) * Utility::pow<2>(c) / las * lamw / Utility::pow<3>(1.0 + lamw);
109 }
110
111 Real
112 d2EffectiveSaturation(Real pressure, Real c, Real sn, Real ss, Real las)
113 {
114  if (pressure >= 0.0)
115  return 0.0;
116  const Real x = (c - 1.0) * std::exp(c - 1.0 - c * pressure / las);
117  const Real lamw = LambertW(x);
118  return -(ss - sn) * Utility::pow<3>(c) / Utility::pow<2>(las) * lamw * (1.0 - 2.0 * lamw) /
119  Utility::pow<5>(1.0 + lamw);
120 }
121
122 Real
123 relativePermeability(Real s, Real c, Real sn, Real ss, Real kn, Real ks)
124 {
125  if (s <= sn)
126  return kn;
127
128  if (s >= ss)
129  return ks;
130
131  const Real coef = (ks - kn) * (c - 1.0);
132  const Real th = (s - sn) / (ss - sn);
133  const Real krel = kn + coef * Utility::pow<2>(th) / (c - th);
134  return krel;
135 }
136
137 Real
138 dRelativePermeability(Real s, Real c, Real sn, Real ss, Real kn, Real ks)
139 {
140  if (s <= sn)
141  return 0.0;
142
143  if (s >= ss)
144  return 0.0;
145
146  const Real coef = (ks - kn) * (c - 1.0);
147  const Real th = (s - sn) / (ss - sn);
148  const Real krelp = coef * (2.0 * th / (c - th) + Utility::pow<2>(th) / Utility::pow<2>(c - th));
149  return krelp / (ss - sn);
150 }
151
152 Real
153 d2RelativePermeability(Real s, Real c, Real sn, Real ss, Real kn, Real ks)
154 {
155  if (s <= sn)
156  return 0.0;
157
158  if (s >= ss)
159  return 0.0;
160
161  const Real coef = (ks - kn) * (c - 1.0);
162  const Real th = (s - sn) / (ss - sn);
163  const Real krelpp = coef * (2.0 / (c - th) + 4.0 * th / Utility::pow<2>(c - th) +
164  2.0 * Utility::pow<2>(th) / Utility::pow<3>(c - th));
165  return krelpp / Utility::pow<2>(ss - sn);
166 }
167 }
Real d2RelativePermeability(Real s, Real c, Real sn, Real ss, Real kn, Real ks)
Second derivative of relative permeability with respect to saturation.
Real dEffectiveSaturation(Real pc, Real c, Real sn, Real ss, Real las)
Derivative of effective saturation wrt capillary pressure.
Real LambertW(Real z)
Provides the Lambert W function, which satisfies W(z) * exp(W(z)) = z.
Real d2EffectiveSaturation(Real pc, Real c, Real sn, Real ss, Real las)
Second derivative of effective saturation wrt capillary pressure.
Real effectiveSaturation(Real pc, Real c, Real sn, Real ss, Real las)
Effective saturation as a function of capillary pressure If pc>=0 this will yield 1...
Real dRelativePermeability(Real s, Real c, Real sn, Real ss, Real kn, Real ks)
Derivative of relative permeability with respect to saturation.
Real relativePermeability(Real s, Real c, Real sn, Real ss, Real kn, Real ks)
Relative permeability as a function of saturation.
const std::string pressure
Definition: NS.h:24
Broadbridge-White version of relative permeability, and effective saturation as a function of capilla...