This is a fully upwinded flux Kernel The Variable of this Kernel should be the porepressure. More...

#include <Q2PPorepressureFlux.h>

Inheritance diagram for Q2PPorepressureFlux:

Public Member Functions
	Q2PPorepressureFlux (const InputParameters &parameters)

Protected Member Functions
virtual Real	computeQpResidual () override
	Note that this is not the complete residual for the quadpoint In computeResidual we sum over the quadpoints and then add the upwind mobility parts. More...

virtual void	computeResidual () override
	This simply calls upwind. More...

virtual void	computeOffDiagJacobian (MooseVariableFEBase &jvar) override
	this simply calls upwind More...

virtual void	computeJacobian () override
	this simply calls upwind More...

Real	computeQpJac (unsigned int dvar)
	the derivative of the flux without the upstream mobility terms More...

void	upwind (bool compute_res, bool compute_jac, unsigned int jvar)
	Do the upwinding for both the residual and jacobian I've put both calculations in the same code to try to reduce code duplication. More...

void	prepareNodalValues ()
	calculates the nodal values of mobility, and derivatives thereof More...

Protected Attributes
const RichardsDensity &	_density
	fluid density More...

const VariableValue &	_sat
	saturation at the nodes More...

unsigned int	_sat_var
	variable number of the saturation variable More...

const RichardsRelPerm &	_relperm
	fluid relative permeability More...

Real	_viscosity
	fluid viscosity More...

const MaterialProperty< RealVectorValue > &	_gravity
	gravity More...

const MaterialProperty< RealTensorValue > &	_permeability
	permeability More...

unsigned int	_num_nodes
	number of nodes in the element More...

std::vector< Real >	_mobility
	nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the residual More...

std::vector< Real >	_dmobility_dp
	d(_mobility)/d(porepressure) These are used in the jacobian calculations More...

std::vector< Real >	_dmobility_ds
	d(_mobility)/d(saturation) These are used in the jacobian calculations More...

Detailed Description

This is a fully upwinded flux Kernel The Variable of this Kernel should be the porepressure.

The residual for the kernel is the darcy flux. This 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, and pot = Porepressure - density*gravity.x

However, in fully-upwind, the first step is to take the mobility outside the integral. R_i = mobility*int{flux_no_mob} = mobility*F_i NOTE: R_i is exactly the mass flux flowing out of node i. Similarly, F_i is a measure of fluid flowing out of node i.

This leads to the definition of upwinding:

If F_i is positive then R_i = mobility_i * F_i That is, we use the upwind value of mobility.

For the F_i<0 nodes we construct their R_i using mass conservation

Definition at line 46 of file Q2PPorepressureFlux.h.

Constructor & Destructor Documentation

◆ Q2PPorepressureFlux()

Q2PPorepressureFlux::Q2PPorepressureFlux ( const InputParameters & parameters )

Definition at line 43 of file Q2PPorepressureFlux.C.

   : Kernel(parameters),
     _density(getUserObject<RichardsDensity>("fluid_density")),
     _sat(coupledDofValues("saturation_variable")),
     _sat_var(coupled("saturation_variable")),
     _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
     _viscosity(getParam<Real>("fluid_viscosity")),
     _gravity(getMaterialProperty<RealVectorValue>("gravity")),
     _permeability(getMaterialProperty<RealTensorValue>("permeability")),
     _num_nodes(0),
     _mobility(0),
     _dmobility_dp(0),
     _dmobility_ds(0)
 {
 }

Member Function Documentation

◆ computeJacobian()

void Q2PPorepressureFlux::computeJacobian ( )

overrideprotectedvirtual

this simply calls upwind

Definition at line 103 of file Q2PPorepressureFlux.C.

 {
   upwind(false, true, _var.number());
 }

◆ computeOffDiagJacobian()

void Q2PPorepressureFlux::computeOffDiagJacobian ( MooseVariableFEBase & jvar )

overrideprotectedvirtual

this simply calls upwind

Definition at line 109 of file Q2PPorepressureFlux.C.

 {
   upwind(false, true, jvar.number());
 }

◆ computeQpJac()

Real Q2PPorepressureFlux::computeQpJac ( unsigned int dvar )

protected

the derivative of the flux without the upstream mobility terms

Definition at line 115 of file Q2PPorepressureFlux.C.

 {
   // this is just the derivative of the flux WITHOUT the upstream mobility terms
   // Those terms get added in during computeJacobian()
   if (dvar == _var.number())
     return _grad_test[_i][_qp] *
            (_permeability[_qp] *
             (_grad_phi[_j][_qp] - _density.ddensity(_u[_qp]) * _gravity[_qp] * _phi[_j][_qp]));
   else
     return 0;
 }

