www.mooseframework.org
Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | List of all members
PorousFlowDarcyBase Class Reference

Darcy advective flux. More...

#include <PorousFlowDarcyBase.h>

Inheritance diagram for PorousFlowDarcyBase:
[legend]

Public Member Functions

 PorousFlowDarcyBase (const InputParameters &parameters)
 

Protected Types

enum  JacRes { JacRes::CALCULATE_RESIDUAL = 0, JacRes::CALCULATE_JACOBIAN = 1 }
 
enum  FallbackEnum { FallbackEnum::QUICK, FallbackEnum::HARMONIC }
 If full upwinding is failing due to nodes swapping between upwind and downwind in successive nonlinear iterations (within one timestep and element) a fallback scheme is used instead. More...
 

Protected Member Functions

virtual void timestepSetup () override
 
virtual Real computeQpResidual () override
 
virtual void computeResidual () override
 
virtual void computeJacobian () override
 
virtual void computeOffDiagJacobian (MooseVariableFEBase &jvar) override
 
virtual Real darcyQp (unsigned int ph) const
 The Darcy part of the flux (this is the non-upwinded part) More...
 
virtual Real darcyQpJacobian (unsigned int jvar, unsigned int ph) const
 Jacobian of the Darcy part of the flux. More...
 
virtual Real mobility (unsigned nodenum, unsigned phase) const
 The mobility of the fluid. More...
 
virtual Real dmobility (unsigned nodenum, unsigned phase, unsigned pvar) const
 The derivative of mobility with respect to PorousFlow variable pvar. More...
 
void computeResidualAndJacobian (JacRes res_or_jac, unsigned int jvar)
 Computation of the residual and Jacobian. More...
 
void fullyUpwind (JacRes res_or_jac, unsigned int ph, unsigned int pvar)
 Calculate the residual or Jacobian using full upwinding. More...
 
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. More...
 
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 element. More...
 

Protected Attributes

const MaterialProperty< RealTensorValue > & _permeability
 Permeability of porous material. More...
 
const MaterialProperty< std::vector< RealTensorValue > > & _dpermeability_dvar
 d(permeabiity)/d(PorousFlow variable) More...
 
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dpermeability_dgradvar
 d(permeabiity)/d(grad(PorousFlow variable)) More...
 
const MaterialProperty< std::vector< Real > > & _fluid_density_node
 Fluid density for each phase (at the node) More...
 
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) More...
 
const MaterialProperty< std::vector< Real > > & _fluid_density_qp
 Fluid density for each phase (at the qp) More...
 
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) More...
 
const MaterialProperty< std::vector< Real > > & _fluid_viscosity
 Viscosity of each component in each phase. More...
 
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_viscosity_dvar
 Derivative of the fluid viscosity for each phase wrt PorousFlow variables. More...
 
const MaterialProperty< std::vector< Real > > & _pp
 Nodal pore pressure in each phase. More...
 
const MaterialProperty< std::vector< RealGradient > > & _grad_p
 Gradient of the pore pressure in each phase. More...
 
const MaterialProperty< std::vector< std::vector< Real > > > & _dgrad_p_dgrad_var
 Derivative of Grad porepressure in each phase wrt grad(PorousFlow variables) More...
 
const MaterialProperty< std::vector< std::vector< RealGradient > > > & _dgrad_p_dvar
 Derivative of Grad porepressure in each phase wrt PorousFlow variables. More...
 
const PorousFlowDictator_dictator
 PorousFlowDictator UserObject. More...
 
const unsigned int _num_phases
 The number of fluid phases. More...
 
const RealVectorValue _gravity
 Gravity. Defaults to 9.81 m/s^2. More...
 
const unsigned _full_upwind_threshold
 If the number of upwind-downwind swaps is less than this amount then full upwinding is used. More...
 
enum PorousFlowDarcyBase::FallbackEnum _fallback_scheme
 
std::vector< std::vector< Real > > _proto_flux
 The Darcy flux. More...
 
std::vector< std::vector< std::vector< Real > > > _jacobian
 Derivative of _proto_flux with respect to nodal variables. More...
 
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 a given fluid phase. More...
 
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 for a given fluid phase. More...
 

Detailed Description

Darcy advective flux.

A fully-updwinded version is implemented, where the mobility of the upstream nodes is used. Alternately, after some number of upwind-downwind swaps, another type of handling of the mobility may be employed (see fallback_scheme, below)

