LCOV - code coverage report
Current view: top level - src/utils - PorousFlowBroadbridgeWhite.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 54 55 98.2 %
Date: 2017-11-21 14:47:27 Functions: 9 9 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       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             : #include "PorousFlowBroadbridgeWhite.h"
       9             : #include "libmesh/utility.h"
      10             : 
      11             : namespace PorousFlowBroadbridgeWhite
      12             : {
      13             : Real
      14     1555454 : 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     1555454 :   if (z < 1.0)
      67             :   {
      68             :     /* series near 0 */
      69      221862 :     p = std::sqrt(2.0 * (2.7182818284590452353602874713526625 * z + 1.0));
      70      221862 :     w = -1.0 + p * (1.0 + p * (-0.333333333333333333333 + p * 0.152777777777777777777777));
      71             :   }
      72             :   else
      73     1333592 :     w = std::log(z); /* asymptotic */
      74     1555454 :   if (z > 3.0)
      75     1288912 :     w -= std::log(w); /* useful? */
      76     8996924 :   for (unsigned int i = 0; i < 10; ++i)
      77             :   {
      78             :     /* Halley iteration */
      79     5276189 :     e = std::exp(w);
      80     5276189 :     t = w * e - z;
      81     5276189 :     p = w + 1.0;
      82     5276189 :     t /= e * p - 0.5 * (p + 1.0) * t / p;
      83     5276189 :     w -= t;
      84     5276189 :     if (std::abs(t) < eps * (1.0 + std::abs(w)))
      85     1555454 :       return w; /* rel-abs error */
      86             :   }
      87             :   /* should never get here */
      88           0 :   mooseError("LambertW: No convergence at z= ", z, "\n");
      89             : }
      90             : 
      91             : Real
      92      630602 : effectiveSaturation(Real pressure, Real c, Real sn, Real ss, Real las)
      93             : {
      94      630602 :   if (pressure >= 0.0)
      95             :     return 1.0;
      96      624654 :   const Real x = (c - 1.0) * std::exp(c - 1.0 - c * pressure / las);
      97      624654 :   const Real th = c / (1.0 + LambertW(x)); // use branch 0 for positive x
      98      624654 :   return sn + (ss - sn) * th;
      99             : }
     100             : 
     101             : Real
     102      626602 : dEffectiveSaturation(Real pressure, Real c, Real sn, Real ss, Real las)
     103             : {
     104      626602 :   if (pressure >= 0.0)
     105             :     return 0.0;
     106      620654 :   const Real x = (c - 1.0) * std::exp(c - 1.0 - c * pressure / las);
     107      620654 :   const Real lamw = LambertW(x);
     108     1861962 :   return (ss - sn) * Utility::pow<2>(c) / las * lamw / Utility::pow<3>(1.0 + lamw);
     109             : }
     110             : 
     111             : Real
     112      313106 : d2EffectiveSaturation(Real pressure, Real c, Real sn, Real ss, Real las)
     113             : {
     114      313106 :   if (pressure >= 0.0)
     115             :     return 0.0;
     116      310146 :   const Real x = (c - 1.0) * std::exp(c - 1.0 - c * pressure / las);
     117      310146 :   const Real lamw = LambertW(x);
     118      930438 :   return -(ss - sn) * Utility::pow<3>(c) / Utility::pow<2>(las) * lamw * (1.0 - 2.0 * lamw) /
     119      620292 :          Utility::pow<5>(1.0 + lamw);
     120             : }
     121             : 
     122             : Real
     123      313502 : relativePermeability(Real s, Real c, Real sn, Real ss, Real kn, Real ks)
     124             : {
     125      313502 :   if (s <= sn)
     126             :     return kn;
     127             : 
     128      313500 :   if (s >= ss)
     129             :     return ks;
     130             : 
     131      310510 :   const Real coef = (ks - kn) * (c - 1.0);
     132      310510 :   const Real th = (s - sn) / (ss - sn);
     133      310510 :   const Real krel = kn + coef * Utility::pow<2>(th) / (c - th);
     134      310510 :   return krel;
     135             : }
     136             : 
     137             : Real
     138      313502 : dRelativePermeability(Real s, Real c, Real sn, Real ss, Real kn, Real ks)
     139             : {
     140      313502 :   if (s <= sn)
     141             :     return 0.0;
     142             : 
     143      313500 :   if (s >= ss)
     144             :     return 0.0;
     145             : 
     146      310510 :   const Real coef = (ks - kn) * (c - 1.0);
     147      310510 :   const Real th = (s - sn) / (ss - sn);
     148      621020 :   const Real krelp = coef * (2.0 * th / (c - th) + Utility::pow<2>(th) / Utility::pow<2>(c - th));
     149      310510 :   return krelp / (ss - sn);
     150             : }
     151             : 
     152             : Real
     153           6 : d2RelativePermeability(Real s, Real c, Real sn, Real ss, Real kn, Real ks)
     154             : {
     155           6 :   if (s <= sn)
     156             :     return 0.0;
     157             : 
     158           4 :   if (s >= ss)
     159             :     return 0.0;
     160             : 
     161           2 :   const Real coef = (ks - kn) * (c - 1.0);
     162           2 :   const Real th = (s - sn) / (ss - sn);
     163           6 :   const Real krelpp = coef * (2.0 / (c - th) + 4.0 * th / Utility::pow<2>(c - th) +
     164           4 :                               2.0 * Utility::pow<2>(th) / Utility::pow<3>(c - th));
     165           2 :   return krelpp / Utility::pow<2>(ss - sn);
     166             : }
     167        2499 : }

Generated by: LCOV version 1.11