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 26 of file RichardsMaterial.h.

Constructor & Destructor Documentation

RichardsMaterial::RichardsMaterial ( const InputParameters &  parameters)

Definition at line 63 of file RichardsMaterial.C.

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

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 387 of file RichardsMaterial.C.

Referenced by computeProperties().

388 {
390 
391  for (unsigned int i = 0; i < _num_p; ++i)
392  {
393  if ((*_material_SUPG_UO[i]).SUPG_trivial())
394  continue; // as the derivatives won't be needed
395 
396  // second derivative of density
397  _d2density[i].resize(_num_p);
398  Real ddens = (*_material_density_UO[i]).ddensity(_pp[qp][i]);
399  Real d2dens = (*_material_density_UO[i]).d2density(_pp[qp][i]);
400  for (unsigned int j = 0; j < _num_p; ++j)
401  {
402  _d2density[i][j].resize(_num_p);
403  for (unsigned int k = 0; k < _num_p; ++k)
404  _d2density[i][j][k] =
405  d2dens * _dpp_dv[qp][i][j] * _dpp_dv[qp][i][k] + ddens * _d2pp_dv[qp][i][j][k];
406  }
407 
408  // second derivative of relative permeability
409  _d2rel_perm_dv[i].resize(_num_p);
410  Real drel = (*_material_relperm_UO[i]).drelperm(_seff[qp][i]);
411  Real d2rel = (*_material_relperm_UO[i]).d2relperm(_seff[qp][i]);
412  for (unsigned int j = 0; j < _num_p; ++j)
413  {
414  _d2rel_perm_dv[i][j].resize(_num_p);
415  for (unsigned int k = 0; k < _num_p; ++k)
416  _d2rel_perm_dv[i][j][k] =
417  d2rel * _dseff_dv[qp][i][j] * _dseff_dv[qp][i][k] + drel * _d2seff_dv[qp][i][j][k];
418  }
419 
420  // now compute the second derivs of the fluxes
421  for (unsigned int j = 0; j < _num_p; ++j)
422  {
423  for (unsigned int k = 0; k < _num_p; ++k)
424  {
425  _d2flux_dvdv[qp][i][j][k] =
426  _d2density[i][j][k] * _rel_perm[qp][i] *
427  (_permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]));
428  _d2flux_dvdv[qp][i][j][k] +=
429  (_ddensity_dv[qp][i][j] * _drel_perm_dv[qp][i][k] +
430  _ddensity_dv[qp][i][k] * _drel_perm_dv[qp][i][j]) *
431  (_permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]));
432  _d2flux_dvdv[qp][i][j][k] +=
433  _density[qp][i] * _d2rel_perm_dv[i][j][k] *
434  (_permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]));
435  _d2flux_dvdv[qp][i][j][k] += (_ddensity_dv[qp][i][j] * _rel_perm[qp][i] +
436  _density[qp][i] * _drel_perm_dv[qp][i][j]) *
437  (_permeability[qp] * (-_ddensity_dv[qp][i][k] * _gravity[qp]));
438  _d2flux_dvdv[qp][i][j][k] += (_ddensity_dv[qp][i][k] * _rel_perm[qp][i] +
439  _density[qp][i] * _drel_perm_dv[qp][i][k]) *
440  (_permeability[qp] * (-_ddensity_dv[qp][i][j] * _gravity[qp]));
441  _d2flux_dvdv[qp][i][j][k] += _density[qp][i] * _rel_perm[qp][i] *
442  (_permeability[qp] * (-_d2density[i][j][k] * _gravity[qp]));
443  }
444  }
445  for (unsigned int j = 0; j < _num_p; ++j)
446  for (unsigned int k = 0; k < _num_p; ++k)
447  _d2flux_dvdv[qp][i][j][k] /= _viscosity[qp][i];
448 
449  for (unsigned int j = 0; j < _num_p; ++j)
450  {
451  for (unsigned int k = 0; k < _num_p; ++k)
452  {
453  _d2flux_dgradvdv[qp][i][j][k] = (_ddensity_dv[qp][i][k] * _rel_perm[qp][i] +
454  _density[qp][i] * _drel_perm_dv[qp][i][k]) *
455  _permeability[qp] * _dpp_dv[qp][i][j] / _viscosity[qp][i];
456  _d2flux_dvdgradv[qp][i][k][j] = _d2flux_dgradvdv[qp][i][j][k];
457  }
458  }
459  }
460 }
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)
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 265 of file RichardsMaterial.C.

Referenced by computeProperties().

