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 1979 : PorousFlowSink::validParams()
22 : {
23 1979 : InputParameters params = IntegratedBC::validParams();
24 3958 : params.addRequiredParam<UserObjectName>(
25 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
26 3958 : 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 3958 : 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 3958 : params.addParam<bool>("use_mobility",
37 3958 : 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 3958 : params.addParam<bool>("use_relperm",
45 3958 : false,
46 : "If true, then fluxes are multiplied by relative "
47 : "permeability. This can be used in conjunction with "
48 : "other use_*");
49 3958 : params.addParam<bool>("use_enthalpy",
50 3958 : 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 3958 : params.addParam<bool>("use_internal_energy",
56 3958 : 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 3958 : params.addParam<bool>("use_thermal_conductivity",
61 3958 : 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 3958 : params.addParam<FunctionName>(
67 : "flux_function",
68 3958 : 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 1979 : params.addClassDescription("Applies a flux sink to a boundary.");
76 1979 : return params;
77 0 : }
78 :
79 1080 : PorousFlowSink::PorousFlowSink(const InputParameters & parameters)
80 : : IntegratedBC(parameters),
81 1080 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
82 2160 : _involves_fluid(isParamValid("fluid_phase")),
83 1927 : _ph(_involves_fluid ? getParam<unsigned int>("fluid_phase") : 0),
84 2160 : _use_mass_fraction(isParamValid("mass_fraction_component")),
85 1080 : _has_mass_fraction(
86 1080 : hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
87 2115 : hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
88 : "dPorousFlow_mass_frac_nodal_dvar")),
89 1645 : _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
90 2160 : _use_mobility(getParam<bool>("use_mobility")),
91 1080 : _has_mobility(
92 2013 : hasMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp") &&
93 2946 : hasMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar") &&
94 2906 : hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
95 1973 : hasMaterialProperty<std::vector<std::vector<Real>>>(
96 893 : "dPorousFlow_fluid_phase_density_nodal_dvar") &&
97 3053 : hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
98 1973 : hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
99 2160 : _use_relperm(getParam<bool>("use_relperm")),
100 1080 : _has_relperm(hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
101 2013 : hasMaterialProperty<std::vector<std::vector<Real>>>(
102 : "dPorousFlow_relative_permeability_nodal_dvar")),
103 2160 : _use_enthalpy(getParam<bool>("use_enthalpy")),
104 1080 : _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
105 2076 : hasMaterialProperty<std::vector<std::vector<Real>>>(
106 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
107 2160 : _use_internal_energy(getParam<bool>("use_internal_energy")),
108 1080 : _has_internal_energy(
109 1080 : hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
110 2076 : hasMaterialProperty<std::vector<std::vector<Real>>>(
111 : "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
112 2160 : _use_thermal_conductivity(getParam<bool>("use_thermal_conductivity")),
113 1080 : _has_thermal_conductivity(
114 1080 : hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
115 1230 : hasMaterialProperty<std::vector<RealTensorValue>>(
116 : "dPorousFlow_thermal_conductivity_qp_dvar")),
117 1080 : _m_func(getFunction("flux_function")),
118 2160 : _permeability(_has_mobility
119 1973 : ? &getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")
120 : : nullptr),
121 1973 : _dpermeability_dvar(_has_mobility ? &getMaterialProperty<std::vector<RealTensorValue>>(
122 : "dPorousFlow_permeability_qp_dvar")
123 : : nullptr),
124 2160 : _dpermeability_dgradvar(_has_mobility
125 1973 : ? &getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
126 : "dPorousFlow_permeability_qp_dgradvar")
127 : : nullptr),
128 1973 : _fluid_density_node(_has_mobility ? &getMaterialProperty<std::vector<Real>>(
129 : "PorousFlow_fluid_phase_density_nodal")
130 : : nullptr),
131 1973 : _dfluid_density_node_dvar(_has_mobility ? &getMaterialProperty<std::vector<std::vector<Real>>>(
132 : "dPorousFlow_fluid_phase_density_nodal_dvar")
133 : : nullptr),
134 2160 : _fluid_viscosity(_has_mobility
135 1973 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
136 : : nullptr),
137 1973 : _dfluid_viscosity_dvar(_has_mobility ? &getMaterialProperty<std::vector<std::vector<Real>>>(
138 : "dPorousFlow_viscosity_nodal_dvar")
139 : : nullptr),
140 2013 : _relative_permeability(_has_relperm ? &getMaterialProperty<std::vector<Real>>(
141 : "PorousFlow_relative_permeability_nodal")
142 : : nullptr),
143 2160 : _drelative_permeability_dvar(_has_relperm
144 2013 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
145 : "dPorousFlow_relative_permeability_nodal_dvar")
146 : : nullptr),
147 2115 : _mass_fractions(_has_mass_fraction ? &getMaterialProperty<std::vector<std::vector<Real>>>(
148 : "PorousFlow_mass_frac_nodal")
149 : : nullptr),
150 2160 : _dmass_fractions_dvar(_has_mass_fraction
151 2115 : ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
152 : "dPorousFlow_mass_frac_nodal_dvar")
153 : : nullptr),
154 2076 : _enthalpy(_has_enthalpy ? &getMaterialPropertyByName<std::vector<Real>>(
155 : "PorousFlow_fluid_phase_enthalpy_nodal")
156 : : nullptr),
157 2076 : _denthalpy_dvar(_has_enthalpy ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
158 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
159 : : nullptr),
160 2076 : _internal_energy(_has_internal_energy ? &getMaterialPropertyByName<std::vector<Real>>(
161 : "PorousFlow_fluid_phase_internal_energy_nodal")
162 : : nullptr),
163 2160 : _dinternal_energy_dvar(_has_internal_energy
164 2076 : ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
165 : "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
166 : : nullptr),
167 1230 : _thermal_conductivity(_has_thermal_conductivity ? &getMaterialProperty<RealTensorValue>(
168 : "PorousFlow_thermal_conductivity_qp")
169 : : nullptr),
170 2160 : _dthermal_conductivity_dvar(_has_thermal_conductivity
171 1230 : ? &getMaterialProperty<std::vector<RealTensorValue>>(
172 : "dPorousFlow_thermal_conductivity_qp_dvar")
173 : : nullptr),
174 2160 : _perm_derivs(_dictator.usePermDerivs())
175 : {
176 1080 : 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 1080 : if (!_involves_fluid && (_use_mass_fraction || _use_mobility || _use_relperm || _use_enthalpy ||
185 233 : _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 1080 : 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 1080 : 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 1080 : 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 1080 : 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 1080 : 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 1080 : 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 1080 : 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 1080 : }
222 :
223 : Real
224 1450572 : PorousFlowSink::computeQpResidual()
225 : {
226 1450572 : Real flux = _test[_i][_qp] * multiplier();
227 1450572 : if (_use_mobility)
228 : {
229 : const Real k =
230 582010 : ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp]; // do not upwind permeability
231 582010 : flux *= (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
232 : }
233 1450572 : if (_use_relperm)
234 613212 : flux *= (*_relative_permeability)[_i][_ph];
235 1450572 : if (_use_mass_fraction)
236 587370 : flux *= (*_mass_fractions)[_i][_ph][_sp];
237 1450572 : if (_use_enthalpy)
238 112928 : flux *= (*_enthalpy)[_i][_ph];
239 1450572 : if (_use_internal_energy)
240 34096 : flux *= (*_internal_energy)[_i][_ph];
241 1450572 : if (_use_thermal_conductivity)
242 33792 : flux *= ((*_thermal_conductivity)[_qp] * _normals[_qp]) *
243 33792 : _normals[_qp]; // do not upwind thermal_conductivity
244 :
245 1450572 : return flux;
246 : }
247 :
248 : Real
249 5472136 : PorousFlowSink::computeQpJacobian()
250 : {
251 5472136 : return jac(_var.number());
252 : }
253 :
254 : Real
255 837276 : PorousFlowSink::computeQpOffDiagJacobian(unsigned int jvar)
256 : {
257 837276 : return jac(jvar);
258 : }
259 :
260 : Real
261 6309412 : PorousFlowSink::jac(unsigned int jvar) const
262 : {
263 6309412 : if (_dictator.notPorousFlowVariable(jvar))
264 : return 0.0;
265 6309412 : 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 6309412 : Real flux = _test[_i][_qp] * multiplier();
274 6309412 : Real deriv = _test[_i][_qp] * (_i != _j ? 0.0 : dmultiplier_dvar(pvar));
275 :
276 6309412 : if (_use_mobility)
277 : {
278 620752 : const Real k = ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp];
279 620752 : const Real mob = (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
280 :
281 : Real mobprime = 0.0;
282 620752 : if (_perm_derivs)
283 : {
284 19200 : RealTensorValue ktprime = (*_dpermeability_dvar)[_qp][pvar] * _phi[_j][_qp];
285 76800 : for (const auto i : make_range(Moose::dim))
286 115200 : ktprime += (*_dpermeability_dgradvar)[_qp][i][pvar] * _grad_phi[_j][_qp](i);
287 19200 : const Real kprime = (ktprime * _normals[_qp]) * _normals[_qp];
288 :
289 19200 : mobprime += (*_fluid_density_node)[_i][_ph] * kprime / (*_fluid_viscosity)[_i][_ph];
290 : }
291 :
292 620752 : mobprime +=
293 620752 : (_i != _j
294 620752 : ? 0.0
295 112232 : : (*_dfluid_density_node_dvar)[_i][_ph][pvar] * k / (*_fluid_viscosity)[_i][_ph] -
296 112232 : (*_fluid_density_node)[_i][_ph] * k * (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
297 : std::pow((*_fluid_viscosity)[_i][_ph], 2));
298 620752 : deriv = mob * deriv + mobprime * flux;
299 620752 : flux *= mob;
300 : }
301 6309412 : if (_use_relperm)
302 : {
303 1041504 : const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
304 1041504 : deriv = (*_relative_permeability)[_i][_ph] * deriv + relperm_prime * flux;
305 1041504 : flux *= (*_relative_permeability)[_i][_ph];
306 : }
307 6309412 : if (_use_mass_fraction)
308 : {
309 841696 : const Real mf_prime = (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
310 841696 : deriv = (*_mass_fractions)[_i][_ph][_sp] * deriv + mf_prime * flux;
311 841696 : flux *= (*_mass_fractions)[_i][_ph][_sp];
312 : }
313 6309412 : if (_use_enthalpy)
314 : {
315 73360 : const Real en_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
316 73360 : deriv = (*_enthalpy)[_i][_ph] * deriv + en_prime * flux;
317 73360 : flux *= (*_enthalpy)[_i][_ph];
318 : }
319 6309412 : if (_use_internal_energy)
320 : {
321 23552 : const Real ie_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
322 23552 : deriv = (*_internal_energy)[_i][_ph] * deriv + ie_prime * flux;
323 23552 : flux *= (*_internal_energy)[_i][_ph];
324 : }
325 6309412 : if (_use_thermal_conductivity)
326 : {
327 21504 : const Real tc = ((*_thermal_conductivity)[_qp] * _normals[_qp]) * _normals[_qp];
328 21504 : const RealTensorValue tctprime = (*_dthermal_conductivity_dvar)[_qp][pvar] * _phi[_j][_qp];
329 21504 : const Real tcprime = (tctprime * _normals[_qp]) * _normals[_qp];
330 21504 : deriv = tc * deriv + tcprime * flux;
331 : // don't need this: flux *= tc;
332 : }
333 : return deriv;
334 : }
335 :
336 : Real
337 7682740 : PorousFlowSink::multiplier() const
338 : {
339 7682740 : return _m_func.value(_t, _q_point[_qp]);
340 : }
341 :
342 : Real
343 831182 : PorousFlowSink::dmultiplier_dvar(unsigned int /*pvar*/) const
344 : {
345 831182 : return 0.0;
346 : }
|
