Flow123d  release_2.2.0-48-gb04af7f
isotherm.hh
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.hh
15  * @brief
16  *
17  * Other possible transformation of coordinates:
18  *
19  * c_l - conc. liquid
20  * c_s - conc. solid
21  * c_t = h_l * c_l + h_s * c_s = X' + Y'
22  * X = c_t
23  * Y = X' - Y' = h_l * c_l - h_s * c_s
24  *
25  * A) make table for function c_t -> Y
26  * 1) for given c_t solve nonlinear eq.
27  * h_l * c_l + h_s * f(c_l) = c_t
28  *
29  * 2) from c_l compute
30  * Y = c_t - 2* h_l * c_l
31  *
32  *
33  * B) calculation of new c_l, c_s from c_t using table:
34  *
35  * 1) use table to get value of Y for given c_t
36  * 2) compute:
37  * c_l = (c_t - Y) / (2 * h_l)
38  * c_s = (c_t + Y) / (2 * h_s)
39  *
40  * ========================
41  * The transformation currently in use transforms
42  * pair (c_l, c_s) directly to (c_t, W) in ortogonal way.
43  * Proposed transformation first scale (c_l, c_s) to (X',Y')
44  * and then transform scaled coordinates in ortogonal way.
45  */
46 
47 #ifndef SORPTION_IMPL_HH_
48 #define SORPTION_IMPL_HH_
49 
50 #include <vector>
52 #include <boost/math/tools/roots.hpp>
53 #include "fields/field.hh"
54 
55 
56 
57 /**
58  * Convergence criteria for interval based nonlinear solver. It is functor, that
59  * returns true if bounds a,b of the solution are close enough.
60  * We use relative criteria.
61  */
62 template <class T>
63 class tolerance
64 {
65 public:
66  tolerance(unsigned bits)
67  {
68  BOOST_MATH_STD_USING
69  eps = T(ldexp(1.0F, 1-bits));
70  }
71  bool operator()(const T& a, const T& b)
72  {
73  BOOST_MATH_STD_USING
74  return fabs(a - b) <= (eps * (std::max)(fabs(a), fabs(b)));
75  }
76 private:
77  T eps;
78 };
79 
80 
81 /**
82  * Complementary functor for isotherm of type "none"
83  */
84 class None {
85 public:
86  /// Constructor to set parameters
87  None(){}
88  /// Isotherm definition.
89  inline double operator()(double x) {
90  return (0.0);
91  }
92 };
93 
94 /**
95  * Functor for linear isotherm
96  */
97 class Linear {
98 public:
99  /// Constructor to set parameters
100  Linear(double mult_coef) : mult_coef_(mult_coef) {}
101  /// Isotherm definition.
102  inline double operator()(double x) {
103  return (mult_coef_*x);
104  }
105 private:
106  /// Parameters of the isotherm.
107  double mult_coef_;
108 };
109 
110 
111 
112 /**
113  * Functor for Langmuir isotherm.
114  */
115 class Langmuir {
116 public:
117  /// Constructor to set parameters
118  Langmuir( double mult_coef, double alpha) : mult_coef_(mult_coef), alpha_(alpha) {}
119  /// Isotherm definition.
120  inline double operator()( double x) {
121  return (mult_coef_*(alpha_ * x)/(alpha_ *x + 1));
122  }
123 
124 private:
125  /// Parameters of the isotherm.
126  double mult_coef_;
127  double alpha_;
128 };
129 
130 
131 
132 /**
133  * Functor for Freundlich isotherm.
134  */
135 class Freundlich {
136 public:
137  /// Constructor to set parameters
138  Freundlich(double mult_coef, double exponent) : mult_coef_(mult_coef), exponent_(exponent){}
139  /// Isotherm definition.
140  inline double operator()(double x) {
141  return (mult_coef_*pow(x, exponent_));
142  }
143 
144 private:
145  /// Parameters of the isotherm.
146  double mult_coef_;
147  double exponent_;
148 };
149 
150 
151 /**
152 * Class describing one isotherm with possibly precalculated interpolation table.
153 */
154 class Isotherm {
155 public:
156  TYPEDEF_ERR_INFO( EI_BoostMessage, std::string);
157  DECLARE_EXCEPTION( ExcBoostSolver, <<
158  "The Boost solver of nonlinear equation did not converge in sorption model.\n"
159  << "Check input at the following address for possible user error: \t" << Input::EI_Address::val << "\n"
160  "Solver message: \n" << EI_BoostMessage::val);
161 
162  TYPEDEF_ERR_INFO( EI_TotalMass, double);
163  DECLARE_EXCEPTION( ExcNegativeTotalMass, <<
164  "The total mass in sorption model became negative during the computation (value: "
165  << EI_TotalMass::val << ").\n"
166  << "Check input at the following address for possible user error: \t" << Input::EI_Address::val << "\n");
167 
168 
169  /// Type of adsorption isotherm.
171  none = 0,
172  linear = 1,
173  freundlich = 2,
174  langmuir = 3
175  };
176 
177  /// Pair of soluted and adsorbed concentration.
178  struct ConcPair {
179  ConcPair(double x, double y) : fluid(x), solid(y) {}
180  double fluid;
181  double solid;
182  };
183 
184  /**
185  * Setting adsorption parameters for general isotherm. These parameters are then used either
186  * for creation of the interpolation table via @p make_table method or just one adsorption is computed
187  * through @p compute method. Provided parameters are:
188  * @param sorption_type - type of isotherm
189  * @param limited_solubility_on - true if @p c_aqua_limit is solubility limit
190  * @param aqua_density - density of the liquid phase
191  * @param scale_aqua - generalized porosity, fraction of the space with liquid phase
192  * @param scale_sorbed - fraction of the space with the solid to which we adsorp
193  * @param c_aqua_limit - limit for interpolation table, possibly solubility limit
194  * @param mult_coef - multiplicative coefficient of the isotherm (all isotherms have one)
195  * @param second_coef - possibly second parameter of the isotherm
196  */
197  inline void reinit(enum SorptionType sorption_type, bool limited_solubility_on,
198  double aqua_density, double scale_aqua, double scale_sorbed,
199  double c_aqua_limit, double mult_coef, double second_coef);
200 
201  /**
202  * Create interpolation table for isotherm in rotated coordinate system with X axes given by total mass in
203  * both phases. Size of the table is the only parameter. Currently we support only linear interpolation.
204  * @p reinit has to be called just before this method.
205  */
206  void make_table(int n_points);
207 
208  /**
209  * Direct calculation of the equilibrium adsorption using a non-linear solver.
210  * @p reinit has to be called just before this method.
211  */
212  inline void compute(double &c_aqua, double &c_sorbed);
213 
214  /**
215  * Use interpolation to determine equilibrium state.
216  * Assumes previous call to @p make_table. If total mass is larger then table limit we either
217  * call @p precipitate (limit_solubility_on) or use direct computation.
218  */
219  inline void interpolate(double &c_aqua, double &c_sorbed);
220 
221  /**
222  * Returns true if interpolation table is created.
223  */
224  inline bool is_precomputed(void) {
225  return interpolation_table.size() != 0;
226  }
227 
228 protected:
229  /**
230  * Implementation of interpolation construction for particular isotherm functor.
231  * Specialized for sorption 'none' - creates empty table.
232  */
233  template<class Func>
234  void make_table(const Func &isotherm, int n_points);
235  /**
236  * Find new values for concentrations in @p c_pair that has same total mass and lies on the
237  * @p isotherm (functor object).
238  * Specialized for sorption 'none' - returns unchanged @p c_pair.
239  */
240  template<class Func>
241  inline ConcPair solve_conc(ConcPair c_pair, const Func &isotherm);
242  /**
243  * Dispatch isotherm type and use appropriate template.
244  */
245  inline ConcPair solve_conc(ConcPair conc);
246  /**
247  * Update concentrations using interopolation.
248  */
249  inline ConcPair compute_projection( ConcPair conc );
250  /**
251  * Modify concentrations after adsorption for limited solubility.
252  */
253  inline ConcPair precipitate( ConcPair conc );
254 
255  /****************************************
256  * Data
257  */
258 
259  /// Type of isotherm
260  enum SorptionType adsorption_type_;
261 
262  /// Multiplication parameter of the isotherm
263  double mult_coef_;
264 
265  /// Optional secod parameter of the isotherm
266  double second_coef_;
267 
268  /* Concentration in liquid phase for limit of the interpolation table, or
269  * solubility limit.
270  */
271  double table_limit_;
272 
273  /// Solubility limit flag
275 
276  /// density of the solvent
277  double rho_aqua_;
278  /// coefficient that convert soluted concentration to mass; porosity = k_W, originally rho_aqua*porosity = k_W
279  double scale_aqua_;
280  /// coefficient that convert adsorbed molar concentration to mass; molar_weight * rho_rock * (1 - porosity) = k_H
282  /// reciprocal values
283  double inv_scale_aqua_, inv_scale_sorbed_;
284  /**
285  * Interpolation table of isotherm in the rotated coordinates.
286  * The X axes of rotated system is total mass, the Y axes is perpendicular.
287  */
289  /**
290  * Step on the rotated X axes (total mass).
291  */
293 
294 };
295 
296 // Sorption none template specializations
297 template<> Isotherm::ConcPair Isotherm::solve_conc(Isotherm::ConcPair c_pair, const None &isotherm);
298 template<> void Isotherm::make_table(const None &isotherm, int n_steps);
299 
300 
301 /**
302  * Functor for solved equation in form F(x) ==0.
303  * Function @p func is an isotherm functor object in concentration based coordinated system.
304  * We solve the equation in modified system (scaled and rotated) for better numerical stability.
305  * The solved equation reads:
306  * F(X) -Y =0, where
307  * X is total mass , Y
308  */
309 template <class Func>
311 {
312 public:
313  CrossFunction(const Func &func_, double total_mass, double scale_aqua, double scale_sorbed, double rho_aqua)
314  : func(func_), total_mass_(total_mass),
315  scale_sorbed_(scale_sorbed), scale_aqua_(scale_aqua), rho_aqua_(rho_aqua)
316  {}
317 
318  double operator()( double conc_aqua)
319  {
320  // that is the selected isotherm // scale_sorbed_ * func( conc_aqua ) + scale_aqua_ * conc_aqua - total_mass_
321  return scale_sorbed_*func( conc_aqua/rho_aqua_) + (scale_aqua_) * conc_aqua - total_mass_;
322  }
323 private:
324  Func func;
325  double total_mass_, scale_sorbed_, scale_aqua_, rho_aqua_;
326 };
327 
328 
329 /*****************************************************************************************************************
330  * IMPLEMENTATION
331  */
332 
333 
334 inline void Isotherm::reinit(enum SorptionType adsorption_type, bool limited_solubility_on,
335  double rho_aqua, double scale_aqua, double scale_sorbed,
336  double c_aqua_limit, double mult_coef, double second_coef)
337 {
338  adsorption_type_ = adsorption_type;
339  rho_aqua_ = rho_aqua;
340  scale_aqua_ = scale_aqua;
341  scale_sorbed_ = scale_sorbed;
342  limited_solubility_on_ = limited_solubility_on;
343  mult_coef_ = mult_coef*rho_aqua;
344  second_coef_ = second_coef;
345 
346  table_limit_ = c_aqua_limit;
347  inv_scale_aqua_ = scale_aqua_/(scale_aqua_*scale_aqua_ + scale_sorbed_*scale_sorbed_);
348  inv_scale_sorbed_ = scale_sorbed_/(scale_aqua_*scale_aqua_ + scale_sorbed_*scale_sorbed_);
349 }
350 
351 
352 inline void Isotherm::compute( double &c_aqua, double &c_sorbed ) {
353  ConcPair c_pair(c_aqua, c_sorbed);
354  ConcPair result(0,0);
355 
356  if (limited_solubility_on_ && (c_pair.fluid > table_limit_)) {
357  result = precipitate( c_pair );
358  } else {
359  result = solve_conc( c_pair );
360  }
361  c_aqua=result.fluid;
362  c_sorbed=result.solid;
363 }
364 
365 
366 inline void Isotherm::interpolate( double &c_aqua, double &c_sorbed ) {
367  ConcPair c_pair(c_aqua, c_sorbed);
368  ConcPair result(0,0);
369 
370  result = compute_projection( c_pair );
371 
372  c_aqua=result.fluid;
373  c_sorbed=result.solid;
374 }
375 
376 
377 
379  double total_mass = (scale_aqua_* c_pair.fluid + scale_sorbed_ * c_pair.solid);
380  if(total_mass < 0.0)
381  THROW( Isotherm::ExcNegativeTotalMass()
382  << EI_TotalMass(total_mass)
383  );
384  // total_mass_step_ is set and checked in make_table
385  double total_mass_steps = total_mass / total_mass_step_;
386  unsigned int total_mass_idx = static_cast <unsigned int>(std::floor(total_mass_steps));
387 
388  if (total_mass_idx < (interpolation_table.size() - 1) ) {
389  double rot_sorbed = interpolation_table[total_mass_idx] + (total_mass_steps - total_mass_idx)*(interpolation_table[total_mass_idx+1] - interpolation_table[total_mass_idx]);
390  return ConcPair( (total_mass * inv_scale_aqua_ - rot_sorbed * inv_scale_sorbed_),
391  (total_mass * inv_scale_sorbed_ + rot_sorbed * inv_scale_aqua_) );
392  } else {
393  if (limited_solubility_on_) {
394  return precipitate( c_pair );
395  } else {
396  return solve_conc( c_pair );
397  }
398  }
399 }
400 
401 
402 
404  double total_mass = (scale_aqua_*c_pair.fluid + scale_sorbed_ * c_pair.solid);
405  return ConcPair( table_limit_,
406  (total_mass - scale_aqua_ * table_limit_)/scale_sorbed_ );
407 }
408 
409 
410 
411 template<class Func>
412 inline Isotherm::ConcPair Isotherm::solve_conc(Isotherm::ConcPair c_pair, const Func &isotherm)
413 {
414  boost::uintmax_t max_iter = 20;
415  tolerance<double> toler(30);
416  double total_mass = (scale_aqua_*c_pair.fluid + scale_sorbed_ * c_pair.solid);
417  double mass_limit = table_limit_*scale_aqua_ + const_cast<Func &>(isotherm)(table_limit_ / this->rho_aqua_)*scale_sorbed_;
418  double upper_solution_bound;
419 
420  if(total_mass >= mass_limit)
421  {
422  mass_limit = total_mass;
423  }
424  upper_solution_bound = mass_limit / scale_aqua_;
425  CrossFunction<Func> eq_func(isotherm, total_mass, scale_aqua_, scale_sorbed_, this->rho_aqua_);
426  pair<double,double> solution;
427  if (total_mass > 0) // here should be probably some kind of tolerance instead of "0"
428  {
429  try {
430  solution = boost::math::tools::toms748_solve(eq_func, 0.0, upper_solution_bound, toler, max_iter);
431  }
432  catch(boost::exception const & e)
433  {
434  THROW( Isotherm::ExcBoostSolver()
435  << EI_BoostMessage(boost::diagnostic_information(e))
436  );
437  }
438  }
439  double difference;
440  difference = (solution.second - solution.first)/2;
441  double c_aqua = solution.first + difference;
442  return ConcPair( c_aqua,
443  (total_mass - scale_aqua_ * c_aqua)/scale_sorbed_);
444 }
445 
446 
448 {
449  switch(adsorption_type_)
450  {
451  case 0: // none
452  {
453  None obj_isotherm;
454  return solve_conc( conc, obj_isotherm);
455  }
456  break;
457  case 1: // linear:
458  {
459  Linear obj_isotherm(mult_coef_);
460  return solve_conc( conc, obj_isotherm);
461  }
462  break;
463  case 2: // freundlich
464  {
465  Freundlich obj_isotherm(mult_coef_, second_coef_);
466  return solve_conc( conc, obj_isotherm);
467  }
468  break;
469  case 3: // langmuir:
470  {
471  Langmuir obj_isotherm(mult_coef_, second_coef_);
472  return solve_conc( conc, obj_isotherm);
473  }
474  break;
475  }
476  return conc;
477 }
478 
479 
480 template<class Func>
481 void Isotherm::make_table(const Func &isotherm, int n_steps)
482 {
483  double mass_limit = scale_aqua_ * table_limit_ + scale_sorbed_ * const_cast<Func &>(isotherm)(table_limit_ / this->rho_aqua_);
484 
485  if(mass_limit < 0.0)
486  THROW( Isotherm::ExcNegativeTotalMass()
487  << EI_TotalMass(mass_limit)
488  );
489  total_mass_step_ = mass_limit / n_steps;
490  double mass = 0.0;
491  for(int i=0; i<= n_steps; i++) {
492  // aqueous concentration (original coordinates c_a) corresponding to i-th total_mass_step_
493  ConcPair c_pair( mass/scale_aqua_, 0.0 );
494 
495  ConcPair result = solve_conc( c_pair, isotherm);
496  double c_sorbed_rot = ( result.solid * scale_aqua_ - result.fluid * scale_sorbed_);
497  interpolation_table.push_back(c_sorbed_rot);
498  mass = mass+total_mass_step_;
499  }
500 
501  return;
502 }
503 
504 #endif /* SORPTION_IMPL_HH_ */
double operator()(double x)
Isotherm definition.
Definition: isotherm.hh:89
tolerance(unsigned bits)
Definition: isotherm.hh:66
double mult_coef_
Multiplication parameter of the isotherm.
Definition: isotherm.hh:263
CrossFunction(const Func &func_, double total_mass, double scale_aqua, double scale_sorbed, double rho_aqua)
Definition: isotherm.hh:313
void compute(double &c_aqua, double &c_sorbed)
Definition: isotherm.hh:352
double second_coef_
Optional secod parameter of the isotherm.
Definition: isotherm.hh:266
IntersectionPoint * interpolate(const IntersectionPoint &A1, const IntersectionPoint &A2, double t)
ConcPair(double x, double y)
Definition: isotherm.hh:179
double operator()(double x)
Isotherm definition.
Definition: isotherm.hh:140
Definition: isotherm.hh:84
ConcPair solve_conc(ConcPair c_pair, const Func &isotherm)
Definition: isotherm.hh:412
Langmuir(double mult_coef, double alpha)
Constructor to set parameters.
Definition: isotherm.hh:118
double total_mass_step_
Definition: isotherm.hh:292
#define DECLARE_EXCEPTION(ExcName, Format)
Macro for simple definition of exceptions.
Definition: exceptions.hh:158
double alpha_
Definition: isotherm.hh:127
double inv_scale_sorbed_
Definition: isotherm.hh:283
double operator()(double x)
Isotherm definition.
Definition: isotherm.hh:102
double mult_coef_
Parameters of the isotherm.
Definition: isotherm.hh:126
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.hh:334
Freundlich(double mult_coef, double exponent)
Constructor to set parameters.
Definition: isotherm.hh:138
double mult_coef_
Parameters of the isotherm.
Definition: isotherm.hh:107
vector< double > interpolation_table
Definition: isotherm.hh:288
double scale_sorbed_
coefficient that convert adsorbed molar concentration to mass; molar_weight * rho_rock * (1 - porosit...
Definition: isotherm.hh:281
ConcPair precipitate(ConcPair conc)
Definition: isotherm.hh:403
Linear(double mult_coef)
Constructor to set parameters.
Definition: isotherm.hh:100
#define TYPEDEF_ERR_INFO(EI_Type, Type)
Macro to simplify declaration of error_info types.
Definition: exceptions.hh:194
bool limited_solubility_on_
Solubility limit flag.
Definition: isotherm.hh:274
Pair of soluted and adsorbed concentration.
Definition: isotherm.hh:178
double operator()(double conc_aqua)
Definition: isotherm.hh:318
double mult_coef_
Parameters of the isotherm.
Definition: isotherm.hh:146
SorptionType
Type of adsorption isotherm.
Definition: isotherm.hh:170
double exponent_
Definition: isotherm.hh:147
bool operator()(const T &a, const T &b)
Definition: isotherm.hh:71
void make_table(int n_points)
Definition: isotherm.cc:22
None()
Constructor to set parameters.
Definition: isotherm.hh:87
void interpolate(double &c_aqua, double &c_sorbed)
Definition: isotherm.hh:366
double total_mass_
Definition: isotherm.hh:325
ConcPair compute_projection(ConcPair conc)
Definition: isotherm.hh:378
double scale_aqua_
coefficient that convert soluted concentration to mass; porosity = k_W, originally rho_aqua*porosity ...
Definition: isotherm.hh:279
double table_limit_
Definition: isotherm.hh:271
bool is_precomputed(void)
Definition: isotherm.hh:224
double operator()(double x)
Isotherm definition.
Definition: isotherm.hh:120
double rho_aqua_
density of the solvent
Definition: isotherm.hh:277
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53