266 {
267  // fluid viscosity
268  _viscosity[qp].resize(_num_p);
269  for (unsigned int i = 0; i < _num_p; ++i)
270  _viscosity[qp][i] = _material_viscosity[i];
271 
272  // fluid saturation
273  _sat_old[qp].resize(_num_p);
274  _sat[qp].resize(_num_p);
275  _dsat_dv[qp].resize(_num_p);
276  for (unsigned int i = 0; i < _num_p; ++i)
277  {
278  _sat_old[qp][i] = (*_material_sat_UO[i]).sat(_seff_old[qp][i]);
279  _sat[qp][i] = (*_material_sat_UO[i]).sat(_seff[qp][i]);
280  _dsat_dv[qp][i].assign(_num_p, (*_material_sat_UO[i]).dsat(_seff[qp][i]));
281  for (unsigned int j = 0; j < _num_p; ++j)
282  _dsat_dv[qp][i][j] *= _dseff_dv[qp][i][j];
283  }
284 
285  // fluid density
286  _density_old[qp].resize(_num_p);
287  _density[qp].resize(_num_p);
288  _ddensity_dv[qp].resize(_num_p);
289  for (unsigned int i = 0; i < _num_p; ++i)
290  {
291  _density_old[qp][i] = (*_material_density_UO[i]).density(_pp_old[qp][i]);
292  _density[qp][i] = (*_material_density_UO[i]).density(_pp[qp][i]);
293  _ddensity_dv[qp][i].assign(_num_p, (*_material_density_UO[i]).ddensity(_pp[qp][i]));
294  for (unsigned int j = 0; j < _num_p; ++j)
295  _ddensity_dv[qp][i][j] *= _dpp_dv[qp][i][j];
296  }
297 
298  // relative permeability
299  _rel_perm[qp].resize(_num_p);
300  _drel_perm_dv[qp].resize(_num_p);
301  for (unsigned int i = 0; i < _num_p; ++i)
302  {
303  _rel_perm[qp][i] = (*_material_relperm_UO[i]).relperm(_seff[qp][i]);
304  _drel_perm_dv[qp][i].assign(_num_p, (*_material_relperm_UO[i]).drelperm(_seff[qp][i]));
305  for (unsigned int j = 0; j < _num_p; ++j)
306  _drel_perm_dv[qp][i][j] *= _dseff_dv[qp][i][j];
307  }
308 
309  // fluid mass
310  _mass_old[qp].resize(_num_p);
311  _mass[qp].resize(_num_p);
312  _dmass[qp].resize(_num_p);
313  for (unsigned int i = 0; i < _num_p; ++i)
314  {
315  _mass_old[qp][i] = _porosity_old[qp] * _density_old[qp][i] * _sat_old[qp][i];
316  _mass[qp][i] = _porosity[qp] * _density[qp][i] * _sat[qp][i];
317  _dmass[qp][i].resize(_num_p);
318  for (unsigned int j = 0; j < _num_p; ++j)
319  _dmass[qp][i][j] = _porosity[qp] * (_ddensity_dv[qp][i][j] * _sat[qp][i] +
320  _density[qp][i] * _dsat_dv[qp][i][j]);
321  }
322 
323  // flux without the mobility part
324  _flux_no_mob[qp].resize(_num_p);
325  _dflux_no_mob_dv[qp].resize(_num_p);
326  _dflux_no_mob_dgradv[qp].resize(_num_p);
327  for (unsigned int i = 0; i < _num_p; ++i)
328  {
329  _flux_no_mob[qp][i] = _permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]);
330 
331  _dflux_no_mob_dv[qp][i].resize(_num_p);
332  for (unsigned int j = 0; j < _num_p; ++j)
333  _dflux_no_mob_dv[qp][i][j] = _permeability[qp] * (-_ddensity_dv[qp][i][j] * _gravity[qp]);
334 
335  _dflux_no_mob_dgradv[qp][i].resize(_num_p);
336  for (unsigned int j = 0; j < _num_p; ++j)
337  _dflux_no_mob_dgradv[qp][i][j] = _permeability[qp] * _dpp_dv[qp][i][j];
338  }
339 
340  // flux
341  _flux[qp].resize(_num_p);
342  _dflux_dv[qp].resize(_num_p);
343  _dflux_dgradv[qp].resize(_num_p);
344  for (unsigned int i = 0; i < _num_p; ++i)
345  {
346  _flux[qp][i] = _density[qp][i] * _rel_perm[qp][i] * _flux_no_mob[qp][i] / _viscosity[qp][i];
347 
348  _dflux_dv[qp][i].resize(_num_p);
349  for (unsigned int j = 0; j < _num_p; ++j)
350  {
351  _dflux_dv[qp][i][j] =
352  _density[qp][i] * _rel_perm[qp][i] * _dflux_no_mob_dv[qp][i][j] / _viscosity[qp][i];
353  _dflux_dv[qp][i][j] +=
354  (_ddensity_dv[qp][i][j] * _rel_perm[qp][i] + _density[qp][i] * _drel_perm_dv[qp][i][j]) *
355  _flux_no_mob[qp][i] / _viscosity[qp][i];
356  }
357 
358  _dflux_dgradv[qp][i].resize(_num_p);
359  for (unsigned int j = 0; j < _num_p; ++j)
360  _dflux_dgradv[qp][i][j] =
361  _density[qp][i] * _rel_perm[qp][i] * _dflux_no_mob_dgradv[qp][i][j] / _viscosity[qp][i];
362  }
363 }
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:15
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)
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 206 of file RichardsMaterial.C.

