Flow123d  master-f44eb46
soil_models.cc
Go to the documentation of this file.
1 /*
2  * soil_models.cc
3  *
4  * Created on: Aug 11, 2016
5  * Author: jb
6  */
7 
8 
9 #include <algorithm>
10 #include <cmath>
11 
12 #include "system/logger.hh"
13 
14 #include "system/asserts.hh"
15 #include "flow/soil_models.hh"
16 #include "tools/include_fadbad.hh" // for "fadbad.h", "badiff.h", "fadiff.h"
17 
18 template <class Model>
20 : cut_fraction_(cut_fraction)
21 {}
22 
23 
24 template <class Model>
26 {
27  data.cut_fraction = cut_fraction_;
28  model_.reset_(data);
29  /*
30  // Check table of hydraulic functions
31  static int checked=0;
32  if (! checked) {
33  checked=1;
34 
35  double factor=pow(10, 0.1);
36  for(double h=0.1; h < 1e4; h*=factor) {
37  LogOut().fmt("h: {:g} sat: {:g} cond: {:g}", -h, water_content(-h), conductivity(-h));
38  }
39  }
40  */
41 }
42 
43 template <class Model>
44 double SoilModelImplBase<Model>::conductivity( const double &p_head) const
45 {
46  return model_.conductivity_(p_head);
47 }
48 
49 template <class Model>
50 auto SoilModelImplBase<Model>::conductivity_diff(const DiffDouble &p_head)->DiffDouble const
51 {
52  return model_.conductivity_(p_head);
53 }
54 
55 template <class Model>
56 double SoilModelImplBase<Model>::water_content( const double &p_head) const
57 {
58  return model_.water_content_(p_head);
59 }
60 
61 template <class Model>
62 auto SoilModelImplBase<Model>::water_content_diff(const DiffDouble &p_head)->DiffDouble const
63 {
64  return model_.water_content_(p_head);
65 }
66 
67 
68 
69 
70 namespace internal {
71 
73 : Bpar(0.5), Ppar(2), K_lower_limit(1.0E-20)
74 {}
75 
77 {
78  // check soil parameters
79  ASSERT_LT( soil.cut_fraction, 1.0);
80  ASSERT_GT( soil.cut_fraction, 0.0);
81  soil_param_ = soil;
82 
83  m = 1-1/soil_param_.n;
84  // soil_param_.cut_fraction = (Qs -Qr)/(Qnc - Qr)
86 
87  // conductivity internal scaling
88  FFQr = 1.0; // pow(1 - pow(Qeer,1/m),m);
89  double Qs_relative = soil_param_.cut_fraction;
90  FFQs = pow(1 - pow(Qs_relative, 1/m),m);
91 
92 
94  /*
95  std::cout << "Hs : " << Hs
96  << " qs: " << soil_param_.Qs
97  << " qsnc: " << Qs_nc
98  << " Q: " << Q_rel(Hs) << std::endl;
99  */
100 }
101 
102 
103 template <class T>
104 T VanGenuchten::Q_rel(const T &h) const
105 {
106  return pow( 1 + pow(-soil_param_.alpha * h, soil_param_.n), -m);
107 }
108 
109 template <class T>
110 T VanGenuchten::Q_rel_inv(const T &q) const
111 {
112  return -pow( pow( q, -1/m ) -1, 1/soil_param_.n)/soil_param_.alpha;
113 }
114 
115 
116 
117 // pri generovani grafu jsem zjistil ze originalni vzorecek pro
118 // vodivost pres theta je numericky nestabilni pro tlaky blizke
119 // nule, zejmena pro velka n
120 // tento vzorec je stabilni:
121 //
122 /// k(h)=t(h)**(0.5)* (1- ((h)**n/(1+(h)**n)) **m)**2
123 
124 template <class T>
125 T VanGenuchten::conductivity_(const T& h) const
126 {
127  T Kr, Q_unscaled, Q_cut_unscaled, FFQ;
128 
129  if (h < Hs) {
130  Q_unscaled = Q_rel(h);
131  Q_cut_unscaled = Q_unscaled / soil_param_.cut_fraction;
132  FFQ = pow(1 - pow(Q_unscaled,1/m),m);
133  Kr = soil_param_.Ks * pow(Q_cut_unscaled,Bpar)*pow((FFQr - FFQ)/(FFQr - FFQs),Ppar);
134  }
135  else Kr = soil_param_.Ks;
136 
137  if (Kr < K_lower_limit) return K_lower_limit;
138  else return Kr;
139 }
140 
141 //template SoilModelBase::DiffDouble VanGenuchten::conductivity_<SoilModelBase::DiffDouble>(const SoilModelBase::DiffDouble &h) const;
142 //template double VanGenuchten::conductivity_<double>(const double &h) const;
143 
144 
145 template <class T>
146 T VanGenuchten::water_content_(const T& h) const
147 {
148  if (h < Hs) return soil_param_.Qr + (Qs_nc - soil_param_.Qr) *Q_rel(h);
149  else return soil_param_.Qs;
150 }
151 
152 //template SoilModelBase::DiffDouble VanGenuchten::water_content_<SoilModelBase::DiffDouble>(const SoilModelBase::DiffDouble &h) const;
153 //template double VanGenuchten::water_content_<double>(const double &h) const;
154 
156 : VanGenuchten()
157 {
158  Ppar=3;
159 }
160 
161 
162 template <class T>
163 T Irmay::conductivity_(const T& h) const
164 {
165  T Kr;
166 
167  if (h < Hs) {
168  T Q = this->Q_rel(h);
169  Kr = soil_param_.Ks * pow(Q, Ppar);
170  }
171  else Kr = soil_param_.Ks;
172 
173  if (Kr < K_lower_limit) return K_lower_limit;
174  else return Kr;
175 }
176 
177 //template SoilModelBase::DiffDouble Irmay::conductivity_<SoilModelBase::DiffDouble>(const SoilModelBase::DiffDouble &h) const;
178 //template double Irmay::conductivity_<double>(const double &h) const;
179 
180 } //close namespace internal
181 
182 
183 
186 
187 
Definitions of ASSERTS.
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
Definition: asserts.hh:317
fadbad::B< double > DiffDouble
Definition: soil_models.hh:60
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
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
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