Flow123d  master-f44eb46
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 ISOTHERM_IMPL_HH_
48 #define ISOTHERM_IMPL_HH_
49 
50 #include <stdint.h> // for uintmax_t
51 #include <algorithm> // for max
52 #include <vector>
53 #include <cmath> // for floor, pow
54 #include <complex> // for fabs
55 #include <ostream> // for operator<<
56 #include <string> // for string
57 #include "input/input_exception.hh" // for EI_Address
58 #include "system/exceptions.hh" // for operator<<
59 #include <boost/math/special_functions/detail/round_fwd.hpp> // for BOOST_M...
60 
61 /**
62  * Convergence criteria for interval based nonlinear solver. It is functor, that
63  * returns true if bounds a,b of the solution are close enough.
64  * We use relative criteria.
65  */
66 template <class T>
67 class tolerance
68 {
69 public:
70  tolerance(unsigned bits)
71  {
72  BOOST_MATH_STD_USING
73  eps = T(ldexp(1.0F, 1-bits));
74  }
75  bool operator()(const T& a, const T& b)
76  {
77  BOOST_MATH_STD_USING
78  return fabs(a - b) <= (eps * (std::max)(fabs(a), fabs(b)));
79  }
80 private:
81  T eps;
82 };
83 
84 
85 /**
86  * Complementary functor for isotherm of type "none"
87  */
88 class None {
89 public:
90  /// Constructor to set parameters
91  None(){}
92  /// Isotherm definition.
93  inline double operator()(double) {
94  return (0.0);
95  }
96 };
97 
98 /**
99  * Functor for linear isotherm
100  */
101 class Linear {
102 public:
103  /// Constructor to set parameters
104  Linear(double mult_coef) : mult_coef_(mult_coef) {}
105  /// Isotherm definition.
106  inline double operator()(double x) {
107  return (mult_coef_*x);
108  }
109 private:
110  /// Parameters of the isotherm.
111  double mult_coef_;
112 };
113 
114 
115 
116 /**
117  * Functor for Langmuir isotherm.
118  */
119 class Langmuir {
120 public:
121  /// Constructor to set parameters
122  Langmuir( double mult_coef, double alpha) : mult_coef_(mult_coef), alpha_(alpha) {}
123  /// Isotherm definition.
124  inline double operator()( double x) {
125  return (mult_coef_*(alpha_ * x)/(alpha_ *x + 1));
126  }
127 
128 private:
129  /// Parameters of the isotherm.
130  double mult_coef_;
131  double alpha_;
132 };
133 
134 
135 
136 /**
137  * Functor for Freundlich isotherm.
138  */
139 class Freundlich {
140 public:
141  /// Constructor to set parameters
142  Freundlich(double mult_coef, double exponent) : mult_coef_(mult_coef), exponent_(exponent){}
143  /// Isotherm definition.
144  inline double operator()(double x) {
145  return (mult_coef_*pow(x, exponent_));
146  }
147 
148 private:
149  /// Parameters of the isotherm.
150  double mult_coef_;
151  double exponent_;
152 };
153 
154 
155 /**
156 * Class describing one isotherm with possibly precalculated interpolation table.
157 */
158 class Isotherm {
159 public:
160  TYPEDEF_ERR_INFO( EI_BoostMessage, std::string);
161  DECLARE_EXCEPTION( ExcBoostSolver, <<
162  "The Boost solver of nonlinear equation did not converge in sorption model.\n"
163  << "Check input at the following address for possible user error: \t" << Input::EI_Address::val << "\n"
164  "Solver message: \n" << EI_BoostMessage::val);
165 
166  TYPEDEF_ERR_INFO( EI_TotalMass, double);
167  DECLARE_EXCEPTION( ExcNegativeTotalMass, <<
168  "The total mass in sorption model became negative during the computation (value: "
169  << EI_TotalMass::val << ").\n"
170  << "Check input at the following address for possible user error: \t" << Input::EI_Address::val << "\n");
171 
172 
173  /// Type of adsorption isotherm.
175  none = 0,
176  linear = 1,
178  langmuir = 3
179  };
180 
181  /// Pair of soluted and adsorbed concentration.
182  struct ConcPair {
183  ConcPair(double x, double y) : fluid(x), solid(y) {}
184  double fluid;
185  double solid;
186  };
187 
188  /// Default constructor.
189  Isotherm();
190 
191  /**
192  * Setting adsorption parameters for general isotherm. These parameters are then used either
193  * for creation of the interpolation table via @p make_table method or just one adsorption is computed
194  * through @p compute method. Provided parameters are:
195  * @param sorption_type - type of isotherm
196  * @param limited_solubility_on - true if @p c_aqua_limit is solubility limit
197  * @param aqua_density - density of the liquid phase
198  * @param scale_aqua - generalized porosity, fraction of the space with liquid phase
199  * @param scale_sorbed - fraction of the space with the solid to which we adsorp
200  * @param c_aqua_limit - solubility limit (limit aqueous concentration)
201  * @param mult_coef - multiplicative coefficient of the isotherm (all isotherms have one)
202  * @param second_coef - possibly second parameter of the isotherm
203  */
204  void reinit(enum SorptionType sorption_type, bool limited_solubility_on,
205  double aqua_density, double scale_aqua, double scale_sorbed,
206  double c_aqua_limit, double mult_coef, double second_coef);
207 
208  /**
209  * Create interpolation table for isotherm in rotated coordinate system with X axes given by total mass in
210  * both phases. Currently we support only linear interpolation.
211  * @p reinit has to be called just before this method.
212  * @param n_points is the size of the table
213  * @param table_limit is the limit value of aqueous concentration
214  */
215  void make_table(unsigned int n_points, double table_limit);
216 
217  /**
218  * Clears the interpolation table and resets the table limit.
219  * (That means, no interpolation takes place in the computation.)
220  */
221  void clear_table();
222 
223  /**
224  * Direct calculation of the equilibrium adsorption using a non-linear solver.
225  * @p reinit has to be called just before this method.
226  */
227  void compute(double &c_aqua, double &c_sorbed);
228 
229  /**
230  * Use interpolation to determine equilibrium state.
231  * Assumes previous call to @p make_table. If total mass is larger then table limit we either
232  * call @p precipitate (limit_solubility_on) or use direct computation.
233  */
234  void interpolate(double &c_aqua, double &c_sorbed);
235 
236  /**
237  * Returns true if interpolation table is created.
238  */
239  inline bool is_precomputed(void) {
240  return interpolation_table.size() != 0;
241  }
242 
243  /// Getter for table limit (limit aqueous concentration).
244  inline double table_limit(void) const
245  { return table_limit_;}
246 
247 protected:
248  /**
249  * Implementation of interpolation construction for particular isotherm functor.
250  * Specialized for sorption 'none' - creates empty table.
251  */
252  template<class Func>
253  void make_table(const Func &isotherm, int n_points);
254  /**
255  * Find new values for concentrations in @p c_pair that has same total mass and lies on the
256  * @p isotherm (functor object).
257  * Specialized for sorption 'none' - returns unchanged @p c_pair.
258  */
259  template<class Func>
260  ConcPair solve_conc(ConcPair c_pair, const Func &isotherm);
261  /**
262  * Dispatch isotherm type and use appropriate template.
263  */
264  ConcPair solve_conc(ConcPair conc);
265  /**
266  * Update concentrations using interopolation.
267  */
268  ConcPair compute_projection( ConcPair conc );
269  /**
270  * Modify concentrations after adsorption for limited solubility.
271  */
272  ConcPair precipitate( ConcPair conc );
273 
274  double get_total_mass( ConcPair conc );
275 
276  /****************************************
277  * Data
278  */
279 
280  /// Type of isotherm
282 
283  /// Multiplication parameter of the isotherm
284  double mult_coef_;
285 
286  /// Optional secod parameter of the isotherm
287  double second_coef_;
288 
289  /// Concentration in liquid phase for limit of the interpolation table.
290  double table_limit_;
291 
292  /// Solubility limit flag
294  /// Concentration limit in liquid phase (solubility limit).
296 
297  /// density of the solvent
298  double rho_aqua_;
299  /// coefficient that convert soluted concentration to mass; porosity = k_W, originally rho_aqua*porosity = k_W
300  double scale_aqua_;
301  /// coefficient that convert adsorbed molar concentration to mass; molar_weight * rho_rock * (1 - porosity) = k_H
303  /// reciprocal values
305  /**
306  * Interpolation table of isotherm in the rotated coordinates.
307  * The X axes of rotated system is total mass, the Y axes is perpendicular.
308  */
310  /**
311  * Step on the rotated X axes (total mass).
312  */
314 
315 };
316 
317 
318 /**
319  * Functor for solved equation in form F(x) ==0.
320  * Function @p func is an isotherm functor object in concentration based coordinated system.
321  * We solve the equation in modified system (scaled and rotated) for better numerical stability.
322  * The solved equation reads:
323  * F(X) -Y =0, where
324  * X is total mass , Y
325  */
326 template <class Func>
328 {
329 public:
330  CrossFunction(const Func &func_, double total_mass, double scale_aqua, double scale_sorbed, double rho_aqua)
331  : func(func_), total_mass_(total_mass),
332  scale_sorbed_(scale_sorbed), scale_aqua_(scale_aqua), rho_aqua_(rho_aqua)
333  {}
334 
335  double operator()( double conc_aqua)
336  {
337  // that is the selected isotherm // scale_sorbed_ * func( conc_aqua ) + scale_aqua_ * conc_aqua - total_mass_
338  return scale_sorbed_*func( conc_aqua/rho_aqua_) + (scale_aqua_) * conc_aqua - total_mass_;
339  }
340 private:
341  Func func;
343 };
344 
345 #endif /* ISOTHERM_IMPL_HH_ */
double rho_aqua_
Definition: isotherm.hh:342
double scale_sorbed_
Definition: isotherm.hh:342
double scale_aqua_
Definition: isotherm.hh:342
CrossFunction(const Func &func_, double total_mass, double scale_aqua, double scale_sorbed, double rho_aqua)
Definition: isotherm.hh:330
double total_mass_
Definition: isotherm.hh:342
double operator()(double conc_aqua)
Definition: isotherm.hh:335
double exponent_
Definition: isotherm.hh:151
double mult_coef_
Parameters of the isotherm.
Definition: isotherm.hh:150
Freundlich(double mult_coef, double exponent)
Constructor to set parameters.
Definition: isotherm.hh:142
double operator()(double x)
Isotherm definition.
Definition: isotherm.hh:144
DECLARE_EXCEPTION(ExcBoostSolver,<< "The Boost solver of nonlinear equation did not converge in sorption model.\n"<< "Check input at the following address for possible user error: \t"<< Input::EI_Address::val<< "\n" "Solver message: \n"<< EI_BoostMessage::val)
TYPEDEF_ERR_INFO(EI_BoostMessage, std::string)
double solubility_limit_
Concentration limit in liquid phase (solubility limit).
Definition: isotherm.hh:295
ConcPair precipitate(ConcPair conc)
Definition: isotherm.cc:127
bool is_precomputed(void)
Definition: isotherm.hh:239
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
TYPEDEF_ERR_INFO(EI_TotalMass, double)
ConcPair solve_conc(ConcPair c_pair, const Func &isotherm)
Definition: isotherm.cc:137
DECLARE_EXCEPTION(ExcNegativeTotalMass,<< "The total mass in sorption model became negative during the computation (value: "<< EI_TotalMass::val<< ").\n"<< "Check input at the following address for possible user error: \t"<< Input::EI_Address::val<< "\n")
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
@ freundlich
Definition: isotherm.hh:177
@ langmuir
Definition: isotherm.hh:178
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
double mult_coef_
Parameters of the isotherm.
Definition: isotherm.hh:130
double operator()(double x)
Isotherm definition.
Definition: isotherm.hh:124
double alpha_
Definition: isotherm.hh:131
Langmuir(double mult_coef, double alpha)
Constructor to set parameters.
Definition: isotherm.hh:122
double mult_coef_
Parameters of the isotherm.
Definition: isotherm.hh:111
Linear(double mult_coef)
Constructor to set parameters.
Definition: isotherm.hh:104
double operator()(double x)
Isotherm definition.
Definition: isotherm.hh:106
Definition: isotherm.hh:88
None()
Constructor to set parameters.
Definition: isotherm.hh:91
double operator()(double)
Isotherm definition.
Definition: isotherm.hh:93
tolerance(unsigned bits)
Definition: isotherm.hh:70
bool operator()(const T &a, const T &b)
Definition: isotherm.hh:75
Pair of soluted and adsorbed concentration.
Definition: isotherm.hh:182
ConcPair(double x, double y)
Definition: isotherm.hh:183