www.mooseframework.org
PorousFlowDarcyBase.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 "PorousFlowDarcyBase.h"
11 
12 #include "Assembly.h"
13 #include "MooseMesh.h"
14 #include "MooseVariable.h"
15 #include "SystemBase.h"
16 
17 #include "libmesh/quadrature.h"
18 
21 {
23  params.addRequiredParam<RealVectorValue>("gravity",
24  "Gravitational acceleration vector downwards (m/s^2)");
25  params.addRequiredParam<UserObjectName>(
26  "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
27  params.addParam<unsigned>("full_upwind_threshold",
28  5,
29  "If, for each timestep, the number of "
30  "upwind-downwind swaps in an element is less than "
31  "this quantity, then full upwinding is used for that element. "
32  "Otherwise the fallback scheme is employed.");
33  MooseEnum fallback_enum("quick harmonic", "quick");
34  params.addParam<MooseEnum>("fallback_scheme",
35  fallback_enum,
36  "quick: use nodal mobility without "
37  "preserving mass. harmonic: use a "
38  "harmonic mean of nodal mobilities "
39  "and preserve fluid mass");
40  params.addClassDescription("Fully-upwinded advective Darcy flux");
41  return params;
42 }
43 
45  : Kernel(parameters),
46  _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
47  _dpermeability_dvar(
48  getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
49  _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
50  "dPorousFlow_permeability_qp_dgradvar")),
51  _fluid_density_node(
52  getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
53  _dfluid_density_node_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
54  "dPorousFlow_fluid_phase_density_nodal_dvar")),
55  _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
56  _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
57  "dPorousFlow_fluid_phase_density_qp_dvar")),
58  _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")),
59  _dfluid_viscosity_dvar(
60  getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
61  _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_nodal")),
62  _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
63  _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
64  "dPorousFlow_grad_porepressure_qp_dgradvar")),
65  _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
66  "dPorousFlow_grad_porepressure_qp_dvar")),
67  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
68  _num_phases(_dictator.numPhases()),
69  _gravity(getParam<RealVectorValue>("gravity")),
70  _perm_derivs(_dictator.usePermDerivs()),
71  _full_upwind_threshold(getParam<unsigned>("full_upwind_threshold")),
72  _fallback_scheme(getParam<MooseEnum>("fallback_scheme").getEnum<FallbackEnum>()),
73  _proto_flux(_num_phases),
74  _jacobian(_num_phases),
75  _num_upwinds(),
76  _num_downwinds()
77 {
78 #ifdef LIBMESH_HAVE_TBB_API
79  if (libMesh::n_threads() > 1)
80  mooseWarning("PorousFlowDarcyBase: num_upwinds and num_downwinds may not be computed "
81  "accurately when using TBB and greater than 1 thread");
82 #endif
83 }
84 
85 void
87 {
89  _num_upwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
90  _num_downwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
91 }
92 
93 Real
94 PorousFlowDarcyBase::darcyQp(unsigned int ph) const
95 {
96  return _grad_test[_i][_qp] *
98 }
99 
100 Real
101 PorousFlowDarcyBase::darcyQpJacobian(unsigned int jvar, unsigned int ph) const
102 {
104  return 0.0;
105 
106  const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
107 
109  _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
110  _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
111 
112  deriv += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
113 
114  if (_perm_derivs)
115  {
116  deriv += _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
117  (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
118  for (const auto i : make_range(Moose::dim))
119  deriv += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
120  (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
121  }
122 
123  return _grad_test[_i][_qp] * deriv;
124 }
125 
126 Real
128 {
129  mooseError("PorousFlowDarcyBase: computeQpResidual called");
130  return 0.0;
131 }
132 
133 void
135 {
137 }
138 
139 void
141 {
143 }
144 
145 void
147 {
149 }
150 
151 void
153 {
154  if ((res_or_jac == JacRes::CALCULATE_JACOBIAN) && _dictator.notPorousFlowVariable(jvar))
155  return;
156 
157  // The PorousFlow variable index corresponding to the variable number jvar
158  const unsigned int pvar =
159  ((res_or_jac == JacRes::CALCULATE_JACOBIAN) ? _dictator.porousFlowVariableNum(jvar) : 0);
160 
162  if ((_local_ke.n() == 0) && (res_or_jac == JacRes::CALCULATE_JACOBIAN)) // this removes a problem
163  // encountered in the
164  // initial timestep when
165  // use_displaced_mesh=true
166  return;
167 
168  // The number of nodes in the element
169  const unsigned int num_nodes = _test.size();
170 
171  // Compute the residual and jacobian without the mobility terms. Even if we are computing the
172  // Jacobian we still need this in order to see which nodes are upwind and which are downwind.
173  for (unsigned ph = 0; ph < _num_phases; ++ph)
174  {
175  _proto_flux[ph].assign(num_nodes, 0);
176  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
177  {
178  for (_i = 0; _i < num_nodes; ++_i)
179  _proto_flux[ph][_i] += _JxW[_qp] * _coord[_qp] * darcyQp(ph);
180  }
181  }
182 
183  // for this element, record whether each node is "upwind" or "downwind" (or neither)
184  const unsigned elem = _current_elem->id();
185  if (_num_upwinds.find(elem) == _num_upwinds.end())
186  {
187  _num_upwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
188  _num_downwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
189  for (unsigned ph = 0; ph < _num_phases; ++ph)
190  {
191  _num_upwinds[elem][ph].assign(num_nodes, 0);
192  _num_downwinds[elem][ph].assign(num_nodes, 0);
193  }
194  }
195  // record the information once per nonlinear iteration
196  if (res_or_jac == JacRes::CALCULATE_JACOBIAN && jvar == _var.number())
197  {
198  for (unsigned ph = 0; ph < _num_phases; ++ph)
199  {
200  for (unsigned nod = 0; nod < num_nodes; ++nod)
201  {
202  if (_proto_flux[ph][nod] > 0)
203  _num_upwinds[elem][ph][nod]++;
204  else if (_proto_flux[ph][nod] < 0)
205  _num_downwinds[elem][ph][nod]++;
206  }
207  }
208  }
209 
210  // based on _num_upwinds and _num_downwinds, calculate the maximum number
211  // of upwind-downwind swaps that have been encountered in this timestep
212  // for this element
213  std::vector<unsigned> max_swaps(_num_phases, 0);
214  for (unsigned ph = 0; ph < _num_phases; ++ph)
215  {
216  for (unsigned nod = 0; nod < num_nodes; ++nod)
217  max_swaps[ph] = std::max(
218  max_swaps[ph], std::min(_num_upwinds[elem][ph][nod], _num_downwinds[elem][ph][nod]));
219  }
220 
221  // size the _jacobian correctly and calculate it for the case residual = _proto_flux
222  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
223  {
224  for (unsigned ph = 0; ph < _num_phases; ++ph)
225  {
226  _jacobian[ph].resize(_local_ke.m());
227  for (_i = 0; _i < _test.size(); _i++)
228  {
229  _jacobian[ph][_i].assign(_local_ke.n(), 0.0);
230  for (_j = 0; _j < _phi.size(); _j++)
231  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
232  _jacobian[ph][_i][_j] += _JxW[_qp] * _coord[_qp] * darcyQpJacobian(jvar, ph);
233  }
234  }
235  }
236 
237  // Loop over all the phases, computing the mass flux, which
238  // gets placed into _proto_flux, and the derivative of this
239  // which gets placed into _jacobian
240  for (unsigned int ph = 0; ph < _num_phases; ++ph)
241  {
242  if (max_swaps[ph] < _full_upwind_threshold)
243  fullyUpwind(res_or_jac, ph, pvar);
244  else
245  {
246  switch (_fallback_scheme)
247  {
248  case FallbackEnum::QUICK:
249  quickUpwind(res_or_jac, ph, pvar);
250  break;
252  harmonicMean(res_or_jac, ph, pvar);
253  break;
254  }
255  }
256  }
257 
258  // Add results to the Residual or Jacobian
259  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
260  {
262  for (_i = 0; _i < _test.size(); _i++)
263  for (unsigned int ph = 0; ph < _num_phases; ++ph)
264  _local_re(_i) += _proto_flux[ph][_i];
266 
267  if (_has_save_in)
268  for (unsigned int i = 0; i < _save_in.size(); i++)
269  _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
270  }
271 
272  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
273  {
274  for (_i = 0; _i < _test.size(); _i++)
275  for (_j = 0; _j < _phi.size(); _j++)
276  for (unsigned int ph = 0; ph < _num_phases; ++ph)
277  _local_ke(_i, _j) += _jacobian[ph][_i][_j];
278 
280 
281  if (_has_diag_save_in && jvar == _var.number())
282  {
283  unsigned int rows = _local_ke.m();
284  DenseVector<Number> diag(rows);
285  for (unsigned int i = 0; i < rows; i++)
286  diag(i) = _local_ke(i, i);
287 
288  for (unsigned int i = 0; i < _diag_save_in.size(); i++)
289  _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
290  }
291  }
292 }
293 
294 void
295 PorousFlowDarcyBase::fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
296 {
326  // The number of nodes in the element
327  const unsigned int num_nodes = _test.size();
328 
329  Real mob;
330  Real dmob;
331  // Define variables used to ensure mass conservation
332  Real total_mass_out = 0.0;
333  Real total_in = 0.0;
334 
335  // The following holds derivatives of these
336  std::vector<Real> dtotal_mass_out;
337  std::vector<Real> dtotal_in;
338  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
339  {
340  dtotal_mass_out.assign(num_nodes, 0.0);
341  dtotal_in.assign(num_nodes, 0.0);
342  }
343 
344  // Perform the upwinding using the mobility
345  std::vector<bool> upwind_node(num_nodes);
346  for (unsigned int n = 0; n < num_nodes; ++n)
347  {
348  if (_proto_flux[ph][n] >= 0.0) // upstream node
349  {
350  upwind_node[n] = true;
351  // The mobility at the upstream node
352  mob = mobility(n, ph);
353  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
354  {
355  // The derivative of the mobility wrt the PorousFlow variable
356  dmob = dmobility(n, ph, pvar);
357 
358  for (_j = 0; _j < _phi.size(); _j++)
359  _jacobian[ph][n][_j] *= mob;
360 
361  if (_test.size() == _phi.size())
362  /* mobility at node=n depends only on the variables at node=n, by construction. For
363  * linear-lagrange variables, this means that Jacobian entries involving the derivative
364  * of mobility will only be nonzero for derivatives wrt variables at node=n. Hence the
365  * [n][n] in the line below. However, for other variable types (eg constant monomials)
366  * I cannot tell what variable number contributes to the derivative. However, in all
367  * cases I can possibly imagine, the derivative is zero anyway, since in the full
368  * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
369  */
370  _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
371 
372  for (_j = 0; _j < _phi.size(); _j++)
373  dtotal_mass_out[_j] += _jacobian[ph][n][_j];
374  }
375  _proto_flux[ph][n] *= mob;
376  total_mass_out += _proto_flux[ph][n];
377  }
378  else
379  {
380  upwind_node[n] = false;
381  total_in -= _proto_flux[ph][n];
382  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
383  for (_j = 0; _j < _phi.size(); _j++)
384  dtotal_in[_j] -= _jacobian[ph][n][_j];
385  }
386  }
387 
388  // Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes,
389  // weighted by their proto_flux values
390  for (unsigned int n = 0; n < num_nodes; ++n)
391  {
392  if (!upwind_node[n]) // downstream node
393  {
394  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
395  for (_j = 0; _j < _phi.size(); _j++)
396  {
397  _jacobian[ph][n][_j] *= total_mass_out / total_in;
398  _jacobian[ph][n][_j] +=
399  _proto_flux[ph][n] * (dtotal_mass_out[_j] / total_in -
400  dtotal_in[_j] * total_mass_out / total_in / total_in);
401  }
402  _proto_flux[ph][n] *= total_mass_out / total_in;
403  }
404  }
405 }
406 
407 void
408 PorousFlowDarcyBase::quickUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
409 {
410  // The number of nodes in the element
411  const unsigned int num_nodes = _test.size();
412 
413  Real mob;
414  Real dmob;
415 
416  // Use the raw nodal mobility
417  for (unsigned int n = 0; n < num_nodes; ++n)
418  {
419  // The mobility at the node
420  mob = mobility(n, ph);
421  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
422  {
423  // The derivative of the mobility wrt the PorousFlow variable
424  dmob = dmobility(n, ph, pvar);
425 
426  for (_j = 0; _j < _phi.size(); _j++)
427  _jacobian[ph][n][_j] *= mob;
428 
429  if (_test.size() == _phi.size())
430  /* mobility at node=n depends only on the variables at node=n, by construction. For
431  * linear-lagrange variables, this means that Jacobian entries involving the derivative
432  * of mobility will only be nonzero for derivatives wrt variables at node=n. Hence the
433  * [n][n] in the line below. However, for other variable types (eg constant monomials)
434  * I cannot tell what variable number contributes to the derivative. However, in all
435  * cases I can possibly imagine, the derivative is zero anyway, since in the full
436  * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
437  */
438  _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
439  }
440  _proto_flux[ph][n] *= mob;
441  }
442 }
443 
444 void
445 PorousFlowDarcyBase::harmonicMean(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
446 {
447  // The number of nodes in the element
448  const unsigned int num_nodes = _test.size();
449 
450  std::vector<Real> mob(num_nodes);
451  unsigned num_zero = 0;
452  unsigned zero_mobility_node = std::numeric_limits<unsigned>::max();
453  Real harmonic_mob = 0;
454  for (unsigned n = 0; n < num_nodes; ++n)
455  {
456  mob[n] = mobility(n, ph);
457  if (mob[n] == 0.0)
458  {
459  zero_mobility_node = n;
460  num_zero++;
461  }
462  else
463  harmonic_mob += 1.0 / mob[n];
464  }
465  if (num_zero > 0)
466  harmonic_mob = 0.0;
467  else
468  harmonic_mob = (1.0 * num_nodes) / harmonic_mob;
469 
470  // d(harmonic_mob)/d(PorousFlow variable at node n)
471  std::vector<Real> dharmonic_mob(num_nodes, 0.0);
472  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
473  {
474  const Real harm2 = std::pow(harmonic_mob, 2) / (1.0 * num_nodes);
475  if (num_zero == 0)
476  for (unsigned n = 0; n < num_nodes; ++n)
477  dharmonic_mob[n] = dmobility(n, ph, pvar) * harm2 / std::pow(mob[n], 2);
478  else if (num_zero == 1)
479  dharmonic_mob[zero_mobility_node] =
480  num_nodes * dmobility(zero_mobility_node, ph, pvar); // other derivs are zero
481  // if num_zero > 1 then all dharmonic_mob = 0.0
482  }
483 
484  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
485  for (unsigned n = 0; n < num_nodes; ++n)
486  for (_j = 0; _j < _phi.size(); _j++)
487  {
488  _jacobian[ph][n][_j] *= harmonic_mob;
489  if (_test.size() == _phi.size())
490  _jacobian[ph][n][_j] += dharmonic_mob[_j] * _proto_flux[ph][n];
491  }
492 
493  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
494  for (unsigned n = 0; n < num_nodes; ++n)
495  _proto_flux[ph][n] *= harmonic_mob;
496 }
497 
498 Real
499 PorousFlowDarcyBase::mobility(unsigned nodenum, unsigned phase) const
500 {
501  return _fluid_density_node[nodenum][phase] / _fluid_viscosity[nodenum][phase];
502 }
503 
504 Real
505 PorousFlowDarcyBase::dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
506 {
507  Real dm = _dfluid_density_node_dvar[nodenum][phase][pvar] / _fluid_viscosity[nodenum][phase];
508  dm -= _fluid_density_node[nodenum][phase] * _dfluid_viscosity_dvar[nodenum][phase][pvar] /
509  std::pow(_fluid_viscosity[nodenum][phase], 2);
510  return dm;
511 }
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_density_node_dvar
Derivative of the fluid density for each phase wrt PorousFlow variables (at the node) ...
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dpermeability_dgradvar
d(permeabiity)/d(grad(PorousFlow variable))
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_density_qp_dvar
Derivative of the fluid density for each phase wrt PorousFlow variables (at the qp) ...
const MaterialProperty< std::vector< RealTensorValue > > & _dpermeability_dvar
d(permeabiity)/d(PorousFlow variable)
unsigned int n_threads()
static InputParameters validParams()
bool notPorousFlowVariable(unsigned int moose_var_num) const
Returns true if moose_var_num is not a porous flow variabe.
void accumulateTaggedLocalResidual()
MooseVariable & _var
std::vector< MooseVariableFEBase *> _save_in
const MooseArray< Real > & _JxW
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
enum PorousFlowDarcyBase::FallbackEnum _fallback_scheme
unsigned int number() const
const VariablePhiGradient & _grad_phi
const MaterialProperty< std::vector< Real > > & _fluid_viscosity
Viscosity of each component in each phase.
const MaterialProperty< RealTensorValue > & _permeability
Permeability of porous material.
FallbackEnum
If full upwinding is failing due to nodes swapping between upwind and downwind in successive nonlinea...
std::unordered_map< unsigned, std::vector< std::vector< unsigned > > > _num_downwinds
Number of nonlinear iterations (in this timestep and this element) that a node is an downwind node fo...
std::vector< std::vector< Real > > _proto_flux
The Darcy flux.
static constexpr std::size_t dim
virtual void computeResidualAndJacobian() override
std::vector< std::vector< std::vector< Real > > > _jacobian
Derivative of _proto_flux with respect to nodal variables.
virtual Real darcyQpJacobian(unsigned int jvar, unsigned int ph) const
Jacobian of the Darcy part of the flux.
unsigned int m() const
const MaterialProperty< std::vector< Real > > & _fluid_density_qp
Fluid density for each phase (at the qp)
bool _has_diag_save_in
DenseMatrix< Number > _local_ke
void mooseWarning(Args &&... args) const
void addRequiredParam(const std::string &name, const std::string &doc_string)
virtual void computeResidual() override
const MaterialProperty< std::vector< RealGradient > > & _grad_p
Gradient of the pore pressure in each phase.
TensorValue< Real > RealTensorValue
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_viscosity_dvar
Derivative of the fluid viscosity for each phase wrt PorousFlow variables.
const VariableTestValue & _test
virtual void timestepSetup()
const MaterialProperty< std::vector< std::vector< RealGradient > > > & _dgrad_p_dvar
Derivative of Grad porepressure in each phase wrt PorousFlow variables.
Real deriv(unsigned n, unsigned alpha, unsigned beta, Real x)
static InputParameters validParams()
std::vector< MooseVariableFEBase *> _diag_save_in
const PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
const unsigned int _num_phases
The number of fluid phases.
const QBase *const & _qrule
bool _has_save_in
unsigned int _i
virtual void computeOffDiagJacobian(unsigned int jvar) override
std::unordered_map< unsigned, std::vector< std::vector< unsigned > > > _num_upwinds
Number of nonlinear iterations (in this timestep and this element) that a node is an upwind node for ...
void accumulateTaggedLocalMatrix()
void harmonicMean(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
Calculate the residual or Jacobian by using the harmonic mean of the nodal mobilities for the entire ...
virtual Real computeQpResidual() override
virtual void timestepSetup() override
const MooseArray< Real > & _coord
const RealVectorValue _gravity
Gravity. Defaults to 9.81 m/s^2.
Assembly & _assembly
unsigned int _j
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual Real mobility(unsigned nodenum, unsigned phase) const
The mobility of the fluid.
PorousFlowDarcyBase(const InputParameters &parameters)
const MaterialProperty< std::vector< std::vector< Real > > > & _dgrad_p_dgrad_var
Derivative of Grad porepressure in each phase wrt grad(PorousFlow variables)
void quickUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
Calculate the residual or Jacobian using the nodal mobilities, but without conserving fluid mass...
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
const VariableTestGradient & _grad_test
virtual Real dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
The derivative of mobility with respect to PorousFlow variable pvar.
const unsigned _full_upwind_threshold
If the number of upwind-downwind swaps is less than this amount then full upwinding is used...
DenseVector< Number > _local_re
virtual void computeJacobian() override
IntRange< T > make_range(T beg, T end)
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
const Elem *const & _current_elem
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
The PorousFlow variable number.
unsigned int n() const
void prepareVectorTag(Assembly &assembly, unsigned int ivar)
void prepareMatrixTag(Assembly &assembly, unsigned int ivar, unsigned int jvar)
virtual Real darcyQp(unsigned int ph) const
The Darcy part of the flux (this is the non-upwinded part)
const VariablePhiValue & _phi
const bool _perm_derivs
Flag to check whether permeabiity derivatives are non-zero.
MooseUnits pow(const MooseUnits &, int)
const MaterialProperty< std::vector< Real > > & _fluid_density_node
Fluid density for each phase (at the node)
unsigned int _qp
void fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
Calculate the residual or Jacobian using full upwinding.