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

Convective flux of component k in fluid phase alpha. More...

#include <PorousFlowAdvectiveFlux.h>

Inheritance diagram for PorousFlowAdvectiveFlux:
[legend]

Public Member Functions

 PorousFlowAdvectiveFlux (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 Real mobility (unsigned nodenum, unsigned phase) const override
 The mobility of the fluid. More...
 
virtual Real dmobility (unsigned nodenum, unsigned phase, unsigned pvar) const override
 The derivative of mobility with respect to PorousFlow variable pvar. More...
 
virtual void timestepSetup () override
 
virtual Real computeQpResidual () override
 
virtual void computeResidual () override
 
virtual void computeJacobian () override
 
virtual void computeOffDiagJacobian (unsigned int 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...
 
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< std::vector< std::vector< Real > > > & _mass_fractions
 Mass fraction of each component in each phase. More...
 
const MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _dmass_fractions_dvar
 Derivative of the mass fraction of each component in each phase wrt PorousFlow variables. More...
 
const MaterialProperty< std::vector< Real > > & _relative_permeability
 Relative permeability of each phase. More...
 
const MaterialProperty< std::vector< std::vector< Real > > > & _drelative_permeability_dvar
 Derivative of relative permeability of each phase wrt PorousFlow variables. More...
 
const unsigned int _fluid_component
 Index of the fluid component that this kernel acts on. More...
 
const MaterialProperty< RealTensorValue > & _permeability
 Permeability of porous material. More...
 
const MaterialProperty< std::vector< RealTensorValue > > & _dpermeability_dvar
 d(permeabiity)/d(porous-flow variable) More...
 
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dpermeability_dgradvar
 d(permeabiity)/d(grad(porous-flow 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_porousflow_dictator
 PorousFlow 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

Convective flux of component k in fluid phase alpha.

A fully-updwinded version is implemented, where the mobility of the upstream nodes is used.

Definition at line 23 of file PorousFlowAdvectiveFlux.h.

Member Enumeration Documentation

enum PorousFlowDarcyBase::FallbackEnum
strongprotectedinherited

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 145 of file PorousFlowDarcyBase.h.

145 { QUICK, HARMONIC } _fallback_scheme;
enum PorousFlowDarcyBase::FallbackEnum _fallback_scheme
enum PorousFlowDarcyBase::JacRes
strongprotectedinherited
Enumerator
CALCULATE_RESIDUAL 
CALCULATE_JACOBIAN 

Definition at line 59 of file PorousFlowDarcyBase.h.

60  {
61  CALCULATE_RESIDUAL = 0,
62  CALCULATE_JACOBIAN = 1
63  };

Constructor & Destructor Documentation

PorousFlowAdvectiveFlux::PorousFlowAdvectiveFlux ( const InputParameters &  parameters)

Definition at line 22 of file PorousFlowAdvectiveFlux.C.

23  : PorousFlowDarcyBase(parameters),
25  getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
26  _dmass_fractions_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
27  "dPorousFlow_mass_frac_nodal_dvar")),
29  getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal")),
30  _drelative_permeability_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
31  "dPorousFlow_relative_permeability_nodal_dvar")),
32  _fluid_component(getParam<unsigned int>("fluid_component"))
33 {
34 }
const unsigned int _fluid_component
Index of the fluid component that this kernel acts on.
const MaterialProperty< std::vector< Real > > & _relative_permeability
Relative permeability of each phase.
const MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _dmass_fractions_dvar
Derivative of the mass fraction of each component in each phase wrt PorousFlow variables.
const MaterialProperty< std::vector< std::vector< Real > > > & _mass_fractions
Mass fraction of each component in each phase.
const MaterialProperty< std::vector< std::vector< Real > > > & _drelative_permeability_dvar
Derivative of relative permeability of each phase wrt PorousFlow variables.
PorousFlowDarcyBase(const InputParameters &parameters)

Member Function Documentation

