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             : //  van-Genuchten effective saturation as a function of pressure (not capillary pressure)
      11             : //
      12             : #include "RichardsSeffVG.h"
      13             : 
      14             : Real
      15    25919422 : RichardsSeffVG::seff(Real p, Real al, Real m)
      16             : {
      17             :   Real n, seff;
      18             : 
      19    25919422 :   if (p >= 0)
      20             :     return 1.0;
      21             :   else
      22             :   {
      23     8442136 :     n = 1.0 / (1.0 - m);
      24     8442136 :     seff = 1 + std::pow(-al * p, n);
      25     8442136 :     return std::pow(seff, -m);
      26             :   }
      27             : }
      28             : 
      29             : Real
      30     9257704 : RichardsSeffVG::dseff(Real p, Real al, Real m)
      31             : {
      32     9257704 :   if (p >= 0)
      33             :     return 0.0;
      34             :   else
      35             :   {
      36     3131403 :     Real n = 1.0 / (1.0 - m);
      37     3131403 :     Real inner = 1 + std::pow(-al * p, n);
      38     3131403 :     Real dinner_dp = -n * al * std::pow(-al * p, n - 1);
      39     3131403 :     Real dseff_dp = -m * std::pow(inner, -m - 1) * dinner_dp;
      40     3131403 :     return dseff_dp;
      41             :   }
      42             : }
      43             : 
      44             : Real
      45     5126892 : RichardsSeffVG::d2seff(Real p, Real al, Real m)
      46             : {
      47     5126892 :   if (p >= 0)
      48             :     return 0.0;
      49             :   else
      50             :   {
      51     1957223 :     Real n = 1.0 / (1.0 - m);
      52     1957223 :     Real inner = 1 + std::pow(-al * p, n);
      53     1957223 :     Real dinner_dp = -n * al * std::pow(-al * p, n - 1);
      54     1957223 :     Real d2inner_dp2 = n * (n - 1) * al * al * std::pow(-al * p, n - 2);
      55     1957223 :     Real d2seff_dp2 = m * (m + 1) * std::pow(inner, -m - 2) * std::pow(dinner_dp, 2) -
      56     1957223 :                       m * std::pow(inner, -m - 1) * d2inner_dp2;
      57     1957223 :     return d2seff_dp2;
      58             :   }
      59             : }