v0.14.0
Loading...
Searching...
No Matches
MaterialUnsaturatedFlow.hpp
Go to the documentation of this file.
1/** \file MaterialUnsaturatedFlow.hpp
2 * \brief Mix implementation of transport element
3 *
4 * \ingroup mofem_mix_transport_elem
5 *
6 */
7
8#ifndef __MATERIALUNSATURATEDFLOW_HPP__
9#define __MATERIALUNSATURATEDFLOW_HPP__
10
11namespace MixTransport {
12
14
16 blockId = 0;
17 matName = "SimpleDarcy";
18 sCale = 1;
19 thetaS = 0.38;
20 thetaM = 0.38;
21 thetaR = 0.068;
22 alpha = 0.8;
23 n = 1.09;
24 Ks = 0.048;
25 hS = 0;
26 Ah = 0.;
27 AhZ = 0.;
28 AhZZ = 0;
29 }
30
31 int blockId; ///< Block Id
32 std::string matName; ///< material name
33
34 double Ks; ///< Saturated hydraulic conductivity [m/day]
35 double hS; ///< minimum capillary height [m]
36 double thetaS; ///< saturated water content
37 double thetaR; ///< residual water contents
38 double thetaM; ///< model parameter
39 double alpha; ///< model parameter
40 double n; ///< model parameter
41
42 // Initial hydraulic head
43 double Ah; ///< Initial hydraulic head coefficient
44 double AhZ; ///< Initial hydraulic head coefficient
45 double AhZZ; ///< Initial hydraulic head coefficient
46
47 double initialPcEval() const { return Ah + AhZ * z + AhZZ * z * z; }
48
49 void addOptions(po::options_description &o, const std::string &prefix) {
50 o.add_options()((prefix + ".material_name").c_str(),
51 po::value<std::string>(&matName)->default_value(matName))(
52 (prefix + ".ePsilon0").c_str(),
53 po::value<double>(&ePsilon0)->default_value(ePsilon0))(
54 (prefix + ".ePsilon1").c_str(),
55 po::value<double>(&ePsilon1)->default_value(ePsilon1))(
56 (prefix + ".sCale").c_str(),
57 po::value<double>(&sCale)->default_value(sCale))(
58 (prefix + ".thetaS").c_str(),
59 po::value<double>(&thetaS)->default_value(thetaS))(
60 (prefix + ".thetaR").c_str(),
61 po::value<double>(&thetaR)->default_value(thetaR))(
62 (prefix + ".thetaM").c_str(),
63 po::value<double>(&thetaM)->default_value(thetaM))(
64 (prefix + ".alpha").c_str(),
65 po::value<double>(&alpha)->default_value(alpha))(
66 (prefix + ".n").c_str(), po::value<double>(&n)->default_value(n))(
67 (prefix + ".hS").c_str(), po::value<double>(&hS)->default_value(hS))(
68 (prefix + ".Ks").c_str(), po::value<double>(&Ks)->default_value(Ks))(
69 (prefix + ".Ah").c_str(), po::value<double>(&Ah)->default_value(Ah))(
70 (prefix + ".AhZ").c_str(), po::value<double>(&AhZ)->default_value(AhZ))(
71 (prefix + ".AhZZ").c_str(),
72 po::value<double>(&AhZZ)->default_value(AhZZ))(
73 (prefix + ".scaleZ").c_str(),
74 po::value<double>(&scaleZ)->default_value(scaleZ));
75 }
76
77 void printMatParameters(const int id, const std::string &prefix) const {
78 PetscPrintf(PETSC_COMM_WORLD, "Mat name %s-%s block id %d\n",
79 prefix.c_str(), matName.c_str(), id);
80 PetscPrintf(PETSC_COMM_WORLD, "Material name: %s\n", matName.c_str());
81 PetscPrintf(PETSC_COMM_WORLD, "thetaS=%6.4g\n", thetaS);
82 PetscPrintf(PETSC_COMM_WORLD, "thetaR=%6.4g\n", thetaR);
83 PetscPrintf(PETSC_COMM_WORLD, "thetaM=%6.4g\n", thetaM);
84 PetscPrintf(PETSC_COMM_WORLD, "alpha=%6.4g\n", alpha);
85 PetscPrintf(PETSC_COMM_WORLD, "n=%6.4g\n", n);
86 PetscPrintf(PETSC_COMM_WORLD, "hS=%6.4g\n", hS);
87 PetscPrintf(PETSC_COMM_WORLD, "Ks=%6.4g\n", Ks);
88 PetscPrintf(PETSC_COMM_WORLD, "Ah=%6.4g\n", Ah);
89 PetscPrintf(PETSC_COMM_WORLD, "AhZ=%6.4g\n", AhZ);
90 PetscPrintf(PETSC_COMM_WORLD, "AhZZ=%6.4g\n", AhZZ);
91 PetscPrintf(PETSC_COMM_WORLD, "ePsilon0=%6.4g\n", ePsilon0);
92 PetscPrintf(PETSC_COMM_WORLD, "ePsilon1=%6.4g\n", ePsilon1);
93 PetscPrintf(PETSC_COMM_WORLD, "sCale=%6.4g\n", sCale);
94 PetscPrintf(PETSC_COMM_WORLD, "scaleZ=%6.4g\n", scaleZ);
95 }
96
97 typedef boost::function<boost::shared_ptr<CommonMaterialData>(
98 const CommonMaterialData &data)>
100};
101
103
104 static boost::shared_ptr<CommonMaterialData>
106 return boost::shared_ptr<CommonMaterialData>(new MaterialDarcy(data));
107 }
108
110
111 MoFEMErrorCode calK() {
113 K = Ks;
115 };
116
117 MoFEMErrorCode calDiffK() {
119 diffK = 0;
121 };
122
123 MoFEMErrorCode calC() {
125 C = ePsilon1;
127 }
128
129 MoFEMErrorCode calDiffC() {
131 diffC = 0;
133 }
134
135 MoFEMErrorCode calTheta() {
137 tHeta = thetaS;
139 }
140
141 virtual MoFEMErrorCode calSe() {
143 Se = 1;
145 }
146};
147
150 : CommonMaterialData(data) {}
151
152 template <typename TYPE> inline TYPE funSe(TYPE &theta) {
153 return (theta - thetaR) / (thetaS - thetaR);
154 }
155
156 virtual void recordTheta() = 0;
157
158 virtual void recordKr() = 0;
159
160 double Kr;
161 MoFEMErrorCode calK() {
163 if (h < hS) {
164 int r = ::function(2 * blockId + 1, 1, 1, &h, &Kr);
165 if (r < 0) {
166 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
167 "ADOL-C function evaluation with error");
168 }
169 K = Ks * Kr;
170 } else {
171 K = Ks;
172 }
173 K += Ks * ePsilon0;
175 };
176
177 double diffKr;
178 MoFEMErrorCode calDiffK() {
180 if (h < hS) {
181 diffK = 0;
182 int r = ::gradient(2 * blockId + 1, 1, &h, &diffKr);
183 if (r < 0) {
184 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
185 "ADOL-C function evaluation with error");
186 }
187 diffK = Ks * diffKr;
188 } else {
189 diffK = 0;
190 }
192 };
193
194 MoFEMErrorCode calC() {
196 if (h < hS) {
197 int r = ::gradient(2 * blockId + 0, 1, &h, &C);
198 if (r < 0) {
199 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
200 "ADOL-C function evaluation with error");
201 }
202 } else {
203 C = 0;
204 }
205 C += ePsilon1;
207 }
208
209 MoFEMErrorCode calDiffC() {
211 if (h < hS) {
212 double v = 1;
213 int r = ::hess_vec(2 * blockId + 0, 1, &h, &v, &diffC);
214 if (r < 0) {
215 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
216 "ADOL-C function evaluation with error");
217 }
218 } else {
219 diffC = 0;
220 }
222 }
223
224 MoFEMErrorCode calTheta() {
226 if (h < hS) {
227 int r = ::function(2 * blockId + 0, 1, 1, &h, &tHeta);
228 if (r < 0) {
229 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
230 "ADOL-C function evaluation with error");
231 }
232 } else {
233 tHeta = thetaS;
234 }
236 }
237
238 MoFEMErrorCode calSe() {
240 if (h < hS) {
241 int r = ::function(2 * blockId + 0, 1, 1, &h, &tHeta);
242 Se = funSe(tHeta);
243 if (r < 0) {
244 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
245 "ADOL-C function evaluation with error");
246 }
247 } else {
248 tHeta = thetaS;
249 }
251 }
252};
253
255
256 static boost::shared_ptr<CommonMaterialData>
258 return boost::shared_ptr<CommonMaterialData>(
259 new MaterialVanGenuchten(data));
260 }
261
264 recordTheta();
265 recordKr();
266 }
267
273
274 template <typename TYPE> inline TYPE funTheta(TYPE &h, const double m) {
275 return thetaR + (thetaM - thetaR) / pow(1 + pow(-alpha * h, n), m);
276 }
277
278 template <typename TYPE>
279 inline TYPE funFunSeStar(TYPE &SeStar, const double m) {
280 return pow(1 - pow(SeStar, 1 / m), m);
281 }
282
284 const double m = 1 - 1 / n;
285 aTheta = funTheta(ah, m);
286 aSe = funSe(aTheta);
287 aSeStar = aSe * (thetaS - thetaR) / (thetaM - thetaR);
288 double one = 1;
289 const double c = funFunSeStar<double>(one, m);
290 return sqrt(aSe) *
291 pow((1 - funFunSeStar<adouble>(aSeStar, m)) / (1 - c), 2);
292 }
293
294 virtual void recordTheta() {
295 h = -1 - hS;
296 trace_on(2 * blockId + 0, true);
297 ah <<= h;
298 const double m = 1 - 1 / n;
299 aTheta = funTheta(ah, m);
300 double r_theta;
301 aTheta >>= r_theta;
302 trace_off();
303 }
304
305 virtual void recordKr() {
306 h = -1 - hS;
307 trace_on(2 * blockId + 1, true);
308 ah <<= h;
309 aKr = funKr(ah);
310 double r_Kr;
311 aKr >>= r_Kr;
312 trace_off();
313 }
314
315 void printTheta(const double b, const double e, double s,
316 const std::string &prefix) {
317 const double m = 1 - 1 / n;
318 h = b;
319 for (; h >= e; h += s) {
320 s = -pow(-s, 0.9);
321 double theta = funTheta(h, m);
322 double Se = (theta - thetaR) / (thetaS - thetaR);
323 PetscPrintf(PETSC_COMM_SELF, "%s %6.4e %6.4e %6.4e\n", prefix.c_str(), h,
324 theta, Se);
325 }
326 }
327
328 void printKappa(const double b, const double e, double s,
329 const std::string &prefix) {
330 h = b;
331 for (; h >= e; h += s) {
332 s = -pow(-s, 0.9);
333 calK();
334 PetscPrintf(PETSC_COMM_SELF, "%s %6.4e %6.4e %6.4e\n", prefix.c_str(), h,
335 Kr, K);
336 }
337 }
338
339 void printC(const double b, const double e, double s,
340 const std::string &prefix) {
341 h = b;
342 for (; h >= e; h += s) {
343 s = -pow(-s, 0.9);
344 calC();
345 PetscPrintf(PETSC_COMM_SELF, "%s %6.4e %6.4e\n", prefix.c_str(), h, C);
346 }
347 }
348};
349
351 static map<std::string, CommonMaterialData::RegisterHook>
353 MoFEMErrorCode operator()() const {
356 mapOfRegistredMaterials["VanGenuchten"] =
359 }
360};
361
362} // namespace MixTransport
363
364#endif //__MATERIALUNSATURATEDFLOW_HPP__
static Index< 'o', 3 > o
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'm', SPACE_DIM > m
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
double initialPcEval() const
Initialize head.
double hS
minimum capillary height [m]
boost::function< boost::shared_ptr< CommonMaterialData >(const CommonMaterialData &data)> RegisterHook
void printMatParameters(const int id, const std::string &prefix) const
double Ks
Saturated hydraulic conductivity [m/day].
double Ah
Initial hydraulic head coefficient.
double AhZZ
Initial hydraulic head coefficient.
double thetaS
saturated water content
double AhZ
Initial hydraulic head coefficient.
void addOptions(po::options_description &o, const std::string &prefix)
double thetaR
residual water contents
Generic material model for unsaturated water transport.
double diffC
Derivative of capacity [S^2/L^2 * L^2/F ].
static double ePsilon0
Regularization parameter.
double sCale
Scale time dependent eq.
double K
Hydraulic conductivity [L/s].
double C
Capacity [S^2/L^2].
double Se
Effective saturation.
static double scaleZ
Scale z direction.
static double ePsilon1
Regularization parameter.
double diffK
Derivative of hydraulic conductivity [L/s * L^2/F].
MaterialDarcy(const CommonMaterialData &data)
static boost::shared_ptr< CommonMaterialData > createMatPtr(const CommonMaterialData &data)
static boost::shared_ptr< CommonMaterialData > createMatPtr(const CommonMaterialData &data)
TYPE funFunSeStar(TYPE &SeStar, const double m)
TYPE funTheta(TYPE &h, const double m)
MaterialVanGenuchten(const CommonMaterialData &data)
void printTheta(const double b, const double e, double s, const std::string &prefix)
void printC(const double b, const double e, double s, const std::string &prefix)
void printKappa(const double b, const double e, double s, const std::string &prefix)
MaterialWithAutomaticDifferentiation(const CommonMaterialData &data)
static map< std::string, CommonMaterialData::RegisterHook > mapOfRegistredMaterials