Flow123d  jenkins-Flow123d-windows32-release-multijob-28
sorption.cc
Go to the documentation of this file.
1 #include <vector>
2 
3 #include "reaction/isotherm.hh"
4 #include "reaction/sorption.hh"
5 //#include "system/system.hh"
6 #include "system/sys_profiler.hh"
7 #include "mesh/accessors.hh"
8 
9 using namespace std;
10 
11 /********************************* *********************************************************/
12 /********************************* SORPTION_SIMPLE *********************************************************/
13 /********************************* *********************************************************/
14 
15 IT::Record SorptionSimple::input_type = SorptionBase::record_factory(SorptionRecord::simple);
16 
17 SorptionSimple::SorptionSimple(Mesh &init_mesh, Input::Record in_rec)
18  : SorptionBase(init_mesh, in_rec)
19 {
20  //DBGMSG("SorptionSimple constructor.\n");
21  data_ = new EqData("conc_solid");
22  this->eq_data_ = data_;
23  output_selection = make_output_selection("conc_solid", "SorptionSimple_Output");
24 }
25 
26 SorptionSimple::~SorptionSimple(void)
27 {}
28 
29 void SorptionSimple::isotherm_reinit(std::vector<Isotherm> &isotherms_vec, const ElementAccessor<3> &elem)
30 {
31  START_TIMER("SorptionSimple::isotherm_reinit");
32 
33  double rock_density = data_->rock_density.value(elem.centre(),elem);
34  double por_m = data_->porosity.value(elem.centre(),elem);
35 
36  // List of types of isotherms in particular regions
37  arma::uvec adsorption_type = data_->sorption_type.value(elem.centre(),elem);
38  arma::Col<double> mult_coef_vec = data_->isotherm_mult.value(elem.centre(),elem);
39  arma::Col<double> second_coef_vec = data_->isotherm_other.value(elem.centre(),elem);
40 
41  for(unsigned int i_subst = 0; i_subst < n_substances_; i_subst++)
42  {
43  double mult_coef = mult_coef_vec[i_subst];
44  double second_coef = second_coef_vec[i_subst];
45  Isotherm & isotherm = isotherms_vec[i_subst];
46 
47  //scales are different for the case of sorption in mobile and immobile pores
48  double scale_aqua = por_m,
49  scale_sorbed = (1 - por_m) * rock_density * molar_masses_[i_subst];
50 
51  //DBGMSG("molar_masses[%d %d]: %f\n",i_subst, substance_id[i_subst], molar_masses[i_subst]);
52  if ( scale_sorbed == 0.0)
53  xprintf(UsrErr, "Parameter scale_sorbed ((1 - por_m) * rock_density * molar_masses[i_subst]) is equal to zero.");
54  bool limited_solubility_on = false;
55  double table_limit;
56  if (solubility_vec_[i_subst] <= 0.0) {
57  limited_solubility_on = false;
58  table_limit=table_limit_[i_subst];
59  } else {
60  limited_solubility_on = true;
61  table_limit=solubility_vec_[i_subst];
62  }
63  isotherm.reinit(Isotherm::SorptionType(adsorption_type[i_subst]), limited_solubility_on,
64  solvent_density_, scale_aqua, scale_sorbed, table_limit, mult_coef, second_coef);
65 
66  }
67 
68  END_TIMER("SorptionSimple::isotherm_reinit");
69 }
70 
71 
72 /*********************************** *********************************************************/
73 /*********************************** SORPTION_DUAL *********************************************************/
74 /*********************************** *********************************************************/
75 
76 SorptionDual::SorptionDual(Mesh &init_mesh, Input::Record in_rec,
77  const string &output_conc_name,
78  const string &output_selection_name)
79  : SorptionBase(init_mesh, in_rec)
80 {
81  data_ = new EqData(output_conc_name);
82  *data_+=immob_porosity_
83  .flags_add(FieldFlag::input_copy)
84  .name("porosity_immobile");
85  this->eq_data_ = data_;
86  output_selection = make_output_selection(output_conc_name, output_selection_name);
87 }
88 
89 SorptionDual::~SorptionDual(void)
90 {}
91 
92 /*
93 void SorptionDual::isotherm_reinit(std::vector<Isotherm> &isotherms_vec, const ElementAccessor<3> &elem)
94 {
95  START_TIMER("SorptionMob::isotherm_reinit");
96 
97  double rock_density = data_.rock_density.value(elem.centre(),elem);
98 
99  double por_m = data_.porosity.value(elem.centre(),elem);
100  double por_imm = immob_porosity_.value(elem.centre(),elem);
101 
102  // List of types of isotherms in particular regions
103  arma::uvec adsorption_type = data_.sorption_types.value(elem.centre(),elem);
104  arma::Col<double> mult_coef_vec = data_.mult_coefs.value(elem.centre(),elem);
105  arma::Col<double> second_coef_vec = data_.second_params.value(elem.centre(),elem);
106 
107  for(int i_subst = 0; i_subst < n_substances_; i_subst++)
108  {
109  double mult_coef = mult_coef_vec[i_subst];
110  double second_coef = second_coef_vec[i_subst];
111  Isotherm & isotherm = isotherms_vec[i_subst];
112 
113  //scales are different for the case of sorption in mobile and immobile pores
114  double scale_aqua, scale_sorbed;
115  scale_aqua = por_imm;
116  scale_sorbed = compute_sorbing_scale(por_m,por_imm) * (1 - por_m - por_imm) * rock_density * molar_masses[i_subst];
117  if(scale_sorbed == 0.0)
118  xprintf(UsrErr, "Parameter scale_sorbed (phi * (1 - por_m - por_imm) * rock_density * molar_masses[i_subst]) is equal to zero.");
119  bool limited_solubility_on;
120  double table_limit;
121  if (solubility_vec_[i_subst] <= 0.0) {
122  limited_solubility_on = false;
123  table_limit=table_limit_[i_subst];
124 
125  } else {
126  limited_solubility_on = true;
127  table_limit=solubility_vec_[i_subst];
128  }
129  isotherm.reinit(Isotherm::SorptionType(adsorption_type[i_subst]), limited_solubility_on,
130  solvent_dens, scale_aqua, scale_sorbed, table_limit, mult_coef, second_coef);
131 
132  }
133 
134  END_TIMER("SorptionMob::isotherm_reinit");
135 
136  return;
137 }
138 */
139 
140 /********************************** *******************************************************/
141 /*********************************** SORPTION_MOBILE *******************************************************/
142 /********************************** *******************************************************/
143 
144 IT::Record SorptionMob::input_type = SorptionBase::record_factory(SorptionRecord::mobile);
145 
146 SorptionMob::SorptionMob(Mesh &init_mesh, Input::Record in_rec)
147  : SorptionDual(init_mesh, in_rec, "conc_solid", "SorptionMobile_Output")
148 {}
149 
150 
151 SorptionMob::~SorptionMob(void)
152 {}
153 
154 /*
155 double SorptionMob::compute_sorbing_scale(double por_m, double por_imm)
156 {
157  double phi = por_m/(por_m + por_imm);
158  return phi;
159 }
160 */
161 
162 void SorptionMob::isotherm_reinit(std::vector<Isotherm> &isotherms_vec, const ElementAccessor<3> &elem)
163 {
164  START_TIMER("SorptionMob::isotherm_reinit");
165 
166  double rock_density = data_->rock_density.value(elem.centre(),elem);
167 
168  double por_m = data_->porosity.value(elem.centre(),elem);
169  double por_imm = immob_porosity_.value(elem.centre(),elem);
170  double phi = por_m/(por_m + por_imm);
171 
172  // List of types of isotherms in particular regions
173  arma::uvec adsorption_type = data_->sorption_type.value(elem.centre(),elem);
174  arma::Col<double> mult_coef_vec = data_->isotherm_mult.value(elem.centre(),elem);
175  arma::Col<double> second_coef_vec = data_->isotherm_other.value(elem.centre(),elem);
176 
177  for(unsigned int i_subst = 0; i_subst < n_substances_; i_subst++)
178  {
179  double mult_coef = mult_coef_vec[i_subst];
180  double second_coef = second_coef_vec[i_subst];
181  Isotherm & isotherm = isotherms_vec[i_subst];
182 
183  //scales are different for the case of sorption in mobile and immobile pores
184  double scale_aqua, scale_sorbed;
185  scale_aqua = por_m;
186  scale_sorbed = phi * (1 - por_m - por_imm) * rock_density * molar_masses_[i_subst];
187  if(scale_sorbed == 0.0)
188  xprintf(UsrErr, "Parameter scale_sorbed (phi * (1 - por_m - por_imm) * rock_density * molar_masses[i_subst]) is equal to zero.");
189  bool limited_solubility_on;
190  double table_limit;
191  if (solubility_vec_[i_subst] <= 0.0) {
192  limited_solubility_on = false;
193  table_limit=table_limit_[i_subst];
194 
195  } else {
196  limited_solubility_on = true;
197  table_limit=solubility_vec_[i_subst];
198  }
199  isotherm.reinit(Isotherm::SorptionType(adsorption_type[i_subst]), limited_solubility_on,
200  solvent_density_, scale_aqua, scale_sorbed, table_limit, mult_coef, second_coef);
201 
202  }
203 
204  END_TIMER("SorptionMob::isotherm_reinit");
205 }
206 
207 
208 /*********************************** *****************************************************/
209 /*********************************** SORPTION_IMMOBILE *****************************************************/
210 /*********************************** *****************************************************/
211 
212 IT::Record SorptionImmob::input_type = SorptionBase::record_factory(SorptionRecord::immobile);
213 
214 SorptionImmob::SorptionImmob(Mesh &init_mesh, Input::Record in_rec)
215 : SorptionDual(init_mesh, in_rec, "conc_immobile_solid", "SorptionImmobile_Output")
216 {}
217 
218 SorptionImmob::~SorptionImmob(void)
219 {}
220 
221 /*
222 double SorptionImmob::compute_sorbing_scale(double por_m, double por_imm)
223 {
224  double phi = por_imm / (por_m + por_imm);
225  return phi;
226 }
227 */
228 
229 void SorptionImmob::isotherm_reinit(std::vector<Isotherm> &isotherms_vec, const ElementAccessor<3> &elem)
230 {
231  START_TIMER("SorptionImmob::isotherm_reinit");
232 
233  double rock_density = data_->rock_density.value(elem.centre(),elem);
234 
235  double por_m = data_->porosity.value(elem.centre(),elem);
236  double por_imm = immob_porosity_.value(elem.centre(),elem);
237  //double phi = por_imm / (por_m + por_imm); // = 1-phi
238  double phi = por_m/(por_m + por_imm);
239 
240  // List of types of isotherms in particular regions
241  arma::uvec adsorption_type = data_->sorption_type.value(elem.centre(),elem);
242  arma::Col<double> mult_coef_vec = data_->isotherm_mult.value(elem.centre(),elem);
243  arma::Col<double> second_coef_vec = data_->isotherm_other.value(elem.centre(),elem);
244 
245  for(unsigned int i_subst = 0; i_subst < n_substances_; i_subst++)
246  {
247  double mult_coef = mult_coef_vec[i_subst];
248  double second_coef = second_coef_vec[i_subst];
249  Isotherm & isotherm = isotherms_vec[i_subst];
250 
251  //scales are different for the case of sorption in mobile and immobile pores
252  double scale_aqua, scale_sorbed;
253  scale_aqua = por_imm;
254  scale_sorbed = (1 - phi) * (1 - por_m - por_imm) * rock_density * molar_masses_[i_subst];
255  if(scale_sorbed == 0.0)
256  xprintf(UsrErr, "Parameter scale_sorbed ((1 - phi) * (1 - por_m - por_imm) * rock_density * molar_masses[i_subst]) is equal to zero.");
257  bool limited_solubility_on;
258  double table_limit;
259  if (solubility_vec_[i_subst] <= 0.0) {
260  limited_solubility_on = false;
261  table_limit=table_limit_[i_subst];
262 
263  } else {
264  limited_solubility_on = true;
265  table_limit=solubility_vec_[i_subst];
266  }
267  isotherm.reinit(Isotherm::SorptionType(adsorption_type[i_subst]), limited_solubility_on,
268  solvent_density_, scale_aqua, scale_sorbed, table_limit, mult_coef, second_coef);
269 
270  }
271 
272  END_TIMER("SorptionImmob::isotherm_reinit");
273 }
SorptionMob::input_type
static Input::Type::Record input_type
Definition: sorption.hh:78
EquationBase::eq_data_
FieldSet * eq_data_
Definition: equation.hh:237
SorptionBase::solvent_density_
double solvent_density_
Definition: sorption_base.hh:173
SorptionDual::~SorptionDual
~SorptionDual(void)
Destructor.
Definition: sorption.cc:89
SorptionBase::EqData::sorption_type
Field< 3, FieldValue< 3 >::EnumVector > sorption_type
Discrete need Selection for initialization.
Definition: sorption_base.hh:58
FieldCommon::flags_add
FieldCommon & flags_add(FieldFlag::Flags::Mask mask)
Definition: field_common.hh:143
SorptionBase::EqData::isotherm_other
Field< 3, FieldValue< 3 >::Vector > isotherm_other
Langmuir sorption coeficients alpha (in fraction c_s = omega * (alpha*c_a)/(1- alpha*c_a)).
Definition: sorption_base.hh:65
SorptionBase::output_selection
Input::Type::Selection output_selection
Definition: sorption_base.hh:202
SorptionMob::~SorptionMob
~SorptionMob(void)
Destructor.
Definition: sorption.cc:151
SorptionDual
Abstract class of sorption model in case dual porosity is considered.
Definition: sorption.hh:45
SorptionImmob::input_type
static Input::Type::Record input_type
Definition: sorption.hh:100
SorptionBase::EqData::rock_density
Field< 3, FieldValue< 3 >::Scalar > rock_density
Rock matrix density.
Definition: sorption_base.hh:59
SorptionMob::isotherm_reinit
void isotherm_reinit(std::vector< Isotherm > &isotherms_vec, const ElementAccessor< 3 > &elem) override
Reinitializes the isotherm.
Definition: sorption.cc:162
Mesh
Definition: mesh.h:108
SorptionBase::molar_masses_
std::vector< double > molar_masses_
Definition: sorption_base.hh:168
SorptionBase::table_limit_
std::vector< double > table_limit_
Definition: sorption_base.hh:181
SorptionBase::EqData::isotherm_mult
Field< 3, FieldValue< 3 >::Vector > isotherm_mult
Multiplication coefficients (k, omega) for all types of isotherms.
Definition: sorption_base.hh:63
Isotherm::reinit
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:295
sys_profiler.hh
SorptionBase::EqData::porosity
Field< 3, FieldValue< 3 >::Scalar > porosity
Porosity field copied from transport.
Definition: sorption_base.hh:68
SorptionDual::immob_porosity_
Field< 3, FieldValue< 3 >::Scalar > immob_porosity_
Definition: sorption.hh:66
FieldFlag::input_copy
static constexpr Mask input_copy
A field that is input of its equation and cna not read from input, thus muzt be set by copy...
Definition: field_flag.hh:29
SorptionBase::solubility_vec_
std::vector< double > solubility_vec_
Definition: sorption_base.hh:177
Input::Record
Accessor to the data with type Type::Record.
Definition: accessors.hh:308
xprintf
#define xprintf(...)
Definition: system.hh:104
START_TIMER
#define START_TIMER(tag)
Starts a timer with specified tag.
Definition: sys_profiler.hh:114
Field::value
virtual Value::return_type const & value(const Point &p, const ElementAccessor< spacedim > &elm) const
Definition: field.hh:399
SorptionBase::record_factory
static Input::Type::Record record_factory(SorptionRecord::Type)
Creates the input record for different cases of sorption model (simple or in dual porosity)...
Definition: sorption_base.cc:86
SorptionBase::make_output_selection
static Input::Type::Selection make_output_selection(const string &output_field_name, const string &selection_name)
Definition: sorption_base.hh:41
Isotherm::SorptionType
SorptionType
Type of adsorption isotherm.
Definition: isotherm.hh:138
UsrErr
Definition: system.hh:75
SorptionImmob::SorptionImmob
SorptionImmob(Mesh &init_mesh, Input::Record in_rec)
Constructor.
Definition: sorption.cc:214
SorptionSimple::input_type
static Input::Type::Record input_type
Definition: sorption.hh:28
SorptionSimple::isotherm_reinit
void isotherm_reinit(std::vector< Isotherm > &isotherms, const ElementAccessor< 3 > &elm) override
Reinitializes the isotherm.
Definition: sorption.cc:29
SorptionSimple::SorptionSimple
SorptionSimple(Mesh &init_mesh, Input::Record in_rec)
Constructor.
Definition: sorption.cc:17
SorptionSimple::~SorptionSimple
~SorptionSimple(void)
Destructor.
Definition: sorption.cc:26
FieldCommon::name
FieldCommon & name(const string &name)
Definition: field_common.hh:71
END_TIMER
#define END_TIMER(tag)
Ends a timer with specified tag.
Definition: sys_profiler.hh:127
ElementAccessor::centre
arma::vec::fixed< spacedim > centre() const
Definition: accessors.hh:81
Input::Type::Record
Record type proxy class.
Definition: type_record.hh:161
SorptionBase::n_substances_
unsigned int n_substances_
Definition: sorption_base.hh:190
SorptionDual::SorptionDual
SorptionDual(Mesh &init_mesh, Input::Record in_rec, const string &output_conc_name, const string &output_selection_name)
Constructor.
Definition: sorption.cc:76
SorptionBase::data_
EqData * data_
Pointer to equation data. The object is constructed in descendants.
Definition: sorption_base.hh:159
SorptionImmob::isotherm_reinit
void isotherm_reinit(std::vector< Isotherm > &isotherms_vec, const ElementAccessor< 3 > &elem) override
Reinitializes the isotherm.
Definition: sorption.cc:229
SorptionMob::SorptionMob
SorptionMob(Mesh &init_mesh, Input::Record in_rec)
Constructor.
Definition: sorption.cc:146
SorptionImmob::~SorptionImmob
~SorptionImmob(void)
Destructor.
Definition: sorption.cc:218
Generated on Tue Aug 19 2014 23:38:30 for Flow123d by   doxygen 1.8.7