www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
RichardsMaterial Class Reference

#include <RichardsMaterial.h>

Inheritance diagram for RichardsMaterial:
[legend]

Public Member Functions

 RichardsMaterial (const InputParameters &parameters)
 

Protected Member Functions

virtual void computeProperties ()
 
void computePandSeff ()
 computes the quadpoint values of the porepressure(s) and effective saturation(s), and their derivatives wrt the variables in the system. More...
 
void computeDerivedQuantities (unsigned int qp)
 Computes the "derived" quantities — those that depend on porepressure(s) and effective saturation(s) — such as density, relative permeability, mass, flux, etc. More...
 
void compute2ndDerivedQuantities (unsigned int qp)
 Computes 2nd derivatives of the flux. More...
 
void zeroSUPG (unsigned int qp)
 Assigns and zeroes the MaterialProperties associated with SUPG. More...
 
void computeSUPG ()
 Computes the tauvel_SUPG and its derivatives. More...
 

Protected Attributes

Real _material_por
 porosity as entered by the user More...
 
const VariableValue & _por_change
 porosity changes. if not entered they default to zero More...
 
const VariableValue & _por_change_old
 
RealTensorValue _material_perm
 permeability as entered by the user More...
 
RealVectorValue _material_gravity
 gravity as entered by user More...
 
MaterialProperty< Real > & _porosity_old
 material properties More...
 
MaterialProperty< Real > & _porosity
 
MaterialProperty< RealTensorValue > & _permeability
 
MaterialProperty< RealVectorValue > & _gravity
 
const RichardsVarNames_richards_name_UO
 The variable names userobject for the Richards variables. More...
 
unsigned int _num_p
 
std::vector< const RichardsRelPerm * > _material_relperm_UO
 
std::vector< const RichardsSeff * > _material_seff_UO
 
std::vector< const RichardsSat * > _material_sat_UO
 
std::vector< const RichardsDensity * > _material_density_UO
 
std::vector< const RichardsSUPG * > _material_SUPG_UO
 
std::vector< const VariableValue * > _perm_change
 

Private Member Functions

void zero2ndDerivedQuantities (unsigned int qp)
 Zeroes 2nd derivatives of the flux. More...
 

Private Attributes

Real _trace_perm
 trace of permeability tensor More...
 
std::vector< Real > _material_viscosity
 
MaterialProperty< std::vector< Real > > & _pp_old
 old values of porepressure(s) More...
 
MaterialProperty< std::vector< Real > > & _pp
 porepressure(s) More...
 
MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
 d(porepressure_i)/d(variable_j) More...
 
MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _d2pp_dv
 d^2(porepressure_i)/d(variable_j)/d(variable_k) More...
 
MaterialProperty< std::vector< Real > > & _viscosity
 fluid viscosity (or viscosities in the multiphase case) More...
 
MaterialProperty< std::vector< Real > > & _density_old
 old fluid density (or densities for multiphase problems) More...
 
MaterialProperty< std::vector< Real > > & _density
 fluid density (or densities for multiphase problems) More...
 
MaterialProperty< std::vector< std::vector< Real > > > & _ddensity_dv
 d(density_i)/d(variable_j) More...
 
MaterialProperty< std::vector< Real > > & _seff_old
 old effective saturation More...
 
MaterialProperty< std::vector< Real > > & _seff
 effective saturation (vector of effective saturations in case of multiphase) More...
 
MaterialProperty< std::vector< std::vector< Real > > > & _dseff_dv
 d(Seff_i)/d(variable_j) More...
 
MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _d2seff_dv
 d^2(Seff_i)/d(variable_j)/d(variable_k) More...
 
MaterialProperty< std::vector< Real > > & _sat_old
 old saturation More...
 
MaterialProperty< std::vector< Real > > & _sat
 saturation (vector of saturations in case of multiphase) More...
 
MaterialProperty< std::vector< std::vector< Real > > > & _dsat_dv
 d(saturation_i)/d(variable_j) More...
 
MaterialProperty< std::vector< Real > > & _rel_perm
 relative permeability (vector of relative permeabilities in case of multiphase) More...
 
MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
 d(relperm_i)/d(variable_j) More...
 
MaterialProperty< std::vector< Real > > & _mass_old
 old value of fluid mass (a vector of masses for multicomponent) More...
 
MaterialProperty< std::vector< Real > > & _mass
 fluid mass (a vector of masses for multicomponent) More...
 
MaterialProperty< std::vector< std::vector< Real > > > & _dmass
 d(fluid mass_i)/dP_j (a vector of masses for multicomponent) More...
 
MaterialProperty< std::vector< RealVectorValue > > & _flux_no_mob
 permeability*(grad(P) - density*gravity) (a vector of these for multicomponent) More...
 
MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dflux_no_mob_dv
 d(_flux_no_mob_i)/d(variable_j) More...
 
MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dflux_no_mob_dgradv
 d(_flux_no_mob_i)/d(grad(variable_j)) More...
 
MaterialProperty< std::vector< RealVectorValue > > & _flux
 fluid flux (a vector of fluxes for multicomponent) More...
 
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 More...
 
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 RealVectorValue More...
 
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 RealVectorValue More...
 
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 More...
 
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. More...
 
MaterialProperty< std::vector< RealVectorValue > > & _tauvel_SUPG
 
MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dtauvel_SUPG_dgradp
 
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 More...
 
std::vector< std::vector< std::vector< Real > > > _d2rel_perm_dv
 d^2(relperm_i)/dP_j/dP_k - used in various derivative calculations More...
 
std::vector< const VariableValue * > _pressure_vals
 
std::vector< const VariableValue * > _pressure_old_vals
 
std::vector< const VariableGradient * > _grad_p
 

Detailed Description

Definition at line 28 of file RichardsMaterial.h.

Constructor & Destructor Documentation

◆ RichardsMaterial()

RichardsMaterial::RichardsMaterial ( const InputParameters &  parameters)

Definition at line 67 of file RichardsMaterial.C.

