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