Applies a flux sink to a boundary The sink is a piecewise linear function of porepressure (the "variable") at the quad points. More...

#include <RichardsPiecewiseLinearSink.h>

Inheritance diagram for RichardsPiecewiseLinearSink:

[legend]

Public Member Functions
	RichardsPiecewiseLinearSink (const InputParameters &parameters)

Protected Member Functions
virtual void	computeResidual () override

virtual Real	computeQpResidual () override

virtual void	computeJacobian () override

virtual Real	computeQpJacobian () override

virtual void	computeJacobianBlock (MooseVariableFEBase &jvar) override

virtual Real	computeQpOffDiagJacobian (unsigned int jvar) override

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

Real	jac (unsigned int wrt_num)
	derivative of residual wrt the wrt_num Richards variable More...

Protected Attributes
bool	_use_mobility
	whether to multiply the sink flux by permeability*density/viscosity More...

bool	_use_relperm
	whether to multiply the sink flux by relative permeability More...

bool	_fully_upwind
	whether to use full upwinding More...

LinearInterpolation	_sink_func
	piecewise-linear function of porepressure (this defines the strength of the sink) More...

const Function &	_m_func
	sink flux gets multiplied by this function More...

const RichardsVarNames &	_richards_name_UO
	holds info about the names and values of richards variable in the simulation More...

unsigned int	_num_p
	number of richards variables More...

unsigned int	_pvar
	the moose internal variable number corresponding to the porepressure of this sink flux More...

const RichardsDensity *	_density_UO
	user object defining the density. Only used if _fully_upwind = true More...

const RichardsSeff *	_seff_UO
	user object defining the effective saturation. Only used if _fully_upwind = true More...

const RichardsRelPerm *	_relperm_UO
	user object defining the relative permeability. Only used if _fully_upwind = true More...

const PostprocessorValue &	_area_pp
	area postprocessor. if given then all bare_fluxes are divided by this quantity More...

unsigned int	_num_nodes
	number of nodes in this element. Only used if _fully_upwind = true More...

std::vector< Real >	_nodal_density
	nodal values of fluid density These are used if _fully_upwind = true More...

std::vector< std::vector< Real > >	_dnodal_density_dv
	d(_nodal_density)/d(variable_ph) (variable_ph is the variable for phase=ph) These are used in the jacobian calculations if _fully_upwind = true More...

std::vector< Real >	_nodal_relperm
	nodal values of relative permeability These are used if _fully_upwind = true More...

std::vector< std::vector< Real > >	_dnodal_relperm_dv
	d(_nodal_relperm)/d(variable_ph) (variable_ph is the variable for phase=ph) These are used in the jacobian calculations if _fully_upwind = true More...

const MaterialProperty< std::vector< Real > > &	_pp
	porepressure values (only the _pvar component is used) More...

const MaterialProperty< std::vector< std::vector< Real > > > &	_dpp_dv
	d(porepressure_i)/d(variable_j) More...

const MaterialProperty< std::vector< Real > > &	_viscosity
	viscosity (only the _pvar component is used) More...

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

const MaterialProperty< std::vector< std::vector< Real > > > &	_dseff_dv
	derivative of effective saturation wrt variables only _dseff_dv[_pvar][i] is used for i being all variables More...

const MaterialProperty< std::vector< Real > > &	_rel_perm
	relative permeability (only the _pvar component is used) More...

const MaterialProperty< std::vector< std::vector< Real > > > &	_drel_perm_dv
	d(relperm_i)/d(variable_j) More...

const MaterialProperty< std::vector< Real > > &	_density
	fluid density (only the _pvar component is used) More...

const MaterialProperty< std::vector< std::vector< Real > > > &	_ddensity_dv
	d(density_i)/d(variable_j) More...

std::vector< const VariableValue * >	_ps_at_nodes
	Holds the values of pressures at all the nodes of the element Only used if _fully_upwind = true Eg: _ps_at_nodes[_pvar] is a pointer to this variable's nodal porepressure values So: (*_ps_at_nodes[_pvar])[i] = _var.dofValues()[i] = porepressure of pressure-variable _pvar at node i. More...

Detailed Description

Applies a flux sink to a boundary The sink is a piecewise linear function of porepressure (the "variable") at the quad points.