68  : Material(parameters),
69 
70  _material_por(getParam<Real>("mat_porosity")),
71  _por_change(coupledValue("por_change")),
72  _por_change_old(coupledValueOld("por_change")),
73 
74  _material_perm(getParam<RealTensorValue>("mat_permeability")),
75 
76  _material_gravity(getParam<RealVectorValue>("gravity")),
77 
78  _porosity_old(declareProperty<Real>("porosity_old")),
79  _porosity(declareProperty<Real>("porosity")),
80  _permeability(declareProperty<RealTensorValue>("permeability")),
81  _gravity(declareProperty<RealVectorValue>("gravity")),
82 
83  _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
85 
87 
88  _material_viscosity(getParam<std::vector<Real>>("viscosity")),
89 
90  _pp_old(declareProperty<std::vector<Real>>("porepressure_old")),
91  _pp(declareProperty<std::vector<Real>>("porepressure")),
92  _dpp_dv(declareProperty<std::vector<std::vector<Real>>>("dporepressure_dv")),
93  _d2pp_dv(declareProperty<std::vector<std::vector<std::vector<Real>>>>("d2porepressure_dvdv")),
94 
95  _viscosity(declareProperty<std::vector<Real>>("viscosity")),
96 
97  _density_old(declareProperty<std::vector<Real>>("density_old")),
98  _density(declareProperty<std::vector<Real>>("density")),
99  _ddensity_dv(declareProperty<std::vector<std::vector<Real>>>("ddensity_dv")),
100 
101  _seff_old(declareProperty<std::vector<Real>>("s_eff_old")),
102  _seff(declareProperty<std::vector<Real>>("s_eff")),
103  _dseff_dv(declareProperty<std::vector<std::vector<Real>>>("ds_eff_dv")),
104  _d2seff_dv(declareProperty<std::vector<std::vector<std::vector<Real>>>>("d2s_eff_dvdv")),
105 
106  _sat_old(declareProperty<std::vector<Real>>("sat_old")),
107  _sat(declareProperty<std::vector<Real>>("sat")),
108  _dsat_dv(declareProperty<std::vector<std::vector<Real>>>("dsat_dv")),
109 
110  _rel_perm(declareProperty<std::vector<Real>>("rel_perm")),
111  _drel_perm_dv(declareProperty<std::vector<std::vector<Real>>>("drel_perm_dv")),
112 
113  _mass_old(declareProperty<std::vector<Real>>("mass_old")),
114  _mass(declareProperty<std::vector<Real>>("mass")),
115  _dmass(declareProperty<std::vector<std::vector<Real>>>("dmass")),
116 
117  _flux_no_mob(declareProperty<std::vector<RealVectorValue>>("flux_no_mob")),
118  _dflux_no_mob_dv(declareProperty<std::vector<std::vector<RealVectorValue>>>("dflux_no_mob_dv")),
120  declareProperty<std::vector<std::vector<RealTensorValue>>>("dflux_no_mob_dgradv")),
121 
122  _flux(declareProperty<std::vector<RealVectorValue>>("flux")),
123  _dflux_dv(declareProperty<std::vector<std::vector<RealVectorValue>>>("dflux_dv")),
124  _dflux_dgradv(declareProperty<std::vector<std::vector<RealTensorValue>>>("dflux_dgradv")),
125  _d2flux_dvdv(
126  declareProperty<std::vector<std::vector<std::vector<RealVectorValue>>>>("d2flux_dvdv")),
128  declareProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>("d2flux_dgradvdv")),
130  declareProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>("d2flux_dvdgradv")),
131 
132  _tauvel_SUPG(declareProperty<std::vector<RealVectorValue>>("tauvel_SUPG")),
134  declareProperty<std::vector<std::vector<RealTensorValue>>>("dtauvel_SUPG_dgradv")),
135  _dtauvel_SUPG_dp(declareProperty<std::vector<std::vector<RealVectorValue>>>("dtauvel_SUPG_dv"))
136 
137 {
138 
139  // Need to add the variables that the user object is coupled to as dependencies so MOOSE will
140  // compute them
141  {
142  const std::vector<MooseVariableFEBase *> & coupled_vars =
143  _richards_name_UO.getCoupledMooseVars();
144  for (unsigned int i = 0; i < coupled_vars.size(); i++)
145  addMooseVariableDependency(coupled_vars[i]);
146  }
147 
148  if (_material_por <= 0 || _material_por >= 1)
149  mooseError("Porosity set to ", _material_por, " but it must be between 0 and 1");
150 
151  if (isCoupled("perm_change") && (coupledComponents("perm_change") != LIBMESH_DIM * LIBMESH_DIM))
152  mooseError(LIBMESH_DIM * LIBMESH_DIM,
153  " components of perm_change must be given to a RichardsMaterial. You supplied ",
154  coupledComponents("perm_change"),
155  "\n");
156 
157  _perm_change.resize(LIBMESH_DIM * LIBMESH_DIM);
158  for (unsigned int i = 0; i < LIBMESH_DIM * LIBMESH_DIM; ++i)
159  _perm_change[i] = (isCoupled("perm_change") ? &coupledValue("perm_change", i)
160  : &_zero); // coupledValue returns a reference (an
161  // alias) to a VariableValue, and the &
162  // turns it into a pointer
163 
164  if (!(_material_viscosity.size() == _num_p &&
165  getParam<std::vector<UserObjectName>>("relperm_UO").size() == _num_p &&
166  getParam<std::vector<UserObjectName>>("seff_UO").size() == _num_p &&
167  getParam<std::vector<UserObjectName>>("sat_UO").size() == _num_p &&
168  getParam<std::vector<UserObjectName>>("density_UO").size() == _num_p &&
169  getParam<std::vector<UserObjectName>>("SUPG_UO").size() == _num_p))
170  mooseError("There are ",
171  _num_p,
172  " Richards fluid variables, so you need to specify this "
173  "number of viscosities, relperm_UO, seff_UO, sat_UO, "
174  "density_UO, SUPG_UO");
175 
176  _d2density.resize(_num_p);
177  _d2rel_perm_dv.resize(_num_p);
178  _pressure_vals.resize(_num_p);
179  _pressure_old_vals.resize(_num_p);
181  _material_seff_UO.resize(_num_p);
182  _material_sat_UO.resize(_num_p);
184  _material_SUPG_UO.resize(_num_p);
185  _grad_p.resize(_num_p);
186 
187  for (unsigned int i = 0; i < _num_p; ++i)
188  {
189  // DON'T WANT "pressure_vars" at all since pp_name_UO contains the same info
190  //_pressure_vals[i] = &coupledValue("pressure_vars", i); // coupled value returns a reference
191  //_pressure_old_vals[i] = (_is_transient ? &coupledValueOld("pressure_vars", i) : &_zero);
192  //_grad_p[i] = &coupledGradient("pressure_vars", i);
193 
194  // in the following. first get the userobject names that were inputted, then get the i_th one
195  // of these, then get the actual userobject that this corresponds to, then finally & gives
196  // pointer to RichardsRelPerm object.
197  _material_relperm_UO[i] = &getUserObjectByName<RichardsRelPerm>(
198  getParam<std::vector<UserObjectName>>("relperm_UO")[i]);
199  _material_seff_UO[i] =
200  &getUserObjectByName<RichardsSeff>(getParam<std::vector<UserObjectName>>("seff_UO")[i]);
201  _material_sat_UO[i] =
202  &getUserObjectByName<RichardsSat>(getParam<std::vector<UserObjectName>>("sat_UO")[i]);
203  _material_density_UO[i] = &getUserObjectByName<RichardsDensity>(
204  getParam<std::vector<UserObjectName>>("density_UO")[i]);
205  _material_SUPG_UO[i] =
206  &getUserObjectByName<RichardsSUPG>(getParam<std::vector<UserObjectName>>("SUPG_UO")[i]);
207  }
208 }
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)
std::vector< const VariableValue * > _perm_change
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
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
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
unsigned int _num_p
MaterialProperty< std::vector< Real > > & _sat
saturation (vector of saturations in case of multiphase)
std::vector< const RichardsSUPG * > _material_SUPG_UO
unsigned int num_v() const
the number of porepressure variables
MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
d(relperm_i)/d(variable_j)
Real _material_por
porosity as entered by the user
const VariableValue & _por_change
porosity changes. if not entered they default to zero
MaterialProperty< std::vector< Real > > & _seff_old
old effective saturation
MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
d(porepressure_i)/d(variable_j)
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)
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
std::vector< const VariableValue * > _pressure_old_vals
MaterialProperty< std::vector< Real > > & _density
fluid density (or densities for multiphase problems)
MaterialProperty< std::vector< Real > > & _sat_old
old saturation
std::vector< const RichardsSeff * > _material_seff_UO
MaterialProperty< std::vector< Real > > & _mass_old
old value of fluid mass (a vector of masses for multicomponent)
MaterialProperty< Real > & _porosity
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)
std::vector< const RichardsSat * > _material_sat_UO
MaterialProperty< std::vector< Real > > & _density_old
old fluid density (or densities for multiphase problems)
MaterialProperty< Real > & _porosity_old
material properties
const VariableValue & _por_change_old
MaterialProperty< std::vector< Real > > & _seff
effective saturation (vector of effective saturations in case of multiphase)
RealTensorValue _material_perm
permeability as entered by the user
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
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)
MaterialProperty< std::vector< std::vector< Real > > > & _dmass
d(fluid mass_i)/dP_j (a vector of masses for multicomponent)
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
MaterialProperty< std::vector< Real > > & _viscosity
fluid viscosity (or viscosities in the multiphase case)

