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