This is specified by _sink_func. In addition, this sink can be multiplied by: (1) the relative permeability of the fluid at the quad point. (2) perm_nn*density/viscosity, where perm_nn is the permeability tensor projected in the normal direction. (3) a Function (which can be time-dependent, for instance) and divided by: (4) an area Postprocessor

Definition at line 39 of file RichardsPiecewiseLinearSink.h.

Constructor & Destructor Documentation

◆ RichardsPiecewiseLinearSink()

RichardsPiecewiseLinearSink::RichardsPiecewiseLinearSink ( const InputParameters & parameters )

Definition at line 78 of file RichardsPiecewiseLinearSink.C.

   : IntegratedBC(parameters),
     _use_mobility(getParam<bool>("use_mobility")),
     _use_relperm(getParam<bool>("use_relperm")),
     _fully_upwind(getParam<bool>("fully_upwind")),
  
     _sink_func(getParam<std::vector<Real>>("pressures"),
                getParam<std::vector<Real>>("bare_fluxes")),
  
     _m_func(getFunction("multiplying_fcn")),
  
     _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
     _num_p(_richards_name_UO.num_v()),
     _pvar(_richards_name_UO.richards_var_num(_var.number())),
  
     // in the following, getUserObjectByName returns a reference (an alias) to a RichardsBLAH user
     // object, and the & turns it into a pointer
     _density_UO(_fully_upwind ? &getUserObjectByName<RichardsDensity>(
                                     getParam<std::vector<UserObjectName>>("density_UO")[_pvar])
                               : NULL),
     _seff_UO(_fully_upwind ? &getUserObjectByName<RichardsSeff>(
                                  getParam<std::vector<UserObjectName>>("seff_UO")[_pvar])
                            : NULL),
     _relperm_UO(_fully_upwind ? &getUserObjectByName<RichardsRelPerm>(
                                     getParam<std::vector<UserObjectName>>("relperm_UO")[_pvar])
                               : NULL),
  
     _area_pp(getPostprocessorValue("area_pp")),
  
     _num_nodes(0),
     _nodal_density(0),
     _dnodal_density_dv(0),
     _nodal_relperm(0),
     _dnodal_relperm_dv(0),
  
     _pp(getMaterialProperty<std::vector<Real>>("porepressure")),
     _dpp_dv(getMaterialProperty<std::vector<std::vector<Real>>>("dporepressure_dv")),
  
     _viscosity(getMaterialProperty<std::vector<Real>>("viscosity")),
     _permeability(getMaterialProperty<RealTensorValue>("permeability")),
  
     _dseff_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ds_eff_dv")),
  
     _rel_perm(getMaterialProperty<std::vector<Real>>("rel_perm")),
     _drel_perm_dv(getMaterialProperty<std::vector<std::vector<Real>>>("drel_perm_dv")),
  
     _density(getMaterialProperty<std::vector<Real>>("density")),
     _ddensity_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ddensity_dv"))
 {
   _ps_at_nodes.resize(_num_p);
   for (unsigned int pnum = 0; pnum < _num_p; ++pnum)
     _ps_at_nodes[pnum] = _richards_name_UO.nodal_var(pnum);
 }

Member Function Documentation

◆ computeJacobian()

void RichardsPiecewiseLinearSink::computeJacobian ( )

overrideprotectedvirtual

Definition at line 230 of file RichardsPiecewiseLinearSink.C.

 {
   if (_fully_upwind)
     prepareNodalValues();
   IntegratedBC::computeJacobian();
 }

◆ computeJacobianBlock()

void RichardsPiecewiseLinearSink::computeJacobianBlock ( MooseVariableFEBase & jvar )

overrideprotectedvirtual

Definition at line 244 of file RichardsPiecewiseLinearSink.C.

 {
   if (_fully_upwind)
     prepareNodalValues();
   IntegratedBC::computeJacobianBlock(jvar);
 }

◆ computeQpJacobian()

Real RichardsPiecewiseLinearSink::computeQpJacobian ( )

overrideprotectedvirtual

Definition at line 238 of file RichardsPiecewiseLinearSink.C.

 {
   return jac(_pvar);
 }

◆ computeQpOffDiagJacobian()

Real RichardsPiecewiseLinearSink::computeQpOffDiagJacobian ( unsigned int jvar )

