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             : //  "cut" van-Genuchten effective saturation as a function of pressure, and its derivs wrt p
      11             : //
      12             : #include "RichardsSeff1VGcut.h"
      13             : 
      14             : registerMooseObject("RichardsApp", RichardsSeff1VGcut);
      15             : 
      16             : InputParameters
      17          13 : RichardsSeff1VGcut::validParams()
      18             : {
      19          13 :   InputParameters params = RichardsSeff1VG::validParams();
      20          26 :   params.addRequiredRangeCheckedParam<Real>(
      21             :       "p_cut",
      22             :       "p_cut < 0",
      23             :       "cutoff in pressure.  Must be negative.  If p>p_cut then use "
      24             :       "van-Genuchten function.  Otherwise use a linear relationship which is "
      25             :       "chosen so the value and derivative match van-Genuchten at p=p_cut");
      26          13 :   params.addClassDescription("cut van-Genuchten effective saturation as a function of capillary "
      27             :                              "pressure.  Single-phase  seff = (1 + (-al*p)^(1/(1-m)))^(-m) for "
      28             :                              "p>p_cut, otherwise user a a linear relationship that is chosen so "
      29             :                              "the value and derivative match van-Genuchten at p=p_cut.");
      30          13 :   return params;
      31           0 : }
      32             : 
      33           6 : RichardsSeff1VGcut::RichardsSeff1VGcut(const InputParameters & parameters)
      34             :   : RichardsSeff1VG(parameters),
      35           6 :     _al(getParam<Real>("al")),
      36          12 :     _m(getParam<Real>("m")),
      37          12 :     _p_cut(getParam<Real>("p_cut")),
      38           6 :     _s_cut(0),
      39           6 :     _ds_cut(0)
      40             : {
      41           6 :   _s_cut = RichardsSeffVG::seff(_p_cut, _al, _m);
      42           6 :   _ds_cut = RichardsSeffVG::dseff(_p_cut, _al, _m);
      43           6 : }
      44             : 
      45             : void
      46           2 : RichardsSeff1VGcut::initialSetup()
      47             : {
      48           2 :   _console << "cut VG Seff has p_cut=" << _p_cut << " so seff_cut=" << _s_cut
      49           2 :            << " and seff=0 at p=" << -_s_cut / _ds_cut + _p_cut << std::endl;
      50           2 : }
      51             : 
      52             : Real
      53         606 : RichardsSeff1VGcut::seff(std::vector<const VariableValue *> p, unsigned int qp) const
      54             : {
      55         606 :   if ((*p[0])[qp] > _p_cut)
      56             :   {
      57         360 :     return RichardsSeff1VG::seff(p, qp);
      58             :   }
      59             :   else
      60             :   {
      61         246 :     Real seff_linear = _s_cut + _ds_cut * ((*p[0])[qp] - _p_cut);
      62             :     // return (seff_linear > 0 ? seff_linear : 0); // andy isn't sure of this - might be useful to
      63             :     // allow negative saturations
      64         246 :     return seff_linear;
      65             :   }
      66             : }
      67             : 
      68             : void
      69         606 : RichardsSeff1VGcut::dseff(std::vector<const VariableValue *> p,
      70             :                           unsigned int qp,
      71             :                           std::vector<Real> & result) const
      72             : {
      73         606 :   if ((*p[0])[qp] > _p_cut)
      74         360 :     return RichardsSeff1VG::dseff(p, qp, result);
      75             :   else
      76         246 :     result[0] = _ds_cut;
      77             : }
      78             : 
      79             : void
      80         606 : RichardsSeff1VGcut::d2seff(std::vector<const VariableValue *> p,
      81             :                            unsigned int qp,
      82             :                            std::vector<std::vector<Real>> & result) const
      83             : {
      84         606 :   if ((*p[0])[qp] > _p_cut)
      85         360 :     return RichardsSeff1VG::d2seff(p, qp, result);
      86             :   else
      87         246 :     result[0][0] = 0;
      88             : }