From c2b8b1d7c8a9ee66f0d32f151804072228edb118 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sun, 20 Jan 2008 20:29:48 +0000 Subject: [PATCH] Split functions in a file into 2 files, in order to clarify roles. --- Cantera/src/equil/Makefile.in | 2 +- Cantera/src/equil/equilibrate.cpp | 4 +- Cantera/src/equil/vcs_MultiPhaseEquil.cpp | 168 +--------- Cantera/src/equil/vcs_MultiPhaseEquil.h | 16 +- Cantera/src/equil/vcs_equilibrate.cpp | 357 ++++++++++++++++++++++ Cantera/src/equil/vcs_solve.h | 2 - 6 files changed, 375 insertions(+), 174 deletions(-) create mode 100644 Cantera/src/equil/vcs_equilibrate.cpp diff --git a/Cantera/src/equil/Makefile.in b/Cantera/src/equil/Makefile.in index c13c582b7..e530375c9 100644 --- a/Cantera/src/equil/Makefile.in +++ b/Cantera/src/equil/Makefile.in @@ -69,7 +69,7 @@ VCSNONIDEAL_OBJ = vcs_solve_TP.o vcs_VolPhase.o vcs_solve.o vcs_prob.o \ vcs_nasa_poly.o vcs_nondim.o vcs_Exception.o \ vcs_funcVtot.o vcs_inest.o vcs_rearrange.o \ vcs_root1d.o vcs_rxnadj.o vcs_timer_generic.o \ - vcs_SpeciesProperties.o \ + vcs_SpeciesProperties.o vcs_equilibrate.o \ vcs_prep.o vcs_species_thermo.o vcs_Gibbs.o \ $(DALT_OBJ) diff --git a/Cantera/src/equil/equilibrate.cpp b/Cantera/src/equil/equilibrate.cpp index a4ca22128..65e2307f8 100644 --- a/Cantera/src/equil/equilibrate.cpp +++ b/Cantera/src/equil/equilibrate.cpp @@ -133,8 +133,8 @@ namespace Cantera { m->addPhase(&s, 1.0); m->init(); nAttempts++; - (void) vcs_equilibrate(*m, XY, estimateEquil, printLvlSub, - rtol, maxsteps, maxiter, loglevel-1); + vcs_equilibrate(*m, XY, estimateEquil, printLvlSub, solver, + rtol, maxsteps, maxiter, loglevel-1); redo = false; if (loglevel > 0) addLogEntry("VCSnonideal solver succeeded."); diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index dd0a17095..037f905ce 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -1,9 +1,8 @@ /** * @file vcs_MultiPhaseEquil.cpp - * Driver routines for equilibrium solvers + * Driver routine for the VCSnonideal equilibrium solver package */ /* - * * $Id$ */ /* @@ -11,6 +10,7 @@ * Contract DE-AC04-94AL85000 with Sandia Corporation, the * U.S. Government retains certain rights in this software. */ + #include "vcs_MultiPhaseEquil.h" #include "vcs_prob.h" #include "vcs_internal.h" @@ -1618,171 +1618,7 @@ namespace Cantera { return VCS_SUCCESS; } - /* - * Set a single-phase chemical solution to chemical equilibrium. - * This is a convenience function that uses one or the other of - * the two chemical equilibrium solvers. - * - * @param s The object to set to an equilibrium state - * - * @param XY An integer specifying the two properties to be held - * constant. - * - * @param estimateEquil Boolean indicating whether the solver - * should estimate its own initial condition. - * If false, the initial mole fraction vector - * in the %ThermoPhase object is used as the - * initial condition. - * - * @param printLvl Determines the amount of printing that - * gets sent to stdout from the vcs package - * (Note, you may have to compile with debug - * flags to get some printing). - * - * @param solver The equilibrium solver to use. If solver = 0, - * the ChemEquil solver will be used, and if - * solver = 1, the vcs_MultiPhaseEquil solver will - * be used (slower than ChemEquil, - * but more stable). If solver < 0 (default, then - * ChemEquil will be tried first, and if it fails - * vcs_MultiPhaseEquil will be tried. - * - * @param maxsteps The maximum number of steps to take to find - * the solution. - * - * @param maxiter For the MultiPhaseEquil solver only, this is - * the maximum number of outer temperature or - * pressure iterations to take when T and/or P is - * not held fixed. - * - * @param loglevel Controls amount of diagnostic output. loglevel - * = 0 suppresses diagnostics, and increasingly-verbose - * messages are written as loglevel increases. The - * messages are written to a file in HTML format for viewing - * in a web browser. @see HTML_logs - */ - int vcs_equilibrate(thermo_t& s, const char* XY, - bool estimateEquil, int printLvl, - int solver, - doublereal rtol, int maxsteps, int maxiter, - int loglevel) { - MultiPhase* m = 0; - bool redo = true; - int retn = 1; - beginLogGroup("equilibrate", loglevel); - addLogEntry("Single-phase equilibrate function"); - { - beginLogGroup("arguments"); - addLogEntry("phase",s.id()); - addLogEntry("XY",XY); - addLogEntry("solver",solver); - addLogEntry("rtol",rtol); - addLogEntry("maxsteps",maxsteps); - addLogEntry("maxiter",maxiter); - addLogEntry("loglevel",loglevel); - endLogGroup("arguments"); - } - - if (solver > 0) { - m = new MultiPhase; - try { - /* - * Set the kmoles of the phase to 1.0, arbitrarily. - * It actually doesn't matter. - */ - m->addPhase(&s, 1.0); - m->init(); - retn = vcs_equilibrate(*m, XY, estimateEquil, printLvl, - rtol, maxsteps, maxiter, loglevel); - redo = false; - addLogEntry("MultiPhaseEquil solver succeeded."); - delete m; - } - catch (CanteraError err) { - addLogEntry("MultiPhaseEquil solver failed."); - delete m; - throw err; - } - } else { - throw CanteraError("vcs_equilibrate", - "ChemEquil not implemented in this interface yet"); - } - - /* - * We are here only for a success - */ - endLogGroup("equilibrate"); - return retn; - } - int vcs_equilibrate(MultiPhase& s, const char* XY, - bool estimateEquil, int printLvl, - doublereal tol, int maxsteps, int maxiter, - int loglevel) { - int ixy = _equilflag(XY); - int retn = vcs_equilibrate_1(s, ixy, estimateEquil, printLvl, - tol, maxsteps, maxiter, loglevel); - return retn; - }; - - /* - * Set a multiphase mixture to a state of chemical equilibrium. - * This is the top-level driver for multiphase equilibrium. It - * doesn't do much more than call the equilibrate method of class - * MultiPhase, except that it adds some messages to the logfile, - * if loglevel is set > 0. - * - * @return Returns the value of 1 if successful. - * - * @ingroup equil - */ - int vcs_equilibrate_1(MultiPhase& s, int ixy, - bool estimateEquil, int printLvl, - doublereal tol, int maxsteps, int maxiter, int loglevel) { - static int counter = 0; - int retn = 1; - beginLogGroup("equilibrate",loglevel); - addLogEntry("multiphase equilibrate function"); - beginLogGroup("arguments"); - addLogEntry("XY",ixy); - addLogEntry("tol",tol); - addLogEntry("maxsteps",maxsteps); - addLogEntry("maxiter",maxiter); - addLogEntry("loglevel",loglevel); - endLogGroup("arguments"); - int printLvlSub = MAX(0, printLvl-1); - - s.init(); - - try { - vcs_MultiPhaseEquil *eqsolve = new vcs_MultiPhaseEquil(&s, printLvlSub); - int err = eqsolve->equilibrate(ixy, estimateEquil, printLvlSub, - tol, maxsteps, maxiter); - if (err != 0) { - retn = 0; - } - addLogEntry("Success. Error", err); - endLogGroup("equilibrate"); - // hard code a csv output file. - if (printLvl > 0) { - string reportFile = "vcs_equilibrate_res.csv"; - if (counter > 0) { - reportFile = "vcs_equilibrate_res_" + int2str(counter) + ".csv"; - } - eqsolve->reportCSV(reportFile); - counter++; - } - delete eqsolve; - } - catch (CanteraError e) { - addLogEntry("Failure.", lastErrorMessage()); - endLogGroup("equilibrate"); - throw e; - } - - return retn; - } - } diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.h b/Cantera/src/equil/vcs_MultiPhaseEquil.h index 5e17528c0..52ececa06 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.h +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.h @@ -111,6 +111,10 @@ namespace Cantera { * (Note, you may have to compile with debug * flags to get some printing). * + * @param solver Determines which solver is used. + * - 1 MultiPhaseEquil solver + * - 2 VCSnonideal Solver (default) + * * @param maxsteps The maximum number of steps to take to find * the solution. * @@ -129,6 +133,7 @@ namespace Cantera { */ int vcs_equilibrate(MultiPhase& s, const char* XY, bool estimateEquil = false, int printLvl = 0, + int solver = 2, doublereal rtol = 1.0e-9, int maxsteps = 1000, int maxiter = 100, int loglevel = -99); @@ -157,6 +162,10 @@ namespace Cantera { * (Note, you may have to compile with debug * flags to get some printing). * + * @param solver Determines which solver is used. + * - 1 MultiPhaseEquil solver + * - 2 VCSnonideal Solver (default) + * * @param maxsteps The maximum number of steps to take to find * the solution. * @@ -174,9 +183,10 @@ namespace Cantera { * @ingroup equilfunctions */ int vcs_equilibrate_1(MultiPhase& s, int ixy, - bool estimateEquil = false, int printLvl = 0, - doublereal rtol = 1.0e-9, int maxsteps = 1000, - int maxiter = 100, int loglevel = -99); + bool estimateEquil = false, int printLvl = 0, + int solver = 2, + doublereal rtol = 1.0e-9, int maxsteps = 1000, + int maxiter = 100, int loglevel = -99); //! Cantera's Interface to the Multiphase chemical equilibrium solver. /*! diff --git a/Cantera/src/equil/vcs_equilibrate.cpp b/Cantera/src/equil/vcs_equilibrate.cpp new file mode 100644 index 000000000..b36563549 --- /dev/null +++ b/Cantera/src/equil/vcs_equilibrate.cpp @@ -0,0 +1,357 @@ +/** + * @file vcs_equilibrate.cpp + * Driver routines for equilibrium solvers + */ +/* + * + * $Id$ + */ +/* + * Copywrite (2006) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ +#include "vcs_MultiPhaseEquil.h" +#include "vcs_prob.h" +#include "vcs_internal.h" +#include "vcs_VolPhase.h" +#include "vcs_species_thermo.h" +#include "vcs_SpeciesProperties.h" +#include "vcs_VolPhase.h" +#include "vcs_nasa_poly.h" +#include "vcs_solve.h" +#include "equil.h" + +#include "ct_defs.h" +#include "mix_defs.h" +#include "speciesThermoTypes.h" +#ifdef WITH_IDEAL_SOLUTIONS +#include "IdealSolidSolnPhase.h" +#endif +#ifdef WITH_ELECTROLYTES +#include "IdealMolalSoln.h" +#endif +#include "ChemEquil.h" + +#include +#include + +using namespace Cantera; +using namespace std; + +namespace Cantera { + + /* + * Set a single-phase chemical solution to chemical equilibrium. + * This is a convenience function that uses one or the other of + * the two chemical equilibrium solvers. + * + * @param s The object to set to an equilibrium state + * + * @param XY An integer specifying the two properties to be held + * constant. + * + * @param estimateEquil Boolean indicating whether the solver + * should estimate its own initial condition. + * If false, the initial mole fraction vector + * in the %ThermoPhase object is used as the + * initial condition. + * + * @param printLvl Determines the amount of printing that + * gets sent to stdout from the vcs package + * (Note, you may have to compile with debug + * flags to get some printing). + * + * @param solver The equilibrium solver to use. If solver = 0, + * the ChemEquil solver will be used, and if + * solver = 1, the vcs_MultiPhaseEquil solver will + * be used (slower than ChemEquil, + * but more stable). If solver < 0 (default, then + * ChemEquil will be tried first, and if it fails + * vcs_MultiPhaseEquil will be tried. + * + * @param maxsteps The maximum number of steps to take to find + * the solution. + * + * @param maxiter For the MultiPhaseEquil solver only, this is + * the maximum number of outer temperature or + * pressure iterations to take when T and/or P is + * not held fixed. + * + * @param loglevel Controls amount of diagnostic output. loglevel + * = 0 suppresses diagnostics, and increasingly-verbose + * messages are written as loglevel increases. The + * messages are written to a file in HTML format for viewing + * in a web browser. @see HTML_logs + */ + int vcs_equilibrate(thermo_t& s, const char* XY, + bool estimateEquil, int printLvl, + int solver, + doublereal rtol, int maxsteps, int maxiter, + int loglevel) { + MultiPhase* m = 0; + int retn = 1; + int retnSub = 0; + + beginLogGroup("equilibrate", loglevel); + // retry: + addLogEntry("Single-phase equilibrate function"); + { + beginLogGroup("arguments"); + addLogEntry("phase",s.id()); + addLogEntry("XY",XY); + addLogEntry("solver",solver); + addLogEntry("rtol",rtol); + addLogEntry("maxsteps",maxsteps); + addLogEntry("maxiter",maxiter); + addLogEntry("loglevel",loglevel); + endLogGroup("arguments"); + } + + if (solver == 2) { + m = new MultiPhase; + try { + /* + * Set the kmoles of the phase to 1.0, arbitrarily. + * It actually doesn't matter. + */ + m->addPhase(&s, 1.0); + m->init(); + + retn = vcs_equilibrate(*m, XY, estimateEquil, printLvl, solver, + rtol, maxsteps, maxiter, loglevel); + addLogEntry("MultiPhaseEquil solver succeeded."); + delete m; + } + catch (CanteraError &err) { + addLogEntry("MultiPhaseEquil solver failed."); + delete m; + throw err; + } + } else if (solver == 1) { + m = new MultiPhase; + try { + m->addPhase(&s, 1.0); + m->init(); + (void) equilibrate(*m, XY, rtol, maxsteps, maxiter, loglevel-1); + if (loglevel > 0) + addLogEntry("MultiPhaseEquil solver succeeded."); + delete m; + retn = 1; + } + catch (CanteraError &err) { + if (loglevel > 0) { + addLogEntry("MultiPhaseEquil solver failed."); + } + delete m; + throw err; + } + } else if (solver == 0) { + ChemEquil *e = new ChemEquil; + try { + e->options.maxIterations = maxsteps; + e->options.relTolerance = rtol; + retnSub = e->equilibrate(s,XY,loglevel-1); + if (retnSub < 0) { + if (loglevel > 0) { + addLogEntry("ChemEquil solver failed."); + } + delete e; + throw CanteraError("equilibrate", + "ChemEquil equilibrium solver failed"); + } + retn = 1; + s.setElementPotentials(e->elementPotentials()); + delete e; + if (loglevel > 0) { + addLogEntry("ChemEquil solver succeeded."); + } + } + catch (CanteraError &err) { + if (loglevel > 0) { + addLogEntry("ChemEquil solver failed."); + } + delete e; + throw err; + } + } else { + throw CanteraError("vcs_equilibrate", + "unknown solver"); + } + + /* + * We are here only for a success + */ + endLogGroup("equilibrate"); + return retn; + } + + // Set a multi-phase chemical solution to chemical equilibrium. + /* + * This function uses the vcs_MultiPhaseEquil interface to the + * vcs solver. + * The function uses the element abundance vector that is + * currently consistent with the composition within the phases + * themselves. Two other thermodynamic quantities, determined by the + * XY string, are held constant during the equilibration. + * + * @param s The object to set to an equilibrium state + * + * @param XY A character string specifying the two properties to + * be held constant + * + * @param estimateEquil Boolean indicating whether the solver + * should estimate its own initial condition. + * If false, the initial mole fraction vector + * in the %ThermoPhase object is used as the + * initial condition. + * + * @param printLvl Determines the amount of printing that + * gets sent to stdout from the vcs package + * (Note, you may have to compile with debug + * flags to get some printing). + * + * @param maxsteps The maximum number of steps to take to find + * the solution. + * + * @param maxiter For the MultiPhaseEquil solver only, this is + * the maximum number of outer temperature or + * pressure iterations to take when T and/or P is + * not held fixed. + * + * @param loglevel Controls amount of diagnostic output. loglevel + * = 0 suppresses diagnostics, and increasingly-verbose + * messages are written as loglevel increases. The + * messages are written to a file in HTML format for viewing + * in a web browser. @see HTML_logs + * + * @ingroup equilfunctions + */ + int vcs_equilibrate(MultiPhase& s, const char* XY, + bool estimateEquil, int printLvl, int solver, + doublereal tol, int maxsteps, int maxiter, + int loglevel) { + int ixy = _equilflag(XY); + int retn = vcs_equilibrate_1(s, ixy, estimateEquil, printLvl, solver, + tol, maxsteps, maxiter, loglevel); + return retn; + }; + + + // Set a multi-phase chemical solution to chemical equilibrium. + /* + * This function uses the vcs_MultiPhaseEquil interface to the + * vcs solver. + * The function uses the element abundance vector that is + * currently consistent with the composition within the phases + * themselves. Two other thermodynamic quantities, determined by the + * XY string, are held constant during the equilibration. + * + * @param s The object to set to an equilibrium state + * + * @param XY An integer specifying the two properties to be held + * constant. + * + * @param estimateEquil Boolean indicating whether the solver + * should estimate its own initial condition. + * If false, the initial mole fraction vector + * in the %ThermoPhase object is used as the + * initial condition. + * + * @param printLvl Determines the amount of printing that + * gets sent to stdout from the vcs package + * (Note, you may have to compile with debug + * flags to get some printing). + * + * @param maxsteps The maximum number of steps to take to find + * the solution. + * + * @param maxiter For the MultiPhaseEquil solver only, this is + * the maximum number of outer temperature or + * pressure iterations to take when T and/or P is + * not held fixed. + * + * @param loglevel Controls amount of diagnostic output. loglevel + * = 0 suppresses diagnostics, and increasingly-verbose + * messages are written as loglevel increases. The + * messages are written to a file in HTML format for viewing + * in a web browser. @see HTML_logs + * + * @ingroup equilfunctions + */ + int vcs_equilibrate_1(MultiPhase& s, int ixy, + bool estimateEquil, int printLvl, int solver, + doublereal tol, int maxsteps, int maxiter, int loglevel) { + static int counter = 0; + int retn = 1; + beginLogGroup("equilibrate",loglevel); + addLogEntry("multiphase equilibrate function"); + beginLogGroup("arguments"); + addLogEntry("XY",ixy); + addLogEntry("tol",tol); + addLogEntry("maxsteps",maxsteps); + addLogEntry("maxiter",maxiter); + addLogEntry("loglevel",loglevel); + endLogGroup("arguments"); + int printLvlSub = MAX(0, printLvl-1); + + s.init(); + + if (solver == 2) { + try { + vcs_MultiPhaseEquil *eqsolve = new vcs_MultiPhaseEquil(&s, printLvlSub); + int err = eqsolve->equilibrate(ixy, estimateEquil, printLvlSub, + tol, maxsteps, maxiter); + if (err != 0) { + retn = 0; + } + addLogEntry("Success. Error", err); + endLogGroup("equilibrate"); + // hard code a csv output file. + if (printLvl > 0) { + string reportFile = "vcs_equilibrate_res.csv"; + if (counter > 0) { + reportFile = "vcs_equilibrate_res_" + int2str(counter) + ".csv"; + } + eqsolve->reportCSV(reportFile); + counter++; + } + delete eqsolve; + } + catch (CanteraError &e) { + addLogEntry("Failure.", lastErrorMessage()); + endLogGroup("equilibrate"); + throw e; + } + } else if (solver == 1) { + if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) { + try { + double err = s.equilibrate(ixy, tol, maxsteps, maxiter, loglevel); + if (loglevel > 0) { + addLogEntry("Success. Error",err); + endLogGroup("equilibrate"); + } + return 0; + } + catch (CanteraError &e) { + if (loglevel > 0) { + addLogEntry("Failure.",lastErrorMessage()); + endLogGroup("equilibrate"); + } + throw e; + } + } + else { + if (loglevel > 0) { + addLogEntry("multiphase equilibrium can be done only for TP, HP, SP, or TV"); + endLogGroup("equilibrate"); + } + throw CanteraError("equilibrate","unsupported option"); + //return -1.0; + } + } else { + throw CanteraError("vcs_equilibrate_1", "unknown solver"); + } + return retn; + } +} diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index f139a96ab..47b57bfa4 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -779,9 +779,7 @@ public: VCS_COUNTERS *m_VCount; -#ifdef DEBUG_MODE int vcs_debug_print_lvl; -#endif //! Units for the chemical potential data: /*!