LCOV - code coverage report
Current view: top level - src/utils - PorousFlowVanGenuchten.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 64 64 100.0 %
Date: 2017-11-21 14:47:27 Functions: 11 11 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 "PorousFlowVanGenuchten.h"
       9             : #include "libmesh/utility.h"
      10             : 
      11             : namespace PorousFlowVanGenuchten
      12             : {
      13             : Real
      14     5989812 : effectiveSaturation(Real p, Real alpha, Real m)
      15             : {
      16             :   Real n, seff;
      17             : 
      18     5989812 :   if (p >= 0.0)
      19             :     return 1.0;
      20             :   else
      21             :   {
      22      635734 :     n = 1.0 / (1.0 - m);
      23      635734 :     seff = 1.0 + std::pow(-alpha * p, n);
      24      635734 :     return std::pow(seff, -m);
      25             :   }
      26             : }
      27             : 
      28             : Real
      29     5872682 : dEffectiveSaturation(Real p, Real alpha, Real m)
      30             : {
      31     5872682 :   if (p >= 0.0)
      32             :     return 0.0;
      33             :   else
      34             :   {
      35      631992 :     Real n = 1.0 / (1.0 - m);
      36      631992 :     Real inner = 1.0 + std::pow(-alpha * p, n);
      37      631992 :     Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
      38      631992 :     Real dseff_dp = -m * std::pow(inner, -m - 1) * dinner_dp;
      39      631992 :     return dseff_dp;
      40             :   }
      41             : }
      42             : 
      43             : Real
      44     2734650 : d2EffectiveSaturation(Real p, Real alpha, Real m)
      45             : {
      46     2734650 :   if (p >= 0.0)
      47             :     return 0.0;
      48             :   else
      49             :   {
      50      300189 :     Real n = 1.0 / (1.0 - m);
      51      300189 :     Real inner = 1.0 + std::pow(-alpha * p, n);
      52      300189 :     Real dinner_dp = -n * alpha * std::pow(-alpha * p, n - 1.0);
      53      300189 :     Real d2inner_dp2 = n * (n - 1.0) * alpha * alpha * std::pow(-alpha * p, n - 2.0);
      54      300189 :     Real d2seff_dp2 = m * (m + 1.0) * std::pow(inner, -m - 2.0) * std::pow(dinner_dp, 2.0) -
      55      300189 :                       m * std::pow(inner, -m - 1.0) * d2inner_dp2;
      56      300189 :     return d2seff_dp2;
      57             :   }
      58             : }
      59             : 
      60             : Real
      61       65161 : capillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
      62             : {
      63       65161 :   if (seff >= 1.0)
      64             :     return 0.0;
      65       60502 :   else if (seff <= 0.0)
      66        5230 :     return pc_max;
      67             :   else
      68             :   {
      69       55272 :     Real a = std::pow(seff, -1.0 / m) - 1.0;
      70      110544 :     return std::min(std::pow(a, 1.0 - m) / alpha, pc_max);
      71             :   }
      72             : }
      73             : 
      74             : Real
      75       61485 : dCapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
      76             : {
      77       61485 :   if (seff <= 0.0 || seff >= 1.0)
      78             :     return 0.0;
      79             :   else
      80             :   {
      81       53509 :     Real a = std::pow(seff, -1.0 / m) - 1.0;
      82             :     // Return 0 if pc > pc_max
      83       53509 :     if (std::pow(a, 1.0 - m) / alpha > pc_max)
      84             :       return 0.0;
      85             :     else
      86       49208 :       return (m - 1.0) * std::pow(a, -m) * std::pow(seff, -1.0 - 1.0 / m) / m / alpha;
      87             :   }
      88             : }
      89             : 
      90             : Real
      91       36403 : d2CapillaryPressure(Real seff, Real alpha, Real m, Real pc_max)
      92             : {
      93       36403 :   if (seff <= 0.0 || seff >= 1.0)
      94             :     return 0.0;
      95             :   else
      96             :   {
      97       31958 :     Real a = std::pow(seff, -1.0 / m) - 1.0;
      98             :     // Return 0 if pc > pc_max
      99       31958 :     if (std::pow(a, 1.0 - m) / alpha > pc_max)
     100             :       return 0.0;
     101             :     else
     102             :     {
     103       29853 :       Real d2pc = -std::pow(a, -1.0 - m) * std::pow(seff, -2.0 - 2.0 / m) +
     104       29853 :                   ((1.0 + m) / m) * std::pow(a, -m) * std::pow(seff, -1.0 / m - 2.0);
     105       29853 :       d2pc *= (1.0 - m) / m / alpha;
     106       29853 :       return d2pc;
     107             :     }
     108             :   }
     109             : }
     110             : 
     111             : Real
     112       96151 : relativePermeability(Real seff, Real m)
     113             : {
     114       96151 :   if (seff <= 0.0)
     115             :     return 0.0;
     116       96149 :   else if (seff >= 1.0)
     117             :     return 1.0;
     118             : 
     119       96120 :   const Real a = 1.0 - std::pow(seff, 1.0 / m);
     120       96120 :   const Real b = 1.0 - std::pow(a, m);
     121             : 
     122      192240 :   return std::sqrt(seff) * Utility::pow<2>(b);
     123             : }
     124             : 
     125             : Real
     126       96153 : dRelativePermeability(Real seff, Real m)
     127             : {
     128             :   // Guard against division by zero
     129       96153 :   if (seff <= 0.0 || seff >= 1.0)
     130             :     return 0.0;
     131             : 
     132       96122 :   const Real a = 1.0 - std::pow(seff, 1.0 / m);
     133       96122 :   const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
     134       96122 :   const Real b = 1.0 - std::pow(a, m);
     135       96122 :   const Real db = -m * std::pow(a, m - 1.0) * da;
     136             : 
     137      192244 :   return 0.5 * std::pow(seff, -0.5) * Utility::pow<2>(b) + 2.0 * std::sqrt(seff) * b * db;
     138             : }
     139             : 
     140             : Real
     141          14 : d2RelativePermeability(Real seff, Real m)
     142             : {
     143             :   // Guard against division by zero
     144          14 :   if (seff <= 0.0 || seff >= 1.0)
     145             :     return 0.0;
     146             : 
     147          10 :   const Real a = 1.0 - std::pow(seff, 1.0 / m);
     148          10 :   const Real da = -1.0 / m * std::pow(seff, 1.0 / m - 1.0);
     149          10 :   const Real d2a = -(1.0 / m) * (1.0 / m - 1.0) * std::pow(seff, 1.0 / m - 2.0);
     150          10 :   const Real b = 1.0 - std::pow(a, m);
     151          10 :   const Real db = -m * std::pow(a, m - 1.0) * da;
     152          10 :   const Real d2b = -m * (m - 1.0) * std::pow(a, m - 2.0) * da * da - m * std::pow(a, m - 1.0) * d2a;
     153             : 
     154          30 :   return -0.25 * std::pow(seff, -1.5) * Utility::pow<2>(b) + 2.0 * std::pow(seff, -0.5) * b * db +
     155          20 :          2.0 * std::sqrt(seff) * db * db + 2.0 * std::sqrt(seff) * b * d2b;
     156             : }
     157        2499 : }

Generated by: LCOV version 1.11