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 "PorousFlowPorosityExponentialBase.h"
      11             : 
      12             : InputParameters
      13       10181 : PorousFlowPorosityExponentialBase::validParams()
      14             : {
      15       10181 :   InputParameters params = PorousFlowPorosityBase::validParams();
      16       20362 :   params.addParam<bool>("strain_at_nearest_qp",
      17       20362 :                         false,
      18             :                         "When calculating nodal porosity that depends on strain, use the strain at "
      19             :                         "the nearest quadpoint.  This adds a small extra computational burden, and "
      20             :                         "is not necessary for simulations involving only linear lagrange elements. "
      21             :                         " If you set this to true, you will also want to set the same parameter to "
      22             :                         "true for related Kernels and Materials");
      23       20362 :   params.addParam<bool>("ensure_positive",
      24       20362 :                         true,
      25             :                         "Modify the usual exponential relationships that "
      26             :                         "governs porosity so that porosity is always "
      27             :                         "positive");
      28       10181 :   params.addClassDescription("Base class Material for porosity that is computed via an exponential "
      29             :                              "relationship with coupled variables (strain, porepressure, "
      30             :                              "temperature, chemistry)");
      31       10181 :   return params;
      32           0 : }
      33             : 
      34        8009 : PorousFlowPorosityExponentialBase::PorousFlowPorosityExponentialBase(
      35        8009 :     const InputParameters & parameters)
      36             :   : PorousFlowPorosityBase(parameters),
      37        8009 :     _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
      38       24027 :     _ensure_positive(getParam<bool>("ensure_positive"))
      39             : {
      40        8009 : }
      41             : 
      42             : void
      43      219092 : PorousFlowPorosityExponentialBase::initQpStatefulProperties()
      44             : {
      45      219092 :   const Real a = atNegInfinityQp();
      46      219092 :   const Real b = atZeroQp();
      47             :   mooseAssert(a > b, "PorousFlowPorosityExponentialBase a must be larger than b");
      48      219092 :   const Real decay = decayQp();
      49             : 
      50      219092 :   if (decay <= 0.0 || !_ensure_positive)
      51      216856 :     _porosity[_qp] = a + (b - a) * std::exp(decay);
      52             :   else
      53             :   {
      54        2236 :     const Real c = std::log(a / (a - b));
      55        2236 :     const Real expx = std::exp(-decay / c);
      56        2236 :     _porosity[_qp] = a + (b - a) * std::exp(c * (1.0 - expx));
      57             :   }
      58      219092 : }
      59             : 
      60             : void
      61     9369357 : PorousFlowPorosityExponentialBase::computeQpProperties()
      62             : {
      63     9369357 :   const Real a = atNegInfinityQp();
      64     9369357 :   const Real b = atZeroQp();
      65     9369357 :   const Real decay = decayQp();
      66             :   Real exp_term = 1.0; // set appropriately below
      67             : 
      68             :   Real deriv = 0.0; // = d(porosity)/d(decay)
      69     9369357 :   if (decay <= 0.0 || !_ensure_positive)
      70             :   {
      71     9322907 :     exp_term = std::exp(decay);
      72     9322907 :     _porosity[_qp] = a + (b - a) * exp_term;
      73     9322907 :     deriv = _porosity[_qp] - a;
      74             :   }
      75             :   else
      76             :   {
      77       46450 :     const Real c = std::log(a / (a - b));
      78       46450 :     const Real expx = std::exp(-decay / c);
      79             :     // note that at decay = 0, we have expx = 1, so porosity = a + b - a = b
      80             :     // and at decay = infinity, expx = 0, so porosity = a + (b - a) * a / (a - b) = 0
      81       46450 :     exp_term = std::exp(c * (1.0 - expx));
      82       46450 :     _porosity[_qp] = a + (b - a) * exp_term;
      83       46450 :     deriv = (_porosity[_qp] - a) * expx;
      84             :   }
      85             : 
      86     9369357 :   (*_dporosity_dvar)[_qp].resize(_num_var);
      87     9369357 :   (*_dporosity_dgradvar)[_qp].resize(_num_var);
      88    32161329 :   for (unsigned int v = 0; v < _num_var; ++v)
      89             :   {
      90    22791972 :     (*_dporosity_dvar)[_qp][v] = ddecayQp_dvar(v) * deriv;
      91    22791972 :     (*_dporosity_dgradvar)[_qp][v] = ddecayQp_dgradvar(v) * deriv;
      92             : 
      93    22791972 :     const Real da = datNegInfinityQp(v);
      94    22791972 :     const Real db = datZeroQp(v);
      95    22791972 :     (*_dporosity_dvar)[_qp][v] += da * (1 - exp_term) + db * exp_term;
      96             : 
      97    22791972 :     if (!(decay <= 0.0 || !_ensure_positive))
      98             :     {
      99      222198 :       const Real c = std::log(a / (a - b));
     100      222198 :       const Real expx = std::exp(-decay / c);
     101      222198 :       const Real dc = (a - b) * (da * b / a - db) / std::pow(a, 2);
     102      222198 :       (*_dporosity_dvar)[_qp][v] += (b - a) * exp_term * dc * (1 - expx - expx / c);
     103             :     }
     104             :   }
     105     9369357 : }