From 7a6b6158430d911d52efdb0c99602dbc27e624af Mon Sep 17 00:00:00 2001
From: Nicholas Malaya <nicholas.malaya@gmail.com>
Date: Fri, 17 Aug 2012 21:25:11 +0000
Subject: [PATCH] adding statmech transport test

---
 test_problems/statmech/Makefile.am            |   7 +-
 test_problems/statmech/statmech_transport.cpp | 113 ++++++++++++++++++
 2 files changed, 118 insertions(+), 2 deletions(-)
 create mode 100644 test_problems/statmech/statmech_transport.cpp

diff --git a/test_problems/statmech/Makefile.am b/test_problems/statmech/Makefile.am
index e72774e5b..ffaf5fa28 100644
--- a/test_problems/statmech/Makefile.am
+++ b/test_problems/statmech/Makefile.am
@@ -10,7 +10,8 @@ LIBS                 =  $(LINK)
 
 check_PROGRAMS         = statmech_test \
                          statmech_test_poly \
-                         statmech_test_Fe
+                         statmech_test_Fe   \
+	                 statmech_transport
 
 library_includedir   = $(INC)
 
@@ -21,10 +22,12 @@ library_includedir   = $(INC)
 statmech_test_SOURCES      = statmech_test.cpp
 statmech_test_poly_SOURCES = statmech_test_poly.cpp
 statmech_test_Fe_SOURCES   = statmech_test_Fe.cpp
+statmech_transport_SOURCES = statmech_transport.cpp
 
 TESTS_ENVIRONMENT = 	
 TESTS             = statmech_test \
-                    statmech_test_Fe
+                    statmech_test_Fe \
+	            statmech_transport
 #                    statmech_test_poly
 
 EXTRA_DIST = runtest_stat 
\ No newline at end of file
diff --git a/test_problems/statmech/statmech_transport.cpp b/test_problems/statmech/statmech_transport.cpp
new file mode 100644
index 000000000..ca9a4ae8f
--- /dev/null
+++ b/test_problems/statmech/statmech_transport.cpp
@@ -0,0 +1,113 @@
+/**
+ *  @file statmech
+ *       test problem for statistical mechanics in cantera
+ */
+
+//  Example 
+//
+// Test case for the statistical mechanics in cantera
+//
+
+#include <iostream>
+#include <string>
+#include <vector>
+#include <string>
+#include <iomanip>
+
+/*****************************************************************/
+/*****************************************************************/
+
+//#include "cantera/Cantera.h"
+#include "cantera/transport.h"
+#include "cantera/IdealGasMix.h"
+#include "cantera/equil/equil.h"
+
+#include "cantera/transport/TransportFactory.h"
+
+using namespace std;
+using namespace Cantera;
+
+int main(int argc, char** argv) 
+{
+
+  try
+    {
+      int k;
+      IdealGasMix g("test_stat.xml");
+      int nsp = g.nSpecies();
+      double pres = 1.0E5;
+
+      vector_fp Xset(nsp, 0.0);
+      Xset[0] =  0.5 ;
+      Xset[1] =  0.5;
+  
+      g.setState_TPX(1500.0, pres, DATA_PTR(Xset));
+      equilibrate(g, "TP", -1);
+
+      vector_fp cp_R(nsp, 0.0);
+      g.getCp_R(DATA_PTR(cp_R));
+
+      //for(int i=0;i<nsp;i++)
+      //{
+      //  std::cout.precision(10);
+      //  std::cout << cp_R[i] << std::endl;
+      //	}  
+
+      // error check-- exactly 2.5 for atoms
+      if(cp_R[0] != 2.5)
+	{
+	  std::cout << "Error for monotomic Species!\n";
+	  return 1;
+	}
+
+      // error check: analytical result is more complicated for 
+      // molecules. One species should suffice, lets try NO2, with
+      // three vibrational modes:
+      /// theta[0]: 1.07900e3
+      /// theta[1]: 1.90000e3
+      /// theta[2]: 2.32700e3
+      // at T = 1500
+      //
+      // This is precisely: 6.655804161 (e.g. 5/2 + 2 + 3.1558..)
+      //
+      double theta[3];
+      theta[0] = 1.07900e3;
+      theta[1] = 1.90000e3;
+      theta[2] = 2.32700e3;
+
+      double T;
+      T = 1500.0;
+
+      double denom;
+      double ctr = 0.0;
+      double GasConstant = 1.0;
+
+      for(int i = 0; i < 3; i++)
+	{
+	  denom = exp(2*theta[i]/T) - 2* exp(theta[i]/T) + 1;
+	  ctr += GasConstant * theta[i] * (theta[i] * exp(theta[i]/T)/(T*T))/ (denom);
+	  //std::cout << "survey says: " << ctr << " and denom is: " << denom << std::endl;
+	}
+      //std::cout << "survey says: " << ctr << " and denom is: " << denom << std::endl;
+      double sol = ctr + 5/2 + 2;
+      double tol = 1e-9;
+
+      if(abs(cp_R[3] - sol) >= tol )
+	{
+	  double diff = cp_R[3]-sol;
+	  std::cout << "Error for Species NO2!\n";
+	  std::cout << "Diff was: " << diff << "\n";
+	  return 1;
+	}
+
+    }
+    catch (CanteraError) 
+      {
+	showErrors(cout);
+	return 1;
+      }
+
+  // Mark it zero!
+  return 0;
+
+}