Flow123d  release_3.0.0-875-gdc24e59
sorption_base.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 sorption_base.hh
15  * @brief Class SorptionBase is abstract class representing model of sorption in transport.
16  *
17  * The sorption is described by several types of isotherms - linear, Freundlich or Langmuir.
18  * Limited solubility can be considered.
19  *
20  * Interpolation tables are used to speed up evaluation of isotherms.
21  *
22  */
23 
24 #ifndef SORPTION_BASE_H
25 #define SORPTION_BASE_H
26 
27 #include <boost/exception/info.hpp> // for operator<<, error_info::~erro...
28 #include <memory> // for shared_ptr
29 #include <string> // for string
30 #include <vector>
31 #include "reaction/reaction_term.hh" // for ReactionTerm
32 #include "fields/field.hh" // for Field
33 #include "fields/field_values.hh" // for FieldValue<>::Scalar, FieldVa...
34 #include "fields/field_set.hh"
35 #include "fields/multi_field.hh"
36 #include "la/vector_mpi.hh"
38 #include "input/input_exception.hh" // for DECLARE_INPUT_EXCEPTION, Exce...
39 #include "input/type_base.hh" // for Array
40 #include "input/type_generic.hh" // for Instance
41 #include "petscvec.h" // for Vec, VecScatter, _p_VecScatter
42 #include "system/exceptions.hh" // for operator<<, ExcStream, EI
43 
44 class Isotherm;
45 class Mesh;
46 namespace Input {
47  class Record;
48  namespace Type {
49  class Record;
50  class Selection;
51  }
52 }
53 template <int spacedim> class ElementAccessor;
54 
55 
56 
57 
59 {
60 public:
61  TYPEDEF_ERR_INFO( EI_ArrayName, std::string);
62  DECLARE_INPUT_EXCEPTION( ExcSubstanceCountMatch, << "The size of the input array " << EI_ArrayName::qval
63  << " does not match the number of substances.");
64 
65  /**
66  * Static variable for new input data types input
67  */
68  static const Input::Type::Record & get_input_type();
69 
70  static Input::Type::Instance make_output_type(const string &equation_name, const string &output_field_name, const string &output_field_desc )
71  {
72  return EqData(output_field_name, output_field_desc).output_fields.make_output_type(equation_name, "");
73  }
74 
75  class EqData : public FieldSet
76  {
77  public:
78  /**
79  * Sorption type specifies a kind of equilibrial description of adsorption.
80  */
81  static const Input::Type::Selection & get_sorption_type_selection();
82 
83  /// Collect all fields
84  EqData(const string &output_field_name, const string &output_field_desc);
85 
86  MultiField<3, FieldValue<3>::Enum > sorption_type; ///< Discrete need Selection for initialization.
87  Field<3, FieldValue<3>::Scalar > rock_density; ///< Rock matrix density.
88 
89  /// Multiplication coefficients (k, omega) for all types of isotherms.
90  /** Langmuir: c_s = omega * (alpha*c_a)/(1- alpha*c_a), Linear: c_s = k*c_a */
92  /// Langmuir sorption coeficients alpha (in fraction c_s = omega * (alpha*c_a)/(1- alpha*c_a)).
94 
95  MultiField<3, FieldValue<3>::Scalar> init_conc_solid; ///< Initial sorbed concentrations.
96  Field<3, FieldValue<3>::Scalar > porosity; ///< Porosity field copied from transport.
97 
98  MultiField<3, FieldValue<3>::Scalar> conc_solid; ///< Calculated sorbed concentrations, for output only.
99 
100  /// Input data set - fields in this set are read from the input file.
102 
103  /// Fields indended for output, i.e. all input fields plus those representing solution.
105 
106  };
107 
108  /**
109  * Constructor with parameter for initialization of a new declared class member
110  */
111  SorptionBase(Mesh &init_mesh, Input::Record in_rec);
112  /**
113  * Destructor.
114  */
115  virtual ~SorptionBase(void);
116 
117  /// Prepares the object to usage.
118  /**
119  * Allocating memory, reading input, initialization of fields.
120  */
121  void initialize() override;
122 
123  /**
124  * Does first computation after initialization process.
125  * The time is set and initial condition is set and output.
126  */
127  void zero_time_step() override;
128 
129  /// Updates the solution.
130  /**
131  * Goes through local distribution of elements and calls @p compute_reaction.
132  */
133  void update_solution(void) override;
134 
135  void output_data(void) override;
136 
137  bool evaluate_time_constraint(double &time_constraint) override { return false; }
138 
139 
140 protected:
141  /**
142  * This method disables to use constructor without parameters.
143  */
144  SorptionBase();
145 
146  /** Initializes possible following reactions from input record.
147  * It should be called after setting mesh, time_governor, distribution and concentration_matrix
148  * if there are some setting methods for reactions called (they are not at the moment, so it could be part of init_from_input).
149  */
150  void make_reactions();
151 
152  /// Reads names of substances from input and creates indexing to global vector of substance.
153  /** Also creates the local vector of molar masses. */
154  void initialize_substance_ids();
155 
156  /// Initializes private members of sorption from the input record.
157  void initialize_from_input();
158 
159  /// Initializes field sets.
160  void initialize_fields();
161 
162  ///Reads and sets initial condition for concentration in solid.
163  void set_initial_condition();
164 
165  /**
166  * For simulation of sorption in just one element either inside of MOBILE or IMMOBILE pores.
167  */
168  double **compute_reaction(double **concentrations, int loc_el) override;
169 
170  /// Reinitializes the isotherm.
171  /**
172  * On data change the isotherm is recomputed, possibly new interpolation table is made.
173  * NOTE: Be sure to update common element data (porosity, rock density etc.)
174  * by @p compute_common_ele_data(), before calling reinitialization!
175  */
176  void isotherm_reinit(unsigned int i_subst, const ElementAccessor<3> &elm);
177 
178  /// Calls @p isotherm_reinit for all isotherms.
179  void isotherm_reinit_all(const ElementAccessor<3> &elm);
180 
181  /**
182  * Creates interpolation table for isotherms.
183  */
184  void make_tables(void);
185 
186  /// Computes maximal aqueous concentration at the current step.
187  void update_max_conc();
188 
189  /// Sets max conc to zeros on all regins.
190  void clear_max_conc();
191 
192  /// Pointer to equation data. The object is constructed in descendants.
194 
195  /**
196  * Temporary nr_of_points can be computed using step_length. Should be |nr_of_region x nr_of_substances| matrix later.
197  */
199  /**
200  * Density of the solvent.
201  * TODO: Could be done region dependent, easily.
202  */
204  /**
205  * Critical concentrations of species dissolved in water.
206  */
208  /**
209  * Concentration table limits of species dissolved in water.
210  */
212  /**
213  * Maximum concentration per region.
214  * It is used for optimization of interpolation table.
215  */
217  /**
218  * Three dimensional array contains intersections between isotherms and mass balance lines.
219  * It describes behaviour of sorbents in mobile pores of various rock matrix enviroments.
220  * Up to |nr_of_region x nr_of_substances x n_points| doubles. Because of equidistant step
221  * lenght in cocidered system of coordinates, just function values are stored.
222  */
224 
225  unsigned int n_substances_; //< number of substances that take part in the sorption mode
226 
227  /// Mapping from local indexing of substances to global.
229 
230  /**
231  * Array for storage infos about sorbed species concentrations.
232  */
233  double** conc_solid;
234 
235  /**
236  * Reaction model that follows the sorption.
237  */
238  std::shared_ptr<ReactionTerm> reaction_liquid;
239  std::shared_ptr<ReactionTerm> reaction_solid;
240 
241  ///@name members used in output routines
242  //@{
243  std::vector<VectorMPI> conc_solid_out; ///< sorbed concentration array output (gathered - sequential)
244  //@}
245 
246  // Temporary objects holding pointers to appropriate FieldFE
247  // TODO remove after final fix of equations
248  /// Fields correspond with \p conc_solid_out.
250 
251  /** Structure for data respectful to element, but indepedent of actual isotherm.
252  * Reads mobile/immobile porosity, rock density and then computes concentration scaling parameters.
253  * Is kind of optimization, so that these data are computed only when necessary.
254  */
256  double scale_aqua;
257  double scale_sorbed;
259  } common_ele_data;
260 
261  /** Computes @p CommonElementData.
262  * Is pure virtual, implemented differently for simple/mobile/immobile sorption class.
263  */
264  virtual void compute_common_ele_data(const ElementAccessor<3> &elem) = 0;
265 };
266 
267 #endif //SORPTION_BASE_H
std::shared_ptr< ReactionTerm > reaction_solid
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:71
MultiField< 3, FieldValue< 3 >::Scalar > conc_solid
Calculated sorbed concentrations, for output only.
double solvent_density_
std::vector< std::vector< Isotherm > > isotherms
TYPEDEF_ERR_INFO(EI_KeyName, const string)
static Input::Type::Instance make_output_type(const string &equation_name, const string &output_field_name, const string &output_field_desc)
Abstract linear system class.
Definition: balance.hh:35
Class template representing a field with values dependent on: point, element, and region...
Definition: field.hh:83
Field< 3, FieldValue< 3 >::Scalar > rock_density
Rock matrix density.
Definition: mesh.h:80
EquationOutput output_fields
Fields indended for output, i.e. all input fields plus those representing solution.
Helper class that stores data of generic types.
Definition: type_generic.hh:89
std::vector< double > table_limit_
Class ReactionTerm is an abstract class representing reaction term in transport.
std::vector< unsigned int > substance_global_idx_
Mapping from local indexing of substances to global.
std::vector< VectorMPI > conc_solid_out
sorbed concentration array output (gathered - sequential)
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.
std::vector< std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > > output_field_ptr
Fields correspond with conc_solid_out.
MultiField< 3, FieldValue< 3 >::Scalar > init_conc_solid
Initial sorbed concentrations.
std::vector< double > solubility_vec_
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
unsigned int n_interpolation_steps_
MultiField< 3, FieldValue< 3 >::Enum > sorption_type
Discrete need Selection for initialization.
double ** conc_solid
Record type proxy class.
Definition: type_record.hh:182
const Input::Type::Instance & make_output_type(const string &equation_name, const string &aditional_description="")
Class for representation of a vector of fields of the same physical quantity.
Definition: multi_field.hh:87
MultiField< 3, FieldValue< 3 >::Scalar > distribution_coefficient
Multiplication coefficients (k, omega) for all types of isotherms.
unsigned int n_substances_
DECLARE_INPUT_EXCEPTION(ExcInputMessage,<< EI_Message::val)
Simple input exception that accepts just string message.
bool evaluate_time_constraint(double &time_constraint) override
Computes a constraint for time step.
EqData * data_
Pointer to equation data. The object is constructed in descendants.
std::vector< std::vector< double > > max_conc
Template for classes storing finite set of named values.
FieldSet input_data_set_
Input data set - fields in this set are read from the input file.
std::shared_ptr< ReactionTerm > reaction_liquid