antioch-0.4.0
Functions
kinetics_regression_vec.C File Reference
#include "antioch_config.h"
#include <valarray>
#include "Eigen/Dense"
#include "metaphysicl/numberarray.h"
#include "vexcl/vexcl.hpp"
#include "antioch/eigen_utils_decl.h"
#include "antioch/metaphysicl_utils_decl.h"
#include "antioch/valarray_utils_decl.h"
#include "antioch/vector_utils_decl.h"
#include "antioch/vexcl_utils_decl.h"
#include "antioch/antioch_asserts.h"
#include "antioch/chemical_species.h"
#include "antioch/chemical_mixture.h"
#include "antioch/reaction_set.h"
#include "antioch/read_reaction_set_data.h"
#include "antioch/cea_thermo.h"
#include "antioch/kinetics_evaluator.h"
#include "antioch/eigen_utils.h"
#include "antioch/metaphysicl_utils.h"
#include "antioch/valarray_utils.h"
#include "antioch/vector_utils.h"
#include "antioch/vexcl_utils.h"
#include <limits>
#include <string>
#include <vector>

Go to the source code of this file.

Functions

template<typename SpeciesVector1 , typename SpeciesVector2 >
int vec_compare (const SpeciesVector1 &a, const SpeciesVector2 &b, const std::string &name)
 
template<typename PairScalars >
int vectester (const std::string &input_name, const PairScalars &example, const std::string &testname)
 
int main (int argc, char *argv[])
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 471 of file kinetics_regression_vec.C.

References antioch_error, and vectester().