Definition at line 29 of file PorousFlowDarcyBase.h.

Member Enumeration Documentation

◆ FallbackEnum

enum PorousFlowDarcyBase::FallbackEnum
strongprotected

If full upwinding is failing due to nodes swapping between upwind and downwind in successive nonlinear iterations (within one timestep and element) a fallback scheme is used instead.

QUICK: use the nodal mobilities, as in full upwinding, but do not attempt to conserve fluid mass. HARMONIC: assign the harmonic mean of the mobilities to the entire element.

Enumerator
QUICK 
HARMONIC 

Definition at line 150 of file PorousFlowDarcyBase.h.

150 { QUICK, HARMONIC } _fallback_scheme;
enum PorousFlowDarcyBase::FallbackEnum _fallback_scheme

◆ JacRes

enum PorousFlowDarcyBase::JacRes
strongprotected
Enumerator
CALCULATE_RESIDUAL 
CALCULATE_JACOBIAN 

Definition at line 64 of file PorousFlowDarcyBase.h.

65  {
66  CALCULATE_RESIDUAL = 0,
67  CALCULATE_JACOBIAN = 1
68  };

Constructor & Destructor Documentation

◆ PorousFlowDarcyBase()

PorousFlowDarcyBase::PorousFlowDarcyBase ( const InputParameters &  parameters)

Definition at line 45 of file PorousFlowDarcyBase.C.

46  : Kernel(parameters),
47  _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
49  getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
50  _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
51  "dPorousFlow_permeability_qp_dgradvar")),
53  getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
54  _dfluid_density_node_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
55  "dPorousFlow_fluid_phase_density_nodal_dvar")),
56  _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
57  _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
58  "dPorousFlow_fluid_phase_density_qp_dvar")),
59  _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")),
61  getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
62  _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_nodal")),
63  _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
64  _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
65  "dPorousFlow_grad_porepressure_qp_dgradvar")),
66  _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
67  "dPorousFlow_grad_porepressure_qp_dvar")),
68  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
70  _gravity(getParam<RealVectorValue>("gravity")),
71  _full_upwind_threshold(getParam<unsigned>("full_upwind_threshold")),
72  _fallback_scheme(getParam<MooseEnum>("fallback_scheme").getEnum<FallbackEnum>()),
75  _num_upwinds(),
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 }
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)
enum PorousFlowDarcyBase::FallbackEnum _fallback_scheme
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)
const MaterialProperty< std::vector< RealGradient > > & _grad_p
Gradient of the pore pressure in each phase.
const MaterialProperty< std::vector< Real > > & _pp
Nodal pore pressure in each phase.
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 PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
const unsigned int _num_phases
The number of fluid phases.
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 ...
unsigned int numPhases() const
The number of fluid phases.
const RealVectorValue _gravity
Gravity. Defaults to 9.81 m/s^2.
const MaterialProperty< std::vector< std::vector< Real > > > & _dgrad_p_dgrad_var
Derivative of Grad porepressure in each phase wrt grad(PorousFlow variables)
const unsigned _full_upwind_threshold
If the number of upwind-downwind swaps is less than this amount then full upwinding is used...
const MaterialProperty< std::vector< Real > > & _fluid_density_node
Fluid density for each phase (at the node)

Member Function Documentation

◆ computeJacobian()

void PorousFlowDarcyBase::computeJacobian ( )
overrideprotectedvirtual

Definition at line 132 of file PorousFlowDarcyBase.C.

133 {
135 }
void computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
Computation of the residual and Jacobian.

◆ computeOffDiagJacobian()

void PorousFlowDarcyBase::computeOffDiagJacobian ( MooseVariableFEBase &  jvar)
overrideprotectedvirtual

Definition at line 138 of file PorousFlowDarcyBase.C.

139 {
141 }
void computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
Computation of the residual and Jacobian.

◆ computeQpResidual()

Real PorousFlowDarcyBase::computeQpResidual ( )
overrideprotectedvirtual

Definition at line 119 of file PorousFlowDarcyBase.C.

120 {
121  mooseError("PorousFlowDarcyBase: computeQpResidual called");
122  return 0.0;
123 }

◆ computeResidual()

void PorousFlowDarcyBase::computeResidual ( )
overrideprotectedvirtual

Definition at line 126 of file PorousFlowDarcyBase.C.

127 {
129 }
void computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
Computation of the residual and Jacobian.

◆ computeResidualAndJacobian()

