Darcy advective flux. More...

#include <PorousFlowDarcyBase.h>

Inheritance diagram for PorousFlowDarcyBase:

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 bool	_perm_derivs
	Flag to check whether permeabiity derivatives are non-zero. 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 28 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 152 of file PorousFlowDarcyBase.h.

152 { QUICK, HARMONIC } _fallback_scheme;

◆ JacRes

enum PorousFlowDarcyBase::JacRes

strongprotected

Enumerator
CALCULATE_RESIDUAL
CALCULATE_JACOBIAN

Definition at line 63 of file PorousFlowDarcyBase.h.

   {
     CALCULATE_RESIDUAL = 0,
     CALCULATE_JACOBIAN = 1
   };

Constructor & Destructor Documentation

◆ PorousFlowDarcyBase()

PorousFlowDarcyBase::PorousFlowDarcyBase ( const InputParameters & parameters )

Definition at line 45 of file PorousFlowDarcyBase.C.

   : Kernel(parameters),
     _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
     _dpermeability_dvar(
         getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
     _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
         "dPorousFlow_permeability_qp_dgradvar")),
     _fluid_density_node(
         getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
     _dfluid_density_node_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
         "dPorousFlow_fluid_phase_density_nodal_dvar")),
     _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
     _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
         "dPorousFlow_fluid_phase_density_qp_dvar")),
     _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")),
     _dfluid_viscosity_dvar(
         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
     _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_nodal")),
     _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
     _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
         "dPorousFlow_grad_porepressure_qp_dgradvar")),
     _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
         "dPorousFlow_grad_porepressure_qp_dvar")),
     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
     _num_phases(_dictator.numPhases()),
     _gravity(getParam<RealVectorValue>("gravity")),
     _perm_derivs(_dictator.usePermDerivs()),
     _full_upwind_threshold(getParam<unsigned>("full_upwind_threshold")),
     _fallback_scheme(getParam<MooseEnum>("fallback_scheme").getEnum<FallbackEnum>()),
     _proto_flux(_num_phases),
     _jacobian(_num_phases),
     _num_upwinds(),
     _num_downwinds()
 {
 #ifdef LIBMESH_HAVE_TBB_API
   if (libMesh::n_threads() > 1)
     mooseWarning("PorousFlowDarcyBase: num_upwinds and num_downwinds may not be computed "
                  "accurately when using TBB and greater than 1 thread");
 #endif
 }

Member Function Documentation

◆ computeJacobian()

void PorousFlowDarcyBase::computeJacobian ( )

overrideprotectedvirtual

Definition at line 141 of file PorousFlowDarcyBase.C.

 {
   computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, _var.number());
 }

◆ computeOffDiagJacobian()

void PorousFlowDarcyBase::computeOffDiagJacobian ( MooseVariableFEBase & jvar )

overrideprotectedvirtual

Definition at line 147 of file PorousFlowDarcyBase.C.

 {
   computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, jvar.number());
 }

◆ computeQpResidual()

Real PorousFlowDarcyBase::computeQpResidual ( )

overrideprotectedvirtual

Definition at line 128 of file PorousFlowDarcyBase.C.

 {
   mooseError("PorousFlowDarcyBase: computeQpResidual called");
   return 0.0;
 }

◆ computeResidual()

void PorousFlowDarcyBase::computeResidual ( )

overrideprotectedvirtual

Definition at line 135 of file PorousFlowDarcyBase.C.

 {
   computeResidualAndJacobian(JacRes::CALCULATE_RESIDUAL, 0.0);
 }

◆ 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_jac	Whether to calculate the residual or the jacobian.
jvar	PorousFlow variable to take derivatives wrt to

