www.mooseframework.org
PorousFlowDarcyBase.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 "PorousFlowDarcyBase.h"
9 
10 // MOOSE includes
11 #include "Assembly.h"
12 #include "MooseMesh.h"
13 #include "MooseVariable.h"
14 #include "SystemBase.h"
15 
16 #include "libmesh/quadrature.h"
17 
18 template <>
19 InputParameters
21 {
22  InputParameters params = validParams<Kernel>();
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 
44 PorousFlowDarcyBase::PorousFlowDarcyBase(const InputParameters & parameters)
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  _porousflow_dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
68  _num_phases(_porousflow_dictator.numPhases()),
69  _gravity(getParam<RealVectorValue>("gravity")),
70  _full_upwind_threshold(getParam<unsigned>("full_upwind_threshold")),
71  _fallback_scheme(getParam<MooseEnum>("fallback_scheme").getEnum<FallbackEnum>()),
72  _proto_flux(_num_phases),
73  _jacobian(_num_phases),
74  _num_upwinds(),
75  _num_downwinds()
76 {
77 #ifdef LIBMESH_HAVE_TBB_API
78  if (libMesh::n_threads() > 1)
79  mooseWarning("PorousFlowDarcyBase: num_upwinds and num_downwinds may not be computed "
80  "accurately when using TBB and greater than 1 thread");
81 #endif
82 }
83 
84 void
86 {
87  Kernel::timestepSetup();
88  _num_upwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
89  _num_downwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
90 }
91 
92 Real
93 PorousFlowDarcyBase::darcyQp(unsigned int ph) const
94 {
95  return _grad_test[_i][_qp] *
96  (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
97 }
98 
99 Real
100 PorousFlowDarcyBase::darcyQpJacobian(unsigned int jvar, unsigned int ph) const
101 {
103  return 0.0;
104 
105  const unsigned int pvar = _porousflow_dictator.porousFlowVariableNum(jvar);
106  RealVectorValue deriv = _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
107  (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
108  for (unsigned i = 0; i < LIBMESH_DIM; ++i)
109  deriv += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
110  (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
111  deriv += _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
112  _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
113  deriv += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
114  return _grad_test[_i][_qp] * deriv;
115 }
116 
117 Real
119 {
120  mooseError("PorousFlowDarcyBase: computeQpResidual called");
121  return 0.0;
122 }
123 
124 void
126 {
128 }
129 
130 void
132 {
134 }
135 
136 void
138 {
140 }
141 
142 void
144 {
145  if ((res_or_jac == JacRes::CALCULATE_JACOBIAN) &&
147  return;
148 
150  const unsigned int pvar =
152  : 0);
153 
154  DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
155  if ((ke.n() == 0) &&
156  (res_or_jac == JacRes::CALCULATE_JACOBIAN)) // this removes a problem encountered in
157  // the initial timestep when
158  // use_displaced_mesh=true
159  return;
160 
162  const unsigned int num_nodes = _test.size();
163 
166  for (unsigned ph = 0; ph < _num_phases; ++ph)
167  {
168  _proto_flux[ph].assign(num_nodes, 0);
169  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
170  {
171  for (_i = 0; _i < num_nodes; ++_i)
172  _proto_flux[ph][_i] += _JxW[_qp] * _coord[_qp] * darcyQp(ph);
173  }
174  }
175 
176  // for this element, record whether each node is "upwind" or "downwind" (or neither)
177  const unsigned elem = _current_elem->id();
178  if (_num_upwinds.find(elem) == _num_upwinds.end())
179  {
180  _num_upwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
181  _num_downwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
182  for (unsigned ph = 0; ph < _num_phases; ++ph)
183  {
184  _num_upwinds[elem][ph].assign(num_nodes, 0);
185  _num_downwinds[elem][ph].assign(num_nodes, 0);
186  }
187  }
188  // record the information once per nonlinear iteration
189  if (res_or_jac == JacRes::CALCULATE_JACOBIAN && jvar == _var.number())
190  {
191  for (unsigned ph = 0; ph < _num_phases; ++ph)
192  {
193  for (unsigned nod = 0; nod < num_nodes; ++nod)
194  {
195  if (_proto_flux[ph][nod] > 0)
196  _num_upwinds[elem][ph][nod]++;
197  else if (_proto_flux[ph][nod] < 0)
198  _num_downwinds[elem][ph][nod]++;
199  }
200  }
201  }
202 
203  // based on _num_upwinds and _num_downwinds, calculate the maximum number
204  // of upwind-downwind swaps that have been encountered in this timestep
205  // for this element
206  std::vector<unsigned> max_swaps(_num_phases, 0);
207  for (unsigned ph = 0; ph < _num_phases; ++ph)
208  {
209  for (unsigned nod = 0; nod < num_nodes; ++nod)
210  max_swaps[ph] = std::max(
211  max_swaps[ph], std::min(_num_upwinds[elem][ph][nod], _num_downwinds[elem][ph][nod]));
212  }
213 
214  // size the _jacobian correctly and calculate it for the case residual = _proto_flux
215  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
216  {
217  for (unsigned ph = 0; ph < _num_phases; ++ph)
218  {
219  _jacobian[ph].resize(ke.m());
220  for (_i = 0; _i < _test.size(); _i++)
221  {
222  _jacobian[ph][_i].assign(ke.n(), 0.0);
223  for (_j = 0; _j < _phi.size(); _j++)
224  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
225  _jacobian[ph][_i][_j] += _JxW[_qp] * _coord[_qp] * darcyQpJacobian(jvar, ph);
226  }
227  }
228  }
229 
230  // Loop over all the phases, computing the mass flux, which
231  // gets placed into _proto_flux, and the derivative of this
232  // which gets placed into _jacobian
233  for (unsigned int ph = 0; ph < _num_phases; ++ph)
234  {
235  if (max_swaps[ph] < _full_upwind_threshold)
236  fullyUpwind(res_or_jac, ph, pvar);
237  else
238  {
239  switch (_fallback_scheme)
240  {
241  case FallbackEnum::QUICK:
242  quickUpwind(res_or_jac, ph, pvar);
243  break;
245  harmonicMean(res_or_jac, ph, pvar);
246  break;
247  }
248  }
249  }
250 
252  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
253  {
254  DenseVector<Number> & re = _assembly.residualBlock(_var.number());
255 
256  _local_re.resize(re.size());
257  _local_re.zero();
258  for (_i = 0; _i < _test.size(); _i++)
259  for (unsigned int ph = 0; ph < _num_phases; ++ph)
260  _local_re(_i) += _proto_flux[ph][_i];
261 
262  re += _local_re;
263 
264  if (_has_save_in)
265  {
266  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
267  for (unsigned int i = 0; i < _save_in.size(); i++)
268  _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
269  }
270  }
271 
272  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
273  {
274  _local_ke.resize(ke.m(), ke.n());
275  _local_ke.zero();
276 
277  for (_i = 0; _i < _test.size(); _i++)
278  for (_j = 0; _j < _phi.size(); _j++)
279  for (unsigned int ph = 0; ph < _num_phases; ++ph)
280  _local_ke(_i, _j) += _jacobian[ph][_i][_j];
281 
282  ke += _local_ke;
283 
284  if (_has_diag_save_in && jvar == _var.number())
285  {
286  unsigned int rows = ke.m();
287  DenseVector<Number> diag(rows);
288  for (unsigned int i = 0; i < rows; i++)
289  diag(i) = _local_ke(i, i);
290 
291  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
292  for (unsigned int i = 0; i < _diag_save_in.size(); i++)
293  _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
294  }
295  }
296 }
297 
298 void
299 PorousFlowDarcyBase::fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
300 {
330  const unsigned int num_nodes = _test.size();
332 
333  Real mob;
334  Real dmob;
336  Real total_mass_out = 0.0;
337  Real total_in = 0.0;
338 
340  std::vector<Real> dtotal_mass_out;
341  std::vector<Real> dtotal_in;
342  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
343  {
344  dtotal_mass_out.assign(num_nodes, 0.0);
345  dtotal_in.assign(num_nodes, 0.0);
346  }
347 
349  std::vector<bool> upwind_node(num_nodes);
350  for (unsigned int n = 0; n < num_nodes; ++n)
351  {
352  if (_proto_flux[ph][n] >= 0.0) // upstream node
353  {
354  upwind_node[n] = true;
356  mob = mobility(n, ph);
357  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
358  {
360  dmob = dmobility(n, ph, pvar);
361 
362  for (_j = 0; _j < _phi.size(); _j++)
363  _jacobian[ph][n][_j] *= mob;
364 
365  if (_test.size() == _phi.size())
366  /* mobility at node=n depends only on the variables at node=n, by construction. For
367  * linear-lagrange variables, this means that Jacobian entries involving the derivative
368  * of mobility will only be nonzero for derivatives wrt variables at node=n. Hence the
369  * [n][n] in the line below. However, for other variable types (eg constant monomials)
370  * I cannot tell what variable number contributes to the derivative. However, in all
371  * cases I can possibly imagine, the derivative is zero anyway, since in the full
372  * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
373  */
374  _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
375 
376  for (_j = 0; _j < _phi.size(); _j++)
377  dtotal_mass_out[_j] += _jacobian[ph][n][_j];
378  }
379  _proto_flux[ph][n] *= mob;
380  total_mass_out += _proto_flux[ph][n];
381  }
382  else
383  {
384  upwind_node[n] = false;
385  total_in -= _proto_flux[ph][n];
386  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
387  for (_j = 0; _j < _phi.size(); _j++)
388  dtotal_in[_j] -= _jacobian[ph][n][_j];
389  }
390  }
391 
393  for (unsigned int n = 0; n < num_nodes; ++n)
394  {
395  if (!upwind_node[n]) // downstream node
396  {
397  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
398  for (_j = 0; _j < _phi.size(); _j++)
399  {
400  _jacobian[ph][n][_j] *= total_mass_out / total_in;
401  _jacobian[ph][n][_j] +=
402  _proto_flux[ph][n] * (dtotal_mass_out[_j] / total_in -
403  dtotal_in[_j] * total_mass_out / total_in / total_in);
404  }
405  _proto_flux[ph][n] *= total_mass_out / total_in;
406  }
407  }
408 }
409 
410 void
411 PorousFlowDarcyBase::quickUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
412 {
414  const unsigned int num_nodes = _test.size();
415 
416  Real mob;
417  Real dmob;
418 
420  for (unsigned int n = 0; n < num_nodes; ++n)
421  {
423  mob = mobility(n, ph);
424  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
425  {
427  dmob = dmobility(n, ph, pvar);
428 
429  for (_j = 0; _j < _phi.size(); _j++)
430  _jacobian[ph][n][_j] *= mob;
431 
432  if (_test.size() == _phi.size())
433  /* mobility at node=n depends only on the variables at node=n, by construction. For
434  * linear-lagrange variables, this means that Jacobian entries involving the derivative
435  * of mobility will only be nonzero for derivatives wrt variables at node=n. Hence the
436  * [n][n] in the line below. However, for other variable types (eg constant monomials)
437  * I cannot tell what variable number contributes to the derivative. However, in all
438  * cases I can possibly imagine, the derivative is zero anyway, since in the full
439  * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
440  */
441  _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
442  }
443  _proto_flux[ph][n] *= mob;
444  }
445 }
446 
447 void
448 PorousFlowDarcyBase::harmonicMean(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
449 {
451  const unsigned int num_nodes = _test.size();
452 
453  std::vector<Real> mob(num_nodes);
454  unsigned num_zero = 0;
455  unsigned zero_mobility_node;
456  Real harmonic_mob = 0;
457  for (unsigned n = 0; n < num_nodes; ++n)
458  {
459  mob[n] = mobility(n, ph);
460  if (mob[n] == 0.0)
461  {
462  zero_mobility_node = n;
463  num_zero++;
464  }
465  else
466  harmonic_mob += 1.0 / mob[n];
467  }
468  if (num_zero > 0)
469  harmonic_mob = 0.0;
470  else
471  harmonic_mob = (1.0 * num_nodes) / harmonic_mob;
472 
473  /*
474  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
475  {
476  Moose::out << "RES ";
477  for (unsigned n = 0; n < num_nodes; ++n)
478  Moose::out << std::setprecision(16) << _pp[n][0] << " ";
479  Moose::out << "\n";
480  for (unsigned n = 0; n < num_nodes; ++n)
481  Moose::out << mob[n] << " ";
482  Moose::out << "harmonic_mob = " << harmonic_mob << "\n";
483  }
484  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
485  {
486  Moose::out << "JAC ";
487  for (unsigned n = 0; n < num_nodes; ++n)
488  Moose::out << mob[n] << " ";
489  Moose::out << "harmonic_mob = " << harmonic_mob << "\n";
490  }
491  */
492 
493  // d(harmonic_mob)/d(PorousFlow variable at node n)
494  std::vector<Real> dharmonic_mob(num_nodes, 0.0);
495  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
496  {
497  const Real harm2 = std::pow(harmonic_mob, 2) / (1.0 * num_nodes);
498  if (num_zero == 0)
499  for (unsigned n = 0; n < num_nodes; ++n)
500  dharmonic_mob[n] = dmobility(n, ph, pvar) * harm2 / std::pow(mob[n], 2);
501  else if (num_zero == 1)
502  dharmonic_mob[zero_mobility_node] =
503  num_nodes * dmobility(zero_mobility_node, ph, pvar); // other derivs are zero
504  // if num_zero > 1 then all dharmonic_mob = 0.0
505  }
506 
507  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
508  for (unsigned n = 0; n < num_nodes; ++n)
509  for (_j = 0; _j < _phi.size(); _j++)
510  {
511  _jacobian[ph][n][_j] *= harmonic_mob;
512  if (_test.size() == _phi.size())
513  _jacobian[ph][n][_j] += dharmonic_mob[_j] * _proto_flux[ph][n];
514  }
515 
516  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
517  for (unsigned n = 0; n < num_nodes; ++n)
518  _proto_flux[ph][n] *= harmonic_mob;
519 }
520 
521 Real
522 PorousFlowDarcyBase::mobility(unsigned nodenum, unsigned phase) const
523 {
524  return _fluid_density_node[nodenum][phase] / _fluid_viscosity[nodenum][phase];
525 }
526 
527 Real
528 PorousFlowDarcyBase::dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
529 {
530  Real dm = _dfluid_density_node_dvar[nodenum][phase][pvar] / _fluid_viscosity[nodenum][phase];
531  dm -= _fluid_density_node[nodenum][phase] * _dfluid_viscosity_dvar[nodenum][phase][pvar] /
532  std::pow(_fluid_viscosity[nodenum][phase], 2);
533  return dm;
534 }
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(porous-flow 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(porous-flow variable)
enum PorousFlowDarcyBase::FallbackEnum _fallback_scheme
void computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
Computation of the residual and Jacobian.
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.
std::vector< std::vector< std::vector< Real > > > _jacobian
Derivative of _proto_flux with respect to nodal variables.
const MaterialProperty< std::vector< Real > > & _fluid_density_qp
Fluid density for each phase (at the qp)
virtual void computeResidual() override
const MaterialProperty< std::vector< RealGradient > > & _grad_p
Gradient of the pore pressure in each phase.
const PorousFlowDictator & _porousflow_dictator
PorousFlow UserObject.
virtual Real darcyQp(unsigned int ph) const
the Darcy part of the flux (this is the non-upwinded part)
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_viscosity_dvar
Derivative of the fluid viscosity for each phase wrt PorousFlow variables.
const MaterialProperty< std::vector< std::vector< RealGradient > > > & _dgrad_p_dvar
Derivative of Grad porepressure in each phase wrt PorousFlow variables.
const unsigned int _num_phases
The number of fluid phases.
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 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 RealVectorValue _gravity
Gravity. Defaults to 9.81 m/s^2.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
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...
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...
virtual Real mobility(unsigned nodenum, unsigned phase) const
The mobility of the fluid.
virtual void computeJacobian() override
InputParameters validParams< PorousFlowDarcyBase >()
bool notPorousFlowVariable(unsigned int moose_var_num) const
returns true if moose_var_num is not a porous flow variabe
virtual Real darcyQpJacobian(unsigned int jvar, unsigned int ph) const
Jacobian of the Darcy part of the flux.
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
the PorousFlow variable number
const MaterialProperty< std::vector< Real > > & _fluid_density_node
Fluid density for each phase (at the node)
void fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
Calculate the residual or Jacobian using full upwinding.