www.mooseframework.org
PorousFlowSink.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 
10 #include "PorousFlowSink.h"
11 
12 #include "MooseVariable.h"
13 
14 #include "libmesh/quadrature.h"
15 
16 #include <iostream>
17 
18 registerMooseObject("PorousFlowApp", PorousFlowSink);
19 
22 {
24  params.addRequiredParam<UserObjectName>(
25  "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
26  params.addParam<unsigned int>("fluid_phase",
27  "If supplied, then this BC will potentially be a function of fluid "
28  "pressure, and you can use mass_fraction_component, use_mobility, "
29  "use_relperm, use_enthalpy and use_energy. If not supplied, then "
30  "this BC can only be a function of temperature");
31  params.addParam<unsigned int>("mass_fraction_component",
32  "The index corresponding to a fluid "
33  "component. If supplied, the flux will "
34  "be multiplied by the nodal mass "
35  "fraction for the component");
36  params.addParam<bool>("use_mobility",
37  false,
38  "If true, then fluxes are multiplied by "
39  "(density*permeability_nn/viscosity), where the "
40  "'_nn' indicates the component normal to the "
41  "boundary. In this case bare_flux is measured in "
42  "Pa.m^-1. This can be used in conjunction with "
43  "other use_*");
44  params.addParam<bool>("use_relperm",
45  false,
46  "If true, then fluxes are multiplied by relative "
47  "permeability. This can be used in conjunction with "
48  "other use_*");
49  params.addParam<bool>("use_enthalpy",
50  false,
51  "If true, then fluxes are multiplied by enthalpy. "
52  "In this case bare_flux is measured in kg.m^-2.s^-1 "
53  "/ (J.kg). This can be used in conjunction with "
54  "other use_*");
55  params.addParam<bool>("use_internal_energy",
56  false,
57  "If true, then fluxes are multiplied by fluid internal energy. "
58  " In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). "
59  " This can be used in conjunction with other use_*");
60  params.addParam<bool>("use_thermal_conductivity",
61  false,
62  "If true, then fluxes are multiplied by "
63  "thermal conductivity projected onto "
64  "the normal direction. This can be "
65  "used in conjunction with other use_*");
66  params.addParam<FunctionName>(
67  "flux_function",
68  1.0,
69  "The flux. The flux is OUT of the medium: hence positive values of "
70  "this function means this BC will act as a SINK, while negative values "
71  "indicate this flux will be a SOURCE. The functional form is useful "
72  "for spatially or temporally varying sinks. Without any use_*, this "
73  "function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case "
74  "with only heat and no fluids)");
75  params.addClassDescription("Applies a flux sink to a boundary.");
76  return params;
77 }
78 
80  : IntegratedBC(parameters),
81  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
82  _involves_fluid(isParamValid("fluid_phase")),
83  _ph(_involves_fluid ? getParam<unsigned int>("fluid_phase") : 0),
84  _use_mass_fraction(isParamValid("mass_fraction_component")),
85  _has_mass_fraction(
86  hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
87  hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
88  "dPorousFlow_mass_frac_nodal_dvar")),
89  _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
90  _use_mobility(getParam<bool>("use_mobility")),
91  _has_mobility(
92  hasMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp") &&
93  hasMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar") &&
94  hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
95  hasMaterialProperty<std::vector<std::vector<Real>>>(
96  "dPorousFlow_fluid_phase_density_nodal_dvar") &&
97  hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
98  hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
99  _use_relperm(getParam<bool>("use_relperm")),
100  _has_relperm(hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
101  hasMaterialProperty<std::vector<std::vector<Real>>>(
102  "dPorousFlow_relative_permeability_nodal_dvar")),
103  _use_enthalpy(getParam<bool>("use_enthalpy")),
104  _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
105  hasMaterialProperty<std::vector<std::vector<Real>>>(
106  "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
107  _use_internal_energy(getParam<bool>("use_internal_energy")),
108  _has_internal_energy(
109  hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
110  hasMaterialProperty<std::vector<std::vector<Real>>>(
111  "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
112  _use_thermal_conductivity(getParam<bool>("use_thermal_conductivity")),
113  _has_thermal_conductivity(
114  hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
115  hasMaterialProperty<std::vector<RealTensorValue>>(
116  "dPorousFlow_thermal_conductivity_qp_dvar")),
117  _m_func(getFunction("flux_function")),
118  _permeability(_has_mobility
119  ? &getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")
120  : nullptr),
121  _dpermeability_dvar(_has_mobility ? &getMaterialProperty<std::vector<RealTensorValue>>(
122  "dPorousFlow_permeability_qp_dvar")
123  : nullptr),
124  _dpermeability_dgradvar(_has_mobility
125  ? &getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
126  "dPorousFlow_permeability_qp_dgradvar")
127  : nullptr),
128  _fluid_density_node(_has_mobility ? &getMaterialProperty<std::vector<Real>>(
129  "PorousFlow_fluid_phase_density_nodal")
130  : nullptr),
131  _dfluid_density_node_dvar(_has_mobility ? &getMaterialProperty<std::vector<std::vector<Real>>>(
132  "dPorousFlow_fluid_phase_density_nodal_dvar")
133  : nullptr),
134  _fluid_viscosity(_has_mobility
135  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
136  : nullptr),
137  _dfluid_viscosity_dvar(_has_mobility ? &getMaterialProperty<std::vector<std::vector<Real>>>(
138  "dPorousFlow_viscosity_nodal_dvar")
139  : nullptr),
140  _relative_permeability(_has_relperm ? &getMaterialProperty<std::vector<Real>>(
141  "PorousFlow_relative_permeability_nodal")
142  : nullptr),
143  _drelative_permeability_dvar(_has_relperm
144  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
145  "dPorousFlow_relative_permeability_nodal_dvar")
146  : nullptr),
147  _mass_fractions(_has_mass_fraction ? &getMaterialProperty<std::vector<std::vector<Real>>>(
148  "PorousFlow_mass_frac_nodal")
149  : nullptr),
150  _dmass_fractions_dvar(_has_mass_fraction
151  ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
152  "dPorousFlow_mass_frac_nodal_dvar")
153  : nullptr),
154  _enthalpy(_has_enthalpy ? &getMaterialPropertyByName<std::vector<Real>>(
155  "PorousFlow_fluid_phase_enthalpy_nodal")
156  : nullptr),
157  _denthalpy_dvar(_has_enthalpy ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
158  "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
159  : nullptr),
160  _internal_energy(_has_internal_energy ? &getMaterialPropertyByName<std::vector<Real>>(
161  "PorousFlow_fluid_phase_internal_energy_nodal")
162  : nullptr),
163  _dinternal_energy_dvar(_has_internal_energy
164  ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
165  "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
166  : nullptr),
167  _thermal_conductivity(_has_thermal_conductivity ? &getMaterialProperty<RealTensorValue>(
168  "PorousFlow_thermal_conductivity_qp")
169  : nullptr),
170  _dthermal_conductivity_dvar(_has_thermal_conductivity
171  ? &getMaterialProperty<std::vector<RealTensorValue>>(
172  "dPorousFlow_thermal_conductivity_qp_dvar")
173  : nullptr),
174  _perm_derivs(_dictator.usePermDerivs())
175 {
177  paramError("fluid_phase",
178  "The Dictator proclaims that the maximum phase index in this simulation is ",
179  _dictator.numPhases() - 1,
180  " whereas you have used ",
181  _ph,
182  ". Remember that indexing starts at 0. You must try harder.");
183 
186  mooseError("PorousFlowSink: To use_mass_fraction, use_mobility, use_relperm, use_enthalpy or "
187  "use_internal_energy, you must provide a fluid phase number");
188 
190  paramError("mass_fraction_component",
191  "The Dictator declares that the maximum fluid component index is ",
192  _dictator.numComponents() - 1,
193  ", but you have set mass_fraction_component to ",
194  _sp,
195  ". Remember that indexing starts at 0. Please be assured that the Dictator has "
196  "noted your error.");
197 
199  mooseError("PorousFlowSink: You have used the use_mass_fraction flag, but you have no "
200  "mass_fraction Material");
201 
203  mooseError("PorousFlowSink: You have used the use_mobility flag, but there are not the "
204  "required Materials for this");
205 
206  if (_use_relperm && !_has_relperm)
207  mooseError(
208  "PorousFlowSink: You have used the use_relperm flag, but you have no relperm Material");
209 
211  mooseError(
212  "PorousFlowSink: You have used the use_enthalpy flag, but you have no enthalpy Material");
213 
215  mooseError("PorousFlowSink: You have used the use_internal_energy flag, but you have no "
216  "internal_energy Material");
217 
219  mooseError("PorousFlowSink: You have used the use_thermal_conductivity flag, but you have no "
220  "thermal_conductivity Material");
221 }
222 
223 Real
225 {
226  Real flux = _test[_i][_qp] * multiplier();
227  if (_use_mobility)
228  {
229  const Real k =
230  ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp]; // do not upwind permeability
231  flux *= (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
232  }
233  if (_use_relperm)
234  flux *= (*_relative_permeability)[_i][_ph];
235  if (_use_mass_fraction)
236  flux *= (*_mass_fractions)[_i][_ph][_sp];
237  if (_use_enthalpy)
238  flux *= (*_enthalpy)[_i][_ph];
240  flux *= (*_internal_energy)[_i][_ph];
242  flux *= ((*_thermal_conductivity)[_qp] * _normals[_qp]) *
243  _normals[_qp]; // do not upwind thermal_conductivity
244 
245  return flux;
246 }
247 
248 Real
250 {
251  return jac(_var.number());
252 }
253 
254 Real
256 {
257  return jac(jvar);
258 }
259 
260 Real
261 PorousFlowSink::jac(unsigned int jvar) const
262 {
264  return 0.0;
265  const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
266 
267  // For _i != _j, note:
268  // since the only non-upwinded contribution to the residual is
269  // from the permeability and thermal_conductivity, the only contribution
270  // of the residual at node _i from changing jvar at node _j is through
271  // the derivative of permeability or thermal_conductivity
272 
273  Real flux = _test[_i][_qp] * multiplier();
274  Real deriv = _test[_i][_qp] * (_i != _j ? 0.0 : dmultiplier_dvar(pvar));
275 
276  if (_use_mobility)
277  {
278  const Real k = ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp];
279  const Real mob = (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
280 
281  Real mobprime = 0.0;
282  if (_perm_derivs)
283  {
284  RealTensorValue ktprime = (*_dpermeability_dvar)[_qp][pvar] * _phi[_j][_qp];
285  for (const auto i : make_range(Moose::dim))
286  ktprime += (*_dpermeability_dgradvar)[_qp][i][pvar] * _grad_phi[_j][_qp](i);
287  const Real kprime = (ktprime * _normals[_qp]) * _normals[_qp];
288 
289  mobprime += (*_fluid_density_node)[_i][_ph] * kprime / (*_fluid_viscosity)[_i][_ph];
290  }
291 
292  mobprime +=
293  (_i != _j
294  ? 0.0
295  : (*_dfluid_density_node_dvar)[_i][_ph][pvar] * k / (*_fluid_viscosity)[_i][_ph] -
296  (*_fluid_density_node)[_i][_ph] * k * (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
297  std::pow((*_fluid_viscosity)[_i][_ph], 2));
298  deriv = mob * deriv + mobprime * flux;
299  flux *= mob;
300  }
301  if (_use_relperm)
302  {
303  const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
304  deriv = (*_relative_permeability)[_i][_ph] * deriv + relperm_prime * flux;
305  flux *= (*_relative_permeability)[_i][_ph];
306  }
307  if (_use_mass_fraction)
308  {
309  const Real mf_prime = (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
310  deriv = (*_mass_fractions)[_i][_ph][_sp] * deriv + mf_prime * flux;
311  flux *= (*_mass_fractions)[_i][_ph][_sp];
312  }
313  if (_use_enthalpy)
314  {
315  const Real en_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
316  deriv = (*_enthalpy)[_i][_ph] * deriv + en_prime * flux;
317  flux *= (*_enthalpy)[_i][_ph];
318  }
320  {
321  const Real ie_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
322  deriv = (*_internal_energy)[_i][_ph] * deriv + ie_prime * flux;
323  flux *= (*_internal_energy)[_i][_ph];
324  }
326  {
327  const Real tc = ((*_thermal_conductivity)[_qp] * _normals[_qp]) * _normals[_qp];
328  const RealTensorValue tctprime = (*_dthermal_conductivity_dvar)[_qp][pvar] * _phi[_j][_qp];
329  const Real tcprime = (tctprime * _normals[_qp]) * _normals[_qp];
330  deriv = tc * deriv + tcprime * flux;
331  // don't need this: flux *= tc;
332  }
333  return deriv;
334 }
335 
336 Real
338 {
339  return _m_func.value(_t, _q_point[_qp]);
340 }
341 
342 Real
343 PorousFlowSink::dmultiplier_dvar(unsigned int /*pvar*/) const
344 {
345  return 0.0;
346 }
const VariableTestValue & _test
const bool _has_enthalpy
Whether there is an "enthalpy" Material. This is just for error checking.
const bool _involves_fluid
Whether this BC involves fluid (whether the user has supplied a fluid phase number) ...
const unsigned int _sp
The component number (only used if _use_mass_fraction==true)
const unsigned int _ph
The phase number.
unsigned int _j
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > *const _dpermeability_dgradvar
d(Permeability)/d(grad(PorousFlow variable))
bool notPorousFlowVariable(unsigned int moose_var_num) const
Returns true if moose_var_num is not a porous flow variabe.
const MooseArray< Point > & _normals
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
unsigned int number() const
const bool _has_internal_energy
Whether there is an "internal_energy" Material. This is just for error checking.
const MaterialProperty< std::vector< Real > > *const _fluid_viscosity
Viscosity of each component in each phase.
const bool _has_thermal_conductivity
Whether there is an "thermal_conductivity" Material. This is just for error checking.
unsigned int numComponents() const
The number of fluid components.
static InputParameters validParams()
const bool _has_mobility
Whether there are Materials that can form "mobility". This is just for error checking.
unsigned int _i
static constexpr std::size_t dim
const Function & _m_func
The flux.
const bool _use_mass_fraction
Whether the flux will be multiplied by the mass fraction.
virtual Real multiplier() const
The flux gets multiplied by this quantity.
const VariablePhiValue & _phi
const bool _has_relperm
Whether there is a "relperm" Material. This is just for error checking.
void addRequiredParam(const std::string &name, const std::string &doc_string)
unsigned int _qp
const VariablePhiGradient & _grad_phi
registerMooseObject("PorousFlowApp", PorousFlowSink)
const MooseArray< Point > & _q_point
TensorValue< Real > RealTensorValue
const bool _use_relperm
Whether to multiply the sink flux by relative permeability.
Applies a flux sink to a boundary.
Real jac(unsigned int jvar) const
Derivative of residual with respect to the jvar variable.
Real deriv(unsigned n, unsigned alpha, unsigned beta, Real x)
const bool _use_mobility
Whether to multiply the sink flux by permeability*density/viscosity.
virtual Real computeQpJacobian() override
virtual Real dmultiplier_dvar(unsigned int pvar) const
d(multiplier)/d(Porous flow variable pvar)
void paramError(const std::string &param, Args... args) const
unsigned int numPhases() const
The number of fluid phases.
const bool _use_enthalpy
Whether to multiply the sink flux by enthalpy.
const bool _perm_derivs
Flag to check whether permeabiity derivatives are non-zero.
MooseVariable & _var
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
const bool _has_mass_fraction
Whether there is a "mass_fraction" Material. This is just for error checking.
const bool _use_internal_energy
Whether to multiply the sink flux by internal_energy.
virtual Real computeQpResidual() override
const bool _use_thermal_conductivity
Whether to multiply the sink flux by thermal_conductivity.
IntRange< T > make_range(T beg, T end)
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
The PorousFlow variable number.
virtual Real value(Real t, const Point &p) const
PorousFlowSink(const InputParameters &parameters)
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
MooseUnits pow(const MooseUnits &, int)
static const std::string k
Definition: NS.h:124
void ErrorVector unsigned int
static InputParameters validParams()