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             : #include <limits>
      13             : 
      14             : InputParameters
      15        5198 : PorousFlowPorosityExponentialBase::validParams()
      16             : {
      17        5198 :   InputParameters params = PorousFlowPorosityBase::validParams();
      18       10396 :   params.addParam<bool>("strain_at_nearest_qp",
      19       10396 :                         false,
      20             :                         "When calculating nodal porosity that depends on strain, use the strain at "
      21             :                         "the nearest quadpoint.  This adds a small extra computational burden, and "
      22             :                         "is not necessary for simulations involving only linear lagrange elements. "
      23             :                         " If you set this to true, you will also want to set the same parameter to "
      24             :                         "true for related Kernels and Materials");
      25       10396 :   params.addParam<bool>("ensure_positive",
      26       10396 :                         true,
      27             :                         "Modify the usual exponential relationships that "
      28             :                         "governs porosity so that porosity is always "
      29             :                         "positive");
      30       10396 :   params.addParam<Real>(
      31             :       "porosity_min",
      32       10396 :       std::numeric_limits<Real>::lowest(),
      33             :       "Minimum allowed value of the porosity: if the computed porosity is less than this value, "
      34             :       "porosity is set to this value instead.  By default no floor is imposed.  The "
      35             :       "ensure_positive "
      36             :       "transform only acts for decay > 0, so chemistry-driven (precipitation) porosity is "
      37             :       "otherwise "
      38             :       "unbounded below and can become negative once pore space is filled by mineral; set this to a "
      39             :       "small positive value in that case.");
      40       10396 :   params.addParam<Real>("zero_modifier",
      41       10396 :                         1E-3,
      42             :                         "If the porosity_min floor is active, the porosity derivatives are set to "
      43             :                         "zero_modifier times their unfloored values (rather than exactly zero) to "
      44             :                         "hint to the Newton-Krylov nonlinear solver that porosity "
      45             :                         "is not strictly constant, which aids convergence");
      46       10396 :   params.addParamNamesToGroup("zero_modifier", "Advanced");
      47        5198 :   params.addClassDescription("Base class Material for porosity that is computed via an exponential "
      48             :                              "relationship with coupled variables (strain, porepressure, "
      49             :                              "temperature, chemistry)");
      50        5198 :   return params;
      51           0 : }
      52             : 
      53        4037 : PorousFlowPorosityExponentialBase::PorousFlowPorosityExponentialBase(
      54        4037 :     const InputParameters & parameters)
      55             :   : PorousFlowPorosityBase(parameters),
      56        4037 :     _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
      57        8074 :     _ensure_positive(getParam<bool>("ensure_positive")),
      58        8074 :     _porosity_min(getParam<Real>("porosity_min")),
      59       12111 :     _zero_modifier(getParam<Real>("zero_modifier"))
      60             : {
      61        4037 : }
      62             : 
      63             : void
      64      154022 : PorousFlowPorosityExponentialBase::initQpStatefulProperties()
      65             : {
      66      154022 :   const Real a = atNegInfinityQp();
      67      154022 :   const Real b = atZeroQp();
      68             :   mooseAssert(a > b, "PorousFlowPorosityExponentialBase a must be larger than b");
      69      154022 :   const Real decay = decayQp();
      70             : 
      71      154022 :   if (decay <= 0.0 || !_ensure_positive)
      72      152770 :     _porosity[_qp] = a + (b - a) * std::exp(decay);
      73             :   else
      74             :   {
      75        1252 :     const Real c = std::log(a / (a - b));
      76        1252 :     const Real expx = std::exp(-decay / c);
      77        1252 :     _porosity[_qp] = a + (b - a) * std::exp(c * (1.0 - expx));
      78             :   }
      79             : 
      80      154022 :   if (_porosity[_qp] < _porosity_min)
      81           0 :     _porosity[_qp] = _porosity_min;
      82      154022 : }
      83             : 
      84             : void
      85     7053656 : PorousFlowPorosityExponentialBase::computeQpProperties()
      86             : {
      87     7053656 :   const Real a = atNegInfinityQp();
      88     7053656 :   const Real b = atZeroQp();
      89     7053656 :   const Real decay = decayQp();
      90             :   Real exp_term = 1.0; // set appropriately below
      91             : 
      92             :   Real deriv = 0.0; // = d(porosity)/d(decay)
      93     7053656 :   if (decay <= 0.0 || !_ensure_positive)
      94             :   {
      95     7016535 :     exp_term = std::exp(decay);
      96     7016535 :     _porosity[_qp] = a + (b - a) * exp_term;
      97     7016535 :     deriv = _porosity[_qp] - a;
      98             :   }
      99             :   else
     100             :   {
     101       37121 :     const Real c = std::log(a / (a - b));
     102       37121 :     const Real expx = std::exp(-decay / c);
     103             :     // note that at decay = 0, we have expx = 1, so porosity = a + b - a = b
     104             :     // and at decay = infinity, expx = 0, so porosity = a + (b - a) * a / (a - b) = 0
     105       37121 :     exp_term = std::exp(c * (1.0 - expx));
     106       37121 :     _porosity[_qp] = a + (b - a) * exp_term;
     107       37121 :     deriv = (_porosity[_qp] - a) * expx;
     108             :   }
     109             : 
     110     7053656 :   (*_dporosity_dvar)[_qp].resize(_num_var);
     111     7053656 :   (*_dporosity_dgradvar)[_qp].resize(_num_var);
     112    23099595 :   for (unsigned int v = 0; v < _num_var; ++v)
     113             :   {
     114    16045939 :     (*_dporosity_dvar)[_qp][v] = ddecayQp_dvar(v) * deriv;
     115    16045939 :     (*_dporosity_dgradvar)[_qp][v] = ddecayQp_dgradvar(v) * deriv;
     116             : 
     117    16045939 :     const Real da = datNegInfinityQp(v);
     118    16045939 :     const Real db = datZeroQp(v);
     119    16045939 :     (*_dporosity_dvar)[_qp][v] += da * (1 - exp_term) + db * exp_term;
     120             : 
     121    16045939 :     if (!(decay <= 0.0 || !_ensure_positive))
     122             :     {
     123      177219 :       const Real c = std::log(a / (a - b));
     124      177219 :       const Real expx = std::exp(-decay / c);
     125      177219 :       const Real dc = (a - b) * (da * b / a - db) / std::pow(a, 2);
     126      177219 :       (*_dporosity_dvar)[_qp][v] += (b - a) * exp_term * dc * (1 - expx - expx / c);
     127             :     }
     128             :   }
     129             : 
     130             :   // Apply the porosity floor last, after the unfloored derivatives above have been
     131             :   // formed.  When floored, soften the derivatives with _zero_modifier so the Newton
     132             :   // process still sees porosity as weakly varying (cf. PorousFlowPorosityLinear).
     133     7053656 :   if (_porosity[_qp] < _porosity_min)
     134             :   {
     135         114 :     _porosity[_qp] = _porosity_min;
     136         228 :     for (unsigned int v = 0; v < _num_var; ++v)
     137             :     {
     138         114 :       (*_dporosity_dvar)[_qp][v] *= _zero_modifier;
     139         114 :       (*_dporosity_dgradvar)[_qp][v] *= _zero_modifier;
     140             :     }
     141             :   }
     142     7053656 : }