LCOV - code coverage report
Current view: top level - src/materials - PorousFlow1PhaseMD_Gaussian.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 75 77 97.4 %
Date: 2017-11-18 13:30:36 Functions: 7 7 100.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11