overrideprotectedvirtual

Definition at line 252 of file RichardsPiecewiseLinearSink.C.

 {
   if (_richards_name_UO.not_richards_var(jvar))
     return 0.0;
   unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
   return jac(dvar);
 }

◆ computeQpResidual()

Real RichardsPiecewiseLinearSink::computeQpResidual ( )

overrideprotectedvirtual

Definition at line 184 of file RichardsPiecewiseLinearSink.C.

 {
   Real flux = 0;
   Real k = 0;
  
   if (!_fully_upwind)
   {
     flux = _test[_i][_qp] * _sink_func.sample(_pp[_qp][_pvar]);
     if (_use_mobility)
     {
       k = (_permeability[_qp] * _normals[_qp]) * _normals[_qp];
       flux *= _density[_qp][_pvar] * k / _viscosity[_qp][_pvar];
     }
     if (_use_relperm)
       flux *= _rel_perm[_qp][_pvar];
   }
   else
   {
     flux = _test[_i][_qp] * _sink_func.sample((*_ps_at_nodes[_pvar])[_i]);
     if (_use_mobility)
     {
       k = (_permeability[0] * _normals[_qp]) * _normals[_qp]; // assume that _permeability is
                                                               // constant throughout element so
                                                               // doesn't need to be upwinded
       flux *= _nodal_density[_i] * k /
               _viscosity[0][_pvar]; // assume that viscosity is constant throughout element
     }
     if (_use_relperm)
       flux *= _nodal_relperm[_i];
   }
  
   flux *= _m_func.value(_t, _q_point[_qp]);
  
   if (_area_pp == 0.0)
   {
     if (flux != 0)
       mooseError("RichardsPiecewiseLinearSink: flux is nonzero, but area is zero!\n");
     // if flux == 0, then leave it as zero.
   }
   else
     flux /= _area_pp;
  
   return flux;
 }

◆ computeResidual()

void RichardsPiecewiseLinearSink::computeResidual ( )

overrideprotectedvirtual

Definition at line 176 of file RichardsPiecewiseLinearSink.C.

 {
   if (_fully_upwind)
     prepareNodalValues();
   IntegratedBC::computeResidual();
 }

◆ jac()

Real RichardsPiecewiseLinearSink::jac ( unsigned int wrt_num )

protected

derivative of residual wrt the wrt_num Richards variable

Definition at line 261 of file RichardsPiecewiseLinearSink.C.

 {
   Real flux = 0;
   Real deriv = 0;
   Real k = 0;
   Real mob = 0;
   Real mobp = 0;
   Real phi = 0;
  
   if (!_fully_upwind)
   {
     flux = _sink_func.sample(_pp[_qp][_pvar]);
     deriv = _sink_func.sampleDerivative(_pp[_qp][_pvar]) * _dpp_dv[_qp][_pvar][wrt_num];
     phi = _phi[_j][_qp];
     if (_use_mobility)
     {
       k = (_permeability[_qp] * _normals[_qp]) * _normals[_qp];
       mob = _density[_qp][_pvar] * k / _viscosity[_qp][_pvar];
       mobp = _ddensity_dv[_qp][_pvar][wrt_num] * k / _viscosity[_qp][_pvar];
       deriv = mob * deriv + mobp * flux;
       flux *= mob;
     }
     if (_use_relperm)
       deriv = _rel_perm[_qp][_pvar] * deriv + _drel_perm_dv[_qp][_pvar][wrt_num] * flux;
   }
   else
   {
     if (_i != _j)
       return 0.0; // residual at node _i only depends on variables at that node
     flux = _sink_func.sample((*_ps_at_nodes[_pvar])[_i]);
     deriv = (_pvar == wrt_num ? _sink_func.sampleDerivative((*_ps_at_nodes[_pvar])[_i])
                               : 0); // NOTE: i'm assuming that the variables are pressure variables
     phi = 1;
     if (_use_mobility)
     {
       k = (_permeability[0] * _normals[_qp]) * _normals[_qp];
       mob = _nodal_density[_i] * k / _viscosity[0][_pvar];
       mobp = _dnodal_density_dv[_i][wrt_num] * k / _viscosity[0][_pvar];
       deriv = mob * deriv + mobp * flux;
       flux *= mob;
     }
     if (_use_relperm)
       deriv = _nodal_relperm[_i] * deriv + _dnodal_relperm_dv[_i][wrt_num] * flux;
   }
  
   deriv *= _m_func.value(_t, _q_point[_qp]);
  
   if (_area_pp == 0.0)
   {
     if (deriv != 0)
       mooseError("RichardsPiecewiseLinearSink: deriv is nonzero, but area is zero!\n");
     // if deriv == 0, then leave it as zero.
   }
   else
     deriv /= _area_pp;
  
   return _test[_i][_qp] * deriv * phi;
 }

