v0.8.23
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  //
44  double Ah; ///< Initial hydraulic head coefficient
45  double AhZ; ///< Initial hydraulic head coefficient
46  double AhZZ; ///< Initial hydraulic head coefficient
47 
48  double initalPcEval() const { return Ah + AhZ * z + AhZZ * z * z; }
49 
50  void addOptions(po::options_description &o, const std::string &prefix) {
51  o.add_options()((prefix + ".material_name").c_str(),
52  po::value<std::string>(&matName)->default_value(matName))(
53  (prefix + ".ePsilon0").c_str(),
54  po::value<double>(&ePsilon0)->default_value(ePsilon0))(
55  (prefix + ".ePsilon1").c_str(),
56  po::value<double>(&ePsilon1)->default_value(ePsilon1))(
57  (prefix + ".sCale").c_str(),
58  po::value<double>(&sCale)->default_value(sCale))(
59  (prefix + ".thetaS").c_str(),
60  po::value<double>(&thetaS)->default_value(thetaS))(
61  (prefix + ".thetaR").c_str(),
62  po::value<double>(&thetaR)->default_value(thetaR))(
63  (prefix + ".thetaM").c_str(),
64  po::value<double>(&thetaM)->default_value(thetaM))(
65  (prefix + ".alpha").c_str(),
66  po::value<double>(&alpha)->default_value(alpha))(
67  (prefix + ".n").c_str(), po::value<double>(&n)->default_value(n))(
68  (prefix + ".hS").c_str(), po::value<double>(&hS)->default_value(hS))(
69  (prefix + ".Ks").c_str(), po::value<double>(&Ks)->default_value(Ks))(
70  (prefix + ".Ah").c_str(), po::value<double>(&Ah)->default_value(Ah))(
71  (prefix + ".AhZ").c_str(), po::value<double>(&AhZ)->default_value(AhZ))(
72  (prefix + ".AhZZ").c_str(),
73  po::value<double>(&AhZZ)->default_value(AhZZ));
74  }
75 
76  void printMatParameters(const int id, const std::string &prefix) const {
77  PetscPrintf(PETSC_COMM_WORLD, "Mat name %s-%s block id %d\n",
78  prefix.c_str(), matName.c_str(), id);
79  PetscPrintf(PETSC_COMM_WORLD, "Material name: %s\n", matName.c_str());
80  PetscPrintf(PETSC_COMM_WORLD, "thetaS=%6.4g\n", thetaS);
81  PetscPrintf(PETSC_COMM_WORLD, "thetaR=%6.4g\n", thetaR);
82  PetscPrintf(PETSC_COMM_WORLD, "thetaM=%6.4g\n", thetaM);
83  PetscPrintf(PETSC_COMM_WORLD, "alpha=%6.4g\n", alpha);
84  PetscPrintf(PETSC_COMM_WORLD, "n=%6.4g\n", n);
85  PetscPrintf(PETSC_COMM_WORLD, "hS=%6.4g\n", hS);
86  PetscPrintf(PETSC_COMM_WORLD, "Ks=%6.4g\n", Ks);
87  PetscPrintf(PETSC_COMM_WORLD, "Ah=%6.4g\n", Ah);
88  PetscPrintf(PETSC_COMM_WORLD, "AhZ=%6.4g\n", AhZ);
89  PetscPrintf(PETSC_COMM_WORLD, "AhZZ=%6.4g\n", AhZZ);
90  PetscPrintf(PETSC_COMM_WORLD, "ePsilon0=%6.4g\n", ePsilon0);
91  PetscPrintf(PETSC_COMM_WORLD, "ePsilon1=%6.4g\n", ePsilon1);
92  PetscPrintf(PETSC_COMM_WORLD, "sCale=%6.4g\n", sCale);
93  }
94 
95  typedef boost::function<boost::shared_ptr<CommonMaterialData>(
96  const CommonMaterialData &data)>
98 };
99 
101 
102  static boost::shared_ptr<CommonMaterialData>
104  return boost::shared_ptr<CommonMaterialData>(new MaterialDarcy(data));
105  }
106 
108 
111  K = Ks;
113  };
114 
117  diffK = 0;
119  };
120 
123  C = ePsilon1;
125  }
126 
129  diffC = 0;
131  }
132 
135  tHeta = thetaS;
137  }
138 
139  virtual MoFEMErrorCode calSe() {
141  Se = 1;
143  }
144 };
145 
148  : CommonMaterialData(data) {}
149 
150  template <typename TYPE> inline TYPE funSe(TYPE &theta) {
151  return (theta - thetaR) / (thetaS - thetaR);
152  }
153 
154  virtual void recordTheta() = 0;
155 
156  virtual void recordKr() = 0;
157 
158  double Kr;
161  if (h < hS) {
162  int r = ::function(2 * blockId + 1, 1, 1, &h, &Kr);
163  if (r < 0) {
164  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
165  "ADOL-C function evaluation with error");
166  }
167  K = Ks * Kr;
168  } else {
169  K = Ks;
170  }
171  K += Ks * ePsilon0;
173  };
174 
175  double diffKr;
178  if (h < hS) {
179  diffK = 0;
180  int r = ::gradient(2 * blockId + 1, 1, &h, &diffKr);
181  if (r < 0) {
182  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
183  "ADOL-C function evaluation with error");
184  }
185  diffK = Ks * diffKr;
186  } else {
187  diffK = 0;
188  }
190  };
191 
194  if (h < hS) {
195  int r = ::gradient(2 * blockId + 0, 1, &h, &C);
196  if (r < 0) {
197  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
198  "ADOL-C function evaluation with error");
199  }
200  } else {
201  C = 0;
202  }
203  C += ePsilon1;
205  }
206 
209  if (h < hS) {
210  double v = 1;
211  int r = ::hess_vec(2 * blockId + 0, 1, &h, &v, &diffC);
212  if (r < 0) {
213  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
214  "ADOL-C function evaluation with error");
215  }
216  } else {
217  diffC = 0;
218  }
220  }
221 
224  if (h < hS) {
225  int r = ::function(2 * blockId + 0, 1, 1, &h, &tHeta);
226  if (r < 0) {
227  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
228  "ADOL-C function evaluation with error");
229  }
230  } else {
231  tHeta = thetaS;
232  }
234  }
235 
238  if (h < hS) {
239  int r = ::function(2 * blockId + 0, 1, 1, &h, &tHeta);
240  Se = funSe(tHeta);
241  if (r < 0) {
242  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
243  "ADOL-C function evaluation with error");
244  }
245  } else {
246  tHeta = thetaS;
247  }
249  }
250 };
251 
253 
254  static boost::shared_ptr<CommonMaterialData>
256  return boost::shared_ptr<CommonMaterialData>(
257  new MaterialVanGenuchten(data));
258  }
259 
262  recordTheta();
263  recordKr();
264  }
265 
271 
272  template <typename TYPE> inline TYPE funTheta(TYPE &h, const double m) {
273  return thetaR + (thetaM - thetaR) / pow(1 + pow(-alpha * h, n), m);
274  }
275 
276  template <typename TYPE>
277  inline TYPE funFunSeStar(TYPE &SeStar, const double m) {
278  return pow(1 - pow(SeStar, 1 / m), m);
279  }
280 
281  inline adouble funKr(adouble &ah) {
282  const double m = 1 - 1 / n;
283  aTheta = funTheta(ah, m);
284  aSe = funSe(aTheta);
285  aSeStar = aSe * (thetaS - thetaR) / (thetaM - thetaR);
286  double one = 1;
287  const double c = funFunSeStar<double>(one, m);
288  return sqrt(aSe) *
289  pow((1 - funFunSeStar<adouble>(aSeStar, m)) / (1 - c), 2);
290  }
291 
292  virtual void recordTheta() {
293  h = -1 - hS;
294  trace_on(2 * blockId + 0, true);
295  ah <<= h;
296  const double m = 1 - 1 / n;
297  aTheta = funTheta(ah, m);
298  double r_theta;
299  aTheta >>= r_theta;
300  trace_off();
301  }
302 
303  virtual void recordKr() {
304  h = -1 - hS;
305  trace_on(2 * blockId + 1, true);
306  ah <<= h;
307  aKr = funKr(ah);
308  double r_Kr;
309  aKr >>= r_Kr;
310  trace_off();
311  }
312 
313  void printTheta(const double b, const double e, double s,
314  const std::string &prefix) {
315  const double m = 1 - 1 / n;
316  h = b;
317  for (; h >= e; h += s) {
318  s = -pow(-s, 0.9);
319  double theta = funTheta(h, m);
320  double Se = (theta - thetaR) / (thetaS - thetaR);
321  PetscPrintf(PETSC_COMM_SELF, "%s %6.4e %6.4e %6.4e\n", prefix.c_str(), h,
322  theta, Se);
323  }
324  }
325 
326  void printKappa(const double b, const double e, double s,
327  const std::string &prefix) {
328  h = b;
329  for (; h >= e; h += s) {
330  s = -pow(-s, 0.9);
331  calK();
332  PetscPrintf(PETSC_COMM_SELF, "%s %6.4e %6.4e %6.4e\n", prefix.c_str(), h,
333  Kr, K);
334  }
335  }
336 
337  void printC(const double b, const double e, double s,
338  const std::string &prefix) {
339  h = b;
340  for (; h >= e; h += s) {
341  s = -pow(-s, 0.9);
342  calC();
343  PetscPrintf(PETSC_COMM_SELF, "%s %6.4e %6.4e\n", prefix.c_str(), h, C);
344  }
345  }
346 };
347 
349  static map<std::string, CommonMaterialData::RegisterHook>
354  mapOfRegistredMaterials["VanGenuchten"] =
357  }
358 };
359 
360 } // namespace MixTransport
361 
362 #endif //__MATERIALUNSATURATEDFLOW_HPP__
double thetaS
saturated water content
double diffK
Derivative of hydraulic conductivity [L/s * L^2/F].
static boost::shared_ptr< CommonMaterialData > createMatPtr(const CommonMaterialData &data)
double initalPcEval() const
Initialize head.
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:500
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:507
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.
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