Referenced by computeProperties().

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

Definition at line 570 of file RichardsMaterial.C.

571 {
572  // compute porepressures and effective saturations
573  // with algorithms depending on the _richards_name_UO.var_types()
574  computePandSeff();
575 
576  // porosity, permeability, and gravity
577  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
578  {
579  _porosity[qp] = _material_por + _por_change[qp];
581 
583  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
584  for (unsigned int j = 0; j < LIBMESH_DIM; j++)
585  _permeability[qp](i, j) *= std::pow(10, (*_perm_change[LIBMESH_DIM * i + j])[qp]);
586 
588  }
589 
590  // compute "derived" quantities -- those that depend on P and Seff --- such as density, relperm
591  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
593 
594  // compute certain second derivatives of the derived quantities
595  // These are needed in Jacobian calculations if doing SUPG
596  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
598 
599  // Now for SUPG itself
600  for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
601  zeroSUPG(qp);
602 
603  // the following saves computational effort if all SUPG is trivial
604  bool trivial_supg = true;
605  for (unsigned int i = 0; i < _num_p; ++i)
606  trivial_supg = trivial_supg && (*_material_SUPG_UO[i]).SUPG_trivial();
607  if (trivial_supg)
608  return;
609  else
610  computeSUPG();
611 }
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
void RichardsMaterial::computeSUPG ( )
protected

Computes the tauvel_SUPG and its derivatives.

Definition at line 476 of file RichardsMaterial.C.

Referenced by computeProperties().

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

Referenced by compute2ndDerivedQuantities().

367 {
368  _d2flux_dvdv[qp].resize(_num_p);
369  _d2flux_dgradvdv[qp].resize(_num_p);
370  _d2flux_dvdgradv[qp].resize(_num_p);
371 
372  for (unsigned int i = 0; i < _num_p; ++i)
373  {
374  _d2flux_dvdv[qp][i].resize(_num_p);
375  _d2flux_dgradvdv[qp][i].resize(_num_p);
376  _d2flux_dvdgradv[qp][i].resize(_num_p);
377  for (unsigned int j = 0; j < _num_p; ++j)
378  {
379  _d2flux_dvdv[qp][i][j].assign(_num_p, RealVectorValue());
380  _d2flux_dgradvdv[qp][i][j].assign(_num_p, RealTensorValue());
381  _d2flux_dvdgradv[qp][i][j].assign(_num_p, RealTensorValue());
382  }
383  }
384 }
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...
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 463 of file RichardsMaterial.C.

Referenced by computeProperties().

464 {
465  _tauvel_SUPG[qp].assign(_num_p, RealVectorValue());
466  _dtauvel_SUPG_dgradp[qp].resize(_num_p);
467  _dtauvel_SUPG_dp[qp].resize(_num_p);
468  for (unsigned int i = 0; i < _num_p; ++i)
469  {
470  _dtauvel_SUPG_dp[qp][i].assign(_num_p, RealVectorValue());
471  _dtauvel_SUPG_dgradp[qp][i].assign(_num_p, RealTensorValue());
472  }
473 }
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

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 199 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and RichardsMaterial().

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 187 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and zero2ndDerivedQuantities().

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 190 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and zero2ndDerivedQuantities().

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 184 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and zero2ndDerivedQuantities().

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 115 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computePandSeff().

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 202 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and RichardsMaterial().

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 139 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computePandSeff().

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

d(density_i)/d(variable_j)

Definition at line 127 of file RichardsMaterial.h.

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

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

fluid density (or densities for multiphase problems)

Definition at line 124 of file RichardsMaterial.h.

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

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

old fluid density (or densities for multiphase problems)

