www.mooseframework.org
RichardsMaterial.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "RichardsMaterial.h"
9 #include "Assembly.h"
10 #include "MooseMesh.h"
11 
12 #include <cmath> // std::sinh and std::cosh
13 
14 #include "libmesh/quadrature.h"
15 
16 template <>
17 InputParameters
19 {
20  InputParameters params = validParams<Material>();
21 
22  params.addRequiredParam<Real>(
23  "mat_porosity", "The porosity of the material. Should be between 0 and 1. Eg, 0.1");
24  params.addCoupledVar("por_change",
25  0,
26  "An auxillary variable describing porosity changes. "
27  "Porosity = mat_porosity + por_change. If this is not "
28  "provided, zero is used.");
29  params.addRequiredParam<RealTensorValue>("mat_permeability", "The permeability tensor (m^2).");
30  params.addCoupledVar("perm_change",
31  "A list of auxillary variable describing permeability "
32  "changes. There must be 9 of these, corresponding to the "
33  "xx, xy, xz, yx, yy, yz, zx, zy, zz components respectively. "
34  " Permeability = mat_permeability*10^(perm_change).");
35  params.addRequiredParam<UserObjectName>(
36  "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
37  params.addRequiredParam<std::vector<UserObjectName>>(
38  "relperm_UO", "List of names of user objects that define relative permeability");
39  params.addRequiredParam<std::vector<UserObjectName>>(
40  "seff_UO",
41  "List of name of user objects that define effective saturation as a function of "
42  "pressure list");
43  params.addRequiredParam<std::vector<UserObjectName>>(
44  "sat_UO",
45  "List of names of user objects that define saturation as a function of effective saturation");
46  params.addRequiredParam<std::vector<UserObjectName>>(
47  "density_UO", "List of names of user objects that define the fluid density");
48  params.addRequiredParam<std::vector<UserObjectName>>(
49  "SUPG_UO", "List of names of user objects that define the SUPG");
50  params.addRequiredParam<std::vector<Real>>(
51  "viscosity", "List of viscosity of fluids (Pa.s). Typical value for water is=1E-3");
52  params.addRequiredParam<RealVectorValue>(
53  "gravity",
54  "Gravitational acceleration (m/s^2) as a vector pointing downwards. Eg (0,0,-10)");
55  // params.addRequiredCoupledVar("pressure_vars", "List of variables that represent the pressure");
56  params.addParam<bool>(
57  "linear_shape_fcns",
58  true,
59  "If you are using second-order Lagrange shape functions you need to set this to false.");
60  return params;
61 }
62 
63 RichardsMaterial::RichardsMaterial(const InputParameters & parameters)
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")),
80  _num_p(_richards_name_UO.num_v()),
81 
82  _trace_perm(_material_perm.tr()),
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")),
115  _dflux_no_mob_dgradv(
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")),
123  _d2flux_dgradvdv(
124  declareProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>("d2flux_dgradvdv")),
125  _d2flux_dvdgradv(
126  declareProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>("d2flux_dvdgradv")),
127 
128  _tauvel_SUPG(declareProperty<std::vector<RealVectorValue>>("tauvel_SUPG")),
129  _dtauvel_SUPG_dgradp(
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 }
204 
205 void
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 }
263 
264 void
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 }
364 
365 void
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 }
385 
386 void
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 }
461 
462 void
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 }
474 
475 void
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 }
568 
569 void
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 }
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
RichardsMaterial(const InputParameters &parameters)
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
const std::string density
Definition: NS.h:15
unsigned int _num_p
This holds maps between pressure_var or pressure_var, sat_var used in RichardsMaterial and kernels...
MaterialProperty< std::vector< Real > > & _sat
saturation (vector of saturations in case of multiphase)
std::vector< const RichardsSUPG * > _material_SUPG_UO
InputParameters validParams< RichardsMaterial >()
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
void zero2ndDerivedQuantities(unsigned int qp)
Zeroes 2nd derivatives of the flux.
MaterialProperty< std::vector< Real > > & _seff_old
old effective saturation
MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
d(porepressure_i)/d(variable_j)
void computePandSeff()
computes the quadpoint values of the porepressure(s) and effective saturation(s), and their derivativ...
MaterialProperty< std::vector< std::vector< std::vector< RealVectorValue > > > > & _d2flux_dvdv
d^2(Richards flux_i)/d(variable_j)/d(variable_k), here flux_i is the i_th flux, which is itself a Rea...
MaterialProperty< RealTensorValue > & _permeability
MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dflux_no_mob_dv
d(_flux_no_mob_i)/d(variable_j)
MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dtauvel_SUPG_dgradp
std::vector< const VariableGradient * > _grad_p
const VariableGradient * grad_var(unsigned int richards_var_num) const
a vector of pointers to grad(Variable)
MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dflux_dv
d(Richards flux_i)/d(variable_j), here flux_i is the i_th flux, which is itself a RealVectorValue ...
std::vector< const VariableValue * > _pressure_vals
MaterialProperty< std::vector< Real > > & _rel_perm
relative permeability (vector of relative permeabilities in case of multiphase)
MaterialProperty< std::vector< Real > > & _mass
fluid mass (a vector of masses for multicomponent)
RealVectorValue _material_gravity
gravity as entered by user
MaterialProperty< std::vector< RealVectorValue > > & _flux_no_mob
permeability*(grad(P) - density*gravity) (a vector of these for multicomponent)
std::vector< const RichardsDensity * > _material_density_UO
void computeDerivedQuantities(unsigned int qp)
Computes the "derived" quantities — those that depend on porepressure(s) and effective saturation(s)...
MaterialProperty< std::vector< Real > > & _density
fluid density (or densities for multiphase problems)
std::vector< const VariableValue * > _pressure_old_vals
MaterialProperty< std::vector< Real > > & _sat_old
old saturation
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
std::vector< const RichardsSeff * > _material_seff_UO
std::string var_types() const
return the _var_types string
MaterialProperty< std::vector< Real > > & _mass_old
old value of fluid mass (a vector of masses for multicomponent)
MaterialProperty< Real > & _porosity
void computeSUPG()
Computes the tauvel_SUPG and its derivatives.
MaterialProperty< std::vector< RealVectorValue > > & _tauvel_SUPG
MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _d2seff_dv
d^2(Seff_i)/d(variable_j)/d(variable_k)
void compute2ndDerivedQuantities(unsigned int qp)
Computes 2nd derivatives of the flux.
std::vector< const RichardsSat * > _material_sat_UO
void zeroSUPG(unsigned int qp)
Assigns and zeroes the MaterialProperties associated with SUPG.
MaterialProperty< std::vector< Real > > & _density_old
old fluid density (or densities for multiphase problems)
MaterialProperty< Real > & _porosity_old
material properties
const VariableValue * richards_vals(unsigned int richards_var_num) const
a vector of pointers to VariableValues
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
const VariableValue * richards_vals_old(unsigned int richards_var_num) const
a vector of pointers to old VariableValues
MaterialProperty< std::vector< Real > > & _viscosity
fluid viscosity (or viscosities in the multiphase case)
virtual void computeProperties()