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 2014 : PorousFlowLineSink::validParams()
15 : {
16 2014 : InputParameters params = PorousFlowLineGeometry::validParams();
17 4028 : MooseEnum p_or_t_choice("pressure=0 temperature=1", "pressure");
18 4028 : 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 4028 : 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 4028 : params.addRequiredParam<UserObjectName>(
29 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
30 4028 : params.addParam<unsigned int>(
31 : "fluid_phase",
32 4028 : 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 4028 : 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 4028 : params.addParam<bool>(
42 4028 : "use_relative_permeability", false, "Multiply the flux by the fluid relative permeability");
43 4028 : params.addParam<bool>("use_mobility", false, "Multiply the flux by the fluid mobility");
44 4028 : params.addParam<bool>("use_enthalpy", false, "Multiply the flux by the fluid enthalpy");
45 4028 : params.addParam<bool>(
46 4028 : "use_internal_energy", false, "Multiply the flux by the fluid internal energy");
47 4028 : params.addCoupledVar("multiplying_var", 1.0, "Fluxes will be moultiplied by this variable");
48 2014 : 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 2014 : return params;
51 2014 : }
52 :
53 1131 : PorousFlowLineSink::PorousFlowLineSink(const InputParameters & parameters)
54 : : PorousFlowLineGeometry(parameters),
55 1131 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
56 1131 : _total_outflow_mass(
57 1131 : const_cast<PorousFlowSumQuantity &>(getUserObject<PorousFlowSumQuantity>("SumQuantityUO"))),
58 :
59 1131 : _has_porepressure(
60 1131 : hasMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp") &&
61 2238 : hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_porepressure_qp_dvar")),
62 1131 : _has_temperature(hasMaterialProperty<Real>("PorousFlow_temperature_qp") &&
63 2260 : hasMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
64 1131 : _has_mass_fraction(
65 1131 : hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
66 2240 : hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
67 : "dPorousFlow_mass_frac_nodal_dvar")),
68 1131 : _has_relative_permeability(
69 1131 : hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
70 2178 : hasMaterialProperty<std::vector<std::vector<Real>>>(
71 : "dPorousFlow_relative_permeability_nodal_dvar")),
72 1131 : _has_mobility(
73 2178 : hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
74 2178 : hasMaterialProperty<std::vector<std::vector<Real>>>(
75 1047 : "dPorousFlow_relative_permeability_nodal_dvar") &&
76 3225 : hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
77 2178 : hasMaterialProperty<std::vector<std::vector<Real>>>(
78 1043 : "dPorousFlow_fluid_phase_density_nodal_dvar") &&
79 3305 : hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
80 2174 : hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
81 1131 : _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
82 2234 : hasMaterialProperty<std::vector<std::vector<Real>>>(
83 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
84 1131 : _has_internal_energy(
85 1131 : hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
86 2234 : hasMaterialProperty<std::vector<std::vector<Real>>>(
87 : "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
88 :
89 2262 : _p_or_t(getParam<MooseEnum>("function_of").getEnum<PorTchoice>()),
90 2262 : _use_mass_fraction(isParamValid("mass_fraction_component")),
91 2262 : _use_relative_permeability(getParam<bool>("use_relative_permeability")),
92 2262 : _use_mobility(getParam<bool>("use_mobility")),
93 2262 : _use_enthalpy(getParam<bool>("use_enthalpy")),
94 2262 : _use_internal_energy(getParam<bool>("use_internal_energy")),
95 :
96 2262 : _ph(getParam<unsigned int>("fluid_phase")),
97 1399 : _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
98 :
99 973 : _pp((_p_or_t == PorTchoice::pressure && _has_porepressure)
100 3075 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")
101 : : nullptr),
102 973 : _dpp_dvar((_p_or_t == PorTchoice::pressure && _has_porepressure)
103 3075 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
104 : "dPorousFlow_porepressure_qp_dvar")
105 : : nullptr),
106 158 : _temperature((_p_or_t == PorTchoice::temperature && _has_temperature)
107 1445 : ? &getMaterialProperty<Real>("PorousFlow_temperature_qp")
108 : : nullptr),
109 1131 : _dtemperature_dvar(
110 158 : (_p_or_t == PorTchoice::temperature && _has_temperature)
111 1445 : ? &getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")
112 : : nullptr),
113 1131 : _fluid_density_node(
114 515 : (_use_mobility && _has_mobility)
115 2155 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")
116 : : nullptr),
117 515 : _dfluid_density_node_dvar((_use_mobility && _has_mobility)
118 2155 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
119 : "dPorousFlow_fluid_phase_density_nodal_dvar")
120 : : nullptr),
121 515 : _fluid_viscosity((_use_mobility && _has_mobility)
122 2155 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
123 : : nullptr),
124 515 : _dfluid_viscosity_dvar((_use_mobility && _has_mobility)
125 2155 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
126 : "dPorousFlow_viscosity_nodal_dvar")
127 : : nullptr),
128 1131 : _relative_permeability(
129 515 : ((_use_mobility && _has_mobility) ||
130 622 : (_use_relative_permeability && _has_relative_permeability))
131 2947 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal")
132 : : nullptr),
133 515 : _drelative_permeability_dvar(((_use_mobility && _has_mobility) ||
134 622 : (_use_relative_permeability && _has_relative_permeability))
135 2947 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
136 : "dPorousFlow_relative_permeability_nodal_dvar")
137 : : nullptr),
138 1131 : _mass_fractions(
139 268 : (_use_mass_fraction && _has_mass_fraction)
140 1665 : ? &getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")
141 : : nullptr),
142 268 : _dmass_fractions_dvar((_use_mass_fraction && _has_mass_fraction)
143 1665 : ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
144 : "dPorousFlow_mass_frac_nodal_dvar")
145 : : nullptr),
146 2234 : _enthalpy(_has_enthalpy ? &getMaterialPropertyByName<std::vector<Real>>(
147 : "PorousFlow_fluid_phase_enthalpy_nodal")
148 : : nullptr),
149 2234 : _denthalpy_dvar(_has_enthalpy ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
150 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
151 : : nullptr),
152 2234 : _internal_energy(_has_internal_energy ? &getMaterialPropertyByName<std::vector<Real>>(
153 : "PorousFlow_fluid_phase_internal_energy_nodal")
154 : : nullptr),
155 2262 : _dinternal_energy_dvar(_has_internal_energy
156 2234 : ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
157 : "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
158 : : nullptr),
159 2262 : _multiplying_var(coupledValue("multiplying_var"))
160 : {
161 : // zero the outflow mass
162 1131 : _total_outflow_mass.zero();
163 :
164 1131 : 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 1129 : 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 1127 : 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 1125 : 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 1123 : 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 1121 : 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 1119 : 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 1113 : 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 1111 : 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 1109 : const std::vector<MooseVariableFEBase *> & coupled_vars = _dictator.getCoupledMooseVars();
213 5386 : for (unsigned int i = 0; i < coupled_vars.size(); i++)
214 8554 : addMooseVariableDependency(coupled_vars[i]);
215 1109 : }
216 :
217 : void
218 62179 : 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 62179 : _total_outflow_mass.zero();
223 :
224 62179 : PorousFlowLineGeometry::addPoints();
225 62179 : }
226 :
227 : Real
228 376870 : PorousFlowLineSink::computeQpResidual()
229 : {
230 : // Get the ID we initially assigned to this point
231 376870 : const unsigned current_dirac_ptid = currentPointCachedID();
232 376870 : Real outflow = computeQpBaseOutflow(current_dirac_ptid);
233 376868 : if (outflow == 0.0)
234 : return 0.0;
235 :
236 302052 : outflow *= _multiplying_var[_qp];
237 :
238 302052 : if (_use_relative_permeability)
239 67840 : outflow *= (*_relative_permeability)[_i][_ph];
240 :
241 302052 : if (_use_mobility)
242 173164 : outflow *= (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
243 173164 : (*_fluid_viscosity)[_i][_ph];
244 :
245 302052 : if (_use_mass_fraction)
246 141600 : outflow *= (*_mass_fractions)[_i][_ph][_sp];
247 :
248 302052 : if (_use_enthalpy)
249 18342 : outflow *= (*_enthalpy)[_i][_ph];
250 :
251 302052 : if (_use_internal_energy)
252 123840 : outflow *= (*_internal_energy)[_i][_ph];
253 :
254 302052 : _total_outflow_mass.add(
255 302052 : outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
256 :
257 302052 : return outflow;
258 : }
259 :
260 : Real
261 714716 : PorousFlowLineSink::computeQpJacobian()
262 : {
263 714716 : return jac(_var.number());
264 : }
265 :
266 : Real
267 261916 : PorousFlowLineSink::computeQpOffDiagJacobian(unsigned int jvar)
268 : {
269 261916 : return jac(jvar);
270 : }
271 :
272 : Real
273 976632 : PorousFlowLineSink::jac(unsigned int jvar)
274 : {
275 976632 : if (_dictator.notPorousFlowVariable(jvar))
276 : return 0.0;
277 976632 : const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
278 :
279 : Real outflow;
280 : Real outflowp;
281 976632 : const unsigned current_dirac_ptid = currentPointCachedID();
282 976632 : computeQpBaseOutflowJacobian(jvar, current_dirac_ptid, outflow, outflowp);
283 976632 : if (outflow == 0.0 && outflowp == 0.0)
284 : return 0.0;
285 :
286 765848 : outflow *= _multiplying_var[_qp];
287 765848 : outflowp *= _multiplying_var[_qp];
288 :
289 765848 : if (_use_relative_permeability)
290 : {
291 47200 : const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
292 47200 : outflowp = (*_relative_permeability)[_i][_ph] * outflowp + relperm_prime * outflow;
293 47200 : outflow *= (*_relative_permeability)[_i][_ph];
294 : }
295 :
296 765848 : if (_use_mobility)
297 : {
298 663552 : const Real mob = (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
299 663552 : (*_fluid_viscosity)[_i][_ph];
300 : const Real mob_prime =
301 663552 : (_i != _j
302 663552 : ? 0.0
303 82944 : : (*_drelative_permeability_dvar)[_i][_ph][pvar] * (*_fluid_density_node)[_i][_ph] /
304 82944 : (*_fluid_viscosity)[_i][_ph] +
305 82944 : (*_relative_permeability)[_i][_ph] *
306 82944 : (*_dfluid_density_node_dvar)[_i][_ph][pvar] / (*_fluid_viscosity)[_i][_ph] -
307 82944 : (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] *
308 82944 : (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
309 : Utility::pow<2>((*_fluid_viscosity)[_i][_ph]));
310 663552 : outflowp = mob * outflowp + mob_prime * outflow;
311 663552 : outflow *= mob;
312 : }
313 :
314 765848 : if (_use_mass_fraction)
315 : {
316 : const Real mass_fractions_prime =
317 100800 : (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
318 100800 : outflowp = (*_mass_fractions)[_i][_ph][_sp] * outflowp + mass_fractions_prime * outflow;
319 100800 : outflow *= (*_mass_fractions)[_i][_ph][_sp];
320 : }
321 :
322 765848 : if (_use_enthalpy)
323 : {
324 19552 : const Real enthalpy_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
325 19552 : outflowp = (*_enthalpy)[_i][_ph] * outflowp + enthalpy_prime * outflow;
326 19552 : outflow *= (*_enthalpy)[_i][_ph];
327 : }
328 :
329 765848 : if (_use_internal_energy)
330 : {
331 87360 : const Real internal_energy_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
332 87360 : 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 765848 : return outflowp;
338 : }
339 :
340 : Real
341 1422078 : PorousFlowLineSink::ptqp() const
342 : {
343 1422078 : return (_p_or_t == PorTchoice::pressure ? (*_pp)[_qp][_ph] : (*_temperature)[_qp]);
344 : }
345 :
346 : Real
347 943032 : PorousFlowLineSink::dptqp(unsigned pvar) const
348 : {
349 943032 : return (_p_or_t == PorTchoice::pressure ? (*_dpp_dvar)[_qp][_ph][pvar]
350 82880 : : (*_dtemperature_dvar)[_qp][pvar]);
351 : }
|
