antioch-0.4.0
Functions
threebody_process_unit.C File Reference
#include <limits>
#include "antioch/vector_utils.h"
#include "antioch/reaction.h"
#include "antioch/threebody_reaction.h"

Go to the source code of this file.

Functions

template<typename Scalar >
int tester ()
 
int main ()
 

Function Documentation

int main ( )

Definition at line 232 of file threebody_process_unit.C.

233 {
234  return (tester<double>() ||
235  tester<long double>() ||
236  tester<float>());
237 }
template<typename Scalar >
int tester ( )

Definition at line 36 of file threebody_process_unit.C.

References Antioch::Reaction< CoeffType >::add_forward_rate(), Antioch::KineticsModel::ARRHENIUS, Antioch::KineticsModel::BERTHELOT, Antioch::KineticsModel::BHE, Antioch::ThreeBodyReaction< CoeffType >::compute_forward_rate_coefficient(), Antioch::ThreeBodyReaction< CoeffType >::compute_forward_rate_coefficient_and_derivatives(), Antioch::KineticsModel::D, Antioch::KineticsModel::HERCOURT_ESSEN, Antioch::KineticsModel::KOOIJ, std::pow(), Antioch::Reaction< CoeffType >::set_efficiency(), Antioch::set_zero(), and Antioch::KineticsModel::VANTHOFF.

37 {
38  using std::abs;
39  using std::exp;
40  using std::pow;
41 
42  const Scalar Cf = 1.4;
43  const Scalar Ea = 5.0;
44  const Scalar beta = 1.2;
45  const Scalar D = 2.5e-2;
46 
47  const std::string equation("A + B -> C + D");
48  const unsigned int n_species(4);
49 
50  int return_flag = 0;
51  std::vector<Scalar> mol_densities;
52  mol_densities.push_back(1e-2);
53  mol_densities.push_back(1e-2);
54  mol_densities.push_back(1e-2);
55  mol_densities.push_back(1e-2);
56  std::vector<Scalar> species_eff;
57  species_eff.push_back(0.4);
58  species_eff.push_back(0.54);
59  species_eff.push_back(0.5);
60  species_eff.push_back(1.0);
61  Scalar TB_term;
62  Antioch::set_zero(TB_term);
63  for(unsigned int i = 0; i < n_species; i++)
64  {
65  TB_term += species_eff[i] * mol_densities[i];
66  }
67 
68 
69  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 100;
70 
71  for(Scalar T = 300.1; T <= 2500.1; T += 10.)
72  {
73 
74  const Antioch::KineticsConditions<Scalar> conditions(T);
75 
76  for(unsigned int ikinmod = 0; ikinmod < 6; ikinmod++)
77  {
78 
79  Antioch::KineticsType<Scalar> *rate_kinetics(NULL);
80  Scalar rate_exact;
81  Scalar derive_exact;
82  std::vector<Scalar> derive_dX_exact;
83  derive_dX_exact.resize(n_species);
85 
86  switch(ikinmod)
87  {
88  case 0:
89  {
91  rate_kinetics = new Antioch::HercourtEssenRate<Scalar>(Cf,beta,1.);
92  rate_exact = Cf * pow(T,beta) * TB_term;
93  derive_exact = Cf * pow (T,beta) * beta/T * TB_term;
94  for(unsigned int i = 0; i < n_species; i++)
95  {
96  derive_dX_exact[i] = Cf * pow(T,beta) * species_eff[i];
97  }
98  break;
99  }
100  case 1:
101  {
103  rate_kinetics = new Antioch::BerthelotRate<Scalar>(Cf,D);
104  rate_exact = Cf * exp(D*T) * TB_term;
105  derive_exact = Cf * exp(D*T) * D * TB_term;
106  for(unsigned int i = 0; i < n_species; i++)
107  {
108  derive_dX_exact[i] = Cf * exp(D*T) * species_eff[i];
109  }
110  break;
111  }
112  case 2:
113  {
115  rate_kinetics = new Antioch::ArrheniusRate<Scalar>(Cf,Ea,1.);
116  rate_exact = Cf * exp(-Ea/T) * TB_term;
117  derive_exact = Cf * exp(-Ea/T) * Ea/pow(T,2) * TB_term;
118  for(unsigned int i = 0; i < n_species; i++)
119  {
120  derive_dX_exact[i] = Cf * exp(-Ea/T) * species_eff[i];
121  }
122  break;
123  }
124  case 3:
125  {
126  kin_mod = Antioch::KineticsModel::BHE;
127  rate_kinetics = new Antioch::BerthelotHercourtEssenRate<Scalar>(Cf,beta,D,1.);
128  rate_exact = Cf * pow(T,beta) * exp(D * T) * TB_term;
129  derive_exact = Cf * pow(T,beta) * exp(D * T) * (beta/T + D) * TB_term;
130  for(unsigned int i = 0; i < n_species; i++)
131  {
132  derive_dX_exact[i] = Cf * pow(T,beta) * exp(D * T) * species_eff[i];
133  }
134  break;
135  }
136  case 4:
137  {
139  rate_kinetics = new Antioch::KooijRate<Scalar>(Cf,beta,Ea,1.,1.);
140  rate_exact = Cf * pow(T,beta) * exp(-Ea/T) * TB_term;
141  derive_exact = Cf * pow(T,beta) * exp(-Ea/T) * (beta/T + Ea/pow(T,2)) * TB_term;
142  for(unsigned int i = 0; i < n_species; i++)
143  {
144  derive_dX_exact[i] = Cf * pow(T,beta) * exp(-Ea/T) * species_eff[i];
145  }
146  break;
147  }
148  case 5:
149  {
151  rate_kinetics = new Antioch::VantHoffRate<Scalar>(Cf,beta,Ea,D,1.,1.);
152  rate_exact = Cf * pow(T,beta) * exp(-Ea/T + D * T) * TB_term;
153  derive_exact = Cf * pow(T,beta) * exp(-Ea/T + D * T) * (D + beta/T + Ea/pow(T,2)) * TB_term;
154  for(unsigned int i = 0; i < n_species; i++)
155  {
156  derive_dX_exact[i] = Cf * pow(T,beta) * exp(-Ea/T + D * T) * species_eff[i];
157  }
158  break;
159  }
160  }
161 
162  Antioch::ThreeBodyReaction<Scalar> * TB_reaction = new Antioch::ThreeBodyReaction<Scalar>(n_species,equation,true,kin_mod);
163  TB_reaction->add_forward_rate(rate_kinetics);
164  for(unsigned int i = 0; i < n_species; i++)
165  {
166  TB_reaction->set_efficiency("",i,species_eff[i]);
167  }
168  Scalar rate1 = TB_reaction->compute_forward_rate_coefficient(mol_densities,conditions);
169  Scalar rate;
170  Scalar drate_dT;
171  std::vector<Scalar> drate_dx;
172  drate_dx.resize(n_species);
173  TB_reaction->compute_forward_rate_coefficient_and_derivatives(mol_densities,conditions,rate,drate_dT,drate_dx);
174 
175 
176  for(unsigned int i = 0; i < n_species; i++)
177  {
178  if( abs( (drate_dx[i] - derive_dX_exact[i])/derive_dX_exact[i] ) > tol )
179  {
180  std::cout << std::scientific << std::setprecision(16)
181  << "Error: Mismatch in derivative values." << std::endl
182  << "Kinetics model (see enum) " << kin_mod << std::endl
183  << "species " << i << std::endl
184  << "species molar densities " << mol_densities[i] << std::endl
185  << "species efficiency " << species_eff[i] << std::endl
186  << "drate_dX = " << drate_dx[i] << std::endl
187  << "drate_dX_exact = " << derive_dX_exact[i] << std::endl;
188  return_flag = 1;
189  }
190  }
191  if( abs( (rate1 - rate_exact)/rate_exact ) > tol )
192  {
193  std::cout << std::scientific << std::setprecision(16)
194  << "Error: Mismatch in rate values." << std::endl
195  << "Kinetics model (see enum) " << kin_mod << std::endl
196  << "T = " << T << " K" << std::endl
197  << "rate(T) = " << rate1 << std::endl
198  << "rate_exact = " << rate_exact << std::endl;
199 
200  return_flag = 1;
201  }
202  if( abs( (rate - rate_exact)/rate_exact ) > tol )
203  {
204  std::cout << std::scientific << std::setprecision(16)
205  << "Error: Mismatch in rate values." << std::endl
206  << "Kinetics model (see enum) " << kin_mod << std::endl
207  << "T = " << T << " K" << std::endl
208  << "rate(T) = " << rate << std::endl
209  << "rate_exact = " << rate_exact << std::endl;
210 
211  return_flag = 1;
212  }
213  if( abs( (drate_dT - derive_exact)/derive_exact ) > tol )
214  {
215  std::cout << std::scientific << std::setprecision(16)
216  << "Error: Mismatch in rate derivative values." << std::endl
217  << "Kinetics model (see enum) " << kin_mod << std::endl
218  << "T = " << T << " K" << std::endl
219  << "drate_dT(T) = " << drate_dT << std::endl
220  << "derive_exact = " << derive_exact << std::endl;
221 
222  return_flag = 1;
223  }
224  delete TB_reaction;
225  if(return_flag)return return_flag;
226  }
227  }
228 
229  return return_flag;
230 }
Berthelot rate equation.
base class for kinetics models
Definition: kinetics_type.h:82
Berthelot Hercourt-Essen rate equation.
A single reaction mechanism.
Definition: reaction.h:64
Antioch::enable_if_c< Antioch::is_valarray< T >::value, typename Antioch::state_type< T >::type >::type pow(const T &in, const T2 &n)
void add_forward_rate(KineticsType< CoeffType, VectorCoeffType > *rate)
Add a forward rate object.
void compute_forward_rate_coefficient_and_derivatives(const VectorStateType &molar_densities, const KineticsConditions< StateType, VectorStateType > &conditions, StateType &kfwd, StateType &dkfwd_dT, VectorStateType &dkfwd_dX) const
void set_zero(_Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &a)
Definition: eigen_utils.h:217
Van't Hoff rate equation.
Definition: kinetics_type.h:62
Hercourt-Essen rate equation.
Kooij rate equation.
Definition: kinetics_type.h:59
StateType compute_forward_rate_coefficient(const VectorStateType &molar_densities, const KineticsConditions< StateType, VectorStateType > &conditions) const
void set_efficiency(const std::string &, const unsigned int s, const CoeffType efficiency)
This class contains the conditions of the chemistry.

Generated on Thu Jul 7 2016 11:09:48 for antioch-0.4.0 by  doxygen 1.8.8