16 #include "libmesh/quadrature.h" 26 "mat_porosity",
"The porosity of the material. Should be between 0 and 1. Eg, 0.1");
29 "An auxillary variable describing porosity changes. " 30 "Porosity = mat_porosity + por_change. If this is not " 31 "provided, zero is used.");
34 "A list of auxillary variable describing permeability " 35 "changes. There must be 9 of these, corresponding to the " 36 "xx, xy, xz, yx, yy, yz, zx, zy, zz components respectively. " 37 " Permeability = mat_permeability*10^(perm_change).");
39 "richardsVarNames_UO",
"The UserObject that holds the list of Richards variable names.");
41 "relperm_UO",
"List of names of user objects that define relative permeability");
44 "List of name of user objects that define effective saturation as a function of " 48 "List of names of user objects that define saturation as a function of effective saturation");
50 "density_UO",
"List of names of user objects that define the fluid density");
52 "SUPG_UO",
"List of names of user objects that define the SUPG");
54 "viscosity",
"List of viscosity of fluids (Pa.s). Typical value for water is=1E-3");
57 "Gravitational acceleration (m/s^2) as a vector pointing downwards. Eg (0,0,-10)");
62 "If you are using second-order Lagrange shape functions you need to set this to false.");
69 _material_por(getParam<
Real>(
"mat_porosity")),
70 _por_change(coupledValue(
"por_change")),
71 _por_change_old(coupledValueOld(
"por_change")),
77 _porosity_old(declareProperty<
Real>(
"porosity_old")),
78 _porosity(declareProperty<
Real>(
"porosity")),
83 _num_p(_richards_name_UO.num_v()),
85 _perm_change(isCoupled(
"perm_change")
86 ? coupledValues(
"perm_change")
89 _trace_perm(_material_perm.tr()),
91 _material_viscosity(getParam<
std::vector<
Real>>(
"viscosity")),
93 _pp_old(declareProperty<
std::vector<
Real>>(
"porepressure_old")),
94 _pp(declareProperty<
std::vector<
Real>>(
"porepressure")),
95 _dpp_dv(declareProperty<
std::vector<
std::vector<
Real>>>(
"dporepressure_dv")),
96 _d2pp_dv(declareProperty<
std::vector<
std::vector<
std::vector<
Real>>>>(
"d2porepressure_dvdv")),
98 _viscosity(declareProperty<
std::vector<
Real>>(
"viscosity")),
100 _density_old(declareProperty<
std::vector<
Real>>(
"density_old")),
101 _density(declareProperty<
std::vector<
Real>>(
"density")),
102 _ddensity_dv(declareProperty<
std::vector<
std::vector<
Real>>>(
"ddensity_dv")),
104 _seff_old(declareProperty<
std::vector<
Real>>(
"s_eff_old")),
105 _seff(declareProperty<
std::vector<
Real>>(
"s_eff")),
106 _dseff_dv(declareProperty<
std::vector<
std::vector<
Real>>>(
"ds_eff_dv")),
107 _d2seff_dv(declareProperty<
std::vector<
std::vector<
std::vector<
Real>>>>(
"d2s_eff_dvdv")),
109 _sat_old(declareProperty<
std::vector<
Real>>(
"sat_old")),
110 _sat(declareProperty<
std::vector<
Real>>(
"sat")),
111 _dsat_dv(declareProperty<
std::vector<
std::vector<
Real>>>(
"dsat_dv")),
113 _rel_perm(declareProperty<
std::vector<
Real>>(
"rel_perm")),
114 _drel_perm_dv(declareProperty<
std::vector<
std::vector<
Real>>>(
"drel_perm_dv")),
116 _mass_old(declareProperty<
std::vector<
Real>>(
"mass_old")),
117 _mass(declareProperty<
std::vector<
Real>>(
"mass")),
118 _dmass(declareProperty<
std::vector<
std::vector<
Real>>>(
"dmass")),
122 _dflux_no_mob_dgradv(
136 _dtauvel_SUPG_dgradp(
145 const std::vector<MooseVariableFEBase *> & coupled_vars =
147 for (
unsigned int i = 0; i < coupled_vars.size(); i++)
151 if (_material_por <= 0 || _material_por >= 1)
156 " components of perm_change must be given to a RichardsMaterial. You supplied ",
161 getParam<std::vector<UserObjectName>>(
"relperm_UO").size() ==
_num_p &&
162 getParam<std::vector<UserObjectName>>(
"seff_UO").size() ==
_num_p &&
163 getParam<std::vector<UserObjectName>>(
"sat_UO").size() ==
_num_p &&
164 getParam<std::vector<UserObjectName>>(
"density_UO").size() ==
_num_p &&
165 getParam<std::vector<UserObjectName>>(
"SUPG_UO").size() ==
_num_p))
168 " Richards fluid variables, so you need to specify this " 169 "number of viscosities, relperm_UO, seff_UO, sat_UO, " 170 "density_UO, SUPG_UO");
183 for (
unsigned int i = 0; i <
_num_p; ++i)
194 getParam<std::vector<UserObjectName>>(
"relperm_UO")[i]);
196 &getUserObjectByName<RichardsSeff>(getParam<std::vector<UserObjectName>>(
"seff_UO")[i]);
198 &getUserObjectByName<RichardsSat>(getParam<std::vector<UserObjectName>>(
"sat_UO")[i]);
200 getParam<std::vector<UserObjectName>>(
"density_UO")[i]);
202 &getUserObjectByName<RichardsSUPG>(getParam<std::vector<UserObjectName>>(
"SUPG_UO")[i]);
213 for (
unsigned int i = 0; i <
_num_p; ++i)
221 for (
unsigned int qp = 0; qp <
_qrule->n_points(); qp++)
235 for (
unsigned int i = 0; i <
_num_p; ++i)
244 for (
unsigned int j = 0;
j <
_num_p; ++
j)
254 for (
unsigned int j = 0;
j <
_num_p; ++
j)
270 for (
unsigned int i = 0; i <
_num_p; ++i)
277 for (
unsigned int i = 0; i <
_num_p; ++i)
282 for (
unsigned int j = 0;
j <
_num_p; ++
j)
290 for (
unsigned int i = 0; i <
_num_p; ++i)
295 for (
unsigned int j = 0;
j <
_num_p; ++
j)
302 for (
unsigned int i = 0; i <
_num_p; ++i)
306 for (
unsigned int j = 0;
j <
_num_p; ++
j)
314 for (
unsigned int i = 0; i <
_num_p; ++i)
319 for (
unsigned int j = 0;
j <
_num_p; ++
j)
328 for (
unsigned int i = 0; i <
_num_p; ++i)
333 for (
unsigned int j = 0;
j <
_num_p; ++
j)
337 for (
unsigned int j = 0;
j <
_num_p; ++
j)
345 for (
unsigned int i = 0; i <
_num_p; ++i)
350 for (
unsigned int j = 0;
j <
_num_p; ++
j)
360 for (
unsigned int j = 0;
j <
_num_p; ++
j)
373 for (
unsigned int i = 0; i <
_num_p; ++i)
378 for (
unsigned int j = 0;
j <
_num_p; ++
j)
392 for (
unsigned int i = 0; i <
_num_p; ++i)
401 for (
unsigned int j = 0;
j <
_num_p; ++
j)
404 for (
unsigned int k = 0;
k <
_num_p; ++
k)
413 for (
unsigned int j = 0;
j <
_num_p; ++
j)
416 for (
unsigned int k = 0;
k <
_num_p; ++
k)
422 for (
unsigned int j = 0;
j <
_num_p; ++
j)
424 for (
unsigned int k = 0;
k <
_num_p; ++
k)
446 for (
unsigned int j = 0;
j <
_num_p; ++
j)
447 for (
unsigned int k = 0;
k <
_num_p; ++
k)
450 for (
unsigned int j = 0;
j <
_num_p; ++
j)
452 for (
unsigned int k = 0;
k <
_num_p; ++
k)
469 for (
unsigned int i = 0; i <
_num_p; ++i)
487 const std::vector<Real> & dxidx(fe->get_dxidx());
488 const std::vector<Real> & dxidy(fe->get_dxidy());
489 const std::vector<Real> & dxidz(fe->get_dxidz());
490 const std::vector<Real> & detadx(fe->get_detadx());
491 const std::vector<Real> & detady(fe->get_detady());
492 const std::vector<Real> & detadz(fe->get_detadz());
493 const std::vector<Real> & dzetadx(fe->get_dzetadx());
494 const std::vector<Real> & dzetady(fe->get_dzetady());
495 const std::vector<Real> & dzetadz(fe->get_dzetadz());
497 for (
unsigned int qp = 0; qp <
_qrule->n_points(); qp++)
501 mooseAssert(qp < dxidx.size(),
"Insufficient data in dxidx array!");
502 mooseAssert(qp < dxidy.size(),
"Insufficient data in dxidy array!");
503 mooseAssert(qp < dxidz.size(),
"Insufficient data in dxidz array!");
506 mooseAssert(qp < detadx.size(),
"Insufficient data in detadx array!");
507 mooseAssert(qp < detady.size(),
"Insufficient data in detady array!");
508 mooseAssert(qp < detadz.size(),
"Insufficient data in detadz array!");
512 mooseAssert(qp < dzetadx.size(),
"Insufficient data in dzetadx array!");
513 mooseAssert(qp < dzetady.size(),
"Insufficient data in dzetady array!");
514 mooseAssert(qp < dzetadz.size(),
"Insufficient data in dzetadz array!");
522 eta_prime(0) = detadx[qp];
523 eta_prime(1) = detady[qp];
527 eta_prime(2) = detadz[qp];
528 zeta_prime(0) = dzetadx[qp];
529 zeta_prime(1) = dzetady[qp];
530 zeta_prime(2) = dzetadz[qp];
534 for (
unsigned int i = 0; i <
_num_p; ++i)
546 (*
_material_SUPG_UO[i]).dbb2_dgradp(vel, dvel_dgradp, xi_prime, eta_prime, zeta_prime);
558 dtauvel_dgradp(
j,
k) +=
559 dtau_dgradp(
j) * vel(
k);
560 for (
unsigned int j = 0;
j <
_num_p; ++
j)
564 for (
unsigned int j = 0;
j <
_num_p; ++
j)
578 for (
unsigned int qp = 0; qp <
_qrule->n_points(); qp++)
592 for (
unsigned int qp = 0; qp <
_qrule->n_points(); qp++)
597 for (
unsigned int qp = 0; qp <
_qrule->n_points(); qp++)
601 for (
unsigned int qp = 0; qp <
_qrule->n_points(); qp++)
605 bool trivial_supg =
true;
606 for (
unsigned int i = 0; i <
_num_p; ++i)
MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dflux_no_mob_dgradv
d(_flux_no_mob_i)/d(grad(variable_j))
MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _d2pp_dv
d^2(porepressure_i)/d(variable_j)/d(variable_k)
MaterialProperty< std::vector< Real > > & _pp
porepressure(s)
virtual bool isCoupled(const std::string &var_name, unsigned int i=0) const
const QBase *const & _qrule
MaterialProperty< RealVectorValue > & _gravity
MaterialProperty< std::vector< std::vector< std::vector< RealTensorValue > > > > & _d2flux_dvdgradv
d^2(Richards flux_i)/d(variable_j)/d(grad(variable_k)), here flux_i is the i_th flux, which is itself a RealVectorValue. We should have _d2flux_dvdgradv[i][j][k] = _d2flux_dgradvdv[i][k][j], but i think it is more clear having both, and hopefully not a blowout on memory/CPU.
std::vector< std::vector< std::vector< Real > > > _d2rel_perm_dv
d^2(relperm_i)/dP_j/dP_k - used in various derivative calculations
RichardsMaterial(const InputParameters ¶meters)
MaterialProperty< std::vector< std::vector< std::vector< RealTensorValue > > > > & _d2flux_dgradvdv
d^2(Richards flux_i)/d(grad(variable_j))/d(variable_k), here flux_i is the i_th flux, which is itself a RealVectorValue
std::vector< const RichardsRelPerm * > _material_relperm_UO
const VariableValue * richards_vals_old(unsigned int richards_var_num) const
a vector of pointers to old VariableValues
MaterialProperty< std::vector< Real > > & _pp_old
old values of porepressure(s)
const RichardsVarNames & _richards_name_UO
The variable names userobject for the Richards variables.
MaterialProperty< std::vector< RealVectorValue > > & _flux
fluid flux (a vector of fluxes for multicomponent)
std::vector< Real > _material_viscosity
const std::vector< const VariableValue * > _perm_change
static const std::string density
static constexpr std::size_t dim
This holds maps between pressure_var or pressure_var, sat_var used in RichardsMaterial and kernels...
MaterialProperty< std::vector< Real > > & _sat
saturation (vector of saturations in case of multiphase)
std::vector< const RichardsSUPG * > _material_SUPG_UO
MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
d(relperm_i)/d(variable_j)
virtual void resize(const std::size_t size) override final
Real _material_por
porosity as entered by the user
const VariableValue & _por_change
porosity changes. if not entered they default to zero
void zero2ndDerivedQuantities(unsigned int qp)
Zeroes 2nd derivatives of the flux.
TensorValue< Real > RealTensorValue
MaterialProperty< std::vector< Real > > & _seff_old
old effective saturation
registerMooseObject("RichardsApp", RichardsMaterial)
MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
d(porepressure_i)/d(variable_j)
void computePandSeff()
computes the quadpoint values of the porepressure(s) and effective saturation(s), and their derivativ...
MaterialProperty< std::vector< std::vector< std::vector< RealVectorValue > > > > & _d2flux_dvdv
d^2(Richards flux_i)/d(variable_j)/d(variable_k), here flux_i is the i_th flux, which is itself a Rea...
MaterialProperty< RealTensorValue > & _permeability
MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dflux_no_mob_dv
d(_flux_no_mob_i)/d(variable_j)
std::string var_types() const
return the _var_types string
static InputParameters validParams()
const VariableGradient * grad_var(unsigned int richards_var_num) const
a vector of pointers to grad(Variable)
virtual unsigned int dimension() const
MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dtauvel_SUPG_dgradp
std::vector< const VariableGradient * > _grad_p
MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dflux_dv
d(Richards flux_i)/d(variable_j), here flux_i is the i_th flux, which is itself a RealVectorValue ...
std::vector< const VariableValue * > _pressure_vals
MaterialProperty< std::vector< Real > > & _rel_perm
relative permeability (vector of relative permeabilities in case of multiphase)
MaterialProperty< std::vector< Real > > & _mass
fluid mass (a vector of masses for multicomponent)
RealVectorValue _material_gravity
gravity as entered by user
MaterialProperty< std::vector< RealVectorValue > > & _flux_no_mob
permeability*(grad(P) - density*gravity) (a vector of these for multicomponent)
std::vector< const RichardsDensity * > _material_density_UO
void computeDerivedQuantities(unsigned int qp)
Computes the "derived" quantities — those that depend on porepressure(s) and effective saturation(s)...
MaterialProperty< std::vector< Real > > & _density
fluid density (or densities for multiphase problems)
std::vector< const VariableValue * > _pressure_old_vals
MaterialProperty< std::vector< Real > > & _sat_old
old saturation
const std::vector< MooseVariableFieldBase *> & getCoupledMooseVars() const
OutputTools< Real >::VariableValue VariableValue
std::vector< const RichardsSeff * > _material_seff_UO
const VariableValue * richards_vals(unsigned int richards_var_num) const
a vector of pointers to VariableValues
unsigned int coupledComponents(const std::string &var_name) const
MaterialProperty< std::vector< Real > > & _mass_old
old value of fluid mass (a vector of masses for multicomponent)
void addMooseVariableDependency(MooseVariableFieldBase *var)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MaterialProperty< Real > & _porosity
void computeSUPG()
Computes the tauvel_SUPG and its derivatives.
MaterialProperty< std::vector< RealVectorValue > > & _tauvel_SUPG
MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _d2seff_dv
d^2(Seff_i)/d(variable_j)/d(variable_k)
void compute2ndDerivedQuantities(unsigned int qp)
Computes 2nd derivatives of the flux.
std::vector< const RichardsSat * > _material_sat_UO
void zeroSUPG(unsigned int qp)
Assigns and zeroes the MaterialProperties associated with SUPG.
MaterialProperty< std::vector< Real > > & _density_old
old fluid density (or densities for multiphase problems)
IntRange< T > make_range(T beg, T end)
MaterialProperty< Real > & _porosity_old
material properties
const VariableValue & _por_change_old
void mooseError(Args &&... args) const
MaterialProperty< std::vector< Real > > & _seff
effective saturation (vector of effective saturations in case of multiphase)
RealTensorValue _material_perm
permeability as entered by the user
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dflux_dgradv
d(Richards flux_i)/d(grad(variable_j)), here flux_i is the i_th flux, which is itself a RealVectorVal...
Real _trace_perm
trace of permeability tensor
const FEBase *const & getFE(FEType type, unsigned int dim) const
MaterialProperty< std::vector< std::vector< Real > > > & _ddensity_dv
d(density_i)/d(variable_j)
MaterialProperty< std::vector< std::vector< Real > > > & _dsat_dv
d(saturation_i)/d(variable_j)
MaterialProperty< std::vector< std::vector< Real > > > & _dseff_dv
d(Seff_i)/d(variable_j)
MooseUnits pow(const MooseUnits &, int)
static const std::string k
MaterialProperty< std::vector< std::vector< Real > > > & _dmass
d(fluid mass_i)/dP_j (a vector of masses for multicomponent)
static InputParameters validParams()
MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dtauvel_SUPG_dp
std::vector< std::vector< std::vector< Real > > > _d2density
d^2(density)/dp_j/dP_k - used in various derivative calculations
const Elem *const & _current_elem
MaterialProperty< std::vector< Real > > & _viscosity
fluid viscosity (or viscosities in the multiphase case)
virtual void computeProperties()