diff --git a/test_problems/cathermo/HMW_graph_GvT/.cvsignore b/test_problems/cathermo/HMW_graph_GvT/.cvsignore new file mode 100644 index 000000000..aa44ea6f1 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/.cvsignore @@ -0,0 +1,11 @@ +Makefile +.cvsignore.swp +.depends +Gex_standalone +HMW_graph_GvT +HMW_graph_GvT.d +diff_test.out +output.txt +outputa.txt +sortAlgorithms.d + diff --git a/test_problems/cathermo/HMW_graph_GvT/Gex_standalone.cpp b/test_problems/cathermo/HMW_graph_GvT/Gex_standalone.cpp new file mode 100644 index 000000000..04e7e3e9e --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/Gex_standalone.cpp @@ -0,0 +1,151 @@ + +#include +#include +#include + +using namespace std; + + +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); + } + 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); + } + 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); + } + return retn; +} + +double calc(double temp, double Iionic) { + + 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); + + double beta0 = Beta0(temp, 0); + printf(" beta0 = %g\n", beta0); + + double beta1 = Beta1(temp, 0); + printf(" beta1 = %g\n", beta1); + + double cphi = Cphi(temp, 0); + printf(" Cphi = %g\n", cphi); + + 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 b = 1.2; + + double osm1 = - zm * zx * Aphi * sqrtI / (1.0 + b * sqrtI) + + m * 2.0 * vm * vx / v * (beta0 + beta1 * exp( - alpha * sqrtI)) + + m * m * 2 * pow((vm*vx), 1.5) / v * cphi; + + double os = osm1 + 1.0; + double a2 = alpha * alpha; + + + printf("osmotic coeff = %20.13g\n", os); + + double lnmeanAct = - zm * zx * Aphi * + (sqrtI / (1.0 + b * sqrtI) + 2 / b * log (1.0 + b * sqrtI)) + + m * 2 * vm * vx / v * (2.0 * beta0 + + 2.0 * beta1/ (a2 * Iionic) + * (1.0- (1.0 + alpha* sqrtI - a2*Iionic/2.0)* exp(- alpha * sqrtI)) + ) + + 3.0 * m * m / 2.0 * (2 * sqrt(vm * vx) * vm * vx / v) * cphi; + + printf("ln(meanac) = %20.13g\n", lnmeanAct); + + double actCoeff = exp(lnmeanAct); + + printf("actCoeff = %20.13g\n", actCoeff); + + /* + * Gas constant in J gmol-1 K-1 + */ + double GasConst = 8.314472; + double gex = v * m * GasConst * 1.0E-3 * temp * (-osm1 + lnmeanAct); + printf(" Gex = %20.13g kJ/kg_water\n", gex); + + 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 DelG = gex + IdealMixing; + printf(" G - G0 = %20.13g kJ/kg_water\n", DelG); + + double mu0[6], mu[6]; + mu0[0] = -307.76256; + double molecWeight = 18.01528; + double diff = - RT * molecWeight / 1000. * 2. * m * os; + mu[0] = mu0[0] + diff; + printf("mus_kJ/gmol - H2O(L) - %20.13g %20.13g\n", mu0[0], mu[0]); + printf(" diff = %20.14g\n", diff); + double xo = 1.0 / (molecWeight/1000. * 2 * m + 1.0); + printf(" no = %g\n", xo); + double tmp = diff / (RT); + double actCoefWater = exp(tmp) / xo; + printf("actCoefWater = %g\n", actCoefWater); + return gex; +} + +main() { + + printf("standalone test of Gibbs excess free energy:\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"); + + res = calc(273.15 + 200., Iionic); +} diff --git a/test_problems/cathermo/HMW_graph_GvT/HMW_NaCl_sp1977_alt.xml b/test_problems/cathermo/HMW_graph_GvT/HMW_NaCl_sp1977_alt.xml new file mode 100644 index 000000000..1ca6c066a --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/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_GvT/HMW_graph_GvT.cpp b/test_problems/cathermo/HMW_graph_GvT/HMW_graph_GvT.cpp new file mode 100644 index 000000000..6f10c9923 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/HMW_graph_GvT.cpp @@ -0,0 +1,319 @@ +/** + * + * @file HMW_graph_1.cpp + */ + +/* + * $Author$ + * $Date$ + * $Revision$ + */ +#include + +#ifdef SRCDIRTREE +#include "ct_defs.h" +#include "logger.h" +#include "ThermoPhase.h" +#include "HMWSoln.h" +#include "importCTML.h" +#else +#include "ThermoPhase.h" + +#include "cantera/Cantera.h" +#include "cantera/kernel/logger.h" +#include "cantera/thermo.h" +#include "cantera/kernel/thermo/HMWSoln.h" +#endif + +#include "TemperatureTable.h" + + +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) { + double acMol[30]; + double mf[30]; + double activities[30]; + double moll[30]; + int nsp = HMW->nSpecies(); + + for (int i = 0; i < 30; i++) { + acMol[i] = 1.0; + mf[i] = 0.0; + activities[i] = 0.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; + int extraCols = 1; + + try { + //Cantera::ThermoPhase *tp = 0; + char iFile[80]; + strcpy(iFile, "HMW_NaCl.xml"); + if (argc > 1) { + strcpy(iFile, argv[1]); + } + + //fileLog *fl = new fileLog("HMW_graph_1.log"); + //setLogger(fl); + + HMWSoln *HMW = new HMWSoln(iFile, "NaCl_electrolyte"); + + + /* + * Load in and initialize the + */ + string nacl_s = "NaCl_Solid.xml"; + string id = "NaCl(S)"; + Cantera::ThermoPhase *solid = Cantera::newPhase(nacl_s, id); + + 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; + moll[i] = 0.0; + } + + HMW->getMoleFractions(mf); + string sName; + + TemperatureTable TTable(29, true, 293.15, 10., 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 using the setState_TPM() function. + */ + 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); + + printf("Fixed Concentration of the System:\n"); + + printf(" Species Mole_Fraction Molality\n"); + printf(" Na+ %g %g\n", Xmol[i1], moll[i1]); + printf(" Cl- %g %g\n", Xmol[i2], moll[i2]); + printf(" H2O(L) %g \n", Xmol[0]); + printf("\n"); + /* + * ThermoUnknowns + */ + double mu0_RT[20], mu[20]; + double mu0_NaCl, mu0_Naplus, mu0_Clminus, Delta_G0; + double mu_NaCl, mu_Naplus, mu_Clminus, Delta_G; + double molarGibbs0, molarGibbs; + + /* + * Create a Table of NaCl Enthalpy Properties as a Function + * of the Temperature + */ + + printf(" Table at fixed molality(Delta_G refers to rxn, NaCl(s) -> Na+ + Cl-)\n"); + printf(" -> pressure follows the saturation pressure above one atmosphere)\n"); + printf(" -> This calculation is meant to test Gibbs_ex -> TODO\n"); + printf("\n"); + + printf(" T, Pres, Aphi, Delta_G0," + " Delta_G," + " molarGibbs0, molarGibbs, Gibbs_ex," + " meanAC_moll, OsmCoeff-1"); + if (extraCols) { + printf(", Gibbs_ex_Formula, IdealMixing"); + } + printf("\n"); + printf(" Kelvin, bars, sqrt(kg/gmol), kJ/gmol," + " kJ/gmol," + " kJ/kgWater, kJ/kgWater, kJ/kgWater," + " , "); + if (extraCols) { + printf(", kJ/kgWater, kJ/kgWater "); + } + printf("\n"); + for (i = 0; i < TTable.NPoints; i++) { + double T = TTable.T[i]; + double RT = GasConstant * T; + + 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->getGibbs_RT(mu0_RT); + mu0_NaCl = mu0_RT[0] * RT * 1.0E-6; + + HMW->getGibbs_RT(mu0_RT); + //double mu0_water = mu0_RT[0] * RT * 1.0E-6; + mu0_Naplus = mu0_RT[i1] * RT * 1.0E-6; + mu0_Clminus = mu0_RT[i2] * RT * 1.0E-6; + Delta_G0 = (mu0_Naplus + mu0_Clminus) - mu0_NaCl; + + HMW->getMolalityActivityCoefficients(acMol); + HMW->getActivities(act); + + double meanAC = sqrt(acMol[i1] * acMol[i2]); + solid->getChemPotentials(mu); + + mu_NaCl = mu[0] * 1.0E-6; + + HMW->getChemPotentials(mu); + for (int k = 0; k < nsp; k++) { + mu[k] *= 1.0E-6; + } + mu_Naplus = mu[i1]; + mu_Clminus = mu[i2]; + Delta_G = (mu_Naplus + mu_Clminus) - mu_NaCl; + + + molarGibbs = HMW->gibbs_mole() * 1.0E-6; + /* + * Now the molarGibbs value is based on a mole of + * solution. This is useless for comparison purposes. + * Change to kg Water + */ + double molecWater = HMW->molecularWeight(0); + double Mo = molecWater / 1000.; + double Gibbs_kgWater = molarGibbs / (Xmol[0] * Mo); + + + double Aphi = HMW->A_Debye_TP() / 3.0; + + for (int k = 0; k < nsp; k++) { + mu0_RT[k] *= RT * 1.0E-6; + } + + molarGibbs0 = 0.0; + for (int k = 0; k < nsp; k++) { + molarGibbs0 += Xmol[k] * mu0_RT[k]; + } + double Gibbs0_kgWater = molarGibbs0 / (Xmol[0] * Mo); + + double osm1 = HMW->osmoticCoefficient(); + osm1 = osm1 - 1.0; + /* + * Need the gas constant in kJ/gmolK + */ + double rgas = 8.314472 * 1.0E-3; + double IdealMixing = moll[i1] * 2.0 * rgas * T * (log(moll[i1]) - 1.0); + + + double G_ex_kgWater = Gibbs_kgWater - Gibbs0_kgWater - IdealMixing; + + /* + * Calcualte excess gibbs free energy from another formula + */ + double G_ex_formula = 2 * Is * rgas * T * ( - osm1 + log(meanAC)); +/* + if (fabs (T-298.15) < 1.0) { + printf("mu0_Naplus = %g\n", mu0_Naplus); + printf("mu0_Clminus = %g\n", mu0_Clminus); + printf("mu0_NaCl(s) = %g, mu_NaCl(s) = %g\n",mu0_NaCl, mu_NaCl); + } +*/ + double pbar = pres * 1.0E-5; + + //if (extraCols && T == 323.15) { + // for (int k = 0; k < nsp; k++) { + // printf("mus_kJ/gmol - %s - %14.8g %14.8g %g\n", + // HMW->speciesName(k).c_str(), mu0_RT[k], mu[k], Xmol[k]); + // } + //} + + printf("%10g, %10g, %12g, %12g, %12g, %12g, %12g, %12g, %14.9g, %14.9g", + T, pbar, Aphi, Delta_G0, Delta_G, Gibbs0_kgWater, Gibbs_kgWater, G_ex_kgWater, + meanAC, osm1); + if (extraCols) { + printf(", %12g, %12g", G_ex_formula, IdealMixing); + } + printf("\n"); + } + + + delete HMW; + HMW = 0; + delete solid; + solid = 0; + Cantera::appdelete(); + + return retn; + + } catch (CanteraError) { + + showErrors(); + Cantera::appdelete(); + return -1; + } +} diff --git a/test_problems/cathermo/HMW_graph_GvT/Makefile.in b/test_problems/cathermo/HMW_graph_GvT/Makefile.in new file mode 100644 index 000000000..83365cf00 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/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_GvT + +# the object files to be linked together. List those generated from Fortran +# and from C/C++ separately +OBJS = HMW_graph_GvT.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 Gex_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) + +Gex_standalone: Gex_standalone.o + $(CXX) -o Gex_standalone Gex_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_GvT/NaCl_Solid.xml b/test_problems/cathermo/HMW_graph_GvT/NaCl_Solid.xml new file mode 100644 index 000000000..d711be8ff --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/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_GvT/README b/test_problems/cathermo/HMW_graph_GvT/README new file mode 100644 index 000000000..4b3f97dd8 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/README @@ -0,0 +1,13 @@ + +Notes on this calculation + +Silverster and Pitzer have the following numbers for equilibrium at 25C + Delta_G0 = -9.0416 kJ gmol-1 + m_sat = 6.146 + meanAC_moll = 1.008 +This table is meant to run through that calculation. at 298.15. +Thus Delta_G(298.15) should be 0.0. It's -0.00487 kJ/gmol -> very close =). + + + + diff --git a/test_problems/cathermo/HMW_graph_GvT/TemperatureTable.h b/test_problems/cathermo/HMW_graph_GvT/TemperatureTable.h new file mode 100644 index 000000000..70cb2bf56 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/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_GvT/output_blessed.txt b/test_problems/cathermo/HMW_graph_GvT/output_blessed.txt new file mode 100644 index 000000000..43baa0e80 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/output_blessed.txt @@ -0,0 +1,42 @@ +Fixed Concentration of the System: + Species Mole_Fraction Molality + Na+ 0.0906484 6.146 + Cl- 0.0906484 6.146 + H2O(L) 0.818703 + + Table at fixed molality(Delta_G refers to rxn, NaCl(s) -> Na+ + Cl-) + -> pressure follows the saturation pressure above one atmosphere) + -> This calculation is meant to test Gibbs_ex -> TODO + + T, Pres, Aphi, Delta_G0, Delta_G, molarGibbs0, molarGibbs, Gibbs_ex, meanAC_moll, OsmCoeff-1, Gibbs_ex_Formula, IdealMixing + Kelvin, bars, sqrt(kg/gmol), kJ/gmol, kJ/gmol, kJ/kgWater, kJ/kgWater, kJ/kgWater, , , kJ/kgWater, kJ/kgWater + 293.15, 1.01325, 0.388171, -8.81994, 0.0281518, -19715.3, -19699.4, -8.55527, 0.999275857, 0.284828391, -8.55527, 24.4417 + 298.15, 1.01325, 0.391447, -9.04161, -0.00455546, -19738.1, -19721.7, -8.48792, 1.00697377, 0.285503623, -8.48793, 24.8586 + 303.15, 1.01325, 0.394888, -9.2523, -0.0401229, -19761.2, -19744.4, -8.47441, 1.01169426, 0.285149994, -8.47441, 25.2755 + 313.15, 1.01325, 0.402262, -9.64279, -0.119784, -19808.3, -19790.8, -8.60278, 1.01304296, 0.281758709, -8.60278, 26.1092 + 323.15, 1.01325, 0.410293, -9.99482, -0.210534, -19856.7, -19838.6, -8.92949, 1.00500239, 0.275364003, -8.92949, 26.943 + 333.15, 1.01325, 0.418983, -10.311, -0.311811, -19906.1, -19887.8, -9.44562, 0.989182754, 0.266541029, -9.44562, 27.7768 + 343.15, 1.01325, 0.428341, -10.593, -0.422831, -19956.8, -19938.3, -10.1442, 0.967055992, 0.255753666, -10.1442, 28.6105 + 353.15, 1.01325, 0.438377, -10.8421, -0.542633, -20008.5, -19990.1, -11.0202, 0.939922258, 0.243373279, -11.0202, 29.4443 + 363.15, 1.01325, 0.44911, -11.0585, -0.670113, -20061.3, -20043.1, -12.0702, 0.908900029, 0.229693982, -12.0702, 30.278 + 373.15, 1.01418, 0.460559, -11.2422, -0.80406, -20115.2, -20097.4, -13.2927, 0.874931349, 0.214945074, -13.2927, 31.1118 + 383.15, 1.43379, 0.472737, -11.3928, -0.943349, -20170, -20152.8, -14.6866, 0.838832235, 0.199310422, -14.6866, 31.9456 + 393.15, 1.98674, 0.485682, -11.5087, -1.08648, -20225.8, -20209.3, -16.2545, 0.801217017, 0.182913805, -16.2545, 32.7793 + 403.15, 2.7028, 0.499432, -11.5883, -1.23212, -20282.5, -20266.9, -18, 0.762601332, 0.165845486, -18, 33.6131 + 413.15, 3.61539, 0.514026, -11.6294, -1.37896, -20340.1, -20325.6, -19.9283, 0.723396901, 0.148161838, -19.9283, 34.4468 + 423.15, 4.76165, 0.529514, -11.6295, -1.52576, -20398.5, -20385.3, -22.0467, 0.683929685, 0.129889688, -22.0467, 35.2806 + 433.15, 6.18235, 0.545955, -11.5858, -1.67144, -20457.7, -20446, -24.3645, 0.64445611, 0.111029683, -24.3645, 36.1144 + 443.15, 7.92187, 0.563414, -11.4949, -1.8151, -20517.7, -20507.6, -26.8932, 0.605177169, 0.0915587819, -26.8932, 36.9481 + 453.15, 10.0281, 0.581971, -11.3534, -1.95611, -20578.3, -20570.2, -29.6471, 0.566250414, 0.0714319193, -29.6471, 37.7819 + 463.15, 12.5524, 0.601719, -11.1573, -2.09422, -20639.6, -20633.6, -32.6429, 0.527799977, 0.0505828468, -32.6429, 38.6157 + 473.15, 15.5493, 0.622769, -10.9025, -2.22964, -20701.5, -20697.9, -35.9013, 0.489924798, 0.0289240925, -35.9013, 39.4494 + 483.15, 19.0767, 0.645254, -10.5846, -2.36325, -20763.9, -20763.1, -39.4466, 0.452705284, 0.00634592246, -39.4466, 40.2832 + 493.15, 23.1959, 0.669336, -10.199, -2.49675, -20826.8, -20829, -43.3084, 0.416208613, -0.0172859171, -43.3084, 41.1169 + 503.15, 27.9709, 0.695212, -9.74073, -2.63297, -20890.2, -20895.8, -47.5225, 0.380492894, -0.0421340312, -47.5225, 41.9507 + 513.15, 33.4693, 0.723127, -9.20474, -2.77622, -20954, -20963.4, -52.1325, 0.345610386, -0.0683953538, -52.1325, 42.7845 + 523.15, 39.7617, 0.753389, -8.58571, -2.93283, -21018.1, -21031.7, -57.1929, 0.311609996, -0.0963101812, -57.1929, 43.6182 + 533.15, 46.9226, 0.786394, -7.87817, -3.11191, -21082.6, -21100.9, -62.7721, 0.278539266, -0.126175738, -62.7721, 44.452 + 543.15, 55.0299, 0.822655, -7.07643, -3.32643, -21147.3, -21170.9, -68.9579, 0.24644614, -0.158366728, -68.9579, 45.2857 + 553.15, 64.1658, 0.862861, -6.17465, -3.59495, -21212.1, -21241.9, -75.8658, 0.215380867, -0.193367067, -75.8658, 46.1195 + 563.15, 74.4178, 0.907952, -5.16684, -3.94426, -21277.1, -21313.8, -83.6516, 0.1853986, -0.231820303, -83.6516, 46.9533 + 573.15, 85.879, 0.959258, -4.04685, -4.41371, -21342.2, -21386.9, -92.5326, 0.156563632, -0.27461274, -92.5326, 47.787 diff --git a/test_problems/cathermo/HMW_graph_GvT/runtest b/test_problems/cathermo/HMW_graph_GvT/runtest new file mode 100755 index 000000000..e08a57210 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/runtest @@ -0,0 +1,42 @@ +#!/bin/sh +# +# +temp_success="1" +/bin/rm -f output.txt outputa.txt + +########################################################################## +prog=HMW_graph_GvT +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_GvT/sortAlgorithms.cpp b/test_problems/cathermo/HMW_graph_GvT/sortAlgorithms.cpp new file mode 100644 index 000000000..d97e51b40 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/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_GvT/sortAlgorithms.h b/test_problems/cathermo/HMW_graph_GvT/sortAlgorithms.h new file mode 100644 index 000000000..72a7fc2a2 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_GvT/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