www.mooseframework.org
PorousFlowLineSink.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "PorousFlowLineSink.h"
9 #include "libmesh/utility.h"
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<PorousFlowLineGeometry>();
16  MooseEnum p_or_t_choice("pressure=0 temperature=1", "pressure");
17  params.addParam<MooseEnum>("function_of",
18  p_or_t_choice,
19  "Modifying functions will be a function of either pressure and "
20  "permeability (eg, for boreholes that pump fluids) or "
21  "temperature and thermal conductivity (eg, for boreholes that "
22  "pump pure heat with no fluid flow)");
23  params.addRequiredParam<UserObjectName>(
24  "SumQuantityUO",
25  "User Object of type=PorousFlowSumQuantity in which to place the total "
26  "outflow from the line sink for each time step.");
27  params.addRequiredParam<UserObjectName>(
28  "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
29  params.addParam<unsigned int>(
30  "fluid_phase",
31  0,
32  "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) "
33  "controls the flux to the line sink. For p_or_t=temperature, and without "
34  "any use_*, this parameter is irrelevant");
35  params.addParam<unsigned int>("mass_fraction_component",
36  "The index corresponding to a fluid "
37  "component. If supplied, the flux will "
38  "be multiplied by the nodal mass "
39  "fraction for the component");
40  params.addParam<bool>(
41  "use_relative_permeability", false, "Multiply the flux by the fluid relative permeability");
42  params.addParam<bool>("use_mobility", false, "Multiply the flux by the fluid mobility");
43  params.addParam<bool>("use_enthalpy", false, "Multiply the flux by the fluid enthalpy");
44  params.addParam<bool>(
45  "use_internal_energy", false, "Multiply the flux by the fluid internal energy");
46  params.addClassDescription("Approximates a line sink in the mesh by a sequence of weighted Dirac "
47  "points whose positions are read from a file");
48  return params;
49 }
50 
51 PorousFlowLineSink::PorousFlowLineSink(const InputParameters & parameters)
52  : PorousFlowLineGeometry(parameters),
53  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
54  _total_outflow_mass(
55  const_cast<PorousFlowSumQuantity &>(getUserObject<PorousFlowSumQuantity>("SumQuantityUO"))),
56 
57  _has_porepressure(
58  hasMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp") &&
59  hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_porepressure_qp_dvar")),
60  _has_temperature(hasMaterialProperty<Real>("PorousFlow_temperature_qp") &&
61  hasMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
62  _has_mass_fraction(
63  hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
64  hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
65  "dPorousFlow_mass_frac_nodal_dvar")),
66  _has_relative_permeability(
67  hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
68  hasMaterialProperty<std::vector<std::vector<Real>>>(
69  "dPorousFlow_relative_permeability_nodal_dvar")),
70  _has_mobility(
71  hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
72  hasMaterialProperty<std::vector<std::vector<Real>>>(
73  "dPorousFlow_relative_permeability_nodal_dvar") &&
74  hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
75  hasMaterialProperty<std::vector<std::vector<Real>>>(
76  "dPorousFlow_fluid_phase_density_nodal_dvar") &&
77  hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
78  hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
79  _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
80  hasMaterialProperty<std::vector<std::vector<Real>>>(
81  "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
82  _has_internal_energy(
83  hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
84  hasMaterialProperty<std::vector<std::vector<Real>>>(
85  "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
86 
87  _p_or_t(getParam<MooseEnum>("function_of").getEnum<PorTchoice>()),
88  _use_mass_fraction(isParamValid("mass_fraction_component")),
89  _use_relative_permeability(getParam<bool>("use_relative_permeability")),
90  _use_mobility(getParam<bool>("use_mobility")),
91  _use_enthalpy(getParam<bool>("use_enthalpy")),
92  _use_internal_energy(getParam<bool>("use_internal_energy")),
93 
94  _ph(getParam<unsigned int>("fluid_phase")),
95  _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
96 
97  _pp((_p_or_t == PorTchoice::pressure && _has_porepressure)
98  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")
99  : nullptr),
100  _dpp_dvar((_p_or_t == PorTchoice::pressure && _has_porepressure)
101  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
102  "dPorousFlow_porepressure_qp_dvar")
103  : nullptr),
104  _temperature((_p_or_t == PorTchoice::temperature && _has_temperature)
105  ? &getMaterialProperty<Real>("PorousFlow_temperature_qp")
106  : nullptr),
107  _dtemperature_dvar(
108  (_p_or_t == PorTchoice::temperature && _has_temperature)
109  ? &getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")
110  : nullptr),
111  _fluid_density_node(
112  (_use_mobility && _has_mobility)
113  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")
114  : nullptr),
115  _dfluid_density_node_dvar((_use_mobility && _has_mobility)
116  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
117  "dPorousFlow_fluid_phase_density_nodal_dvar")
118  : nullptr),
119  _fluid_viscosity((_use_mobility && _has_mobility)
120  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
121  : nullptr),
122  _dfluid_viscosity_dvar((_use_mobility && _has_mobility)
123  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
124  "dPorousFlow_viscosity_nodal_dvar")
125  : nullptr),
126  _relative_permeability(
127  ((_use_mobility && _has_mobility) ||
128  (_use_relative_permeability && _has_relative_permeability))
129  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal")
130  : nullptr),
131  _drelative_permeability_dvar(((_use_mobility && _has_mobility) ||
132  (_use_relative_permeability && _has_relative_permeability))
133  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
134  "dPorousFlow_relative_permeability_nodal_dvar")
135  : nullptr),
136  _mass_fractions(
137  (_use_mass_fraction && _has_mass_fraction)
138  ? &getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")
139  : nullptr),
140  _dmass_fractions_dvar((_use_mass_fraction && _has_mass_fraction)
141  ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
142  "dPorousFlow_mass_frac_nodal_dvar")
143  : nullptr),
144  _enthalpy(
145  _has_enthalpy
146  ? &getMaterialPropertyByName<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal")
147  : nullptr),
148  _denthalpy_dvar(_has_enthalpy
149  ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
150  "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
151  : nullptr),
152  _internal_energy(_has_internal_energy
153  ? &getMaterialPropertyByName<std::vector<Real>>(
154  "PorousFlow_fluid_phase_internal_energy_nodal")
155  : nullptr),
156  _dinternal_energy_dvar(_has_internal_energy
157  ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
158  "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
159  : nullptr)
160 {
161  // zero the outflow mass
163 
164  if (_ph >= _dictator.numPhases())
165  mooseError("PorousFlowLineSink: The Dictator declares that the number of fluid phases is ",
167  ", but you have set the fluid_phase to ",
168  _ph,
169  ". You must try harder.");
171  mooseError("PorousFlowLineSink: The Dictator declares that the number of fluid components is ",
173  ", but you have set the mass_fraction_component to ",
174  _sp,
175  ". Please be assured that the Dictator has noted your error.");
177  mooseError("PorousFlowLineSink: You have specified function_of=porepressure, but you do not "
178  "have a quadpoint porepressure material");
180  mooseError("PorousFlowLineSink: You have specified function_of=temperature, but you do not "
181  "have a quadpoint temperature material");
183  mooseError("PorousFlowLineSink: You have specified a fluid component, but do not have a nodal "
184  "mass-fraction material");
186  mooseError("PorousFlowLineSink: You have set use_relative_permeability=true, but do not have a "
187  "nodal relative permeability material");
189  mooseError("PorousFlowLineSink: You have set use_mobility=true, but do not have nodal density, "
190  "relative permeability or viscosity material");
192  mooseError("PorousFlowLineSink: You have set use_enthalpy=true, but do not have a nodal "
193  "enthalpy material");
195  mooseError("PorousFlowLineSink: You have set use_internal_energy=true, but do not have a nodal "
196  "internal-energy material");
197 
198  // To correctly compute the Jacobian terms,
199  // tell MOOSE that this DiracKernel depends on all the PorousFlow Variables
200  const std::vector<MooseVariable *> & coupled_vars = _dictator.getCoupledMooseVars();
201  for (unsigned int i = 0; i < coupled_vars.size(); i++)
202  addMooseVariableDependency(coupled_vars[i]);
203 }
204 
205 void
207 {
208  // This function gets called just before the DiracKernel is evaluated
209  // so this is a handy place to zero this out.
211 
213 }
214 
215 Real
217 {
218  // Get the ID we initially assigned to this point
219  const unsigned current_dirac_ptid = currentPointCachedID();
220  Real outflow = computeQpBaseOutflow(current_dirac_ptid);
221  if (outflow == 0.0)
222  return 0.0;
223 
225  outflow *= (*_relative_permeability)[_i][_ph];
226 
227  if (_use_mobility)
228  outflow *= (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
229  (*_fluid_viscosity)[_i][_ph];
230 
231  if (_use_mass_fraction)
232  outflow *= (*_mass_fractions)[_i][_ph][_sp];
233 
234  if (_use_enthalpy)
235  outflow *= (*_enthalpy)[_i][_ph];
236 
238  outflow *= (*_internal_energy)[_i][_ph];
239 
241  outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
242 
243  return outflow;
244 }
245 
246 Real
248 {
249  return jac(_var.number());
250 }
251 
252 Real
254 {
255  return jac(jvar);
256 }
257 
258 Real
259 PorousFlowLineSink::jac(unsigned int jvar)
260 {
262  return 0.0;
263  const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
264 
265  Real outflow;
266  Real outflowp;
267  const unsigned current_dirac_ptid = currentPointCachedID();
268  computeQpBaseOutflowJacobian(jvar, current_dirac_ptid, outflow, outflowp);
269  if (outflow == 0.0 && outflowp == 0.0)
270  return 0.0;
271 
273  {
274  const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
275  outflowp = (*_relative_permeability)[_i][_ph] * outflowp + relperm_prime * outflow;
276  outflow *= (*_relative_permeability)[_i][_ph];
277  }
278 
279  if (_use_mobility)
280  {
281  const Real mob = (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
282  (*_fluid_viscosity)[_i][_ph];
283  const Real mob_prime =
284  (_i != _j
285  ? 0.0
286  : (*_drelative_permeability_dvar)[_i][_ph][pvar] * (*_fluid_density_node)[_i][_ph] /
287  (*_fluid_viscosity)[_i][_ph] +
288  (*_relative_permeability)[_i][_ph] *
289  (*_dfluid_density_node_dvar)[_i][_ph][pvar] / (*_fluid_viscosity)[_i][_ph] -
290  (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] *
291  (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
292  Utility::pow<2>((*_fluid_viscosity)[_i][_ph]));
293  outflowp = mob * outflowp + mob_prime * outflow;
294  outflow *= mob;
295  }
296 
297  if (_use_mass_fraction)
298  {
299  const Real mass_fractions_prime =
300  (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
301  outflowp = (*_mass_fractions)[_i][_ph][_sp] * outflowp + mass_fractions_prime * outflow;
302  outflow *= (*_mass_fractions)[_i][_ph][_sp];
303  }
304 
305  if (_use_enthalpy)
306  {
307  const Real enthalpy_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
308  outflowp = (*_enthalpy)[_i][_ph] * outflowp + enthalpy_prime * outflow;
309  outflow *= (*_enthalpy)[_i][_ph];
310  }
311 
313  {
314  const Real internal_energy_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
315  outflowp = (*_internal_energy)[_i][_ph] * outflowp + internal_energy_prime * outflow;
316  // this multiplication was performed, but the code does not need to know: outflow *=
317  // (*_internal_energy)[_i][_ph];
318  }
319 
320  return outflowp;
321 }
322 
323 Real
325 {
326  return (_p_or_t == PorTchoice::pressure ? (*_pp)[_qp][_ph] : (*_temperature)[_qp]);
327 }
328 
329 Real
330 PorousFlowLineSink::dptqp(unsigned pvar) const
331 {
332  return (_p_or_t == PorTchoice::pressure ? (*_dpp_dvar)[_qp][_ph][pvar]
333  : (*_dtemperature_dvar)[_qp][pvar]);
334 }
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.
const PorousFlowDictator & _dictator
PorousFlow UserObject.
const bool _has_enthalpy
Whether an enthalpy material exists (for error checking)
const bool _has_temperature
Whether a quadpoint temperature material exists (for error checking)
const bool _has_relative_permeability
Whether a relative permeability material exists (for error checking)
InputParameters validParams< PorousFlowLineSink >()
Approximates a borehole by a sequence of Dirac Points.
const bool _use_relative_permeability
Whether the flux will be multiplied by the relative permeability.
InputParameters validParams< PorousFlowLineGeometry >()
virtual Real computeQpResidual() override
unsigned int numComponents() const
the number of fluid components
const std::string temperature
Definition: NS.h:25
void zero()
sets _total = 0
const bool _has_porepressure
Whether a quadpoint porepressure material exists (for error checking)
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)
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 _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...
Sums into _total This is used, for instance, to record the total mass flowing into a borehole...
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
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.
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
unsigned int numPhases() const
the number of fluid phases
virtual void addPoints() override
Add Dirac Points to the borehole.
const std::string pressure
Definition: NS.h:24
const bool _use_internal_energy
Whether the flux will be multiplied by the internal-energy.
PorousFlowLineSink(const InputParameters &parameters)
const bool _has_mass_fraction
Whether a mass_fraction material exists (for error checking)
Real dptqp(unsigned pvar) const
If _p_or_t==0, then returns d(quadpoint porepressure)/d(PorousFlow variable), else returns d(quadpoin...
const unsigned int _sp
The component number (only used if _use_mass_fraction==true)
bool notPorousFlowVariable(unsigned int moose_var_num) const
returns true if moose_var_num is not a porous flow variabe
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
the PorousFlow variable number
enum PorousFlowLineSink::PorTchoice _p_or_t
Real jac(unsigned int jvar)
Jacobian contribution for the derivative wrt the variable jvar.
const bool _use_mobility
Whether the flux will be multiplied by the mobility.
const MaterialProperty< Real > *const _temperature
Quadpoint temperature.