Member Function Documentation

◆ compute2ndDerivedQuantities()

void RichardsMaterial::compute2ndDerivedQuantities ( unsigned int  qp)
protected

Computes 2nd derivatives of the flux.

These are needed by kernels if doing SUPG

Parameters
qpThe quadpoint to evaluate the quantites at

Definition at line 392 of file RichardsMaterial.C.

Referenced by computeProperties().

393 {
395 
396  for (unsigned int i = 0; i < _num_p; ++i)
397  {
398  if ((*_material_SUPG_UO[i]).SUPG_trivial())
399  continue; // as the derivatives won't be needed
400 
401  // second derivative of density
402  _d2density[i].resize(_num_p);
403  Real ddens = (*_material_density_UO[i]).ddensity(_pp[qp][i]);
404  Real d2dens = (*_material_density_UO[i]).d2density(_pp[qp][i]);
405  for (unsigned int j = 0; j < _num_p; ++j)
406  {
407  _d2density[i][j].resize(_num_p);
408  for (unsigned int k = 0; k < _num_p; ++k)
409  _d2density[i][j][k] =
410  d2dens * _dpp_dv[qp][i][j] * _dpp_dv[qp][i][k] + ddens * _d2pp_dv[qp][i][j][k];
411  }
412 
413  // second derivative of relative permeability
414  _d2rel_perm_dv[i].resize(_num_p);
415  Real drel = (*_material_relperm_UO[i]).drelperm(_seff[qp][i]);
416  Real d2rel = (*_material_relperm_UO[i]).d2relperm(_seff[qp][i]);
417  for (unsigned int j = 0; j < _num_p; ++j)
418  {
419  _d2rel_perm_dv[i][j].resize(_num_p);
420  for (unsigned int k = 0; k < _num_p; ++k)
421  _d2rel_perm_dv[i][j][k] =
422  d2rel * _dseff_dv[qp][i][j] * _dseff_dv[qp][i][k] + drel * _d2seff_dv[qp][i][j][k];
423  }
424 
425  // now compute the second derivs of the fluxes
426  for (unsigned int j = 0; j < _num_p; ++j)
427  {
428  for (unsigned int k = 0; k < _num_p; ++k)
429  {
430  _d2flux_dvdv[qp][i][j][k] =
431  _d2density[i][j][k] * _rel_perm[qp][i] *
432  (_permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]));
433  _d2flux_dvdv[qp][i][j][k] +=
434  (_ddensity_dv[qp][i][j] * _drel_perm_dv[qp][i][k] +
435  _ddensity_dv[qp][i][k] * _drel_perm_dv[qp][i][j]) *
436  (_permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]));
437  _d2flux_dvdv[qp][i][j][k] +=
438  _density[qp][i] * _d2rel_perm_dv[i][j][k] *
439  (_permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]));
440  _d2flux_dvdv[qp][i][j][k] += (_ddensity_dv[qp][i][j] * _rel_perm[qp][i] +
441  _density[qp][i] * _drel_perm_dv[qp][i][j]) *
442  (_permeability[qp] * (-_ddensity_dv[qp][i][k] * _gravity[qp]));
443  _d2flux_dvdv[qp][i][j][k] += (_ddensity_dv[qp][i][k] * _rel_perm[qp][i] +
444  _density[qp][i] * _drel_perm_dv[qp][i][k]) *
445  (_permeability[qp] * (-_ddensity_dv[qp][i][j] * _gravity[qp]));
446  _d2flux_dvdv[qp][i][j][k] += _density[qp][i] * _rel_perm[qp][i] *
447  (_permeability[qp] * (-_d2density[i][j][k] * _gravity[qp]));
448  }
449  }
450  for (unsigned int j = 0; j < _num_p; ++j)
451  for (unsigned int k = 0; k < _num_p; ++k)
452  _d2flux_dvdv[qp][i][j][k] /= _viscosity[qp][i];
453 
454  for (unsigned int j = 0; j < _num_p; ++j)
455  {
456  for (unsigned int k = 0; k < _num_p; ++k)
457  {
458  _d2flux_dgradvdv[qp][i][j][k] = (_ddensity_dv[qp][i][k] * _rel_perm[qp][i] +
459  _density[qp][i] * _drel_perm_dv[qp][i][k]) *
460  _permeability[qp] * _dpp_dv[qp][i][j] / _viscosity[qp][i];
461  _d2flux_dvdgradv[qp][i][k][j] = _d2flux_dgradvdv[qp][i][j][k];
462  }
463  }
464  }
465 }
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)
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
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
unsigned int _num_p
std::vector< const RichardsSUPG * > _material_SUPG_UO
MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
d(relperm_i)/d(variable_j)
void zero2ndDerivedQuantities(unsigned int qp)
Zeroes 2nd derivatives of the flux.
MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
d(porepressure_i)/d(variable_j)
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
std::vector< const VariableGradient * > _grad_p
MaterialProperty< std::vector< Real > > & _rel_perm
relative permeability (vector of relative permeabilities in case of multiphase)
std::vector< const RichardsDensity * > _material_density_UO
MaterialProperty< std::vector< Real > > & _density
fluid density (or densities for multiphase problems)
MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _d2seff_dv
d^2(Seff_i)/d(variable_j)/d(variable_k)
MaterialProperty< std::vector< Real > > & _seff
effective saturation (vector of effective saturations in case of multiphase)
MaterialProperty< std::vector< std::vector< Real > > > & _ddensity_dv
d(density_i)/d(variable_j)
MaterialProperty< std::vector< std::vector< Real > > > & _dseff_dv
d(Seff_i)/d(variable_j)
std::vector< std::vector< std::vector< Real > > > _d2density
d^2(density)/dp_j/dP_k - used in various derivative calculations
MaterialProperty< std::vector< Real > > & _viscosity
fluid viscosity (or viscosities in the multiphase case)