Referenced by upwind().

◆ computeQpResidual()

Real Q2PPorepressureFlux::computeQpResidual ( )

overrideprotectedvirtual

Note that this is not the complete residual for the quadpoint In computeResidual we sum over the quadpoints and then add the upwind mobility parts.

Definition at line 88 of file Q2PPorepressureFlux.C.

 {
   // note this is not the complete residual:
   // the upwind mobility parts get added in computeResidual
   return _grad_test[_i][_qp] *
          (_permeability[_qp] * (_grad_u[_qp] - _density.density(_u[_qp]) * _gravity[_qp]));
 }

Referenced by upwind().

◆ computeResidual()

void Q2PPorepressureFlux::computeResidual ( )

overrideprotectedvirtual

This simply calls upwind.

Definition at line 97 of file Q2PPorepressureFlux.C.

 {
   upwind(true, false, 0);
 }

◆ prepareNodalValues()

void Q2PPorepressureFlux::prepareNodalValues ( )

protected

calculates the nodal values of mobility, and derivatives thereof

Definition at line 60 of file Q2PPorepressureFlux.C.

 {
   _num_nodes = _sat.size();
  
   Real density;
   Real ddensity_dp;
   Real relperm;
   Real drelperm_ds;
  
   _mobility.resize(_num_nodes);
   _dmobility_dp.resize(_num_nodes);
   _dmobility_ds.resize(_num_nodes);
  
   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
   {
     density = _density.density(_var.dofValues()[nodenum]);      // fluid density at the node
     ddensity_dp = _density.ddensity(_var.dofValues()[nodenum]); // d(fluid density)/dP at the node
     relperm = _relperm.relperm(_sat[nodenum]); // relative permeability of the fluid at node nodenum
     drelperm_ds = _relperm.drelperm(_sat[nodenum]); // d(relperm)/dsat
  
     // calculate the mobility and its derivatives wrt P and S
     _mobility[nodenum] = density * relperm / _viscosity;
     _dmobility_dp[nodenum] = ddensity_dp * relperm / _viscosity;
     _dmobility_ds[nodenum] = density * drelperm_ds / _viscosity;
   }
 }

Referenced by upwind().

◆ upwind()

void Q2PPorepressureFlux::upwind	(	bool	compute_res,
		bool	compute_jac,
		unsigned int	jvar
	)

protected

Do the upwinding for both the residual and jacobian I've put both calculations in the same code to try to reduce code duplication.

This is because when calculating the jacobian we need to calculate the residual to see which nodes are upwind and which are downwind

