From 2d2c16ec3869bff519b06761cb39d5d77886896c Mon Sep 17 00:00:00 2001
From: Nicholas Malaya <nicholas.malaya@gmail.com>
Date: Tue, 17 Jul 2012 21:22:21 +0000
Subject: [PATCH] adding statmech testing suite

---
 configure.ac                                  |   1 +
 test_problems/Makefile.am                     |   2 +-
 test_problems/statmech/Makefile.am            |  30 +++++
 test_problems/statmech/runtest_stat           |  15 +++
 .../statmech/statmech_properties.cpp          |  64 ++++++++++
 test_problems/statmech/statmech_test.cpp      | 115 ++++++++++++++++++
 test_problems/statmech/statmech_test_Fe.cpp   |  75 ++++++++++++
 test_problems/statmech/statmech_test_poly.cpp |  77 ++++++++++++
 test_problems/statmech/test.xml               |  37 ++++++
 test_problems/statmech/test_stat.xml          |  66 ++++++++++
 test_problems/statmech/test_stat_Fe.xml       |  44 +++++++
 test_problems/statmech/test_stat_err.xml      |  47 +++++++
 12 files changed, 572 insertions(+), 1 deletion(-)
 create mode 100644 test_problems/statmech/Makefile.am
 create mode 100755 test_problems/statmech/runtest_stat
 create mode 100644 test_problems/statmech/statmech_properties.cpp
 create mode 100644 test_problems/statmech/statmech_test.cpp
 create mode 100644 test_problems/statmech/statmech_test_Fe.cpp
 create mode 100644 test_problems/statmech/statmech_test_poly.cpp
 create mode 100644 test_problems/statmech/test.xml
 create mode 100644 test_problems/statmech/test_stat.xml
 create mode 100644 test_problems/statmech/test_stat_Fe.xml
 create mode 100644 test_problems/statmech/test_stat_err.xml