◆ computeDerivedQuantities()

void RichardsMaterial::computeDerivedQuantities ( unsigned int  qp)
protected

Computes the "derived" quantities — those that depend on porepressure(s) and effective saturation(s) — such as density, relative permeability, mass, flux, etc.

Parameters
qpThe quadpoint to evaluate the quantites at

Definition at line 270 of file RichardsMaterial.C.

Referenced by computeProperties().

271 {
272  // fluid viscosity
273  _viscosity[qp].resize(_num_p);
274  for (unsigned int i = 0; i < _num_p; ++i)
275  _viscosity[qp][i] = _material_viscosity[i];
276 
277  // fluid saturation
278  _sat_old[qp].resize(_num_p);
279  _sat[qp].resize(_num_p);
280  _dsat_dv[qp].resize(_num_p);
281  for (unsigned int i = 0; i < _num_p; ++i)
282  {
283  _sat_old[qp][i] = (*_material_sat_UO[i]).sat(_seff_old[qp][i]);
284  _sat[qp][i] = (*_material_sat_UO[i]).sat(_seff[qp][i]);
285  _dsat_dv[qp][i].assign(_num_p, (*_material_sat_UO[i]).dsat(_seff[qp][i]));
286  for (unsigned int j = 0; j < _num_p; ++j)
287  _dsat_dv[qp][i][j] *= _dseff_dv[qp][i][j];
288  }
289 
290  // fluid density
291  _density_old[qp].resize(_num_p);
292  _density[qp].resize(_num_p);
293  _ddensity_dv[qp].resize(_num_p);
294  for (unsigned int i = 0; i < _num_p; ++i)
295  {
296  _density_old[qp][i] = (*_material_density_UO[i]).density(_pp_old[qp][i]);
297  _density[qp][i] = (*_material_density_UO[i]).density(_pp[qp][i]);
298  _ddensity_dv[qp][i].assign(_num_p, (*_material_density_UO[i]).ddensity(_pp[qp][i]));
299  for (unsigned int j = 0; j < _num_p; ++j)
300  _ddensity_dv[qp][i][j] *= _dpp_dv[qp][i][j];
301  }
302 
303  // relative permeability
304  _rel_perm[qp].resize(_num_p);
305  _drel_perm_dv[qp].resize(_num_p);
306  for (unsigned int i = 0; i < _num_p; ++i)
307  {
308  _rel_perm[qp][i] = (*_material_relperm_UO[i]).relperm(_seff[qp][i]);
309  _drel_perm_dv[qp][i].assign(_num_p, (*_material_relperm_UO[i]).drelperm(_seff[qp][i]));
310  for (unsigned int j = 0; j < _num_p; ++j)
311  _drel_perm_dv[qp][i][j] *= _dseff_dv[qp][i][j];
312  }
313 
314  // fluid mass
315  _mass_old[qp].resize(_num_p);
316  _mass[qp].resize(_num_p);
317  _dmass[qp].resize(_num_p);
318  for (unsigned int i = 0; i < _num_p; ++i)
319  {
320  _mass_old[qp][i] = _porosity_old[qp] * _density_old[qp][i] * _sat_old[qp][i];
321  _mass[qp][i] = _porosity[qp] * _density[qp][i] * _sat[qp][i];
322  _dmass[qp][i].resize(_num_p);
323  for (unsigned int j = 0; j < _num_p; ++j)
324  _dmass[qp][i][j] = _porosity[qp] * (_ddensity_dv[qp][i][j] * _sat[qp][i] +
325  _density[qp][i] * _dsat_dv[qp][i][j]);
326  }
327 
328  // flux without the mobility part
329  _flux_no_mob[qp].resize(_num_p);
330  _dflux_no_mob_dv[qp].resize(_num_p);
331  _dflux_no_mob_dgradv[qp].resize(_num_p);
332  for (unsigned int i = 0; i < _num_p; ++i)
333  {
334  _flux_no_mob[qp][i] = _permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]);
335 
336  _dflux_no_mob_dv[qp][i].resize(_num_p);
337  for (unsigned int j = 0; j < _num_p; ++j)
338  _dflux_no_mob_dv[qp][i][j] = _permeability[qp] * (-_ddensity_dv[qp][i][j] * _gravity[qp]);
339 
340  _dflux_no_mob_dgradv[qp][i].resize(_num_p);
341  for (unsigned int j = 0; j < _num_p; ++j)
342  _dflux_no_mob_dgradv[qp][i][j] = _permeability[qp] * _dpp_dv[qp][i][j];
343  }
344 
345  // flux
346  _flux[qp].resize(_num_p);
347  _dflux_dv[qp].resize(_num_p);
348  _dflux_dgradv[qp].resize(_num_p);
349  for (unsigned int i = 0; i < _num_p; ++i)
350  {
351  _flux[qp][i] = _density[qp][i] * _rel_perm[qp][i] * _flux_no_mob[qp][i] / _viscosity[qp][i];
352 
353  _dflux_dv[qp][i].resize(_num_p);
354  for (unsigned int j = 0; j < _num_p; ++j)
355  {
356  _dflux_dv[qp][i][j] =
357  _density[qp][i] * _rel_perm[qp][i] * _dflux_no_mob_dv[qp][i][j] / _viscosity[qp][i];
358  _dflux_dv[qp][i][j] +=
359  (_ddensity_dv[qp][i][j] * _rel_perm[qp][i] + _density[qp][i] * _drel_perm_dv[qp][i][j]) *
360  _flux_no_mob[qp][i] / _viscosity[qp][i];
361  }
362 
363  _dflux_dgradv[qp][i].resize(_num_p);
364  for (unsigned int j = 0; j < _num_p; ++j)
365  _dflux_dgradv[qp][i][j] =
366  _density[qp][i] * _rel_perm[qp][i] * _dflux_no_mob_dgradv[qp][i][j] / _viscosity[qp][i];
367  }
368 }
MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dflux_no_mob_dgradv
d(_flux_no_mob_i)/d(grad(variable_j))
MaterialProperty< std::vector< Real > > & _pp
porepressure(s)
MaterialProperty< RealVectorValue > & _gravity
std::vector< const RichardsRelPerm * > _material_relperm_UO
MaterialProperty< std::vector< Real > > & _pp_old
old values of porepressure(s)
MaterialProperty< std::vector< RealVectorValue > > & _flux
fluid flux (a vector of fluxes for multicomponent)
std::vector< Real > _material_viscosity
const std::string density
Definition: NS.h:17
unsigned int _num_p
MaterialProperty< std::vector< Real > > & _sat
saturation (vector of saturations in case of multiphase)
MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
d(relperm_i)/d(variable_j)
MaterialProperty< std::vector< Real > > & _seff_old
old effective saturation
MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
d(porepressure_i)/d(variable_j)
MaterialProperty< RealTensorValue > & _permeability
MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dflux_no_mob_dv
d(_flux_no_mob_i)/d(variable_j)
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 ...
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)
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
MaterialProperty< std::vector< Real > > & _density
fluid density (or densities for multiphase problems)
MaterialProperty< std::vector< Real > > & _sat_old
old saturation
MaterialProperty< std::vector< Real > > & _mass_old
old value of fluid mass (a vector of masses for multicomponent)
MaterialProperty< Real > & _porosity
std::vector< const RichardsSat * > _material_sat_UO
MaterialProperty< std::vector< Real > > & _density_old
old fluid density (or densities for multiphase problems)
MaterialProperty< Real > & _porosity_old
material properties
MaterialProperty< std::vector< Real > > & _seff
effective saturation (vector of effective saturations in case of multiphase)
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...
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)
MaterialProperty< std::vector< std::vector< Real > > > & _dmass
d(fluid mass_i)/dP_j (a vector of masses for multicomponent)
MaterialProperty< std::vector< Real > > & _viscosity
fluid viscosity (or viscosities in the multiphase case)

