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