www.mooseframework.org
PorousFlowLineSink.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 "PorousFlowLineSink.h"
11 #include "libmesh/utility.h"
12 
15 {
17  MooseEnum p_or_t_choice("pressure=0 temperature=1", "pressure");
18  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  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  params.addRequiredParam<UserObjectName>(
29  "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
30  params.addParam<unsigned int>(
31  "fluid_phase",
32  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  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  params.addParam<bool>(
42  "use_relative_permeability", false, "Multiply the flux by the fluid relative permeability");
43  params.addParam<bool>("use_mobility", false, "Multiply the flux by the fluid mobility");
44  params.addParam<bool>("use_enthalpy", false, "Multiply the flux by the fluid enthalpy");
45  params.addParam<bool>(
46  "use_internal_energy", false, "Multiply the flux by the fluid internal energy");
47  params.addCoupledVar("multiplying_var", 1.0, "Fluxes will be moultiplied by this variable");
48  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  return params;
51 }
52 
54  : PorousFlowLineGeometry(parameters),
55  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
56  _total_outflow_mass(
57  const_cast<PorousFlowSumQuantity &>(getUserObject<PorousFlowSumQuantity>("SumQuantityUO"))),
58 
59  _has_porepressure(
60  hasMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp") &&
61  hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_porepressure_qp_dvar")),
62  _has_temperature(hasMaterialProperty<Real>("PorousFlow_temperature_qp") &&
63  hasMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
64  _has_mass_fraction(
65  hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
66  hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
67  "dPorousFlow_mass_frac_nodal_dvar")),
68  _has_relative_permeability(
69  hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
70  hasMaterialProperty<std::vector<std::vector<Real>>>(
71  "dPorousFlow_relative_permeability_nodal_dvar")),
72  _has_mobility(
73  hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
74  hasMaterialProperty<std::vector<std::vector<Real>>>(
75  "dPorousFlow_relative_permeability_nodal_dvar") &&
76  hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
77  hasMaterialProperty<std::vector<std::vector<Real>>>(
78  "dPorousFlow_fluid_phase_density_nodal_dvar") &&
79  hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
80  hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
81  _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
82  hasMaterialProperty<std::vector<std::vector<Real>>>(
83  "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
84  _has_internal_energy(
85  hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
86  hasMaterialProperty<std::vector<std::vector<Real>>>(
87  "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
88 
89  _p_or_t(getParam<MooseEnum>("function_of").getEnum<PorTchoice>()),
90  _use_mass_fraction(isParamValid("mass_fraction_component")),
91  _use_relative_permeability(getParam<bool>("use_relative_permeability")),
92  _use_mobility(getParam<bool>("use_mobility")),
93  _use_enthalpy(getParam<bool>("use_enthalpy")),
94  _use_internal_energy(getParam<bool>("use_internal_energy")),
95 
96  _ph(getParam<unsigned int>("fluid_phase")),
97  _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
98 
99  _pp((_p_or_t == PorTchoice::pressure && _has_porepressure)
100  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")
101  : nullptr),
102  _dpp_dvar((_p_or_t == PorTchoice::pressure && _has_porepressure)
103  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
104  "dPorousFlow_porepressure_qp_dvar")
105  : nullptr),
106  _temperature((_p_or_t == PorTchoice::temperature && _has_temperature)
107  ? &getMaterialProperty<Real>("PorousFlow_temperature_qp")
108  : nullptr),
109  _dtemperature_dvar(
110  (_p_or_t == PorTchoice::temperature && _has_temperature)
111  ? &getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")
112  : nullptr),
113  _fluid_density_node(
114  (_use_mobility && _has_mobility)
115  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")
116  : nullptr),
117  _dfluid_density_node_dvar((_use_mobility && _has_mobility)
118  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
119  "dPorousFlow_fluid_phase_density_nodal_dvar")
120  : nullptr),
121  _fluid_viscosity((_use_mobility && _has_mobility)
122  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
123  : nullptr),
124  _dfluid_viscosity_dvar((_use_mobility && _has_mobility)
125  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
126  "dPorousFlow_viscosity_nodal_dvar")
127  : nullptr),
128  _relative_permeability(
129  ((_use_mobility && _has_mobility) ||
130  (_use_relative_permeability && _has_relative_permeability))
131  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal")
132  : nullptr),
133  _drelative_permeability_dvar(((_use_mobility && _has_mobility) ||
134  (_use_relative_permeability && _has_relative_permeability))
135  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
136  "dPorousFlow_relative_permeability_nodal_dvar")
137  : nullptr),
138  _mass_fractions(
139  (_use_mass_fraction && _has_mass_fraction)
140  ? &getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")
141  : nullptr),
142  _dmass_fractions_dvar((_use_mass_fraction && _has_mass_fraction)
143  ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
144  "dPorousFlow_mass_frac_nodal_dvar")
145  : nullptr),
146  _enthalpy(_has_enthalpy ? &getMaterialPropertyByName<std::vector<Real>>(
147  "PorousFlow_fluid_phase_enthalpy_nodal")
148  : nullptr),
149  _denthalpy_dvar(_has_enthalpy ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
150  "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
151  : nullptr),
152  _internal_energy(_has_internal_energy ? &getMaterialPropertyByName<std::vector<Real>>(
153  "PorousFlow_fluid_phase_internal_energy_nodal")
154  : nullptr),
155  _dinternal_energy_dvar(_has_internal_energy
156  ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
157  "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
158  : nullptr),
159  _multiplying_var(coupledValue("multiplying_var"))
160 {
161  // zero the outflow mass
163 
164  if (_ph >= _dictator.numPhases())
165  paramError("fluid_phase",
166  "The Dictator proclaims that the maximum phase index in this simulation is ",
167  _dictator.numPhases() - 1,
168  " whereas you have used ",
169  _ph,
170  ". Remember that indexing starts at 0. You must try harder.");
171 
173  paramError(
174  "mass_fraction_component",
175  "The Dictator proclaims that the maximum fluid component index in this simulation is ",
176  _dictator.numComponents() - 1,
177  " whereas you have used ",
178  _sp,
179  ". Remember that indexing starts at 0. Please be assured that the Dictator has noted your "
180  "error.");
181 
183  mooseError("PorousFlowLineSink: You have specified function_of=porepressure, but you do not "
184  "have a quadpoint porepressure material");
185 
187  mooseError("PorousFlowLineSink: You have specified function_of=temperature, but you do not "
188  "have a quadpoint temperature material");
189 
191  mooseError("PorousFlowLineSink: You have specified a fluid component, but do not have a nodal "
192  "mass-fraction material");
193 
195  mooseError("PorousFlowLineSink: You have set use_relative_permeability=true, but do not have a "
196  "nodal relative permeability material");
197 
199  mooseError("PorousFlowLineSink: You have set use_mobility=true, but do not have nodal density, "
200  "relative permeability or viscosity material");
201 
203  mooseError("PorousFlowLineSink: You have set use_enthalpy=true, but do not have a nodal "
204  "enthalpy material");
205 
207  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  const std::vector<MooseVariableFEBase *> & coupled_vars = _dictator.getCoupledMooseVars();
213  for (unsigned int i = 0; i < coupled_vars.size(); i++)
214  addMooseVariableDependency(coupled_vars[i]);
215 }
216 
217 void
219 {
220  // This function gets called just before the DiracKernel is evaluated
221  // so this is a handy place to zero this out.
223 
225 }
226 
227 Real
229 {
230  // Get the ID we initially assigned to this point
231  const unsigned current_dirac_ptid = currentPointCachedID();
232  Real outflow = computeQpBaseOutflow(current_dirac_ptid);
233  if (outflow == 0.0)
234  return 0.0;
235 
236  outflow *= _multiplying_var[_qp];
237 
239  outflow *= (*_relative_permeability)[_i][_ph];
240 
241  if (_use_mobility)
242  outflow *= (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
243  (*_fluid_viscosity)[_i][_ph];
244 
245  if (_use_mass_fraction)
246  outflow *= (*_mass_fractions)[_i][_ph][_sp];
247 
248  if (_use_enthalpy)
249  outflow *= (*_enthalpy)[_i][_ph];
250 
252  outflow *= (*_internal_energy)[_i][_ph];
253 
255  outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
256 
257  return outflow;
258 }
259 
260 Real
262 {
263  return jac(_var.number());
264 }
265 
266 Real
268 {
269  return jac(jvar);
270 }
271 
272 Real
273 PorousFlowLineSink::jac(unsigned int jvar)
274 {
276  return 0.0;
277  const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
278 
279  Real outflow;
280  Real outflowp;
281  const unsigned current_dirac_ptid = currentPointCachedID();
282  computeQpBaseOutflowJacobian(jvar, current_dirac_ptid, outflow, outflowp);
283  if (outflow == 0.0 && outflowp == 0.0)
284  return 0.0;
285 
286  outflow *= _multiplying_var[_qp];
287  outflowp *= _multiplying_var[_qp];
288 
290  {
291  const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
292  outflowp = (*_relative_permeability)[_i][_ph] * outflowp + relperm_prime * outflow;
293  outflow *= (*_relative_permeability)[_i][_ph];
294  }
295 
296  if (_use_mobility)
297  {
298  const Real mob = (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
299  (*_fluid_viscosity)[_i][_ph];
300  const Real mob_prime =
301  (_i != _j
302  ? 0.0
303  : (*_drelative_permeability_dvar)[_i][_ph][pvar] * (*_fluid_density_node)[_i][_ph] /
304  (*_fluid_viscosity)[_i][_ph] +
305  (*_relative_permeability)[_i][_ph] *
306  (*_dfluid_density_node_dvar)[_i][_ph][pvar] / (*_fluid_viscosity)[_i][_ph] -
307  (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] *
308  (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
309  Utility::pow<2>((*_fluid_viscosity)[_i][_ph]));
310  outflowp = mob * outflowp + mob_prime * outflow;
311  outflow *= mob;
312  }
313 
314  if (_use_mass_fraction)
315  {
316  const Real mass_fractions_prime =
317  (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
318  outflowp = (*_mass_fractions)[_i][_ph][_sp] * outflowp + mass_fractions_prime * outflow;
319  outflow *= (*_mass_fractions)[_i][_ph][_sp];
320  }
321 
322  if (_use_enthalpy)
323  {
324  const Real enthalpy_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
325  outflowp = (*_enthalpy)[_i][_ph] * outflowp + enthalpy_prime * outflow;
326  outflow *= (*_enthalpy)[_i][_ph];
327  }
328 
330  {
331  const Real internal_energy_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
332  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  return outflowp;
338 }
339 
340 Real
342 {
343  return (_p_or_t == PorTchoice::pressure ? (*_pp)[_qp][_ph] : (*_temperature)[_qp]);
344 }
345 
346 Real
347 PorousFlowLineSink::dptqp(unsigned pvar) const
348 {
349  return (_p_or_t == PorTchoice::pressure ? (*_dpp_dvar)[_qp][_ph][pvar]
350  : (*_dtemperature_dvar)[_qp][pvar]);
351 }
virtual void addPoints() override
Add Dirac Points to the line sink.
const bool _has_internal_energy
Whether an internal-energy material exists (for error checking)
void add(Real contrib)
Adds contrib to _total.
const bool _use_enthalpy
Whether the flux will be multiplied by the enthalpy.
const unsigned int _ph
The phase number.
bool notPorousFlowVariable(unsigned int moose_var_num) const
Returns true if moose_var_num is not a porous flow variabe.
const PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
const bool _has_enthalpy
Whether an enthalpy material exists (for error checking)
static InputParameters validParams()
const bool _has_temperature
Whether a quadpoint temperature material exists (for error checking)
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const bool _has_relative_permeability
Whether a relative permeability material exists (for error checking)
unsigned int number() const
Approximates a borehole by a sequence of Dirac Points.
unsigned int numComponents() const
The number of fluid components.
const bool _use_relative_permeability
Whether the flux will be multiplied by the relative permeability.
unsigned currentPointCachedID()
virtual Real computeQpResidual() override
static const std::string temperature
Definition: NS.h:57
void zero()
Sets _total = 0.
const bool _has_porepressure
Whether a quadpoint porepressure material exists (for error checking)
void addRequiredParam(const std::string &name, const std::string &doc_string)
virtual Real computeQpJacobian() override
const MaterialProperty< std::vector< Real > > *const _dtemperature_dvar
d(quadpoint temperature)/d(PorousFlow variable)
const MaterialProperty< std::vector< std::vector< Real > > > *const _dpp_dvar
d(quadpoint pore pressure in each phase)/d(PorousFlow variable)
const bool _has_mobility
Whether enough materials exist to form the mobility (for error checking)
virtual Real computeQpBaseOutflow(unsigned current_dirac_ptid) const =0
Returns the flux from the line sink (before modification by mobility, etc). Derived classes should ov...
Real ptqp() const
If _p_or_t==0, then returns the quadpoint porepressure, else returns the quadpoint temperature...
MooseVariableField< T > & _var
void paramError(const std::string &param, Args... args) const
unsigned int numPhases() const
The number of fluid phases.
Sums into _total This is used, for instance, to record the total mass flowing into a borehole...
static InputParameters validParams()
Creates a new PorousFlowLineGeometry This reads the file containing the lines of the form weight x y ...
const MaterialProperty< std::vector< Real > > *const _pp
Quadpoint pore pressure in each phase.
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
PorTchoice
whether the flux is a function of pressure or temperature
void addCoupledVar(const std::string &name, const std::string &doc_string)
const std::vector< MooseVariableFieldBase *> & getCoupledMooseVars() const
PorousFlowSumQuantity & _total_outflow_mass
This is used to hold the total fluid flowing into the line sink for each time step.
const bool _use_mass_fraction
Whether the flux will be multiplied by the mass fraction.
void addMooseVariableDependency(MooseVariableFieldBase *var)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int _j
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
unsigned int _qp
virtual void addPoints() override
Add Dirac Points to the borehole.
static const std::string pressure
Definition: NS.h:56
const bool _use_internal_energy
Whether the flux will be multiplied by the internal-energy.
void mooseError(Args &&... args) const
PorousFlowLineSink(const InputParameters &parameters)
void addClassDescription(const std::string &doc_string)
const bool _has_mass_fraction
Whether a mass_fraction material exists (for error checking)
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
The PorousFlow variable number.
const VariableValue & _multiplying_var
mass flux is multiplied by this variable evaluated at quadpoints
const unsigned int _sp
The component number (only used if _use_mass_fraction==true)
unsigned int _i
enum PorousFlowLineSink::PorTchoice _p_or_t
Real jac(unsigned int jvar)
Jacobian contribution for the derivative wrt the variable jvar.
virtual void computeQpBaseOutflowJacobian(unsigned jvar, unsigned current_dirac_ptid, Real &outflow, Real &outflowp) const =0
Calculates the BaseOutflow as well as its derivative wrt jvar. Derived classes should override this...
const bool _use_mobility
Whether the flux will be multiplied by the mobility.
void ErrorVector unsigned int
const MaterialProperty< Real > *const _temperature
Quadpoint temperature.
Real dptqp(unsigned pvar) const
If _p_or_t==0, then returns d(quadpoint porepressure)/d(PorousFlow variable), else returns d(quadpoin...