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