◆ computePandSeff()

void RichardsMaterial::computePandSeff ( )
protected

computes the quadpoint values of the porepressure(s) and effective saturation(s), and their derivatives wrt the variables in the system.

These are then used in computeProperties to compute relative permeability, density, and so on

Definition at line 211 of file RichardsMaterial.C.

Referenced by computeProperties().

212 {
213  // Get the pressure and effective saturation at each quadpoint
214  // From these we will build the relative permeability, density, flux, etc
215  if (_richards_name_UO.var_types() == "pppp")
216  {
217  for (unsigned int i = 0; i < _num_p; ++i)
218  {
222  }
223  }
224 
225  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
226  {
227  _pp_old[qp].resize(_num_p);
228  _pp[qp].resize(_num_p);
229  _dpp_dv[qp].resize(_num_p);
230  _d2pp_dv[qp].resize(_num_p);
231 
232  _seff_old[qp].resize(_num_p);
233  _seff[qp].resize(_num_p);
234  _dseff_dv[qp].resize(_num_p);
235  _d2seff_dv[qp].resize(_num_p);
236 
237  if (_richards_name_UO.var_types() == "pppp")
238  {
239  for (unsigned int i = 0; i < _num_p; ++i)
240  {
241  _pp_old[qp][i] = (*_pressure_old_vals[i])[qp];
242  _pp[qp][i] = (*_pressure_vals[i])[qp];
243 
244  _dpp_dv[qp][i].assign(_num_p, 0);
245  _dpp_dv[qp][i][i] = 1;
246 
247  _d2pp_dv[qp][i].resize(_num_p);
248  for (unsigned int j = 0; j < _num_p; ++j)
249  _d2pp_dv[qp][i][j].assign(_num_p, 0);
250 
251  _seff_old[qp][i] = (*_material_seff_UO[i]).seff(_pressure_old_vals, qp);
252  _seff[qp][i] = (*_material_seff_UO[i]).seff(_pressure_vals, qp);
253 
254  _dseff_dv[qp][i].resize(_num_p);
255  (*_material_seff_UO[i]).dseff(_pressure_vals, qp, _dseff_dv[qp][i]);
256 
257  _d2seff_dv[qp][i].resize(_num_p);
258  for (unsigned int j = 0; j < _num_p; ++j)
259  _d2seff_dv[qp][i][j].resize(_num_p);
260  (*_material_seff_UO[i]).d2seff(_pressure_vals, qp, _d2seff_dv[qp][i]);
261  }
262  }
263  // the above lines of code are only valid for "pppp"
264  // if you decide to code other RichardsVariables (eg "psss")
265  // you will need to add some lines here
266  }
267 }
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)
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.
unsigned int _num_p
MaterialProperty< std::vector< Real > > & _seff_old
old effective saturation
MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
d(porepressure_i)/d(variable_j)
std::string var_types() const
return the _var_types string
const VariableGradient * grad_var(unsigned int richards_var_num) const
a vector of pointers to grad(Variable)
std::vector< const VariableGradient * > _grad_p
std::vector< const VariableValue * > _pressure_vals
std::vector< const VariableValue * > _pressure_old_vals
std::vector< const RichardsSeff * > _material_seff_UO
const VariableValue * richards_vals(unsigned int richards_var_num) const
a vector of pointers to VariableValues
MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _d2seff_dv
d^2(Seff_i)/d(variable_j)/d(variable_k)
MaterialProperty< std::vector< Real > > & _seff
effective saturation (vector of effective saturations in case of multiphase)
MaterialProperty< std::vector< std::vector< Real > > > & _dseff_dv
d(Seff_i)/d(variable_j)

◆ computeProperties()

void RichardsMaterial::computeProperties ( )
protectedvirtual

Definition at line 575 of file RichardsMaterial.C.

576 {
577  // compute porepressures and effective saturations
578  // with algorithms depending on the _richards_name_UO.var_types()
579  computePandSeff();
580 
581  // porosity, permeability, and gravity
582  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
583  {
584  _porosity[qp] = _material_por + _por_change[qp];
586 
588  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
589  for (unsigned int j = 0; j < LIBMESH_DIM; j++)
590  _permeability[qp](i, j) *= std::pow(10, (*_perm_change[LIBMESH_DIM * i + j])[qp]);
591 
593  }
594 
595  // compute "derived" quantities -- those that depend on P and Seff --- such as density, relperm
596  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
598 
599  // compute certain second derivatives of the derived quantities
600  // These are needed in Jacobian calculations if doing SUPG
601  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
603 
604  // Now for SUPG itself
605  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
606  zeroSUPG(qp);
607 
608  // the following saves computational effort if all SUPG is trivial
609  bool trivial_supg = true;
610  for (unsigned int i = 0; i < _num_p; ++i)
611  trivial_supg = trivial_supg && (*_material_SUPG_UO[i]).SUPG_trivial();
612  if (trivial_supg)
613  return;
614  else
615  computeSUPG();
616 }
std::vector< const VariableValue * > _perm_change
MaterialProperty< RealVectorValue > & _gravity
unsigned int _num_p
std::vector< const RichardsSUPG * > _material_SUPG_UO
Real _material_por
porosity as entered by the user
const VariableValue & _por_change
porosity changes. if not entered they default to zero
void computePandSeff()
computes the quadpoint values of the porepressure(s) and effective saturation(s), and their derivativ...
MaterialProperty< RealTensorValue > & _permeability
RealVectorValue _material_gravity
gravity as entered by user
void computeDerivedQuantities(unsigned int qp)
Computes the "derived" quantities — those that depend on porepressure(s) and effective saturation(s)...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
MaterialProperty< Real > & _porosity
void computeSUPG()
Computes the tauvel_SUPG and its derivatives.
void compute2ndDerivedQuantities(unsigned int qp)
Computes 2nd derivatives of the flux.
void zeroSUPG(unsigned int qp)
Assigns and zeroes the MaterialProperties associated with SUPG.
MaterialProperty< Real > & _porosity_old
material properties
const VariableValue & _por_change_old
RealTensorValue _material_perm
permeability as entered by the user

