www.mooseframework.org
PorousFlowMassTimeDerivative.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 #include "MooseVariable.h"
13 
14 #include "libmesh/quadrature.h"
15 
16 #include <limits>
17 
19 
22 {
24  params.addParam<bool>("strain_at_nearest_qp",
25  false,
26  "When calculating nodal porosity that depends on strain, use the strain at "
27  "the nearest quadpoint. This adds a small extra computational burden, and "
28  "is not necessary for simulations involving only linear lagrange elements. "
29  " If you set this to true, you will also want to set the same parameter to "
30  "true for related Kernels and Materials");
31  params.addParam<std::string>(
32  "base_name",
33  "For mechanically-coupled systems, this Kernel will depend on the volumetric strain. "
34  "base_name should almost always be the same base_name as given to the TensorMechanics object "
35  "that computes strain. Supplying a base_name to this Kernel but not defining an associated "
36  "TensorMechanics strain calculator means that this Kernel will not depend on volumetric "
37  "strain. That could be useful when models contain solid mechanics that is not coupled to "
38  "porous flow, for example");
39  params.addParam<bool>(
40  "multiply_by_density",
41  true,
42  "If true, then this Kernel represents the time derivative of the fluid mass. If false, then "
43  "this Kernel represents the time derivative of the fluid volume (care must then be taken "
44  "when using other PorousFlow objects, such as the PorousFlowFluidMass postprocessor).");
45  params.addParam<unsigned int>(
46  "fluid_component", 0, "The index corresponding to the component for this kernel");
47  params.addRequiredParam<UserObjectName>(
48  "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
49  params.set<bool>("use_displaced_mesh") = false;
50  params.suppressParameter<bool>("use_displaced_mesh");
51  params.addClassDescription("Derivative of fluid-component mass with respect to time. Mass "
52  "lumping to the nodes is used.");
53  return params;
54 }
55 
57  : TimeKernel(parameters),
58  _fluid_component(getParam<unsigned int>("fluid_component")),
59  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
60  _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
61  _num_phases(_dictator.numPhases()),
62  _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
63  _multiply_by_density(getParam<bool>("multiply_by_density")),
64  _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
65  _has_total_strain(hasMaterialProperty<RankTwoTensor>(_base_name + "total_strain")),
66  _total_strain_old(_has_total_strain
67  ? &getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")
68  : nullptr),
69  _porosity(getMaterialProperty<Real>("PorousFlow_porosity_nodal")),
70  _porosity_old(getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")),
71  _dporosity_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_nodal_dvar")),
72  _dporosity_dgradvar(
73  getMaterialProperty<std::vector<RealGradient>>("dPorousFlow_porosity_nodal_dgradvar")),
74  _nearest_qp(_strain_at_nearest_qp
75  ? &getMaterialProperty<unsigned int>("PorousFlow_nearestqp_nodal")
76  : nullptr),
77  _fluid_density(_multiply_by_density ? &getMaterialProperty<std::vector<Real>>(
78  "PorousFlow_fluid_phase_density_nodal")
79  : nullptr),
80  _fluid_density_old(_multiply_by_density ? &getMaterialPropertyOld<std::vector<Real>>(
81  "PorousFlow_fluid_phase_density_nodal")
82  : nullptr),
83  _dfluid_density_dvar(_multiply_by_density
84  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
85  "dPorousFlow_fluid_phase_density_nodal_dvar")
86  : nullptr),
87  _fluid_saturation_nodal(getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")),
88  _fluid_saturation_nodal_old(
89  getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")),
90  _dfluid_saturation_nodal_dvar(
91  getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_saturation_nodal_dvar")),
92  _mass_frac(getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
93  _mass_frac_old(
94  getMaterialPropertyOld<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
95  _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
96  "dPorousFlow_mass_frac_nodal_dvar"))
97 {
99  paramError(
100  "fluid_component",
101  "The Dictator proclaims that the maximum fluid component index in this simulation is ",
102  _dictator.numComponents() - 1,
103  " whereas you have used ",
105  ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
106 }
107 
108 Real
110 {
111  Real mass = 0.0;
112  Real mass_old = 0.0;
113  for (unsigned ph = 0; ph < _num_phases; ++ph)
114  {
115  const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
116  mass += dens * _fluid_saturation_nodal[_i][ph] * _mass_frac[_i][ph][_fluid_component];
117  const Real dens_old = (_multiply_by_density ? (*_fluid_density_old)[_i][ph] : 1.0);
118  mass_old +=
120  }
121  const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
122 
123  return _test[_i][_qp] * (1.0 + strain) * (_porosity[_i] * mass - _porosity_old[_i] * mass_old) /
124  _dt;
125 }
126 
127 Real
129 {
130  // If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
131  if (!_var_is_porflow_var)
132  return 0.0;
134 }
135 
136 Real
138 {
139  // If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
141  return 0.0;
143 }
144 
145 Real
147 {
148  const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
149 
150  const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
151 
152  // porosity is dependent on variables that are lumped to the nodes,
153  // but it can depend on the gradient
154  // of variables, which are NOT lumped to the nodes, hence:
155  Real dmass = 0.0;
156  for (unsigned ph = 0; ph < _num_phases; ++ph)
157  {
158  const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
159  dmass += dens * _fluid_saturation_nodal[_i][ph] * _mass_frac[_i][ph][_fluid_component] *
160  _dporosity_dgradvar[_i][pvar] * _grad_phi[_j][nearest_qp];
161  }
162 
163  if (_i != _j)
164  return _test[_i][_qp] * (1.0 + strain) * dmass / _dt;
165 
166  // As the fluid mass is lumped to the nodes, only non-zero terms are for _i==_j
167  for (unsigned ph = 0; ph < _num_phases; ++ph)
168  {
170  dmass += (*_dfluid_density_dvar)[_i][ph][pvar] * _fluid_saturation_nodal[_i][ph] *
172  const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
173  dmass += dens * _dfluid_saturation_nodal_dvar[_i][ph][pvar] *
175  dmass += dens * _fluid_saturation_nodal[_i][ph] *
177  dmass += dens * _fluid_saturation_nodal[_i][ph] * _mass_frac[_i][ph][_fluid_component] *
178  _dporosity_dvar[_i][pvar];
179  }
180  return _test[_i][_qp] * (1.0 + strain) * dmass / _dt;
181 }
registerMooseObject("PorousFlowApp", PorousFlowMassTimeDerivative)
const bool _strain_at_nearest_qp
Whether the porosity uses the volumetric strain at the closest quadpoint.
bool notPorousFlowVariable(unsigned int moose_var_num) const
Returns true if moose_var_num is not a porous flow variabe.
MooseVariable & _var
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const MaterialProperty< Real > & _porosity_old
Old value of porosity.
unsigned int number() const
const VariablePhiGradient & _grad_phi
T & set(const std::string &name, bool quiet_mode=false)
const MaterialProperty< Real > & _porosity
Porosity at the nodes, but it can depend on grad(variables) which are actually evaluated at the qps...
unsigned int numComponents() const
The number of fluid components.
Real computeQpJac(unsigned int pvar)
Derivative of residual with respect to PorousFlow variable number pvar This is used by both computeQp...
Real & _dt
const MaterialProperty< std::vector< Real > > & _fluid_saturation_nodal
Nodal fluid saturation.
void addRequiredParam(const std::string &name, const std::string &doc_string)
void suppressParameter(const std::string &name)
const MaterialProperty< std::vector< RealGradient > > & _dporosity_dgradvar
d(porosity)/d(grad PorousFlow variable) - remember these derivatives will be wrt grad(vars) at qps ...
const VariableTestValue & _test
const bool _var_is_porflow_var
Whether the Variable for this Kernel is a PorousFlow variable according to the Dictator.
unsigned int _i
PorousFlowMassTimeDerivative(const InputParameters &parameters)
const MaterialProperty< std::vector< Real > > & _fluid_saturation_nodal_old
Old value of fluid saturation.
void paramError(const std::string &param, Args... args) const
const MaterialProperty< std::vector< std::vector< Real > > > & _mass_frac
Nodal mass fraction.
const PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
const bool _has_total_strain
Whether there is a Material called _base_name_total_strain.
const MaterialProperty< std::vector< std::vector< Real > > > & _mass_frac_old
Old value of nodal mass fraction.
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_saturation_nodal_dvar
d(nodal fluid saturation)/d(PorousFlow variable)
unsigned int _j
const MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _dmass_frac_dvar
d(nodal mass fraction)/d(PorousFlow variable)
const bool _multiply_by_density
Whether to multiply by density: if true then this Kernel is computing the time-derivative of fluid ma...
Kernel = (mass_component - mass_component_old)/dt where mass_component = porosity*sum_phases(density_...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
const unsigned int _num_phases
Number of fluid phases.
const unsigned int _fluid_component
The fluid component index.
void addClassDescription(const std::string &doc_string)
const MaterialProperty< std::vector< Real > > & _dporosity_dvar
d(porosity)/d(PorousFlow variable) - these derivatives will be wrt variables at the nodes ...
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
The PorousFlow variable number.
static InputParameters validParams()
void ErrorVector unsigned int
unsigned int _qp