diff --git a/configure.ac b/configure.ac
index e0f0a2ab4..e198538c9 100644
--- a/configure.ac
+++ b/configure.ac
@@ -151,6 +151,7 @@ AC_OUTPUT(Makefile                                                         \
                    test_problems/cathermo/VPissp/Makefile                  \
                    test_problems/cathermo/wtWater/Makefile                 \
                    test_problems/PecosTransport/Makefile                   \
+                   test_problems/statmech/Makefile                         \
 		   doc/Makefile                                            \
 		   doc/doxygen/Makefile                                    \
 		   src/Makefile                                            \
diff --git a/test_problems/Makefile.am b/test_problems/Makefile.am
index 62144887d..17e52406f 100644
--- a/test_problems/Makefile.am
+++ b/test_problems/Makefile.am
@@ -4,7 +4,7 @@ SUBDIRS += ChemEquil_ionizedGas ChemEquil_red1 CpJump cxx_ex diamondSurf
 SUBDIRS += diamondSurf_dupl fracCoeff multiGasTransport NASA9poly_test
 SUBDIRS += negATest printUtilUnitTest pureFluidTest silane_equil 
 SUBDIRS += spectroscopy surfkin surfSolverTest VPsilane_test 
-SUBDIRS += VCSnonideal PecosTransport 
+SUBDIRS += VCSnonideal PecosTransport statmech
 
 #f77test
 
diff --git a/test_problems/statmech/Makefile.am b/test_problems/statmech/Makefile.am
new file mode 100644
index 000000000..a42826f01
--- /dev/null
+++ b/test_problems/statmech/Makefile.am
@@ -0,0 +1,30 @@
+INC                  = -I. -I$(top_builddir)/build/include/ -g
+AM_CPPFLAGS          = $(INC)
+AM_CXXFLAGS          = $(AM_CPPFLAGS)
+
+LINK                 = -lctcxx -luser -loneD -lzeroD -lequil -lkinetics -ltransport -lthermo 
+LINK                += -lctnumerics -lctmath -ltpx -lctspectra -lconverters -lctbase -lcvode 
+LINK                += -lctlapack -lctblas -lctf2c -lm -lstdc++
+AM_LDFLAGS           = -L$(top_builddir)/build/lib/ 
+LIBS                 =  $(LINK)
+
+check_PROGRAMS         = statmech_test \
+                         statmech_test_poly \
+                         statmech_test_Fe
+
+library_includedir   = $(INC)
+
+#-----------------------
+# Cantera DH graph test
+#-----------------------
+
+statmech_test_SOURCES      = statmech_test.cpp
+statmech_test_poly_SOURCES = statmech_test_poly.cpp
+statmech_test_Fe_SOURCES   = statmech_test_Fe.cpp
+
+TESTS_ENVIRONMENT = 	
+TESTS             = statmech_test \
+                    statmech_test_poly \
+                    statmech_test_Fe
+
+EXTRA_DIST = runtest_stat 
\ No newline at end of file
diff --git a/test_problems/statmech/runtest_stat b/test_problems/statmech/runtest_stat
new file mode 100755
index 000000000..48da88d0c
--- /dev/null
+++ b/test_problems/statmech/runtest_stat
@@ -0,0 +1,15 @@
+#!/bin/sh
+#
+#
+temp_success="1"
+/bin/rm  -f output.txt outputa.txt 
+tname="mixGasTransport"
+#################################################################
+#
+#################################################################
+CANTERA_DATA=${CANTERA_DATA:=../../data/inputs}; export CANTERA_DATA
+
+CANTERA_BIN=${CANTERA_BIN:=../../bin}
+./statmech_test > output.txt
+
+exit $?
diff --git a/test_problems/statmech/statmech_properties.cpp b/test_problems/statmech/statmech_properties.cpp
new file mode 100644
index 000000000..38fa68052
--- /dev/null
+++ b/test_problems/statmech/statmech_properties.cpp
@@ -0,0 +1,64 @@
+/**
+ *  @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>
+
+using namespace std;
+
+/*****************************************************************/
+/*****************************************************************/
+
+#include "Cantera.h"
+#include "transport.h"
+#include "IdealGasMix.h"
+#include "equil.h"
+
+#include "TransportFactory.h"
+
+using namespace Cantera;
+using namespace Cantera_CXX;
+
+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));
+
+      
+
+    }
+    catch (CanteraError) 
+      {
+	showErrors(cout);
+	return 1;
+      }
+
+  // Mark it zero!
+  return 0;
+
+}
diff --git a/test_problems/statmech/statmech_test.cpp b/test_problems/statmech/statmech_test.cpp
new file mode 100644
index 000000000..ab3f2a68d
--- /dev/null
+++ b/test_problems/statmech/statmech_test.cpp
@@ -0,0 +1,115 @@
+/**
+ *  @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>
+
+using namespace std;
+
+/*****************************************************************/
+/*****************************************************************/
+
+#include "Cantera.h"
+#include "transport.h"
+#include "IdealGasMix.h"
+#include "equil.h"
+
+#include "TransportFactory.h"
+
+using namespace Cantera;
+using namespace Cantera_CXX;
+
+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;
+
+}
diff --git a/test_problems/statmech/statmech_test_Fe.cpp b/test_problems/statmech/statmech_test_Fe.cpp
new file mode 100644
index 000000000..86e76b733
--- /dev/null
+++ b/test_problems/statmech/statmech_test_Fe.cpp
@@ -0,0 +1,75 @@
+/**
+ *  @file statmech
+ *       test problem for statistical mechanics in cantera
+ */
+
+//  Example 
+//
+// Test case to check error thrown if using Fe (not supported species)
+//
+
+#include <iostream>
+#include <string>
+#include <vector>
+#include <string>
+#include <iomanip>
+
+using namespace std;
+
+/*****************************************************************/
+/*****************************************************************/
+
+#include "Cantera.h"
+#include "transport.h"
+#include "IdealGasMix.h"
+#include "equil.h"
+
+#include "TransportFactory.h"
+
+using namespace Cantera;
+using namespace Cantera_CXX;
+
+int main(int argc, char** argv) 
+{
+
+  try
+    {
+      int k;
+      IdealGasMix g("test_stat_Fe.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 << cp_R[i] << std::endl;
+	}  
+
+      // error check
+      if(cp_R[0] != 0)
+	{
+	  std::cout << "Error for monotomic Species!\n";
+	  return 1;
+	}
+    }
+    catch (CanteraError) 
+      {
+	// need to get error here because of loading Fe in input file, when
+	// no Fe exists in the species information table, in statmech.cpp	
+	//showErrors(cout);
+	return 0;
+      }
+
+  // Mark it zero!
+  return 1;
+
+}
diff --git a/test_problems/statmech/statmech_test_poly.cpp b/test_problems/statmech/statmech_test_poly.cpp
new file mode 100644
index 000000000..b13721b74
--- /dev/null
+++ b/test_problems/statmech/statmech_test_poly.cpp
@@ -0,0 +1,77 @@
+/**
+ *  @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>
+
+using namespace std;
+
+/*****************************************************************/
+/*****************************************************************/
+
+#include "Cantera.h"
+#include "transport.h"
+#include "IdealGasMix.h"
+#include "equil.h"
+
+#include "TransportFactory.h"
+
+using namespace Cantera;
+using namespace Cantera_CXX;
+
+int main(int argc, char** argv) 
+{
+
+  try
+    {
+      int k;
+      IdealGasMix g("test_stat_err.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 << cp_R[i] << std::endl;
+	}  
+
+      // error check
+      if(cp_R[0] != 0)
+	{
+	  std::cout << "Error for monotomic Species!\n";
+	  return 1;
+	}
+    }
+    catch (CanteraError) 
+      {
+	// we wanted to catch an error here for incorrectly trying to use poly methods
+	// for the statmech species data, so exit successfully
+	//showErrors(cout);
+
+	// Mark it zero!
+	return 0;
+      }
+
+  // something is wrong: we were suppose to catch an error here: paradox!
+  return 1;
+
+}
diff --git a/test_problems/statmech/test.xml b/test_problems/statmech/test.xml
new file mode 100644
index 000000000..5eb783af8
--- /dev/null
+++ b/test_problems/statmech/test.xml
@@ -0,0 +1,37 @@
+
+<?xml version="1.0"?>
+<ctml>
+  <validate reactions="yes" species="yes"/>
+                   
+ <!-- phase Na    -->
+  <phase dim="3" id="H">
+    <elementArray datasrc="elements.xml">
+        O H C Fe Ca N Na Cl
+    </elementArray>
+    <speciesArray datasrc="#species_test"> H </speciesArray>
+    <thermo model="IdealGas">
+      <density units="g/cm3">2.165</density>
+    </thermo>
+    <transport model="None"/>
+    <kinetics model="none"/>
+  </phase>
+
+ <!-- species definitions     -->
+  <speciesData id="species_test">
+
+    <!-- species H   -->
+    <species name="H">
+      <atomArray> H:1 </atomArray>
+      <thermo>
+        <NASA P0="100000.0" Tmax="1000.0" Tmin="10.0">
+          <floatArray name="coeffs" size="7">
+            2.344331120E+00,   7.980520750E-03,  -1.947815100E-05,   2.015720940E-08, 
+           -7.376117610E-12,  -9.179351730E+02,   6.830102380E-01</floatArray>
+        </NASA>
+      </thermo>
+      <density units="g/cm3">2.165</density>
+    </species>
+
+  </speciesData>
+
+</ctml>
diff --git a/test_problems/statmech/test_stat.xml b/test_problems/statmech/test_stat.xml
new file mode 100644
index 000000000..2938c4153
--- /dev/null
+++ b/test_problems/statmech/test_stat.xml
@@ -0,0 +1,66 @@
+
+<?xml version="1.0"?>
+<ctml>
+  <validate reactions="yes" species="yes"/>
+                   
+ <!-- phase H    -->
+  <phase dim="3" id="H">
+    <elementArray datasrc="elements.xml">
+        O H C N Na
+    </elementArray>
+    <speciesArray datasrc="#species_test"> H O N NO2</speciesArray>
+    <thermo model="IdealGas">
+      <density units="g/cm3">2.165</density>
+    </thermo>
+    <transport model="None"/>
+    <kinetics model="none"/>
+  </phase>
+
+ <!-- species definitions     -->
+  <speciesData id="species_test">
+
+    <!-- species H   -->
+    <species name="H">
+      <atomArray> H:1 </atomArray>
+      <thermo>
+        <StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
+        </StatMech>
+      </thermo>
+      <density units="g/cm3">2.165</density>
+    </species>
+
+    <!-- species O   -->
+    <species name="O">
+      <atomArray>O:1</atomArray>
+      <thermo>
+        <StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
+        </StatMech>
+      </thermo>
+      <density units="g/cm3">2.165</density>
+    </species>
+
+    <!-- species N   -->
+    <species name="N">
+      <atomArray>N:1</atomArray>
+      <thermo>
+        <StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
+        </StatMech>
+      </thermo>
+      <density units="g/cm3">2.165</density>
+    </species>
+
+    <!-- species NO2   -->
+    <species name="NO2">
+      <atomArray>O:2 N:1</atomArray>
+      <thermo>
+        <StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
+        </StatMech>
+        <StatMech P0="100000.0" Tmax="3000.0" Tmin="1.0">
+        </StatMech>
+      </thermo>
+      <density units="g/cm3">2.165</density>
+    </species>
+
+  </speciesData>
+
+</ctml>
diff --git a/test_problems/statmech/test_stat_Fe.xml b/test_problems/statmech/test_stat_Fe.xml
new file mode 100644
index 000000000..f1573815f
--- /dev/null
+++ b/test_problems/statmech/test_stat_Fe.xml
@@ -0,0 +1,44 @@
+
+<?xml version="1.0"?>
+<ctml>
+  <validate reactions="yes" species="yes"/>
+                   
+ <!-- phase H    -->
+  <phase dim="3" id="H">
+    <elementArray datasrc="elements.xml">
+        O H C N Na Cl Fe
+    </elementArray>
+    <speciesArray datasrc="#species_test"> H Fe</speciesArray>
+    <thermo model="IdealGas">
+      <density units="g/cm3">2.165</density>
+    </thermo>
+    <transport model="None"/>
+    <kinetics model="none"/>
+  </phase>
+
+ <!-- species definitions     -->
+  <speciesData id="species_test">
+
+    <!-- species H   -->
+    <species name="H">
+      <atomArray> H:1 </atomArray>
+      <thermo>
+        <StatMech P0="100000.0" Tmax="2000.0" Tmin="1.0">
+        </StatMech>
+      </thermo>
+      <density units="g/cm3">2.165</density>
+    </species>
+
+    <!-- species Fe   -->
+    <species name="Fe">
+      <atomArray> Fe:1 </atomArray>
+      <thermo>
+        <StatMech P0="100000.0" Tmax="2000.0" Tmin="1.0">
+        </StatMech>
+      </thermo>
+      <density units="g/cm3">2.165</density>
+    </species>
+
+  </speciesData>
+
+</ctml>
diff --git a/test_problems/statmech/test_stat_err.xml b/test_problems/statmech/test_stat_err.xml
new file mode 100644
index 000000000..ae8634a6b
--- /dev/null
+++ b/test_problems/statmech/test_stat_err.xml
@@ -0,0 +1,47 @@
+
+<?xml version="1.0"?>
+<ctml>
+  <validate reactions="yes" species="yes"/>
+                   
+ <!-- phase H    -->
+  <phase dim="3" id="H">
+    <elementArray datasrc="elements.xml">
+        O H C N Na Cl
+    </elementArray>
+    <speciesArray datasrc="#species_test"> H O</speciesArray>
+    <thermo model="IdealGas">
+      <density units="g/cm3">2.165</density>
+    </thermo>
+    <transport model="None"/>
+    <kinetics model="none"/>
+  </phase>
+
+ <!-- species definitions     -->
+  <speciesData id="species_test">
+
+    <!-- species H   -->
+    <species name="H">
+      <atomArray> H:1 </atomArray>
+      <thermo>
+        <StatMech P0="100000.0" Tmax="1000.0" Tmin="0.1">
+        </StatMech>
+      </thermo>
+      <density units="g/cm3">2.165</density>
+    </species>
+
+    <!-- species O   -->
+    <species name="O">
+      <atomArray> O:1 </atomArray>
+      <thermo>
+        <NASA9 P0="100000.0" Tmax="1000.0" Tmin="0.1">
+          <floatArray name="coeffs" size="7">
+            2.344331120E+00,   7.980520750E-03,  -1.947815100E-05,   2.015720940E-08, 
+           -7.376117610E-12,  -9.179351730E+02,   6.830102380E-01</floatArray>
+        </NASA9>
+      </thermo>
+      <density units="g/cm3">2.165</density>
+    </species>
+
+  </speciesData>
+
+</ctml>