diff --git a/test_problems/cathermo/HMW_graph_HvT/.cvsignore b/test_problems/cathermo/HMW_graph_HvT/.cvsignore new file mode 100644 index 000000000..72f5b5172 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/.cvsignore @@ -0,0 +1,11 @@ +Makefile +output.txt +outputa.txt +.depends +HMW_NaCl_sp1977_alt.out +HMW_graph_HvT +HMW_graph_HvT.d +L_standalone +diff_test.out +sortAlgorithms.d +table.csv diff --git a/test_problems/cathermo/HMW_graph_HvT/HMW_NaCl_sp1977_alt.xml b/test_problems/cathermo/HMW_graph_HvT/HMW_NaCl_sp1977_alt.xml new file mode 100644 index 000000000..1ca6c066a --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/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_HvT/HMW_graph_HvT.cpp b/test_problems/cathermo/HMW_graph_HvT/HMW_graph_HvT.cpp new file mode 100644 index 000000000..f9c6fd009 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/HMW_graph_HvT.cpp @@ -0,0 +1,347 @@ +/** + * + * @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/thermo.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[30]; + double mf[30]; + + double activities[30]; + double moll[30]; + for (int i = 0; i < 30; i++) { + acMol[i] = 1.0; + mf[i] = 0.0; + activities[i] = 1.0; + moll[i] = 0.0; + } + + 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 Enth0_RT[20], pmEnth[20], molarEnth; + + //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; + mf[i] = 0.0; + act[i] = 1.0; + moll[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 H0_NaCl, H0_Naplus, H0_Clminus, Delta_H0s, H0_H2O; + double H_NaCl, H_Naplus, H_Clminus, H_H2O; + double molarEnth0; + + printf("A_L/RT: Comparison to Pitzer's book, p. 99, can be made.\n"); + printf(" Agreement to 3-4 sig digits for Aphi and A_L/RT\n"); + printf("\n"); + + printf("Delta_H0: Heat 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_Hs: Heat 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("phiL: phiL, calculated from the program, is checked\n"); + printf(" against analytical formula in L_standalone program.\n"); + printf(" (comparison against Eq. 12, Silvester and Pitzer)\n"); + +#ifdef DEBUG_HKM + FILE *ttt = fopen("table.csv","w"); +#endif + /* + * Create a Table of NaCl Enthalpy Properties as a Function + * of the Temperature + */ + printf("\n\n"); + printf(" T, Pres, Aphi, A_L/RT," + " Delta_H0," + " Delta_Hs, L, phiL," + " L_rel_molal," + " MolarEnth, MolarEnth0\n"); + printf(" Kelvin, bar, sqrt(kg/gmol), sqrt(kg/gmol)," + " kJ/gmolSalt," + " kJ/gmolSalt, kJ/gmolSoln, kJ/gmolSalt," + " kJ/gmolSalt, kJ/gmol, kJ/gmol\n"); +#ifdef DEBUG_HKM + fprintf(ttt,"T, Pres, A_L/RT, Delta_H0, Delta_Hs, phiL\n"); + fprintf(ttt,"Kelvin, bar, sqrt(kg/gmol), kJ/gmolSalt, kJ/gmolSalt, 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->getEnthalpy_RT(Enth0_RT); + H0_NaCl = Enth0_RT[0] * RT * 1.0E-6; + HMW->getEnthalpy_RT(Enth0_RT); + H0_H2O = Enth0_RT[0] * RT * 1.0E-6; + H0_Naplus = Enth0_RT[i1] * RT * 1.0E-6; + H0_Clminus = Enth0_RT[i2] * RT * 1.0E-6; + + /* + * Calculate the standard state heat of solution + * for NaCl(s) -> Na+ + Cl- + * units: kJ/gmolSalt + */ + Delta_H0s = H0_Naplus + H0_Clminus - H0_NaCl; + + + solid->getPartialMolarEnthalpies(pmEnth); + H_NaCl = pmEnth[0] * 1.0E-6; + HMW->getPartialMolarEnthalpies(pmEnth); + H_H2O = pmEnth[0] * 1.0E-6; + H_Naplus = pmEnth[i1] * 1.0E-6; + H_Clminus = pmEnth[i2] * 1.0E-6; + //double Delta_H_Salt = H_NaCl - (H_Naplus + H_Clminus); + //double Lfunc = HMW->relative_enthalpy() * 1.0E-6; + molarEnth = HMW->enthalpy_mole() * 1.0E-6; + + + + /* + * + */ + double Delta_Hs = (Xmol[0] * H_H2O + + Xmol[i1] * H_Naplus + + Xmol[i2] * H_Clminus + - Xmol[0] * H0_H2O + - Xmol[i1] * H_NaCl); + Delta_Hs /= Xmol[i1]; + + /* + * Calculate the relative enthalpy, L, from the + * partial molar quantities. units kJ/gmolSolution + */ + double L = (Xmol[0] * (H_H2O - H0_H2O) + + Xmol[i1] * (H_Naplus - H0_Naplus) + + Xmol[i2] * (H_Clminus - H0_Clminus)); + + /* + * Calculate the apparent relative molal enthalpy, phiL, + * units of kJ/gmolSaltAdded + */ + double phiL = L / Xmol[i1]; + + + double Aphi = HMW->A_Debye_TP() / 3.0; + double AL = HMW->ADebye_L(); + double LrelMol = HMW->relative_molal_enthalpy() * 1.0E-6; + + + + for (int k = 0; k < nsp; k++) { + Enth0_RT[k] *= RT * 1.0E-6; + } + + molarEnth0 = 0.0; + for (int k = 0; k < nsp; k++) { + molarEnth0 += Xmol[k] * Enth0_RT[k]; + } + + if (i != TTable.NPoints+1) { + printf("%13g, %13g, %13g, %13g, %13g, %13g, %13g, " + "%13g, %13g, %13g, %13g\n", + T, pres*1.0E-5, Aphi, AL/RT, Delta_H0s, Delta_Hs, + L, phiL, + LrelMol, + molarEnth , molarEnth0 ); + +#ifdef DEBUG_HKM + fprintf(ttt,"%g, %g, %g, %g, %g, %g\n", + T, pres*1.0E-5, AL/RT, Delta_H0s, Delta_Hs, phiL); +#endif + } + + } + + printf("Breakdown of Enthalpy Calculation at 323.15 K, 1atm:\n"); + + printf(" Species MoleFrac Molal H0 " + " partH (partH - H0)\n"); + printf(" H2O(L)"); + printf("%13g %13g %13g %13g %13g\n", Xmol[0], moll[0], H0_H2O , H_H2O, H_H2O-H0_H2O); + printf(" Na+ "); + printf("%13g %13g %13g %13g %13g\n", Xmol[i1], moll[i1], + H0_Naplus , H_Naplus, H_Naplus -H0_Naplus); + printf(" Cl- "); + printf("%13g %13g %13g %13g %13g\n", Xmol[i2], moll[i2], + H0_Clminus , H_Clminus, H_Clminus - H0_Clminus); + + + + 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_HvT/L_standalone.cpp b/test_problems/cathermo/HMW_graph_HvT/L_standalone.cpp new file mode 100644 index 000000000..876e16c1c --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/L_standalone.cpp @@ -0,0 +1,168 @@ + +#include +#include +#include + +using namespace std; + + +/* + * Values of A_L/RT : tabular form + * units sqrt(kg/gmol), + */ +double A_LdRT(double temp) { + double retn; + if (temp == 323.15) { + retn = 1.08059671; + } else if (temp == 473.15) { + retn = 4.17423374; + } else { + printf("A_LdRT 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 { + 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 { + 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 { + 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_L = A_LdRT(temp); + A_L *= GasConst * temp; + + double beta0prime = Beta0(temp, 1); + printf(" beta0prime = %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 Cpmx = 0.5 * sqrt(vm * vx) * cphiprime; + + double phiL = v * zm * zx * (A_L/(2.*b)) * log(1 + 1.2 * sqrtI) - + 2 * vm * vx * GasConst * temp * temp * ( m * Bpmx + m * m * Cpmx); + phiL *= 1.0E-3; + printf(" phiL = %15.8g kJ/gmolSalt\n", phiL); + + double molecWeight = 18.01528; + + double RT = GasConst * temp * 1.0E-3; + //double IdealMixing = 2.0 * RT * m * (log(m) - 1.0); + //printf(" IdealMixing = %20.13g kJ/kg_water\n", IdealMixing); + + + double xo = 1.0 / (molecWeight/1000. * 2 * m + 1.0); + printf(" no = %g\n", xo); + + return phiL; +} + +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_HvT/Makefile.in b/test_problems/cathermo/HMW_graph_HvT/Makefile.in new file mode 100644 index 000000000..fceafb1b8 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/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_HvT + +# the object files to be linked together. List those generated from Fortran +# and from C/C++ separately +OBJS = HMW_graph_HvT.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 L_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) + +L_standalone: L_standalone.o + $(CXX) -o L_standalone L_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_HvT/NaCl_Solid.xml b/test_problems/cathermo/HMW_graph_HvT/NaCl_Solid.xml new file mode 100644 index 000000000..d711be8ff --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/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_HvT/README b/test_problems/cathermo/HMW_graph_HvT/README new file mode 100644 index 000000000..f22bf8f5e --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/README @@ -0,0 +1 @@ +Check on the enthalpy routines for HMWSoln diff --git a/test_problems/cathermo/HMW_graph_HvT/TemperatureTable.h b/test_problems/cathermo/HMW_graph_HvT/TemperatureTable.h new file mode 100644 index 000000000..70cb2bf56 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/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_HvT/output_blessed.txt b/test_problems/cathermo/HMW_graph_HvT/output_blessed.txt new file mode 100644 index 000000000..3f980dc4d --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/output_blessed.txt @@ -0,0 +1,38 @@ +A_L/RT: Comparison to Pitzer's book, p. 99, can be made. + Agreement to 3-4 sig digits for Aphi and A_L/RT + +Delta_H0: Heat of Solution per mole of salt (standard states) + rxn for the ss heat of soln: NaCl(s) -> Na+(aq) + Cl-(aq) + +Delta_Hs: Heat 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 + +phiL: phiL, calculated from the program, is checked + against analytical formula in L_standalone program. + (comparison against Eq. 12, Silvester and Pitzer) + + + T, Pres, Aphi, A_L/RT, Delta_H0, Delta_Hs, L, phiL, L_rel_molal, MolarEnth, MolarEnth0 + Kelvin, bar, sqrt(kg/gmol), sqrt(kg/gmol), kJ/gmolSalt, kJ/gmolSalt, kJ/gmolSoln, kJ/gmolSalt, kJ/gmolSalt, kJ/gmol, kJ/gmol + 273.15, 1.01325, 0.376717, 0.553555, 7.43878, 2.95802, -0.406173, -4.48075, -4.48075, -272.676, -272.27 + 298.15, 1.01325, 0.391447, 0.801126, 3.84389, 2.07326, -0.160505, -1.77063, -1.77063, -271.096, -270.936 + 323.15, 1.01325, 0.410293, 1.08058, 0.789801, 1.55964, 0.0697849, 0.769841, 0.769841, -269.487, -269.556 + 348.15, 1.01325, 0.433273, 1.39763, -2.05146, 1.19304, 0.294109, 3.24451, 3.24451, -267.86, -268.154 + 373.15, 1.01418, 0.460559, 1.76368, -5.00826, 0.774675, 0.524213, 5.78293, 5.78293, -266.231, -266.755 + 398.15, 2.32238, 0.492454, 2.1941, -8.41039, 0.130309, 0.774201, 8.5407, 8.5407, -264.608, -265.383 + 423.15, 4.76165, 0.529514, 2.71189, -12.5861, -0.867087, 1.06231, 11.719, 11.719, -262.999, -264.062 + 448.15, 8.92602, 0.572549, 3.35237, -17.8642, -2.27038, 1.41355, 15.5938, 15.5938, -261.4, -262.814 + 473.15, 15.5493, 0.622769, 4.17418, -24.5734, -3.99315, 1.86557, 20.5802, 20.5802, -259.793, -261.659 + 498.15, 25.4972, 0.682036, 5.28228, -33.0422, -5.67195, 2.48107, 27.3703, 27.3703, -258.132, -260.613 + 523.15, 39.7617, 0.753389, 6.88208, -43.599, -6.34036, 3.37744, 37.2587, 37.2587, -256.312, -259.69 + 548.15, 59.4639, 0.842213, 9.42265, -56.5718, -3.55817, 4.8056, 53.0137, 53.0137, -254.086, -258.892 + 573.15, 85.879, 0.959258, 14.046, -72.2884, 9.3421, 7.39968, 81.6305, 81.6305, -250.808, -258.208 + 598.15, 120.51, 1.13023, 24.3993, -91.0764, 54.4473, 13.1915, 145.524, 145.524, -244.395, -257.587 + 623.15, 165.294, 1.43872, 57.186, -113.263, 235.127, 31.581, 348.39, 348.39, -225.264, -256.845 + 323.15, 1.01325, 0.410293, 1.08058, 0.789801, 1.55964, 0.0697849, 0.769841, 0.769841, -269.487, -269.556 +Breakdown of Enthalpy Calculation at 323.15 K, 1atm: + Species MoleFrac Molal H0 partH (partH - H0) + H2O(L) 0.818703 0 -283.956 -284.104 -0.147946 + Na+ 0.0906484 6.146 -239.506 -238.453 1.05302 + Cl- 0.0906484 6.146 -169.555 -168.502 1.05302 diff --git a/test_problems/cathermo/HMW_graph_HvT/runtest b/test_problems/cathermo/HMW_graph_HvT/runtest new file mode 100755 index 000000000..c0f7ca57f --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/runtest @@ -0,0 +1,42 @@ +#!/bin/sh +# +# +temp_success="1" +/bin/rm -f output.txt outputa.txt + +########################################################################## +prog=HMW_graph_HvT +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_HvT/sortAlgorithms.cpp b/test_problems/cathermo/HMW_graph_HvT/sortAlgorithms.cpp new file mode 100644 index 000000000..d97e51b40 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/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_HvT/sortAlgorithms.h b/test_problems/cathermo/HMW_graph_HvT/sortAlgorithms.h new file mode 100644 index 000000000..72a7fc2a2 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_HvT/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