www.mooseframework.org
PorousFlowMassRadioactiveDecay.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 // MOOSE includes
11 #include "MooseVariable.h"
12 
13 #include "libmesh/quadrature.h"
14 
15 // C++ includes
16 #include <limits>
17 
18 template <>
19 InputParameters
21 {
22  InputParameters params = validParams<TimeKernel>();
23  params.addParam<bool>("strain_at_nearest_qp",
24  false,
25  "When calculating nodal porosity that depends on strain, use the strain at "
26  "the nearest quadpoint. This adds a small extra computational burden, and "
27  "is not necessary for simulations involving only linear lagrange elements. "
28  " If you set this to true, you will also want to set the same parameter to "
29  "true for related Kernels and Materials");
30  params.addParam<unsigned int>(
31  "fluid_component", 0, "The index corresponding to the fluid component for this kernel");
32  params.addRequiredParam<Real>("decay_rate",
33  "The decay rate (units 1/time) for the fluid component");
34  params.addRequiredParam<UserObjectName>(
35  "PorousFlowDictator", "The UserObject that holds the list of Porous-Flow variable names.");
36  params.addClassDescription("Radioactive decay of a fluid component");
37  return params;
38 }
39 
41  : TimeKernel(parameters),
42  _decay_rate(getParam<Real>("decay_rate")),
43  _fluid_component(getParam<unsigned int>("fluid_component")),
44  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
45  _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
46  _num_phases(_dictator.numPhases()),
47  _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
48  _porosity(getMaterialProperty<Real>("PorousFlow_porosity_nodal")),
49  _dporosity_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_nodal_dvar")),
50  _dporosity_dgradvar(
51  getMaterialProperty<std::vector<RealGradient>>("dPorousFlow_porosity_nodal_dgradvar")),
52  _nearest_qp(_strain_at_nearest_qp
53  ? &getMaterialProperty<unsigned int>("PorousFlow_nearestqp_nodal")
54  : nullptr),
55  _fluid_density(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
56  _dfluid_density_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
57  "dPorousFlow_fluid_phase_density_nodal_dvar")),
58  _fluid_saturation_nodal(getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")),
59  _dfluid_saturation_nodal_dvar(
60  getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_saturation_nodal_dvar")),
61  _mass_frac(getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
62  _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
63  "dPorousFlow_mass_frac_nodal_dvar"))
64 {
66  mooseError("The Dictator proclaims that the number of components in this simulation is ",
68  " whereas you have used the Kernel PorousFlowComponetMassRadioactiveDecay with "
69  "component = ",
71  ". The Dictator does not take such mistakes lightly");
72 }
73 
74 Real
76 {
77  Real mass = 0.0;
78  for (unsigned ph = 0; ph < _num_phases; ++ph)
79  mass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
81 
82  return _test[_i][_qp] * _decay_rate * _porosity[_i] * mass;
83 }
84 
85 Real
87 {
90  return 0.0;
91  return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
92 }
93 
94 Real
96 {
99  return 0.0;
101 }
102 
103 Real
105 {
106  const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
107 
108  // porosity is dependent on variables that are lumped to the nodes,
109  // but it can depend on the gradient
110  // of variables, which are NOT lumped to the nodes, hence:
111  Real dmass = 0.0;
112  for (unsigned ph = 0; ph < _num_phases; ++ph)
113  dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
114  _mass_frac[_i][ph][_fluid_component] * _dporosity_dgradvar[_i][pvar] *
115  _grad_phi[_j][nearest_qp];
116 
117  if (_i != _j)
118  return _test[_i][_qp] * _decay_rate * dmass;
119 
121  for (unsigned ph = 0; ph < _num_phases; ++ph)
122  {
123  dmass += _dfluid_density_dvar[_i][ph][pvar] * _fluid_saturation_nodal[_i][ph] *
124  _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
125  dmass += _fluid_density[_i][ph] * _dfluid_saturation_nodal_dvar[_i][ph][pvar] *
126  _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
127  dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
128  _dmass_frac_dvar[_i][ph][_fluid_component][pvar] * _porosity[_i];
129  dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
130  _mass_frac[_i][ph][_fluid_component] * _dporosity_dvar[_i][pvar];
131  }
132  return _test[_i][_qp] * _decay_rate * dmass;
133 }
const MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _dmass_frac_dvar
d(nodal mass fraction)/d(porous-flow variable)
const MaterialProperty< std::vector< Real > > & _dporosity_dvar
d(porosity)/d(porous-flow variable) - these derivatives will be wrt variables at the nodes ...
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
const unsigned int _fluid_component
the fluid component index
const unsigned int _num_phases
number of fluid phases
const MaterialProperty< std::vector< std::vector< Real > > > & _mass_frac
nodal mass fraction
InputParameters validParams< PorousFlowMassRadioactiveDecay >()
unsigned int numComponents() const
the number of fluid components
const bool _strain_at_nearest_qp
whether the porosity uses the volumetric strain at the closest quadpoint
const MaterialProperty< std::vector< RealGradient > > & _dporosity_dgradvar
d(porosity)/d(grad porous-flow variable) - remember these derivatives will be wrt grad(vars) at qps ...
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_density_dvar
d(nodal fluid density)/d(porous-flow variable)
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_saturation_nodal_dvar
d(nodal fluid saturation)/d(porous-flow variable)
const PorousFlowDictator & _dictator
holds info on the PorousFlow variables
const MaterialProperty< std::vector< Real > > & _fluid_density
nodal fluid density
PorousFlowMassRadioactiveDecay(const InputParameters &parameters)
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
const MaterialProperty< std::vector< Real > > & _fluid_saturation_nodal
nodal fluid saturation
bool notPorousFlowVariable(unsigned int moose_var_num) const
returns true if moose_var_num is not a porous flow variabe
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
the PorousFlow variable number
Real computeQpJac(unsigned int pvar)
Derivative of residual with respect to PorousFlow variable number pvar This is used by both computeQp...
const bool _var_is_porflow_var
whether the Variable for this Kernel is a porous-flow variable according to the Dictator ...
const MaterialProperty< Real > & _porosity
porosity at the nodes, but it can depend on grad(variables) which are actually evaluated at the qps ...