Definition at line 121 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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 181 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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 178 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

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

Definition at line 172 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

d(_flux_no_mob_i)/d(variable_j)

Definition at line 169 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

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

Definition at line 163 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

d(porepressure_i)/d(variable_j)

Definition at line 112 of file RichardsMaterial.h.

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

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

d(relperm_i)/d(variable_j)

Definition at line 154 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computeDerivedQuantities().

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

d(saturation_i)/d(variable_j)

Definition at line 148 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

d(Seff_i)/d(variable_j)

Definition at line 136 of file RichardsMaterial.h.

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

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

Definition at line 194 of file RichardsMaterial.h.

Referenced by computeSUPG(), and zeroSUPG().

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

Definition at line 196 of file RichardsMaterial.h.

Referenced by computeSUPG(), and zeroSUPG().

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

fluid flux (a vector of fluxes for multicomponent)

Definition at line 175 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

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

Definition at line 166 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

std::vector<const VariableGradient *> RichardsMaterial::_grad_p
private
MaterialProperty<RealVectorValue>& RichardsMaterial::_gravity
protected
MaterialProperty<std::vector<Real> >& RichardsMaterial::_mass
private

fluid mass (a vector of masses for multicomponent)

Definition at line 160 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

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

Definition at line 157 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

gravity as entered by user

Definition at line 43 of file RichardsMaterial.h.

Referenced by computeProperties().

RealTensorValue RichardsMaterial::_material_perm
protected

permeability as entered by the user

Definition at line 40 of file RichardsMaterial.h.

Referenced by computeProperties().

Real RichardsMaterial::_material_por
protected

porosity as entered by the user

Definition at line 33 of file RichardsMaterial.h.

Referenced by computeProperties(), and RichardsMaterial().

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

Definition at line 57 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and RichardsMaterial().

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

Definition at line 56 of file RichardsMaterial.h.

Referenced by computePandSeff(), and RichardsMaterial().

std::vector<const RichardsSUPG *> RichardsMaterial::_material_SUPG_UO
protected
std::vector<Real> RichardsMaterial::_material_viscosity
private

Definition at line 103 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and RichardsMaterial().

unsigned int RichardsMaterial::_num_p
protected
std::vector<const VariableValue *> RichardsMaterial::_perm_change
protected

Definition at line 61 of file RichardsMaterial.h.

Referenced by computeProperties(), and RichardsMaterial().

MaterialProperty<RealTensorValue>& RichardsMaterial::_permeability
protected
const VariableValue& RichardsMaterial::_por_change
protected

porosity changes. if not entered they default to zero

Definition at line 36 of file RichardsMaterial.h.

Referenced by computeProperties().

const VariableValue& RichardsMaterial::_por_change_old
protected

Definition at line 37 of file RichardsMaterial.h.

Referenced by computeProperties().

MaterialProperty<Real>& RichardsMaterial::_porosity
protected

Definition at line 47 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and computeProperties().

MaterialProperty<Real>& RichardsMaterial::_porosity_old
protected

material properties

Definition at line 46 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and computeProperties().

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

porepressure(s)

Definition at line 109 of file RichardsMaterial.h.

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

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

old values of porepressure(s)

Definition at line 106 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and computePandSeff().

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

Definition at line 205 of file RichardsMaterial.h.

Referenced by computePandSeff(), and RichardsMaterial().

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

Definition at line 204 of file RichardsMaterial.h.

Referenced by computePandSeff(), and RichardsMaterial().

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

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

Definition at line 151 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computeDerivedQuantities().

const RichardsVarNames& RichardsMaterial::_richards_name_UO
protected

The variable names userobject for the Richards variables.

Definition at line 52 of file RichardsMaterial.h.

Referenced by computePandSeff(), and RichardsMaterial().

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

saturation (vector of saturations in case of multiphase)

Definition at line 145 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

old saturation

Definition at line 142 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities().

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

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

Definition at line 133 of file RichardsMaterial.h.

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

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

old effective saturation

Definition at line 130 of file RichardsMaterial.h.

Referenced by computeDerivedQuantities(), and computePandSeff().

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

Definition at line 192 of file RichardsMaterial.h.

Referenced by computeSUPG(), and zeroSUPG().

Real RichardsMaterial::_trace_perm
private

trace of permeability tensor

Definition at line 101 of file RichardsMaterial.h.

Referenced by computeSUPG().

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

fluid viscosity (or viscosities in the multiphase case)

Definition at line 118 of file RichardsMaterial.h.

Referenced by compute2ndDerivedQuantities(), and computeDerivedQuantities().


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