17 #include "libmesh/quadrature.h" 30 "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
32 "The variable representing the porepressure");
35 "A RichardsRelPerm UserObject (eg RichardsRelPermPowerGas) that defines the " 36 "fluid relative permeability as a function of the saturation Variable.");
39 "Flux according to Darcy-Richards flow. The Variable of this Kernel must be the saturation");
46 _pp(coupledValue(
"porepressure_variable")),
47 _grad_pp(coupledGradient(
"porepressure_variable")),
48 _pp_nodal(coupledDofValues(
"porepressure_variable")),
49 _pp_var(coupled(
"porepressure_variable")),
51 _viscosity(getParam<
Real>(
"fluid_viscosity")),
75 for (
unsigned int nodenum = 0; nodenum <
_num_nodes; ++nodenum)
114 upwind(
false,
true, jvar);
187 bool reached_steady =
true;
188 for (
unsigned int nodenum = 0; nodenum <
_num_nodes; ++nodenum)
192 reached_steady =
false;
196 reached_steady =
false;
200 Real total_mass_out = 0;
205 std::vector<Real> dtotal_mass_out;
206 std::vector<Real> dtotal_in;
214 for (
unsigned int nodenum = 0; nodenum <
_num_nodes; ++nodenum)
216 if (
_local_re(nodenum) >= 0 || reached_steady)
246 for (
unsigned int nodenum = 0; nodenum <
_num_nodes; ++nodenum)
252 _local_ke(nodenum,
_j) *= total_mass_out / total_in;
254 _local_re(nodenum) * (dtotal_mass_out[
_j] / total_in -
255 dtotal_in[
_j] * total_mass_out / total_in / total_in);
257 _local_re(nodenum) *= total_mass_out / total_in;
266 for (
unsigned int i = 0; i <
_save_in.size(); i++)
277 for (
unsigned int i = 0; i < rows; i++)
virtual Real computeQpResidual() override
Note that this is not the complete residual for the quadpoint In computeResidual we sum over the quad...
void prepareNodalValues()
calculates the nodal values of mobility, and derivatives thereof
const RichardsRelPerm & _relperm
fluid relative permeability
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the...
virtual void computeOffDiagJacobian(unsigned int jvar) override
this simply calls upwind
static InputParameters validParams()
registerMooseObject("RichardsApp", Q2PSaturationFlux)
void accumulateTaggedLocalResidual()
virtual Real drelperm(Real seff) const =0
derivative of relative permeability wrt effective saturation This must be over-ridden in your derived...
std::vector< MooseVariableFEBase *> _save_in
virtual Real ddensity(Real p) const =0
derivative of fluid density wrt porepressure This must be over-ridden in derived classes to provide a...
const MooseArray< Real > & _JxW
unsigned int number() const
unsigned int _num_nodes
number of nodes in the element
Base class for Richards relative permeability classes that provide relative permeability as a functio...
std::vector< Real > _dmobility_dp
d(_mobility)/d(porepressure) These are used in the jacobian calculations
Real _viscosity
fluid viscosity
std::vector< Real > _dmobility_ds
d(_mobility)/d(saturation) These are used in the jacobian calculations
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 tr...
const VariablePhiGradient & _grad_phi
const VariableGradient & _grad_pp
grad(porepressure) at the quadpoints
static const std::string density
const VariableValue & _pp
porepressure at the quadpoints
unsigned int _pp_var
variable number of the porepressure variable
virtual void computeJacobian() override
this simply calls upwind
DenseMatrix< Number > _local_ke
This is a fully upwinded flux Kernel The Variable of this Kernel should be the saturation.
virtual Real density(Real p) const =0
fluid density as a function of porepressure This must be over-ridden in derived classes to provide an...
TensorValue< Real > RealTensorValue
const RichardsDensity & _density
fluid density
const MaterialProperty< RealVectorValue > & _gravity
gravity
const VariableTestValue & _test
Q2PSaturationFlux(const InputParameters ¶meters)
std::vector< MooseVariableFEBase *> _diag_save_in
const QBase *const & _qrule
virtual Real relperm(Real seff) const =0
relative permeability as a function of effective saturation This must be over-ridden in your derived ...
void accumulateTaggedLocalMatrix()
static InputParameters validParams()
const MooseArray< Real > & _coord
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const VariableTestGradient & _grad_test
const VariableValue & _pp_nodal
porepressure at the nodes
const DoFValue & dofValues() const override
virtual void computeResidual() override
This simply calls upwind.
DenseVector< Number > _local_re
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
Real computeQpJac(unsigned int dvar)
the derivative of the flux without the upstream mobility terms
void prepareVectorTag(Assembly &assembly, unsigned int ivar)
void prepareMatrixTag(Assembly &assembly, unsigned int ivar, unsigned int jvar)
const MaterialProperty< RealTensorValue > & _permeability
permeability
const VariablePhiValue & _phi