void PorousFlowDarcyBase::computeResidualAndJacobian ( JacRes  res_or_jac,
unsigned int  jvar 
)
protected

Computation of the residual and Jacobian.

If res_or_jac=CALCULATE_JACOBIAN then the residual gets calculated anyway (becuase we have to know which nodes are upwind and which are downwind) except we don't actually need to store the residual in _assembly.residualBlock and we don't need to save in save_in.

If res_or_jac=CALCULATE_RESIDUAL then we don't have to do all the Jacobian calculations.

Parameters
res_or_jacWhether to calculate the residual or the jacobian.
jvarPorousFlow variable to take derivatives wrt to

Definition at line 144 of file PorousFlowDarcyBase.C.

Referenced by computeJacobian(), computeOffDiagJacobian(), and computeResidual().

145 {
146  if ((res_or_jac == JacRes::CALCULATE_JACOBIAN) && _dictator.notPorousFlowVariable(jvar))
147  return;
148 
149  // The PorousFlow variable index corresponding to the variable number jvar
150  const unsigned int pvar =
151  ((res_or_jac == JacRes::CALCULATE_JACOBIAN) ? _dictator.porousFlowVariableNum(jvar) : 0);
152 
153  DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
154  if ((ke.n() == 0) && (res_or_jac == JacRes::CALCULATE_JACOBIAN)) // this removes a problem
155  // encountered in the initial
156  // timestep when
157  // use_displaced_mesh=true
158  return;
159 
160  // The number of nodes in the element
161  const unsigned int num_nodes = _test.size();
162 
163  // Compute the residual and jacobian without the mobility terms. Even if we are computing the
164  // Jacobian we still need this in order to see which nodes are upwind and which are downwind.
165  for (unsigned ph = 0; ph < _num_phases; ++ph)
166  {
167  _proto_flux[ph].assign(num_nodes, 0);
168  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
169  {
170  for (_i = 0; _i < num_nodes; ++_i)
171  _proto_flux[ph][_i] += _JxW[_qp] * _coord[_qp] * darcyQp(ph);
172  }
173  }
174 
175  // for this element, record whether each node is "upwind" or "downwind" (or neither)
176  const unsigned elem = _current_elem->id();
177  if (_num_upwinds.find(elem) == _num_upwinds.end())
178  {
179  _num_upwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
180  _num_downwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
181  for (unsigned ph = 0; ph < _num_phases; ++ph)
182  {
183  _num_upwinds[elem][ph].assign(num_nodes, 0);
184  _num_downwinds[elem][ph].assign(num_nodes, 0);
185  }
186  }
187  // record the information once per nonlinear iteration
188  if (res_or_jac == JacRes::CALCULATE_JACOBIAN && jvar == _var.number())
189  {
190  for (unsigned ph = 0; ph < _num_phases; ++ph)
191  {
192  for (unsigned nod = 0; nod < num_nodes; ++nod)
193  {
194  if (_proto_flux[ph][nod] > 0)
195  _num_upwinds[elem][ph][nod]++;
196  else if (_proto_flux[ph][nod] < 0)
197  _num_downwinds[elem][ph][nod]++;
198  }
199  }
200  }
201 
202  // based on _num_upwinds and _num_downwinds, calculate the maximum number
203  // of upwind-downwind swaps that have been encountered in this timestep
204  // for this element
205  std::vector<unsigned> max_swaps(_num_phases, 0);
206  for (unsigned ph = 0; ph < _num_phases; ++ph)
207  {
208  for (unsigned nod = 0; nod < num_nodes; ++nod)
209  max_swaps[ph] = std::max(
210  max_swaps[ph], std::min(_num_upwinds[elem][ph][nod], _num_downwinds[elem][ph][nod]));
211  }
212 
213  // size the _jacobian correctly and calculate it for the case residual = _proto_flux
214  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
215  {
216  for (unsigned ph = 0; ph < _num_phases; ++ph)
217  {
218  _jacobian[ph].resize(ke.m());
219  for (_i = 0; _i < _test.size(); _i++)
220  {
221  _jacobian[ph][_i].assign(ke.n(), 0.0);
222  for (_j = 0; _j < _phi.size(); _j++)
223  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
224  _jacobian[ph][_i][_j] += _JxW[_qp] * _coord[_qp] * darcyQpJacobian(jvar, ph);
225  }
226  }
227  }
228 
229  // Loop over all the phases, computing the mass flux, which
230  // gets placed into _proto_flux, and the derivative of this
231  // which gets placed into _jacobian
232  for (unsigned int ph = 0; ph < _num_phases; ++ph)
233  {
234  if (max_swaps[ph] < _full_upwind_threshold)
235  fullyUpwind(res_or_jac, ph, pvar);
236  else
237  {
238  switch (_fallback_scheme)
239  {
240  case FallbackEnum::QUICK:
241  quickUpwind(res_or_jac, ph, pvar);
242  break;
244  harmonicMean(res_or_jac, ph, pvar);
245  break;
246  }
247  }
248  }
249 
250  // Add results to the Residual or Jacobian
251  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
252  {
253  DenseVector<Number> & re = _assembly.residualBlock(_var.number());
254 
255  _local_re.resize(re.size());
256  _local_re.zero();
257  for (_i = 0; _i < _test.size(); _i++)
258  for (unsigned int ph = 0; ph < _num_phases; ++ph)
259  _local_re(_i) += _proto_flux[ph][_i];
260 
261  re += _local_re;
262 
263  if (_has_save_in)
264  {
265  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
266  for (unsigned int i = 0; i < _save_in.size(); i++)
267  _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
268  }
269  }
270 
271  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
272  {
273  _local_ke.resize(ke.m(), ke.n());
274  _local_ke.zero();
275 
276  for (_i = 0; _i < _test.size(); _i++)
277  for (_j = 0; _j < _phi.size(); _j++)
278  for (unsigned int ph = 0; ph < _num_phases; ++ph)
279  _local_ke(_i, _j) += _jacobian[ph][_i][_j];
280 
281  ke += _local_ke;
282 
283  if (_has_diag_save_in && jvar == _var.number())
284  {
285  unsigned int rows = ke.m();
286  DenseVector<Number> diag(rows);
287  for (unsigned int i = 0; i < rows; i++)
288  diag(i) = _local_ke(i, i);
289 
290  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
291  for (unsigned int i = 0; i < _diag_save_in.size(); i++)
292  _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
293  }
294  }
295 }
bool notPorousFlowVariable(unsigned int moose_var_num) const
Returns true if moose_var_num is not a porous flow variabe.
enum PorousFlowDarcyBase::FallbackEnum _fallback_scheme
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.
virtual Real darcyQpJacobian(unsigned int jvar, unsigned int ph) const
Jacobian of the Darcy part of the flux.
const PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
const unsigned int _num_phases
The number of fluid phases.
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 ...
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...
const unsigned _full_upwind_threshold
If the number of upwind-downwind swaps is less than this amount then full upwinding is used...
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
The PorousFlow variable number.
virtual Real darcyQp(unsigned int ph) const
The Darcy part of the flux (this is the non-upwinded part)
void fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
Calculate the residual or Jacobian using full upwinding.

