From 25ba149aab46eb53f712d34ff25c0900a1e79e0f Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 5 Apr 2012 00:24:31 +0000 Subject: [PATCH] Sorry for monolithic commit. Will break it up in the future. Moved the external libraries to separate library files so that libcantera.a just contains its own namespace externals. Fixed several errors in the equilibrium program that occurred during the port. (int to size_t issues). Moved some equilibrium program headers to the include file system, so that it can link with equilibrium program. Worked on Cantera.mak. Needs more work. Fixed an issue with the Residual virtual base classes within numerics. They didn't inherit due to int to size_t migration. This caused numerous test problems to fail (issue with backwards compatibility - do we want it and how much do we want?). Added csvdiff back so that it's available for shell environment runtests. --- SConstruct | 76 ++- ext/SConscript | 50 +- ext/f2c_math/ddaspk.c | 7 +- include/cantera/base/ct_defs.h | 17 + include/cantera/kinetics/RxnRates.h | 6 +- include/cantera/numerics/ResidJacEval.h | 2 +- include/cantera/thermo/Elements.h | 309 ++++++++++ .../legacy/{Cantera.h => Cantera_legacy.h} | 0 platform/posix/Cantera.mak.in | 29 +- platform/posix/SConscript | 37 +- platform/posix/bin/install_tsc | 59 ++ src/SConscript | 18 +- src/apps/SConscript | 13 +- src/apps/csvdiff.cpp | 22 +- src/apps/mdp_allo.cpp | 21 +- src/equil/BasisOptimize.cpp | 2 +- src/equil/ChemEquil.cpp | 4 +- src/equil/vcs_DoubleStarStar.cpp | 2 +- src/equil/vcs_Gibbs.cpp | 6 +- src/equil/vcs_MultiPhaseEquil.cpp | 9 +- src/equil/vcs_SpeciesProperties.cpp | 4 +- src/equil/vcs_TP.cpp | 6 +- src/equil/vcs_VolPhase.cpp | 7 +- src/equil/vcs_dbocls.c | 5 +- src/equil/vcs_dbolsm.c | 14 +- src/equil/vcs_elem.cpp | 13 +- src/equil/vcs_elem_rearrange.cpp | 8 +- src/equil/vcs_equilibrate.cpp | 10 +- src/equil/vcs_inest.cpp | 6 +- src/equil/vcs_nondim.cpp | 7 +- src/equil/vcs_phaseStability.cpp | 13 +- src/equil/vcs_prep.cpp | 18 +- src/equil/vcs_prob.cpp | 6 +- src/equil/vcs_rearrange.cpp | 4 +- src/equil/vcs_report.cpp | 6 +- src/equil/vcs_root1d.cpp | 2 +- src/equil/vcs_rxnadj.cpp | 10 +- src/equil/vcs_setMolesLinProg.cpp | 10 +- src/equil/vcs_solve.cpp | 22 +- src/equil/vcs_solve_TP.cpp | 37 +- src/equil/vcs_solve_phaseStability.cpp | 8 +- src/equil/vcs_species_thermo.cpp | 6 +- src/equil/vcs_species_thermo.h | 1 + src/equil/vcs_util.cpp | 4 +- src/fortran/SConscript | 3 + src/kinetics/Kinetics.cpp | 4 +- src/kinetics/ReactionStoichMgr.cpp | 2 + src/matlab/SConscript | 9 +- src/numerics/BEulerInt.cpp | 3 +- src/numerics/NonlinearSolver.cpp | 11 +- src/numerics/ResidJacEval.cpp | 2 +- src/python/SConscript | 6 +- src/thermo/Elements.cpp | 549 +++++++++++++++++- src/thermo/IonsFromNeutralVPSSTP.cpp | 2 +- src/transport/TransportFactory.cpp | 1 + 55 files changed, 1323 insertions(+), 185 deletions(-) rename platform/legacy/{Cantera.h => Cantera_legacy.h} (100%) create mode 100755 platform/posix/bin/install_tsc diff --git a/SConstruct b/SConstruct index ef93dd37b..fc829554e 100644 --- a/SConstruct +++ b/SConstruct @@ -7,7 +7,7 @@ Basic usage: 'scons build' - Compile Cantera and the language interfaces using default options. - 'scons clean' - Delete files created while building Cantera. + 'cons clean' - Delete files created while building Cantera. '[sudo] scons install' - Install Cantera. @@ -112,6 +112,12 @@ env = Environment(tools=toolchain+['textfile', 'subst', 'recursiveInstall', 'wix toolchain=toolchain, **extraEnvArgs) +# +# To print the current environment +# +# print env.Dump() + + env['OS'] = platform.system() env['OS_BITS'] = int(platform.architecture()[0][:2]) @@ -138,8 +144,10 @@ class defaults: pass if os.name == 'posix': defaults.prefix = '/usr/local' - defaults.boostIncDir = '/usr/include' - defaults.boostLibDir = '/usr/lib' + # defaults.boostIncDir = '/usr/include' + # defaults.boostLibDir = '/usr/lib' + defaults.boostIncDir = '' + defaults.boostLibDir = '' env['INSTALL_MANPAGES'] = True elif os.name == 'nt': defaults.prefix = pjoin(os.environ['ProgramFiles'], 'Cantera') @@ -568,6 +576,15 @@ opts.AddVariables( """Create symbolic links for headers that were installed to the 'kernel' subdirectory in previous versions of Cantera.""", False), + BoolVariable( + 'renamed_shared_libraries', + """If this option is turned on, which is the default, the shared libraries that are created + will be renamed to have a "_shared" extension added to their base name. + If not, the base names will be the same as the static libraries. + In some cases this simplifies subsequent linking environments with + static libaries and avoids a bug with using valgrind with + the -static linking flag.""", + True), PathVariable( 'graphvizdir', """The directory location of the graphviz program, "dot". dot is @@ -931,7 +948,7 @@ cdefine('WITH_METAL', 'with_metal') cdefine('WITH_STOICH_SUBSTANCE', 'with_stoich_substance') cdefine('WITH_SEMICONDUCTOR', 'with_semiconductor') cdefine('WITH_PRIME', 'with_prime') -cdefine('H298MODIFY_CAPABILITY', 'with_n298modify_capability') +cdefine('H298MODIFY_CAPABILITY', 'with_h298modify_capability') cdefine('WITH_PURE_FLUIDS', 'with_pure_fluids') cdefine('WITH_HTML_LOGS', 'with_html_log_files') cdefine('WITH_VCSNONIDEAL', 'with_vcsnonideal') @@ -1018,7 +1035,7 @@ if 'install' in COMMAND_LINE_TARGETS: # Make symlinks to replicate old header directory structure if env['legacy_headers']: install(env.Command, pjoin('$inst_incdir', 'kernel'), [], Mkdir("$TARGET")) - install('$inst_incdir', 'platform/legacy/Cantera.h') + install('$inst_incdir', 'platform/legacy/Cantera_legacy.h') if env['OS'] == 'Windows': cmd = Copy("$TARGET", "$SOURCE") @@ -1050,19 +1067,60 @@ if 'install' in COMMAND_LINE_TARGETS: ### List of libraries needed to link to Cantera ### linkLibs = ['cantera'] +linkSharedLibs = ['cantera_shared'] + if env['use_sundials'] == 'y': env['sundials_libs'] = ['sundials_cvodes', 'sundials_ida', 'sundials_nvecserial'] linkLibs.extend(('sundials_cvodes', 'sundials_ida', 'sundials_nvecserial')) + linkSharedLibs.extend(('sundials_cvodes', 'sundials_ida', 'sundials_nvecserial')) else: env['sundials_libs'] = [] + linkLibs.extend(['cvode']) + linkSharedLibs.extend(['cvode_shared']) + #print 'linkLibs = ', linkLibs -linkLibs.extend(env['blas_lapack_libs']) +linkLibs.append('ctmath') +linkSharedLibs.append('ctmath_shared') + +# +# Add lapack and blas to the link line +# +if env['blas_lapack_libs'] == []: + linkLibs.extend(('ctlapack', 'ctblas')) + linkSharedLibs.extend(('ctlapack_shared', 'ctblas_shared')) +else: + linkLibs.extend(env['blas_lapack_libs']) +# +# Add execstream to the link line +# +linkLibs.append('execstream') +linkSharedLibs.append('execstream_shared') +# +# Add the f2c library if it is necessary to link fortran libraries +# if not env['build_with_f2c']: + # + # Early gcc compilers will fail on this line as they use g77 and not gfortran + # linkLibs.append('gfortran') + linkSharedLibs.append('gfortran') +else: + # Add the f2c library when f2c is requested + # + linkLibs.append('ctf2c') + linkSharedLibs.append('ctf2c_shared') +# +# Store the list of needed static link libraries in the environment +# env['cantera_libs'] = linkLibs +env['cantera_shared_libs'] = linkSharedLibs +if env['renamed_shared_libraries'] == False : + env['cantera_shared_libs'] = linkLibs + + # Add targets from the SConscript files in the various subdirectories Export('env', 'build', 'libraryTargets', 'install', 'buildSample') @@ -1082,15 +1140,15 @@ SConscript('build/src/SConscript') if env['python_package'] in ('full','minimal'): SConscript('src/python/SConscript') -if env['matlab_toolbox'] == 'y': - SConscript('build/src/matlab/SConscript') - SConscript('build/src/apps/SConscript') if env['OS'] != 'Windows': VariantDir('build/platform', 'platform/posix', duplicate=0) SConscript('build/platform/SConscript') +if env['matlab_toolbox'] == 'y': + SConscript('build/src/matlab/SConscript') + if env['doxygen_docs'] or env['sphinx_docs']: SConscript('doc/SConscript') diff --git a/ext/SConscript b/ext/SConscript index 22f84e1a0..75ef9e52b 100644 --- a/ext/SConscript +++ b/ext/SConscript @@ -96,13 +96,57 @@ if env['use_sundials'] == 'n': for subdir, extensions, prepFunction in libs: localenv = prepFunction(env) - objects = localenv.SharedObject(mglob(localenv, subdir, *extensions)) - libraryTargets.extend(objects) + # objects = localenv.SharedObject(mglob(localenv, subdir, *extensions)) + # libraryTargets.extend(objects) + libName = 'ct' + subdir + if libName == 'ctf2c_blas': + libName = 'ctblas' + if libName == 'ctf2c_lapack': + libName = 'ctlapack' + if libName == 'ctf2c_math' : + libName = 'ctmath' + if libName == 'ctcvode/source' : + libName = 'cvode' + if libName == 'ctlibexecstream' : + libName = 'execstream' + if libName == 'ctf2c_libs': + libName = 'ctf2c' + + + if localenv['renamed_shared_libraries'] : + sharedLibName = libName + '_shared' + else: + sharedLibName = libName + + # + # Build the static library + # + slTargets = [] + + objects = localenv.SharedObject(mglob(localenv, subdir, *extensions)) + slTargets.extend(objects) + # ttarget = pjoin('..', 'lib', libName) + ss = pjoin('..', 'lib', libName) + lib = build(localenv.StaticLibrary(ss, slTargets, SPAWN=getSpawn(localenv))) + #lib = build(localenv.StaticLibrary('../lib/cantera', libraryTargets, + # SPAWN=getSpawn(localenv))) + + install('$inst_libdir', lib) + # + # Build the shared library + # + liby = build(localenv.SharedLibrary(target = pjoin('..', 'lib', sharedLibName), source = mglob(localenv, subdir, *extensions))) + #liby = build(localenv.SharedLibrary(target = pjoin('..', 'lib', sharedLibName), source = mglob(localenv, subdir, *extensions), + # SPAWN=getSpawn(localenv))) + install('$inst_libdir', liby) + -# Google Test localenv = env.Clone() localenv.Append(CPPPATH=[Dir('#ext/gtest'), Dir('#ext/gtest/include')], CPPDEFINES={'GTEST_HAS_PTHREAD': 0}) build(localenv.Library(pjoin('../lib', 'gtest'), source=['gtest/src/gtest-all.cc'])) + + + diff --git a/ext/f2c_math/ddaspk.c b/ext/f2c_math/ddaspk.c index 052eb72cd..299eb91ae 100644 --- a/ext/f2c_math/ddaspk.c +++ b/ext/f2c_math/ddaspk.c @@ -7,6 +7,7 @@ /* Table of constant values */ +static long int lc__4 = 4; static integer c__49 = 49; static integer c__201 = 201; static integer c__0 = 0; @@ -164,7 +165,7 @@ static integer c__926 = 926; integer lenwp, lenrw, mxord, nwarn; doublereal rtoli; integer lsavr; - extern doublereal d1mach_(integer *); + extern doublereal d1mach_(long int *); doublereal tdist, tnext, fmaxl; extern /* Subroutine */ int ddstp_(doublereal *, doublereal *, doublereal *, integer *, U_fp, U_fp, U_fp, doublereal *, doublereal *, @@ -2016,8 +2017,8 @@ L200: } /* Compute unit roundoff and HMIN. */ - - uround = d1mach_(&c__4); + + uround = d1mach_(&lc__4); rwork[9] = uround; /* Computing MAX */ d__1 = abs(*t), d__2 = abs(*tout); diff --git a/include/cantera/base/ct_defs.h b/include/cantera/base/ct_defs.h index 8e7d72418..83717fc45 100644 --- a/include/cantera/base/ct_defs.h +++ b/include/cantera/base/ct_defs.h @@ -19,6 +19,7 @@ #include // STL includes +#include #include #include #include @@ -166,6 +167,22 @@ const doublereal Undef = -999.1234; //! Small number to compare differences of mole fractions against. const doublereal Tiny = 1.e-20; +//! inline function to return the max value of two doubles. + /*! + * @param x double value + * @param y second double value + */ +inline doublereal fmaxx(doublereal x, doublereal y) +{ return (x > y) ? x : y; } + +//! inline function to return the min value of two doubles. +/*! + * @param x double value + * @param y second double value + */ +inline doublereal fminn(doublereal x, doublereal y) +{ return (x < y) ? x : y; } + //! Map connecting a string name with a double. /*! * This is used mostly to assign concentrations and mole fractions diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h index ee8b48ba6..8d520efed 100644 --- a/include/cantera/kinetics/RxnRates.h +++ b/include/cantera/kinetics/RxnRates.h @@ -548,9 +548,9 @@ public: pressureIter iter = pressures_.upper_bound(c[0]); AssertThrowMsg(iter != pressures_.end(), "Plog::update_C", - "Pressure out of range: " + fp2str(logP)); - AssertThrowMsg(iter != pressures.begin(), "Plog::update_C", - "Pressure out of range: " + fp2str(logP)); + "Pressure out of range: " + fp2str(logP_)); + AssertThrowMsg(iter != pressures_.begin(), "Plog::update_C", + "Pressure out of range: " + fp2str(logP_)); // upper interpolation pressure logP2_ = iter->first; diff --git a/include/cantera/numerics/ResidJacEval.h b/include/cantera/numerics/ResidJacEval.h index 7b7ebe0bc..dfd09b985 100644 --- a/include/cantera/numerics/ResidJacEval.h +++ b/include/cantera/numerics/ResidJacEval.h @@ -122,7 +122,7 @@ public: const doublereal* const ydot, doublereal* const resid, const ResidEval_Type_Enum evalType = Base_ResidEval, - const size_t id_x = npos, + const int id_x = -1, const doublereal delta_x = 0.0); diff --git a/include/cantera/thermo/Elements.h b/include/cantera/thermo/Elements.h index 1e38797f3..7a2503f7c 100644 --- a/include/cantera/thermo/Elements.h +++ b/include/cantera/thermo/Elements.h @@ -95,6 +95,315 @@ namespace Cantera //! @exception CanteraError If a match is not found, throws a CanteraError double LookupWtElements(const std::string& ename); + class XML_Node; + + //! Object containing the elements that make up species in a phase. + /*! + * Class %Elements manages the elements that are part of a + * chemistry specification. This class may support calculations + * employing Multiple phases. In this case, a single Elements object may + * be shared by more than one Constituents class. Reactions between + * the phases may then be described using stoichiometry base on the + * same Elements class object. + * + * The member functions return information about the elements described + * in a particular instantiation of the class. + * + * @ingroup phases + */ + class Elements { + + public: + + //! Default constructor for the elements class + Elements(); + + //! Default destructor for the elements class + ~Elements(); + + + //! copy constructor + /*! + * This copy constructor just calls the assignment operator for this + * class. It sets the number of subscribers to zer0. + * + * @param right Reference to the object to be copied. + */ + Elements(const Elements& right); + + //! Assigntment operator + /*! + * This is the assignment operator for the Elements class. + * Right now we pretty much do a straight uncomplicated + * assignment. However, subscribers are not mucked with, as they + * have to do with the address of the object to be subscribed to + * + * @param right Reference to the object to be copied. + */ + Elements& operator=(const Elements& right); + + + //! Static function to look up an atomic weight + /*! + * This static function looks up the argument string in the + * database above and returns the associated molecular weight. + * The data are from the periodic table. + * + * Note: The idea behind this function is to provide a unified + * source for the element atomic weights. This helps to + * ensure that mass is conserved. + * + * @param ename String, Only the first 3 characters are significant + * + * @return + * Return value contains the atomic weight of the element + * If a match for the string is not found, a value of -1.0 is + * returned. + * + * @exception CanteraError + * If a match is not found, a CanteraError is thrown as well + */ + static double LookupWtElements(const std::string &ename); + /// Atomic weight of element m. + /*! + * @param m element index + */ + doublereal atomicWeight(int m) const { return m_atomicWeights[m]; } + + /// Atomic number of element m. + /*! + * @param m element index + */ + int atomicNumber(int m) const { return m_atomicNumbers[m]; } + + //! Entropy at 298.15 K and 1 bar of stable state + //! of the element + /*! + * units J kmol-1 K-1 + * + * @param m Element index + */ + doublereal entropyElement298(int m) const; + + //! Return the element constraint type + /*! + * Possible types include: + * + * CT_ELEM_TYPE_ABSPOS 0 + * CT_ELEM_TYPE_ELECTRONCHARGE 1 + * CT_ELEM_TYPE_CHARGENEUTRALITY 2 + * CT_ELEM_TYPE_LATTICERATIO 3 + * CT_ELEM_TYPE_KINETICFROZEN 4 + * CT_ELEM_TYPE_SURFACECONSTRAINT 5 + * CT_ELEM_TYPE_OTHERCONSTRAINT 6 + * + * The default is CT_ELEM_TYPE_ABSPOS + * + * @param m Element index + * + * @return Returns the element type + */ + int elementType(int m) const; + + //! Change the element type of the mth constraint + /*! + * Reassigns an element type + * + * @param m Element index + * @param elem_type New elem type to be assigned + * + * @return Returns the old element type + */ + int changeElementType(int m, int elem_type); + + /// vector of element atomic weights + const vector_fp& atomicWeights() const { return m_atomicWeights; } + /** + * Inline function that returns the number of elements in the object. + * + * @return + * \c int: The number of elements in the object. + */ + int nElements() const { return m_mm; } + + //! Function that returns the index of an element. + /*! + * Index of element named \c name. The index is an integer + * assigned to each element in the order it was added, + * beginning with 0 for the first element. If \c name is not + * the name of an element in the set, then the value -1 is + * returned. + * + * @param name String containing the index. + */ + int elementIndex(std::string name) const; + + //! Name of the element with index \c m. + /*! + * @param m Element index. If m < 0 or m >= nElements() an exception is thrown. + */ + std::string elementName(int m) const; + + //! Returns a string vector containing the element names + /*! + * Returns a read-only reference to the vector of element names. + * @return const vector& : The vector contains + * the element names in their indexed order. + */ + const std::vector& elementNames() const { + return m_elementNames; + } + + //! Add an element to the current set of elements in the current object. + /*! + * The default weight is a special value, which will cause the + * routine to look up the actual weight via a string lookup. + * + * There are two interfaces to this routine. The XML interface + * looks up the required parameters for the regular interface + * and then calls the base routine. + * + * @param symbol string symbol for the element. + * @param weight Atomic weight of the element. If no argument + * is provided, a lookup is attempted. + */ + void addElement(const std::string& symbol, + doublereal weight = -12345.0); + //! Add an element to the current set of elements in the current object. + /*! + * @param e Reference to the XML_Node containing the element information + * The node name is the element symbol and the atomWt attribute + * is used as the atomic weight. + */ + void addElement(const XML_Node& e); + + //! Add an element only if the element hasn't been added before. + /*! + * This is accomplished via a string match on symbol. + * + * @param symbol string symbol for the element. + * @param weight Atomic weight of the element. If no argument + * is provided, a lookup is attempted. + * @param atomicNumber defaults to 0 + * @param entropy298 Value of the entropy at 298 and 1 bar of the + * element in its most stable form. + * The default is to specify an ENTROPY298_UNKNOWN value, + * which will cause a throw error if its ever + * needed. + * @param elem_type New elem type to be assigned. + * The default is a regular element, CT_ELEM_TYPE_ABSPOS + */ + void addUniqueElement(const std::string& symbol, + doublereal weight = -12345.0, int atomicNumber = 0, + doublereal entropy298 = ENTROPY298_UNKNOWN, int elem_type = CT_ELEM_TYPE_ABSPOS); + + //! Add an element to the current set of elements in the current object. + /*! + * @param e Reference to the XML_Node containing the element information + * The node name is the element symbol and the atomWt attribute + * is used as the atomic weight. + */ + void addUniqueElement(const XML_Node& e); + + //! Add multiple elements from a XML_Node phase description + /*! + * @param phase XML_Node reference to a phase + */ + void addElementsFromXML(const XML_Node& phase); + + //! Prohibit addition of more elements, and prepare to add species. + void freezeElements(); + + //! True if freezeElements has been called. + bool elementsFrozen() const; + + /// Remove all elements + void clear(); + + /// True if both elements and species have been frozen + bool ready() const; + + //! subscribe to this object + /*! + * Increment by one the number of subscriptions to this object. + */ + void subscribe(); + + //! unsubscribe to this object + /*! + * decrement by one the number of subscriptions to this object. + */ + int unsubscribe(); + + //! report the number of subscriptions + int reportSubscriptions() const; + + protected: + + /******************************************************************/ + /* Description of DATA in the Object */ + /******************************************************************/ + + //! Number of elements. + int m_mm; + + /* m_elementsFrozen: */ + /** boolean indicating completion of object + * + * If this is true, then no elements may be added to the + * object. + */ + bool m_elementsFrozen; + + /** + * Vector of element atomic weights: + * + * units = kg / kmol + */ + vector_fp m_atomicWeights; + + /** + * Vector of element atomic numbers: + * + */ + vector_int m_atomicNumbers; + + /** Vector of strings containing the names of the elements + * + * Note, a string search is the primary way to identify elements. + */ + std::vector m_elementNames; + + //! Entropy at 298.15 K and 1 bar of stable state + /*! + * units J kmol-1 + */ + vector_fp m_entropy298; + + //! Vector of element types + vector_int m_elem_type; + /** + * Number of Constituents Objects that use this object + * + * Number of Constituents Objects that require this Elements object + * to complete its definition. + * The destructor checks to see that this is equal to zero. + * when the element object is released. + */ + int numSubscribers; + + /********* GLOBAL STATIC SECTION *************/ + + public: + /** Vector of pointers to Elements Objects + * + */ + static std::vector Global_Elements_List; + + friend class Constituents; + }; + + } // namespace #endif diff --git a/platform/legacy/Cantera.h b/platform/legacy/Cantera_legacy.h similarity index 100% rename from platform/legacy/Cantera.h rename to platform/legacy/Cantera_legacy.h diff --git a/platform/posix/Cantera.mak.in b/platform/posix/Cantera.mak.in index 3a2c6c95a..ec2d08246 100644 --- a/platform/posix/Cantera.mak.in +++ b/platform/posix/Cantera.mak.in @@ -28,7 +28,9 @@ CANTERA_VERSION=@cantera_version@ CANTERA_CORE_INCLUDES=-I@ct_incroot@ # Required Cantera libraries -CANTERA_CORE_LIBS=@mak_threadflags@ -L@ct_libdir@ -lcantera +CANTERA_CORE_LIBS=@mak_threadflags@ -L@ct_libdir@ -lcantera -lctmath -lexecstream + +CANTERA_CORE_LIBS_DEP = @ct_libdir@/libcantera.a CANTERA_CORE_FTN=-L@ct_libdir@ -lcantera_fortran -lcantera @@ -62,6 +64,19 @@ else CANTERA_BLAS_LAPACK_LIBS=@mak_blas_lapack_libs@ endif +############################################################################### +# F2C USAGE +############################################################################### + +CANTERA_mak_have_f2c_lib=@mak_have_f2c_lib@ + +ifeq ($(CANTERA_mak_have_f2c_lib), 1) + CANTERA_F2C_LIBS= -lctf2c +else + CANTERA_F2C_LIBS= +endif + + ############################################################################### # COMBINATIONS OF INCLUDES AND LIBS ############################################################################### @@ -69,11 +84,21 @@ endif CANTERA_INCLUDES=$(CANTERA_CORE_INCLUDES) $(CANTERA_SUNDIALS_INCLUDE) \ $(CANTERA_BOOST_INCLUDES) +CANTERA_TOTAL_INCLUDES = $(CANTERA_INCLUDES) + # Add this into the compilation environment to identify the version number CANTERA_DEFINES = -DCANTERA_VERSION=@cantera_version@ CANTERA_LIBS=$(CANTERA_CORE_LIBS) $(CANTERA_SUNDIALS_LIBS) \ - $(CANTERA_BLAS_LAPACK_LIBS) $(CANTERA_BOOST_LIBS) + $(CANTERA_BLAS_LAPACK_LIBS) $(CANTERA_BOOST_LIBS) \ + $(CANTERA_F2C_LIBS) + +CANTERA_TOTAL_LIBS=$(CANTERA_LIBS) + +CANTERA_TOTAL_LIBS_DEP= $(CANTERA_CORE_LIBS_DEP) \ + $(CANTERA_SUNDIALS_LIBS_DEP) \ + $(CANTERA_BLAS_LAPACK_LIBS_DEP) + CANTERA_FORTRAN_LIBS=$(CANTERA_CORE_FTN) $(CANTERA_SUNDIALS_LIBS) \ $(CANTERA_BLAS_LAPACK_LIBS) $(CANTERA_BOOST_LIBS) \ diff --git a/platform/posix/SConscript b/platform/posix/SConscript index 2ed5fe8a8..b96dce218 100644 --- a/platform/posix/SConscript +++ b/platform/posix/SConscript @@ -21,11 +21,20 @@ inst = install('$inst_bindir', target) install(localenv.AddPostAction, inst, Chmod('$TARGET', 0755)) # Cantera.mak include file for Makefile projects -localenv['mak_sundials_libs'] = ' '.join('-l%s' % s + +if env['use_sundials'] == 'n': + localenv['mak_sundials_libs'] = '-lcvode' + localenv['mak_sundials_libdir'] = '' + localenv['mak_sundials_include'] = '' +else: + # + # Add links to the sundials environment + # + localenv['mak_sundials_libs'] = ' '.join('-l%s' % s for s in localenv['sundials_libs']) -localenv['mak_sundials_libdir'] = ('-L' + localenv['sundials_libdir'] + localenv['mak_sundials_libdir'] = ('-L' + localenv['sundials_libdir'] if localenv['sundials_libdir'] else '') -localenv['mak_sundials_include'] = ('-I' + localenv['sundials_include'] + localenv['mak_sundials_include'] = ('-I' + localenv['sundials_include'] if localenv['sundials_include'] else '') localenv['mak_boost_include'] = ('-I' + localenv['boost_inc_dir'] @@ -36,15 +45,29 @@ localenv['mak_boost_libdir'] = ('-L' + localenv['boost_lib_dir'] localenv['mak_boost_libs'] = ' '.join('-l%s' % s for s in localenv['boost_libs']) - +# +# Handle blas lapack linkage +# localenv['mak_have_blas_lapack_dir'] = '1' if localenv['blas_lapack_dir'] else '0' -localenv['mak_blas_lapack_libs'] = ' '.join('-l%s' % s - for s in localenv['blas_lapack_libs']) +if localenv['blas_lapack_dir']: + localenv['mak_blas_lapack_libs'] = ' '.join('-l%s' % s + for s in localenv['blas_lapack_libs']) +else: + localenv['mak_blas_lapack_libs'] = ('-L' + '$inst_libdir' + ' -lctlapack -lctblas') + localenv['mak_threadflags'] = localenv['thread_flags'] if '-pthread' in localenv['thread_flags']: localenv['mak_fort_threadflags'] = '-lpthread' else: localenv['mak_fort_threadflags'] = '' +# +# Handle f2c Linkage +# +localenv['mak_have_f2c_lib'] = '1' if localenv['build_with_f2c'] else '0' +# print 'make_have-f2c_lib = ', localenv['mak_have_f2c_lib'] + mak = build(localenv.SubstFile('Cantera.mak', 'Cantera.mak.in')) -install('$inst_sampledir', mak) +install('$inst_incdir', mak) + +install('$inst_bindir', '#platform/posix/bin/install_tsc') diff --git a/platform/posix/bin/install_tsc b/platform/posix/bin/install_tsc new file mode 100755 index 000000000..a452a136f --- /dev/null +++ b/platform/posix/bin/install_tsc @@ -0,0 +1,59 @@ +#!/bin/sh +# +# install_ts source DEST_DIR' +# +# Install a file checking the time step first +# +# HKM -> openBSD doesn't have a -v option +# +if test $# -ne 2 ; then + echo 'install_ts ERROR: This program requires 2 and only 2 parameters:' + echo ' install_ts source DEST_DIR' + exit -1 +fi +lh=$1 +INCDIR=$2 +# +VERBOSE=-v +verbose='y' +if test "x$VERBOSE" = "x-v" ; then + verbose= +fi +INSTALL_XV="/usr/bin/install -c -m 755 -c $VERBOSE " +INSTALL_V="/usr/bin/install -c -m 644 -c $VERBOSE " +# +# If the destination directory doesn't exist yet, create it +# +if test ! -d $INCDIR ; then + /usr/bin/install -c -d $INCDIR +fi +th="${INCDIR}"/"${lh}" ; +if test ! -f "${lh}" ; then + echo 'install_ts ERROR: the file, ' ${lh} ', does not exist' + exit -1 +fi +if test ! -f "${th}" ; then + if test -x "${lh}" ; then + $INSTALL_XV "${lh}" "${th}" + else + $INSTALL_V "${lh}" "${th}" + fi + if test x"$verbose" = xy ; then + echo "${lh}" ' -> ' "${th}" + fi +else +# if test "${lh}" -nt "${th}" ; then + /ascldap/users/hkmoffa/Cantera/gc/canteraLiquidTransportDevelop/bin/tscompare "${lh}" "${th}" + res=$? + if test "$res" = "0" ; then + if test -x "${lh}" ; then + $INSTALL_XV "${lh}" "${th}" + else + $INSTALL_V "${lh}" "${th}" + fi + if test x"$verbose" = xy ; then + echo "${lh}" ' -> ' "${th}" + fi + fi +fi +exit 0 diff --git a/src/SConscript b/src/SConscript index 62dbcb073..9753cfdd5 100644 --- a/src/SConscript +++ b/src/SConscript @@ -63,9 +63,21 @@ if (localenv['use_sundials'] == 'y' and if localenv['toolchain'] == 'mingw': localenv.Append(LINKFLAGS=['-static-libgcc', '-static-libstdc++']) -# Build the Cantera shared library -lib = build(localenv.SharedLibrary('../lib/cantera_shared', libraryTargets, - SPAWN=getSpawn(localenv))) +# +# Define the name according to the environmental variable +# +if localenv['renamed_shared_libraries'] == True : + sharedName = '../lib/cantera_shared' +else: + sharedName = '../lib/cantera' + +# Build the Cantera shared library using the correct name +lib = build(localenv.SharedLibrary(sharedName, libraryTargets, SPAWN=getSpawn(localenv))) + env['cantera_shlib'] = lib localenv.Depends(lib, localenv['config_h_target']) +# +# Create an install rule for the shared library +# install('$inst_libdir', lib) + diff --git a/src/apps/SConscript b/src/apps/SConscript index 82fe75184..ed99e8db4 100644 --- a/src/apps/SConscript +++ b/src/apps/SConscript @@ -4,14 +4,21 @@ Import('env', 'build', 'install') localenv = env.Clone() programs = [('cti2ctml', ['cti2ctml.cpp']), - ('ck2cti', ['ck2cti.cpp'])] + ('ck2cti', ['ck2cti.cpp']), + ('csvdiff', ['csvdiff.cpp', 'tok_input_util.cpp', 'mdp_allo.cpp'])] -localenv.Append(CPPPATH=['#src', '#include']) +localenv.Append(CPPPATH=['#src', '#include', '#src/apps']) + +llibs = env['cantera_libs'] +print 'llibs = ', llibs for name, src in programs: prog = build(localenv.Program(target=pjoin('#build/bin', name), source=src, - LIBS=['cantera'])) + LIBS=llibs)) + #prog = build(localenv.Program(target=pjoin('#build/bin', name), + # source=src, + # LIBS=['cantera'])) install('$inst_bindir', prog) # Copy man pages diff --git a/src/apps/csvdiff.cpp b/src/apps/csvdiff.cpp index b6b1b935b..2031faaaf 100644 --- a/src/apps/csvdiff.cpp +++ b/src/apps/csvdiff.cpp @@ -21,7 +21,7 @@ */ #include -#include +#include #include #include #include @@ -38,8 +38,16 @@ using namespace std; #include #endif -#include "cantera/base/mdp_allo.h" +#include "mdp_allo.h" +//#include "cantera/base/mdp_allo.h" #include "tok_input_util.h" +#ifndef MAX +# define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif +#ifndef MIN +# define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) +#endif + int Debug_Flag = true; double grtol = 1.0E-3; @@ -790,8 +798,8 @@ int main(int argc, char* argv[]) * Right now, if the number of data rows differ, we will punt. * Maybe later we can do something more significant */ - int nDataRowsMIN = std::min(nDataRows1, nDataRows2); - int nDataRowsMAX = std::max(nDataRows1, nDataRows2); + int nDataRowsMIN = MIN(nDataRows1, nDataRows2); + int nDataRowsMAX = MAX(nDataRows1, nDataRows2); if (nDataRows1 != nDataRows2) { printf("Number of Data rows in file1, %d, is different than file2, %d\n", nDataRows1, nDataRows2); @@ -805,7 +813,7 @@ int main(int argc, char* argv[]) read_title(fp2, &title2, nTitleLines2); if (nTitleLines1 > 0 || nTitleLines2 > 0) { - int n = std::min(nTitleLines1, nTitleLines2); + int n = MIN(nTitleLines1, nTitleLines2); for (i = 0; i < n; i++) { if (strcmp(title1[i], title2[i]) != 0) { printf("Title Line %d differ:\n\t\"%s\"\n\t\"%s\"\n", i, title1[i], title2[i]); @@ -858,7 +866,7 @@ int main(int argc, char* argv[]) * Do a Comparison of the names to find the maximum number * of matches. */ - nColMAX = std::max(nCol1, nCol2); + nColMAX = MAX(nCol1, nCol2); compColList = mdp_alloc_int_2(nColMAX, 2, -1); nColcomparisons = 0; @@ -939,7 +947,7 @@ int main(int argc, char* argv[]) curVarValues1 = NVValues1[i1]; curVarValues2 = NVValues2[i2]; atol_j = get_atol(curVarValues1, nDataRows1, gatol); - atol_j = std::min(atol_j, get_atol(curVarValues2, nDataRows2, gatol)); + atol_j = MIN(atol_j, get_atol(curVarValues2, nDataRows2, gatol)); for (j = 0; j < nDataRowsMIN; j++) { slope1 = 0.0; diff --git a/src/apps/mdp_allo.cpp b/src/apps/mdp_allo.cpp index f488f04bf..51a7c1f7b 100644 --- a/src/apps/mdp_allo.cpp +++ b/src/apps/mdp_allo.cpp @@ -5,7 +5,7 @@ #include #include -#include "cantera/base/mdp_allo.h" +#include "mdp_allo.h" /* * Allocate global storage for 2 debugging ints that are used in IO of @@ -32,6 +32,13 @@ int MDP_MP_myproc = 0; int MDP_ALLO_errorOption = 3; #define MDP_ALLOC_INTERFACE_ERROR 230346 +#ifndef MAX +# define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif +#ifndef MIN +# define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) +#endif + /****************************************************************************/ /****************************************************************************/ @@ -982,8 +989,8 @@ void mdp_realloc_dbl_2(double** *array_hdl, int ndim1, int ndim2, if (ndim2 <= 0) { ndim2 = 1; } - ndim1Old = std::max(ndim1Old, 0); - ndim2Old = std::max(ndim2Old, 0); + ndim1Old = MAX(ndim1Old, 0); + ndim2Old = MAX(ndim2Old, 0); /* * One way to do it, if old information isn't needed. In this algorithm * the arrays are never malloced at the same time. @@ -1007,8 +1014,8 @@ void mdp_realloc_dbl_2(double** *array_hdl, int ndim1, int ndim2, /* * Now, let's initialize the arrays */ - int ndim1Min = std::min(ndim1, ndim1Old); - int ndim2Min = std::min(ndim2, ndim2Old); + int ndim1Min = MIN(ndim1, ndim1Old); + int ndim2Min = MIN(ndim2, ndim2Old); double** array_new = *array_hdl; /* * When the second dimensions are equal, we can copy blocks @@ -1143,7 +1150,7 @@ void mdp_realloc_VecFixedStrings(char** *array_hdl, int numStrings, } array = (char**) mdp_array_alloc(2, numStrings, lenString, sizeof(char)); if (array != NULL) { - int len = std::min(numStrings, numOldStrings); + int len = MIN(numStrings, numOldStrings); ao = *array_hdl; if (ao) { for (i = 0; i < len; i++) { @@ -1407,7 +1414,7 @@ void mdp_realloc_ptr_1(void** *array_hdl, int numLen, int numOldLen) size_t bytenum = sizeof(void*) * numLen; void** array = (void**) smalloc(bytenum); if (array != NULL) { - int len = std::min(numLen, numOldLen); + int len = MIN(numLen, numOldLen); if (*array_hdl) { void** ao = *array_hdl; for (int i = 0; i < len; i++) { diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp index eeecb99bf..68f288283 100644 --- a/src/equil/BasisOptimize.cpp +++ b/src/equil/BasisOptimize.cpp @@ -8,7 +8,7 @@ #include "cantera/thermo/ThermoPhase.h" #include "cantera/equil/MultiPhase.h" -#include +#include using namespace Cantera; using namespace std; diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index 1fd059e63..0fdd6f1fe 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -1112,7 +1112,7 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx, if (ChemEquil_print_lvl > 0) { writelogf("Solution Unknowns: damp = %g\n", damp); writelog(" X_new X_old Step\n"); - for (m = 0; m < nvar; m++) { + for (size_t m = 0; m < m_mm; m++) { writelogf(" % -10.5g % -10.5g % -10.5g\n", x[m], oldx[m], step[m]); } } @@ -1165,7 +1165,7 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, if (ChemEquil_print_lvl > 0 && !m_doResPerturb) { PrintCtrl pc(std::cout, -14, PrintCtrl::CT_OFF_GLOBALOBEY); writelog("Residual: ElFracGoal ElFracCurrent Resid\n"); - for (n = 0; n < m_mm; n++) { + for (int n = 0; n < m_mm; n++) { double rrr = pc.cropAbs10(resid[n], -14); writelogf(" % -14.7E % -14.7E % -10.5E\n", elmFracGoal[n], elmFrac[n], rrr); diff --git a/src/equil/vcs_DoubleStarStar.cpp b/src/equil/vcs_DoubleStarStar.cpp index c044a6780..2e0606d27 100644 --- a/src/equil/vcs_DoubleStarStar.cpp +++ b/src/equil/vcs_DoubleStarStar.cpp @@ -3,7 +3,7 @@ * * Header file for class DoubleStarStar */ -#include "vcs_DoubleStarStar.h" +#include "cantera/equil/vcs_DoubleStarStar.h" namespace VCSnonideal { diff --git a/src/equil/vcs_Gibbs.cpp b/src/equil/vcs_Gibbs.cpp index b169e040b..62735d0ad 100644 --- a/src/equil/vcs_Gibbs.cpp +++ b/src/equil/vcs_Gibbs.cpp @@ -8,9 +8,9 @@ * U.S. Government retains certain rights in this software. */ -#include "vcs_solve.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include "math.h" namespace VCSnonideal diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index 0641aa249..a03e9739e 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -9,14 +9,13 @@ */ #include "cantera/equil/vcs_MultiPhaseEquil.h" -#include "vcs_prob.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_prob.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_species_thermo.h" #include "vcs_SpeciesProperties.h" -#include "vcs_VolPhase.h" -#include "vcs_solve.h" +#include "cantera/equil/vcs_solve.h" #include "cantera/base/ct_defs.h" #include "cantera/thermo/mix_defs.h" diff --git a/src/equil/vcs_SpeciesProperties.cpp b/src/equil/vcs_SpeciesProperties.cpp index b6b32e7c6..386d7aad7 100644 --- a/src/equil/vcs_SpeciesProperties.cpp +++ b/src/equil/vcs_SpeciesProperties.cpp @@ -3,9 +3,9 @@ */ #include "cantera/equil/vcs_defs.h" #include "vcs_SpeciesProperties.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_species_thermo.h" -#include "vcs_internal.h" +#include "cantera/equil/vcs_internal.h" #include #include diff --git a/src/equil/vcs_TP.cpp b/src/equil/vcs_TP.cpp index cc1c7de8c..eef7ded8f 100644 --- a/src/equil/vcs_TP.cpp +++ b/src/equil/vcs_TP.cpp @@ -1,7 +1,7 @@ -#include "vcs_solve.h" -#include "vcs_internal.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" #include "vcs_species_thermo.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_VolPhase.h" #include #include diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 270fb2acb..ecca50f80 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -6,14 +6,15 @@ * Contract DE-AC04-94AL85000 with Sandia Corporation, the * U.S. Government retains certain rights in this software. */ -#include "vcs_VolPhase.h" -#include "vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" +#include "cantera/equil/vcs_internal.h" #include "vcs_SpeciesProperties.h" #include "vcs_species_thermo.h" -#include "vcs_solve.h" +#include "cantera/equil/vcs_solve.h" #include "cantera/thermo/ThermoPhase.h" #include "cantera/thermo/mix_defs.h" +#include "vcs_Exception.h" #include #include diff --git a/src/equil/vcs_dbocls.c b/src/equil/vcs_dbocls.c index a7880f8f8..33fdf9127 100644 --- a/src/equil/vcs_dbocls.c +++ b/src/equil/vcs_dbocls.c @@ -18,6 +18,7 @@ static integer c__1 = 1; static integer c__0 = 0; static integer c__2 = 2; static integer c__4 = 4; +static long int lc__4 = 4; /* DECK DBOCLS */ /* Subroutine */ @@ -67,7 +68,7 @@ int dbocls_(doublereal* w, integer* mdw, integer* mcon, extern /* Subroutine */ int dcopy_(integer*, doublereal*, integer*, doublereal*, integer*); static integer liopt; - extern doublereal d1mach_(integer*); + extern doublereal d1mach_(long int *); static integer locacc; static logical checkl; static integer iscale, locdim; @@ -717,7 +718,7 @@ int dbocls_(doublereal* w, integer* mdw, integer* mcon, /* PROCEDURE(PROCESS OPTION ARRAY) */ zero = 0.; one = 1.; - drelpr = d1mach_(&c__4); + drelpr = d1mach_(&lc__4); checkl = FALSE_; filter = TRUE_; /*lenx = (*ncols + *mcon << 1) + 2;*/ diff --git a/src/equil/vcs_dbolsm.c b/src/equil/vcs_dbolsm.c index 019c33385..5e8f03033 100644 --- a/src/equil/vcs_dbolsm.c +++ b/src/equil/vcs_dbolsm.c @@ -17,8 +17,10 @@ static integer c__1 = 1; static integer c__0 = 0; static integer c__2 = 2; +static long int lc__2 = 2; static integer c_n4 = -4; static integer c__4 = 4; +static long int lc__4 = 4; /* DECK DBOLSM */ @@ -74,7 +76,7 @@ int dbolsm_(doublereal* w, integer* mdw, integer* minput, extern /* Subroutine */ int dmout_(integer*, integer*, integer*, doublereal*, char*, integer*, ftnlen); static integer jdrop; - extern doublereal d1mach_(integer*); + extern doublereal d1mach_(long int*); extern /* Subroutine */ int dvout_(integer*, doublereal*, char*, integer*, ftnlen), ivout_(integer*, integer*, char*, integer * , ftnlen); @@ -1315,7 +1317,7 @@ L470: /* PROGRAM UNIT. */ /* THE COL. SCALING IS DEFINED IN THE CALLING PROGRAM UNIT. */ /* 'BIG' IS PLUS INFINITY ON THIS MACHINE. */ - big = d1mach_(&c__2); + big = d1mach_(&lc__2); i__1 = *ncols; for (j = 1; j <= i__1; ++j) { icase = ind[j]; @@ -1388,8 +1390,8 @@ L580: fac = .75; one = 1.; two = 2.; - tolind = sqrt(d1mach_(&c__4)); - tolsze = sqrt(d1mach_(&c__4)); + tolind = sqrt(d1mach_(&lc__4)); + tolsze = sqrt(d1mach_(&lc__4)); itmax = max(mrows,*ncols) * 5; wt = one; mval = 0; @@ -1437,12 +1439,12 @@ L590: goto L610; } tolind = x[*ncols + ioff]; - if (tolind < d1mach_(&c__4)) { + if (tolind < d1mach_(&lc__4)) { nerr = 25; nlevel = 0; nchar = 88; rdum2 = (real) tolind; - rdum = (real) d1mach_(&c__4); + rdum = (real) d1mach_(&lc__4); xerrwv_("DBOLSM(). THE TOLERANCE FOR RANK DETERMINATION=(R1)" " IS LESS THAN MACHINE PRECISION=(R2).", &nchar, &nerr, &nlevel, &c__0, &idum, &idum, &c__2, &rdum2, &rdum, ( diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index 3b36af3b8..be09de131 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -4,8 +4,9 @@ * element abundances constraints and the algorithm for fixing violations * of the element abundances constraints. */ -#include "vcs_solve.h" -#include "vcs_internal.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "vcs_Exception.h" #include "math.h" namespace VCSnonideal @@ -217,11 +218,11 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) plogf("\n"); } - for (i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_numElemConstraints; ++i) { x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i]; } l2before = 0.0; - for (i = 0; i < m_numElemConstraints; ++i) { + for (size_t i = 0; i < m_numElemConstraints; ++i) { l2before += x[i] * x[i]; } l2before = sqrt(l2before/m_numElemConstraints); @@ -585,14 +586,14 @@ L_CLEANUP: vcs_tmoles(); #ifdef DEBUG_MODE l2after = 0.0; - for (i = 0; i < m_numElemConstraints; ++i) { + for (int i = 0; i < m_numElemConstraints; ++i) { l2after += SQUARE(m_elemAbundances[i] - m_elemAbundancesGoal[i]); } l2after = sqrt(l2after/m_numElemConstraints); if (m_debug_print_lvl >= 2) { plogf(" --- Elem_Abund: Correct Initial " " Final\n"); - for (i = 0; i < m_numElemConstraints; ++i) { + for (int i = 0; i < m_numElemConstraints; ++i) { plogf(" --- "); plogf("%-2.2s", m_elementName[i].c_str()); plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]); diff --git a/src/equil/vcs_elem_rearrange.cpp b/src/equil/vcs_elem_rearrange.cpp index 3be6c86f9..29733da09 100644 --- a/src/equil/vcs_elem_rearrange.cpp +++ b/src/equil/vcs_elem_rearrange.cpp @@ -9,9 +9,9 @@ * U.S. Government retains certain rights in this software. */ -#include "vcs_solve.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include #include @@ -258,7 +258,7 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos) } } std::swap(m_elemAbundancesGoal[ipos], m_elemAbundancesGoal[jpos]); - std::swap(m_elemAbundances[ipos], m_elemAbundancesGoal[jpos]); + std::swap(m_elemAbundances[ipos], m_elemAbundances[jpos]); std::swap(m_elementMapIndex[ipos], m_elementMapIndex[jpos]); std::swap(m_elType[ipos], m_elType[jpos]); std::swap(m_elementActive[ipos], m_elementActive[jpos]); diff --git a/src/equil/vcs_equilibrate.cpp b/src/equil/vcs_equilibrate.cpp index 7ba1886cb..17211f9ca 100644 --- a/src/equil/vcs_equilibrate.cpp +++ b/src/equil/vcs_equilibrate.cpp @@ -8,13 +8,13 @@ * U.S. Government retains certain rights in this software. */ #include "cantera/equil/vcs_MultiPhaseEquil.h" -#include "vcs_prob.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_prob.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_species_thermo.h" #include "vcs_SpeciesProperties.h" -#include "vcs_VolPhase.h" -#include "vcs_solve.h" +#include "cantera/equil/vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" #include "cantera/equil/equil.h" #include "cantera/base/ct_defs.h" diff --git a/src/equil/vcs_inest.cpp b/src/equil/vcs_inest.cpp index dfd3d7434..3076a22f2 100644 --- a/src/equil/vcs_inest.cpp +++ b/src/equil/vcs_inest.cpp @@ -8,9 +8,9 @@ * U.S. Government retains certain rights in this software. */ -#include "vcs_solve.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include "cantera/base/clockWC.h" diff --git a/src/equil/vcs_nondim.cpp b/src/equil/vcs_nondim.cpp index 4d935224f..5b6029dac 100644 --- a/src/equil/vcs_nondim.cpp +++ b/src/equil/vcs_nondim.cpp @@ -8,10 +8,11 @@ * U.S. Government retains certain rights in this software. */ -#include "vcs_solve.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include "cantera/base/stringUtils.h" +#include "vcs_Exception.h" #include #include diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index 2abb3987e..682aee9af 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -3,10 +3,11 @@ * Implementation class for functions associated with determining the stability of a phase * (see Class \link Cantera::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ). */ -#include "vcs_solve.h" -#include "vcs_internal.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" #include "vcs_species_thermo.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_VolPhase.h" +#include "vcs_Exception.h" #include #include @@ -58,7 +59,7 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const size_t kspec = Vphase->spGlobalIndexVCS(k); #ifdef DEBUG_MODE if (m_molNumSpecies_old[kspec] > 0.0) { - printf("ERROR vcs_popPhasePossible we shouldn't be here %d %g > 0.0", + printf("ERROR vcs_popPhasePossible we shouldn't be here %lu %g > 0.0", kspec, m_molNumSpecies_old[kspec]); exit(-1); } @@ -505,14 +506,14 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) if (m_molNumSpecies_old[j] > 0.0) { #ifdef DEBUG_MODE sprintf(anote, "Delta damped from %g " - "to %g due to component %d (%10s) going neg", m_deltaMolNumSpecies[kspec], + "to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec], -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str()); #endif m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC; } else { #ifdef DEBUG_MODE sprintf(anote, "Delta damped from %g " - "to %g due to component %d (%10s) zero", m_deltaMolNumSpecies[kspec], + "to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec], -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str()); #endif m_deltaMolNumSpecies[kspec] = 0.0; diff --git a/src/equil/vcs_prep.cpp b/src/equil/vcs_prep.cpp index 6dde7318a..770a1da2f 100644 --- a/src/equil/vcs_prep.cpp +++ b/src/equil/vcs_prep.cpp @@ -9,10 +9,10 @@ * U.S. Government retains certain rights in this software. */ -#include "vcs_solve.h" -#include "vcs_internal.h" -#include "vcs_prob.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_prob.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_SpeciesProperties.h" #include @@ -125,7 +125,11 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) * Set an initial estimate for the number of noncomponent species * equal to nspecies - nelements. This may be changed below */ - m_numRxnTot = m_numSpeciesTot - m_numElemConstraints; + if (m_numElemConstraints > m_numSpeciesTot) { + m_numRxnTot = 0; + } else { + m_numRxnTot = m_numSpeciesTot - m_numElemConstraints; + } m_numRxnRdc = m_numRxnTot; m_numSpeciesRdc = m_numSpeciesTot; for (i = 0; i < m_numRxnRdc; ++i) { @@ -218,11 +222,13 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) return retn; } - if (m_numElemConstraints != m_numComponents) { + if (m_numSpeciesTot >= m_numComponents) { m_numRxnTot = m_numRxnRdc = m_numSpeciesTot - m_numComponents; for (i = 0; i < m_numRxnRdc; ++i) { m_indexRxnToSpecies[i] = m_numComponents + i; } + } else { + m_numRxnTot = m_numRxnRdc = 0; } /* diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index e2205fcef..53416749b 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -9,10 +9,10 @@ * U.S. Government retains certain rights in this software. */ -#include "vcs_prob.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_prob.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_species_thermo.h" -#include "vcs_internal.h" +#include "cantera/equil/vcs_internal.h" #include "cantera/thermo/ThermoPhase.h" #include "cantera/thermo/MolalityVPSSTP.h" diff --git a/src/equil/vcs_rearrange.cpp b/src/equil/vcs_rearrange.cpp index 6c34fcfc6..18777d1f1 100644 --- a/src/equil/vcs_rearrange.cpp +++ b/src/equil/vcs_rearrange.cpp @@ -8,8 +8,8 @@ * U.S. Government retains certain rights in this software. */ -#include "vcs_solve.h" -#include "vcs_internal.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" #include #include diff --git a/src/equil/vcs_report.cpp b/src/equil/vcs_report.cpp index 0ceacfb79..6410b2553 100644 --- a/src/equil/vcs_report.cpp +++ b/src/equil/vcs_report.cpp @@ -6,9 +6,9 @@ */ -#include "vcs_solve.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include #include diff --git a/src/equil/vcs_root1d.cpp b/src/equil/vcs_root1d.cpp index 2ce141dc7..848ce999a 100644 --- a/src/equil/vcs_root1d.cpp +++ b/src/equil/vcs_root1d.cpp @@ -9,7 +9,7 @@ */ -#include "vcs_internal.h" +#include "cantera/equil/vcs_internal.h" #include #include diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index 692a9acd7..971cdb729 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -8,9 +8,9 @@ * U.S. Government retains certain rights in this software. */ -#include "vcs_solve.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include #include @@ -53,13 +53,13 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) char ANOTE[128]; if (m_debug_print_lvl >= 2) { plogf(" "); - for (j = 0; j < 82; j++) { + for (int j = 0; j < 82; j++) { plogf("-"); } plogf("\n"); plogf(" --- Subroutine vcs_RxnStepSizes called - Details:\n"); plogf(" "); - for (j = 0; j < 82; j++) { + for (int j = 0; j < 82; j++) { plogf("-"); } plogf("\n"); diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index 248c3f864..10e14e190 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -8,10 +8,10 @@ * U.S. Government retains certain rights in this software. */ -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_species_thermo.h" -#include "vcs_solve.h" +#include "cantera/equil/vcs_solve.h" #include #include @@ -24,6 +24,10 @@ #ifdef DEBUG_MODE //extern int vcs_debug_print_lvl; #endif +#ifndef MAX +#define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif + extern "C" void dbocls_(double* W, int* MDW, int* MCON, int* MROWS, int* NCOLS, diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 771831e44..24952e5a3 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -9,12 +9,12 @@ */ -#include "vcs_solve.h" +#include "cantera/equil/vcs_solve.h" #include "vcs_Exception.h" -#include "vcs_internal.h" -#include "vcs_prob.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_prob.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_SpeciesProperties.h" #include "vcs_species_thermo.h" @@ -39,6 +39,7 @@ VCS_SOLVE::VCS_SOLVE() : m_numComponents(0), m_numRxnTot(0), m_numSpeciesRdc(0), + m_numRxnRdc(0), m_numRxnMinorZeroed(0), m_numPhases(0), m_doEstimateEquil(0), @@ -485,19 +486,26 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) /* * OK, We have room. Now, transfer the integer numbers */ - m_numElemConstraints = nelements; + m_numElemConstraints = nelements; m_numSpeciesTot = nspecies; m_numSpeciesRdc = m_numSpeciesTot; /* * nc = number of components -> will be determined later. * but set it to its maximum possible value here. */ - m_numComponents = nelements; + m_numComponents = nelements; /* * m_numRxnTot = number of noncomponents, also equal to the * number of reactions + * Note, it's possible that the number of elements is greater than + * the number of species. In that case set the number of reactions + * to zero. */ - m_numRxnTot = std::max(nspecies - nelements, 0); + if (nelements > nspecies) { + m_numRxnTot = 0; + } else { + m_numRxnTot = nspecies - nelements; + } m_numRxnRdc = m_numRxnTot; /* * number of minor species rxn -> all species rxn are major at the start. diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 6d0b3857a..3a280a995 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -13,16 +13,21 @@ #include #include -#include "vcs_solve.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_species_thermo.h" +#include "vcs_Exception.h" #include "cantera/base/clockWC.h" #include "cantera/base/stringUtils.h" using namespace std; +#ifndef MAX +#define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#endif + namespace VCSnonideal { @@ -112,7 +117,6 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) bool justDeletedMultiPhase = false; bool usedZeroedSpecies; /* return flag from basopt indicating that one of the components had a zero concentration */ - size_t doPhaseDeleteIph = npos; vcs_VolPhase* Vphase; double* sc_irxn = NULL; /* Stoichiometric coefficients for cur rxn */ double* dnPhase_irxn; @@ -121,6 +125,9 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) int forceComponentCalc = 1; size_t iphaseDelete; /* integer that determines which phase is being deleted */ std::vector phasePopPhaseIDs(0); + int doPhaseDeleteIph = -1; + int doPhaseDeleteKspec = -1; + #ifdef DEBUG_MODE char ANOTE[128]; /* @@ -158,7 +165,9 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) solveFail = false; - +#if DEBUG_MODE + int ll; +#endif /* ****************************************************** */ /* **** Evaluate the elemental composition ****** */ /* ****************************************************** */ @@ -454,7 +463,8 @@ L_MAINLOOP_ALL_SPECIES: } #endif lec = false; - doPhaseDeleteIph = npos; + doPhaseDeleteIph = -1; + doPhaseDeleteKspec = -1; /* * Zero out the net change in moles of multispecies phases */ @@ -901,6 +911,7 @@ L_MAINLOOP_ALL_SPECIES: */ m_molNumSpecies_new[kspec] = 0.0; doPhaseDeleteIph = iph; + doPhaseDeleteKspec = kspec; #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { @@ -1016,7 +1027,7 @@ L_MAIN_LOOP_END: L_MAIN_LOOP_END_NO_PRINT: ; #endif - if (doPhaseDeleteIph != npos) { + if (doPhaseDeleteIph != -1) { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { plogf(" --- "); @@ -3670,7 +3681,7 @@ L_END_LOOP: double sumMax = -1.0; int iMax = -1; int jMax = -1; - int n; + size_t n; for (i = 0; i < m_numRxnTot; ++i) { k = m_indexRxnToSpecies[i]; for (j = 0; j < ncTrial; ++j) { @@ -4336,7 +4347,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old); actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old); molNum = VCS_DATA_PTR(m_molNumSpecies_old); - } else { // stateCalc == VCS_STATECALC_NEW + } else if (stateCalc == VCS_STATECALC_NEW) { feSpecies = VCS_DATA_PTR(m_feSpecies_new); tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new); actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new); @@ -4680,12 +4691,12 @@ void VCS_SOLVE::vcs_printSpeciesChemPot(const int stateCalc) const //! Print out and check the elemental abundance vector void VCS_SOLVE::prneav() const { - int j; + size_t j; bool kerr; std::vector eav(m_numElemConstraints, 0.0); for (j = 0; j < m_numElemConstraints; ++j) { - for (int i = 0; i < m_numSpeciesTot; ++i) { + for (size_t i = 0; i < m_numSpeciesTot; ++i) { if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { eav[j] += m_formulaMatrix[j][i] * m_molNumSpecies_old[i]; } @@ -4781,7 +4792,7 @@ double VCS_SOLVE::vcs_tmoles() #ifdef DEBUG_MODE void VCS_SOLVE::check_tmoles() const { - int i; + size_t i; double sum = 0.0; for (i = 0; i < m_numPhases; i++) { double m_tPhaseMoles_old_a = TPhInertMoles[i]; @@ -5400,7 +5411,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, feSpecies = VCS_DATA_PTR(m_feSpecies_new); deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_new); actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_new); - } else { // stateCalc == VCS_STATECALC_OLD + } else if (stateCalc == VCS_STATECALC_OLD) { feSpecies = VCS_DATA_PTR(m_feSpecies_old); deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old); actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_old); diff --git a/src/equil/vcs_solve_phaseStability.cpp b/src/equil/vcs_solve_phaseStability.cpp index dc321ac36..0f73f1900 100644 --- a/src/equil/vcs_solve_phaseStability.cpp +++ b/src/equil/vcs_solve_phaseStability.cpp @@ -16,11 +16,11 @@ #include #include -#include "vcs_solve.h" -#include "vcs_internal.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_solve.h" +#include "cantera/equil/vcs_internal.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_species_thermo.h" -#include "vcs_prob.h" +#include "cantera/equil/vcs_prob.h" #include "cantera/base/clockWC.h" diff --git a/src/equil/vcs_species_thermo.cpp b/src/equil/vcs_species_thermo.cpp index d1d61cff7..cc4926df1 100644 --- a/src/equil/vcs_species_thermo.cpp +++ b/src/equil/vcs_species_thermo.cpp @@ -10,13 +10,13 @@ -#include "vcs_solve.h" +#include "cantera/equil/vcs_solve.h" #include "vcs_species_thermo.h" #include "cantera/equil/vcs_defs.h" -#include "vcs_VolPhase.h" +#include "cantera/equil/vcs_VolPhase.h" #include "vcs_Exception.h" -#include "vcs_internal.h" +#include "cantera/equil/vcs_internal.h" #include #include diff --git a/src/equil/vcs_species_thermo.h b/src/equil/vcs_species_thermo.h index 4726a8a0d..24f7364d5 100644 --- a/src/equil/vcs_species_thermo.h +++ b/src/equil/vcs_species_thermo.h @@ -9,6 +9,7 @@ #define VCS_SPECIES_THERMO_H //#include +#include namespace VCSnonideal { diff --git a/src/equil/vcs_util.cpp b/src/equil/vcs_util.cpp index 33528e756..db3296966 100644 --- a/src/equil/vcs_util.cpp +++ b/src/equil/vcs_util.cpp @@ -11,8 +11,8 @@ #include #include -#include "vcs_internal.h" -#include +#include "cantera/equil/vcs_internal.h" +#include #include using namespace std; diff --git a/src/fortran/SConscript b/src/fortran/SConscript index 7d5f39a8e..c2c3b91d0 100644 --- a/src/fortran/SConscript +++ b/src/fortran/SConscript @@ -17,3 +17,6 @@ lib = build(localenv.Library(target=pjoin('..','..','lib','cantera_fortran'), install('$inst_libdir', lib) install('$inst_incdir', mods) +# +# We are only building the static fortran library here. This is probably ok +# diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 077c7b0db..052b4cb84 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -8,7 +8,9 @@ */ // Copyright 2001-2004 California Institute of Technology -#include "cantera/kinetics/InterfaceKinetics.h" +// Why InterfaceKinetics.h and not Kinetics.h ?? + +#include "cantera/kinetics/Kinetics.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/kinetics/StoichManager.h" #include "cantera/kinetics/RateCoeffMgr.h" diff --git a/src/kinetics/ReactionStoichMgr.cpp b/src/kinetics/ReactionStoichMgr.cpp index 2b8a05d0f..0dcf815a1 100644 --- a/src/kinetics/ReactionStoichMgr.cpp +++ b/src/kinetics/ReactionStoichMgr.cpp @@ -6,6 +6,8 @@ //------------------------------------------------ #include "cantera/kinetics/ReactionStoichMgr.h" + +#include "cantera/kinetics/StoichManager.h" #include "cantera/base/ctexceptions.h" #include "cantera/kinetics/ReactionData.h" diff --git a/src/matlab/SConscript b/src/matlab/SConscript index dbbf214bc..f8ffd328a 100644 --- a/src/matlab/SConscript +++ b/src/matlab/SConscript @@ -9,7 +9,8 @@ matlab_include = pjoin(localenv['matlab_path'], 'extern', 'include') if localenv['OS'] == 'Windows': mexPlatform = 'w' - linklibs = ['cantera_shared'] + linklibs = [] + linkLibs += env['cantera_shared_libs'] linklibs += ['libmx', 'libmex', 'libmat'] if localenv['OS_BITS'] == 32: matlab_libs = pjoin(localenv['matlab_path'], 'extern', @@ -26,7 +27,8 @@ if localenv['OS'] == 'Windows': linkflags.extend(['-static-libgcc', '-static-libstdc++']) elif localenv['OS'] == 'Darwin': - linklibs = ['cantera'] + linklibs = [] + linkLibs += env['cantera_libs'] linklibs += ['mx', 'mex', 'mat'] + env['LIBM'] mexPlatform = 'maci' linkflags.extend(['-Wl,-exported_symbol,_mexFunction']) @@ -37,7 +39,8 @@ elif localenv['OS'] == 'Darwin': matlab_libs = pjoin(localenv['matlab_path'], 'bin', 'macx86') elif os.name == 'posix': - linklibs = ['cantera'] + linklibs = [] + linklibs += env['cantera_libs'] linklibs += ['mx', 'mex', 'mat'] + env['LIBM'] mexPlatform = 'a' diff --git a/src/numerics/BEulerInt.cpp b/src/numerics/BEulerInt.cpp index b7f12e02a..e71642e29 100644 --- a/src/numerics/BEulerInt.cpp +++ b/src/numerics/BEulerInt.cpp @@ -10,8 +10,7 @@ * See file License.txt for licensing information. */ -#include "BEulerInt.h" - +#include "cantera/numerics/BEulerInt.h" #include "cantera/base/mdp_allo.h" #include diff --git a/src/numerics/NonlinearSolver.cpp b/src/numerics/NonlinearSolver.cpp index eb016aa08..8925523a2 100644 --- a/src/numerics/NonlinearSolver.cpp +++ b/src/numerics/NonlinearSolver.cpp @@ -3744,11 +3744,12 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, */ doublereal* dyVector = mdp::mdp_alloc_dbl_1((int) neq_, MDP_DBL_NOINIT); retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt)); - - - if (s_print_NumJac) { if (m_print_flag >= 7) { + if (retn != 1) { + printf("\t\tbeuler_jac ERROR: calcDeltaSolnVariables() returned an error condition.\n"); + printf("\t\t We will bail after calculating the Jacobian\n"); + } if (neq_ < 20) { printf("\t\tUnk m_ewt y dyVector ResN\n"); for (size_t iii = 0; iii < neq_; iii++) { @@ -3841,6 +3842,10 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt)); if (s_print_NumJac) { if (m_print_flag >= 7) { + if (retn != 1) { + printf("\t\tbeuler_jac ERROR: calcDeltaSolnVariables() returned an error condition.\n"); + printf("\t\t We will bail after calculating the Jacobian\n"); + } if (neq_ < 20) { printf("\t\tUnk m_ewt y dyVector ResN\n"); for (size_t iii = 0; iii < neq_; iii++) { diff --git a/src/numerics/ResidJacEval.cpp b/src/numerics/ResidJacEval.cpp index 29194cc04..9d80ff011 100644 --- a/src/numerics/ResidJacEval.cpp +++ b/src/numerics/ResidJacEval.cpp @@ -312,7 +312,7 @@ matrixConditioning(doublereal* const matrix, const int nrows, doublereal* const int ResidJacEval:: evalResidNJ(const doublereal t, const doublereal deltaT, const doublereal* y, const doublereal* ydot, doublereal* const resid, const ResidEval_Type_Enum evalType, - const size_t id_x, const doublereal delta_x) + const int id_x, const doublereal delta_x) { throw CanteraError("ResidJacEval::evalResidNJ()", "Not implemented\n"); return 1; diff --git a/src/python/SConscript b/src/python/SConscript index 81ddeb011..ec188f9de 100644 --- a/src/python/SConscript +++ b/src/python/SConscript @@ -31,7 +31,11 @@ if localenv['python_package'] == 'full': localenv.Append(CPPPATH=['#src', '#include']) cantera_libname = 'cantera_shared' if os.name=='nt' else 'cantera' - pylinklibs = [cantera_libname] + #pylinklibs = [cantera_libname, 'cvode', 'ctmath', 'ctlapack', 'ctblas', 'ctf2c', 'execstream'] + if os.name == 'nt': + pylinklibs = env['cantera_shared_libs'] + else: + pylinklibs = env['cantera_libs'] if localenv['toolchain'] == 'mingw': # On Windows, we need to link against the Python "import" library. diff --git a/src/thermo/Elements.cpp b/src/thermo/Elements.cpp index c8a82b1aa..4f58e184f 100644 --- a/src/thermo/Elements.cpp +++ b/src/thermo/Elements.cpp @@ -2,17 +2,38 @@ * @file Elements.cpp * This file contains a database of atomic weights. */ + +/**************************************************************************** + * $RCSfile: Elements.cpp,v $ + * $Author$ + * $Date$ + * $Revision$ + * + * + ****************************************************************************/ // Copyright 2003 California Institute of Technology -#include "cantera/thermo/Elements.h" -#include "cantera/base/ctexceptions.h" +#ifdef WIN32 +#pragma warning(disable:4786) +#endif + +#include "cantera/thermo/Elements.h" +#include "cantera/base/xml.h" +#include "cantera/base/ctml.h" +#include "cantera/base/ctexceptions.h" +#include "cantera/base/stringUtils.h" + +using namespace ctml; using namespace std; +#ifdef USE_DGG_CODE +#include +#endif #include -namespace Cantera -{ +namespace Cantera { + /*! Database for atomic molecular weights * Values are taken from the 1989 Standard Atomic Weights, CRC @@ -29,13 +50,13 @@ struct awData { double atomicWeight; //!< atomic weight in kg / kg-mol }; -/*! - * @var static struct awData aWTable[] - * \brief aWTable is a vector containing the atomic weights database. - * - * The size of the table is given by the initial instantiation. - */ -static struct awData aWTable[] = { + /*! + * @var static struct awData aWTable[] + * \brief aWTable is a vector containing the atomic weights database. + * + * The size of the table is given by the initial instantiation. + */ + static struct awData aWTable[] = { {"H", 1.00794}, {"D", 2.0 }, {"Tr", 3.0 }, @@ -132,19 +153,511 @@ static struct awData aWTable[] = { {"U", 238.0508 }, {"Np", 237.0482 }, {"Pu", 244.0482 } -}; + }; + -doublereal LookupWtElements(const std::string& ename) -{ + // Static function to look up an atomic weight + /* + * This static function looks up the argument string in the + * database above and returns the associated molecular weight. + * The data are from the periodic table. + * + * Note: The idea behind this function is to provide a unified + * source for the element atomic weights. This helps to + * ensure that mass is conserved. + * + * @param s String, Only the first 3 characters are significant + * + * @return + * Return value contains the atomic weight of the element + * If a match for the string is not found, a value of -1.0 is + * returned. + * + * @exception CanteraError + * If a match is not found, a CanteraError is thrown as well + */ + doublereal Elements::LookupWtElements(const std::string& ename) { int num = sizeof(aWTable) / sizeof(struct awData); string s3 = ename.substr(0,3); for (int i = 0; i < num; i++) { - if (s3 == aWTable[i].name) { - return (aWTable[i].atomicWeight); - } + if (s3 == aWTable[i].name) { + return (aWTable[i].atomicWeight); + } } throw CanteraError("LookupWtElements", "element not found"); return -1.0; -} + } + doublereal LookupWtElements(const std::string& ename) { + int num = sizeof(aWTable) / sizeof(struct awData); + string s3 = ename.substr(0,3); + for (int i = 0; i < num; i++) { + if (s3 == aWTable[i].name) { + return (aWTable[i].atomicWeight); + } + } + throw CanteraError("LookupWtElements", "element not found"); + return -1.0; + } + + + + + + //! Exception class to indicate a fixed set of elements. + /*! + * This class is used to warn the user when the number of elements + * are changed after at least one species is defined. + */ + class ElementsFrozen : public CanteraError { + public: + //! Constructor for class + /*! + * @param func Function where the error occurred. + */ + ElementsFrozen(string func) + : CanteraError(func, + "elements cannot be added after species.") {} + }; + + /* + * Elements Class Constructor + * We initialize all internal variables to zero here. + */ + Elements::Elements() : + m_mm(0), + m_elementsFrozen(false), + m_elem_type(0), + numSubscribers(0) + { + } + + /* + * Elements Class Destructor + * If the number of subscribers is not zero, through an error. + * A logic problem has occurred. + * + * @exception CanteraError + */ + Elements::~Elements() { + if (numSubscribers != 0) { + throw CanteraError("~Elements", "numSubscribers not zero"); + } + } + + Elements::Elements(const Elements &right) : + m_mm(0), + m_elementsFrozen(false), + numSubscribers(0) + { + *this = operator=(right); + } + + Elements& Elements::operator=(const Elements &right) { + if (&right == this) return *this; + + m_mm = right.m_mm; + m_elementsFrozen = right.m_elementsFrozen; + m_atomicWeights = right.m_atomicWeights; + m_atomicNumbers = right.m_atomicNumbers; + m_elementNames = right.m_elementNames; + m_entropy298 = right.m_entropy298; + m_elem_type = right.m_elem_type; + numSubscribers = 0; + + return *this; + } + + + + /* + * freezeElements(): + * + * Set the freeze flag. This is a prerequesite to other + * activivities, i.e., this is done before species are defined. + */ + void Elements::freezeElements() { + m_elementsFrozen = true; + } + +#ifdef INCL_DEPRECATED_METHODS + /* + * + * Returns an ElementData struct that contains the parameters + * for element index m. + */ + ElementData Elements::element(int m) const { + ElementData e; + e.name = m_elementNames[m]; + e.atomicWeight = m_atomicWeights[m]; + return e; + } +#endif + /* + * elementIndex(): + * + * Index of element named \c name. The index is an integer + * assigned to each element in the order it was added, + * beginning with 0 for the first element. If \c name is not + * the name of an element in the set, then the value -1 is + * returned. + * + */ +#ifdef USE_DGG_CODE + int Elements::elementIndex(std::string name) const{ + map::const_iterator it; + it = m_definedElements.find(name); + if (it != m_definedElements.end()) { + return it->second; + } + return -1; + } +#else + int Elements::elementIndex(std::string name) const { + for (int i = 0; i < m_mm; i++) { + if (m_elementNames[i] == name) return i; + } + return -1; + } +#endif + + /* + * + * Name of the element with index \c m. @param m Element + * index. If m < 0 or m >= nElements() an exception is thrown. + */ + string Elements::elementName(int m) const { + if (m < 0 || m >= nElements()) { + throw CanteraError("Elements::elementName()", "out of bounds: " + int2str(m) + " " + int2str(nElements())); + } + return m_elementNames[m]; + } + + + doublereal Elements::entropyElement298(int m) const { + AssertThrowMsg(m_entropy298[m] != ENTROPY298_UNKNOWN, + "Elements::entropy298", + "Entropy at 298 K of element is unknown"); + AssertTrace(m >= 0 && m < m_mm); + return (m_entropy298[m]); + } + //==================================================================================================================== + //! Return the element constraint type + /*! + * Possible types include: + * + * CT_ELEM_TYPE_TURNEDOFF -1 + * CT_ELEM_TYPE_ABSPOS 0 + * CT_ELEM_TYPE_ELECTRONCHARGE 1 + * CT_ELEM_TYPE_CHARGENEUTRALITY 2 + * CT_ELEM_TYPE_LATTICERATIO 3 + * CT_ELEM_TYPE_KINETICFROZEN 4 + * CT_ELEM_TYPE_SURFACECONSTRAINT 5 + * CT_ELEM_TYPE_OTHERCONSTRAINT 6 + * + * The default is CT_ELEM_TYPE_ABSPOS + */ + int Elements::elementType(int m) const + { + return m_elem_type[m]; + } + //==================================================================================================================== + // Change the element type of the mth constraint + /* + * Reassigns an element type + * + * @param m Element index + * @param elem_type New elem type to be assigned + * + * @return Returns the old element type + */ + int Elements::changeElementType(int m, int elem_type) + { + int old = m_elem_type[m]; + m_elem_type[m] = elem_type; + return old; + } + //==================================================================================================================== + /* + * + * Add an element to the current set of elements in the current object. + * @param symbol symbol string + * @param weight atomic weight in kg/kmol. + * + * The default weight is a special value, which will cause the + * routine to look up the actual weight via a string lookup. + * + * There are two interfaces to this routine. The XML interface + * looks up the required parameters for the regular interface + * and then calls the base routine. + */ + void Elements:: + addElement(const std::string& symbol, doublereal weight) + { + if (weight == -12345.0) { + weight = LookupWtElements(symbol); + if (weight < 0.0) { + throw ElementsFrozen("addElement"); + } + } + if (m_elementsFrozen) { + throw ElementsFrozen("addElement"); + return; + } + m_atomicWeights.push_back(weight); + m_elementNames.push_back(symbol); +#ifdef USE_DGG_CODE + m_definedElements[symbol] = nElements() + 1; +#endif + if (symbol == "E") { + m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE); + } else { + m_elem_type.push_back(CT_ELEM_TYPE_ABSPOS); + } + + m_mm++; + } + //=========================================================================================================== + void Elements:: + addElement(const XML_Node& e) { + doublereal weight = atof(e["atomicWt"].c_str()); + string symbol = e["name"]; + addElement(symbol, weight); + } + //=========================================================================================================== + /* + * addUniqueElement(): + * + * Add a unique element to the set. This routine will not allow + * duplicate elements to be input. + * + * @param symbol symbol string + * @param weight atomic weight in kg/kmol. + * + * + * The default weight is a special value, which will cause the + * routine to look up the actual weight via a string lookup. + */ +#ifdef USE_DGG_CODE + void Elements:: + addUniqueElement(const std::string& symbol, doublereal weight, int atomicNumber, + doublereal entropy298, int elem_type) + { + if (m_elementsFrozen) + throw ElementsFrozen("addElement"); + + if (weight == -12345.0) { + weight = LookupWtElements(symbol); + } + + /* + * First decide if this element has been previously added. + * If it unique, add it to the list. + */ + + int i = m_definedElements[symbol] - 1; + if (i < 0) { + m_atomicWeights.push_back(weight); + m_elementNames.push_back(symbol); + m_atomicNumbers.push_back(atomicNumber); + m_entropy298.push_back(entropy298); + if (symbol == "E") { + m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE); + } else { + m_elem_type.push_back(elem_type); + } + m_mm++; + } + else { + if (m_atomicWeights[i] != weight) { + throw CanteraError("AddUniqueElement", + "Duplicate Elements (" + symbol + ") have different weights"); + } + } + } + +#else + void Elements:: + addUniqueElement(const std::string& symbol, + doublereal weight, int atomicNumber, doublereal entropy298, + int elem_type) + { + if (weight == -12345.0) { + weight = LookupWtElements(symbol); + if (weight < 0.0) { + throw ElementsFrozen("addElement"); + } + } + /* + * First decide if this element has been previously added + * by conducting a string search. If it unique, add it to + * the list. + */ + int ifound = 0; + int i = 0; + for (vector::const_iterator it = m_elementNames.begin(); + it < m_elementNames.end(); ++it, ++i) { + if (*it == symbol) { + ifound = 1; + break; + } + } + if (!ifound) { + if (m_elementsFrozen) { + throw ElementsFrozen("addElement"); + return; + } + m_atomicWeights.push_back(weight); + m_elementNames.push_back(symbol); + m_atomicNumbers.push_back(atomicNumber); + m_entropy298.push_back(entropy298); + if (symbol == "E") { + m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE); + } else { + m_elem_type.push_back(elem_type); + } + m_mm++; + } else { + if (m_atomicWeights[i] != weight) { + throw CanteraError("AddUniqueElement", + "Duplicate Elements (" + symbol + ") have different weights"); + } + } + } +#endif + + + /* + * @todo call addUniqueElement(symbol, weight) instead of + * addElement. + */ + void Elements:: + addUniqueElement(const XML_Node& e) { + doublereal weight = 0.0; + if (e.hasAttrib("atomicWt")) + weight = atof(stripws(e["atomicWt"]).c_str()); + int anum = 0; + if (e.hasAttrib("atomicNumber")) + anum = atoi(stripws(e["atomicNumber"]).c_str()); + string symbol = e["name"]; + doublereal entropy298 = ENTROPY298_UNKNOWN; + if (e.hasChild("entropy298")) { + XML_Node& e298Node = e.child("entropy298"); + if (e298Node.hasAttrib("value")) { + entropy298 = atofCheck(stripws(e298Node["value"]).c_str()); + } + } + if (weight != 0.0) { + addUniqueElement(symbol, weight, anum, entropy298); + } else { + addUniqueElement(symbol); + } + } + + // True if freezeElements has been called. + bool Elements::elementsFrozen() const + { + return m_elementsFrozen; + } + + /* + * clear() + * + * Remove all elements from the structure. + */ + void Elements::clear() { + m_mm = 0; + m_atomicWeights.resize(0); + m_elementNames.resize(0); + m_entropy298.resize(0); + m_elem_type.resize(0); + m_elementsFrozen = false; + } + + /* + * ready(): + * + * True if the elements have been frozen + */ + bool Elements::ready() const { + return (m_elementsFrozen); + } + + + void Elements::addElementsFromXML(const XML_Node& phase) { + + // get the declared element names + if (! phase.hasChild("elementArray")) { + throw CanteraError("Elements::addElementsFromXML", + "phase xml node doesn't have \"elementArray\" XML Node"); + } + XML_Node& elements = phase.child("elementArray"); + vector enames; + getStringArray(elements, enames); + + // // element database defaults to elements.xml + string element_database = "elements.xml"; + if (elements.hasAttrib("datasrc")) + element_database = elements["datasrc"]; + + XML_Node* doc = get_XML_File(element_database); + XML_Node* dbe = &doc->child("ctml/elementData"); + + XML_Node& root = phase.root(); + XML_Node* local_db = 0; + if (root.hasChild("ctml")) { + if (root.child("ctml").hasChild("elementData")) { + local_db = &root.child("ctml/elementData"); + } + } + + int nel = static_cast(enames.size()); + int i; + string enm; + XML_Node* e = 0; + for (i = 0; i < nel; i++) { + e = 0; + if (local_db) { + //writelog("looking in local database."); + e = local_db->findByAttr("name",enames[i]); + //if (!e) writelog(enames[i]+" not found."); + } + if (!e) + e = dbe->findByAttr("name",enames[i]); + if (e) { + addUniqueElement(*e); + } + else { + throw CanteraError("addElementsFromXML","no data for element " + +enames[i]); + } + } + + } + + /* + * subscribe(), unsubscribe(), and reportSubscriptions(): + * + * Handles setting and reporting the number of subscriptions to this + * object. + */ + void Elements::subscribe() { + ++numSubscribers; + } + int Elements::unsubscribe() { + --numSubscribers; + return numSubscribers; + } + int Elements::reportSubscriptions() const { + return numSubscribers; + } + + /********************* GLOBAL STATIC SECTION **************************/ + /* + * We keep track of a vector of pointers to element objects. + * Initially there are no Elements objects. Whenever one is created, + * the pointer to that object is added onto this list. + */ + vector Elements::Global_Elements_List; + /***********************************************************************/ } diff --git a/src/thermo/IonsFromNeutralVPSSTP.cpp b/src/thermo/IonsFromNeutralVPSSTP.cpp index 39a9ce0ec..30668b6a9 100644 --- a/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -701,7 +701,7 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const } #ifdef DEBUG_MODE sum = -1.0; - for (k = 0; k < m_kk; k++) { + for (int k = 0; k < m_kk; k++) { sum += moleFractions_[k]; } if (fabs(sum) > 1.0E-11) { diff --git a/src/transport/TransportFactory.cpp b/src/transport/TransportFactory.cpp index 0aa13a481..8c65c7e89 100644 --- a/src/transport/TransportFactory.cpp +++ b/src/transport/TransportFactory.cpp @@ -36,6 +36,7 @@ #include #include +#include using namespace std;