◆ computeSUPG()

void RichardsMaterial::computeSUPG ( )
protected

Computes the tauvel_SUPG and its derivatives.

Definition at line 481 of file RichardsMaterial.C.

Referenced by computeProperties().

482 {
483  // Grab reference to linear Lagrange finite element object pointer,
484  // currently this is always a linear Lagrange element, so this might need to
485  // be generalized if we start working with higher-order elements...
486  FEBase *& fe(_assembly.getFE(getParam<bool>("linear_shape_fcns") ? FEType(FIRST, LAGRANGE)
487  : FEType(SECOND, LAGRANGE),
488  _current_elem->dim()));
489 
490  // Grab references to FE object's mapping data from the _subproblem's FE object
491  const std::vector<Real> & dxidx(fe->get_dxidx());
492  const std::vector<Real> & dxidy(fe->get_dxidy());
493  const std::vector<Real> & dxidz(fe->get_dxidz());
494  const std::vector<Real> & detadx(fe->get_detadx());
495  const std::vector<Real> & detady(fe->get_detady());
496  const std::vector<Real> & detadz(fe->get_detadz());
497  const std::vector<Real> & dzetadx(fe->get_dzetadx());
498  const std::vector<Real> & dzetady(fe->get_dzetady());
499  const std::vector<Real> & dzetadz(fe->get_dzetadz());
500 
501  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
502  {
503 
504  // Bounds checking on element data and putting into vector form
505  mooseAssert(qp < dxidx.size(), "Insufficient data in dxidx array!");
506  mooseAssert(qp < dxidy.size(), "Insufficient data in dxidy array!");
507  mooseAssert(qp < dxidz.size(), "Insufficient data in dxidz array!");
508  if (_mesh.dimension() >= 2)
509  {
510  mooseAssert(qp < detadx.size(), "Insufficient data in detadx array!");
511  mooseAssert(qp < detady.size(), "Insufficient data in detady array!");
512  mooseAssert(qp < detadz.size(), "Insufficient data in detadz array!");
513  }
514  if (_mesh.dimension() >= 3)
515  {
516  mooseAssert(qp < dzetadx.size(), "Insufficient data in dzetadx array!");
517  mooseAssert(qp < dzetady.size(), "Insufficient data in dzetady array!");
518  mooseAssert(qp < dzetadz.size(), "Insufficient data in dzetadz array!");
519  }
520 
521  // CHECK : Does this work spherical, cylindrical, etc?
522  RealVectorValue xi_prime(dxidx[qp], dxidy[qp], dxidz[qp]);
523  RealVectorValue eta_prime, zeta_prime;
524  if (_mesh.dimension() >= 2)
525  {
526  eta_prime(0) = detadx[qp];
527  eta_prime(1) = detady[qp];
528  }
529  if (_mesh.dimension() == 3)
530  {
531  eta_prime(2) = detadz[qp];
532  zeta_prime(0) = dzetadx[qp];
533  zeta_prime(1) = dzetady[qp];
534  zeta_prime(2) = dzetadz[qp];
535  }
536 
537  _trace_perm = _permeability[qp].tr();
538  for (unsigned int i = 0; i < _num_p; ++i)
539  {
540  RealVectorValue vel =
541  (*_material_SUPG_UO[i])
542  .velSUPG(_permeability[qp], (*_grad_p[i])[qp], _density[qp][i], _gravity[qp]);
543  RealTensorValue dvel_dgradp = (*_material_SUPG_UO[i]).dvelSUPG_dgradp(_permeability[qp]);
544  RealVectorValue dvel_dp =
545  (*_material_SUPG_UO[i])
546  .dvelSUPG_dp(_permeability[qp], _ddensity_dv[qp][i][i], _gravity[qp]);
547  RealVectorValue bb =
548  (*_material_SUPG_UO[i]).bb(vel, _mesh.dimension(), xi_prime, eta_prime, zeta_prime);
549  RealVectorValue dbb2_dgradp =
550  (*_material_SUPG_UO[i]).dbb2_dgradp(vel, dvel_dgradp, xi_prime, eta_prime, zeta_prime);
551  Real dbb2_dp = (*_material_SUPG_UO[i]).dbb2_dp(vel, dvel_dp, xi_prime, eta_prime, zeta_prime);
552  Real tau = (*_material_SUPG_UO[i]).tauSUPG(vel, _trace_perm, bb);
553  RealVectorValue dtau_dgradp =
554  (*_material_SUPG_UO[i]).dtauSUPG_dgradp(vel, dvel_dgradp, _trace_perm, bb, dbb2_dgradp);
555  Real dtau_dp = (*_material_SUPG_UO[i]).dtauSUPG_dp(vel, dvel_dp, _trace_perm, bb, dbb2_dp);
556 
557  _tauvel_SUPG[qp][i] = tau * vel;
558 
559  RealTensorValue dtauvel_dgradp = tau * dvel_dgradp;
560  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
561  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
562  dtauvel_dgradp(j, k) +=
563  dtau_dgradp(j) * vel(k); // this is outerproduct - maybe libmesh can do it better?
564  for (unsigned int j = 0; j < _num_p; ++j)
565  _dtauvel_SUPG_dgradp[qp][i][j] = dtauvel_dgradp * _dpp_dv[qp][i][j];
566 
567  RealVectorValue dtauvel_dp = dtau_dp * vel + tau * dvel_dp;
568  for (unsigned int j = 0; j < _num_p; ++j)
569  _dtauvel_SUPG_dp[qp][i][j] = dtauvel_dp * _dpp_dv[qp][i][j];
570  }
571  }
572 }
MaterialProperty< RealVectorValue > & _gravity
unsigned int _num_p
std::vector< const RichardsSUPG * > _material_SUPG_UO
MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
d(porepressure_i)/d(variable_j)
MaterialProperty< RealTensorValue > & _permeability
MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dtauvel_SUPG_dgradp
std::vector< const VariableGradient * > _grad_p
MaterialProperty< std::vector< Real > > & _density
fluid density (or densities for multiphase problems)
MaterialProperty< std::vector< RealVectorValue > > & _tauvel_SUPG
Real _trace_perm
trace of permeability tensor
MaterialProperty< std::vector< std::vector< Real > > > & _ddensity_dv
d(density_i)/d(variable_j)
MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dtauvel_SUPG_dp

◆ zero2ndDerivedQuantities()

void RichardsMaterial::zero2ndDerivedQuantities ( unsigned int  qp)
private

