Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 : InputParameters
21 3905 : PorousFlowSink::validParams()
22 : {
23 3905 : InputParameters params = IntegratedBC::validParams();
24 7810 : params.addRequiredParam<UserObjectName>(
25 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
26 7810 : 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 7810 : 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 7810 : params.addParam<bool>("use_mobility",
37 7810 : 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 7810 : params.addParam<bool>("use_relperm",
45 7810 : false,
46 : "If true, then fluxes are multiplied by relative "
47 : "permeability. This can be used in conjunction with "
48 : "other use_*");
49 7810 : params.addParam<bool>("use_enthalpy",
50 7810 : 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 7810 : params.addParam<bool>("use_internal_energy",
56 7810 : 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 7810 : params.addParam<bool>("use_thermal_conductivity",
61 7810 : 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 7810 : params.addParam<FunctionName>(
67 : "flux_function",
68 7810 : 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 3905 : params.addClassDescription("Applies a flux sink to a boundary.");
76 3905 : return params;
77 0 : }
78 :
79 2182 : PorousFlowSink::PorousFlowSink(const InputParameters & parameters)
80 : : IntegratedBC(parameters),
81 2182 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
82 4364 : _involves_fluid(isParamValid("fluid_phase")),
83 3883 : _ph(_involves_fluid ? getParam<unsigned int>("fluid_phase") : 0),
84 4364 : _use_mass_fraction(isParamValid("mass_fraction_component")),
85 2182 : _has_mass_fraction(
86 2182 : hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
87 4287 : hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
88 : "dPorousFlow_mass_frac_nodal_dvar")),
89 3309 : _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
90 4364 : _use_mobility(getParam<bool>("use_mobility")),
91 2182 : _has_mobility(
92 4069 : hasMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp") &&
93 5956 : hasMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar") &&
94 5868 : hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
95 3981 : hasMaterialProperty<std::vector<std::vector<Real>>>(
96 1799 : "dPorousFlow_fluid_phase_density_nodal_dvar") &&
97 6163 : hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
98 3981 : hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
99 4364 : _use_relperm(getParam<bool>("use_relperm")),
100 2182 : _has_relperm(hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
101 4069 : hasMaterialProperty<std::vector<std::vector<Real>>>(
102 : "dPorousFlow_relative_permeability_nodal_dvar")),
103 4364 : _use_enthalpy(getParam<bool>("use_enthalpy")),
104 2182 : _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
105 4208 : hasMaterialProperty<std::vector<std::vector<Real>>>(
106 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
107 4364 : _use_internal_energy(getParam<bool>("use_internal_energy")),
108 2182 : _has_internal_energy(
109 2182 : hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
110 4208 : hasMaterialProperty<std::vector<std::vector<Real>>>(
111 : "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
112 4364 : _use_thermal_conductivity(getParam<bool>("use_thermal_conductivity")),
113 2182 : _has_thermal_conductivity(
114 2182 : hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
115 2452 : hasMaterialProperty<std::vector<RealTensorValue>>(
116 : "dPorousFlow_thermal_conductivity_qp_dvar")),
117 2182 : _m_func(getFunction("flux_function")),
118 4364 : _permeability(_has_mobility
119 3981 : ? &getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")
120 : : nullptr),
121 3981 : _dpermeability_dvar(_has_mobility ? &getMaterialProperty<std::vector<RealTensorValue>>(
122 : "dPorousFlow_permeability_qp_dvar")
123 : : nullptr),
124 4364 : _dpermeability_dgradvar(_has_mobility
125 3981 : ? &getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
126 : "dPorousFlow_permeability_qp_dgradvar")
127 : : nullptr),
128 3981 : _fluid_density_node(_has_mobility ? &getMaterialProperty<std::vector<Real>>(
129 : "PorousFlow_fluid_phase_density_nodal")
130 : : nullptr),
131 3981 : _dfluid_density_node_dvar(_has_mobility ? &getMaterialProperty<std::vector<std::vector<Real>>>(
132 : "dPorousFlow_fluid_phase_density_nodal_dvar")
133 : : nullptr),
134 4364 : _fluid_viscosity(_has_mobility
135 3981 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
136 : : nullptr),
137 3981 : _dfluid_viscosity_dvar(_has_mobility ? &getMaterialProperty<std::vector<std::vector<Real>>>(
138 : "dPorousFlow_viscosity_nodal_dvar")
139 : : nullptr),
140 4069 : _relative_permeability(_has_relperm ? &getMaterialProperty<std::vector<Real>>(
141 : "PorousFlow_relative_permeability_nodal")
142 : : nullptr),
143 4364 : _drelative_permeability_dvar(_has_relperm
144 4069 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
145 : "dPorousFlow_relative_permeability_nodal_dvar")
146 : : nullptr),
147 4287 : _mass_fractions(_has_mass_fraction ? &getMaterialProperty<std::vector<std::vector<Real>>>(
148 : "PorousFlow_mass_frac_nodal")
149 : : nullptr),
150 4364 : _dmass_fractions_dvar(_has_mass_fraction
151 4287 : ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
152 : "dPorousFlow_mass_frac_nodal_dvar")
153 : : nullptr),
154 4208 : _enthalpy(_has_enthalpy ? &getMaterialPropertyByName<std::vector<Real>>(
155 : "PorousFlow_fluid_phase_enthalpy_nodal")
156 : : nullptr),
157 4208 : _denthalpy_dvar(_has_enthalpy ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
158 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
159 : : nullptr),
160 4208 : _internal_energy(_has_internal_energy ? &getMaterialPropertyByName<std::vector<Real>>(
161 : "PorousFlow_fluid_phase_internal_energy_nodal")
162 : : nullptr),
163 4364 : _dinternal_energy_dvar(_has_internal_energy
164 4208 : ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
165 : "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
166 : : nullptr),
167 2452 : _thermal_conductivity(_has_thermal_conductivity ? &getMaterialProperty<RealTensorValue>(
168 : "PorousFlow_thermal_conductivity_qp")
169 : : nullptr),
170 4364 : _dthermal_conductivity_dvar(_has_thermal_conductivity
171 2452 : ? &getMaterialProperty<std::vector<RealTensorValue>>(
172 : "dPorousFlow_thermal_conductivity_qp_dvar")
173 : : nullptr),
174 4364 : _perm_derivs(_dictator.usePermDerivs())
175 : {
176 2182 : if (_involves_fluid && _ph >= _dictator.numPhases())
177 0 : paramError("fluid_phase",
178 : "The Dictator proclaims that the maximum phase index in this simulation is ",
179 0 : _dictator.numPhases() - 1,
180 : " whereas you have used ",
181 0 : _ph,
182 : ". Remember that indexing starts at 0. You must try harder.");
183 :
184 2182 : if (!_involves_fluid && (_use_mass_fraction || _use_mobility || _use_relperm || _use_enthalpy ||
185 481 : _use_internal_energy))
186 0 : 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 :
189 2182 : if (_use_mass_fraction && _sp >= _dictator.numComponents())
190 0 : paramError("mass_fraction_component",
191 : "The Dictator declares that the maximum fluid component index is ",
192 0 : _dictator.numComponents() - 1,
193 : ", but you have set mass_fraction_component to ",
194 0 : _sp,
195 : ". Remember that indexing starts at 0. Please be assured that the Dictator has "
196 : "noted your error.");
197 :
198 2182 : if (_use_mass_fraction && !_has_mass_fraction)
199 0 : mooseError("PorousFlowSink: You have used the use_mass_fraction flag, but you have no "
200 : "mass_fraction Material");
201 :
202 2182 : if (_use_mobility && !_has_mobility)
203 0 : mooseError("PorousFlowSink: You have used the use_mobility flag, but there are not the "
204 : "required Materials for this");
205 :
206 2182 : if (_use_relperm && !_has_relperm)
207 0 : mooseError(
208 : "PorousFlowSink: You have used the use_relperm flag, but you have no relperm Material");
209 :
210 2182 : if (_use_enthalpy && !_has_enthalpy)
211 0 : mooseError(
212 : "PorousFlowSink: You have used the use_enthalpy flag, but you have no enthalpy Material");
213 :
214 2182 : if (_use_internal_energy && !_has_internal_energy)
215 0 : mooseError("PorousFlowSink: You have used the use_internal_energy flag, but you have no "
216 : "internal_energy Material");
217 :
218 2182 : if (_use_thermal_conductivity && !_has_thermal_conductivity)
219 0 : mooseError("PorousFlowSink: You have used the use_thermal_conductivity flag, but you have no "
220 : "thermal_conductivity Material");
221 2182 : }
222 :
223 : Real
224 2010136 : PorousFlowSink::computeQpResidual()
225 : {
226 2010136 : Real flux = _test[_i][_qp] * multiplier();
227 2010136 : if (_use_mobility)
228 : {
229 : const Real k =
230 745610 : ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp]; // do not upwind permeability
231 745610 : flux *= (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
232 : }
233 2010136 : if (_use_relperm)
234 769928 : flux *= (*_relative_permeability)[_i][_ph];
235 2010136 : if (_use_mass_fraction)
236 754686 : flux *= (*_mass_fractions)[_i][_ph][_sp];
237 2010136 : if (_use_enthalpy)
238 141216 : flux *= (*_enthalpy)[_i][_ph];
239 2010136 : if (_use_internal_energy)
240 42712 : flux *= (*_internal_energy)[_i][_ph];
241 2010136 : if (_use_thermal_conductivity)
242 42240 : flux *= ((*_thermal_conductivity)[_qp] * _normals[_qp]) *
243 42240 : _normals[_qp]; // do not upwind thermal_conductivity
244 :
245 2010136 : return flux;
246 : }
247 :
248 : Real
249 7979156 : PorousFlowSink::computeQpJacobian()
250 : {
251 7979156 : return jac(_var.number());
252 : }
253 :
254 : Real
255 1185508 : PorousFlowSink::computeQpOffDiagJacobian(unsigned int jvar)
256 : {
257 1185508 : return jac(jvar);
258 : }
259 :
260 : Real
261 9164664 : PorousFlowSink::jac(unsigned int jvar) const
262 : {
263 9164664 : if (_dictator.notPorousFlowVariable(jvar))
264 : return 0.0;
265 9164664 : 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 9164664 : Real flux = _test[_i][_qp] * multiplier();
274 9164664 : Real deriv = _test[_i][_qp] * (_i != _j ? 0.0 : dmultiplier_dvar(pvar));
275 :
276 9164664 : if (_use_mobility)
277 : {
278 853468 : const Real k = ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp];
279 853468 : const Real mob = (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
280 :
281 : Real mobprime = 0.0;
282 853468 : if (_perm_derivs)
283 : {
284 23040 : RealTensorValue ktprime = (*_dpermeability_dvar)[_qp][pvar] * _phi[_j][_qp];
285 92160 : for (const auto i : make_range(Moose::dim))
286 138240 : ktprime += (*_dpermeability_dgradvar)[_qp][i][pvar] * _grad_phi[_j][_qp](i);
287 23040 : const Real kprime = (ktprime * _normals[_qp]) * _normals[_qp];
288 :
289 23040 : mobprime += (*_fluid_density_node)[_i][_ph] * kprime / (*_fluid_viscosity)[_i][_ph];
290 : }
291 :
292 853468 : mobprime +=
293 853468 : (_i != _j
294 853468 : ? 0.0
295 160134 : : (*_dfluid_density_node_dvar)[_i][_ph][pvar] * k / (*_fluid_viscosity)[_i][_ph] -
296 160134 : (*_fluid_density_node)[_i][_ph] * k * (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
297 : std::pow((*_fluid_viscosity)[_i][_ph], 2));
298 853468 : deriv = mob * deriv + mobprime * flux;
299 853468 : flux *= mob;
300 : }
301 9164664 : if (_use_relperm)
302 : {
303 1348832 : const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
304 1348832 : deriv = (*_relative_permeability)[_i][_ph] * deriv + relperm_prime * flux;
305 1348832 : flux *= (*_relative_permeability)[_i][_ph];
306 : }
307 9164664 : if (_use_mass_fraction)
308 : {
309 1190396 : const Real mf_prime = (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
310 1190396 : deriv = (*_mass_fractions)[_i][_ph][_sp] * deriv + mf_prime * flux;
311 1190396 : flux *= (*_mass_fractions)[_i][_ph][_sp];
312 : }
313 9164664 : if (_use_enthalpy)
314 : {
315 92072 : const Real en_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
316 92072 : deriv = (*_enthalpy)[_i][_ph] * deriv + en_prime * flux;
317 92072 : flux *= (*_enthalpy)[_i][_ph];
318 : }
319 9164664 : if (_use_internal_energy)
320 : {
321 30080 : const Real ie_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
322 30080 : deriv = (*_internal_energy)[_i][_ph] * deriv + ie_prime * flux;
323 30080 : flux *= (*_internal_energy)[_i][_ph];
324 : }
325 9164664 : if (_use_thermal_conductivity)
326 : {
327 26880 : const Real tc = ((*_thermal_conductivity)[_qp] * _normals[_qp]) * _normals[_qp];
328 26880 : const RealTensorValue tctprime = (*_dthermal_conductivity_dvar)[_qp][pvar] * _phi[_j][_qp];
329 26880 : const Real tcprime = (tctprime * _normals[_qp]) * _normals[_qp];
330 26880 : deriv = tc * deriv + tcprime * flux;
331 : // don't need this: flux *= tc;
332 : }
333 : return deriv;
334 : }
335 :
336 : Real
337 11062170 : PorousFlowSink::multiplier() const
338 : {
339 11062170 : return _m_func.value(_t, _q_point[_qp]);
340 : }
341 :
342 : Real
343 1211788 : PorousFlowSink::dmultiplier_dvar(unsigned int /*pvar*/) const
344 : {
345 1211788 : return 0.0;
346 : }
|