Definition at line 153 of file PorousFlowDarcyBase.C.

 {
   if ((res_or_jac == JacRes::CALCULATE_JACOBIAN) && _dictator.notPorousFlowVariable(jvar))
     return;
  
   // The PorousFlow variable index corresponding to the variable number jvar
   const unsigned int pvar =
       ((res_or_jac == JacRes::CALCULATE_JACOBIAN) ? _dictator.porousFlowVariableNum(jvar) : 0);
  
   DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
   if ((ke.n() == 0) && (res_or_jac == JacRes::CALCULATE_JACOBIAN)) // this removes a problem
                                                                    // encountered in the initial
                                                                    // timestep when
                                                                    // use_displaced_mesh=true
     return;
  
   // The number of nodes in the element
   const unsigned int num_nodes = _test.size();
  
   // 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.
   for (unsigned ph = 0; ph < _num_phases; ++ph)
   {
     _proto_flux[ph].assign(num_nodes, 0);
     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     {
       for (_i = 0; _i < num_nodes; ++_i)
         _proto_flux[ph][_i] += _JxW[_qp] * _coord[_qp] * darcyQp(ph);
     }
   }
  
   // for this element, record whether each node is "upwind" or "downwind" (or neither)
   const unsigned elem = _current_elem->id();
   if (_num_upwinds.find(elem) == _num_upwinds.end())
   {
     _num_upwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
     _num_downwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
     for (unsigned ph = 0; ph < _num_phases; ++ph)
     {
       _num_upwinds[elem][ph].assign(num_nodes, 0);
       _num_downwinds[elem][ph].assign(num_nodes, 0);
     }
   }
   // record the information once per nonlinear iteration
   if (res_or_jac == JacRes::CALCULATE_JACOBIAN && jvar == _var.number())
   {
     for (unsigned ph = 0; ph < _num_phases; ++ph)
     {
       for (unsigned nod = 0; nod < num_nodes; ++nod)
       {
         if (_proto_flux[ph][nod] > 0)
           _num_upwinds[elem][ph][nod]++;
         else if (_proto_flux[ph][nod] < 0)
           _num_downwinds[elem][ph][nod]++;
       }
     }
   }
  
   // based on _num_upwinds and _num_downwinds, calculate the maximum number
   // of upwind-downwind swaps that have been encountered in this timestep
   // for this element
   std::vector<unsigned> max_swaps(_num_phases, 0);
   for (unsigned ph = 0; ph < _num_phases; ++ph)
   {
     for (unsigned nod = 0; nod < num_nodes; ++nod)
       max_swaps[ph] = std::max(
           max_swaps[ph], std::min(_num_upwinds[elem][ph][nod], _num_downwinds[elem][ph][nod]));
   }
  
   // size the _jacobian correctly and calculate it for the case residual = _proto_flux
   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
   {
     for (unsigned ph = 0; ph < _num_phases; ++ph)
     {
       _jacobian[ph].resize(ke.m());
       for (_i = 0; _i < _test.size(); _i++)
       {
         _jacobian[ph][_i].assign(ke.n(), 0.0);
         for (_j = 0; _j < _phi.size(); _j++)
           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
             _jacobian[ph][_i][_j] += _JxW[_qp] * _coord[_qp] * darcyQpJacobian(jvar, ph);
       }
     }
   }
  
   // Loop over all the phases, computing the mass flux, which
   // gets placed into _proto_flux, and the derivative of this
   // which gets placed into _jacobian
   for (unsigned int ph = 0; ph < _num_phases; ++ph)
   {
     if (max_swaps[ph] < _full_upwind_threshold)
       fullyUpwind(res_or_jac, ph, pvar);
     else
     {
       switch (_fallback_scheme)
       {
         case FallbackEnum::QUICK:
           quickUpwind(res_or_jac, ph, pvar);
           break;
         case FallbackEnum::HARMONIC:
           harmonicMean(res_or_jac, ph, pvar);
           break;
       }
     }
   }
  
   // Add results to the Residual or Jacobian
   if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
   {
     DenseVector<Number> & re = _assembly.residualBlock(_var.number());
  
     _local_re.resize(re.size());
     _local_re.zero();
     for (_i = 0; _i < _test.size(); _i++)
       for (unsigned int ph = 0; ph < _num_phases; ++ph)
         _local_re(_i) += _proto_flux[ph][_i];
  
     re += _local_re;
  
     if (_has_save_in)
     {
       Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
       for (unsigned int i = 0; i < _save_in.size(); i++)
         _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
     }
   }
  
   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
   {
     _local_ke.resize(ke.m(), ke.n());
     _local_ke.zero();
  
     for (_i = 0; _i < _test.size(); _i++)
       for (_j = 0; _j < _phi.size(); _j++)
         for (unsigned int ph = 0; ph < _num_phases; ++ph)
           _local_ke(_i, _j) += _jacobian[ph][_i][_j];
  
     ke += _local_ke;
  
     if (_has_diag_save_in && jvar == _var.number())
     {
       unsigned int rows = ke.m();
       DenseVector<Number> diag(rows);
       for (unsigned int i = 0; i < rows; i++)
         diag(i) = _local_ke(i, i);
  
       Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
       for (unsigned int i = 0; i < _diag_save_in.size(); i++)
         _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
     }
   }
 }

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