void PorousFlowDarcyBase::computeJacobian ( )
overrideprotectedvirtualinherited

Definition at line 131 of file PorousFlowDarcyBase.C.

132 {
134 }
void computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
Computation of the residual and Jacobian.
void PorousFlowDarcyBase::computeOffDiagJacobian ( unsigned int  jvar)
overrideprotectedvirtualinherited

Definition at line 137 of file PorousFlowDarcyBase.C.

138 {
140 }
void computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
Computation of the residual and Jacobian.
Real PorousFlowDarcyBase::computeQpResidual ( )
overrideprotectedvirtualinherited

Definition at line 118 of file PorousFlowDarcyBase.C.

119 {
120  mooseError("PorousFlowDarcyBase: computeQpResidual called");
121  return 0.0;
122 }
void PorousFlowDarcyBase::computeResidual ( )
overrideprotectedvirtualinherited

Definition at line 125 of file PorousFlowDarcyBase.C.

126 {
128 }
void computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
Computation of the residual and Jacobian.
void PorousFlowDarcyBase::computeResidualAndJacobian ( JacRes  res_or_jac,
unsigned int  jvar 
)
protectedinherited

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

The PorousFlow variable index corresponding to the variable number jvar

The number of nodes in the element

Compute the residual and jacobian without the mobility terms. Even if we are computing the Jacobian we still need this in order to see which nodes are upwind and which are downwind.

Add results to the Residual or Jacobian

Definition at line 143 of file PorousFlowDarcyBase.C.

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

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 }
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.
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 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...
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
void fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
Calculate the residual or Jacobian using full upwinding.
Real PorousFlowDarcyBase::darcyQp ( unsigned int  ph) const
protectedvirtualinherited

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

Definition at line 93 of file PorousFlowDarcyBase.C.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian().

94 {
95  return _grad_test[_i][_qp] *
96  (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
97 }
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.
Real PorousFlowDarcyBase::darcyQpJacobian ( unsigned int  jvar,
unsigned int  ph 
) const
protectedvirtualinherited

Jacobian of the Darcy part of the flux.

Definition at line 100 of file PorousFlowDarcyBase.C.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian().

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 }
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)
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 PorousFlowDictator & _porousflow_dictator
PorousFlow UserObject.
const MaterialProperty< std::vector< std::vector< RealGradient > > > & _dgrad_p_dvar
Derivative of Grad porepressure in each phase wrt PorousFlow variables.
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)
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
Real PorousFlowAdvectiveFlux::dmobility ( unsigned  nodenum,
unsigned  phase,
unsigned  pvar 
) const
overrideprotectedvirtual

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 from PorousFlowDarcyBase.

Definition at line 44 of file PorousFlowAdvectiveFlux.C.

45 {
46  Real dm = _dmass_fractions_dvar[nodenum][phase][_fluid_component][pvar] *
47  _fluid_density_node[nodenum][phase] * _relative_permeability[nodenum][phase] /
48  _fluid_viscosity[nodenum][phase];
49  dm += _mass_fractions[nodenum][phase][_fluid_component] *
50  _dfluid_density_node_dvar[nodenum][phase][pvar] * _relative_permeability[nodenum][phase] /
51  _fluid_viscosity[nodenum][phase];
52  dm += _mass_fractions[nodenum][phase][_fluid_component] * _fluid_density_node[nodenum][phase] *
53  _drelative_permeability_dvar[nodenum][phase][pvar] / _fluid_viscosity[nodenum][phase];
54  dm -= _mass_fractions[nodenum][phase][_fluid_component] * _fluid_density_node[nodenum][phase] *
55  _relative_permeability[nodenum][phase] * _dfluid_viscosity_dvar[nodenum][phase][pvar] /
56  std::pow(_fluid_viscosity[nodenum][phase], 2);
57  return dm;
58 }
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 unsigned int _fluid_component
Index of the fluid component that this kernel acts on.
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.
const MaterialProperty< std::vector< Real > > & _relative_permeability
Relative permeability of each phase.
const MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _dmass_fractions_dvar
Derivative of the mass fraction of each component in each phase wrt PorousFlow variables.
const MaterialProperty< std::vector< std::vector< Real > > > & _mass_fractions
Mass fraction of each component in each phase.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const MaterialProperty< std::vector< std::vector< Real > > > & _drelative_permeability_dvar
Derivative of relative permeability of each phase wrt PorousFlow variables.
const MaterialProperty< std::vector< Real > > & _fluid_density_node
Fluid density for each phase (at the node)
void PorousFlowDarcyBase::fullyUpwind ( JacRes  res_or_jac,
unsigned int  ph,
unsigned int  pvar 
)
protectedinherited

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.

