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 "PorousFlowRelativePermeabilityVG.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13             : registerMooseObject("PorousFlowApp", PorousFlowRelativePermeabilityVG);
      14             : registerMooseObject("PorousFlowApp", ADPorousFlowRelativePermeabilityVG);
      15             : 
      16             : template <bool is_ad>
      17             : InputParameters
      18        1021 : PorousFlowRelativePermeabilityVGTempl<is_ad>::validParams()
      19             : {
      20        1021 :   InputParameters params = PorousFlowRelativePermeabilityBaseTempl<is_ad>::validParams();
      21        2042 :   params.addRequiredRangeCheckedParam<Real>(
      22             :       "m", "m > 0 & m < 1", "The van Genuchten exponent of the phase");
      23        3063 :   params.addRangeCheckedParam<Real>("seff_turnover",
      24        2042 :                                     1.0,
      25             :                                     "seff_turnover > 0 & seff_turnover <= 1",
      26             :                                     "The relative permeability will be a cubic for seff > "
      27             :                                     "seff_turnover.  The cubic is chosen so that its derivative "
      28             :                                     "and value matche the VG function at seff=seff_turnover");
      29        2042 :   params.addParam<bool>("zero_derivative",
      30        2042 :                         false,
      31             :                         "Employ a cubic for seff>seff_turnover that has zero derivative at seff=1");
      32        2042 :   params.addParam<bool>("wetting",
      33        2042 :                         true,
      34             :                         "If true, use the van Genuchten form appropriate for a wetting (liquid) "
      35             :                         "phase.  If false, use the non-wetting (gas) expression.");
      36        1021 :   params.addClassDescription("This Material calculates relative permeability of a phase "
      37             :                              "using the van Genuchten model");
      38        1021 :   return params;
      39           0 : }
      40             : 
      41             : template <bool is_ad>
      42         795 : PorousFlowRelativePermeabilityVGTempl<is_ad>::PorousFlowRelativePermeabilityVGTempl(
      43             :     const InputParameters & parameters)
      44             :   : PorousFlowRelativePermeabilityBaseTempl<is_ad>(parameters),
      45         795 :     _m(this->template getParam<Real>("m")),
      46        1590 :     _wetting(this->template getParam<bool>("wetting")),
      47        1590 :     _cut(this->template getParam<Real>("seff_turnover")),
      48         795 :     _cub0(_wetting ? PorousFlowVanGenuchten::relativePermeability(_cut, _m)
      49         132 :                    : PorousFlowVanGenuchten::relativePermeabilityNW(_cut, _m)),
      50         795 :     _cub1(_wetting ? PorousFlowVanGenuchten::dRelativePermeability(_cut, _m)
      51         132 :                    : PorousFlowVanGenuchten::dRelativePermeabilityNW(_cut, _m)),
      52        1590 :     _cub2(_cut < 1.0
      53         894 :               ? (this->template getParam<bool>("zero_derivative")
      54          99 :                      ? 3.0 * (1.0 - _cub0 - _cub1 * (1.0 - _cut)) / Utility::pow<2>(1.0 - _cut) +
      55           0 :                            _cub1 / (1.0 - _cut)
      56          99 :                      : (_wetting ? PorousFlowVanGenuchten::d2RelativePermeability(_cut, _m)
      57           0 :                                  : PorousFlowVanGenuchten::d2RelativePermeabilityNW(_cut, _m)))
      58             :               : 0.0),
      59        1590 :     _cub3(_cut < 1.0
      60         894 :               ? (this->template getParam<bool>("zero_derivative")
      61          99 :                      ? -2.0 * (1.0 - _cub0 - _cub1 * (1.0 - _cut)) / Utility::pow<3>(1.0 - _cut) -
      62           0 :                            _cub1 / Utility::pow<2>(1.0 - _cut)
      63          99 :                      : (1.0 - _cub0 - _cub1 * (1.0 - _cut) - _cub2 * Utility::pow<2>(1.0 - _cut)) /
      64             :                            Utility::pow<3>(1.0 - _cut))
      65         795 :               : 0.0)
      66             : {
      67         795 : }
      68             : 
      69             : template <bool is_ad>
      70             : GenericReal<is_ad>
      71      767720 : PorousFlowRelativePermeabilityVGTempl<is_ad>::relativePermeability(GenericReal<is_ad> seff) const
      72             : {
      73             :   if (MetaPhysicL::raw_value(seff) < _cut)
      74             :   {
      75      759523 :     if (_wetting)
      76      750343 :       return PorousFlowVanGenuchten::relativePermeability(seff, _m);
      77             :     else
      78        9180 :       return PorousFlowVanGenuchten::relativePermeabilityNW(seff, _m);
      79             :   }
      80             : 
      81        8197 :   return _cub0 + _cub1 * (seff - _cut) + _cub2 * Utility::pow<2>(seff - _cut) +
      82        8197 :          _cub3 * Utility::pow<3>(seff - _cut);
      83             : }
      84             : 
      85             : template <bool is_ad>
      86             : Real
      87      767720 : PorousFlowRelativePermeabilityVGTempl<is_ad>::dRelativePermeability(Real seff) const
      88             : {
      89      767720 :   if (seff < _cut)
      90             :   {
      91      759523 :     if (_wetting)
      92      750343 :       return PorousFlowVanGenuchten::dRelativePermeability(seff, _m);
      93             :     else
      94        9180 :       return PorousFlowVanGenuchten::dRelativePermeabilityNW(seff, _m);
      95             :   }
      96             : 
      97        8197 :   return _cub1 + 2.0 * _cub2 * (seff - _cut) + 3.0 * _cub3 * Utility::pow<2>(seff - _cut);
      98             : }
      99             : 
     100             : template class PorousFlowRelativePermeabilityVGTempl<false>;
     101             : template class PorousFlowRelativePermeabilityVGTempl<true>;