From 38d41d8be2c3a0e288143cc2c1aca6ce54d75a14 Mon Sep 17 00:00:00 2001 From: Sylvain Plessis Date: Wed, 15 Apr 2015 16:12:21 -0500 Subject: [PATCH 1/4] Some forgotten overload with the KineticsConditions<> object --- src/kinetics/include/antioch/arrhenius_rate.h | 36 +++++++++++++++++++ src/kinetics/include/antioch/berthelot_rate.h | 35 ++++++++++++++++++ src/kinetics/include/antioch/constant_rate.h | 35 ++++++++++++++++++ .../include/antioch/hercourtessen_rate.h | 36 +++++++++++++++++++ 4 files changed, 142 insertions(+) diff --git a/src/kinetics/include/antioch/arrhenius_rate.h b/src/kinetics/include/antioch/arrhenius_rate.h index 0ed8a98d..d9a74020 100644 --- a/src/kinetics/include/antioch/arrhenius_rate.h +++ b/src/kinetics/include/antioch/arrhenius_rate.h @@ -136,6 +136,30 @@ namespace Antioch template void rate_and_derivative(const StateType& T, StateType& rate, StateType& drate_dT) const; +// KineticsConditions overloads + + //! \return the rate evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + rate(const KineticsConditions& T) const + ANTIOCH_AUTOFUNC(StateType, _Cf * ant_exp(- _Ea/T.T())) + + //! \return the rate evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + operator()(const KineticsConditions& T) const + ANTIOCH_AUTOFUNC(StateType, this->rate(T)) + + //! \return the derivative with respect to temperature evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + derivative( const KineticsConditions& T ) const + ANTIOCH_AUTOFUNC(StateType, (*this)(T) * (_Ea/T.temp_cache().T2) ) + + //! Simultaneously evaluate the rate and its derivative at \p T. + template + void rate_and_derivative(const KineticsConditions& T, StateType& rate, StateType& drate_dT) const; + //! print equation const std::string numeric() const; @@ -307,6 +331,18 @@ namespace Antioch return; } + template + template + inline + void ArrheniusRate::rate_and_derivative(const KineticsConditions& T, + StateType& rate, + StateType& drate_dT) const + { + rate = (*this)(T); + drate_dT = rate * _Ea/(T.temp_cache().T2); + return; + } + } // end namespace Antioch #endif // ANTIOCH_ARRHENIUS_RATE_H diff --git a/src/kinetics/include/antioch/berthelot_rate.h b/src/kinetics/include/antioch/berthelot_rate.h index 46166bc0..37c9e383 100644 --- a/src/kinetics/include/antioch/berthelot_rate.h +++ b/src/kinetics/include/antioch/berthelot_rate.h @@ -114,6 +114,30 @@ namespace Antioch template void rate_and_derivative(const StateType& T, StateType& rate, StateType& drate_dT) const; +// KineticsConditions overloads + + //! \return the rate evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + rate(const KineticsConditions& T) const + ANTIOCH_AUTOFUNC(StateType, _Cf * ant_exp(_D * T.T()) ) + + //! \return the rate evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + operator()(const KineticsConditions& T) const + ANTIOCH_AUTOFUNC(StateType, this->rate(T)) + + //! \return the derivative with respect to temperature evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + derivative( const KineticsConditions& T ) const + ANTIOCH_AUTOFUNC(StateType, (*this)(T) * _D ) + + //! Simultaneously evaluate the rate and its derivative at \p T. + template + void rate_and_derivative(const KineticsConditions& T, StateType& rate, StateType& drate_dT) const; + //! print equation const std::string numeric() const; }; @@ -249,6 +273,17 @@ namespace Antioch return; } + template + template + inline + void BerthelotRate::rate_and_derivative(const KineticsConditions& T, + StateType& rate, StateType& drate_dT) const + { + rate = (*this)(T); + drate_dT = rate*_D; + return; + } + } // end namespace Antioch #endif // ANTIOCH_BERTHELOT_RATE_H diff --git a/src/kinetics/include/antioch/constant_rate.h b/src/kinetics/include/antioch/constant_rate.h index 0b680558..da88a484 100644 --- a/src/kinetics/include/antioch/constant_rate.h +++ b/src/kinetics/include/antioch/constant_rate.h @@ -108,6 +108,29 @@ namespace Antioch template void rate_and_derivative(const StateType& T, StateType& rate, StateType& drate_dT) const; + // KineticsConditions overloads + //! \return the rate evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + rate(const KineticsConditions& cond) const + ANTIOCH_AUTOFUNC(StateType, constant_clone(cond.T(),_Cf)) + + //! \return the rate evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + operator()(const KineticsConditions& cond) const + ANTIOCH_AUTOFUNC(StateType, this->rate(cond)) + + //! \return the derivative with respect to temperature evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + derivative( const KineticsConditions& cond ) const + ANTIOCH_AUTOFUNC(StateType, zero_clone(cond.T())) + + //! Simultaneously evaluate the rate and its derivative at \p T. + template + void rate_and_derivative(const KineticsConditions& cond, StateType& rate, StateType& drate_dT) const; + //! print equation const std::string numeric() const; @@ -193,6 +216,18 @@ namespace Antioch return; } + template + template + inline + void ConstantRate::rate_and_derivative(const KineticsConditions& /*cond*/, + StateType& rate, + StateType& drate_dT) const + { + Antioch::constant_fill(rate, _Cf); + Antioch::set_zero(drate_dT); + return; + } + } // end namespace Antioch #endif // ANTIOCH_HERCOURT_ESSEN_RATE_H diff --git a/src/kinetics/include/antioch/hercourtessen_rate.h b/src/kinetics/include/antioch/hercourtessen_rate.h index 96440e05..3753c416 100644 --- a/src/kinetics/include/antioch/hercourtessen_rate.h +++ b/src/kinetics/include/antioch/hercourtessen_rate.h @@ -123,6 +123,30 @@ namespace Antioch template void rate_and_derivative(const StateType& T, StateType& rate, StateType& drate_dT) const; +// KineticsConditions overloads + + //! \return the rate evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + rate(const KineticsConditions& T) const + ANTIOCH_AUTOFUNC(StateType, _Cf * ant_exp(_eta * T.temp_cache().lnT)) + + //! \return the rate evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + operator()(const KineticsConditions& T) const + ANTIOCH_AUTOFUNC(StateType, this->rate(T)) + + //! \return the derivative with respect to temperature evaluated at \p T. + template + ANTIOCH_AUTO(StateType) + derivative( const KineticsConditions& T ) const + ANTIOCH_AUTOFUNC(StateType, (*this)(T)/T.T()*(_eta)) + + //! Simultaneously evaluate the rate and its derivative at \p T. + template + void rate_and_derivative(const KineticsConditions& T, StateType& rate, StateType& drate_dT) const; + //! print equation const std::string numeric() const; @@ -298,6 +322,18 @@ namespace Antioch return; } + template + template + inline + void HercourtEssenRate::rate_and_derivative(const KineticsConditions& cond, + StateType& rate, + StateType& drate_dT) const + { + rate = (*this)(cond); + drate_dT = rate/cond.T() * _eta; + return; + } + } // end namespace Antioch #endif // ANTIOCH_HERCOURT_ESSEN_RATE_H From 1d617ccde651a90685bbfb232b151362cd557cb0 Mon Sep 17 00:00:00 2001 From: Sylvain Plessis Date: Wed, 15 Apr 2015 16:13:52 -0500 Subject: [PATCH 2/4] Update of all the kinetics model tests. Now all the tests are standardized, testing both KineticsConditions and StateType for the temperature. All the vectorizations are also tested. --- test/arrhenius_rate_unit.C | 94 ++++---- test/arrhenius_rate_vec_unit.C | 161 +++++++++----- test/berthelot_rate_unit.C | 92 ++++---- test/berthelot_rate_vec_unit.C | 221 +++++++++++++------ test/berthelothercourtessen_rate_unit.C | 91 ++++---- test/berthelothercourtessen_rate_vec_unit.C | 222 +++++++++++++------ test/constant_rate_unit.C | 87 ++++---- test/constant_rate_vec_unit.C | 224 ++++++++++++++------ test/hercourtessen_rate_unit.C | 99 ++++----- test/hercourtessen_rate_vec_unit.C | 222 +++++++++++++------ test/kooij_rate_unit.C | 92 ++++---- test/kooij_rate_vec_unit.C | 204 +++++++++++++----- test/vanthoff_rate_vec_unit.C | 208 +++++++++++++----- 13 files changed, 1356 insertions(+), 661 deletions(-) diff --git a/test/arrhenius_rate_unit.C b/test/arrhenius_rate_unit.C index 80792c25..849f3a90 100644 --- a/test/arrhenius_rate_unit.C +++ b/test/arrhenius_rate_unit.C @@ -37,67 +37,73 @@ #include "antioch/units.h" template -int test_values(const Scalar & Cf, const Scalar & Ea, const Scalar & R, const Antioch::ArrheniusRate & arrhenius_rate) +int check_rate_and_derivative(const Scalar & rate_exact, const Scalar & derive_exact, + const Scalar & rate, const Scalar & derive, const Scalar & T) { - using std::abs; - using std::exp; - int return_flag = 0; - - const Scalar tol = std::numeric_limits::epsilon() * 100; - - for(Scalar T = 300.1; T <= 2500.1; T += 10.) - { - const Scalar rate_exact = Cf*exp(-Ea/(R*T)); - const Scalar derive_exact = Ea/(R*T*T) * Cf * exp(-Ea/(R*T)); - - Scalar rate1 = arrhenius_rate(T); - Scalar deriveRate1 = arrhenius_rate.derivative(T); - Scalar rate; - Scalar deriveRate; - - arrhenius_rate.rate_and_derivative(T,rate,deriveRate); - - if( abs( (rate1 - rate_exact)/rate_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate values." << std::endl - << "T = " << T << " K" << std::endl - << "rate(T) = " << rate1 << std::endl - << "rate_exact = " << rate_exact << std::endl; - - return_flag = 1; - } + const Scalar tol = std::numeric_limits::epsilon() * 2; + int return_flag(0); if( abs( (rate - rate_exact)/rate_exact ) > tol ) { std::cout << std::scientific << std::setprecision(16) << "Error: Mismatch in rate values." << std::endl << "T = " << T << " K" << std::endl << "rate(T) = " << rate << std::endl - << "rate_exact = " << rate_exact << std::endl; - + << "rate_exact = " << rate_exact << std::endl + << "relative difference = " << abs( (rate - rate_exact)/rate_exact ) << std::endl + << "tolerance = " << tol << std::endl; + return_flag = 1; } - if( abs( (deriveRate1 - derive_exact)/derive_exact ) > tol ) + if( abs( (derive - derive_exact)/derive_exact ) > tol ) { std::cout << std::scientific << std::setprecision(16) << "Error: Mismatch in rate derivative values." << std::endl << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate1 << std::endl - << "derive_exact = " << derive_exact << std::endl; + << "drate_dT(T) = " << derive << std::endl + << "derive_exact = " << derive_exact << std::endl + << "relative difference = " << abs( (derive - derive_exact)/derive_exact ) << std::endl + << "tolerance = " << tol << std::endl; return_flag = 1; } - if( abs( (deriveRate - derive_exact)/derive_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate << std::endl - << "derive_exact = " << derive_exact << std::endl; - return_flag = 1; - } - if(return_flag)break; + return return_flag; +} + +template +int test_values(const Scalar & Cf, const Scalar & Ea, const Scalar & R, const Antioch::ArrheniusRate & arrhenius_rate) +{ + using std::abs; + using std::exp; + int return_flag = 0; + + for(Scalar T = 300.1; T <= 2500.1; T += 10.) + { + const Scalar rate_exact = Cf*exp(-Ea/(R*T)); + const Scalar derive_exact = Ea/(R*T*T) * Cf * exp(-Ea/(R*T)); + Antioch::KineticsConditions cond(T); + +// KineticsConditions + + Scalar rate = arrhenius_rate(cond); + Scalar deriveRate = arrhenius_rate.derivative(cond); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + + arrhenius_rate.rate_and_derivative(cond,rate,deriveRate); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + +// T + + rate = arrhenius_rate(T); + deriveRate = arrhenius_rate.derivative(T); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + + arrhenius_rate.rate_and_derivative(T,rate,deriveRate); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; } return return_flag; diff --git a/test/arrhenius_rate_vec_unit.C b/test/arrhenius_rate_vec_unit.C index 6f0b9361..262e2b50 100644 --- a/test/arrhenius_rate_vec_unit.C +++ b/test/arrhenius_rate_vec_unit.C @@ -68,6 +68,69 @@ GRVY::GRVY_Timer_Class gt; #include #include +template +int check_rate_and_derivative(const PairScalars & rate_exact, const PairScalars & derive_exact, + const PairScalars & rate, const PairScalars & derive, const PairScalars & T) +{ + typedef typename Antioch::value_type::type Scalar; + const Scalar tol = std::numeric_limits::epsilon() * 2; + + int return_flag(0); + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + if( abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (rate[2*tuple+1] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple] << std::endl + << "derive_exact = " << derive_exact[2*tuple] << std::endl + << "relative difference = " << abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple+1] << std::endl + << "derive_exact = " << derive_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + + } + return return_flag; +} + template int vectester(const PairScalars& example, const std::string& testname) { @@ -83,74 +146,72 @@ int vectester(const PairScalars& example, const std::string& testname) // Construct from example to avoid resizing issues PairScalars T = example; + PairScalars rate_exact = example; + PairScalars derive_exact = example; for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) { T[2*tuple] = 1500.1; T[2*tuple+1] = 1600.1; + rate_exact[2*tuple] = Cf*exp(-Ea/1500.1); + rate_exact[2*tuple+1] = Cf*exp(-Ea/1600.1); + derive_exact[2*tuple] = Ea/(Scalar(1500.1)*Scalar(1500.1)) * Cf * exp(-Ea/Scalar(1500.1)); + derive_exact[2*tuple+1] = Ea/(Scalar(1600.1)*Scalar(1600.1)) * Cf * exp(-Ea/Scalar(1600.1)); } - - const Scalar rate_exact0 = Cf*exp(-Ea/1500.1); - const Scalar rate_exact1 = Cf*exp(-Ea/1600.1); - const Scalar derive_exact0 = Ea/(Scalar(1500.1)*Scalar(1500.1)) * Cf * exp(-Ea/Scalar(1500.1)); - const Scalar derive_exact1 = Ea/(Scalar(1600.1)*Scalar(1600.1)) * Cf * exp(-Ea/Scalar(1600.1)); + Antioch::KineticsConditions cond(T); int return_flag = 0; +// KineticsConditions #ifdef ANTIOCH_HAVE_GRVY gt.BeginTimer(testname); #endif - const PairScalars rate = arrhenius_rate(T); - const PairScalars deriveRate = arrhenius_rate.derivative(T); + PairScalars rate = arrhenius_rate(cond); + PairScalars derive = arrhenius_rate.derivative(cond); #ifdef ANTIOCH_HAVE_GRVY gt.EndTimer(testname); #endif - const Scalar tol = std::numeric_limits::epsilon()*10; + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; - for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) - { - if( abs( (rate[2*tuple] - rate_exact0)/rate_exact0 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T0) = " << rate[2*tuple] << std::endl - << "rate_exact = " << rate_exact0 << std::endl - << "difference = " << rate[2*tuple] - rate_exact0 << std::endl; - - return_flag = 1; - break; - } - - if( abs( (rate[2*tuple+1] - rate_exact1)/rate_exact1 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T1) = " << rate[2*tuple+1] << std::endl - << "rate_exact = " << rate_exact1 << std::endl - << "difference = " << rate[2*tuple+1] - rate_exact1 << std::endl; - - return_flag = 1; - break; - } - } - if( abs( (deriveRate[0] - derive_exact0)/derive_exact0 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T0) = " << deriveRate[0] << std::endl - << "derive_exact = " << derive_exact0 << std::endl; +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif - return_flag = 1; - } - if( abs( (deriveRate[1] - derive_exact1)/derive_exact1 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T1) = " << deriveRate[1] << std::endl - << "derive_exact = " << derive_exact1 << std::endl; + arrhenius_rate.rate_and_derivative(cond,rate,derive); - return_flag = 1; - } +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +// T +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + rate = arrhenius_rate(T); + derive = arrhenius_rate.derivative(T); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + arrhenius_rate.rate_and_derivative(T,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; std::cout << "Arrhenius rate: " << arrhenius_rate << std::endl; @@ -185,13 +246,11 @@ int main() vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray"); #endif #ifdef ANTIOCH_HAVE_VEXCL -std::cout << "vexcl start" << std::endl; vex::Context ctx_f (vex::Filter::All); -std::cout << "vexcl start" << std::endl; if (!ctx_f.empty()) returnval = returnval || vectester (vex::vector (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector"); -std::cout << "vexcl float ok" << std::endl; + vex::Context ctx_d (vex::Filter::DoublePrecision); if (!ctx_d.empty()) returnval = returnval || diff --git a/test/berthelot_rate_unit.C b/test/berthelot_rate_unit.C index 7e679920..2862c40c 100644 --- a/test/berthelot_rate_unit.C +++ b/test/berthelot_rate_unit.C @@ -35,66 +35,72 @@ #include "antioch/berthelot_rate.h" template -int test_values(const Scalar & Cf, const Scalar & D, const Antioch::BerthelotRate & berthelot_rate) +int check_rate_and_derivative(const Scalar & rate_exact, const Scalar & derive_exact, + const Scalar & rate, const Scalar & derive, const Scalar & T) { - using std::abs; - using std::exp; - int return_flag = 0; - const Scalar tol = std::numeric_limits::epsilon() * 100; - - for(Scalar T = 300.1; T <= 2500.1; T += 10.) - { - const Scalar rate_exact = Cf*exp(D*T); - const Scalar derive_exact = D * Cf * exp(D*T); - - Scalar rate1 = berthelot_rate(T); - Scalar deriveRate1 = berthelot_rate.derivative(T); - Scalar rate; - Scalar deriveRate; - - berthelot_rate.rate_and_derivative(T,rate,deriveRate); - - if( abs( (rate1 - rate_exact)/rate_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate values." << std::endl - << "T = " << T << " K" << std::endl - << "rate(T) = " << rate1 << std::endl - << "rate_exact = " << rate_exact << std::endl; - - return_flag = 1; - } + const Scalar tol = std::numeric_limits::epsilon() * 2; + int return_flag(0); if( abs( (rate - rate_exact)/rate_exact ) > tol ) { std::cout << std::scientific << std::setprecision(16) << "Error: Mismatch in rate values." << std::endl << "T = " << T << " K" << std::endl << "rate(T) = " << rate << std::endl - << "rate_exact = " << rate_exact << std::endl; - + << "rate_exact = " << rate_exact << std::endl + << "relative difference = " << abs( (rate - rate_exact)/rate_exact ) << std::endl + << "tolerance = " << tol << std::endl; + return_flag = 1; } - if( abs( (deriveRate1 - derive_exact)/derive_exact ) > tol ) + if( abs( (derive - derive_exact)/derive_exact ) > tol ) { std::cout << std::scientific << std::setprecision(16) << "Error: Mismatch in rate derivative values." << std::endl << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate1 << std::endl - << "derive_exact = " << derive_exact << std::endl; + << "drate_dT(T) = " << derive << std::endl + << "derive_exact = " << derive_exact << std::endl + << "relative difference = " << abs( (derive - derive_exact)/derive_exact ) << std::endl + << "tolerance = " << tol << std::endl; return_flag = 1; } - if( abs( (deriveRate - derive_exact)/derive_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate << std::endl - << "derive_exact = " << derive_exact << std::endl; - return_flag = 1; - } - if(return_flag)break; + return return_flag; +} + +template +int test_values(const Scalar & Cf, const Scalar & D, const Antioch::BerthelotRate & berthelot_rate) +{ + using std::abs; + using std::exp; + int return_flag = 0; + + for(Scalar T = 300.1; T <= 2500.1; T += 10.) + { + const Scalar rate_exact = Cf*exp(D*T); + const Scalar derive_exact = D * Cf * exp(D*T); + Antioch::KineticsConditions cond(T); + +//KineticsConditions + Scalar rate = berthelot_rate(cond); + Scalar deriveRate = berthelot_rate.derivative(cond); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + + berthelot_rate.rate_and_derivative(cond,rate,deriveRate); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + +// T + rate = berthelot_rate(T); + deriveRate = berthelot_rate.derivative(T); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + + berthelot_rate.rate_and_derivative(T,rate,deriveRate); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + } return return_flag; } diff --git a/test/berthelot_rate_vec_unit.C b/test/berthelot_rate_vec_unit.C index c2cc7e78..e037062b 100644 --- a/test/berthelot_rate_vec_unit.C +++ b/test/berthelot_rate_vec_unit.C @@ -43,6 +43,10 @@ #include "metaphysicl/numberarray.h" #endif +#ifdef ANTIOCH_HAVE_VEXCL +#include "vexcl/vexcl.hpp" +#endif + #include "antioch/eigen_utils_decl.h" #include "antioch/metaphysicl_utils_decl.h" #include "antioch/valarray_utils_decl.h" @@ -56,8 +60,77 @@ #include #include +#ifdef ANTIOCH_HAVE_GRVY +#include "grvy.h" + +GRVY::GRVY_Timer_Class gt; +#endif + +template +int check_rate_and_derivative(const PairScalars & rate_exact, const PairScalars & derive_exact, + const PairScalars & rate, const PairScalars & derive, const PairScalars & T) +{ + typedef typename Antioch::value_type::type Scalar; + const Scalar tol = std::numeric_limits::epsilon() * 2; + + int return_flag(0); + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + if( abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (rate[2*tuple+1] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple] << std::endl + << "derive_exact = " << derive_exact[2*tuple] << std::endl + << "relative difference = " << abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple+1] << std::endl + << "derive_exact = " << derive_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + + } + return return_flag; +} + template -int vectester(const PairScalars& example) +int vectester(const PairScalars& example, const std::string & testname) { using std::abs; using std::exp; @@ -71,58 +144,72 @@ int vectester(const PairScalars& example) // Construct from example to avoid resizing issues PairScalars T = example; - T[0] = 1500.1; - T[1] = 1600.1; - - const Scalar rate_exact0 = Cf*exp(D*1500.1); - const Scalar rate_exact1 = Cf*exp(D*1600.1); - const Scalar derive_exact0 = D * Cf * exp(D*Scalar(1500.1)); - const Scalar derive_exact1 = D * Cf * exp(D*Scalar(1600.1)); + PairScalars rate_exact = example; + PairScalars derive_exact = example; + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + T[2*tuple] = 1500.1; + T[2*tuple+1] = 1600.1; + rate_exact[2*tuple] = Cf*exp(D*1500.1); + rate_exact[2*tuple+1] = Cf*exp(D*1600.1); + derive_exact[2*tuple] = D * Cf * exp(D*Scalar(1500.1)); + derive_exact[2*tuple+1] = D * Cf * exp(D*Scalar(1600.1)); + } + Antioch::KineticsConditions cond(T); int return_flag = 0; - const PairScalars rate = berthelot_rate(T); - const PairScalars deriveRate = berthelot_rate.derivative(T); - - const Scalar tol = std::numeric_limits::epsilon()*10; - - if( abs( (rate[0] - rate_exact0)/rate_exact0 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T0) = " << rate[0] << std::endl - << "rate_exact = " << rate_exact0 << std::endl - << "difference = " << rate[0] - rate_exact0 << std::endl; - - return_flag = 1; - } - - if( abs( (rate[1] - rate_exact1)/rate_exact1 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T1) = " << rate[1] << std::endl - << "rate_exact = " << rate_exact1 << std::endl - << "difference = " << rate[1] - rate_exact1 << std::endl; - - return_flag = 1; - } - if( abs( (deriveRate[0] - derive_exact0)/derive_exact0 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T0) = " << deriveRate[0] << std::endl - << "derive_exact = " << derive_exact0 << std::endl; - - return_flag = 1; - } - if( abs( (deriveRate[1] - derive_exact1)/derive_exact1 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T1) = " << deriveRate[1] << std::endl - << "derive_exact = " << derive_exact1 << std::endl; - - return_flag = 1; - } +// KineticsConditions +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + PairScalars rate = berthelot_rate(cond); + PairScalars derive = berthelot_rate.derivative(cond); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + berthelot_rate.rate_and_derivative(cond,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +// T +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + rate = berthelot_rate(T); + derive = berthelot_rate.derivative(T); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + berthelot_rate.rate_and_derivative(T,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; std::cout << "Berthelot rate: " << berthelot_rate << std::endl; @@ -135,26 +222,42 @@ int main() int returnval = 0; returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); #ifdef ANTIOCH_HAVE_EIGEN returnval = returnval || - vectester (Eigen::Array2f()); + vectester (Eigen::Array(), "Eigen::ArrayXf"); returnval = returnval || - vectester (Eigen::Array2d()); + vectester (Eigen::Array(), "Eigen::ArrayXd"); returnval = returnval || - vectester (Eigen::Array()); + vectester (Eigen::Array(), "Eigen::ArrayXld"); #endif #ifdef ANTIOCH_HAVE_METAPHYSICL returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, float> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, long double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray"); +#endif +#ifdef ANTIOCH_HAVE_VEXCL + vex::Context ctx_f (vex::Filter::All); + if (!ctx_f.empty()) + returnval = returnval || + vectester (vex::vector (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector"); + + vex::Context ctx_d (vex::Filter::DoublePrecision); + if (!ctx_d.empty()) + returnval = returnval || + vectester (vex::vector (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector"); +#endif + +#ifdef ANTIOCH_HAVE_GRVY + gt.Finalize(); + gt.Summarize(); #endif return returnval; diff --git a/test/berthelothercourtessen_rate_unit.C b/test/berthelothercourtessen_rate_unit.C index c370aa2f..4b0c014d 100644 --- a/test/berthelothercourtessen_rate_unit.C +++ b/test/berthelothercourtessen_rate_unit.C @@ -35,6 +35,40 @@ #include "antioch/berthelothercourtessen_rate.h" +template +int check_rate_and_derivative(const Scalar & rate_exact, const Scalar & derive_exact, + const Scalar & rate, const Scalar & derive, const Scalar & T) +{ + const Scalar tol = std::numeric_limits::epsilon() * 2; + int return_flag(0); + if( abs( (rate - rate_exact)/rate_exact ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate << std::endl + << "rate_exact = " << rate_exact << std::endl + << "relative difference = " << abs( (rate - rate_exact)/rate_exact ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive - derive_exact)/derive_exact ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive << std::endl + << "derive_exact = " << derive_exact << std::endl + << "relative difference = " << abs( (derive - derive_exact)/derive_exact ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + + return return_flag; +} + template int test_values(const Scalar & Cf, const Scalar & eta, const Scalar & D, const Scalar & Tref, const Antioch::BerthelotHercourtEssenRate & berthelothercourtessen_rate) { @@ -42,62 +76,33 @@ int test_values(const Scalar & Cf, const Scalar & eta, const Scalar & D, const S using std::exp; using std::pow; int return_flag = 0; - const Scalar tol = std::numeric_limits::epsilon() * 100; for(Scalar T = 300.1L; T <= 2500.1L; T += 10.L) { const Scalar rate_exact = Cf*pow(T/Tref,eta)*exp(D*T); const Scalar derive_exact = Cf * pow(T/Tref,eta) * exp(D*T) * (D + eta/T); + Antioch::KineticsConditions cond(T); - Scalar rate1 = berthelothercourtessen_rate(T); - Scalar deriveRate1 = berthelothercourtessen_rate.derivative(T); - Scalar rate; - Scalar deriveRate; +// KineticsConditions + Scalar rate = berthelothercourtessen_rate(cond); + Scalar deriveRate = berthelothercourtessen_rate.derivative(cond); - berthelothercourtessen_rate.rate_and_derivative(T,rate,deriveRate); + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; - if( abs( (rate1 - rate_exact)/rate_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate values." << std::endl - << "T = " << T << " K" << std::endl - << "rate(T) = " << rate1 << std::endl - << "rate_exact = " << rate_exact << std::endl; + berthelothercourtessen_rate.rate_and_derivative(cond,rate,deriveRate); - return_flag = 1; - } - if( abs( (rate - rate_exact)/rate_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate values." << std::endl - << "T = " << T << " K" << std::endl - << "rate(T) = " << rate << std::endl - << "rate_exact = " << rate_exact << std::endl; + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; - return_flag = 1; - } - if( abs( (deriveRate1 - derive_exact)/derive_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate1 << std::endl - << "derive_exact = " << derive_exact << std::endl; +// T + rate = berthelothercourtessen_rate(T); + deriveRate = berthelothercourtessen_rate.derivative(T); - return_flag = 1; - } - if( abs( (deriveRate - derive_exact)/derive_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate << std::endl - << "derive_exact = " << derive_exact << std::endl; + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; - return_flag = 1; - } + berthelothercourtessen_rate.rate_and_derivative(T,rate,deriveRate); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; - if(return_flag)break; } return return_flag; diff --git a/test/berthelothercourtessen_rate_vec_unit.C b/test/berthelothercourtessen_rate_vec_unit.C index d954b03a..f75a398a 100644 --- a/test/berthelothercourtessen_rate_vec_unit.C +++ b/test/berthelothercourtessen_rate_vec_unit.C @@ -43,6 +43,10 @@ #include "metaphysicl/numberarray.h" #endif +#ifdef ANTIOCH_HAVE_VEXCL +#include "vexcl/vexcl.hpp" +#endif + #include "antioch/eigen_utils_decl.h" #include "antioch/metaphysicl_utils_decl.h" #include "antioch/valarray_utils_decl.h" @@ -56,8 +60,77 @@ #include #include +#ifdef ANTIOCH_HAVE_GRVY +#include "grvy.h" + +GRVY::GRVY_Timer_Class gt; +#endif + template -int vectester(const PairScalars& example) +int check_rate_and_derivative(const PairScalars & rate_exact, const PairScalars & derive_exact, + const PairScalars & rate, const PairScalars & derive, const PairScalars & T) +{ + typedef typename Antioch::value_type::type Scalar; + const Scalar tol = std::numeric_limits::epsilon() * 2; + + int return_flag(0); + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + if( abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (rate[2*tuple+1] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple] << std::endl + << "derive_exact = " << derive_exact[2*tuple] << std::endl + << "relative difference = " << abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple+1] << std::endl + << "derive_exact = " << derive_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + + } + return return_flag; +} + +template +int vectester(const PairScalars& example, const std::string & testname) { using std::abs; using std::exp; @@ -73,58 +146,73 @@ int vectester(const PairScalars& example) // Construct from example to avoid resizing issues PairScalars T = example; - T[0] = 1500.1; - T[1] = 1600.1; - - const Scalar rate_exact0 = Cf*pow(Scalar(1500.1),eta)*exp(D*1500.1); - const Scalar rate_exact1 = Cf*pow(Scalar(1600.1),eta)*exp(D*1600.1); - const Scalar derive_exact0 = Cf * pow(Scalar(1500.1),eta) * exp(D*Scalar(1500.1)) * (D + eta/Scalar(1500.1)); - const Scalar derive_exact1 = Cf * pow(Scalar(1600.1),eta) * exp(D*Scalar(1600.1)) * (D + eta/Scalar(1600.1)); + PairScalars rate_exact = example; + PairScalars derive_exact = example; + + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + T[2*tuple] = 1500.1; + T[2*tuple+1] = 1600.1; + rate_exact[2*tuple] = Cf*pow(Scalar(1500.1),eta)*exp(D*1500.1); + rate_exact[2*tuple+1] = Cf*pow(Scalar(1600.1),eta)*exp(D*1600.1); + derive_exact[2*tuple] = Cf * pow(Scalar(1500.1),eta) * exp(D*Scalar(1500.1)) * (D + eta/Scalar(1500.1)); + derive_exact[2*tuple+1] = Cf * pow(Scalar(1600.1),eta) * exp(D*Scalar(1600.1)) * (D + eta/Scalar(1600.1)); + } + Antioch::KineticsConditions cond(T); int return_flag = 0; - const PairScalars rate = berthelothercourtessen_rate(T); - const PairScalars deriveRate = berthelothercourtessen_rate.derivative(T); - - const Scalar tol = std::numeric_limits::epsilon()*10; - - if( abs( (rate[0] - rate_exact0)/rate_exact0 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T0) = " << rate[0] << std::endl - << "rate_exact = " << rate_exact0 << std::endl - << "difference = " << rate[0] - rate_exact0 << std::endl; - - return_flag = 1; - } - - if( abs( (rate[1] - rate_exact1)/rate_exact1 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T1) = " << rate[1] << std::endl - << "rate_exact = " << rate_exact1 << std::endl - << "difference = " << rate[1] - rate_exact1 << std::endl; - - return_flag = 1; - } - if( abs( (deriveRate[0] - derive_exact0)/derive_exact0 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T0) = " << deriveRate[0] << std::endl - << "derive_exact = " << derive_exact0 << std::endl; - - return_flag = 1; - } - if( abs( (deriveRate[1] - derive_exact1)/derive_exact1 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T1) = " << deriveRate[1] << std::endl - << "derive_exact = " << derive_exact1 << std::endl; - - return_flag = 1; - } +// KineticsConditions +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + PairScalars rate = berthelothercourtessen_rate(cond); + PairScalars derive = berthelothercourtessen_rate.derivative(cond); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + berthelothercourtessen_rate.rate_and_derivative(cond,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + +// T +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + rate = berthelothercourtessen_rate(T); + derive = berthelothercourtessen_rate.derivative(T); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + berthelothercourtessen_rate.rate_and_derivative(T,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + + + std::cout << "Berthelot Hercourt Essen rate: " << berthelothercourtessen_rate << std::endl; @@ -137,26 +225,42 @@ int main() int returnval = 0; returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); #ifdef ANTIOCH_HAVE_EIGEN returnval = returnval || - vectester (Eigen::Array2f()); + vectester (Eigen::Array(), "Eigen::ArrayXf"); returnval = returnval || - vectester (Eigen::Array2d()); + vectester (Eigen::Array(), "Eigen::ArrayXd"); returnval = returnval || - vectester (Eigen::Array()); + vectester (Eigen::Array(), "Eigen::ArrayXld"); #endif #ifdef ANTIOCH_HAVE_METAPHYSICL returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, float> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, long double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray"); +#endif +#ifdef ANTIOCH_HAVE_VEXCL + vex::Context ctx_f (vex::Filter::All); + if (!ctx_f.empty()) + returnval = returnval || + vectester (vex::vector (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector"); + + vex::Context ctx_d (vex::Filter::DoublePrecision); + if (!ctx_d.empty()) + returnval = returnval || + vectester (vex::vector (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector"); +#endif + +#ifdef ANTIOCH_HAVE_GRVY + gt.Finalize(); + gt.Summarize(); #endif return returnval; diff --git a/test/constant_rate_unit.C b/test/constant_rate_unit.C index 7bdb77d5..e26fcec4 100644 --- a/test/constant_rate_unit.C +++ b/test/constant_rate_unit.C @@ -35,64 +35,71 @@ #include "antioch/constant_rate.h" template -int test_values(const Scalar & Cf, const Antioch::ConstantRate & constant_rate) +int check_rate_and_derivative(const Scalar & rate_exact, const Scalar & derive_exact, + const Scalar & rate, const Scalar & derive, const Scalar & T) { - const Scalar tol = std::numeric_limits::epsilon() * 100; - int return_flag = 0; - for(Scalar T = 300.1; T <= 2500.1; T += 10.) - { - - const Scalar rate_exact = Cf; - const Scalar derive_exact = 0.L; - - Scalar rate1 = constant_rate(T); - Scalar deriveRate1 = constant_rate.derivative(T); - Scalar rate; - Scalar deriveRate; - - constant_rate.rate_and_derivative(T,rate,deriveRate); - - if( abs( (rate1 - rate_exact)/rate_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate values." << std::endl - << "T = " << T << " K" << std::endl - << "rate(T) = " << rate1 << std::endl - << "rate_exact = " << rate_exact << std::endl; - - return_flag = 1; - } + const Scalar tol = std::numeric_limits::epsilon() * 2; + int return_flag(0); if( abs( (rate - rate_exact)/rate_exact ) > tol ) { std::cout << std::scientific << std::setprecision(16) << "Error: Mismatch in rate values." << std::endl << "T = " << T << " K" << std::endl << "rate(T) = " << rate << std::endl - << "rate_exact = " << rate_exact << std::endl; + << "rate_exact = " << rate_exact << std::endl + << "relative difference = " << abs( (rate - rate_exact)/rate_exact ) << std::endl + << "tolerance = " << tol << std::endl; return_flag = 1; } - if( abs( (deriveRate1 - derive_exact)/derive_exact ) > tol ) + if( abs( (derive - derive_exact)/derive_exact ) > tol ) { std::cout << std::scientific << std::setprecision(16) << "Error: Mismatch in rate derivative values." << std::endl << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate1 << std::endl - << "derive_exact = " << derive_exact << std::endl; + << "drate_dT(T) = " << derive << std::endl + << "derive_exact = " << derive_exact << std::endl + << "relative difference = " << abs( (derive - derive_exact)/derive_exact ) << std::endl + << "tolerance = " << tol << std::endl; return_flag = 1; } - if( abs( (deriveRate - derive_exact)/derive_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate << std::endl - << "derive_exact = " << derive_exact << std::endl; - return_flag = 1; - } - if(return_flag)break; + return return_flag; +} + + +template +int test_values(const Scalar & Cf, const Antioch::ConstantRate & constant_rate) +{ + int return_flag = 0; + for(Scalar T = 300.1; T <= 2500.1; T += 10.) + { + + const Scalar rate_exact = Cf; + const Scalar derive_exact = 0.L; + Antioch::KineticsConditions cond(T); + +// KineticsConditions + Scalar rate = constant_rate(cond); + Scalar deriveRate = constant_rate.derivative(cond); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + + constant_rate.rate_and_derivative(cond,rate,deriveRate); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + +// T + rate = constant_rate(T); + deriveRate = constant_rate.derivative(T); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + + constant_rate.rate_and_derivative(T,rate,deriveRate); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + } return return_flag; } diff --git a/test/constant_rate_vec_unit.C b/test/constant_rate_vec_unit.C index fe1de0bd..ee0fc468 100644 --- a/test/constant_rate_vec_unit.C +++ b/test/constant_rate_vec_unit.C @@ -44,6 +44,10 @@ #include "metaphysicl/numberarray.h" #endif +#ifdef ANTIOCH_HAVE_VEXCL +#include "vexcl/vexcl.hpp" +#endif + #include "antioch/eigen_utils_decl.h" #include "antioch/metaphysicl_utils_decl.h" #include "antioch/valarray_utils_decl.h" @@ -57,73 +61,153 @@ #include #include +#ifdef ANTIOCH_HAVE_GRVY +#include "grvy.h" + +GRVY::GRVY_Timer_Class gt; +#endif + template -int vectester(const PairScalars& example) +int check_rate_and_derivative(const PairScalars & rate_exact, const PairScalars & derive_exact, + const PairScalars & rate, const PairScalars & derive, const PairScalars & T) +{ + typedef typename Antioch::value_type::type Scalar; + const Scalar tol = std::numeric_limits::epsilon() * 2; + + int return_flag(0); + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + if( abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T[2*tuple] << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (rate[2*tuple+1] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T[2*tuple+1] << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T[2*tuple] << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple] << std::endl + << "derive_exact = " << derive_exact[2*tuple] << std::endl + << "relative difference = " << abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T[2*tuple+1] << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple+1] << std::endl + << "derive_exact = " << derive_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + + } + return return_flag; +} + +template +int vectester(const PairScalars& example, const std::string & testname) { using std::abs; typedef typename Antioch::value_type::type Scalar; - const Scalar Cf = 1.4; + const Scalar Cf = 1.4L; Antioch::ConstantRate constant_rate(Cf); // Construct from example to avoid resizing issues PairScalars T = example; - T[0] = 1500.1; - T[1] = 1600.1; + PairScalars rate_exact = example; + PairScalars derive_exact = Antioch::zero_clone(example); + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + T[2*tuple] = 1500.1L; + T[2*tuple+1] = 1600.1L; + rate_exact[2*tuple] = Cf; + rate_exact[2*tuple+1] = Cf; + } + Antioch::KineticsConditions cond(T); - const Scalar rate_exact0 = Cf; - const Scalar rate_exact1 = Cf; - const Scalar derive_exact0 = 0.L; - const Scalar derive_exact1 = 0.L; - int return_flag = 0; - const PairScalars rate = constant_rate(T);//Antioch::zero_clone(T); - const PairScalars deriveRate = constant_rate.derivative(T);//Antioch::zero_clone(T); - - const Scalar tol = std::numeric_limits::epsilon()*10; - -// hercourtessen_rate.rate_and_derivative(T,rate,deriveRate); - - if( abs( (rate[0] - rate_exact0)/rate_exact0 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T0) = " << rate[0] << std::endl - << "rate_exact = " << rate_exact0 << std::endl - << "difference = " << rate[0] - rate_exact0 << std::endl; - - return_flag = 1; - } - - if( abs( (rate[1] - rate_exact1)/rate_exact1 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T1) = " << rate[1] << std::endl - << "rate_exact = " << rate_exact1 << std::endl - << "difference = " << rate[1] - rate_exact1 << std::endl; - - return_flag = 1; - } - if( abs( (deriveRate[0] - derive_exact0)/derive_exact0 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T0) = " << deriveRate[0] << std::endl - << "derive_exact = " << derive_exact0 << std::endl; - - return_flag = 1; - } - if( abs( (deriveRate[1] - derive_exact1)/derive_exact1 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T1) = " << deriveRate[1] << std::endl - << "derive_exact = " << derive_exact1 << std::endl; - - return_flag = 1; - } +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + +// KineticsConditions + PairScalars rate = constant_rate(cond);//Antioch::zero_clone(T.T()); + PairScalars derive = constant_rate.derivative(cond);//Antioch::zero_clone(T.T()); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + constant_rate.rate_and_derivative(cond,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +// T + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + rate = constant_rate(T); + derive = constant_rate.derivative(T); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + constant_rate.rate_and_derivative(T,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; return return_flag; } @@ -134,26 +218,42 @@ int main() int returnval = 0; returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); #ifdef ANTIOCH_HAVE_EIGEN returnval = returnval || - vectester (Eigen::Array2f()); + vectester (Eigen::Array(), "Eigen::ArrayXf"); returnval = returnval || - vectester (Eigen::Array2d()); + vectester (Eigen::Array(), "Eigen::ArrayXd"); returnval = returnval || - vectester (Eigen::Array()); + vectester (Eigen::Array(), "Eigen::ArrayXld"); #endif #ifdef ANTIOCH_HAVE_METAPHYSICL returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, float> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, long double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray"); +#endif +#ifdef ANTIOCH_HAVE_VEXCL + vex::Context ctx_f (vex::Filter::All); + if (!ctx_f.empty()) + returnval = returnval || + vectester (vex::vector (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector"); + + vex::Context ctx_d (vex::Filter::DoublePrecision); + if (!ctx_d.empty()) + returnval = returnval || + vectester (vex::vector (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector"); +#endif + +#ifdef ANTIOCH_HAVE_GRVY + gt.Finalize(); + gt.Summarize(); #endif return returnval; diff --git a/test/hercourtessen_rate_unit.C b/test/hercourtessen_rate_unit.C index 10b8c145..aefcecc9 100644 --- a/test/hercourtessen_rate_unit.C +++ b/test/hercourtessen_rate_unit.C @@ -35,38 +35,11 @@ #include "antioch/hercourtessen_rate.h" template -int test_values(const Scalar & Cf, const Scalar & eta, const Scalar & Tref, const Antioch::HercourtEssenRate & hercourtessen_rate) +int check_rate_and_derivative(const Scalar & rate_exact, const Scalar & derive_exact, + const Scalar & rate, const Scalar & derive, const Scalar & T) { - using std::abs; - using std::pow; - int return_flag = 0; - const Scalar tol = std::numeric_limits::epsilon() * 100; - - for(Scalar T = 300.1L; T <= 2500.1L; T += 10.L) - { - const Scalar rate_exact = Cf * pow(T/Tref,eta); - const Scalar derive_exact = Cf * eta * pow(T/Tref,eta)/T; - - Scalar rate1 = hercourtessen_rate(T); - Scalar deriveRate1 = hercourtessen_rate.derivative(T); - Scalar rate; - Scalar deriveRate; - - hercourtessen_rate.rate_and_derivative(T,rate,deriveRate); - - if( abs( (rate1 - rate_exact)/rate_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate values." << std::endl - << "T = " << T << " K" << std::endl - << "rate(T) = " << rate1 << std::endl - << "rate_exact = " << rate_exact << std::endl - << "tol = " << tol << std::endl - << "relative error = " << abs( (rate1 - rate_exact)/rate_exact ) << std::endl - << "object: " << hercourtessen_rate << std::endl; - - return_flag = 1; - } + const Scalar tol = std::numeric_limits::epsilon() * 2; + int return_flag(0); if( abs( (rate - rate_exact)/rate_exact ) > tol ) { std::cout << std::scientific << std::setprecision(16) @@ -74,39 +47,61 @@ int test_values(const Scalar & Cf, const Scalar & eta, const Scalar & Tref, cons << "T = " << T << " K" << std::endl << "rate(T) = " << rate << std::endl << "rate_exact = " << rate_exact << std::endl - << "tol = " << tol << std::endl - << "relative error = " << abs( (rate1 - rate_exact)/rate_exact ) << std::endl - << "object: " << hercourtessen_rate << std::endl; + << "relative difference = " << abs( (rate - rate_exact)/rate_exact ) << std::endl + << "tolerance = " << tol << std::endl; return_flag = 1; } - if( abs( (deriveRate1 - derive_exact)/derive_exact ) > tol ) + if( abs( (derive - derive_exact)/derive_exact ) > tol ) { std::cout << std::scientific << std::setprecision(16) << "Error: Mismatch in rate derivative values." << std::endl << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate1 << std::endl + << "drate_dT(T) = " << derive << std::endl << "derive_exact = " << derive_exact << std::endl - << "tol = " << tol << std::endl - << "relative error = " << abs( (rate1 - rate_exact)/rate_exact ) << std::endl - << "object: " << hercourtessen_rate << std::endl; + << "relative difference = " << abs( (derive - derive_exact)/derive_exact ) << std::endl + << "tolerance = " << tol << std::endl; return_flag = 1; } - if( abs( (deriveRate - derive_exact)/derive_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate << std::endl - << "derive_exact = " << derive_exact << std::endl - << "tol = " << tol << std::endl - << "relative error = " << abs( (rate1 - rate_exact)/rate_exact ) << std::endl - << "object: " << hercourtessen_rate << std::endl; - return_flag = 1; - } - if(return_flag)break; + return return_flag; +} + +template +int test_values(const Scalar & Cf, const Scalar & eta, const Scalar & Tref, const Antioch::HercourtEssenRate & hercourtessen_rate) +{ + using std::abs; + using std::pow; + int return_flag = 0; + + for(Scalar T = 300.1L; T <= 2500.1L; T += 10.L) + { + const Scalar rate_exact = Cf * pow(T/Tref,eta); + const Scalar derive_exact = Cf * eta * pow(T/Tref,eta)/T; + Antioch::KineticsConditions cond(T); + + +// KineticsConditions + Scalar rate = hercourtessen_rate(cond); + Scalar deriveRate = hercourtessen_rate.derivative(cond); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + + hercourtessen_rate.rate_and_derivative(cond,rate,deriveRate); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + +// T + rate = hercourtessen_rate(T); + deriveRate = hercourtessen_rate.derivative(T); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + + hercourtessen_rate.rate_and_derivative(T,rate,deriveRate); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + } return return_flag; } diff --git a/test/hercourtessen_rate_vec_unit.C b/test/hercourtessen_rate_vec_unit.C index 8af89c79..e457b0da 100644 --- a/test/hercourtessen_rate_vec_unit.C +++ b/test/hercourtessen_rate_vec_unit.C @@ -43,6 +43,10 @@ #include "metaphysicl/numberarray.h" #endif +#ifdef ANTIOCH_HAVE_VEXCL +#include "vexcl/vexcl.hpp" +#endif + #include "antioch/eigen_utils_decl.h" #include "antioch/metaphysicl_utils_decl.h" #include "antioch/valarray_utils_decl.h" @@ -56,8 +60,78 @@ #include #include +#ifdef ANTIOCH_HAVE_GRVY +#include "grvy.h" + +GRVY::GRVY_Timer_Class gt; +#endif + + +template +int check_rate_and_derivative(const PairScalars & rate_exact, const PairScalars & derive_exact, + const PairScalars & rate, const PairScalars & derive, const PairScalars & T) +{ + typedef typename Antioch::value_type::type Scalar; + const Scalar tol = std::numeric_limits::epsilon() * 2; + + int return_flag(0); + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + if( abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (rate[2*tuple+1] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple] << std::endl + << "derive_exact = " << derive_exact[2*tuple] << std::endl + << "relative difference = " << abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple+1] << std::endl + << "derive_exact = " << derive_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + + } + return return_flag; +} + template -int vectester(const PairScalars& example) +int vectester(const PairScalars& example, const std::string & testname) { using std::abs; using std::pow; @@ -71,60 +145,70 @@ int vectester(const PairScalars& example) // Construct from example to avoid resizing issues PairScalars T = example; - T[0] = 1500.1; - T[1] = 1600.1; + PairScalars rate_exact = example; + PairScalars derive_exact = example; + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + T[2*tuple] = 1500.1L; + T[2*tuple+1] = 1600.1L; + rate_exact[2*tuple] = Cf*pow(Scalar(1500.1),eta); + rate_exact[2*tuple+1] = Cf*pow(Scalar(1600.1),eta); + derive_exact[2*tuple] = eta * Cf * pow(Scalar(1500.1),eta)/Scalar(1500.1); + derive_exact[2*tuple+1] = eta * Cf * pow(Scalar(1600.1),eta)/Scalar(1600.1); + } + Antioch::KineticsConditions cond(T); - const Scalar rate_exact0 = Cf*pow(Scalar(1500.1),eta); - const Scalar rate_exact1 = Cf*pow(Scalar(1600.1),eta); - const Scalar derive_exact0 = eta * Cf * pow(Scalar(1500.1),eta)/Scalar(1500.1); - const Scalar derive_exact1 = eta * Cf * pow(Scalar(1600.1),eta)/Scalar(1600.1); - int return_flag = 0; - const PairScalars rate = hercourtessen_rate(T);//Antioch::zero_clone(T); - const PairScalars deriveRate = hercourtessen_rate.derivative(T);//Antioch::zero_clone(T); - - const Scalar tol = std::numeric_limits::epsilon()*10; - -// hercourtessen_rate.rate_and_derivative(T,rate,deriveRate); - - if( abs( (rate[0] - rate_exact0)/rate_exact0 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T0) = " << rate[0] << std::endl - << "rate_exact = " << rate_exact0 << std::endl - << "difference = " << rate[0] - rate_exact0 << std::endl; - - return_flag = 1; - } - - if( abs( (rate[1] - rate_exact1)/rate_exact1 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T1) = " << rate[1] << std::endl - << "rate_exact = " << rate_exact1 << std::endl - << "difference = " << rate[1] - rate_exact1 << std::endl; - - return_flag = 1; - } - if( abs( (deriveRate[0] - derive_exact0)/derive_exact0 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T0) = " << deriveRate[0] << std::endl - << "derive_exact = " << derive_exact0 << std::endl; - - return_flag = 1; - } - if( abs( (deriveRate[1] - derive_exact1)/derive_exact1 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T1) = " << deriveRate[1] << std::endl - << "derive_exact = " << derive_exact1 << std::endl; - - return_flag = 1; - } +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + +// KineticsConditions + PairScalars rate = hercourtessen_rate(cond); + PairScalars derive = hercourtessen_rate.derivative(cond); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + hercourtessen_rate.rate_and_derivative(cond,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +// T +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + rate = hercourtessen_rate(T); + derive = hercourtessen_rate.derivative(T); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + hercourtessen_rate.rate_and_derivative(T,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; std::cout << "Hercourt Essen rate: " << hercourtessen_rate << std::endl; @@ -137,26 +221,42 @@ int main() int returnval = 0; returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); #ifdef ANTIOCH_HAVE_EIGEN returnval = returnval || - vectester (Eigen::Array2f()); + vectester (Eigen::Array(), "Eigen::ArrayXf"); returnval = returnval || - vectester (Eigen::Array2d()); + vectester (Eigen::Array(), "Eigen::ArrayXd"); returnval = returnval || - vectester (Eigen::Array()); + vectester (Eigen::Array(), "Eigen::ArrayXld"); #endif #ifdef ANTIOCH_HAVE_METAPHYSICL returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, float> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, long double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray"); +#endif +#ifdef ANTIOCH_HAVE_VEXCL + vex::Context ctx_f (vex::Filter::All); + if (!ctx_f.empty()) + returnval = returnval || + vectester (vex::vector (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector"); + + vex::Context ctx_d (vex::Filter::DoublePrecision); + if (!ctx_d.empty()) + returnval = returnval || + vectester (vex::vector (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector"); +#endif + +#ifdef ANTIOCH_HAVE_GRVY + gt.Finalize(); + gt.Summarize(); #endif return returnval; diff --git a/test/kooij_rate_unit.C b/test/kooij_rate_unit.C index cfce015b..4a395194 100644 --- a/test/kooij_rate_unit.C +++ b/test/kooij_rate_unit.C @@ -36,6 +36,40 @@ #include "antioch/physical_constants.h" #include "antioch/units.h" +template +int check_rate_and_derivative(const Scalar & rate_exact, const Scalar & derive_exact, + const Scalar & rate, const Scalar & derive, const Scalar & T) +{ + const Scalar tol = std::numeric_limits::epsilon() * 2; + int return_flag(0); + if( abs( (rate - rate_exact)/rate_exact ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate << std::endl + << "rate_exact = " << rate_exact << std::endl + << "relative difference = " << abs( (rate - rate_exact)/rate_exact ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive - derive_exact)/derive_exact ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive << std::endl + << "derive_exact = " << derive_exact << std::endl + << "relative difference = " << abs( (derive - derive_exact)/derive_exact ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + + return return_flag; +} + template int test_values(const Scalar & Cf, const Scalar & eta, const Scalar & Ea, const Scalar & Tref, const Scalar & R, const Antioch::KooijRate & kooij_rate) { @@ -52,56 +86,28 @@ int test_values(const Scalar & Cf, const Scalar & eta, const Scalar & Ea, const const Scalar rate_exact = Cf*pow(T/Tref,eta)*exp(-Ea/(R*T)); const Scalar derive_exact = exp(-Ea/(R*T)) * pow(T/Tref,eta) * Cf * (Ea/(R*T*T) + eta/T ); + Antioch::KineticsConditions cond(T); - Scalar rate1 = kooij_rate(T); - Scalar deriveRate1 = kooij_rate.derivative(T); - Scalar rate; - Scalar deriveRate; +// KineticsConditions + Scalar rate = kooij_rate(cond); + Scalar deriveRate = kooij_rate.derivative(cond); - kooij_rate.rate_and_derivative(T,rate,deriveRate); + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; - if( abs( (rate1 - rate_exact)/rate_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate values." << std::endl - << "T = " << T << " K" << std::endl - << "rate(T) = " << rate1 << std::endl - << "rate_exact = " << rate_exact << std::endl; + kooij_rate.rate_and_derivative(cond,rate,deriveRate); - return_flag = 1; - } - if( abs( (rate - rate_exact)/rate_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate values." << std::endl - << "T = " << T << " K" << std::endl - << "rate(T) = " << rate << std::endl - << "rate_exact = " << rate_exact << std::endl; + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; - return_flag = 1; - } - if( abs( (deriveRate1 - derive_exact)/derive_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate1 << std::endl - << "derive_exact = " << derive_exact << std::endl; - - return_flag = 1; - } - if( abs( (deriveRate - derive_exact)/derive_exact ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "T = " << T << " K" << std::endl - << "drate_dT(T) = " << deriveRate << std::endl - << "derive_exact = " << derive_exact << std::endl; +// T + rate = kooij_rate(T); + deriveRate = kooij_rate.derivative(T); + + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; + + kooij_rate.rate_and_derivative(T,rate,deriveRate); - return_flag = 1; - } + return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag; - if(return_flag)break; } return return_flag; } diff --git a/test/kooij_rate_vec_unit.C b/test/kooij_rate_vec_unit.C index 583fd56f..287289a9 100644 --- a/test/kooij_rate_vec_unit.C +++ b/test/kooij_rate_vec_unit.C @@ -43,6 +43,10 @@ #include "metaphysicl/numberarray.h" #endif +#ifdef ANTIOCH_HAVE_VEXCL +#include "vexcl/vexcl.hpp" +#endif + #include "antioch/eigen_utils_decl.h" #include "antioch/metaphysicl_utils_decl.h" #include "antioch/valarray_utils_decl.h" @@ -56,8 +60,77 @@ #include #include +#ifdef ANTIOCH_HAVE_GRVY +#include "grvy.h" + +GRVY::GRVY_Timer_Class gt; +#endif + +template +int check_rate_and_derivative(const PairScalars & rate_exact, const PairScalars & derive_exact, + const PairScalars & rate, const PairScalars & derive, const PairScalars & T) +{ + typedef typename Antioch::value_type::type Scalar; + const Scalar tol = std::numeric_limits::epsilon() * 2; + + int return_flag(0); + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + if( abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (rate[2*tuple+1] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple] << std::endl + << "derive_exact = " << derive_exact[2*tuple] << std::endl + << "relative difference = " << abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple+1] << std::endl + << "derive_exact = " << derive_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + + } + return return_flag; +} + template -int vectester(const PairScalars& example) +int vectester(const PairScalars& example, const std::string & testname) { using std::abs; using std::exp; @@ -73,60 +146,71 @@ int vectester(const PairScalars& example) // Construct from example to avoid resizing issues PairScalars T = example; - T[0] = 1500.1; - T[1] = 1600.1; + PairScalars rate_exact = example; + PairScalars derive_exact = example; - const Scalar rate_exact0 = Cf*pow(Scalar(1500.1),eta)*exp(-Ea/1500.1); - const Scalar rate_exact1 = Cf*pow(Scalar(1600.1),eta)*exp(-Ea/1600.1); - const Scalar derive_exact0 = exp(-Ea/Scalar(1500.1)) * pow(Scalar(1500.1),eta) * Cf/Scalar(1500.1) * (Ea/Scalar(1500.1) + eta ); - const Scalar derive_exact1 = exp(-Ea/Scalar(1600.1)) * pow(Scalar(1600.1),eta) * Cf/Scalar(1600.1) * (Ea/Scalar(1600.1) + eta ); + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + T[2*tuple] = 1500.1; + T[2*tuple+1] = 1600.1; + rate_exact[2*tuple] = Cf*pow(Scalar(1500.1),eta)*exp(-Ea/1500.1); + rate_exact[2*tuple+1] = Cf*pow(Scalar(1600.1),eta)*exp(-Ea/1600.1); + derive_exact[2*tuple] = exp(-Ea/Scalar(1500.1)) * pow(Scalar(1500.1),eta) * Cf/Scalar(1500.1) * (Ea/Scalar(1500.1) + eta ); + derive_exact[2*tuple+1] = exp(-Ea/Scalar(1600.1)) * pow(Scalar(1600.1),eta) * Cf/Scalar(1600.1) * (Ea/Scalar(1600.1) + eta ); + } + Antioch::KineticsConditions cond(T); int return_flag = 0; - const PairScalars rate = kooij_rate(T); - const PairScalars deriveRate = kooij_rate.derivative(T); +// KineticsConditions +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + PairScalars rate = kooij_rate(cond); + PairScalars derive = kooij_rate.derivative(cond); - const Scalar tol = std::numeric_limits::epsilon()*10; +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif -// kooij_rate.rate_and_derivative(T,rate,deriveRate); + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; - if( abs( (rate[0] - rate_exact0)/rate_exact0 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T0) = " << rate[0] << std::endl - << "rate_exact = " << rate_exact0 << std::endl - << "difference = " << rate[0] - rate_exact0 << std::endl; +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif - return_flag = 1; - } + kooij_rate.rate_and_derivative(cond,rate,derive); - if( abs( (rate[1] - rate_exact1)/rate_exact1 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T1) = " << rate[1] << std::endl - << "rate_exact = " << rate_exact1 << std::endl - << "difference = " << rate[1] - rate_exact1 << std::endl; +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif - return_flag = 1; - } - if( abs( (deriveRate[0] - derive_exact0)/derive_exact0 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T0) = " << deriveRate[0] << std::endl - << "derive_exact = " << derive_exact0 << std::endl; + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; - return_flag = 1; - } - if( abs( (deriveRate[1] - derive_exact1)/derive_exact1 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T1) = " << deriveRate[1] << std::endl - << "derive_exact = " << derive_exact1 << std::endl; +// T +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + rate = kooij_rate(T); + derive = kooij_rate.derivative(T); - return_flag = 1; - } +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + kooij_rate.rate_and_derivative(T,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; std::cout << "Kooij rate: " << kooij_rate << std::endl; @@ -139,26 +223,42 @@ int main() int returnval = 0; returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); #ifdef ANTIOCH_HAVE_EIGEN returnval = returnval || - vectester (Eigen::Array2f()); + vectester (Eigen::Array(), "Eigen::ArrayXf"); returnval = returnval || - vectester (Eigen::Array2d()); + vectester (Eigen::Array(), "Eigen::ArrayXd"); returnval = returnval || - vectester (Eigen::Array()); + vectester (Eigen::Array(), "Eigen::ArrayXld"); #endif #ifdef ANTIOCH_HAVE_METAPHYSICL returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, float> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, long double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray"); +#endif +#ifdef ANTIOCH_HAVE_VEXCL + vex::Context ctx_f (vex::Filter::All); + if (!ctx_f.empty()) + returnval = returnval || + vectester (vex::vector (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector"); + + vex::Context ctx_d (vex::Filter::DoublePrecision); + if (!ctx_d.empty()) + returnval = returnval || + vectester (vex::vector (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector"); +#endif + +#ifdef ANTIOCH_HAVE_GRVY + gt.Finalize(); + gt.Summarize(); #endif return returnval; diff --git a/test/vanthoff_rate_vec_unit.C b/test/vanthoff_rate_vec_unit.C index 363892db..cc5c8ecf 100644 --- a/test/vanthoff_rate_vec_unit.C +++ b/test/vanthoff_rate_vec_unit.C @@ -43,6 +43,10 @@ #include "metaphysicl/numberarray.h" #endif +#ifdef ANTIOCH_HAVE_VEXCL +#include "vexcl/vexcl.hpp" +#endif + #include "antioch/eigen_utils_decl.h" #include "antioch/metaphysicl_utils_decl.h" #include "antioch/valarray_utils_decl.h" @@ -56,8 +60,78 @@ #include #include +#ifdef ANTIOCH_HAVE_GRVY +#include "grvy.h" + +GRVY::GRVY_Timer_Class gt; +#endif + + +template +int check_rate_and_derivative(const PairScalars & rate_exact, const PairScalars & derive_exact, + const PairScalars & rate, const PairScalars & derive, const PairScalars & T) +{ + typedef typename Antioch::value_type::type Scalar; + const Scalar tol = std::numeric_limits::epsilon() * 2; + + int return_flag(0); + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + if( abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (rate[2*tuple+1] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate values." << std::endl + << "T = " << T << " K" << std::endl + << "rate(T) = " << rate[2*tuple] << std::endl + << "rate_exact = " << rate_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple] << std::endl + << "derive_exact = " << derive_exact[2*tuple] << std::endl + << "relative difference = " << abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + if( abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) > tol ) + { + std::cout << std::scientific << std::setprecision(16) + << "Error: Mismatch in rate derivative values." << std::endl + << "T = " << T << " K" << std::endl + << "drate_dT(T) = " << derive[2*tuple+1] << std::endl + << "derive_exact = " << derive_exact[2*tuple+1] << std::endl + << "relative difference = " << abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) << std::endl + << "tolerance = " << tol << std::endl; + + return_flag = 1; + } + + } + return return_flag; +} + template -int vectester(const PairScalars& example) +int vectester(const PairScalars& example, const std::string & testname) { using std::abs; using std::exp; @@ -74,58 +148,72 @@ int vectester(const PairScalars& example) // Construct from example to avoid resizing issues PairScalars T = example; - T[0] = 1500.1; - T[1] = 1600.1; - - const Scalar rate_exact0 = Cf*pow(Scalar(1500.1),eta)*exp(-Ea/1500.1+D*1500.1); - const Scalar rate_exact1 = Cf*pow(Scalar(1600.1),eta)*exp(-Ea/1600.1+D*1600.1); - const Scalar derive_exact0 = Cf * pow(Scalar(1500.1),eta) * exp(-Ea/Scalar(1500.1) + D*Scalar(1500.1)) * (D + eta/Scalar(1500.1) + Ea/(Scalar(1500.1)*Scalar(1500.1))); - const Scalar derive_exact1 = Cf * pow(Scalar(1600.1),eta) * exp(-Ea/Scalar(1600.1) + D*Scalar(1600.1)) * (D + eta/Scalar(1600.1) + Ea/(Scalar(1600.1)*Scalar(1600.1))); + PairScalars rate_exact = example; + PairScalars derive_exact = example; + for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple) + { + T[2*tuple] = 1500.1; + T[2*tuple+1] = 1600.1; + rate_exact[2*tuple] = Cf*pow(Scalar(1500.1),eta)*exp(-Ea/1500.1+D*1500.1); + rate_exact[2*tuple+1] = Cf*pow(Scalar(1600.1),eta)*exp(-Ea/1600.1+D*1600.1); + derive_exact[2*tuple] = Cf * pow(Scalar(1500.1),eta) * exp(-Ea/Scalar(1500.1) + D*Scalar(1500.1)) * (D + eta/Scalar(1500.1) + Ea/(Scalar(1500.1)*Scalar(1500.1))); + derive_exact[2*tuple+1] = Cf * pow(Scalar(1600.1),eta) * exp(-Ea/Scalar(1600.1) + D*Scalar(1600.1)) * (D + eta/Scalar(1600.1) + Ea/(Scalar(1600.1)*Scalar(1600.1))); + } + Antioch::KineticsConditions cond(T); int return_flag = 0; - const PairScalars rate = vanthoff_rate(T); - const PairScalars deriveRate = vanthoff_rate.derivative(T); +// KineticsConditions +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif - const Scalar tol = std::numeric_limits::epsilon()*10; + PairScalars rate = vanthoff_rate(cond); + PairScalars derive = vanthoff_rate.derivative(cond); - if( abs( (rate[0] - rate_exact0)/rate_exact0 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T0) = " << rate[0] << std::endl - << "rate_exact = " << rate_exact0 << std::endl - << "difference = " << rate[0] - rate_exact0 << std::endl; +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif - return_flag = 1; - } + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; - if( abs( (rate[1] - rate_exact1)/rate_exact1 ) > tol ) - { - std::cout << "Error: Mismatch in rate values." << std::endl - << "rate(T1) = " << rate[1] << std::endl - << "rate_exact = " << rate_exact1 << std::endl - << "difference = " << rate[1] - rate_exact1 << std::endl; +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif - return_flag = 1; - } - if( abs( (deriveRate[0] - derive_exact0)/derive_exact0 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T0) = " << deriveRate[0] << std::endl - << "derive_exact = " << derive_exact0 << std::endl; + vanthoff_rate.rate_and_derivative(cond,rate,derive); - return_flag = 1; - } - if( abs( (deriveRate[1] - derive_exact1)/derive_exact1 ) > tol ) - { - std::cout << std::scientific << std::setprecision(16) - << "Error: Mismatch in rate derivative values." << std::endl - << "drate_dT(T1) = " << deriveRate[1] << std::endl - << "derive_exact = " << derive_exact1 << std::endl; +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif - return_flag = 1; - } + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +// T +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + rate = vanthoff_rate(T); + derive = vanthoff_rate.derivative(T); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; + +#ifdef ANTIOCH_HAVE_GRVY + gt.BeginTimer(testname); +#endif + + vanthoff_rate.rate_and_derivative(cond,rate,derive); + +#ifdef ANTIOCH_HAVE_GRVY + gt.EndTimer(testname); +#endif + + return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag; std::cout << "Van't Hoff rate: " << vanthoff_rate << std::endl; @@ -138,26 +226,42 @@ int main() int returnval = 0; returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); returnval = returnval || - vectester (std::valarray(2)); + vectester (std::valarray(2*ANTIOCH_N_TUPLES), "valarray"); #ifdef ANTIOCH_HAVE_EIGEN returnval = returnval || - vectester (Eigen::Array2f()); + vectester (Eigen::Array(), "Eigen::ArrayXf"); returnval = returnval || - vectester (Eigen::Array2d()); + vectester (Eigen::Array(), "Eigen::ArrayXd"); returnval = returnval || - vectester (Eigen::Array()); + vectester (Eigen::Array(), "Eigen::ArrayXld"); #endif #ifdef ANTIOCH_HAVE_METAPHYSICL returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, float> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray"); returnval = returnval || - vectester (MetaPhysicL::NumberArray<2, long double> (0)); + vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray"); +#endif +#ifdef ANTIOCH_HAVE_VEXCL + vex::Context ctx_f (vex::Filter::All); + if (!ctx_f.empty()) + returnval = returnval || + vectester (vex::vector (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector"); + + vex::Context ctx_d (vex::Filter::DoublePrecision); + if (!ctx_d.empty()) + returnval = returnval || + vectester (vex::vector (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector"); +#endif + +#ifdef ANTIOCH_HAVE_GRVY + gt.Finalize(); + gt.Summarize(); #endif return returnval; From 484731e0762b85b11b583fb46cd608fb85a0997e Mon Sep 17 00:00:00 2001 From: Sylvain Plessis Date: Fri, 1 May 2015 10:54:34 -0500 Subject: [PATCH 3/4] unused variable removed --- test/kooij_rate_unit.C | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/kooij_rate_unit.C b/test/kooij_rate_unit.C index 4a395194..758e2a79 100644 --- a/test/kooij_rate_unit.C +++ b/test/kooij_rate_unit.C @@ -79,8 +79,6 @@ int test_values(const Scalar & Cf, const Scalar & eta, const Scalar & Ea, const int return_flag = 0; - const Scalar tol = std::numeric_limits::epsilon() * 100; - for(Scalar T = 300.1L; T <= 2500.1L; T += 10.L) { From 679546f0ce62ca39c41ad531df5a93f871b74816 Mon Sep 17 00:00:00 2001 From: Sylvain Plessis Date: Fri, 1 May 2015 22:29:05 -0500 Subject: [PATCH 4/4] fork function of base object now passing KineticsConditions to all kinetics models --- src/kinetics/include/antioch/kinetics_type.h | 24 ++++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/kinetics/include/antioch/kinetics_type.h b/src/kinetics/include/antioch/kinetics_type.h index 2db8b107..2d062481 100644 --- a/src/kinetics/include/antioch/kinetics_type.h +++ b/src/kinetics/include/antioch/kinetics_type.h @@ -184,25 +184,25 @@ namespace Antioch{ { case(KineticsModel::CONSTANT): { - return (static_cast*>(this))->rate(conditions.T()); + return (static_cast*>(this))->rate(conditions); } break; case(KineticsModel::HERCOURT_ESSEN): { - return (static_cast*>(this))->rate(conditions.T()); + return (static_cast*>(this))->rate(conditions); } break; case(KineticsModel::BERTHELOT): { - return (static_cast*>(this))->rate(conditions.T()); + return (static_cast*>(this))->rate(conditions); } break; case(KineticsModel::ARRHENIUS): { - return (static_cast*>(this))->rate(conditions.T()); + return (static_cast*>(this))->rate(conditions); } break; @@ -261,25 +261,25 @@ namespace Antioch{ { case(KineticsModel::CONSTANT): { - return (static_cast*>(this))->derivative(conditions.T()); + return (static_cast*>(this))->derivative(conditions); } break; case(KineticsModel::HERCOURT_ESSEN): { - return (static_cast*>(this))->derivative(conditions.T()); + return (static_cast*>(this))->derivative(conditions); } break; case(KineticsModel::BERTHELOT): { - return (static_cast*>(this))->derivative(conditions.T()); + return (static_cast*>(this))->derivative(conditions); } break; case(KineticsModel::ARRHENIUS): { - return (static_cast*>(this))->derivative(conditions.T()); + return (static_cast*>(this))->derivative(conditions); } break; @@ -339,25 +339,25 @@ namespace Antioch{ { case(KineticsModel::CONSTANT): { - (static_cast*>(this))->rate_and_derivative(conditions.T(),rate,drate_dT); + (static_cast*>(this))->rate_and_derivative(conditions,rate,drate_dT); } break; case(KineticsModel::HERCOURT_ESSEN): { - (static_cast*>(this))->rate_and_derivative(conditions.T(),rate,drate_dT); + (static_cast*>(this))->rate_and_derivative(conditions,rate,drate_dT); } break; case(KineticsModel::BERTHELOT): { - (static_cast*>(this))->rate_and_derivative(conditions.T(),rate,drate_dT); + (static_cast*>(this))->rate_and_derivative(conditions,rate,drate_dT); } break; case(KineticsModel::ARRHENIUS): { - (static_cast*>(this))->rate_and_derivative(conditions.T(),rate,drate_dT); + (static_cast*>(this))->rate_and_derivative(conditions,rate,drate_dT); } break;