Flow123d  master-f44eb46
radioactive_decay.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 radioactive_decay.cc
15  * @brief
16  */
17 
22 
23 #include "system/global_defs.h"
24 #include "mesh/mesh.h"
25 #include "input/factory.hh"
26 #include "input/accessors.hh"
27 
28 #include "armadillo"
29 
30 FLOW123D_FORCE_LINK_IN_CHILD(radioactiveDecay)
31 
32 
33 using namespace Input::Type;
34 
36  return Record("RadioactiveDecayProduct", "A record describing a product of a radioactive decay.")
37  .allow_auto_conversion("name")
39  "The name of the product.")
40  .declare_key("energy", Double(0.0), Default("0.0"),
41  "Not used at the moment! The released energy in MeV from the decay of the radionuclide into the product.")
42  .declare_key("branching_ratio", Double(0.0), Default("1.0"),
43  "The branching ratio of the product when there is more than one."
44  "Considering only one product, the default ratio 1.0 is used."
45  "Its value must be positive. Further, the branching ratios of all products are normalized"
46  "by their sum, so the sum then gives 1.0 (this also resolves possible rounding errors).")
47  .close();
48 }
49 
51  return Record("Decay", "A model of a radioactive decay.")
52  .declare_key("radionuclide", String(), Default::obligatory(),
53  "The name of the parent radionuclide.")
54  .declare_key("half_life", Double(0.0), Default::obligatory(),
55  "The half life of the parent radionuclide in seconds.")
57  "An array of the decay products (daughters).")
58  .close();
59 }
60 
62  return Record("RadioactiveDecay", "A model of a radioactive decay and possibly of a decay chain.")
68  "An array of radioactive decays.")
69  .close();
70 }
71 
73  Input::register_class< RadioactiveDecay, Mesh &, Input::Record >("RadioactiveDecay") +
75 
76 
77 
79  : FirstOrderReactionBase(init_mesh, in_rec)
80 {
81 }
82 
84 {
85 }
86 
87 
89 {
90  Input::Array decay_array = input_record_.val<Input::Array>("decays");
91 
92  substance_ids_.resize( decay_array.size() );
93  half_lives_.resize( decay_array.size() );
94  bifurcation_.resize( decay_array.size() );
95 
96  unsigned int idx;
97  int i_decay = 0;
98  for (Input::Iterator<Input::Record> dec_it = decay_array.begin<Input::Record>();
99  dec_it != decay_array.end(); ++dec_it, ++i_decay)
100  {
101  //read half-life
102  half_lives_[i_decay] = dec_it->val<double>("half_life");
103 
104  //read radionuclide name and products array
105  string radionuclide = dec_it->val<string>("radionuclide");
106  Input::Array product_array = dec_it->val<Input::Array>("products");
107 
108  //resizing according to product count
109  substance_ids_[i_decay].resize(product_array.size()+1);
110  bifurcation_[i_decay].resize(product_array.size());
111 
112  //Reading products - setting substance ids, branching ratios.
113  //iterating over products
114  unsigned int i_product = 0;
115  for (Input::Iterator<Input::Record> prod_it = product_array.begin<Input::Record>();
116  prod_it != product_array.end(); ++prod_it, ++i_product)
117  {
118  string product_name = prod_it->val<string>("name");
119  idx = find_subst_name(product_name);
120  if (idx < eq_data_base_->substances_.size())
121  substance_ids_[i_decay][i_product+1] = idx;
122  else THROW(ReactionTerm::ExcUnknownSubstance()
123  << ReactionTerm::EI_Substance(product_name)
124  << (*prod_it).ei_address());
125 
126  bifurcation_[i_decay][i_product] = prod_it->val<double>("branching_ratio");
127  }
128 
129  // set radionuclide substance index
130  idx = find_subst_name(radionuclide);
131  if (idx < eq_data_base_->substances_.size())
132  substance_ids_[i_decay][0] = idx;
133  else THROW(ReactionTerm::ExcUnknownSubstance()
134  << ReactionTerm::EI_Substance(radionuclide)
135  << (*dec_it).ei_address());
136 
137  //Normalization of branching ratios in bifurcation vector by its norm.
138  //sum:
139  double sum=0;
140  for(auto &b : bifurcation_[i_decay])
141  sum += b;
142  //Normalization:
143  for(auto &b : bifurcation_[i_decay])
144  b = b / sum;
145  }
146 }
147 
148 
150 {
151  // create decay matrix
153  unsigned int reactant_index, product_index; //global indices of the substances
154  double exponent; //temporary variable k = ln(2)/t_half
155  for (unsigned int i_decay = 0; i_decay < half_lives_.size(); i_decay++) {
156  reactant_index = substance_ids_[i_decay][0];
157  exponent = log(2) / half_lives_[i_decay];
158  reaction_matrix_(reactant_index, reactant_index) = -exponent;
159 
160  for (unsigned int i_product = 1; i_product < substance_ids_[i_decay].size(); ++i_product){
161  product_index = substance_ids_[i_decay][i_product];
162  reaction_matrix_(product_index, reactant_index) = exponent * bifurcation_[i_decay][i_product-1];
163  }
164  }
165 }
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)
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 const int registrar
Registrar of class to factory.
static const Input::Type::Record & get_input_type_single_decay()
Input record which defines particular decay step.
void initialize_from_input() override
Initializes private members of sorption from the input record.
static const Input::Type::Record & get_input_type()
Input record for class RadioactiveDecay.
void assemble_ode_matrix(void) override
Implements the assembly of the system matrix of the ODEs.
std::vector< double > half_lives_
Half-lives of the substances.
static const Input::Type::Record & get_input_type_product()
Input record for a product of a radioactive decay.
~RadioactiveDecay(void)
Destructor.
RadioactiveDecay(Mesh &mesh, Input::Record in_rec)
Constructor.
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.