◆ darcyQp()

Real PorousFlowDarcyBase::darcyQp ( unsigned int ph ) const

protectedvirtual

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

Definition at line 95 of file PorousFlowDarcyBase.C.

 {
   return _grad_test[_i][_qp] *
          (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
 }

Referenced by computeResidualAndJacobian().

◆ darcyQpJacobian()

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

protectedvirtual

Jacobian of the Darcy part of the flux.

Definition at line 102 of file PorousFlowDarcyBase.C.

 {
   if (_dictator.notPorousFlowVariable(jvar))
     return 0.0;
  
   const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
  
   RealVectorValue deriv =
       _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
                             _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
  
   deriv += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
  
   if (_perm_derivs)
   {
     deriv += _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
              (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
     for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
       deriv += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
                (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
   }
  
   return _grad_test[_i][_qp] * deriv;
 }

Referenced by computeResidualAndJacobian().

◆ dmobility()

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

protectedvirtual

The derivative of mobility with respect to PorousFlow variable pvar.

Parameters

nodenum	The node-number to evaluate the mobility for
phase	the fluid phase number
pvar	the PorousFlow variable pvar

Reimplemented in PorousFlowAdvectiveFlux, and PorousFlowHeatAdvection.

Definition at line 517 of file PorousFlowDarcyBase.C.

 {
   Real dm = _dfluid_density_node_dvar[nodenum][phase][pvar] / _fluid_viscosity[nodenum][phase];
   dm -= _fluid_density_node[nodenum][phase] * _dfluid_viscosity_dvar[nodenum][phase][pvar] /
         std::pow(_fluid_viscosity[nodenum][phase], 2);
   return dm;
 }

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

◆ 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_jac	whether to compute the residual or jacobian
ph	fluid phase number
pvar	differentiate 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 307 of file PorousFlowDarcyBase.C.

 {
   // The number of nodes in the element
   const unsigned int num_nodes = _test.size();
  
   Real mob;
   Real dmob;
   // Define variables used to ensure mass conservation
   Real total_mass_out = 0.0;
   Real total_in = 0.0;
  
   // The following holds derivatives of these
   std::vector<Real> dtotal_mass_out;
   std::vector<Real> dtotal_in;
   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
   {
     dtotal_mass_out.assign(num_nodes, 0.0);
     dtotal_in.assign(num_nodes, 0.0);
   }
  
   // Perform the upwinding using the mobility
   std::vector<bool> upwind_node(num_nodes);
   for (unsigned int n = 0; n < num_nodes; ++n)
   {
     if (_proto_flux[ph][n] >= 0.0) // upstream node
     {
       upwind_node[n] = true;
       // The mobility at the upstream node
       mob = mobility(n, ph);
       if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
       {
         // The derivative of the mobility wrt the PorousFlow variable
         dmob = dmobility(n, ph, pvar);
  
         for (_j = 0; _j < _phi.size(); _j++)
           _jacobian[ph][n][_j] *= mob;
  
         if (_test.size() == _phi.size())
           /* mobility at node=n depends only on the variables at node=n, by construction.  For
            * linear-lagrange variables, this means that Jacobian entries involving the derivative
            * of mobility will only be nonzero for derivatives wrt variables at node=n.  Hence the
            * [n][n] in the line below.  However, for other variable types (eg constant monomials)
            * I cannot tell what variable number contributes to the derivative.  However, in all
            * cases I can possibly imagine, the derivative is zero anyway, since in the full
            * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
            */
           _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
  
         for (_j = 0; _j < _phi.size(); _j++)
           dtotal_mass_out[_j] += _jacobian[ph][n][_j];
       }
       _proto_flux[ph][n] *= mob;
       total_mass_out += _proto_flux[ph][n];
     }
     else
     {
       upwind_node[n] = false;
       total_in -= _proto_flux[ph][n]; 
       if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
         for (_j = 0; _j < _phi.size(); _j++)
           dtotal_in[_j] -= _jacobian[ph][n][_j];
     }
   }
  
   // Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes,
   // weighted by their proto_flux values
   for (unsigned int n = 0; n < num_nodes; ++n)
   {
     if (!upwind_node[n]) // downstream node
     {
       if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
         for (_j = 0; _j < _phi.size(); _j++)
         {
           _jacobian[ph][n][_j] *= total_mass_out / total_in;
           _jacobian[ph][n][_j] +=
               _proto_flux[ph][n] * (dtotal_mass_out[_j] / total_in -
                                     dtotal_in[_j] * total_mass_out / total_in / total_in);
         }
       _proto_flux[ph][n] *= total_mass_out / total_in;
     }
   }
 }

Referenced by computeResidualAndJacobian().

◆ 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_jac	whether to compute the residual or jacobian
ph	fluid phase number
pvar	differentiate wrt to this PorousFlow variable (when computing the jacobian)

Definition at line 457 of file PorousFlowDarcyBase.C.

 {
   // The number of nodes in the element
   const unsigned int num_nodes = _test.size();
  
   std::vector<Real> mob(num_nodes);
   unsigned num_zero = 0;
   unsigned zero_mobility_node = std::numeric_limits<unsigned>::max();
   Real harmonic_mob = 0;
   for (unsigned n = 0; n < num_nodes; ++n)
   {
     mob[n] = mobility(n, ph);
     if (mob[n] == 0.0)
     {
       zero_mobility_node = n;
       num_zero++;
     }
     else
       harmonic_mob += 1.0 / mob[n];
   }
   if (num_zero > 0)
     harmonic_mob = 0.0;
   else
     harmonic_mob = (1.0 * num_nodes) / harmonic_mob;
  
   // d(harmonic_mob)/d(PorousFlow variable at node n)
   std::vector<Real> dharmonic_mob(num_nodes, 0.0);
   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
   {
     const Real harm2 = std::pow(harmonic_mob, 2) / (1.0 * num_nodes);
     if (num_zero == 0)
       for (unsigned n = 0; n < num_nodes; ++n)
         dharmonic_mob[n] = dmobility(n, ph, pvar) * harm2 / std::pow(mob[n], 2);
     else if (num_zero == 1)
       dharmonic_mob[zero_mobility_node] =
           num_nodes * dmobility(zero_mobility_node, ph, pvar); // other derivs are zero
     // if num_zero > 1 then all dharmonic_mob = 0.0
   }
  
   if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     for (unsigned n = 0; n < num_nodes; ++n)
       for (_j = 0; _j < _phi.size(); _j++)
       {
         _jacobian[ph][n][_j] *= harmonic_mob;
         if (_test.size() == _phi.size())
           _jacobian[ph][n][_j] += dharmonic_mob[_j] * _proto_flux[ph][n];
       }
  
   if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
     for (unsigned n = 0; n < num_nodes; ++n)
       _proto_flux[ph][n] *= harmonic_mob;
 }

Referenced by computeResidualAndJacobian().

◆ 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

nodenum	The node-number to evaluate the mobility for
phase	the fluid phase number

Reimplemented in PorousFlowAdvectiveFlux, and PorousFlowHeatAdvection.

Definition at line 511 of file PorousFlowDarcyBase.C.

 {
   return _fluid_density_node[nodenum][phase] / _fluid_viscosity[nodenum][phase];
 }

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

◆ 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_jac	whether to compute the residual or jacobian
ph	fluid phase number
pvar	differentiate wrt to this PorousFlow variable (when computing the jacobian)

Definition at line 420 of file PorousFlowDarcyBase.C.

 {
   // The number of nodes in the element
   const unsigned int num_nodes = _test.size();
  
   Real mob;
   Real dmob;
  
   // Use the raw nodal mobility
   for (unsigned int n = 0; n < num_nodes; ++n)
   {
     // The mobility at the node
     mob = mobility(n, ph);
     if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
     {
       // The derivative of the mobility wrt the PorousFlow variable
       dmob = dmobility(n, ph, pvar);
  
       for (_j = 0; _j < _phi.size(); _j++)
         _jacobian[ph][n][_j] *= mob;
  
       if (_test.size() == _phi.size())
         /* mobility at node=n depends only on the variables at node=n, by construction.  For
          * linear-lagrange variables, this means that Jacobian entries involving the derivative
          * of mobility will only be nonzero for derivatives wrt variables at node=n.  Hence the
          * [n][n] in the line below.  However, for other variable types (eg constant monomials)
          * I cannot tell what variable number contributes to the derivative.  However, in all
          * cases I can possibly imagine, the derivative is zero anyway, since in the full
          * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
          */
         _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
     }
     _proto_flux[ph][n] *= mob;
   }
 }

Referenced by computeResidualAndJacobian().

◆ timestepSetup()

void PorousFlowDarcyBase::timestepSetup ( )

overrideprotectedvirtual

Definition at line 87 of file PorousFlowDarcyBase.C.

 {
   Kernel::timestepSetup();
   _num_upwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
   _num_downwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
 }

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 99 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 105 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 111 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 120 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 123 of file PorousFlowDarcyBase.h.

Referenced by darcyQpJacobian().

◆ _dictator

const PorousFlowDictator& PorousFlowDarcyBase::_dictator

protected

PorousFlowDictator UserObject.

Definition at line 126 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 93 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 90 of file PorousFlowDarcyBase.h.

Referenced by darcyQpJacobian().

◆ _fallback_scheme

enum PorousFlowDarcyBase::FallbackEnum PorousFlowDarcyBase::_fallback_scheme

protected

Referenced by computeResidualAndJacobian().

◆ _fluid_density_node

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

protected

Fluid density for each phase (at the node)

Definition at line 96 of file PorousFlowDarcyBase.h.

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

◆ _fluid_density_qp

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

protected

Fluid density for each phase (at the qp)

Definition at line 102 of file PorousFlowDarcyBase.h.

Referenced by darcyQp(), and darcyQpJacobian().

◆ _fluid_viscosity

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

protected

Viscosity of each component in each phase.

Definition at line 108 of file PorousFlowDarcyBase.h.

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

◆ _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 141 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 117 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 132 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 167 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 181 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 129 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 174 of file PorousFlowDarcyBase.h.

Referenced by computeResidualAndJacobian(), and timestepSetup().

◆ _perm_derivs

const bool PorousFlowDarcyBase::_perm_derivs

protected

Flag to check whether permeabiity derivatives are non-zero.

Definition at line 135 of file PorousFlowDarcyBase.h.

Referenced by darcyQpJacobian().

◆ _permeability

const MaterialProperty<RealTensorValue>& PorousFlowDarcyBase::_permeability

protected

Permeability of porous material.

Definition at line 87 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 114 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 160 of file PorousFlowDarcyBase.h.

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

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

porous_flow/include/kernels/PorousFlowDarcyBase.h
porous_flow/src/kernels/PorousFlowDarcyBase.C

Public Member Functions

Protected Types

Protected Member Functions

Protected Attributes

Detailed Description

Member Enumeration Documentation

◆ FallbackEnum

◆ JacRes

Constructor & Destructor Documentation

◆ PorousFlowDarcyBase()

Member Function Documentation

◆ computeJacobian()

◆ computeOffDiagJacobian()

◆ computeQpResidual()

◆ computeResidual()

◆ computeResidualAndJacobian()

◆ darcyQp()

◆ darcyQpJacobian()

◆ dmobility()

◆ fullyUpwind()

◆ harmonicMean()

◆ mobility()

◆ quickUpwind()

◆ timestepSetup()

Member Data Documentation

◆ _dfluid_density_node_dvar

◆ _dfluid_density_qp_dvar

◆ _dfluid_viscosity_dvar

◆ _dgrad_p_dgrad_var

◆ _dgrad_p_dvar

◆ _dictator

◆ _dpermeability_dgradvar

◆ _dpermeability_dvar

◆ _fallback_scheme

◆ _fluid_density_node

◆ _fluid_density_qp

◆ _fluid_viscosity

◆ _full_upwind_threshold

◆ _grad_p

◆ _gravity

◆ _jacobian

◆ _num_downwinds

◆ _num_phases

◆ _num_upwinds

◆ _perm_derivs

◆ _permeability

◆ _pp

◆ _proto_flux