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         240 : PorousFlow1PhaseMD_Gaussian::validParams()
      16             : {
      17         240 :   InputParameters params = PorousFlowVariableBase::validParams();
      18         480 :   params.addRequiredCoupledVar("mass_density",
      19             :                                "Variable that represents log(mass-density) of the single phase");
      20         480 :   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         480 :   params.addRequiredRangeCheckedParam<Real>(
      26             :       "density_P0", "density_P0>0", "The density of the fluid phase at zero porepressure");
      27         480 :   params.addRequiredRangeCheckedParam<Real>(
      28             :       "bulk_modulus", "bulk_modulus>0", "The constant bulk modulus of the fluid phase");
      29         240 :   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         240 :   return params;
      35           0 : }
      36             : 
      37         186 : PorousFlow1PhaseMD_Gaussian::PorousFlow1PhaseMD_Gaussian(const InputParameters & parameters)
      38             :   : PorousFlowVariableBase(parameters),
      39             : 
      40         186 :     _al(getParam<Real>("al")),
      41         186 :     _al2(std::pow(_al, 2)),
      42         372 :     _logdens0(std::log(getParam<Real>("density_P0"))),
      43         372 :     _bulk(getParam<Real>("bulk_modulus")),
      44         186 :     _recip_bulk(1.0 / _al / _bulk),
      45         186 :     _recip_bulk2(std::pow(_recip_bulk, 2)),
      46             : 
      47         258 :     _md_var(_nodal_material ? coupledDofValues("mass_density") : coupledValue("mass_density")),
      48         186 :     _gradmd_qp_var(coupledGradient("mass_density")),
      49         186 :     _md_varnum(coupled("mass_density")),
      50         186 :     _pvar(_dictator.isPorousFlowVariable(_md_varnum) ? _dictator.porousFlowVariableNum(_md_varnum)
      51         186 :                                                      : 0)
      52             : {
      53         186 :   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         186 : }
      59             : 
      60             : void
      61         584 : PorousFlow1PhaseMD_Gaussian::initQpStatefulProperties()
      62             : {
      63         584 :   PorousFlowVariableBase::initQpStatefulProperties();
      64         584 :   buildPS();
      65         584 : }
      66             : 
      67             : void
      68       13896 : PorousFlow1PhaseMD_Gaussian::computeQpProperties()
      69             : {
      70             :   // size stuff correctly and prepare the derivative matrices with zeroes
      71       13896 :   PorousFlowVariableBase::computeQpProperties();
      72             : 
      73       13896 :   buildPS();
      74             : 
      75       13896 :   if (!_nodal_material)
      76             :   {
      77        6804 :     if (_md_var[_qp] >= _logdens0)
      78             :     {
      79        6444 :       (*_gradp_qp)[_qp][0] = _gradmd_qp_var[_qp] * _bulk;
      80        6444 :       (*_grads_qp)[_qp][0] = 0.0;
      81             :     }
      82             :     else
      83             :     {
      84         360 :       (*_gradp_qp)[_qp][0] =
      85         360 :           _gradmd_qp_var[_qp] / (_recip_bulk - 2.0 * _al * _porepressure[_qp][0]) / _al;
      86         360 :       (*_grads_qp)[_qp][0] =
      87         360 :           -2.0 * _al2 * _porepressure[_qp][0] * _saturation[_qp][0] * (*_gradp_qp)[_qp][0];
      88             :     }
      89             :   }
      90             : 
      91       13896 :   if (_dictator.notPorousFlowVariable(_md_varnum))
      92             :     return;
      93             : 
      94       13896 :   if (_md_var[_qp] >= _logdens0)
      95             :   {
      96             :     // fully saturated at the node or quadpoint
      97       13030 :     (*_dporepressure_dvar)[_qp][0][_pvar] = _bulk;
      98       13030 :     (*_dsaturation_dvar)[_qp][0][_pvar] = 0.0;
      99             :   }
     100             :   else
     101             :   {
     102         866 :     const Real pp = _porepressure[_qp][0];
     103         866 :     (*_dporepressure_dvar)[_qp][0][_pvar] = 1.0 / (_recip_bulk - 2.0 * _al * pp) / _al;
     104         866 :     const Real sat = _saturation[_qp][0];
     105         866 :     (*_dsaturation_dvar)[_qp][0][_pvar] =
     106         866 :         -2.0 * _al2 * pp * sat * (*_dporepressure_dvar)[_qp][0][_pvar];
     107             :   }
     108             : 
     109       13896 :   if (!_nodal_material)
     110             :   {
     111        6804 :     if (_md_var[_qp] >= _logdens0)
     112             :     {
     113             :       // fully saturated at the quadpoint
     114        6444 :       (*_dgradp_qp_dgradv)[_qp][0][_pvar] = _bulk;
     115        6444 :       (*_dgradp_qp_dv)[_qp][0][_pvar] = 0.0;
     116        6444 :       (*_dgrads_qp_dgradv)[_qp][0][_pvar] = 0.0;
     117        6444 :       (*_dgrads_qp_dv)[_qp][0][_pvar] = 0.0;
     118             :     }
     119             :     else
     120             :     {
     121         360 :       const Real pp = _porepressure[_qp][0];
     122         360 :       (*_dgradp_qp_dgradv)[_qp][0][_pvar] = 1.0 / (_recip_bulk - 2.0 * _al * pp) / _al;
     123         360 :       (*_dgradp_qp_dv)[_qp][0][_pvar] = _gradmd_qp_var[_qp] * 2.0 * _al *
     124         360 :                                         (*_dporepressure_dvar)[_qp][0][_pvar] /
     125         360 :                                         std::pow(_recip_bulk - 2.0 * _al * pp, 2.0) / _al;
     126         360 :       const Real sat = _saturation[_qp][0];
     127         360 :       (*_dgrads_qp_dgradv)[_qp][0][_pvar] =
     128         360 :           -2.0 * _al2 * pp * sat * (*_dgradp_qp_dgradv)[_qp][0][_pvar];
     129         360 :       (*_dgrads_qp_dv)[_qp][0][_pvar] =
     130         360 :           -2.0 * _al2 * (*_dporepressure_dvar)[_qp][0][_pvar] * sat * (*_gradp_qp)[_qp][0];
     131             :       (*_dgrads_qp_dv)[_qp][0][_pvar] +=
     132         360 :           -2.0 * _al2 * pp * (*_dsaturation_dvar)[_qp][0][_pvar] * (*_gradp_qp)[_qp][0];
     133         360 :       (*_dgrads_qp_dv)[_qp][0][_pvar] += -2.0 * _al2 * pp * sat * (*_dgradp_qp_dv)[_qp][0][_pvar];
     134             :     }
     135             :   }
     136             : }
     137             : 
     138             : void
     139       14480 : PorousFlow1PhaseMD_Gaussian::buildPS()
     140             : {
     141       14480 :   if (_md_var[_qp] >= _logdens0)
     142             :   {
     143             :     // full saturation
     144       13582 :     _porepressure[_qp][0] = (_md_var[_qp] - _logdens0) * _bulk;
     145       13582 :     _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         898 :     _porepressure[_qp][0] =
     154         898 :         (_recip_bulk - std::sqrt(_recip_bulk2 + 4.0 * (_logdens0 - _md_var[_qp]))) / (2.0 * _al);
     155         898 :     _saturation[_qp][0] = std::exp(-std::pow(_al * _porepressure[_qp][0], 2.0));
     156             :   }
     157       14480 : }