◆ darcyQp()

Real PorousFlowDarcyBase::darcyQp ( unsigned int  ph) const
protectedvirtual

The Darcy part of the flux (this is the non-upwinded part)

Definition at line 94 of file PorousFlowDarcyBase.C.

Referenced by computeResidualAndJacobian().

95 {
96  return _grad_test[_i][_qp] *
97  (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
98 }
const MaterialProperty< RealTensorValue > & _permeability
Permeability of porous material.
const MaterialProperty< std::vector< Real > > & _fluid_density_qp
Fluid density for each phase (at the qp)
const MaterialProperty< std::vector< RealGradient > > & _grad_p
Gradient of the pore pressure in each phase.
const RealVectorValue _gravity
Gravity. Defaults to 9.81 m/s^2.

◆ darcyQpJacobian()

Real PorousFlowDarcyBase::darcyQpJacobian ( unsigned int  jvar,
unsigned int  ph 
) const
protectedvirtual

Jacobian of the Darcy part of the flux.

Definition at line 101 of file PorousFlowDarcyBase.C.

Referenced by computeResidualAndJacobian().

102 {
104  return 0.0;
105 
106  const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
107  RealVectorValue deriv = _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
108  (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
109  for (unsigned i = 0; i < LIBMESH_DIM; ++i)
110  deriv += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
111  (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
112  deriv += _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
113  _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
114  deriv += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
115  return _grad_test[_i][_qp] * deriv;
116 }
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)
bool notPorousFlowVariable(unsigned int moose_var_num) const
Returns true if moose_var_num is not a porous flow variabe.
const MaterialProperty< RealTensorValue > & _permeability
Permeability of porous material.
const MaterialProperty< std::vector< Real > > & _fluid_density_qp
Fluid density for each phase (at the qp)
const MaterialProperty< std::vector< RealGradient > > & _grad_p
Gradient of the pore pressure in each phase.
const MaterialProperty< std::vector< std::vector< RealGradient > > > & _dgrad_p_dvar
Derivative of Grad porepressure in each phase wrt PorousFlow variables.
const PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
const RealVectorValue _gravity
Gravity. Defaults to 9.81 m/s^2.
const MaterialProperty< std::vector< std::vector< Real > > > & _dgrad_p_dgrad_var
Derivative of Grad porepressure in each phase wrt grad(PorousFlow variables)
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
The PorousFlow variable number.

◆ dmobility()

Real PorousFlowDarcyBase::dmobility ( unsigned  nodenum,
unsigned  phase,
unsigned  pvar 
) const
protectedvirtual

The derivative of mobility with respect to PorousFlow variable pvar.

Parameters
nodenumThe node-number to evaluate the mobility for
phasethe fluid phase number
pvarthe PorousFlow variable pvar

Reimplemented in PorousFlowAdvectiveFlux, and PorousFlowHeatAdvection.

Definition at line 508 of file PorousFlowDarcyBase.C.

Referenced by fullyUpwind(), harmonicMean(), and quickUpwind().

509 {
510  Real dm = _dfluid_density_node_dvar[nodenum][phase][pvar] / _fluid_viscosity[nodenum][phase];
511  dm -= _fluid_density_node[nodenum][phase] * _dfluid_viscosity_dvar[nodenum][phase][pvar] /
512  std::pow(_fluid_viscosity[nodenum][phase], 2);
513  return dm;
514 }
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< Real > > & _fluid_viscosity
Viscosity of each component in each phase.
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_viscosity_dvar
Derivative of the fluid viscosity for each phase wrt PorousFlow variables.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const MaterialProperty< std::vector< Real > > & _fluid_density_node
Fluid density for each phase (at the node)

◆ fullyUpwind()

void PorousFlowDarcyBase::fullyUpwind ( JacRes  res_or_jac,
unsigned int  ph,
unsigned int  pvar 
)
protected

Calculate the residual or Jacobian using full upwinding.

Parameters
res_or_jacwhether to compute the residual or jacobian
phfluid phase number
pvardifferentiate wrt to this PorousFlow variable (when computing the jacobian)

Perform the full upwinding by multiplying the residuals at the upstream nodes by their mobilities. Mobility is different for each phase, and in each situation: mobility = density / viscosity for single-component Darcy flow mobility = mass_fraction * density * relative_perm / viscosity for multi-component, multiphase flow mobility = enthalpy * density * relative_perm / viscosity for heat convection

The residual for the kernel is the sum over Darcy fluxes for each phase. The Darcy flux for a particular phase is R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)} for node i. where int is the integral over the element. However, in fully-upwind, the first step is to take the mobility outside the integral, which was done in the _proto_flux calculation above.