Referenced by computeQpJacobian(), and computeQpOffDiagJacobian().

◆ prepareNodalValues()

void RichardsPiecewiseLinearSink::prepareNodalValues ( )

protected

calculates the nodal values of pressure, mobility, and derivatives thereof

Definition at line 133 of file RichardsPiecewiseLinearSink.C.

 {
   // NOTE: i'm assuming that all the richards variables are pressure values
  
   _num_nodes = (*_ps_at_nodes[_pvar]).size();
  
   Real p;
   Real seff;
   std::vector<Real> dseff_dp;
   Real drelperm_ds;
  
   _nodal_density.resize(_num_nodes);
   _dnodal_density_dv.resize(_num_nodes);
   _nodal_relperm.resize(_num_nodes);
   _dnodal_relperm_dv.resize(_num_nodes);
   dseff_dp.resize(_num_p);
   for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
   {
     // retrieve and calculate basic things at the node
     p = (*_ps_at_nodes[_pvar])[nodenum]; // pressure of fluid _pvar at node nodenum
  
     _nodal_density[nodenum] = _density_UO->density(p); // density of fluid _pvar at node nodenum
     _dnodal_density_dv[nodenum].resize(_num_p);
     for (unsigned int ph = 0; ph < _num_p; ++ph)
       _dnodal_density_dv[nodenum][ph] = 0;
     _dnodal_density_dv[nodenum][_pvar] = _density_UO->ddensity(p); // d(density)/dP
  
     seff = _seff_UO->seff(_ps_at_nodes,
                           nodenum); // effective saturation of fluid _pvar at node nodenum
     _seff_UO->dseff(
         _ps_at_nodes, nodenum, dseff_dp); // d(seff)/d(P_ph), for ph = 0, ..., _num_p - 1
  
     _nodal_relperm[nodenum] =
         _relperm_UO->relperm(seff); // relative permeability of fluid _pvar at node nodenum
     drelperm_ds = _relperm_UO->drelperm(seff); // d(relperm)/dseff
  
     _dnodal_relperm_dv[nodenum].resize(_num_p);
     for (unsigned int ph = 0; ph < _num_p; ++ph)
       _dnodal_relperm_dv[nodenum][ph] = drelperm_ds * dseff_dp[ph];
   }
 }

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

Member Data Documentation

◆ _area_pp

const PostprocessorValue& RichardsPiecewiseLinearSink::_area_pp

protected

area postprocessor. if given then all bare_fluxes are divided by this quantity

Definition at line 92 of file RichardsPiecewiseLinearSink.h.

Referenced by computeQpResidual(), and jac().

◆ _ddensity_dv

const MaterialProperty<std::vector<std::vector<Real> > >& RichardsPiecewiseLinearSink::_ddensity_dv

protected

d(density_i)/d(variable_j)

Definition at line 149 of file RichardsPiecewiseLinearSink.h.

Referenced by jac().

◆ _density

const MaterialProperty<std::vector<Real> >& RichardsPiecewiseLinearSink::_density

protected

fluid density (only the _pvar component is used)

Definition at line 146 of file RichardsPiecewiseLinearSink.h.

Referenced by computeQpResidual(), and jac().

◆ _density_UO

const RichardsDensity* RichardsPiecewiseLinearSink::_density_UO

protected

user object defining the density. Only used if _fully_upwind = true

Definition at line 83 of file RichardsPiecewiseLinearSink.h.

Referenced by prepareNodalValues().

◆ _dnodal_density_dv

std::vector<std::vector<Real> > RichardsPiecewiseLinearSink::_dnodal_density_dv

protected

d(_nodal_density)/d(variable_ph) (variable_ph is the variable for phase=ph) These are used in the jacobian calculations if _fully_upwind = true

