From ba5928083f6bc15c09258a5b65cb8783366e5671 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Thu, 10 May 2007 03:28:30 +0000 Subject: [PATCH] added thread safety option --- Cantera/clib/src/Makefile.in | 2 +- Cantera/python/setup.py.in | 10 +- Cantera/src/base/Makefile.in | 8 +- Cantera/src/base/global.h | 120 +- Cantera/src/base/misc.cpp | 740 +++++++--- Cantera/src/base/units.h | 121 +- Cantera/src/equil/BasisOptimize.cpp | 284 ++-- Cantera/src/equil/ChemEquil.cpp | 1268 +++++++++-------- Cantera/src/equil/ChemEquil.h | 12 +- Cantera/src/equil/MultiPhase.cpp | 129 +- Cantera/src/equil/MultiPhaseEquil.cpp | 128 +- Cantera/src/equil/MultiPhaseEquil.h | 16 +- Cantera/src/equil/equilibrate.cpp | 133 +- Cantera/src/kinetics/FalloffFactory.cpp | 4 + Cantera/src/kinetics/FalloffFactory.h | 37 +- Cantera/src/kinetics/FalloffMgr.h | 21 +- Cantera/src/kinetics/KineticsFactory.cpp | 46 +- Cantera/src/kinetics/KineticsFactory.h | 38 +- Cantera/src/kinetics/importKinetics.cpp | 91 +- Cantera/src/thermo/SpeciesThermoFactory.cpp | 118 +- Cantera/src/thermo/SpeciesThermoFactory.h | 16 +- Cantera/src/thermo/ThermoFactory.cpp | 3 + Cantera/src/thermo/ThermoFactory.h | 47 +- Cantera/src/transport/TransportFactory.cpp | 25 +- Cantera/src/transport/TransportFactory.h | 66 +- Cantera/src/zeroD/Makefile.in | 2 +- Cantera/src/zeroD/ReactorFactory.cpp | 4 + Cantera/src/zeroD/ReactorFactory.h | 31 +- config.h.in | 4 +- configure | 51 +- configure.in | 29 + ext/f2c_libs/arith.h | 5 +- preconfig | 25 + test_problems/diamondSurf/Makefile.in | 2 +- .../min_python/minDiamond/Makefile.in | 2 +- test_problems/min_python/negATest/Makefile.in | 2 +- test_problems/negATest/Makefile.in | 2 +- test_problems/pureFluidTest/Makefile.in | 2 +- test_problems/silane_equil/Makefile.in | 2 +- test_problems/surfkin/Makefile.in | 2 +- 40 files changed, 2138 insertions(+), 1510 deletions(-) diff --git a/Cantera/clib/src/Makefile.in b/Cantera/clib/src/Makefile.in index 089e2fd4e..20325c92d 100755 --- a/Cantera/clib/src/Makefile.in +++ b/Cantera/clib/src/Makefile.in @@ -50,7 +50,7 @@ LIB_DEPS = $(CANTERA_LIBDIR)/libctbase.a \ $(CANTERA_LIBDIR)/liboneD.a # the directory where Cantera include files may be found. -CXX_INCLUDES = -I../../src/base -I../../src/thermo -I../../src/kinetics -I../../src/transport -I../../src/numerics -I../../src/oneD -I../../src/zeroD -I../../src/equil -I../../src/converters +CXX_INCLUDES = -I../../src/base -I../../src/thermo -I../../src/kinetics -I../../src/transport -I../../src/numerics -I../../src/oneD -I../../src/zeroD -I../../src/equil -I../../src/converters @CXX_INCLUDES@ # flags passed to the C++ compiler/linker for the linking step LCXX_FLAGS = -L$(CANTERA_LIBDIR) @CXXFLAGS@ diff --git a/Cantera/python/setup.py.in b/Cantera/python/setup.py.in index 94792edb4..3371ffa59 100644 --- a/Cantera/python/setup.py.in +++ b/Cantera/python/setup.py.in @@ -45,13 +45,17 @@ endlibstr1 = "@LCXX_END_LIBS@" endlib1 = endlibstr1.replace('-l', ' ') endlib = endlib1.split() +locallibstr = "@LOCAL_LIBS@" +local1 = locallibstr.replace('-l', ' ') +locallibs = local1.split() + if platform == "win32": libs = ["clib", "zeroD","oneD","kinetics", "transport", "cantera"] + bllist + cvlist + ["ctbase", "ctmath", "tpx"] else: - - libs = ["clib", "zeroD","oneD", "kinetics", "transport", - "equil", "ctnumerics", "converters"] + thermolib + bllist + cvlist + ["ctbase", "ctmath", "tpx"] + libs = ["clib"] + locallibs + #libs = ["clib", "zeroD","oneD", "kinetics", "transport", + # "equil", "ctnumerics", "converters"] + thermolib + bllist + cvlist + ["ctbase", "ctmath", "tpx", "boost_thread-mt-1_34"] if @build_with_f2c@ == 1: libs.append("ctf2c") diff --git a/Cantera/src/base/Makefile.in b/Cantera/src/base/Makefile.in index 29faf35ea..f849666bc 100644 --- a/Cantera/src/base/Makefile.in +++ b/Cantera/src/base/Makefile.in @@ -26,13 +26,13 @@ PIC_FLAG=@PIC@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG) -BASE_OBJ = misc.o ct2ctml.o ctml.o plots.o stringUtils.o xml.o +BASE_OBJ = ct2ctml.o ctml.o misc.o plots.o stringUtils.o xml.o -BASE_H = ct_defs.h ctexceptions.h global.h logger.h XML_Writer.h \ +BASE_H = ct_defs.h ctexceptions.h logger.h XML_Writer.h \ ctml.h plots.h stringUtils.h xml.h config.h utilities.h \ - Array.h vec_functions.h + Array.h vec_functions.h global.h FactoryBase.h -CXX_INCLUDES = -I. +CXX_INCLUDES = -I. @CXX_INCLUDES@ LIB = @buildlib@/libctbase.a DEPENDS = $(BASE_OBJ:.o=.d) diff --git a/Cantera/src/base/global.h b/Cantera/src/base/global.h index f1c695769..e631ac8a0 100755 --- a/Cantera/src/base/global.h +++ b/Cantera/src/base/global.h @@ -191,6 +191,14 @@ namespace Cantera { */ void appdelete(); + //! Delete and free memory allocated per thread in multithreaded applications + /*! + * Delete the memory allocated per thread by Cantera. It should be called from + * within the thread just before the thread terminates. If your version of Cantera has not + * been specifically compiled for thread safety this function does nothing. + */ + void thread_complete() ; + //! Returns root directory where %Cantera where installed /*! * @return @@ -279,6 +287,16 @@ namespace Cantera { * @ingroup textlogs */ void writelog(const char* msg); + + //! Write a formated message to the screen + /*! + * Using the printf formatting of C write a message to the screen + * with variable values. + * + * @param fmt c format string for the following arguments + * @ingroup textlogs + */ + void writelogf(const char* fmt,...); //! Write an error message and terminate execution. /*! @@ -452,104 +470,12 @@ namespace Cantera { inline void write_logfile(std::string file = "log.html") {} #endif + XML_Node* get_XML_Node(const std::string& file_ID, XML_Node* root); - //! Search for an XML_Node either wiithin an existing XML tree structure, or in another file, - //! based on the file name or the XML id attribute. - /*! - * This routine will locate an XML node in either the input - * XML tree or in another input file specified by the file - * part of the file_ID string. Searches are based on the - * ID attribute of the XML element only. - * - * @param file_ID This is a concatenation of two strings seperated - * by the "#" character. The string before the - * pound character is the file name of an xml - * file to carry out the search. The string after - * the # character is the ID attribute - * of the xml element to search for. - * The string is interpreted as a file string if - * no # character is in the string. - * - * @param root If the file string is empty, searches for the - * xml element with matching ID attribute are - * carried out from this XML node. - * - * @return - * This routine will process the XML file, creating an XML - * tree structure. It returns a pointer to the top of the tree. - * - * - * For example, - * @code - * - * XML_Node* xn = get_XML_Node("phase", "gri30.xml#gri30_mix", 0); - * - * @endcode - * - * will search in the file gri30.xml for an XML element of the following form, where - * the XML element name, phase, is an optional hit: - * @verbatim - +#include #ifdef WIN32 #include #include #endif using namespace std; +// RFB : If running multiple threads in a cpp application, +// the Application class is the only internal object +// that is single instance with static data. Synchronize access to those data +// structures. +// Using macros to avoid polluting code with alot of ifdef's + +#ifdef THREAD_SAFE_CANTERA + #include + #include + static boost::mutex dir_mutex; // For input directory access + static boost::mutex msg_mutex; // For access to string messages + static boost::mutex app_mutex; // Application state including creating singleton + //static boost::mutex log_mutex; // Logger pointer + static boost::mutex xml_mutex; // XML file storage + + #define DIR_LOCK() boost::mutex::scoped_lock d_lock(dir_mutex) + #define MSG_LOCK() boost::mutex::scoped_lock m_lock(msg_mutex) + #define APP_LOCK() boost::mutex::scoped_lock a_lock(app_mutex) + //#define LOG_LOCK() boost::mutex::scoped_lock l_lock(log_mutex) + #define XML_LOCK() boost::mutex::scoped_lock x_lock(xml_mutex) + + #ifdef WITH_HTML_LOGS + //static boost::mutex html_mutex; // html logs + //#define HTML_LOCK() boost::mutex::scoped_lock h_lock(html_mutex) + #endif + + #if defined(BOOST_HAS_WINTHREADS) + #include + typedef unsigned int cthreadId_t ; + class thread_equal + { + public: + bool operator()(cthreadId_t L, cthreadId_t R) + { + return L == R ; + } + } ; + cthreadId_t getThisThreadId() + { + return ::GetCurrentThreadId() ; + } + #elif defined(BOOST_HAS_PTHREADS) + + typedef pthread_t cthreadId_t ; + class thread_equal + { + public: + bool operator()(cthreadId_t L, cthreadId_t R) + { + return pthread_equal(L, R) ; + } + } ; + cthreadId_t getThisThreadId() + { + return pthread_self() ; + } + #elif defined(BOOST_HAS_MPTASKS) + typedef MPTaskID cthreadId_t ; + class thread_equal + { + public: + bool operator()(cthreadId_t L, cthreadId_t R) + { + return L == R ; + } + } ; + cthreadId_t getThisThreadId() + { + return MPCurrentTaskID() ; + } + #endif + +#else + #define DIR_LOCK() + #define MSG_LOCK() + #define APP_LOCK() + //#define LOG_LOCK() + #define XML_LOCK() + + #ifdef WITH_HTML_LOGS + //#define HTML_LOCK() + #endif + +#endif + namespace Cantera { /*! @@ -71,32 +153,172 @@ namespace Cantera { * @ingroup globalData */ class Application { + + protected: + class Messages + { public: + Messages() + { + // install a default logwriter that writes to standard + // output / standard error + logwriter = new Logger(); + #ifdef WITH_HTML_LOGS + // HTML log files + xmllog = 0; + current = 0; + loglevel = 0; + #endif + } + ~Messages() + { + delete logwriter; + #ifdef WITH_HTML_LOGS + if (xmllog) { + write_logfile("orphan"); + //delete xmllog; + } + #endif + } + + void addError(std::string r, std::string msg) ; + int getErrorCount() ; + void popError() ; + std::string lastErrorMessage() ; + void getErrors( std::ostream& f) ; + void logErrors() ; + + + void writelog(const std::string& msg) ; + void writelog(const char* pszmsg) ; + void logerror(const std::string& msg) ; + int getUserEnv() ; + void setLogger(Logger* logwriter) ; + + #ifdef WITH_HTML_LOGS + void beginLogGroup(std::string title, int loglevel) ; + void addLogEntry(std::string tag, std::string value) ; + void addLogEntry(std::string tag, doublereal value) ; + void addLogEntry(std::string tag, int value) ; + void addLogEntry(std::string msg) ; + void endLogGroup(std::string title) ; + void write_logfile(std::string file) ; + #endif + + protected: + //! Current list of error messages + vector errorMessage; + //! Current list of warning messages + //vector warning; + //! Current error Routine + vector errorRoutine; + //! Last error message + //string msglog; + //! Current line length + // size_t linelen; + //! Current pointer to the logwriter + Logger* logwriter; + #ifdef WITH_HTML_LOGS + //! Current pointer to the top of the XML_Node tree for the current HTML log + XML_Node *xmllog; + //! Pointer to the last current position in the XML_Node tree for the current HTML log + XML_Node *current; + //! Current value of loglevel + int loglevel; + //! Vector of loglevels for loggroups that are open + vector loglevels; + //! Current vector of loggroups that are open + vector loggroups; + #endif + } ; + + #ifdef THREAD_SAFE_CANTERA + typedef boost::shared_ptr< Messages > pMessages_t ; + typedef std::map< cthreadId_t, pMessages_t > threadMsgMap_t ; + class ThreadMessages + { + public: + ThreadMessages() + { + } + + Messages* operator->() + { + MSG_LOCK() ; + cthreadId_t curId = getThisThreadId() ; + threadMsgMap_t::iterator iter = m_threadMsgMap.find( curId ) ; + if ( iter != m_threadMsgMap.end() ) + { + return (iter->second.get()) ; + } + pMessages_t pMsgs( new Messages() ) ; + m_threadMsgMap.insert( std::pair< cthreadId_t, pMessages_t >( curId, pMsgs ) ) ; + return pMsgs.get() ; + } + + void removeThreadMessages() { + MSG_LOCK() ; + cthreadId_t curId = getThisThreadId() ; + threadMsgMap_t::iterator iter = m_threadMsgMap.find( curId ) ; + if ( iter != m_threadMsgMap.end() ) { + m_threadMsgMap.erase( iter ) ; + } + } + + private: + threadMsgMap_t m_threadMsgMap ; + } ; + #endif + + + protected: //RFB Protected ctor access thru static member function Instance //! Constructor for class sets up the initial conditions - Application() : linelen(0), stop_on_error(false), - tmp_dir("."), sleep("1") + Application() : /*linelen(0),*/ stop_on_error(false), + tmp_dir("."), m_sleep("1") +#if !defined( THREAD_SAFE_CANTERA ) + , pMessenger( new Messages() ) +#endif { // if TMP or TEMP is set, use it for the temporary // directory char* tmpdir = getenv("TMP"); if (tmpdir == 0) - tmpdir = getenv("TEMP"); + tmpdir = getenv("TEMP"); if (tmpdir != 0) - tmp_dir = string(tmpdir); + tmp_dir = string(tmpdir); // if SLEEP is set, use it as the sleep time char* sleepstr = getenv("SLEEP"); if (sleepstr != 0) { - sleep = string(sleepstr); + m_sleep = string(sleepstr); } // install a default logwriter that writes to standard // output / standard error - logwriter = new Logger(); - // HTML log files - xmllog = 0; - current = 0; - loglevel = 0; +// logwriter = new Logger(); +//#ifdef WITH_HTML_LOGS +// // HTML log files +// xmllog = 0; +// current = 0; +// loglevel = 0; +//#endif + setDefaultDirectories(); +#if defined(THREAD_SAFE_CANTERA) + Unit::units() ; +#endif + } + + public: + /// Return a pointer to the one and only instance of class Application + /* + * If the an Application object has not yet been created it is created + */ + static Application* Instance() { + APP_LOCK(); + if ( Application::s_app == 0 ) { + Application::s_app = new Application() ; + } + return s_app ; } //! Destructor for class deletes global data @@ -107,28 +329,74 @@ namespace Cantera { virtual ~Application() { map::iterator pos; for (pos = xmlfiles.begin(); pos != xmlfiles.end(); ++pos) { - pos->second->unlock(); - delete pos->second; - pos->second = 0; + pos->second->unlock(); + delete pos->second; + pos->second = 0; } - delete logwriter; - if (xmllog) { - write_logfile("orphan"); - //delete xmllog; +// delete logwriter; +//#ifdef WITH_HTML_LOGS +// if (xmllog) { +// write_logfile("orphan"); +// //delete xmllog; +// } +//#endif + } + + static void ApplicationDestroy() { + APP_LOCK() ; + if ( Application::s_app != 0 ) { + delete Application::s_app ; + Application::s_app = 0 ; } } + + void setTmpDir(std::string tmp) ; + std::string getTmpDir() ; + std::string sleep() ; + void addError(std::string r, std::string msg) { pMessenger->addError(r, msg) ; } + int getErrorCount() { return pMessenger->getErrorCount() ; } + void popError() { pMessenger->popError() ; } + std::string lastErrorMessage() { return pMessenger->lastErrorMessage() ; } + void getErrors( std::ostream& f) { pMessenger->getErrors(f) ; } + void logErrors() { pMessenger->logErrors() ; } + void addDataDirectory( std::string dir ) ; + std::string findInputFile(std::string name) ; + XML_Node* get_XML_File(std::string file, int debug=0) ; + void close_XML_File(std::string file) ; + + void writelog(const std::string& msg) { pMessenger->writelog(msg); } + void writelog(const char* pszmsg) { pMessenger->writelog(pszmsg); } + void logerror(const std::string& msg) {pMessenger->logerror(msg); } + int getUserEnv() { return pMessenger->getUserEnv() ; } + void setLogger(Logger* logwriter) {pMessenger->setLogger(logwriter);} + + void thread_complete() ; + +#ifdef WITH_HTML_LOGS + void beginLogGroup(std::string title, int loglevel) { pMessenger->beginLogGroup(title,loglevel);} + void addLogEntry(std::string tag, std::string value) { pMessenger->addLogEntry(tag, value);} + void addLogEntry(std::string tag, doublereal value) { pMessenger->addLogEntry(tag, value);} + void addLogEntry(std::string tag, int value) { pMessenger->addLogEntry(tag, value);} + void addLogEntry(std::string msg) { pMessenger->addLogEntry(msg); } + void endLogGroup(std::string title) { pMessenger->endLogGroup(title) ;} + void write_logfile(std::string file) { pMessenger->write_logfile(file) ; } +#endif + +protected: + void setDefaultDirectories(); + //! Current vector of input directories to search for input files vector inputDirs; //! Current list of error messages - vector errorMessage; + //vector errorMessage; //! Current list of warning messages - vector warning; + //vector warning; //! Current error Routine - vector errorRoutine; + //vector errorRoutine; //! Last error message - string msglog; + //string msglog; //! Current line length - size_t linelen; + // size_t linelen; //! Current value of stop_on_error bool stop_on_error; //! Current map of options @@ -138,69 +406,79 @@ namespace Cantera { //! Current vector of xml file trees that have been previously parsed map xmlfiles; //! Current sleep command. - string sleep; + string m_sleep; //! Current pointer to the logwriter - Logger* logwriter; + //Logger* logwriter; + #ifdef WITH_HTML_LOGS //! Current pointer to the top of the XML_Node tree for the current HTML log - XML_Node *xmllog; + //XML_Node *xmllog; //! Pointer to the last current position in the XML_Node tree for the current HTML log - XML_Node *current; + //XML_Node *current; //! Current value of loglevel - int loglevel; + //int loglevel; //! Vector of loglevels for loggroups that are open - vector loglevels; + //vector loglevels; //! Current vector of loggroups that are open - vector loggroups; - }; + //vector loggroups; + #endif +#if defined(THREAD_SAFE_CANTERA) + ThreadMessages pMessenger ; +#else + auto_ptr< Messages > pMessenger ; +#endif + private: + static Application* s_app ; + }; - /// Return a pointer to the one and only instance of class Application - Application* app(); + //Application* app(); + inline Application* app() { + return Application::Instance() ; + } - static void setDefaultDirectories(); /// Pointer to the single Application instance - static Application* s_app = 0; + Application* Application::s_app = 0; /** * Definition of the static member of the Unit class. */ Unit* Unit::s_u = 0; + #if defined(THREAD_SAFE_CANTERA) + boost::mutex Unit::units_mutex ; + #endif - static void appinit() { - if (s_app == 0) { - s_app = new Application; - } - } /* * Delete all global data. It should be called at the end of the * application if leak checking is to be done. */ void appdelete() { - if (s_app) { - delete s_app; - s_app = 0; - } -#ifndef DGG_MOD - SpeciesThermoFactory::deleteFactory(); - ThermoFactory::deleteFactory(); - FalloffFactory::deleteFalloffFactory(); + Application::ApplicationDestroy() ; + FactoryBase::deleteFactories(); Unit::deleteUnit(); + } + + /* + * Delete Messenger object allocated per thread. + */ + void thread_complete() { + app()->thread_complete() ; + } + + void Application::thread_complete() { +#if defined(THREAD_SAFE_CANTERA) + pMessenger.removeThreadMessages() ; #endif } - Application* app() { - if (s_app == 0) { - s_app = new Application; - setDefaultDirectories(); - } - return s_app; + XML_Node* get_XML_File(std::string file, int debug) { + return app()->get_XML_File(file, debug) ; } - - XML_Node* get_XML_File(std::string file, int debug) { + XML_Node* Application::get_XML_File(std::string file, int debug) { + XML_LOCK(); string path = ""; /* try { @@ -235,8 +513,8 @@ namespace Cantera { #endif string ff = path; - if (app()->xmlfiles.find(path) - == app()->xmlfiles.end()) { + if (xmlfiles.find(path) + == xmlfiles.end()) { /* * Check whether or not the file is XML. If not, it will * be first processed with the preprocessor. We determine @@ -273,12 +551,12 @@ namespace Cantera { * already processed this file. If we have, return a pointer to * the processed xml tree. */ - if (app()->xmlfiles.find(ff) != app()->xmlfiles.end()) { + if (xmlfiles.find(ff) != xmlfiles.end()) { #ifdef DEBUG_PATHS cout << "get_XML_File(): File, " << ff << ", was previously read." << " Retrieving the storred xml tree." << endl; #endif - return s_app->xmlfiles[ff]; + return xmlfiles[ff]; } /* * Ok, we didn't find the processed XML tree. Do the conversion @@ -299,7 +577,7 @@ namespace Cantera { if (s) { x->build(s); x->lock(); - s_app->xmlfiles[ff] = x; + xmlfiles[ff] = x; } else { string estring = "cannot open "+ff+" for reading."; @@ -312,30 +590,48 @@ namespace Cantera { * lookup operation in the return statement will return a valid * pointer. */ - return s_app->xmlfiles[ff]; + return xmlfiles[ff]; } void close_XML_File(std::string file) { + app()->close_XML_File(file) ; + } + void Application::close_XML_File(std::string file) { + XML_LOCK(); if (file == "all") { map::iterator - b = app()->xmlfiles.begin(), e = app()->xmlfiles.end(); + b = xmlfiles.begin(), + e = xmlfiles.end(); for(; b != e; ++b) { b->second->unlock(); delete b->second; - s_app->xmlfiles.erase(b->first); + xmlfiles.erase(b->first); } } - else if (app()->xmlfiles.find(file) - != app()->xmlfiles.end()) { - s_app->xmlfiles[file]->unlock(); - delete s_app->xmlfiles[file]; - s_app->xmlfiles.erase(file); + else if (xmlfiles.find(file) != xmlfiles.end()) { + xmlfiles[file]->unlock(); + delete xmlfiles[file]; + xmlfiles.erase(file); } } - void setTmpDir(std::string tmp) { app()->tmp_dir = tmp; } - string tmpDir() { appinit(); return app()->tmp_dir; } - string sleep() { appinit(); return app()->sleep; } + void setTmpDir(std::string tmp) { app()->setTmpDir(tmp); } + void Application::setTmpDir(std::string tmp) { + APP_LOCK(); + tmp_dir = tmp ; + } + + string tmpDir() { return app()->getTmpDir(); } + std::string Application::getTmpDir() { + APP_LOCK(); + return tmp_dir ; + } + + string sleep() { return app()->sleep(); } + std::string Application::sleep() { + APP_LOCK(); + return m_sleep ; + } /* @@ -343,7 +639,10 @@ namespace Cantera { * \ingroup errorhandling */ int nErrors() { - return static_cast(app()->errorMessage.size()); + return app()->getErrorCount(); + } + int Application::Messages::getErrorCount() { + return static_cast(errorMessage.size()) ; } /* @@ -354,10 +653,12 @@ namespace Cantera { * \ingroup errorhandling */ void popError() { - appinit(); - if (nErrors() > 0) { - s_app->errorMessage.pop_back(); - s_app->errorRoutine.pop_back(); + app()->popError() ; + } + void Application::Messages::popError() { + if ( static_cast(errorMessage.size()) > 0) { + errorRoutine.pop_back() ; + errorMessage.pop_back() ; } } @@ -368,18 +669,19 @@ namespace Cantera { * \ingroup errorhandling */ string lastErrorMessage() { - appinit(); - if (nErrors() > 0) { + return app()->lastErrorMessage() ; + } + std::string Application::Messages::lastErrorMessage() { + if ( static_cast(errorMessage.size()) > 0) { string head = "\n\n************************************************\n" " Cantera Error! \n" "************************************************\n\n"; - return head+string("\nProcedure: ")+s_app->errorRoutine.back() - +string("\nError: ")+s_app->errorMessage.back(); - } - else { - return ""; - } + return head+string("\nProcedure: ")+errorRoutine.back() + +string("\nError: ")+errorMessage.back(); + } else { + return ""; + } } /* @@ -392,23 +694,25 @@ namespace Cantera { * \ingroup errorhandling */ void showErrors(std::ostream& f) { - appinit(); - int i = static_cast(s_app->errorMessage.size()); + app()->getErrors(f) ; + } + void Application::Messages::getErrors( std::ostream& f) { + int i = static_cast(errorMessage.size()); if (i == 0) return; f << endl << endl; f << "************************************************" << endl; f << " Cantera Error! " << endl; f << "************************************************" << endl - << endl; + << endl; int j; for (j = 0; j < i; j++) { f << endl; - f << "Procedure: " << s_app->errorRoutine[j] << endl; - f << "Error: " << s_app->errorMessage[j] << endl; + f << "Procedure: " << errorRoutine[j] << endl; + f << "Error: " << errorMessage[j] << endl; } f << endl << endl; - s_app->errorMessage.clear(); - s_app->errorRoutine.clear(); + errorMessage.clear(); + errorRoutine.clear(); } /* @@ -421,8 +725,10 @@ namespace Cantera { * \ingroup errorhandling */ void showErrors() { - appinit(); - int i = static_cast(s_app->errorMessage.size()); + app()->logErrors() ; + } + void Application::Messages::logErrors() { + int i = static_cast(errorMessage.size()); if (i == 0) return; writelog("\n\n"); writelog("************************************************\n"); @@ -431,12 +737,12 @@ namespace Cantera { int j; for (j = 0; j < i; j++) { writelog("\n"); - writelog(string("Procedure: ")+ s_app->errorRoutine[j]+" \n"); - writelog(string("Error: ")+s_app->errorMessage[j]+" \n"); + writelog(string("Procedure: ")+ errorRoutine[j]+" \n"); + writelog(string("Error: ")+ errorMessage[j]+" \n"); } writelog("\n\n"); - s_app->errorMessage.clear(); - s_app->errorRoutine.clear(); + errorMessage.clear(); + errorRoutine.clear(); } /* @@ -448,9 +754,11 @@ namespace Cantera { * \ingroup errorhandling */ void setError(std::string r, std::string msg) { - appinit(); - s_app->errorMessage.push_back(msg); - s_app->errorRoutine.push_back(r); + app()->addError(r, msg) ; + } + void Application::Messages::addError(std::string r, std::string msg) { + errorMessage.push_back(msg); + errorRoutine.push_back(r); } @@ -482,9 +790,8 @@ namespace Cantera { * Additional directories may be added by calling function addDirectory. * @ingroup inputfiles */ - static void setDefaultDirectories() { - appinit(); - vector& dirs = s_app->inputDirs; + void Application::setDefaultDirectories() { + vector& dirs = inputDirs; // always look in the local directory first dirs.push_back("."); @@ -546,16 +853,20 @@ namespace Cantera { // Add a directory to the input file search path. // @ingroup inputfiles void addDirectory(std::string dir) { - appinit(); - if (s_app->inputDirs.size() == 0) setDefaultDirectories(); + app()->addDataDirectory( dir ) ; + } + void Application::addDataDirectory( std::string dir ) { + DIR_LOCK() ; + if (inputDirs.size() == 0) setDefaultDirectories(); string d = stripnonprint(dir); - size_t m, n = s_app->inputDirs.size(); + size_t m, n = inputDirs.size(); // don't add if already present - for (m = 0; m < n; m++) - if (d == s_app->inputDirs[m]) return; + for (m = 0; m < n; m++) { + if (d == inputDirs[m]) return; + } - s_app->inputDirs.push_back(stripnonprint(dir)); + inputDirs.push_back(d); } /* @@ -580,13 +891,17 @@ namespace Cantera { * If the file is not found, a message is written to * stdout and a CanteraError exception is thrown. */ - string findInputFile(std::string name) { - appinit(); + std::string findInputFile(std::string name) { + return app()->findInputFile(name) ; + } + + std::string Application::findInputFile(std::string name) { + DIR_LOCK() ; string::size_type islash = name.find('/'); string::size_type ibslash = name.find('\\'); string inname; - vector& dirs = s_app->inputDirs; - if (dirs.size() == 0) setDefaultDirectories(); + vector& dirs = inputDirs; + //if (dirs.size() == 0) setDefaultDirectories(); int nd; if (islash == string::npos && ibslash == string::npos) { @@ -650,7 +965,7 @@ namespace Cantera { // exceptions CanteraError::CanteraError(std::string proc, std::string msg) { - setError(proc, msg); + app()->addError(proc, msg); } ArraySizeError::ArraySizeError(std::string proc, int sz, int reqd) : @@ -716,22 +1031,54 @@ namespace Cantera { // is used from MATLAB or other application that do not have a // standard output stream. @ingroup textlogs void writelog(const std::string& msg) { - app()->logwriter->write(msg); + app()->writelog(msg); + } + void Application::Messages::writelog(const std::string& msg) { + logwriter->write(msg) ; } // Write a message to the screen. - void writelog(const char* msg) {writelog(string(msg));} + void writelog(const char* msg) { + app()->writelog(msg); + } + void Application::Messages::writelog(const char* pszmsg) { + logwriter->write( pszmsg ) ; + } + + void writelogf(const char* fmt,...) { + enum { BUFSIZE = 2048 } ; + char sbuf[BUFSIZE] ; + + va_list args ; + + va_start( args, fmt ) ; + + #if defined(WIN32) && defined(MSC_VER) + _vsnprintf( sbuf, BUFSIZE, fmt, args ) ; + #else + vsprintf( sbuf, fmt, args ) ; + #endif + + writelog( sbuf ) ; + + va_end(args) ; + } // Write an error message and terminate execution. test. // @ingroup textlogs void error(const std::string& msg) { - app()->logwriter->error(msg); + app()->logerror(msg); + } + void Application::Messages::logerror(const std::string& msg) { + logwriter->error(msg) ; } // @ingroup textlogs int userInterface() { - appinit(); - return app()->logwriter->env(); + return app()->getUserEnv(); + } + int Application::Messages::getUserEnv() { + return logwriter->env() ; } // Install a logger. Called by the language interfaces to install an @@ -739,9 +1086,15 @@ namespace Cantera { // @see Logger. // @ingroup textlogs void setLogger(Logger* logwriter) { - appinit(); - delete s_app->logwriter; - s_app->logwriter = logwriter; + app()->setLogger(logwriter) ; + } + void Application::Messages::setLogger(Logger* _logwriter) { + if ( logwriter == _logwriter ) return ; + if (logwriter != 0 ) { + delete logwriter; + logwriter = 0 ; + } + logwriter = _logwriter; } #ifdef WITH_HTML_LOGS @@ -778,49 +1131,68 @@ namespace Cantera { // heading, until endLogGroup() is called. // @ingroup HTML_logs void beginLogGroup(std::string title, int loglevel) { - appinit(); - if (loglevel != -99) s_app->loglevel = loglevel; - else s_app->loglevel--; - s_app->loglevels.push_back(s_app->loglevel); - s_app->loggroups.push_back(title); - if (s_app->loglevel <= 0) return; - if (s_app->xmllog == 0) { - s_app->xmllog = new XML_Node("html"); - s_app->current = &s_app->xmllog->addChild("ul"); + app()->beginLogGroup(title, loglevel) ; } - s_app->current = &s_app->current->addChild("li",""+title+""); - s_app->current = &s_app->current->addChild("ul"); + + void Application::Messages::beginLogGroup(std::string title, int _loglevel /*=-99*/) { + if (_loglevel != -99) loglevel = _loglevel; + else loglevel--; + if (loglevel <= 0) return; + loglevels.push_back(loglevel); + loggroups.push_back(title); + if (xmllog == 0) { + xmllog = new XML_Node("html"); + current = &xmllog->addChild("ul"); + } + current = ¤t->addChild("li",""+title+""); + current = ¤t->addChild("ul"); } // Add an entry to the log file. Entries appear in the form "tag: // value". // @ingroup HTML_logs void addLogEntry(std::string tag, std::string value) { - if (s_app->loglevel > 0 && s_app->current) - s_app->current->addChild("li",tag+": "+value); + app()->addLogEntry(tag, value) ; + } + + void Application::Messages::addLogEntry(std::string tag, std::string value) { + if (loglevel > 0 && current) + current->addChild("li",tag+": "+value); } // Add an entry to the log file. Entries appear in the form "tag: // value". // @ingroup HTML_logs void addLogEntry(std::string tag, doublereal value) { - if (s_app->loglevel > 0 && s_app->current) - s_app->current->addChild("li",tag+": "+fp2str(value)); + app()->addLogEntry( tag, value ) ; + } + + void Application::Messages::addLogEntry(std::string tag, doublereal value) { + if (loglevel > 0 && current) + current->addChild("li",tag+": "+fp2str(value)); } // Add an entry to the log file. Entries appear in the form "tag: // value". // @ingroup HTML_logs void addLogEntry(std::string tag, int value) { - if (s_app->loglevel > 0 && s_app->current) - s_app->current->addChild("li",tag+": "+int2str(value)); + app()->addLogEntry( tag, value ) ; + } + + void Application::Messages::addLogEntry(std::string tag, int value) { + if (loglevel > 0 && current) + current->addChild("li",tag+": "+int2str(value)); } // Add an entry to the log file. // @ingroup HTML_logs void addLogEntry(std::string msg) { - if (s_app->loglevel > 0 && s_app->current) - s_app->current->addChild("li",msg); + app()->addLogEntry(msg) ; + } + + void Application::Messages::addLogEntry(std::string msg) { + if (loglevel > 0 && current) + current->addChild("li",msg); } // Close the current group of log messages. This is typically @@ -831,27 +1203,32 @@ namespace Cantera { // beginLogGroup is called first to create a new group. // @ingroup HTML_logs void endLogGroup(std::string title) { - if (s_app->loglevel > 0) { - s_app->current = s_app->current->parent(); - s_app->current = s_app->current->parent(); - } - s_app->loglevel = s_app->loglevels.back(); - s_app->loglevels.pop_back(); - if (title != "" && title != s_app->loggroups.back()) { + app()->endLogGroup(title) ; + } + + void Application::Messages::endLogGroup(std::string title) { + if (loglevel <= 0) return; + //if (loglevel > 0) { + current = current->parent(); + current = current->parent(); + //} + loglevel = loglevels.back(); + loglevels.pop_back(); + if (title != "" && title != loggroups.back()) { writelog("Logfile error." - "\n beginLogGroup: "+ s_app->loggroups.back()+ + "\n beginLogGroup: "+ loggroups.back()+ "\n endLogGroup: "+title+"\n"); write_logfile("logerror"); //s_app->loggroups.clear(); //s_app->loglevels.clear(); } - else if (s_app->loggroups.size() == 1) { - write_logfile(s_app->loggroups.back()+"_log"); - s_app->loggroups.clear(); - s_app->loglevels.clear(); + else if (loggroups.size() == 1) { + write_logfile(loggroups.back()+"_log"); + loggroups.clear(); + loglevels.clear(); } else - s_app->loggroups.pop_back(); + loggroups.pop_back(); } // Write the HTML log file. Log entries are stored in memory in @@ -864,7 +1241,11 @@ namespace Cantera { // file will be overwritten. will be appended to the name. // @ingroup HTML_logs void write_logfile(std::string file) { - if (!s_app->xmllog) { + app()->write_logfile(file) ; + } + + void Application::Messages::write_logfile(std::string file) { + if (!xmllog) { return; } string::size_type idot = file.rfind('.'); @@ -898,35 +1279,24 @@ namespace Cantera { // existing file. Open it as an output stream, and dump the // XML (HTML) tree to it. - if (s_app->xmllog) { + if (xmllog) { ofstream f(fname.c_str()); // go to the top of the tree, and write it all. - s_app->xmllog->root().write(f); + xmllog->root().write(f); f.close(); writelog("Log file " + fname + " written.\n"); - delete s_app->xmllog; - s_app->xmllog = 0; - s_app->current = 0; + delete xmllog; + xmllog = 0; + current = 0; } } #endif // WITH_HTML_LOGS - - /* - * First we define a couple of typedefs that will - * be used throught this file - */ - //! typedef for a pointer to an XML_Node - typedef const vector nodeset_t; - //! typedef for an XML_Node - typedef XML_Node node_t; - - /// split a string at a '#' sign. Used to separate a file name /// from an id string. - static void split(const std::string& src, std::string& file, std::string& id) { + static void split_at_pound(const std::string& src, std::string& file, std::string& id) { string::size_type ipound = src.find('#'); if (ipound != string::npos) { id = src.substr(ipound+1,src.size()); @@ -937,7 +1307,6 @@ namespace Cantera { file = src; } } - /* * This routine will locate an XML node in either the input * XML tree or in another input file specified by the file @@ -958,9 +1327,9 @@ namespace Cantera { * carried out from this XML node. */ XML_Node* get_XML_Node(const std::string& file_ID, XML_Node* root) { - string fname, idstr; + std::string fname, idstr; XML_Node *db, *doc; - split(file_ID, fname, idstr); + split_at_pound(file_ID, fname, idstr); if (fname == "") { if (!root) throw CanteraError("get_XML_Node", "no file name given. file_ID = "+file_ID); @@ -1009,7 +1378,7 @@ namespace Cantera { XML_Node* root) { string fname, idTarget; XML_Node *db, *doc; - split(file_ID, fname, idTarget); + split_at_pound(file_ID, fname, idTarget); if (fname == "") { if (!root) return 0; db = root->findNameID(nameTarget, idTarget); @@ -1021,6 +1390,7 @@ namespace Cantera { return db; } - + std::vector FactoryBase::s_vFactoryRegistry; } + diff --git a/Cantera/src/base/units.h b/Cantera/src/base/units.h index 094988d61..681e1381c 100644 --- a/Cantera/src/base/units.h +++ b/Cantera/src/base/units.h @@ -17,6 +17,12 @@ #include "ct_defs.h" #include "ctexceptions.h" +#include + +#if defined(THREAD_SAFE_CANTERA) +#include +#endif + namespace Cantera { //! Unit conversion utility @@ -29,6 +35,9 @@ namespace Cantera { //! Initialize the static Unit class. static Unit* units() { + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(units_mutex) ; + #endif if (!s_u) s_u = new Unit; return s_u; } @@ -38,9 +47,12 @@ namespace Cantera { * Note this can't be done in a destructor. */ static void deleteUnit() { + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(units_mutex) ; + #endif if (s_u) { - delete s_u; - s_u = 0; + delete s_u; + s_u = 0; } } @@ -54,10 +66,10 @@ namespace Cantera { */ doublereal actEnergyToSI(std::string units) { if (m_act_u.find(units) != m_act_u.end()) { - return m_act_u[units]; + return m_act_u[units]; } else { - return toSI(units); + return toSI(units); } } @@ -83,55 +95,55 @@ namespace Cantera { while (1 > 0) { - // get token consisting of all characters up to the next - // dash, slash, or the end of the string - k = u.find_first_of("/-"); - if (k != std::string::npos) - tok = u.substr(0,k); - else - tok = u; - tsize = static_cast(tok.size()); - if (tsize == 0) - fctr = 1.0; - else if (tok[tsize - 1] == '2') { - tsub = tok.substr(0,tsize-1); - fctr = m_u[tsub]; - fctr *= fctr; - } - else if (tok[tsize - 1] == '3') { - tsub = tok.substr(0,tsize-1); - fctr = m_u[tsub]; - fctr *= fctr*fctr; - } - else if (tok[tsize - 1] == '4') { - tsub = tok.substr(0,tsize-1); - fctr = m_u[tsub]; - fctr *= fctr*fctr*fctr; - } - else if (tok[tsize - 1] == '5') { - tsub = tok.substr(0,tsize-1); - fctr = m_u[tsub]; - fctr *= fctr*fctr*fctr*fctr; - } - else if (tok[tsize - 1] == '6') { - tsub = tok.substr(0,tsize-1); - fctr = m_u[tsub]; - fctr *= fctr*fctr*fctr*fctr*fctr; - } - else { - tsub = tok; - fctr = m_u[tok]; - } + // get token consisting of all characters up to the next + // dash, slash, or the end of the string + k = u.find_first_of("/-"); + if (k != std::string::npos) + tok = u.substr(0,k); + else + tok = u; + tsize = static_cast(tok.size()); + if (tsize == 0) + fctr = 1.0; + else if (tok[tsize - 1] == '2') { + tsub = tok.substr(0,tsize-1); + fctr = m_u[tsub]; + fctr *= fctr; + } + else if (tok[tsize - 1] == '3') { + tsub = tok.substr(0,tsize-1); + fctr = m_u[tsub]; + fctr *= fctr*fctr; + } + else if (tok[tsize - 1] == '4') { + tsub = tok.substr(0,tsize-1); + fctr = m_u[tsub]; + fctr *= fctr*fctr*fctr; + } + else if (tok[tsize - 1] == '5') { + tsub = tok.substr(0,tsize-1); + fctr = m_u[tsub]; + fctr *= fctr*fctr*fctr*fctr; + } + else if (tok[tsize - 1] == '6') { + tsub = tok.substr(0,tsize-1); + fctr = m_u[tsub]; + fctr *= fctr*fctr*fctr*fctr*fctr; + } + else { + tsub = tok; + fctr = m_u[tok]; + } - // tok is not one of the entries in map m_u, then - // m_u[tok] returns 0.0. Check for this. - if (fctr == 0) - throw CanteraError("toSI","unknown unit: "+tsub); - if (action == '-') f *= fctr; - else if (action == '/') f /= fctr; - if (k == std::string::npos) break; - action = u[k]; - u = u.substr(k+1,u.size()); + // tok is not one of the entries in map m_u, then + // m_u[tok] returns 0.0. Check for this. + if (fctr == 0) + throw CanteraError("toSI","unknown unit: "+tsub); + if (action == '-') f *= fctr; + else if (action == '/') f /= fctr; + if (k == std::string::npos) break; + action = u[k]; + u = u.substr(k+1,u.size()); } return f; } @@ -156,6 +168,10 @@ namespace Cantera { */ std::map m_act_u; +#if defined(THREAD_SAFE_CANTERA) + static boost::mutex units_mutex; +#endif + /*! * Units class constructor, containing the default mappings between * strings and units. @@ -215,3 +231,4 @@ namespace Cantera { } #endif + diff --git a/Cantera/src/equil/BasisOptimize.cpp b/Cantera/src/equil/BasisOptimize.cpp index 43dc00163..75debe911 100644 --- a/Cantera/src/equil/BasisOptimize.cpp +++ b/Cantera/src/equil/BasisOptimize.cpp @@ -5,6 +5,7 @@ * stoichiometric coefficient matrix (see /ref equil functions) */ /* + /* * $Author$ * $Date$ * $Revision$ @@ -19,7 +20,7 @@ using namespace std; #ifdef DEBUG_HKM namespace Cantera { int Cantera::BasisOptimize_print_lvl = 0; -static char sbuf[1024]; +//static char sbuf[1024]; } static void print_stringTrunc(const char *str, int space, int alignment); #endif @@ -79,9 +80,9 @@ static int mlequ(double *c, int idem, int n, double *b, int m); * */ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, - MultiPhase *mphase, vector_int & orderVectorSpecies, - vector_int & orderVectorElements, - vector_fp & formRxnMatrix) { + MultiPhase *mphase, vector_int & orderVectorSpecies, + vector_int & orderVectorElements, + vector_fp & formRxnMatrix) { int j, jj, k=0, kk, l, i, jl, ml; bool lindep; @@ -131,25 +132,25 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, writelog(" --- Formula Matrix used in BASOPT calculation\n"); writelog(" --- Species | Order | "); for (j = 0; j < ne; j++) { - jj = orderVectorElements[j]; - writelog(" "); - ename = mphase->elementName(jj); - print_stringTrunc(ename.c_str(), 4, 1); - sprintf(sbuf,"(%1d)", j); writelog(sbuf); + jj = orderVectorElements[j]; + writelog(" "); + ename = mphase->elementName(jj); + print_stringTrunc(ename.c_str(), 4, 1); + writelogf("(%1d)", j); } writelog("\n"); for (k = 0; k < nspecies; k++) { - kk = orderVectorSpecies[k]; - writelog(" --- "); - sname = mphase->speciesName(kk); - print_stringTrunc(sname.c_str(), 11, 1); - sprintf(sbuf," | %4d |", k); writelog(sbuf); - for (j = 0; j < ne; j++) { - jj = orderVectorElements[j]; - double num = mphase->nAtoms(kk,jj); - sprintf(sbuf,"%6.1g ", num); writelog(sbuf); - } - writelog("\n"); + kk = orderVectorSpecies[k]; + writelog(" --- "); + sname = mphase->speciesName(kk); + print_stringTrunc(sname.c_str(), 11, 1); + writelogf(" | %4d |", k); + for (j = 0; j < ne; j++) { + jj = orderVectorElements[j]; + double num = mphase->nAtoms(kk,jj); + writelogf("%6.1g ", num); + } + writelog("\n"); } writelog(" --- \n"); } @@ -209,13 +210,13 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, */ kk = amax(DATA_PTR(molNum), 0, nspecies); for (j = 0; j < nspecies; j++) { - if (orderVectorSpecies[j] == kk) { - k = j; - break; - } + if (orderVectorSpecies[j] == kk) { + k = j; + break; + } } if (j == nspecies) { - throw CanteraError("BasisOptimize", "orderVectorSpecies contains an error"); + throw CanteraError("BasisOptimize", "orderVectorSpecies contains an error"); } if (molNum[kk] == 0.0) *usedZeroedSpecies = true; @@ -223,9 +224,9 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, * If the largest molNum is negative, then we are done. */ if (molNum[kk] == USEDBEFORE) { - nComponents = jr; - nNonComponents = nspecies - nComponents; - goto L_END_LOOP; + nComponents = jr; + nNonComponents = nspecies - nComponents; + goto L_END_LOOP; } /* * Assign a small negative number to the component that we have @@ -245,32 +246,32 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, */ jl = jr; for (j = 0; j < ne; ++j) { - jj = orderVectorElements[j]; - sm[j + jr*ne] = mphase->nAtoms(kk,jj); + jj = orderVectorElements[j]; + sm[j + jr*ne] = mphase->nAtoms(kk,jj); } if (jl > 0) { - /* - * Compute the coefficients of JA column of the - * the upper triangular R matrix, SS(J) = R_J_JR - * (this is slightly different than Dalquist) - * R_JA_JA = 1 - */ - for (j = 0; j < jl; ++j) { - ss[j] = 0.0; - for (i = 0; i < ne; ++i) { - ss[j] += sm[i + jr*ne] * sm[i + j*ne]; - } - ss[j] /= sa[j]; - } - /* - * Now make the new column, (*,JR), orthogonal to the - * previous columns - */ - for (j = 0; j < jl; ++j) { - for (l = 0; l < ne; ++l) { - sm[l + jr*ne] -= ss[j] * sm[l + j*ne]; - } - } + /* + * Compute the coefficients of JA column of the + * the upper triangular R matrix, SS(J) = R_J_JR + * (this is slightly different than Dalquist) + * R_JA_JA = 1 + */ + for (j = 0; j < jl; ++j) { + ss[j] = 0.0; + for (i = 0; i < ne; ++i) { + ss[j] += sm[i + jr*ne] * sm[i + j*ne]; + } + ss[j] /= sa[j]; + } + /* + * Now make the new column, (*,JR), orthogonal to the + * previous columns + */ + for (j = 0; j < jl; ++j) { + for (l = 0; l < ne; ++l) { + sm[l + jr*ne] -= ss[j] * sm[l + j*ne]; + } + } } /* * Find the new length of the new column in Q. @@ -278,8 +279,8 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, */ sa[jr] = 0.0; for (ml = 0; ml < ne; ++ml) { - tmp = sm[ml + jr*ne]; - sa[jr] += tmp * tmp; + tmp = sm[ml + jr*ne]; + sa[jr] += tmp * tmp; } /* **************************************************** */ /* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */ @@ -293,15 +294,13 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, if (jr != k) { #ifdef DEBUG_HKM if (BasisOptimize_print_lvl >= 1) { - kk = orderVectorSpecies[k]; - sname = mphase->speciesName(kk); - sprintf(sbuf," --- %-12.12s", sname.c_str()); writelog(sbuf); - jj = orderVectorSpecies[jr]; - ename = mphase->speciesName(jj); - sprintf(sbuf,"(%9.2g) replaces %-12.12s", molSave, ename.c_str()); - writelog(sbuf); - sprintf(sbuf,"(%9.2g) as component %3d\n", molNum[jj], jr); - writelog(sbuf); + kk = orderVectorSpecies[k]; + sname = mphase->speciesName(kk); + writelogf(" --- %-12.12s", sname.c_str()); + jj = orderVectorSpecies[jr]; + ename = mphase->speciesName(jj); + writelogf("(%9.2g) replaces %-12.12s", molSave, ename.c_str()); + writelogf("(%9.2g) as component %3d\n", molNum[jj], jr); } #endif switch_pos(orderVectorSpecies, jr, k); @@ -381,40 +380,38 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, #ifdef DEBUG_HKM if (Cantera::BasisOptimize_print_lvl >= 1) { writelog(" ---\n"); - sprintf(sbuf," --- Number of Components = %d\n", nComponents); - writelog(sbuf); + writelogf(" --- Number of Components = %d\n", nComponents); writelog(" --- Formula Matrix:\n"); writelog(" --- Components: "); for (k = 0; k < nComponents; k++) { kk = orderVectorSpecies[k]; - sprintf(sbuf," %3d (%3d) ", k, kk); writelog(sbuf); + writelogf(" %3d (%3d) ", k, kk); } writelog("\n --- Components Moles: "); for (k = 0; k < nComponents; k++) { kk = orderVectorSpecies[k]; - sprintf(sbuf,"%-11.3g", molNumBase[kk]); writelog(sbuf); + writelogf("%-11.3g", molNumBase[kk]); } writelog("\n --- NonComponent | Moles | "); for (i = 0; i < nComponents; i++) { kk = orderVectorSpecies[i]; sname = mphase->speciesName(kk); - sprintf(sbuf,"%-11.10s", sname.c_str()); writelog(sbuf); + writelogf("%-11.10s", sname.c_str()); } writelog("\n"); for (i = 0; i < nNonComponents; i++) { k = i + nComponents; kk = orderVectorSpecies[k]; - sprintf(sbuf," --- %3d (%3d) ", k, kk); writelog(sbuf); + writelogf(" --- %3d (%3d) ", k, kk); sname = mphase->speciesName(kk); - sprintf(sbuf,"%-10.10s", sname.c_str()); writelog(sbuf); - sprintf(sbuf,"|%10.3g|", molNumBase[kk]); writelog(sbuf); + writelogf("%-10.10s", sname.c_str()); + writelogf("|%10.3g|", molNumBase[kk]); /* * Print the negative of formRxnMatrix[]; it's easier to interpret. */ for (j = 0; j < nComponents; j++) { - sprintf(sbuf," %6.2f", - formRxnMatrix[j + i * ne]); - writelog(sbuf); + writelogf(" %6.2f", - formRxnMatrix[j + i * ne]); } writelog("\n"); } @@ -448,7 +445,7 @@ static void print_stringTrunc(const char *str, int space, int alignment) int len = strlen(str); if ((len) >= space) { for (i = 0; i < space; i++) { - sprintf(sbuf,"%c", str[i]); writelog(sbuf); + writelogf("%c", str[i]); } } else { if (alignment == 1) { @@ -462,7 +459,7 @@ static void print_stringTrunc(const char *str, int space, int alignment) if (ls != 0) { for (i = 0; i < ls; i++) writelog(" "); } - sprintf(sbuf,"%s", str); writelog(sbuf); + writelogf("%s", str); if (rs != 0) { for (i = 0; i < rs; i++) writelog(" "); } @@ -540,14 +537,13 @@ static int amax(double *x, int j, int n) { for (i = 0; i < n; ++i) { if (c[i + i * idem] == 0.0) { /* - * Do a simple form of row pivoting to find a non-zero pivot - */ + * Do a simple form of row pivoting to find a non-zero pivot + */ for (k = i + 1; k < n; ++k) { - if (c[k + i * idem] != 0.0) goto FOUND_PIVOT; + if (c[k + i * idem] != 0.0) goto FOUND_PIVOT; } #ifdef DEBUG_HKM - sprintf(sbuf,"vcs_mlequ ERROR: Encountered a zero column: %d\n", i); - writelog(sbuf); + writelogf("vcs_mlequ ERROR: Encountered a zero column: %d\n", i); #endif return 1; FOUND_PIVOT: ; @@ -557,10 +553,10 @@ static int amax(double *x, int j, int n) { for (l = 0; l < n; ++l) { if (l != i && c[l + i * idem] != 0.0) { - R = c[l + i * idem] / c[i + i * idem]; - c[l + i * idem] = 0.0; - for (j = i+1; j < n; ++j) c[l + j * idem] -= c[i + j * idem] * R; - for (j = 0; j < m; ++j) b[l + j * idem] -= b[i + j * idem] * R; + R = c[l + i * idem] / c[i + i * idem]; + c[l + i * idem] = 0.0; + for (j = i+1; j < n; ++j) c[l + j * idem] -= c[i + j * idem] * R; + for (j = 0; j < m; ++j) b[l + j * idem] -= b[i + j * idem] * R; } } } @@ -612,9 +608,9 @@ static int amax(double *x, int j, int n) { * nonsingular matrix to invert. */ int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances, - MultiPhase *mphase, - vector_int & orderVectorSpecies, - vector_int & orderVectorElements) { + MultiPhase *mphase, + vector_int & orderVectorSpecies, + vector_int & orderVectorElements) { int j, k, l, i, jl, ml, jr, ielem, jj, kk=0; @@ -669,12 +665,12 @@ int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances, for (j = 0; j < nelements; j++) { eAbund[j] = 0.0; for (k = 0; k < nspecies; k++) { - eAbund[j] += fabs(mphase->nAtoms(k, j)); + eAbund[j] += fabs(mphase->nAtoms(k, j)); } } } else { copy(elementAbundances.begin(), elementAbundances.end(), - eAbund.begin()); + eAbund.begin()); } vector_fp sa(nelements,0.0); @@ -701,39 +697,38 @@ int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances, */ k = nelements; for (ielem = jr; ielem < nelements; ielem++) { - kk = orderVectorElements[ielem]; - if (eAbund[kk] != test && eAbund[kk] > 0.0) { - k = ielem; - break; - } + kk = orderVectorElements[ielem]; + if (eAbund[kk] != test && eAbund[kk] > 0.0) { + k = ielem; + break; + } } for (ielem = jr; ielem < nelements; ielem++) { - kk = orderVectorElements[ielem]; - if (eAbund[kk] != test) { - k = ielem; - break; - } + kk = orderVectorElements[ielem]; + if (eAbund[kk] != test) { + k = ielem; + break; + } } if (k == nelements) { - // When we are here, there is an error usually. - // We haven't found the number of elements necessary. - // This is signalled by returning jr != nComponents. + // When we are here, there is an error usually. + // We haven't found the number of elements necessary. + // This is signalled by returning jr != nComponents. #ifdef DEBUG_HKM if (BasisOptimize_print_lvl > 0) { - sprintf(sbuf,"Error exit: returning with nComponents = %d\n", jr); - writelog(sbuf); + writelogf("Error exit: returning with nComponents = %d\n", jr); } #endif - return jr; + return jr; } - + /* * Assign a large negative number to the element that we have * just found, in order to take it out of further consideration. */ eAbund[kk] = test; - + /* *********************************************************** */ /* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX */ /* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ****** */ @@ -750,43 +745,43 @@ int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances, * (note j and k indecises are flipped compared to the previous routine) */ for (j = 0; j < nComponents; ++j) { - jj = orderVectorSpecies[j]; - kk = orderVectorElements[k]; - sm[j + jr*nComponents] = mphase->nAtoms(jj,kk); + jj = orderVectorSpecies[j]; + kk = orderVectorElements[k]; + sm[j + jr*nComponents] = mphase->nAtoms(jj,kk); } if (jl > 0) { - /* - * Compute the coefficients of JA column of the - * the upper triangular R matrix, SS(J) = R_J_JR - * (this is slightly different than Dalquist) - * R_JA_JA = 1 - */ - for (j = 0; j < jl; ++j) { - ss[j] = 0.0; - for (i = 0; i < nComponents; ++i) { - ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents]; - } - ss[j] /= sa[j]; - } - /* - * Now make the new column, (*,JR), orthogonal to the - * previous columns - */ - for (j = 0; j < jl; ++j) { - for (l = 0; l < nComponents; ++l) { - sm[l + jr*nComponents] -= ss[j] * sm[l + j*nComponents]; - } - } + /* + * Compute the coefficients of JA column of the + * the upper triangular R matrix, SS(J) = R_J_JR + * (this is slightly different than Dalquist) + * R_JA_JA = 1 + */ + for (j = 0; j < jl; ++j) { + ss[j] = 0.0; + for (i = 0; i < nComponents; ++i) { + ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents]; } - + ss[j] /= sa[j]; + } + /* + * Now make the new column, (*,JR), orthogonal to the + * previous columns + */ + for (j = 0; j < jl; ++j) { + for (l = 0; l < nComponents; ++l) { + sm[l + jr*nComponents] -= ss[j] * sm[l + j*nComponents]; + } + } + } + /* * Find the new length of the new column in Q. * It will be used in the denominator in future row calcs. */ sa[jr] = 0.0; for (ml = 0; ml < nComponents; ++ml) { - double tmp = sm[ml + jr*nComponents]; - sa[jr] += tmp * tmp; + double tmp = sm[ml + jr*nComponents]; + sa[jr] += tmp * tmp; } /* **************************************************** */ /* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */ @@ -800,15 +795,15 @@ int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances, if (jr != k) { #ifdef DEBUG_HKM if (BasisOptimize_print_lvl > 0) { - kk = orderVectorElements[k]; - ename = mphase->elementName(kk); - writelog(" --- "); - sprintf(sbuf,"%-2.2s", ename.c_str()); writelog(sbuf); - writelog("replaces "); - kk = orderVectorElements[jr]; - ename = mphase->elementName(kk); - sprintf(sbuf,"%-2.2s", ename.c_str()); writelog(sbuf); - sprintf(sbuf," as element %3d\n", jr); writelog(sbuf); + kk = orderVectorElements[k]; + ename = mphase->elementName(kk); + writelog(" --- "); + writelogf("%-2.2s", ename.c_str()); + writelog("replaces "); + kk = orderVectorElements[jr]; + ename = mphase->elementName(kk); + writelogf("%-2.2s", ename.c_str()); + writelogf(" as element %3d\n", jr); } #endif switch_pos(orderVectorElements, jr, k); @@ -822,3 +817,4 @@ int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances, } while (jr < (nComponents-1)); return nComponents; } /* vcs_elem_rearrange() ****************************************************/ + diff --git a/Cantera/src/equil/ChemEquil.cpp b/Cantera/src/equil/ChemEquil.cpp index 86ee84243..58397c24d 100755 --- a/Cantera/src/equil/ChemEquil.cpp +++ b/Cantera/src/equil/ChemEquil.cpp @@ -33,7 +33,7 @@ using namespace std; #ifdef DEBUG_HKM #include "stdio.h" int Cantera::ChemEquil_print_lvl = 0; -static char sbuf[1024]; +//static char sbuf[1024]; #endif #ifndef MIN #define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) @@ -61,9 +61,9 @@ namespace Cantera { /// Default Constructor. ChemEquil::ChemEquil() : m_skip(-1), m_p1(0), m_p2(0), m_elementTotalSum(1.0), - m_p0(OneAtm), m_eloc(-1), - m_elemFracCutoff(1.0E-100), - m_doResPerturb(false) + m_p0(OneAtm), m_eloc(-1), + m_elemFracCutoff(1.0E-100), + m_doResPerturb(false) {} //! Constructor combined with the initialization function @@ -107,7 +107,7 @@ namespace Cantera { m_nComponents = m_mm; //if (m_kk < m_mm) { //throw CanteraError("ChemEquil::initialize", - // "number of species cannot be less than the number of elements."); + // "number of species cannot be less than the number of elements."); //} // allocate space in internal work arrays within the ChemEquil object @@ -137,32 +137,32 @@ namespace Cantera { doublereal na, ewt; for (m = 0; m < m_mm; m++) { for (k = 0; k < m_kk; k++) { - na = s.nAtoms(k,m); + na = s.nAtoms(k,m); - // handle the case of negative atom numbers (used to - // represent positive ions, where the 'element' is an - // electron - if (na < 0.0) { + // handle the case of negative atom numbers (used to + // represent positive ions, where the 'element' is an + // electron + if (na < 0.0) { - // if negative atom numbers have already been specified - // for some element other than this one, throw - // an exception - if (mneg >= 0 && mneg != m) - throw CanteraError("ChemEquil::initialize", - "negative atom numbers allowed for only one element"); - mneg = m; - ewt = s.atomicWeight(m); + // if negative atom numbers have already been specified + // for some element other than this one, throw + // an exception + if (mneg >= 0 && mneg != m) + throw CanteraError("ChemEquil::initialize", + "negative atom numbers allowed for only one element"); + mneg = m; + ewt = s.atomicWeight(m); - // the element should be an electron... if it isn't - // print a warning. - if (ewt > 1.0e-3) - writelog(string("WARNING: species " - +s.speciesName(k) - +" has "+fp2str(s.nAtoms(k,m)) - +" atoms of element " - +s.elementName(m)+ - ", but this element is not an electron.\n")); - } + // the element should be an electron... if it isn't + // print a warning. + if (ewt > 1.0e-3) + writelog(string("WARNING: species " + +s.speciesName(k) + +" has "+fp2str(s.nAtoms(k,m)) + +" atoms of element " + +s.elementName(m)+ + ", but this element is not an electron.\n")); + } } } m_eloc = mneg; @@ -170,7 +170,7 @@ namespace Cantera { // set up the elemental composition matrix for (k = 0; k < m_kk; k++) { for (m = 0; m < m_mm; m++) { - m_comp[k*m_mm + m] = s.nAtoms(k,m); + m_comp[k*m_mm + m] = s.nAtoms(k,m); } } } @@ -186,13 +186,13 @@ namespace Cantera { * */ void ChemEquil::setToEquilState(thermo_t& s, - const vector_fp& lambda_RT, doublereal t) + const vector_fp& lambda_RT, doublereal t) { // Construct the chemical potentials by summing element potentials fill(m_mu_RT.begin(), m_mu_RT.end(), 0.0); for (int k = 0; k < m_kk; k++) for (int m = 0; m < m_mm; m++) - m_mu_RT[k] += lambda_RT[m]*nAtoms(k,m); + m_mu_RT[k] += lambda_RT[m]*nAtoms(k,m); // Set the temperature s.setTemperature(t); @@ -222,12 +222,12 @@ namespace Cantera { for (m = 0; m < m_mm; m++) { m_elementmolefracs[m] = 0.0; for (k = 0; k < m_kk; k++) { - m_elementmolefracs[m] += nAtoms(k,m) * m_molefractions[k]; - if (m_molefractions[k] < 0.0) { - throw CanteraError("update", - "negative mole fraction for "+s.speciesName(k)+ - ": "+fp2str(m_molefractions[k])); - } + m_elementmolefracs[m] += nAtoms(k,m) * m_molefractions[k]; + if (m_molefractions[k] < 0.0) { + throw CanteraError("update", + "negative mole fraction for "+s.speciesName(k)+ + ": "+fp2str(m_molefractions[k])); + } } sum += m_elementmolefracs[m]; } @@ -239,30 +239,33 @@ namespace Cantera { /// Estimate the initial mole numbers. This version borrows from the /// MultiPhaseEquil solver. - int ChemEquil::setInitialMoles(thermo_t& s, vector_fp & elMoleGoal) { + int ChemEquil::setInitialMoles(thermo_t& s, vector_fp & elMoleGoal, + int loglevel) { MultiPhase* mp = 0; MultiPhaseEquil* e = 0; int iok = 0; - beginLogGroup("ChemEquil::setInitialMoles"); + if (loglevel > 0) + beginLogGroup("ChemEquil::setInitialMoles"); try { mp = new MultiPhase; mp->addPhase(&s, 1.0); mp->init(); - e = new MultiPhaseEquil(mp, true); - e->setInitialMixMoles(); + e = new MultiPhaseEquil(mp, true, loglevel-1); + e->setInitialMixMoles(loglevel-1); // store component indices if (m_nComponents > m_kk) { m_nComponents = m_kk; } for (int m = 0; m < m_nComponents; m++) { - m_component[m] = e->componentIndex(m); + m_component[m] = e->componentIndex(m); } for (int k = 0; k < m_kk; k++) { - if (s.moleFraction(k) > 0.0) { - addLogEntry(s.speciesName(k), - s.moleFraction(k)); - } + if (s.moleFraction(k) > 0.0) { + if (loglevel > 0) + addLogEntry(s.speciesName(k), + s.moleFraction(k)); + } } /* * Update the current values of the temp, density, and @@ -273,20 +276,20 @@ namespace Cantera { #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - writelog("setInitialMoles: Estimated Mole Fractions\n"); - sprintf(sbuf," Temperature = %g\n", s.temperature()); writelog(sbuf); - sprintf(sbuf," Pressure = %g\n", s.pressure()); writelog(sbuf); - for (int k = 0; k < m_kk; k++) { - string nnn = s.speciesName(k); - double mf = s.moleFraction(k); - sprintf(sbuf," %-12s % -10.5g\n", nnn.c_str(), mf); writelog(sbuf); - } - writelog(" Element_Name ElementGoal ElementMF\n"); - for (int m = 0; m < m_mm; m++) { - string nnn = s.elementName(m); - sprintf(sbuf," %-12s % -10.5g% -10.5g\n", - nnn.c_str(), elMoleGoal[m], m_elementmolefracs[m]); writelog(sbuf); - } + writelog("setInitialMoles: Estimated Mole Fractions\n"); + writelogf(" Temperature = %g\n", s.temperature()); + writelogf(" Pressure = %g\n", s.pressure()); + for (int k = 0; k < m_kk; k++) { + string nnn = s.speciesName(k); + double mf = s.moleFraction(k); + writelogf(" %-12s % -10.5g\n", nnn.c_str(), mf); + } + writelog(" Element_Name ElementGoal ElementMF\n"); + for (int m = 0; m < m_mm; m++) { + string nnn = s.elementName(m); + writelogf(" %-12s % -10.5g% -10.5g\n", + nnn.c_str(), elMoleGoal[m], m_elementmolefracs[m]); + } } #endif @@ -299,7 +302,8 @@ namespace Cantera { delete mp; iok = -1; } - endLogGroup(); + if (loglevel > 0) + endLogGroup(); return iok; } @@ -308,10 +312,11 @@ namespace Cantera { * Generate a starting estimate for the element potentials. */ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, - vector_fp& elMolesGoal) + vector_fp& elMolesGoal, int loglevel) { int m, n; - beginLogGroup("estimateElementPotentials"); + if (loglevel > 0) + beginLogGroup("estimateElementPotentials"); //for (k = 0; k < m_kk; k++) { // if (m_molefractions[k] > 0.0) { // m_molefractions[k] = fmaxx(m_molefractions[k], 0.05); @@ -327,7 +332,7 @@ namespace Cantera { s.getMoleFractions(DATA_PTR(xMF_est)); for (n = 0; n < s.nSpecies(); n++) { if (xMF_est[n] < 1.0E-20) { - xMF_est[n] = 1.0E-20; + xMF_est[n] = 1.0E-20; } } s.setMoleFractions(DATA_PTR(xMF_est)); @@ -339,24 +344,24 @@ namespace Cantera { int usedZeroedSpecies = 0; vector_fp formRxnMatrix; m_nComponents = BasisOptimize(&usedZeroedSpecies, false, - mp, m_orderVectorSpecies, - m_orderVectorElements, formRxnMatrix); + mp, m_orderVectorSpecies, + m_orderVectorElements, formRxnMatrix); for (m = 0; m < m_nComponents; m++) { int k = m_orderVectorSpecies[m]; m_component[m] = k; if (xMF_est[k] < 1.0E-8) { - xMF_est[k] = 1.0E-8; + xMF_est[k] = 1.0E-8; } } s.setMoleFractions(DATA_PTR(xMF_est)); s.getMoleFractions(DATA_PTR(xMF_est)); int nct = Cantera::ElemRearrange(m_nComponents, elMolesGoal, mp, - m_orderVectorSpecies, m_orderVectorElements); + m_orderVectorSpecies, m_orderVectorElements); if (nct != m_nComponents) { throw CanteraError("ChemEquil::estimateElementPotentials", - "confused"); + "confused"); } delete mp; @@ -369,28 +374,26 @@ namespace Cantera { #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { for (m = 0; m < m_nComponents; m++) { - int isp = m_component[m]; - string nnn = s.speciesName(isp); - sprintf(sbuf,"isp = %d, %s\n", isp, nnn.c_str()); - writelog(sbuf); + int isp = m_component[m]; + string nnn = s.speciesName(isp); + writelogf("isp = %d, %s\n", isp, nnn.c_str()); } double pres = s.pressure(); double temp = s.temperature(); - sprintf(sbuf,"Pressure = %g\n", pres); writelog(sbuf); - sprintf(sbuf,"Temperature = %g\n", temp); writelog(sbuf); + writelogf("Pressure = %g\n", pres); + writelogf("Temperature = %g\n", temp); writelog(" id Name MF mu/RT \n"); for (n = 0; n < s.nSpecies(); n++) { - string nnn = s.speciesName(n); - sprintf(sbuf,"%10d %15s %10.5g %10.5g\n", - n, nnn.c_str(), xMF_est[n], mu_RT[n]); - writelog(sbuf); + string nnn = s.speciesName(n); + writelogf("%10d %15s %10.5g %10.5g\n", + n, nnn.c_str(), xMF_est[n], mu_RT[n]); } } #endif DenseMatrix aa(m_nComponents, m_nComponents, 0.0); for (m = 0; m < m_nComponents; m++) { for (n = 0; n < m_nComponents; n++) { - aa(m,n) = nAtoms(m_component[m], m_orderVectorElements[n]); + aa(m,n) = nAtoms(m_component[m], m_orderVectorElements[n]); } b[m] = mu_RT[m_component[m]]; } @@ -400,8 +403,9 @@ namespace Cantera { info = solve(aa, DATA_PTR(b)); } catch (CanteraError) { - addLogEntry("failed to estimate initial element potentials."); - info = -2; + if (loglevel > 0) + addLogEntry("failed to estimate initial element potentials."); + info = -2; } for (m = 0; m < m_nComponents; m++) { lambda_RT[m_orderVectorElements[m]] = b[m]; @@ -410,35 +414,36 @@ namespace Cantera { lambda_RT[m_orderVectorElements[m]] = 0.0; } if (info == 0) { - for (m = 0; m < m_mm; m++) { - addLogEntry(s.elementName(m),lambda_RT[m]); - } + if (loglevel > 0) { + for (m = 0; m < m_mm; m++) { + addLogEntry(s.elementName(m),lambda_RT[m]); + } + } } #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { writelog(" id CompSpecies ChemPot EstChemPot Diff\n"); for (m = 0; m < m_nComponents; m++) { - int isp = m_component[m]; - double tmp = 0.0; - string sname = s.speciesName(isp); - for (n = 0; n < m_mm; n++) { - tmp += nAtoms(isp, n) * lambda_RT[n]; - } - sprintf(sbuf,"%3d %16s %10.5g %10.5g %10.5g\n", - m, sname.c_str(), mu_RT[isp], tmp, tmp - mu_RT[isp]); - writelog(sbuf); + int isp = m_component[m]; + double tmp = 0.0; + string sname = s.speciesName(isp); + for (n = 0; n < m_mm; n++) { + tmp += nAtoms(isp, n) * lambda_RT[n]; + } + writelogf("%3d %16s %10.5g %10.5g %10.5g\n", + m, sname.c_str(), mu_RT[isp], tmp, tmp - mu_RT[isp]); } writelog(" id ElName Lambda_RT\n"); for (m = 0; m < m_mm; m++) { - string ename = s.elementName(m); - sprintf(sbuf," %3d %6s %10.5g\n", m, ename.c_str(), lambda_RT[m]); - writelog(sbuf); + string ename = s.elementName(m); + writelogf(" %3d %6s %10.5g\n", m, ename.c_str(), lambda_RT[m]); } } #endif - endLogGroup(); + if (loglevel > 0) + endLogGroup(); return info; } @@ -452,13 +457,14 @@ namespace Cantera { * within the current thermodynamic state of the system. */ int ChemEquil::equilibrate(thermo_t& s, const char* XY, - bool useThermoPhaseElementPotentials) { + bool useThermoPhaseElementPotentials, int loglevel) { vector_fp elMolesGoal(s.nElements()); initialize(s); update(s); copy(m_elementmolefracs.begin(), m_elementmolefracs.end(), - elMolesGoal.begin()); - return equilibrate(s, XY, elMolesGoal, useThermoPhaseElementPotentials); + elMolesGoal.begin()); + return equilibrate(s, XY, elMolesGoal, useThermoPhaseElementPotentials, + loglevel-1); } @@ -477,8 +483,10 @@ namespace Cantera { * Unsuccessful returns are indicated by a return value of -1 for * lack of convergence or -3 for a singular jacobian. */ - int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, vector_fp& elMolesGoal, - bool useThermoPhaseElementPotentials) + int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, + vector_fp& elMolesGoal, + bool useThermoPhaseElementPotentials, + int loglevel) { doublereal xval, yval, tmp; int fail = 0; @@ -497,14 +505,15 @@ namespace Cantera { */ if (m_mm != s.nElements() || m_kk != s.nSpecies()) { throw CanteraError("ChemEquil::equilibrate ERROR", - "Input ThermoPhase is incompatible with initialization"); + "Input ThermoPhase is incompatible with initialization"); } #ifdef DEBUG_HKM int n; const vector& eNames = s.elementNames(); #endif - beginLogGroup("ChemEquil::equilibrate"); + if (loglevel > 0) + beginLogGroup("ChemEquil::equilibrate"); initialize(s); update(s); switch (XY) { @@ -537,24 +546,26 @@ namespace Cantera { m_p2 = new DensityCalculator; break; default: - endLogGroup("ChemEquil::equilibrate"); + if (loglevel > 0) + endLogGroup("ChemEquil::equilibrate"); throw CanteraError("equilibrate","illegal property pair."); } - - addLogEntry("Problem type","fixed "+m_p1->symbol()+", "+m_p2->symbol()); - addLogEntry(m_p1->symbol(), m_p1->value(s)); - addLogEntry(m_p2->symbol(), m_p2->value(s)); - + if (loglevel > 0) { + addLogEntry("Problem type","fixed "+m_p1->symbol()+", "+m_p2->symbol()); + addLogEntry(m_p1->symbol(), m_p1->value(s)); + addLogEntry(m_p2->symbol(), m_p2->value(s)); + } // If the temperature is one of the specified variables, and // it is outside the valid range, throw an exception. if (tempFixed) { double tfixed = s.temperature(); if (tfixed > s.maxTemp() + 1.0 || tfixed < s.minTemp() - 1.0) { - endLogGroup("ChemEquil::equilibrate"); - throw CanteraError("ChemEquil","Specified temperature (" - +fp2str(s.temperature())+" K) outside " - "valid range of "+fp2str(s.minTemp())+" K to " - +fp2str(s.maxTemp())+" K\n"); + if (loglevel > 0) + endLogGroup("ChemEquil::equilibrate"); + throw CanteraError("ChemEquil","Specified temperature (" + +fp2str(s.temperature())+" K) outside " + "valid range of "+fp2str(s.minTemp())+" K to " + +fp2str(s.maxTemp())+" K\n"); } } @@ -583,13 +594,13 @@ namespace Cantera { for (im = 0; im < m_nComponents; im++) { m = m_orderVectorElements[im]; if (elMolesGoal[m] > tmp ) { - m_skip = m; - tmp = elMolesGoal[m]; + m_skip = m; + tmp = elMolesGoal[m]; } } if (tmp <= 0.0) { throw CanteraError("ChemEquil", - "Element Abundance Vector is zeroed"); + "Element Abundance Vector is zeroed"); } // start with a composition with everything non-zero. Note @@ -610,8 +621,8 @@ namespace Cantera { // loop to estimate T if (!tempFixed) { - - beginLogGroup("Initial T Estimate"); + if (loglevel > 0) + beginLogGroup("Initial T Estimate"); doublereal tmax = s.maxTemp(); doublereal tmin = s.minTemp(); @@ -623,11 +634,11 @@ namespace Cantera { // bounnds (phigh, plow) for p1. s.setTemperature(tmax); - setInitialMoles(s, elMolesGoal); + setInitialMoles(s, elMolesGoal, loglevel - 1); phigh = m_p1->value(s); s.setTemperature(tmin); - setInitialMoles(s, elMolesGoal); + setInitialMoles(s, elMolesGoal, loglevel - 1); plow = m_p1->value(s); // start with T at the midpoint of the range @@ -637,43 +648,45 @@ namespace Cantera { // loop up to 5 times for (int it = 0; it < 5; it++) { - // set the composition and get p1 - setInitialMoles(s, elMolesGoal); - pval = m_p1->value(s); + // set the composition and get p1 + setInitialMoles(s, elMolesGoal, loglevel - 1); + pval = m_p1->value(s); - // If this value of p1 is greater than the specified - // property value, then the current temperature is too - // high. Use it as the new upper bound. Otherwise, it - // is too low, so use it as the new lower bound. - if (pval > xval) { - tmax = t0; - phigh = pval; - } - else { - tmin = t0; - plow = pval; - } - - // Determine the new T estimate by linearly intepolation - // between the upper and lower bounds - slope = (phigh - plow)/(tmax - tmin); - dt = (xval - plow)/slope; - - // If within 100 K, terminate the search - if (fabs(dt) < 100.0) break; - - // update the T estimate - t0 = tmin + dt; - addLogEntry("new T estimate", t0); - - s.setTemperature(t0); - } - endLogGroup("Initial T Estimate"); // initial T estimate + // If this value of p1 is greater than the specified + // property value, then the current temperature is too + // high. Use it as the new upper bound. Otherwise, it + // is too low, so use it as the new lower bound. + if (pval > xval) { + tmax = t0; + phigh = pval; + } + else { + tmin = t0; + plow = pval; } - - setInitialMoles(s, elMolesGoal); + // Determine the new T estimate by linearly intepolation + // between the upper and lower bounds + slope = (phigh - plow)/(tmax - tmin); + dt = (xval - plow)/slope; + + // If within 100 K, terminate the search + if (fabs(dt) < 100.0) break; + + // update the T estimate + t0 = tmin + dt; + if (loglevel > 0) + addLogEntry("new T estimate", t0); + + s.setTemperature(t0); + } + if (loglevel > 0) + endLogGroup("Initial T Estimate"); // initial T estimate + } + + + setInitialMoles(s, elMolesGoal,loglevel); /* * If requested, get the initial estimate for the @@ -683,12 +696,12 @@ namespace Cantera { if (useThermoPhaseElementPotentials) { bool haveEm = s.getElementPotentials(DATA_PTR(x)); if (haveEm) { - doublereal rt = GasConstant * s.temperature(); - for (m = 0; m < m_mm; m++) { - x[m] /= rt; - } + doublereal rt = GasConstant * s.temperature(); + for (m = 0; m < m_mm; m++) { + x[m] /= rt; + } } else { - estimateElementPotentials(s, x, elMolesGoal); + estimateElementPotentials(s, x, elMolesGoal); } } else { /* @@ -716,9 +729,9 @@ namespace Cantera { int info = estimateEP_Brinkley(s, x, elMolesGoal); if (info != 0) { if (info == 1) { - addLogEntry("estimateEP_Brinkley didn't converge in given max interations"); + addLogEntry("estimateEP_Brinkley didn't converge in given max interations"); } else if (info == -3) { - addLogEntry("estimateEP_Brinkley had a singular Jacobian. Continuing anyway"); + addLogEntry("estimateEP_Brinkley had a singular Jacobian. Continuing anyway"); } } else { setToEquilState(s, x, s.temperature()); @@ -765,7 +778,8 @@ namespace Cantera { iter++; if (iter > 1) endLogGroup("Iteration "+int2str(iter-1)); // iteration - beginLogGroup("Iteration "+int2str(iter)); + if (loglevel > 0) + beginLogGroup("Iteration "+int2str(iter)); // compute the residual and the jacobian using the current // solution vector @@ -778,29 +792,28 @@ namespace Cantera { #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - sprintf(sbuf,"Jacobian matrix %d:\n", iter); writelog(sbuf); + writelogf("Jacobian matrix %d:\n", iter); for (m = 0; m <= m_mm; m++) { - writelog(" [ "); - for (n = 0; n <= m_mm; n++) { - sprintf(sbuf,"%10.5g ", jac(m,n)); writelog(sbuf); - } - writelog(" ]"); - char xName[32]; - if (m < m_mm) { - string nnn = eNames[m]; - sprintf(xName, "x_%-10s", nnn.c_str()); - } else { - sprintf(xName, "x_XX"); - } - if (m_eloc == m) { - sprintf(xName, "x_ELOC"); - } - if (m == m_skip) { - sprintf(xName, "x_YY"); - } - sprintf(sbuf,"%-12s", xName); writelog(sbuf); - sprintf(sbuf, " = - (%10.5g)\n", res_trial[m]); - writelog(sbuf); + writelog(" [ "); + for (n = 0; n <= m_mm; n++) { + writelogf("%10.5g ", jac(m,n)); + } + writelog(" ]"); + char xName[32]; + if (m < m_mm) { + string nnn = eNames[m]; + sprintf(xName, "x_%-10s", nnn.c_str()); + } else { + sprintf(xName, "x_XX"); + } + if (m_eloc == m) { + sprintf(xName, "x_ELOC"); + } + if (m == m_skip) { + sprintf(xName, "x_YY"); + } + writelogf("%-12s", xName); + writelogf(" = - (%10.5g)\n", res_trial[m]); } } #endif @@ -825,9 +838,9 @@ namespace Cantera { s.restoreState(state); throw CanteraError("equilibrate", - "Jacobian is singular. \nTry adding more species, " - "changing the elemental composition slightly, \nor removing " - "unused elements."); + "Jacobian is singular. \nTry adding more species, " + "changing the elemental composition slightly, \nor removing " + "unused elements."); return -3; } @@ -837,32 +850,32 @@ namespace Cantera { for (m = 0; m < nvar; m++) { newval = x[m] + res_trial[m]; if (newval > above[m]) { - fctr = fmaxx( 0.0, fminn( fctr, - 0.8*(above[m] - x[m])/(newval - x[m]))); + fctr = fmaxx( 0.0, fminn( fctr, + 0.8*(above[m] - x[m])/(newval - x[m]))); } else if (newval < below[m]) { - fctr = fminn(fctr, 0.8*(x[m] - below[m])/(x[m] - newval)); + fctr = fminn(fctr, 0.8*(x[m] - below[m])/(x[m] - newval)); } } if (fctr != 1.0) addLogEntry("factor to keep solution in bounds", - fctr); + fctr); // multiply the step by the scaling factor scale(res_trial.begin(), res_trial.end(), res_trial.begin(), fctr); if (!dampStep(s, oldx, oldf, grad, res_trial, - x, f, elMolesGoal , xval, yval)) + x, f, elMolesGoal , xval, yval)) { - fail++; - if (fail > 3) { - addLogEntry("dampStep","Failed 3 times. Giving up."); - endLogGroup(); // iteration - endLogGroup(); // equilibrate - s.restoreState(state); - throw CanteraError("equilibrate", - "Cannot find an acceptable Newton damping coefficient."); - return -4; - } + fail++; + if (fail > 3) { + addLogEntry("dampStep","Failed 3 times. Giving up."); + endLogGroup(); // iteration + endLogGroup(); // equilibrate + s.restoreState(state); + throw CanteraError("equilibrate", + "Cannot find an acceptable Newton damping coefficient."); + return -4; + } } else fail = 0; @@ -881,67 +894,75 @@ namespace Cantera { for (m = 0; m < nvar; m++) { double tval = options.relTolerance; if (m < mm) { - if (m == m_eloc) { - tval = elMolesGoal[m] * options.relTolerance + options.absElemTol - + 1.0E-15; - } else { - tval = elMolesGoal[m] * options.relTolerance + options.absElemTol; - } + if (m == m_eloc) { + tval = elMolesGoal[m] * options.relTolerance + options.absElemTol + + 1.0E-15; + } else { + tval = elMolesGoal[m] * options.relTolerance + options.absElemTol; + } } if (fabs(res_trial[m]) > tval) { - passThis = false; + passThis = false; } } if (iter > 0 && passThis - && fabs(deltax) < options.relTolerance - && fabs(deltay) < options.relTolerance) { + && fabs(deltax) < options.relTolerance + && fabs(deltay) < options.relTolerance) { options.iterations = iter; - endLogGroup("Iteration "+int2str(iter)); // iteration - beginLogGroup("Converged solution"); - addLogEntry("Iterations",iter); - addLogEntry("Relative error in "+m_p1->symbol(),deltax); - addLogEntry("Relative error in "+m_p2->symbol(),deltay); - addLogEntry("Max residual",rmax); - beginLogGroup("Element potentials"); + if (loglevel > 0) { + endLogGroup("Iteration "+int2str(iter)); // iteration + beginLogGroup("Converged solution"); + addLogEntry("Iterations",iter); + addLogEntry("Relative error in "+m_p1->symbol(),deltax); + addLogEntry("Relative error in "+m_p2->symbol(),deltay); + addLogEntry("Max residual",rmax); + beginLogGroup("Element potentials"); + } doublereal rt = GasConstant* s.temperature(); for (m = 0; m < m_mm; m++) { - m_lambda[m] = x[m]*rt; - addLogEntry("element "+ s.elementName(m), fp2str(x[m])); + m_lambda[m] = x[m]*rt; + if (loglevel > 0) + addLogEntry("element "+ s.elementName(m), fp2str(x[m])); } if (m_eloc >= 0) { - adjustEloc(s, elMolesGoal); + adjustEloc(s, elMolesGoal); } /* * Save the calculated and converged element potentials * to the original ThermoPhase object. */ s.setElementPotentials(m_lambda); - addLogEntry("Saving Element Potentials to ThermoPhase Object"); - endLogGroup("Element potentials"); - - if (s.temperature() > s.maxTemp() + 1.0 || - s.temperature() < s.minTemp() - 1.0 ) { - writelog("Warning: Temperature (" - +fp2str(s.temperature())+" K) outside " - "valid range of "+fp2str(s.minTemp())+" K to " - +fp2str(s.maxTemp())+" K\n"); + if (loglevel > 0) { + addLogEntry("Saving Element Potentials to ThermoPhase Object"); + endLogGroup("Element potentials"); + } + if (s.temperature() > s.maxTemp() + 1.0 || + s.temperature() < s.minTemp() - 1.0 ) { + writelog("Warning: Temperature (" + +fp2str(s.temperature())+" K) outside " + "valid range of "+fp2str(s.minTemp())+" K to " + +fp2str(s.maxTemp())+" K\n"); + } + if (loglevel > 0) { + endLogGroup("Converged solution"); + endLogGroup("ChemEquil::equilibrate"); } - endLogGroup("Converged solution"); - endLogGroup("ChemEquil::equilibrate"); return 0; } // no convergence if (iter > options.maxIterations) { - addLogEntry("equilibrate","no convergence"); - endLogGroup("Iteration "+int2str(iter)); - endLogGroup("ChemEquil::equilibrate"); + if (loglevel > 0) { + addLogEntry("equilibrate","no convergence"); + endLogGroup("Iteration "+int2str(iter)); + endLogGroup("ChemEquil::equilibrate"); + } s.restoreState(state); throw CanteraError("equilibrate", - "no convergence in "+int2str(options.maxIterations) - +" iterations."); + "no convergence in "+int2str(options.maxIterations) + +" iterations."); return -1; } goto next; @@ -962,8 +983,8 @@ namespace Cantera { * Near convergence, the delta damping gets out of the way. */ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx, - double oldf, vector_fp& grad, vector_fp& step, vector_fp& x, - double& f, vector_fp& elmols, double xval, double yval ) + double oldf, vector_fp& grad, vector_fp& step, vector_fp& x, + double& f, vector_fp& elmols, double xval, double yval ) { int nvar = x.size(); int m; @@ -975,19 +996,19 @@ namespace Cantera { damp = 1.0; for (m = 0; m < m_mm; m++) { if (m == m_eloc) { - if (step[m] > 1.25) { - damp = MIN(damp, 1.25 /step[m]); - } - if (step[m] < -1.25) { - damp = MIN(damp, -1.25 / step[m]); - } + if (step[m] > 1.25) { + damp = MIN(damp, 1.25 /step[m]); + } + if (step[m] < -1.25) { + damp = MIN(damp, -1.25 / step[m]); + } } else { - if (step[m] > 0.75) { - damp = MIN(damp, 0.75 /step[m]); - } - if (step[m] < -0.75) { - damp = MIN(damp, -0.75 / step[m]); - } + if (step[m] > 0.75) { + damp = MIN(damp, 0.75 /step[m]); + } + if (step[m] < -0.75) { + damp = MIN(damp, -0.75 / step[m]); + } } } @@ -999,11 +1020,10 @@ namespace Cantera { } #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - sprintf(sbuf, "Solution Unknowns: damp = %g\n", damp); writelog(sbuf); + writelogf("Solution Unknowns: damp = %g\n", damp); writelog(" X_new X_old Step\n"); for (m = 0; m < nvar; m++) { - sprintf(sbuf," % -10.5g % -10.5g % -10.5g\n", x[m], oldx[m], step[m]); - writelog(sbuf); + writelogf(" % -10.5g % -10.5g % -10.5g\n", x[m], oldx[m], step[m]); } } #endif @@ -1015,10 +1035,11 @@ namespace Cantera { * Evaluates the residual vector F, of length mm */ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, - const vector_fp& elmFracGoal, vector_fp& resid, - doublereal xval, doublereal yval) + const vector_fp& elmFracGoal, vector_fp& resid, + doublereal xval, doublereal yval, int loglevel) { - beginLogGroup("ChemEquil::equilResidual"); + if (loglevel > 0) + beginLogGroup("ChemEquil::equilResidual"); int n, m; doublereal xx, yy; doublereal temp = exp(x[m_mm]); @@ -1030,32 +1051,32 @@ namespace Cantera { m = m_orderVectorElements[n]; // drive element potential for absent elements to -1000 if (elmFracGoal[m] < m_elemFracCutoff && m != m_eloc) { - resid[m] = x[m] + 1000.0; + resid[m] = x[m] + 1000.0; } else if (n >= m_nComponents) { - resid[m] = x[m]; + resid[m] = x[m]; } else { - /* - * Change the calculation for small element number, using - * L'Hopital's rule. - * The log formulation is unstable. - */ - if (elmFracGoal[m] < 1.0E-10 || elmFrac[m] < 1.0E-10 || m == m_eloc) { - resid[m] = elmFracGoal[m] - elmFrac[m]; - } else { - resid[m] = log( (1.0 + elmFracGoal[m]) / (1.0 + elmFrac[m]) ); - } + /* + * Change the calculation for small element number, using + * L'Hopital's rule. + * The log formulation is unstable. + */ + if (elmFracGoal[m] < 1.0E-10 || elmFrac[m] < 1.0E-10 || m == m_eloc) { + resid[m] = elmFracGoal[m] - elmFrac[m]; + } else { + resid[m] = log( (1.0 + elmFracGoal[m]) / (1.0 + elmFrac[m]) ); + } } - addLogEntry(s.elementName(m),fp2str(elmFrac[m])+" (" - +fp2str(elmFracGoal[m])+")"); + if (loglevel > 0) + addLogEntry(s.elementName(m),fp2str(elmFrac[m])+" (" + +fp2str(elmFracGoal[m])+")"); } #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0 && !m_doResPerturb) { writelog("Residual: ElFracGoal ElFracCurrent Resid\n"); for (n = 0; n < m_mm; n++) { - sprintf(sbuf," % -14.7E % -14.7E % -10.5E\n", - elmFracGoal[n], elmFrac[n], resid[n]); - writelog(sbuf); + writelogf(" % -14.7E % -14.7E % -10.5E\n", + elmFracGoal[n], elmFrac[n], resid[n]); } } #endif @@ -1064,19 +1085,19 @@ namespace Cantera { yy = m_p2->value(s); resid[m_mm] = xx/xval - 1.0; resid[m_skip] = yy/yval - 1.0; - string xstr = fp2str(xx)+" ("+fp2str(xval)+")"; - addLogEntry(m_p1->symbol(), xstr); - string ystr = fp2str(yy)+" ("+fp2str(yval)+")"; - addLogEntry(m_p2->symbol(), ystr); - endLogGroup("ChemEquil::equilResidual"); + if (loglevel > 0) { + string xstr = fp2str(xx)+" ("+fp2str(xval)+")"; + addLogEntry(m_p1->symbol(), xstr); + string ystr = fp2str(yy)+" ("+fp2str(yval)+")"; + addLogEntry(m_p2->symbol(), ystr); + endLogGroup("ChemEquil::equilResidual"); + } #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0 && !m_doResPerturb) { writelog(" Goal Xvalue Resid\n"); - sprintf(sbuf," XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[m_mm]); - writelog(sbuf); - sprintf(sbuf," YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]); - writelog(sbuf); + writelogf(" XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[m_mm]); + writelogf(" YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]); } #endif } @@ -1085,10 +1106,11 @@ namespace Cantera { //-------------------- Jacobian evaluation --------------------------- void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x, - const vector_fp& elmols, DenseMatrix& jac, - doublereal xval, doublereal yval) + const vector_fp& elmols, DenseMatrix& jac, + doublereal xval, doublereal yval, int loglevel) { - beginLogGroup("equilJacobian"); + if (loglevel > 0) + beginLogGroup("equilJacobian"); int len = x.size(); vector_fp& r0 = m_jwork1; vector_fp& r1 = m_jwork2; @@ -1098,7 +1120,7 @@ namespace Cantera { doublereal rdx, dx, xsave, dx2; doublereal atol = 1.e-10; - equilResidual(s, x, elmols, r0, xval, yval); + equilResidual(s, x, elmols, r0, xval, yval, loglevel-1); m_doResPerturb = false; for (n = 0; n < len; n++) { @@ -1112,19 +1134,20 @@ namespace Cantera { // calculate perturbed residual - equilResidual(s, x, elmols, r1, xval, yval); + equilResidual(s, x, elmols, r1, xval, yval, loglevel-1); // compute nth column of Jacobian for (m = 0; m < len; m++) { - jac(m, n) = (r1[m] - r0[m])*rdx; + jac(m, n) = (r1[m] - r0[m])*rdx; } x[n] = xsave; } m_doResPerturb = false; - endLogGroup("equilJacobian"); + if (loglevel > 0) + endLogGroup("equilJacobian"); } - + /** * Given a vector of dimensionless element abundances, * this routine calculates the moles of the elements and @@ -1133,10 +1156,10 @@ namespace Cantera { * -------- * x[m] = current dimensionless element potentials.. */ - double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double & n_t, - const vector_fp & Xmol_i_calc, - vector_fp& eMolesCalc, vector_fp& n_i_calc, - double pressureConst) { + double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double & n_t, + const vector_fp & Xmol_i_calc, + vector_fp& eMolesCalc, vector_fp& n_i_calc, + double pressureConst) { int k, m; double n_t_calc = 0.0; double tmp; @@ -1152,20 +1175,20 @@ namespace Cantera { for (k = 0; k < m_kk; k++) { tmp = - (m_muSS_RT[k] + log(actCoeff[k])); for (m = 0; m < m_mm; m++) { - tmp += nAtoms(k,m) * x[m]; + tmp += nAtoms(k,m) * x[m]; } if (tmp > 100.) tmp = 100.; if (tmp < -300.) { - n_i_calc[k] = 0.0; + n_i_calc[k] = 0.0; } else { - n_i_calc[k] = n_t * exp(tmp); + n_i_calc[k] = n_t * exp(tmp); } n_t_calc += n_i_calc[k]; } for (m = 0; m < m_mm; m++) { eMolesCalc[m] = 0.0; for (k = 0; k < m_kk; k++) { - eMolesCalc[m] += nAtoms(k,m) * n_i_calc[k]; + eMolesCalc[m] += nAtoms(k,m) * n_i_calc[k]; } } return n_t_calc; @@ -1208,7 +1231,7 @@ namespace Cantera { * NOTE: update for activity coefficients. */ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, - vector_fp& elMoles) { + vector_fp& elMoles) { /* * Before we do anything, we will save the state of the solution. * Then, if things go drastically wrong, we will restore the @@ -1253,22 +1276,22 @@ namespace Cantera { for (m = 0; m < m_mm; m++) { elMolesTotal += elMoles[m]; for (k = 0; k < m_kk; k++) { - eMolesFix[m] += nAtoms(k,m) * n_i[k]; + eMolesFix[m] += nAtoms(k,m) * n_i[k]; } } for (m = 0; m < m_mm; m++) { if (x[m] > 50.0) { - x[m] = 50.; + x[m] = 50.; } if (elMoles[m] > 1.0E-70) { - if (x[m] < -100) { - x[m] = -100.; - } + if (x[m] < -100) { + x[m] = -100.; + } } else { - if (x[m] < -1000.) { - x[m] = -1000.; - } + if (x[m] < -1000.) { + x[m] = -1000.; + } } } @@ -1283,17 +1306,17 @@ namespace Cantera { tmp = - (m_muSS_RT[k] + log(actCoeff[k])); sum2 = 0.0; for (m = 0; m < m_mm; m++) { - sum = nAtoms(k,m); - tmp += sum * x[m]; - sum2 += sum; - if (sum2 > nAtomsMax) { - nAtomsMax = sum2; - } + sum = nAtoms(k,m); + tmp += sum * x[m]; + sum2 += sum; + if (sum2 > nAtomsMax) { + nAtomsMax = sum2; + } } if (tmp > 100.) { - n_t += 2.8E43; + n_t += 2.8E43; } else { - n_t += exp(tmp); + n_t += exp(tmp); } } @@ -1304,31 +1327,29 @@ namespace Cantera { writelog("estimateEP_Brinkley::\n\n"); double temp = s.temperature(); double pres = s.pressure(); - sprintf(sbuf, "temp = %g\n", temp); writelog(sbuf); - sprintf(sbuf, "pres = %g\n", pres); writelog(sbuf); + writelogf("temp = %g\n", temp); + writelogf("pres = %g\n", pres); writelog("Initial mole numbers and mu_SS:\n"); writelog(" Name MoleNum mu_SS actCoeff\n"); for (k = 0; k < m_kk; k++) { - string nnn = s.speciesName(k); - sprintf(sbuf,"%15s %13.5g %13.5g %13.5g\n", - nnn.c_str(), n_i[k], m_muSS_RT[k], actCoeff[k]); - writelog(sbuf); + string nnn = s.speciesName(k); + writelogf("%15s %13.5g %13.5g %13.5g\n", + nnn.c_str(), n_i[k], m_muSS_RT[k], actCoeff[k]); } - sprintf(sbuf,"Initial n_t = %10.5g\n", n_t); writelog(sbuf); + writelogf("Initial n_t = %10.5g\n", n_t); writelog("Comparison of Goal Element Abundance with Initial Guess:\n"); writelog(" eName eCurrent eGoal\n"); for (m = 0; m < m_mm; m++) { - string nnn = s.elementName(m); - sprintf(sbuf,"%5s %13.5g %13.5g\n",nnn.c_str(), eMolesFix[m], elMoles[m]); - writelog(sbuf); + string nnn = s.elementName(m); + writelogf("%5s %13.5g %13.5g\n",nnn.c_str(), eMolesFix[m], elMoles[m]); } } #endif for (m = 0; m < m_mm; m++) { if (m != m_eloc) { - if (elMoles[m] <= options.absElemTol) { - x[m] = -200.; - } + if (elMoles[m] <= options.absElemTol) { + x[m] = -200.; + } } } /* @@ -1340,7 +1361,7 @@ namespace Cantera { * Save the old solution */ for (m = 0; m < m_mm; m++) { - x_old[m] = x[m]; + x_old[m] = x[m]; } x_old[m_mm] = n_t; /* @@ -1348,35 +1369,31 @@ namespace Cantera { */ #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - sprintf(sbuf, "START ITERATION %d:\n", iter); writelog(sbuf); + writelogf("START ITERATION %d:\n", iter); } #endif /* * Calculate the mole numbers of species and elements. */ double n_t_calc = calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc, - pressureConst); + pressureConst); for (k = 0; k < m_kk; k++) { - Xmol_i_calc[k] = n_i_calc[k]/n_t_calc; + Xmol_i_calc[k] = n_i_calc[k]/n_t_calc; } #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n"); - for (k = 0; k < m_kk; k++) { - string nnn = s.speciesName(k); - sprintf(sbuf,"%15s: %10.5g %10.5g\n", nnn.c_str(), n_i_calc[k], Xmol_i_calc[k]); - writelog(sbuf); - } - sprintf(sbuf,"%15s: %10.5g\n", "Total Molar Sum", n_t_calc); - writelog(sbuf); - sprintf(sbuf,"(iter %d) element moles bal: Goal Calculated\n", iter); - writelog(sbuf); - for (m = 0; m < m_mm; m++) { - string nnn = eNames[m]; - sprintf(sbuf," %8s: %10.5g %10.5g \n", nnn.c_str(), elMoles[m], eMolesCalc[m]); - writelog(sbuf); - } + writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n"); + for (k = 0; k < m_kk; k++) { + string nnn = s.speciesName(k); + writelogf("%15s: %10.5g %10.5g\n", nnn.c_str(), n_i_calc[k], Xmol_i_calc[k]); + } + writelogf("%15s: %10.5g\n", "Total Molar Sum", n_t_calc); + writelogf("(iter %d) element moles bal: Goal Calculated\n", iter); + for (m = 0; m < m_mm; m++) { + string nnn = eNames[m]; + writelogf(" %8s: %10.5g %10.5g \n", nnn.c_str(), elMoles[m], eMolesCalc[m]); + } } #endif @@ -1388,54 +1405,53 @@ namespace Cantera { */ int iM = -1; for (m = 0; m < m_mm; m++) { - if (elMoles[m] > 0.001 * elMolesTotal) { - if (eMolesCalc[m] > 1000. * elMoles[m]) { - normalStep = false; - iM = m; - } - if (1000 * eMolesCalc[m] < elMoles[m]) { - normalStep = false; - iM = m; - } - } + if (elMoles[m] > 0.001 * elMolesTotal) { + if (eMolesCalc[m] > 1000. * elMoles[m]) { + normalStep = false; + iM = m; + } + if (1000 * eMolesCalc[m] < elMoles[m]) { + normalStep = false; + iM = m; + } + } } #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - if (!normalStep) { - sprintf(sbuf," NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM); - writelog(sbuf); - } + if (!normalStep) { + writelogf(" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM); + } } #endif if (!normalStep) { - beta = 1.0; - resid[m_mm] = 0.0; - for (im = 0; im < m_mm; im++) { - m = m_orderVectorElements[im]; - resid[m] = 0.0; - if (im < m_nComponents) { - if (elMoles[m] > 0.001 * elMolesTotal) { - if (eMolesCalc[m] > 1000. * elMoles[m]) { - resid[m] = -0.5; - resid[m_mm] -= 0.5; - } - if (1000 * eMolesCalc[m] < elMoles[m]) { - resid[m] = 0.5; - resid[m_mm] += 0.5; - } - } - } - } - if (n_t < (elMolesTotal / nAtomsMax)) { - if (resid[m_mm] < 0.0) { - resid[m_mm] = 0.1; - } - } else if (n_t > elMolesTotal) { - if (resid[m_mm] > 0.0) { - resid[m_mm] = 0.0; - } - } - goto updateSolnVector; + beta = 1.0; + resid[m_mm] = 0.0; + for (im = 0; im < m_mm; im++) { + m = m_orderVectorElements[im]; + resid[m] = 0.0; + if (im < m_nComponents) { + if (elMoles[m] > 0.001 * elMolesTotal) { + if (eMolesCalc[m] > 1000. * elMoles[m]) { + resid[m] = -0.5; + resid[m_mm] -= 0.5; + } + if (1000 * eMolesCalc[m] < elMoles[m]) { + resid[m] = 0.5; + resid[m_mm] += 0.5; + } + } + } + } + if (n_t < (elMolesTotal / nAtomsMax)) { + if (resid[m_mm] < 0.0) { + resid[m_mm] = 0.1; + } + } else if (n_t > elMolesTotal) { + if (resid[m_mm] > 0.0) { + resid[m_mm] = 0.0; + } + } + goto updateSolnVector; } @@ -1464,7 +1480,7 @@ namespace Cantera { * Hopefully, it's caught by the equal rows logic below. */ for (m = 0; m < m_mm; m++) { - lumpSum[m] = 1; + lumpSum[m] = 1; } nCutoff = 1.0E-9 * n_t_calc; @@ -1474,34 +1490,33 @@ namespace Cantera { } #endif for (m = 0; m < m_mm; m++) { - int kMSp = -1; - int kMSp2 = -1; - int nSpeciesWithElem = 0; - for (k = 0; k < m_kk; k++) { - if (n_i_calc[k] > nCutoff) { - if (fabs(nAtoms(k,m)) > 0.001) { - nSpeciesWithElem++; - if (kMSp != -1) { - kMSp2 = k; - double factor = fabs(nAtoms(kMSp,m) / nAtoms(kMSp2,m)); - for (n = 0; n < m_mm; n++) { - if (fabs(factor * nAtoms(kMSp2,n) - nAtoms(kMSp,n)) > 1.0E-8) { - lumpSum[m] = 0; - break; - } - } - } else { - kMSp = k; - } - } - } - } + int kMSp = -1; + int kMSp2 = -1; + int nSpeciesWithElem = 0; + for (k = 0; k < m_kk; k++) { + if (n_i_calc[k] > nCutoff) { + if (fabs(nAtoms(k,m)) > 0.001) { + nSpeciesWithElem++; + if (kMSp != -1) { + kMSp2 = k; + double factor = fabs(nAtoms(kMSp,m) / nAtoms(kMSp2,m)); + for (n = 0; n < m_mm; n++) { + if (fabs(factor * nAtoms(kMSp2,n) - nAtoms(kMSp,n)) > 1.0E-8) { + lumpSum[m] = 0; + break; + } + } + } else { + kMSp = k; + } + } + } + } #ifdef DEBUG_HKM - if (ChemEquil_print_lvl > 0) { - string nnn = eNames[m]; - sprintf(sbuf," %5s %3d : %5d %5d\n",nnn.c_str(), lumpSum[m], kMSp, kMSp2); - writelog(sbuf); - } + if (ChemEquil_print_lvl > 0) { + string nnn = eNames[m]; + writelogf(" %5s %3d : %5d %5d\n",nnn.c_str(), lumpSum[m], kMSp, kMSp2); + } #endif } @@ -1509,22 +1524,22 @@ namespace Cantera { * Formulate the matrix. */ for (im = 0; im < m_mm; im++) { - m = m_orderVectorElements[im]; - if (im < m_nComponents) { - for (n = 0; n < m_mm; n++) { - a1(m,n) = 0.0; - for (k = 0; k < m_kk; k++) { - a1(m,n) += nAtoms(k,m) * nAtoms(k,n) * n_i_calc[k]; - } - } - a1(m,m_mm) = eMolesCalc[m]; - a1(m_mm, m) = eMolesCalc[m]; - } else { - for (n = 0; n <= m_mm; n++) { - a1(m,n) = 0.0; - } - a1(m,m) = 1.0; - } + m = m_orderVectorElements[im]; + if (im < m_nComponents) { + for (n = 0; n < m_mm; n++) { + a1(m,n) = 0.0; + for (k = 0; k < m_kk; k++) { + a1(m,n) += nAtoms(k,m) * nAtoms(k,n) * n_i_calc[k]; + } + } + a1(m,m_mm) = eMolesCalc[m]; + a1(m_mm, m) = eMolesCalc[m]; + } else { + for (n = 0; n <= m_mm; n++) { + a1(m,n) = 0.0; + } + a1(m,m) = 1.0; + } } a1(m_mm, m_mm) = 0.0; @@ -1533,46 +1548,45 @@ namespace Cantera { */ sum = 0.0; for (im = 0; im < m_mm; im++) { - m = m_orderVectorElements[im]; - if (im < m_nComponents) { - resid[m] = elMoles[m] - eMolesCalc[m]; - } else { - resid[m] = 0.0; - } - /* - * For equations with positive and negative coefficients, (electronic charge), - * we must mitigate the convergence criteria by a condition limited by - * finite precision of inverting a matrix. - * Other equations with just positive coefficients aren't limited by this. - */ - if (m == m_eloc) { - tmp = resid[m] / (elMoles[m] + elMolesTotal*1.0E-6 + options.absElemTol); - } else { - tmp = resid[m] / (elMoles[m] + options.absElemTol); - } - sum += tmp * tmp; + m = m_orderVectorElements[im]; + if (im < m_nComponents) { + resid[m] = elMoles[m] - eMolesCalc[m]; + } else { + resid[m] = 0.0; + } + /* + * For equations with positive and negative coefficients, (electronic charge), + * we must mitigate the convergence criteria by a condition limited by + * finite precision of inverting a matrix. + * Other equations with just positive coefficients aren't limited by this. + */ + if (m == m_eloc) { + tmp = resid[m] / (elMoles[m] + elMolesTotal*1.0E-6 + options.absElemTol); + } else { + tmp = resid[m] / (elMoles[m] + options.absElemTol); + } + sum += tmp * tmp; } for (m = 0; m < m_mm; m++) { - if (a1(m,m) < 1.0E-50) { + if (a1(m,m) < 1.0E-50) { #ifdef DEBUG_HKM - if (ChemEquil_print_lvl > 0) { - sprintf(sbuf," NOTE: Diagonalizing the analytical Jac row %d\n", m); - writelog(sbuf); - } + if (ChemEquil_print_lvl > 0) { + writelogf(" NOTE: Diagonalizing the analytical Jac row %d\n", m); + } #endif - for (n = 0; n < m_mm; n++) { - a1(m,n) = 0.0; - } - a1(m,m) = 1.0; - if (resid[m] > 0.0) { - resid[m] = 1.0; - } else if (resid[m] < 0.0) { - resid[m] = -1.0; - } else { - resid[m] = 0.0; - } - } + for (n = 0; n < m_mm; n++) { + a1(m,n) = 0.0; + } + a1(m,m) = 1.0; + if (resid[m] > 0.0) { + resid[m] = 1.0; + } else if (resid[m] < 0.0) { + resid[m] = -1.0; + } else { + resid[m] = 0.0; + } + } } @@ -1580,14 +1594,14 @@ namespace Cantera { #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - writelog("Matrix:\n"); - for (m = 0; m <= m_mm; m++) { - writelog(" ["); - for (n = 0; n <= m_mm; n++) { - sprintf(sbuf," %10.5g", a1(m,n)); writelog(sbuf); - } - sprintf(sbuf,"] = %10.5g\n", resid[m]); writelog(sbuf); - } + writelog("Matrix:\n"); + for (m = 0; m <= m_mm; m++) { + writelog(" ["); + for (n = 0; n <= m_mm; n++) { + writelogf(" %10.5g", a1(m,n)); + } + writelogf("] = %10.5g\n", resid[m]); + } } #endif @@ -1595,8 +1609,7 @@ namespace Cantera { sum += tmp * tmp; #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - sprintf(sbuf,"(it %d) Convergence = %g\n", iter, sum); - writelog(sbuf); + writelogf("(it %d) Convergence = %g\n", iter, sum); } #endif /* @@ -1607,49 +1620,48 @@ namespace Cantera { * singular. */ if (sum < 0.05 * options.relTolerance) { - retn = 0; - goto exit; + retn = 0; + goto exit; } /* * Row Sum scaling */ for (m = 0; m <= m_mm; m++) { - tmp = 0.0; - for (n = 0; n <= m_mm; n++) { - tmp += fabs(a1(m,n)); - } - if (m < m_mm && tmp < 1.0E-30) { + tmp = 0.0; + for (n = 0; n <= m_mm; n++) { + tmp += fabs(a1(m,n)); + } + if (m < m_mm && tmp < 1.0E-30) { #ifdef DEBUG_HKM - if (ChemEquil_print_lvl > 0) { - sprintf(sbuf," NOTE: Diagonalizing row %d\n", m); - writelog(sbuf); - } + if (ChemEquil_print_lvl > 0) { + writelogf(" NOTE: Diagonalizing row %d\n", m); + } #endif - for (n = 0; n <= m_mm; n++) { - if (n != m) { - a1(m,n) = 0.0; - a1(n,m) = 0.0; - } - } - } - tmp = 1.0/tmp; - for (n = 0; n <= m_mm; n++) { - a1(m,n) *= tmp; - } - resid[m] *= tmp; + for (n = 0; n <= m_mm; n++) { + if (n != m) { + a1(m,n) = 0.0; + a1(n,m) = 0.0; + } + } + } + tmp = 1.0/tmp; + for (n = 0; n <= m_mm; n++) { + a1(m,n) *= tmp; + } + resid[m] *= tmp; } #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - writelog("Row Summed Matrix:\n"); - for (m = 0; m <= m_mm; m++) { - writelog(" ["); - for (n = 0; n <= m_mm; n++) { - sprintf(sbuf," %10.5g", a1(m,n)); writelog(sbuf); - } - sprintf(sbuf,"] = %10.5g\n", resid[m]); writelog(sbuf); - } + writelog("Row Summed Matrix:\n"); + for (m = 0; m <= m_mm; m++) { + writelog(" ["); + for (n = 0; n <= m_mm; n++) { + writelogf(" %10.5g", a1(m,n)); + } + writelogf("] = %10.5g\n", resid[m]); + } } #endif @@ -1676,70 +1688,68 @@ namespace Cantera { */ modifiedMatrix = false; for (m = 0; m < m_mm; m++) { - int sameAsRow = -1; - for (int im = 0; im < m; im++) { - bool theSame = true; - for (n = 0; n < m_mm; n++) { - if (fabs(a1(m,n) - a1(im,n)) > 1.0E-7) { - theSame = false; - break; - } - } - if (theSame) { - sameAsRow = im; - } - } - if (sameAsRow >= 0 || lumpSum[m]) { + int sameAsRow = -1; + for (int im = 0; im < m; im++) { + bool theSame = true; + for (n = 0; n < m_mm; n++) { + if (fabs(a1(m,n) - a1(im,n)) > 1.0E-7) { + theSame = false; + break; + } + } + if (theSame) { + sameAsRow = im; + } + } + if (sameAsRow >= 0 || lumpSum[m]) { #ifdef DEBUG_HKM - if (ChemEquil_print_lvl > 0) { - if (lumpSum[m]) { - sprintf(sbuf,"Lump summing row %d, due to rank deficiency analysis\n", m); - writelog(sbuf); - } else if (sameAsRow >= 0) { - sprintf(sbuf,"Identified that rows %d and %d are the same\n", m, sameAsRow); - writelog(sbuf); - } - } + if (ChemEquil_print_lvl > 0) { + if (lumpSum[m]) { + writelogf("Lump summing row %d, due to rank deficiency analysis\n", m); + } else if (sameAsRow >= 0) { + writelogf("Identified that rows %d and %d are the same\n", m, sameAsRow); + } + } #endif - modifiedMatrix = true; - for (n = 0; n < m_mm; n++) { - if (n != m) { - a1(m,m) += fabs(a1(m,n)); - a1(m,n) = 0.0; - } - } - } + modifiedMatrix = true; + for (n = 0; n < m_mm; n++) { + if (n != m) { + a1(m,m) += fabs(a1(m,n)); + a1(m,n) = 0.0; + } + } + } } #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0 && modifiedMatrix) { - writelog("Row Summed, MODIFIED Matrix:\n"); - for (m = 0; m <= m_mm; m++) { - writelog(" ["); - for (n = 0; n <= m_mm; n++) { - sprintf(sbuf," %10.5g", a1(m,n)); writelog(sbuf); - } - sprintf(sbuf,"] = %10.5g\n", resid[m]); writelog(sbuf); - } + writelog("Row Summed, MODIFIED Matrix:\n"); + for (m = 0; m <= m_mm; m++) { + writelog(" ["); + for (n = 0; n <= m_mm; n++) { + writelogf(" %10.5g", a1(m,n)); + } + writelogf("] = %10.5g\n", resid[m]); + } } #endif try { - info = solve(a1, DATA_PTR(resid)); + info = solve(a1, DATA_PTR(resid)); } catch (CanteraError) { - addLogEntry("estimateEP_Brinkley:Jacobian is singular."); + addLogEntry("estimateEP_Brinkley:Jacobian is singular."); #ifdef DEBUG_HKM - if (ChemEquil_print_lvl > 0) { - writelog("Matrix is SINGULAR.ERROR\n"); - } + if (ChemEquil_print_lvl > 0) { + writelog("Matrix is SINGULAR.ERROR\n"); + } #endif - s.restoreState(state); - throw CanteraError("equilibrate:estimateEP_Brinkley()", - "Jacobian is singular. \nTry adding more species, " - "changing the elemental composition slightly, \nor removing " - "unused elements."); - return -3; + s.restoreState(state); + throw CanteraError("equilibrate:estimateEP_Brinkley()", + "Jacobian is singular. \nTry adding more species, " + "changing the elemental composition slightly, \nor removing " + "unused elements."); + return -3; } /* @@ -1749,24 +1759,24 @@ namespace Cantera { */ beta = 1.0; for (m = 0; m < m_mm; m++) { - if (resid[m] > 1.0) { - double betat = 1.0 / resid[m]; - if (betat < beta) { - beta = betat; - } - } - if (resid[m] < -1.0) { - double betat = -1.0 / resid[m]; - if (betat < beta) { - beta = betat; - } - } + if (resid[m] > 1.0) { + double betat = 1.0 / resid[m]; + if (betat < beta) { + beta = betat; + } + } + if (resid[m] < -1.0) { + double betat = -1.0 / resid[m]; + if (betat < beta) { + beta = betat; + } + } } #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - if (beta != 1.0) { - sprintf(sbuf,"(it %d) Beta = %g\n", iter, beta); writelog(sbuf); - } + if (beta != 1.0) { + writelogf("(it %d) Beta = %g\n", iter, beta); + } } #endif @@ -1775,40 +1785,35 @@ namespace Cantera { */ updateSolnVector: for (m = 0; m < m_mm; m++) { - x[m] += beta * resid[m]; + x[m] += beta * resid[m]; } n_t *= exp(beta * resid[m_mm]); #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - sprintf(sbuf,"(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter); - writelog(sbuf); - for (m = 0; m < m_mm; m++) { - string eee = eNames[m]; - sprintf(sbuf," %5s %10.5g %10.5g %10.5g\n", eee.c_str(), x_old[m], x[m], resid[m]); - writelog(sbuf); - } - sprintf(sbuf," n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm])); - writelog(sbuf); + writelogf("(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter); + for (m = 0; m < m_mm; m++) { + string eee = eNames[m]; + writelogf(" %5s %10.5g %10.5g %10.5g\n", eee.c_str(), x_old[m], x[m], resid[m]); + } + writelogf(" n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm])); } #endif } exit: #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - double temp = s.temperature(); - double pres = s.pressure(); + double temp = s.temperature(); + double pres = s.pressure(); - if (retn == 0) { - sprintf(sbuf," ChemEquil::estimateEP_Brinkley() SUCCESS: equilibrium found at T = %g, Pres = %g\n", - temp, pres); - writelog(sbuf); - } else { - sprintf(sbuf," ChemEquil::estimateEP_Brinkley() FAILURE: equilibrium not found at T = %g, Pres = %g\n", - temp, pres); - writelog(sbuf); - } + if (retn == 0) { + writelogf(" ChemEquil::estimateEP_Brinkley() SUCCESS: equilibrium found at T = %g, Pres = %g\n", + temp, pres); + } else { + writelogf(" ChemEquil::estimateEP_Brinkley() FAILURE: equilibrium not found at T = %g, Pres = %g\n", + temp, pres); + } } #endif return retn; @@ -1831,18 +1836,18 @@ namespace Cantera { double maxNegVal = -1.0; if (ChemEquil_print_lvl > 0) { for (k = 0; k < m_kk; k++) { - if (nAtoms(k,m_eloc) > 0.0) { - if (m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) { - maxPosVal = m_molefractions[k]; - maxPosEloc = k; - } - } - if (nAtoms(k,m_eloc) < 0.0) { - if (m_molefractions[k] > maxNegVal && m_molefractions[k] > 0.0) { - maxNegVal = m_molefractions[k]; - maxNegEloc = k; - } - } + if (nAtoms(k,m_eloc) > 0.0) { + if (m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) { + maxPosVal = m_molefractions[k]; + maxPosEloc = k; + } + } + if (nAtoms(k,m_eloc) < 0.0) { + if (m_molefractions[k] > maxNegVal && m_molefractions[k] > 0.0) { + maxNegVal = m_molefractions[k]; + maxNegEloc = k; + } + } } } #endif @@ -1851,10 +1856,10 @@ namespace Cantera { double sumNeg = 0.0; for (k = 0; k < m_kk; k++) { if (nAtoms(k,m_eloc) > 0.0) { - sumPos += nAtoms(k,m_eloc) * m_molefractions[k]; + sumPos += nAtoms(k,m_eloc) * m_molefractions[k]; } if (nAtoms(k,m_eloc) < 0.0) { - sumNeg += nAtoms(k,m_eloc) * m_molefractions[k]; + sumNeg += nAtoms(k,m_eloc) * m_molefractions[k]; } } sumNeg = - sumNeg; @@ -1864,37 +1869,35 @@ namespace Cantera { double factor = (elMolesGoal[m_eloc] + sumNeg) / sumPos; #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - if (factor < 0.9999999999) { - string nnn = s.speciesName(maxPosEloc); - sprintf(sbuf,"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", - nnn.c_str(), - m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor); - writelog(sbuf); - } + if (factor < 0.9999999999) { + string nnn = s.speciesName(maxPosEloc); + writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", + nnn.c_str(), + m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor); + } } #endif for (k = 0; k < m_kk; k++) { - if (nAtoms(k,m_eloc) > 0.0) { - m_molefractions[k] *= factor; - } + if (nAtoms(k,m_eloc) > 0.0) { + m_molefractions[k] *= factor; + } } } else { double factor = (-elMolesGoal[m_eloc] + sumPos) / sumNeg; #ifdef DEBUG_HKM if (ChemEquil_print_lvl > 0) { - if (factor < 0.9999999999) { - string nnn = s.speciesName(maxNegEloc); - sprintf(sbuf,"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", - nnn.c_str(), - m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor); - writelog(sbuf); - } + if (factor < 0.9999999999) { + string nnn = s.speciesName(maxNegEloc); + writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", + nnn.c_str(), + m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor); + } } #endif for (k = 0; k < m_kk; k++) { - if (nAtoms(k,m_eloc) < 0.0) { - m_molefractions[k] *= factor; - } + if (nAtoms(k,m_eloc) < 0.0) { + m_molefractions[k] *= factor; + } } } @@ -1904,3 +1907,4 @@ namespace Cantera { } } // namespace + diff --git a/Cantera/src/equil/ChemEquil.h b/Cantera/src/equil/ChemEquil.h index 8ce080eb4..4e8b4c74d 100755 --- a/Cantera/src/equil/ChemEquil.h +++ b/Cantera/src/equil/ChemEquil.h @@ -112,9 +112,9 @@ namespace Cantera { virtual ~ChemEquil(); int equilibrate(thermo_t& s, const char* XY, - bool useThermoPhaseElementPotentials = false); + bool useThermoPhaseElementPotentials = false, int loglevel = 0); int equilibrate(thermo_t& s, const char* XY, vector_fp& elMoles, - bool useThermoPhaseElementPotentials = false); + bool useThermoPhaseElementPotentials = false, int loglevel = 0); const vector_fp& elementPotentials() const { return m_lambda; } /** @@ -144,10 +144,10 @@ namespace Cantera { void setToEquilState(thermo_t& s, const vector_fp& x, doublereal t); - int setInitialMoles(thermo_t& s, vector_fp& elMoleGoal); + int setInitialMoles(thermo_t& s, vector_fp& elMoleGoal, int loglevel = 0); int estimateElementPotentials(thermo_t& s, vector_fp& lambda, - vector_fp& elMolesGoal); + vector_fp& elMolesGoal, int loglevel = 0); int estimateEP_Brinkley(thermo_t&s, vector_fp& lambda, vector_fp& elMoles); @@ -157,11 +157,11 @@ namespace Cantera { void equilResidual(thermo_t& s, const vector_fp& x, const vector_fp& elmtotal, vector_fp& resid, - double xval, double yval); + double xval, double yval, int loglevel = 0); void equilJacobian(thermo_t& s, vector_fp& x, const vector_fp& elmols, DenseMatrix& jac, - double xval, double yval); + double xval, double yval, int loglevel = 0); void adjustEloc(thermo_t& s, vector_fp & elMolesGoal); diff --git a/Cantera/src/equil/MultiPhase.cpp b/Cantera/src/equil/MultiPhase.cpp index 145016881..1d4ef9b2c 100644 --- a/Cantera/src/equil/MultiPhase.cpp +++ b/Cantera/src/equil/MultiPhase.cpp @@ -491,13 +491,15 @@ namespace Cantera { MultiPhaseEquil* e = 0; if (!m_init) init(); - beginLogGroup("MultiPhase::equilibrate", loglevel); + if (loglevel > 0) + beginLogGroup("MultiPhase::equilibrate", loglevel); if (XY == TP) { - addLogEntry("problem type","fixed T,P"); - addLogEntry("Temperature",temperature()); - addLogEntry("Pressure", pressure()); - + if (loglevel > 0) { + addLogEntry("problem type","fixed T,P"); + addLogEntry("Temperature",temperature()); + addLogEntry("Pressure", pressure()); + } // create an equilibrium manager e = new MultiPhaseEquil(this); @@ -505,7 +507,8 @@ namespace Cantera { error = e->equilibrate(XY, err, maxsteps); } catch (CanteraError err) { - endLogGroup(); + if (loglevel > 0) + endLogGroup(); delete e; e = 0; throw err; @@ -517,9 +520,10 @@ namespace Cantera { h0 = enthalpy(); Tlow = 0.5*m_Tmin; // lower bound on T Thigh = 2.0*m_Tmax; // upper bound on T - addLogEntry("problem type","fixed H,P"); - addLogEntry("H target",fp2str(h0)); - + if (loglevel > 0) { + addLogEntry("problem type","fixed H,P"); + addLogEntry("H target",fp2str(h0)); + } for (n = 0; n < maxiter; n++) { // if 'strt' is false, the current composition will be used as @@ -531,7 +535,8 @@ namespace Cantera { e = new MultiPhaseEquil(this, strt); // start with a loose error tolerance, but tighten it as we get // close to the final temperature - beginLogGroup("iteration "+int2str(n)); + if (loglevel > 0) + beginLogGroup("iteration "+int2str(n)); try { error = e->equilibrate(TP, err, maxsteps); @@ -567,18 +572,22 @@ namespace Cantera { } herr = fabs((h0 - hnow)/h0); - addLogEntry("T",fp2str(temperature())); - addLogEntry("H",fp2str(hnow)); - addLogEntry("H rel error",fp2str(herr)); - addLogEntry("lower T bound",fp2str(Tlow)); - addLogEntry("upper T bound",fp2str(Thigh)); - endLogGroup(); // iteration + if (loglevel > 0) { + addLogEntry("T",fp2str(temperature())); + addLogEntry("H",fp2str(hnow)); + addLogEntry("H rel error",fp2str(herr)); + addLogEntry("lower T bound",fp2str(Tlow)); + addLogEntry("upper T bound",fp2str(Thigh)); + endLogGroup(); // iteration + } if (herr < err) { // || dta < 1.0e-4) { - addLogEntry("T iterations",int2str(n)); - addLogEntry("Final T",fp2str(temperature())); - addLogEntry("H rel error",fp2str(herr)); + if (loglevel > 0) { + addLogEntry("T iterations",int2str(n)); + addLogEntry("Final T",fp2str(temperature())); + addLogEntry("H rel error",fp2str(herr)); + } goto done; } tnew = m_temp + dt; @@ -593,25 +602,30 @@ namespace Cantera { } catch (CanteraError err) { - if (!strt) { - addLogEntry("no convergence", - "try estimating starting composition"); + if (!strt) { + if (loglevel > 0) + addLogEntry("no convergence", + "try estimating starting composition"); strt = true; } else { tnew = 0.5*(m_temp + Thigh); if (fabs(tnew - m_temp) < 1.0) tnew = m_temp + 1.0; setTemperature(tnew); - addLogEntry("no convergence", - "trying T = "+fp2str(m_temp)); + if (loglevel > 0) + addLogEntry("no convergence", + "trying T = "+fp2str(m_temp)); } - endLogGroup(); + if (loglevel > 0) + endLogGroup(); } delete e; e = 0; } - addLogEntry("reached max number of T iterations",int2str(maxiter)); - endLogGroup(); + if (loglevel > 0) { + addLogEntry("reached max number of T iterations",int2str(maxiter)); + endLogGroup(); + } throw CanteraError("MultiPhase::equilibrate", "No convergence for T"); } @@ -620,18 +634,20 @@ namespace Cantera { start = true; Tlow = 1.0; // m_Tmin; // lower bound on T Thigh = 1.0e6; // m_Tmax; // upper bound on T - addLogEntry("problem type","fixed S,P"); - addLogEntry("S target",fp2str(s0)); - addLogEntry("min T",fp2str(Tlow)); - addLogEntry("max T",fp2str(Thigh)); - + if (loglevel > 0) { + addLogEntry("problem type","fixed S,P"); + addLogEntry("S target",fp2str(s0)); + addLogEntry("min T",fp2str(Tlow)); + addLogEntry("max T",fp2str(Thigh)); + } for (n = 0; n < maxiter; n++) { if (e) delete e; e = new MultiPhaseEquil(this, strt); ferr = 0.1; if (fabs(dt) < 1.0) ferr = err; //start = false; - beginLogGroup("iteration "+int2str(n)); + if (loglevel > 0) + beginLogGroup("iteration "+int2str(n)); try { error = e->equilibrate(TP, err, maxsteps); @@ -643,20 +659,23 @@ namespace Cantera { if (m_temp < Thigh) Thigh = m_temp; } serr = fabs((s0 - snow)/s0); - addLogEntry("T",fp2str(temperature())); - addLogEntry("S",fp2str(snow)); - addLogEntry("S rel error",fp2str(serr)); - endLogGroup(); - + if (loglevel > 0) { + addLogEntry("T",fp2str(temperature())); + addLogEntry("S",fp2str(snow)); + addLogEntry("S rel error",fp2str(serr)); + endLogGroup(); + } dt = (s0 - snow)*m_temp/cp(); dtmax = 0.5*fabs(Thigh - Tlow); dtmax = (dtmax > 500.0 ? 500.0 : dtmax); dta = fabs(dt); if (dta > dtmax) dt *= dtmax/dta; if (herr < err || dta < 1.0e-4) { - addLogEntry("T iterations",int2str(n)); - addLogEntry("Final T",fp2str(temperature())); - addLogEntry("S rel error",fp2str(serr)); + if (loglevel > 0) { + addLogEntry("T iterations",int2str(n)); + addLogEntry("Final T",fp2str(temperature())); + addLogEntry("S rel error",fp2str(serr)); + } goto done; } tnew = m_temp + dt; @@ -669,24 +688,30 @@ namespace Cantera { catch (CanteraError err) { if (!strt) { - addLogEntry("no convergence", - "setting strt to True"); + if (loglevel > 0) { + addLogEntry("no convergence", + "setting strt to True"); + } strt = true; } else { tnew = 0.5*(m_temp + Thigh); setTemperature(tnew); - addLogEntry("no convergence", - "trying T = "+fp2str(m_temp)); - + if (loglevel > 0) { + addLogEntry("no convergence", + "trying T = "+fp2str(m_temp)); + } } - endLogGroup(); + if (loglevel > 0) + endLogGroup(); } delete e; e = 0; } - addLogEntry("reached max number of T iterations",int2str(maxiter)); - endLogGroup(); + if (loglevel > 0) { + addLogEntry("reached max number of T iterations",int2str(maxiter)); + endLogGroup(); + } throw CanteraError("MultiPhase::equilibrate", "No convergence for T"); } @@ -785,14 +810,16 @@ namespace Cantera { } else { - endLogGroup(); + if (loglevel > 0) + endLogGroup(); throw CanteraError("MultiPhase::equilibrate","unknown option"); } return -1.0; done: delete e; e = 0; - endLogGroup(); + if (loglevel > 0) + endLogGroup(); return err; } diff --git a/Cantera/src/equil/MultiPhaseEquil.cpp b/Cantera/src/equil/MultiPhaseEquil.cpp index 171d51567..0ec39e84b 100644 --- a/Cantera/src/equil/MultiPhaseEquil.cpp +++ b/Cantera/src/equil/MultiPhaseEquil.cpp @@ -39,7 +39,7 @@ namespace Cantera { /// @param start If true, the initial composition will be /// determined by a linear Gibbs minimization, otherwise the /// initial mixture composition will be used. - MultiPhaseEquil::MultiPhaseEquil(mix_t* mix, bool start) : m_mix(mix) + MultiPhaseEquil::MultiPhaseEquil(mix_t* mix, bool start, int loglevel) : m_mix(mix) { // the multi-phase mixture // m_mix = mix; @@ -165,7 +165,7 @@ namespace Cantera { // only the elemental composition of the initial mixture state // matters. if (start) { - setInitialMoles(); + setInitialMoles(loglevel-1); } computeN(); @@ -198,28 +198,37 @@ namespace Cantera { int i; m_iter = 0; string iterstr; - beginLogGroup("MultiPhaseEquil::equilibrate", loglevel); + if (loglevel > 0) + beginLogGroup("MultiPhaseEquil::equilibrate", loglevel); for (i = 0; i < maxsteps; i++) { - iterstr = "iteration "+int2str(i); - beginLogGroup(iterstr); - stepComposition(); - addLogEntry("error",fp2str(error())); - endLogGroup(iterstr); + if (loglevel > 0) { + iterstr = "iteration "+int2str(i); + beginLogGroup(iterstr); + } + stepComposition(loglevel-1); + if (loglevel > 0) { + addLogEntry("error",fp2str(error())); + endLogGroup(iterstr); + } if (error() < err) break; } if (i >= maxsteps) { - addLogEntry("Error","no convergence in "+int2str(maxsteps) - +" iterations"); - endLogGroup("MultiPhaseEquil::equilibrate"); + if (loglevel > 0) { + addLogEntry("Error","no convergence in "+int2str(maxsteps) + +" iterations"); + endLogGroup("MultiPhaseEquil::equilibrate"); + } throw CanteraError("MultiPhaseEquil::equilibrate", "no convergence in " + int2str(maxsteps) + " iterations. Error = " + fp2str(error())); } - addLogEntry("iterations",int2str(iterations())); - addLogEntry("error tolerance",fp2str(err)); - addLogEntry("error",fp2str(error())); - endLogGroup("MultiPhaseEquil::equilibrate"); + if (loglevel > 0) { + addLogEntry("iterations",int2str(iterations())); + addLogEntry("error tolerance",fp2str(err)); + addLogEntry("error",fp2str(error())); + endLogGroup("MultiPhaseEquil::equilibrate"); + } finish(); return error(); } @@ -256,11 +265,12 @@ namespace Cantera { /// to solving the linear programming problem of minimizing the /// linear Gibbs function subject to the element and /// non-negativity constraints. - int MultiPhaseEquil::setInitialMoles() { + int MultiPhaseEquil::setInitialMoles(int loglevel) { index_t ik, j; double not_mu = 1.0e12; - beginLogGroup("MultiPhaseEquil::setInitialMoles"); + if (loglevel > 0) + beginLogGroup("MultiPhaseEquil::setInitialMoles"); m_mix->getValidChemPotentials(not_mu, DATA_PTR(m_mu), true); doublereal dg_rt; @@ -275,7 +285,8 @@ namespace Cantera { // choose a set of components based on the current // composition computeN(); - addLogEntry("iteration",iter); + if (loglevel > 0) + addLogEntry("iteration",iter); redo = false; iter++; if (iter > 4) break; @@ -300,7 +311,8 @@ namespace Cantera { // if a component has nearly zero moles, redo // with a new set of components if (!redo && delta_xi < 1.0e-10 && ik < m_nel) { - addLogEntry("component too small",speciesName(ik)); + if (loglevel > 0) + addLogEntry("component too small",speciesName(ik)); redo = true; } if (delta_xi < dxi_min) dxi_min = delta_xi; @@ -316,8 +328,8 @@ namespace Cantera { } for (ik = 0; ik < m_nsp; ik++) if (moles(ik) != 0.0) addLogEntry(speciesName(ik), moles(ik)); - - endLogGroup("MultiPhaseEquil::setInitialMoles"); + if (loglevel > 0) + endLogGroup("MultiPhaseEquil::setInitialMoles"); return 0; } @@ -477,30 +489,41 @@ namespace Cantera { } #if defined(WITH_HTML_LOGS) - void MultiPhaseEquil::printInfo() { + void MultiPhaseEquil::printInfo(int loglevel) { index_t m, ik, k; - beginLogGroup("info"); - beginLogGroup("components"); + if (loglevel > 0) { + beginLogGroup("info"); + beginLogGroup("components"); + } for (m = 0; m < m_nel; m++) { ik = m_order[m]; k = m_species[ik]; - addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik])); + if (loglevel > 0) + addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik])); + } + if (loglevel > 0) { + endLogGroup("components"); + beginLogGroup("non-components"); } - endLogGroup("components"); - beginLogGroup("non-components"); for (m = m_nel; m < m_nsp; m++) { ik = m_order[m]; k = m_species[ik]; - addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik])); + if (loglevel > 0) + addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik])); + } + if (loglevel > 0) { + endLogGroup("non-components"); + addLogEntry("Error",fp2str(error())); + beginLogGroup("Delta G / RT"); } - endLogGroup("non-components"); - addLogEntry("Error",fp2str(error())); - beginLogGroup("Delta G / RT"); for (k = 0; k < m_nsp - m_nel; k++) { - addLogEntry(reactionString(k), fp2str(m_deltaG_RT[k])); + if (loglevel > 0) + addLogEntry(reactionString(k), fp2str(m_deltaG_RT[k])); + } + if (loglevel > 0) { + endLogGroup("Delta G / RT"); + endLogGroup("info"); } - endLogGroup("Delta G / RT"); - endLogGroup("info"); } /// Return a string specifying the jth reaction. @@ -526,19 +549,23 @@ namespace Cantera { } #endif - void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN) { + void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN, + int loglevel) { index_t k, ik; - beginLogGroup("MultiPhaseEquil::step"); + if (loglevel > 0) + beginLogGroup("MultiPhaseEquil::step"); if (omega < 0.0) throw CanteraError("step","negative omega"); for (ik = 0; ik < m_nel; ik++) { k = m_order[ik]; m_lastmoles[k] = m_moles[k]; - addLogEntry("component "+m_mix->speciesName(m_species[k])+" moles", + if (loglevel > 0) { + addLogEntry("component "+m_mix->speciesName(m_species[k])+" moles", m_moles[k]); - addLogEntry("component "+m_mix->speciesName(m_species[k])+" step", + addLogEntry("component "+m_mix->speciesName(m_species[k])+" step", omega*deltaN[k]); + } m_moles[k] += omega * deltaN[k]; } @@ -554,16 +581,17 @@ namespace Cantera { } } updateMixMoles(); - endLogGroup("MultiPhaseEquil::step"); + if (loglevel > 0) + endLogGroup("MultiPhaseEquil::step"); } /// Take one step in composition, given the gradient of G at the /// starting point, and a vector of reaction steps dxi. doublereal MultiPhaseEquil:: - stepComposition() { - - beginLogGroup("MultiPhaseEquil::stepComposition"); + stepComposition(int loglevel) { + if (loglevel > 0) + beginLogGroup("MultiPhaseEquil::stepComposition"); m_iter++; index_t ik, k = 0; @@ -619,15 +647,17 @@ namespace Cantera { } } if (m_moles[k] < -Tiny) { - addLogEntry("Negative moles for " - +m_mix->speciesName(m_species[k]), fp2str(m_moles[k])); + if (loglevel > 0) + addLogEntry("Negative moles for " + +m_mix->speciesName(m_species[k]), fp2str(m_moles[k])); } m_majorsp[k] = true; } } // now take a step with this scaled omega - addLogEntry("Stepping by ", fp2str(omegamax)); + if (loglevel > 0) + addLogEntry("Stepping by ", fp2str(omegamax)); step(omegamax, m_work); // compute the gradient of G at this new position in the // current direction. If it is positive, then we have overshot @@ -643,11 +673,13 @@ namespace Cantera { if (grad1 > 0.0) { omega *= fabs(grad0) / (grad1 + fabs(grad0)); for (k = 0; k < m_nsp; k++) m_moles[k] = m_lastmoles[k]; - addLogEntry("Stepped over minimum. Take smaller step ", fp2str(omega)); + if (loglevel > 0) + addLogEntry("Stepped over minimum. Take smaller step ", fp2str(omega)); step(omega, m_work); } - printInfo(); - endLogGroup("MultiPhaseEquil::stepComposition"); + printInfo(loglevel); + if (loglevel > 0) + endLogGroup("MultiPhaseEquil::stepComposition"); return omega; } diff --git a/Cantera/src/equil/MultiPhaseEquil.h b/Cantera/src/equil/MultiPhaseEquil.h index 75984935c..52f4fba4e 100644 --- a/Cantera/src/equil/MultiPhaseEquil.h +++ b/Cantera/src/equil/MultiPhaseEquil.h @@ -36,7 +36,7 @@ namespace Cantera { typedef size_t index_t; typedef DenseMatrix matrix_t; - MultiPhaseEquil(mix_t* mix, bool start=true); + MultiPhaseEquil(mix_t* mix, bool start=true, int loglevel = 0); virtual ~MultiPhaseEquil() {} @@ -62,14 +62,14 @@ namespace Cantera { #if defined(WITH_HTML_LOGS) std::string reactionString(index_t j); - void printInfo(); + void printInfo(int loglevel); #else inline std::string reactionString(index_t j) { return std::string(""); } - inline void printInfo() {} + inline void printInfo(int loglevel) {} #endif - void setInitialMixMoles() { - setInitialMoles(); + void setInitialMixMoles(int loglevel = 0) { + setInitialMoles(loglevel); finish(); } @@ -78,12 +78,12 @@ namespace Cantera { protected: void getComponents(const vector_int& order); - int setInitialMoles(); + int setInitialMoles(int loglevel = 0); void computeN(); - doublereal stepComposition(); + doublereal stepComposition(int loglevel = 0); //void sort(vector_fp& x); void unsort(vector_fp& x); - void step(doublereal omega, vector_fp& deltaN); + void step(doublereal omega, vector_fp& deltaN, int loglevel = 0); doublereal computeReactionSteps(vector_fp& dxi); void updateMixMoles(); void finish(); diff --git a/Cantera/src/equil/equilibrate.cpp b/Cantera/src/equil/equilibrate.cpp index e0769d809..55274725d 100644 --- a/Cantera/src/equil/equilibrate.cpp +++ b/Cantera/src/equil/equilibrate.cpp @@ -23,35 +23,42 @@ namespace Cantera { doublereal equilibrate(MultiPhase& s, const char* XY, doublereal tol, int maxsteps, int maxiter, int loglevel) { + if (loglevel > 0) { + beginLogGroup("equilibrate",loglevel); + addLogEntry("multiphase equilibrate function"); + beginLogGroup("arguments"); + addLogEntry("XY",XY); + addLogEntry("tol",tol); + addLogEntry("maxsteps",maxsteps); + addLogEntry("maxiter",maxiter); + addLogEntry("loglevel",loglevel); + endLogGroup("arguments"); + } + s.init(); + int ixy = _equilflag(XY); + 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"); - beginLogGroup("equilibrate",loglevel); - addLogEntry("multiphase equilibrate function"); - beginLogGroup("arguments"); - addLogEntry("XY",XY); - addLogEntry("tol",tol); - addLogEntry("maxsteps",maxsteps); - addLogEntry("maxiter",maxiter); - addLogEntry("loglevel",loglevel); - endLogGroup("arguments"); - - s.init(); - int ixy = _equilflag(XY); - if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) { - try { - double err = s.equilibrate(ixy, tol, maxsteps, maxiter); - addLogEntry("Success. Error",err); - endLogGroup("equilibrate"); - return err; + } + return err; + } + catch (CanteraError e) { + if (loglevel > 0) { + addLogEntry("Failure.",lastErrorMessage()); + endLogGroup("equilibrate"); + } + throw e; + } } - catch (CanteraError e) { - addLogEntry("Failure.",lastErrorMessage()); - endLogGroup("equilibrate"); - throw e; - } - } - else { - addLogEntry("multiphase equilibrium can be done only for TP, HP, SP, or TV"); - endLogGroup("equilibrate"); + 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; } @@ -97,18 +104,20 @@ namespace Cantera { int nAttempts = 0; int retnSub = 0; - 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 (loglevel > 0) { + 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"); + } } while (redo) { if (solver > 0) { @@ -117,22 +126,27 @@ namespace Cantera { m->addPhase(&s, 1.0); m->init(); nAttempts++; - (void) equilibrate(*m, XY, rtol, maxsteps, maxiter, loglevel); + (void) equilibrate(*m, XY, rtol, maxsteps, maxiter, loglevel-1); redo = false; - addLogEntry("MultiPhaseEquil solver succeeded."); + if (loglevel > 0) + addLogEntry("MultiPhaseEquil solver succeeded."); delete m; retn = nAttempts; } catch (CanteraError err) { - addLogEntry("MultiPhaseEquil solver failed."); + if (loglevel > 0) + addLogEntry("MultiPhaseEquil solver failed."); delete m; if (nAttempts < 2) { - addLogEntry("Trying single phase ChemEquil solver."); + if (loglevel > 0) + addLogEntry("Trying single phase ChemEquil solver."); solver = -1; - } else { - endLogGroup("equilibrate"); - throw err; - } + } + else { + if (loglevel > 0) + endLogGroup("equilibrate"); + throw err; + } } } else { // solver <= 0 @@ -144,12 +158,14 @@ namespace Cantera { e->options.maxIterations = maxsteps; e->options.relTolerance = rtol; nAttempts++; - retnSub = e->equilibrate(s,XY); + retnSub = e->equilibrate(s,XY,loglevel-1); if (retnSub < 0) { - addLogEntry("ChemEquil solver failed."); + if (loglevel > 0) + addLogEntry("ChemEquil solver failed."); if (nAttempts < 2) { - addLogEntry("Trying MultiPhaseEquil solver."); - solver = 1; + if (loglevel > 0) + addLogEntry("Trying MultiPhaseEquil solver."); + solver = 1; } else { throw CanteraError("equilibrate", "Both equilibrium solvers failed"); @@ -159,20 +175,24 @@ namespace Cantera { s.setElementPotentials(e->elementPotentials()); redo = false; delete e; - addLogEntry("ChemEquil solver succeeded."); + if (loglevel > 0) + addLogEntry("ChemEquil solver succeeded."); } catch (CanteraError err) { delete e; - addLogEntry("ChemEquil solver failed."); + if (loglevel > 0) + addLogEntry("ChemEquil solver failed."); // If ChemEquil fails, try the MultiPhase solver if (solver < 0) { - addLogEntry("Trying MultiPhaseEquil solver."); + if (loglevel > 0) + addLogEntry("Trying MultiPhaseEquil solver."); solver = 1; } else { redo = false; - endLogGroup("equilibrate"); + if (loglevel > 0) + endLogGroup("equilibrate"); throw err; } } @@ -181,7 +201,8 @@ namespace Cantera { /* * We are here only for a success */ - endLogGroup("equilibrate"); + if (loglevel > 0) + endLogGroup("equilibrate"); return retn; } } diff --git a/Cantera/src/kinetics/FalloffFactory.cpp b/Cantera/src/kinetics/FalloffFactory.cpp index f42f5b6ba..002b710d8 100755 --- a/Cantera/src/kinetics/FalloffFactory.cpp +++ b/Cantera/src/kinetics/FalloffFactory.cpp @@ -21,6 +21,9 @@ namespace Cantera { FalloffFactory* FalloffFactory::s_factory = 0; + #if defined(THREAD_SAFE_CANTERA) + boost::mutex FalloffFactory::falloff_mutex ; + #endif /** * The 3-parameter Troe falloff parameterization. @@ -286,3 +289,4 @@ namespace Cantera { } } + diff --git a/Cantera/src/kinetics/FalloffFactory.h b/Cantera/src/kinetics/FalloffFactory.h index f014df598..adfcd53ba 100755 --- a/Cantera/src/kinetics/FalloffFactory.h +++ b/Cantera/src/kinetics/FalloffFactory.h @@ -20,6 +20,11 @@ #include "ct_defs.h" #include "reaction_defs.h" +#include "FactoryBase.h" + +#if defined(THREAD_SAFE_CANTERA) +#include +#endif namespace Cantera { @@ -87,7 +92,7 @@ namespace Cantera { * @endcode * @ingroup falloffGroup */ - class FalloffFactory { + class FalloffFactory : public FactoryBase { public: /** @@ -97,23 +102,29 @@ namespace Cantera { * to the existing factory is returned. */ static FalloffFactory* factory() { + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(falloff_mutex) ; + #endif if (!s_factory) s_factory = new FalloffFactory; return s_factory; } - static void deleteFalloffFactory() { - if (s_factory) { - delete s_factory; - s_factory = 0; - } - } + virtual void deleteFactory() { + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(falloff_mutex) ; + #endif + if (s_factory) { + delete s_factory; + s_factory = 0; + } + } /** * Destructor doesn't do anything. We do not delete statically - * created single instance of this class here, because it would - * create an infinite loop if destructor is called for that - * single instance. Instead, to delete single instance, we - * call delete[] from FalloffMng's destructor. + * created single instance of this class here, because it would + * create an infinite loop if destructor is called for that + * single instance. Instead, to delete single instance, we + * call delete[] from FalloffMng's destructor. */ virtual ~FalloffFactory() { } @@ -129,8 +140,12 @@ namespace Cantera { private: static FalloffFactory* s_factory; FalloffFactory(){} + #if defined(THREAD_SAFE_CANTERA) + static boost::mutex falloff_mutex ; + #endif }; } #endif + diff --git a/Cantera/src/kinetics/FalloffMgr.h b/Cantera/src/kinetics/FalloffMgr.h index 4b8c82026..a36d7bcce 100755 --- a/Cantera/src/kinetics/FalloffMgr.h +++ b/Cantera/src/kinetics/FalloffMgr.h @@ -30,10 +30,12 @@ namespace Cantera { * falloff function calculators. If omitted, the standard factory * will be used. */ - FalloffMgr(FalloffFactory* f = 0) : - m_n(0), m_n0(0), m_worksize(0) { - if (f == 0) m_factory = FalloffFactory::factory(); - else m_factory = f; + FalloffMgr(/*FalloffFactory* f = 0*/) : + m_n(0), m_n0(0), m_worksize(0) { + //if (f == 0) + m_factory = FalloffFactory::factory(); // RFB:TODO This raw pointer should be encapsulated + // because accessing a 'Singleton Factory' + //else m_factory = f; } /** @@ -43,10 +45,10 @@ namespace Cantera { virtual ~FalloffMgr(){ int i; for (i = 0; i < m_n; i++) delete m_falloff[i]; - if (m_factory) { - FalloffFactory::deleteFalloffFactory(); - m_factory = 0; - } + //if (m_factory) { + //FalloffFactory::deleteFalloffFactory(); + //m_factory = 0; + //} } /** @@ -58,7 +60,7 @@ namespace Cantera { * @param c vector of coefficients for the falloff function. */ void install(int rxn, int type, - const vector_fp& c) { + const vector_fp& c) { if (type != SIMPLE_FALLOFF) { m_rxn.push_back(rxn); Falloff* f = m_factory->newFalloff(type,c); @@ -121,3 +123,4 @@ namespace Cantera { } #endif + diff --git a/Cantera/src/kinetics/KineticsFactory.cpp b/Cantera/src/kinetics/KineticsFactory.cpp index 8e66fdfe9..770eef374 100644 --- a/Cantera/src/kinetics/KineticsFactory.cpp +++ b/Cantera/src/kinetics/KineticsFactory.cpp @@ -28,6 +28,9 @@ using namespace std; namespace Cantera { KineticsFactory* KineticsFactory::s_factory = 0; + #if defined(THREAD_SAFE_CANTERA) + boost::mutex KineticsFactory::kinetics_mutex ; + #endif static int ntypes = 5; static string _types[] = {"none", "GasKinetics", "GRI30", "Interface", "Edge"}; @@ -58,28 +61,28 @@ namespace Cantera { Kinetics* KineticsFactory:: newKinetics(XML_Node& phaseData, vector th) { - /* - * Look for a child of the xml element phase called - * "kinetics". It has an attribute name "model". - * Store the value of that attribute in the variable kintype - */ + /* + * Look for a child of the xml element phase called + * "kinetics". It has an attribute name "model". + * Store the value of that attribute in the variable kintype + */ string kintype = phaseData.child("kinetics")["model"]; - /* - * look up the string kintype in the list of known - * kinetics managers (list is kept at the top of this file). - * Translate it to an integer value, ikin. - */ + /* + * look up the string kintype in the list of known + * kinetics managers (list is kept at the top of this file). + * Translate it to an integer value, ikin. + */ int ikin=-1; int n; for (n = 0; n < ntypes; n++) { - if (kintype == _types[n]) ikin = _itypes[n]; + if (kintype == _types[n]) ikin = _itypes[n]; } - /* - * Assign the kinetics manager based on the value of ikin. - * Kinetics managers are classes derived from the base - * Kinetics class. Unknown kinetics managers will throw a - * CanteraError here. - */ + /* + * Assign the kinetics manager based on the value of ikin. + * Kinetics managers are classes derived from the base + * Kinetics class. Unknown kinetics managers will throw a + * CanteraError here. + */ Kinetics* k=0; switch (ikin) { @@ -104,8 +107,8 @@ namespace Cantera { break; default: - throw UnknownKineticsModel("KineticsFactory::newKinetics", - kintype); + throw UnknownKineticsModel("KineticsFactory::newKinetics", + kintype); } // Now that we have the kinetics manager, we can @@ -143,10 +146,11 @@ namespace Cantera { break; default: - throw UnknownKineticsModel("KineticsFactory::newKinetics", - model); + throw UnknownKineticsModel("KineticsFactory::newKinetics", + model); } return k; } } + diff --git a/Cantera/src/kinetics/KineticsFactory.h b/Cantera/src/kinetics/KineticsFactory.h index 9270dbc1a..488b9255b 100644 --- a/Cantera/src/kinetics/KineticsFactory.h +++ b/Cantera/src/kinetics/KineticsFactory.h @@ -16,35 +16,53 @@ #include "Kinetics.h" #include "xml.h" +#include "FactoryBase.h" + +#if defined(THREAD_SAFE_CANTERA) +#include +#endif namespace Cantera { class UnknownKineticsModel : public CanteraError { public: - UnknownKineticsModel(std::string proc, std::string kineticsModel) : - CanteraError(proc, "Specified Kinetics model " - + kineticsModel + - " does not match any known type.") {} - virtual ~UnknownKineticsModel() {} + UnknownKineticsModel(std::string proc, std::string kineticsModel) : + CanteraError(proc, "Specified Kinetics model " + + kineticsModel + + " does not match any known type.") {} + virtual ~UnknownKineticsModel() {} }; /** * Factory for kinetics managers. */ - class KineticsFactory { + class KineticsFactory : public FactoryBase { public: static KineticsFactory* factory() { + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(kinetics_mutex) ; + #endif if (!s_factory) s_factory = new KineticsFactory; return s_factory; } virtual ~KineticsFactory() { - delete s_factory; - s_factory = 0; + //delete s_factory; + //s_factory = 0; + } + + virtual void deleteFactory() { + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(kinetics_mutex) ; + #endif + if ( s_factory ) { + delete s_factory ; + s_factory = 0 ; + } } /** @@ -59,6 +77,9 @@ namespace Cantera { static KineticsFactory* s_factory; KineticsFactory(){} + #if defined(THREAD_SAFE_CANTERA) + static boost::mutex kinetics_mutex ; + #endif }; @@ -89,3 +110,4 @@ namespace Cantera { #endif + diff --git a/Cantera/src/kinetics/importKinetics.cpp b/Cantera/src/kinetics/importKinetics.cpp index d4d387541..3ca21cbd5 100644 --- a/Cantera/src/kinetics/importKinetics.cpp +++ b/Cantera/src/kinetics/importKinetics.cpp @@ -25,6 +25,7 @@ #include "importKinetics.h" #include "mix_defs.h" #include +#include // Cantera includes #include "speciesThermoTypes.h" @@ -45,36 +46,47 @@ using namespace ctml; using namespace std; +namespace Cantera { + //! these are all used to check for duplicate reactions class rxninfo { public: //! rdata - std::vector< std::map > rdata; + std::vector< std::map > m_rdata; //! string name - std::vector eqn; + std::vector m_eqn; //! string vector of ints - std::vector dup; + std::vector m_dup; //! string vector of ints - std::vector nr; + std::vector m_nr; //! string vector of ints - std::vector typ; + std::vector m_typ; //! vector of bools. - std::vector rev; + std::vector m_rev; + ~rxninfo() { + m_eqn.clear(); + m_dup.clear(); + m_nr.clear(); + m_typ.clear(); + m_rdata.clear(); + } + bool installReaction(int i, const XML_Node& r, Kinetics* k, + std::string default_phase, int rule, + bool validate_rxn) ; }; //! Temporary storage of rxninfo -rxninfo* _rxns = 0; +//rxninfo* _rxns = 0; //! @name utilitydefines. //@{ -#define _reactiondata _rxns->rdata -#define _eqn _rxns->eqn -#define _dup _rxns->dup -#define _nr _rxns->nr -#define _typ _rxns->typ -#define _rev _rxns->rev +#define _reactiondata m_rdata +#define _eqn m_eqn +#define _dup m_dup +#define _nr m_nr +#define _typ m_typ +#define _rev m_rev //@} -namespace Cantera { /* * First we define a couple of typedefs that will @@ -588,7 +600,7 @@ namespace Cantera { * * @ingroup kineticsmgr */ - static bool installReaction(int i, const XML_Node& r, Kinetics* k, + bool rxninfo::installReaction(int i, const XML_Node& r, Kinetics* k, string default_phase, int rule, bool validate_rxn) { @@ -732,12 +744,12 @@ namespace Cantera { string msg = string("Undeclared duplicate reactions detected: \n") +"Reaction "+int2str(nn+1)+": "+_eqn[nn] +"\nReaction "+int2str(i+1)+": "+eqn+"\n"; - _reactiondata.clear(); - _eqn.clear(); - _rev.clear(); - _nr.clear(); - _typ.clear(); - _dup.clear(); + //_reactiondata.clear(); + //_eqn.clear(); + //_rev.clear(); + //_nr.clear(); + //_typ.clear(); + //_dup.clear(); throw CanteraError("installReaction",msg); } } @@ -795,15 +807,16 @@ namespace Cantera { bool installReactionArrays(const XML_Node& p, Kinetics& kin, std::string default_phase, bool check_for_duplicates) { - if (_rxns == 0) { - _rxns = new rxninfo; - } - _eqn.clear(); - _dup.clear(); - _nr.clear(); - _typ.clear(); - _reactiondata.clear(); - _rev.clear(); + std::auto_ptr< rxninfo > _rxns( new rxninfo ) ; + //if (_rxns == 0) { + //_rxns = new rxninfo; + //} + //_eqn.clear(); + //_dup.clear(); + //_nr.clear(); + //_typ.clear(); + //_reactiondata.clear(); + //_rev.clear(); vector rarrays; int itot = 0; @@ -870,7 +883,7 @@ namespace Cantera { for (i = 0; i < nrxns; i++) { const XML_Node* r = allrxns[i]; if (r) { - if (installReaction(itot, *r, &kin, + if (_rxns->installReaction(itot, *r, &kin, default_phase, rxnrule, check_for_duplicates)) ++itot; } } @@ -904,7 +917,7 @@ namespace Cantera { * sometimes has surprising results. */ if ((rxid >= imin) && (rxid <= imax)) { - if (installReaction(itot, *r, &kin, + if (_rxns->installReaction(itot, *r, &kin, default_phase, rxnrule, check_for_duplicates)) ++itot; } } @@ -919,13 +932,13 @@ namespace Cantera { */ kin.finalize(); //writer = 0; - _eqn.clear(); - _dup.clear(); - _nr.clear(); - _typ.clear(); - _reactiondata.clear(); - delete _rxns; - _rxns = 0; + //_eqn.clear(); + //_dup.clear(); + //_nr.clear(); + //_typ.clear(); + //_reactiondata.clear(); + //delete _rxns; + //_rxns = 0; return true; } diff --git a/Cantera/src/thermo/SpeciesThermoFactory.cpp b/Cantera/src/thermo/SpeciesThermoFactory.cpp index ccc7eac47..25279635e 100755 --- a/Cantera/src/thermo/SpeciesThermoFactory.cpp +++ b/Cantera/src/thermo/SpeciesThermoFactory.cpp @@ -36,6 +36,10 @@ using namespace ctml; namespace Cantera { SpeciesThermoFactory* SpeciesThermoFactory::s_factory = 0; + #if defined(THREAD_SAFE_CANTERA) + boost::mutex SpeciesThermoFactory::species_thermo_mutex ; + #endif + /** * Examine the types of species thermo parameterizations, @@ -48,8 +52,8 @@ namespace Cantera { * @todo Make sure that spDadta_node is species Data XML node by checking its name is speciesData */ static void getSpeciesThermoTypes(XML_Node* spData_node, - int& has_nasa, int& has_shomate, int& has_simple, - int &has_other) { + int& has_nasa, int& has_shomate, int& has_simple, + int &has_other) { const XML_Node& sparray = *spData_node; std::vector sp; @@ -59,19 +63,19 @@ namespace Cantera { for (n = 0; n < ns; n++) { XML_Node* spNode = sp[n]; if (spNode->hasChild("thermo")) { - const XML_Node& th = sp[n]->child("thermo"); - if (th.hasChild("NASA")) has_nasa = 1; - if (th.hasChild("Shomate")) has_shomate = 1; - if (th.hasChild("const_cp")) has_simple = 1; - if (th.hasChild("poly")) { - if (th.child("poly")["order"] == "1") has_simple = 1; - else throw CanteraError("newSpeciesThermo", - "poly with order > 1 not yet supported"); - } - if (th.hasChild("Mu0")) has_other = 1; + const XML_Node& th = sp[n]->child("thermo"); + if (th.hasChild("NASA")) has_nasa = 1; + if (th.hasChild("Shomate")) has_shomate = 1; + if (th.hasChild("const_cp")) has_simple = 1; + if (th.hasChild("poly")) { + if (th.child("poly")["order"] == "1") has_simple = 1; + else throw CanteraError("newSpeciesThermo", + "poly with order > 1 not yet supported"); + } + if (th.hasChild("Mu0")) has_other = 1; } else { - throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:", - spNode->attrib("name"), "missing"); + throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:", + spNode->attrib("name"), "missing"); } } } @@ -94,7 +98,7 @@ namespace Cantera { return new GeneralSpeciesThermo(); } return newSpeciesThermo(NASA*inasa - + SHOMATE*ishomate + SIMPLE*isimple); + + SHOMATE*ishomate + SIMPLE*isimple); } SpeciesThermo* SpeciesThermoFactory:: @@ -103,17 +107,17 @@ namespace Cantera { int inasa = 0, ishomate = 0, isimple = 0, iother = 0; for (int j = 0; j < n; j++) { try { - getSpeciesThermoTypes(spData_nodes[j], inasa, ishomate, isimple, iother); + getSpeciesThermoTypes(spData_nodes[j], inasa, ishomate, isimple, iother); } catch (UnknownSpeciesThermoModel) { - iother = 1; - popError(); + iother = 1; + popError(); } } if (iother) { return new GeneralSpeciesThermo(); } return newSpeciesThermo(NASA*inasa - + SHOMATE*ishomate + SIMPLE*isimple); + + SHOMATE*ishomate + SIMPLE*isimple); } @@ -126,17 +130,17 @@ namespace Cantera { int inasa = 0, ishomate = 0, isimple = 0, iother = 0; for (int j = 0; j < n; j++) { try { - getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple, iother); + getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple, iother); } catch (UnknownSpeciesThermoModel) { - iother = 1; - popError(); + iother = 1; + popError(); } } if (iother) { return new GeneralSpeciesThermo(); } return newSpeciesThermo(NASA*inasa - + SHOMATE*ishomate + SIMPLE*isimple); + + SHOMATE*ishomate + SIMPLE*isimple); } @@ -158,7 +162,7 @@ namespace Cantera { return new SpeciesThermoDuo; default: throw UnknownSpeciesThermo( - "SpeciesThermoFactory::newSpeciesThermo",type); + "SpeciesThermoFactory::newSpeciesThermo",type); return 0; } } @@ -169,7 +173,7 @@ namespace Cantera { * temperature. */ void NasaThermo::checkContinuity(std::string name, double tmid, const doublereal* clow, - doublereal* chigh) { + doublereal* chigh) { // heat capacity doublereal cplow = poly4(tmid, clow); @@ -269,7 +273,7 @@ namespace Cantera { } else { throw CanteraError("installNasaThermo", - "non-continuous temperature ranges."); + "non-continuous temperature ranges."); } // The NasaThermo species property manager expects the @@ -331,7 +335,7 @@ namespace Cantera { } else { throw CanteraError("installNasaThermo", - "non-continuous temperature ranges."); + "non-continuous temperature ranges."); } array_fp c(15); c[0] = tmid; @@ -378,9 +382,9 @@ namespace Cantera { if (dualRange) getFloatArray(f1ptr->child("floatArray"), c1, false); else { - c1.resize(7,0.0); - copy(c0.begin(), c0.begin()+7, c1.begin()); - } + c1.resize(7,0.0); + copy(c0.begin(), c0.begin()+7, c1.begin()); + } } else if (fabs(tmax1 - tmin0) < 0.01) { tmin = tmin1; @@ -391,7 +395,7 @@ namespace Cantera { } else { throw CanteraError("installShomateThermo", - "non-continuous temperature ranges."); + "non-continuous temperature ranges."); } array_fp c(15); c[0] = tmid; @@ -434,18 +438,18 @@ namespace Cantera { void SpeciesThermoFactory:: installThermoForSpecies(int k, const XML_Node& s, SpeciesThermo& spthermo) { - /* - * Check to see that the species block has a thermo block - * before processing. Throw an error if not there. - */ - if (!(s.hasChild("thermo"))) { - throw UnknownSpeciesThermoModel("installSpecies", - s["name"], ""); - } - const XML_Node& thermo = s.child("thermo"); - const std::vector& tp = thermo.children(); - int nc = static_cast(tp.size()); - if (nc == 1) { + /* + * Check to see that the species block has a thermo block + * before processing. Throw an error if not there. + */ + if (!(s.hasChild("thermo"))) { + throw UnknownSpeciesThermoModel("installSpecies", + s["name"], ""); + } + const XML_Node& thermo = s.child("thermo"); + const std::vector& tp = thermo.children(); + int nc = static_cast(tp.size()); + if (nc == 1) { const XML_Node* f = tp[0]; if (f->name() == "Shomate") { installShomateThermoFromXML(s["name"], spthermo, k, f, 0); @@ -456,15 +460,15 @@ namespace Cantera { else if (f->name() == "NASA") { installNasaThermoFromXML(s["name"], spthermo, k, f, 0); } - else if (f->name() == "Mu0") { - installMu0ThermoFromXML(s["name"], spthermo, k, f); - } + else if (f->name() == "Mu0") { + installMu0ThermoFromXML(s["name"], spthermo, k, f); + } else { throw UnknownSpeciesThermoModel("installSpecies", - s["name"], f->name()); + s["name"], f->name()); } - } - else if (nc == 2) { + } + else if (nc == 2) { const XML_Node* f0 = tp[0]; const XML_Node* f1 = tp[1]; if (f0->name() == "NASA" && f1->name() == "NASA") { @@ -475,14 +479,16 @@ namespace Cantera { } else { throw UnknownSpeciesThermoModel("installSpecies", s["name"], - f0->name() + " and " - + f1->name()); + f0->name() + " and " + + f1->name()); } - } - else { - throw UnknownSpeciesThermoModel("installSpecies", s["name"], - "multiple"); - } + } + else { + throw UnknownSpeciesThermoModel("installSpecies", s["name"], + "multiple"); + } } } + + diff --git a/Cantera/src/thermo/SpeciesThermoFactory.h b/Cantera/src/thermo/SpeciesThermoFactory.h index 576de2259..5e0b4c587 100755 --- a/Cantera/src/thermo/SpeciesThermoFactory.h +++ b/Cantera/src/thermo/SpeciesThermoFactory.h @@ -19,6 +19,8 @@ #include "SpeciesThermo.h" #include "ctexceptions.h" +#include "FactoryBase.h" + namespace Cantera { @@ -59,7 +61,7 @@ namespace Cantera { * * @ingroup thermoprops */ - class SpeciesThermoFactory { + class SpeciesThermoFactory : public FactoryBase { public: @@ -73,6 +75,9 @@ namespace Cantera { * instance. */ static SpeciesThermoFactory* factory() { +#if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(species_thermo_mutex); +#endif if (!s_factory) s_factory = new SpeciesThermoFactory; return s_factory; } @@ -83,7 +88,10 @@ namespace Cantera { * the process terminates (for example, when checking for * memory leaks) then this method can be called to delete it. */ - static void deleteFactory() { + virtual void deleteFactory() { +#if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(species_thermo_mutex); +#endif if (s_factory) { delete s_factory; s_factory = 0; @@ -154,6 +162,10 @@ namespace Cantera { //! pointer to the sole instance of this class static SpeciesThermoFactory* s_factory; +#if defined(THREAD_SAFE_CANTERA) + static boost::mutex species_thermo_mutex; +#endif + //! Constructor. This is made private, so that only the static //! method factory() can instantiate the class. SpeciesThermoFactory(){} diff --git a/Cantera/src/thermo/ThermoFactory.cpp b/Cantera/src/thermo/ThermoFactory.cpp index 705b4e763..a30978b31 100644 --- a/Cantera/src/thermo/ThermoFactory.cpp +++ b/Cantera/src/thermo/ThermoFactory.cpp @@ -62,6 +62,9 @@ using namespace std; namespace Cantera { ThermoFactory* ThermoFactory::s_factory = 0; +#if defined(THREAD_SAFE_CANTERA) + boost::mutex ThermoFactory::thermo_mutex; +#endif static int ntypes = 10; static string _types[] = {"IdealGas", "Incompressible", diff --git a/Cantera/src/thermo/ThermoFactory.h b/Cantera/src/thermo/ThermoFactory.h index fa465c02c..2e5857f2b 100644 --- a/Cantera/src/thermo/ThermoFactory.h +++ b/Cantera/src/thermo/ThermoFactory.h @@ -19,10 +19,17 @@ #include "ThermoPhase.h" #include "xml.h" -#include "SpeciesThermoFactory.h" + +#if defined(THREAD_SAFE_CANTERA) +#include +#endif + +//#include "SpeciesThermoFactory.h" +#include "FactoryBase.h" namespace Cantera { + class SpeciesThermoFactory; /*! * @addtogroup thermoprops @@ -55,33 +62,38 @@ namespace Cantera { //! Factory class for thermodynamic property managers. /*! - * This class keeps a list of the known ThermoPhase classes, and is used - * to create new instances of these classes. + * This class keeps a list of the known ThermoPhase classes, and is + * used to create new instances of these classes. */ - class ThermoFactory { + class ThermoFactory : public FactoryBase { public: - //! Static function that creates a static instance of the factor. + //! Static function that creates a static instance of the factory. static ThermoFactory* factory() { +#if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(thermo_mutex); +#endif if (!s_factory) s_factory = new ThermoFactory; return s_factory; } - //! delete the static instance of this factory - static void deleteFactory() { - if (s_factory) { - delete s_factory; - s_factory = 0; - } + //! delete the static instance of this factory + virtual void deleteFactory() { +#if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(thermo_mutex); +#endif + if (s_factory) { + delete s_factory; + s_factory = 0; + } } //! Destructor doesn't do anything. /*! - * We do not delete statically - * created single instance of this class here, because it would - * create an infinite loop if destructor is called for that - * single instance. + * We do not delete statically created single instance of this + * class here, because it would create an infinite loop if + * destructor is called for that single instance. */ virtual ~ThermoFactory() { } @@ -103,6 +115,11 @@ namespace Cantera { //! Private constructor prevents usage ThermoFactory(){} + +#if defined(THREAD_SAFE_CANTERA) + static boost::mutex thermo_mutex; +#endif + }; //! Create a new thermo manager instance. diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp index 24d3119ff..7803f4a49 100755 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -46,6 +46,9 @@ using namespace std; namespace Cantera { TransportFactory* TransportFactory::s_factory = 0; +#if defined(THREAD_SAFE_CANTERA) + boost::mutex TransportFactory::transport_mutex; +#endif ////////////////////////// exceptions ///////////////////////// @@ -259,20 +262,23 @@ namespace Cantera { * explicitly deleted. */ TransportFactory::~TransportFactory() { - if (m_integrals) { + if (m_integrals) { delete m_integrals; m_integrals = 0; - } } - + } + /** * This static function deletes the statically allocated instance. */ - void TransportFactory::deleteTransportFactory() { - if (s_factory) { - delete s_factory; - s_factory = 0; - } + void TransportFactory::deleteFactory() { + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(transport_mutex) ; + #endif + if (s_factory) { + delete s_factory; + s_factory = 0; + } } /** @@ -527,7 +533,7 @@ namespace Cantera { tr.astar_poly.push_back(ca); tr.bstar_poly.push_back(cb); tr.cstar_poly.push_back(cc); - tr.poly[i][j] = static_cast(tr.astar_poly.size()) - 1; + tr.poly[i][j] = static_cast(tr.astar_poly.size()) - 1; tr.fitlist.push_back(dstar); } @@ -990,3 +996,4 @@ namespace Cantera { tr.xml->XML_close(logfile, "binary_diffusion_coefficients"); } } + diff --git a/Cantera/src/transport/TransportFactory.h b/Cantera/src/transport/TransportFactory.h index 7cf92e087..74277e813 100755 --- a/Cantera/src/transport/TransportFactory.h +++ b/Cantera/src/transport/TransportFactory.h @@ -35,6 +35,11 @@ using namespace std; // Cantera includes #include "ct_defs.h" #include "TransportBase.h" +#include "FactoryBase.h" + +#if defined(THREAD_SAFE_CANTERA) +#include +#endif namespace Cantera { @@ -43,11 +48,11 @@ namespace Cantera { */ struct GasTransportData { GasTransportData() : speciesName("-"), - geometry(-1), wellDepth(-1.0), - diameter(-1.0), - dipoleMoment(-1.0), - polarizability(-1.0), - rotRelaxNumber(-1.0) {} + geometry(-1), wellDepth(-1.0), + diameter(-1.0), + dipoleMoment(-1.0), + polarizability(-1.0), + rotRelaxNumber(-1.0) {} string speciesName; int geometry; @@ -73,7 +78,7 @@ namespace Cantera { * created in other ways. @ingroup transportgroup * @ingroup transportProps */ - class TransportFactory { + class TransportFactory : FactoryBase { public: @@ -91,27 +96,30 @@ namespace Cantera { * @endcode */ static TransportFactory* factory() { - if (!s_factory) { - s_factory = new TransportFactory(); - } + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(transport_mutex) ; + #endif + if (!s_factory) { + s_factory = new TransportFactory(); + } return s_factory; } - /** - * Deletes the statically malloced instance. - */ - static void deleteTransportFactory(); + /** + * Deletes the statically malloced instance. + */ + virtual void deleteFactory(); /** * Destructor - * - * We do not delete statically - * created single instance of this class here, because it would - * create an infinite loop if destructor is called for that - * single instance. + * + * We do not delete statically + * created single instance of this class here, because it would + * create an infinite loop if destructor is called for that + * single instance. */ virtual ~TransportFactory(); - + /// Build a new transport manager virtual Transport* @@ -125,6 +133,9 @@ namespace Cantera { private: static TransportFactory* s_factory; + #if defined(THREAD_SAFE_CANTERA) + static boost::mutex transport_mutex ; + #endif // The constructor is private; use static method factory() to // get a pointer to a factory instance @@ -190,14 +201,14 @@ namespace Cantera { f = TransportFactory::factory(); } Transport* ptr = f->newTransport(transportModel, thermo, loglevel); - /* - * Note: We delete the static s_factory instance here, instead of in - * appdelete() in misc.cpp, to avoid linking problems involving - * the need for multiple cantera and transport library statements - * for applications that don't have transport in them. - */ - TransportFactory::deleteTransportFactory(); - return ptr; + /* + * Note: We delete the static s_factory instance here, instead of in + * appdelete() in misc.cpp, to avoid linking problems involving + * the need for multiple cantera and transport library statements + * for applications that don't have transport in them. + */ + //TransportFactory::deleteFactory(); + return ptr; } @@ -209,3 +220,4 @@ namespace Cantera { + diff --git a/Cantera/src/zeroD/Makefile.in b/Cantera/src/zeroD/Makefile.in index d610b433f..cf8e23de1 100644 --- a/Cantera/src/zeroD/Makefile.in +++ b/Cantera/src/zeroD/Makefile.in @@ -26,7 +26,7 @@ ZEROD_H = Reactor.h ReactorBase.h FlowDevice.h Wall.h ReactorNet.h \ flowControllers.h PID_Controller.h Reservoir.h FlowReactor.h \ ConstPressureReactor.h ReactorFactory.h -CXX_INCLUDES = -I../base -I../thermo -I../kinetics -I../numerics +CXX_INCLUDES = -I../base -I../thermo -I../kinetics -I../numerics @CXX_INCLUDES@ ZEROD_LIB = @buildlib@/libzeroD.a DEPENDS = $(OBJS:.o=.d) diff --git a/Cantera/src/zeroD/ReactorFactory.cpp b/Cantera/src/zeroD/ReactorFactory.cpp index 0bd165a70..ee07a21f3 100644 --- a/Cantera/src/zeroD/ReactorFactory.cpp +++ b/Cantera/src/zeroD/ReactorFactory.cpp @@ -26,6 +26,9 @@ using namespace std; namespace CanteraZeroD { ReactorFactory* ReactorFactory::s_factory = 0; + #ifdef THREAD_SAFE_CANTERA + boost::mutex ReactorFactory::reactor_mutex ; + #endif static int ntypes = 4; static string _types[] = {"Reservoir", "Reactor", "ConstPressureReactor", @@ -67,3 +70,4 @@ namespace CanteraZeroD { } } + diff --git a/Cantera/src/zeroD/ReactorFactory.h b/Cantera/src/zeroD/ReactorFactory.h index 6dae2ad20..d1da8024f 100644 --- a/Cantera/src/zeroD/ReactorFactory.h +++ b/Cantera/src/zeroD/ReactorFactory.h @@ -15,27 +15,38 @@ #define REACTOR_FACTORY_H #include "ReactorBase.h" +#include "FactoryBase.h" + +#if defined(THREAD_SAFE_CANTERA) +#include +#endif namespace CanteraZeroD { - class ReactorFactory { + class ReactorFactory : FactoryBase { public: static ReactorFactory* factory() { + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(reactor_mutex) ; + #endif if (!s_factory) s_factory = new ReactorFactory; return s_factory; } - static void deleteFactory() { - if (s_factory) { - delete s_factory; - s_factory = 0; - } - } + virtual void deleteFactory() { + #if defined(THREAD_SAFE_CANTERA) + boost::mutex::scoped_lock lock(reactor_mutex) ; + #endif + if (s_factory) { + delete s_factory; + s_factory = 0; + } + } - /** + /** * Destructor doesn't do anything. */ virtual ~ReactorFactory() {} @@ -50,6 +61,9 @@ namespace CanteraZeroD { private: static ReactorFactory* s_factory; + #if defined(THREAD_SAFE_CANTERA) + static boost::mutex reactor_mutex ; + #endif ReactorFactory(){} }; @@ -66,3 +80,4 @@ namespace CanteraZeroD { #endif + diff --git a/config.h.in b/config.h.in index 75ca1cf1a..7d8383107 100755 --- a/config.h.in +++ b/config.h.in @@ -100,7 +100,7 @@ typedef int ftnlen; // Fortran hidden string length type #undef CANTERA_DATA -#define WITH_HTML_LOGS +#undef WITH_HTML_LOGS // define STORE_MOLE_FRACTIONS if you want Cantera to internally // represent the composition of a mixture as mole fractions. Usually @@ -116,6 +116,8 @@ typedef int ftnlen; // Fortran hidden string length type //--------------------- compile options ---------------------------- #undef USE_PCH +#define THREAD_SAFE_CANTERA + //--------------------- optional phase models ---------------------- // This define indicates the enabling of the inclusion of // accurate liquid/vapor equations diff --git a/configure b/configure index 4f450ebcd..1203723d2 100755 --- a/configure +++ b/configure @@ -308,7 +308,7 @@ ac_includes_default="\ # include #endif" -ac_subst_vars='SHELL PATH_SEPARATOR PACKAGE_NAME PACKAGE_TARNAME PACKAGE_VERSION PACKAGE_STRING PACKAGE_BUGREPORT exec_prefix prefix program_transform_name bindir sbindir libexecdir datadir sysconfdir sharedstatedir localstatedir libdir includedir oldincludedir infodir mandir build_alias host_alias target_alias DEFS ECHO_C ECHO_N ECHO_T LIBS BITCOMPILE BITHARDWARE BITCHANGE ldemulationarg CVF_LIBDIR USE_CLIB_DLL local_inst local_python_inst python_prefix python_win_prefix ctversion homedir ct_libdir ct_bindir ct_incdir ct_incroot ct_datadir ct_demodir ct_templdir ct_tutdir ct_docdir ct_dir ct_mandir COMPACT_INSTALL build build_cpu build_vendor build_os host host_cpu host_vendor host_os target target_cpu target_vendor target_os username ctroot buildinc buildlib buildbin MAKE GRAPHVIZDIR ARCHIVE DO_RANLIB RANLIB CXX_DEPENDS USERDIR INCL_USER_CODE CXX CXXFLAGS LDFLAGS CPPFLAGS ac_ct_CXX EXEEXT OBJEXT use_sundials CVODE_LIBS IDA_LIBS sundials_include CANTERA_DEBUG_MODE COMPILE_PURE_FLUIDS phase_object_files phase_header_files COMPILE_IDEAL_SOLUTIONS COMPILE_ELECTROLYTES NEED_CATHERMO COMPILE_KINETICS COMPILE_HETEROKIN COMPILE_RXNPATH WITH_REACTORS KERNEL KERNEL_OBJ BUILD_CK LIB_DIR build_lapack build_blas BLAS_LAPACK_LIBS BLAS_LAPACK_DIR build_with_f2c build_f2c_lib LOCAL_LIB_DIRS LOCAL_LIBS CT_SHARED_LIB PYTHON_CMD BUILD_PYTHON NUMARRAY_INC_DIR NUMARRAY_HOME CANTERA_PYTHON_HOME CVSTAG MATLAB_CMD BUILD_MATLAB BUILD_CLIB export_name INSTALL_PROGRAM INSTALL_SCRIPT INSTALL_DATA CC CFLAGS ac_ct_CC CXXCPP EGREP SOEXT SHARED PIC CXX_INCLUDES LCXX_FLAGS LCXX_END_LIBS HAVE_STRIPSYMBOLS F77 FFLAGS ac_ct_F77 F90 BUILD_F90 F90FLAGS F90BUILDFLAGS precompile_headers FLIBS OS_IS_DARWIN OS_IS_WIN OS_IS_CYGWIN SHARED_CTLIB mex_ext F77_EXT CXX_EXT OBJ_EXT EXE_EXT math_libs SO LDSHARED EXTRA_LINK TSCOMPARE_abs INSTALL_abs INSTALL_VERBOSE LIBOBJS LTLIBOBJS' +ac_subst_vars='SHELL PATH_SEPARATOR PACKAGE_NAME PACKAGE_TARNAME PACKAGE_VERSION PACKAGE_STRING PACKAGE_BUGREPORT exec_prefix prefix program_transform_name bindir sbindir libexecdir datadir sysconfdir sharedstatedir localstatedir libdir includedir oldincludedir infodir mandir build_alias host_alias target_alias DEFS ECHO_C ECHO_N ECHO_T LIBS BITCOMPILE BITHARDWARE BITCHANGE ldemulationarg CVF_LIBDIR USE_CLIB_DLL local_inst local_python_inst python_prefix python_win_prefix ctversion homedir ct_libdir ct_bindir ct_incdir ct_incroot ct_datadir ct_demodir ct_templdir ct_tutdir ct_docdir ct_dir ct_mandir COMPACT_INSTALL build build_cpu build_vendor build_os host host_cpu host_vendor host_os target target_cpu target_vendor target_os username ctroot buildinc buildlib buildbin MAKE GRAPHVIZDIR ARCHIVE DO_RANLIB RANLIB CXX_DEPENDS USERDIR INCL_USER_CODE CXX CXXFLAGS LDFLAGS CPPFLAGS ac_ct_CXX EXEEXT OBJEXT use_sundials CVODE_LIBS IDA_LIBS sundials_include CANTERA_DEBUG_MODE COMPILE_PURE_FLUIDS phase_object_files phase_header_files COMPILE_IDEAL_SOLUTIONS COMPILE_ELECTROLYTES NEED_CATHERMO COMPILE_KINETICS COMPILE_HETEROKIN COMPILE_RXNPATH WITH_REACTORS KERNEL KERNEL_OBJ BUILD_CK LIB_DIR BOOST_INCLUDE BOOST_LIB build_lapack build_blas BLAS_LAPACK_LIBS BLAS_LAPACK_DIR build_with_f2c build_f2c_lib LOCAL_LIB_DIRS LOCAL_LIBS CT_SHARED_LIB PYTHON_CMD BUILD_PYTHON NUMARRAY_INC_DIR NUMARRAY_HOME CANTERA_PYTHON_HOME CVSTAG MATLAB_CMD BUILD_MATLAB BUILD_CLIB export_name INSTALL_PROGRAM INSTALL_SCRIPT INSTALL_DATA CC CFLAGS ac_ct_CC CXXCPP EGREP SOEXT SHARED PIC CXX_INCLUDES LCXX_FLAGS LCXX_END_LIBS HAVE_STRIPSYMBOLS F77 FFLAGS ac_ct_F77 F90 BUILD_F90 F90FLAGS F90BUILDFLAGS precompile_headers FLIBS OS_IS_DARWIN OS_IS_WIN OS_IS_CYGWIN SHARED_CTLIB mex_ext F77_EXT CXX_EXT OBJ_EXT EXE_EXT math_libs SO LDSHARED EXTRA_LINK TSCOMPARE_abs INSTALL_abs INSTALL_VERBOSE LIBOBJS LTLIBOBJS' ac_subst_files='' # Initialize some variables set by options. @@ -1705,6 +1705,7 @@ GRAPHVIZDIR=" " fi + #----------- ARCHIVE -------------------- if test "x${OS_IS_DARWIN}" = "x1"; then @@ -2644,6 +2645,27 @@ fi +if test "$WITH_HTML_LOG_FILES" = "y"; then + cat >>confdefs.h <<\_ACEOF +#define WITH_HTML_LOGS 1 +_ACEOF + +fi + +BOOST_INCLUDE= +BOOST_LIB= +if test "$BUILD_THREAD_SAFE" = "y" ; then + cat >>confdefs.h <<\_ACEOF +#define THREAD_SAFE_CANTERA 1 +_ACEOF + + BOOST_INCLUDE=-I$BOOST_INC_DIR + BOOST_LIB=$BOOST_THREAD_LIB +fi + + + + # # Report to the user what has been included/excluded from the compilation # @@ -2655,6 +2677,11 @@ if test "$CANTERA_DEBUG_MODE" = "1" ; then else echo " Cantera Debug Mode = OFF" fi +if test "$BUILD_THREAD_SAFE" = "y" ; then + echo " Thread Safe = YES" +else + echo " Thread Safe = NO" +fi if test "$NEED_TRANSPORT" = "1" ; then echo " Calculation of Transport Props = ON" else @@ -2805,10 +2832,6 @@ fi LOCAL_LIBS=$LOCAL_LIBS' '-lequil -if test -n "$NEED_CATHERMO" -then LOCAL_LIBS=$LOCAL_LIBS' '-lthermo -fi - if test -n "$COMPILE_KINETICS" then LOCAL_LIBS=$LOCAL_LIBS' '-lkinetics fi @@ -2817,6 +2840,10 @@ if test -n "$NEED_TRANSPORT" then LOCAL_LIBS=$LOCAL_LIBS' '-ltransport fi +if test -n "$NEED_CATHERMO" +then LOCAL_LIBS=$LOCAL_LIBS' '-lthermo +fi + LOCAL_LIBS=$LOCAL_LIBS' '-lctnumerics if test -n "$NEED_CVODE"; then @@ -2859,6 +2886,12 @@ if test ${use_sundials} = 1; then LOCAL_LIB_DIRS=$LOCAL_LIB_DIRS' -L'$SUNDIALS_LIB_DIR fi +if test ${BUILD_THREAD_SAFE} = "y"; then +LOCAL_LIBS=$LOCAL_LIBS' '-l$BOOST_LIB +LOCAL_LIB_DIRS=$LOCAL_LIB_DIRS' -L'$BOOST_LIB_DIR +fi + + @@ -8178,6 +8211,7 @@ if test -z "$SHARED"; then SHARED='-shared'; fi if test -z "$PIC"; then PIC='-fPIC'; fi +CXX_INCLUDES=$BOOST_INCLUDE # # LCXX_FLAGS: Flags that get attached to the CXX loader @@ -8339,7 +8373,7 @@ fi # Provide some information about the compiler. -echo "$as_me:8342:" \ +echo "$as_me:8376:" \ "checking for Fortran 77 compiler version" >&5 ac_compiler=`set X $ac_compile; echo $2` { (eval echo "$as_me:$LINENO: \"$ac_compiler --version &5\"") >&5 @@ -8806,7 +8840,7 @@ fi - ac_config_files="$ac_config_files Makefile Cantera/Makefile Cantera/src/Makefile Cantera/src/base/Makefile Cantera/src/zeroD/Makefile Cantera/src/oneD/Makefile Cantera/src/converters/Makefile Cantera/src/transport/Makefile Cantera/src/thermo/Makefile Cantera/src/kinetics/Makefile Cantera/src/numerics/Makefile Cantera/src/equil/Makefile Cantera/clib/src/Makefile Cantera/fortran/src/Makefile Cantera/fortran/f77demos/f77demos.mak Cantera/fortran/f77demos/isentropic.dsp Cantera/matlab/Makefile Cantera/matlab/setup_matlab.py Cantera/matlab/setup_winmatlab.py Cantera/python/Makefile Cantera/python/setup.py Cantera/cxx/Makefile Cantera/cxx/src/Makefile Cantera/cxx/demos/Makefile Cantera/user/Makefile Cantera/python/src/Makefile ext/lapack/Makefile ext/blas/Makefile ext/cvode/Makefile ext/math/Makefile ext/recipes/Makefile ext/tpx/Makefile ext/Makefile ext/f2c_libs/Makefile ext/f2c_blas/Makefile ext/f2c_lapack/Makefile ext/f2c_math/Makefile examples/Makefile examples/cxx/Makefile tools/Makefile tools/doc/Cantera.cfg tools/doc/Makefile tools/src/Makefile tools/src/sample.mak tools/src/finish_install.py tools/src/package4mac tools/templates/f77/demo.mak tools/templates/f90/demo.mak tools/templates/cxx/demo.mak tools/testtools/Makefile data/inputs/Makefile test_problems/Makefile test_problems/cxx_ex/Makefile test_problems/silane_equil/Makefile test_problems/surfkin/Makefile test_problems/diamondSurf/Makefile test_problems/ChemEquil_gri_matrix/Makefile test_problems/ChemEquil_gri_pairs/Makefile test_problems/ChemEquil_ionizedGas/Makefile test_problems/ChemEquil_red1/Makefile test_problems/fracCoeff/Makefile test_problems/negATest/Makefile test_problems/ck2cti_test/Makefile test_problems/ck2cti_test/runtest test_problems/min_python/Makefile test_problems/min_python/minDiamond/Makefile test_problems/min_python/negATest/Makefile test_problems/pureFluidTest/Makefile test_problems/python/Makefile test_problems/cathermo/Makefile test_problems/cathermo/issp/Makefile test_problems/cathermo/ims/Makefile test_problems/cathermo/stoichSubSSTP/Makefile test_problems/cathermo/testIAPWS/Makefile test_problems/cathermo/testIAPWSPres/Makefile test_problems/cathermo/testIAPWSTripP/Makefile test_problems/cathermo/testWaterPDSS/Makefile test_problems/cathermo/testWaterTP/Makefile test_problems/cathermo/HMW_test_1/Makefile test_problems/cathermo/HMW_test_3/Makefile test_problems/cathermo/HMW_graph_GvT/Makefile test_problems/cathermo/HMW_graph_GvI/Makefile test_problems/cathermo/HMW_graph_HvT/Makefile test_problems/cathermo/HMW_graph_CpvT/Makefile test_problems/cathermo/HMW_graph_VvT/Makefile test_problems/cathermo/DH_graph_1/Makefile test_problems/cathermo/DH_graph_acommon/Makefile test_problems/cathermo/DH_graph_NM/Makefile test_problems/cathermo/DH_graph_Pitzer/Makefile test_problems/cathermo/DH_graph_bdotak/Makefile bin/install_tsc" + ac_config_files="$ac_config_files Makefile Cantera/Makefile Cantera/src/Makefile Cantera/src/app/Makefile Cantera/src/base/Makefile Cantera/src/zeroD/Makefile Cantera/src/oneD/Makefile Cantera/src/converters/Makefile Cantera/src/transport/Makefile Cantera/src/thermo/Makefile Cantera/src/kinetics/Makefile Cantera/src/numerics/Makefile Cantera/src/equil/Makefile Cantera/clib/src/Makefile Cantera/fortran/src/Makefile Cantera/fortran/f77demos/f77demos.mak Cantera/fortran/f77demos/isentropic.dsp Cantera/matlab/Makefile Cantera/matlab/setup_matlab.py Cantera/matlab/setup_winmatlab.py Cantera/python/Makefile Cantera/python/setup.py Cantera/cxx/Makefile Cantera/cxx/src/Makefile Cantera/cxx/demos/Makefile Cantera/user/Makefile Cantera/python/src/Makefile ext/lapack/Makefile ext/blas/Makefile ext/cvode/Makefile ext/math/Makefile ext/recipes/Makefile ext/tpx/Makefile ext/Makefile ext/f2c_libs/Makefile ext/f2c_blas/Makefile ext/f2c_lapack/Makefile ext/f2c_math/Makefile examples/Makefile examples/cxx/Makefile tools/Makefile tools/doc/Cantera.cfg tools/doc/Makefile tools/src/Makefile tools/src/sample.mak tools/src/finish_install.py tools/src/package4mac tools/templates/f77/demo.mak tools/templates/f90/demo.mak tools/templates/cxx/demo.mak tools/testtools/Makefile data/inputs/Makefile test_problems/Makefile test_problems/cxx_ex/Makefile test_problems/silane_equil/Makefile test_problems/surfkin/Makefile test_problems/diamondSurf/Makefile test_problems/ChemEquil_gri_matrix/Makefile test_problems/ChemEquil_gri_pairs/Makefile test_problems/ChemEquil_ionizedGas/Makefile test_problems/ChemEquil_red1/Makefile test_problems/fracCoeff/Makefile test_problems/negATest/Makefile test_problems/ck2cti_test/Makefile test_problems/ck2cti_test/runtest test_problems/min_python/Makefile test_problems/min_python/minDiamond/Makefile test_problems/min_python/negATest/Makefile test_problems/pureFluidTest/Makefile test_problems/python/Makefile test_problems/cathermo/Makefile test_problems/cathermo/issp/Makefile test_problems/cathermo/ims/Makefile test_problems/cathermo/stoichSubSSTP/Makefile test_problems/cathermo/testIAPWS/Makefile test_problems/cathermo/testIAPWSPres/Makefile test_problems/cathermo/testIAPWSTripP/Makefile test_problems/cathermo/testWaterPDSS/Makefile test_problems/cathermo/testWaterTP/Makefile test_problems/cathermo/HMW_test_1/Makefile test_problems/cathermo/HMW_test_3/Makefile test_problems/cathermo/HMW_graph_GvT/Makefile test_problems/cathermo/HMW_graph_GvI/Makefile test_problems/cathermo/HMW_graph_HvT/Makefile test_problems/cathermo/HMW_graph_CpvT/Makefile test_problems/cathermo/HMW_graph_VvT/Makefile test_problems/cathermo/DH_graph_1/Makefile test_problems/cathermo/DH_graph_acommon/Makefile test_problems/cathermo/DH_graph_NM/Makefile test_problems/cathermo/DH_graph_Pitzer/Makefile test_problems/cathermo/DH_graph_bdotak/Makefile bin/install_tsc" test "x$prefix" = xNONE && prefix=$ac_default_prefix @@ -9282,6 +9316,7 @@ do "Makefile" ) CONFIG_FILES="$CONFIG_FILES Makefile" ;; "Cantera/Makefile" ) CONFIG_FILES="$CONFIG_FILES Cantera/Makefile" ;; "Cantera/src/Makefile" ) CONFIG_FILES="$CONFIG_FILES Cantera/src/Makefile" ;; + "Cantera/src/app/Makefile" ) CONFIG_FILES="$CONFIG_FILES Cantera/src/app/Makefile" ;; "Cantera/src/base/Makefile" ) CONFIG_FILES="$CONFIG_FILES Cantera/src/base/Makefile" ;; "Cantera/src/zeroD/Makefile" ) CONFIG_FILES="$CONFIG_FILES Cantera/src/zeroD/Makefile" ;; "Cantera/src/oneD/Makefile" ) CONFIG_FILES="$CONFIG_FILES Cantera/src/oneD/Makefile" ;; @@ -9530,6 +9565,8 @@ s,@KERNEL@,$KERNEL,;t t s,@KERNEL_OBJ@,$KERNEL_OBJ,;t t s,@BUILD_CK@,$BUILD_CK,;t t s,@LIB_DIR@,$LIB_DIR,;t t +s,@BOOST_INCLUDE@,$BOOST_INCLUDE,;t t +s,@BOOST_LIB@,$BOOST_LIB,;t t s,@build_lapack@,$build_lapack,;t t s,@build_blas@,$build_blas,;t t s,@BLAS_LAPACK_LIBS@,$BLAS_LAPACK_LIBS,;t t diff --git a/configure.in b/configure.in index 381c980e0..47a377854 100755 --- a/configure.in +++ b/configure.in @@ -276,6 +276,7 @@ GRAPHVIZDIR=" " fi AC_SUBST(GRAPHVIZDIR) + #----------- ARCHIVE -------------------- if test "x${OS_IS_DARWIN}" = "x1"; then @@ -539,6 +540,21 @@ AC_SUBST(KERNEL_OBJ) AC_SUBST(BUILD_CK) AC_SUBST(LIB_DIR) +if test "$WITH_HTML_LOG_FILES" = "y"; then + AC_DEFINE(WITH_HTML_LOGS) +fi + +BOOST_INCLUDE= +BOOST_LIB= +if test "$BUILD_THREAD_SAFE" = "y" ; then + AC_DEFINE(THREAD_SAFE_CANTERA) + BOOST_INCLUDE=-I$BOOST_INC_DIR + BOOST_LIB=$BOOST_THREAD_LIB +fi +AC_SUBST(BOOST_INCLUDE) +AC_SUBST(BOOST_LIB) + + # # Report to the user what has been included/excluded from the compilation # @@ -550,6 +566,11 @@ if test "$CANTERA_DEBUG_MODE" = "1" ; then else echo " Cantera Debug Mode = OFF" fi +if test "$BUILD_THREAD_SAFE" = "y" ; then + echo " Thread Safe = YES" +else + echo " Thread Safe = NO" +fi if test "$NEED_TRANSPORT" = "1" ; then echo " Calculation of Transport Props = ON" else @@ -754,6 +775,12 @@ if test ${use_sundials} = 1; then LOCAL_LIB_DIRS=$LOCAL_LIB_DIRS' -L'$SUNDIALS_LIB_DIR fi +if test ${BUILD_THREAD_SAFE} = "y"; then +LOCAL_LIBS=$LOCAL_LIBS' '-l$BOOST_LIB +LOCAL_LIB_DIRS=$LOCAL_LIB_DIRS' -L'$BOOST_LIB_DIR +fi + + AC_SUBST(LOCAL_LIB_DIRS) AC_SUBST(LOCAL_LIBS) @@ -1028,6 +1055,7 @@ AC_SUBST(SHARED) if test -z "$PIC"; then PIC='-fPIC'; fi AC_SUBST(PIC) +CXX_INCLUDES=$BOOST_INCLUDE AC_SUBST(CXX_INCLUDES) # # LCXX_FLAGS: Flags that get attached to the CXX loader @@ -1358,6 +1386,7 @@ dnl Checks for library functions. AC_OUTPUT(Makefile \ Cantera/Makefile \ Cantera/src/Makefile \ + Cantera/src/app/Makefile \ Cantera/src/base/Makefile \ Cantera/src/zeroD/Makefile \ Cantera/src/oneD/Makefile \ diff --git a/ext/f2c_libs/arith.h b/ext/f2c_libs/arith.h index 9841db38a..508eb414f 100644 --- a/ext/f2c_libs/arith.h +++ b/ext/f2c_libs/arith.h @@ -1,3 +1,4 @@ -#define IEEE_8087 -#define Arith_Kind_ASL 1 +#define IEEE_MC68k +#define Arith_Kind_ASL 2 +#define Double_Align #define NANCHECK diff --git a/preconfig b/preconfig index a4f7a0c1f..73e4fa7c6 100755 --- a/preconfig +++ b/preconfig @@ -244,6 +244,12 @@ ENABLE_RXNPATH='y' # two phase pure fluids ENABLE_TPX='y' +# write HTMl log files. Some multiphase equilibrium procedures can +# write copious diagnostic log messages. Set this to anything but 'y' +# to remove this capability (results in slightly faster equilibrium +# calculations) +WITH_HTML_LOG_FILES='y' + #----------------------------------------------------------------- # CVODE / CVODES #----------------------------------------------------------------- @@ -343,6 +349,18 @@ PIC=${PIC:=-fPIC} # the compiler option to create a shared library from object files SHARED=${SHARED:="-dynamic"} +#------------------------------------------------------------------- +# Thread Safety +#------------------------------------------------------------------- + +# Cantera can be built so that it is thread safe. Doing so requires +# using procedures from the boost library, so if you want thread +# safety then you need to get and install boost if you don't have it. + +BUILD_THREAD_SAFE=${BUILD_THREAD_SAFE:="y"} +BOOST_INC_DIR=${BOOST_INC_DIR:="/usr/local/include/boost-1_34"} +BOOST_LIB_DIR=${BOOST_LIB_DIR:="/usr/local/lib"} +BOOST_THREAD_LIB=${BOOST_THREAD_LIB:="boost_thread-mt-1_34"} #------------------------------------------------------------------- # External procedures @@ -497,6 +515,13 @@ export WITH_IDEAL_SOLUTIONS export WITH_ELECTROLYTES export WITH_PRIME +export BUILD_THREAD_SAFE +export BOOST_INC_DIR +export BOOST_LIB_DIR +export BOOST_THREAD_LIB + +export WITH_HTML_LOG_FILES + #cd config chmod +x ./configure chmod +x config/config.guess diff --git a/test_problems/diamondSurf/Makefile.in b/test_problems/diamondSurf/Makefile.in index ed6fdeadf..27476d25b 100644 --- a/test_problems/diamondSurf/Makefile.in +++ b/test_problems/diamondSurf/Makefile.in @@ -64,7 +64,7 @@ LCXX_FLAGS = -L$(CANTERA_LIBDIR) @LOCAL_LIB_DIRS@ @CXXFLAGS@ # How to compile C++ source files to object files .@CXX_EXT@.@OBJ_EXT@: Interface.h - $(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS) + $(CXX) -c $< -I$(CANTERA_INCDIR) @CXX_INCLUDES@ $(CXX_FLAGS) # How to compile the dependency file .cpp.d: diff --git a/test_problems/min_python/minDiamond/Makefile.in b/test_problems/min_python/minDiamond/Makefile.in index f74d26286..cf0c13445 100644 --- a/test_problems/min_python/minDiamond/Makefile.in +++ b/test_problems/min_python/minDiamond/Makefile.in @@ -64,7 +64,7 @@ LCXX_FLAGS = -L$(CANTERA_LIBDIR) @LOCAL_LIB_DIRS@ @CXXFLAGS@ # How to compile C++ source files to object files .@CXX_EXT@.@OBJ_EXT@: Interface.h - $(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS) + $(CXX) -c $< -I$(CANTERA_INCDIR) @CXX_INCLUDES@ $(CXX_FLAGS) # How to compile the dependency file .cpp.d: diff --git a/test_problems/min_python/negATest/Makefile.in b/test_problems/min_python/negATest/Makefile.in index 42b43746b..5e2625b2c 100644 --- a/test_problems/min_python/negATest/Makefile.in +++ b/test_problems/min_python/negATest/Makefile.in @@ -64,7 +64,7 @@ LCXX_FLAGS = -L$(CANTERA_LIBDIR) @LOCAL_LIB_DIRS@ @CXXFLAGS@ # How to compile C++ source files to object files .@CXX_EXT@.@OBJ_EXT@: Interface.h - $(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS) + $(CXX) -c $< -I$(CANTERA_INCDIR) @CXX_INCLUDES@ $(CXX_FLAGS) # How to compile the dependency file .cpp.d: diff --git a/test_problems/negATest/Makefile.in b/test_problems/negATest/Makefile.in index 9d98fdc50..fa67f8003 100644 --- a/test_problems/negATest/Makefile.in +++ b/test_problems/negATest/Makefile.in @@ -64,7 +64,7 @@ LCXX_FLAGS = -L$(CANTERA_LIBDIR) @LOCAL_LIB_DIRS@ @CXXFLAGS@ # How to compile C++ source files to object files .@CXX_EXT@.@OBJ_EXT@: Interface.h - $(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS) + $(CXX) -c $< -I$(CANTERA_INCDIR) @CXX_INCLUDES@ $(CXX_FLAGS) # How to compile the dependency file .cpp.d: diff --git a/test_problems/pureFluidTest/Makefile.in b/test_problems/pureFluidTest/Makefile.in index 744c11a6b..7e58be34b 100644 --- a/test_problems/pureFluidTest/Makefile.in +++ b/test_problems/pureFluidTest/Makefile.in @@ -58,7 +58,7 @@ 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 +INCLUDES=-I$(CANTERA_INCDIR) -I$(CANTERA_INCDIR)/kernel @CXX_INCLUDES@ endif # flags passed to the C++ compiler/linker for the linking step diff --git a/test_problems/silane_equil/Makefile.in b/test_problems/silane_equil/Makefile.in index 43ba09ea0..7f69ec1e3 100644 --- a/test_problems/silane_equil/Makefile.in +++ b/test_problems/silane_equil/Makefile.in @@ -64,7 +64,7 @@ LCXX_FLAGS = -L$(CANTERA_LIBDIR) @LOCAL_LIB_DIRS@ @CXXFLAGS@ # How to compile C++ source files to object files .@CXX_EXT@.@OBJ_EXT@: - $(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS) + $(CXX) -c $< -I$(CANTERA_INCDIR) @CXX_INCLUDES@ $(CXX_FLAGS) # How to compile the dependency file .cpp.d: diff --git a/test_problems/surfkin/Makefile.in b/test_problems/surfkin/Makefile.in index af52ab18e..d0871cebe 100644 --- a/test_problems/surfkin/Makefile.in +++ b/test_problems/surfkin/Makefile.in @@ -65,7 +65,7 @@ LCXX_FLAGS = -L$(CANTERA_LIBDIR) @LOCAL_LIB_DIRS@ @CXXFLAGS@ # how to compile C++ source files to object files .@CXX_EXT@.@OBJ_EXT@: Interface.h - $(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS) + $(CXX) -c $< -I$(CANTERA_INCDIR) @CXX_INCLUDES@ $(CXX_FLAGS) DEPENDS=$(OBJS:.o=.d)