Zeroes 2nd derivatives of the flux.

The values are only calculated in compute2ndDerivedQuantites if the SUPG UserObject says they are need. This is for computational efficiency as these are very expensive to calculate

Parameters
qpThe quadpoint to evaluate the quantites at

Definition at line 371 of file RichardsMaterial.C.

Referenced by compute2ndDerivedQuantities().

372 {
373  _d2flux_dvdv[qp].resize(_num_p);
374  _d2flux_dgradvdv[qp].resize(_num_p);
375  _d2flux_dvdgradv[qp].resize(_num_p);
376 
377  for (unsigned int i = 0; i < _num_p; ++i)
378  {
379  _d2flux_dvdv[qp][i].resize(_num_p);
380  _d2flux_dgradvdv[qp][i].resize(_num_p);
381  _d2flux_dvdgradv[qp][i].resize(_num_p);
382  for (unsigned int j = 0; j < _num_p; ++j)
383  {
384  _d2flux_dvdv[qp][i][j].assign(_num_p, RealVectorValue());
385  _d2flux_dgradvdv[qp][i][j].assign(_num_p, RealTensorValue());
386  _d2flux_dvdgradv[qp][i][j].assign(_num_p, RealTensorValue());
387  }
388  }
389 }
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.
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
unsigned int _num_p
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...

◆ zeroSUPG()

void RichardsMaterial::zeroSUPG ( unsigned int  qp)
protected

Assigns and zeroes the MaterialProperties associated with SUPG.

Parameters
qpThe quadpoint to assign+zero at

Definition at line 468 of file RichardsMaterial.C.

Referenced by computeProperties().

469 {
470  _tauvel_SUPG[qp].assign(_num_p, RealVectorValue());
471  _dtauvel_SUPG_dgradp[qp].resize(_num_p);
472  _dtauvel_SUPG_dp[qp].resize(_num_p);
473  for (unsigned int i = 0; i < _num_p; ++i)
474  {
475  _dtauvel_SUPG_dp[qp][i].assign(_num_p, RealVectorValue());
476  _dtauvel_SUPG_dgradp[qp][i].assign(_num_p, RealTensorValue());
477  }
478 }
unsigned int _num_p
MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dtauvel_SUPG_dgradp
MaterialProperty< std::vector< RealVectorValue > > & _tauvel_SUPG
MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dtauvel_SUPG_dp

Member Data Documentation

◆ _d2density

std::vector<std::vector<std::vector<Real> > > RichardsMaterial::_d2density
private

d^2(density)/dp_j/dP_k - used in various derivative calculations

Definition at line 201 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and RichardsMaterial().

◆ _d2flux_dgradvdv

MaterialProperty<std::vector<std::vector<std::vector<RealTensorValue> > > >& RichardsMaterial::_d2flux_dgradvdv
private

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

Definition at line 189 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and zero2ndDerivedQuantities().

◆ _d2flux_dvdgradv

MaterialProperty<std::vector<std::vector<std::vector<RealTensorValue> > > >& RichardsMaterial::_d2flux_dvdgradv
private

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.

Definition at line 192 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and zero2ndDerivedQuantities().

◆ _d2flux_dvdv

MaterialProperty<std::vector<std::vector<std::vector<RealVectorValue> > > >& RichardsMaterial::_d2flux_dvdv
private

d^2(Richards flux_i)/d(variable_j)/d(variable_k), here flux_i is the i_th flux, which is itself a RealVectorValue

Definition at line 186 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and zero2ndDerivedQuantities().

◆ _d2pp_dv

MaterialProperty<std::vector<std::vector<std::vector<Real> > > >& RichardsMaterial::_d2pp_dv
private

d^2(porepressure_i)/d(variable_j)/d(variable_k)

Definition at line 117 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computePandSeff().

◆ _d2rel_perm_dv

std::vector<std::vector<std::vector<Real> > > RichardsMaterial::_d2rel_perm_dv
private

d^2(relperm_i)/dP_j/dP_k - used in various derivative calculations

Definition at line 204 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and RichardsMaterial().

◆ _d2seff_dv

MaterialProperty<std::vector<std::vector<std::vector<Real> > > >& RichardsMaterial::_d2seff_dv
private

d^2(Seff_i)/d(variable_j)/d(variable_k)

Definition at line 141 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computePandSeff().

◆ _ddensity_dv

MaterialProperty<std::vector<std::vector<Real> > >& RichardsMaterial::_ddensity_dv
private

d(density_i)/d(variable_j)

Definition at line 129 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), computeDerivedQuantities(), and computeSUPG().

◆ _density

MaterialProperty<std::vector<Real> >& RichardsMaterial::_density
private

fluid density (or densities for multiphase problems)

Definition at line 126 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), computeDerivedQuantities(), and computeSUPG().

◆ _density_old

MaterialProperty<std::vector<Real> >& RichardsMaterial::_density_old
private

old fluid density (or densities for multiphase problems)

Definition at line 123 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _dflux_dgradv

MaterialProperty<std::vector<std::vector<RealTensorValue> > >& RichardsMaterial::_dflux_dgradv
private

d(Richards flux_i)/d(grad(variable_j)), here flux_i is the i_th flux, which is itself a RealVectorValue

Definition at line 183 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _dflux_dv

MaterialProperty<std::vector<std::vector<RealVectorValue> > >& RichardsMaterial::_dflux_dv
private

d(Richards flux_i)/d(variable_j), here flux_i is the i_th flux, which is itself a RealVectorValue

Definition at line 180 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _dflux_no_mob_dgradv

MaterialProperty<std::vector<std::vector<RealTensorValue> > >& RichardsMaterial::_dflux_no_mob_dgradv
private

d(_flux_no_mob_i)/d(grad(variable_j))

Definition at line 174 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _dflux_no_mob_dv

MaterialProperty<std::vector<std::vector<RealVectorValue> > >& RichardsMaterial::_dflux_no_mob_dv
private

d(_flux_no_mob_i)/d(variable_j)

Definition at line 171 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _dmass

MaterialProperty<std::vector<std::vector<Real> > >& RichardsMaterial::_dmass
private

d(fluid mass_i)/dP_j (a vector of masses for multicomponent)

Definition at line 165 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _dpp_dv

MaterialProperty<std::vector<std::vector<Real> > >& RichardsMaterial::_dpp_dv
private

d(porepressure_i)/d(variable_j)

Definition at line 114 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), computeDerivedQuantities(), computePandSeff(), and computeSUPG().

◆ _drel_perm_dv

MaterialProperty<std::vector<std::vector<Real> > >& RichardsMaterial::_drel_perm_dv
private

d(relperm_i)/d(variable_j)

Definition at line 156 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computeDerivedQuantities().

◆ _dsat_dv

MaterialProperty<std::vector<std::vector<Real> > >& RichardsMaterial::_dsat_dv
private

d(saturation_i)/d(variable_j)

Definition at line 150 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _dseff_dv

MaterialProperty<std::vector<std::vector<Real> > >& RichardsMaterial::_dseff_dv
private

