Flow123d  DF_profiler_memory_monitor-0108f36
isotherm.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file isotherm.cc
15  * @brief
16  */
17 
18 #include "reaction/isotherm.hh"
19 #include "system/sys_profiler.hh"
20 #include "system/logger.hh"
21 
22 #include <boost/core/explicit_operator_bool.hpp> // for optiona...
23 
24 #include <boost/exception/exception.hpp> // for exception
25 
26 #include <boost/format/alt_sstream.hpp> // for basic_a...
27 #include <boost/format/alt_sstream_impl.hpp> // for basic_a...
28 #include <boost/optional/optional.hpp> // for get_poi...
29 #include <boost/exception/diagnostic_information.hpp> // for diagnos...
30 #include <boost/math/tools/toms748_solve.hpp> // for toms748...
31 
33 : table_limit_(0.0)
34 {
35 }
36 
38 {
39  table_limit_ = 0.0;
40  interpolation_table.clear();
41 }
42 
43 
44 void Isotherm::reinit(enum SorptionType adsorption_type, bool limited_solubility_on,
45  double rho_aqua, double scale_aqua, double scale_sorbed,
46  double c_aqua_limit, double mult_coef, double second_coef)
47 {
48  adsorption_type_ = adsorption_type;
49  rho_aqua_ = rho_aqua;
50  scale_aqua_ = scale_aqua;
51  scale_sorbed_ = scale_sorbed;
52  limited_solubility_on_ = limited_solubility_on;
53  mult_coef_ = mult_coef*rho_aqua;
54  second_coef_ = second_coef;
55 
56  solubility_limit_ = c_aqua_limit;
59 }
60 
61 
63 {
64  return scale_aqua_* conc.fluid + scale_sorbed_ * conc.solid;
65 }
66 
67 
68 void Isotherm::compute( double &c_aqua, double &c_sorbed )
69 {
70  // if sorption is switched off, do not compute anything
72  return;
73 
74  ConcPair c_pair(c_aqua, c_sorbed);
75  ConcPair result(0,0);
76 
77  result = solve_conc( c_pair );
78 
79  c_aqua=result.fluid;
80  c_sorbed=result.solid;
81 }
82 
83 
84 void Isotherm::interpolate( double &c_aqua, double &c_sorbed )
85 {
86  // if sorption is switched off, do not compute anything
88  return;
89 
90  ConcPair c_pair(c_aqua, c_sorbed);
91  ConcPair result(0,0);
92 
93  result = compute_projection( c_pair );
94 
95  c_aqua=result.fluid;
96  c_sorbed=result.solid;
97 }
98 
99 
101 {
102  double total_mass = get_total_mass(c_pair);
103 // DebugOut().fmt("compute_projection: total mass = {}, c_aqua = {}\n", total_mass, c_pair.fluid);
104  if(total_mass < 0.0)
105  THROW( Isotherm::ExcNegativeTotalMass()
106  << EI_TotalMass(total_mass)
107  );
108  // total_mass_step_ is set and checked in make_table
109  double total_mass_steps = total_mass / total_mass_step_;
110  unsigned int total_mass_idx = static_cast <unsigned int>(std::floor(total_mass_steps));
111 
112  if (total_mass_idx < (interpolation_table.size() - 1) ) {
113  double rot_sorbed = interpolation_table[total_mass_idx]
114  + (total_mass_steps - total_mass_idx)*(interpolation_table[total_mass_idx+1]
115  - interpolation_table[total_mass_idx]);
116  return ConcPair( (total_mass * inv_scale_aqua_ - rot_sorbed * inv_scale_sorbed_),
117  (total_mass * inv_scale_sorbed_ + rot_sorbed * inv_scale_aqua_) );
118  } else {
120  return precipitate( c_pair );
121  else
122  return solve_conc( c_pair );
123  }
124 }
125 
126 
128 {
129  double total_mass = get_total_mass(c_pair);
130 // DebugOut().fmt("precipitate: total mass = {}, c_aqua = {}\n", total_mass, c_pair.fluid);
132  (total_mass - scale_aqua_ * solubility_limit_) / scale_sorbed_ );
133 }
134 
135 
136 template<class Func>
138 {
139  double total_mass = get_total_mass(c_pair);
140  double mass_limit = get_total_mass(Isotherm::ConcPair(solubility_limit_, const_cast<Func &>(isotherm)(solubility_limit_ / this->rho_aqua_)));
141 
142  // condition on limited solubility in the rotated coordinate system (total mass)
143 // DebugOut().fmt("total_mass {}, mass_limit {} \n",total_mass, mass_limit);
144  if (total_mass > mass_limit){
146  return precipitate( c_pair );
147  else
148  // if solubility is not limited, increase mass limit
149  mass_limit = total_mass;
150  }
151 
152  double upper_solution_bound = mass_limit / scale_aqua_;
153  CrossFunction<Func> eq_func(isotherm, total_mass, scale_aqua_, scale_sorbed_, this->rho_aqua_);
154  std::pair<double,double> solution;
155  if (total_mass > 0) // here should be probably some kind of tolerance instead of "0"
156  {
157  try {
158  boost::uintmax_t max_iter = 20;
159  tolerance<double> toler(30);
160  solution = boost::math::tools::toms748_solve(eq_func, 0.0, upper_solution_bound, toler, max_iter);
161  }
162  catch(boost::exception const & e)
163  {
164  THROW( Isotherm::ExcBoostSolver()
165  << EI_BoostMessage(boost::diagnostic_information(e))
166  );
167  }
168  }
169  double difference;
170  difference = (solution.second - solution.first)/2;
171  double c_aqua = solution.first + difference;
172  return ConcPair(c_aqua,
173  (total_mass - scale_aqua_ * c_aqua)/scale_sorbed_);
174 }
175 
176 // Isotherm None specialization
178 {
179  return c_pair;
180 }
181 
182 
184 {
185  switch(adsorption_type_)
186  {
187  case 0: // none
188 // {
189 // None obj_isotherm;
190 // return solve_conc( conc, obj_isotherm);
191 // }
192  break;
193  case 1: // linear:
194  {
195  Linear obj_isotherm(mult_coef_);
196  return solve_conc( conc, obj_isotherm);
197  }
198  break;
199  case 2: // freundlich
200  {
201  Freundlich obj_isotherm(mult_coef_, second_coef_);
202  return solve_conc( conc, obj_isotherm);
203  }
204  break;
205  case 3: // langmuir:
206  {
207  Langmuir obj_isotherm(mult_coef_, second_coef_);
208  return solve_conc( conc, obj_isotherm);
209  }
210  break;
211  }
212  return conc;
213 }
214 
215 
216 template<class Func>
217 void Isotherm::make_table( const Func &isotherm, int n_steps )
218 {
219  // limit aqueous concentration for the interpolation table; cannot be higher than solubility limit
220  double aqua_limit = table_limit_;
222  aqua_limit = solubility_limit_;
223 
224  double mass_limit = scale_aqua_ * aqua_limit + scale_sorbed_ * const_cast<Func &>(isotherm)(aqua_limit / this->rho_aqua_);
225 // DebugOut().fmt("make_table: mass_limit = {}, aqua_limit = {}\n", mass_limit, aqua_limit);
226 
227  if(mass_limit < 0.0)
228  THROW( Isotherm::ExcNegativeTotalMass()
229  << EI_TotalMass(mass_limit)
230  );
231  total_mass_step_ = mass_limit / n_steps;
232  double mass = 0.0;
233  interpolation_table.clear();
234  interpolation_table.reserve(n_steps);
235  for(int i=0; i<= n_steps; i++) {
236  // aqueous concentration (original coordinates c_a) corresponding to i-th total_mass_step_
237  ConcPair c_pair( mass/scale_aqua_, 0.0 );
238 
239  ConcPair result = solve_conc( c_pair, isotherm);
240  double c_sorbed_rot = ( result.solid * scale_aqua_ - result.fluid * scale_sorbed_);
241  interpolation_table.push_back(c_sorbed_rot);
242  mass = mass+total_mass_step_;
243  }
244 }
245 
246 // Isotherm None specialization
247 template<> void Isotherm::make_table( const None &, int )
248 {
249  // Solve_conc returns the same, so we need to do that also in compute_projection.
250  // We set size of the table to 1, so it follow the conditions into solve_conc again.
251 
252  limited_solubility_on_ = false; // so it cannot go in precipitate function
253 
254  total_mass_step_ = 1; // set just one step in the table, so we void zero division
255  interpolation_table.resize(1,0); // set one value in the table so the condition in compute_projection fails
256  return;
257 }
258 
259 void Isotherm::make_table( unsigned int n_points, double table_limit )
260 {
261  START_TIMER("Isotherm::make_table");
263  if(table_limit_ > 0.0)
264  switch(adsorption_type_)
265  {
266  case 0: // none
267  {
268  None obj_isotherm;
269  make_table(obj_isotherm, 1);
270  }
271  break;
272  case 1: // linear:
273  {
274  Linear obj_isotherm(mult_coef_);
275  make_table(obj_isotherm, n_points);
276  }
277  break;
278  case 2: // freundlich:
279  {
280  Freundlich obj_isotherm(mult_coef_, second_coef_);
281  make_table(obj_isotherm, n_points);
282  }
283  break;
284  case 3: // langmuir:
285  {
286  Langmuir obj_isotherm(mult_coef_, second_coef_);
287  make_table(obj_isotherm, n_points);
288  }
289  break;
290  default:
291  break;
292  }
293  else
294  clear_table();
295 }
double solubility_limit_
Concentration limit in liquid phase (solubility limit).
Definition: isotherm.hh:295
ConcPair precipitate(ConcPair conc)
Definition: isotherm.cc:127
void make_table(unsigned int n_points, double table_limit)
Definition: isotherm.cc:259
double get_total_mass(ConcPair conc)
Definition: isotherm.cc:62
ConcPair compute_projection(ConcPair conc)
Definition: isotherm.cc:100
double inv_scale_sorbed_
Definition: isotherm.hh:304
ConcPair solve_conc(ConcPair c_pair, const Func &isotherm)
Definition: isotherm.cc:137
void compute(double &c_aqua, double &c_sorbed)
Definition: isotherm.cc:68
bool limited_solubility_on_
Solubility limit flag.
Definition: isotherm.hh:293
double mult_coef_
Multiplication parameter of the isotherm.
Definition: isotherm.hh:284
double table_limit(void) const
Getter for table limit (limit aqueous concentration).
Definition: isotherm.hh:244
Isotherm()
Default constructor.
Definition: isotherm.cc:32
SorptionType
Type of adsorption isotherm.
Definition: isotherm.hh:174
double table_limit_
Concentration in liquid phase for limit of the interpolation table.
Definition: isotherm.hh:290
double second_coef_
Optional secod parameter of the isotherm.
Definition: isotherm.hh:287
void interpolate(double &c_aqua, double &c_sorbed)
Definition: isotherm.cc:84
void reinit(enum SorptionType sorption_type, bool limited_solubility_on, double aqua_density, double scale_aqua, double scale_sorbed, double c_aqua_limit, double mult_coef, double second_coef)
Definition: isotherm.cc:44
double rho_aqua_
density of the solvent
Definition: isotherm.hh:298
double total_mass_step_
Definition: isotherm.hh:313
std::vector< double > interpolation_table
Definition: isotherm.hh:309
double scale_aqua_
coefficient that convert soluted concentration to mass; porosity = k_W, originally rho_aqua*porosity ...
Definition: isotherm.hh:300
enum SorptionType adsorption_type_
Type of isotherm.
Definition: isotherm.hh:281
double inv_scale_aqua_
reciprocal values
Definition: isotherm.hh:304
double scale_sorbed_
coefficient that convert adsorbed molar concentration to mass; molar_weight * rho_rock * (1 - porosit...
Definition: isotherm.hh:302
void clear_table()
Definition: isotherm.cc:37
Definition: isotherm.hh:88
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
Pair of soluted and adsorbed concentration.
Definition: isotherm.hh:182
#define START_TIMER(tag)
Starts a timer with specified tag.