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 "PorousFlow1PhaseMD_Gaussian.h"
      11             : 
      12             : registerMooseObject("PorousFlowApp", PorousFlow1PhaseMD_Gaussian);
      13             : 
      14             : InputParameters
      15         416 : PorousFlow1PhaseMD_Gaussian::validParams()
      16             : {
      17         416 :   InputParameters params = PorousFlowVariableBase::validParams();
      18         832 :   params.addRequiredCoupledVar("mass_density",
      19             :                                "Variable that represents log(mass-density) of the single phase");
      20         832 :   params.addRequiredRangeCheckedParam<Real>(
      21             :       "al",
      22             :       "al>0",
      23             :       "For this class, the capillary function is assumed to be saturation = "
      24             :       "exp(-(al*porepressure)^2) for porepressure<0.");
      25         832 :   params.addRequiredRangeCheckedParam<Real>(
      26             :       "density_P0", "density_P0>0", "The density of the fluid phase at zero porepressure");
      27         832 :   params.addRequiredRangeCheckedParam<Real>(
      28             :       "bulk_modulus", "bulk_modulus>0", "The constant bulk modulus of the fluid phase");
      29         416 :   params.addClassDescription("This Material is used for the single-phase situation where "
      30             :                              "log(mass-density) is the primary variable.  calculates the 1 "
      31             :                              "porepressure and the 1 saturation in a 1-phase situation, "
      32             :                              "and derivatives of these with respect to the PorousFlowVariables.  A "
      33             :                              "gaussian capillary function is assumed");
      34         416 :   return params;
      35           0 : }
      36             : 
      37         330 : PorousFlow1PhaseMD_Gaussian::PorousFlow1PhaseMD_Gaussian(const InputParameters & parameters)
      38             :   : PorousFlowVariableBase(parameters),
      39             : 
      40         330 :     _al(getParam<Real>("al")),
      41         330 :     _al2(std::pow(_al, 2)),
      42         660 :     _logdens0(std::log(getParam<Real>("density_P0"))),
      43         660 :     _bulk(getParam<Real>("bulk_modulus")),
      44         330 :     _recip_bulk(1.0 / _al / _bulk),
      45         330 :     _recip_bulk2(std::pow(_recip_bulk, 2)),
      46             : 
      47         462 :     _md_var(_nodal_material ? coupledDofValues("mass_density") : coupledValue("mass_density")),
      48         330 :     _gradmd_qp_var(coupledGradient("mass_density")),
      49         330 :     _md_varnum(coupled("mass_density")),
      50         330 :     _pvar(_dictator.isPorousFlowVariable(_md_varnum) ? _dictator.porousFlowVariableNum(_md_varnum)
      51         330 :                                                      : 0)
      52             : {
      53         330 :   if (_num_phases != 1)
      54           0 :     mooseError("The Dictator proclaims that the number of phases is ",
      55           0 :                _dictator.numPhases(),
      56             :                " whereas PorousFlow1PhaseMD_Gaussian can only be used for 1-phase simulations.  Be "
      57             :                "aware that the Dictator has noted your mistake.");
      58         330 : }
      59             : 
      60             : void
      61         920 : PorousFlow1PhaseMD_Gaussian::initQpStatefulProperties()
      62             : {
      63         920 :   PorousFlowVariableBase::initQpStatefulProperties();
      64         920 :   buildPS();
      65         920 : }
      66             : 
      67             : void
      68       20320 : PorousFlow1PhaseMD_Gaussian::computeQpProperties()
      69             : {
      70             :   // size stuff correctly and prepare the derivative matrices with zeroes
      71       20320 :   PorousFlowVariableBase::computeQpProperties();
      72             : 
      73       20320 :   buildPS();
      74             : 
      75       20320 :   if (!_nodal_material)
      76             :   {
      77        9980 :     if (_md_var[_qp] >= _logdens0)
      78             :     {
      79        9530 :       (*_gradp_qp)[_qp][0] = _gradmd_qp_var[_qp] * _bulk;
      80        9530 :       (*_grads_qp)[_qp][0] = 0.0;
      81             :     }
      82             :     else
      83             :     {
      84         450 :       (*_gradp_qp)[_qp][0] =
      85         450 :           _gradmd_qp_var[_qp] / (_recip_bulk - 2.0 * _al * _porepressure[_qp][0]) / _al;
      86         450 :       (*_grads_qp)[_qp][0] =
      87         450 :           -2.0 * _al2 * _porepressure[_qp][0] * _saturation[_qp][0] * (*_gradp_qp)[_qp][0];
      88             :     }
      89             :   }
      90             : 
      91       20320 :   if (_dictator.notPorousFlowVariable(_md_varnum))
      92             :     return;
      93             : 
      94       20320 :   if (_md_var[_qp] >= _logdens0)
      95             :   {
      96             :     // fully saturated at the node or quadpoint
      97       19236 :     (*_dporepressure_dvar)[_qp][0][_pvar] = _bulk;
      98       19236 :     (*_dsaturation_dvar)[_qp][0][_pvar] = 0.0;
      99             :   }
     100             :   else
     101             :   {
     102        1084 :     const Real pp = _porepressure[_qp][0];
     103        1084 :     (*_dporepressure_dvar)[_qp][0][_pvar] = 1.0 / (_recip_bulk - 2.0 * _al * pp) / _al;
     104        1084 :     const Real sat = _saturation[_qp][0];
     105        1084 :     (*_dsaturation_dvar)[_qp][0][_pvar] =
     106        1084 :         -2.0 * _al2 * pp * sat * (*_dporepressure_dvar)[_qp][0][_pvar];
     107             :   }
     108             : 
     109       20320 :   if (!_nodal_material)
     110             :   {
     111        9980 :     if (_md_var[_qp] >= _logdens0)
     112             :     {
     113             :       // fully saturated at the quadpoint
     114        9530 :       (*_dgradp_qp_dgradv)[_qp][0][_pvar] = _bulk;
     115        9530 :       (*_dgradp_qp_dv)[_qp][0][_pvar] = 0.0;
     116        9530 :       (*_dgrads_qp_dgradv)[_qp][0][_pvar] = 0.0;
     117        9530 :       (*_dgrads_qp_dv)[_qp][0][_pvar] = 0.0;
     118             :     }
     119             :     else
     120             :     {
     121         450 :       const Real pp = _porepressure[_qp][0];
     122         450 :       (*_dgradp_qp_dgradv)[_qp][0][_pvar] = 1.0 / (_recip_bulk - 2.0 * _al * pp) / _al;
     123         450 :       (*_dgradp_qp_dv)[_qp][0][_pvar] = _gradmd_qp_var[_qp] * 2.0 * _al *
     124         450 :                                         (*_dporepressure_dvar)[_qp][0][_pvar] /
     125         450 :                                         std::pow(_recip_bulk - 2.0 * _al * pp, 2.0) / _al;
     126         450 :       const Real sat = _saturation[_qp][0];
     127         450 :       (*_dgrads_qp_dgradv)[_qp][0][_pvar] =
     128         450 :           -2.0 * _al2 * pp * sat * (*_dgradp_qp_dgradv)[_qp][0][_pvar];
     129         450 :       (*_dgrads_qp_dv)[_qp][0][_pvar] =
     130         450 :           -2.0 * _al2 * (*_dporepressure_dvar)[_qp][0][_pvar] * sat * (*_gradp_qp)[_qp][0];
     131             :       (*_dgrads_qp_dv)[_qp][0][_pvar] +=
     132         450 :           -2.0 * _al2 * pp * (*_dsaturation_dvar)[_qp][0][_pvar] * (*_gradp_qp)[_qp][0];
     133         450 :       (*_dgrads_qp_dv)[_qp][0][_pvar] += -2.0 * _al2 * pp * sat * (*_dgradp_qp_dv)[_qp][0][_pvar];
     134             :     }
     135             :   }
     136             : }
     137             : 
     138             : void
     139       21240 : PorousFlow1PhaseMD_Gaussian::buildPS()
     140             : {
     141       21240 :   if (_md_var[_qp] >= _logdens0)
     142             :   {
     143             :     // full saturation
     144       20116 :     _porepressure[_qp][0] = (_md_var[_qp] - _logdens0) * _bulk;
     145       20116 :     _saturation[_qp][0] = 1.0;
     146             :   }
     147             :   else
     148             :   {
     149             :     // v = logdens0 + p/bulk - (al p)^2
     150             :     // 0 = (v-logdens0) - p/bulk + (al p)^2
     151             :     // 2 al p = (1/al/bulk) +/- sqrt((1/al/bulk)^2 - 4(v-logdens0))  (the "minus" sign is chosen)
     152             :     // s = exp(-(al*p)^2)
     153        1124 :     _porepressure[_qp][0] =
     154        1124 :         (_recip_bulk - std::sqrt(_recip_bulk2 + 4.0 * (_logdens0 - _md_var[_qp]))) / (2.0 * _al);
     155        1124 :     _saturation[_qp][0] = std::exp(-std::pow(_al * _porepressure[_qp][0], 2.0));
     156             :   }
     157       21240 : }