NOTE: Physically _proto_flux[i][ph] is a measure of fluid of phase ph flowing out of node i. If we had left in mobility, it would be exactly the component mass flux flowing out of node i.

This leads to the definition of upwinding:

If _proto_flux(i)[ph] is positive then we use mobility_i. That is we use the upwind value of mobility.

The final subtle thing is we must also conserve fluid mass: the total component mass flowing out of node i must be the sum of the masses flowing into the other nodes.

note the -= means the result is positive

Definition at line 298 of file PorousFlowDarcyBase.C.

Referenced by computeResidualAndJacobian().

299 {
329  // The number of nodes in the element
330  const unsigned int num_nodes = _test.size();
331 
332  Real mob;
333  Real dmob;
334  // Define variables used to ensure mass conservation
335  Real total_mass_out = 0.0;
336  Real total_in = 0.0;
337 
338  // The following holds derivatives of these
339  std::vector<Real> dtotal_mass_out;
340  std::vector<Real> dtotal_in;
341  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
342  {
343  dtotal_mass_out.assign(num_nodes, 0.0);
344  dtotal_in.assign(num_nodes, 0.0);
345  }
346 
347  // Perform the upwinding using the mobility
348  std::vector<bool> upwind_node(num_nodes);
349  for (unsigned int n = 0; n < num_nodes; ++n)
350  {
351  if (_proto_flux[ph][n] >= 0.0) // upstream node
352  {
353  upwind_node[n] = true;
354  // The mobility at the upstream node
355  mob = mobility(n, ph);
356  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
357  {
358  // The derivative of the mobility wrt the PorousFlow variable
359  dmob = dmobility(n, ph, pvar);
360 
361  for (_j = 0; _j < _phi.size(); _j++)
362  _jacobian[ph][n][_j] *= mob;
363 
364  if (_test.size() == _phi.size())
365  /* mobility at node=n depends only on the variables at node=n, by construction. For
366  * linear-lagrange variables, this means that Jacobian entries involving the derivative
367  * of mobility will only be nonzero for derivatives wrt variables at node=n. Hence the
368  * [n][n] in the line below. However, for other variable types (eg constant monomials)
369  * I cannot tell what variable number contributes to the derivative. However, in all
370  * cases I can possibly imagine, the derivative is zero anyway, since in the full
371  * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
372  */
373  _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
374 
375  for (_j = 0; _j < _phi.size(); _j++)
376  dtotal_mass_out[_j] += _jacobian[ph][n][_j];
377  }
378  _proto_flux[ph][n] *= mob;
379  total_mass_out += _proto_flux[ph][n];
380  }
381  else
382  {
383  upwind_node[n] = false;
384  total_in -= _proto_flux[ph][n];
385  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
386  for (_j = 0; _j < _phi.size(); _j++)
387  dtotal_in[_j] -= _jacobian[ph][n][_j];
388  }
389  }
390 
391  // Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes,
392  // weighted by their proto_flux values
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 }
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.
virtual Real mobility(unsigned nodenum, unsigned phase) const
The mobility of the fluid.
virtual Real dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
The derivative of mobility with respect to PorousFlow variable pvar.