The number of nodes in the element

Define variables used to ensure mass conservation

The following holds derivatives of these

Perform the upwinding using the mobility

The mobility at the upstream node

The derivative of the mobility wrt the PorousFlow variable

note the -= means the result is positive

Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes, weighted by their proto_flux values

Definition at line 299 of file PorousFlowDarcyBase.C.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian().

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 }
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 dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
The derivative of mobility with respect to PorousFlow variable pvar.
virtual Real mobility(unsigned nodenum, unsigned phase) const
The mobility of the fluid.
void PorousFlowDarcyBase::harmonicMean ( JacRes  res_or_jac,
unsigned int  ph,
unsigned int  pvar 
)
protectedinherited

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)

The number of nodes in the element

Definition at line 448 of file PorousFlowDarcyBase.C.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian().

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 }
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 dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
The derivative of mobility with respect to PorousFlow variable pvar.
virtual Real mobility(unsigned nodenum, unsigned phase) const
The mobility of the fluid.
Real PorousFlowAdvectiveFlux::mobility ( unsigned  nodenum,
unsigned  phase 
) const
overrideprotectedvirtual

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 from PorousFlowDarcyBase.

Definition at line 37 of file PorousFlowAdvectiveFlux.C.

38 {
39  return _mass_fractions[nodenum][phase][_fluid_component] * _fluid_density_node[nodenum][phase] *
40  _relative_permeability[nodenum][phase] / _fluid_viscosity[nodenum][phase];
41 }
const unsigned int _fluid_component
Index of the fluid component that this kernel acts on.
const MaterialProperty< std::vector< Real > > & _fluid_viscosity
Viscosity of each component in each phase.
const MaterialProperty< std::vector< Real > > & _relative_permeability
Relative permeability of each phase.
const MaterialProperty< std::vector< std::vector< Real > > > & _mass_fractions
Mass fraction of each component in each phase.
const MaterialProperty< std::vector< Real > > & _fluid_density_node
Fluid density for each phase (at the node)
void PorousFlowDarcyBase::quickUpwind ( JacRes  res_or_jac,
unsigned int  ph,
unsigned int  pvar 
)
protectedinherited

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)

The number of nodes in the element

Use the raw nodal mobility

The mobility at the node

The derivative of the mobility wrt the PorousFlow variable

Definition at line 411 of file PorousFlowDarcyBase.C.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian().

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 }
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 dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
The derivative of mobility with respect to PorousFlow variable pvar.
virtual Real mobility(unsigned nodenum, unsigned phase) const
The mobility of the fluid.
void PorousFlowDarcyBase::timestepSetup ( )
overrideprotectedvirtualinherited

Definition at line 85 of file PorousFlowDarcyBase.C.

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 }
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

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

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

Definition at line 95 of file PorousFlowDarcyBase.h.

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

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

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

Definition at line 101 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::darcyQpJacobian().

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

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

Definition at line 107 of file PorousFlowDarcyBase.h.

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

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

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

Definition at line 116 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::darcyQpJacobian().

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

Derivative of Grad porepressure in each phase wrt PorousFlow variables.

Definition at line 119 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::darcyQpJacobian().

