26 #include "grins_config.h"
28 #ifdef GRINS_HAVE_ANTIOCH
41 #include "libmesh/getpot.h"
44 #include "antioch/cea_evaluator.h"
45 #include "antioch/stat_mech_thermo.h"
48 #include "boost/math/special_functions/fpclassify.hpp"
50 int test_generic(
const libMesh::Real value,
const libMesh::Real value_reg,
const std::string& name )
54 const double tol = std::numeric_limits<double>::epsilon()*40;
56 const double rel_error = std::fabs( (value - value_reg)/value_reg );
61 std::cout <<
"Mismatch in "+name << std::endl
62 << name+
" = " << value << std::endl
63 << name+
"_reg = " << value_reg << std::endl
64 <<
"rel_error = " << rel_error << std::endl;
70 template<
typename Thermo>
71 int test_cp(
const libMesh::Real cp );
73 template<
typename Thermo>
74 int test_cv(
const libMesh::Real cv );
76 template<
typename Thermo>
77 int test_h_s(
const libMesh::Real h_s );
80 int test_cp<Antioch::CEAEvaluator<libMesh::Real> >(
const libMesh::Real cp )
82 double cp_reg = 1.2361869971209990e+03;
88 int test_cv<Antioch::CEAEvaluator<libMesh::Real> >(
const libMesh::Real cv )
90 double cv_reg = 8.4681056933423179e+02;
96 int test_h_s<Antioch::CEAEvaluator<libMesh::Real> >(
const libMesh::Real h_s )
98 double h_s_reg = 7.6606764036494098e+05;
104 int test_cp<Antioch::StatMechThermodynamics<libMesh::Real> >(
const libMesh::Real cp )
106 double cp_reg = 1.2309250447693457e+03;
112 int test_cv<Antioch::StatMechThermodynamics<libMesh::Real> >(
const libMesh::Real cv )
114 double cv_reg = 8.4154861698257866e+02;
120 int test_h_s<Antioch::StatMechThermodynamics<libMesh::Real> >(
const libMesh::Real h_s )
122 double h_s_reg = 7.6397682090011111e+05;
128 template <
typename Thermo>
133 libMesh::Real T = 1000;
135 libMesh::Real rho = 1.0e-3;
137 const unsigned int n_species = 5;
139 std::vector<libMesh::Real> Y(n_species,0.2);
141 libMesh::Real R_mix = antioch_mixture.
R_mix(Y);
146 std::vector<double> Tqp(1,T);
150 std::vector<double> rhoqp(1,rho);
154 std::vector<double> Rqp(1,R_mix);
158 std::vector<std::vector<double> > Yqp(1,Y);
161 std::vector<libMesh::Real> omega_dot(n_species,0.0);
163 libMesh::Real cp = antioch_evaluator.
cp( cache, 0 );
165 std::cout << std::scientific << std::setprecision(16)
166 <<
"cp = " << cp << std::endl;
168 libMesh::Real cv = antioch_evaluator.
cv( cache, 0 );
170 std::cout << std::scientific << std::setprecision(16)
171 <<
"cv = " << cv << std::endl;
173 libMesh::Real h_s = antioch_evaluator.
h_s( cache, 0, 0 );
175 std::cout << std::scientific << std::setprecision(16)
176 <<
"h_s = " << h_s << std::endl;
179 int return_flag_temp = 0;
181 return_flag_temp = test_cp<Thermo>( cp );
182 if( return_flag_temp != 0 ) return_flag = 1;
184 return_flag_temp = test_cv<Thermo>( cv );
185 if( return_flag_temp != 0 ) return_flag = 1;
187 return_flag_temp = test_h_s<Thermo>( h_s );
188 if( return_flag_temp != 0 ) return_flag = 1;
190 antioch_evaluator.
omega_dot( cache, 0, omega_dot );
192 for(
unsigned int i = 0; i < n_species; i++ )
194 std::cout << std::scientific << std::setprecision(16)
195 <<
"omega_dot(" << i <<
") = " << omega_dot[i] << std::endl;
199 double tol = std::numeric_limits<double>::epsilon()*40;
202 libMesh::Real sum = 0.0;
203 for(
unsigned int s = 0; s < n_species; s++ )
208 if( std::fabs(sum) > tol )
210 std::cout <<
"Error: Sum of mass sources not equal to zero!" << std::endl
211 <<
"sum = " << sum << std::endl;
215 bool omega_is_nan =
false;
216 for(
unsigned int s = 0; s < n_species; s++ )
218 if( boost::math::isnan(omega_dot[s]) )
227 std::cerr <<
"Error: Detected NAN's for omega_dot!" << std::endl;
231 std::vector<libMesh::Real> omega_dot_reg(n_species,0.0);
232 omega_dot_reg[0] = 4.3561014567820650e-01;
233 omega_dot_reg[1] = -3.6797207156306633e+00;
234 omega_dot_reg[2] = 2.9970582875850726e+00;
235 omega_dot_reg[3] = -1.8346634812051223e+00;
236 omega_dot_reg[4] = 2.0817157635725065e+00;
238 for(
unsigned int s = 0; s < n_species; s++ )
240 if( std::fabs( (omega_dot[s] - omega_dot_reg[s])/omega_dot_reg[s] ) > tol )
242 std::cerr <<
"Error: Mismatch in omega_dot." << std::endl
243 << std::setprecision(16) << std::scientific
244 <<
"s = " << s << std::endl
245 <<
"omega_dot = " << omega_dot[s] << std::endl
246 <<
"omega_dot_reg = " << omega_dot_reg[s] << std::endl;
255 int main(
int argc,
char* argv[] )
261 std::cerr <<
"Error: Must specify input file." << std::endl;
265 GetPot input( argv[1] );
271 std::cout <<
"Running AntiochCEAThermo regression test." << std::endl;
272 return_flag = test_evaluator<Antioch::CEAEvaluator<libMesh::Real> >( antioch_mixture );
274 std::cout << std::endl <<
"Running AntiochStatMechThermo regression test." << std::endl;
275 return_flag = test_evaluator<Antioch::StatMechThermodynamics<libMesh::Real> >( antioch_mixture );
280 #else //GRINS_HAVE_ANTIOCH
libMesh::Real cv(const CachedValues &cache, unsigned int qp)
Wrapper class for evaluating chemistry and thermo properties using Antioch.
libMesh::Real R_mix(const std::vector< libMesh::Real > &mass_fractions) const
Mixture gas constant, [J/kg-K].
void set_values(unsigned int quantity, std::vector< libMesh::Number > &values)
int main(int argc, char *argv[])
libMesh::Real cp(const CachedValues &cache, unsigned int qp)
int test_generic(const libMesh::Real value, const libMesh::Real value_reg, const std::string &name)
int test_h_s(const libMesh::Real h_s)
void omega_dot(const CachedValues &cache, unsigned int qp, std::vector< libMesh::Real > &omega_dot)
libMesh::Real h_s(const CachedValues &cache, unsigned int qp, unsigned int species)
int test_cv(const libMesh::Real cv)
int test_cp(const libMesh::Real cp)
int test_evaluator(const GRINS::AntiochMixture &antioch_mixture)
Wrapper class for storing state for Antioch thermo and kinetics.
void add_quantity(unsigned int quantity)
void set_vector_values(unsigned int quantity, std::vector< std::vector< libMesh::Number > > &values)