◆ harmonicMean()

void PorousFlowDarcyBase::harmonicMean ( JacRes  res_or_jac,
unsigned int  ph,
unsigned int  pvar 
)
protected

Calculate the residual or Jacobian by using the harmonic mean of the nodal mobilities for the entire element.

Parameters
res_or_jacwhether to compute the residual or jacobian
phfluid phase number
pvardifferentiate wrt to this PorousFlow variable (when computing the jacobian)

Definition at line 448 of file PorousFlowDarcyBase.C.

Referenced by computeResidualAndJacobian().

449 {
450  // The number of nodes in the element
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  // d(harmonic_mob)/d(PorousFlow variable at node n)
474  std::vector<Real> dharmonic_mob(num_nodes, 0.0);
475  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
476  {
477  const Real harm2 = std::pow(harmonic_mob, 2) / (1.0 * num_nodes);
478  if (num_zero == 0)
479  for (unsigned n = 0; n < num_nodes; ++n)
480  dharmonic_mob[n] = dmobility(n, ph, pvar) * harm2 / std::pow(mob[n], 2);
481  else if (num_zero == 1)
482  dharmonic_mob[zero_mobility_node] =
483  num_nodes * dmobility(zero_mobility_node, ph, pvar); // other derivs are zero
484  // if num_zero > 1 then all dharmonic_mob = 0.0
485  }
486 
487  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
488  for (unsigned n = 0; n < num_nodes; ++n)
489  for (_j = 0; _j < _phi.size(); _j++)
490  {
491  _jacobian[ph][n][_j] *= harmonic_mob;
492  if (_test.size() == _phi.size())
493  _jacobian[ph][n][_j] += dharmonic_mob[_j] * _proto_flux[ph][n];
494  }
495 
496  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
497  for (unsigned n = 0; n < num_nodes; ++n)
498  _proto_flux[ph][n] *= harmonic_mob;
499 }
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.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
virtual Real mobility(unsigned nodenum, unsigned phase) const
The mobility of the fluid.
virtual Real dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
The derivative of mobility with respect to PorousFlow variable pvar.

◆ mobility()

Real PorousFlowDarcyBase::mobility ( unsigned  nodenum,
unsigned  phase 
) const
protectedvirtual

The mobility of the fluid.

For multi-component Darcy flow this is mass_fraction * fluid_density * relative_permeability / fluid_viscosity

Parameters
nodenumThe node-number to evaluate the mobility for
phasethe fluid phase number

Reimplemented in PorousFlowAdvectiveFlux, and PorousFlowHeatAdvection.

Definition at line 502 of file PorousFlowDarcyBase.C.

Referenced by fullyUpwind(), harmonicMean(), and quickUpwind().