Definition at line 107 of file RichardsPiecewiseLinearSink.h.

Referenced by jac(), and prepareNodalValues().

◆ _dnodal_relperm_dv

std::vector<std::vector<Real> > RichardsPiecewiseLinearSink::_dnodal_relperm_dv

protected

d(_nodal_relperm)/d(variable_ph) (variable_ph is the variable for phase=ph) These are used in the jacobian calculations if _fully_upwind = true

Definition at line 119 of file RichardsPiecewiseLinearSink.h.

Referenced by jac(), and prepareNodalValues().

◆ _dpp_dv

const MaterialProperty<std::vector<std::vector<Real> > >& RichardsPiecewiseLinearSink::_dpp_dv

protected

d(porepressure_i)/d(variable_j)

Definition at line 125 of file RichardsPiecewiseLinearSink.h.

Referenced by jac().

◆ _drel_perm_dv

const MaterialProperty<std::vector<std::vector<Real> > >& RichardsPiecewiseLinearSink::_drel_perm_dv

protected

d(relperm_i)/d(variable_j)

Definition at line 143 of file RichardsPiecewiseLinearSink.h.

Referenced by jac().

◆ _dseff_dv

const MaterialProperty<std::vector<std::vector<Real> > >& RichardsPiecewiseLinearSink::_dseff_dv

protected

derivative of effective saturation wrt variables only _dseff_dv[_pvar][i] is used for i being all variables

Definition at line 137 of file RichardsPiecewiseLinearSink.h.

◆ _fully_upwind

bool RichardsPiecewiseLinearSink::_fully_upwind

protected

whether to use full upwinding

Definition at line 65 of file RichardsPiecewiseLinearSink.h.

Referenced by computeJacobian(), computeJacobianBlock(), computeQpResidual(), computeResidual(), and jac().

◆ _m_func

const Function& RichardsPiecewiseLinearSink::_m_func

protected

sink flux gets multiplied by this function

Definition at line 71 of file RichardsPiecewiseLinearSink.h.

Referenced by computeQpResidual(), and jac().

◆ _nodal_density

std::vector<Real> RichardsPiecewiseLinearSink::_nodal_density

protected

nodal values of fluid density These are used if _fully_upwind = true

Definition at line 101 of file RichardsPiecewiseLinearSink.h.

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

◆ _nodal_relperm

std::vector<Real> RichardsPiecewiseLinearSink::_nodal_relperm

protected

nodal values of relative permeability These are used if _fully_upwind = true

Definition at line 113 of file RichardsPiecewiseLinearSink.h.

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

◆ _num_nodes

unsigned int RichardsPiecewiseLinearSink::_num_nodes

protected

number of nodes in this element. Only used if _fully_upwind = true

Definition at line 95 of file RichardsPiecewiseLinearSink.h.

Referenced by prepareNodalValues().

◆ _num_p

unsigned int RichardsPiecewiseLinearSink::_num_p

protected

number of richards variables

Definition at line 77 of file RichardsPiecewiseLinearSink.h.

Referenced by prepareNodalValues(), and RichardsPiecewiseLinearSink().

◆ _permeability

const MaterialProperty<RealTensorValue>& RichardsPiecewiseLinearSink::_permeability

protected

permeability

Definition at line 131 of file RichardsPiecewiseLinearSink.h.

Referenced by computeQpResidual(), and jac().

◆ _pp

const MaterialProperty<std::vector<Real> >& RichardsPiecewiseLinearSink::_pp

protected

porepressure values (only the _pvar component is used)

Definition at line 122 of file RichardsPiecewiseLinearSink.h.

Referenced by computeQpResidual(), and jac().

◆ _ps_at_nodes

std::vector<const VariableValue *> RichardsPiecewiseLinearSink::_ps_at_nodes

protected

Holds the values of pressures at all the nodes of the element Only used if _fully_upwind = true Eg: _ps_at_nodes[_pvar] is a pointer to this variable's nodal porepressure values So: (*_ps_at_nodes[_pvar])[i] = _var.dofValues()[i] = porepressure of pressure-variable _pvar at node i.

Definition at line 159 of file RichardsPiecewiseLinearSink.h.

Referenced by computeQpResidual(), jac(), prepareNodalValues(), and RichardsPiecewiseLinearSink().