Flow123d  master-f44eb46
soil_models.hh
Go to the documentation of this file.
1 /*
2  * File: hydro_functions.hh
3  * Author: jb
4  *
5  * Created on May 17, 2010, 1:52 PM
6  */
7 
8 #ifndef _HYDRO_FUNCTIONS_HH
9 #define _HYDRO_FUNCTIONS_HH
10 
11 
12 // #include "badiff.h" // for B::d, B::deriv, B::diff, B::getBTypeNameHV, B::o...
13 // #include "fadbad.h" // for B
14 #include "tools/include_fadbad.hh" // for "fadbad.h", "badiff.h"
15 namespace internal { class Irmay; }
16 namespace internal { class VanGenuchten; }
17 
18 
19 // todo:
20 // objekt por SoilModel - jedna sada primarnich i pomocnych parametru
21 // konzistentni hydrologicke funkce
22 //
23 // potreboval bych mit vice modelu, z orezavanim i bez nej,
24 // orezavani zpusobuje velke residuum v miste nespojitosti derivace saturace
25 
26 // problem - jak pouzit funktory, FADBAD, aproximace a
27 // inverzi (inverze tlak v zav. na vlhkosti)
28 //
29 // podrobnejsi vyzkum ohledne hydrologickych funkci
30 // vodivost vyplyva s kapilarnich modelu (jak presne?)
31 // co vlhkost?
32 
33 /**
34  * One of possibly more classes specifying dependence of the relative conductivity and the relative water content on the
35  * pressure head. These functions (@p conductivity and @p water_content, respectively) are templates to allow evaluation of
36  * the derivatives through fadbad lib.
37  */
38 
39 struct SoilData {
40  double n; // power parameter
41  double alpha; // pressure head scaling
42  double Qr; // residual water content
43  double Qs; // saturated water content
44  double cut_fraction; // in (0,1), fraction of relative Q where to cut and rescale the retention curve
45  double Ks; // saturated conductivity
46 // double Kk; // conductivity at the cut point
47 };
48 
49 
50 /**
51  * Pure virtual interface, of all models.
52  */
54 public:
57  irmay=1
58  };
59 
60  typedef fadbad::B<double> DiffDouble;
61 
62  virtual void reset(SoilData soil)=0;
63 
64  virtual double conductivity( const double &phead) const =0;
65  virtual auto conductivity_diff(const DiffDouble &p_head)->DiffDouble const =0;
66 
67  virtual double water_content( const double &phead) const =0;
68  virtual auto water_content_diff(const DiffDouble &p_head)->DiffDouble const =0;
69 
70  virtual ~SoilModelBase() {};
71 };
72 
73 
74 template <class Model>
76 public:
77  // We assume that Model have method templates:
78  // template <class T> T conductivity_(const T &h) const;
79  // template <class T> T water_content_(const T &h) const;
80 
81 
82 
84 
85  SoilModelImplBase(double cut_fraction = 0.999);
86 
87  void reset(SoilData data) override;
88 
89  double conductivity( const double &p_head) const override;
90  auto conductivity_diff(const DiffDouble &p_head)->DiffDouble const override;
91 
92  double water_content( const double &p_head) const override;
93  auto water_content_diff(const DiffDouble &p_head)->DiffDouble const override;
94 
96 
97 private:
99  double cut_fraction_;
100 };
101 
102 
103 
104 // ************************************** Internal model definitions
105 namespace internal {
106 
108 public:
109  VanGenuchten();
110 
111  void reset_(SoilData soil);
112 
113  template <class T>
114  T conductivity_(const T &h) const;
115 
116  template <class T>
117  T water_content_(const T &h) const;
118 
119 protected:
120 
121  template <class T> T Q_rel(const T &h) const;
122  template <class T> T Q_rel_inv(const T &q) const;
123 
124  // input parameters
126 
127  // constant conductivity parameters
128  double Bpar;
129  double Ppar;
131 
132  // aux values
133  double m;
134  double Qs_nc; // saturated value of continuation of cut water content function
135 
136 
137  double FFQr, FFQs;
138  double Hs;
139 
140 };
141 
142 
143 
144 class Irmay : public VanGenuchten {
145 public:
146  Irmay();
147 
148  template <class T>
149  T conductivity_(const T &h) const;
150 };
151 
152 } // close namespace internal
153 
154 
157 
158 
159 
160 
161 
162 
163 
164 
165 
166 
167 /*
168 
169 //FK-------------------------------------------------------------
170 template <class T>
171 class FK //FK - hydraulic conductivity function
172 {
173 private:
174 // function parameters
175 static const double Bpar; // = 0.5;
176 static const double PPar; // = 2;
177 
178 HydrologyParams par;
179 
180 // auxiliary parameters
181 double Qa,Qm,Qk,m,Hr,Hk,Hs,C1Qee,C2Qee,Qeer,Qees,Qeek;
182 public:
183  FK(const HydrologyParams &par);
184  T operator()(const T& h);
185 };
186 
187 template <class T> const double FK<T>::Bpar =0.5;
188 template <class T> const double FK<T>::PPar =2;
189 
190 template <class T>
191 FK<T> :: FK(const HydrologyParams &par)
192 : par(par)
193 {
194 
195  m = 1-1/par.n;
196 
197  Qa=par.Qr;
198  Qm=par.Qs / par.cut_fraction;
199  //Qk=0.99 * par.Qs;
200  Qk=par.Qs;
201  C1Qee = 1/(Qm - Qa);
202  C2Qee = -Qa/(Qm - Qa);
203 // if (par.cut_fraction >0.9999) {
204 // Hr=-1000;
205 // Hk=Hs=0;
206 // return;
207 // }
208  // fraction =1
209  Qees = C1Qee * par.Qs + C2Qee; // 1
210  Qeek = min(C1Qee * Qk + C2Qee, Qees); // 1
211  Qeer = min(C1Qee * par.Qr + C2Qee, Qeek); // 0
212 
213  Hr = -1/par.alfa*pow(pow(Qeer,-1/m)-1,1/par.n); //
214  Hs = -1/par.alfa*pow(pow(Qees,-1/m)-1,1/par.n);
215  Hk = -1/par.alfa*pow(pow(Qeek,-1/m)-1,1/par.n);
216 
217  cout << "K: Hr, Hs, Hk:" << Hr << " " << Hs << " " << Hk << endl;
218 }
219 
220 template <class T>
221 T FK<T> :: operator()(const T& h)
222 {
223  T Kr,FFQr,FFQ,FFQk,Qee,Qe,Qek,C1Qe,C2Qe,Q;
224 
225  if (h < Hr) return par.Ks*(1E-200); // numericaly zero Kr
226  else
227  if (h < Hk) {
228  Q = Qa + (Qm - Qa)*pow((1 + pow(-par.alfa*h,par.n)),-m);
229  Qee = C1Qee*Q + C2Qee;
230  FFQr = pow(1 - pow(Qeer,1/m),m);
231  FFQ = pow(1 - pow(Qee,1/m),m);
232  FFQk = pow(1 - pow(Qeek,1/m),m);
233  C1Qe = 1/(par.Qs - par.Qr);
234  C2Qe = -par.Qr/(par.Qs - par.Qr);
235  Qe = C1Qe*Q + C2Qe;
236  Qek = C1Qe*Qk + C2Qe;
237  Kr = pow(Qe/Qek,Bpar)*pow((FFQr - FFQ)/(FFQr - FFQk),PPar) * par.Kk/par.Ks;
238  return ( max<T>(par.Ks*Kr,par.Ks*(1E-10)) );
239  }
240  else if(h <= Hs)
241  {
242  Kr = (1-par.Kk/par.Ks)/(Hs-Hk)*(h-Hs) + 1;
243  return ( par.Ks*Kr );
244  }
245  else return par.Ks;
246 }
247 
248 
249 //FQ--------------------------------------------------------------
250 template <class T>
251 class FQ //FQ - water saturation function
252 {
253 HydrologyParams par;
254 // auxiliary parameters
255 double m, Qeer,Qees,Hr,Hs, Qa, Qm;
256 public:
257  FQ(const HydrologyParams &par);
258  T operator()(const T& h);
259 };
260 
261 template <class T>
262 FQ<T>::FQ(const HydrologyParams &par)
263 : par(par)
264 {
265  m = 1 - 1 / par.n;
266 
267 // if (par.cut_fraction >0.9999) {
268 // Hr=-1000;
269 // Hs=0;
270 // return;
271 // }
272  Qm=par.Qs / par.cut_fraction;
273  Qa=par.Qr;
274 
275  Qeer = (par.Qr - Qa)/ (Qm - Qa);
276  Qees = (par.Qs - Qa)/ (Qm - Qa);
277  Hr = -1 / par.alfa * pow( pow(Qeer,-1/m)-1, 1/par.n);
278  Hs = -1 / par.alfa * pow( pow(Qees,-1/m)-1, 1/par.n);
279 
280  cout << "Q: Hr, Hs:" << Hr << " " << Hs << " " << endl;
281 }
282 
283 
284 template <class T>
285 T FQ<T>::operator()(const T& h)
286 {
287  T Qee;
288 
289  if (h < Hr) return par.Qr;
290  else if (h < Hs) {
291  Qee = pow( 1 + pow(-par.alfa * h, par.n), -m);
292  return Qa + (Qm - Qa) * Qee;
293  }
294  else return par.Qs;
295 
296 }
297 
298 //FC--------------------------------------------------------------
299 template <class T>
300 class FC //FC - water capacity function FQ derivative
301 {
302 HydrologyParams par;
303 // auxiliary parameters
304 double m, C1Qee,C2Qee,Qeer,Qees,Hr,Hs, Qa, Qm;
305 public:
306  FC(const HydrologyParams &par);
307  T operator()(const T& h);
308 };
309 
310 template <class T>
311 FC<T>::FC(const HydrologyParams &par)
312 : par(par)
313 {
314  m = 1 - 1 / par.n;
315  Qm=par.Qs / par.cut_fraction;
316  Qa=par.Qr;
317 
318  C1Qee = 1/(Qm - Qa);
319  C2Qee = -Qa/(Qm - Qa);
320  Qeer = C1Qee * par.Qr + C2Qee;
321  Qees = C1Qee * par.Qs + C2Qee;
322  Hr = -1 / par.alfa * pow( pow(Qeer,-1/m)-1, 1/par.n);
323  Hs = -1 / par.alfa * pow( pow(Qees,-1/m)-1, 1/par.n);
324 }
325 
326 
327 template <class T>
328 T FC<T>::operator()(const T& h)
329 {
330  T C1;
331 
332  if (h <= Hr)
333  return 0.0;
334  else if ( h < Hs ) {
335  C1=pow( 1 + pow( -par.alfa * h, par.n), -m-1);
336  return (Qm - Qa) * m * par.n * pow(par.alfa, par.n)*pow(-h, par.n-1)*C1;
337  }
338  else
339  return 0.0;
340 }
341 
342 
343 ! ************************************************************************
344 ! FH - inverse water saturation function
345 ! **************************************
346 real pure function FH_4(h,Matr)
347 implicit none
348 integer,INTENT(IN) :: Matr
349 real, INTENT(IN) :: h
350  FH_4=FH_8(DBLE(h),Matr)
351 end function FH_4
352 
353 double precision pure function FH_8(Qe,Matr)
354 implicit none
355 integer, INTENT(IN) :: Matr
356 double precision,INTENT(IN) :: Qe
357 double precision :: n,m,Qr,Qs,Qa,Qm,Alfa,Q,Qee
358 
359  Qr=par(1,Matr)
360  Qs=par(2,Matr)
361  Qa=par(3,Matr)
362  Qm=par(4,Matr)
363  Alfa=par(5,Matr)
364  n=par(6,Matr)
365  m=1-1/n
366  Q=Qr+(Qs-Qr)*Qe
367  Qee=dmax1((Q-Qa)/(Qm-Qa),1d-3)
368  Qee=dmin1(Qee,1-1d-15)
369  FH_8=-1/Alfa*(Qee**(-1/m)-1)**(1/n)
370 end function FH_8
371 */
372 
373 // Conductivity for analytical solution
374 /*
375 class HydroModel_analytical
376 {
377 public:
378  HydroModel_analytical(ParameterHandler &prm) {};
379 
380  /// Maximum point of capacity. (computed numericaly from eq. atan(x)*2x=1 )
381  inline double cap_arg_max() const { return -0.765378926665788882857; }
382 
383  template <class T>
384  T FK(const T &h) const
385  {
386  if (h>=0.0) return ( 2.0 );
387  return ( 2.0 / (1+h*h) );
388  }
389 
390  template <class T>
391  T FQ(const T &h) const
392  {
393  static T pi_half =std::atan(1.0)*2;
394  if (h>=0.0) return ( 2.0*pi_half*pi_half );
395  T a_tan = atan(h);
396  return ( 2*pi_half*pi_half - a_tan*a_tan );
397  }
398 private:
399  double arg_cap_max;
400 };
401 
402 
403 class HydroModel_linear
404 {
405  HydroModel_linear(ParameterHandler &prm) {};
406 
407 public:
408  template <class T>
409  T FK(const T &h) const {
410  return 1.0;
411  }
412 
413  template <class T>
414  T FQ(const T &h) const {
415  return h;
416  }
417 };
418 
419 */
420 #endif /* _HYDRO_FUNCTIONS_HH */
421 
virtual ~SoilModelBase()
Definition: soil_models.hh:70
virtual auto water_content_diff(const DiffDouble &p_head) -> DiffDouble const =0
virtual auto conductivity_diff(const DiffDouble &p_head) -> DiffDouble const =0
virtual double water_content(const double &phead) const =0
fadbad::B< double > DiffDouble
Definition: soil_models.hh:60
virtual void reset(SoilData soil)=0
virtual double conductivity(const double &phead) const =0
auto water_content_diff(const DiffDouble &p_head) -> DiffDouble const override
Definition: soil_models.cc:62
auto conductivity_diff(const DiffDouble &p_head) -> DiffDouble const override
Definition: soil_models.cc:50
double conductivity(const double &p_head) const override
Definition: soil_models.cc:44
double water_content(const double &p_head) const override
Definition: soil_models.cc:56
void reset(SoilData data) override
Definition: soil_models.cc:25
SoilModelBase::DiffDouble DiffDouble
Definition: soil_models.hh:83
SoilModelImplBase(double cut_fraction=0.999)
Definition: soil_models.cc:19
T conductivity_(const T &h) const
Definition: soil_models.cc:163
T water_content_(const T &h) const
Definition: soil_models.cc:146
void reset_(SoilData soil)
Definition: soil_models.cc:76
T Q_rel_inv(const T &q) const
Definition: soil_models.cc:110
T Q_rel(const T &h) const
Definition: soil_models.cc:104
T conductivity_(const T &h) const
k(h)=t(h)**(0.5)* (1- ((h)**n/(1+(h)**n)) **m)**2
Definition: soil_models.cc:125
SoilModelImplBase< internal::VanGenuchten > SoilModel_VanGenuchten
Definition: soil_models.hh:155
SoilModelImplBase< internal::Irmay > SoilModel_Irmay
Definition: soil_models.hh:156
double Ks
Definition: soil_models.hh:45
double Qr
Definition: soil_models.hh:42
double n
Definition: soil_models.hh:40
double alpha
Definition: soil_models.hh:41
double cut_fraction
Definition: soil_models.hh:44
double Qs
Definition: soil_models.hh:43