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" );
159 Antioch::read_reaction_set_data_xml<Scalar>( input_name,
true, reaction_set );
161 PairScalars T = example;
162 PairScalars P = example;
165 PairScalars massfrac = example;
167 for (
unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
169 T[2*tuple ] = 1500.0;
170 T[2*tuple+1] = 1500.0;
173 P[2*tuple+1] = 1.0e5;
175 massfrac[2*tuple ] = 0.2;
176 massfrac[2*tuple+1] = 0.2;
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);
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);
191 std::vector<std::vector<PairScalars> > domega_dot_drhos
192 (n_species, omega_dot);
195 #ifdef ANTIOCH_HAVE_GRVY
196 const std::string testnormal = testname +
"-normal";
197 gt.BeginTimer(testnormal);
200 const PairScalars R_mix = chem_mixture.R(Y);
202 const PairScalars rho = P/(R_mix*T);
204 chem_mixture.molar_densities(rho,Y,molar_densities);
207 template Cache<PairScalars> Cache;
209 thermo.dh_RT_minus_s_R_dT(Cache(T),dh_RT_minus_s_R_dT);
211 kinetics.compute_mass_sources( conditions, molar_densities, h_RT_minus_s_R, omega_dot );
213 kinetics.compute_mass_sources_and_derivs ( conditions,
221 #ifdef ANTIOCH_HAVE_GRVY
222 gt.EndTimer(testnormal);
227 #ifdef ANTIOCH_HAVE_EIGEN
229 typedef Eigen::Array<PairScalars,n_species,1> SpeciesVecEigenType;
231 #ifdef ANTIOCH_HAVE_GRVY
232 const std::string testeigenA = testname +
"-eigenA";
233 gt.BeginTimer(testeigenA);
238 SpeciesVecEigenType eigen_Y;
241 SpeciesVecEigenType eigen_molar_densities;
244 SpeciesVecEigenType eigen_h_RT_minus_s_R;
247 SpeciesVecEigenType eigen_dh_RT_minus_s_R_dT;
250 SpeciesVecEigenType eigen_omega_dot;
253 SpeciesVecEigenType eigen_domega_dot_dT;
258 const PairScalars eigen_R = chem_mixture.R(eigen_Y);
259 chem_mixture.molar_densities(rho,eigen_Y,eigen_molar_densities);
261 thermo.h_RT_minus_s_R(Cache(T),eigen_h_RT_minus_s_R);
263 thermo.dh_RT_minus_s_R_dT(Cache(T),eigen_dh_RT_minus_s_R_dT);
265 kinetics.compute_mass_sources( eigen_conditions, eigen_molar_densities, eigen_h_RT_minus_s_R, eigen_omega_dot );
267 #ifdef ANTIOCH_HAVE_GRVY
268 gt.EndTimer(testeigenA);
275 vec_compare(eigen_molar_densities, molar_densities,
276 "eigen_molar_densities");
280 "eigen_h_RT_minus_s_R");
283 vec_compare(eigen_dh_RT_minus_s_R_dT, dh_RT_minus_s_R_dT,
284 "eigen_dh_RT_minus_s_R_dT");
287 vec_compare(eigen_omega_dot,omega_dot,
"eigen_omega_dot");
291 typedef Eigen::Matrix<PairScalars,Eigen::Dynamic,1> SpeciesVecEigenType;
294 SpeciesVecEigenType eigen_Y(n_species,1);
297 SpeciesVecEigenType eigen_molar_densities(n_species,1);
300 SpeciesVecEigenType eigen_h_RT_minus_s_R(n_species,1);
303 SpeciesVecEigenType eigen_dh_RT_minus_s_R_dT(n_species,1);
306 SpeciesVecEigenType eigen_omega_dot(n_species,1);
309 SpeciesVecEigenType eigen_domega_dot_dT(n_species,1);
314 #ifdef ANTIOCH_HAVE_GRVY
315 const std::string testeigenV = testname +
"-eigenV";
316 gt.BeginTimer(testeigenV);
319 const PairScalars eigen_R = chem_mixture.R(eigen_Y);
321 chem_mixture.molar_densities(rho,eigen_Y,eigen_molar_densities);
323 thermo.h_RT_minus_s_R(Cache(T),eigen_h_RT_minus_s_R);
325 thermo.dh_RT_minus_s_R_dT(Cache(T),eigen_dh_RT_minus_s_R_dT);
327 kinetics.compute_mass_sources( eigen_conditions, eigen_molar_densities, eigen_h_RT_minus_s_R, eigen_omega_dot );
329 #ifdef ANTIOCH_HAVE_GRVY
330 gt.EndTimer(testeigenV);
337 vec_compare(eigen_molar_densities, molar_densities,
338 "eigen_molar_densities");
342 "eigen_h_RT_minus_s_R");
345 vec_compare(eigen_dh_RT_minus_s_R_dT, dh_RT_minus_s_R_dT,
346 "eigen_dh_RT_minus_s_R_dT");
349 vec_compare(eigen_omega_dot,omega_dot,
"eigen_omega_dot");
352 #endif // ANTIOCH_HAVE_EIGEN
354 for(
unsigned int s = 0; s < n_species; s++)
356 std::cout << std::scientific << std::setprecision(16)
357 <<
"omega_dot(" << chem_mixture.chemical_species()[s]->species() <<
") = "
358 << omega_dot[s] << std::endl;
361 for(
unsigned int s = 0; s < n_species; s++)
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;
368 for(
unsigned int s = 0; s < n_species; s++)
370 for(
unsigned int t = 0; t < n_species; t++)
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;
380 std::vector<PairScalars> omega_dot_reg(n_species,example);
383 std::vector<std::vector<PairScalars> > domega_dot_drhos_reg
384 (n_species, omega_dot);
387 std::vector<PairScalars> domega_dot_dT_reg(n_species, example);
389 for (
unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
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]);
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]);
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;
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;
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;
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;
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;
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 ]);
450 vec_compare(omega_dot, omega_dot_reg,
"omega_dot");
453 vec_compare(omega_dot_2, omega_dot_reg,
"omega_dot_2");
456 vec_compare(domega_dot_dT, domega_dot_dT_reg,
"domega_dot_dT");
458 char vecname[] =
"domega_dot_drho_0";
459 for (
unsigned int i = 0; i != 5; ++i)
462 vec_compare(domega_dot_drhos[i], domega_dot_drhos_reg[i], vecname);
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.
This class contains the conditions of the chemistry.
int vec_compare(const SpeciesVector1 &a, const SpeciesVector2 &b, const std::string &name)