d(Seff_i)/d(variable_j)

Definition at line 138 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), computeDerivedQuantities(), and computePandSeff().

◆ _dtauvel_SUPG_dgradp

MaterialProperty<std::vector<std::vector<RealTensorValue> > >& RichardsMaterial::_dtauvel_SUPG_dgradp
private

Definition at line 196 of file RichardsMaterial.h.

Referenced by computeSUPG(), and zeroSUPG().

◆ _dtauvel_SUPG_dp

MaterialProperty<std::vector<std::vector<RealVectorValue> > >& RichardsMaterial::_dtauvel_SUPG_dp
private

Definition at line 198 of file RichardsMaterial.h.

Referenced by computeSUPG(), and zeroSUPG().

◆ _flux

MaterialProperty<std::vector<RealVectorValue> >& RichardsMaterial::_flux
private

fluid flux (a vector of fluxes for multicomponent)

Definition at line 177 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _flux_no_mob

MaterialProperty<std::vector<RealVectorValue> >& RichardsMaterial::_flux_no_mob
private

permeability*(grad(P) - density*gravity) (a vector of these for multicomponent)

Definition at line 168 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _grad_p

std::vector<const VariableGradient *> RichardsMaterial::_grad_p
private

◆ _gravity

MaterialProperty<RealVectorValue>& RichardsMaterial::_gravity
protected

◆ _mass

MaterialProperty<std::vector<Real> >& RichardsMaterial::_mass
private

fluid mass (a vector of masses for multicomponent)

Definition at line 162 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _mass_old

MaterialProperty<std::vector<Real> >& RichardsMaterial::_mass_old
private

old value of fluid mass (a vector of masses for multicomponent)

Definition at line 159 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _material_density_UO

std::vector<const RichardsDensity *> RichardsMaterial::_material_density_UO
protected

◆ _material_gravity

RealVectorValue RichardsMaterial::_material_gravity
protected

gravity as entered by user

Definition at line 45 of file RichardsMaterial.h.

Referenced by computeProperties().

◆ _material_perm

RealTensorValue RichardsMaterial::_material_perm
protected

permeability as entered by the user

Definition at line 42 of file RichardsMaterial.h.

Referenced by computeProperties().

◆ _material_por

Real RichardsMaterial::_material_por
protected

porosity as entered by the user

Definition at line 35 of file RichardsMaterial.h.

Referenced by computeProperties(), and RichardsMaterial().

◆ _material_relperm_UO

std::vector<const RichardsRelPerm *> RichardsMaterial::_material_relperm_UO
protected

◆ _material_sat_UO

std::vector<const RichardsSat *> RichardsMaterial::_material_sat_UO
protected

Definition at line 59 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and RichardsMaterial().

◆ _material_seff_UO

std::vector<const RichardsSeff *> RichardsMaterial::_material_seff_UO
protected

Definition at line 58 of file RichardsMaterial.h.

Referenced by computePandSeff(), and RichardsMaterial().

◆ _material_SUPG_UO

std::vector<const RichardsSUPG *> RichardsMaterial::_material_SUPG_UO
protected

◆ _material_viscosity

std::vector<Real> RichardsMaterial::_material_viscosity
private

Definition at line 105 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and RichardsMaterial().

◆ _num_p

unsigned int RichardsMaterial::_num_p
protected

◆ _perm_change

std::vector<const VariableValue *> RichardsMaterial::_perm_change
protected

Definition at line 63 of file RichardsMaterial.h.

Referenced by computeProperties(), and RichardsMaterial().

◆ _permeability

MaterialProperty<RealTensorValue>& RichardsMaterial::_permeability
protected

◆ _por_change

const VariableValue& RichardsMaterial::_por_change
protected

porosity changes. if not entered they default to zero

Definition at line 38 of file RichardsMaterial.h.

Referenced by computeProperties().

◆ _por_change_old

const VariableValue& RichardsMaterial::_por_change_old
protected

Definition at line 39 of file RichardsMaterial.h.

Referenced by computeProperties().

◆ _porosity

MaterialProperty<Real>& RichardsMaterial::_porosity
protected

Definition at line 49 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and computeProperties().

◆ _porosity_old

MaterialProperty<Real>& RichardsMaterial::_porosity_old
protected

material properties

Definition at line 48 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and computeProperties().

◆ _pp

MaterialProperty<std::vector<Real> >& RichardsMaterial::_pp
private

porepressure(s)

Definition at line 111 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), computeDerivedQuantities(), and computePandSeff().

◆ _pp_old

MaterialProperty<std::vector<Real> >& RichardsMaterial::_pp_old
private

old values of porepressure(s)

Definition at line 108 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and computePandSeff().

◆ _pressure_old_vals

std::vector<const VariableValue *> RichardsMaterial::_pressure_old_vals
private

Definition at line 207 of file RichardsMaterial.h.

Referenced by computePandSeff(), and RichardsMaterial().

◆ _pressure_vals

std::vector<const VariableValue *> RichardsMaterial::_pressure_vals
private

Definition at line 206 of file RichardsMaterial.h.

Referenced by computePandSeff(), and RichardsMaterial().

◆ _rel_perm

MaterialProperty<std::vector<Real> >& RichardsMaterial::_rel_perm
private

relative permeability (vector of relative permeabilities in case of multiphase)

Definition at line 153 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computeDerivedQuantities().

◆ _richards_name_UO

const RichardsVarNames& RichardsMaterial::_richards_name_UO
protected

The variable names userobject for the Richards variables.

Definition at line 54 of file RichardsMaterial.h.

Referenced by computePandSeff(), and RichardsMaterial().

◆ _sat

MaterialProperty<std::vector<Real> >& RichardsMaterial::_sat
private

saturation (vector of saturations in case of multiphase)

Definition at line 147 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _sat_old

MaterialProperty<std::vector<Real> >& RichardsMaterial::_sat_old
private

old saturation

Definition at line 144 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

◆ _seff

MaterialProperty<std::vector<Real> >& RichardsMaterial::_seff
private

effective saturation (vector of effective saturations in case of multiphase)

Definition at line 135 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), computeDerivedQuantities(), and computePandSeff().

◆ _seff_old

MaterialProperty<std::vector<Real> >& RichardsMaterial::_seff_old
private

old effective saturation

Definition at line 132 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and computePandSeff().

◆ _tauvel_SUPG

MaterialProperty<std::vector<RealVectorValue> >& RichardsMaterial::_tauvel_SUPG
private

Definition at line 194 of file RichardsMaterial.h.

Referenced by computeSUPG(), and zeroSUPG().

◆ _trace_perm

Real RichardsMaterial::_trace_perm
private

trace of permeability tensor

Definition at line 103 of file RichardsMaterial.h.

Referenced by computeSUPG().

◆ _viscosity

MaterialProperty<std::vector<Real> >& RichardsMaterial::_viscosity
private

fluid viscosity (or viscosities in the multiphase case)

Definition at line 120 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computeDerivedQuantities().


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