rotational_relaxation_vec_unit.C File Reference

#include "antioch_config.h"
#include <valarray>
#include "Eigen/Dense"
#include "metaphysicl/numberarray.h"
#include "vexcl/vexcl.hpp"
#include "antioch/eigen_utils_decl.h"
#include "antioch/metaphysicl_utils_decl.h"
#include "antioch/valarray_utils_decl.h"
#include "antioch/vexcl_utils_decl.h"
#include "antioch/rotational_relaxation.h"
#include "antioch/eigen_utils.h"
#include "antioch/metaphysicl_utils.h"
#include "antioch/valarray_utils.h"
#include "antioch/vexcl_utils.h"
#include "antioch/cmath_shims.h"
#include <cmath>
#include <limits>

Go to the source code of this file.

Functions
const long double	pi (3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825)
template<typename Scalar >
Scalar	F (const Scalar &x)
template<typename Scalar >
Scalar	z (const Scalar &T, const Scalar &eps_kb, const Scalar &z_298)
template<typename PairScalars >
int	vectester (const PairScalars &example, const std::string &testname)
int	main ()

Function Documentation

template<typename Scalar >

Scalar F

(

const Scalar &

x

)

Definition at line 72 of file rotational_relaxation_vec_unit.C.

References pi().

Referenced by z().

{
  return 1.L + Antioch::ant_pow(pi,(float)1.5)/2.L  * Antioch::ant_sqrt(x)
             + (pi*pi/4.L + 2.L) * x
             + Antioch::ant_pow(pi * x,(float)1.5);
}

int main

(

)

Definition at line 160 of file rotational_relaxation_vec_unit.C.

References vectester().

{
  int returnval = 0;

  returnval = returnval ||
    vectester (std::valarray<float>(2*ANTIOCH_N_TUPLES), "valarray<float>");
  returnval = returnval ||
    vectester (std::valarray<double>(2*ANTIOCH_N_TUPLES), "valarray<double>");
  returnval = returnval ||
    vectester (std::valarray<long double>(2*ANTIOCH_N_TUPLES), "valarray<ld>");
#ifdef ANTIOCH_HAVE_EIGEN
  returnval = returnval ||
    vectester (Eigen::Array<float, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXf");
  returnval = returnval ||
    vectester (Eigen::Array<double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXd");
  returnval = returnval ||
    vectester (Eigen::Array<long double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXld");
#endif
#ifdef ANTIOCH_HAVE_METAPHYSICL
  returnval = returnval ||
    vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray<float>");
  returnval = returnval ||
    vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray<double>");
  returnval = returnval ||
    vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray<ld>");
#endif
#ifdef ANTIOCH_HAVE_VEXCL
  vex::Context ctx_f (vex::Filter::All);
  if (!ctx_f.empty())
    returnval = returnval ||
      vectester (vex::vector<float> (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector<float>");

  vex::Context ctx_d (vex::Filter::DoublePrecision);
  if (!ctx_d.empty())
    returnval = returnval ||
      vectester (vex::vector<double> (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector<double>");
#endif

#ifdef ANTIOCH_HAVE_GRVY
  gt.Finalize();
  gt.Summarize();
#endif

  return returnval;
}

const long double pi

(

3.

141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825

)

Referenced by F().

template<typename PairScalars >

int vectester	(	const PairScalars &	example,
		const std::string &	testname
	)

Definition at line 87 of file rotational_relaxation_vec_unit.C.

References Z(), and z().

Referenced by main().

{

std::cout << "entering " << testname << std::endl;
  typedef typename Antioch::value_type<PairScalars>::type Scalar;

  const Scalar eps_kb = 97.53L; // N2 value
  const Scalar z_298 = 4.0L; // N2 value

  Antioch::RotationalRelaxation<Scalar> rot(z_298,eps_kb);

  // Construct from example to avoid resizing issues
  PairScalars T = example;
  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
    {
      T[2*tuple]   = 1500.1L;
      T[2*tuple+1] = 1600.1L;
    }
  
  int return_flag = 0;

#ifdef ANTIOCH_HAVE_GRVY
  gt.BeginTimer(testname);
#endif

  const PairScalars Z = rot(T);

#ifdef ANTIOCH_HAVE_GRVY
  gt.EndTimer(testname);
#endif

  const Scalar tol = (std::numeric_limits<Scalar>::epsilon()*10 < 5e-17)?
                        5e-17:
                        std::numeric_limits<Scalar>::epsilon()*10;

  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
    {
      const Scalar Z_exact0 = z( (Scalar)1500.1, eps_kb, z_298);
      const Scalar Z_exact1 = z( (Scalar)1600.1, eps_kb, z_298);

      if( abs( (Z[2*tuple] - Z_exact0)/Z_exact0 ) > tol )
        {
          std::cout << std::scientific << std::setprecision(20)
                    << "Error: Mismatch in Z values in test " << testname << std::endl
                    << "Z(T0)   = "             << Z[2*tuple]                                 << std::endl
                    << "Z_exact = "             << Z_exact0                                   << std::endl
                    << "absolute difference = " << Z[2*tuple] - Z_exact0                      << std::endl
                    << "relative difference = " << std::abs(Z[2*tuple] - Z_exact0) / Z_exact0 << std::endl
                    << "tolerance = "           << tol << std::endl;

          return_flag = 1;
          break;
        }

      if( abs( (Z[2*tuple+1] - Z_exact1)/Z_exact1 ) > tol )
        {
          std::cout << std::scientific << std::setprecision(20)
                    << "Error: Mismatch in Z values in test " << testname << std::endl
                    << "Z(T1)   = "             << Z[2*tuple+1]                                 << std::endl
                    << "Z_exact = "             << Z_exact1                                     << std::endl
                    << "absolute difference = " << Z[2*tuple+1] - Z_exact1                      << std::endl
                    << "relative difference = " << std::abs(Z[2*tuple+1] - Z_exact1) / Z_exact1 << std::endl
                    << "tolerance = "           << tol << std::endl;

          return_flag = 1;
          break;
        }
    }

  return return_flag;
}

template<typename Scalar >

Scalar z	(	const Scalar &	T,
		const Scalar &	eps_kb,
		const Scalar &	z_298
	)

Definition at line 80 of file rotational_relaxation_vec_unit.C.

References F().

Referenced by tester(), and vectester().

{
  return z_298 * F(eps_kb / 298.L) / F(eps_kb / T);
}