472 {
473  // Check command line count.
474  if( argc < 2 )
475  {
476  // TODO: Need more consistent error handling.
477  std::cerr << "Error: Must specify reaction set XML input file." << std::endl;
478  antioch_error();
479  }
480 
481  int returnval = 0;
482 
483  returnval +=
484  vectester (argv[1], std::valarray<float>(2*ANTIOCH_N_TUPLES), "valarray<float>");
485  returnval +=
486  vectester (argv[1], std::valarray<double>(2*ANTIOCH_N_TUPLES), "valarray<double>");
487 // We're not getting the full long double precision yet?
488 // returnval = returnval ||
489 // vectester (argv[1], std::valarray<long double>(2*ANTIOCH_N_TUPLES), "valarray<ld>");
490 #ifdef ANTIOCH_HAVE_EIGEN
491  returnval +=
492  vectester (argv[1], Eigen::Array<float, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXf");
493  returnval +=
494  vectester (argv[1], Eigen::Array<double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXd");
495 // We're not getting the full long double precision yet?
496 // returnval = returnval ||
497 // vectester (argv[1], Eigen::Array<long double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXld");
498 #endif
499 #ifdef ANTIOCH_HAVE_METAPHYSICL
500  returnval +=
501  vectester (argv[1], MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray<float>");
502  returnval +=
503  vectester (argv[1], MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray<double>");
504 // returnval = returnval ||
505 // vectester (argv[1], MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray<ld>");
506 #endif
507 #ifdef ANTIOCH_HAVE_VEXCL
508  vex::Context ctx_f (vex::Filter::All);
509  if (!ctx_f.empty())
510  returnval = returnval ||
511  vectester (argv[1], vex::vector<float> (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector<float>");
512 
513  vex::Context ctx_d (vex::Filter::DoublePrecision);
514  if (!ctx_d.empty())
515  returnval = returnval ||
516  vectester (argv[1], vex::vector<double> (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector<double>");
517 #endif
518 
519  std::cout << "Found " << returnval << " errors" << std::endl;
520 
521 #ifdef ANTIOCH_HAVE_GRVY
522  gt.Finalize();
523  gt.Summarize();
524 #endif
525 
526  return (returnval != 0) ? 1 : 0;
527 }
#define antioch_error()
int vectester(const std::string &input_name, const PairScalars &example, const std::string &testname)
template<typename SpeciesVector1 , typename SpeciesVector2 >
int vec_compare ( const SpeciesVector1 &  a,
const SpeciesVector2 &  b,
const std::string &  name 
)

Definition at line 85 of file kinetics_regression_vec.C.

References std::max(), and Antioch::max().

Referenced by vectester().

86 {
87  int found_errors = 0;
88 
89  typedef typename Antioch::value_type<SpeciesVector1>::type StateType;
90  typedef typename Antioch::value_type<StateType>::type Scalar;
91 
92  if (static_cast<std::size_t>(a.size()) !=
93  static_cast<std::size_t>(b.size()))
94  {
95  std::cerr << "Error: Mismatch in vector sizes " << name << std::endl;
96  }
97 
98  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 200;
99 
100  for (unsigned int s=0; s != a.size(); s++)
101  {
102  using std::abs;
103  using std::max;
104 
105  // Break this expression up to workaround bugs in my Eigen
106  // and VexCL versions - RHS
107  const StateType as = a[s], bs = b[s];
108  const StateType rel_error = (as - bs)/max(as,bs);
109  const StateType abs_rel_error = abs(rel_error);
110 
111  if( Antioch::max(abs_rel_error) > tol )
112  {
113  found_errors++;
114  }
115  }
116 
117  if (found_errors)
118  {
119  using std::max;
120 
121  std::cerr << "Error: Mismatch in vectors " << name << std::endl;
122  for( unsigned int s = 0; s < n_species; s++)
123  {
124  const StateType as = a[s], bs = b[s];
125  const StateType rel_error = (as - bs)/max(as,bs);
126  std::cout << std::scientific << std::setprecision(16)
127  << "a(" << s << ") = " << as << std::endl
128  << "b(" << s << ") = " << bs << std::endl
129  << "a-b(" << s << ") = " << StateType(as-bs)
130  << std::endl <<
131  "a-b/(max(a,b)) = " << rel_error << std::endl
132  << "tol =" << tol
133  << std::endl;
134  }
135  }
136 
137  return found_errors;
138 }
Antioch::enable_if_c< is_eigen< T >::value, typename value_type< T >::type >::type max(const T &in)
Definition: eigen_utils.h:88
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 >
template<typename PairScalars >
int vectester ( const std::string &  input_name,
const PairScalars &  example,
const std::string &  testname 
)

Definition at line 141 of file kinetics_regression_vec.C.

References Antioch::ChemicalMixture< CoeffType >::chemical_species(), Antioch::KineticsEvaluator< CoeffType, StateType >::compute_mass_sources(), Antioch::KineticsEvaluator< CoeffType, StateType >::compute_mass_sources_and_derivs(), Antioch::CEAThermodynamics< CoeffType >::dh_RT_minus_s_R_dT(), Antioch::CEAThermodynamics< CoeffType >::h_RT_minus_s_R(), Antioch::init_constant(), Antioch::ChemicalMixture< CoeffType >::molar_densities, Antioch::ChemicalMixture< CoeffType >::R(), and vec_compare().

Referenced by main().

144 {
145  typedef typename Antioch::value_type<PairScalars>::type Scalar;
146 
147  std::vector<std::string> species_str_list;
148  species_str_list.reserve(n_species);
149  species_str_list.push_back( "N2" );
150  species_str_list.push_back( "O2" );
151  species_str_list.push_back( "N" );
152  species_str_list.push_back( "O" );
153  species_str_list.push_back( "NO" );
154 
155  Antioch::ChemicalMixture<Scalar> chem_mixture( species_str_list );
156  Antioch::ReactionSet<Scalar> reaction_set( chem_mixture );
157  Antioch::CEAThermodynamics<Scalar> thermo( chem_mixture );
158 
159  Antioch::read_reaction_set_data_xml<Scalar>( input_name, true, reaction_set );
160 
161  PairScalars T = example;
162  PairScalars P = example;
163 
164  // Mass fractions
165  PairScalars massfrac = example;
166 
167  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
168  {
169  T[2*tuple ] = 1500.0;
170  T[2*tuple+1] = 1500.0;
171 
172  P[2*tuple ] = 1.0e5;
173  P[2*tuple+1] = 1.0e5;
174 
175  massfrac[2*tuple ] = 0.2;
176  massfrac[2*tuple+1] = 0.2;
177  }
178 
179  const Antioch::KineticsConditions<PairScalars> conditions(T);
180  const std::vector<PairScalars> Y(n_species,massfrac);
181  std::vector<PairScalars> molar_densities(n_species, example);
182  std::vector<PairScalars> h_RT_minus_s_R(n_species, example);
183  std::vector<PairScalars> dh_RT_minus_s_R_dT(n_species, example);
184 
185  Antioch::KineticsEvaluator<Scalar,PairScalars> kinetics( reaction_set, example );
186 
187  std::vector<PairScalars> omega_dot(n_species, example);
188  std::vector<PairScalars> omega_dot_2(n_species, example);
189  std::vector<PairScalars> domega_dot_dT(n_species, example);
190 
191  std::vector<std::vector<PairScalars> > domega_dot_drhos
192  (n_species, omega_dot); // omega_dot is a good example vec<Pair>
193 
194 
195 #ifdef ANTIOCH_HAVE_GRVY
196  const std::string testnormal = testname + "-normal";
197  gt.BeginTimer(testnormal);
198 #endif
199 
200  const PairScalars R_mix = chem_mixture.R(Y);
201 
202  const PairScalars rho = P/(R_mix*T);
203 
204  chem_mixture.molar_densities(rho,Y,molar_densities);
205 
206  typedef typename Antioch::CEAThermodynamics<Scalar>::
207  template Cache<PairScalars> Cache;
208  thermo.h_RT_minus_s_R(Cache(T),h_RT_minus_s_R);
209  thermo.dh_RT_minus_s_R_dT(Cache(T),dh_RT_minus_s_R_dT);
210 
211  kinetics.compute_mass_sources( conditions, molar_densities, h_RT_minus_s_R, omega_dot );
212 
213  kinetics.compute_mass_sources_and_derivs ( conditions,
214  molar_densities,
215  h_RT_minus_s_R,
216  dh_RT_minus_s_R_dT,
217  omega_dot_2,
218  domega_dot_dT,
219  domega_dot_drhos );
220 
221 #ifdef ANTIOCH_HAVE_GRVY
222  gt.EndTimer(testnormal);
223 #endif
224 
225  int return_flag = 0;
226 
227 #ifdef ANTIOCH_HAVE_EIGEN
228  {
229  typedef Eigen::Array<PairScalars,n_species,1> SpeciesVecEigenType;
230 
231 #ifdef ANTIOCH_HAVE_GRVY
232  const std::string testeigenA = testname + "-eigenA";
233  gt.BeginTimer(testeigenA);
234 #endif
235 
237 
238  SpeciesVecEigenType eigen_Y;
239  Antioch::init_constant(eigen_Y, massfrac);
240 
241  SpeciesVecEigenType eigen_molar_densities;
242  Antioch::init_constant(eigen_molar_densities, example);
243 
244  SpeciesVecEigenType eigen_h_RT_minus_s_R;
245  Antioch::init_constant(eigen_h_RT_minus_s_R, example);
246 
247  SpeciesVecEigenType eigen_dh_RT_minus_s_R_dT;
248  Antioch::init_constant(eigen_dh_RT_minus_s_R_dT, example);
249 
250  SpeciesVecEigenType eigen_omega_dot;
251  Antioch::init_constant(eigen_omega_dot, example);
252 
253  SpeciesVecEigenType eigen_domega_dot_dT;
254  Antioch::init_constant(eigen_domega_dot_dT, example);
255 
256  // FIXME: What to do for domega_dot_drhos type?
257 
258  const PairScalars eigen_R = chem_mixture.R(eigen_Y);
259  chem_mixture.molar_densities(rho,eigen_Y,eigen_molar_densities);
260 
261  thermo.h_RT_minus_s_R(Cache(T),eigen_h_RT_minus_s_R);
262 
263  thermo.dh_RT_minus_s_R_dT(Cache(T),eigen_dh_RT_minus_s_R_dT);
264 
265  kinetics.compute_mass_sources( eigen_conditions, eigen_molar_densities, eigen_h_RT_minus_s_R, eigen_omega_dot );
266 
267 #ifdef ANTIOCH_HAVE_GRVY
268  gt.EndTimer(testeigenA);
269 #endif
270 
271  return_flag +=
272  vec_compare(eigen_R,R_mix,"eigen_R");
273 
274  return_flag +=
275  vec_compare(eigen_molar_densities, molar_densities,
276  "eigen_molar_densities");
277 
278  return_flag +=
279  vec_compare(eigen_h_RT_minus_s_R, h_RT_minus_s_R,
280  "eigen_h_RT_minus_s_R");
281 
282  return_flag +=
283  vec_compare(eigen_dh_RT_minus_s_R_dT, dh_RT_minus_s_R_dT,
284  "eigen_dh_RT_minus_s_R_dT");
285 
286  return_flag +=
287  vec_compare(eigen_omega_dot,omega_dot,"eigen_omega_dot");
288  }
289 
290  {
291  typedef Eigen::Matrix<PairScalars,Eigen::Dynamic,1> SpeciesVecEigenType;
293 
294  SpeciesVecEigenType eigen_Y(n_species,1);
295  Antioch::init_constant(eigen_Y, massfrac);
296 
297  SpeciesVecEigenType eigen_molar_densities(n_species,1);
298  Antioch::init_constant(eigen_molar_densities, example);
299 
300  SpeciesVecEigenType eigen_h_RT_minus_s_R(n_species,1);
301  Antioch::init_constant(eigen_h_RT_minus_s_R, example);
302 
303  SpeciesVecEigenType eigen_dh_RT_minus_s_R_dT(n_species,1);
304  Antioch::init_constant(eigen_dh_RT_minus_s_R_dT, example);
305 
306  SpeciesVecEigenType eigen_omega_dot(n_species,1);
307  Antioch::init_constant(eigen_omega_dot, example);
308 
309  SpeciesVecEigenType eigen_domega_dot_dT(n_species,1);
310  Antioch::init_constant(eigen_domega_dot_dT, example);
311 
312  // FIXME: What to do for domega_dot_drhos type?
313 
314 #ifdef ANTIOCH_HAVE_GRVY
315  const std::string testeigenV = testname + "-eigenV";
316  gt.BeginTimer(testeigenV);
317 #endif
318 
319  const PairScalars eigen_R = chem_mixture.R(eigen_Y);
320 
321  chem_mixture.molar_densities(rho,eigen_Y,eigen_molar_densities);
322 
323  thermo.h_RT_minus_s_R(Cache(T),eigen_h_RT_minus_s_R);
324 
325  thermo.dh_RT_minus_s_R_dT(Cache(T),eigen_dh_RT_minus_s_R_dT);
326 
327  kinetics.compute_mass_sources( eigen_conditions, eigen_molar_densities, eigen_h_RT_minus_s_R, eigen_omega_dot );
328 
329 #ifdef ANTIOCH_HAVE_GRVY
330  gt.EndTimer(testeigenV);
331 #endif
332 
333  return_flag +=
334  vec_compare(eigen_R,R_mix,"eigen_R");
335 
336  return_flag +=
337  vec_compare(eigen_molar_densities, molar_densities,
338  "eigen_molar_densities");
339 
340  return_flag +=
341  vec_compare(eigen_h_RT_minus_s_R, h_RT_minus_s_R,
342  "eigen_h_RT_minus_s_R");
343 
344  return_flag +=
345  vec_compare(eigen_dh_RT_minus_s_R_dT, dh_RT_minus_s_R_dT,
346  "eigen_dh_RT_minus_s_R_dT");
347 
348  return_flag +=
349  vec_compare(eigen_omega_dot,omega_dot,"eigen_omega_dot");
350  }
351 
352 #endif // ANTIOCH_HAVE_EIGEN
353 
354  for( unsigned int s = 0; s < n_species; s++)
355  {
356  std::cout << std::scientific << std::setprecision(16)
357  << "omega_dot(" << chem_mixture.chemical_species()[s]->species() << ") = "
358  << omega_dot[s] << std::endl;
359  }
360 
361  for( unsigned int s = 0; s < n_species; s++)
362  {
363  std::cout << std::scientific << std::setprecision(16)
364  << "domega_dot_dT(" << chem_mixture.chemical_species()[s]->species() << ") = "
365  << domega_dot_dT[s] << std::endl;
366  }
367 
368  for( unsigned int s = 0; s < n_species; s++)
369  {
370  for( unsigned int t = 0; t < n_species; t++)
371  {
372  std::cout << std::scientific << std::setprecision(16)
373  << "domega_dot_drhos(" << chem_mixture.chemical_species()[s]->species()
374  << ", " << chem_mixture.chemical_species()[t]->species() << ") = "
375  << domega_dot_drhos[s][t] << std::endl;
376  }
377  }
378 
379  // Regression values for omega_dot
380  std::vector<PairScalars> omega_dot_reg(n_species,example);
381 
382  // Regression values for domega_dot_drhos
383  std::vector<std::vector<PairScalars> > domega_dot_drhos_reg
384  (n_species, omega_dot); // omega_dot is the right size for an example
385 
386  // Regression values for domega_dot_dT
387  std::vector<PairScalars> domega_dot_dT_reg(n_species, example);
388 
389  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
390  {
391  omega_dot_reg[0][2*tuple ] = 9.1623705357123753e+04;
392  omega_dot_reg[1][2*tuple ] = -3.3462025680272243e+05;
393  omega_dot_reg[2][2*tuple ] = -2.1139216712069495e+05;
394  omega_dot_reg[3][2*tuple ] = 1.9782018625609628e+05;
395  omega_dot_reg[4][2*tuple ] = 2.5656853231019735e+05;
396  omega_dot_reg[0][2*tuple+1] = Scalar(omega_dot_reg[0][0]);
397  omega_dot_reg[1][2*tuple+1] = Scalar(omega_dot_reg[1][0]);
398  omega_dot_reg[2][2*tuple+1] = Scalar(omega_dot_reg[2][0]);
399  omega_dot_reg[3][2*tuple+1] = Scalar(omega_dot_reg[3][0]);
400  omega_dot_reg[4][2*tuple+1] = Scalar(omega_dot_reg[4][0]);
401 
402  domega_dot_dT_reg[0][2*tuple ] = 1.8014990183270937e+02;
403  domega_dot_dT_reg[1][2*tuple ] = -5.2724437115534380e+02;
404  domega_dot_dT_reg[2][2*tuple ] = -3.0930094476883017e+02;
405  domega_dot_dT_reg[3][2*tuple ] = 3.7972747459781005e+02;
406  domega_dot_dT_reg[4][2*tuple ] = 2.7666793949365456e+02;
407  domega_dot_dT_reg[0][2*tuple+1] = Scalar(domega_dot_dT_reg[0][0]);
408  domega_dot_dT_reg[1][2*tuple+1] = Scalar(domega_dot_dT_reg[1][0]);
409  domega_dot_dT_reg[2][2*tuple+1] = Scalar(domega_dot_dT_reg[2][0]);
410  domega_dot_dT_reg[3][2*tuple+1] = Scalar(domega_dot_dT_reg[3][0]);
411  domega_dot_dT_reg[4][2*tuple+1] = Scalar(domega_dot_dT_reg[4][0]);
412 
413  domega_dot_drhos_reg[0][0][2*tuple ] = 1.9675775188085109e+04;
414  domega_dot_drhos_reg[0][1][2*tuple ] = 1.7226141262419737e+04;
415  domega_dot_drhos_reg[0][2][2*tuple ] = 3.2159299284723610e+06;
416  domega_dot_drhos_reg[0][3][2*tuple ] = 1.4765214711933021e+05;
417  domega_dot_drhos_reg[0][4][2*tuple ] = 2.3225053279918131e+06;
418 
419  domega_dot_drhos_reg[1][0][2*tuple ] = 8.8927385505978492e+03;
420  domega_dot_drhos_reg[1][1][2*tuple ] = -9.9560178070099482e+06;
421  domega_dot_drhos_reg[1][2][2*tuple ] = -9.8748760140991123e+06;
422  domega_dot_drhos_reg[1][3][2*tuple ] = 4.6143036700500813e+05;
423  domega_dot_drhos_reg[1][4][2*tuple ] = 8.3487375168772399e+03;
424 
425  domega_dot_drhos_reg[2][0][2*tuple ] = -2.2420842426881281e+04;
426  domega_dot_drhos_reg[2][1][2*tuple ] = -4.3812843857644886e+06;
427  domega_dot_drhos_reg[2][2][2*tuple ] = -6.8343593463263955e+06;
428  domega_dot_drhos_reg[2][3][2*tuple ] = -5.4143671040862988e+05;
429  domega_dot_drhos_reg[2][4][2*tuple ] = -1.2267997668149246e+06;
430 
431  domega_dot_drhos_reg[3][0][2*tuple ] = -1.2028166578920147e+04;
432  domega_dot_drhos_reg[3][1][2*tuple ] = 4.9713710400172938e+06;
433  domega_dot_drhos_reg[3][2][2*tuple ] = 5.7418898143800552e+06;
434  domega_dot_drhos_reg[3][3][2*tuple ] = -9.1121284934572734e+05;
435  domega_dot_drhos_reg[3][4][2*tuple ] = 1.2431710353864791e+06;
436 
437  domega_dot_drhos_reg[4][0][2*tuple ] = 5.8804952671184686e+03;
438  domega_dot_drhos_reg[4][1][2*tuple ] = 9.3487050114947233e+06;
439  domega_dot_drhos_reg[4][2][2*tuple ] = 7.7514156175730915e+06;
440  domega_dot_drhos_reg[4][3][2*tuple ] = 8.4356704563001888e+05;
441  domega_dot_drhos_reg[4][4][2*tuple ] = -2.3472253340802449e+06;
442 
443  for (unsigned int i = 0; i != 5; ++i)
444  for (unsigned int j = 0; j != 5; ++j)
445  domega_dot_drhos_reg[i][j][2*tuple+1] =
446  Scalar(domega_dot_drhos_reg[i][j][2*tuple ]);
447  }
448 
449  return_flag +=
450  vec_compare(omega_dot, omega_dot_reg, "omega_dot");
451 
452  return_flag +=
453  vec_compare(omega_dot_2, omega_dot_reg, "omega_dot_2");
454 
455  return_flag +=
456  vec_compare(domega_dot_dT, domega_dot_dT_reg, "domega_dot_dT");
457 
458  char vecname[] = "domega_dot_drho_0";
459  for (unsigned int i = 0; i != 5; ++i)
460  {
461  return_flag +=
462  vec_compare(domega_dot_drhos[i], domega_dot_drhos_reg[i], vecname);
463 
464  // Should b=e safe to assume "1"+1 = "2" etc.
465  ++vecname[16];
466  }
467 
468  return return_flag;
469 }
This class encapsulates all the reaction mechanisms considered in a chemical nonequilibrium simulatio...
Class to handle computing mass source terms for a given ReactionSet.
Class storing chemical mixture properties.
void init_constant(Vector &output, const Scalar &example)
StateType h_RT_minus_s_R(const Cache< StateType > &cache, unsigned int species) const
We currently need different specializations for scalar vs vector inputs here.
Definition: cea_thermo.h:420
This class contains the conditions of the chemistry.
int vec_compare(const SpeciesVector1 &a, const SpeciesVector2 &b, const std::string &name)

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