Definition at line 128 of file Q2PPorepressureFlux.C.

 {
   if (compute_jac && !(jvar == _var.number() || jvar == _sat_var))
     return;
  
   // calculate the mobility values and their derivatives
   prepareNodalValues();
  
   // compute the residual 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
   DenseVector<Number> & re = _assembly.residualBlock(_var.number());
   _local_re.resize(re.size());
   _local_re.zero();
  
   for (_i = 0; _i < _test.size(); _i++)
     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
       _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
  
   DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
   if (compute_jac)
   {
     _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 (_qp = 0; _qp < _qrule->n_points(); _qp++)
           _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJac(jvar);
   }
  
   // Now perform the upwinding by multiplying the residuals at the
   // upstream nodes by their mobilities
   //
   // The residual for the kernel is the darcy flux.
   // This 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 _local_re calculation above.
   //
   // NOTE: Physically _local_re(_i) is a measure of fluid flowing out of node i
   // If we had left in mobility, it would be exactly the mass flux flowing out of node i.
   //
   // This leads to the definition of upwinding:
   // ***
   // If _local_re(i) 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 mass
   // flowing out of node i must be the sum of the masses flowing
   // into the other nodes.
  
   // FIRST:
   // this is a dirty way of getting around precision loss problems
   // and problems at steadystate where upwinding oscillates from
   // node-to-node causing nonconvergence.
   // I'm not sure if i actually need to do this in moose.  Certainly
   // in cosflow it is necessary.
   // I will code a better algorithm if necessary
   bool reached_steady = true;
   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
   {
     if (_local_re(nodenum) >= 1E-20)
     {
       reached_steady = false;
       break;
     }
   }
   reached_steady = false;
  
   // DEFINE VARIABLES USED TO ENSURE MASS CONSERVATION
   // total mass out - used for mass conservation
   Real total_mass_out = 0;
   // total flux in
   Real total_in = 0;
  
   // the following holds derivatives of these
   std::vector<Real> dtotal_mass_out;
   std::vector<Real> dtotal_in;
   if (compute_jac)
   {
     dtotal_mass_out.assign(_num_nodes, 0);
     dtotal_in.assign(_num_nodes, 0);
   }
  
   // PERFORM THE UPWINDING!
   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
   {
     if (_local_re(nodenum) >= 0 || reached_steady) // upstream node
     {
       if (compute_jac)
       {
         for (_j = 0; _j < _phi.size(); _j++)
           _local_ke(nodenum, _j) *= _mobility[nodenum];
         if (jvar == _var.number())
           // deriv wrt P
           _local_ke(nodenum, nodenum) += _dmobility_dp[nodenum] * _local_re(nodenum);
         else
           // deriv wrt S
           _local_ke(nodenum, nodenum) += _dmobility_ds[nodenum] * _local_re(nodenum);
         for (_j = 0; _j < _phi.size(); _j++)
           dtotal_mass_out[_j] += _local_ke(nodenum, _j);
       }
       _local_re(nodenum) *= _mobility[nodenum];
       total_mass_out += _local_re(nodenum);
     }
     else
     {
       total_in -= _local_re(nodenum); // note the -= means the result is positive
       if (compute_jac)
         for (_j = 0; _j < _phi.size(); _j++)
           dtotal_in[_j] -= _local_ke(nodenum, _j);
     }
   }
  
   // CONSERVE MASS
   // proportion the total_mass_out mass to the inflow nodes, weighting by their _local_re values
   if (!reached_steady)
     for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
       if (_local_re(nodenum) < 0)
       {
         if (compute_jac)
           for (_j = 0; _j < _phi.size(); _j++)
           {
             _local_ke(nodenum, _j) *= total_mass_out / total_in;
             _local_ke(nodenum, _j) +=
                 _local_re(nodenum) * (dtotal_mass_out[_j] / total_in -
                                       dtotal_in[_j] * total_mass_out / total_in / total_in);
           }
         _local_re(nodenum) *= total_mass_out / total_in;
       }
  
   // ADD RESULTS TO RESIDUAL OR JACOBIAN
   if (compute_res)
   {
     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 (compute_jac)
   {
     ke += _local_ke;
  
     if (_has_diag_save_in && jvar == _var.number())
     {
       const 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().

Member Data Documentation

◆ _density

const RichardsDensity& Q2PPorepressureFlux::_density

protected

fluid density

Definition at line 85 of file Q2PPorepressureFlux.h.

Referenced by computeQpJac(), computeQpResidual(), and prepareNodalValues().

◆ _dmobility_dp

std::vector<Real> Q2PPorepressureFlux::_dmobility_dp

protected

d(_mobility)/d(porepressure) These are used in the jacobian calculations

Definition at line 118 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _dmobility_ds

std::vector<Real> Q2PPorepressureFlux::_dmobility_ds

protected

d(_mobility)/d(saturation) These are used in the jacobian calculations

Definition at line 124 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _gravity

const MaterialProperty<RealVectorValue>& Q2PPorepressureFlux::_gravity

protected

gravity

Definition at line 100 of file Q2PPorepressureFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

◆ _mobility

std::vector<Real> Q2PPorepressureFlux::_mobility

protected

nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the residual

Definition at line 112 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _num_nodes

unsigned int Q2PPorepressureFlux::_num_nodes

protected

number of nodes in the element

Definition at line 106 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _permeability

const MaterialProperty<RealTensorValue>& Q2PPorepressureFlux::_permeability

protected

permeability

Definition at line 103 of file Q2PPorepressureFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

◆ _relperm

const RichardsRelPerm& Q2PPorepressureFlux::_relperm

protected

fluid relative permeability

Definition at line 94 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues().

◆ _sat

const VariableValue& Q2PPorepressureFlux::_sat

protected

saturation at the nodes

Definition at line 88 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues().

◆ _sat_var

unsigned int Q2PPorepressureFlux::_sat_var

protected

variable number of the saturation variable

Definition at line 91 of file Q2PPorepressureFlux.h.

Referenced by upwind().

◆ _viscosity

Real Q2PPorepressureFlux::_viscosity

protected

fluid viscosity

Definition at line 97 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues().

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

richards/include/kernels/Q2PPorepressureFlux.h
richards/src/kernels/Q2PPorepressureFlux.C

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ Q2PPorepressureFlux()

Member Function Documentation

◆ computeJacobian()

◆ computeOffDiagJacobian()

◆ computeQpJac()

◆ computeQpResidual()

◆ computeResidual()

◆ prepareNodalValues()

◆ upwind()

Member Data Documentation

◆ _density

◆ _dmobility_dp

◆ _dmobility_ds

◆ _gravity

◆ _mobility

◆ _num_nodes

◆ _permeability

◆ _relperm

◆ _sat

◆ _sat_var

◆ _viscosity