32 #include "antioch_config.h"
36 #ifdef ANTIOCH_HAVE_EIGEN
37 #include "Eigen/Dense"
40 #ifdef ANTIOCH_HAVE_METAPHYSICL
41 #include "metaphysicl/numberarray.h"
44 #ifdef ANTIOCH_HAVE_VEXCL
45 #include "vexcl/vexcl.hpp"
72 #ifdef ANTIOCH_HAVE_GRVY
75 GRVY::GRVY_Timer_Class gt;
83 template <
typename Scalar>
86 const std::string& species_name,
87 Scalar molar_mass, Scalar gas_constant, Scalar formation_enthalpy,
88 Scalar n_tr_dofs,
int charge )
94 const Scalar tol = (std::numeric_limits<Scalar>::epsilon() * 10 < 5e-17)?5e-17:
95 std::numeric_limits<Scalar>::epsilon() * 10;
97 if( chem_species.
species() != species_name )
99 std::cerr <<
"Error: Name mismatch for "<< species_name << std::endl
100 <<
"name = " << chem_species.
species() << std::endl;
104 if( std::abs(chem_species.
molar_mass() - molar_mass)/molar_mass > tol )
106 std::cerr <<
"Error: Molar mass mismatch for "<< species_name << std::endl
107 <<
"molar mass = " << chem_species.
molar_mass() << std::endl;
111 if( std::abs(chem_species.
gas_constant() - gas_constant)/gas_constant > tol )
113 std::cerr <<
"Error: Gas constant mismatch for "<< species_name << std::endl
114 <<
"gas constant = " << chem_species.
gas_constant() << std::endl;
118 if( std::abs(chem_species.
formation_enthalpy() - formation_enthalpy)/formation_enthalpy )
120 std::cerr <<
"Error: Formation enthalpy mismatch for "<< species_name << std::endl
125 if( chem_species.
n_tr_dofs() != n_tr_dofs )
127 std::cerr <<
"Error: Number translational DoFs mismatch for "<< species_name << std::endl
128 <<
"n_tr_dofs = " << chem_species.
n_tr_dofs() << std::endl;
132 if( chem_species.
charge() != charge )
134 std::cerr <<
"Error: Charge mismatch for "<< species_name << std::endl
135 <<
"charge = " << chem_species.
charge() << std::endl;
143 template <
typename PairScalars>
144 int vectester(
const PairScalars& example,
const std::string& testname)
150 const Scalar Mm_N = 14.008e-3L;
151 const Scalar Mm_O = 16.000e-3L;
152 const Scalar Mm_N2 = 2.L * Mm_N;
153 const Scalar Mm_O2 = 2.L * Mm_O;
154 const Scalar Mm_NO = Mm_N + Mm_O;
156 std::vector<std::string> species_str_list;
157 const unsigned int n_species = 5;
158 species_str_list.reserve(n_species);
159 species_str_list.push_back(
"N2" );
160 species_str_list.push_back(
"O2" );
161 species_str_list.push_back(
"N" );
162 species_str_list.push_back(
"O" );
163 species_str_list.push_back(
"NO" );
167 const std::map<std::string,Antioch::Species>& species_name_map = chem_mixture.
species_name_map();
169 const std::vector<Antioch::ChemicalSpecies<Scalar>*>& chemical_species = chem_mixture.
chemical_species();
170 const std::vector<Antioch::Species> species_list = chem_mixture.
species_list();
174 const Scalar tol = (std::numeric_limits<Scalar>::epsilon() * 10 < 5e-17)?5e-17:
175 std::numeric_limits<Scalar>::epsilon() * 10;
178 for(
unsigned int i = 0; i < n_species; i++ )
180 if( species_name_map.find( species_str_list[i] )->second != species_list[i] )
182 std::cerr <<
"Error: species name map and species list ordering mismatch" << std::endl
183 <<
"species_name_map = " << species_name_map.find( species_str_list[i] )->second
184 <<
", species_list = " << species_list[i] << std::endl;
190 for(
unsigned int i = 0; i < n_species; i++ )
192 if( species_inverse_name_map.find( species_list[i] )->second != species_str_list[i] )
194 std::cerr <<
"Error: species inverse name map and species list ordering mismatch" << std::endl
195 <<
"species_inverse_name_map = " << species_inverse_name_map.find( species_list[i] )->second
196 <<
", species_str_list = " << species_str_list[i] << std::endl;
203 unsigned int index = 0;
204 Scalar molar_mass = Mm_N2;
205 if( std::abs(molar_mass - chem_mixture.
M(index))/molar_mass > tol )
207 std::cerr <<
"Error: Molar mass inconsistency in mixture" << std::endl
208 <<
"molar mass = " << chem_mixture.
M(index) << std::endl;
211 return_flag =
test_species( index, chemical_species,
"N2", molar_mass,
212 Scalar(Antioch::Constants::R_universal<Scalar>()/molar_mass),
213 Scalar(0.0), Scalar(2.5), Scalar(0));
218 unsigned int index = 1;
219 Scalar molar_mass = Mm_O2;
220 if( std::abs(molar_mass - chem_mixture.
M(index))/molar_mass > tol )
222 std::cerr <<
"Error: Molar mass inconsistency in mixture" << std::endl
223 <<
"molar mass = " << chem_mixture.
M(index) << std::endl;
226 return_flag =
test_species( index, chemical_species,
"O2", molar_mass,
227 Scalar(Antioch::Constants::R_universal<Scalar>()/molar_mass),
228 Scalar(0.0), Scalar(2.5), Scalar(0));
233 unsigned int index = 2;
234 Scalar molar_mass = Mm_N;
235 if( std::abs(molar_mass - chem_mixture.
M(index))/molar_mass > tol )
237 std::cerr <<
"Error: Molar mass inconsistency in mixture" << std::endl
238 <<
"molar mass = " << chem_mixture.
M(index) << std::endl;
241 return_flag =
test_species( index, chemical_species,
"N", molar_mass,
242 Scalar(Antioch::Constants::R_universal<Scalar>()/molar_mass),
243 Scalar(3.3621610000e7), Scalar(1.5), Scalar(0));
248 unsigned int index = 3;
249 Scalar molar_mass = Mm_O;
250 if( std::abs(molar_mass - chem_mixture.
M(index))/molar_mass > tol )
252 std::cerr <<
"Error: Molar mass inconsistency in mixture" << std::endl
253 <<
"molar mass = " << chem_mixture.
M(index) << std::endl;
256 return_flag =
test_species( index, chemical_species,
"O", molar_mass,
257 Scalar(Antioch::Constants::R_universal<Scalar>()/molar_mass),
258 Scalar(1.5420000000e7), Scalar(1.5), Scalar(0));
263 unsigned int index = 4;
264 Scalar molar_mass = Mm_NO;
265 if( std::abs(molar_mass - chem_mixture.
M(index))/molar_mass > tol )
267 std::cerr <<
"Error: Molar mass inconsistency in mixture" << std::endl
268 <<
"molar mass = " << chem_mixture.
M(index) << std::endl;
271 return_flag =
test_species( index, chemical_species,
"NO", molar_mass,
272 Scalar(Antioch::Constants::R_universal<Scalar>()/molar_mass),
273 Scalar(2.9961230000e6), Scalar(2.5), Scalar(0));
276 std::vector<PairScalars> mass_fractions( n_species, example );
277 PairScalars R_exact = example;
278 PairScalars M_exact = example;
279 std::vector<PairScalars> X_exact(n_species, example);
281 for (
unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
283 mass_fractions[0][2*tuple ] = 0.25L;
284 mass_fractions[1][2*tuple ] = 0.25L;
285 mass_fractions[2][2*tuple ] = 0.25L;
286 mass_fractions[3][2*tuple ] = 0.25L;
287 mass_fractions[4][2*tuple ] = 0L;
288 mass_fractions[0][2*tuple+1] = 0.2L;
289 mass_fractions[1][2*tuple+1] = 0.2L;
290 mass_fractions[2][2*tuple+1] = 0.2L;
291 mass_fractions[3][2*tuple+1] = 0.2L;
292 mass_fractions[4][2*tuple+1] = 0.2L;
294 R_exact[2*tuple ] = Antioch::Constants::R_universal<Scalar>()*( 0.25L/Mm_N2 + 0.25L/Mm_O2 + 0.25L/Mm_N + 0.25L/Mm_O);
295 R_exact[2*tuple+1] = Antioch::Constants::R_universal<Scalar>()*( 0.2L/Mm_N2 + 0.2L/Mm_O2 + 0.2L/Mm_N + 0.2L/Mm_O + 0.2L/Mm_NO );
297 M_exact[2*tuple ] = 1.0L/( 0.25L*( 1.0L/Mm_N2 + 1.0L/Mm_O2 + 1.0L/Mm_N + 1.0L/Mm_O) );
298 M_exact[2*tuple+1] = 1.0L/( 0.2L*( 1.0L/Mm_N2 + 1.0L/Mm_O2 + 1.0L/Mm_N + 1.0L/Mm_O + 1.0L/Mm_NO) );
300 X_exact[0][2*tuple ] = 0.25L*M_exact[0]/Mm_N2;
301 X_exact[1][2*tuple ] = 0.25L*M_exact[0]/Mm_O2;
302 X_exact[2][2*tuple ] = 0.25L*M_exact[0]/Mm_N;
303 X_exact[3][2*tuple ] = 0.25L*M_exact[0]/Mm_O;
304 X_exact[4][2*tuple ] = 0L;
305 X_exact[0][2*tuple+1] = 0.2L*M_exact[1]/Mm_N2;
306 X_exact[1][2*tuple+1] = 0.2L*M_exact[1]/Mm_O2;
307 X_exact[2][2*tuple+1] = 0.2L*M_exact[1]/Mm_N;
308 X_exact[3][2*tuple+1] = 0.2L*M_exact[1]/Mm_O;
309 X_exact[4][2*tuple+1] = 0.2L*M_exact[1]/Mm_NO;
312 #ifdef ANTIOCH_HAVE_GRVY
313 const std::string testnormal = testname +
"-normal";
314 gt.BeginTimer(testnormal);
317 std::vector<PairScalars> X(n_species, example);
319 const PairScalars R = chem_mixture.
R(mass_fractions);
320 const PairScalars M = chem_mixture.
M(mass_fractions);
321 chem_mixture.X( M, mass_fractions, X );
323 #ifdef ANTIOCH_HAVE_GRVY
324 gt.EndTimer(testnormal);
327 const PairScalars rel_R_error =
328 abs( (R - R_exact)/R_exact);
331 std::cerr <<
"Error: Mismatch in mixture gas constant." << std::endl
332 << std::setprecision(16) << std::scientific
333 <<
"R = " << R << std::endl
334 <<
"R_exact = " << R_exact << std::endl;
338 const PairScalars rel_M_error =
339 abs( (M - M_exact)/M_exact);
342 std::cerr <<
"Error: Mismatch in mixture molar mass." << std::endl
343 << std::setprecision(16) << std::scientific
344 <<
"M = " << M << std::endl
345 <<
"M_exact = " << M_exact << std::endl;
349 for(
unsigned int s = 0; s < n_species; s++ )
351 const PairScalars rel_X_error =
352 abs( (X[s] - X_exact[s])/X_exact[s]);
355 std::cerr <<
"Error: Mismatch in mole fraction for species " << s << std::endl
356 << std::setprecision(16) << std::scientific
357 <<
"X = " << X[s] << std::endl
358 <<
"X_exact = " << X_exact[s] << std::endl;
364 #ifdef ANTIOCH_HAVE_EIGEN
366 typedef Eigen::Array<PairScalars,n_species,1> SpeciesVecEigenType;
368 SpeciesVecEigenType eigen_mass_fractions;
371 SpeciesVecEigenType eigen_X_exact;
373 for (
unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
375 eigen_mass_fractions[4][2*tuple] = 0L;
377 eigen_X_exact[0][2*tuple ] = 0.25L*M_exact[0]/28.016L;
378 eigen_X_exact[1][2*tuple ] = 0.25L*M_exact[0]/32.0L;
379 eigen_X_exact[2][2*tuple ] = 0.25L*M_exact[0]/14.008L;
380 eigen_X_exact[3][2*tuple ] = 0.25L*M_exact[0]/16.0L;
381 eigen_X_exact[4][2*tuple ] = 0L;
382 eigen_X_exact[0][2*tuple+1] = 0.2L*M_exact[1]/28.016L;
383 eigen_X_exact[1][2*tuple+1] = 0.2L*M_exact[1]/32.0L;
384 eigen_X_exact[2][2*tuple+1] = 0.2L*M_exact[1]/14.008L;
385 eigen_X_exact[3][2*tuple+1] = 0.2L*M_exact[1]/16.0L;
386 eigen_X_exact[4][2*tuple+1] = 0.2L*M_exact[1]/30.008L;
389 #ifdef ANTIOCH_HAVE_GRVY
390 const std::string testeigenA = testname +
"-eigenA";
391 gt.BeginTimer(testeigenA);
394 const PairScalars R_eigen = chem_mixture.
R(eigen_mass_fractions);
395 const PairScalars M_eigen = chem_mixture.
M(eigen_mass_fractions);
396 SpeciesVecEigenType eigen_X;
398 chem_mixture.X( M_eigen, eigen_mass_fractions, eigen_X );
400 #ifdef ANTIOCH_HAVE_GRVY
401 gt.EndTimer(testeigenA);
404 const Scalar tol = (std::numeric_limits<Scalar>::epsilon() * 10 < 5e-17)?5e-17:
405 std::numeric_limits<Scalar>::epsilon() * 10;
407 const PairScalars eigen_rel_R_error =
408 abs( (R_eigen - R_exact)/R_exact);
411 std::cerr <<
"Error: Mismatch in mixture gas constant." << std::endl
412 << std::setprecision(16) << std::scientific
413 <<
"R_eigen = " << R_eigen << std::endl
414 <<
"R_exact = " << R_exact << std::endl;
418 const PairScalars eigen_rel_M_error =
419 abs( (M_eigen - M_exact)/M_exact);
422 std::cerr <<
"Error: Mismatch in mixture molar mass." << std::endl
423 << std::setprecision(16) << std::scientific
424 <<
"M_eigen = " << M_eigen << std::endl
425 <<
"M_exact = " << M_exact << std::endl;
429 for(
unsigned int s = 0; s < n_species; s++ )
431 const PairScalars eigen_rel_X_error =
432 abs( (eigen_X[s] - X_exact[s])/X_exact[s]);
435 std::cerr <<
"Error: Mismatch in mole fraction for species " << s << std::endl
436 << std::setprecision(16) << std::scientific
437 <<
"X_eigen = " << eigen_X[s] << std::endl
438 <<
"X_exact = " << X_exact[s] << std::endl;
445 typedef Eigen::Matrix<PairScalars,Eigen::Dynamic,1> SpeciesVecEigenType;
447 SpeciesVecEigenType eigen_mass_fractions(n_species, 1);
450 SpeciesVecEigenType eigen_X_exact(n_species, 1);
453 for (
unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
455 eigen_mass_fractions[4][2*tuple] = 0L;
457 eigen_X_exact[0][2*tuple ] = 0.25L*M_exact[0]/28.016L;
458 eigen_X_exact[1][2*tuple ] = 0.25L*M_exact[0]/32.0L;
459 eigen_X_exact[2][2*tuple ] = 0.25L*M_exact[0]/14.008L;
460 eigen_X_exact[3][2*tuple ] = 0.25L*M_exact[0]/16.0L;
461 eigen_X_exact[4][2*tuple ] = 0L;
462 eigen_X_exact[0][2*tuple+1] = 0.2L*M_exact[1]/28.016L;
463 eigen_X_exact[1][2*tuple+1] = 0.2L*M_exact[1]/32.0L;
464 eigen_X_exact[2][2*tuple+1] = 0.2L*M_exact[1]/14.008L;
465 eigen_X_exact[3][2*tuple+1] = 0.2L*M_exact[1]/16.0L;
466 eigen_X_exact[4][2*tuple+1] = 0.2L*M_exact[1]/30.008L;
469 #ifdef ANTIOCH_HAVE_GRVY
470 const std::string testeigenV = testname +
"-eigenV";
471 gt.BeginTimer(testeigenV);
474 const PairScalars R_eigen = chem_mixture.
R(eigen_mass_fractions);
475 const PairScalars M_eigen = chem_mixture.
M(eigen_mass_fractions);
476 SpeciesVecEigenType eigen_X(n_species, 1);
478 chem_mixture.X( M_eigen, eigen_mass_fractions, eigen_X );
480 #ifdef ANTIOCH_HAVE_GRVY
481 gt.EndTimer(testeigenV);
484 const Scalar tol = (std::numeric_limits<Scalar>::epsilon() * 10 < 5e-17)?5e-17:
485 std::numeric_limits<Scalar>::epsilon() * 10;
487 const PairScalars eigen_rel_R_error =
488 abs( (R_eigen - R_exact)/R_exact);
491 std::cerr <<
"Error: Mismatch in mixture gas constant." << std::endl
492 << std::setprecision(16) << std::scientific
493 <<
"R_eigen = " << R_eigen << std::endl
494 <<
"R_exact = " << R_exact << std::endl;
498 const PairScalars eigen_rel_M_error =
499 abs( (M_eigen - M_exact)/M_exact);
502 std::cerr <<
"Error: Mismatch in mixture molar mass." << std::endl
503 << std::setprecision(16) << std::scientific
504 <<
"M_eigen = " << M_eigen << std::endl
505 <<
"M_exact = " << M_exact << std::endl;
509 for(
unsigned int s = 0; s < n_species; s++ )
511 const PairScalars eigen_rel_X_error =
512 abs( (eigen_X[s] - X_exact[s])/X_exact[s]);
515 std::cerr <<
"Error: Mismatch in mole fraction for species " << s << std::endl
516 << std::setprecision(16) << std::scientific
517 <<
"X_eigen = " << eigen_X[s] << std::endl
518 <<
"X_exact = " << X_exact[s] << std::endl;
523 #endif // ANTIOCH_HAVE_EIGEN
533 returnval = returnval ||
534 vectester (std::valarray<float>(2*ANTIOCH_N_TUPLES),
"valarray<float>");
535 returnval = returnval ||
536 vectester (std::valarray<double>(2*ANTIOCH_N_TUPLES),
"valarray<double>");
537 returnval = returnval ||
538 vectester (std::valarray<long double>(2*ANTIOCH_N_TUPLES),
"valarray<ld>");
539 #ifdef ANTIOCH_HAVE_EIGEN
540 returnval = returnval ||
541 vectester (Eigen::Array<float, 2*ANTIOCH_N_TUPLES, 1>(),
"Eigen::ArrayXf");
542 returnval = returnval ||
543 vectester (Eigen::Array<double, 2*ANTIOCH_N_TUPLES, 1>(),
"Eigen::ArrayXd");
544 returnval = returnval ||
545 vectester (Eigen::Array<long double, 2*ANTIOCH_N_TUPLES, 1>(),
"Eigen::ArrayXld");
547 #ifdef ANTIOCH_HAVE_METAPHYSICL
548 returnval = returnval ||
549 vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0),
"NumberArray<float>");
550 returnval = returnval ||
551 vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0),
"NumberArray<double>");
552 returnval = returnval ||
553 vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0),
"NumberArray<ld>");
555 #ifdef ANTIOCH_HAVE_VEXCL
556 vex::Context ctx_f (vex::Filter::All);
558 returnval = returnval ||
559 vectester (vex::vector<float> (ctx_f, 2*ANTIOCH_N_TUPLES),
"vex::vector<float>");
561 vex::Context ctx_d (vex::Filter::DoublePrecision);
563 returnval = returnval ||
564 vectester (vex::vector<double> (ctx_d, 2*ANTIOCH_N_TUPLES),
"vex::vector<double>");
567 #ifdef ANTIOCH_HAVE_GRVY
Antioch::enable_if_c< is_eigen< T >::value, typename value_type< T >::type >::type max(const T &in)
CoeffType formation_enthalpy() const
Returns formation enthalpy in units of [J/kg].
CoeffType R(const unsigned int s) const
Gas constant for species s in [J/kg-K].
const std::vector< Species > & species_list() const
const std::string & species() const
Returns a descriptive name for this species.
CoeffType M(const unsigned int s) const
Molecular weight (molar mass) for species s in kg/mol.
const std::map< std::string, Species > & species_name_map() const
int vectester(const PairScalars &example, const std::string &testname)
Class to encapsulate data for each chemical species.
const std::map< Species, std::string > & species_inverse_name_map() const
int charge() const
Returns electrical charge number.
int test_species(const unsigned int species, const std::vector< Antioch::ChemicalSpecies< Scalar > * > &chemical_species, const std::string &species_name, Scalar molar_mass, Scalar gas_constant, Scalar formation_enthalpy, Scalar n_tr_dofs, int charge)
CoeffType n_tr_dofs() const
Returns number of translational degrees of freedom.
CoeffType gas_constant() const
Returns the species ideal gas constant in [J/kg-K].
Class storing chemical mixture properties.
void init_constant(Vector &output, const Scalar &example)
CoeffType molar_mass() const
Returns the molar mass in (kg/mol)
const std::vector< ChemicalSpecies< CoeffType > * > & chemical_species() const