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 "PorousFlowLineSink.h"
11 : #include "libmesh/utility.h"
12 :
13 : InputParameters
14 1159 : PorousFlowLineSink::validParams()
15 : {
16 1159 : InputParameters params = PorousFlowLineGeometry::validParams();
17 2318 : MooseEnum p_or_t_choice("pressure=0 temperature=1", "pressure");
18 2318 : params.addParam<MooseEnum>("function_of",
19 : p_or_t_choice,
20 : "Modifying functions will be a function of either pressure and "
21 : "permeability (eg, for boreholes that pump fluids) or "
22 : "temperature and thermal conductivity (eg, for boreholes that "
23 : "pump pure heat with no fluid flow)");
24 2318 : params.addRequiredParam<UserObjectName>(
25 : "SumQuantityUO",
26 : "User Object of type=PorousFlowSumQuantity in which to place the total "
27 : "outflow from the line sink for each time step.");
28 2318 : params.addRequiredParam<UserObjectName>(
29 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
30 2318 : params.addParam<unsigned int>(
31 : "fluid_phase",
32 2318 : 0,
33 : "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) "
34 : "controls the flux to the line sink. For p_or_t=temperature, and without "
35 : "any use_*, this parameter is irrelevant");
36 2318 : params.addParam<unsigned int>("mass_fraction_component",
37 : "The index corresponding to a fluid "
38 : "component. If supplied, the flux will "
39 : "be multiplied by the nodal mass "
40 : "fraction for the component");
41 2318 : params.addParam<bool>(
42 2318 : "use_relative_permeability", false, "Multiply the flux by the fluid relative permeability");
43 2318 : params.addParam<bool>("use_mobility", false, "Multiply the flux by the fluid mobility");
44 2318 : params.addParam<bool>("use_enthalpy", false, "Multiply the flux by the fluid enthalpy");
45 2318 : params.addParam<bool>(
46 2318 : "use_internal_energy", false, "Multiply the flux by the fluid internal energy");
47 2318 : params.addCoupledVar("multiplying_var", 1.0, "Fluxes will be moultiplied by this variable");
48 1159 : params.addClassDescription("Approximates a line sink in the mesh by a sequence of weighted Dirac "
49 : "points whose positions are read from a file");
50 1159 : return params;
51 1159 : }
52 :
53 633 : PorousFlowLineSink::PorousFlowLineSink(const InputParameters & parameters)
54 : : PorousFlowLineGeometry(parameters),
55 633 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
56 633 : _total_outflow_mass(
57 633 : const_cast<PorousFlowSumQuantity &>(getUserObject<PorousFlowSumQuantity>("SumQuantityUO"))),
58 :
59 633 : _has_porepressure(
60 633 : hasMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp") &&
61 1254 : hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_porepressure_qp_dvar")),
62 633 : _has_temperature(hasMaterialProperty<Real>("PorousFlow_temperature_qp") &&
63 1264 : hasMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
64 633 : _has_mass_fraction(
65 633 : hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
66 1256 : hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
67 : "dPorousFlow_mass_frac_nodal_dvar")),
68 633 : _has_relative_permeability(
69 633 : hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
70 1210 : hasMaterialProperty<std::vector<std::vector<Real>>>(
71 : "dPorousFlow_relative_permeability_nodal_dvar")),
72 633 : _has_mobility(
73 1210 : hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
74 1210 : hasMaterialProperty<std::vector<std::vector<Real>>>(
75 577 : "dPorousFlow_relative_permeability_nodal_dvar") &&
76 1787 : hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
77 1210 : hasMaterialProperty<std::vector<std::vector<Real>>>(
78 573 : "dPorousFlow_fluid_phase_density_nodal_dvar") &&
79 1839 : hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
80 1206 : hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
81 633 : _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
82 1250 : hasMaterialProperty<std::vector<std::vector<Real>>>(
83 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
84 633 : _has_internal_energy(
85 633 : hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
86 1250 : hasMaterialProperty<std::vector<std::vector<Real>>>(
87 : "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
88 :
89 1266 : _p_or_t(getParam<MooseEnum>("function_of").getEnum<PorTchoice>()),
90 1266 : _use_mass_fraction(isParamValid("mass_fraction_component")),
91 1266 : _use_relative_permeability(getParam<bool>("use_relative_permeability")),
92 1266 : _use_mobility(getParam<bool>("use_mobility")),
93 1266 : _use_enthalpy(getParam<bool>("use_enthalpy")),
94 1266 : _use_internal_energy(getParam<bool>("use_internal_energy")),
95 :
96 1266 : _ph(getParam<unsigned int>("fluid_phase")),
97 805 : _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
98 :
99 535 : _pp((_p_or_t == PorTchoice::pressure && _has_porepressure)
100 1701 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")
101 : : nullptr),
102 535 : _dpp_dvar((_p_or_t == PorTchoice::pressure && _has_porepressure)
103 1701 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
104 : "dPorousFlow_porepressure_qp_dvar")
105 : : nullptr),
106 98 : _temperature((_p_or_t == PorTchoice::temperature && _has_temperature)
107 827 : ? &getMaterialProperty<Real>("PorousFlow_temperature_qp")
108 : : nullptr),
109 633 : _dtemperature_dvar(
110 98 : (_p_or_t == PorTchoice::temperature && _has_temperature)
111 827 : ? &getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")
112 : : nullptr),
113 633 : _fluid_density_node(
114 271 : (_use_mobility && _has_mobility)
115 1169 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")
116 : : nullptr),
117 271 : _dfluid_density_node_dvar((_use_mobility && _has_mobility)
118 1169 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
119 : "dPorousFlow_fluid_phase_density_nodal_dvar")
120 : : nullptr),
121 271 : _fluid_viscosity((_use_mobility && _has_mobility)
122 1169 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
123 : : nullptr),
124 271 : _dfluid_viscosity_dvar((_use_mobility && _has_mobility)
125 1169 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
126 : "dPorousFlow_viscosity_nodal_dvar")
127 : : nullptr),
128 633 : _relative_permeability(
129 271 : ((_use_mobility && _has_mobility) ||
130 368 : (_use_relative_permeability && _has_relative_permeability))
131 1613 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal")
132 : : nullptr),
133 271 : _drelative_permeability_dvar(((_use_mobility && _has_mobility) ||
134 368 : (_use_relative_permeability && _has_relative_permeability))
135 1613 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
136 : "dPorousFlow_relative_permeability_nodal_dvar")
137 : : nullptr),
138 633 : _mass_fractions(
139 172 : (_use_mass_fraction && _has_mass_fraction)
140 975 : ? &getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")
141 : : nullptr),
142 172 : _dmass_fractions_dvar((_use_mass_fraction && _has_mass_fraction)
143 975 : ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
144 : "dPorousFlow_mass_frac_nodal_dvar")
145 : : nullptr),
146 1250 : _enthalpy(_has_enthalpy ? &getMaterialPropertyByName<std::vector<Real>>(
147 : "PorousFlow_fluid_phase_enthalpy_nodal")
148 : : nullptr),
149 1250 : _denthalpy_dvar(_has_enthalpy ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
150 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
151 : : nullptr),
152 1250 : _internal_energy(_has_internal_energy ? &getMaterialPropertyByName<std::vector<Real>>(
153 : "PorousFlow_fluid_phase_internal_energy_nodal")
154 : : nullptr),
155 1266 : _dinternal_energy_dvar(_has_internal_energy
156 1250 : ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
157 : "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
158 : : nullptr),
159 1266 : _multiplying_var(coupledValue("multiplying_var"))
160 : {
161 : // zero the outflow mass
162 633 : _total_outflow_mass.zero();
163 :
164 633 : if (_ph >= _dictator.numPhases())
165 2 : paramError("fluid_phase",
166 : "The Dictator proclaims that the maximum phase index in this simulation is ",
167 2 : _dictator.numPhases() - 1,
168 : " whereas you have used ",
169 2 : _ph,
170 : ". Remember that indexing starts at 0. You must try harder.");
171 :
172 631 : if (_use_mass_fraction && _sp >= _dictator.numComponents())
173 2 : paramError(
174 : "mass_fraction_component",
175 : "The Dictator proclaims that the maximum fluid component index in this simulation is ",
176 2 : _dictator.numComponents() - 1,
177 : " whereas you have used ",
178 2 : _sp,
179 : ". Remember that indexing starts at 0. Please be assured that the Dictator has noted your "
180 : "error.");
181 :
182 629 : if (_p_or_t == PorTchoice::pressure && !_has_porepressure)
183 2 : mooseError("PorousFlowLineSink: You have specified function_of=porepressure, but you do not "
184 : "have a quadpoint porepressure material");
185 :
186 627 : if (_p_or_t == PorTchoice::temperature && !_has_temperature)
187 2 : mooseError("PorousFlowLineSink: You have specified function_of=temperature, but you do not "
188 : "have a quadpoint temperature material");
189 :
190 625 : if (_use_mass_fraction && !_has_mass_fraction)
191 2 : mooseError("PorousFlowLineSink: You have specified a fluid component, but do not have a nodal "
192 : "mass-fraction material");
193 :
194 623 : if (_use_relative_permeability && !_has_relative_permeability)
195 2 : mooseError("PorousFlowLineSink: You have set use_relative_permeability=true, but do not have a "
196 : "nodal relative permeability material");
197 :
198 621 : if (_use_mobility && !_has_mobility)
199 6 : mooseError("PorousFlowLineSink: You have set use_mobility=true, but do not have nodal density, "
200 : "relative permeability or viscosity material");
201 :
202 615 : if (_use_enthalpy && !_has_enthalpy)
203 2 : mooseError("PorousFlowLineSink: You have set use_enthalpy=true, but do not have a nodal "
204 : "enthalpy material");
205 :
206 613 : if (_use_internal_energy && !_has_internal_energy)
207 2 : mooseError("PorousFlowLineSink: You have set use_internal_energy=true, but do not have a nodal "
208 : "internal-energy material");
209 :
210 : // To correctly compute the Jacobian terms,
211 : // tell MOOSE that this DiracKernel depends on all the PorousFlow Variables
212 611 : const std::vector<MooseVariableFEBase *> & coupled_vars = _dictator.getCoupledMooseVars();
213 3202 : for (unsigned int i = 0; i < coupled_vars.size(); i++)
214 5182 : addMooseVariableDependency(coupled_vars[i]);
215 611 : }
216 :
217 : void
218 37031 : PorousFlowLineSink::addPoints()
219 : {
220 : // This function gets called just before the DiracKernel is evaluated
221 : // so this is a handy place to zero this out.
222 37031 : _total_outflow_mass.zero();
223 :
224 37031 : PorousFlowLineGeometry::addPoints();
225 37031 : }
226 :
227 : Real
228 285888 : PorousFlowLineSink::computeQpResidual()
229 : {
230 : // Get the ID we initially assigned to this point
231 285888 : const unsigned current_dirac_ptid = currentPointCachedID();
232 285888 : Real outflow = computeQpBaseOutflow(current_dirac_ptid);
233 285886 : if (outflow == 0.0)
234 : return 0.0;
235 :
236 227974 : outflow *= _multiplying_var[_qp];
237 :
238 227974 : if (_use_relative_permeability)
239 54080 : outflow *= (*_relative_permeability)[_i][_ph];
240 :
241 227974 : if (_use_mobility)
242 125078 : outflow *= (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
243 125078 : (*_fluid_viscosity)[_i][_ph];
244 :
245 227974 : if (_use_mass_fraction)
246 113280 : outflow *= (*_mass_fractions)[_i][_ph][_sp];
247 :
248 227974 : if (_use_enthalpy)
249 14598 : outflow *= (*_enthalpy)[_i][_ph];
250 :
251 227974 : if (_use_internal_energy)
252 99072 : outflow *= (*_internal_energy)[_i][_ph];
253 :
254 227974 : _total_outflow_mass.add(
255 227974 : outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
256 :
257 227974 : return outflow;
258 : }
259 :
260 : Real
261 484954 : PorousFlowLineSink::computeQpJacobian()
262 : {
263 484954 : return jac(_var.number());
264 : }
265 :
266 : Real
267 197242 : PorousFlowLineSink::computeQpOffDiagJacobian(unsigned int jvar)
268 : {
269 197242 : return jac(jvar);
270 : }
271 :
272 : Real
273 682196 : PorousFlowLineSink::jac(unsigned int jvar)
274 : {
275 682196 : if (_dictator.notPorousFlowVariable(jvar))
276 : return 0.0;
277 682196 : const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
278 :
279 : Real outflow;
280 : Real outflowp;
281 682196 : const unsigned current_dirac_ptid = currentPointCachedID();
282 682196 : computeQpBaseOutflowJacobian(jvar, current_dirac_ptid, outflow, outflowp);
283 682196 : if (outflow == 0.0 && outflowp == 0.0)
284 : return 0.0;
285 :
286 533368 : outflow *= _multiplying_var[_qp];
287 533368 : outflowp *= _multiplying_var[_qp];
288 :
289 533368 : if (_use_relative_permeability)
290 : {
291 37184 : const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
292 37184 : outflowp = (*_relative_permeability)[_i][_ph] * outflowp + relperm_prime * outflow;
293 37184 : outflow *= (*_relative_permeability)[_i][_ph];
294 : }
295 :
296 533368 : if (_use_mobility)
297 : {
298 452416 : const Real mob = (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
299 452416 : (*_fluid_viscosity)[_i][_ph];
300 : const Real mob_prime =
301 452416 : (_i != _j
302 452416 : ? 0.0
303 56552 : : (*_drelative_permeability_dvar)[_i][_ph][pvar] * (*_fluid_density_node)[_i][_ph] /
304 56552 : (*_fluid_viscosity)[_i][_ph] +
305 56552 : (*_relative_permeability)[_i][_ph] *
306 56552 : (*_dfluid_density_node_dvar)[_i][_ph][pvar] / (*_fluid_viscosity)[_i][_ph] -
307 56552 : (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] *
308 56552 : (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
309 : Utility::pow<2>((*_fluid_viscosity)[_i][_ph]));
310 452416 : outflowp = mob * outflowp + mob_prime * outflow;
311 452416 : outflow *= mob;
312 : }
313 :
314 533368 : if (_use_mass_fraction)
315 : {
316 : const Real mass_fractions_prime =
317 80640 : (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
318 80640 : outflowp = (*_mass_fractions)[_i][_ph][_sp] * outflowp + mass_fractions_prime * outflow;
319 80640 : outflow *= (*_mass_fractions)[_i][_ph][_sp];
320 : }
321 :
322 533368 : if (_use_enthalpy)
323 : {
324 14848 : const Real enthalpy_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
325 14848 : outflowp = (*_enthalpy)[_i][_ph] * outflowp + enthalpy_prime * outflow;
326 14848 : outflow *= (*_enthalpy)[_i][_ph];
327 : }
328 :
329 533368 : if (_use_internal_energy)
330 : {
331 69888 : const Real internal_energy_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
332 69888 : outflowp = (*_internal_energy)[_i][_ph] * outflowp + internal_energy_prime * outflow;
333 : // this multiplication was performed, but the code does not need to know: outflow *=
334 : // (*_internal_energy)[_i][_ph];
335 : }
336 :
337 533368 : return outflowp;
338 : }
339 :
340 : Real
341 1020788 : PorousFlowLineSink::ptqp() const
342 : {
343 1020788 : return (_p_or_t == PorTchoice::pressure ? (*_pp)[_qp][_ph] : (*_temperature)[_qp]);
344 : }
345 :
346 : Real
347 655316 : PorousFlowLineSink::dptqp(unsigned pvar) const
348 : {
349 655316 : return (_p_or_t == PorTchoice::pressure ? (*_dpp_dvar)[_qp][_ph][pvar]
350 66304 : : (*_dtemperature_dvar)[_qp][pvar]);
351 : }
|
