Q2PSaturationFlux Class Reference

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

#include <Q2PSaturationFlux.h>

Inheritance diagram for Q2PSaturationFlux:

Public Member Functions
	Q2PSaturationFlux (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 &	_pp
	porepressure at the quadpoints More...
const VariableGradient &	_grad_pp
	grad(porepressure) at the quadpoints More...
const VariableValue &	_pp_nodal
	porepressure at the nodes More...
unsigned int	_pp_var
	variable number of the porepressure 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 saturation.

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 47 of file Q2PSaturationFlux.h.

Constructor & Destructor Documentation

Q2PSaturationFlux()

Q2PSaturationFlux::Q2PSaturationFlux

(

const InputParameters &

parameters

)

Definition at line 44 of file Q2PSaturationFlux.C.

  : Kernel(parameters),
    _density(getUserObject<RichardsDensity>("fluid_density")),
    _pp(coupledValue("porepressure_variable")),
    _grad_pp(coupledGradient("porepressure_variable")),
    _pp_nodal(coupledDofValues("porepressure_variable")),
    _pp_var(coupled("porepressure_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 Q2PSaturationFlux::computeJacobian ( )

overrideprotectedvirtual

this simply calls upwind

Definition at line 107 of file Q2PSaturationFlux.C.

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

computeOffDiagJacobian()

void Q2PSaturationFlux::computeOffDiagJacobian ( MooseVariableFEBase &  jvar )

overrideprotectedvirtual

this simply calls upwind

Definition at line 113 of file Q2PSaturationFlux.C.

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

computeQpJac()

Real Q2PSaturationFlux::computeQpJac ( unsigned int  dvar )

protected

the derivative of the flux without the upstream mobility terms

Definition at line 119 of file Q2PSaturationFlux.C.

Referenced by upwind().

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

computeQpResidual()

Real Q2PSaturationFlux::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 92 of file Q2PSaturationFlux.C.

Referenced by upwind().

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

computeResidual()

void Q2PSaturationFlux::computeResidual ( )

overrideprotectedvirtual

This simply calls upwind.

Definition at line 101 of file Q2PSaturationFlux.C.

{
  upwind(true, false, 0);
}

prepareNodalValues()

void Q2PSaturationFlux::prepareNodalValues ( )

protected

calculates the nodal values of mobility, and derivatives thereof

Definition at line 63 of file Q2PSaturationFlux.C.

Referenced by upwind().

{
  _num_nodes = _pp_nodal.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(_pp_nodal[nodenum]);      // fluid density at the node
    ddensity_dp = _density.ddensity(_pp_nodal[nodenum]); // d(fluid density)/dP at the node
    relperm = _relperm.relperm(
        _var.dofValues()[nodenum]); // relative permeability of the fluid at node nodenum
    drelperm_ds = _relperm.drelperm(_var.dofValues()[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;
  }
}

upwind()

void Q2PSaturationFlux::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 132 of file Q2PSaturationFlux.C.

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

{
  if (compute_jac && !(jvar == _var.number() || jvar == _pp_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 S
          _local_ke(nodenum, nodenum) += _dmobility_ds[nodenum] * _local_re(nodenum);
        else
          // deriv wrt P
          _local_ke(nodenum, nodenum) += _dmobility_dp[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());
    }
  }
}

Member Data Documentation

_density

const RichardsDensity& Q2PSaturationFlux::_density

protected

fluid density

Definition at line 86 of file Q2PSaturationFlux.h.

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

_dmobility_dp

std::vector<Real> Q2PSaturationFlux::_dmobility_dp

protected

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

Definition at line 125 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

_dmobility_ds

std::vector<Real> Q2PSaturationFlux::_dmobility_ds

protected

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

Definition at line 131 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

_grad_pp

const VariableGradient& Q2PSaturationFlux::_grad_pp

protected

grad(porepressure) at the quadpoints

Definition at line 92 of file Q2PSaturationFlux.h.

Referenced by computeQpResidual().

_gravity

const MaterialProperty<RealVectorValue>& Q2PSaturationFlux::_gravity

protected

gravity

Definition at line 107 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

_mobility

std::vector<Real> Q2PSaturationFlux::_mobility

protected

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

Definition at line 119 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

_num_nodes

unsigned int Q2PSaturationFlux::_num_nodes

protected

number of nodes in the element

Definition at line 113 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

_permeability

const MaterialProperty<RealTensorValue>& Q2PSaturationFlux::_permeability

protected

permeability

Definition at line 110 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

_pp

const VariableValue& Q2PSaturationFlux::_pp

protected

porepressure at the quadpoints

Definition at line 89 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

_pp_nodal

const VariableValue& Q2PSaturationFlux::_pp_nodal

protected

porepressure at the nodes

Definition at line 95 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues().

_pp_var

unsigned int Q2PSaturationFlux::_pp_var

protected

variable number of the porepressure variable

Definition at line 98 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and upwind().

_relperm

const RichardsRelPerm& Q2PSaturationFlux::_relperm

protected

fluid relative permeability

Definition at line 101 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues().

_viscosity

Real Q2PSaturationFlux::_viscosity

protected

fluid viscosity

Definition at line 104 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues().

---

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

• richards/include/kernels/Q2PSaturationFlux.h
• richards/src/kernels/Q2PSaturationFlux.C