Flow123d  release_2.2.0-26-ge868538
sorption.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 sorption.cc
15  * @brief
16  */
17 
18 #include <vector>
19 #include <limits>
20 
21 #include "reaction/isotherm.hh"
22 #include "reaction/sorption.hh"
23 #include "system/sys_profiler.hh"
24 #include "mesh/accessors.hh"
25 #include "input/factory.hh"
26 
27 FLOW123D_FORCE_LINK_IN_CHILD(sorptionMobile)
28 FLOW123D_FORCE_LINK_IN_CHILD(sorptionImmobile)
30 
31 
32 /********************************* SORPTION_SIMPLE *********************************************************/
33 /********************************* *********************************************************/
34 
35 const IT::Record & SorptionSimple::get_input_type() {
36  return IT::Record("Sorption", "Sorption model in the reaction term of transport.")
39  //.declare_key("output_fields", IT::Array(make_output_selection("conc_solid", "Sorption_Output")),
40  // IT::Default("\"conc_solid\""), "List of fields to write to output stream.")
41  .declare_key("output", make_output_type("Sorption", "conc_solid"),
42  IT::Default("{ \"fields\": [ \"conc_solid\" ] }"),
43  "Setting of the fields output.")
44 
45  .close();
46 }
47 
49  : SorptionBase(init_mesh, in_rec)
50 {
51  data_ = new EqData("conc_solid");
52  this->eq_data_ = data_;
53  //output_selection = make_output_selection(
54  // "SorptionSimple_output_fields",
55  // "Selection of field names of Simple Sorption model available for output.");
56 }
57 
58 const int SorptionSimple::registrar =
59  Input::register_class< SorptionSimple, Mesh &, Input::Record >("Sorption") +
61 
63 {}
64 
66 {
67  START_TIMER("SorptionSimple::isotherm_reinit");
68 
69  double rock_density = data_->rock_density.value(elem.centre(),elem);
70  double por_m = data_->porosity.value(elem.centre(),elem);
71 
72  for(unsigned int i_subst = 0; i_subst < n_substances_; i_subst++)
73  {
74  double mult_coef = data_->distribution_coefficient[i_subst].value(elem.centre(),elem);
75  double second_coef = data_->isotherm_other[i_subst].value(elem.centre(),elem);
76  Isotherm & isotherm = isotherms_vec[i_subst];
77 
78  //scales are different for the case of sorption in mobile and immobile pores
79  double scale_aqua = por_m,
80  scale_sorbed = (1 - por_m) * rock_density;
81 
82  bool limited_solubility_on = false;
83  double table_limit;
84  if (solubility_vec_[i_subst] <= 0.0) {
85  limited_solubility_on = false;
86  table_limit=table_limit_[i_subst];
87  } else {
88  limited_solubility_on = true;
89  table_limit=solubility_vec_[i_subst];
90  }
91 
92  if( (1-por_m) <= std::numeric_limits<double>::epsilon()) //means there is no sorbing surface
93  {
94  isotherm.reinit(Isotherm::none, false, solvent_density_, scale_aqua, scale_sorbed,table_limit,0,0);
95  continue;
96  }
97 
98  if ( scale_sorbed <= 0.0)
99  xprintf(UsrErr, "Scaling parameter in sorption is not positive. Check the input for rock density and molar mass of %d. substance.",i_subst);
100 
101  isotherm.reinit(Isotherm::SorptionType(data_->sorption_type[i_subst].value(elem.centre(),elem)), limited_solubility_on,
102  solvent_density_, scale_aqua, scale_sorbed, table_limit, mult_coef, second_coef);
103 
104  }
105 
106  END_TIMER("SorptionSimple::isotherm_reinit");
107 }
108 
109 
110 /*********************************** *********************************************************/
111 /*********************************** SORPTION_DUAL *********************************************************/
112 /*********************************** *********************************************************/
113 
115  const string &output_conc_name,
116  const string &output_selection_name)
117  : SorptionBase(init_mesh, in_rec)
118 {
119  data_ = new EqData(output_conc_name);
122  .name("porosity_immobile")
123  .set_limits(0.0);
124  this->eq_data_ = data_;
125  //output_selection = make_output_selection(output_conc_name, output_selection_name);
126 }
127 
129 {}
130 
131 /********************************** *******************************************************/
132 /*********************************** SORPTION_MOBILE *******************************************************/
133 /********************************** *******************************************************/
134 
136  return IT::Record("SorptionMobile", "Sorption model in the mobile zone, following the dual porosity model.")
139  //.declare_key("output_fields", IT::Array(make_output_selection("conc_solid", "SorptionMobile_Output")),
140  // IT::Default("\"conc_solid\""), "List of fields to write to output stream.")
141  .declare_key("output", make_output_type("SorptionMobile", "conc_solid"),
142  IT::Default("{ \"fields\": [ \"conc_solid\" ] }"),
143  "Setting of the fields output.")
144 
145  .close();
146 }
147 
148 
149 const int SorptionMob::registrar =
150  Input::register_class< SorptionMob, Mesh &, Input::Record >("SorptionMobile") +
152 
153 
155  : SorptionDual(init_mesh, in_rec, "conc_solid", "SorptionMobile_Output")
156 {}
157 
158 
160 {}
161 
162 /*
163 double SorptionMob::compute_sorbing_scale(double por_m, double por_imm)
164 {
165  double phi = por_m/(por_m + por_imm);
166  return phi;
167 }
168 */
169 
171 {
172  START_TIMER("SorptionMob::isotherm_reinit");
173 
174  double rock_density = data_->rock_density.value(elem.centre(),elem);
175 
176  double por_m = data_->porosity.value(elem.centre(),elem);
177  double por_imm = immob_porosity_.value(elem.centre(),elem);
178  double phi = por_m/(por_m + por_imm);
179 
180  for(unsigned int i_subst = 0; i_subst < n_substances_; i_subst++)
181  {
182  double mult_coef = data_->distribution_coefficient[i_subst].value(elem.centre(),elem);
183  double second_coef = data_->isotherm_other[i_subst].value(elem.centre(),elem);
184  Isotherm & isotherm = isotherms_vec[i_subst];
185 
186  //scales are different for the case of sorption in mobile and immobile pores
187  double scale_aqua = por_m,
188  scale_sorbed = phi * (1 - por_m - por_imm) * rock_density;
189 
190  bool limited_solubility_on;
191  double table_limit;
192  if (solubility_vec_[i_subst] <= 0.0) {
193  limited_solubility_on = false;
194  table_limit=table_limit_[i_subst];
195 
196  } else {
197  limited_solubility_on = true;
198  table_limit=solubility_vec_[i_subst];
199  }
200 
201  if( (1-por_m-por_imm) <= std::numeric_limits<double>::epsilon()) //means there is no sorbing surface
202  {
203  isotherm.reinit(Isotherm::none, false, solvent_density_, scale_aqua, scale_sorbed,table_limit,0,0);
204  continue;
205  }
206 
207  if ( scale_sorbed <= 0.0)
208  xprintf(UsrErr, "Scaling parameter in sorption is not positive. Check the input for rock density and molar mass of %d. substance.",i_subst);
209 
210  isotherm.reinit(Isotherm::SorptionType(data_->sorption_type[i_subst].value(elem.centre(),elem)), limited_solubility_on,
211  solvent_density_, scale_aqua, scale_sorbed, table_limit, mult_coef, second_coef);
212 
213  }
214 
215  END_TIMER("SorptionMob::isotherm_reinit");
216 }
217 
218 
219 /*********************************** *****************************************************/
220 /*********************************** SORPTION_IMMOBILE *****************************************************/
221 /*********************************** *****************************************************/
222 
224  return IT::Record("SorptionImmobile", "Sorption model in the immobile zone, following the dual porosity model.")
227  //.declare_key("output_fields", IT::Array(make_output_selection("conc_immobile_solid", "SorptionImmobile_Output")),
228  // IT::Default("\"conc_immobile_solid\""), "List of fields to write to output stream.")
229  .declare_key("output", make_output_type("SorptionImmobile", "conc_immobile_solid"),
230  IT::Default("{ \"fields\": [ \"conc_immobile_solid\" ] }"),
231  "Setting of the fields output.")
232 
233  .close();
234 }
235 
236 const int SorptionImmob::registrar =
237  Input::register_class< SorptionImmob, Mesh &, Input::Record >("SorptionImmobile") +
239 
241 : SorptionDual(init_mesh, in_rec, "conc_immobile_solid", "SorptionImmobile_Output")
242 {}
243 
245 {}
246 
247 /*
248 double SorptionImmob::compute_sorbing_scale(double por_m, double por_imm)
249 {
250  double phi = por_imm / (por_m + por_imm);
251  return phi;
252 }
253 */
254 
256 {
257  START_TIMER("SorptionImmob::isotherm_reinit");
258 
259  double rock_density = data_->rock_density.value(elem.centre(),elem);
260 
261  double por_m = data_->porosity.value(elem.centre(),elem);
262  double por_imm = immob_porosity_.value(elem.centre(),elem);
263  double phi = por_m/(por_m + por_imm);
264 
265  for(unsigned int i_subst = 0; i_subst < n_substances_; i_subst++)
266  {
267  double mult_coef = data_->distribution_coefficient[i_subst].value(elem.centre(),elem);
268  double second_coef = data_->isotherm_other[i_subst].value(elem.centre(),elem);
269  Isotherm & isotherm = isotherms_vec[i_subst];
270 
271  //scales are different for the case of sorption in mobile and immobile pores
272  double scale_aqua = por_imm,
273  scale_sorbed = (1 - phi) * (1 - por_m - por_imm) * rock_density;
274 
275  bool limited_solubility_on;
276  double table_limit;
277  if (solubility_vec_[i_subst] <= 0.0) {
278  limited_solubility_on = false;
279  table_limit=table_limit_[i_subst];
280 
281  } else {
282  limited_solubility_on = true;
283  table_limit=solubility_vec_[i_subst];
284  }
285 
286  if( (1-por_m-por_imm) <= std::numeric_limits<double>::epsilon()) //means there is no sorbing surface
287  {
288  isotherm.reinit(Isotherm::none, false, solvent_density_, scale_aqua, scale_sorbed,0,0,0);
289  continue;
290  }
291 
292  if ( scale_sorbed <= 0.0)
293  xprintf(UsrErr, "Scaling parameter in sorption is not positive. Check the input for rock density and molar mass of %d. substance.",i_subst);
294 
295  isotherm.reinit(Isotherm::SorptionType(data_->sorption_type[i_subst].value(elem.centre(),elem)), limited_solubility_on,
296  solvent_density_, scale_aqua, scale_sorbed, table_limit, mult_coef, second_coef);
297 
298  }
299 
300  END_TIMER("SorptionImmob::isotherm_reinit");
301 }
FieldSet * eq_data_
Definition: equation.hh:232
double solvent_density_
~SorptionDual(void)
Destructor.
Definition: sorption.cc:128
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:598
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
FieldCommon & flags_add(FieldFlag::Flags::Mask mask)
~SorptionMob(void)
Destructor.
Definition: sorption.cc:159
Abstract class of sorption model in case dual porosity is considered.
Definition: sorption.hh:65
Field< 3, FieldValue< 3 >::Scalar > rock_density
Rock matrix density.
void isotherm_reinit(std::vector< Isotherm > &isotherms_vec, const ElementAccessor< 3 > &elem) override
Reinitializes the isotherm.
Definition: sorption.cc:170
Definition: mesh.h:97
static const Input::Type::Record & get_input_type()
Definition: sorption.cc:35
static Input::Type::Abstract & it_abstract_mobile_term()
std::vector< double > table_limit_
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
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:303
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:195
MultiField< 3, FieldValue< 3 >::Scalar > isotherm_other
Langmuir sorption coeficients alpha (in fraction c_s = omega * (alpha*c_a)/(1- alpha*c_a)).
Field< 3, FieldValue< 3 >::Scalar > porosity
Porosity field copied from transport.
Field< 3, FieldValue< 3 >::Scalar > immob_porosity_
Definition: sorption.hh:86
static constexpr Mask input_copy
Definition: field_flag.hh:44
std::vector< double > solubility_vec_
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
#define xprintf(...)
Definition: system.hh:92
#define START_TIMER(tag)
Starts a timer with specified tag.
static const Input::Type::Record & get_input_type()
Definition: sorption.cc:135
static Input::Type::Abstract & it_abstract_immobile_term()
virtual Value::return_type const & value(const Point &p, const ElementAccessor< spacedim > &elm) const
Definition: field.hh:362
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:490
MultiField< 3, FieldValue< 3 >::Enum > sorption_type
Discrete need Selection for initialization.
static Input::Type::Instance make_output_type(const string &equation_name, const string &output_field_name)
SorptionType
Type of adsorption isotherm.
Definition: isotherm.hh:170
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:215
static Input::Type::Abstract & it_abstract_term()
Definition: system.hh:64
const double epsilon
Definition: mathfce.h:23
SorptionImmob(Mesh &init_mesh, Input::Record in_rec)
Constructor.
Definition: sorption.cc:240
static const int registrar
Registrar of class to factory.
Definition: sorption.hh:58
void isotherm_reinit(std::vector< Isotherm > &isotherms, const ElementAccessor< 3 > &elm) override
Reinitializes the isotherm.
Definition: sorption.cc:65
SorptionSimple(Mesh &init_mesh, Input::Record in_rec)
Constructor.
Definition: sorption.cc:48
static const int registrar
Registrar of class to factory.
Definition: sorption.hh:144
~SorptionSimple(void)
Destructor.
Definition: sorption.cc:62
FieldCommon & name(const string &name)
Definition: field_common.hh:97
static const int registrar
Registrar of class to factory.
Definition: sorption.hh:116
#define END_TIMER(tag)
Ends a timer with specified tag.
arma::vec::fixed< spacedim > centre() const
Definition: accessors.hh:91
Simple sorption model without dual porosity.
Definition: sorption.hh:39
Record type proxy class.
Definition: type_record.hh:182
static const Input::Type::Record & get_input_type()
FieldCommon & set_limits(double min, double max=std::numeric_limits< double >::max())
MultiField< 3, FieldValue< 3 >::Scalar > distribution_coefficient
Multiplication coefficients (k, omega) for all types of isotherms.
unsigned int n_substances_
SorptionDual(Mesh &init_mesh, Input::Record in_rec, const string &output_conc_name, const string &output_selection_name)
Constructor.
Definition: sorption.cc:114
EqData * data_
Pointer to equation data. The object is constructed in descendants.
void isotherm_reinit(std::vector< Isotherm > &isotherms_vec, const ElementAccessor< 3 > &elem) override
Reinitializes the isotherm.
Definition: sorption.cc:255
Other possible transformation of coordinates:
static const Input::Type::Record & get_input_type()
Definition: sorption.cc:223
This file contains classes representing sorption model. Sorption model can be computed both in case t...
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:180
SorptionMob(Mesh &init_mesh, Input::Record in_rec)
Constructor.
Definition: sorption.cc:154
~SorptionImmob(void)
Destructor.
Definition: sorption.cc:244