Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "RichardsMaterial.h"
11 : #include "Assembly.h"
12 : #include "MooseMesh.h"
13 :
14 : #include <cmath> // std::sinh and std::cosh
15 :
16 : #include "libmesh/quadrature.h"
17 :
18 : registerMooseObject("RichardsApp", RichardsMaterial);
19 :
20 : InputParameters
21 1205 : RichardsMaterial::validParams()
22 : {
23 1205 : InputParameters params = Material::validParams();
24 :
25 2410 : params.addRequiredParam<Real>(
26 : "mat_porosity", "The porosity of the material. Should be between 0 and 1. Eg, 0.1");
27 2410 : params.addCoupledVar("por_change",
28 : 0,
29 : "An auxillary variable describing porosity changes. "
30 : "Porosity = mat_porosity + por_change. If this is not "
31 : "provided, zero is used.");
32 2410 : params.addRequiredParam<RealTensorValue>("mat_permeability", "The permeability tensor (m^2).");
33 2410 : params.addCoupledVar("perm_change",
34 : "A list of auxillary variable describing permeability "
35 : "changes. There must be 9 of these, corresponding to the "
36 : "xx, xy, xz, yx, yy, yz, zx, zy, zz components respectively. "
37 : " Permeability = mat_permeability*10^(perm_change).");
38 2410 : params.addRequiredParam<UserObjectName>(
39 : "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
40 2410 : params.addRequiredParam<std::vector<UserObjectName>>(
41 : "relperm_UO", "List of names of user objects that define relative permeability");
42 2410 : params.addRequiredParam<std::vector<UserObjectName>>(
43 : "seff_UO",
44 : "List of name of user objects that define effective saturation as a function of "
45 : "pressure list");
46 2410 : params.addRequiredParam<std::vector<UserObjectName>>(
47 : "sat_UO",
48 : "List of names of user objects that define saturation as a function of effective saturation");
49 2410 : params.addRequiredParam<std::vector<UserObjectName>>(
50 : "density_UO", "List of names of user objects that define the fluid density");
51 2410 : params.addRequiredParam<std::vector<UserObjectName>>(
52 : "SUPG_UO", "List of names of user objects that define the SUPG");
53 2410 : params.addRequiredParam<std::vector<Real>>(
54 : "viscosity", "List of viscosity of fluids (Pa.s). Typical value for water is=1E-3");
55 2410 : params.addRequiredParam<RealVectorValue>(
56 : "gravity",
57 : "Gravitational acceleration (m/s^2) as a vector pointing downwards. Eg (0,0,-10)");
58 : // params.addRequiredCoupledVar("pressure_vars", "List of variables that represent the pressure");
59 2410 : params.addParam<bool>(
60 : "linear_shape_fcns",
61 2410 : true,
62 : "If you are using second-order Lagrange shape functions you need to set this to false.");
63 1205 : return params;
64 0 : }
65 :
66 900 : RichardsMaterial::RichardsMaterial(const InputParameters & parameters)
67 : : Material(parameters),
68 :
69 900 : _material_por(getParam<Real>("mat_porosity")),
70 900 : _por_change(coupledValue("por_change")),
71 900 : _por_change_old(coupledValueOld("por_change")),
72 :
73 1800 : _material_perm(getParam<RealTensorValue>("mat_permeability")),
74 :
75 1800 : _material_gravity(getParam<RealVectorValue>("gravity")),
76 :
77 900 : _porosity_old(declareProperty<Real>("porosity_old")),
78 900 : _porosity(declareProperty<Real>("porosity")),
79 900 : _permeability(declareProperty<RealTensorValue>("permeability")),
80 900 : _gravity(declareProperty<RealVectorValue>("gravity")),
81 :
82 900 : _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
83 900 : _num_p(_richards_name_UO.num_v()),
84 :
85 1807 : _perm_change(isCoupled("perm_change")
86 900 : ? coupledValues("perm_change")
87 900 : : std::vector<const VariableValue *>(LIBMESH_DIM * LIBMESH_DIM, &_zero)),
88 :
89 900 : _trace_perm(_material_perm.tr()),
90 :
91 1800 : _material_viscosity(getParam<std::vector<Real>>("viscosity")),
92 :
93 900 : _pp_old(declareProperty<std::vector<Real>>("porepressure_old")),
94 900 : _pp(declareProperty<std::vector<Real>>("porepressure")),
95 900 : _dpp_dv(declareProperty<std::vector<std::vector<Real>>>("dporepressure_dv")),
96 900 : _d2pp_dv(declareProperty<std::vector<std::vector<std::vector<Real>>>>("d2porepressure_dvdv")),
97 :
98 900 : _viscosity(declareProperty<std::vector<Real>>("viscosity")),
99 :
100 900 : _density_old(declareProperty<std::vector<Real>>("density_old")),
101 900 : _density(declareProperty<std::vector<Real>>("density")),
102 900 : _ddensity_dv(declareProperty<std::vector<std::vector<Real>>>("ddensity_dv")),
103 :
104 900 : _seff_old(declareProperty<std::vector<Real>>("s_eff_old")),
105 900 : _seff(declareProperty<std::vector<Real>>("s_eff")),
106 900 : _dseff_dv(declareProperty<std::vector<std::vector<Real>>>("ds_eff_dv")),
107 900 : _d2seff_dv(declareProperty<std::vector<std::vector<std::vector<Real>>>>("d2s_eff_dvdv")),
108 :
109 900 : _sat_old(declareProperty<std::vector<Real>>("sat_old")),
110 900 : _sat(declareProperty<std::vector<Real>>("sat")),
111 900 : _dsat_dv(declareProperty<std::vector<std::vector<Real>>>("dsat_dv")),
112 :
113 900 : _rel_perm(declareProperty<std::vector<Real>>("rel_perm")),
114 900 : _drel_perm_dv(declareProperty<std::vector<std::vector<Real>>>("drel_perm_dv")),
115 :
116 900 : _mass_old(declareProperty<std::vector<Real>>("mass_old")),
117 900 : _mass(declareProperty<std::vector<Real>>("mass")),
118 900 : _dmass(declareProperty<std::vector<std::vector<Real>>>("dmass")),
119 :
120 900 : _flux_no_mob(declareProperty<std::vector<RealVectorValue>>("flux_no_mob")),
121 900 : _dflux_no_mob_dv(declareProperty<std::vector<std::vector<RealVectorValue>>>("dflux_no_mob_dv")),
122 900 : _dflux_no_mob_dgradv(
123 900 : declareProperty<std::vector<std::vector<RealTensorValue>>>("dflux_no_mob_dgradv")),
124 :
125 900 : _flux(declareProperty<std::vector<RealVectorValue>>("flux")),
126 900 : _dflux_dv(declareProperty<std::vector<std::vector<RealVectorValue>>>("dflux_dv")),
127 900 : _dflux_dgradv(declareProperty<std::vector<std::vector<RealTensorValue>>>("dflux_dgradv")),
128 900 : _d2flux_dvdv(
129 900 : declareProperty<std::vector<std::vector<std::vector<RealVectorValue>>>>("d2flux_dvdv")),
130 900 : _d2flux_dgradvdv(
131 900 : declareProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>("d2flux_dgradvdv")),
132 900 : _d2flux_dvdgradv(
133 900 : declareProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>("d2flux_dvdgradv")),
134 :
135 900 : _tauvel_SUPG(declareProperty<std::vector<RealVectorValue>>("tauvel_SUPG")),
136 900 : _dtauvel_SUPG_dgradp(
137 900 : declareProperty<std::vector<std::vector<RealTensorValue>>>("dtauvel_SUPG_dgradv")),
138 1800 : _dtauvel_SUPG_dp(declareProperty<std::vector<std::vector<RealVectorValue>>>("dtauvel_SUPG_dv"))
139 :
140 : {
141 :
142 : // Need to add the variables that the user object is coupled to as dependencies so MOOSE will
143 : // compute them
144 : {
145 : const std::vector<MooseVariableFEBase *> & coupled_vars =
146 900 : _richards_name_UO.getCoupledMooseVars();
147 2085 : for (unsigned int i = 0; i < coupled_vars.size(); i++)
148 2370 : addMooseVariableDependency(coupled_vars[i]);
149 : }
150 :
151 900 : if (_material_por <= 0 || _material_por >= 1)
152 2 : mooseError("Porosity set to ", _material_por, " but it must be between 0 and 1");
153 :
154 911 : if (isCoupled("perm_change") && (coupledComponents("perm_change") != LIBMESH_DIM * LIBMESH_DIM))
155 0 : mooseError(LIBMESH_DIM * LIBMESH_DIM,
156 : " components of perm_change must be given to a RichardsMaterial. You supplied ",
157 1 : coupledComponents("perm_change"),
158 : "\n");
159 :
160 897 : if (!(_material_viscosity.size() == _num_p &&
161 2689 : getParam<std::vector<UserObjectName>>("relperm_UO").size() == _num_p &&
162 2687 : getParam<std::vector<UserObjectName>>("seff_UO").size() == _num_p &&
163 2685 : getParam<std::vector<UserObjectName>>("sat_UO").size() == _num_p &&
164 2683 : getParam<std::vector<UserObjectName>>("density_UO").size() == _num_p &&
165 2682 : getParam<std::vector<UserObjectName>>("SUPG_UO").size() == _num_p))
166 6 : mooseError("There are ",
167 6 : _num_p,
168 : " Richards fluid variables, so you need to specify this "
169 : "number of viscosities, relperm_UO, seff_UO, sat_UO, "
170 : "density_UO, SUPG_UO");
171 :
172 891 : _d2density.resize(_num_p);
173 891 : _d2rel_perm_dv.resize(_num_p);
174 891 : _pressure_vals.resize(_num_p);
175 891 : _pressure_old_vals.resize(_num_p);
176 891 : _material_relperm_UO.resize(_num_p);
177 891 : _material_seff_UO.resize(_num_p);
178 891 : _material_sat_UO.resize(_num_p);
179 891 : _material_density_UO.resize(_num_p);
180 891 : _material_SUPG_UO.resize(_num_p);
181 891 : _grad_p.resize(_num_p);
182 :
183 2061 : for (unsigned int i = 0; i < _num_p; ++i)
184 : {
185 : // DON'T WANT "pressure_vars" at all since pp_name_UO contains the same info
186 : //_pressure_vals[i] = &coupledValue("pressure_vars", i); // coupled value returns a reference
187 : //_pressure_old_vals[i] = (_is_transient ? &coupledValueOld("pressure_vars", i) : &_zero);
188 : //_grad_p[i] = &coupledGradient("pressure_vars", i);
189 :
190 : // in the following. first get the userobject names that were inputted, then get the i_th one
191 : // of these, then get the actual userobject that this corresponds to, then finally & gives
192 : // pointer to RichardsRelPerm object.
193 1170 : _material_relperm_UO[i] = &getUserObjectByName<RichardsRelPerm>(
194 2340 : getParam<std::vector<UserObjectName>>("relperm_UO")[i]);
195 1170 : _material_seff_UO[i] =
196 2340 : &getUserObjectByName<RichardsSeff>(getParam<std::vector<UserObjectName>>("seff_UO")[i]);
197 1170 : _material_sat_UO[i] =
198 2340 : &getUserObjectByName<RichardsSat>(getParam<std::vector<UserObjectName>>("sat_UO")[i]);
199 1170 : _material_density_UO[i] = &getUserObjectByName<RichardsDensity>(
200 1170 : getParam<std::vector<UserObjectName>>("density_UO")[i]);
201 1170 : _material_SUPG_UO[i] =
202 3510 : &getUserObjectByName<RichardsSUPG>(getParam<std::vector<UserObjectName>>("SUPG_UO")[i]);
203 : }
204 891 : }
205 :
206 : void
207 1899373 : RichardsMaterial::computePandSeff()
208 : {
209 : // Get the pressure and effective saturation at each quadpoint
210 : // From these we will build the relative permeability, density, flux, etc
211 1899373 : if (_richards_name_UO.var_types() == "pppp")
212 : {
213 4108964 : for (unsigned int i = 0; i < _num_p; ++i)
214 : {
215 2209591 : _pressure_vals[i] = _richards_name_UO.richards_vals(i);
216 2209591 : _pressure_old_vals[i] = _richards_name_UO.richards_vals_old(i);
217 2209591 : _grad_p[i] = _richards_name_UO.grad_var(i);
218 : }
219 : }
220 :
221 7784090 : for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
222 : {
223 5884717 : _pp_old[qp].resize(_num_p);
224 5884717 : _pp[qp].resize(_num_p);
225 5884717 : _dpp_dv[qp].resize(_num_p);
226 5884717 : _d2pp_dv[qp].resize(_num_p);
227 :
228 5884717 : _seff_old[qp].resize(_num_p);
229 5884717 : _seff[qp].resize(_num_p);
230 5884717 : _dseff_dv[qp].resize(_num_p);
231 5884717 : _d2seff_dv[qp].resize(_num_p);
232 :
233 5884717 : if (_richards_name_UO.var_types() == "pppp")
234 : {
235 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
236 : {
237 6753142 : _pp_old[qp][i] = (*_pressure_old_vals[i])[qp];
238 6753142 : _pp[qp][i] = (*_pressure_vals[i])[qp];
239 :
240 6753142 : _dpp_dv[qp][i].assign(_num_p, 0);
241 6753142 : _dpp_dv[qp][i][i] = 1;
242 :
243 6753142 : _d2pp_dv[qp][i].resize(_num_p);
244 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
245 8489992 : _d2pp_dv[qp][i][j].assign(_num_p, 0);
246 :
247 6753142 : _seff_old[qp][i] = (*_material_seff_UO[i]).seff(_pressure_old_vals, qp);
248 6753142 : _seff[qp][i] = (*_material_seff_UO[i]).seff(_pressure_vals, qp);
249 :
250 6753142 : _dseff_dv[qp][i].resize(_num_p);
251 6753142 : (*_material_seff_UO[i]).dseff(_pressure_vals, qp, _dseff_dv[qp][i]);
252 :
253 6753142 : _d2seff_dv[qp][i].resize(_num_p);
254 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
255 8489992 : _d2seff_dv[qp][i][j].resize(_num_p);
256 6753142 : (*_material_seff_UO[i]).d2seff(_pressure_vals, qp, _d2seff_dv[qp][i]);
257 : }
258 : }
259 : // the above lines of code are only valid for "pppp"
260 : // if you decide to code other RichardsVariables (eg "psss")
261 : // you will need to add some lines here
262 : }
263 1899373 : }
264 :
265 : void
266 5884717 : RichardsMaterial::computeDerivedQuantities(unsigned int qp)
267 : {
268 : // fluid viscosity
269 5884717 : _viscosity[qp].resize(_num_p);
270 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
271 6753142 : _viscosity[qp][i] = _material_viscosity[i];
272 :
273 : // fluid saturation
274 5884717 : _sat_old[qp].resize(_num_p);
275 5884717 : _sat[qp].resize(_num_p);
276 5884717 : _dsat_dv[qp].resize(_num_p);
277 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
278 : {
279 6753142 : _sat_old[qp][i] = (*_material_sat_UO[i]).sat(_seff_old[qp][i]);
280 6753142 : _sat[qp][i] = (*_material_sat_UO[i]).sat(_seff[qp][i]);
281 6753142 : _dsat_dv[qp][i].assign(_num_p, (*_material_sat_UO[i]).dsat(_seff[qp][i]));
282 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
283 8489992 : _dsat_dv[qp][i][j] *= _dseff_dv[qp][i][j];
284 : }
285 :
286 : // fluid density
287 5884717 : _density_old[qp].resize(_num_p);
288 5884717 : _density[qp].resize(_num_p);
289 5884717 : _ddensity_dv[qp].resize(_num_p);
290 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
291 : {
292 6753142 : _density_old[qp][i] = (*_material_density_UO[i]).density(_pp_old[qp][i]);
293 6753142 : _density[qp][i] = (*_material_density_UO[i]).density(_pp[qp][i]);
294 6753142 : _ddensity_dv[qp][i].assign(_num_p, (*_material_density_UO[i]).ddensity(_pp[qp][i]));
295 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
296 8489992 : _ddensity_dv[qp][i][j] *= _dpp_dv[qp][i][j];
297 : }
298 :
299 : // relative permeability
300 5884717 : _rel_perm[qp].resize(_num_p);
301 5884717 : _drel_perm_dv[qp].resize(_num_p);
302 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
303 : {
304 6753142 : _rel_perm[qp][i] = (*_material_relperm_UO[i]).relperm(_seff[qp][i]);
305 6753142 : _drel_perm_dv[qp][i].assign(_num_p, (*_material_relperm_UO[i]).drelperm(_seff[qp][i]));
306 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
307 8489992 : _drel_perm_dv[qp][i][j] *= _dseff_dv[qp][i][j];
308 : }
309 :
310 : // fluid mass
311 5884717 : _mass_old[qp].resize(_num_p);
312 5884717 : _mass[qp].resize(_num_p);
313 5884717 : _dmass[qp].resize(_num_p);
314 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
315 : {
316 6753142 : _mass_old[qp][i] = _porosity_old[qp] * _density_old[qp][i] * _sat_old[qp][i];
317 6753142 : _mass[qp][i] = _porosity[qp] * _density[qp][i] * _sat[qp][i];
318 6753142 : _dmass[qp][i].resize(_num_p);
319 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
320 8489992 : _dmass[qp][i][j] = _porosity[qp] * (_ddensity_dv[qp][i][j] * _sat[qp][i] +
321 8489992 : _density[qp][i] * _dsat_dv[qp][i][j]);
322 : }
323 :
324 : // flux without the mobility part
325 5884717 : _flux_no_mob[qp].resize(_num_p);
326 5884717 : _dflux_no_mob_dv[qp].resize(_num_p);
327 5884717 : _dflux_no_mob_dgradv[qp].resize(_num_p);
328 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
329 : {
330 6753142 : _flux_no_mob[qp][i] = _permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]);
331 :
332 6753142 : _dflux_no_mob_dv[qp][i].resize(_num_p);
333 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
334 8489992 : _dflux_no_mob_dv[qp][i][j] = _permeability[qp] * (-_ddensity_dv[qp][i][j] * _gravity[qp]);
335 :
336 6753142 : _dflux_no_mob_dgradv[qp][i].resize(_num_p);
337 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
338 8489992 : _dflux_no_mob_dgradv[qp][i][j] = _permeability[qp] * _dpp_dv[qp][i][j];
339 : }
340 :
341 : // flux
342 5884717 : _flux[qp].resize(_num_p);
343 5884717 : _dflux_dv[qp].resize(_num_p);
344 5884717 : _dflux_dgradv[qp].resize(_num_p);
345 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
346 : {
347 6753142 : _flux[qp][i] = _density[qp][i] * _rel_perm[qp][i] * _flux_no_mob[qp][i] / _viscosity[qp][i];
348 :
349 6753142 : _dflux_dv[qp][i].resize(_num_p);
350 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
351 : {
352 8489992 : _dflux_dv[qp][i][j] =
353 8489992 : _density[qp][i] * _rel_perm[qp][i] * _dflux_no_mob_dv[qp][i][j] / _viscosity[qp][i];
354 : _dflux_dv[qp][i][j] +=
355 8489992 : (_ddensity_dv[qp][i][j] * _rel_perm[qp][i] + _density[qp][i] * _drel_perm_dv[qp][i][j]) *
356 8489992 : _flux_no_mob[qp][i] / _viscosity[qp][i];
357 : }
358 :
359 6753142 : _dflux_dgradv[qp][i].resize(_num_p);
360 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
361 8489992 : _dflux_dgradv[qp][i][j] =
362 8489992 : _density[qp][i] * _rel_perm[qp][i] * _dflux_no_mob_dgradv[qp][i][j] / _viscosity[qp][i];
363 : }
364 5884717 : }
365 :
366 : void
367 5884717 : RichardsMaterial::zero2ndDerivedQuantities(unsigned int qp)
368 : {
369 5884717 : _d2flux_dvdv[qp].resize(_num_p);
370 5884717 : _d2flux_dgradvdv[qp].resize(_num_p);
371 5884717 : _d2flux_dvdgradv[qp].resize(_num_p);
372 :
373 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
374 : {
375 6753142 : _d2flux_dvdv[qp][i].resize(_num_p);
376 6753142 : _d2flux_dgradvdv[qp][i].resize(_num_p);
377 6753142 : _d2flux_dvdgradv[qp][i].resize(_num_p);
378 15243134 : for (unsigned int j = 0; j < _num_p; ++j)
379 : {
380 8489992 : _d2flux_dvdv[qp][i][j].assign(_num_p, RealVectorValue());
381 8489992 : _d2flux_dgradvdv[qp][i][j].assign(_num_p, RealTensorValue());
382 8489992 : _d2flux_dvdgradv[qp][i][j].assign(_num_p, RealTensorValue());
383 : }
384 : }
385 5884717 : }
386 :
387 : void
388 5884717 : RichardsMaterial::compute2ndDerivedQuantities(unsigned int qp)
389 : {
390 5884717 : zero2ndDerivedQuantities(qp);
391 :
392 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
393 : {
394 6753142 : if ((*_material_SUPG_UO[i]).SUPG_trivial())
395 475124 : continue; // as the derivatives won't be needed
396 :
397 : // second derivative of density
398 6278018 : _d2density[i].resize(_num_p);
399 6278018 : Real ddens = (*_material_density_UO[i]).ddensity(_pp[qp][i]);
400 6278018 : Real d2dens = (*_material_density_UO[i]).d2density(_pp[qp][i]);
401 14199510 : for (unsigned int j = 0; j < _num_p; ++j)
402 : {
403 7921492 : _d2density[i][j].resize(_num_p);
404 19129932 : for (unsigned int k = 0; k < _num_p; ++k)
405 11208440 : _d2density[i][j][k] =
406 11208440 : d2dens * _dpp_dv[qp][i][j] * _dpp_dv[qp][i][k] + ddens * _d2pp_dv[qp][i][j][k];
407 : }
408 :
409 : // second derivative of relative permeability
410 6278018 : _d2rel_perm_dv[i].resize(_num_p);
411 6278018 : Real drel = (*_material_relperm_UO[i]).drelperm(_seff[qp][i]);
412 6278018 : Real d2rel = (*_material_relperm_UO[i]).d2relperm(_seff[qp][i]);
413 14199510 : for (unsigned int j = 0; j < _num_p; ++j)
414 : {
415 7921492 : _d2rel_perm_dv[i][j].resize(_num_p);
416 19129932 : for (unsigned int k = 0; k < _num_p; ++k)
417 11208440 : _d2rel_perm_dv[i][j][k] =
418 11208440 : d2rel * _dseff_dv[qp][i][j] * _dseff_dv[qp][i][k] + drel * _d2seff_dv[qp][i][j][k];
419 : }
420 :
421 : // now compute the second derivs of the fluxes
422 14199510 : for (unsigned int j = 0; j < _num_p; ++j)
423 : {
424 19129932 : for (unsigned int k = 0; k < _num_p; ++k)
425 : {
426 11208440 : _d2flux_dvdv[qp][i][j][k] =
427 11208440 : _d2density[i][j][k] * _rel_perm[qp][i] *
428 11208440 : (_permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]));
429 : _d2flux_dvdv[qp][i][j][k] +=
430 11208440 : (_ddensity_dv[qp][i][j] * _drel_perm_dv[qp][i][k] +
431 11208440 : _ddensity_dv[qp][i][k] * _drel_perm_dv[qp][i][j]) *
432 11208440 : (_permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]));
433 : _d2flux_dvdv[qp][i][j][k] +=
434 11208440 : _density[qp][i] * _d2rel_perm_dv[i][j][k] *
435 11208440 : (_permeability[qp] * ((*_grad_p[i])[qp] - _density[qp][i] * _gravity[qp]));
436 11208440 : _d2flux_dvdv[qp][i][j][k] += (_ddensity_dv[qp][i][j] * _rel_perm[qp][i] +
437 11208440 : _density[qp][i] * _drel_perm_dv[qp][i][j]) *
438 11208440 : (_permeability[qp] * (-_ddensity_dv[qp][i][k] * _gravity[qp]));
439 11208440 : _d2flux_dvdv[qp][i][j][k] += (_ddensity_dv[qp][i][k] * _rel_perm[qp][i] +
440 11208440 : _density[qp][i] * _drel_perm_dv[qp][i][k]) *
441 11208440 : (_permeability[qp] * (-_ddensity_dv[qp][i][j] * _gravity[qp]));
442 11208440 : _d2flux_dvdv[qp][i][j][k] += _density[qp][i] * _rel_perm[qp][i] *
443 11208440 : (_permeability[qp] * (-_d2density[i][j][k] * _gravity[qp]));
444 : }
445 : }
446 14199510 : for (unsigned int j = 0; j < _num_p; ++j)
447 19129932 : for (unsigned int k = 0; k < _num_p; ++k)
448 11208440 : _d2flux_dvdv[qp][i][j][k] /= _viscosity[qp][i];
449 :
450 14199510 : for (unsigned int j = 0; j < _num_p; ++j)
451 : {
452 19129932 : for (unsigned int k = 0; k < _num_p; ++k)
453 : {
454 11208440 : _d2flux_dgradvdv[qp][i][j][k] = (_ddensity_dv[qp][i][k] * _rel_perm[qp][i] +
455 11208440 : _density[qp][i] * _drel_perm_dv[qp][i][k]) *
456 11208440 : _permeability[qp] * _dpp_dv[qp][i][j] / _viscosity[qp][i];
457 11208440 : _d2flux_dvdgradv[qp][i][k][j] = _d2flux_dgradvdv[qp][i][j][k];
458 : }
459 : }
460 : }
461 5884717 : }
462 :
463 : void
464 5884717 : RichardsMaterial::zeroSUPG(unsigned int qp)
465 : {
466 5884717 : _tauvel_SUPG[qp].assign(_num_p, RealVectorValue());
467 5884717 : _dtauvel_SUPG_dgradp[qp].resize(_num_p);
468 5884717 : _dtauvel_SUPG_dp[qp].resize(_num_p);
469 12637859 : for (unsigned int i = 0; i < _num_p; ++i)
470 : {
471 6753142 : _dtauvel_SUPG_dp[qp][i].assign(_num_p, RealVectorValue());
472 6753142 : _dtauvel_SUPG_dgradp[qp][i].assign(_num_p, RealTensorValue());
473 : }
474 5884717 : }
475 :
476 : void
477 1756151 : RichardsMaterial::computeSUPG()
478 : {
479 : // Grab reference to linear Lagrange finite element object pointer,
480 : // currently this is always a linear Lagrange element, so this might need to
481 : // be generalized if we start working with higher-order elements...
482 5268453 : auto && fe = _assembly.getFE(getParam<bool>("linear_shape_fcns") ? FEType(FIRST, LAGRANGE)
483 : : FEType(SECOND, LAGRANGE),
484 1756151 : _current_elem->dim());
485 :
486 : // Grab references to FE object's mapping data from the _subproblem's FE object
487 1756151 : const std::vector<Real> & dxidx(fe->get_dxidx());
488 : const std::vector<Real> & dxidy(fe->get_dxidy());
489 : const std::vector<Real> & dxidz(fe->get_dxidz());
490 : const std::vector<Real> & detadx(fe->get_detadx());
491 : const std::vector<Real> & detady(fe->get_detady());
492 : const std::vector<Real> & detadz(fe->get_detadz());
493 : const std::vector<Real> & dzetadx(fe->get_dzetadx());
494 : const std::vector<Real> & dzetady(fe->get_dzetady());
495 : const std::vector<Real> & dzetadz(fe->get_dzetadz());
496 :
497 7212432 : for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
498 : {
499 :
500 : // Bounds checking on element data and putting into vector form
501 : mooseAssert(qp < dxidx.size(), "Insufficient data in dxidx array!");
502 : mooseAssert(qp < dxidy.size(), "Insufficient data in dxidy array!");
503 : mooseAssert(qp < dxidz.size(), "Insufficient data in dxidz array!");
504 5456281 : if (_mesh.dimension() >= 2)
505 : {
506 : mooseAssert(qp < detadx.size(), "Insufficient data in detadx array!");
507 : mooseAssert(qp < detady.size(), "Insufficient data in detady array!");
508 : mooseAssert(qp < detadz.size(), "Insufficient data in detadz array!");
509 : }
510 5456281 : if (_mesh.dimension() >= 3)
511 : {
512 : mooseAssert(qp < dzetadx.size(), "Insufficient data in dzetadx array!");
513 : mooseAssert(qp < dzetady.size(), "Insufficient data in dzetady array!");
514 : mooseAssert(qp < dzetadz.size(), "Insufficient data in dzetadz array!");
515 : }
516 :
517 : // CHECK : Does this work spherical, cylindrical, etc?
518 5456281 : RealVectorValue xi_prime(dxidx[qp], dxidy[qp], dxidz[qp]);
519 : RealVectorValue eta_prime, zeta_prime;
520 5456281 : if (_mesh.dimension() >= 2)
521 : {
522 3122603 : eta_prime(0) = detadx[qp];
523 3122603 : eta_prime(1) = detady[qp];
524 : }
525 5456281 : if (_mesh.dimension() == 3)
526 : {
527 1725471 : eta_prime(2) = detadz[qp];
528 1725471 : zeta_prime(0) = dzetadx[qp];
529 1725471 : zeta_prime(1) = dzetady[qp];
530 1725471 : zeta_prime(2) = dzetadz[qp];
531 : }
532 :
533 5456281 : _trace_perm = _permeability[qp].tr();
534 11734299 : for (unsigned int i = 0; i < _num_p; ++i)
535 : {
536 : RealVectorValue vel =
537 6278018 : (*_material_SUPG_UO[i])
538 6278018 : .velSUPG(_permeability[qp], (*_grad_p[i])[qp], _density[qp][i], _gravity[qp]);
539 6278018 : RealTensorValue dvel_dgradp = (*_material_SUPG_UO[i]).dvelSUPG_dgradp(_permeability[qp]);
540 : RealVectorValue dvel_dp =
541 : (*_material_SUPG_UO[i])
542 6278018 : .dvelSUPG_dp(_permeability[qp], _ddensity_dv[qp][i][i], _gravity[qp]);
543 : RealVectorValue bb =
544 6278018 : (*_material_SUPG_UO[i]).bb(vel, _mesh.dimension(), xi_prime, eta_prime, zeta_prime);
545 : RealVectorValue dbb2_dgradp =
546 6278018 : (*_material_SUPG_UO[i]).dbb2_dgradp(vel, dvel_dgradp, xi_prime, eta_prime, zeta_prime);
547 6278018 : Real dbb2_dp = (*_material_SUPG_UO[i]).dbb2_dp(vel, dvel_dp, xi_prime, eta_prime, zeta_prime);
548 6278018 : Real tau = (*_material_SUPG_UO[i]).tauSUPG(vel, _trace_perm, bb);
549 : RealVectorValue dtau_dgradp =
550 6278018 : (*_material_SUPG_UO[i]).dtauSUPG_dgradp(vel, dvel_dgradp, _trace_perm, bb, dbb2_dgradp);
551 6278018 : Real dtau_dp = (*_material_SUPG_UO[i]).dtauSUPG_dp(vel, dvel_dp, _trace_perm, bb, dbb2_dp);
552 :
553 6278018 : _tauvel_SUPG[qp][i] = tau * vel;
554 :
555 6278018 : RealTensorValue dtauvel_dgradp = tau * dvel_dgradp;
556 25112072 : for (const auto j : make_range(Moose::dim))
557 75336216 : for (const auto k : make_range(Moose::dim))
558 56502162 : dtauvel_dgradp(j, k) +=
559 56502162 : dtau_dgradp(j) * vel(k); // this is outerproduct - maybe libmesh can do it better?
560 14199510 : for (unsigned int j = 0; j < _num_p; ++j)
561 7921492 : _dtauvel_SUPG_dgradp[qp][i][j] = dtauvel_dgradp * _dpp_dv[qp][i][j];
562 :
563 : RealVectorValue dtauvel_dp = dtau_dp * vel + tau * dvel_dp;
564 14199510 : for (unsigned int j = 0; j < _num_p; ++j)
565 7921492 : _dtauvel_SUPG_dp[qp][i][j] = dtauvel_dp * _dpp_dv[qp][i][j];
566 : }
567 : }
568 1756151 : }
569 :
570 : void
571 1899373 : RichardsMaterial::computeProperties()
572 : {
573 : // compute porepressures and effective saturations
574 : // with algorithms depending on the _richards_name_UO.var_types()
575 1899373 : computePandSeff();
576 :
577 : // porosity, permeability, and gravity
578 7784090 : for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
579 : {
580 5884717 : _porosity[qp] = _material_por + _por_change[qp];
581 5884717 : _porosity_old[qp] = _material_por + _por_change_old[qp];
582 :
583 5884717 : _permeability[qp] = _material_perm;
584 23538868 : for (const auto i : make_range(Moose::dim))
585 70616604 : for (const auto j : make_range(Moose::dim))
586 52962453 : _permeability[qp](i, j) *= std::pow(10, (*_perm_change[LIBMESH_DIM * i + j])[qp]);
587 :
588 5884717 : _gravity[qp] = _material_gravity;
589 : }
590 :
591 : // compute "derived" quantities -- those that depend on P and Seff --- such as density, relperm
592 7784090 : for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
593 5884717 : computeDerivedQuantities(qp);
594 :
595 : // compute certain second derivatives of the derived quantities
596 : // These are needed in Jacobian calculations if doing SUPG
597 7784090 : for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
598 5884717 : compute2ndDerivedQuantities(qp);
599 :
600 : // Now for SUPG itself
601 7784090 : for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
602 5884717 : zeroSUPG(qp);
603 :
604 : // the following saves computational effort if all SUPG is trivial
605 : bool trivial_supg = true;
606 4108964 : for (unsigned int i = 0; i < _num_p; ++i)
607 2209591 : trivial_supg = trivial_supg && (*_material_SUPG_UO[i]).SUPG_trivial();
608 1899373 : if (trivial_supg)
609 : return;
610 : else
611 1756151 : computeSUPG();
612 : }
|