const MaterialProperty<std::vector<std::vector<std::vector<Real> > > >& PorousFlowAdvectiveFlux::_dmass_fractions_dvar
protected

Derivative of the mass fraction of each component in each phase wrt PorousFlow variables.

Definition at line 36 of file PorousFlowAdvectiveFlux.h.

Referenced by dmobility().

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

d(permeabiity)/d(grad(porous-flow variable))

Definition at line 89 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::darcyQpJacobian().

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

d(permeabiity)/d(porous-flow variable)

Definition at line 86 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::darcyQpJacobian().

const MaterialProperty<std::vector<std::vector<Real> > >& PorousFlowAdvectiveFlux::_drelative_permeability_dvar
protected

Derivative of relative permeability of each phase wrt PorousFlow variables.

Definition at line 42 of file PorousFlowAdvectiveFlux.h.

Referenced by dmobility().

enum PorousFlowDarcyBase::FallbackEnum PorousFlowDarcyBase::_fallback_scheme
protectedinherited
const unsigned int PorousFlowAdvectiveFlux::_fluid_component
protected

Index of the fluid component that this kernel acts on.

Definition at line 45 of file PorousFlowAdvectiveFlux.h.

Referenced by dmobility(), and mobility().

const MaterialProperty<std::vector<Real> >& PorousFlowDarcyBase::_fluid_density_node
protectedinherited
const MaterialProperty<std::vector<Real> >& PorousFlowDarcyBase::_fluid_density_qp
protectedinherited

Fluid density for each phase (at the qp)

Definition at line 98 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::darcyQp(), and PorousFlowDarcyBase::darcyQpJacobian().

const MaterialProperty<std::vector<Real> >& PorousFlowDarcyBase::_fluid_viscosity
protectedinherited
const unsigned PorousFlowDarcyBase::_full_upwind_threshold
protectedinherited

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 134 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian().

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

Gradient of the pore pressure in each phase.

Definition at line 113 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::darcyQp(), and PorousFlowDarcyBase::darcyQpJacobian().

const RealVectorValue PorousFlowDarcyBase::_gravity
protectedinherited

Gravity. Defaults to 9.81 m/s^2.

Definition at line 128 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::darcyQp(), and PorousFlowDarcyBase::darcyQpJacobian().

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

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 160 of file PorousFlowDarcyBase.h.

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

const MaterialProperty<std::vector<std::vector<Real> > >& PorousFlowAdvectiveFlux::_mass_fractions
protected

Mass fraction of each component in each phase.

Definition at line 33 of file PorousFlowAdvectiveFlux.h.

Referenced by dmobility(), and mobility().

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

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 174 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian(), and PorousFlowDarcyBase::timestepSetup().

const unsigned int PorousFlowDarcyBase::_num_phases
protectedinherited

The number of fluid phases.

Definition at line 125 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian().

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

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 167 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian(), and PorousFlowDarcyBase::timestepSetup().

const MaterialProperty<RealTensorValue>& PorousFlowDarcyBase::_permeability
protectedinherited

Permeability of porous material.

Definition at line 83 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::darcyQp(), and PorousFlowDarcyBase::darcyQpJacobian().

const PorousFlowDictator& PorousFlowDarcyBase::_porousflow_dictator
protectedinherited

PorousFlow UserObject.

Definition at line 122 of file PorousFlowDarcyBase.h.

Referenced by PorousFlowDarcyBase::computeResidualAndJacobian(), and PorousFlowDarcyBase::darcyQpJacobian().

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

Nodal pore pressure in each phase.

Definition at line 110 of file PorousFlowDarcyBase.h.

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

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 153 of file PorousFlowDarcyBase.h.

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

const MaterialProperty<std::vector<Real> >& PorousFlowAdvectiveFlux::_relative_permeability
protected

Relative permeability of each phase.

Definition at line 39 of file PorousFlowAdvectiveFlux.h.

Referenced by dmobility(), and mobility().


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