antioch-0.4.0
reaction.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_H
29 #define ANTIOCH_REACTION_H
30 
31 // Antioch
32 #include "antioch/antioch_asserts.h"
33 #include "antioch/cmath_shims.h"
34 #include "antioch/kinetics_type.h"
35 #include "antioch/constant_rate.h"
36 #include "antioch/hercourtessen_rate.h"
37 #include "antioch/berthelot_rate.h"
38 #include "antioch/arrhenius_rate.h"
39 #include "antioch/berthelothercourtessen_rate.h"
40 #include "antioch/kooij_rate.h"
41 #include "antioch/vanthoff_rate.h"
42 #include "antioch/photochemical_rate.h"
43 #include "antioch/reaction_enum.h"
44 #include "antioch/chemical_mixture.h"
45 #include "antioch/kinetics_conditions.h"
46 #include "antioch/kinetics_parsing.h" // reset_parameter_of_rate
47 
48 //C++
49 #include <string>
50 #include <vector>
51 #include <iostream>
52 #include <limits>
53 
54 namespace Antioch
55 {
56  // Forward declarations
57  template <typename CoeffType>
58  class ElementaryReaction;
59 
60  template <typename CoeffType>
61  class DuplicateReaction;
62 
63  template <typename CoeffType>
64  class ThreeBodyReaction;
65 
66  template <typename CoeffType,typename FalloffType>
67  class FalloffReaction;
68 
69  template <typename CoeffType,typename FalloffType>
70  class FalloffThreeBodyReaction;
71 
72  template <typename CoeffType>
73  class LindemannFalloff;
74 
75  template <typename CoeffType>
76  class TroeFalloff;
77 
79 
107  template<typename CoeffType=double, typename VectorCoeffType = std::vector<CoeffType> >
108  class Reaction
109  {
110  public:
111 
113  Reaction( const unsigned int n_species, const std::string &equation,
114  const bool &reversible = true,
115  const ReactionType::ReactionType type = ReactionType::ELEMENTARY,
116  const KineticsModel::KineticsModel kin = KineticsModel::KOOIJ);
117 
118  virtual ~Reaction();
119 
120  unsigned int n_species() const;
121 
123  const std::string & equation() const;
124 
126  const std::string & id() const;
127 
129  void set_id(const std::string & id);
130 
134  ReactionType::ReactionType type() const;
135 
139  void set_type( const ReactionType::ReactionType type);
140 
143  void set_reversibility( const bool reversible);
144 
146  KineticsModel::KineticsModel kinetics_model() const;
147 
149  void set_kinetics_model( const KineticsModel::KineticsModel kin);
150 
151  bool initialized() const;
152 
154  void set_parameter_of_rate(KineticsModel::Parameters parameter, CoeffType new_value, unsigned int n_reaction = 0, const std::string & unit = "SI");
155 
157  void set_parameter_of_rate(KineticsModel::Parameters parameter, CoeffType new_value, unsigned int n_reaction, int l, const std::string & unit = "SI");
158 
159 
161  CoeffType get_parameter_of_rate(KineticsModel::Parameters parameter, unsigned int n_reaction = 0, const std::string & unit = "SI") const;
162 
164  CoeffType get_parameter_of_rate(KineticsModel::Parameters parameter, unsigned int n_reaction, const std::string & unit, int l) const;
165 
167  void set_parameter_of_chemical_process(ReactionType::Parameters parameter, CoeffType new_value, unsigned int species = std::numeric_limits<unsigned int>::max() );
168 
170  CoeffType get_parameter_of_chemical_process(ReactionType::Parameters parameter, unsigned int species = std::numeric_limits<unsigned int>::max() ) const;
171 
174  bool reversible() const;
175 
177  unsigned int n_reactants() const;
178 
180  unsigned int n_products() const;
181 
183  const std::string& reactant_name(const unsigned int r) const;
184 
186  const std::string& product_name(const unsigned int p) const;
187 
189  unsigned int reactant_id(const unsigned int r) const;
190 
192  unsigned int product_id(const unsigned int p) const;
193 
195  unsigned int reactant_stoichiometric_coefficient(const unsigned int r) const;
196 
198  unsigned int product_stoichiometric_coefficient(const unsigned int p) const;
199 
201  CoeffType reactant_partial_order(const unsigned int r) const;
202 
204  CoeffType product_partial_order(const unsigned int p) const;
205 
207  void add_reactant( const std::string &name,
208  const unsigned int r_id,
209  const unsigned int stoichiometric_coeff,
210  const CoeffType partial_order = std::numeric_limits<CoeffType>::infinity());// what test could be reliable?
211 
213  void add_product( const std::string &name,
214  const unsigned int p_id,
215  const unsigned int stoichiometric_coeff,
216  const CoeffType partial_order = std::numeric_limits<CoeffType>::infinity()); // what test could be reliable?
217 
219  void clear_reactant();
220 
222  void clear_product();
223 
225  void set_efficiency( const std::string &,
226  const unsigned int s,
227  const CoeffType efficiency);
228 
230  CoeffType get_efficiency( const unsigned int s) const;
231 
233  CoeffType efficiency( const unsigned int s) const;
234 
236  void initialize(unsigned int index = 0);
237 
239  int gamma() const;
240 
242  template <typename StateType, typename VectorStateType>
243  StateType equilibrium_constant( const StateType& P0_RT,
244  const VectorStateType& h_RT_minus_s_R ) const;
245 
247  template <typename StateType, typename VectorStateType>
248  void equilibrium_constant_and_derivative( const StateType& T,
249  const StateType& P0_RT,
250  const VectorStateType& h_RT_minus_s_R,
251  const VectorStateType& ddT_h_RT_minus_s_R,
252  StateType& keq,
253  StateType& dkeq_dT) const;
254 
256  template <typename StateType, typename VectorStateType>
257  StateType compute_forward_rate_coefficient( const VectorStateType& molar_densities,
258  const KineticsConditions<StateType,VectorStateType>& conditions) const;
259 
260  // Deprecated API for backwards compatibility
261  template <typename StateType, typename VectorStateType>
262  StateType compute_forward_rate_coefficient( const VectorStateType& molar_densities,
263  const StateType& temp)
264  const;
265 
268  template <typename StateType, typename VectorStateType>
269  void compute_forward_rate_coefficient_and_derivatives( const VectorStateType& molar_densities,
270  const KineticsConditions<StateType,VectorStateType>& conditions,
271  StateType& kfwd,
272  StateType& dkfwd_dT,
273  VectorStateType& dkfwd_dX) const;
274 
275  // Deprecated API for backwards compatibility
276  template <typename StateType, typename VectorStateType>
277  void compute_forward_rate_coefficient_and_derivatives( const VectorStateType& molar_densities,
278  const
279  StateType& temp,
280  StateType& kfwd,
281  StateType& dkfwd_dT,
282  VectorStateType& dkfwd_dX) const;
283 
285  template <typename StateType, typename VectorStateType>
286  StateType compute_rate_of_progress( const VectorStateType& molar_densities,
287  const KineticsConditions<StateType,VectorStateType>& conditions,
288  const StateType& P0_RT,
289  const VectorStateType& h_RT_minus_s_R) const;
290 
291  // Deprecated API for backwards compatibility
292  template <typename StateType, typename VectorStateType>
293  StateType compute_rate_of_progress( const VectorStateType& molar_densities,
294  const StateType& temp,
295  const StateType& P0_RT,
296  const VectorStateType& h_RT_minus_s_R) const;
297 
298  template <typename StateType, typename VectorStateType>
299  void compute_rate_of_progress_and_derivatives( const VectorStateType &molar_densities,
300  const ChemicalMixture<CoeffType>& /*chem_mixture*/, // fully useless, why is it here?
301  const KineticsConditions<StateType,VectorStateType>& conditions,
302  const StateType &P0_RT,
303  const VectorStateType &h_RT_minus_s_R,
304  const VectorStateType &dh_RT_minus_s_R_dT,
305  StateType& net_reaction_rate,
306  StateType& dnet_rate_dT,
307  VectorStateType& dnet_rate_dX_s ) const;
308 
309  // Deprecated API for backwards compatibility
310  template <typename StateType, typename VectorStateType>
311  void compute_rate_of_progress_and_derivatives( const VectorStateType &molar_densities,
312  const ChemicalMixture<CoeffType>& chem_mixture,
313  const StateType& temp,
314  const StateType &P0_RT,
315  const VectorStateType &h_RT_minus_s_R,
316  const VectorStateType &dh_RT_minus_s_R_dT,
317  StateType& net_reaction_rate,
318  StateType& dnet_rate_dT,
319  VectorStateType& dnet_rate_dX_s ) const;
320 
322  const KineticsType<CoeffType,VectorCoeffType>& forward_rate(unsigned int ir = 0) const;
323 
325  KineticsType<CoeffType,VectorCoeffType>& forward_rate(unsigned int ir = 0);
326 
328  //
329  // just a wrapper around the FalloffType & F() method
330  // of the FalloffReaction object, test consistency of type
331  template <typename FalloffType>
332  FalloffType & falloff();
333 
335  void add_forward_rate(KineticsType<CoeffType,VectorCoeffType> *rate);
336 
338  void swap_forward_rates(unsigned int irate, unsigned int jrate);
339 
341  unsigned int n_rate_constants() const;
342 
344  void print(std::ostream& os = std::cout) const;
345 
347  friend std::ostream& operator << (std::ostream& os, const Reaction &rxn)
348  {
349  rxn.print(os);
350  return os;
351  }
352 
353  protected:
354 
355  unsigned int _n_species;
356  std::string _id;
357  std::string _equation;
358  std::vector<std::string> _reactant_names;
359  std::vector<std::string> _product_names;
360  std::vector<unsigned int> _reactant_ids;
361  std::vector<unsigned int> _product_ids;
362  std::vector<unsigned int> _reactant_stoichiometry;
363  std::vector<unsigned int> _product_stoichiometry;
364  std::vector<unsigned int> _species_reactant_stoichiometry;
365  std::vector<unsigned int> _species_product_stoichiometry;
366  std::vector<CoeffType> _species_reactant_partial_order;
367  std::vector<CoeffType> _species_product_partial_order;
368  std::vector<int> _species_delta_stoichiometry;
369  int _gamma;
370  bool _initialized;
371  bool _reversible;
372  ReactionType::ReactionType _type;
373  KineticsModel::KineticsModel _kintype;
374 
376  std::vector<KineticsType<CoeffType,VectorCoeffType>* > _forward_rate;
377 
379  std::vector<CoeffType> _efficiencies;
380 
381  private:
382  Reaction();
383 
384  };
385 
386  /* ------------------------- Inline Functions -------------------------*/
387  template<typename CoeffType, typename VectorCoeffType>
388  inline
389  unsigned int Reaction<CoeffType,VectorCoeffType>::n_species() const
390  {
391  return _n_species;
392  }
393 
394  template<typename CoeffType, typename VectorCoeffType>
395  inline
396  const std::string & Reaction<CoeffType,VectorCoeffType>::id() const
397  {
398  return _id;
399  }
400 
401  template<typename CoeffType, typename VectorCoeffType>
402  inline
403  void Reaction<CoeffType,VectorCoeffType>::set_id(const std::string & id)
404  {
405  _id = id;
406  }
407 
408  template<typename CoeffType, typename VectorCoeffType>
409  inline
410  const std::string & Reaction<CoeffType,VectorCoeffType>::equation() const
411  {
412  return _equation;
413  }
414 
415  template<typename CoeffType, typename VectorCoeffType>
416  inline
417  ReactionType::ReactionType Reaction<CoeffType,VectorCoeffType>::type() const
418  {
419  return _type;
420  }
421 
422  template<typename CoeffType, typename VectorCoeffType>
423  inline
424  void Reaction<CoeffType,VectorCoeffType>::set_type( const ReactionType::ReactionType type)
425  {
426  _type = type;
427  return;
428  }
429 
430  template<typename CoeffType, typename VectorCoeffType>
431  inline
432  void Reaction<CoeffType,VectorCoeffType>::set_reversibility( const bool reversible)
433  {
434  _reversible = reversible;
435  return;
436  }
437 
438  template<typename CoeffType,typename VectorCoeffType>
439  inline
440  KineticsModel::KineticsModel Reaction<CoeffType,VectorCoeffType>::kinetics_model() const
441  {
442  return _kintype;
443  }
444 
445  template<typename CoeffType, typename VectorCoeffType>
446  inline
447  void Reaction<CoeffType,VectorCoeffType>::set_kinetics_model( const KineticsModel::KineticsModel kin)
448  {
449  _kintype = kin;
450  return;
451  }
452 
453  template<typename CoeffType, typename VectorCoeffType>
454  inline
455  bool Reaction<CoeffType,VectorCoeffType>::initialized() const
456  {
457  return _initialized;
458  }
459 
460  template<typename CoeffType, typename VectorCoeffType>
461  inline
462  bool Reaction<CoeffType,VectorCoeffType>::reversible() const
463  {
464  return _reversible;
465  }
466 
467  template<typename CoeffType, typename VectorCoeffType>
468  inline
469  void Reaction<CoeffType,VectorCoeffType>::add_forward_rate(KineticsType<CoeffType,VectorCoeffType> *rate)
470  {
471  _forward_rate.push_back(rate);
472  return;
473  }
474 
475  template<typename CoeffType, typename VectorCoeffType>
476  inline
477  unsigned int Reaction<CoeffType,VectorCoeffType>::n_rate_constants() const
478  {
479  return _forward_rate.size();
480  }
481 
482  template<typename CoeffType, typename VectorCoeffType>
483  inline
484  void Reaction<CoeffType,VectorCoeffType>::swap_forward_rates(unsigned int irate, unsigned int jrate)
485  {
486  antioch_assert_less(irate,_forward_rate.size());
487  antioch_assert_less(jrate,_forward_rate.size());
488 
489  KineticsType<CoeffType>* rate_tmp = _forward_rate[jrate];
490  _forward_rate[jrate] = _forward_rate[irate];
491  _forward_rate[irate] = rate_tmp;
492 
493  return;
494  }
495 
496  template<typename CoeffType, typename VectorCoeffType>
497  inline
498  unsigned int Reaction<CoeffType,VectorCoeffType>::n_reactants () const
499  {
500  antioch_assert_less(_reactant_ids.size(), this->n_species());
501  antioch_assert_equal_to(_reactant_ids.size(), _reactant_stoichiometry.size());
502  antioch_assert_equal_to(_reactant_ids.size(), _reactant_names.size());
503  return _reactant_ids.size();
504  }
505 
506  template<typename CoeffType, typename VectorCoeffType>
507  inline
508  unsigned int Reaction<CoeffType,VectorCoeffType>::n_products() const
509  {
510  antioch_assert_less(_product_ids.size(), this->n_species());
511  antioch_assert_equal_to(_product_ids.size(), _product_stoichiometry.size());
512  antioch_assert_equal_to(_product_ids.size(), _product_names.size());
513  return _product_ids.size();
514  }
515 
516  template<typename CoeffType, typename VectorCoeffType>
517  inline
518  const std::string& Reaction<CoeffType,VectorCoeffType>::reactant_name(const unsigned int r) const
519  {
520  antioch_assert_less(r, _reactant_names.size());
521  return _reactant_names[r];
522  }
523 
524  template<typename CoeffType, typename VectorCoeffType>
525  inline
526  const std::string& Reaction<CoeffType,VectorCoeffType>::product_name(const unsigned int p) const
527  {
528  antioch_assert_less(p, _product_names.size());
529  return _product_names[p];
530  }
531 
532  template<typename CoeffType, typename VectorCoeffType>
533  inline
534  unsigned int Reaction<CoeffType,VectorCoeffType>::reactant_id(const unsigned int r) const
535  {
536  antioch_assert_less(r, _reactant_ids.size());
537  antioch_assert_less(_reactant_ids[r], this->n_species());
538  return _reactant_ids[r];
539  }
540 
541  template<typename CoeffType, typename VectorCoeffType>
542  inline
543  unsigned int Reaction<CoeffType,VectorCoeffType>::product_id(const unsigned int p) const
544  {
545  antioch_assert_less(p, _product_ids.size());
546  antioch_assert_less(_product_ids[p], this->n_species());
547  return _product_ids[p];
548  }
549 
550  template<typename CoeffType, typename VectorCoeffType>
551  inline
552  unsigned int Reaction<CoeffType,VectorCoeffType>::reactant_stoichiometric_coefficient(const unsigned int r) const
553  {
554  antioch_assert_less(r, _reactant_stoichiometry.size());
555  antioch_assert_less(_reactant_ids[r], this->n_species());
556  return _reactant_stoichiometry[r];
557  }
558 
559  template<typename CoeffType, typename VectorCoeffType>
560  inline
561  unsigned int Reaction<CoeffType,VectorCoeffType>::product_stoichiometric_coefficient(const unsigned int p) const
562  {
563  antioch_assert_less(p, _product_stoichiometry.size());
564  antioch_assert_less(_product_ids[p], this->n_species());
565  return _product_stoichiometry[p];
566  }
567 
568  template<typename CoeffType, typename VectorCoeffType>
569  inline
570  CoeffType Reaction<CoeffType,VectorCoeffType>::reactant_partial_order(const unsigned int r) const
571  {
572  antioch_assert_less(r, _species_reactant_partial_order.size());
573  antioch_assert_less(_reactant_ids[r], this->n_species());
574  return _species_reactant_partial_order[r];
575  }
576 
577  template<typename CoeffType, typename VectorCoeffType>
578  inline
579  CoeffType Reaction<CoeffType,VectorCoeffType>::product_partial_order(const unsigned int p) const
580  {
581  antioch_assert_less(p, _species_product_partial_order.size());
582  antioch_assert_less(_product_ids[p], this->n_species());
583  return _species_product_partial_order[p];
584  }
585 
586  template<typename CoeffType, typename VectorCoeffType>
587  inline
588  void Reaction<CoeffType,VectorCoeffType>::add_reactant (const std::string &name,
589  const unsigned int r_id,
590  const unsigned int stoichiometric_coeff,
591  const CoeffType partial_order)
592  {
593  antioch_assert_less(r_id, this->n_species());
594  _reactant_names.push_back(name);
595  _reactant_ids.push_back(r_id);
596  _reactant_stoichiometry.push_back(stoichiometric_coeff);
597 
598  CoeffType order = (partial_order == std::numeric_limits<CoeffType>::infinity() )?static_cast<CoeffType>(stoichiometric_coeff):partial_order;
599  _species_reactant_partial_order.push_back(order);
600  return;
601  }
602 
603  template<typename CoeffType, typename VectorCoeffType>
604  inline
605  void Reaction<CoeffType,VectorCoeffType>::add_product (const std::string &name,
606  const unsigned int p_id,
607  const unsigned int stoichiometric_coeff,
608  const CoeffType partial_order)
609  {
610  antioch_assert_less(p_id, this->n_species());
611  _product_names.push_back(name);
612  _product_ids.push_back(p_id);
613  _product_stoichiometry.push_back(stoichiometric_coeff);
614 
615  CoeffType order = (partial_order == std::numeric_limits<CoeffType>::infinity() )?static_cast<CoeffType>(stoichiometric_coeff):partial_order;
616  _species_product_partial_order.push_back(order);
617  return;
618  }
619 
620  template<typename CoeffType, typename VectorCoeffType>
621  inline
622  void Reaction<CoeffType,VectorCoeffType>::clear_reactant()
623  {
624  _reactant_names.clear();
625  _reactant_ids.clear();
626  _reactant_stoichiometry.clear();
627  }
628 
629  template<typename CoeffType, typename VectorCoeffType>
630  inline
631  void Reaction<CoeffType,VectorCoeffType>::clear_product()
632  {
633  _product_names.clear();
634  _product_ids.clear();
635  _product_stoichiometry.clear();
636  }
637 
638  template<typename CoeffType, typename VectorCoeffType>
639  inline
640  void Reaction<CoeffType,VectorCoeffType>::set_efficiency (const std::string & , //? where does that come from?
641  const unsigned int s,
642  const CoeffType efficiency)
643  {
644  antioch_assert_less(s, _efficiencies.size());
645  _efficiencies[s] = efficiency;
646  }
647 
648  template<typename CoeffType, typename VectorCoeffType>
649  inline
650  CoeffType Reaction<CoeffType,VectorCoeffType>::get_efficiency (const unsigned int s) const
651  {
652  antioch_assert_less(s, _efficiencies.size());
653  return _efficiencies[s];
654  }
655 
656  template<typename CoeffType, typename VectorCoeffType>
657  inline
658  CoeffType Reaction<CoeffType,VectorCoeffType>::efficiency( const unsigned int s ) const
659  {
660 
661  antioch_assert(s < _efficiencies.size());
662  return _efficiencies[s];
663  }
664 
665  template<typename CoeffType, typename VectorCoeffType>
666  inline
667  int Reaction<CoeffType,VectorCoeffType>::gamma () const
668  {
669  return _gamma;
670  }
671 
672  template<typename CoeffType, typename VectorCoeffType>
673  inline
674  const KineticsType<CoeffType,VectorCoeffType>& Reaction<CoeffType,VectorCoeffType>::forward_rate(unsigned int ir) const
675  {
676  antioch_assert_less(ir,_forward_rate.size());
677  return *_forward_rate[ir];
678  }
679 
680  template<typename CoeffType, typename VectorCoeffType>
681  inline
682  KineticsType<CoeffType,VectorCoeffType>& Reaction<CoeffType,VectorCoeffType>::forward_rate(unsigned int ir)
683  {
684  antioch_assert_less(ir,_forward_rate.size());
685  return *_forward_rate[ir];
686  }
687 
688 
689  template<typename CoeffType, typename VectorCoeffType>
690  inline
691  Reaction<CoeffType,VectorCoeffType>::Reaction( const unsigned int n_species,
692  const std::string &equation,
693  const bool &reversible,
694  const ReactionType::ReactionType type,
695  const KineticsModel::KineticsModel kin)
696  : _n_species(n_species),
697  _equation(equation),
698  _gamma(0),
699  _initialized(false),
700  _reversible(reversible),
701  _type(type),
702  _kintype(kin)
703  {
704  return;
705  }
706 
707 
708  template<typename CoeffType, typename VectorCoeffType>
709  inline
710  Reaction<CoeffType,VectorCoeffType>::~Reaction()
711  {
712  for(unsigned int ir = 0; ir < _forward_rate.size(); ir++)
713  {
714  delete _forward_rate[ir];
715  }
716 
717  return;
718  }
719 
720 
721  template<typename CoeffType, typename VectorCoeffType>
722  inline
723  void Reaction<CoeffType,VectorCoeffType>::initialize(unsigned int index)
724  {
725  _gamma = 0;
726  for (unsigned int r=0; r< this->n_reactants(); r++)
727  {
728  _gamma -= this->reactant_stoichiometric_coefficient(r);
729  }
730 
731  for (unsigned int p=0; p < this->n_products(); p++)
732  {
733  _gamma += this->product_stoichiometric_coefficient(p);
734  }
735 
736  // gives kinetics object index in reaction set
737  for(typename std::vector<KineticsType<CoeffType,VectorCoeffType>* >::iterator it = _forward_rate.begin();
738  it != _forward_rate.end(); it++)
739  {
740  (*it)->set_index(index);
741  }
742 
743  // set initialization flag
744  _initialized = true;
745  }
746 
747 
748  template<typename CoeffType, typename VectorCoeffType>
749  template<typename StateType, typename VectorStateType>
750  inline
751  StateType Reaction<CoeffType,VectorCoeffType>::equilibrium_constant( const StateType& P0_RT,
752  const VectorStateType& h_RT_minus_s_R ) const
753  {
754  antioch_assert( this->initialized() );
756  // antioch_assert_greater( P0_RT, 0.0 );
757  antioch_assert_greater( h_RT_minus_s_R.size(), 0 );
758  antioch_assert_equal_to( h_RT_minus_s_R.size(), this->n_species() );
759 
760 
761 // DrG0 = - reactants + product
762 // K = (P0/(RT))^gamma exp(-DrG0)
763 // exppower = -DrG0 = reactants - products
764  StateType exppower = ( static_cast<CoeffType>(_reactant_stoichiometry[0])*
765  h_RT_minus_s_R[_reactant_ids[0]] );
766 
767  for (unsigned int s=1; s < this->n_reactants(); s++)
768  {
769  exppower += ( static_cast<CoeffType>(_reactant_stoichiometry[s])*
770  h_RT_minus_s_R[_reactant_ids[s]] );
771  }
772 
773  for (unsigned int s=0; s < this->n_products(); s++)
774  {
775  exppower -= ( static_cast<CoeffType>(_product_stoichiometry[s])*
776  h_RT_minus_s_R[_product_ids[s]] );
777  }
778  return ant_pow( P0_RT, static_cast<CoeffType>(this->gamma()) ) * ant_exp(exppower);
779  }
780 
781 
782  template<typename CoeffType, typename VectorCoeffType>
783  template<typename StateType, typename VectorStateType>
784  inline
785  void Reaction<CoeffType,VectorCoeffType>::equilibrium_constant_and_derivative( const StateType& T,
786  const StateType& P0_RT,
787  const VectorStateType& h_RT_minus_s_R,
788  const VectorStateType& ddT_h_RT_minus_s_R,
789  StateType& keq,
790  StateType& dkeq_dT) const
791  {
792  antioch_assert(this->initialized());
793 
795  // antioch_assert_greater( P0_RT, 0.0 );
796  // antioch_assert_greater( T, 0.0 );
797  antioch_assert_greater( h_RT_minus_s_R.size(), 0 );
798  antioch_assert_equal_to( h_RT_minus_s_R.size(), this->n_species() );
799  antioch_assert_equal_to( ddT_h_RT_minus_s_R.size(), this->n_species() );
800 
801  // get the equilibrium constant
802  keq = this->equilibrium_constant( P0_RT, h_RT_minus_s_R );
803 
804  StateType ddT_exppower = ( static_cast<CoeffType>(_reactant_stoichiometry[0])*
805  ddT_h_RT_minus_s_R[_reactant_ids[0]] );
806 
807  for (unsigned int s=1; s<this->n_reactants(); s++)
808  ddT_exppower += ( static_cast<CoeffType>(_reactant_stoichiometry[s])*
809  ddT_h_RT_minus_s_R[_reactant_ids[s]] );
810 
811  for (unsigned int s=0; s<this->n_products(); s++)
812  ddT_exppower -= ( static_cast<CoeffType>(_product_stoichiometry[s])*
813  ddT_h_RT_minus_s_R[_product_ids[s]] );
814 
815  // compute its derivative
816  dkeq_dT = keq*(-static_cast<CoeffType>(this->gamma())/T + ddT_exppower);
817 
818  return;
819  }
820 
821 
822  template<typename CoeffType, typename VectorCoeffType>
823  inline
824  void Reaction<CoeffType,VectorCoeffType>::print( std::ostream& os ) const
825  {
826  os << "# Gas-Phase Reaction \"" << _equation << "\":\n";
827  if (this->n_species())
828  {
829  os << "# reactants: ";
830  for (unsigned int r=0; r<this->n_reactants(); r++)
831  os << this->reactant_name(r) << ":"
832  << this->reactant_stoichiometric_coefficient(r) << ","
833  << this->reactant_partial_order(r) << " ";
834  os << "\n"
835  << "# products: ";
836  for (unsigned int p=0; p<this->n_products(); p++)
837  os << this->product_name(p) << ":"
838  << this->product_stoichiometric_coefficient(p) << ","
839  << this->product_partial_order(p) << " ";
840  }
841  os << "\n# Chemical process: " << _type;
842  os << "\n# Kinetics model: " << _kintype;
843  (_reversible)?os << "\n# reversible":
844  os << "\n# irreversible";
845  for(unsigned int ir = 0; ir < _forward_rate.size(); ir++)
846  {
847  os << "\n# forward rate eqn: " << *_forward_rate[ir];
848  }
849 
850  if (_type == ReactionType::THREE_BODY ||
851  _type == ReactionType::LINDEMANN_FALLOFF_THREE_BODY ||
852  _type == ReactionType::TROE_FALLOFF_THREE_BODY)
853  {
854  os << "\n# efficiencies: ";
855  for (unsigned int s=0; s<this->n_species(); s++)
856  os << s << ":" << this->efficiency(s) << " ";
857  }
858  os << "\n#";
859  return;
860  }
861 
862 
863  template<typename CoeffType, typename VectorCoeffType>
864  template <typename StateType, typename VectorStateType>
865  inline
866  StateType Reaction<CoeffType,VectorCoeffType>::compute_forward_rate_coefficient( const VectorStateType& molar_densities,
867  const KineticsConditions<StateType,VectorStateType>& conditions) const
868  {
869  switch(_type)
870  {
871  case(ReactionType::ELEMENTARY):
872  {
873  return (static_cast<const ElementaryReaction<CoeffType>*>(this))->compute_forward_rate_coefficient(molar_densities,conditions);
874  }
875  break;
876 
877  case(ReactionType::DUPLICATE):
878  {
879  return (static_cast<const DuplicateReaction<CoeffType>*>(this))->compute_forward_rate_coefficient(molar_densities,conditions);
880  }
881  break;
882 
883  case(ReactionType::THREE_BODY):
884  {
885  return (static_cast<const ThreeBodyReaction<CoeffType>*>(this))->compute_forward_rate_coefficient(molar_densities,conditions);
886  }
887  break;
888 
889  case(ReactionType::LINDEMANN_FALLOFF):
890  {
891  return (static_cast<const FalloffReaction<CoeffType,LindemannFalloff<CoeffType> >*>(this))->compute_forward_rate_coefficient(molar_densities,conditions);
892  }
893  break;
894 
895  case(ReactionType::TROE_FALLOFF):
896  {
897  return (static_cast<const FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->compute_forward_rate_coefficient(molar_densities,conditions);
898  }
899  break;
900 
901  case(ReactionType::LINDEMANN_FALLOFF_THREE_BODY):
902  {
903  return (static_cast<const FalloffThreeBodyReaction<CoeffType,LindemannFalloff<CoeffType> >*>(this))->compute_forward_rate_coefficient(molar_densities,conditions);
904  }
905  break;
906 
907  case(ReactionType::TROE_FALLOFF_THREE_BODY):
908  {
909  return (static_cast<const FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->compute_forward_rate_coefficient(molar_densities,conditions);
910  }
911  break;
912 
913  default:
914  {
915  antioch_error();
916  }
917  } // switch(_type)
918 
919  // Dummy
920  return zero_clone(conditions.T());
921  }
922 
923  template<typename CoeffType, typename VectorCoeffType>
924  template <typename StateType, typename VectorStateType>
925  inline
926  StateType Reaction<CoeffType,VectorCoeffType>::compute_forward_rate_coefficient
927  ( const VectorStateType& molar_densities,
928  const StateType& temp) const
929  {
930  antioch_deprecated();
931  return compute_forward_rate_coefficient
932  (molar_densities,
933  KineticsConditions<StateType,VectorStateType>(temp));
934  }
935 
936  template<typename CoeffType, typename VectorCoeffType>
937  template <typename StateType, typename VectorStateType>
938  inline
939  void Reaction<CoeffType,VectorCoeffType>::compute_forward_rate_coefficient_and_derivatives( const VectorStateType& molar_densities,
940  const KineticsConditions<StateType,VectorStateType>& conditions,
941  StateType& kfwd,
942  StateType& dkfwd_dT,
943  VectorStateType& dkfwd_dX) const
944  {
945  switch(_type)
946  {
947  case(ReactionType::ELEMENTARY):
948  {
949  (static_cast<const ElementaryReaction<CoeffType>*>(this))->compute_forward_rate_coefficient_and_derivatives(molar_densities,conditions,kfwd,dkfwd_dT,dkfwd_dX);
950  }
951  break;
952 
953  case(ReactionType::DUPLICATE):
954  {
955  (static_cast<const DuplicateReaction<CoeffType>*>(this))->compute_forward_rate_coefficient_and_derivatives(molar_densities,conditions,kfwd,dkfwd_dT,dkfwd_dX);
956  }
957  break;
958 
959  case(ReactionType::THREE_BODY):
960  {
961  (static_cast<const ThreeBodyReaction<CoeffType>*>(this))->compute_forward_rate_coefficient_and_derivatives(molar_densities,conditions,kfwd,dkfwd_dT,dkfwd_dX);
962  }
963  break;
964 
965  case(ReactionType::LINDEMANN_FALLOFF):
966  {
967  (static_cast<const FalloffReaction<CoeffType,LindemannFalloff<CoeffType> >*>(this))->compute_forward_rate_coefficient_and_derivatives(molar_densities,conditions,kfwd,dkfwd_dT,dkfwd_dX);
968  }
969  break;
970 
971  case(ReactionType::TROE_FALLOFF):
972  {
973  (static_cast<const FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->compute_forward_rate_coefficient_and_derivatives(molar_densities,conditions,kfwd,dkfwd_dT,dkfwd_dX);
974  }
975 
976  case(ReactionType::LINDEMANN_FALLOFF_THREE_BODY):
977  {
978  (static_cast<const FalloffThreeBodyReaction<CoeffType,LindemannFalloff<CoeffType> >*>(this))->compute_forward_rate_coefficient_and_derivatives(molar_densities,conditions,kfwd,dkfwd_dT,dkfwd_dX);
979  }
980  break;
981 
982  case(ReactionType::TROE_FALLOFF_THREE_BODY):
983  {
984  (static_cast<const FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->compute_forward_rate_coefficient_and_derivatives(molar_densities,conditions,kfwd,dkfwd_dT,dkfwd_dX);
985  }
986  break;
987 
988  default:
989  {
990  antioch_error();
991  }
992 
993  } // switch(type)
994 
995  return;
996  }
997 
998  template<typename CoeffType, typename VectorCoeffType>
999  template <typename StateType, typename VectorStateType>
1000  inline
1001  void
1002  Reaction<CoeffType,VectorCoeffType>::compute_forward_rate_coefficient_and_derivatives
1003  ( const VectorStateType& molar_densities,
1004  const StateType& temp,
1005  StateType& kfwd,
1006  StateType& dkfwd_dT,
1007  VectorStateType& dkfwd_dX) const
1008  {
1009  antioch_deprecated();
1010  return compute_forward_rate_coefficient_and_derivatives
1011  (molar_densities,
1012  KineticsConditions<StateType,VectorStateType>(temp),
1013  kfwd, dkfwd_dT, dkfwd_dX);
1014  }
1015 
1016  //kfwd *prod_r [R]^nu_r - kbkwd * prod_p [P]^nu_p ( = - 1/nu_r d[R]/dt)
1017  template<typename CoeffType, typename VectorCoeffType>
1018  template <typename StateType, typename VectorStateType>
1019  inline
1020  StateType Reaction<CoeffType,VectorCoeffType>::compute_rate_of_progress( const VectorStateType& molar_densities,
1021  const KineticsConditions<StateType,VectorStateType>& conditions,
1022  const StateType& P0_RT,
1023  const VectorStateType& h_RT_minus_s_R) const
1024  {
1025  StateType kfwd = this->compute_forward_rate_coefficient(molar_densities,conditions);
1026  StateType kfwd_times_reactants = kfwd;
1027 
1028  // Rfwd
1029  for (unsigned int ro=0; ro < this->n_reactants(); ro++)
1030  {
1031  kfwd_times_reactants *=
1032  ant_pow( molar_densities[this->reactant_id(ro)],
1033  this->reactant_partial_order(ro));
1034  }
1035 
1036  StateType kbkwd_times_products = Antioch::zero_clone(kfwd_times_reactants);
1037  if(_reversible)
1038  {
1039 
1040  StateType Keq = this->equilibrium_constant( P0_RT, h_RT_minus_s_R );
1041  kbkwd_times_products = kfwd/Keq;
1042 
1043  // Rbkwd
1044  for (unsigned int po=0; po< this->n_products(); po++)
1045  {
1046  kbkwd_times_products *=
1047  ant_pow( molar_densities[this->product_id(po)],
1048  this->product_partial_order(po));
1049 
1050  }
1051  }
1052 
1053  return kfwd_times_reactants - kbkwd_times_products;
1054 
1055  }
1056 
1057  template<typename CoeffType, typename VectorCoeffType>
1058  template <typename StateType, typename VectorStateType>
1059  inline
1060  StateType Reaction<CoeffType,VectorCoeffType>::compute_rate_of_progress( const VectorStateType& molar_densities,
1061  const StateType& temp,
1062  const StateType& P0_RT,
1063  const VectorStateType& h_RT_minus_s_R) const
1064  {
1065  antioch_deprecated();
1066  return compute_rate_of_progress
1067  (molar_densities,
1068  KineticsConditions<StateType,VectorStateType>(temp),
1069  P0_RT, h_RT_minus_s_R);
1070  }
1071 
1072  template<typename CoeffType, typename VectorCoeffType>
1073  template <typename StateType, typename VectorStateType>
1074  inline
1075  void Reaction<CoeffType,VectorCoeffType>::compute_rate_of_progress_and_derivatives( const VectorStateType &molar_densities,
1076  const ChemicalMixture<CoeffType>& /*chem_mixture*/,
1077  const KineticsConditions<StateType,VectorStateType>& conditions,
1078  const StateType &P0_RT,
1079  const VectorStateType &h_RT_minus_s_R,
1080  const VectorStateType &dh_RT_minus_s_R_dT,
1081  StateType& net_reaction_rate,
1082  StateType& dnet_rate_dT,
1083  VectorStateType& dnet_rate_dX_s ) const
1084  {
1085  antioch_assert_equal_to (molar_densities.size(), this->n_species());
1086 
1087 // First the forward component, if reversible, compute and add the backward component
1088 
1089  StateType kfwd = Antioch::zero_clone(conditions.T());
1090  StateType dkfwd_dT = Antioch::zero_clone(conditions.T());
1091  VectorStateType dkfwd_dX_s = Antioch::zero_clone(molar_densities);
1092  VectorStateType dkbkwd_dX_s = Antioch::zero_clone(molar_densities);
1093 
1094  this->compute_forward_rate_coefficient_and_derivatives(molar_densities, conditions, kfwd, dkfwd_dT ,dkfwd_dX_s);
1095 
1096  // If users want to use valarrays, then the output reference sizes
1097  // had better already match the input value sizes...
1098  StateType dRfwd_dT = Antioch::zero_clone(kfwd);
1099 
1100  // We need to construct using an input StateType argument if we
1101  // want StateType==valarray to have the right sizes
1102  // valarray compatibility makes this a bit redundant, but not much
1103  // worse than the previous version
1105  VectorStateType dRfwd_dX_s(this->n_species(), kfwd);
1106 
1107  Antioch::set_zero(dRfwd_dX_s);
1108 
1109  // pre-fill the participating species partials with the rates
1110  for (unsigned int r=0; r< this->n_reactants(); r++)
1111  {
1112  dRfwd_dX_s[this->reactant_id(r)] = kfwd;
1113  }
1114 
1115  //init
1116  StateType facfwd = constant_clone(conditions.T(),1);
1117  dRfwd_dT = dkfwd_dT;
1118 
1119  // Rfwd & derivatives
1120  for (unsigned int ro=0; ro < this->n_reactants(); ro++)
1121  {
1122  const StateType val =
1123  ant_pow( molar_densities[this->reactant_id(ro)],
1124  this->reactant_partial_order(ro));
1125 
1126  const StateType dval =
1127  ( static_cast<CoeffType>(this->reactant_stoichiometric_coefficient(ro))*
1128  ant_pow( molar_densities[this->reactant_id(ro)],
1129  this->reactant_partial_order(ro) - 1)
1130  );
1131 
1132  facfwd *= val;
1133  dRfwd_dT *= val;
1134 
1135  for (unsigned int ri=0; ri<this->n_reactants(); ri++)
1136  {
1137  dRfwd_dX_s[this->reactant_id(ri)] *= (ri == ro) ? dval : val;
1138  }
1139  }
1140 
1141  for (unsigned int s = 0; s < this->n_species(); s++)
1142  {
1143  dRfwd_dX_s[s] += facfwd * dkfwd_dX_s[s];
1144  }
1145 
1146  net_reaction_rate = facfwd * kfwd;
1147 
1148  dnet_rate_dT = dRfwd_dT;
1149 
1150  for (unsigned int s = 0; s < this->n_species(); s++)
1151  {
1152  dnet_rate_dX_s[s] = dRfwd_dX_s[s];
1153  }
1154 
1155  if(_reversible)
1156  {
1157 
1158  //backward to be computed and added
1159 
1160  StateType keq = Antioch::zero_clone(conditions.T());
1161  StateType dkeq_dT = Antioch::zero_clone(conditions.T());
1162 
1163  equilibrium_constant_and_derivative( conditions.T(), P0_RT, h_RT_minus_s_R,
1164  dh_RT_minus_s_R_dT,
1165  keq, dkeq_dT );
1166 
1167  const StateType kbkwd = kfwd/keq;
1168  const StateType dkbkwd_dT = (dkfwd_dT - kbkwd*dkeq_dT)/keq;
1169  for(unsigned int s = 0; s < this->n_species(); s++)
1170  {
1171  dkbkwd_dX_s[s] = dkfwd_dX_s[s]/keq;
1172  }
1173 
1174  // If users want to use valarrays, then the output reference sizes
1175  // had better already match the input value sizes...
1176  StateType dRbkwd_dT = Antioch::zero_clone(dkbkwd_dT);
1177 
1178  // We need to construct using an input StateType argument if we
1179  // want StateType==valarray to have the right sizes
1180  // valarray compatibility makes this a bit redundant, but not much
1181  // worse than the previous version
1183  VectorStateType dRbkwd_dX_s(this->n_species(), kbkwd);
1184 
1185  Antioch::set_zero(dRbkwd_dX_s);
1186 
1187  // pre-fill the participating species partials with the rates
1188  for (unsigned int p=0; p < this->n_products(); p++)
1189  {
1190  dRbkwd_dX_s[this->product_id(p)] = kbkwd;
1191  }
1192 
1193  //init
1194  StateType facbkwd = constant_clone(conditions.T(),1);
1195  dRbkwd_dT = dkbkwd_dT;
1196 
1197  // Rbkwd & derivatives
1198  for (unsigned int po=0; po< this->n_products(); po++)
1199  {
1200  const StateType val =
1201  ant_pow( molar_densities[this->product_id(po)],
1202  this->product_partial_order(po));
1203 
1204  const StateType dval =
1205  ( static_cast<CoeffType>(this->product_stoichiometric_coefficient(po))*
1206  ant_pow( molar_densities[this->product_id(po)],
1207  this->product_partial_order(po) - 1)
1208  );
1209 
1210  facbkwd *= val;
1211  dRbkwd_dT *= val;
1212 
1213  for (unsigned int pi=0; pi<this->n_products(); pi++)
1214  {
1215  dRbkwd_dX_s[this->product_id(pi)] *= (pi == po) ? dval : val;
1216  }
1217 
1218  }
1219 
1220  for (unsigned int s = 0; s < this->n_species(); s++)
1221  {
1222  dRbkwd_dX_s[s] += facbkwd * dkbkwd_dX_s[s];
1223  }
1224 
1225  net_reaction_rate -= facbkwd * kbkwd;
1226 
1227  dnet_rate_dT -= dRbkwd_dT;
1228 
1229  for (unsigned int s = 0; s < this->n_species(); s++)
1230  {
1231  dnet_rate_dX_s[s] -= dRbkwd_dX_s[s];
1232  }
1233 
1234  } //end of if(_reversible) condition
1235 
1236  return;
1237  }
1238 
1239  template<typename CoeffType, typename VectorCoeffType>
1240  template <typename StateType, typename VectorStateType>
1241  inline
1242  void Reaction<CoeffType,VectorCoeffType>::compute_rate_of_progress_and_derivatives( const VectorStateType &molar_densities,
1243  const ChemicalMixture<CoeffType>& chem_mixture,
1244  const StateType& temp,
1245  const StateType &P0_RT,
1246  const VectorStateType &h_RT_minus_s_R,
1247  const VectorStateType &dh_RT_minus_s_R_dT,
1248  StateType& net_reaction_rate,
1249  StateType& dnet_rate_dT,
1250  VectorStateType& dnet_rate_dX_s ) const
1251  {
1252  antioch_deprecated();
1253  return compute_rate_of_progress_and_derivatives
1254  (molar_densities, chem_mixture,
1255  KineticsConditions<StateType,VectorStateType>(temp),
1256  P0_RT, h_RT_minus_s_R, dh_RT_minus_s_R_dT, net_reaction_rate,
1257  dnet_rate_dT, dnet_rate_dX_s);
1258  }
1259 
1260  template<typename CoeffType, typename VectorCoeffType>
1261  template <typename FalloffType>
1262  inline
1263  FalloffType & Reaction<CoeffType,VectorCoeffType>::falloff()
1264  {
1265  switch(_type)
1266  {
1267  case(ReactionType::LINDEMANN_FALLOFF):
1268  {
1269  return (static_cast<const FalloffReaction<CoeffType,LindemannFalloff<CoeffType> >*>(this))->F();
1270  }
1271  break;
1272 
1273  case(ReactionType::TROE_FALLOFF):
1274  {
1275  return (static_cast<const FalloffReaction<CoeffType,LindemannFalloff<CoeffType> >*>(this))->F();
1276  }
1277  break;
1278 
1279  case(ReactionType::LINDEMANN_FALLOFF_THREE_BODY):
1280  {
1281  return (static_cast<const FalloffThreeBodyReaction<CoeffType,LindemannFalloff<CoeffType> >*>(this))->F();
1282  }
1283  break;
1284 
1285  case(ReactionType::TROE_FALLOFF_THREE_BODY):
1286  {
1287  return (static_cast<const FalloffThreeBodyReaction<CoeffType,LindemannFalloff<CoeffType> >*>(this))->F();
1288  }
1289  break;
1290 
1291  default:
1292  {
1293  std::cerr << "You are trying to retrieve a Falloff object in a reaction that is not a falloff.\n"
1294  << "The reaction is " << *this << std::endl;
1295  antioch_error();
1296  return NULL;
1297  }
1298 
1299  }
1300  }
1301 
1302  template<typename CoeffType, typename VectorCoeffType>
1303  inline
1304  void Reaction<CoeffType,VectorCoeffType>::set_parameter_of_rate(KineticsModel::Parameters parameter, CoeffType new_value, unsigned int n_kin, const std::string & unit)
1305  {
1306  antioch_assert_less(n_kin,_forward_rate.size());
1307  reset_parameter_of_rate(*_forward_rate[n_kin], parameter, new_value, unit);
1308  }
1309 
1310  template<typename CoeffType, typename VectorCoeffType>
1311  inline
1312  void Reaction<CoeffType,VectorCoeffType>::set_parameter_of_rate(KineticsModel::Parameters parameter, CoeffType new_value, unsigned int n_kin, int l, const std::string & unit)
1313  {
1314  antioch_assert_less(n_kin,_forward_rate.size());
1315  reset_parameter_of_rate(*_forward_rate[n_kin], parameter, new_value, l, unit);
1316  }
1317 
1318  template<typename CoeffType, typename VectorCoeffType>
1319  inline
1320  CoeffType Reaction<CoeffType,VectorCoeffType>::get_parameter_of_rate(KineticsModel::Parameters parameter, unsigned int n_reaction, const std::string & /* unit*/) const
1321  {
1322  antioch_assert_less(n_reaction,_forward_rate.size());
1323  return _forward_rate[n_reaction]->get_parameter(parameter);
1324  }
1325 
1326  template<typename CoeffType, typename VectorCoeffType>
1327  inline
1328  CoeffType Reaction<CoeffType,VectorCoeffType>::get_parameter_of_rate(KineticsModel::Parameters parameter, unsigned int n_reaction, const std::string & /* unit*/, int l) const
1329  {
1330  antioch_assert_less(n_reaction,_forward_rate.size());
1331  return _forward_rate[n_reaction]->get_parameter(parameter,l);
1332  }
1333 
1334  template<typename CoeffType, typename VectorCoeffType>
1335  inline
1336  void Reaction<CoeffType,VectorCoeffType>::set_parameter_of_chemical_process(ReactionType::Parameters parameter, CoeffType new_value, unsigned int species)
1337  {
1338  switch(parameter)
1339  {
1340  case ReactionType::EFFICIENCIES:
1341  {
1342  this->set_efficiency(std::string(),species,new_value); // tests inside, if not three-body or species wrong
1343  }
1344  break;
1345  case ReactionType::TROE_ALPHA:
1346  {
1347  if(this->type() == ReactionType::TROE_FALLOFF)
1348  {
1349  (static_cast<FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().set_alpha(new_value);
1350  }
1351  else if(this->type() == ReactionType::TROE_FALLOFF_THREE_BODY)
1352  {
1353  (static_cast<FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().set_alpha(new_value);
1354  }else
1355  {
1356  antioch_error();
1357  }
1358  }
1359  break;
1360  case ReactionType::TROE_T1:
1361  {
1362  if(this->type() == ReactionType::TROE_FALLOFF)
1363  {
1364  (static_cast<FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().set_T1(new_value);
1365  }
1366  else if(this->type() == ReactionType::TROE_FALLOFF_THREE_BODY)
1367  {
1368  (static_cast<FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().set_T1(new_value);
1369  }else
1370  {
1371  antioch_error();
1372  }
1373  }
1374  break;
1375  case ReactionType::TROE_T2:
1376  {
1377  if(this->type() == ReactionType::TROE_FALLOFF)
1378  {
1379  (static_cast<FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().set_T2(new_value);
1380  }
1381  else if(this->type() == ReactionType::TROE_FALLOFF_THREE_BODY)
1382  {
1383  (static_cast<FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().set_T2(new_value);
1384  }else
1385  {
1386  antioch_error();
1387  }
1388  }
1389  break;
1390  case ReactionType::TROE_T3:
1391  {
1392  if(this->type() == ReactionType::TROE_FALLOFF)
1393  {
1394  (static_cast<FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().set_T3(new_value);
1395  }
1396  else if(this->type() == ReactionType::TROE_FALLOFF_THREE_BODY)
1397  {
1398  (static_cast<FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().set_T3(new_value);
1399  }else
1400  {
1401  antioch_error();
1402  }
1403  }
1404  break;
1405  default:
1406  {
1407  antioch_error();
1408  }
1409  break;
1410  }
1411  }
1412 
1413  template<typename CoeffType, typename VectorCoeffType>
1414  inline
1415  CoeffType Reaction<CoeffType,VectorCoeffType>::get_parameter_of_chemical_process(ReactionType::Parameters parameter, unsigned int species) const
1416  {
1417  CoeffType value(0);
1418  switch(parameter)
1419  {
1420  case ReactionType::EFFICIENCIES:
1421  {
1422  value = this->get_efficiency(species); // tests inside, if not three-body or species wrong
1423  }
1424  break;
1425  case ReactionType::TROE_ALPHA:
1426  {
1427  if(this->type() == ReactionType::TROE_FALLOFF)
1428  {
1429  value = (static_cast<const FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().get_alpha();
1430  }
1431  else if(this->type() == ReactionType::TROE_FALLOFF_THREE_BODY)
1432  {
1433  value = (static_cast<const FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().get_alpha();
1434  }else
1435  {
1436  antioch_error();
1437  }
1438  }
1439  break;
1440  case ReactionType::TROE_T1:
1441  {
1442  if(this->type() == ReactionType::TROE_FALLOFF)
1443  {
1444  value = (static_cast<const FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().get_T1();
1445  }
1446  else if(this->type() == ReactionType::TROE_FALLOFF_THREE_BODY)
1447  {
1448  value = (static_cast<const FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().get_T1();
1449  }else
1450  {
1451  antioch_error();
1452  }
1453  }
1454  break;
1455  case ReactionType::TROE_T2:
1456  {
1457  if(this->type() == ReactionType::TROE_FALLOFF)
1458  {
1459  value = (static_cast<const FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().get_T2();
1460  }
1461  else if(this->type() == ReactionType::TROE_FALLOFF_THREE_BODY)
1462  {
1463  value = (static_cast<const FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().get_T2();
1464  }else
1465  {
1466  antioch_error();
1467  }
1468  }
1469  break;
1470  case ReactionType::TROE_T3:
1471  {
1472  if(this->type() == ReactionType::TROE_FALLOFF)
1473  {
1474  value = (static_cast<const FalloffReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().get_T3();
1475  }
1476  else if(this->type() == ReactionType::TROE_FALLOFF_THREE_BODY)
1477  {
1478  value = (static_cast<const FalloffThreeBodyReaction<CoeffType,TroeFalloff<CoeffType> >*>(this))->F().get_T3();
1479  }else
1480  {
1481  antioch_error();
1482  }
1483  }
1484  break;
1485  default:
1486  {
1487  antioch_error();
1488  }
1489  break;
1490  }
1491 
1492  return value;
1493  }
1494 
1495 } // namespace Antioch
1496 
1497 #endif // ANTIOCH_REACTION_H
antioch_assert
#define antioch_assert(asserted)
Definition: antioch_asserts.h:58
Antioch::Reaction::gamma
int gamma() const
Definition: reaction.h:667
Antioch::Reaction::equilibrium_constant
StateType equilibrium_constant(const StateType &P0_RT, const VectorStateType &h_RT_minus_s_R) const
Definition: reaction.h:751
Antioch::Reaction::add_reactant
void add_reactant(const std::string &name, const unsigned int r_id, const unsigned int stoichiometric_coeff, const CoeffType partial_order=std::numeric_limits< CoeffType >::infinity())
Definition: reaction.h:588
Antioch::Reaction::initialized
bool initialized() const
Definition: reaction.h:455
Antioch::Reaction::falloff
FalloffType & falloff()
Return writeable reference to the falloff object, test type.
Definition: reaction.h:1263
Antioch::Reaction::_id
std::string _id
Definition: reaction.h:356
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::Reaction
A single reaction mechanism.
Definition: reaction.h:108
Antioch::LindemannFalloff
Simplest falloff model ( )
Definition: lindemann_falloff.h:52
Antioch::Reaction::get_parameter_of_chemical_process
CoeffType get_parameter_of_chemical_process(ReactionType::Parameters parameter, unsigned int species=std::numeric_limits< unsigned int >::max()) const
get a parameter from the chemical process
Definition: reaction.h:1415
Antioch::Reaction::initialize
void initialize(unsigned int index=0)
Computes derived quantities.
Definition: reaction.h:723
Antioch::Reaction::_n_species
unsigned int _n_species
Definition: reaction.h:355
Antioch::Reaction::_product_ids
std::vector< unsigned int > _product_ids
Definition: reaction.h:361
Antioch::Reaction::equation
const std::string & equation() const
Definition: reaction.h:410
Antioch::KineticsType
base class for kinetics models
Definition: kinetics_type.h:82
antioch_assert_equal_to
#define antioch_assert_equal_to(expr1, expr2)
Definition: antioch_asserts.h:74
Antioch::Reaction::_type
ReactionType::ReactionType _type
Definition: reaction.h:372
Antioch::Reaction::n_species
unsigned int n_species() const
Definition: reaction.h:389
antioch_assert_greater
#define antioch_assert_greater(expr1, expr2)
Definition: antioch_asserts.h:77
Antioch::FalloffThreeBodyReaction
Base class for falloff processes coupled with efficiencies.
Definition: falloff_threebody_reaction.h:102
Antioch::Reaction::_kintype
KineticsModel::KineticsModel _kintype
Definition: reaction.h:373
Antioch::reset_parameter_of_rate
void reset_parameter_of_rate(KineticsType< CoeffType, VectorCoeffType > &rate, KineticsModel::Parameters parameter, const CoeffType new_value, const std::string &unit="SI")
The rate constant models derived from the Arrhenius model have an activation energy which is stored a...
Definition: kinetics_parsing.h:235
Antioch::Reaction::set_type
void set_type(const ReactionType::ReactionType type)
Set the type of reaction.
Definition: reaction.h:424
Antioch::Reaction::_product_stoichiometry
std::vector< unsigned int > _product_stoichiometry
Definition: reaction.h:363
Antioch::Reaction::set_id
void set_id(const std::string &id)
set the reaction id.
Definition: reaction.h:403
antioch_assert_less
#define antioch_assert_less(expr1, expr2)
Definition: antioch_asserts.h:76
Antioch::Reaction::id
const std::string & id() const
Definition: reaction.h:396
Antioch::Reaction::_species_reactant_partial_order
std::vector< CoeffType > _species_reactant_partial_order
Definition: reaction.h:366
Antioch::Reaction::reversible
bool reversible() const
Definition: reaction.h:462
Antioch::Reaction::get_efficiency
CoeffType get_efficiency(const unsigned int s) const
Definition: reaction.h:650
Antioch::Reaction::set_parameter_of_rate
void set_parameter_of_rate(KineticsModel::Parameters parameter, CoeffType new_value, unsigned int n_reaction=0, const std::string &unit="SI")
reset a parameter from the rate constant
Definition: reaction.h:1304
Antioch::Reaction::_species_delta_stoichiometry
std::vector< int > _species_delta_stoichiometry
Definition: reaction.h:368
Antioch::Reaction::reactant_name
const std::string & reactant_name(const unsigned int r) const
Definition: reaction.h:518
Antioch::ThreeBodyReaction
A single reaction mechanism.
Definition: reaction.h:64
Antioch::Reaction::_reactant_stoichiometry
std::vector< unsigned int > _reactant_stoichiometry
Definition: reaction.h:362
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::Reaction::equilibrium_constant_and_derivative
void equilibrium_constant_and_derivative(const StateType &T, const StateType &P0_RT, const VectorStateType &h_RT_minus_s_R, const VectorStateType &ddT_h_RT_minus_s_R, StateType &keq, StateType &dkeq_dT) const
Definition: reaction.h:785
Antioch::Reaction::clear_reactant
void clear_reactant()
Definition: reaction.h:622
Antioch::FalloffReaction
Base class for falloff processes.
Definition: falloff_reaction.h:98
Antioch::Reaction::set_reversibility
void set_reversibility(const bool reversible)
Set the reversibility of reaction.
Definition: reaction.h:432
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
Antioch::constant_clone
_Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > constant_clone(const _Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &ex, const Scalar &value)
Definition: eigen_utils.h:181
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::Reaction::n_rate_constants
unsigned int n_rate_constants() const
Return the number of rate constant objects.
Definition: reaction.h:477
Antioch::Reaction::set_parameter_of_chemical_process
void set_parameter_of_chemical_process(ReactionType::Parameters parameter, CoeffType new_value, unsigned int species=std::numeric_limits< unsigned int >::max())
reset a parameter from the chemical process
Definition: reaction.h:1336
Antioch::Reaction::compute_rate_of_progress
StateType compute_rate_of_progress(const VectorStateType &molar_densities, const KineticsConditions< StateType, VectorStateType > &conditions, const StateType &P0_RT, const VectorStateType &h_RT_minus_s_R) const
Definition: reaction.h:1020
Antioch::Reaction::_product_names
std::vector< std::string > _product_names
Definition: reaction.h:359
Antioch::Reaction::clear_product
void clear_product()
Definition: reaction.h:631
Antioch::Reaction::compute_forward_rate_coefficient_and_derivatives
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
Definition: reaction.h:939
Antioch::Reaction::add_forward_rate
void add_forward_rate(KineticsType< CoeffType, VectorCoeffType > *rate)
Add a forward rate object.
Definition: reaction.h:469
Antioch::Reaction::compute_rate_of_progress_and_derivatives
void compute_rate_of_progress_and_derivatives(const VectorStateType &molar_densities, const ChemicalMixture< CoeffType > &, const KineticsConditions< StateType, VectorStateType > &conditions, const StateType &P0_RT, const VectorStateType &h_RT_minus_s_R, const VectorStateType &dh_RT_minus_s_R_dT, StateType &net_reaction_rate, StateType &dnet_rate_dT, VectorStateType &dnet_rate_dX_s) const
Definition: reaction.h:1075
Antioch::Reaction::_forward_rate
std::vector< KineticsType< CoeffType, VectorCoeffType > * > _forward_rate
The forward reaction rate modified Arrhenius form.
Definition: reaction.h:376
Antioch::Reaction::_species_product_partial_order
std::vector< CoeffType > _species_product_partial_order
Definition: reaction.h:367
antioch_deprecated
#define antioch_deprecated()
Definition: antioch_asserts.h:117
Antioch::DuplicateReaction
A single reaction mechanism.
Definition: duplicate_reaction.h:64
Antioch::Reaction::add_product
void add_product(const std::string &name, const unsigned int p_id, const unsigned int stoichiometric_coeff, const CoeffType partial_order=std::numeric_limits< CoeffType >::infinity())
Definition: reaction.h:605
Antioch::Reaction::get_parameter_of_rate
CoeffType get_parameter_of_rate(KineticsModel::Parameters parameter, unsigned int n_reaction=0, const std::string &unit="SI") const
get a parameter from the rate constant
Definition: reaction.h:1320
Antioch::Reaction::forward_rate
const KineticsType< CoeffType, VectorCoeffType > & forward_rate(unsigned int ir=0) const
Return const reference to the forward rate object.
Definition: reaction.h:674
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::Reaction::~Reaction
virtual ~Reaction()
Definition: reaction.h:710
Antioch::Reaction::_reactant_names
std::vector< std::string > _reactant_names
Definition: reaction.h:358
Antioch::Reaction::_species_reactant_stoichiometry
std::vector< unsigned int > _species_reactant_stoichiometry
Definition: reaction.h:364
Antioch::Reaction::set_kinetics_model
void set_kinetics_model(const KineticsModel::KineticsModel kin)
Set the model of kinetics.
Definition: reaction.h:447
Antioch::ElementaryReaction
A single reaction mechanism.
Definition: elementary_reaction.h:64
Antioch::Reaction::_equation
std::string _equation
Definition: reaction.h:357
Antioch::ChemicalMixture
Class storing chemical mixture properties.
Definition: chemical_mixture.h:75
Antioch::Reaction::_reactant_ids
std::vector< unsigned int > _reactant_ids
Definition: reaction.h:360
Antioch::Reaction::type
ReactionType::ReactionType type() const
Type of reaction.
Definition: reaction.h:417
Antioch
The parameters are reduced parameters.
Definition: chemical_mixture.h:47
Antioch::Reaction::_species_product_stoichiometry
std::vector< unsigned int > _species_product_stoichiometry
Definition: reaction.h:365
Antioch::Reaction::product_id
unsigned int product_id(const unsigned int p) const
Definition: reaction.h:543
Antioch::Reaction::kinetics_model
KineticsModel::KineticsModel kinetics_model() const
Model of kinetics.
Definition: reaction.h:440
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::Reaction::set_efficiency
void set_efficiency(const std::string &, const unsigned int s, const CoeffType efficiency)
Definition: reaction.h:640
Antioch::zero_clone
_Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > zero_clone(const _Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &ex)
Definition: eigen_utils.h:145
Antioch::KineticsConditions::T
const StateType & T() const
Definition: kinetics_conditions.h:120
Antioch::Reaction::product_name
const std::string & product_name(const unsigned int p) const
Definition: reaction.h:526
Antioch::Reaction::print
void print(std::ostream &os=std::cout) const
Formatted print, by default to std::cout.
Definition: reaction.h:824
Antioch::Reaction::swap_forward_rates
void swap_forward_rates(unsigned int irate, unsigned int jrate)
Swap two forward rates object.
Definition: reaction.h:484
Antioch::Reaction::efficiency
CoeffType efficiency(const unsigned int s) const
Definition: reaction.h:658
Antioch::Reaction::n_reactants
unsigned int n_reactants() const
Definition: reaction.h:498
Antioch::Constants::pi
CoeffType pi()
pi
Definition: math_constants.h:51
Antioch::KineticsConditions
This class contains the conditions of the chemistry.
Definition: kinetics_conditions.h:64
Antioch::TroeFalloff
The Troe falloff model is defined by: with The derivatives are therefore: .
Definition: reaction.h:76
Antioch::Reaction::operator<<
friend std::ostream & operator<<(std::ostream &os, const Reaction &rxn)
Formatted print.
Definition: reaction.h:347
Antioch::Reaction::_efficiencies
std::vector< CoeffType > _efficiencies
efficiencies for three body reactions
Definition: reaction.h:379
Generated on Thu Jul 7 2016 11:09:45 for antioch-0.4.0 by  doxygen 1.8.8