Flow123d  master-f44eb46
first_order_reaction.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 first_order_reaction.cc
15  * @brief
16  */
17 
21 
22 #include "system/global_defs.h"
23 #include "mesh/mesh.h"
24 #include "input/factory.hh"
25 #include "input/accessors.hh"
26 
27 FLOW123D_FORCE_LINK_IN_CHILD(firstOrderReaction)
28 
29 
30 using namespace Input::Type;
31 
33  return Record("FirstOrderReactionReactant", "A record describing a reactant of a reaction.")
34  .allow_auto_conversion("name")
36  "The name of the reactant.")
37  //.declare_key("stoichiometric_coefficient", Integer(0.0), Default::optional(1.0)) //in future
38  .close();
39 }
40 
42  return Record("FirstOrderReactionProduct", "A record describing a product of a reaction.")
43  .allow_auto_conversion("name")
45  "The name of the product.")
46  //.declare_key("stoichiometric_coefficient", Integer(0.0), Default::optional(1.0)) //in future
47  .declare_key("branching_ratio", Double(0.0), Default("1.0"),
48  "The branching ratio of the product when there are more products.\n"
49  "The value must be positive. Further, the branching ratios of all products are normalized "
50  "in order to sum to one.\n"
51  "The default value 1.0, should only be used in the case of single product.")
52  .close();
53 }
54 
56  return Record("Reaction", "Describes a single first order chemical reaction.")
58  "An array of reactants. Do not use array, reactions with only one reactant (decays) are implemented at the moment!")
59  .declare_key("reaction_rate", Double(0.0), Default::obligatory(),
60  "The reaction rate coefficient of the first order reaction.")
62  "An array of products.")
63  .close();
64 }
65 
66 
68  return Record("FirstOrderReaction", "A model of first order chemical reactions (decompositions of a reactant into products).")
74  "An array of first order chemical reactions.")
75  .close();
76 }
77 
79  Input::register_class< FirstOrderReaction, Mesh &, Input::Record >("FirstOrderReaction") +
81 
82 
84  : FirstOrderReactionBase(init_mesh, in_rec)
85 {
86 }
87 
89 {
90 }
91 
93 {
94  // create decay matrix
96  unsigned int reactant_index, product_index; //global indices of the substances
97  double exponent; //temporary variable for k
98  for (unsigned int i_reaction = 0; i_reaction < reaction_rates_.size(); i_reaction++) {
99  reactant_index = substance_ids_[i_reaction][0];
100  exponent = reaction_rates_[i_reaction];
101  reaction_matrix_(reactant_index, reactant_index) = -exponent;
102 
103  for (unsigned int i_product = 1; i_product < substance_ids_[i_reaction].size(); ++i_product){
104  product_index = substance_ids_[i_reaction][i_product];
105  reaction_matrix_(product_index, reactant_index) = exponent * bifurcation_[i_reaction][i_product-1];
106  }
107  }
108 }
109 
110 
112 {
113  unsigned int idx; //temporary variable, indexing substances
114 
115  Input::Array reactions_array = input_record_.val<Input::Array>("reactions");
116 
117  substance_ids_.resize( reactions_array.size() );
118  reaction_rates_.resize( reactions_array.size() );
119  bifurcation_.resize( reactions_array.size() );
120 
121  int i_reaction=0;
122  for (Input::Iterator<Input::Record> react_it = reactions_array.begin<Input::Record>();
123  react_it != reactions_array.end(); ++react_it, ++i_reaction)
124  {
125  //read reaction rate
126  reaction_rates_[i_reaction] = react_it->val<double>("reaction_rate");
127 
128  //read reactant name, product names and branching ratios
129  Input::Array reactant_array = react_it->val<Input::Array>("reactants");
130  Input::Array product_array = react_it->val<Input::Array>("products");
131 
132  //resize substance_ids array
133  substance_ids_[i_reaction].resize( product_array.size()+1 );
134  bifurcation_[i_reaction].resize(product_array.size());
135 
136  if(reactant_array.size() != 1)
137  THROW( ExcTooManyReactants() );
138 
139  //take only one reactant
140  Input::Iterator<Input::Record> reactant_it = reactant_array.begin<Input::Record>();
141  {
142  string reactant_name = reactant_it->val<string>("name");
143  idx = find_subst_name(reactant_name);
144  if (idx < eq_data_base_->substances_.size())
145  substance_ids_[i_reaction][0] = idx;
146  else THROW(ReactionTerm::ExcUnknownSubstance()
147  << ReactionTerm::EI_Substance(reactant_name)
148  << (*reactant_it).ei_address());
149  }
150 
151  // set products
152  unsigned int i_product = 0;
153  for(Input::Iterator<Input::Record> product_it = product_array.begin<Input::Record>();
154  product_it != product_array.end(); ++product_it, i_product++)
155  {
156  string product_name = product_it->val<string>("name");
157  idx = find_subst_name(product_name);
158  if (idx < eq_data_base_->substances_.size())
159  substance_ids_[i_reaction][i_product+1] = idx;
160  else THROW(ReactionTerm::ExcUnknownSubstance()
161  << ReactionTerm::EI_Substance(product_name)
162  << product_array.ei_address());
163 
164  bifurcation_[i_reaction][i_product] = product_it->val<double>("branching_ratio");
165  }
166 
167 
168  //Normalization of branching ratios in bifurcation vector by its norm.
169  //sum:
170  double sum=0;
171  for(auto &b : bifurcation_[i_reaction])
172  sum += b;
173  //Normalization:
174  for(auto &b : bifurcation_[i_reaction])
175  b = b / sum;
176  }
177 }
Input::Record input_record_
Definition: equation.hh:242
Base class for linear reactions and decay chain.
unsigned int n_substances_
Number of all transported substances. It is the dimension of the reaction matrix.
std::vector< std::vector< unsigned int > > substance_ids_
arma::mat reaction_matrix_
Reaction matrix.
std::vector< std::vector< double > > bifurcation_
unsigned int find_subst_name(const std::string &name)
static const Input::Type::Record & get_input_type_product()
Input record for a product of a reaction.
static const Input::Type::Record & get_input_type_single_reaction()
Input record which defines particular reaction.
static const Input::Type::Record & get_input_type()
Input record for class FirstOrderReaction.
~FirstOrderReaction(void)
Destructor.
void initialize_from_input() override
Initializes private members of sorption from the input record.
static const int registrar
Registrar of class to factory.
std::vector< double > reaction_rates_
Vector of reaction rates of the transported substances.
FirstOrderReaction(Mesh &init_mesh, Input::Record in_rec)
Constructor.
void assemble_ode_matrix(void) override
Implements the assembly of the system matrix of the ODEs.
static const Input::Type::Record & get_input_type_reactant()
Input record for a reactant of a reaction.
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
Iterator< ValueType > begin() const
unsigned int size() const
IteratorBase end() const
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
const Ret val(const string &key) const
Class for declaration of inputs sequences.
Definition: type_base.hh:339
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Record type proxy class.
Definition: type_record.hh:182
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
virtual Record & allow_auto_conversion(const string &from_key)
Allows shorter input of the Record providing only value of the from_key given as the parameter.
Definition: type_record.cc:133
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
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:503
Class for declaration of the input data that are in string format.
Definition: type_base.hh:582
Definition: mesh.h:362
static Input::Type::Abstract & it_abstract_immobile_term()
static Input::Type::Abstract & it_abstract_term()
static Input::Type::Abstract & it_abstract_reaction()
static Input::Type::Abstract & it_abstract_mobile_term()
Global macros to enhance readability and debugging, general constants.
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:104
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
Class ReactionTerm is an abstract class representing reaction term in transport.