Flow123d  master-f44eb46
assembly_reaction.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 assembly_reaction.hh
15  * @brief
16  */
17 
18 #ifndef ASSEMBLY_REACTION_HH_
19 #define ASSEMBLY_REACTION_HH_
20 
24 #include "reaction/isotherm.hh"
25 #include "fem/fe_p.hh"
26 #include "fem/fe_values.hh"
28 //#include "coupling/balance.hh"
30 
31 
32 
33 template <unsigned int dim>
35 {
36 public:
37  typedef typename DualPorosity::EqFields EqFields;
38  typedef typename DualPorosity::EqData EqData;
39 
40  static constexpr const char * name() { return "InitConditionAssemblyDp"; }
41 
42  /// Constructor.
43  InitConditionAssemblyDp(EqFields *eq_fields, EqData *eq_data)
44  : AssemblyBase<dim>(0), eq_fields_(eq_fields), eq_data_(eq_data) {
47  }
48 
49  /// Destructor.
51 
52  /// Initialize auxiliary vectors and other data members
53  void initialize(ElementCacheMap *element_cache_map) {
54  //this->balance_ = eq_data_->balance_;
55  this->element_cache_map_ = element_cache_map;
56  }
57 
58 
59  /// Assemble integral over element
60  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
61  {
62  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
63 
64  dof_p0_ = cell.get_loc_dof_indices()[0];
65  auto p = *( this->bulk_points(element_patch_idx).begin() );
66 
67  //setting initial solid concentration for substances involved in adsorption
68  for (unsigned int sbi = 0; sbi < eq_data_->substances_.size(); sbi++)
69  {
70  eq_fields_->conc_immobile_fe[sbi]->vec().set( dof_p0_, eq_fields_->init_conc_immobile[sbi](p) );
71  }
72  }
73 
74 
75 private:
76  /// Data objects shared with Elasticity
79 
80  /// Sub field set contains fields used in calculation.
82 
83  IntIdx dof_p0_; ///< Index of local DOF
84 
85  template < template<IntDim...> class DimAssembly>
86  friend class GenericAssembly;
87 };
88 
89 template <unsigned int dim>
90 class ReactionAssemblyDp : public AssemblyBase<dim>
91 {
92 public:
93  typedef typename DualPorosity::EqFields EqFields;
94  typedef typename DualPorosity::EqData EqData;
95 
96  static constexpr const char * name() { return "InitConditionAssemblyDp"; }
97 
98  /// Constructor.
99  ReactionAssemblyDp(EqFields *eq_fields, EqData *eq_data)
100  : AssemblyBase<dim>(0), eq_fields_(eq_fields), eq_data_(eq_data) {
102  this->used_fields_ += eq_fields_->porosity;
105  }
106 
107  /// Destructor.
109 
110  /// Initialize auxiliary vectors and other data members
111  void initialize(ElementCacheMap *element_cache_map) {
112  //this->balance_ = eq_data_->balance_;
113  this->element_cache_map_ = element_cache_map;
114  }
115 
116 
117  /// Assemble integral over element
118  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
119  {
120  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
121 
122  auto p = *( this->bulk_points(element_patch_idx).begin() );
123 
124  // if porosity_immobile == 0 then mobile concentration stays the same
125  // and immobile concentration cannot change
127  if (por_immob_ == 0.0) return;
128 
129  // get data from fields
130  dof_p0_ = cell.get_loc_dof_indices()[0];
132  arma::Col<double> diff_vec(eq_data_->substances_.size());
133  for (sbi_=0; sbi_<eq_data_->substances_.size(); sbi_++) // Optimize: SWAP LOOPS
134  diff_vec[sbi_] = eq_fields_->diffusion_rate_immobile[sbi_](p);
135 
137 
138  for (sbi_ = 0; sbi_ < eq_data_->substances_.size(); sbi_++) //over all substances
139  {
140  exponent_ = diff_vec[sbi_] * temp_exponent_;
141  //previous values
144 
145  // ---compute average concentration------------------------------------------
147  / (por_mob_ + por_immob_);
148 
150 
151  // the following 2 conditions guarantee:
152  // 1) stability of forward Euler's method
153  // 2) the error of forward Euler's method will not be large
154  if(eq_data_->time_->dt() <= por_mob_*por_immob_/(max(diff_vec)*(por_mob_+por_immob_)) &&
156  {
158  // ---compute concentration in mobile area
160 
161  // ---compute concentration in immobile area
163  }
164  else //analytic solution
165  {
166  temp_ = exp(-exponent_);
167  // ---compute concentration in mobile area
169 
170  // ---compute concentration in immobile area
172  }
173 
176  }
177  }
178 
179 
180 private:
181  /// Data objects shared with Elasticity
184 
185  /// Sub field set contains fields used in calculation.
187 
188  unsigned int sbi_; ///< Index of substance
189  IntIdx dof_p0_; ///< Index of local DOF
190  double conc_average_; ///< weighted (by porosity) average of concentration
191  double conc_mob_, conc_immob_; ///< new mobile and immobile concentration
192  double previous_conc_mob_, previous_conc_immob_; ///< mobile and immobile concentration in previous time step
193  double conc_max_; ///< difference between concentration and average concentration
194  double por_mob_, por_immob_; ///< mobile and immobile porosity
195  double exponent_, temp_exponent_, temp_; ///< Precomputed values
196 
197  template < template<IntDim...> class DimAssembly>
198  friend class GenericAssembly;
199 };
200 
201 template <unsigned int dim>
203 {
204 public:
206  typedef typename SorptionBase::EqData EqData;
207 
208  static constexpr const char * name() { return "InitConditionAssemblySorp"; }
209 
210  /// Constructor.
212  : AssemblyBase<dim>(0), eq_fields_(eq_fields), eq_data_(eq_data) {
215  }
216 
217  /// Destructor.
219 
220  /// Initialize auxiliary vectors and other data members
221  void initialize(ElementCacheMap *element_cache_map) {
222  //this->balance_ = eq_data_->balance_;
223  this->element_cache_map_ = element_cache_map;
224  }
225 
226 
227  /// Assemble integral over element
228  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
229  {
230  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
231 
232  dof_p0_ = cell.get_loc_dof_indices()[0];
233  auto p = *( this->bulk_points(element_patch_idx).begin() );
234 
235  //setting initial solid concentration for substances involved in adsorption
236  for (unsigned int sbi = 0; sbi < eq_data_->n_substances_; sbi++)
237  {
239  }
240  }
241 
242 
243 private:
244  /// Data objects shared with Elasticity
247 
248  /// Sub field set contains fields used in calculation.
250 
251  IntIdx dof_p0_; ///< Index of local DOF
252 
253  template < template<IntDim...> class DimAssembly>
254  friend class GenericAssembly;
255 };
256 
257 template <unsigned int dim>
259 {
260 public:
262  typedef typename SorptionBase::EqData EqData;
263 
264  static constexpr const char * name() { return "ReactionAssemblySorp"; }
265 
266  /// Constructor.
267  ReactionAssemblySorp(EqFields *eq_fields, EqData *eq_data)
268  : AssemblyBase<dim>(0), eq_fields_(eq_fields), eq_data_(eq_data) {
276  }
277 
278  /// Destructor.
280 
281  /// Initialize auxiliary vectors and other data members
282  void initialize(ElementCacheMap *element_cache_map) {
283  //this->balance_ = eq_data_->balance_;
284  this->element_cache_map_ = element_cache_map;
285  }
286 
287 
288  /// Assemble integral over element
289  inline void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
290  {
291  ASSERT_EQ(cell.dim(), dim).error("Dimension of element mismatch!");
292 
293  unsigned int i_subst, subst_id;
294  auto p = *( this->bulk_points(element_patch_idx).begin() );
295 
296  reg_idx_ = cell.elm().region().bulk_idx();
297  dof_p0_ = cell.get_loc_dof_indices()[0];
298 
299  double scale_aqua = eq_fields_->scale_aqua(p);
300  double scale_sorbed = eq_fields_->scale_sorbed(p);
301  double no_sorbing_surface_cond = eq_fields_->no_sorbing_surface_cond(p);
302 
303  for(i_subst = 0; i_subst < eq_data_->n_substances_; i_subst++)
304  {
305  subst_id = eq_data_->substance_global_idx_[i_subst];
306 
307  Isotherm & isotherm = eq_data_->isotherms[reg_idx_][i_subst];
308 
309  bool limited_solubility_on = eq_data_->solubility_vec_[i_subst] > 0.0;
310 
311  // in case of no sorbing surface, set isotherm type None
312  if( no_sorbing_surface_cond <= std::numeric_limits<double>::epsilon())
313  {
315  scale_aqua, scale_sorbed,
316  0,0,0);
317  } else {
318  if ( scale_sorbed <= 0.0)
319  THROW( SorptionBase::ExcNotPositiveScaling() << SorptionBase::EI_Subst(i_subst) );
320 
321  isotherm.reinit(Isotherm::SorptionType(eq_fields_->sorption_type[i_subst](p)),
322  limited_solubility_on, eq_data_->solvent_density_,
323  scale_aqua, scale_sorbed,
324  eq_data_->solubility_vec_[i_subst],
326  eq_fields_->isotherm_other[i_subst](p));
327  }
328 
329  double c_aqua = eq_fields_->conc_mobile_fe[subst_id]->vec().get(dof_p0_);
330  double c_sorbed = eq_fields_->conc_solid_fe[subst_id]->vec().get(dof_p0_);
331  isotherm.compute(c_aqua, c_sorbed);
332  eq_fields_->conc_mobile_fe[subst_id]->vec().set(dof_p0_, c_aqua);
333  eq_fields_->conc_solid_fe[subst_id]->vec().set(dof_p0_, c_sorbed);
334 
335  // update maximal concentration per region (optimization for interpolation)
336  if(eq_data_->table_limit_[i_subst] < 0)
337  eq_data_->max_conc[reg_idx_][i_subst] = std::max(eq_data_->max_conc[reg_idx_][i_subst],
338  eq_fields_->conc_mobile_fe[subst_id]->vec().get(dof_p0_));
339  }
340  }
341 
342 
343 private:
344  /// Data objects shared with Elasticity
347 
348  /// Sub field set contains fields used in calculation.
350 
351  IntIdx dof_p0_; ///< Index of local DOF
352  int reg_idx_; ///< Bulk region idx
353  //unsigned int sbi_; ///< Index of substance
354 
355  template < template<IntDim...> class DimAssembly>
356  friend class GenericAssembly;
357 };
358 
359 #endif /* ASSEMBLY_REACTION_HH_ */
360 
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
Range< BulkPoint > bulk_points(unsigned int element_patch_idx) const
Return BulkPoint range of appropriate dimension.
int active_integrals_
Holds mask of active integrals.
ElementCacheMap * element_cache_map_
ElementCacheMap shared with GenericAssembly object.
Cell accessor allow iterate over DOF handler cells.
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
unsigned int dim() const
Return dimension of element appropriate to cell.
DualPorosity data.
TimeGovernor * time_
TimeGovernor object shared with assembly classes.
double scheme_tolerance_
Dual porosity computational scheme tolerance.
DualPorosity fields.
FieldFEScalarVec conc_immobile_fe
Underlaying FieldFE for each substance of conc_immobile.
MultiField< 3, FieldValue< 3 >::Scalar > init_conc_immobile
Initial concentrations in the immobile zone.
MultiField< 3, FieldValue< 3 >::Scalar > diffusion_rate_immobile
Mass transfer coefficients between mobile and immobile pores.
Field< 3, FieldValue< 3 >::Scalar > porosity_immobile
Immobile porosity field.
Field< 3, FieldValue< 3 >::Scalar > porosity
Porosity field.
Region region() const
Definition: accessors.hh:198
Directing class of FieldValueCache.
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
Generic class of assemblation.
FieldSet used_fields_
Sub field set contains fields used in calculation.
static constexpr const char * name()
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
InitConditionAssemblyDp(EqFields *eq_fields, EqData *eq_data)
Constructor.
DualPorosity::EqFields EqFields
IntIdx dof_p0_
Index of local DOF.
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
DualPorosity::EqData EqData
~InitConditionAssemblyDp()
Destructor.
EqFields * eq_fields_
Data objects shared with Elasticity.
SorptionBase::EqData EqData
IntIdx dof_p0_
Index of local DOF.
static constexpr const char * name()
EqFields * eq_fields_
Data objects shared with Elasticity.
InitConditionAssemblySorp(EqFields *eq_fields, EqData *eq_data)
Constructor.
FieldSet used_fields_
Sub field set contains fields used in calculation.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
SorptionBase::EqFields EqFields
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
void compute(double &c_aqua, double &c_sorbed)
Definition: isotherm.cc:68
SorptionType
Type of adsorption isotherm.
Definition: isotherm.hh:174
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 conc_max_
difference between concentration and average concentration
double conc_immob_
new mobile and immobile concentration
DualPorosity::EqData EqData
double por_immob_
mobile and immobile porosity
unsigned int sbi_
Index of substance.
ReactionAssemblyDp(EqFields *eq_fields, EqData *eq_data)
Constructor.
DualPorosity::EqFields EqFields
EqFields * eq_fields_
Data objects shared with Elasticity.
static constexpr const char * name()
double previous_conc_immob_
mobile and immobile concentration in previous time step
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
double temp_
Precomputed values.
IntIdx dof_p0_
Index of local DOF.
double conc_average_
weighted (by porosity) average of concentration
FieldSet used_fields_
Sub field set contains fields used in calculation.
~ReactionAssemblyDp()
Destructor.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
int reg_idx_
Bulk region idx.
EqFields * eq_fields_
Data objects shared with Elasticity.
SorptionBase::EqData EqData
void cell_integral(DHCellAccessor cell, unsigned int element_patch_idx)
Assemble integral over element.
IntIdx dof_p0_
Index of local DOF.
FieldSet used_fields_
Sub field set contains fields used in calculation.
~ReactionAssemblySorp()
Destructor.
ReactionAssemblySorp(EqFields *eq_fields, EqData *eq_data)
Constructor.
void initialize(ElementCacheMap *element_cache_map)
Initialize auxiliary vectors and other data members.
static constexpr const char * name()
SorptionBase::EqFields EqFields
SubstanceList substances_
FieldFEScalarVec conc_mobile_fe
FieldFEs representing P0 interpolation of mobile concentration (passed from transport).
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
Definition: region.hh:90
DualPorosity data.
std::vector< unsigned int > substance_global_idx_
Mapping from local indexing of substances to global.
std::vector< double > table_limit_
unsigned int n_substances_
number of substances that take part in the sorption mode
std::vector< std::vector< double > > max_conc
std::vector< std::vector< Isotherm > > isotherms
std::vector< double > solubility_vec_
MultiField< 3, FieldValue< 3 >::Scalar > init_conc_solid
Initial sorbed concentrations.
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 > scale_sorbed
FieldFEScalarVec conc_solid_fe
Underlaying FieldFE for each substance of conc_solid.
MultiField< 3, FieldValue< 3 >::Enum > sorption_type
Discrete need Selection for initialization.
Field< 3, FieldValue< 3 >::Scalar > scale_aqua
Instances of FieldModel used in assembly methods.
Field< 3, FieldValue< 3 >::Scalar > no_sorbing_surface_cond
MultiField< 3, FieldValue< 3 >::Scalar > distribution_coefficient
Multiplication coefficients (k, omega) for all types of isotherms.
unsigned int size() const
Definition: substance.hh:87
double dt() const
Class Dual_por_exchange implements the model of dual porosity.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
@ bulk
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
int IntIdx
Definition: index_types.hh:25
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
Definitions of particular quadrature rules on simplices.