v0.9.1
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 
11 namespace 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 
113  K = Ks;
115  };
116 
119  diffK = 0;
121  };
122 
125  C = ePsilon1;
127  }
128 
131  diffC = 0;
133  }
134 
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;
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;
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 
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 
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 
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 
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 
283  inline adouble funKr(adouble &ah) {
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>
356  mapOfRegistredMaterials["VanGenuchten"] =
359  }
360 };
361 
362 } // namespace MixTransport
363 
364 #endif //__MATERIALUNSATURATEDFLOW_HPP__
static double scaleZ
Scale z direction.
double thetaS
saturated water content
double initialPcEval() const
Initialize head.
double diffK
Derivative of hydraulic conductivity [L/s * L^2/F].
static boost::shared_ptr< CommonMaterialData > createMatPtr(const CommonMaterialData &data)
double diffC
Derivative of capacity [S^2/L^2 * L^2/F ].
static boost::shared_ptr< CommonMaterialData > createMatPtr(const CommonMaterialData &data)
boost::function< boost::shared_ptr< CommonMaterialData > const CommonMaterialData &data)> RegisterHook
double C
Capacity [S^2/L^2].
static double ePsilon0
Regularization parameter.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
double Ah
Initial hydraulic head coefficient.
void addOptions(po::options_description &o, const std::string &prefix)
double Ks
Saturated hydraulic conductivity [m/day].
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
MaterialVanGenuchten(const CommonMaterialData &data)
double hS
minimum capillary height [m]
double AhZZ
Initial hydraulic head coefficient.
TYPE funFunSeStar(TYPE &SeStar, const double m)
MaterialDarcy(const CommonMaterialData &data)
void printKappa(const double b, const double e, double s, const std::string &prefix)
static double ePsilon1
Regularization parameter.
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
void printTheta(const double b, const double e, double s, const std::string &prefix)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
double K
Hydraulic conductivity [L/s].
MaterialWithAutomaticDifferentiation(const CommonMaterialData &data)
double Se
Effective saturation.
TYPE funTheta(TYPE &h, const double m)
double sCale
Scale time dependent eq.
void printC(const double b, const double e, double s, const std::string &prefix)
static map< std::string, CommonMaterialData::RegisterHook > mapOfRegistredMaterials
void printMatParameters(const int id, const std::string &prefix) const
Generic material model for unsaturated water transport.
double AhZ
Initial hydraulic head coefficient.
double thetaR
residual water contents