diff --git a/test_problems/cathermo/HMW_graph_CpvT/.cvsignore b/test_problems/cathermo/HMW_graph_CpvT/.cvsignore new file mode 100644 index 000000000..c26dfc44d --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/.cvsignore @@ -0,0 +1,11 @@ +Makefile +diff_test.out +output.txt +outputa.txt +sortAlgorithms.d +table.csv +.depends +Cp_standalone +HMW_graph_CpvT +HMW_graph_CpvT.d + diff --git a/test_problems/cathermo/HMW_graph_CpvT/Cp_standalone.cpp b/test_problems/cathermo/HMW_graph_CpvT/Cp_standalone.cpp new file mode 100644 index 000000000..2b11c0a61 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/Cp_standalone.cpp @@ -0,0 +1,192 @@ + +#include +#include +#include + +using namespace std; + + +/* + * Values of A_J/R : tabular form + * units sqrt(kg/gmol), + */ +double A_JdR(double temp) { + double retn; + if (temp == 323.15) { + retn = 5.50274; + } else if (temp == 473.15) { + retn = 24.8263; + } else { + printf("A_JdR unknown temp value %g\n", temp); + exit(-1); + } + return retn; +} + +double Beta0(double temp, int ifunc) { + double q1 = 0.0765; + double q2 = -777.03; + double q3 = -4.4706; + double q4 = 0.008946; + double q5 = -3.3158E-6; + double retn; + double tref = 298.15; + if (ifunc == 0) { + retn = q1 + q2 * (1.0/temp - 1.0/tref) + + q3 * (log(temp/tref)) + q4 * (temp - tref) + + q5 * (temp * temp - tref * tref); + } else if (ifunc == 1) { + retn = (- q2 * 1.0/(temp* temp) + + q3 / temp + + q4 + + 2.0 * temp * q5); + } else if (ifunc == 2) { + retn = ( 2.0 * q2 * 1.0/(temp* temp*temp) + - q3 / (temp*temp) + + 2.0 * q5); + } else { + exit(-1); + } + return retn; +} + +double Beta1(double temp, int ifunc) { + double q6 = 0.2664; + double q9 = 6.1608E-5; + double q10 = 1.0715E-6; + double retn; + double tref = 298.15; + if (ifunc == 0) { + retn = q6 + q9 * (temp - tref) + + q10 * (temp * temp - tref * tref); + } else if (ifunc == 1) { + retn = q9 + 2.0 * q10 * temp; + } else if (ifunc == 2) { + retn = 2.0 * q10; + } else { + exit(-1); + } + return retn; +} + +double Cphi(double temp, int ifunc) { + double q11 = 0.00127; + double q12 = 33.317; + double q13 = 0.09421; + double q14 = -4.655E-5; + double retn; + double tref = 298.15; + if (ifunc == 0) { + retn = q11 + q12 * (1.0/temp - 1.0/tref) + + q13 * (log(temp/tref)) + q14 * (temp - tref); + } else if (ifunc == 1) { + retn = - q12 / (temp * temp) + + q13 / temp + q14; + } else if (ifunc == 2) { + retn = + 2.0 * q12 / (temp * temp * temp) + - q13 / (temp * temp) ; + } else { + exit(-1); + } + return retn; +} + +double calc(double temp, double Iionic) { + /* + * Gas Constant in J gmol-1 K-1 + */ + double GasConst = 8.314472; + + double Aphi = 0.0; + if (temp == 323.15) { + Aphi = 0.4102995331359; + } else if (temp == 473.15) { + Aphi = 0.622777; + } else { + printf("ERROR: unknown temp\n"); + exit(-1); + } + //printf(" Aphi = %g\n", Aphi); + + /* + * Calculate A_H in J gmol-1 sqrt(kg/gmol) + */ + double A_J = A_JdR(temp); + A_J *= GasConst; + + double beta0prime2 = Beta0(temp, 2); + printf(" beta0prime2 = %g\n", beta0prime2); + + double beta1prime2 = Beta1(temp, 2); + printf(" beta1prime2 = %g\n", beta1prime2); + + double cphiprime2 = Cphi(temp, 2); + printf(" Cphiprime2 = %g\n", cphiprime2); + + double beta0prime = Beta0(temp, 1); + printf(" beta0prime2 = %g\n", beta0prime); + + double beta1prime = Beta1(temp, 1); + printf(" beta1prime = %g\n", beta1prime); + + double cphiprime = Cphi(temp, 1); + printf(" Cphiprime = %g\n", cphiprime); + + double vm = 1.0; + double vx = 1.0; + double v = vm + vx; + double m = Iionic; + double zm = 1.; + double zx = 1.0; + + double sqrtI = sqrt(Iionic); + + double alpha = 2.0; + double a2 = alpha * alpha; + double b = 1.2; + + double Bpmx = beta0prime + 2.0 * beta1prime / (a2* Iionic) * + (1.0 - (1.0 + alpha * sqrtI) * exp(-alpha*sqrtI) ); + + double Bppmx = beta0prime2 + 2.0 * beta1prime2 / (a2* Iionic) * + (1.0 - (1.0 + alpha * sqrtI) * exp(-alpha*sqrtI) ); + + double Cpmx = 0.5 * sqrt(vm * vx) * cphiprime; + + double Cppmx = 0.5 * sqrt(vm * vx) * cphiprime2; + + double Bmx = Bppmx + 2.0 / temp * Bpmx; + double Cmx = Cppmx + 2.0 / temp * Cpmx; + + double phiJ = v * zm * zx * (A_J/(2.*b)) * log(1 + 1.2 * sqrtI) - + 2 * vm * vx * GasConst * temp * temp * ( m * Bmx + m * m * Cmx); + phiJ *= 1.0E-3; + printf(" phiJ = %15.8g kJ/gmolSalt\n", phiJ); + + double molecWeight = 18.01528; + + double RT = GasConst * temp * 1.0E-3; + + + double xo = 1.0 / (molecWeight/1000. * 2 * m + 1.0); + printf(" no = %g\n", xo); + + return phiJ; +} + +main() { + + printf("Standalone test of the apparent relative molal enthalpy, phiL:\n"); + printf(" (Check against simple formula in Silvester&Pitzer, J. Phys. Chem. 81, 1822 (1977)\n"); + + printf("T = 50C\n"); + double Iionic = 6.146; + printf("Ionic Strength = %g\n", Iionic); + + double res = calc(273.15 + 50., Iionic); + printf("T = 200C\n"); + printf("Ionic Strength = %g\n", Iionic); + + res = calc(273.15 + 200., Iionic); + +} diff --git a/test_problems/cathermo/HMW_graph_CpvT/HMW_NaCl_sp1977_alt.xml b/test_problems/cathermo/HMW_graph_CpvT/HMW_NaCl_sp1977_alt.xml new file mode 100644 index 000000000..1ca6c066a --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/HMW_NaCl_sp1977_alt.xml @@ -0,0 +1,243 @@ + + + + + + H2O(L) Cl- H+ Na+ OH- + + + 298.15 + 101325.0 + + Na+:6.0954 + Cl-:6.0954 + H+:2.1628E-9 + OH-:1.3977E-6 + + + + + + + + + + + + 0.0765, 0.008946, -3.3158E-6, + -777.03, -4.4706 + + 0.2664, 6.1608E-5, 1.0715E-6 + 0.0 + 0.00127, -4.655E-5, 0.0, + 33.317, 0.09421 + + 2.0 + + + + 0.1775, 0.0, 0.0, 0.0, 0.0 + 0.2945, 0.0, 0.0 + 0.0 + 0.0008, 0.0, 0.0, 0.0, 0.0 + 2.0 + + + + 0.0864, 0.0, 0.0, 0.0, 0.0 + 0.253, 0.0, 0.0 + 0.0 + 0.0044, 0.0, 0.0, 0.0, 0.0 + 2.0 + + + + -0.05 + + + + -0.05 + -0.006 + + + + 0.036 + + + + 0.036 + -0.004 + + + + H2O(L) + + O H C Fe Si N Na Cl + + + + + + + + + + H:2 O:1 + + + + 7.255750050E+01, -6.624454020E-01, 2.561987460E-03, -4.365919230E-06, + 2.781789810E-09, -4.188654990E+04, -2.882801370E+02 + + + + + + 0.018068 + + + + + + Na:1 + +1 + + + + -57993.47558 , 305112.6040 , -592222.1591 , + 401977.9827 , 804.4195980 , 10625.24901 , + -133796.2298 + + + + + + + 0.00834 + + + + + + Cl:1 + -1 + + + + 0.00834 + + + + + 56696.2042 , -297835.978 , 581426.549 , + -401759.991 , -804.301136 , -10873.8257 , + 130650.697 + + + + + + + + H:1 + +1 + + 0.0 + + + + 0.0 + 3 + + 0.0 , 0.0, 0.0 + + + 273.15, 298.15 , 623.15 + + + + + + + + O:1 H:1 + -1 + + + 0.00834 + + + + + 44674.99961 , -234943.0414 , 460522.8260 , + -320695.1836 , -638.5044716 , -8683.955813 , + 102874.2667 + + + + + + + + diff --git a/test_problems/cathermo/HMW_graph_CpvT/HMW_graph_CpvT.cpp b/test_problems/cathermo/HMW_graph_CpvT/HMW_graph_CpvT.cpp new file mode 100755 index 000000000..5133b2f81 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/HMW_graph_CpvT.cpp @@ -0,0 +1,351 @@ +/** + * + * @file HMW_graph_1.cpp + */ + +/* + * $Author$ + * $Date$ + * $Revision$ + */ +#include + +#ifdef SRCDIRTREE +#include "ct_defs.h" +#include "logger.h" +#include "TemperatureTable.h" +#include "ThermoPhase.h" +#include "HMWSoln.h" +#include "importCTML.h" +#else +#include "cantera/Cantera.h" +#include "cantera/kernel/logger.h" +#include "cantera/thermo.h" +#include "TemperatureTable.h" +#include "ThermoPhase.h" +#include "HMWSoln.h" +#include "importCTML.h" +#endif + +using namespace Cantera; + +class fileLog: public Logger { +public: + fileLog(string fName) { + m_fName = fName; + m_fs.open(fName.c_str()); + } + + virtual void write(const string& msg) { + m_fs << msg; + m_fs.flush(); + } + + virtual ~fileLog() { + m_fs.close(); + } + + string m_fName; + ofstream m_fs; + +}; + +void printUsage() { + cout << "usage: HMW_test " << endl; + cout <<" -> Everything is hardwired" << endl; +} + +void pAtable(HMWSoln *HMW) { + int nsp = HMW->nSpecies(); + double acMol[100]; + double mf[100]; + double activities[100]; + double moll[100]; + + HMW->getMolalityActivityCoefficients(acMol); + HMW->getMoleFractions(mf); + HMW->getActivities(activities); + HMW->getMolalities(moll); + string sName; + printf(" Name Activity ActCoeffMolal " + " MoleFract Molality\n"); + for (int k = 0; k < nsp; k++) { + sName = HMW->speciesName(k); + printf("%16s %13g %13g %13g %13g\n", + sName.c_str(), activities[k], acMol[k], mf[k], moll[k]); + } + +} + +int main(int argc, char **argv) +{ + + int retn = 0; + int i; + + try { + + char iFile[80]; + strcpy(iFile, "HMW_NaCl.xml"); + if (argc > 1) { + strcpy(iFile, argv[1]); + } + double Cp0_R[20], pmCp[20]; + + //fileLog *fl = new fileLog("HMW_graph_1.log"); + //setLogger(fl); + + HMWSoln *HMW = new HMWSoln(iFile, "NaCl_electrolyte"); + + + /* + * Load in and initialize the + */ + Cantera::ThermoPhase *solid = newPhase("NaCl_Solid.xml","NaCl(S)"); + + + int nsp = HMW->nSpecies(); + double acMol[100]; + double act[100]; + double mf[100]; + double moll[100]; + for (i = 0; i < 100; i++) { + acMol[i] = 1.0; + act[i] = 1.0; + mf[i] = 0.0; + act[i] = 0.0; + } + + HMW->getMoleFractions(mf); + string sName; + + TemperatureTable TTable(15, false, 273.15, 25., 0, 0); + + + HMW->setState_TP(298.15, 1.01325E5); + + int i1 = HMW->speciesIndex("Na+"); + int i2 = HMW->speciesIndex("Cl-"); + //int i3 = HMW->speciesIndex("H2O(L)"); + for (i = 1; i < nsp; i++) { + moll[i] = 0.0; + } + HMW->setMolalities(moll); + + double ISQRT; + double Is = 0.0; + + /* + * Set the Pressure + */ + double pres = OneAtm; + + /* + * Fix the molality + */ + Is = 6.146; + ISQRT = sqrt(Is); + moll[i1] = Is; + moll[i2] = Is; + HMW->setState_TPM(298.15, pres, moll); + double Xmol[30]; + HMW->getMoleFractions(Xmol); + + /* + * ThermoUnknowns + */ + double T; + + double Cp0_NaCl, Cp0_Naplus, Cp0_Clminus, Delta_Cp0s, Cp0_H2O; + double Cp_NaCl, Cp_Naplus, Cp_Clminus, Cp_H2O; + double molarCp0; +#ifdef DEBUG_HKM + FILE *ttt = fopen("table.csv","w"); +#endif + printf("A_J/R: Comparison to Pitzer's book, p. 99, can be made.\n"); + printf(" Agreement is within 12 pc \n"); + printf("\n"); + + printf("Delta_Cp0: Heat Capacity of Solution per mole of salt (standard states)\n"); + printf(" rxn for the ss heat of soln: " + "NaCl(s) -> Na+(aq) + Cl-(aq)\n"); + + printf("\n"); + printf("Delta_Cps: Delta heat Capacity of Solution per mole of salt\n"); + printf(" rxn for heat of soln: " + " n1 H2O(l,pure) + n2 NaCl(s) -> n2 MX(aq) + n1 H2O(l) \n"); + printf(" Delta_Hs = (n1 h_H2O_bar + n2 h_MX_bar " + "- n1 h_H2O_0 - n2 h_MX_0)/n2\n"); + printf("\n"); + printf("phiJ: phiJ, calculated from the program, is checked\n"); + printf(" against analytical formula in J_standalone program.\n"); + printf(" (comparison against Eq. 12, Silvester and Pitzer)\n"); + + /* + * Create a Table of NaCl Enthalpy Properties as a Function + * of the Temperature + */ + printf("\n\n"); + printf(" T, Pres, Aphi, A_J/R," + " Delta_Cp0," + " Delta_Cps, J, phiJ," + " MolarCp, MolarCp0\n"); + printf(" Kelvin, bar, sqrt(kg/gmol), sqrt(kg/gmol)," + " kJ/gmolSalt," + " kJ/gmolSalt, kJ/gmolSoln, kJ/gmolSalt," + " kJ/gmol, kJ/gmol\n"); +#ifdef DEBUG_HKM + fprintf(ttt,"T, Pres, A_J/R, Delta_Cp0, Delta_Cps, J, phiJ\n"); + fprintf(ttt,"Kelvin, bar, sqrt(kg/gmol), kJ/gmolSalt, kJ/gmolSalt, kJ/gmolSoln," + "kJ/gmolSalt\n"); +#endif + for (i = 0; i < TTable.NPoints + 1; i++) { + if (i == TTable.NPoints) { + T = 323.15; + } else { + T = TTable.T[i]; + } + /* + * RT is in units of J/kmolK + */ + //double RT = GasConstant * T; + + /* + * Make sure we are at the saturation pressure or above. + */ + + double psat = HMW->satPressure(T); + + pres = OneAtm; + if (psat > pres) pres = psat; + + + HMW->setState_TPM(T, pres, moll); + + solid->setState_TP(T, pres); + + /* + * Get the Standard State DeltaH + */ + + solid->getCp_R(Cp0_R); + Cp0_NaCl = Cp0_R[0] * GasConstant * 1.0E-6; + + + HMW->getCp_R(Cp0_R); + Cp0_H2O = Cp0_R[0] * GasConstant * 1.0E-6; + Cp0_Naplus = Cp0_R[i1] * GasConstant * 1.0E-6; + Cp0_Clminus = Cp0_R[i2] * GasConstant * 1.0E-6; + /* + * Calculate the standard state heat of solution + * for NaCl(s) -> Na+ + Cl- + * units: kJ/gmolSalt + */ + + Delta_Cp0s = Cp0_Naplus + Cp0_Clminus - Cp0_NaCl; + + pmCp[0] = solid->cp_mole(); + + Cp_NaCl = pmCp[0] * 1.0E-6; + + + HMW->getPartialMolarCp(pmCp); + Cp_H2O = pmCp[0] * 1.0E-6; + Cp_Naplus = pmCp[i1] * 1.0E-6; + Cp_Clminus = pmCp[i2] * 1.0E-6; + + //double Delta_Cp_Salt = Cp_NaCl - (Cp_Naplus + Cp_Clminus); + + double molarCp = HMW->cp_mole() * 1.0E-6; + + /* + * Calculate the heat capacity of solution for the reaction + * NaCl(s) -> Na+ + Cl- + */ + double Delta_Cps = (Xmol[0] * Cp_H2O + + Xmol[i1] * Cp_Naplus + + Xmol[i2] * Cp_Clminus + - Xmol[0] * Cp0_H2O + - Xmol[i1] * Cp_NaCl); + Delta_Cps /= Xmol[i1]; + + + /* + * Calculate the relative heat capacity, J, from the + * partial molar quantities, units J/gmolSolutionK + */ + double J = (Xmol[0] * (Cp_H2O - Cp0_H2O) + + Xmol[i1] * (Cp_Naplus - Cp0_Naplus) + + Xmol[i2] * (Cp_Clminus - Cp0_Clminus)); + + /* + * Calculate the apparent relative molal heat capacity, phiJ, + * units of J/gmolSaltAddedK + */ + double phiJ = J / Xmol[i1]; + + + double Aphi = HMW->A_Debye_TP(T, pres) / 3.0; + //double AL = HMW->ADebye_L(T,pres); + double AJ = HMW->ADebye_J(T, pres); + + + + for (int k = 0; k < nsp; k++) { + Cp0_R[k] *= GasConstant * 1.0E-6; + } + + molarCp0 = 0.0; + for (int k = 0; k < nsp; k++) { + molarCp0 += Xmol[k] * Cp0_R[k]; + } + + if (i != TTable.NPoints+1) { + printf("%13g, %13g, %13g, %13g, %13g, %13g, " + "%13g, %13g, %13g, %13g\n", + T, pres*1.0E-5, Aphi, AJ/GasConstant, Delta_Cp0s, Delta_Cps, + J, phiJ, molarCp , molarCp0 ); +#ifdef DEBUG_HKM + fprintf(ttt,"%g, %g, %g, %g, %g, %g, %g\n", + T, pres*1.0E-5, AJ/GasConstant, Delta_Cp0s, Delta_Cps, J, phiJ); +#endif + } + + } + + printf("Breakdown of Heat Capacity Calculation at 323.15 K, 1atm:\n"); + + printf(" Species MoleFrac Molal Cp0 " + " partCp (partCp - Cp0)\n"); + printf(" H2O(L)"); + printf("%13g %13g %13g %13g %13g\n", Xmol[0], moll[0], Cp0_H2O , Cp_H2O, Cp_H2O-Cp0_H2O); + printf(" Na+ "); + printf("%13g %13g %13g %13g %13g\n", Xmol[i1], moll[i1], + Cp0_Naplus , Cp_Naplus, Cp_Naplus -Cp0_Naplus); + printf(" Cl- "); + printf("%13g %13g %13g %13g %13g\n", Xmol[i2], moll[i2], + Cp0_Clminus , Cp_Clminus, Cp_Clminus - Cp0_Clminus); + + printf(" NaCl(s)"); + printf("%13g %13g %13g %13g\n", 1.0, + Cp0_NaCl , Cp_NaCl, Cp_NaCl - Cp0_NaCl); + + + delete HMW; + HMW = 0; + delete solid; + solid = 0; + Cantera::appdelete(); + +#ifdef DEBUG_HKM + fclose(ttt); +#endif + return retn; + + } catch (CanteraError) { + printf("caught error\n"); + showErrors(); + Cantera::appdelete(); + return -1; + } +} diff --git a/test_problems/cathermo/HMW_graph_CpvT/Makefile.in b/test_problems/cathermo/HMW_graph_CpvT/Makefile.in new file mode 100644 index 000000000..9bac62d28 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/Makefile.in @@ -0,0 +1,116 @@ +#!/bin/sh + +############################################################################ +# +# Makefile to compile and link a C++ application to +# Cantera. +# +############################################################################# + +# addition to suffixes +.SUFFIXES : .d + +# the name of the executable program to be created +PROG_NAME = HMW_graph_CpvT + +# the object files to be linked together. List those generated from Fortran +# and from C/C++ separately +OBJS = HMW_graph_CpvT.o sortAlgorithms.o + +# Location of the current build. Will assume that tests are run +# in the source directory tree location +src_dir_tree = 1 + +# additional flags to be passed to the linker. If your program +# requires other external libraries, put them here +LINK_OPTIONS = @EXTRA_LINK@ + +############################################################################# + +# Check to see whether we are in the msvc++ environment +os_is_win = @OS_IS_WIN@ + +# Fortran libraries +FORT_LIBS = @FLIBS@ + +# the C++ compiler +CXX = @CXX@ + +# C++ compile flags +ifeq ($(src_dir_tree), 1) +CXX_FLAGS = -DSRCDIRTREE @CXXFLAGS@ +else +CXX_FLAGS = @CXXFLAGS@ +endif + +# Ending C++ linking libraries +LCXX_END_LIBS = @LCXX_END_LIBS@ + +# the directory where the Cantera libraries are located +CANTERA_LIBDIR=@buildlib@ + +# required Cantera libraries +CANTERA_LIBS = @LOCAL_LIBS@ -lctcxx + +# the directory where Cantera include files may be found. +ifeq ($(src_dir_tree), 1) +CANTERA_INCDIR=../../../Cantera/src +INCLUDES=-I$(CANTERA_INCDIR) -I$(CANTERA_INCDIR)/thermo +else +CANTERA_INCDIR=@ctroot@/build/include/cantera +INCLUDES=-I$(CANTERA_INCDIR) -I$(CANTERA_INCDIR)/kernel +endif + +# flags passed to the C++ compiler/linker for the linking step +LCXX_FLAGS = -L$(CANTERA_LIBDIR) @LOCAL_LIB_DIRS@ @CXXFLAGS@ + +# How to compile C++ source files to object files +.@CXX_EXT@.@OBJ_EXT@: + $(CXX) -c $< $(INCLUDES) $(CXX_FLAGS) + +# How to compile the dependency file +.cpp.d: + g++ -MM $(INCLUDES) $(CXX_FLAGS) $*.cpp > $*.d + +# List of dependency files to be created +DEPENDS=$(OBJS:.o=.d) + +# Program Name +PROGRAM = $(PROG_NAME)$(EXE_EXT) + +all: $(PROGRAM) .depends Cp_standalone + +$(PROGRAM): $(OBJS) $(CANTERA_LIBDIR)/libcantera.a \ + $(CANTERA_LIBDIR)/libcaThermo.a + $(CXX) -o $(PROGRAM) $(OBJS) $(LCXX_FLAGS) $(LINK_OPTIONS) \ + $(CANTERA_LIBS) @LIBS@ $(FORT_LIBS) \ + $(LCXX_END_LIBS) + +Cp_standalone: Cp_standalone.o + $(CXX) -o Cp_standalone Cp_standalone.o \ + $(LCXX_FLAGS) $(LINK_OPTIONS) $(LCXX_END_LIBS) + +# depends target -> forces recalculation of dependencies +depends: + @MAKE@ .depends + +.depends: $(DEPENDS) + cat $(DEPENDS) > .depends + +# Do the test -> For the windows vc++ environment, we have to skip checking on +# whether the program is uptodate, because we don't utilize make +# in that environment to build programs. +test: +ifeq ($(os_is_win), 1) +else + @MAKE@ $(PROGRAM) +endif + ./runtest + +clean: + $(RM) $(OBJS) $(PROGRAM) $(DEPENDS) .depends + ../../../bin/rm_cvsignore + (if test -d SunWS_cache ; then \ + $(RM) -rf SunWS_cache ; \ + fi ) + diff --git a/test_problems/cathermo/HMW_graph_CpvT/NaCl_Solid.xml b/test_problems/cathermo/HMW_graph_CpvT/NaCl_Solid.xml new file mode 100644 index 000000000..d711be8ff --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/NaCl_Solid.xml @@ -0,0 +1,39 @@ + + + + + + + + + O H C Fe Ca N Na Cl + + NaCl(S) + + 2.165 + + + + + + + + + + + Na:1 Cl:1 + + + + 50.72389, 6.672267, -2.517167, + 10.15934, -0.200675, -427.2115, + 130.3973 + + + + 2.165 + + + + + diff --git a/test_problems/cathermo/HMW_graph_CpvT/README b/test_problems/cathermo/HMW_graph_CpvT/README new file mode 100644 index 000000000..f22bf8f5e --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/README @@ -0,0 +1 @@ +Check on the enthalpy routines for HMWSoln diff --git a/test_problems/cathermo/HMW_graph_CpvT/TemperatureTable.h b/test_problems/cathermo/HMW_graph_CpvT/TemperatureTable.h new file mode 100644 index 000000000..70cb2bf56 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/TemperatureTable.h @@ -0,0 +1,129 @@ +/* + * $Id$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000, there is a non-exclusive license for use of this + * work by or on behalf of the U.S. Government. Export of this program + * may require a license from the United States Government. + */ + +#ifndef TEMPERATURE_TABLE_H +#define TEMPERATURE_TABLE_H +#include "sortAlgorithms.h" +//#include "mdp_allo.h" +#include +using std::vector; + +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ +/** + * This Class constructs a vector of temperature from which to make + * a table. + */ +class TemperatureTable { + +public: + int NPoints; + bool Include298; + double Tlow; //!< Min temperature for thermo data fit + double Thigh; //!< Max temperature for thermo table + double DeltaT; + vector T; + int numAddedTs; + vector AddedTempVector; +public: + /* + * Default constructor for TemperatureTable() + */ + TemperatureTable(const int nPts = 14, + const bool inc298 = true, + const double tlow = 300., + const double deltaT = 100., + const int numAdded = 0, + const double *addedTempVector = 0) : + NPoints(nPts), + Include298(inc298), + Tlow(tlow), + DeltaT(deltaT), + T(0), + numAddedTs(numAdded) { + /****************************/ + int i; + // AddedTempVector = mdp_alloc_dbl_1(numAdded, 0.0); + AddedTempVector.resize(numAdded, 0.0); + for (int i = 0; i < numAdded; i++) { + AddedTempVector[i] = addedTempVector[i]; + } + //mdp_copy_dbl_1(AddedTempVector, addedTempVector, numAdded); + // T = mdp_alloc_dbl_1(NPoints, 0.0); + T.resize(NPoints, 0.0); + double TCurrent = Tlow; + for (i = 0; i < NPoints; i++) { + T[i] = TCurrent; + TCurrent += DeltaT; + } + if (Include298) { + T.push_back(298.15); + //mdp_realloc_dbl_1(&T, NPoints+1, NPoints, 298.15); + NPoints++; + } + if (numAdded > 0) { + //mdp_realloc_dbl_1(&T, NPoints+numAdded, NPoints, 0.0); + T.resize( NPoints+numAdded, 0.0); + for (i = 0; i < numAdded; i++) { + T[i+NPoints] = addedTempVector[i]; + } + NPoints += numAdded; + } + + sort_dbl_1(DATA_PTR(T), NPoints); + + + } + /***********************************************************************/ + /***********************************************************************/ + /***********************************************************************/ + /* + * Destructor + */ + ~TemperatureTable() { + //mdp_safe_free((void **) &AddedTempVector); + // mdp_safe_free((void **) &T); + } + + /***********************************************************************/ + /***********************************************************************/ + /***********************************************************************/ + /* + * Overloaded operator[] + * + * return the array value in the vector + */ + double operator[](const int i) { + return T[i]; + } + /***********************************************************************/ + /***********************************************************************/ + /***********************************************************************/ + /* + * size() + */ + int size() { + return NPoints; + } +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ + /* + * Block assignment and copy constructors: not needed. + */ +private: + TemperatureTable(const TemperatureTable &); + TemperatureTable& operator=(const TemperatureTable&); +}; +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ +#endif diff --git a/test_problems/cathermo/HMW_graph_CpvT/output_blessed.txt b/test_problems/cathermo/HMW_graph_CpvT/output_blessed.txt new file mode 100644 index 000000000..29401e696 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/output_blessed.txt @@ -0,0 +1,39 @@ +A_J/R: Comparison to Pitzer's book, p. 99, can be made. + Agreement is within 12 pc + +Delta_Cp0: Heat Capacity of Solution per mole of salt (standard states) + rxn for the ss heat of soln: NaCl(s) -> Na+(aq) + Cl-(aq) + +Delta_Cps: Delta heat Capacity of Solution per mole of salt + rxn for heat of soln: n1 H2O(l,pure) + n2 NaCl(s) -> n2 MX(aq) + n1 H2O(l) + Delta_Hs = (n1 h_H2O_bar + n2 h_MX_bar - n1 h_H2O_0 - n2 h_MX_0)/n2 + +phiJ: phiJ, calculated from the program, is checked + against analytical formula in J_standalone program. + (comparison against Eq. 12, Silvester and Pitzer) + + + T, Pres, Aphi, A_J/R, Delta_Cp0, Delta_Cps, J, phiJ, MolarCp, MolarCp0 + Kelvin, bar, sqrt(kg/gmol), sqrt(kg/gmol), kJ/gmolSalt, kJ/gmolSalt, kJ/gmolSoln, kJ/gmolSalt, kJ/gmol, kJ/gmol + 273.15, 1.01325, 0.376717, 4.08576, -0.158977, -0.0364733, 0.0111048, 0.122504, 0.0634482, 0.0523434 + 298.15, 1.01325, 0.391447, 4.61226, -0.130795, -0.0199566, 0.0100473, 0.110838, 0.0644397, 0.0543924 + 323.15, 1.01325, 0.410293, 5.50264, -0.115719, -0.0105888, 0.00952988, 0.10513, 0.065338, 0.0558081 + 348.15, 1.01325, 0.433273, 6.72554, -0.113771, -0.00869263, 0.00952522, 0.105079, 0.065728, 0.0562028 + 373.15, 1.01418, 0.460559, 8.35444, -0.124962, -0.0135747, 0.0100971, 0.111387, 0.0656563, 0.0555592 + 398.15, 2.32238, 0.492454, 10.5438, -0.14929, -0.0237587, 0.0113792, 0.125532, 0.0653078, 0.0539285 + 423.15, 4.76165, 0.529514, 13.5812, -0.186752, -0.0364892, 0.0136211, 0.150263, 0.0649994, 0.0513783 + 448.15, 8.92602, 0.572549, 17.9961, -0.237338, -0.0466955, 0.0172814, 0.190643, 0.0652692, 0.0479877 + 473.15, 15.5493, 0.622769, 24.826, -0.301036, -0.0444443, 0.0232596, 0.256592, 0.0671305, 0.0438709 + 498.15, 25.4972, 0.682036, 36.2735, -0.377833, -0.00866191, 0.0334648, 0.369171, 0.0726787, 0.0392139 + 523.15, 39.7617, 0.753389, 57.5198, -0.467713, 0.110207, 0.0523875, 0.57792, 0.0867379, 0.0343503 + 548.15, 59.4639, 0.842213, 102.542, -0.570661, 0.445824, 0.0921428, 1.01649, 0.122089, 0.0299461 + 573.15, 85.879, 0.959258, 217.131, -0.68666, 1.43604, 0.19242, 2.1227, 0.219976, 0.0275567 + 598.15, 120.51, 1.13023, 604.241, -0.815692, 5.02147, 0.52913, 5.83716, 0.560947, 0.0318178 + 623.15, 165.294, 1.43872, 2813.78, -0.957741, 26.0219, 2.44566, 26.9797, 2.51311, 0.0674475 + 323.15, 1.01325, 0.410293, 5.50264, -0.115719, -0.0105888, 0.00952988, 0.10513, 0.065338, 0.0558081 +Breakdown of Heat Capacity Calculation at 323.15 K, 1atm: + Species MoleFrac Molal Cp0 partCp (partCp - Cp0) + H2O(L) 0.818703 0 0.075328 0.0695928 -0.00573522 + Na+ 0.0906484 6.146 0.0284152 0.106879 0.0784643 + Cl- 0.0906484 6.146 -0.0930958 -0.0146316 0.0784643 + NaCl(s) 1 0.0510383 0.0510383 0 diff --git a/test_problems/cathermo/HMW_graph_CpvT/runtest b/test_problems/cathermo/HMW_graph_CpvT/runtest new file mode 100755 index 000000000..cee966321 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/runtest @@ -0,0 +1,42 @@ +#!/bin/sh +# +# +temp_success="1" +/bin/rm -f output.txt outputa.txt + +########################################################################## +prog=HMW_graph_CpvT +if test ! -x $prog ; then + echo $prog ' does not exist' + exit -1 +fi +########################################################################## +/bin/rm -f test.out test.diff output.txt + +################################################################# +# +CANTERA_DATA=${CANTERA_DATA:=../../../data/inputs}; export CANTERA_DATA +CANTERA_BIN=${CANTERA_BIN:=../../../bin} + +################################################################# + +$prog HMW_NaCl_sp1977_alt.xml > output.txt +retnStat=$? +if [ $retnStat != "0" ] +then + temp_success="0" + echo "$prog returned with bad status, $retnStat, check output" +fi + +$CANTERA_BIN/exp3to2.sh output.txt > outputa.txt +diff -w outputa.txt output_blessed.txt > diff_test.out +retnStat=$? +if [ $retnStat = "0" ] +then + echo "successful diff comparison on $prog test" +else + echo "unsuccessful diff comparison on $prog test" + echo "FAILED" > csvCode.txt + temp_success="0" +fi + diff --git a/test_problems/cathermo/HMW_graph_CpvT/sortAlgorithms.cpp b/test_problems/cathermo/HMW_graph_CpvT/sortAlgorithms.cpp new file mode 100644 index 000000000..d97e51b40 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/sortAlgorithms.cpp @@ -0,0 +1,54 @@ +/* + * @file sortAlgorithms.h + * + * $Author$ + * $Revision$ + * $Date$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + +#include "sortAlgorithms.h" + +/**************************************************************/ + +void sort_dbl_1(double * const x, const int n) { + double rra; + int ll = n/2; + int iret = n - 1; + while (1 > 0) { + if (ll > 0) { + ll--; + rra = x[ll]; + } else { + rra = x[iret]; + x[iret] = x[0]; + iret--; + if (iret == 0) { + x[0] = rra; + return; + } + } + int i = ll; + int j = ll + ll + 1; + while (j <= iret) { + if (j < iret) { + if (x[j] < x[j+1]) + j++; + } + if (rra < x[j]) { + x[i] = x[j]; + i = j; + j = j + j + 1; + } else { + j = iret + 1; + } + } + x[i] = rra; + } +} +/*****************************************************/ diff --git a/test_problems/cathermo/HMW_graph_CpvT/sortAlgorithms.h b/test_problems/cathermo/HMW_graph_CpvT/sortAlgorithms.h new file mode 100644 index 000000000..72a7fc2a2 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_CpvT/sortAlgorithms.h @@ -0,0 +1,21 @@ +/* + * @file sortAlgorithms.h + * + * $Author$ + * $Revision$ + * $Date$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + +#ifndef SORTALGORITHMS_H +#define SORTALGORITHMS_H + + +void sort_dbl_1(double * const x, const int n); + +#endif