503 {
504  return _fluid_density_node[nodenum][phase] / _fluid_viscosity[nodenum][phase];
505 }
const MaterialProperty< std::vector< Real > > & _fluid_viscosity
Viscosity of each component in each phase.
const MaterialProperty< std::vector< Real > > & _fluid_density_node
Fluid density for each phase (at the node)

◆ quickUpwind()

void PorousFlowDarcyBase::quickUpwind ( JacRes  res_or_jac,
unsigned int  ph,
unsigned int  pvar 
)
protected

Calculate the residual or Jacobian using the nodal mobilities, but without conserving fluid mass.

Parameters
res_or_jacwhether to compute the residual or jacobian
phfluid phase number
pvardifferentiate wrt to this PorousFlow variable (when computing the jacobian)

Definition at line 411 of file PorousFlowDarcyBase.C.

Referenced by computeResidualAndJacobian().

412 {
413  // The number of nodes in the element
414  const unsigned int num_nodes = _test.size();
415 
416  Real mob;
417  Real dmob;
418 
419  // Use the raw nodal mobility
420  for (unsigned int n = 0; n < num_nodes; ++n)
421  {
422  // The mobility at the node
423  mob = mobility(n, ph);
424  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
425  {
426  // The derivative of the mobility wrt the PorousFlow variable
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 }
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.
virtual Real mobility(unsigned nodenum, unsigned phase) const
The mobility of the fluid.
virtual Real dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
The derivative of mobility with respect to PorousFlow variable pvar.

◆ timestepSetup()

void PorousFlowDarcyBase::timestepSetup ( )
overrideprotectedvirtual

Definition at line 86 of file PorousFlowDarcyBase.C.

87 {
88  Kernel::timestepSetup();
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 }
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::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 ...

Member Data Documentation

◆ _dfluid_density_node_dvar

const MaterialProperty<std::vector<std::vector<Real> > >& PorousFlowDarcyBase::_dfluid_density_node_dvar
protected

Derivative of the fluid density for each phase wrt PorousFlow variables (at the node)

Definition at line 100 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowAdvectiveFlux::dmobility(), PorousFlowHeatAdvection::dmobility(), and dmobility().

◆ _dfluid_density_qp_dvar

const MaterialProperty<std::vector<std::vector<Real> > >& PorousFlowDarcyBase::_dfluid_density_qp_dvar
protected

Derivative of the fluid density for each phase wrt PorousFlow variables (at the qp)

Definition at line 106 of file PorousFlowDarcyBase.h.

Referenced by darcyQpJacobian().

◆ _dfluid_viscosity_dvar

const MaterialProperty<std::vector<std::vector<Real> > >& PorousFlowDarcyBase::_dfluid_viscosity_dvar
protected

Derivative of the fluid viscosity for each phase wrt PorousFlow variables.

Definition at line 112 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowAdvectiveFlux::dmobility(), PorousFlowHeatAdvection::dmobility(), and dmobility().

◆ _dgrad_p_dgrad_var

const MaterialProperty<std::vector<std::vector<Real> > >& PorousFlowDarcyBase::_dgrad_p_dgrad_var
protected

Derivative of Grad porepressure in each phase wrt grad(PorousFlow variables)

Definition at line 121 of file PorousFlowDarcyBase.h.

Referenced by darcyQpJacobian().

◆ _dgrad_p_dvar

const MaterialProperty<std::vector<std::vector<RealGradient> > >& PorousFlowDarcyBase::_dgrad_p_dvar
protected

Derivative of Grad porepressure in each phase wrt PorousFlow variables.

Definition at line 124 of file PorousFlowDarcyBase.h.

Referenced by darcyQpJacobian().

◆ _dictator

const PorousFlowDictator& PorousFlowDarcyBase::_dictator
protected

PorousFlowDictator UserObject.

Definition at line 127 of file PorousFlowDarcyBase.h.

Referenced by computeResidualAndJacobian(), and darcyQpJacobian().

◆ _dpermeability_dgradvar

const MaterialProperty<std::vector<std::vector<RealTensorValue> > >& PorousFlowDarcyBase::_dpermeability_dgradvar
protected

d(permeabiity)/d(grad(PorousFlow variable))

Definition at line 94 of file PorousFlowDarcyBase.h.

Referenced by darcyQpJacobian().

◆ _dpermeability_dvar

const MaterialProperty<std::vector<RealTensorValue> >& PorousFlowDarcyBase::_dpermeability_dvar
protected

d(permeabiity)/d(PorousFlow variable)

Definition at line 91 of file PorousFlowDarcyBase.h.

Referenced by darcyQpJacobian().

◆ _fallback_scheme

enum PorousFlowDarcyBase::FallbackEnum PorousFlowDarcyBase::_fallback_scheme
protected

◆ _fluid_density_node

const MaterialProperty<std::vector<Real> >& PorousFlowDarcyBase::_fluid_density_node
protected

◆ _fluid_density_qp

const MaterialProperty<std::vector<Real> >& PorousFlowDarcyBase::_fluid_density_qp
protected

Fluid density for each phase (at the qp)

Definition at line 103 of file PorousFlowDarcyBase.h.

Referenced by darcyQp(), and darcyQpJacobian().

◆ _fluid_viscosity

const MaterialProperty<std::vector<Real> >& PorousFlowDarcyBase::_fluid_viscosity
protected

◆ _full_upwind_threshold

const unsigned PorousFlowDarcyBase::_full_upwind_threshold
protected

If the number of upwind-downwind swaps is less than this amount then full upwinding is used.

Otherwise the fallback scheme is employed

Definition at line 139 of file PorousFlowDarcyBase.h.

Referenced by computeResidualAndJacobian().

◆ _grad_p

const MaterialProperty<std::vector<RealGradient> >& PorousFlowDarcyBase::_grad_p
protected

Gradient of the pore pressure in each phase.

Definition at line 118 of file PorousFlowDarcyBase.h.

Referenced by darcyQp(), and darcyQpJacobian().

◆ _gravity

const RealVectorValue PorousFlowDarcyBase::_gravity
protected

Gravity. Defaults to 9.81 m/s^2.

Definition at line 133 of file PorousFlowDarcyBase.h.

Referenced by darcyQp(), and darcyQpJacobian().

◆ _jacobian

std::vector<std::vector<std::vector<Real> > > PorousFlowDarcyBase::_jacobian
protected

Derivative of _proto_flux with respect to nodal variables.

This gets modified with the mobility and its derivative during computation of MOOSE's Jacobian

Definition at line 165 of file PorousFlowDarcyBase.h.

Referenced by computeResidualAndJacobian(), fullyUpwind(), harmonicMean(), and quickUpwind().

◆ _num_downwinds

std::unordered_map<unsigned, std::vector<std::vector<unsigned> > > PorousFlowDarcyBase::_num_downwinds
protected

Number of nonlinear iterations (in this timestep and this element) that a node is an downwind node for a given fluid phase.

_num_upwinds[element_number][phase][node_number_in_element]

Definition at line 179 of file PorousFlowDarcyBase.h.

Referenced by computeResidualAndJacobian(), and timestepSetup().

◆ _num_phases

const unsigned int PorousFlowDarcyBase::_num_phases
protected

The number of fluid phases.

Definition at line 130 of file PorousFlowDarcyBase.h.

Referenced by computeResidualAndJacobian().

◆ _num_upwinds

std::unordered_map<unsigned, std::vector<std::vector<unsigned> > > PorousFlowDarcyBase::_num_upwinds
protected

Number of nonlinear iterations (in this timestep and this element) that a node is an upwind node for a given fluid phase.

_num_upwinds[element_number][phase][node_number_in_element]

Definition at line 172 of file PorousFlowDarcyBase.h.

Referenced by computeResidualAndJacobian(), and timestepSetup().

◆ _permeability

const MaterialProperty<RealTensorValue>& PorousFlowDarcyBase::_permeability
protected

Permeability of porous material.

Definition at line 88 of file PorousFlowDarcyBase.h.

Referenced by darcyQp(), and darcyQpJacobian().

◆ _pp

const MaterialProperty<std::vector<Real> >& PorousFlowDarcyBase::_pp
protected

Nodal pore pressure in each phase.

Definition at line 115 of file PorousFlowDarcyBase.h.

◆ _proto_flux

std::vector<std::vector<Real> > PorousFlowDarcyBase::_proto_flux
protected

The Darcy flux.

When multiplied by the mobility, this is the mass flux of fluid moving out of a node. This multiplication occurs during the formation of the residual and Jacobian. _proto_flux[num_phases][num_nodes]

Definition at line 158 of file PorousFlowDarcyBase.h.

Referenced by computeResidualAndJacobian(), fullyUpwind(), harmonicMean(), and quickUpwind().


The documentation for this class was generated from the following files: