antioch-0.4.0
reaction_set.h
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // Antioch - A Gas Dynamics Thermochemistry Library
5 //
6 // Copyright (C) 2014-2016 Paul T. Bauman, Benjamin S. Kirk,
7 // Sylvain Plessis, Roy H. Stonger
8 //
9 // Copyright (C) 2013 The PECOS Development Team
10 //
11 // This library is free software; you can redistribute it and/or
12 // modify it under the terms of the Version 2.1 GNU Lesser General
13 // Public License as published by the Free Software Foundation.
14 //
15 // This library is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
23 // Boston, MA 02110-1301 USA
24 //
25 //-----------------------------------------------------------------------el-
26 
27 
28 #ifndef ANTIOCH_REACTION_SET_H
29 #define ANTIOCH_REACTION_SET_H
30 
31 // Antioch
32 #include "antioch/chemical_mixture.h"
33 #include "antioch/reaction.h"
34 #include "antioch/kinetics_conditions.h"
35 #include "antioch/elementary_reaction.h"
36 #include "antioch/duplicate_reaction.h"
37 #include "antioch/threebody_reaction.h"
38 #include "antioch/falloff_reaction.h"
39 #include "antioch/falloff_threebody_reaction.h"
40 #include "antioch/lindemann_falloff.h"
41 #include "antioch/troe_falloff.h"
42 #include "antioch/string_utils.h"
43 
44 // C++
45 #include <iostream>
46 #include <fstream>
47 #include <iomanip>
48 #include <vector>
49 #include <limits>
50 
51 namespace Antioch
52 {
53 
58  template<typename CoeffType=double>
59  class ReactionSet
60  {
61 
62  public:
63 
65  ReactionSet( const ChemicalMixture<CoeffType>& chem_mixture );
66 
67  ~ReactionSet();
68 
70  unsigned int n_species() const;
71 
73  unsigned int n_reactions() const;
74 
76  //
77  // The ownership is transfered to the ReactionSet
78  // object, they are deleted in the destructor.
79  void add_reaction(Reaction<CoeffType>* reaction);
80 
82  //
83  // The corresponding pointer is deleted
84  void remove_reaction(unsigned int nr);
85 
87  const Reaction<CoeffType>& reaction(const unsigned int r) const;
88 
90  Reaction<CoeffType>& reaction(const unsigned int r);
91 
93  unsigned int reaction_by_id(const std::string & reaction_id) const;
94 
96  //
97  // in charge of the human-to-antioch translation
98  template <typename ParamType>
99  void set_parameter_of_reaction(const std::string & reaction_id, const std::vector<std::string> & keywords, ParamType value);
100 
102  //
103  // in charge of the human-to-antioch translation
104  CoeffType get_parameter_of_reaction(const std::string & reaction_id, const std::vector<std::string> & keywords) const;
105 
106  const ChemicalMixture<CoeffType>& chemical_mixture() const;
107 
109  template <typename StateType, typename VectorStateType, typename VectorReactionsType>
110  void compute_reaction_rates( const KineticsConditions<StateType,VectorStateType>& conditions,
111  const VectorStateType& molar_densities,
112  const VectorStateType& h_RT_minus_s_R,
113  VectorReactionsType& net_reaction_rates ) const;
114 
116  template <typename StateType, typename VectorStateType, typename VectorReactionsType, typename MatrixReactionsType>
117  void compute_reaction_rates_and_derivs( const KineticsConditions<StateType,VectorStateType>& conditions,
118  const VectorStateType& molar_densities,
119  const VectorStateType& h_RT_minus_s_R,
120  const VectorStateType& dh_RT_minus_s_R_dT,
121  VectorReactionsType& net_reaction_rates,
122  VectorReactionsType& dnet_rate_dT,
123  MatrixReactionsType& dnet_rate_dX_s ) const;
124 
126  template <typename StateType, typename VectorStateType>
127  void print_chemical_scheme( std::ostream& output,
128  const KineticsConditions<StateType,VectorStateType>& conditions,
129  const VectorStateType& molar_densities,
130  const VectorStateType& h_RT_minus_s_R ,
131  std::vector<VectorStateType> &lossMatrix,
132  std::vector<VectorStateType> &prodMatrix,
133  std::vector<VectorStateType> &netMatrix) const;
134 
136  template<typename StateType, typename VectorStateType>
137  void get_reactive_scheme( const KineticsConditions<StateType,VectorStateType>& conditions,
138  const VectorStateType& molar_densities,
139  const VectorStateType& h_RT_minus_s_R,
140  VectorStateType& net_rates,
141  VectorStateType& kfwd_const,
142  VectorStateType& kbkwd_const,
143  VectorStateType& kfwd,
144  VectorStateType& kbkwd,
145  VectorStateType& fwd_conc,
146  VectorStateType& bkwd_conc) const;
147 
148 
150  void print( std::ostream& os = std::cout ) const;
151 
153  friend std::ostream& operator<<( std::ostream& os, const ReactionSet<CoeffType>& rset )
154  {
155  rset.print(os);
156  return os;
157  }
158 
159  private:
160 
161  ReactionSet();
162 
164  //
165  // It finds the kinetics model parameter given keywords and reaction:
166  // - nr is the index of the rate to change
167  // - unit is the unit if provided
168  //
169  // This function is used for both getter and setter.
170  void find_kinetics_model_parameter(const unsigned int r,const std::vector<std::string> & keywords, unsigned int & nr, std::string & unit, int & l) const;
171 
173  //
174  // It finds the chemical process parameter given keywords and reaction type:
175  // - species is the species whose efficiency we want to change
176  //
177  // If we want to change anything other than an efficiency, well nothing happens here.
178  // This function is used for both getter and setter.
179  void find_chemical_process_parameter(ReactionType::Parameters paramChem ,const std::vector<std::string> & keywords, unsigned int & species) const;
180 
181  const ChemicalMixture<CoeffType>& _chem_mixture;
182 
183  std::vector<Reaction<CoeffType>* > _reactions;
184 
186  const CoeffType _P0_R;
187 
188  };
189 
190  /* ------------------------- Inline Functions -------------------------*/
191  template<typename CoeffType>
192  inline
193  unsigned int ReactionSet<CoeffType>::n_species() const
194  {
195  return _chem_mixture.n_species();
196  }
197 
198  template<typename CoeffType>
199  inline
200  unsigned int ReactionSet<CoeffType>::n_reactions() const
201  {
202  return _reactions.size();
203  }
204 
205  template<typename CoeffType>
206  inline
207  void ReactionSet<CoeffType>::add_reaction(Reaction<CoeffType>* reaction)
208  {
209  _reactions.push_back(reaction);
210 
211  // and make sure it is initialized!
212  _reactions.back()->initialize(_reactions.size() - 1);
213 
214  return;
215  }
216 
217  template<typename CoeffType>
218  inline
219  void ReactionSet<CoeffType>::remove_reaction(unsigned int nr)
220  {
221  antioch_assert_less(nr,_reactions.size());
222 
223  //first clear the memory
224  delete _reactions[nr];
225 
226  //second, release the spot
227  _reactions.erase(_reactions.begin() + nr);
228  }
229 
230  template<typename CoeffType>
231  inline
232  const Reaction<CoeffType>& ReactionSet<CoeffType>::reaction(const unsigned int r) const
233  {
234  antioch_assert_less(r, this->n_reactions());
235  return *_reactions[r];
236  }
237 
238  template<typename CoeffType>
239  inline
240  Reaction<CoeffType>& ReactionSet<CoeffType>::reaction(const unsigned int r)
241  {
242  antioch_assert_less(r, this->n_reactions());
243  return *_reactions[r];
244  }
245 
246  template<typename CoeffType>
247  inline
248  const ChemicalMixture<CoeffType>& ReactionSet<CoeffType>::chemical_mixture() const
249  {
250  return _chem_mixture;
251  }
252 
253 
254  template<typename CoeffType>
255  inline
256  ReactionSet<CoeffType>::ReactionSet( const ChemicalMixture<CoeffType>& chem_mixture )
257  : _chem_mixture(chem_mixture),
258  _P0_R(1.0e5/Constants::R_universal<CoeffType>()) //SI
259  {
260  return;
261  }
262 
263 
264  template<typename CoeffType>
265  inline
266  ReactionSet<CoeffType>::~ReactionSet()
267  {
268  for(unsigned int ir = 0; ir < _reactions.size(); ir++)delete _reactions[ir];
269  return;
270  }
271 
272  template<typename CoeffType>
273  template<typename StateType, typename VectorStateType, typename VectorReactionsType>
274  inline
275  void ReactionSet<CoeffType>::compute_reaction_rates ( const KineticsConditions<StateType,VectorStateType>& conditions,
276  const VectorStateType& molar_densities,
277  const VectorStateType& h_RT_minus_s_R,
278  VectorReactionsType& net_reaction_rates ) const
279  {
280  antioch_assert_equal_to( net_reaction_rates.size(), this->n_reactions() );
281 
283  // antioch_assert_greater(T, 0.0);
284 
285  antioch_assert_equal_to( molar_densities.size(), this->n_species() );
286  antioch_assert_equal_to( h_RT_minus_s_R.size(), this->n_species() );
287 
288  // useful constants
289  const StateType P0_RT = _P0_R/conditions.T(); // used to transform equilibrium constant from pressure units
290 
291  // compute reaction forward rates & other reaction-sized arrays
292  for (unsigned int rxn=0; rxn<this->n_reactions(); rxn++)
293  {
294  net_reaction_rates[rxn] = this->reaction(rxn).compute_rate_of_progress(molar_densities, conditions, P0_RT, h_RT_minus_s_R);
295  }
296 
297  return;
298  }
299 
300  template<typename CoeffType>
301  template<typename StateType, typename VectorStateType, typename VectorReactionsType, typename MatrixReactionsType>
302  inline
303  void ReactionSet<CoeffType>::compute_reaction_rates_and_derivs( const KineticsConditions<StateType,VectorStateType>& conditions,
304  const VectorStateType& molar_densities,
305  const VectorStateType& h_RT_minus_s_R,
306  const VectorStateType& dh_RT_minus_s_R_dT,
307  VectorReactionsType& net_reaction_rates,
308  VectorReactionsType& dnet_rate_dT,
309  MatrixReactionsType& dnet_rate_dX_s ) const
310  {
311  antioch_assert_equal_to( net_reaction_rates.size(), this->n_reactions() );
312  antioch_assert_equal_to( dnet_rate_dT.size(), this->n_reactions() );
313  antioch_assert_equal_to( dnet_rate_dX_s.size(), this->n_reactions() );
314 #ifdef NDEBUG
315 #else
316  for (unsigned int r=0; r < this->n_reactions(); r++)
317  {
318  antioch_assert_equal_to( dnet_rate_dX_s[r].size(), this->n_species() );
319  }
320 #endif
321 
323  // antioch_assert_greater(T, 0.0);
324 
325  antioch_assert_equal_to( molar_densities.size(), this->n_species() );
326  antioch_assert_equal_to( h_RT_minus_s_R.size(), this->n_species() );
327 
328  // useful constants
329  const StateType P0_RT = _P0_R/conditions.T(); // used to transform equilibrium constant from pressure units
330 
331  // compute reaction forward rates & other reaction-sized arrays
332  for (unsigned int rxn=0; rxn<this->n_reactions(); rxn++)
333  {
334  this->reaction(rxn).compute_rate_of_progress_and_derivatives( molar_densities, _chem_mixture,
335  conditions, P0_RT, h_RT_minus_s_R, dh_RT_minus_s_R_dT,
336  net_reaction_rates[rxn],
337  dnet_rate_dT[rxn],
338  dnet_rate_dX_s[rxn] );
339  }
340 
341  return;
342  }
343 
344 
345  template<typename CoeffType>
346  inline
347  unsigned int ReactionSet<CoeffType>::reaction_by_id(const std::string & reaction_id) const
348  {
349  unsigned int r(0);
350  for(r = 0; r < this->n_reactions(); r++)
351  {
352  if(this->reaction(r).id() == reaction_id)
353  break;
354  }
355  if(r >= this->n_reactions())
356  {
357  std::string errmsg = "Error: did not find reaction \"" + reaction_id + "\"\nIds are: ";
358  for(r = 0; r < this->n_reactions(); r++)
359  {
360  errmsg += this->reaction(r).id() + ", ";
361  if(r%10 == 0)errmsg += "\n"; // a few formatting is nice
362  }
363  antioch_warning(errmsg);
364  antioch_error();
365  }
366 
367  return r;
368  }
369 
370 
371 
372  template<typename CoeffType>
373  inline
374  void ReactionSet<CoeffType>::find_kinetics_model_parameter(const unsigned int r,const std::vector<std::string> & keywords, unsigned int & nr, std::string & unit, int & l) const
375  {
376 // now we need to know a few things:
377 // if we are in a falloff or duplicate, next keyword is which rate we want, then unit
378 // if we are in an elementary photochemistry, next keyword is which index of sigma or lamba we want, then unit
379 // else we may have a unit
380  switch(this->reaction(r).type())
381  {
382  case ReactionType::DUPLICATE:
383  {
384  nr = std::stoi(keywords[1]); // C++11, throws an exception on error
385  if(keywords.size() > 2)unit = keywords[2];
386  }
387  break;
388  case ReactionType::LINDEMANN_FALLOFF:
389  case ReactionType::TROE_FALLOFF:
390  case ReactionType::LINDEMANN_FALLOFF_THREE_BODY:
391  case ReactionType::TROE_FALLOFF_THREE_BODY:
392  {
393  switch(string_to_kin_enum(keywords[1])) // falloff, if here by error, NOT_FOUND is get (0)
394  {
395  // This is hard-coded (because we want a vector and not a map)
396  // low pressure is the first, high pressure the second
397  // see FalloffReaction object
398  case KineticsModel::Parameters::LOW_PRESSURE:
399  {
400  nr = 0;
401  }
402  break;
403  case KineticsModel::Parameters::HIGH_PRESSURE:
404  {
405  nr = 1;
406  }
407  break;
408  default: // ?
409  {
410  antioch_error();
411  }
412  break;
413  }
414  if(keywords.size() > 2)unit = keywords[2]; // unit baby!
415  }
416  break;
417  case ReactionType::ELEMENTARY: // photochem only
418  {
419  if(this->reaction(r).kinetics_model() == KineticsModel::PHOTOCHEM)
420  {
421  l = std::stoi(keywords[1]); // C++11, throws an exception on error
422  if(keywords.size() > 2)unit = keywords[2];
423  }
424  }
425  break;
426  default:
427  {
428  if(keywords.size() > 1)unit = keywords[1]; // unit baby!
429  }
430  break;
431  }
432 
433  return;
434  }
435 
436  template<typename CoeffType>
437  inline
438  void ReactionSet<CoeffType>::find_chemical_process_parameter(ReactionType::Parameters paramChem ,const std::vector<std::string> & keywords, unsigned int & species) const
439  {
440 // there's not really a unit issue here:
441 // * efficiencies: unitless
442 // * Troe alpha: unitless
443 // * Troe T*: temperature, if it's other than K, you're a bad (like really really bad) scientist, get off my library
444 // * Troe T**: temperature, see above
445 // * Troe T***: temperature, do I really need to say anything?
446 
447  if(paramChem == ReactionType::Parameters::EFFICIENCIES) // who?
448  {
449  antioch_assert_greater(keywords.size(),1); // we need a name
450 
451  if(!this->chemical_mixture().species_name_map().count(keywords[1]))
452  antioch_error(); //who's this?
453 
454  species = this->chemical_mixture().species_name_map().at(keywords[1]);
455  }
456 
457  return;
458  }
459 
460  template<typename CoeffType>
461  inline
462  CoeffType ReactionSet<CoeffType>::get_parameter_of_reaction(const std::string & reaction_id, const std::vector<std::string> & keywords) const
463  {
464  antioch_assert(keywords.size()); // not zero
465 
466  CoeffType parameter;
467  // when CoeffType is vectorizable, we want
468  // to have as little rewriting as possible
469  set_zero(parameter);
470 
471  unsigned int r = this->reaction_by_id(reaction_id);
472 
473  // 1 parse high level
474  KineticsModel::Parameters paramKin = string_to_kin_enum(keywords[0]);
475  ReactionType::Parameters paramChem = string_to_chem_enum(keywords[0]);
476 
477 // provide the necessary enum,
478 // index of reaction rate if kinetics
479 // index of species if chemical
480  if(paramKin != KineticsModel::Parameters::NOT_FOUND)
481  {
482  // which rate? Duplicate want an unsigned int, falloff a keyword
483  unsigned int nr(0); // default is one rate
484  std::string unit("SI"); // default internal parameter unit system
485  int l(-1); // if vectorized parameter
486 
487  this->find_kinetics_model_parameter(r,keywords,nr,unit,l);
488 
489  // let's get this
490  parameter = reaction(r).get_parameter_of_rate(paramKin, nr, unit);
491 
492  }else if(paramChem != ReactionType::Parameters::NOT_FOUND)
493  {
494  unsigned int species = std::numeric_limits<unsigned int>::max(); // sensible default
495 
496  this->find_chemical_process_parameter(paramChem, keywords, species);
497 
498  // let's get this
499  parameter = reaction(r).get_parameter_of_chemical_process(paramChem, species);
500  }else
501  {
502  antioch_error();
503  }
504 
505  return parameter;
506 
507  }
508 
509  template<typename CoeffType>
510  template <typename ParamType>
511  inline
512  void ReactionSet<CoeffType>::set_parameter_of_reaction(const std::string & reaction_id, const std::vector<std::string> & keywords, ParamType value)
513  {
514  antioch_assert(keywords.size()); // not zero
515 
516  // 1 find the reaction
517  unsigned int r = this->reaction_by_id(reaction_id);
518 
519  // 1 parse high level
520  KineticsModel::Parameters paramKin = string_to_kin_enum(keywords[0]);
521  ReactionType::Parameters paramChem = string_to_chem_enum(keywords[0]);
522 // provide the necessary enum,
523 // index of reaction rate if kinetics
524 // index of species if chemical
525  if(paramKin != KineticsModel::Parameters::NOT_FOUND)
526  {
527  // which rate? Duplicate want an unsigned int, falloff a keyword
528  unsigned int nr(0); // default is one rate
529  std::string unit("SI"); // default internal parameter unit system
530  int l(-1); // for vectorized parameter
531 
532  this->find_kinetics_model_parameter(r,keywords,nr,unit,l);
533 
534  // let's set this
535  (l < 0)?this->reaction(r).set_parameter_of_rate(paramKin, value, nr, unit):
536  this->reaction(r).set_parameter_of_rate(paramKin, value, nr, l, unit);
537 
538  }else if(paramChem != ReactionType::Parameters::NOT_FOUND)
539  {
540  unsigned int species = std::numeric_limits<unsigned int>::max(); // sensible default
541 
542  this->find_chemical_process_parameter(paramChem, keywords, species);
543 
544  // let's set this
545  this->reaction(r).set_parameter_of_chemical_process(paramChem, value, species);
546  }else
547  {
548  antioch_error();
549  }
550 
551  }
552 
553  template<typename CoeffType>
554  template<typename StateType, typename VectorStateType>
555  inline
556  void ReactionSet<CoeffType>::print_chemical_scheme( std::ostream& output,
557  const KineticsConditions<StateType,VectorStateType>& conditions,
558  const VectorStateType& molar_densities,
559  const VectorStateType& h_RT_minus_s_R,
560  std::vector<VectorStateType>& lossMatrix,
561  std::vector<VectorStateType>& prodMatrix,
562  std::vector<VectorStateType>& netMatrix ) const
563  {
564 
565  //filling matrixes
566  VectorStateType netRate,kfwd_const,kbkwd_const,kfwd,kbkwd,fwd_conc,bkwd_conc;
567 
568  //getting reaction infos
569  get_reactive_scheme(conditions,molar_densities,h_RT_minus_s_R,netRate,kfwd_const,kbkwd_const,kfwd,kbkwd,fwd_conc,bkwd_conc);
570 
571  lossMatrix.resize(this->n_species());
572  prodMatrix.resize(this->n_species());
573  netMatrix.resize(this->n_species());
574 
575  for(unsigned int s = 0; s < this->n_species(); s++)
576  {
577  lossMatrix[s].resize(this->n_reactions(),0.);
578  prodMatrix[s].resize(this->n_reactions(),0.);
579  netMatrix[s].resize(this->n_reactions(),0.);
580  }
581 
582 
583  //filling matrixes
584  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
585  {
586  const Reaction<CoeffType>& reaction = this->reaction(rxn);
587  for (unsigned int r=0; r<reaction.n_reactants(); r++)
588  {
589  lossMatrix[reaction.reactant_id(r)][rxn] += - static_cast<CoeffType>(reaction.reactant_stoichiometric_coefficient(r)) * kfwd[rxn];
590  prodMatrix[reaction.reactant_id(r)][rxn] += static_cast<CoeffType>(reaction.reactant_stoichiometric_coefficient(r)) * kbkwd[rxn];
591  netMatrix[reaction.reactant_id(r)][rxn] += lossMatrix[reaction.reactant_id(r)][rxn] + prodMatrix[reaction.reactant_id(r)][rxn];
592  }
593  for (unsigned int p=0; p<reaction.n_products(); p++)
594  {
595  lossMatrix[reaction.product_id(p)][rxn] += - static_cast<CoeffType>(reaction.product_stoichiometric_coefficient(p)) * kbkwd[rxn];
596  prodMatrix[reaction.product_id(p)][rxn] += static_cast<CoeffType>(reaction.product_stoichiometric_coefficient(p)) * kfwd[rxn];
597  netMatrix[reaction.product_id(p)][rxn] += lossMatrix[reaction.product_id(p)][rxn] + prodMatrix[reaction.product_id(p)][rxn];
598  }
599  }
600 
601  //explanation of header
602  output << "# molar units considered for those budgets (kmol, m3, s)" << std::endl;
603  output << "# formatted as follow:" << std::endl;
604  output << "# rows: species, columns: reactions" << std::endl;
605  output << "#" << std::endl;
606  output << "# | equation" << std::endl;
607  output << "# | rate coefficient forward , rate coefficient backward" << std::endl;
608  output << "# | forward concentrations , backward concentrations" << std::endl;
609  output << "# | rate forward , rate backward" << std::endl;
610  output << "# ------------------------------------------------------------------" << std::endl;
611  output << "# species | production term , loss term" << std::endl;
612  output << "# species | net production/loss term" << std::endl;
613  output << "#" << std::endl << std::endl;
614 
615 
616  //header
617  // // equation
618  output << std::setw(10) << " ";
619  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
620  {
621  output << std::setw(31) << this->_reactions[rxn]->equation();
622  }
623  output << std::endl;
624  // // coeffs
625  output << std::setw(10) << " ";
626  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
627  {
628  output << std::setw(15) << kfwd_const[rxn] << "," << std::setw(15) << kbkwd_const[rxn];
629  }
630  output << std::endl;
631  // // conc
632  output << std::setw(10) << " ";
633  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
634  {
635  output << std::setw(15) << fwd_conc[rxn] << "," << std::setw(15) << bkwd_conc[rxn];
636  }
637  output << std::endl;
638  // // rates
639  output << std::setw(10) << " ";
640  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
641  {
642  output << std::setw(15) << kfwd[rxn] << "," << std::setw(15) << kbkwd[rxn];
643  }
644  output << std::endl;
645  output << "----------------------------" << std::endl;
646 
647  //species
648 
649  for(unsigned int isp = 0; isp < this->n_species(); isp++)
650  {
651  output << std::setw(10) << _chem_mixture.chemical_species()[isp]->species();
652  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
653  {
654  output << std::setw(15) << std::scientific << std::setprecision(6) << prodMatrix[isp][rxn] << ","
655  << std::setw(15) << std::scientific << std::setprecision(6) << lossMatrix[isp][rxn];
656  }
657  output << std::endl;
658  // // net
659  output << std::setw(10) << "";
660  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
661  {
662  output << std::setw(31) << std::scientific << std::setprecision(16) << netMatrix[isp][rxn];
663  }
664  output << std::endl;
665  }
666  output << std::endl;
667  output << std::endl;
668  output << std::endl;
669 
670  output << "# production table" << std::endl << std::endl;
671  //header
672  // // equation
673  output << std::setw(10) << " ";
674  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
675  {
676  output << std::setw(31) << this->_reactions[rxn]->equation();
677  }
678  output << std::setw(31) << "Total" << std::endl;
679  for(unsigned int isp = 0; isp < this->n_species(); isp++)
680  {
681  output << std::setw(10) << _chem_mixture.chemical_species()[isp]->species();
682  CoeffType rowTotal(0.);
683  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
684  {
685  output << std::setw(31) << std::scientific << std::setprecision(16) << prodMatrix[isp][rxn];
686  rowTotal += prodMatrix[isp][rxn];
687  }
688  output << std::setw(31) << std::scientific << std::setprecision(16) << rowTotal << std::endl;
689  }
690  output << std::endl << std::endl;
691 
692  output << "# loss table" << std::endl << std::endl;
693  //header
694  // // equation
695  output << std::setw(10) << " ";
696  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
697  {
698  output << std::setw(31) << this->_reactions[rxn]->equation();
699  }
700  output << std::setw(31) << "Total" << std::endl;
701  for(unsigned int isp = 0; isp < this->n_species(); isp++)
702  {
703  output << std::setw(10) << _chem_mixture.chemical_species()[isp]->species();
704  CoeffType rowTotal(0.);
705  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
706  {
707  output << std::setw(31) << std::scientific << std::setprecision(16) << lossMatrix[isp][rxn];
708  rowTotal += lossMatrix[isp][rxn];
709  }
710  output << std::setw(31) << std::scientific << std::setprecision(16) << rowTotal << std::endl;
711  }
712  output << std::endl << std::endl;
713 
714  output << "# net table" << std::endl << std::endl;
715  //header
716  // // equation
717  output << std::setw(10) << " ";
718  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
719  {
720  output << std::setw(31) << this->_reactions[rxn]->equation();
721  }
722  output << std::setw(31) << "Total" << std::endl;
723 
724  CoeffType columnTotal(0.);
725  std::vector<CoeffType> columnSum;
726  columnSum.resize(this->n_reactions(),0.);
727  CoeffType netTotalRow(0.);
728 
729  for(unsigned int isp = 0; isp < this->n_species(); isp++)
730  {
731  output << std::setw(10) << _chem_mixture.chemical_species()[isp]->species();
732  CoeffType rowTotal(0.);
733  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
734  {
735  output << std::setw(31) << std::scientific << std::setprecision(16) << netMatrix[isp][rxn];
736  rowTotal += netMatrix[isp][rxn];
737  columnSum[rxn] += netMatrix[isp][rxn];
738  }
739  output << std::setw(31) << std::scientific << std::setprecision(16) << rowTotal << std::endl;
740  netTotalRow += rowTotal;
741  }
742  output << std::setw(10) << "Total";
743  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
744  {
745  output << std::setw(31) << std::scientific << std::setprecision(16) << columnSum[rxn];
746  columnTotal += columnSum[rxn];
747  }
748  output << std::endl << std::endl;
749  output << "sum of row sums: " << std::scientific << std::setprecision(16) << netTotalRow << std::endl;
750  output << "sum of column sums: " << std::scientific << std::setprecision(16) << columnTotal << std::endl << std::endl;
751  output << "# net table, mass here (kg, m3, s)" << std::endl << std::endl;
752  //header
753  // // equation
754  output << std::setw(10) << " ";
755  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
756  {
757  output << std::setw(31) << this->_reactions[rxn]->equation();
758  }
759  output << std::setw(31) << "Total" << std::endl;
760 
761  columnTotal = 0.;
762  columnSum.resize(this->n_reactions(),0.);
763  netTotalRow = 0.;
764 
765  for(unsigned int isp = 0; isp < this->n_species(); isp++)
766  {
767  output << std::setw(10) << _chem_mixture.chemical_species()[isp]->species();
768  CoeffType rowTotal(0.);
769  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
770  {
771  output << std::setw(31) << std::scientific << std::setprecision(16) << netMatrix[isp][rxn];
772  rowTotal += netMatrix[isp][rxn] * _chem_mixture.M(isp);
773  columnSum[rxn] += netMatrix[isp][rxn] * _chem_mixture.M(isp);
774  }
775  output << std::setw(31) << std::scientific << std::setprecision(16) << rowTotal << std::endl;
776  netTotalRow += rowTotal;
777  }
778  output << std::setw(10) << "Total";
779  for(unsigned int rxn = 0; rxn < this->n_reactions(); rxn++)
780  {
781  output << std::setw(31) << std::scientific << std::setprecision(16) << columnSum[rxn];
782  columnTotal += columnSum[rxn];
783  }
784  output << std::endl << std::endl;
785 
786  output << std::endl << std::endl;
787  output << "# Mixture informations" << std::endl;
788  output << "# species / concentration / molar mass" << std::endl;
789  for(unsigned int isp = 0; isp < this->n_species(); isp++)
790  {
791  output << std::setw(10) << _chem_mixture.chemical_species()[isp]->species()
792  << std::setw(31) << molar_densities[isp]
793  << std::setw(31) << _chem_mixture.M(isp) << std::endl;
794  }
795 
796  return;
797  }
798 
799  template<typename CoeffType>
800  template<typename StateType, typename VectorStateType>
801  inline
802  void ReactionSet<CoeffType>::get_reactive_scheme( const KineticsConditions<StateType,VectorStateType>& conditions,
803  const VectorStateType& molar_densities,
804  const VectorStateType& h_RT_minus_s_R,
805  VectorStateType& net_rates,
806  VectorStateType& kfwd_const,
807  VectorStateType& kbkwd_const,
808  VectorStateType& kfwd,
809  VectorStateType& kbkwd,
810  VectorStateType& fwd_conc,
811  VectorStateType& bkwd_conc) const
812  {
814  // antioch_assert_greater(T, 0.0);
815 
816  antioch_assert_equal_to( molar_densities.size(), this->n_species() );
817  antioch_assert_equal_to( h_RT_minus_s_R.size(), this->n_species() );
818 
819  // useful constants
820  const StateType P0_RT = _P0_R/conditions.T(); // used to transform equilibrium constant from pressure units
821 
822  net_rates.resize(this->n_reactions(),0);
823  kfwd_const.resize(this->n_reactions(),0);
824  kfwd.resize(this->n_reactions(),0);
825  kbkwd_const.resize(this->n_reactions(),0);
826  kbkwd.resize(this->n_reactions(),0);
827  fwd_conc.resize(this->n_reactions(),1);
828  bkwd_conc.resize(this->n_reactions(),1);
829 
830  // compute reaction forward rates & other reaction-sized arrays
831  for (unsigned int rxn=0; rxn<this->n_reactions(); rxn++)
832  {
833  const Reaction<CoeffType>& reaction = this->reaction(rxn);
834  kfwd_const[rxn] = reaction.compute_forward_rate_coefficient(molar_densities,conditions);
835  kfwd[rxn] = kfwd_const[rxn];
836 
837  for (unsigned int r=0; r<reaction.n_reactants(); r++)
838  {
839  fwd_conc[rxn] *= pow( molar_densities[reaction.reactant_id(r)],
840  reaction.reactant_partial_order(r));
841  }
842  kfwd[rxn] *= fwd_conc[rxn];
843 
844  if(reaction.reversible())
845  {
846  StateType keq = reaction.equilibrium_constant( P0_RT, h_RT_minus_s_R );
847 
848  kbkwd_const[rxn] = kfwd_const[rxn]/keq;
849  kbkwd[rxn] = kbkwd_const[rxn];
850  for (unsigned int p=0; p<reaction.n_products(); p++)
851  {
852  bkwd_conc[rxn] *= pow( molar_densities[reaction.product_id(p)],
853  reaction.product_partial_order(p));
854  }
855  kbkwd[rxn] *= bkwd_conc[rxn];
856  }
857 
858  net_rates[rxn] = kfwd[rxn] - kbkwd[rxn];
859 
860  }
861 
862  return;
863  }
864 
865 
866  template<typename CoeffType>
867  inline
868  void ReactionSet<CoeffType>::print(std::ostream& os) const
869  {
870  os << "# Number of reactions: " << this->n_reactions() << "\n";
871 
872  for (unsigned int r=0; r < this->n_reactions(); r++)
873  {
874  os << "# " << r << '\n'
875  << (this->reaction(r)) << "\n";
876  }
877 
878  return;
879  }
880 
881 } // end namespace Antioch
882 
883 #endif // ANTIOCH_REACTION_SET_H
Antioch::ReactionSet::_reactions
std::vector< Reaction< CoeffType > * > _reactions
Definition: reaction_set.h:183
antioch_assert
#define antioch_assert(asserted)
Definition: antioch_asserts.h:58
Antioch::Constants::R_universal
CoeffType R_universal()
Universal Gas Constant, R, expressed in J/(mol-K)
Definition: physical_constants.h:42
Antioch::Reaction::equilibrium_constant
StateType equilibrium_constant(const StateType &P0_RT, const VectorStateType &h_RT_minus_s_R) const
Definition: reaction.h:751
antioch_warning
#define antioch_warning(message)
Definition: antioch_asserts.h:99
Antioch::ReactionSet::print_chemical_scheme
void print_chemical_scheme(std::ostream &output, const KineticsConditions< StateType, VectorStateType > &conditions, const VectorStateType &molar_densities, const VectorStateType &h_RT_minus_s_R, std::vector< VectorStateType > &lossMatrix, std::vector< VectorStateType > &prodMatrix, std::vector< VectorStateType > &netMatrix) const
Definition: reaction_set.h:556
Antioch::Reaction::compute_forward_rate_coefficient
StateType compute_forward_rate_coefficient(const VectorStateType &molar_densities, const KineticsConditions< StateType, VectorStateType > &conditions) const
Definition: reaction.h:866
antioch_assert_equal_to
#define antioch_assert_equal_to(expr1, expr2)
Definition: antioch_asserts.h:74
Antioch::ReactionSet::reaction
const Reaction< CoeffType > & reaction(const unsigned int r) const
Definition: reaction_set.h:232
antioch_assert_greater
#define antioch_assert_greater(expr1, expr2)
Definition: antioch_asserts.h:77
Antioch::ReactionSet::set_parameter_of_reaction
void set_parameter_of_reaction(const std::string &reaction_id, const std::vector< std::string > &keywords, ParamType value)
change a parameter of a reaction
Definition: reaction_set.h:512
antioch_assert_less
#define antioch_assert_less(expr1, expr2)
Definition: antioch_asserts.h:76
Antioch::Reaction::reversible
bool reversible() const
Definition: reaction.h:462
Antioch::Reaction::product_stoichiometric_coefficient
unsigned int product_stoichiometric_coefficient(const unsigned int p) const
Definition: reaction.h:561
antioch_error
#define antioch_error()
Definition: antioch_asserts.h:104
Antioch::ReactionSet::remove_reaction
void remove_reaction(unsigned int nr)
remove a reaction from the system.
Definition: reaction_set.h:219
std::pow
Antioch::enable_if_c< Antioch::is_valarray< T >::value, typename Antioch::state_type< T >::type >::type pow(const T &in, const T2 &n)
Definition: valarray_utils.h:71
Antioch::ReactionSet::n_species
unsigned int n_species() const
Definition: reaction_set.h:193
Antioch::Reaction::reactant_id
unsigned int reactant_id(const unsigned int r) const
Definition: reaction.h:534
Antioch::Reaction::product_partial_order
CoeffType product_partial_order(const unsigned int p) const
Definition: reaction.h:579
std::max
max(const _Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &a, const _Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &b) ANTIOCH_AUTOFUNC(_Matrix< _Scalar ANTIOCH_COMMA _Rows ANTIOCH_COMMA _Cols ANTIOCH_COMMA _Options ANTIOCH_COMMA _MaxRows ANTIOCH_COMMA _MaxCols >
Antioch::ReactionSet::n_reactions
unsigned int n_reactions() const
Definition: reaction_set.h:200
Antioch::ReactionSet::print
void print(std::ostream &os=std::cout) const
Formatted print, by default to std::cout.
Definition: reaction_set.h:868
Antioch::ReactionSet::compute_reaction_rates_and_derivs
void compute_reaction_rates_and_derivs(const KineticsConditions< StateType, VectorStateType > &conditions, const VectorStateType &molar_densities, const VectorStateType &h_RT_minus_s_R, const VectorStateType &dh_RT_minus_s_R_dT, VectorReactionsType &net_reaction_rates, VectorReactionsType &dnet_rate_dT, MatrixReactionsType &dnet_rate_dX_s) const
Compute the rates of progress and derivatives for each reaction.
Definition: reaction_set.h:303
Antioch::string_to_kin_enum
KineticsModel::Parameters string_to_kin_enum(const std::string &str)
Definition: string_utils.h:192
Antioch::ReactionSet::reaction_by_id
unsigned int reaction_by_id(const std::string &reaction_id) const
Definition: reaction_set.h:347
Antioch::ReactionSet::chemical_mixture
const ChemicalMixture< CoeffType > & chemical_mixture() const
Definition: reaction_set.h:248
Antioch::ReactionSet::get_parameter_of_reaction
CoeffType get_parameter_of_reaction(const std::string &reaction_id, const std::vector< std::string > &keywords) const
Definition: reaction_set.h:462
Antioch::ReactionSet::find_chemical_process_parameter
void find_chemical_process_parameter(ReactionType::Parameters paramChem, const std::vector< std::string > &keywords, unsigned int &species) const
helper function
Definition: reaction_set.h:438
Antioch::Reaction::n_products
unsigned int n_products() const
Definition: reaction.h:508
Antioch::set_zero
void set_zero(_Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &a)
Definition: eigen_utils.h:217
Antioch::ReactionSet::_P0_R
const CoeffType _P0_R
Scaling for equilibrium constant.
Definition: reaction_set.h:186
Antioch::ChemicalMixture
Class storing chemical mixture properties.
Definition: chemical_mixture.h:75
Antioch
The parameters are reduced parameters.
Definition: chemical_mixture.h:47
Antioch::Reaction::product_id
unsigned int product_id(const unsigned int p) const
Definition: reaction.h:543
Antioch::Reaction::reactant_stoichiometric_coefficient
unsigned int reactant_stoichiometric_coefficient(const unsigned int r) const
Definition: reaction.h:552
Antioch::Reaction::reactant_partial_order
CoeffType reactant_partial_order(const unsigned int r) const
Definition: reaction.h:570
Antioch::ReactionSet::_chem_mixture
const ChemicalMixture< CoeffType > & _chem_mixture
Definition: reaction_set.h:181
Antioch::ReactionSet::add_reaction
void add_reaction(Reaction< CoeffType > *reaction)
Add a reaction to the system.
Definition: reaction_set.h:207
Antioch::ReactionSet::compute_reaction_rates
void compute_reaction_rates(const KineticsConditions< StateType, VectorStateType > &conditions, const VectorStateType &molar_densities, const VectorStateType &h_RT_minus_s_R, VectorReactionsType &net_reaction_rates) const
Compute the rates of progress for each reaction.
Definition: reaction_set.h:275
Antioch::KineticsConditions::T
const StateType & T() const
Definition: kinetics_conditions.h:120
Antioch::Reaction::n_reactants
unsigned int n_reactants() const
Definition: reaction.h:498
Antioch::ReactionSet::find_kinetics_model_parameter
void find_kinetics_model_parameter(const unsigned int r, const std::vector< std::string > &keywords, unsigned int &nr, std::string &unit, int &l) const
helper function
Definition: reaction_set.h:374
Antioch::KineticsConditions
This class contains the conditions of the chemistry.
Definition: kinetics_conditions.h:64
Antioch::ReactionSet::get_reactive_scheme
void get_reactive_scheme(const KineticsConditions< StateType, VectorStateType > &conditions, const VectorStateType &molar_densities, const VectorStateType &h_RT_minus_s_R, VectorStateType &net_rates, VectorStateType &kfwd_const, VectorStateType &kbkwd_const, VectorStateType &kfwd, VectorStateType &kbkwd, VectorStateType &fwd_conc, VectorStateType &bkwd_conc) const
Definition: reaction_set.h:802
Antioch::string_to_chem_enum
ReactionType::Parameters string_to_chem_enum(const std::string &str)
Definition: string_utils.h:232
Generated on Thu Jul 7 2016 11:09:45 for antioch-0.4.0 by  doxygen 1.8.8