From 2a00814a9d25ae0541d3c0abbd625252483c453e Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Sat, 7 Jun 2003 23:53:49 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/clib/src/Makefile.in | 2 +- Cantera/clib/src/ctbdry.cpp | 26 +- Cantera/clib/src/ctonedim.cpp | 461 ++++++++++++++++++ Cantera/clib/src/ctonedim.h | 67 +++ Cantera/clib/src/ctstagn.cpp | 2 +- Cantera/cxx/Makefile.in | 4 +- Cantera/matlab/Makefile.in | 17 +- .../matlab/cantera/1D/@Domain1D/Domain1D.m | 21 + .../cantera/1D/@Domain1D/componentIndex.m | 9 + .../cantera/1D/@Domain1D/componentName.m | 7 + .../cantera/1D/@Domain1D/disableEnergy.m | 4 + .../matlab/cantera/1D/@Domain1D/domainIndex.m | 8 + .../matlab/cantera/1D/@Domain1D/domainType.m | 7 + .../matlab/cantera/1D/@Domain1D/domain_hndl.m | 5 + .../cantera/1D/@Domain1D/enableEnergy.m | 7 + .../matlab/cantera/1D/@Domain1D/gridPoints.m | 15 + Cantera/matlab/cantera/1D/@Domain1D/isFlow.m | 10 + Cantera/matlab/cantera/1D/@Domain1D/isInlet.m | 10 + .../matlab/cantera/1D/@Domain1D/isSurface.m | 10 + .../matlab/cantera/1D/@Domain1D/massFlux.m | 5 + .../cantera/1D/@Domain1D/massFraction.m | 17 + .../matlab/cantera/1D/@Domain1D/nComponents.m | 4 + Cantera/matlab/cantera/1D/@Domain1D/nPoints.m | 5 + .../1D/@Domain1D/private/domain_methods.m | 13 + Cantera/matlab/cantera/1D/@Domain1D/set.m | 88 ++++ .../matlab/cantera/1D/@Domain1D/setBounds.m | 4 + .../1D/@Domain1D/setFixedTempProfile.m | 12 + Cantera/matlab/cantera/1D/@Domain1D/setID.m | 4 + Cantera/matlab/cantera/1D/@Domain1D/setMdot.m | 4 + .../cantera/1D/@Domain1D/setMoleFractions.m | 4 + .../matlab/cantera/1D/@Domain1D/setPressure.m | 5 + .../matlab/cantera/1D/@Domain1D/setProfile.m | 9 + .../cantera/1D/@Domain1D/setTemperature.m | 4 + .../cantera/1D/@Domain1D/setTolerances.m | 17 + .../matlab/cantera/1D/@Domain1D/setupGrid.m | 5 + .../matlab/cantera/1D/@Domain1D/temperature.m | 5 + Cantera/matlab/cantera/1D/@Domain1D/z.m | 15 + Cantera/matlab/cantera/1D/@Stack/Stack.m | 25 + Cantera/matlab/cantera/1D/@Stack/display.m | 8 + .../matlab/cantera/1D/@Stack/domainIndex.m | 7 + Cantera/matlab/cantera/1D/@Stack/grid.m | 6 + .../matlab/cantera/1D/@Stack/plotSolution.m | 13 + .../cantera/1D/@Stack/private/stack_methods.m | 22 + Cantera/matlab/cantera/1D/@Stack/resid.m | 16 + Cantera/matlab/cantera/1D/@Stack/restore.m | 6 + Cantera/matlab/cantera/1D/@Stack/save.m | 10 + Cantera/matlab/cantera/1D/@Stack/saveSoln.m | 14 + .../matlab/cantera/1D/@Stack/setFlatProfile.m | 5 + Cantera/matlab/cantera/1D/@Stack/setProfile.m | 52 ++ .../cantera/1D/@Stack/setRefineCriteria.m | 5 + Cantera/matlab/cantera/1D/@Stack/setValue.m | 20 + Cantera/matlab/cantera/1D/@Stack/solution.m | 23 + Cantera/matlab/cantera/1D/@Stack/solve.m | 5 + Cantera/matlab/cantera/1D/@Stack/subsref.m | 11 + Cantera/matlab/cantera/1D/@Stack/writeStats.m | 10 + Cantera/matlab/cantera/1D/AxiStagnFlow.m | 6 + Cantera/matlab/cantera/1D/AxisymmetricFlow.m | 11 + Cantera/matlab/cantera/1D/Inlet.m | 13 + Cantera/matlab/cantera/1D/Outlet.m | 11 + Cantera/matlab/cantera/1D/Surface.m | 9 + Cantera/matlab/cantera/1D/SymmPlane.m | 11 + Cantera/matlab/cantera/private/ctmethods.cpp | 63 ++- .../matlab/cantera/private/onedimmethods.cpp | 368 ++++++++++++++ Cantera/matlab/setup_matlab.py | 1 + Cantera/python/Cantera/Flow.py | 16 +- Cantera/python/Cantera/OneDim.py | 61 ++- Cantera/python/Cantera/flame.py | 16 +- Cantera/python/Cantera/rxnpath.py | 1 + Cantera/python/examples/mix1.py | 86 ++++ Cantera/python/examples/mix2.py | 93 ++++ Cantera/python/examples/stflame1.py | 1 + Cantera/python/setup.py.in | 2 +- Cantera/src/Flow1D.h | 82 ---- Cantera/src/Jac2.h | 198 -------- Cantera/src/Resid1D.h | 237 --------- Cantera/src/funcs.cpp | 6 +- Cantera/src/importCTML.cpp | 8 - Cantera/src/misc.cpp | 25 +- Cantera/src/oneD/Domain1D.h | 433 ++++++++++++++++ Cantera/src/oneD/Inlet1D.h | 71 ++- Cantera/src/oneD/Jac1D.h | 8 +- Cantera/src/oneD/MultiNewton.cpp | 34 +- Cantera/src/oneD/Newton1D.h | 14 +- Cantera/src/oneD/OneDim.cpp | 44 +- Cantera/src/oneD/OneDim.h | 25 +- Cantera/src/oneD/Sim1D.cpp | 170 +++++-- Cantera/src/oneD/Sim1D.h | 54 +- Cantera/src/oneD/StFlow.cpp | 441 +++++++++++------ Cantera/src/oneD/StFlow.h | 24 +- Cantera/src/oneD/Surf1D.h | 10 +- Cantera/src/oneD/newton_utils.cpp | 55 ++- Cantera/src/oneD/refine.cpp | 97 ++-- Cantera/src/oneD/refine.h | 14 +- Cantera/src/stringUtils.h | 10 + Cantera/src/{ => transport}/L_matrix.h | 0 Makefile.in | 2 +- apps/MixMaster/KineticsFrame.py | 2 + config/configure | 390 +++++++-------- config/configure.in | 29 +- configure | 9 +- examples/cxx/Makefile.in | 2 +- include/onedim.h | 1 + test_problems/silane_equil/Makefile.in | 2 +- tools/Makefile.in | 3 + tools/src/Makefile.in | 2 +- tools/templates/cxx/demo.mak.in | 2 +- tools/templates/f77/demo.mak.in | 2 +- tools/testtools/Makefile.in | 2 +- 108 files changed, 3305 insertions(+), 1142 deletions(-) create mode 100644 Cantera/clib/src/ctonedim.cpp create mode 100644 Cantera/clib/src/ctonedim.h create mode 100755 Cantera/matlab/cantera/1D/@Domain1D/Domain1D.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/componentIndex.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/componentName.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/disableEnergy.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/domainIndex.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/domainType.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/domain_hndl.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/enableEnergy.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/gridPoints.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/isFlow.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/isInlet.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/isSurface.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/massFlux.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/massFraction.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/nComponents.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/nPoints.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/private/domain_methods.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/set.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setBounds.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setFixedTempProfile.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setID.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setMdot.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setMoleFractions.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setPressure.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setProfile.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setTemperature.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setTolerances.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/setupGrid.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/temperature.m create mode 100644 Cantera/matlab/cantera/1D/@Domain1D/z.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/Stack.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/display.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/domainIndex.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/grid.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/plotSolution.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/private/stack_methods.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/resid.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/restore.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/save.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/saveSoln.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/setFlatProfile.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/setProfile.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/setRefineCriteria.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/setValue.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/solution.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/solve.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/subsref.m create mode 100644 Cantera/matlab/cantera/1D/@Stack/writeStats.m create mode 100644 Cantera/matlab/cantera/1D/AxiStagnFlow.m create mode 100644 Cantera/matlab/cantera/1D/AxisymmetricFlow.m create mode 100644 Cantera/matlab/cantera/1D/Inlet.m create mode 100644 Cantera/matlab/cantera/1D/Outlet.m create mode 100644 Cantera/matlab/cantera/1D/Surface.m create mode 100644 Cantera/matlab/cantera/1D/SymmPlane.m create mode 100644 Cantera/matlab/cantera/private/onedimmethods.cpp create mode 100644 Cantera/python/examples/mix1.py create mode 100644 Cantera/python/examples/mix2.py delete mode 100755 Cantera/src/Flow1D.h delete mode 100755 Cantera/src/Jac2.h delete mode 100755 Cantera/src/Resid1D.h create mode 100644 Cantera/src/oneD/Domain1D.h rename Cantera/src/{ => transport}/L_matrix.h (100%) diff --git a/Cantera/clib/src/Makefile.in b/Cantera/clib/src/Makefile.in index 77435b4b4..b59e2904c 100755 --- a/Cantera/clib/src/Makefile.in +++ b/Cantera/clib/src/Makefile.in @@ -14,7 +14,7 @@ SUFFIXES= .cpp .d .o CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) -OBJS = ct.o Storage.o ctstagn.o ctsurf.o ctrpath.o ctbdry.o ctreactor.o ctfunc.o ctxml.o +OBJS = ct.o Storage.o ctstagn.o ctsurf.o ctrpath.o ctbdry.o ctreactor.o ctfunc.o ctxml.o ctonedim.o DEPENDS = $(OBJS:.o=.d) diff --git a/Cantera/clib/src/ctbdry.cpp b/Cantera/clib/src/ctbdry.cpp index cb5e2a0a9..669aa38a4 100755 --- a/Cantera/clib/src/ctbdry.cpp +++ b/Cantera/clib/src/ctbdry.cpp @@ -62,16 +62,26 @@ extern "C" { } int DLL_EXPORT bndry_settemperature(int i, double t) { - _bndry(i)->setTemperature(t); + try { + _bndry(i)->setTemperature(t); + } + catch (CanteraError) {return -1;} return 0; } double DLL_EXPORT bndry_spreadrate(int i) { - return ((Inlet1D*)_bndry(i))->spreadRate(); + try { + return ((Inlet1D*)_bndry(i))->spreadRate(); + } + catch (CanteraError) {return -1;} + return 0; } int DLL_EXPORT bndry_setSpreadRate(int i, double v) { - ((Inlet1D*)_bndry(i))->setSpreadRate(v); + try { + ((Inlet1D*)_bndry(i))->setSpreadRate(v); + } + catch (CanteraError) {return -1;} return 0; } @@ -90,12 +100,18 @@ extern "C" { } int DLL_EXPORT bndry_setxin(int i, double* xin) { - _bndry(i)->setMoleFractions(xin); + try { + _bndry(i)->setMoleFractions(xin); + } + catch (CanteraError) {return -1;} return 0; } int DLL_EXPORT bndry_setxinbyname(int i, char* xin) { - _bndry(i)->setMoleFractions(string(xin)); + try { + _bndry(i)->setMoleFractions(string(xin)); + } + catch (CanteraError) {return -1;} return 0; } } diff --git a/Cantera/clib/src/ctonedim.cpp b/Cantera/clib/src/ctonedim.cpp new file mode 100644 index 000000000..387eef95e --- /dev/null +++ b/Cantera/clib/src/ctonedim.cpp @@ -0,0 +1,461 @@ + +// Cantera includes +#include "oneD/Sim1D.h" +#include "oneD/StFlow.h" +#include "oneD/Inlet1D.h" +#include "DenseMatrix.h" +#include "Cabinet.h" +#include "Storage.h" + +// Build as a DLL under Windows +#ifdef WIN32 +#define DLL_EXPORT __declspec(dllexport) +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#else +#define DLL_EXPORT +#endif + +// Values returned for error conditions +#define ERR -999 +#define DERR -999.999 + + +Cabinet* Cabinet::__storage = 0; +Cabinet* Cabinet::__storage = 0; + + +inline Sim1D* _sim1D(int i) { + return Cabinet::cabinet()->item(i); +} + +inline Domain1D* _domain(int i) { + return Cabinet::cabinet()->item(i); +} + +static StFlow* _stflow(int i) { + Domain1D* d = _domain(i); + if (d->domainType() == cFlowType) return (StFlow*)d; + else + throw CanteraError("_stflow","wrong domain type"); +} + +static Bdry1D* _bdry(int i) { + Domain1D* d = _domain(i); + if (d->isConnector()) return (Bdry1D*)d; + else + throw CanteraError("_bdry","wrong domain type: " + +int2str(d->domainType())); +} + +inline ThermoPhase* _phase(int n) { + return Storage::__storage->__thtable[n]; +} + +inline Kinetics* _kinetics(int n) { + return Storage::__storage->__ktable[n]; +} + +inline ThermoPhase* _thermo(int n) { + return Storage::__storage->__thtable[n]; +} + +inline Transport* _transport(int n) { + return Storage::__storage->__trtable[n]; +} + + +extern "C" { + + int DLL_EXPORT domain_del(int i) { + Cabinet::cabinet()->del(i); + return 0; + } + + int DLL_EXPORT domain_type(int i) { + return _domain(i)->domainType(); + } + + int DLL_EXPORT domain_index(int i) { + return _domain(i)->domainIndex(); + } + + int DLL_EXPORT domain_nComponents(int i) { + return _domain(i)->nComponents(); + } + + int DLL_EXPORT domain_nPoints(int i) { + return _domain(i)->nPoints(); + } + + + int DLL_EXPORT domain_componentName(int i, int n, int sz, char* buf) { + try { + string nm = _domain(i)->componentName(n); + int lout = min(sz, nm.size()); + copy(nm.c_str(), nm.c_str() + lout, buf); + buf[lout] = '\0'; + return nm.size(); + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT domain_componentIndex(int i, char* name) { + try { + int n = _domain(i)->componentIndex(string(name)); + return n; + } + catch (CanteraError) { return -1; } + } + + double DLL_EXPORT domain_grid(int i, int n) { + try { + return _domain(i)->grid(n); + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT domain_setBounds(int i, int nl, double* lower, + int nu, double* upper) { + try { + _domain(i)->setBounds(nl, lower, nu, upper); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT domain_setTolerances(int i, int nr, double* rtol, + int na, double* atol, int itime) { + try { + _domain(i)->setTolerances(nr, rtol, na, atol, itime); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT domain_setupGrid(int i, int npts, double* grid) { + try { + _domain(i)->setupGrid(npts, grid); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT domain_setID(int i, char* id) { + try { + string s = string(id); + _domain(i)->setID(s); + return 0; + } + catch (CanteraError) { return -1; } + } + + + int DLL_EXPORT domain_setDesc(int i, char* desc) { + try { + string s = string(desc); + _domain(i)->setDesc(s); + return 0; + } + catch (CanteraError) { return -1; } + } + + + //------------------ inlet domains ------------------------------ + + int DLL_EXPORT inlet_new() { + try { + Inlet1D* i = new Inlet1D(); + return Cabinet::cabinet()->add(i); + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT surf_new() { + try { + Surf1D* i = new Surf1D(); + return Cabinet::cabinet()->add(i); + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT symm_new() { + try { + Symm1D* i = new Symm1D(); + return Cabinet::cabinet()->add(i); + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT outlet_new() { + try { + Outlet1D* i = new Outlet1D(); + return Cabinet::cabinet()->add(i); + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT bdry_setMdot(int i, double mdot) { + try { + _bdry(i)->setMdot(mdot); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT bdry_setTemperature(int i, double t) { + try { + _bdry(i)->setTemperature(t); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT bdry_setMoleFractions(int i, char* x) { + try { + _bdry(i)->setMoleFractions(string(x)); + return 0; + } + catch (CanteraError) { return -1; } + } + + double DLL_EXPORT bdry_temperature(int i) { + try { + return _bdry(i)->temperature(); + } + catch (CanteraError) { return -1; } + } + + double DLL_EXPORT bdry_massFraction(int i, int k) { + try { + return _bdry(i)->massFraction(k); + } + catch (CanteraError) { return -1; } + } + + double DLL_EXPORT bdry_mdot(int i) { + try { + return _bdry(i)->mdot(); + } + catch (CanteraError) { return -1; } + } + + + //------------------ stagnation flow domains -------------------- + + int DLL_EXPORT stflow_new(int iph, int ikin, int itr) { + try { + IdealGasPhase* ph = (IdealGasPhase*)_thermo(iph); + AxiStagnFlow* x = new AxiStagnFlow(ph, ph->nSpecies(), 2); + x->setKinetics(*_kinetics(ikin)); + x->setTransport(*_transport(ikin)); + + return Cabinet::cabinet()->add(x); + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT stflow_setPressure(int i, double p) { + try { + _stflow(i)->setPressure(p); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT stflow_setFixedTempProfile(int i, int n, double* pos, + double* temp) { + try { + int j; + vector_fp vpos(n), vtemp(n); + for (j = 0; j < n; j++) { + vpos[j] = pos[j]; + vtemp[j] = temp[j]; + } + _stflow(i)->setFixedTempProfile(vpos, vtemp); + return 0; + } + catch (CanteraError) { return -1; } + } + + + int DLL_EXPORT stflow_solveSpeciesEqs(int i, int flag) { + try { + if (flag > 0) + _stflow(i)->solveSpecies(-1); + else + _stflow(i)->fixSpecies(-1); + return 0; + } + catch (CanteraError) { return -1; } + } + + + int DLL_EXPORT stflow_solveEnergyEqn(int i, int flag) { + try { + if (flag > 0) + _stflow(i)->solveEnergyEqn(-1); + else + _stflow(i)->fixTemperature(-1); + return 0; + } + catch (CanteraError) { return -1; } + } + + + //------------------- Sim1D -------------------------------------- + + int DLL_EXPORT sim1D_new(int nd, int* domains) { + vector d; + try { + for (int n = 0; n < nd; n++) { + d.push_back(_domain(domains[n])); + } + Sim1D* s = new Sim1D(d); + return Cabinet::cabinet()->add(s); + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_del(int i) { + Cabinet::cabinet()->del(i); + return 0; + } + + int DLL_EXPORT sim1D_setValue(int i, int dom, int comp, + int localPoint, double value) { + try { + _sim1D(i)->setValue(dom, comp, localPoint, value); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_setProfile(int i, int dom, int comp, + int np, double* pos, double* v) { + try { + vector_fp vv, pv; + for (int n = 0; n < np; n++) { + vv.push_back(v[n]); + pv.push_back(pos[n]); + } + _sim1D(i)->setProfile(dom, comp, pv, vv); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_setFlatProfile(int i, int dom, int comp, double v) { + try { + _sim1D(i)->setFlatProfile(dom, comp, v); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_showSolution(int i, char* fname) { + string fn = string(fname); + if (fn == "-") + _sim1D(i)->showSolution(); + else { + ofstream fout(fname); + _sim1D(i)->showSolution(fout); + fout.close(); + } + return 0; + } + + int DLL_EXPORT sim1D_setTimeStep(int i, double stepsize, int ns, int* nsteps) { + try { + _sim1D(i)->setTimeStep(stepsize, ns, nsteps); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_solve(int i, int loglevel, int refine_grid) { + try { + bool r = (refine_grid == 0 ? false : true); + _sim1D(i)->solve(loglevel, r); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_refine(int i, int loglevel) { + try { + _sim1D(i)->refine(loglevel); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_setRefineCriteria(int i, int dom, double ratio, + double slope, double curve) { + try { + _sim1D(i)->setRefineCriteria(dom, ratio, slope, curve); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_save(int i, char* fname, char* id, + char* desc) { + try { + string sname = string(fname); + string sid = string(id); + string sdesc = string(desc); + _sim1D(i)->save(sname, sid, sdesc); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_restore(int i, char* fname, char* id) { + try { + string sname = string(fname); + string sid = string(id); + _sim1D(i)->restore(sname, sid); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_writeStats(int i) { + try { + _sim1D(i)->writeStats(); + return 0; + } + catch (CanteraError) { return -1; } + } + + int DLL_EXPORT sim1D_domainIndex(int i, char* name) { + try { + return _sim1D(i)->domainIndex(string(name)); + } + catch (CanteraError) { return -1; } + } + + double DLL_EXPORT sim1D_value(int i, int idom, int icomp, int localPoint) { + try { + return _sim1D(i)->value(idom, icomp, localPoint); + } + catch (CanteraError) { return -1.0; } + } + + double DLL_EXPORT sim1D_workValue(int i, int idom, int icomp, int localPoint) { + try { + return _sim1D(i)->workValue(idom, icomp, localPoint); + } + catch (CanteraError) { return -1.0; } + } + + int DLL_EXPORT sim1D_eval(int i, double rdt, int count) { + try { + _sim1D(i)->eval(rdt, count); + return 0; + } + catch (CanteraError) { return -1; } + } + +} diff --git a/Cantera/clib/src/ctonedim.h b/Cantera/clib/src/ctonedim.h new file mode 100644 index 000000000..1eab0aa77 --- /dev/null +++ b/Cantera/clib/src/ctonedim.h @@ -0,0 +1,67 @@ +#ifndef CTC_ONEDIM_H +#define CTC_ONEDIM_H + +#include "clib_defs.h" + +extern "C" { + + int DLL_IMPORT domain_del(int i); + int DLL_IMPORT domain_type(int i); + int DLL_IMPORT domain_index(int i); + int DLL_IMPORT domain_nComponents(int i); + int DLL_IMPORT domain_nPoints(int i); + int DLL_IMPORT domain_componentName(int i, int n, int sz, char* buf); + int DLL_IMPORT domain_componentIndex(int i, char* name); + int DLL_IMPORT domain_setBounds(int i, int nl, double* lower, + int nu, double* upper); + int DLL_IMPORT domain_setTolerances(int i, int nr, double* rtol, + int na, double* atol, int itime); + int DLL_IMPORT domain_setupGrid(int i, int npts, double* grid); + int DLL_IMPORT domain_setID(int, char* id); + int DLL_IMPORT domain_setDesc(int, char* desc); + double DLL_IMPORT domain_grid(int i, int n); + + int DLL_IMPORT bdry_setMdot(int i, double mdot); + int DLL_IMPORT bdry_setTemperature(int i, double t); + int DLL_IMPORT bdry_setMoleFractions(int i, char* x); + double DLL_IMPORT bdry_temperature(int i); + double DLL_IMPORT bdry_massFraction(int i, int k); + double DLL_IMPORT bdry_mdot(int i); + + int DLL_IMPORT inlet_new(); + int DLL_IMPORT outlet_new(); + int DLL_IMPORT symm_new(); + int DLL_IMPORT surf_new(); + + int DLL_IMPORT stflow_new(int iph, int ikin, int itr); + int DLL_IMPORT stflow_setPressure(int i, double p); + int DLL_IMPORT stflow_setFixedTempProfile(int i, int n, double* pos, + double* temp); + int DLL_IMPORT stflow_solveSpeciesEqs(int i, int flag); + int DLL_IMPORT stflow_solveEnergyEqn(int i, int flag); + + int DLL_IMPORT sim1D_new(int nd, int* domains); + int DLL_IMPORT sim1D_del(int i); + int DLL_IMPORT sim1D_setValue(int i, int dom, int comp, int localPoint, double value); + int DLL_IMPORT sim1D_setProfile(int i, int dom, int comp, + int np, double* pos, double* v); + int DLL_IMPORT sim1D_setFlatProfile(int i, int dom, int comp, double v); + int DLL_IMPORT sim1D_showSolution(int i, char* fname); + int DLL_IMPORT sim1D_setTimeStep(int i, double stepsize, int ns, int* nsteps); + int DLL_IMPORT sim1D_solve(int i, int loglevel, int refine_grid); + int DLL_IMPORT sim1D_refine(int i, int loglevel); + int DLL_IMPORT sim1D_setRefineCriteria(int i, int dom, double ratio, + double slope, double curve); + int DLL_IMPORT sim1D_save(int i, char* fname, char* id, + char* desc); + int DLL_IMPORT sim1D_restore(int i, char* fname, char* id); + int DLL_IMPORT sim1D_writeStats(int i); + int DLL_IMPORT sim1D_domainIndex(int i, char* name); + double DLL_IMPORT sim1D_value(int i, int idom, int icomp, int localPoint); + double DLL_IMPORT sim1D_workValue(int i, int idom, + int icomp, int localPoint); + int DLL_IMPORT sim1D_eval(int i, double rdt, int count); +} + + +#endif diff --git a/Cantera/clib/src/ctstagn.cpp b/Cantera/clib/src/ctstagn.cpp index 68bc7d0d0..c2e5d9734 100755 --- a/Cantera/clib/src/ctstagn.cpp +++ b/Cantera/clib/src/ctstagn.cpp @@ -330,7 +330,7 @@ extern "C" { int DLL_EXPORT onedim_new(int nd, int* domains, int* types) { int i; - vector doms; + vector doms; for (i = 0; i < nd; i++) { switch (types[i]) { case 0: diff --git a/Cantera/cxx/Makefile.in b/Cantera/cxx/Makefile.in index 9e413dd4b..e7cb89e5a 100644 --- a/Cantera/cxx/Makefile.in +++ b/Cantera/cxx/Makefile.in @@ -14,7 +14,7 @@ SUFFIXES= .cpp .d .o CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) -OBJS = cxxutils.o +OBJS = cxxutils.o writelog.o DEPENDS = $(OBJS:.o=.d) @@ -34,7 +34,7 @@ LCXX_FLAGS = -L$(CANTERA_LIBDIR) @CXXFLAGS@ $(CXX) -c $< $(CXX_INCLUDES) $(CXX_FLAGS) LIB_NAME=libctcxx -CXXLIB=$(LIB_NAME).a +CXXLIB=@buildlib@/$(LIB_NAME).a all: $(OBJS) @ARCHIVE@ $(CXXLIB) $(OBJS) diff --git a/Cantera/matlab/Makefile.in b/Cantera/matlab/Makefile.in index d9e069993..aceb3ab9d 100644 --- a/Cantera/matlab/Makefile.in +++ b/Cantera/matlab/Makefile.in @@ -19,7 +19,8 @@ SRCS = cantera/private/ctmethods.cpp \ cantera/private/transportmethods.cpp \ cantera/private/reactormethods.cpp \ cantera/private/wallmethods.cpp \ - cantera/private/flowdevicemethods.cpp + cantera/private/flowdevicemethods.cpp \ + cantera/private/onedimmethods.cpp CANTERA_LIBDIR=@buildlib@ LIB_DEPS = $(CANTERA_LIBDIR)/libcantera.a $(CANTERA_LIBDIR)/libzeroD.a \ @@ -44,6 +45,11 @@ install: @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Reactor/private @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Wall/private @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@FlowDevice/private + @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/1D + @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/1D/@Domain1D + @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/1D/@Domain1D/private + @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/1D/@Stack + @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/1D/@Stack/private cd cantera; @INSTALL@ *.m *.@mex_ext@ @prefix@/matlab/toolbox/cantera/cantera cd cantera/private; @INSTALL@ *.m \ @prefix@/matlab/toolbox/cantera/cantera/private @@ -76,6 +82,15 @@ install: @prefix@/matlab/toolbox/cantera/cantera/@FlowDevice cd cantera/@FlowDevice/private; @INSTALL@ *.m \ @prefix@/matlab/toolbox/cantera/cantera/@FlowDevice/private + cd cantera/1D/@Domain1D; @INSTALL@ *.m \ + @prefix@/matlab/toolbox/cantera/cantera/1D/@Domain1D + cd cantera/1D/@Domain1D/private; @INSTALL@ *.m \ + @prefix@/matlab/toolbox/cantera/cantera/1D/@Domain1D/private + cd cantera/1D/@Stack; @INSTALL@ *.m \ + @prefix@/matlab/toolbox/cantera/cantera/1D/@Stack + cd cantera/1D/@Stack/private; @INSTALL@ *.m \ + @prefix@/matlab/toolbox/cantera/cantera/1D/@Stack/private + clean: echo '-' diff --git a/Cantera/matlab/cantera/1D/@Domain1D/Domain1D.m b/Cantera/matlab/cantera/1D/@Domain1D/Domain1D.m new file mode 100755 index 000000000..d9e54ca55 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/Domain1D.m @@ -0,0 +1,21 @@ +function d = Domain1D(a, b, c, d, e) +% DOMAIN1D - Create a new one-dimensional domain. +% +d.dom_id = -1; + +if nargin == 1 + d.dom_id = domain_methods(0, a); +elseif nargin == 2 + if a == 1 + if isa(b,'Solution') + d.dom_id = domain_methods(0, 1, thermo_hndl(b), kinetics_hndl(b), ... + trans_hndl(b)); + else + error('Wrong argument type. Expecting instance of class Solution.') + end + end +end +if d.dom_id < 0 + error(geterr); +end +d = class(d, 'Domain1D'); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/componentIndex.m b/Cantera/matlab/cantera/1D/@Domain1D/componentIndex.m new file mode 100644 index 000000000..cfc223ed6 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/componentIndex.m @@ -0,0 +1,9 @@ +function n = componentIndex(d, name) +% COMPONENTINDEX - +% +if isa(name,'double') + n = name; +else + n = domain_methods(d.dom_id, 18, name); +end + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/componentName.m b/Cantera/matlab/cantera/1D/@Domain1D/componentName.m new file mode 100644 index 000000000..91eb86d6d --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/componentName.m @@ -0,0 +1,7 @@ +function s = componentName(d, n) +% COMPONENTNAME - Name of component n. +% +m = length(n); +for i = 1:m + s{i} = domain_methods(d.dom_id, 40, n(i)); +end diff --git a/Cantera/matlab/cantera/1D/@Domain1D/disableEnergy.m b/Cantera/matlab/cantera/1D/@Domain1D/disableEnergy.m new file mode 100644 index 000000000..f60bede4d --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/disableEnergy.m @@ -0,0 +1,4 @@ +function d = disableEnergy(d) +% ENABLEENERGY - enable the energy equation +% +domain_methods(d.dom_id, 66, 0); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/domainIndex.m b/Cantera/matlab/cantera/1D/@Domain1D/domainIndex.m new file mode 100644 index 000000000..687c29008 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/domainIndex.m @@ -0,0 +1,8 @@ +function i = domainIndex(d) +% DOMAININDEX - domain index. +% +% This function returns an integer flag denoting the location +% of the domain, beginning with 1 at the left. +% +i = domain_methods(d.dom_id, 13) + 1; + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/domainType.m b/Cantera/matlab/cantera/1D/@Domain1D/domainType.m new file mode 100644 index 000000000..6902ec7af --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/domainType.m @@ -0,0 +1,7 @@ +function i = domainType(d) +% DOMAINTYPE - Type of domain. +% +% This function returns an integer flag denoting the domain +% type. +i = domain_methods(d.dom_id, 12); + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/domain_hndl.m b/Cantera/matlab/cantera/1D/@Domain1D/domain_hndl.m new file mode 100644 index 000000000..17ca3d452 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/domain_hndl.m @@ -0,0 +1,5 @@ +function n = domain_hndl(d) +% DOMAIN_HNDL - Integer used to access kernel object. +% +n = d.dom_id; + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/enableEnergy.m b/Cantera/matlab/cantera/1D/@Domain1D/enableEnergy.m new file mode 100644 index 000000000..d8e93f6b2 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/enableEnergy.m @@ -0,0 +1,7 @@ +function d = enableEnergy(d) +% ENABLEENERGY - enable the energy equation +% +disp(' '); +disp('Enabling the energy equation...'); + +domain_methods(d.dom_id, 66, 1); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/gridPoints.m b/Cantera/matlab/cantera/1D/@Domain1D/gridPoints.m new file mode 100644 index 000000000..83385236d --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/gridPoints.m @@ -0,0 +1,15 @@ +function zz = gridPoints(d, n) +% GRID - +% +if nargin == 1 + for i = 1:nPoints(d) + zz(i) = domain_methods(d.dom_id, 19, i); + end +else + m = length(n); + for i = 1:m + zz(i) = domain_methods(d.dom_id, 19, n(i)); + end +end + + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/isFlow.m b/Cantera/matlab/cantera/1D/@Domain1D/isFlow.m new file mode 100644 index 000000000..fc4c6fd29 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/isFlow.m @@ -0,0 +1,10 @@ +function a = isFlow(d) +% ISFLOW - Returns 1 if the domain is a flow domain, and 0 otherwise. +% +t = domainType(d); +if t == 50 + a = 1; +else + a = 0; +end + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/isInlet.m b/Cantera/matlab/cantera/1D/@Domain1D/isInlet.m new file mode 100644 index 000000000..9daac3580 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/isInlet.m @@ -0,0 +1,10 @@ +function a = isInlet(d) +% ISINLET - Returns 1 if the domain is an inlet, and 0 otherwise. +% +t = domainType(d); +if t == 104 + a = 1; +else + a = 0; +end + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/isSurface.m b/Cantera/matlab/cantera/1D/@Domain1D/isSurface.m new file mode 100644 index 000000000..4c82a94d6 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/isSurface.m @@ -0,0 +1,10 @@ +function a = isSurface(d) +% ISSURFACE - Returns 1 if the domain is a surface, and 0 otherwise. +% +t = domainType(d); +if t == 102 + a = 1; +else + a = 0; +end + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/massFlux.m b/Cantera/matlab/cantera/1D/@Domain1D/massFlux.m new file mode 100644 index 000000000..c6e391cf4 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/massFlux.m @@ -0,0 +1,5 @@ +function mdot = massFlux(d) +% MASSFLUX - +% +mdot = domain_methods(d.dom_id, 17); + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/massFraction.m b/Cantera/matlab/cantera/1D/@Domain1D/massFraction.m new file mode 100644 index 000000000..f5864dd78 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/massFraction.m @@ -0,0 +1,17 @@ +function y = massFraction(d, k) +% MASSFRACTION - Mass fraction of species k. +% +% This method returns the mass fraction of species k, where +% k is the integer index of the species in the flow domain +% to which the boundary domain is attached. +% +if domainIndex(d) == 0 + error('no flow domain attached!') +end + +if isInlet(d) + y = domain_methods(d.dom_id,16,k-1); +else + error('not yet...'); +end + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/nComponents.m b/Cantera/matlab/cantera/1D/@Domain1D/nComponents.m new file mode 100644 index 000000000..b1ec5d366 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/nComponents.m @@ -0,0 +1,4 @@ +function n = nComponents(d) +% NCOMPONENTS - number of components +% +n = domain_methods(d.dom_id, 11); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/nPoints.m b/Cantera/matlab/cantera/1D/@Domain1D/nPoints.m new file mode 100644 index 000000000..d2c7da4e6 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/nPoints.m @@ -0,0 +1,5 @@ +function npts = nPoints(d) +% NPOINTS - Number of grid points. +% +npts = domain_methods(d.dom_id, 14); + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/private/domain_methods.m b/Cantera/matlab/cantera/1D/@Domain1D/private/domain_methods.m new file mode 100644 index 000000000..7ae071844 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/private/domain_methods.m @@ -0,0 +1,13 @@ +function v = domain_methods(n, job, a, b, c, d) +% +if nargin == 2 + v = ctmethods(90, n, job); +elseif nargin == 3 + v = ctmethods(90, n, job, a); +elseif nargin == 4 + v = ctmethods(90, n, job, a, b); +elseif nargin == 5 + v = ctmethods(90, n, job, a, b, c); +elseif nargin == 6 + v = ctmethods(90, n, job, a, b, c, d); +end diff --git a/Cantera/matlab/cantera/1D/@Domain1D/set.m b/Cantera/matlab/cantera/1D/@Domain1D/set.m new file mode 100644 index 000000000..1e51bfb6c --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/set.m @@ -0,0 +1,88 @@ +function a = set(a,varargin) +% SET - Set properties. +% +% The properties that may be set are +% +% Either the full property name or the symbol may be +% specified. For the extensive properties (V,H,U,S), the values +% must be given per unit mass. H, U, and S must be set in +% conjunction with pressure (for H,S) or volume (for U,S). Either +% (specific) volume or density may be specified. Mole and mass +% fractions must be input as vectors (either row or column) with +% length equal to the number of species. +% +% Examples: +% +% set(gas,'Temperature',600.0); +% set(gas,'T',600.0); +% set(gas,'T',600.0,'P',2*oneatm,'Y',massfracs); +% set(gas,'H',0.5*enthalpy_mass(gas),'P',pressure(gas)); +% set(gas,'S',entropy_mass(gas),'P',0.5*pressure(gas)); +% set(gas,'X',ones(nSpecies(gas),1)); +% +% Alternatively, individual methods to set properties may be +% called (setTemperature, setMoleFractions, etc.) +% + +property_argin = varargin; + +while length(property_argin) >= 2, + prop = property_argin{1}; + val = property_argin{2}; + property_argin = property_argin(3:end); + switch prop + case 'Temperature' + setTemperature(a,val); + case 'T' + setTemperature(a,val); + case 'MassFractions' + setMassFractions(a,val); + case 'Y' + setMassFractions(a,val); + case 'mdot' + setMdot(a,val); + case 'MassFlux' + setMdot(a,val); + case 'P' + setPressure(a,val); + case 'Pressure' + setPressure(a,val); + case 'tol' + sz = size(val); + if sz == nComponents(a) + setTolerances(a, val(1,:), val(2,:)); + elseif length(val) == 2 + rt = val(1)*ones(1,nComponents(a)); + at = val(2)*ones(1,nComponents(a)); + setTolerances(a, rt, at); + else + error('wrong array size for error tolerances'); + end + case 'tol-time' + sz = size(val); + if sz == nComponents(a) + setTolerances(a, val(1,:), val(2,:)); + elseif length(val) == 2 + rt = val(1)*ones(1,nComponents(a)); + at = val(2)*ones(1,nComponents(a)); + setTolerances(a, rt, at, 'ts'); + else + error('wrong array size for error tolerances'); + end + case 'grid' + setupGrid(a, val); + case 'bounds' + setBounds(a, val(1,:), val(2,:)); + case 'X' + setMoleFractions(a, val); + case 'MoleFractions' + setMoleFractions(a, val); + case 'T_fixed' + setFixedTempProfile(a, val); + case 'ID' + setID(a, val); + otherwise + error(['unknown property ' char(prop)]); + end +end + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setBounds.m b/Cantera/matlab/cantera/1D/@Domain1D/setBounds.m new file mode 100644 index 000000000..451fe3d2a --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setBounds.m @@ -0,0 +1,4 @@ +function d = setBounds(d, lower, upper) +% SETBOUNDS - +% +domain_methods(d.dom_id, 51, lower, upper); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setFixedTempProfile.m b/Cantera/matlab/cantera/1D/@Domain1D/setFixedTempProfile.m new file mode 100644 index 000000000..ece525165 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setFixedTempProfile.m @@ -0,0 +1,12 @@ +function d = setFixedTempProfile(d, profile) +% SETFIXEDTEMPPROFILE - +% +sz = size(profile); +if sz(1) == 2 + domain_methods(d.dom_id, 64, profile(1,:), profile(2,:)); +elseif sz(2) == 2 + domain_methods(d.dom_id, 64, profile(:,1), profile(:,2)); +else + error('wrong temperature profile array shape'); +end + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setID.m b/Cantera/matlab/cantera/1D/@Domain1D/setID.m new file mode 100644 index 000000000..4a325d0f0 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setID.m @@ -0,0 +1,4 @@ +function d = setID(d, id) +% SETID - Set the ID tag for the domain. +% +domain_methods(d.dom_id, 54, id); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setMdot.m b/Cantera/matlab/cantera/1D/@Domain1D/setMdot.m new file mode 100644 index 000000000..98647488a --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setMdot.m @@ -0,0 +1,4 @@ +function d = setMdot(d, mdot) +% SETMDOT - +% +domain_methods(d.dom_id, 60, mdot); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setMoleFractions.m b/Cantera/matlab/cantera/1D/@Domain1D/setMoleFractions.m new file mode 100644 index 000000000..ac9d7b5b7 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setMoleFractions.m @@ -0,0 +1,4 @@ +function d = setMoleFractions(d, x) +% SETMOLEFRACTIONS - +% +domain_methods(d.dom_id, 62, x); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setPressure.m b/Cantera/matlab/cantera/1D/@Domain1D/setPressure.m new file mode 100644 index 000000000..02df91ed4 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setPressure.m @@ -0,0 +1,5 @@ +function d = setPressure(d, p) +% SETPRESSURE - +% +domain_methods(d.dom_id, 63, p); + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setProfile.m b/Cantera/matlab/cantera/1D/@Domain1D/setProfile.m new file mode 100644 index 000000000..9abdfbec3 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setProfile.m @@ -0,0 +1,9 @@ +function d = setProfile(d, n, p) +% SETPROFILE - +% +if d.stack == 0 + error('install domain in stack before calling setProfile.'); +end + +setProfile(d.stack,domainIndex(d),n,p); + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setTemperature.m b/Cantera/matlab/cantera/1D/@Domain1D/setTemperature.m new file mode 100644 index 000000000..1a23e7058 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setTemperature.m @@ -0,0 +1,4 @@ +function d = setTemperature(d, t) +% SETTEMPERATURE - Set the temperature [K]. +% +domain_methods(d.dom_id, 61, t); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setTolerances.m b/Cantera/matlab/cantera/1D/@Domain1D/setTolerances.m new file mode 100644 index 000000000..10a5510f5 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setTolerances.m @@ -0,0 +1,17 @@ +function d = setTolerances(d, rtol, atol, typ) +% SETTOLERANCES - +% +ityp = 0; +if nargin == 4 + switch typ + case 'ts' + itype = -1; + case 'time' + itype = -1; + case 'ss' + itype = 1; + case 'steady' + itype = 1; + end +end +domain_methods(d.dom_id, 52, rtol, atol, ityp); diff --git a/Cantera/matlab/cantera/1D/@Domain1D/setupGrid.m b/Cantera/matlab/cantera/1D/@Domain1D/setupGrid.m new file mode 100644 index 000000000..0922f24dc --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/setupGrid.m @@ -0,0 +1,5 @@ +function d = setupGrid(d, grid) +% SETUPGRID - +% +domain_methods(d.dom_id, 53, grid); + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/temperature.m b/Cantera/matlab/cantera/1D/@Domain1D/temperature.m new file mode 100644 index 000000000..d91e5e332 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/temperature.m @@ -0,0 +1,5 @@ +function t = temperature(d) +% TEMPERATURE - Temperature [K]. +% +t = domain_methods(d.dom_id, 15); + diff --git a/Cantera/matlab/cantera/1D/@Domain1D/z.m b/Cantera/matlab/cantera/1D/@Domain1D/z.m new file mode 100644 index 000000000..3da29ffe2 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Domain1D/z.m @@ -0,0 +1,15 @@ +function zz = z(d, n) +% GRID - +% +if nargin == 1 + for i = 1:nPoints(d) + zz(i) = domain_methods(d.dom_id, 19, i); + end +else + m = length(n); + for i = 1:m + zz(i) = domain_methods(d.dom_id, 19, n(i)); + end +end + + diff --git a/Cantera/matlab/cantera/1D/@Stack/Stack.m b/Cantera/matlab/cantera/1D/@Stack/Stack.m new file mode 100644 index 000000000..e19c49134 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/Stack.m @@ -0,0 +1,25 @@ +function s = Stack(domains) +% +% STACK - A one-dimensional 'stack' of domains. +% +% A stack object is a container for one-dimensional domains, +% which are instances of class Domain1D. The domains are of two +% types - extended domains, and connector domains. +% +s.stack_id = -1; +s.domains = domains; +if nargin == 1 + nd = length(domains); + for n=1:nd + ids(n) = domain_hndl(domains(n)); + end + s.stack_id = stack_methods(0, 8, nd, ids); +else + help(Stack); + error('wrong number of parameters'); +end +if s.stack_id < 0 + error(geterr); +end +s = class(s, 'Stack'); + diff --git a/Cantera/matlab/cantera/1D/@Stack/display.m b/Cantera/matlab/cantera/1D/@Stack/display.m new file mode 100644 index 000000000..4a93e77b2 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/display.m @@ -0,0 +1,8 @@ +function display(s, fname) +% DISPLAY - show all domains. +% +if nargin == 1 + fname = '-'; +end +stack_methods(s.stack_id, 103, fname); + diff --git a/Cantera/matlab/cantera/1D/@Stack/domainIndex.m b/Cantera/matlab/cantera/1D/@Stack/domainIndex.m new file mode 100644 index 000000000..7beff60e0 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/domainIndex.m @@ -0,0 +1,7 @@ +function n = domainIndex(d, name) +% DOMAININDEX - Index of the domain with a specified name. +if isa(name,'double') + n = name +else + n = stack_methods(d.stack_id, 109, name); +end diff --git a/Cantera/matlab/cantera/1D/@Stack/grid.m b/Cantera/matlab/cantera/1D/@Stack/grid.m new file mode 100644 index 000000000..0ffde0e29 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/grid.m @@ -0,0 +1,6 @@ +function z = grid(s, d) +% GRID - the grid in one domain. +% +n = domainIndex(s,d); +d = s.domains(n); +z = gridPoints(d); diff --git a/Cantera/matlab/cantera/1D/@Stack/plotSolution.m b/Cantera/matlab/cantera/1D/@Stack/plotSolution.m new file mode 100644 index 000000000..856278c3f --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/plotSolution.m @@ -0,0 +1,13 @@ +function plotSolution(s, domain, component) +% PLOTSOLUTION - plot a specified solution component +% +% plotSolution(s, 'flow', 'T') plots component 'T' in domain 'flow' +% +n = domainIndex(s,domain); +d = s.domains(n); +z = gridPoints(d); +x = solution(s, domain, component); +plot(z, x); +xlabel('z (m)'); +ylabel(component); + diff --git a/Cantera/matlab/cantera/1D/@Stack/private/stack_methods.m b/Cantera/matlab/cantera/1D/@Stack/private/stack_methods.m new file mode 100644 index 000000000..caf710b04 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/private/stack_methods.m @@ -0,0 +1,22 @@ +function v = stack_methods(n, job, a, b, c, d, e) +% STACK_METHODS - converter function for methods of class Stack +% +% All Cantera functions and methods are handled by the single MEX +% file 'ctmethods.' This function is provided only for convenience, +% and simply calls ctmethods with a flag associated with this class +% as the first parameter, followed by the input arguments. +if nargin == 2 + v = ctmethods(90, n, job); +elseif nargin == 3 + v = ctmethods(90, n, job, a); +elseif nargin == 4 + v = ctmethods(90, n, job, a, b); +elseif nargin == 5 + v = ctmethods(90, n, job, a, b, c); +elseif nargin == 6 + v = ctmethods(90, n, job, a, b, c, d); +elseif nargin == 7 + v = ctmethods(90, n, job, a, b, c, d, e); +else + error('wrong number of arguments'); +end diff --git a/Cantera/matlab/cantera/1D/@Stack/resid.m b/Cantera/matlab/cantera/1D/@Stack/resid.m new file mode 100644 index 000000000..b6764e42d --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/resid.m @@ -0,0 +1,16 @@ +function r = resid(s, domain, rdt, count) +if nargin == 2 + rdt = 0.0; + count = 0; +end + +idom = domainIndex(s, domain); +d = s.domains(idom); + +r = zeros(nComponents(d), nPoints(d)); +stack_methods(s.stack_id, 113, rdt, count); +for m = 1:nComponents(d) + for n = 1:nPoints(d) + r(m,n) = stack_methods(s.stack_id, 31, idom, m, n); + end +end diff --git a/Cantera/matlab/cantera/1D/@Stack/restore.m b/Cantera/matlab/cantera/1D/@Stack/restore.m new file mode 100644 index 000000000..90e6718fa --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/restore.m @@ -0,0 +1,6 @@ +function restore(s, fname, id) +% RESTORE - Restore a previously-saved solution. +% +% This method can be used to provide an initial guess for the +% solution. +stack_methods(s.stack_id, 111, fname, id); diff --git a/Cantera/matlab/cantera/1D/@Stack/save.m b/Cantera/matlab/cantera/1D/@Stack/save.m new file mode 100644 index 000000000..7dfd2d5f4 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/save.m @@ -0,0 +1,10 @@ +function saveSoln(s, fname, id, desc) +% SAVE - +% +if nargin == 2 + id = 'solution'; + desc = '-'; +elseif nargin == 3 + desc = '-'; +end +stack_methods(s.stack_id, 107, fname, id, desc); diff --git a/Cantera/matlab/cantera/1D/@Stack/saveSoln.m b/Cantera/matlab/cantera/1D/@Stack/saveSoln.m new file mode 100644 index 000000000..edad98dec --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/saveSoln.m @@ -0,0 +1,14 @@ +function saveSoln(s, fname, id, desc) +% SAVE - Save solution. +% +if nargin == 1 + fname = 'soln.xml'; + id = 'solution'; + desc = '--'; +elseif nargin == 2 + id = 'solution'; + desc = '--'; +elseif nargin == 3 + desc = '--'; +end +stack_methods(s.stack_id, 107, fname, id, desc); diff --git a/Cantera/matlab/cantera/1D/@Stack/setFlatProfile.m b/Cantera/matlab/cantera/1D/@Stack/setFlatProfile.m new file mode 100644 index 000000000..8c58e46ad --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/setFlatProfile.m @@ -0,0 +1,5 @@ +function setFlatProfile(s, n, comp, v) +% SETFLATPROFILE - +% +stack_methods(s.stack_id, 102, n, comp, v); + diff --git a/Cantera/matlab/cantera/1D/@Stack/setProfile.m b/Cantera/matlab/cantera/1D/@Stack/setProfile.m new file mode 100644 index 000000000..5459ca12c --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/setProfile.m @@ -0,0 +1,52 @@ +function setProfile(s, name, comp, p) +% SETPROFILE - Specify a profile for one component. +% +% name -- domain name +% comp -- component number +% zr -- array of relative positions (0.0 to 1.0) +% v -- array of values +% +% The solution vector values for this component will be linearly +% interpolated from the discrete function defined by v vs. zr. +% Note that zr = 0.0 corresponds to the leftmost grid point in +% the specified domain, and zr = 1.0 corresponds to the rightmost +% grid point. This method can be called at any time, but is +% usually used to set the initial guess for the solution. +% +% Example: +% +% zr = [0 0.1 0.2 0.4 0.8 1]; +% v = [500 650 700 730 800 900]; +% setProfile(1, 2, zr, v); +% +if isa(name,'double') + n = name; +else + n = domainIndex(s, name); +end + +d = s.domains(n); + +if isa(comp,'double') | isa(comp,'cell') + c = comp; +elseif isa(comp,'char') + c = {comp}; +else + error('wrong type'); +end + +np = length(c); +sz = size(p); +if sz(1) == np + 1; + for j = 1:np + ic = componentIndex(d,c{j}); + stack_methods(s.stack_id, 101, n, ic, p(1,:), p(j+1,:)); + end +elseif sz(2) == np + 1; + ic = componentIndex(d,c{j}); + stack_methods(s.stack_id, 101, n, ic, p(:,1), p(:,j+1)); +else + error('wrong profile shape'); +end + + diff --git a/Cantera/matlab/cantera/1D/@Stack/setRefineCriteria.m b/Cantera/matlab/cantera/1D/@Stack/setRefineCriteria.m new file mode 100644 index 000000000..897c92119 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/setRefineCriteria.m @@ -0,0 +1,5 @@ +function d = setRefineCriteria(d, n, ratio, slope, curve) +% SETREFINECRITERIA - +% +stack_methods(d.stack_id, 106, n, ratio, slope, curve); + diff --git a/Cantera/matlab/cantera/1D/@Stack/setValue.m b/Cantera/matlab/cantera/1D/@Stack/setValue.m new file mode 100644 index 000000000..4c95d238a --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/setValue.m @@ -0,0 +1,20 @@ +function setValue(s, n, comp, localPoint, v) +% SETVALUE - Set the value of a single entry in the solution vector. +% +% n -- domain number +% comp -- component number +% localPoint -- local index of the grid point in the domain +% v -- value +% +% Example: +% +% setValue(s, 3, 5, 1, 5.6) +% +% This sets component 5 at the leftmost point (local point 1) in domain 3 +% to the value 5.6. Note that the local index always begins at 1 +% at the left of each domain, independent of the global index of +% the point, which depends on the location of this domain in the +% stack. +% +stack_methods(s.stack_id, 100, n, comp, localPoint, v); + diff --git a/Cantera/matlab/cantera/1D/@Stack/solution.m b/Cantera/matlab/cantera/1D/@Stack/solution.m new file mode 100644 index 000000000..a7f3a4129 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/solution.m @@ -0,0 +1,23 @@ +function x = solution(s, domain, component) +% SOLUTION - get a solution component in one domain. +% +% x = solution(s, 'flow', 'T') returns in vector x the values of +% solution component 'T' in domain 'flow'. +% +idom = domainIndex(s, domain); +d = s.domains(idom); + +if nargin == 3 + icomp = componentIndex(d, component); + for n = 1:nPoints(d) + x(n) = stack_methods(s.stack_id, 30, idom, icomp, n); + end +else + nc = nComponents(d); + np = nPoints(d); + for m = 1:nc + for n = 1:np + x(m,n) = stack_methods(s.stack_id, 30, idom, m, n); + end + end +end diff --git a/Cantera/matlab/cantera/1D/@Stack/solve.m b/Cantera/matlab/cantera/1D/@Stack/solve.m new file mode 100644 index 000000000..0481024e3 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/solve.m @@ -0,0 +1,5 @@ +function solve(s, loglevel, refine_grid) +% SOLVE - +% +stack_methods(s.stack_id, 104, loglevel, refine_grid); + diff --git a/Cantera/matlab/cantera/1D/@Stack/subsref.m b/Cantera/matlab/cantera/1D/@Stack/subsref.m new file mode 100644 index 000000000..edb54ce5a --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/subsref.m @@ -0,0 +1,11 @@ +function b = subsref(s,index) +% SUBSREF - +switch index.type + case '()' + b = s.domains(index.subs{:}); + case '.' + n = domainIndex(s, index.subs); + b = s.domains(n); + otherwise + error('syntax error'); +end diff --git a/Cantera/matlab/cantera/1D/@Stack/writeStats.m b/Cantera/matlab/cantera/1D/@Stack/writeStats.m new file mode 100644 index 000000000..683ddcef1 --- /dev/null +++ b/Cantera/matlab/cantera/1D/@Stack/writeStats.m @@ -0,0 +1,10 @@ +function writeStats(s) +% WRITESTATS - Print statistics for the current solution. +% +% writeStats(s) prints a summary of the number of function and +% Jacobian evaluations for each grid, and the CPU time spent on +% each one. +% +stack_methods(s.stack_id, 108); + + diff --git a/Cantera/matlab/cantera/1D/AxiStagnFlow.m b/Cantera/matlab/cantera/1D/AxiStagnFlow.m new file mode 100644 index 000000000..1f675538b --- /dev/null +++ b/Cantera/matlab/cantera/1D/AxiStagnFlow.m @@ -0,0 +1,6 @@ +function m = AxiStagnFlow(gas) +% AXISTAGNFLOW - Axisymmetric stagnation flow. +% +% Return a Domain1D instance representing an axisymmetric +% stagnation flow. +m = Domain1D(1, gas); diff --git a/Cantera/matlab/cantera/1D/AxisymmetricFlow.m b/Cantera/matlab/cantera/1D/AxisymmetricFlow.m new file mode 100644 index 000000000..8c423e9d6 --- /dev/null +++ b/Cantera/matlab/cantera/1D/AxisymmetricFlow.m @@ -0,0 +1,11 @@ +function m = AxisymmetricFlow(gas, id) +% AXISYMMETRICFLOW - Axisymmetric flow. +% +% Return a Domain1D instance representing an axisymmetric flow. +% +m = Domain1D(1, gas); +if nargin == 1 + setID(m,'flow'); +else + setID(m,id); +end \ No newline at end of file diff --git a/Cantera/matlab/cantera/1D/Inlet.m b/Cantera/matlab/cantera/1D/Inlet.m new file mode 100644 index 000000000..07ae84fc7 --- /dev/null +++ b/Cantera/matlab/cantera/1D/Inlet.m @@ -0,0 +1,13 @@ +function m = Inlet(id) +% INLET - Return a Domain1D instance representing an inlet. +% +% Note that an inlet can only be a terminal domain - it must be +% either the leftmost or rightmost domain in a stack. +m = Domain1D(2); +if nargin == 0 + setID(m,'inlet'); +else + setID(m,id); +end + + diff --git a/Cantera/matlab/cantera/1D/Outlet.m b/Cantera/matlab/cantera/1D/Outlet.m new file mode 100644 index 000000000..4fbe4a967 --- /dev/null +++ b/Cantera/matlab/cantera/1D/Outlet.m @@ -0,0 +1,11 @@ +function m = Outlet(id) +% OUTLET - Return a Domain1D instance representing an outlet. +% +m = Domain1D(5); +if nargin == 0 + setID(m,'outlet'); +else + setID(m,id); +end + + diff --git a/Cantera/matlab/cantera/1D/Surface.m b/Cantera/matlab/cantera/1D/Surface.m new file mode 100644 index 000000000..fc3493426 --- /dev/null +++ b/Cantera/matlab/cantera/1D/Surface.m @@ -0,0 +1,9 @@ +function m = Surface(id) +% SURFACE - Return a Domain1D instance representing a non-reacting +% surface. +m = Domain1D(3); +if nargin == 0 + setID(m,'surface'); +else + setID(m,id); +end diff --git a/Cantera/matlab/cantera/1D/SymmPlane.m b/Cantera/matlab/cantera/1D/SymmPlane.m new file mode 100644 index 000000000..8e164e5d5 --- /dev/null +++ b/Cantera/matlab/cantera/1D/SymmPlane.m @@ -0,0 +1,11 @@ +function m = SymmPlane(id) +% SYMMPLANE - Return a Domain1D instance representing a symmetry plane. +% +m = Domain1D(4); +if nargin == 0 + setID(m,'symmetry_plane'); +else + setID(m,id); +end + + diff --git a/Cantera/matlab/cantera/private/ctmethods.cpp b/Cantera/matlab/cantera/private/ctmethods.cpp index 99392c3d6..02399690b 100644 --- a/Cantera/matlab/cantera/private/ctmethods.cpp +++ b/Cantera/matlab/cantera/private/ctmethods.cpp @@ -1,8 +1,22 @@ +/** + * @file ctmethods.cpp + * + * The interface between the MATLAB environment and the C++ Cantera + * kernel is through a single MEX file. This is top-level driver for + * the MEX file. + * + * This file handles the methods of all Cantera MATLAB classes. The + * class is indicated by the first parameter in the call from MATLAB. + */ #include "mex.h" #include "../../../clib/src/ct.h" #include "ctmatutils.h" +namespace Cantera { + void setMatlabMode(bool m); +} + const int NO_CLASS = 0; const int XML_CLASS = 10; const int THERMO_CLASS = 20; @@ -12,24 +26,53 @@ const int TRANSPORT_CLASS = 50; const int REACTOR_CLASS = 60; const int WALL_CLASS = 70; const int FLOWDEVICE_CLASS = 80; +const int ONEDIM_CLASS = 90; + +void ctfunctions( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + +void xmlmethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + +void thermomethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + +void phasemethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + +void kineticsmethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + +void transportmethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + +void reactormethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + +void wallmethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + +void flowdevicemethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + +void onedimmethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + -void ctfunctions( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); -void xmlmethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); -void thermomethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); -void phasemethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); -void kineticsmethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); -void transportmethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); -void reactormethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); -void wallmethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); -void flowdevicemethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); extern "C" { void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) { + // flag specifying the class int iclass = getInt(prhs[0]); + + // specifies that function writelog should write to MATLAB + Cantera::setMatlabMode(true); + // Hand off to the appropriate routine, based on the + // value of the first parameter switch (iclass) { case NO_CLASS: ctfunctions(nlhs, plhs, nrhs, prhs); break; @@ -49,6 +92,8 @@ extern "C" { wallmethods(nlhs, plhs, nrhs, prhs); break; case FLOWDEVICE_CLASS: flowdevicemethods(nlhs, plhs, nrhs, prhs); break; + case ONEDIM_CLASS: + onedimmethods(nlhs, plhs, nrhs, prhs); break; default: mexErrMsgTxt("unknown class"); } diff --git a/Cantera/matlab/cantera/private/onedimmethods.cpp b/Cantera/matlab/cantera/private/onedimmethods.cpp new file mode 100644 index 000000000..9563892da --- /dev/null +++ b/Cantera/matlab/cantera/private/onedimmethods.cpp @@ -0,0 +1,368 @@ + +#include "mex.h" +#include "ctmatutils.h" +#include "../../../clib/src/ctonedim.h" + +#include +using namespace std; + +void onedimmethods( int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] ) { + double vv; + int job = getInt(prhs[2]); + int n, m; + double *dom_ids, *h; + int indx; + char *nm; + + int dom; + dom = getInt(prhs[1]); + + int idom, icomp, localPoint; + if (job < 10) { + int ph, kin, tr, nd, sz, k, *ptrs; + + switch (job) { + + // construct a new stagnation flow instance + case 1: + checkNArgs(6, nrhs); + ph = getInt(prhs[3]); + kin = getInt(prhs[4]); + tr = getInt(prhs[5]); + indx = stflow_new(ph, kin, tr); + break; + + // construct a new Inlet1D instance + case 2: + checkNArgs(3, nrhs); + indx = inlet_new(); + break; + + // construct a new Surf1D instance + case 3: + checkNArgs(3, nrhs); + indx = surf_new(); + break; + + // construct a new Symm1D instance + case 4: + checkNArgs(3, nrhs); + indx = symm_new(); + break; + + // construct a new Outlet1D instance + case 5: + checkNArgs(3, nrhs); + indx = outlet_new(); + break; + + // construct a new Sim1D instance + case 8: + checkNArgs(5, nrhs); + nd = getInt(prhs[3]); + dom_ids = mxGetPr(prhs[4]); + m = mxGetM(prhs[4]); + n = mxGetN(prhs[4]); + if (m == 1) + sz = n; + else + sz = m; + if (sz != nd) + mexErrMsgTxt("wrong size for domain array"); + + ptrs = new int[sz]; + for (k = 0; k < sz; k++) { + ptrs[k] = int(dom_ids[k]); + } + indx = sim1D_new(sz, ptrs); + delete[] ptrs; + break; + + default: + mexErrMsgTxt("onedimmethods: unknown object type"); + } + + plhs[0] = mxCreateNumericMatrix(1,1,mxDOUBLE_CLASS,mxREAL); + h = mxGetPr(plhs[0]); + *h = double(indx); + if (indx < 0) reportError(); + return; + } + + + // methods + + else if (job < 40) { + + int k; + + switch (job) { + + case 10: + checkNArgs(3, nrhs); + vv = domain_del(dom); break; + case 11: + checkNArgs(3, nrhs); + vv = domain_nComponents(dom); break; + case 12: + checkNArgs(3, nrhs); + vv = domain_type(dom); break; + case 13: + checkNArgs(3, nrhs); + vv = domain_index(dom); + if (vv >= 0.0) vv += 1.0; break; + case 14: + checkNArgs(3, nrhs); + vv = domain_nPoints(dom); break; + case 15: + checkNArgs(3, nrhs); + vv = bdry_temperature(dom); break; + case 16: + checkNArgs(4, nrhs); + k = getInt(prhs[3]); + vv = bdry_massFraction(dom, k); break; + case 17: + checkNArgs(3, nrhs); + vv = bdry_mdot(dom); break; + case 18: + checkNArgs(4, nrhs); + nm = getString(prhs[3]); + vv = domain_componentIndex(dom, nm) ; + if (vv >= 0.0) vv += 1.0; break; + case 19: + checkNArgs(4, nrhs); + localPoint = getInt(prhs[3]) - 1; + vv = domain_grid(dom, localPoint); break; + case 30: + checkNArgs(6, nrhs); + idom = getInt(prhs[3]) - 1; + icomp = getInt(prhs[4]) - 1; + localPoint = getInt(prhs[5]) - 1; + vv = sim1D_value(dom, idom, icomp, localPoint); + break; + case 31: + checkNArgs(6, nrhs); + idom = getInt(prhs[3]) - 1; + icomp = getInt(prhs[4]) - 1; + localPoint = getInt(prhs[5]) - 1; + vv = sim1D_workValue(dom, idom, icomp, localPoint); + break; + default: + mexErrMsgTxt("unknown job"); + } + plhs[0] = mxCreateNumericMatrix(1,1,mxDOUBLE_CLASS,mxREAL); + double *h = mxGetPr(plhs[0]); + *h = vv; + if ((job != 30) && (vv == -1.0)) reportError(); + return; + } + + else if (job < 50) { + int iok = -1; + int buflen, icomp; + char* output_buf; + switch (job) { + case 40: + icomp = getInt(prhs[3]) - 1; + buflen = 40; + output_buf = (char*)mxCalloc(buflen, sizeof(char)); + iok = domain_componentName(dom, icomp, buflen, output_buf); + break; + default: + iok = -1; + } + if (iok >= 0) { + plhs[0] = mxCreateString(output_buf); + return; + } + else { + mexErrMsgTxt("error or unknown method."); + return; + } + } + + + // set parameters + + else { + + int iok = -1; + double *lower, *upper, *rtol, *atol, *grid, *pos, *values, + mdot, t, p, val, *temp, ratio, slope, curve, tstep, *dts, + rdt; + int nlower, nupper, nr, na, npts, np, comp, localPoint, idom, + loglevel, refine_grid, n, flag, itime, ns, *nsteps, icount; + char *xstr, *fname, *id, *desc, *name; + switch (job) { + case 51: + checkNArgs(5, nrhs); + lower = mxGetPr(prhs[3]); + nlower = mxGetM(prhs[3]) * mxGetN(prhs[3]); + upper = mxGetPr(prhs[4]); + nupper = mxGetM(prhs[4]) * mxGetN(prhs[4]); + iok = domain_setBounds(dom, nlower, lower, nupper, upper); + break; + case 52: + checkNArgs(6, nrhs); + rtol = mxGetPr(prhs[3]); + nr = mxGetM(prhs[3]) * mxGetN(prhs[3]); + atol = mxGetPr(prhs[4]); + na = mxGetM(prhs[4]) * mxGetN(prhs[4]); + itime = getInt(prhs[5]); + iok = domain_setTolerances(dom, nr, rtol, na, atol, itime); + break; + case 53: + checkNArgs(4, nrhs); + grid = mxGetPr(prhs[3]); + npts = mxGetM(prhs[3]) * mxGetN(prhs[3]); + iok = domain_setupGrid(dom, npts, grid); + break; + case 54: + id = getString(prhs[3]); + iok = domain_setID(dom, id); + break; + case 60: + checkNArgs(4, nrhs); + mdot = getDouble(prhs[3]); + iok = bdry_setMdot(dom, mdot); + break; + case 61: + checkNArgs(4, nrhs); + t = getDouble(prhs[3]); + iok = bdry_setTemperature(dom, t); + break; + case 62: + checkNArgs(4, nrhs); + xstr = getString(prhs[3]); + iok = bdry_setMoleFractions(dom, xstr); + break; + + case 63: + checkNArgs(4, nrhs); + p = getDouble(prhs[3]); + iok = stflow_setPressure(dom, p); + break; + case 64: + checkNArgs(5, nrhs); + pos = mxGetPr(prhs[3]); + temp = mxGetPr(prhs[4]); + n = mxGetM(prhs[3])*mxGetN(prhs[3]); + iok = stflow_setFixedTempProfile(dom, n, pos, temp); + break; + case 65: + checkNArgs(4, nrhs); + flag = getInt(prhs[3]); + iok = stflow_solveSpeciesEqs(dom, flag); + break; + case 66: + checkNArgs(4, nrhs); + flag = getInt(prhs[3]); + iok = stflow_solveEnergyEqn(dom, flag); + break; + + case 100: + checkNArgs(7, nrhs); + idom = getInt(prhs[3]) - 1; + comp = getInt(prhs[4]) - 1; + localPoint = getInt(prhs[5]) -1; + val = getDouble(prhs[6]); + iok = sim1D_setValue(dom, idom, comp, localPoint, val); + break; + case 101: + checkNArgs(7, nrhs); + idom = getInt(prhs[3]) - 1; + comp = getInt(prhs[4]) - 1; + pos = mxGetPr(prhs[5]); + values = mxGetPr(prhs[6]); + np = mxGetM(prhs[5])*mxGetN(prhs[5]); + iok = sim1D_setProfile(dom, idom, comp, np, pos, values); + break; + case 102: + checkNArgs(6, nrhs); + idom = getInt(prhs[3]) - 1; + comp = getInt(prhs[4]) - 1; + val = getDouble(prhs[5]); + iok = sim1D_setFlatProfile(dom, idom, comp, val); + break; + case 103: + checkNArgs(4, nrhs); + fname = getString(prhs[3]); + iok = sim1D_showSolution(dom, fname); + break; + case 104: + checkNArgs(5, nrhs); + loglevel = getInt(prhs[3]); + refine_grid = getInt(prhs[4]); + iok = sim1D_solve(dom, loglevel, refine_grid); + break; + case 105: + checkNArgs(4, nrhs); + loglevel = getInt(prhs[3]); + iok = sim1D_refine(dom, loglevel); + break; + case 106: + checkNArgs(7, nrhs); + idom = getInt(prhs[3]) - 1; + ratio = getDouble(prhs[4]); + slope = getDouble(prhs[5]); + curve = getDouble(prhs[6]); + iok = sim1D_setRefineCriteria(dom, idom, ratio, slope, curve); + break; + case 107: + iok = 0; + checkNArgs(6, nrhs); + fname = getString(prhs[3]); + id = getString(prhs[4]); + desc = getString(prhs[5]); + iok = sim1D_save(dom, fname, id, desc); + break; + case 108: + checkNArgs(3, nrhs); + iok = sim1D_writeStats(dom); + break; + case 109: + checkNArgs(4, nrhs); + name = getString(prhs[3]); + iok = sim1D_domainIndex(dom, name); + if (iok >= 0) iok++; + break; + case 110: + checkNArgs(3, nrhs); + iok = sim1D_del(dom); + break; + case 111: + iok = 0; + checkNArgs(5, nrhs); + fname = getString(prhs[3]); + id = getString(prhs[4]); + iok = sim1D_restore(dom, fname, id); + break; + case 112: + tstep = getDouble(prhs[3]); + ns = getInt(prhs[4]); + dts = mxGetPr(prhs[5]); + nsteps = new int[ns]; + for (n = 0; n < ns; n++) { + nsteps[n] = int(dts[n]); + } + iok = sim1D_setTimeStep(dom, tstep, ns, nsteps); + delete[] nsteps; + break; + case 113: + checkNArgs(5, nrhs); + rdt = getDouble(prhs[3]); + icount = getInt(prhs[4]); + iok = sim1D_eval(dom, rdt, icount); + break; + default: + mexPrintf(" job = %d ",job); + mexErrMsgTxt("unknown parameter"); + } + if (iok < 0) reportError(); + plhs[0] = mxCreateNumericMatrix(1,1,mxDOUBLE_CLASS,mxREAL); + double *h = mxGetPr(plhs[0]); + *h = double(iok); + return; + } +} diff --git a/Cantera/matlab/setup_matlab.py b/Cantera/matlab/setup_matlab.py index 7852ed652..780b4dd9d 100644 --- a/Cantera/matlab/setup_matlab.py +++ b/Cantera/matlab/setup_matlab.py @@ -20,6 +20,7 @@ mex private/ctmethods.cpp private/ctfunctions.cpp ... private/thermomethods.cpp private/kineticsmethods.cpp ... private/transportmethods.cpp private/reactormethods.cpp ... private/wallmethods.cpp private/flowdevicemethods.cpp ... + private/onedimmethods.cpp private/write.cpp ... """+' -L'+libdir+' '+libs+'\n'+"""disp('done.'); """) fb.close() diff --git a/Cantera/python/Cantera/Flow.py b/Cantera/python/Cantera/Flow.py index 82e0bebad..fdeac6e38 100755 --- a/Cantera/python/Cantera/Flow.py +++ b/Cantera/python/Cantera/Flow.py @@ -25,7 +25,12 @@ _geom = {'Axisymmetric':0, 'Axi':0, 'Planar':1} class Flow1D: """ One-dimensional reacting flows. - Class Flow1D models several types of steady 'one dimensional' + Class Flow1D simulates a one-dimensional flow domain. To use + Flow1D objects, they must be installed in a container, which is an + object of class OneDim. Each Flow1D domain must be terminated by + boundary domains. + + Class Flow1D can model several types of steady 'one dimensional' reacting flows. The flows are one-dimensional in the sense that the governing equations for the steady-state solution can be cast in the form of a set of ordinary differential equations in one @@ -140,7 +145,7 @@ class Flow1D: # fraction and velocity profiles until the energy equation is # enabled. self.refine_components = range(4,4+self.nsp) - if self.type == 0: self.refine_components.append([0,1]) + self.refiner = refine.Refiner(components = self.refine_components, names = self.names) @@ -331,9 +336,9 @@ class Flow1D: def __repr__(self): return self.show() - def show(self): + def show(self, x = None): fname = tempfile.mktemp('.dat') - x = self.x + x = self.x _cantera.flow_showsolution(self.__flow_id, fname, x) f = open(fname,'r') y = f.readlines() @@ -519,7 +524,8 @@ class Flow1D: r.components = range(4,self.nsp+4) if self.energy: r.components.append(2) - + if self.type == 0: + r.components += [0,1] #dsave = r.delta #while 1 > 0: diff --git a/Cantera/python/Cantera/OneDim.py b/Cantera/python/Cantera/OneDim.py index 4ac6b10e1..f836d3d73 100755 --- a/Cantera/python/Cantera/OneDim.py +++ b/Cantera/python/Cantera/OneDim.py @@ -49,55 +49,74 @@ class OneDim: def __init__(self, domains): + """Create a new one-didmensional model from a list of domains. """ + + # instance variables self._size = [] self._start = [] self._end = [] - self._domain = [] - self._flow = [] + self._domain = [] # all domains + self._flow = [] # extended domains self._shape = [] self._loc = 0 self._opt = {} self.time = 0.0 self.x = array([0.0,],'d') self._surf = [] - dtype = [] - dlist = [] self.npts = [] + + # local variables + dtype = [] # list of integer domain types + dlist = [] # list of integer domain ids + + # add each domain for d in domains: + if d.domainType == 0: self.addFlow(d) - dtype.append(0) dlist.append(d.flow_id()) self.npts.append(d.nPoints()) + elif d.domainType == 1: self.addSurface(d) - dtype.append(1) dlist.append(d.surf_id()) self.npts.append(1) + elif d.domainType == 2: self.addBoundary(d) - dtype.append(2) dlist.append(d.bndry_id()) self.npts.append(1) + else: raise 'unknown domain type' + dtype.append(d.domainType) + self.__onedim_id = _cantera.onedim_new(len(dlist), - array(dlist,'i'), - array(dtype,'i')) + array(dlist,'i'), + array(dtype,'i')) self.collect() self.restoreDefaults(); self.ienergy = 0 self.ts_jac_age = 50 + def __del__(self): + """Delete the kernel object. + + This does not delete the individual domains.""" _cantera.onedim_del(self.__onedim_id) + def addFlow(self, flow): - + + # add the domain to the list of all domains and to the list of + # extended domains self._domain.append(flow) self._flow.append(flow) + + # set the index of this domain flow.index = len(self._domain) - 1 - + np, nv = flow.shape() self._shape.append((np,nv)) self._size.append(np*nv) @@ -129,6 +148,7 @@ class OneDim: """ for i in range(len(self._domain)): self._domain[i].x = self.solution(i) + def solution(self, i): """ Return the solution array for domain i. @@ -167,6 +187,7 @@ class OneDim: self._loc += np*nv self._end.append(self._loc) + def addBoundary(self, b): """Add a boundary domain.""" #self._surf.append(surf) @@ -196,7 +217,7 @@ class OneDim: iok = _cantera.onedim_solve(self.__onedim_id, self.x, self.xnew, loglevel) - if loglevel > 0: print _cantera.readlog() + #if loglevel > 0: print _cantera.readlog() if iok >= 0: _cantera.copy(size(self.x),self.xnew,self.x) elif iok > -10: @@ -258,6 +279,7 @@ class OneDim: self._opt = {} self.setOptions( max_jac_age = 20, + ts_jac_age = 30, timestep = 1.e-6, min_timestep = 1.e-12, max_timestep = 0.1, @@ -369,7 +391,7 @@ class OneDim: def c_timeStep(self, nsteps, dt, loglevel = 0): dtnew = _cantera.onedim_timestep(self.__onedim_id, nsteps, dt, self.x, self.xnew, loglevel) - print _cantera.readlog() + #print _cantera.readlog() return dtnew @@ -381,8 +403,9 @@ class OneDim: loglevel -- controls amount of printed diagnostics """ - self.setNewtonOptions(max_jac_age = self.ts_jac_age) - + self.setNewtonOptions(max_jac_age = self._opt['ts_jac_age']) + print 'max jac age = ',self._opt['ts_jac_age'] + if loglevel > 0: print_heading('Begin time integration.\n\n') print(' step size (s) log10(ss) ') @@ -400,8 +423,8 @@ class OneDim: m = self.newton_solve(loglevel-1) self.time += dt n += 1 - if m == 100: dt *= 1.5 - #if m > 0: dt *= 1.5 + if m == 100: + dt *= 1.5 if dt > maxdt: dt = maxdt if loglevel > 0: print @@ -425,7 +448,7 @@ class OneDim: def showStatistics(self): _cantera.onedim_writestats(self.__onedim_id) - print _cantera.readlog() + #print _cantera.readlog() def save(self, filename, id, desc=""): @@ -449,6 +472,6 @@ class OneDim: fn = filename + '.xml' _cantera.onedim_save(self.__onedim_id, fn, id, desc, self.x) - print _cantera.readlog() + #print _cantera.readlog() diff --git a/Cantera/python/Cantera/flame.py b/Cantera/python/Cantera/flame.py index a0ef89413..2365496c0 100755 --- a/Cantera/python/Cantera/flame.py +++ b/Cantera/python/Cantera/flame.py @@ -1,5 +1,6 @@ from Cantera import OneAtm +from Cantera.exceptions import CanteraError from Cantera.Flow import Flow1D from Cantera.boundaries1D import Inlet1D, Surf1D, Symm1D from Numeric import array, zeros, arrayrange @@ -306,15 +307,16 @@ class StagnationFlame: if grid == None: grid = dx * array([0.0, 0.01, 0.03, 0.1, 0.3, 0.6, 1.0]) - + self.__flow = Flow1D(flow_type = 'Stag', gas = gas, grid = grid, pressure = self.p) self.__left = Inlet1D() + self.__right = Surf1D() + self.__container = OneDim([self.__left, self.__flow, self.__right]) self.start = 0 - # get the compositions of the fuel and oxidizer streams, and # calculate the fuel/oxidizer ratio for stoichiometric @@ -449,8 +451,13 @@ class StagnationFlame: self.__flow.setTolerances(u = v, V = v, T = v, Y = v) elif o == 'max_jac_age': self.__container.setOptions(max_jac_age = v) + elif o == 'jac_age': + self.__container.setOptions(max_jac_age = v[0]) + self.__container.setOptions(ts_jac_age = v[1]) elif o == 'timesteps': self.__container.setOptions(nsteps = v[0], timestep = v[1]) + else: + raise CanteraError("unknown option: "+o) def solve(self, loglevel = 0): if not self.start: @@ -458,11 +465,6 @@ class StagnationFlame: self.start = 1 solve(self.__container, loglevel = loglevel, refine_grid = 1) -## def esolve(self, loglevel = 0, efactor = 1.0e4): -## if not self.start: -## self.setEquilProducts() -## self.start = 1 -## esolve(self.__container, efactor = efactor, loglevel = loglevel, refine_grid = 1) def save(self, soln, desc, file = 'flame.xml'): self.__container.save(file, soln, desc) diff --git a/Cantera/python/Cantera/rxnpath.py b/Cantera/python/Cantera/rxnpath.py index 8b281963b..4d2a252bc 100755 --- a/Cantera/python/Cantera/rxnpath.py +++ b/Cantera/python/Cantera/rxnpath.py @@ -138,6 +138,7 @@ class PathBuilder: "buildlog", diagram.id(), 1) if format == "dot": diagram.write(0, dotfile) + diagram.write(1, "rp.txt") elif format == "plain": diagram.write(1, dotfile) diff --git a/Cantera/python/examples/mix1.py b/Cantera/python/examples/mix1.py new file mode 100644 index 000000000..495019a53 --- /dev/null +++ b/Cantera/python/examples/mix1.py @@ -0,0 +1,86 @@ +# Mixing two streams. + +# Since reactors can have multiple inlets and outlets, they can be +# used to implement mixers, splitters, etc. In this example, air and +# methane are mixed in stoichiometric proportions. Due to the low +# temperature, no reactions occur. Note that the air stream and the +# methane stream use *different* reaction mechanisms, with different +# numbers of species and reactions. When gas flows from one reactor or +# reservoir to another one with a different reaction mechanism, +# species are matched by name. If the upstream reactor contains a +# species that is not present in the downstream reaction mechanism, it +# will be ignored. In general, reaction mechanisms for downstream +# reactors should contain all species that might be present in any +# upstream reactor. +# +#----------------------------------------------------------------------- + +from Cantera import * +from Cantera.Reactor import Reactor, Reservoir, MassFlowController, Valve + + +# Use air for stream a. Note that the Air() function does not set the +# composition correctly; thus, we need to explicitly set the +# composition to that of air. +gas_a = Air() +gas_a.setState_TPX(300.0, OneAtm, 'O2:0.21, N2:0.78, AR:0.01') +rho_a = gas_a.density() + + +# Use GRI-Mech 3.0 for stream b (methane) and for the mixer. If it is +# desired to have a pure mixer, with no chemistry, use instead a +# reaction mechanism for gas_b that has no reactions. +gas_b = GRI30() +gas_b.setState_TPX(300.0, OneAtm, 'CH4:1') +rho_b = gas_b.density() + + +# Create reservoirs for the two inlet streams and for the outlet +# stream. The upsteam reservoirs could be replaced by reactors, which +# might themselves be connected to reactors further upstream. The +# outlet reservoir could be replaced with a reactor with no outlet, if +# it is desired to integrate the composition leaving the mixer in +# time, or by an arbitrary network of downstream reactors. +res_a = Reservoir(gas_a) +res_b = Reservoir(gas_b) +downstream = Reservoir(gas_b) + + +# Create a reactor for the mixer. A reactor is required instead of a +# reservoir, since the state will change with time if the inlet mass +# flow rates change or if there is chemistry occurring. +mixer = Reactor(gas_b) + + +# create two mass flow controllers connecting the upstream reservoirs +# to the mixer, and set their mass flow rates to values corresponding +# to stoichiometric combustion. +mfc1 = MassFlowController(res_a, mixer) +mfc1.setMassFlowRate(rho_a*2.5/0.21) + +mfc2 = MassFlowController(res_b, mixer) +mfc2.setMassFlowRate(rho_b*1.0) + + +# connect the mixer to the downstream reservoir with a valve. +outlet = Valve(mixer, downstream) +outlet.setValveCoeff(1.0) + + +# Since the mixer is a reactor, we need to integrate in time to reach +# steady state. A few residence times should be enough. +t = 0.0 +for n in range(30): + tres = mixer.mass()/(mfc1.massFlowRate() + mfc2.massFlowRate()) + t += 0.5*tres + mixer.advance(t) + print '%14.5g %14.5g %14.5g %14.5g %14.5g' % (t, mixer.temperature(), + mixer.enthalpy_mass(), + mixer.pressure(), + mixer.massFraction('CH4')) + +# view the state of the gas in the mixer +gas_b.setState_TPY(mixer.temperature(), mixer.pressure(), + mixer.massFractions()) +print gas_b + diff --git a/Cantera/python/examples/mix2.py b/Cantera/python/examples/mix2.py new file mode 100644 index 000000000..7e98dd39d --- /dev/null +++ b/Cantera/python/examples/mix2.py @@ -0,0 +1,93 @@ +# Mixing two streams with reaction. This is the same as mix1.py, +# except that a source of H atoms is added to ignite the fuel/air +# mixture. Once ignited, the flow of H atoms is stopped. + +import math +from Cantera import * +from Cantera.Reactor import Reactor, Reservoir, MassFlowController, Valve + + +# Use air for stream a. Note that the Air() function does not set the +# composition correctly; thus, we need to explicitly set the +# composition to that of air. +gas_a = Air() +gas_a.setState_TPX(300.0, OneAtm, 'O2:0.21, N2:0.78, AR:0.01') +rho_a = gas_a.density() + + +# Use GRI-Mech 3.0 for stream b (methane) and for the mixer. If it is +# desired to have a pure mixer, with no chemistry, use instead a +# reaction mechanism for gas_b that has no reactions. +gas_b = GRI30() +gas_b.setState_TPX(300.0, OneAtm, 'CH4:1') +rho_b = gas_b.density() + + +# Create reservoirs for the two inlet streams and for the outlet +# stream. The upsteam reservoirs could be replaced by reactors, which +# might themselves be connected to reactors further upstream. The +# outlet reservoir could be replaced with a reactor with no outlet, if +# it is desired to integrate the composition leaving the mixer in +# time, or by an arbitrary network of downstream reactors. +res_a = Reservoir(gas_a) +res_b = Reservoir(gas_b) +downstream = Reservoir(gas_b) + + +# Create a reactor for the mixer. A reactor is required instead of a +# reservoir, since the state will change with time if the inlet mass +# flow rates change or if there is chemistry occurring. +mixer = Reactor(gas_b) + + +# create two mass flow controllers connecting the upstream reservoirs +# to the mixer, and set their mass flow rates to values corresponding +# to stoichiometric combustion. +mfc1 = MassFlowController(res_a, mixer) +mfc1.setMassFlowRate(rho_a*2.5/0.21) + +mfc2 = MassFlowController(res_b, mixer) +mfc2.setMassFlowRate(rho_b*1.0) + + +# connect the mixer to the downstream reservoir with a valve. +outlet = Valve(mixer, downstream) +outlet.setValveCoeff(1.0) + + +# add an igniter to ignite the mixture. The 'igniter' consists of a +# stream of pure H. +gas_c = IdealGasMix('h2o2.xml') +gas_c.setState_TPX(300.0, OneAtm, 'H:1') +igniter = Reactor(gas_c) + +mfc3 = MassFlowController(igniter, mixer) +mfc3.setMassFlowRate(0.05) + + + +# Since the mixer is a reactor, we need to integrate in time to reach +# steady state. A few residence times should be enough. +t = 0.0 +for n in range(30): + tres = mixer.mass()/(mfc1.massFlowRate() + mfc2.massFlowRate()) + tnow = t + t += 0.5*tres + mixer.advance(t) + + # if ignited, turn the igniter off. + # We also need to restart the integration in this case. + if mixer.temperature() > 1200.0: + mfc3.setMassFlowRate(0.0) + mixer.setInitialTime(t) + + print '%14.5g %14.5g %14.5g %14.5g %14.5g' % (t, mixer.temperature(), + mixer.enthalpy_mass(), + mixer.pressure(), + mixer.massFraction('CH4')) +gas_b.setState_TPY(mixer.temperature(), mixer.pressure(), mixer.massFractions()) + +# view the state of the gas in the mixer +gas_b.setState_TPY(mixer.temperature(), mixer.pressure(), + mixer.massFractions()) +print gas_b diff --git a/Cantera/python/examples/stflame1.py b/Cantera/python/examples/stflame1.py index b79501a69..9053d6431 100644 --- a/Cantera/python/examples/stflame1.py +++ b/Cantera/python/examples/stflame1.py @@ -56,6 +56,7 @@ flame.set(mdot = 0.1, # turn the energy equation off (default) flame.set(energy = 'off') +flame.show() # solve the flame, with output level 1 flame.solve(1) diff --git a/Cantera/python/setup.py.in b/Cantera/python/setup.py.in index 0cfdbdd2f..4eda42107 100644 --- a/Cantera/python/setup.py.in +++ b/Cantera/python/setup.py.in @@ -36,7 +36,7 @@ try: packages = ["","Cantera","MixMaster","MixMaster.Units"], ext_modules=[ Extension("Cantera._cantera", - ["src/pycantera.cpp"], + ["src/pycantera.cpp", "src/writelog.cpp"], include_dirs=["../../build/include", "src", "../clib/src"], library_dirs = ["@buildlib@"], libraries = libs) diff --git a/Cantera/src/Flow1D.h b/Cantera/src/Flow1D.h deleted file mode 100755 index 5ac26e174..000000000 --- a/Cantera/src/Flow1D.h +++ /dev/null @@ -1,82 +0,0 @@ - -deprecated - - -#ifndef CT_FLOW1D_H -#define CT_FLOW1D_H - -#include "MultiDomain.h" -#include "MultiJac.h" -#include "MultiNewton.h" - -#include "Jac1D.h" -#include "Newton1D.h" - -#include "Surf1D.h" - -namespace Cantera { - - /** - * Container class for multiple-domain 1D problems. - */ - class OneDim : public MultiDomain { - public: - OneDim() : m_jac(0), m_newt(0) { - m_newt = new MultiNewton(1); - } - - virtual void addDomain(Resid1D* d) { - MultiDomain::addDomain(d); - m_newt->resize(size()); - delete m_jac; - m_jac = 0; - m_jac = new MultiJac(*this); - m_jac_ok = false; - int nd = m_dom.size(); - for (int i = 0; i < nd; i++) - m_dom[i]->setJac(m_jac); - } - - virtual ~OneDim() { - delete m_jac; - delete m_newt; - } - - MultiJac& jacobian() { return *m_jac; } - MultiNewton& newton() { return *m_newt; } - - - virtual int solve(doublereal* x, doublereal* xnew, int loglevel) { - if (!m_jac) { - cout << "creating new jac.." << endl; - m_jac = new MultiJac(*this); - m_jac_ok = false; - int nd = m_dom.size(); - for (int i = 0; i < nd; i++) - m_dom[i]->setJac(m_jac); - cout << "done" << endl; - } - - if (!m_jac_ok) { - cout << "eval jac" << endl; - eval(-1, x, xnew, m_rdt); - m_jac->eval(x, xnew, m_rdt); - m_jac_ok = true; - } - return m_newt->solve(x, xnew, *this, *m_jac, loglevel); - } - - - - protected: - MultiJac* m_jac; - MultiNewton* m_newt; - Jac1D* m_jac1; - Newton1D* m_newt1; - }; - -} - -#endif - - diff --git a/Cantera/src/Jac2.h b/Cantera/src/Jac2.h deleted file mode 100755 index 65bedf19f..000000000 --- a/Cantera/src/Jac2.h +++ /dev/null @@ -1,198 +0,0 @@ -/** - * - * @file Jac2.h - * - * >>>>> Under construction! <<<<< - * - * $Author$ - * $Date$ - * $Revision$ - * - * Copyright 2002 California Institute of Technology - * - */ - -#ifndef CT_JAC2_H -#define CT_JAC2_H - - -#include "BandMatrix.h" -#include "ctlapack.h" -#include "../ext/math/gmres.h" -#include "stringUtils.h" -#include "Array.h" -#include "time.h" - -namespace Cantera { - - - /** - * Residual function evaluator for a one-dimensional problem. - */ - class ResidFunc2 { - public: - - /** - * Constructor. - * @param nv Number of variables at each grid point. - * @param points Number of grid points. - */ - ResidFunc2(int nv, int points) { - m_nv = nv; - m_points = points; - m_max.resize(m_nv, 0.0); - m_min.resize(m_nv, 0.0); - m_rtol.resize(m_nv, 0.0); - m_atol.resize(m_nv, 0.0); - m_slast.resize(m_nv, m_points); - m_rdt = 0.0; - - // m_soln.resize(m_nv, m_points, 0.0); - //m_resid.resize(m_nv, m_points, 0.0); - } - - /// Destructor. - virtual ~ResidFunc2(){} - - /// Number of components at each grid point. - int nComponents() const { return m_nv; } - - /// Number of grid points. - int nPoints() const { return m_points; } - - /// Name of the nth component. - virtual string componentName(int n) const { - return "component " + int2str(n); } - - void setBounds(const vector_fp& lower, const vector_fp& upper) { - if (lower.size() != m_nv || upper.size() != m_nv) - throw CanteraError("ResidFunc2::setBounds", - "wrong array size for solution bounds"); - m_max = upper; - m_min = lower; - } - - void setTolerances(vector_fp& rtol, vector_fp& atol) { - m_rtol = rtol; - m_atol = atol; - } - - doublereal rtol(int n) { return m_rtol[n]; } - doublereal atol(int n) { return m_atol[n]; } - - doublereal upperBound(int n) const { - return m_max[n]; - } - - doublereal lowerBound(int n) const { - return m_min[n]; - } - - - void initTimeInteg(doublereal dt, const Array2D& x0) { - m_slast = x0; - m_rdt = 1.0/dt; - } - - void setSteadyMode() { - m_rdt = 0.0; - } - - bool steady() { return (m_rdt == 0.0); } - bool transient() { return (m_rdt != 0.0); } - - /// Evaluate the residual function at point j. - virtual void eval(int j, Array2D& x, Array2D& r) { - throw CanteraError("ResidFunc2::eval", - "residual function not defined."); - } - - void evalss(Array2D& x, Array2D& r) { - doublereal rdt_save = m_rdt; - m_rdt = 0.0; - eval(-1, x, r); - m_rdt = rdt_save; - } - - protected: - int m_nv; - int m_points; - // Array2D m_soln; - //Array2D m_resid; - Array2D m_slast; - doublereal m_rdt; - vector_fp m_max; - vector_fp m_min; - vector_fp m_rtol; - vector_fp m_atol; - - private: - - }; - - - /////////////////////////////////////////////////////////////// - - - /** - * Class Jac2 evaluates the Jacobian of a system of equations - * defined by a residual function of class ResidFunc2. It is - * assumed that the Jacobian is banded. - */ - class Jac2 : public BandMatrix { - - public: - - /** - * Constructor. The residual function defining the system of - * equations must be supplied. - */ - Jac2(ResidFunc2& r); - - /// Destructor. Does nothing. - virtual ~Jac2(){} - - /** - * Evaluate the Jacobian. - */ - void eval(Array2D& x0, Array2D& resid0); - - /** - * Returns the matrix element describing the influence of the - * nth component at point j on the mth equation at point - * i. Due to the assumption of a banded Jacobian, this will be - * zero unless |i - j| <= 1. - */ - doublereal& v(int m, int i, int n, int j) { - return value(i*m_nv + m, j*m_nv + n); - } - - doublereal elapsedTime() const { - return m_elapsed; - } - - int nEvals() const { return m_nevals; } - - int age() const { return m_age; } - - void incrementAge() { m_age++; } - void setAge(int age) { m_age = age; } - - protected: - - ResidFunc2* m_resid; - Array2D m_r1; - int m_nv, m_points; - doublereal m_atol; - doublereal m_elapsed; - int m_nevals; - int m_age; - - private: - - }; -} - -#endif - - diff --git a/Cantera/src/Resid1D.h b/Cantera/src/Resid1D.h deleted file mode 100755 index 676a9babc..000000000 --- a/Cantera/src/Resid1D.h +++ /dev/null @@ -1,237 +0,0 @@ -/** - * @file Resid1D.h - * - * $Author$ - * $Date$ - * $Revision$ - * - * Copyright 2002 California Institute of Technology - * - */ - -#ifndef CT_RESID1D_H -#define CT_RESID1D_H - - -//#include "stringUtils.h" -#include "ctexceptions.h" -#include "xml.h" - -namespace Cantera { - - // domain types - const int cFlowType = 101; - const int cSurfType = 102; - const int cConnectorType = 103; - const int cInletType = 104; - const int cSymmType = 105; - const int cOutletType = 106; - - class MultiJac; - class OneDim; - - - /** - * Base class for single-domain, one-dimensional residual function - * evaluators. - */ - class Resid1D { - public: - - /** - * Constructor. - * @param nv Number of variables at each grid point. - * @param points Number of grid points. - */ - Resid1D(int nv=1, int points=1, - doublereal time = 0.0) : - m_time(time), - m_container(0), - m_index(-1), - m_type(0), - m_iloc(0), - m_jstart(0), - m_left(0), - m_right(0) { - resize(nv, points); - } - - /// Destructor. - virtual ~Resid1D(){} - - /// Domain type flag. - const int domainType() { return m_type; } - - const OneDim& container() const{ return *m_container; } - - /** - * Specify the container object for this domain, and the - * position of this domain in the list. - */ - void setContainer(OneDim* c, int index){ - m_container = c; - m_index = index; - } - - /** Initialize. Base class method does nothing, but may be - * overloaded. - */ - virtual void init(){} - - /** - * Resize the domain to have nv components and np grid points. - */ - virtual void resize(int nv, int np) { - m_nv = nv; - m_max.resize(m_nv, 0.0); - m_min.resize(m_nv, 0.0); - m_rtol.resize(m_nv, 0.0); - m_atol.resize(m_nv, 0.0); - m_points = np; - m_slast.resize(m_nv * m_points, 0.0); - locate(); - } - - /// Number of components at each grid point. - int nComponents() const { return m_nv; } - - /// Number of grid points in this domain. - int nPoints() const { return m_points; } - - /// Name of the nth component. May be overloaded. - virtual string componentName(int n) const { - return "component " + int2str(n); } - - /** - * Set the lower and upper bounds for each solution component. - */ - void setBounds(int nl, const doublereal* lower, - int nu, const doublereal* upper) { - if (nl != m_nv || nu != m_nv) - throw CanteraError("Resid1D::setBounds", - "wrong array size for solution bounds"); - copy(upper, upper + m_nv, m_max.begin()); - copy(lower, lower + m_nv, m_min.begin()); - } - - void setTolerances(int nr, const doublereal* rtol, - int na, const doublereal* atol) { - if (nr != m_nv || na != m_nv) - throw CanteraError("Resid1D::setTolerances", - "wrong array size for solution error tolerances. Size should be "+int2str(m_nv)); - copy(rtol, rtol + m_nv, m_rtol.begin()); - copy(atol, atol + m_nv, m_atol.begin()); - } - - doublereal rtol(int n) { return m_rtol[n]; } - doublereal atol(int n) { return m_atol[n]; } - - doublereal upperBound(int n) const { return m_max[n]; } - doublereal lowerBound(int n) const { return m_min[n]; } - - void initTimeInteg(doublereal dt, const doublereal* x0) { - copy(x0 + loc(), x0 + loc() + size(), m_slast.begin()); - m_rdt = 1.0/dt; - } - - void setSteadyMode() { m_rdt = 0.0; } - - bool steady() { return (m_rdt == 0.0); } - bool transient() { return (m_rdt != 0.0); } - - void needJacUpdate(); - - void evalss(doublereal* x, doublereal* r, integer* mask) { - eval(-1,x,r,mask,0.0); - } - - /** - * Evaluate the residual function at point j. If j < 0, - * evaluate the residual function at all points. - */ - virtual void eval(int j, doublereal* x, doublereal* r, - integer* mask, doublereal rdt=0.0) { - throw CanteraError("Resid1D::eval", - "residual function not defined."); - } - - virtual void update(doublereal* x) {} - - doublereal time() { return m_time;} - void incrementTime(doublereal dt) { m_time += dt; } - size_t index(int n, int j) const { return m_nv*j + n; } - - virtual void setJac(MultiJac* jac){} - virtual void save(XML_Node& o, doublereal* sol) { - throw CanteraError("Resid1D::save","base class method called"); - } - - int size() { return m_nv*m_points; } - - void locate() { - if (m_left) { - m_jstart = m_left->lastPoint() + 1; - m_iloc = m_left->loc() + m_left->size(); - } - else { - m_jstart = 0; - m_iloc = 0; - } - if (m_right) m_right->locate(); - } - - virtual int loc(int j = 0) { return m_iloc; } - - int firstPoint() { return m_jstart; } - int lastPoint() { return m_jstart + m_points - 1; } - - void append(Resid1D* right) { - linkRight(right); - right->linkLeft(this); - } - - void linkLeft(Resid1D* left) { - m_left = left; - locate(); - } - void linkRight(Resid1D* right) { m_right = right; } - - Resid1D* left() { return m_left; } - Resid1D* right() { return m_right; } - - double prevSoln(int n, int j) const{ - return m_slast[m_nv*j + n]; - } - - void setID(const string& s) {m_id = s;} - void setDesc(const string& s) {m_desc = s;} - - virtual void getTransientMask(integer* mask){} - - protected: - - doublereal m_rdt; - int m_nv; - int m_points; - vector_fp m_slast; - doublereal m_time; - vector_fp m_max; - vector_fp m_min; - vector_fp m_rtol; - vector_fp m_atol; - OneDim* m_container; - int m_index; - int m_type; - int m_iloc; - int m_jstart; - Resid1D *m_left, *m_right; - string m_id, m_desc; - - private: - - }; -} - -#endif - - diff --git a/Cantera/src/funcs.cpp b/Cantera/src/funcs.cpp index 4e6712da4..3a33e94f5 100755 --- a/Cantera/src/funcs.cpp +++ b/Cantera/src/funcs.cpp @@ -34,8 +34,10 @@ namespace Cantera { doublereal linearInterp(doublereal x, const vector_fp& xpts, const vector_fp& fpts) { - if (x <= xpts[0]) return fpts[0]; - if (x >= xpts.back()) return fpts.back(); + if (x <= xpts[0]) + return fpts[0]; + if (x >= xpts.back()) + return fpts.back(); const doublereal* loc = lower_bound(xpts.begin(), xpts.end(), x); int iloc = int(loc - xpts.begin()) - 1; doublereal ff = fpts[iloc] + diff --git a/Cantera/src/importCTML.cpp b/Cantera/src/importCTML.cpp index 4b346d217..be852c26f 100755 --- a/Cantera/src/importCTML.cpp +++ b/Cantera/src/importCTML.cpp @@ -45,14 +45,6 @@ namespace Cantera { typedef vector nodeset_t; typedef XML_Node node_t; - static int intValue(string val) { - return atoi(stripws(val).c_str()); - } - - static doublereal fpValue(string val) { - return atof(stripws(val).c_str()); - } - /// Number of reactant molecules static int nReacMolecules(ReactionData& r) { return accumulate(r.rstoich.begin(), r.rstoich.end(), 0); diff --git a/Cantera/src/misc.cpp b/Cantera/src/misc.cpp index f9e2398c3..8617b1d16 100755 --- a/Cantera/src/misc.cpp +++ b/Cantera/src/misc.cpp @@ -18,6 +18,8 @@ #include using namespace std; +//int writeToMatlab(const char* buf); + namespace Cantera { /** @@ -30,7 +32,8 @@ namespace Cantera { */ class Application { public: - Application() : linelen(0), stop_on_error(false), write_log_to_cout(true) {} + Application() : linelen(0), stop_on_error(false), + write_log_to_cout(true), matlab(false) {} virtual ~Application(){} vector inputDirs; vector errorMessage; @@ -40,6 +43,7 @@ namespace Cantera { size_t linelen; bool stop_on_error; bool write_log_to_cout; + bool matlab; map options; }; @@ -249,22 +253,13 @@ namespace Cantera { } } + void setMatlabMode(bool m) { + appinit(); + __app->matlab = m; + } + void write(const string& msg) {cout << msg;} void write(const char* msg) {cout << msg;} - void writelog(const string& msg) { - appinit(); - __app->msglog += msg; - __app->linelen += msg.size(); - if (msg[msg.size()-1] == '\n') __app->linelen = 0; - if (__app->linelen > 70) { - __app->msglog += "\n"; - __app->linelen = 0; - } - if (__app->write_log_to_cout) { - cout << __app->msglog; - clearlog(); - } - } void writelog(const char* msg) {writelog(string(msg));} void getlog(string& s) { appinit(); diff --git a/Cantera/src/oneD/Domain1D.h b/Cantera/src/oneD/Domain1D.h new file mode 100644 index 000000000..6707b77e1 --- /dev/null +++ b/Cantera/src/oneD/Domain1D.h @@ -0,0 +1,433 @@ +/** + * @file Domain1D.h + * + * $Author$ + * $Date$ + * $Revision$ + * + * Copyright 2002 California Institute of Technology + * + */ + +#ifndef CT_DOMAIN1D_H +#define CT_DOMAIN1D_H + + +#include "../ctexceptions.h" +#include "../xml.h" +#include "refine.h" + + +namespace Cantera { + + // domain types + const int cFlowType = 50; + const int cConnectorType = 100; + const int cSurfType = 102; + const int cInletType = 104; + const int cSymmType = 105; + const int cOutletType = 106; + + class MultiJac; + class OneDim; + + + /** + * Base class for one-dimensional domains. + */ + class Domain1D { + public: + + /** + * Constructor. + * @param nv Number of variables at each grid point. + * @param points Number of grid points. + */ + Domain1D(int nv=1, int points=1, + doublereal time = 0.0) : + m_time(time), + m_container(0), + m_index(-1), + m_type(0), + m_iloc(0), + m_jstart(0), + m_left(0), + m_right(0), + m_refiner(0) { + resize(nv, points); + } + + /// Destructor. Does nothing + virtual ~Domain1D(){ delete m_refiner; } + + /// Domain type flag. + const int domainType() { return m_type; } + + /** + * The left-to-right location of this domain. + */ + const int domainIndex() { return m_index; } + + bool isConnector() { return (m_type >= cConnectorType); } + + /** + * The container holding this domain. + */ + const OneDim& container() const { return *m_container; } + + /** + * Specify the container object for this domain, and the + * position of this domain in the list. + */ + void setContainer(OneDim* c, int index){ + m_container = c; + m_index = index; + } + + /** + * Initialize. Base class method does nothing, but may be + * overloaded. + */ + virtual void init(){} + + virtual void setInitialState(doublereal* xlocal = 0){} + virtual void setState(int point, const doublereal* state, doublereal* x) {} + + /** + * Resize the domain to have nv components and np grid points. + * This method is virtual so that subclasses can perform other + * actions required to resize the domain. + */ + virtual void resize(int nv, int np) { + if (nv != m_nv || !m_refiner) { + m_nv = nv; + delete m_refiner; + m_refiner = new Refiner(*this); + } + m_nv = nv; + m_max.resize(m_nv, 0.0); + m_min.resize(m_nv, 0.0); + m_rtol_ss.resize(m_nv, 0.0); + m_atol_ss.resize(m_nv, 0.0); + m_rtol_ts.resize(m_nv, 0.0); + m_atol_ts.resize(m_nv, 0.0); + m_points = np; + m_z.resize(np, 0.0); + m_slast.resize(m_nv * m_points, 0.0); + locate(); + } + + Refiner& refiner() { return *m_refiner; } + + /// Number of components at each grid point. + int nComponents() const { return m_nv; } + + /// Number of grid points in this domain. + int nPoints() const { return m_points; } + + /// Name of the nth component. May be overloaded. + virtual string componentName(int n) const { + return "component " + int2str(n); + } + + int componentIndex(string name) { + int nc = nComponents(); + for (int n = 0; n < nc; n++) { + if (name == componentName(n)) return n; + } + throw CanteraError("Domain1D::componentIndex", + "no component named "+name); + } + + /** + * Set the lower and upper bounds for each solution component. + */ + void setBounds(int nl, const doublereal* lower, + int nu, const doublereal* upper) { + if (nl < m_nv || nu < m_nv) + throw CanteraError("Domain1D::setBounds", + "wrong array size for solution bounds. " + "Size should be at least "+int2str(m_nv)); + copy(upper, upper + m_nv, m_max.begin()); + copy(lower, lower + m_nv, m_min.begin()); + } + + void setTolerances(int nr, const doublereal* rtol, + int na, const doublereal* atol, int ts = 0) { + if (nr < m_nv || na < m_nv) + throw CanteraError("Domain1D::setTolerances", + "wrong array size for solution error tolerances. " + "Size should be at least "+int2str(m_nv)); + if (ts >= 0) { + copy(rtol, rtol + m_nv, m_rtol_ss.begin()); + copy(atol, atol + m_nv, m_atol_ss.begin()); + } + if (ts <= 0) { + copy(rtol, rtol + m_nv, m_rtol_ts.begin()); + copy(atol, atol + m_nv, m_atol_ts.begin()); + } + } + + /// Relative tolerance of the nth component. + doublereal rtol(int n) { return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]); } + + /// Absolute tolerance of the nth component. + doublereal atol(int n) { return (m_rdt == 0.0 ? m_atol_ss[n] : m_atol_ts[n]); } + + /// Upper bound on the nth component. + doublereal upperBound(int n) const { return m_max[n]; } + + /// Lower bound on the nth component + doublereal lowerBound(int n) const { return m_min[n]; } + + + /** + * Prepare to do time stepping with time step dt. Copy the + * internally-stored solution at the last time step to array + * x0. + */ + void initTimeInteg(doublereal dt, const doublereal* x0) { + copy(x0 + loc(), x0 + loc() + size(), m_slast.begin()); + m_rdt = 1.0/dt; + } + + /** + * Prepare to solve the steady-state problem. + * Set the internally-stored reciprocal of the time step to 0,0 + */ + void setSteadyMode() { m_rdt = 0.0; } + + /// True if in steady-state mode + bool steady() { return (m_rdt == 0.0); } + + /// True if not in steady-state mode + bool transient() { return (m_rdt != 0.0); } + + /** + * Set this if something has changed in the governing + * equations (e.g. the value of a constant has been changed, + * so that the last-computed Jacobian is no longer valid. + * Note: see file OneDim.cpp for the implementation of this method. + */ + void needJacUpdate(); + + /** + * Evaluate the steady-state residual at all points, even if in + * transient mode. Used only to print diagnostic output. + */ + void evalss(doublereal* x, doublereal* r, integer* mask) { + eval(-1,x,r,mask,0.0); + } + + /** + * Evaluate the residual function at point j. If j < 0, + * evaluate the residual function at all points. + */ + virtual void eval(int j, doublereal* x, doublereal* r, + integer* mask, doublereal rdt=0.0) { + throw CanteraError("Domain1D::eval", + "residual function not defined."); + } + + /** + * Does nothing. + */ + virtual void update(doublereal* x) {} + + doublereal time() const { return m_time;} + void incrementTime(doublereal dt) { m_time += dt; } + size_t index(int n, int j) const { return m_nv*j + n; } + doublereal value(doublereal* x, int n, int j) const { + return x[index(n,j)]; + } + + virtual void setJac(MultiJac* jac){} + virtual void save(XML_Node& o, doublereal* sol) { + throw CanteraError("Domain1D::save","base class method called"); + } + + int size() const { return m_nv*m_points; } + + /** + * Find the index of the first grid point in this domain, and + * the start of its variables in the global solution vector. + */ + void locate() { + + if (m_left) { + // there is a domain on the left, so the first grid point + // in this domain is one more than the last one on the left + m_jstart = m_left->lastPoint() + 1; + + // the starting location in the solution vector + m_iloc = m_left->loc() + m_left->size(); + } + else { + // this is the left-most domain + m_jstart = 0; + m_iloc = 0; + } + // if there is a domain to the right of this one, then + // repeat this for it + if (m_right) m_right->locate(); + } + + /** + * Location of the start of the local solution vector in the global + * solution vector, + */ + virtual int loc(int j = 0) const { return m_iloc; } + + /** + * The index of the first (i.e., left-most) grid point + * belonging to this domain. + */ + int firstPoint() const { return m_jstart; } + + /** + * The index of the last (i.e., right-most) grid point + * belonging to this domain. + */ + int lastPoint() const { return m_jstart + m_points - 1; } + + /** + * Set the left neighbor to domain 'left.' Method 'locate' is + * called to update the global positions of this domain and + * all those to its right. + */ + void linkLeft(Domain1D* left) { + m_left = left; + locate(); + } + + /** + * Set the right neighbor to domain 'right.' + */ + void linkRight(Domain1D* right) { m_right = right; } + + /** + * Append domain 'right' to this one, and update all links. + */ + void append(Domain1D* right) { + linkRight(right); + right->linkLeft(this); + } + + /** + * Return a pointer to the left neighbor. + */ + Domain1D* left() const { return m_left; } + + /** + * Return a pointer to the right neighbor. + */ + Domain1D* right() const { return m_right; } + + /** + * Value of component n at point j in the previous solution. + */ + double prevSoln(int n, int j) const { + return m_slast[m_nv*j + n]; + } + + /** + * Specify an identifying tag for this domain. + */ + void setID(const string& s) {m_id = s;} + + string id() { + if (m_id != "") return m_id; + else return string("domain ") + int2str(m_index); + } + + /** + * Specify descriptive text for this domain. + */ + void setDesc(const string& s) {m_desc = s;} + const string& desc() { return m_desc; } + + virtual void getTransientMask(integer* mask){} + + virtual void showSolution(ostream& s, const doublereal* x) {} + virtual void showSolution(const doublereal* x) {} + + virtual void restore(XML_Node& dom, doublereal* soln) {} + + doublereal z(int jlocal) const { + return m_z[jlocal]; + } + doublereal zmin() const { return m_z[0]; } + doublereal zmax() const { return m_z[m_points - 1]; } + + + void setProfile(string name, doublereal* values, doublereal* soln) { + int n, j; + for (n = 0; n < m_nv; n++) { + if (name == componentName(n)) { + for (j = 0; j < m_points; j++) { + soln[index(n, j) + m_iloc] = values[j]; + } + return; + } + } + throw CanteraError("Domain1D::setProfile", + "unknown component: "+name); + } + + vector_fp& grid() { return m_z; } + const vector_fp& grid() const { return m_z; } + doublereal grid(int point) { return m_z[point]; } + + virtual void setupGrid(int n, const doublereal* z) {} + + /** + * Writes some or all initial solution values into array x, + * which is the solution vector for this domain. This allows + * initial values that have been set prior to installing this + * domain into the container to be written to the global + * solution vector. + */ + virtual void _getInitialSoln(doublereal* x) { + throw CanteraError("Domain1D::_getInitialSoln", + "base class method _getInitialSoln called!"); + } + + /** + * Perform any necessary domain-specific initialization using + * local solution vector x. + */ + virtual void _finalize(const doublereal* x) { + throw CanteraError("Domain1D::_finalize", + "base class method _finalize called!"); + } + + + protected: + + doublereal m_rdt; + int m_nv; + int m_points; + vector_fp m_slast; + doublereal m_time; + vector_fp m_max; + vector_fp m_min; + vector_fp m_rtol_ss, m_rtol_ts; + vector_fp m_atol_ss, m_atol_ts; + vector_fp m_z; + OneDim* m_container; + int m_index; + int m_type; + int m_iloc; + int m_jstart; + Domain1D *m_left, *m_right; + string m_id, m_desc; + Refiner* m_refiner; + + private: + + }; +} + +#endif + + diff --git a/Cantera/src/oneD/Inlet1D.h b/Cantera/src/oneD/Inlet1D.h index 6b53e2ee5..7c9854117 100644 --- a/Cantera/src/oneD/Inlet1D.h +++ b/Cantera/src/oneD/Inlet1D.h @@ -17,7 +17,7 @@ #ifndef CT_BDRY1D_H #define CT_BDRY1D_H -#include "Resid1D.h" +#include "Domain1D.h" #include "../SurfPhase.h" #include "../InterfaceKinetics.h" #include "StFlow.h" @@ -40,7 +40,7 @@ namespace Cantera { * The public methods are all virtual, and the base class * implementations throw exceptions. */ - class Bdry1D : public Resid1D { + class Bdry1D : public Domain1D { public: Bdry1D(); @@ -71,6 +71,10 @@ namespace Cantera { /// The total mass flow rate [kg/m2/s]. virtual doublereal mdot() {return m_mdot;} + virtual void _getInitialSoln(doublereal* x) { + cout << "Bdry1D::_getInitialSoln called! " << m_index << endl; + } + protected: void _init(int n); @@ -138,13 +142,32 @@ namespace Cantera { s << endl; } + virtual void showSolution(const doublereal* x) { + char buf[80]; + sprintf(buf, " Mass Flux: %10.4g kg/m^2/s \n", x[0]); + writelog(buf); + sprintf(buf, " Temperature: %10.4g K \n", x[1]); + writelog(buf); + if (m_flow) { + writelog(" Mass Fractions: \n"); + for (int k = 0; k < m_flow->phase().nSpecies(); k++) { + if (m_yin[k] != 0.0) { + sprintf(buf, " %16s %10.4g \n", + m_flow->phase().speciesName(k).c_str(), m_yin[k]); + writelog(buf); + } + } + } + writelog("\n"); + } + virtual void _getInitialSoln(doublereal* x) { x[0] = m_mdot; x[1] = m_temp; } virtual void _finalize(const doublereal* x) { - ; //m_mdot = x[0]; + ;//m_mdot = x[0]; //m_temp = x[1]; } @@ -156,7 +179,7 @@ namespace Cantera { virtual void eval(int jg, doublereal* xg, doublereal* rg, integer* diagg, doublereal rdt); virtual void save(XML_Node& o, doublereal* soln); - + virtual void restore(XML_Node& dom, doublereal* soln); protected: @@ -190,7 +213,38 @@ namespace Cantera { integer* diagg, doublereal rdt); virtual void save(XML_Node& o, doublereal* soln); + virtual void restore(XML_Node& dom, doublereal* soln); + virtual void _finalize(const doublereal* x) { + ; //m_temp = x[0]; + } + protected: + }; + + + /** + */ + class Outlet1D : public Bdry1D { + + public: + + Outlet1D() { + m_type = cOutletType; + } + virtual ~Outlet1D(){} + + virtual string componentName(int n) const; + + virtual void init(); + + virtual void eval(int jg, doublereal* xg, doublereal* rg, + integer* diagg, doublereal rdt); + + virtual void save(XML_Node& o, doublereal* soln); + virtual void restore(XML_Node& dom, doublereal* soln); + virtual void _finalize(const doublereal* x) { + ; //m_temp = x[0]; + } protected: }; @@ -219,6 +273,7 @@ namespace Cantera { integer* diagg, doublereal rdt); virtual void save(XML_Node& o, doublereal* soln); + virtual void restore(XML_Node& dom, doublereal* soln); virtual void _getInitialSoln(doublereal* x) { x[0] = m_temp; @@ -232,6 +287,14 @@ namespace Cantera { s << "------------------- Surface " << domainIndex() << " ------------------- " << endl; s << " temperature: " << m_temp << " K" << " " << x[0] << endl; } + + virtual void showSolution(const doublereal* x) { + char buf[80]; + sprintf(buf, " Temperature: %10.4g K \n", x[0]); + writelog(buf); + writelog("\n"); + } + protected: }; diff --git a/Cantera/src/oneD/Jac1D.h b/Cantera/src/oneD/Jac1D.h index 03dbb88a6..a68b22b2a 100644 --- a/Cantera/src/oneD/Jac1D.h +++ b/Cantera/src/oneD/Jac1D.h @@ -15,7 +15,7 @@ #ifndef CT_JAC1D_H #define CT_JAC1D_H -#include "Resid1D.h" +#include "Domain1D.h" #include "BandMatrix.h" //#include "ArrayViewer.h" #include "Array.h" @@ -25,7 +25,7 @@ namespace Cantera { /** * Class Jac1D evaluates the Jacobian of a system of equations - * defined by a residual function of class Resid1D. It is + * defined by a residual function of class Domain1D. It is * assumed that the Jacobian is banded. */ class Jac1D : public BandMatrix { @@ -36,7 +36,7 @@ namespace Cantera { * Constructor. The residual function defining the system of * equations must be supplied. */ - Jac1D(Resid1D& r); + Jac1D(Domain1D& r); /// Destructor. Does nothing. virtual ~Jac1D(){} @@ -69,7 +69,7 @@ namespace Cantera { protected: - Resid1D* m_resid; + Domain1D* m_resid; Array2D m_r1; // ArrayViewer m_x0, m_r0; int m_nv, m_points; diff --git a/Cantera/src/oneD/MultiNewton.cpp b/Cantera/src/oneD/MultiNewton.cpp index 7e89da014..637410393 100644 --- a/Cantera/src/oneD/MultiNewton.cpp +++ b/Cantera/src/oneD/MultiNewton.cpp @@ -39,9 +39,9 @@ namespace Cantera { // declarations for functions in newton_utils.h doublereal bound_step(const doublereal* x, - const doublereal* step, Resid1D& r, int loglevel=0); + const doublereal* step, Domain1D& r, int loglevel=0); doublereal norm_square(const doublereal* x, - const doublereal* step, Resid1D& r); + const doublereal* step, Domain1D& r); @@ -122,7 +122,7 @@ namespace Cantera { step[n] = -step[n]; } #ifdef DEBUG_STEP - Resid1D* d; + Domain1D* d; if (!ok) { for (n = 0; n < sz; n++) { d = r.pointDomain(n); @@ -176,9 +176,9 @@ namespace Cantera { writelog("\n\nDamped Newton iteration:\n"); writelog(dashedline); - sprintf(m_buf,"\n%s %9s %9s %9s %9s %9s %5s\n", + sprintf(m_buf,"\n%s %9s %9s %9s %9s %9s %5s %5s\n", "m","F_damp","F_bound","log10(ss)", - "log10(s0)","log10(s1)","N_jac"); + "log10(s0)","log10(s1)","N_jac","Age"); writelog(m_buf); writelog(dashedline+"\n"); } @@ -226,11 +226,11 @@ namespace Cantera { // write log information if (loglevel > 0) { doublereal ss = r.ssnorm(x1,step1); - sprintf(m_buf,"\n%d %9.5f %9.5f %9.5f %9.5f %9.5f %5d ", + sprintf(m_buf,"\n%d %9.5f %9.5f %9.5f %9.5f %9.5f %4d %d/%d", m,damp,fbound,log10(ss+SmallNumber), log10(s0+SmallNumber), log10(s1+SmallNumber), - jac.nEvals()); + jac.nEvals(), jac.age(), m_maxAge); writelog(m_buf); } @@ -284,6 +284,8 @@ namespace Cantera { // Check whether the Jacobian should be re-evaluated. if (jac.age() > m_maxAge) { + if (loglevel > 0) + writelog("\nMaximum Jacobian age reached ("+int2str(m_maxAge)+")\n"); forceNewJac = true; } @@ -302,13 +304,20 @@ namespace Cantera { // damp the Newton step m = dampStep(x, stp, x1, stp1, s1, r, jac, loglevel-1, frst); - frst = false; if (loglevel == 1 && m >= 0) { + if (frst) { + sprintf(m_buf,"\n\n %10s %10s %5s ", + "log10(ss)","log10(s1)","N_jac"); + writelog(m_buf); + sprintf(m_buf,"\n ------------------------------------"); + writelog(m_buf); + } doublereal ss = r.ssnorm(x, stp); - sprintf(m_buf,"\n %10.4f %10.4f %d ", + sprintf(m_buf,"\n %10.4f %10.4f %d ", log10(ss),log10(s1),jac.nEvals()); writelog(m_buf); } + frst = false; // Successful step, but not converged yet. Take the damped // step, and try again. @@ -323,7 +332,12 @@ namespace Cantera { // one was being used. If it was a new Jacobian, then // return -1 to signify failure. else if (m < 0) { - if (jac.age() > 1) forceNewJac = true; + if (jac.age() > 1) { + forceNewJac = true; + if (loglevel > 0) + writelog("Re-evaluating Jacobian, since no damping " + "coefficient\ncould be found with this Jacobian."); + } else goto done; } } diff --git a/Cantera/src/oneD/Newton1D.h b/Cantera/src/oneD/Newton1D.h index 2756a900c..9d030045a 100644 --- a/Cantera/src/oneD/Newton1D.h +++ b/Cantera/src/oneD/Newton1D.h @@ -28,21 +28,21 @@ namespace Cantera { doublereal norm(const doublereal* step); void step(doublereal* x, doublereal* step, - Resid1D& r, Jac1D& jac, int loglevel, int update=1); + Domain1D& r, Jac1D& jac, int loglevel, int update=1); doublereal boundStep(const doublereal* x0, const doublereal* step0, - const Resid1D& r, int loglevel); + const Domain1D& r, int loglevel); int dampStep(const doublereal* x0, const doublereal* step0, doublereal* x1, doublereal* step1, doublereal& s1, - Resid1D& r, Jac1D& jac, int loglevel, bool writetitle); - void getErrorWeights(const doublereal* x, doublereal* ewt, Resid1D& r); + Domain1D& r, Jac1D& jac, int loglevel, bool writetitle); + void getErrorWeights(const doublereal* x, doublereal* ewt, Domain1D& r); doublereal norm2(const doublereal* step, doublereal* ewt); doublereal norm_infty(const doublereal* step, doublereal* ewt); - int solve(doublereal* x0, doublereal* x1, Resid1D& r, Jac1D& jac, + int solve(doublereal* x0, doublereal* x1, Domain1D& r, Jac1D& jac, int loglevel); int timeIntegrate(int n, doublereal dt, doublereal* x0, doublereal* x1, - Resid1D& r, Jac1D& jac, int loglevel); - doublereal ssnorm(doublereal* x, doublereal* resid, Resid1D& r); + Domain1D& r, Jac1D& jac, int loglevel); + doublereal ssnorm(doublereal* x, doublereal* resid, Domain1D& r); void setOptions(int maxJacAge = 5, doublereal maxNormRatio = 0.001) { m_maxAge = maxJacAge; diff --git a/Cantera/src/oneD/OneDim.cpp b/Cantera/src/oneD/OneDim.cpp index 1a324a129..97ba02ae6 100644 --- a/Cantera/src/oneD/OneDim.cpp +++ b/Cantera/src/oneD/OneDim.cpp @@ -34,7 +34,7 @@ namespace Cantera { * Construct a OneDim container for the domains pointed at by the * input vector of pointers. */ - OneDim::OneDim(vector domains) : + OneDim::OneDim(vector domains) : m_tmin(1.0e-16), m_tmax(0.1), m_tfactor(0.5), m_jac(0), m_newt(0), m_rdt(0.0), m_jac_ok(false), @@ -56,10 +56,18 @@ namespace Cantera { } + int OneDim::domainIndex(string name) { + for (int n = 0; n < m_nd; n++) { + if (domain(n).id() == name) return n; + } + throw CanteraError("OneDim::domainIndex","no domain named >>"+name+"<<"); + } + + /** * Domains are added left-to-right. */ - void OneDim::addDomain(Resid1D* d) { + void OneDim::addDomain(Domain1D* d) { // if 'd' is not the first domain, link it to the last domain // added (the rightmost one) @@ -103,6 +111,7 @@ namespace Cantera { } } + /** * Save statistics on function and Jacobiab evaulation, and reset * the counters. Statistics are saved only if the number of @@ -118,7 +127,7 @@ namespace Cantera { void OneDim::saveStats() { if (m_jac) { int nev = m_jac->nEvals(); - if (nev > 0) { + if (nev > 0 && m_nevals > 0) { m_gridpts.push_back(m_pts); m_jacEvals.push_back(m_jac->nEvals()); m_jacElapsed.push_back(m_jac->elapsedTime()); @@ -144,7 +153,7 @@ namespace Cantera { saveStats(); m_pts = 0; for (i = 0; i < m_nd; i++) { - Resid1D* d = m_dom[i]; + Domain1D* d = m_dom[i]; int np = d->nPoints(); int nv = d->nComponents(); @@ -207,8 +216,8 @@ namespace Cantera { * 8/26/02 changed '<' to '<=' DGG * */ - Resid1D* OneDim::pointDomain(int i) { - Resid1D* d = right(); + Domain1D* OneDim::pointDomain(int i) { + Domain1D* d = right(); while (d) { if (d->loc() <= i) return d; d = d->left(); @@ -227,7 +236,7 @@ namespace Cantera { fill(m_mask.begin(), m_mask.end(), 0); if (rdt < 0.0) rdt = m_rdt; - vector::iterator d; + vector::iterator d; // iterate over the bulk domains first for (d = m_bulk.begin(); d != m_bulk.end(); ++d) @@ -242,6 +251,8 @@ namespace Cantera { clock_t t1 = clock(); m_evaltime += double(t1 - t0)/CLOCKS_PER_SEC; m_nevals++; + //string ne = string("evals = ")+int2str(m_nevals)+"\n"; + //writelog(ne.c_str()); } } @@ -251,7 +262,7 @@ namespace Cantera { * residual. Used only for diagnostic output. */ doublereal OneDim::ssnorm(doublereal* x, doublereal* r) { - eval(-1, x, r, 0.0); + eval(-1, x, r, 0.0, 0); doublereal ss = 0.0; for (int i = 0; i < m_size; i++) { ss = fmaxx(fabs(r[i]),ss); @@ -275,7 +286,7 @@ namespace Cantera { // iterate over all domains, preparing each one to begin // time stepping - Resid1D* d = left(); + Domain1D* d = left(); while (d) { d->initTimeInteg(dt, x); d = d->right(); @@ -300,7 +311,7 @@ namespace Cantera { */ void OneDim::init() { if (!m_init) { - Resid1D* d = left(); + Domain1D* d = left(); while (d) { d->init(); d = d->right(); @@ -313,9 +324,11 @@ namespace Cantera { /** * Signal that the current Jacobian is no longer valid. */ - void Resid1D::needJacUpdate() { - if (m_container) + void Domain1D::needJacUpdate() { + if (m_container) { m_container->jacobian().setAge(10000); + m_container->saveStats(); + } } @@ -362,8 +375,8 @@ namespace Cantera { copy(r, r + m_size, x); if (m == 100) { dt *= 1.5; - cout << "m = 100, dt = " << dt << endl; } + // else dt /= 1.5; if (dt > m_tmax) dt = m_tmax; } @@ -372,7 +385,6 @@ namespace Cantera { else { if (loglevel > 0) writelog("...failure.\n"); dt *= m_tfactor; - cout << "halved dt = " << dt << endl; if (dt < m_tmin) throw CanteraError("OneDim::timeStep", "Time integration failed."); @@ -421,7 +433,7 @@ namespace Cantera { addString(sim,"timestamp",asctime(newtime)); if (desc != "") addString(sim,"description",desc); - Resid1D* d = left(); + Domain1D* d = left(); while (d) { d->save(sim, sol); d = d->right(); @@ -431,6 +443,6 @@ namespace Cantera { throw CanteraError("save","could not open file "+fname); ct->write(s); s.close(); - writelog("Solution saved to file "+fname+" as solution '"+id+"'.\n"); + writelog("Solution saved to file "+fname+" as solution "+id+".\n"); } } diff --git a/Cantera/src/oneD/OneDim.h b/Cantera/src/oneD/OneDim.h index 7b2ff0849..6bff3256f 100644 --- a/Cantera/src/oneD/OneDim.h +++ b/Cantera/src/oneD/OneDim.h @@ -5,7 +5,7 @@ #ifndef CT_ONEDIM_H #define CT_ONEDIM_H -#include "Resid1D.h" +#include "Domain1D.h" namespace Cantera { @@ -14,7 +14,7 @@ namespace Cantera { /** * Container class for multiple-domain 1D problems. Each domain is - * represented by an instance of Resid1D. + * represented by an instance of Domain1D. */ class OneDim { @@ -24,13 +24,13 @@ namespace Cantera { OneDim(); // Constructor. - OneDim(vector domains); + OneDim(vector domains); /// Destructor. virtual ~OneDim(); /// Add a domain. - void addDomain(Resid1D* d); + void addDomain(Domain1D* d); /// Return a reference to the Jacobian evaluator. MultiJac& jacobian(); @@ -50,7 +50,9 @@ namespace Cantera { int nDomains() const { return m_nd; } /// Return a reference to domain i. - Resid1D& domain(int i) const { return *m_dom[i]; } + Domain1D& domain(int i) const { return *m_dom[i]; } + + int domainIndex(string name); /// The index of the start of domain i in the solution vector. int start(int i) const { return m_dom[i]->loc(); } @@ -59,10 +61,10 @@ namespace Cantera { int size() const { return m_size; } /// Pointer to left-most domain (first added). - Resid1D* left() { return m_dom[0]; } + Domain1D* left() { return m_dom[0]; } /// Pointer to right-most domain (last added). - Resid1D* right() { return m_dom.back(); } + Domain1D* right() { return m_dom.back(); } /// Number of solution components at global point jg. int nVars(int jg) { return m_nvars[jg]; } @@ -123,7 +125,7 @@ namespace Cantera { int count = 1); /// Pointer to the domain global point i belongs to. - Resid1D* pointDomain(int i); + Domain1D* pointDomain(int i); void resize(); @@ -136,7 +138,6 @@ namespace Cantera { double* r, int loglevel); void writeStats(); - void saveStats(); void save(string fname, string id, string desc, doublereal* sol); @@ -151,6 +152,7 @@ namespace Cantera { else m_ts_jac_age = m_ss_jac_age; } + void saveStats(); protected: @@ -167,7 +169,7 @@ namespace Cantera { int m_bw; // Jacobian bandwidth int m_size; // solution vector size - vector m_dom, m_connect, m_bulk; + vector m_dom, m_connect, m_bulk; bool m_init; vector_int m_nvars; @@ -179,6 +181,8 @@ namespace Cantera { // options int m_ss_jac_age, m_ts_jac_age; + private: + // statistics int m_nevals; doublereal m_evaltime; @@ -188,7 +192,6 @@ namespace Cantera { vector_int m_funcEvals; vector_fp m_funcElapsed; - private: }; diff --git a/Cantera/src/oneD/Sim1D.cpp b/Cantera/src/oneD/Sim1D.cpp index e7ffc7fbb..83551e37e 100644 --- a/Cantera/src/oneD/Sim1D.cpp +++ b/Cantera/src/oneD/Sim1D.cpp @@ -7,7 +7,15 @@ namespace Cantera { - Sim1D::Sim1D(vector& domains) : OneDim(domains) { + static void drawline() { + string s(78,'.'); + s += '\n'; + writelog(s.c_str()); + } + + Sim1D::Sim1D() : OneDim() {} + + Sim1D::Sim1D(vector& domains) : OneDim(domains) { // resize the internal solution vector and the wprk array, // and perform domain-specific initialization of the @@ -54,6 +62,11 @@ namespace Cantera { return m_x[iloc]; } + doublereal Sim1D::workValue(int dom, int comp, int localPoint) const { + int iloc = domain(dom).loc() + domain(dom).index(comp, localPoint); + return m_xnew[iloc]; + } + /** * @param dom domain number, beginning with 0 for the leftmost domain. @@ -72,7 +85,7 @@ namespace Cantera { void Sim1D::setProfile(int dom, int comp, const vector_fp& pos, const vector_fp& values) { - Resid1D& d = domain(dom); + Domain1D& d = domain(dom); int np = d.nPoints(); int n; doublereal z0 = d.zmin(); @@ -87,6 +100,57 @@ namespace Cantera { } + void Sim1D::save(string fname, string id, string desc) { + OneDim::save(fname, id, desc, m_x.begin()); + } + + /** + * Initialize the solution with a previously-saved solution. + */ + void Sim1D::restore(string fname, string id) { + ifstream s(fname.c_str()); + char buf[100]; + if (!s) + throw CanteraError("Sim1D::restore", + "could not open input file "+fname); + + XML_Node root; + root.build(s); + s.close(); + + XML_Node* f = root.findID(id); + if (!f) { + throw CanteraError("Sim1D::restore","No solution with id = "+id); + } + + vector xd; + int sz = 0, np, nv, m; + for (m = 0; m < m_nd; m++) { + XML_Node* d = f->findID(domain(m).id()); + if (!d) { + writelog("No data for domain "+domain(m).id()); + xd.push_back(0); + sz += domain(m).nComponents(); + } + else { + XML_Node& node = *d; + xd.push_back(d); + np = intValue(node["points"]); + nv = intValue(node["components"]); + sz += np*domain(m).nComponents(); + } + } + m_x.resize(sz); + m_xnew.resize(sz); + for (m = 0; m < m_nd; m++) { + if (xd[m]) { + domain(m).restore(*xd[m], m_x.begin() + domain(m).loc()); + } + } + resize(); + finalize(); + } + void Sim1D::setFlatProfile(int dom, int comp, doublereal v) { int np = domain(dom).nPoints(); @@ -101,6 +165,14 @@ namespace Cantera { } } + void Sim1D::showSolution() { + for (int n = 0; n < m_nd; n++) { + writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id() + +" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n"); + domain(n).showSolution(m_x.begin() + start(n)); + } + } + void Sim1D::finalize() { for (int n = 0; n < m_nd; n++) { @@ -123,7 +195,7 @@ namespace Cantera { else if (m > -10) throw CanteraError("Sim1D::newtonSolve","no solution found"); else { - cout << "ERROR: solve returned m = " << m << endl; + writelog(string("ERROR: solve returned m = ") + int2str(m) + "\n"); exit(-1); } } @@ -148,13 +220,13 @@ namespace Cantera { try { if (loglevel > 0) { - writelog("Attempt Newton solution of steady-state problem..."); + drawline(); + writelog("\nAttempt Newton solution of steady-state problem..."); } newtonSolve(loglevel-1); if (loglevel > 0) { - writelog("success.\n\n"); - //writelog("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n"); + writelog(" success.\n\n"); writelog("Problem solved on ["); for (int mm = 1; mm < nDomains(); mm+=2) { writelog(int2str(domain(mm).nPoints())); @@ -162,7 +234,6 @@ namespace Cantera { } writelog("]"); writelog(" point grid(s).\n\n"); - //writelog("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"); } ok = true; soln_number++; @@ -172,9 +243,16 @@ namespace Cantera { catch (CanteraError) { char buf[100]; - if (loglevel > 0) writelog("failure. \n\n"); - if (loglevel == 1) writelog("Take "+int2str(nsteps)+" timesteps "); - dt = timeStep(nsteps, dt, m_x.begin(), m_xnew.begin(), loglevel-1); + if (loglevel > 0) { + writelog(" failure. \n\n"); + drawline(); + // } + //if (loglevel == 1) + writelog("Take "+int2str(nsteps)+ + " timesteps "); + } + dt = timeStep(nsteps, dt, m_x.begin(), m_xnew.begin(), + loglevel-1); if (loglevel == 1) { sprintf(buf, " %10.4g %10.4g \n", dt, log10(ssnorm(m_x.begin(), m_xnew.begin()))); @@ -183,8 +261,7 @@ namespace Cantera { istep++; if (istep >= int(m_steps.size())) { nsteps = m_steps.back(); - dt *= 2.0; - cout << " doubled dt = " << dt << endl; + // dt *= 2.0; } else { nsteps = m_steps[istep]; @@ -192,7 +269,7 @@ namespace Cantera { if (dt > m_tmax) dt = m_tmax; } } - if (loglevel > 2) showSolution(cout); + if (loglevel > 2) showSolution(); if (refine_grid) { new_points = refine(loglevel); @@ -203,6 +280,7 @@ namespace Cantera { } } + /** * Refine the grid in all domains. */ @@ -211,10 +289,11 @@ namespace Cantera { vector_fp znew, xnew; doublereal xmid, zmid; int strt, n, m, i; + vector_int dsize; for (n = 0; n < m_nd; n++) { strt = znew.size(); - Resid1D& d = domain(n); + Domain1D& d = domain(n); Refiner& r = d.refiner(); // determine where new points are needed @@ -226,45 +305,54 @@ namespace Cantera { // loop over points in the current grid int npnow = d.nPoints(); + int nstart = znew.size(); for (m = 0; m < npnow; m++) { - // add the current grid point to the new grid - znew.push_back(d.grid(m)); - - // do the same for the solution at this point - for (i = 0; i < comp; i++) { - xnew.push_back(value(n, i, m)); - } - - // now check whether a new point is needed in the interval to the - // right of point m, and if so, add entries to znew and xnew for - // this new point - - if (r.newPointNeeded(m)) { - - // add new point at midpoint - zmid = 0.5*(d.grid(m) + d.grid(m+1)); - znew.push_back(zmid); - - // for each component, linearly interpolate the solution to - // this point + if (r.keepPoint(m)) { + // add the current grid point to the new grid + znew.push_back(d.grid(m)); + + // do the same for the solution at this point for (i = 0; i < comp; i++) { - xmid = 0.5*(value(n, i, m) + value(n, i, m+1)); - xnew.push_back(xmid); + xnew.push_back(value(n, i, m)); + } + + // now check whether a new point is needed in the + // interval to the right of point m, and if so, add + // entries to znew and xnew for this new point + + if (r.newPointNeeded(m) && m < npnow - 1) { + + // add new point at midpoint + zmid = 0.5*(d.grid(m) + d.grid(m+1)); + znew.push_back(zmid); + np++; + + // for each component, linearly interpolate + // the solution to this point + for (i = 0; i < comp; i++) { + xmid = 0.5*(value(n, i, m) + value(n, i, m+1)); + xnew.push_back(xmid); + } } } + //else { + // throw CanteraError("refine","keepPoint is false at m = "+int2str(m)); + //} } + dsize.push_back(znew.size() - nstart); } - // At this point, the new grid znew and the new solution vector xnew have - // been constructed, but the domains themselves have not yet been modified. - // Now update each domain with the new grid. + // At this point, the new grid znew and the new solution + // vector xnew have been constructed, but the domains + // themselves have not yet been modified. Now update each + // domain with the new grid. int gridstart = 0, gridsize; for (n = 0; n < m_nd; n++) { - Resid1D& d = domain(n); + Domain1D& d = domain(n); Refiner& r = d.refiner(); - gridsize = d.nPoints() + r.nNewPoints(); + gridsize = dsize[n]; // d.nPoints() + r.nNewPoints(); d.setupGrid(gridsize, znew.begin() + gridstart); gridstart += gridsize; } diff --git a/Cantera/src/oneD/Sim1D.h b/Cantera/src/oneD/Sim1D.h index 7760b6fe1..a1981cc2e 100644 --- a/Cantera/src/oneD/Sim1D.h +++ b/Cantera/src/oneD/Sim1D.h @@ -11,16 +11,19 @@ namespace Cantera { /** - * One-dimensional simulations. + * One-dimensional simulations. Class Sim1D extends class OneDim + * by storing the solution vector, and by adding a hybrid + * Newton/time-stepping solver. */ class Sim1D : public OneDim { public: /** - * Default constructor. This constructor can be used to create - * a dummy object if necessary, but is not usually called in - * user programs. Use the next constructor instead. + * Default constructor. This constructor is provided to make + * the class default-constructible, but is not meant to be + * used in most applications. Use the next constructor + * instead. */ Sim1D(); @@ -28,20 +31,32 @@ namespace Cantera { /** * Standard constructor. * @param domains A vector of pointers to the domains to be linked together. - * The domains must appear in left-to-right order. + * The domain pointers must be entered in left-to-right order --- i.e., + * the pointer to the leftmost domain is domain[0], the pointer to the + * domain to its right is domain[1], etc. */ - Sim1D(vector& domains); + Sim1D(vector& domains); /// Destructor. Does nothing. virtual ~Sim1D(){} + /** + * @name Setting initial values + * + * These methods are used to set the initial values of + * solution components. + */ + //@{ + /// Set one entry in the solution vector. void setValue(int dom, int comp, int localPoint, doublereal value); /// Get one entry in the solution vector. doublereal value(int dom, int comp, int localPoint) const; + doublereal workValue(int dom, int comp, int localPoint) const; + /// Specify a profile for one component of one domain. void setProfile(int dom, int comp, const vector_fp& pos, const vector_fp& values); @@ -49,11 +64,13 @@ namespace Cantera { /// Set component 'comp' of domain 'dom' to value 'v' at all points. void setFlatProfile(int dom, int comp, doublereal v); + //@} + + void save(string fname, string id, string desc); + /// Print to stream s the current solution for all domains. void showSolution(ostream& s); - - /// Calls method _finalize in each domain. - void finalize(); + void showSolution(); void setTimeStep(doublereal stepsize, int n, integer* tsteps); @@ -61,29 +78,40 @@ namespace Cantera { void solve(int loglevel = 0, bool refine_grid = true); + void eval(doublereal rdt=-1.0, int count = 1) { + OneDim::eval(-1, m_x.begin(), m_xnew.begin(), rdt, count); + } + + /// Refine the grid in all domains. int refine(int loglevel=0); + /// Set the criteria for grid refinement. void setRefineCriteria(int dom = -1, doublereal ratio = 10.0, doublereal slope = 0.8, doublereal curve = 0.8); + void restore(string fname, string id); + protected: - vector_fp m_x; // the solution vector + vector_fp m_x; // the solution vector vector_fp m_xnew; // a work array used to hold the residual - // or the new solution + // or the new solution doublereal m_tstep; // timestep vector_int m_steps; // array of number of steps to take before - // re-attempting the steady-state solution + // re-attempting the steady-state solution + private: + /// Calls method _finalize in each domain. + void finalize(); + void newtonSolve(int loglevel); }; } - #endif diff --git a/Cantera/src/oneD/StFlow.cpp b/Cantera/src/oneD/StFlow.cpp index 75d05d435..ffe227002 100644 --- a/Cantera/src/oneD/StFlow.cpp +++ b/Cantera/src/oneD/StFlow.cpp @@ -93,16 +93,21 @@ namespace Cantera { //---------------------- drawline ---------------------------------- - inline void drawline(ostream& s) { + static void drawline(ostream& s) { s << "\n-------------------------------------" << "------------------------------------------"; } + static void drawline() { + writelog("\n-------------------------------------" + "------------------------------------------"); + } + //--------------------- linear interp ------------------------------ StFlow::StFlow(igthermo_t* ph, int nsp, int points) : - Resid1D(nsp+4, points), + Domain1D(nsp+4, points), m_inlet_u(0.0), m_inlet_V(0.0), m_inlet_T(-1.0), @@ -130,7 +135,7 @@ namespace Cantera { int nsp2 = m_thermo->nSpecies(); if (nsp2 != m_nsp) { m_nsp = nsp2; - Resid1D::resize(m_nsp+4, points); + Domain1D::resize(m_nsp+4, points); } @@ -148,7 +153,7 @@ namespace Cantera { // but turn off the energy equation at all points m_do_energy.resize(m_points,false); - m_diff.resize(m_nsp,m_points); + m_diff.resize(m_nsp*m_points); m_flux.resize(m_nsp,m_points); m_wdot.resize(m_nsp,m_points, 0.0); m_surfdot.resize(m_nsp, 0.0); @@ -163,8 +168,8 @@ namespace Cantera { vmin[0] = -1.e20; vmax[0] = 1.e20; - // no negative V - vmin[1] = -0.1; + // V + vmin[1] = -1.e20; vmax[1] = 1.e20; // temperature bounds @@ -173,13 +178,13 @@ namespace Cantera { // lamda should be negative vmin[3] = -1.e20; - vmax[3] = 1.0; + vmax[3] = 1.e20; // mass fraction bounds int k; for (k = 0; k < m_nsp; k++) { vmin[4+k] = -1.0e-5; - vmax[4+k] = 1.1; + vmax[4+k] = 1.0e5; } setBounds(vmin.size(), vmin.begin(), vmax.size(), vmax.begin()); @@ -187,13 +192,19 @@ namespace Cantera { //-------------------- default error tolerances ---------------- vector_fp rtol(m_nv, 1.0e-8); vector_fp atol(m_nv, 1.0e-15); - setTolerances(rtol.size(), rtol.begin(), atol.size(), atol.begin()); + setTolerances(rtol.size(), rtol.begin(), atol.size(), atol.begin(),false); + setTolerances(rtol.size(), rtol.begin(), atol.size(), atol.begin(),true); //-------------------- grid refinement ------------------------- m_refiner->setActive(0, false); m_refiner->setActive(1, false); m_refiner->setActive(2, false); m_refiner->setActive(3, false); + + vector_fp gr; + for (int ng = 0; ng < m_points; ng++) gr.push_back(1.0*ng/m_points); + setupGrid(m_points, gr.begin()); + setID("stagnation flow"); } @@ -201,7 +212,7 @@ namespace Cantera { * Change the grid size. Called after grid refinement. */ void StFlow::resize(int points) { - Resid1D::resize(m_nv, points); + Domain1D::resize(m_nv, points); m_rho.resize(m_points, 0.0); m_wtm.resize(m_points, 0.0); @@ -210,7 +221,12 @@ namespace Cantera { m_visc.resize(m_points, 0.0); m_tcon.resize(m_points, 0.0); - m_diff.resize(m_nsp,m_points); + if (m_transport_option == c_Mixav_Transport) { + m_diff.resize(m_nsp*m_points); + } + else { + m_diff.resize(m_nsp*m_nsp*m_points); + } m_flux.resize(m_nsp,m_points); m_wdot.resize(m_nsp,m_points, 0.0); m_do_energy.resize(m_points,false); @@ -338,11 +354,13 @@ namespace Cantera { // thermodynamic properties only if a Jacobian is // not being evaluated - if (jpt < 0) updateThermo(x, j0, j1); + if (jpt < 0) + updateThermo(x, j0, j1); // update transport properties only if a Jacobian is // not being evaluated - if (jpt < 0) updateTransport(x, j0, j1); + if (jpt < 0) + updateTransport(x, j0, j1); // update the species diffusive mass fluxes whether or not a // Jacobian is being evaluated @@ -367,8 +385,6 @@ namespace Cantera { // these may be modified by a boundary object -#define NEW_INLET -#ifdef NEW_INLET // Continuity. This propagates information right-to-left, // since rho_u at point 0 is dependent on rho_u at point 1, @@ -385,8 +401,6 @@ namespace Cantera { rsd[index(c_offset_V,0)] = V(x,0); rsd[index(c_offset_T,0)] = T(x,0); rsd[index(c_offset_L,0)] = -rho_u(x,0); - //cout << "density = " << density(0) << " " << u(x,0) - // << " " << rho_u(x,0) << endl; // The default boundary condition for species is zero // flux. However, the boundary object may modify @@ -395,28 +409,6 @@ namespace Cantera { rsd[index(c_offset_Y + k, 0)] = -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); } -#else - // first, call the left boundary object to evaluate - // the residual - - m_boundary[0]->eval(x + index(0,0), m_rho[0], m_flux.begin(), - rsd + index(0,0)); - - - // Now modify the left boundary conditions to allow - // specifying the mass flux at both boundaries. The - // right mass flux is specified directly as a boundary - // condition on the continuity equation; the left mass - // flux is matched by adjusting lambda. - - // Shift the left continuity boundary condition to lambda, - rsd[index(c_offset_L, 0)] = rsd[index(c_offset_U, 0)]; - - // and replace it with the continuity equation. - rsd[index(c_offset_U,0)] = - -(rho_u(x,1) - rho_u(x,0))/m_dz[0] - -(density(1)*V(x,1) + density(0)*V(x,0)); -#endif } @@ -429,22 +421,18 @@ namespace Cantera { else if (j == m_points - 1) { // the boundary object connected to the right of this - // one may modify these equations by subtracting its - // values for V, T, and mdot. As a result, these - // residual equations will force the solution - // variables to the values for the boundary object + // one may modify or replace these equations. The + // default boundary conditions are zero u, V, and T, + // and zero diffusive flux for all species. + rsd[index(0,j)] = rho_u(x,j); rsd[index(1,j)] = V(x,j); rsd[index(2,j)] = T(x,j); - doublereal sum = 0.0; for (k = 0; k < m_nsp; k++) { sum += Y(x,k,j); - rsd[index(k+4,j)] = rho_u(x,j)*Y(x,k,j) + m_flux(k,j-1); + rsd[index(k+4,j)] = m_flux(k,j-1); } - - // TODO: why is this done here, but not for the left - // boundary or interior? rsd[index(4,j)] = 1.0 - sum; diag[index(4,j)] = 0; } @@ -481,7 +469,9 @@ namespace Cantera { //------------------------------------------------- rsd[index(c_offset_V,j)] = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j) - - m_rho[j]*V(x,j)*V(x,j))/m_rho[j]; + - m_rho[j]*V(x,j)*V(x,j))/m_rho[j] + - rdt*(V(x,j) - V_prev(j)); + diag[index(c_offset_V, j)] = 1; //------------------------------------------------- @@ -494,16 +484,16 @@ namespace Cantera { doublereal convec, diffus; for (k = 0; k < m_nsp; k++) { - if (m_do_species[k]) { - convec = rho_u(x,j)*dYdz(x,k,j); - diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) - /(z(j+1) - z(j-1)); - rsd[index(c_offset_Y + k, j)] - = (m_wt[k]*(wdot(k,j) ) - - convec - diffus)/m_rho[j] - - rdt*(Y(x,k,j) - Y_prev(k,j)); - diag[index(c_offset_Y + k, j)] = 1; - } + //if (m_do_species[k]) { + convec = rho_u(x,j)*dYdz(x,k,j); + diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) + /(z(j+1) - z(j-1)); + rsd[index(c_offset_Y + k, j)] + = (m_wt[k]*(wdot(k,j) ) + - convec - diffus)/m_rho[j] + - rdt*(Y(x,k,j) - Y_prev(k,j)); + diag[index(c_offset_Y + k, j)] = 1; + //} } @@ -540,8 +530,6 @@ namespace Cantera { rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); diag[index(c_offset_T, j)] = 1; - - } } @@ -554,12 +542,13 @@ namespace Cantera { diag[index(c_offset_Y+k, j)] = 0; } } + if (!m_do_energy[j]) { rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); diag[index(c_offset_T, j)] = 0; } - // lambda + // Propagate lambda from left to right if (j > 0) { rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); diag[index(c_offset_L, j)] = 0; @@ -575,12 +564,22 @@ namespace Cantera { */ void AxiStagnFlow::updateTransport(doublereal* x,int j0, int j1) { int j; - //for (j = j0; j <= j1; j++) { - for (j = j0; j < j1; j++) { - setGasAtMidpoint(x,j); - m_visc[j] = m_trans->viscosity(); - m_trans->getMixDiffCoeffs(&m_diff(0,j)); - m_tcon[j] = m_trans->thermalConductivity(); + if (m_transport_option == c_Mixav_Transport) { + for (j = j0; j < j1; j++) { + setGasAtMidpoint(x,j); + m_visc[j] = m_trans->viscosity(); + m_trans->getMixDiffCoeffs(m_diff.begin() + j*m_nsp); + m_tcon[j] = m_trans->thermalConductivity(); + } + } + else if (m_transport_option == c_Multi_Transport) { + for (j = j0; j < j1; j++) { + setGasAtMidpoint(x,j); + m_visc[j] = m_trans->viscosity(); + m_trans->getMultiDiffCoeffs(m_nsp, + m_diff.begin() + mindex(0,0,j)); + m_tcon[j] = m_trans->thermalConductivity(); + } } } @@ -762,10 +761,19 @@ namespace Cantera { */ void OneDFlow::updateTransport(doublereal* x,int j0, int j1) { int j; - for (j = j0; j < j1; j++) { - setGasAtMidpoint(x,j); - m_trans->getMixDiffCoeffs(&m_diff(0,j)); - m_tcon[j] = m_trans->thermalConductivity(); + if (m_transport_option == c_Mixav_Transport) { + for (j = j0; j < j1; j++) { + setGasAtMidpoint(x,j); + m_trans->getMixDiffCoeffs(m_diff.begin() + j*m_nsp); + m_tcon[j] = m_trans->thermalConductivity(); + } + } + else if (m_transport_option == c_Multi_Transport) { + for (j = j0; j < j1; j++) { + setGasAtMidpoint(x,j); + m_trans->getMultiDiffCoeffs(m_nsp, m_diff.begin() + mindex(0,0,j)); + m_tcon[j] = m_trans->thermalConductivity(); + } } } @@ -821,12 +829,66 @@ namespace Cantera { } + /** + * Print the solution. + */ + void StFlow::showSolution(const doublereal* x) { + int nn = m_nv/5; + int i, j, n; + char* buf = new char[100]; + + // The mean molecular weight is needed to convert + updateThermo(x, 0, m_points-1); + + sprintf(buf, " Pressure: %10.4g Pa \n", m_press); + writelog(buf); + for (i = 0; i < nn; i++) { + drawline(); + sprintf(buf, "\n z "); + writelog(buf); + for (n = 0; n < 5; n++) { + sprintf(buf, " %10s ",componentName(i*5 + n).c_str()); + writelog(buf); + } + drawline(); + for (j = 0; j < m_points; j++) { + sprintf(buf, "\n %10.4g ",m_z[j]); + writelog(buf); + for (n = 0; n < 5; n++) { + sprintf(buf, " %10.4g ",component(x, i*5+n,j)); + writelog(buf); + } + } + writelog("\n"); + } + int nrem = m_nv - 5*nn; + drawline(); + sprintf(buf, "\n z "); + writelog(buf); + for (n = 0; n < nrem; n++) { + sprintf(buf, " %10s ", componentName(nn*5 + n).c_str()); + writelog(buf); + } + drawline(); + for (j = 0; j < m_points; j++) { + sprintf(buf, "\n %10.4g ",m_z[j]); + writelog(buf); + for (n = 0; n < nrem; n++) { + sprintf(buf, " %10.4g ",component(x, nn*5+n,j)); + writelog(buf); + } + } + writelog("\n"); + } + + /** * Update the diffusive mass fluxes. */ void StFlow::updateDiffFluxes(const doublereal* x, int j0, int j1) { - int j, k; - double sum, wtm, rho, dz; + int j, k, m; + doublereal sum, wtm, rho, dz, gradlogT, s; + switch (m_transport_option) { case c_Mixav_Transport: @@ -837,20 +899,44 @@ namespace Cantera { dz = z(j+1) - z(j); for (k = 0; k < m_nsp; k++) { - m_flux(k,j) = m_wt[k]*(rho*m_diff(k,j)/wtm); + m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm); m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz; sum -= m_flux(k,j); } - // correction flux to insure that \sum_k Y_k j_k = 0. + // correction flux to insure that \sum_k Y_k V_k = 0. for (k = 0; k < m_nsp; k++) m_flux(k,j) += sum*Y(x,k,j); } break; case c_Multi_Transport: - cout << " not yet implemented... " << endl; + for (m = j0; m < j1; m++) { + wtm = m_wtm[m]; + rho = density(m); + dz = z(m+1) - z(m); + + for (k = 0; k < m_nsp; k++) { + sum = 0.0; + for (j = 0; j < m_nsp; j++) { + s = m_wt[j]*m_diff[mindex(k,j,m)]; + s *= (X(x,k,m+1) - X(x,k,m))/dz; + sum += s; + } + m_flux(k,m) = sum*rho*m_wt[k]/(wtm*wtm); + } + } + break; + default: + throw CanteraError("updateDiffFluxes","unknown transport model"); } + if (m_do_soret) { - cout << " net yet implemented... " << endl; + throw CanteraError("updateDiffFluxes","not yet"); + for (m = j0; m < j1; m++) { + gradlogT = 2.0*(T(x,m+1) - T(x,m))/(T(x,m+1) + T(x,m)); + for (k = 0; k < m_nsp; k++) { + m_flux(k,m) -= m_dthermal(k,m)*gradlogT; + } + } } } @@ -1015,7 +1101,6 @@ namespace Cantera { "solution contains no grid points."); } - cout << "importing...." << endl; writelog("Importing datasets:\n"); for (n = 0; n < nd; n++) { XML_Node& fa = *d[n]; @@ -1025,13 +1110,11 @@ namespace Cantera { writelog("axial velocity "); if ((int) x.size() == np) { for (j = 0; j < np; j++) { - cout << j << " " << x[j] << " " << np << endl; - cout << index(0,j) << " " << size_soln << endl; soln[index(0,j)] = x[j]; } } else { - cout << "error..." << endl; + //cout << "error..." << endl; goto error; } } @@ -1101,96 +1184,152 @@ namespace Cantera { } - void StFlow::save(string fname, string id, string desc, doublereal* sol) { - int k; - struct tm *newtime; - time_t aclock; - ::time( &aclock ); /* Get time in seconds */ - newtime = localtime( &aclock ); /* Convert time to struct tm form */ + void StFlow::restore(XML_Node& dom, doublereal* soln) { - ArrayViewer soln(m_nv, m_points, sol); + vector ignored; + int nsp = m_thermo->nSpecies(); + vector_int did_species(nsp, 0); - XML_Node root("doc"); - ifstream fin(fname.c_str()); - XML_Node* ct; - if (fin) { - root.build(fin); - XML_Node* same_ID = root.findID(id); - int jid = 1; - string idnew = id; - while (same_ID != 0) { - idnew = id + "_" + int2str(jid); - jid++; - same_ID = root.findID(idnew); + vector str; + dom.getChildren("string",str); + int nstr = str.size(); + for (int istr = 0; istr < nstr; istr++) { + XML_Node& nd = *str[istr]; + writelog(nd["title"]+": "+nd.value()+"\n"); + } + + map params; + getFloats(dom, params); + setPressure(params["pressure"]); + + + vector d; + dom.child("grid_data").getChildren("floatArray",d); + int nd = d.size(); + + vector_fp x; + int n, np, j, ks, k; + string nm; + bool readgrid = false, wrote_header = false; + for (n = 0; n < nd; n++) { + XML_Node& fa = *d[n]; + nm = fa["title"]; + if (nm == "z") { + getFloatArray(fa,x,false); + np = x.size(); + writelog("Grid contains "+int2str(np)+ + " points.\n"); + readgrid = true; + + // note that setupGrid also resizes the domain. + setupGrid(np, x.begin()); } - id = idnew; - fin.close(); - ct = &root.child("ctml"); } - else { - ct = &root.addChild("ctml"); + if (!readgrid) { + throw CanteraError("StFlow::restore", + "domain contains no grid points."); } - XML_Node& flow = (XML_Node&)ct->addChild("flowfield"); - flow.addAttribute("type",flowType()); - flow.addAttribute("id",id); - addString(flow,"timestamp",asctime(newtime)); - addFloat(flow, "pressure", m_press, "Pa", "pressure"); - // addString(flow,"solve_time",fp2str(m_container->solveTime())); - if (desc != "") addString(flow,"description",desc); - XML_Node& gv = flow.addChild("grid_data"); - addFloatArray(gv,"z",m_z.size(),m_z.begin(), - "m","length"); - vector_fp x(soln.nColumns()); + writelog("Importing datasets:\n"); + for (n = 0; n < nd; n++) { + XML_Node& fa = *d[n]; + nm = fa["title"]; + getFloatArray(fa,x,false); + if (nm == "u") { + writelog("axial velocity "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) { + soln[index(0,j)] = x[j]; + } + } + else { + goto error; + } + } + else if (nm == "z") { + ; // already read grid + } + else if (nm == "V") { + writelog("radial velocity "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) + soln[index(1,j)] = x[j]; + } + else goto error; + } + else if (nm == "T") { + writelog("temperature "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) + soln[index(2,j)] = x[j]; - soln.getRow(0,x.begin()); - addFloatArray(gv,"u",x.size(),x.begin(),"m/s","velocity"); - - soln.getRow(1,x.begin()); - addFloatArray(gv,"V", - x.size(),x.begin(),"1/s","strainrate"); - - soln.getRow(2,x.begin()); - addFloatArray(gv,"T",x.size(),x.begin(),"K","temperature",0.0); - - soln.getRow(3,x.begin()); - addFloatArray(gv,"L",x.size(),x.begin(),"N/m^4"); - - for (k = 0; k < m_nsp; k++) { - soln.getRow(4+k,x.begin()); - addFloatArray(gv,m_thermo->speciesName(k), - x.size(),x.begin(),"","massFraction",0.0,1.0); + // For fixed-temperature simulations, use the imported temperature profile by default. + // If this is not desired, call setFixedTempProfile *after* restoring the solution. + vector_fp zz(np); + for (int jj = 0; jj < np; jj++) zz[jj] = (grid(jj) - zmin())/(zmax() - zmin()); + setFixedTempProfile(zz, x); + } + else goto error; + } + else if (nm == "L") { + writelog("lambda "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) + soln[index(3,j)] = x[j]; + } + else goto error; + } + else if (m_thermo->speciesIndex(nm) >= 0) { + writelog(nm+" "); + if ((int) x.size() == np) { + k = m_thermo->speciesIndex(nm); + did_species[k] = 1; + for (j = 0; j < np; j++) + soln[index(k+4,j)] = x[j]; + } + } + else + ignored.push_back(nm); } -// XML_Node& inlt = flow.addChild("inlet"); -// addFloat(inlt,"T",m_inlet_T,"K","temperature",0.0); -// addFloat(inlt,"P",m_press,"Pa","pressure",0.0); -// for (k = 0; k < m_nsp; k++) { -// if (m_yin[k] != 0.0) -// addFloat(inlt, m_thermo->speciesName(k), m_yin[k], -// "", "massFraction",0.0,1.0); -// } + if (ignored.size() != 0) { + writelog("\n\n"); + writelog("Ignoring datasets:\n"); + int nn = ignored.size(); + for (int n = 0; n < nn; n++) { + writelog(ignored[n]+" "); + } + } - ofstream s(fname.c_str()); - if (!s) - throw CanteraError("save","could not open file "+fname); - ct->writeHeader(s); - ct->write(s); - s.close(); - writelog("Solution saved to file "+fname+" as solution '"+id+"'.\n"); - m_container->writeStats(); + for (ks = 0; ks < nsp; ks++) { + if (did_species[ks] == 0) { + if (!wrote_header) { + writelog("Missing data for species:\n"); + wrote_header = true; + } + writelog(m_thermo->speciesName(ks)+" "); + } + } + + return; + error: + throw CanteraError("StFlow::restore","Data size error"); } + void StFlow::save(XML_Node& o, doublereal* sol) { int k; ArrayViewer soln(m_nv, m_points, sol + loc()); - XML_Node& flow = (XML_Node&)o.addChild("flowfield"); + XML_Node& flow = (XML_Node&)o.addChild("domain"); flow.addAttribute("type",flowType()); flow.addAttribute("id",m_id); + flow.addAttribute("points",m_points); + flow.addAttribute("components",m_nv); + if (m_desc != "") addString(flow,"description",m_desc); XML_Node& gv = flow.addChild("grid_data"); addFloat(flow, "pressure", m_press, "Pa", "pressure"); diff --git a/Cantera/src/oneD/StFlow.h b/Cantera/src/oneD/StFlow.h index 63d5f494d..e3f19d7fc 100644 --- a/Cantera/src/oneD/StFlow.h +++ b/Cantera/src/oneD/StFlow.h @@ -15,7 +15,7 @@ #define CT_STFLOW_H #include "../transport/TransportBase.h" -#include "Resid1D.h" +#include "Domain1D.h" #include "../Array.h" #include "../sort.h" #include "../IdealGasPhase.h" @@ -61,7 +61,7 @@ namespace Cantera { * solution for a chemically-reacting, axisymmetric, * stagnation-point flow. */ - class StFlow : public Resid1D { + class StFlow : public Domain1D { public: @@ -191,13 +191,15 @@ namespace Cantera { string title, int zone); virtual void showSolution(ostream& s, const doublereal* x); + virtual void showSolution(const doublereal* x); - void save(string fname, string id, string desc, doublereal* soln); + //void save(string fname, string id, string desc, doublereal* soln); virtual void save(XML_Node& o, doublereal* sol); void restore(int job, string fname, string id, int& size_z, doublereal* z, int& size_soln, doublereal* soln); + virtual void restore(XML_Node& dom, doublereal* soln); // overloaded in subclasses virtual string flowType() { return ""; } @@ -381,15 +383,18 @@ namespace Cantera { // differencing, assuming u(z) is negative doublereal dVdz(const doublereal* x,int j) const { - return (V(x,j) - V(x,j-1))/m_dz[j-1]; + int jloc = (u(x,j) > 0.0 ? j : j + 1); + return (V(x,jloc) - V(x,jloc-1))/m_dz[jloc-1]; } doublereal dYdz(const doublereal* x,int k, int j) const { - return (Y(x,k,j) - Y(x,k,j-1))/m_dz[j-1]; + int jloc = (u(x,j) > 0.0 ? j : j + 1); + return (Y(x,k,jloc) - Y(x,k,jloc-1))/m_dz[jloc-1]; } doublereal dTdz(const doublereal* x,int j) const { - return (T(x,j) - T(x,j-1))/m_dz[j-1]; + int jloc = (u(x,j) > 0.0 ? j : j + 1); + return (T(x,jloc) - T(x,jloc-1))/m_dz[jloc-1]; } doublereal shear(const doublereal* x,int j) const { @@ -404,6 +409,10 @@ namespace Cantera { return -2.0*(c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/(z(j+1) - z(j-1)); } + int mindex(int k, int j, int m) { + return m*m_nsp*m_nsp + m_nsp*j + k; + } + void updateDiffFluxes(const doublereal* x, int j0, int j1); //--------------------------------------------------------- @@ -440,7 +449,8 @@ namespace Cantera { // transport properties vector_fp m_visc; vector_fp m_tcon; - Array2D m_diff; + vector_fp m_diff; + Array2D m_dthermal; Array2D m_flux; // production rates diff --git a/Cantera/src/oneD/Surf1D.h b/Cantera/src/oneD/Surf1D.h index a2feb1400..eb0456d20 100644 --- a/Cantera/src/oneD/Surf1D.h +++ b/Cantera/src/oneD/Surf1D.h @@ -3,7 +3,7 @@ dep #ifndef CT_SURF1D_H #define CT_SURF1D_H -#include "Resid1D.h" +#include "Domain1D.h" #include "SurfPhase.h" #include "InterfaceKinetics.h" #include "StFlow.h" @@ -16,11 +16,11 @@ namespace Cantera { // surface is zero-dimensional, and defined by a set of surface // species coverages. - class Surf1D : public Resid1D { + class Surf1D : public Domain1D { public: - Surf1D(InterfaceKinetics* skin = 0) : Resid1D(1, 1, 0.0) { + Surf1D(InterfaceKinetics* skin = 0) : Domain1D(1, 1, 0.0) { m_type = cSurfType; m_flow_left = 0; m_flow_right = 0; @@ -119,7 +119,7 @@ namespace Cantera { // check for left and right flow objects if (m_index > 0) { - Resid1D& r = container().domain(m_index-1); + Domain1D& r = container().domain(m_index-1); if (r.domainType() == cFlowType) { m_flow_left = (StFlow*)&r; m_left_nv = m_flow_left->nComponents(); @@ -143,7 +143,7 @@ namespace Cantera { } if (m_index < container().nDomains() - 1) { - Resid1D& r = container().domain(m_index+1); + Domain1D& r = container().domain(m_index+1); if (r.domainType() == cFlowType) { m_flow_right = (StFlow*)&r; m_right_nv = m_flow_right->nComponents(); diff --git a/Cantera/src/oneD/newton_utils.cpp b/Cantera/src/oneD/newton_utils.cpp index 8b20750c8..ab25326cd 100644 --- a/Cantera/src/oneD/newton_utils.cpp +++ b/Cantera/src/oneD/newton_utils.cpp @@ -8,7 +8,7 @@ #endif #include "../ct_defs.h" -#include "Resid1D.h" +#include "Domain1D.h" namespace Cantera { @@ -20,8 +20,13 @@ namespace Cantera { }; + /** + * Return a damping coefficient that keeps the solution after taking one + * Newton step between specified lower and upper bounds. This function only + * considers one domain. + */ doublereal bound_step(const doublereal* x, const doublereal* step, - Resid1D& r, int loglevel=0) { + Domain1D& r, int loglevel=0) { char buf[100]; int np = r.nPoints(); @@ -74,11 +79,38 @@ namespace Cantera { } + + /** + * This function computes the square of a weighted norm of a step + * vector for one domain. + * + * @param x Solution vector for this domain. + * @param step Newton step vector for this domain. + * @param r Object representing the domain. Used to get tolerances, number of components, + * and number of points. + * + * The return value is + * \f[ + * \sum_{n,j} \left(\frac{s_{n,j}}{w_n}\right)^2 + * \f] + * where the error weight for solution component \f$n\f$ is given by + * \f[ + * w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}. + * \f] + * Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for + * component \f$ n \f$, and multiplies the average magnitude of + * solution component n in the domain. The second term, \f$ + * \epsilon_{a,n}$, is the absolute error tolerance for component + * \f$ n \f$. + * + */ doublereal norm_square(const doublereal* x, - const doublereal* step, Resid1D& r) { + const doublereal* step, Domain1D& r) { doublereal f, ewt, esum, sum = 0.0; int n, j; - + doublereal f2max = 0.0; + int nmax = 0; + int jmax = 0; int nv = r.nComponents(); int np = r.nPoints(); @@ -89,13 +121,18 @@ namespace Cantera { for (j = 0; j < np; j++) { f = step[nv*j + n]/ewt; sum += f*f; -// if (fabs(f) > fmx) { -// fmx = fabs(f); -// jmx = j; -// nmx = n; -// } + if (f*f > f2max) { + jmax = j; + nmax = n; + f2max = f*f; + } } } +#undef DEBUG_NORM +#ifdef DEBUG_NORM + cout << "max step in domain " << r.id() << ": " << f2max << endl << + " for component " << r.componentName(nmax) << " at point " << jmax << endl; +#endif return sum; } } diff --git a/Cantera/src/oneD/refine.cpp b/Cantera/src/oneD/refine.cpp index 44a7018f0..4362e5bde 100644 --- a/Cantera/src/oneD/refine.cpp +++ b/Cantera/src/oneD/refine.cpp @@ -1,7 +1,7 @@ #include #include -#include "Resid1D.h" +#include "Domain1D.h" #include "refine.h" @@ -16,6 +16,11 @@ namespace Cantera { return false; } + static void drawline() { + string s(78,'#'); + s += '\n'; + writelog(s.c_str()); + } /** * Return the square root of machine precision. @@ -27,9 +32,9 @@ namespace Cantera { } - Refiner::Refiner(Resid1D& domain) : - m_ratio(10.0), m_slope(0.8), m_curve(0.8), m_min_range(0.01), - m_domain(&domain) + Refiner::Refiner(Domain1D& domain) : + m_ratio(10.0), m_slope(0.8), m_curve(0.8), m_min_range(0.001), + m_domain(&domain), m_npmax(200) { m_nv = m_domain->nComponents(); m_active.resize(m_nv, true); @@ -40,39 +45,48 @@ namespace Cantera { int Refiner::analyze(int n, const doublereal* z, const doublereal* x) { - if (m_domain->nPoints() <= 1) return 0; - m_nv = m_domain->nComponents(); + m_loc.clear(); + m_c.clear(); + m_keep.clear(); - //m_ok = false; + m_keep[0] = 1; + m_keep[n-1] = 1; + + + if (m_domain->nPoints() <= 1) return 0; + + m_nv = m_domain->nComponents(); // check consistency if (n != m_domain->nPoints()) return -1; - m_loc.clear(); - m_c.clear(); + + if (n >= m_npmax) return 0; /** * find locations where cell size ratio is too large. */ int j; vector_fp dz(n-1, 0.0); - dz[0] = z[1] - z[0]; - for (j = 1; j < n-1; j++) { - dz[j] = z[j+1] - z[j]; - if (dz[j] > m_ratio*dz[j-1]) { - m_loc[j] = 1; - m_c["point "+int2str(j)] = 1; - } - if (dz[j] < dz[j-1]/m_ratio) { - m_loc[j-1] = 1; - m_c["point "+int2str(j-1)] = 1; - } - } - string name; doublereal vmin, vmax, smin, smax, aa, ss; doublereal dmax, r; vector_fp v(n), s(n-1); + + dz[0] = z[1] - z[0]; + // for (j = 1; j < n-1; j++) { +// dz[j] = z[j+1] - z[j]; +// if (dz[j] > m_ratio*dz[j-1]) { +// m_loc[j] = 1; +// m_c["point "+int2str(j)] = 1; +// } +// if (dz[j] < dz[j-1]/m_ratio) { +// m_loc[j-1] = 1; +// m_c["point "+int2str(j-1)] = 1; +// } +// if (m_loc.size() + n > m_npmax) goto done; +// } + for (int i = 0; i < m_nv; i++) { //cout << i << " " << m_nv << " " << m_active[i] << endl; if (m_active[i]) { @@ -97,7 +111,6 @@ namespace Cantera { aa = fmaxx(abs(vmax), abs(vmin)); ss = fmaxx(abs(smax), abs(smin)); - // refine based on component i only if the range of v is // greater than a fraction 'min_range' of max |v|. This // eliminates components that consist of small fluctuations @@ -114,6 +127,11 @@ namespace Cantera { if (r > 1.0) { m_loc[j] = 1; m_c[name] = 1; + if (m_loc.size() + n > m_npmax) goto done; + } + if (r >= 0.0) { + m_keep[j] = 1; + m_keep[j+1] = 1; } } } @@ -128,22 +146,27 @@ namespace Cantera { // maximum allowable difference in slope between // adjacent points. - dmax = m_curve*(smax - smin); + dmax = m_curve*(smax - smin); // + 0.5*m_curve*(smax + smin); for (j = 0; j < n-2; j++) { r = abs(s[j+1] - s[j]) / (dmax + m_thresh/dz[j]); if (r > 1.0) { m_c[name] = 1; m_loc[j] = 1; m_loc[j+1] = 1; + if (m_loc.size() + n > m_npmax) goto done; + } + if (r >= 0.0) { + m_keep[j+1] = 1; } //cout << "at point " << j << " slope r = " // << r << " for " << name << endl // << " threshold = " << m_thresh << endl; } } - //cout << name << " " << m_curve << " " << smax << " " << smin << " " << ss << " " << m_min_range << endl; + } } +done: return m_loc.size(); } @@ -154,14 +177,16 @@ namespace Cantera { void Refiner::show() { int nnew = m_loc.size(); if (nnew > 0) { - writelog("Refining grid. " - "New points inserted after grid points "); + drawline(); + writelog(string("Refining grid in ") + + m_domain->id()+".\n" + +" New points inserted after grid points "); map::const_iterator b = m_loc.begin(); for (; b != m_loc.end(); ++b) { writelog(int2str(b->first)+" "); } writelog("\n"); - writelog("to resolve "); + writelog(" to resolve "); map::const_iterator bb = m_c.begin(); for (; bb != m_c.end(); ++bb) { writelog(string(bb->first)+" "); @@ -197,19 +222,5 @@ namespace Cantera { } zn[jn] = z[n-1]; return 0; - } - -// int npts = znew.size(); -// newsoln.resize(npts*ncomp); -// newsoln = Numeric.zeros((npts, ncomp),'d') -// for i in range(ncomp): -// for j in range(npts): -// newsoln[j,i] = interp.interp(znew[j],grid,solution[:,i]) - -// return (Numeric.array(znew), Numeric.array(znew), newsoln, self.ok) - - - - - + } } diff --git a/Cantera/src/oneD/refine.h b/Cantera/src/oneD/refine.h index 363a2b5ec..ec114a25d 100644 --- a/Cantera/src/oneD/refine.h +++ b/Cantera/src/oneD/refine.h @@ -3,13 +3,13 @@ namespace Cantera { - class Resid1D; + class Domain1D; class Refiner { public: - Refiner(Resid1D& domain); + Refiner(Domain1D& domain); virtual ~Refiner(){} void setCriteria(doublereal ratio = 10.0, @@ -18,7 +18,7 @@ namespace Cantera { m_ratio = ratio; m_slope = slope; m_curve = curve; } void setActive(int comp, bool state = true) { m_active[comp] = state; } - + void setMaxPoints(int npmax) { m_npmax = npmax; } int analyze(int n, const doublereal* z, const doublereal* x); int getNewGrid(int n, const doublereal* z, int nn, doublereal* znew); //int getNewSoln(int n, const doublereal* x, doublereal* xnew); @@ -27,17 +27,21 @@ namespace Cantera { bool newPointNeeded(int j) { return m_loc.find(j) != m_loc.end(); } + bool keepPoint(int j) { + return m_keep.find(j) != m_keep.end(); + } double value(const double* x, int i, int j); protected: map m_loc; + map m_keep; map m_c; vector m_active; doublereal m_ratio, m_slope, m_curve; doublereal m_min_range; - Resid1D* m_domain; - int m_nv; + Domain1D* m_domain; + int m_nv, m_npmax; doublereal m_thresh; }; diff --git a/Cantera/src/stringUtils.h b/Cantera/src/stringUtils.h index e80efb5f7..5e8935d73 100755 --- a/Cantera/src/stringUtils.h +++ b/Cantera/src/stringUtils.h @@ -28,6 +28,16 @@ namespace Cantera { string formatCompList(const Phase& mix, int xyc); string logfileName(const string& infile); string getFileName(const string& path); + + inline int intValue(string val) { + return atoi(stripws(val).c_str()); + } + + inline doublereal fpValue(string val) { + return atof(stripws(val).c_str()); + } + + } #endif diff --git a/Cantera/src/L_matrix.h b/Cantera/src/transport/L_matrix.h similarity index 100% rename from Cantera/src/L_matrix.h rename to Cantera/src/transport/L_matrix.h diff --git a/Makefile.in b/Makefile.in index cd0274f3c..e4a5a2fd7 100755 --- a/Makefile.in +++ b/Makefile.in @@ -17,7 +17,7 @@ build_matlab = @BUILD_MATLAB@ LIBDIR=@LIB_DIR@ # removed utils temporarily -all: kernel hdr-collect clib python matlab utils +all: kernel cxxlib hdr-collect clib python matlab utils install: hdr-install kernel-install data-install python-install matlab-install tools-install finish-install diff --git a/apps/MixMaster/KineticsFrame.py b/apps/MixMaster/KineticsFrame.py index 09b1068d8..378894a28 100644 --- a/apps/MixMaster/KineticsFrame.py +++ b/apps/MixMaster/KineticsFrame.py @@ -332,6 +332,8 @@ class ReactionPathFrame(Frame): self.b.build(element = el, diagram = self.d, dotfile = 'rxnpath.dot', format = 'dot') + self.b.build(element = el, diagram = self.d, + dotfile = 'rxnpath.txt', format = 'plain') if self.browser.get() == 1: fmt = self.fmt.get() diff --git a/config/configure b/config/configure index 4a164c617..7d8386f2e 100755 --- a/config/configure +++ b/config/configure @@ -1223,12 +1223,13 @@ if test "$BUILD_FORTRAN_90_INTERFACE" = "y"; then BUILD_F90=1; BUILD_CLIB=1; fi # Python Interface # BUILD_PYTHON=0 -if test -z "$PYTHON_CMD"; then - for ac_prog in python2 python +if test "$BUILD_PYTHON_INTERFACE" != "n"; then + if test -z "$PYTHON_CMD"; then + for ac_prog in python2 python do # Extract the first word of "$ac_prog", so it can be a program name with args. set dummy $ac_prog; ac_word=$2 -echo "$as_me:1231: checking for $ac_word" >&5 +echo "$as_me:1232: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_path_PYTHON_CMD+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -1245,7 +1246,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. if $as_executable_p "$ac_dir/$ac_word"; then ac_cv_path_PYTHON_CMD="$ac_dir/$ac_word" - echo "$as_me:1248: found $ac_dir/$ac_word" >&5 + echo "$as_me:1249: found $ac_dir/$ac_word" >&5 break fi done @@ -1256,10 +1257,10 @@ fi PYTHON_CMD=$ac_cv_path_PYTHON_CMD if test -n "$PYTHON_CMD"; then - echo "$as_me:1259: result: $PYTHON_CMD" >&5 + echo "$as_me:1260: result: $PYTHON_CMD" >&5 echo "${ECHO_T}$PYTHON_CMD" >&6 else - echo "$as_me:1262: result: no" >&5 + echo "$as_me:1263: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -1267,21 +1268,21 @@ fi done test -n "$PYTHON_CMD" || PYTHON_CMD=""none"" - if test "$PYTHON_CMD" != "none"; then BUILD_PYTHON=1; BUILD_CLIB=1; fi -else - BUILD_PYTHON=1 + if test "$PYTHON_CMD" != "none"; then BUILD_PYTHON=1; BUILD_CLIB=1; fi + else + BUILD_PYTHON=1 + fi fi -#AC_SUBST(PYTHON_CMD) - # # Matlab Interface # BUILD_MATLAB=0 -if test -z "$MATLAB_CMD"; then - # Extract the first word of "matlab", so it can be a program name with args. +if test "$BUILD_MATLAB_TOOLBOX" != "n"; then + if test -z "$MATLAB_CMD"; then + # Extract the first word of "matlab", so it can be a program name with args. set dummy matlab; ac_word=$2 -echo "$as_me:1284: checking for $ac_word" >&5 +echo "$as_me:1285: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_path_MATLAB_CMD+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -1298,7 +1299,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. if $as_executable_p "$ac_dir/$ac_word"; then ac_cv_path_MATLAB_CMD="$ac_dir/$ac_word" - echo "$as_me:1301: found $ac_dir/$ac_word" >&5 + echo "$as_me:1302: found $ac_dir/$ac_word" >&5 break fi done @@ -1310,18 +1311,16 @@ fi MATLAB_CMD=$ac_cv_path_MATLAB_CMD if test -n "$MATLAB_CMD"; then - echo "$as_me:1313: result: $MATLAB_CMD" >&5 + echo "$as_me:1314: result: $MATLAB_CMD" >&5 echo "${ECHO_T}$MATLAB_CMD" >&6 else - echo "$as_me:1316: result: no" >&5 + echo "$as_me:1317: result: no" >&5 echo "${ECHO_T}no" >&6 fi - if test "$MATLAB_CMD" != "none"; then BUILD_MATLAB=1; BUILD_CLIB=1; fi + if test "$MATLAB_CMD" != "none"; then BUILD_MATLAB=1; BUILD_CLIB=1; fi + fi fi -#if test "$BUILD_MATLAB_TOOLBOX" = "y"; then BUILD_MATLAB=1; BUILD_CLIB=1; fi - -# AC_SUBST(MATLAB_CMD) #------------------------------------------------- @@ -1337,7 +1336,7 @@ if test -n "$ac_tool_prefix"; then do # Extract the first word of "$ac_tool_prefix$ac_prog", so it can be a program name with args. set dummy $ac_tool_prefix$ac_prog; ac_word=$2 -echo "$as_me:1340: checking for $ac_word" >&5 +echo "$as_me:1339: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_CXX+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -1352,7 +1351,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_CXX="$ac_tool_prefix$ac_prog" -echo "$as_me:1355: found $ac_dir/$ac_word" >&5 +echo "$as_me:1354: found $ac_dir/$ac_word" >&5 break done @@ -1360,10 +1359,10 @@ fi fi CXX=$ac_cv_prog_CXX if test -n "$CXX"; then - echo "$as_me:1363: result: $CXX" >&5 + echo "$as_me:1362: result: $CXX" >&5 echo "${ECHO_T}$CXX" >&6 else - echo "$as_me:1366: result: no" >&5 + echo "$as_me:1365: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -1376,7 +1375,7 @@ if test -z "$CXX"; then do # Extract the first word of "$ac_prog", so it can be a program name with args. set dummy $ac_prog; ac_word=$2 -echo "$as_me:1379: checking for $ac_word" >&5 +echo "$as_me:1378: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_ac_ct_CXX+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -1391,7 +1390,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_ac_ct_CXX="$ac_prog" -echo "$as_me:1394: found $ac_dir/$ac_word" >&5 +echo "$as_me:1393: found $ac_dir/$ac_word" >&5 break done @@ -1399,10 +1398,10 @@ fi fi ac_ct_CXX=$ac_cv_prog_ac_ct_CXX if test -n "$ac_ct_CXX"; then - echo "$as_me:1402: result: $ac_ct_CXX" >&5 + echo "$as_me:1401: result: $ac_ct_CXX" >&5 echo "${ECHO_T}$ac_ct_CXX" >&6 else - echo "$as_me:1405: result: no" >&5 + echo "$as_me:1404: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -1414,27 +1413,27 @@ test -n "$ac_ct_CXX" || ac_ct_CXX="g++" fi # Provide some information about the compiler. -echo "$as_me:1417:" \ +echo "$as_me:1416:" \ "checking for C++ compiler version" >&5 ac_compiler=`set X $ac_compile; echo $2` -{ (eval echo "$as_me:1420: \"$ac_compiler --version &5\"") >&5 +{ (eval echo "$as_me:1419: \"$ac_compiler --version &5\"") >&5 (eval $ac_compiler --version &5) 2>&5 ac_status=$? - echo "$as_me:1423: \$? = $ac_status" >&5 + echo "$as_me:1422: \$? = $ac_status" >&5 (exit $ac_status); } -{ (eval echo "$as_me:1425: \"$ac_compiler -v &5\"") >&5 +{ (eval echo "$as_me:1424: \"$ac_compiler -v &5\"") >&5 (eval $ac_compiler -v &5) 2>&5 ac_status=$? - echo "$as_me:1428: \$? = $ac_status" >&5 + echo "$as_me:1427: \$? = $ac_status" >&5 (exit $ac_status); } -{ (eval echo "$as_me:1430: \"$ac_compiler -V &5\"") >&5 +{ (eval echo "$as_me:1429: \"$ac_compiler -V &5\"") >&5 (eval $ac_compiler -V &5) 2>&5 ac_status=$? - echo "$as_me:1433: \$? = $ac_status" >&5 + echo "$as_me:1432: \$? = $ac_status" >&5 (exit $ac_status); } cat >conftest.$ac_ext <<_ACEOF -#line 1437 "configure" +#line 1436 "configure" #include "confdefs.h" int @@ -1450,13 +1449,13 @@ ac_clean_files="$ac_clean_files a.out a.exe" # Try to create an executable without -o first, disregard a.out. # It will help us diagnose broken compilers, and finding out an intuition # of exeext. -echo "$as_me:1453: checking for C++ compiler default output" >&5 +echo "$as_me:1452: checking for C++ compiler default output" >&5 echo $ECHO_N "checking for C++ compiler default output... $ECHO_C" >&6 ac_link_default=`echo "$ac_link" | sed 's/ -o *conftest[^ ]*//'` -if { (eval echo "$as_me:1456: \"$ac_link_default\"") >&5 +if { (eval echo "$as_me:1455: \"$ac_link_default\"") >&5 (eval $ac_link_default) 2>&5 ac_status=$? - echo "$as_me:1459: \$? = $ac_status" >&5 + echo "$as_me:1458: \$? = $ac_status" >&5 (exit $ac_status); }; then # Find the output, starting from the most likely. This scheme is # not robust to junk in `.', hence go to wildcards (a.*) only as a last @@ -1479,34 +1478,34 @@ done else echo "$as_me: failed program was:" >&5 cat conftest.$ac_ext >&5 -{ { echo "$as_me:1482: error: C++ compiler cannot create executables" >&5 +{ { echo "$as_me:1481: error: C++ compiler cannot create executables" >&5 echo "$as_me: error: C++ compiler cannot create executables" >&2;} { (exit 77); exit 77; }; } fi ac_exeext=$ac_cv_exeext -echo "$as_me:1488: result: $ac_file" >&5 +echo "$as_me:1487: result: $ac_file" >&5 echo "${ECHO_T}$ac_file" >&6 # Check the compiler produces executables we can run. If not, either # the compiler is broken, or we cross compile. -echo "$as_me:1493: checking whether the C++ compiler works" >&5 +echo "$as_me:1492: checking whether the C++ compiler works" >&5 echo $ECHO_N "checking whether the C++ compiler works... $ECHO_C" >&6 # FIXME: These cross compiler hacks should be removed for Autoconf 3.0 # If not cross compiling, check that we can run a simple program. if test "$cross_compiling" != yes; then if { ac_try='./$ac_file' - { (eval echo "$as_me:1499: \"$ac_try\"") >&5 + { (eval echo "$as_me:1498: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:1502: \$? = $ac_status" >&5 + echo "$as_me:1501: \$? = $ac_status" >&5 (exit $ac_status); }; }; then cross_compiling=no else if test "$cross_compiling" = maybe; then cross_compiling=yes else - { { echo "$as_me:1509: error: cannot run C++ compiled programs. + { { echo "$as_me:1508: error: cannot run C++ compiled programs. If you meant to cross compile, use \`--host'." >&5 echo "$as_me: error: cannot run C++ compiled programs. If you meant to cross compile, use \`--host'." >&2;} @@ -1514,24 +1513,24 @@ If you meant to cross compile, use \`--host'." >&2;} fi fi fi -echo "$as_me:1517: result: yes" >&5 +echo "$as_me:1516: result: yes" >&5 echo "${ECHO_T}yes" >&6 rm -f a.out a.exe conftest$ac_cv_exeext ac_clean_files=$ac_clean_files_save # Check the compiler produces executables we can run. If not, either # the compiler is broken, or we cross compile. -echo "$as_me:1524: checking whether we are cross compiling" >&5 +echo "$as_me:1523: checking whether we are cross compiling" >&5 echo $ECHO_N "checking whether we are cross compiling... $ECHO_C" >&6 -echo "$as_me:1526: result: $cross_compiling" >&5 +echo "$as_me:1525: result: $cross_compiling" >&5 echo "${ECHO_T}$cross_compiling" >&6 -echo "$as_me:1529: checking for executable suffix" >&5 +echo "$as_me:1528: checking for executable suffix" >&5 echo $ECHO_N "checking for executable suffix... $ECHO_C" >&6 -if { (eval echo "$as_me:1531: \"$ac_link\"") >&5 +if { (eval echo "$as_me:1530: \"$ac_link\"") >&5 (eval $ac_link) 2>&5 ac_status=$? - echo "$as_me:1534: \$? = $ac_status" >&5 + echo "$as_me:1533: \$? = $ac_status" >&5 (exit $ac_status); }; then # If both `conftest.exe' and `conftest' are `present' (well, observable) # catch `conftest.exe'. For instance with Cygwin, `ls conftest' will @@ -1547,25 +1546,25 @@ for ac_file in `(ls conftest.exe; ls conftest; ls conftest.*) 2>/dev/null`; do esac done else - { { echo "$as_me:1550: error: cannot compute EXEEXT: cannot compile and link" >&5 + { { echo "$as_me:1549: error: cannot compute EXEEXT: cannot compile and link" >&5 echo "$as_me: error: cannot compute EXEEXT: cannot compile and link" >&2;} { (exit 1); exit 1; }; } fi rm -f conftest$ac_cv_exeext -echo "$as_me:1556: result: $ac_cv_exeext" >&5 +echo "$as_me:1555: result: $ac_cv_exeext" >&5 echo "${ECHO_T}$ac_cv_exeext" >&6 rm -f conftest.$ac_ext EXEEXT=$ac_cv_exeext ac_exeext=$EXEEXT -echo "$as_me:1562: checking for object suffix" >&5 +echo "$as_me:1561: checking for object suffix" >&5 echo $ECHO_N "checking for object suffix... $ECHO_C" >&6 if test "${ac_cv_objext+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 else cat >conftest.$ac_ext <<_ACEOF -#line 1568 "configure" +#line 1567 "configure" #include "confdefs.h" int @@ -1577,10 +1576,10 @@ main () } _ACEOF rm -f conftest.o conftest.obj -if { (eval echo "$as_me:1580: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:1579: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:1583: \$? = $ac_status" >&5 + echo "$as_me:1582: \$? = $ac_status" >&5 (exit $ac_status); }; then for ac_file in `(ls conftest.o conftest.obj; ls conftest.*) 2>/dev/null`; do case $ac_file in @@ -1592,24 +1591,24 @@ done else echo "$as_me: failed program was:" >&5 cat conftest.$ac_ext >&5 -{ { echo "$as_me:1595: error: cannot compute OBJEXT: cannot compile" >&5 +{ { echo "$as_me:1594: error: cannot compute OBJEXT: cannot compile" >&5 echo "$as_me: error: cannot compute OBJEXT: cannot compile" >&2;} { (exit 1); exit 1; }; } fi rm -f conftest.$ac_cv_objext conftest.$ac_ext fi -echo "$as_me:1602: result: $ac_cv_objext" >&5 +echo "$as_me:1601: result: $ac_cv_objext" >&5 echo "${ECHO_T}$ac_cv_objext" >&6 OBJEXT=$ac_cv_objext ac_objext=$OBJEXT -echo "$as_me:1606: checking whether we are using the GNU C++ compiler" >&5 +echo "$as_me:1605: checking whether we are using the GNU C++ compiler" >&5 echo $ECHO_N "checking whether we are using the GNU C++ compiler... $ECHO_C" >&6 if test "${ac_cv_cxx_compiler_gnu+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 else cat >conftest.$ac_ext <<_ACEOF -#line 1612 "configure" +#line 1611 "configure" #include "confdefs.h" int @@ -1624,16 +1623,16 @@ main () } _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:1627: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:1626: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:1630: \$? = $ac_status" >&5 + echo "$as_me:1629: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:1633: \"$ac_try\"") >&5 + { (eval echo "$as_me:1632: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:1636: \$? = $ac_status" >&5 + echo "$as_me:1635: \$? = $ac_status" >&5 (exit $ac_status); }; }; then ac_compiler_gnu=yes else @@ -1645,19 +1644,19 @@ rm -f conftest.$ac_objext conftest.$ac_ext ac_cv_cxx_compiler_gnu=$ac_compiler_gnu fi -echo "$as_me:1648: result: $ac_cv_cxx_compiler_gnu" >&5 +echo "$as_me:1647: result: $ac_cv_cxx_compiler_gnu" >&5 echo "${ECHO_T}$ac_cv_cxx_compiler_gnu" >&6 GXX=`test $ac_compiler_gnu = yes && echo yes` ac_test_CXXFLAGS=${CXXFLAGS+set} ac_save_CXXFLAGS=$CXXFLAGS CXXFLAGS="-g" -echo "$as_me:1654: checking whether $CXX accepts -g" >&5 +echo "$as_me:1653: checking whether $CXX accepts -g" >&5 echo $ECHO_N "checking whether $CXX accepts -g... $ECHO_C" >&6 if test "${ac_cv_prog_cxx_g+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 else cat >conftest.$ac_ext <<_ACEOF -#line 1660 "configure" +#line 1659 "configure" #include "confdefs.h" int @@ -1669,16 +1668,16 @@ main () } _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:1672: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:1671: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:1675: \$? = $ac_status" >&5 + echo "$as_me:1674: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:1678: \"$ac_try\"") >&5 + { (eval echo "$as_me:1677: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:1681: \$? = $ac_status" >&5 + echo "$as_me:1680: \$? = $ac_status" >&5 (exit $ac_status); }; }; then ac_cv_prog_cxx_g=yes else @@ -1688,7 +1687,7 @@ ac_cv_prog_cxx_g=no fi rm -f conftest.$ac_objext conftest.$ac_ext fi -echo "$as_me:1691: result: $ac_cv_prog_cxx_g" >&5 +echo "$as_me:1690: result: $ac_cv_prog_cxx_g" >&5 echo "${ECHO_T}$ac_cv_prog_cxx_g" >&6 if test "$ac_test_CXXFLAGS" = set; then CXXFLAGS=$ac_save_CXXFLAGS @@ -1715,7 +1714,7 @@ for ac_declaration in \ 'void exit (int);' do cat >conftest.$ac_ext <<_ACEOF -#line 1718 "configure" +#line 1717 "configure" #include "confdefs.h" #include $ac_declaration @@ -1728,16 +1727,16 @@ exit (42); } _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:1731: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:1730: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:1734: \$? = $ac_status" >&5 + echo "$as_me:1733: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:1737: \"$ac_try\"") >&5 + { (eval echo "$as_me:1736: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:1740: \$? = $ac_status" >&5 + echo "$as_me:1739: \$? = $ac_status" >&5 (exit $ac_status); }; }; then : else @@ -1747,7 +1746,7 @@ continue fi rm -f conftest.$ac_objext conftest.$ac_ext cat >conftest.$ac_ext <<_ACEOF -#line 1750 "configure" +#line 1749 "configure" #include "confdefs.h" $ac_declaration int @@ -1759,16 +1758,16 @@ exit (42); } _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:1762: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:1761: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:1765: \$? = $ac_status" >&5 + echo "$as_me:1764: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:1768: \"$ac_try\"") >&5 + { (eval echo "$as_me:1767: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:1771: \$? = $ac_status" >&5 + echo "$as_me:1770: \$? = $ac_status" >&5 (exit $ac_status); }; }; then break else @@ -1798,7 +1797,7 @@ ac_compiler_gnu=$ac_cv_c_compiler_gnu if test -n "$ac_tool_prefix"; then # Extract the first word of "${ac_tool_prefix}gcc", so it can be a program name with args. set dummy ${ac_tool_prefix}gcc; ac_word=$2 -echo "$as_me:1801: checking for $ac_word" >&5 +echo "$as_me:1800: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_CC+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -1813,7 +1812,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_CC="${ac_tool_prefix}gcc" -echo "$as_me:1816: found $ac_dir/$ac_word" >&5 +echo "$as_me:1815: found $ac_dir/$ac_word" >&5 break done @@ -1821,10 +1820,10 @@ fi fi CC=$ac_cv_prog_CC if test -n "$CC"; then - echo "$as_me:1824: result: $CC" >&5 + echo "$as_me:1823: result: $CC" >&5 echo "${ECHO_T}$CC" >&6 else - echo "$as_me:1827: result: no" >&5 + echo "$as_me:1826: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -1833,7 +1832,7 @@ if test -z "$ac_cv_prog_CC"; then ac_ct_CC=$CC # Extract the first word of "gcc", so it can be a program name with args. set dummy gcc; ac_word=$2 -echo "$as_me:1836: checking for $ac_word" >&5 +echo "$as_me:1835: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_ac_ct_CC+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -1848,7 +1847,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_ac_ct_CC="gcc" -echo "$as_me:1851: found $ac_dir/$ac_word" >&5 +echo "$as_me:1850: found $ac_dir/$ac_word" >&5 break done @@ -1856,10 +1855,10 @@ fi fi ac_ct_CC=$ac_cv_prog_ac_ct_CC if test -n "$ac_ct_CC"; then - echo "$as_me:1859: result: $ac_ct_CC" >&5 + echo "$as_me:1858: result: $ac_ct_CC" >&5 echo "${ECHO_T}$ac_ct_CC" >&6 else - echo "$as_me:1862: result: no" >&5 + echo "$as_me:1861: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -1872,7 +1871,7 @@ if test -z "$CC"; then if test -n "$ac_tool_prefix"; then # Extract the first word of "${ac_tool_prefix}cc", so it can be a program name with args. set dummy ${ac_tool_prefix}cc; ac_word=$2 -echo "$as_me:1875: checking for $ac_word" >&5 +echo "$as_me:1874: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_CC+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -1887,7 +1886,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_CC="${ac_tool_prefix}cc" -echo "$as_me:1890: found $ac_dir/$ac_word" >&5 +echo "$as_me:1889: found $ac_dir/$ac_word" >&5 break done @@ -1895,10 +1894,10 @@ fi fi CC=$ac_cv_prog_CC if test -n "$CC"; then - echo "$as_me:1898: result: $CC" >&5 + echo "$as_me:1897: result: $CC" >&5 echo "${ECHO_T}$CC" >&6 else - echo "$as_me:1901: result: no" >&5 + echo "$as_me:1900: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -1907,7 +1906,7 @@ if test -z "$ac_cv_prog_CC"; then ac_ct_CC=$CC # Extract the first word of "cc", so it can be a program name with args. set dummy cc; ac_word=$2 -echo "$as_me:1910: checking for $ac_word" >&5 +echo "$as_me:1909: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_ac_ct_CC+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -1922,7 +1921,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_ac_ct_CC="cc" -echo "$as_me:1925: found $ac_dir/$ac_word" >&5 +echo "$as_me:1924: found $ac_dir/$ac_word" >&5 break done @@ -1930,10 +1929,10 @@ fi fi ac_ct_CC=$ac_cv_prog_ac_ct_CC if test -n "$ac_ct_CC"; then - echo "$as_me:1933: result: $ac_ct_CC" >&5 + echo "$as_me:1932: result: $ac_ct_CC" >&5 echo "${ECHO_T}$ac_ct_CC" >&6 else - echo "$as_me:1936: result: no" >&5 + echo "$as_me:1935: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -1946,7 +1945,7 @@ fi if test -z "$CC"; then # Extract the first word of "cc", so it can be a program name with args. set dummy cc; ac_word=$2 -echo "$as_me:1949: checking for $ac_word" >&5 +echo "$as_me:1948: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_CC+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -1966,7 +1965,7 @@ if test "$ac_dir/$ac_word" = "/usr/ucb/cc"; then continue fi ac_cv_prog_CC="cc" -echo "$as_me:1969: found $ac_dir/$ac_word" >&5 +echo "$as_me:1968: found $ac_dir/$ac_word" >&5 break done @@ -1988,10 +1987,10 @@ fi fi CC=$ac_cv_prog_CC if test -n "$CC"; then - echo "$as_me:1991: result: $CC" >&5 + echo "$as_me:1990: result: $CC" >&5 echo "${ECHO_T}$CC" >&6 else - echo "$as_me:1994: result: no" >&5 + echo "$as_me:1993: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -2002,7 +2001,7 @@ if test -z "$CC"; then do # Extract the first word of "$ac_tool_prefix$ac_prog", so it can be a program name with args. set dummy $ac_tool_prefix$ac_prog; ac_word=$2 -echo "$as_me:2005: checking for $ac_word" >&5 +echo "$as_me:2004: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_CC+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -2017,7 +2016,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_CC="$ac_tool_prefix$ac_prog" -echo "$as_me:2020: found $ac_dir/$ac_word" >&5 +echo "$as_me:2019: found $ac_dir/$ac_word" >&5 break done @@ -2025,10 +2024,10 @@ fi fi CC=$ac_cv_prog_CC if test -n "$CC"; then - echo "$as_me:2028: result: $CC" >&5 + echo "$as_me:2027: result: $CC" >&5 echo "${ECHO_T}$CC" >&6 else - echo "$as_me:2031: result: no" >&5 + echo "$as_me:2030: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -2041,7 +2040,7 @@ if test -z "$CC"; then do # Extract the first word of "$ac_prog", so it can be a program name with args. set dummy $ac_prog; ac_word=$2 -echo "$as_me:2044: checking for $ac_word" >&5 +echo "$as_me:2043: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_ac_ct_CC+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -2056,7 +2055,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_ac_ct_CC="$ac_prog" -echo "$as_me:2059: found $ac_dir/$ac_word" >&5 +echo "$as_me:2058: found $ac_dir/$ac_word" >&5 break done @@ -2064,10 +2063,10 @@ fi fi ac_ct_CC=$ac_cv_prog_ac_ct_CC if test -n "$ac_ct_CC"; then - echo "$as_me:2067: result: $ac_ct_CC" >&5 + echo "$as_me:2066: result: $ac_ct_CC" >&5 echo "${ECHO_T}$ac_ct_CC" >&6 else - echo "$as_me:2070: result: no" >&5 + echo "$as_me:2069: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -2079,37 +2078,37 @@ fi fi -test -z "$CC" && { { echo "$as_me:2082: error: no acceptable cc found in \$PATH" >&5 +test -z "$CC" && { { echo "$as_me:2081: error: no acceptable cc found in \$PATH" >&5 echo "$as_me: error: no acceptable cc found in \$PATH" >&2;} { (exit 1); exit 1; }; } # Provide some information about the compiler. -echo "$as_me:2087:" \ +echo "$as_me:2086:" \ "checking for C compiler version" >&5 ac_compiler=`set X $ac_compile; echo $2` -{ (eval echo "$as_me:2090: \"$ac_compiler --version &5\"") >&5 +{ (eval echo "$as_me:2089: \"$ac_compiler --version &5\"") >&5 (eval $ac_compiler --version &5) 2>&5 ac_status=$? - echo "$as_me:2093: \$? = $ac_status" >&5 + echo "$as_me:2092: \$? = $ac_status" >&5 (exit $ac_status); } -{ (eval echo "$as_me:2095: \"$ac_compiler -v &5\"") >&5 +{ (eval echo "$as_me:2094: \"$ac_compiler -v &5\"") >&5 (eval $ac_compiler -v &5) 2>&5 ac_status=$? - echo "$as_me:2098: \$? = $ac_status" >&5 + echo "$as_me:2097: \$? = $ac_status" >&5 (exit $ac_status); } -{ (eval echo "$as_me:2100: \"$ac_compiler -V &5\"") >&5 +{ (eval echo "$as_me:2099: \"$ac_compiler -V &5\"") >&5 (eval $ac_compiler -V &5) 2>&5 ac_status=$? - echo "$as_me:2103: \$? = $ac_status" >&5 + echo "$as_me:2102: \$? = $ac_status" >&5 (exit $ac_status); } -echo "$as_me:2106: checking whether we are using the GNU C compiler" >&5 +echo "$as_me:2105: checking whether we are using the GNU C compiler" >&5 echo $ECHO_N "checking whether we are using the GNU C compiler... $ECHO_C" >&6 if test "${ac_cv_c_compiler_gnu+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 else cat >conftest.$ac_ext <<_ACEOF -#line 2112 "configure" +#line 2111 "configure" #include "confdefs.h" int @@ -2124,16 +2123,16 @@ main () } _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:2127: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:2126: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:2130: \$? = $ac_status" >&5 + echo "$as_me:2129: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:2133: \"$ac_try\"") >&5 + { (eval echo "$as_me:2132: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:2136: \$? = $ac_status" >&5 + echo "$as_me:2135: \$? = $ac_status" >&5 (exit $ac_status); }; }; then ac_compiler_gnu=yes else @@ -2145,19 +2144,19 @@ rm -f conftest.$ac_objext conftest.$ac_ext ac_cv_c_compiler_gnu=$ac_compiler_gnu fi -echo "$as_me:2148: result: $ac_cv_c_compiler_gnu" >&5 +echo "$as_me:2147: result: $ac_cv_c_compiler_gnu" >&5 echo "${ECHO_T}$ac_cv_c_compiler_gnu" >&6 GCC=`test $ac_compiler_gnu = yes && echo yes` ac_test_CFLAGS=${CFLAGS+set} ac_save_CFLAGS=$CFLAGS CFLAGS="-g" -echo "$as_me:2154: checking whether $CC accepts -g" >&5 +echo "$as_me:2153: checking whether $CC accepts -g" >&5 echo $ECHO_N "checking whether $CC accepts -g... $ECHO_C" >&6 if test "${ac_cv_prog_cc_g+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 else cat >conftest.$ac_ext <<_ACEOF -#line 2160 "configure" +#line 2159 "configure" #include "confdefs.h" int @@ -2169,16 +2168,16 @@ main () } _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:2172: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:2171: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:2175: \$? = $ac_status" >&5 + echo "$as_me:2174: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:2178: \"$ac_try\"") >&5 + { (eval echo "$as_me:2177: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:2181: \$? = $ac_status" >&5 + echo "$as_me:2180: \$? = $ac_status" >&5 (exit $ac_status); }; }; then ac_cv_prog_cc_g=yes else @@ -2188,7 +2187,7 @@ ac_cv_prog_cc_g=no fi rm -f conftest.$ac_objext conftest.$ac_ext fi -echo "$as_me:2191: result: $ac_cv_prog_cc_g" >&5 +echo "$as_me:2190: result: $ac_cv_prog_cc_g" >&5 echo "${ECHO_T}$ac_cv_prog_cc_g" >&6 if test "$ac_test_CFLAGS" = set; then CFLAGS=$ac_save_CFLAGS @@ -2215,16 +2214,16 @@ cat >conftest.$ac_ext <<_ACEOF #endif _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:2218: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:2217: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:2221: \$? = $ac_status" >&5 + echo "$as_me:2220: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:2224: \"$ac_try\"") >&5 + { (eval echo "$as_me:2223: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:2227: \$? = $ac_status" >&5 + echo "$as_me:2226: \$? = $ac_status" >&5 (exit $ac_status); }; }; then for ac_declaration in \ ''\ @@ -2236,7 +2235,7 @@ if { (eval echo "$as_me:2218: \"$ac_compile\"") >&5 'void exit (int);' do cat >conftest.$ac_ext <<_ACEOF -#line 2239 "configure" +#line 2238 "configure" #include "confdefs.h" #include $ac_declaration @@ -2249,16 +2248,16 @@ exit (42); } _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:2252: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:2251: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:2255: \$? = $ac_status" >&5 + echo "$as_me:2254: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:2258: \"$ac_try\"") >&5 + { (eval echo "$as_me:2257: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:2261: \$? = $ac_status" >&5 + echo "$as_me:2260: \$? = $ac_status" >&5 (exit $ac_status); }; }; then : else @@ -2268,7 +2267,7 @@ continue fi rm -f conftest.$ac_objext conftest.$ac_ext cat >conftest.$ac_ext <<_ACEOF -#line 2271 "configure" +#line 2270 "configure" #include "confdefs.h" $ac_declaration int @@ -2280,16 +2279,16 @@ exit (42); } _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:2283: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:2282: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:2286: \$? = $ac_status" >&5 + echo "$as_me:2285: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:2289: \"$ac_try\"") >&5 + { (eval echo "$as_me:2288: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:2292: \$? = $ac_status" >&5 + echo "$as_me:2291: \$? = $ac_status" >&5 (exit $ac_status); }; }; then break else @@ -2328,7 +2327,7 @@ ac_compiler_gnu=$ac_cv_c_compiler_gnu # AFS /usr/afsws/bin/install, which mishandles nonexistent args # SVR4 /usr/ucb/install, which tries to use the nonexistent group "staff" # ./install, which can be erroneously created by make from ./install.sh. -echo "$as_me:2331: checking for a BSD compatible install" >&5 +echo "$as_me:2330: checking for a BSD compatible install" >&5 echo $ECHO_N "checking for a BSD compatible install... $ECHO_C" >&6 if test -z "$INSTALL"; then if test "${ac_cv_path_install+set}" = set; then @@ -2377,7 +2376,7 @@ fi INSTALL=$ac_install_sh fi fi -echo "$as_me:2380: result: $INSTALL" >&5 +echo "$as_me:2379: result: $INSTALL" >&5 echo "${ECHO_T}$INSTALL" >&6 # Use test -z because SunOS4 sh mishandles braces in ${var-val}. @@ -2401,7 +2400,7 @@ if test -n "$ac_tool_prefix"; then do # Extract the first word of "$ac_tool_prefix$ac_prog", so it can be a program name with args. set dummy $ac_tool_prefix$ac_prog; ac_word=$2 -echo "$as_me:2404: checking for $ac_word" >&5 +echo "$as_me:2403: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_F77+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -2416,7 +2415,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_F77="$ac_tool_prefix$ac_prog" -echo "$as_me:2419: found $ac_dir/$ac_word" >&5 +echo "$as_me:2418: found $ac_dir/$ac_word" >&5 break done @@ -2424,10 +2423,10 @@ fi fi F77=$ac_cv_prog_F77 if test -n "$F77"; then - echo "$as_me:2427: result: $F77" >&5 + echo "$as_me:2426: result: $F77" >&5 echo "${ECHO_T}$F77" >&6 else - echo "$as_me:2430: result: no" >&5 + echo "$as_me:2429: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -2440,7 +2439,7 @@ if test -z "$F77"; then do # Extract the first word of "$ac_prog", so it can be a program name with args. set dummy $ac_prog; ac_word=$2 -echo "$as_me:2443: checking for $ac_word" >&5 +echo "$as_me:2442: checking for $ac_word" >&5 echo $ECHO_N "checking for $ac_word... $ECHO_C" >&6 if test "${ac_cv_prog_ac_ct_F77+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -2455,7 +2454,7 @@ for ac_dir in $ac_dummy; do test -z "$ac_dir" && ac_dir=. $as_executable_p "$ac_dir/$ac_word" || continue ac_cv_prog_ac_ct_F77="$ac_prog" -echo "$as_me:2458: found $ac_dir/$ac_word" >&5 +echo "$as_me:2457: found $ac_dir/$ac_word" >&5 break done @@ -2463,10 +2462,10 @@ fi fi ac_ct_F77=$ac_cv_prog_ac_ct_F77 if test -n "$ac_ct_F77"; then - echo "$as_me:2466: result: $ac_ct_F77" >&5 + echo "$as_me:2465: result: $ac_ct_F77" >&5 echo "${ECHO_T}$ac_ct_F77" >&6 else - echo "$as_me:2469: result: no" >&5 + echo "$as_me:2468: result: no" >&5 echo "${ECHO_T}no" >&6 fi @@ -2477,30 +2476,30 @@ done fi # Provide some information about the compiler. -echo "$as_me:2480:" \ +echo "$as_me:2479:" \ "checking for Fortran 77 compiler version" >&5 ac_compiler=`set X $ac_compile; echo $2` -{ (eval echo "$as_me:2483: \"$ac_compiler --version &5\"") >&5 +{ (eval echo "$as_me:2482: \"$ac_compiler --version &5\"") >&5 (eval $ac_compiler --version &5) 2>&5 ac_status=$? - echo "$as_me:2486: \$? = $ac_status" >&5 + echo "$as_me:2485: \$? = $ac_status" >&5 (exit $ac_status); } -{ (eval echo "$as_me:2488: \"$ac_compiler -v &5\"") >&5 +{ (eval echo "$as_me:2487: \"$ac_compiler -v &5\"") >&5 (eval $ac_compiler -v &5) 2>&5 ac_status=$? - echo "$as_me:2491: \$? = $ac_status" >&5 + echo "$as_me:2490: \$? = $ac_status" >&5 (exit $ac_status); } -{ (eval echo "$as_me:2493: \"$ac_compiler -V &5\"") >&5 +{ (eval echo "$as_me:2492: \"$ac_compiler -V &5\"") >&5 (eval $ac_compiler -V &5) 2>&5 ac_status=$? - echo "$as_me:2496: \$? = $ac_status" >&5 + echo "$as_me:2495: \$? = $ac_status" >&5 (exit $ac_status); } # If we don't use `.F' as extension, the preprocessor is not run on the # input file. ac_save_ext=$ac_ext ac_ext=F -echo "$as_me:2503: checking whether we are using the GNU Fortran 77 compiler" >&5 +echo "$as_me:2502: checking whether we are using the GNU Fortran 77 compiler" >&5 echo $ECHO_N "checking whether we are using the GNU Fortran 77 compiler... $ECHO_C" >&6 if test "${ac_cv_f77_compiler_gnu+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -2514,16 +2513,16 @@ else end _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:2517: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:2516: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:2520: \$? = $ac_status" >&5 + echo "$as_me:2519: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:2523: \"$ac_try\"") >&5 + { (eval echo "$as_me:2522: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:2526: \$? = $ac_status" >&5 + echo "$as_me:2525: \$? = $ac_status" >&5 (exit $ac_status); }; }; then ac_compiler_gnu=yes else @@ -2535,14 +2534,14 @@ rm -f conftest.$ac_objext conftest.$ac_ext ac_cv_f77_compiler_gnu=$ac_compiler_gnu fi -echo "$as_me:2538: result: $ac_cv_f77_compiler_gnu" >&5 +echo "$as_me:2537: result: $ac_cv_f77_compiler_gnu" >&5 echo "${ECHO_T}$ac_cv_f77_compiler_gnu" >&6 ac_ext=$ac_save_ext G77=`test $ac_compiler_gnu = yes && echo yes` ac_test_FFLAGS=${FFLAGS+set} ac_save_FFLAGS=$FFLAGS FFLAGS= -echo "$as_me:2545: checking whether $F77 accepts -g" >&5 +echo "$as_me:2544: checking whether $F77 accepts -g" >&5 echo $ECHO_N "checking whether $F77 accepts -g... $ECHO_C" >&6 if test "${ac_cv_prog_f77_g+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -2554,16 +2553,16 @@ cat >conftest.$ac_ext <<_ACEOF end _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:2557: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:2556: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:2560: \$? = $ac_status" >&5 + echo "$as_me:2559: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:2563: \"$ac_try\"") >&5 + { (eval echo "$as_me:2562: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:2566: \$? = $ac_status" >&5 + echo "$as_me:2565: \$? = $ac_status" >&5 (exit $ac_status); }; }; then ac_cv_prog_f77_g=yes else @@ -2574,7 +2573,7 @@ fi rm -f conftest.$ac_objext conftest.$ac_ext fi -echo "$as_me:2577: result: $ac_cv_prog_f77_g" >&5 +echo "$as_me:2576: result: $ac_cv_prog_f77_g" >&5 echo "${ECHO_T}$ac_cv_prog_f77_g" >&6 if test "$ac_test_FFLAGS" = set; then FFLAGS=$ac_save_FFLAGS @@ -2610,7 +2609,7 @@ ac_ext=f ac_compile='$F77 -c $FFLAGS conftest.$ac_ext >&5' ac_link='$F77 -o conftest$ac_exeext $FFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' ac_compiler_gnu=$ac_cv_f77_compiler_gnu -echo "$as_me:2613: checking how to get verbose linking output from $F77" >&5 +echo "$as_me:2612: checking how to get verbose linking output from $F77" >&5 echo $ECHO_N "checking how to get verbose linking output from $F77... $ECHO_C" >&6 if test "${ac_cv_prog_f77_v+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -2622,16 +2621,16 @@ cat >conftest.$ac_ext <<_ACEOF end _ACEOF rm -f conftest.$ac_objext -if { (eval echo "$as_me:2625: \"$ac_compile\"") >&5 +if { (eval echo "$as_me:2624: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? - echo "$as_me:2628: \$? = $ac_status" >&5 + echo "$as_me:2627: \$? = $ac_status" >&5 (exit $ac_status); } && { ac_try='test -s conftest.$ac_objext' - { (eval echo "$as_me:2631: \"$ac_try\"") >&5 + { (eval echo "$as_me:2630: \"$ac_try\"") >&5 (eval $ac_try) 2>&5 ac_status=$? - echo "$as_me:2634: \$? = $ac_status" >&5 + echo "$as_me:2633: \$? = $ac_status" >&5 (exit $ac_status); }; }; then ac_cv_prog_f77_v= # Try some options frequently used verbose output @@ -2653,7 +2652,7 @@ _ACEOF # flags. ac_save_FFLAGS=$FFLAGS FFLAGS="$FFLAGS $ac_verb" -(eval echo $as_me:2656: \"$ac_link\") >&5 +(eval echo $as_me:2655: \"$ac_link\") >&5 ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` echo "$ac_f77_v_output" >&5 FFLAGS=$ac_save_FFLAGS @@ -2685,21 +2684,21 @@ fi done done if test -z "$ac_cv_prog_f77_v"; then - { echo "$as_me:2688: WARNING: cannot determine how to obtain linking information from $F77" >&5 + { echo "$as_me:2687: WARNING: cannot determine how to obtain linking information from $F77" >&5 echo "$as_me: WARNING: cannot determine how to obtain linking information from $F77" >&2;} fi else echo "$as_me: failed program was:" >&5 cat conftest.$ac_ext >&5 -{ echo "$as_me:2694: WARNING: compilation failed" >&5 +{ echo "$as_me:2693: WARNING: compilation failed" >&5 echo "$as_me: WARNING: compilation failed" >&2;} fi rm -f conftest.$ac_objext conftest.$ac_ext fi -echo "$as_me:2700: result: $ac_cv_prog_f77_v" >&5 +echo "$as_me:2699: result: $ac_cv_prog_f77_v" >&5 echo "${ECHO_T}$ac_cv_prog_f77_v" >&6 -echo "$as_me:2702: checking for Fortran 77 libraries" >&5 +echo "$as_me:2701: checking for Fortran 77 libraries" >&5 echo $ECHO_N "checking for Fortran 77 libraries... $ECHO_C" >&6 if test "${ac_cv_flibs+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 @@ -2725,7 +2724,7 @@ _ACEOF # flags. ac_save_FFLAGS=$FFLAGS FFLAGS="$FFLAGS $ac_cv_prog_f77_v" -(eval echo $as_me:2728: \"$ac_link\") >&5 +(eval echo $as_me:2727: \"$ac_link\") >&5 ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` echo "$ac_f77_v_output" >&5 FFLAGS=$ac_save_FFLAGS @@ -2880,7 +2879,7 @@ esac fi # test "x$FLIBS" = "x" fi -echo "$as_me:2883: result: $ac_cv_flibs" >&5 +echo "$as_me:2882: result: $ac_cv_flibs" >&5 echo "${ECHO_T}$ac_cv_flibs" >&6 FLIBS="$ac_cv_flibs" @@ -2947,7 +2946,7 @@ fi # SO is the extension of shared libraries `(including the dot!) # -- usually .so, .sl on HP-UX, .dll on Cygwin -echo "$as_me:2950: checking SO" >&5 +echo "$as_me:2949: checking SO" >&5 echo $ECHO_N "checking SO... $ECHO_C" >&6 if test -z "$SO" then @@ -2958,7 +2957,7 @@ then *) SO=.so;; esac fi -echo "$as_me:2961: result: $SO" >&5 +echo "$as_me:2960: result: $SO" >&5 echo "${ECHO_T}$SO" >&6 ac_ext=cc @@ -2967,7 +2966,7 @@ ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' ac_compiler_gnu=$ac_cv_cxx_compiler_gnu -ac_config_files="$ac_config_files ../Cantera/Makefile ../Cantera/src/Makefile ../Cantera/src/zeroD/Makefile ../Cantera/src/oneD/Makefile ../Cantera/src/converters/Makefile ../Cantera/src/transport/Makefile ../Cantera/clib/src/Makefile ../Cantera/matlab/Makefile ../Cantera/python/Makefile ../Cantera/python/setup.py ../Cantera/cxx/Makefile ../Cantera/python/src/Makefile ../ext/lapack/Makefile ../ext/blas/Makefile ../ext/cvode/Makefile ../ext/math/Makefile ../ext/tpx/Makefile ../ext/Makefile ../ext/recipes/Makefile ../examples/Makefile ../examples/cxx/Makefile ../Makefile ../tools/Makefile ../tools/src/Makefile ../tools/src/sample.mak ../tools/templates/f77/demo.mak ../tools/testtools/Makefile ../data/inputs/Makefile ../test_problems/Makefile ../test_problems/cxx_ex/Makefile ../test_problems/silane_equil/Makefile" +ac_config_files="$ac_config_files ../Cantera/Makefile ../Cantera/src/Makefile ../Cantera/src/zeroD/Makefile ../Cantera/src/oneD/Makefile ../Cantera/src/converters/Makefile ../Cantera/src/transport/Makefile ../Cantera/clib/src/Makefile ../Cantera/matlab/Makefile ../Cantera/python/Makefile ../Cantera/python/setup.py ../Cantera/cxx/Makefile ../Cantera/python/src/Makefile ../ext/lapack/Makefile ../ext/blas/Makefile ../ext/cvode/Makefile ../ext/math/Makefile ../ext/tpx/Makefile ../ext/Makefile ../ext/recipes/Makefile ../examples/Makefile ../examples/cxx/Makefile ../Makefile ../tools/Makefile ../tools/src/Makefile ../tools/src/sample.mak ../tools/templates/f77/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 "x$prefix" = xNONE && prefix=$ac_default_prefix # Let make expand exec_prefix. @@ -2993,7 +2992,7 @@ DEFS=-DHAVE_CONFIG_H : ${CONFIG_STATUS=./config.status} ac_clean_files_save=$ac_clean_files ac_clean_files="$ac_clean_files $CONFIG_STATUS" -{ echo "$as_me:2996: creating $CONFIG_STATUS" >&5 +{ echo "$as_me:2995: creating $CONFIG_STATUS" >&5 echo "$as_me: creating $CONFIG_STATUS" >&6;} cat >$CONFIG_STATUS <<_ACEOF #! $SHELL @@ -3166,7 +3165,7 @@ cat >>$CONFIG_STATUS <<\EOF echo "$ac_cs_version"; exit 0 ;; --he | --h) # Conflict between --help and --header - { { echo "$as_me:3169: error: ambiguous option: $1 + { { echo "$as_me:3168: error: ambiguous option: $1 Try \`$0 --help' for more information." >&5 echo "$as_me: error: ambiguous option: $1 Try \`$0 --help' for more information." >&2;} @@ -3185,7 +3184,7 @@ Try \`$0 --help' for more information." >&2;} ac_need_defaults=false;; # This is an error. - -*) { { echo "$as_me:3188: error: unrecognized option: $1 + -*) { { echo "$as_me:3187: error: unrecognized option: $1 Try \`$0 --help' for more information." >&5 echo "$as_me: error: unrecognized option: $1 Try \`$0 --help' for more information." >&2;} @@ -3246,6 +3245,7 @@ do "../tools/src/Makefile" ) CONFIG_FILES="$CONFIG_FILES ../tools/src/Makefile" ;; "../tools/src/sample.mak" ) CONFIG_FILES="$CONFIG_FILES ../tools/src/sample.mak" ;; "../tools/templates/f77/demo.mak" ) CONFIG_FILES="$CONFIG_FILES ../tools/templates/f77/demo.mak" ;; + "../tools/templates/cxx/demo.mak" ) CONFIG_FILES="$CONFIG_FILES ../tools/templates/cxx/demo.mak" ;; "../tools/testtools/Makefile" ) CONFIG_FILES="$CONFIG_FILES ../tools/testtools/Makefile" ;; "../data/inputs/Makefile" ) CONFIG_FILES="$CONFIG_FILES ../data/inputs/Makefile" ;; "../test_problems/Makefile" ) CONFIG_FILES="$CONFIG_FILES ../test_problems/Makefile" ;; diff --git a/config/configure.in b/config/configure.in index 17b9134ae..c674c22fa 100755 --- a/config/configure.in +++ b/config/configure.in @@ -147,7 +147,7 @@ AC_SUBST(BLAS_LIBRARY) AC_SUBST(build_blas) -LOCAL_LIBS= +LOCAL_LIBS= if test -n "$NEED_ONED" then LOCAL_LIBS=$LOCAL_LIBS' '-loneD @@ -232,28 +232,28 @@ AC_SUBST(BUILD_F90) # Python Interface # BUILD_PYTHON=0 -if test -z "$PYTHON_CMD"; then - AC_PATH_PROGS(PYTHON_CMD, python2 python, "none") - if test "$PYTHON_CMD" != "none"; then BUILD_PYTHON=1; BUILD_CLIB=1; fi -else - BUILD_PYTHON=1 +if test "$BUILD_PYTHON_INTERFACE" != "n"; then + if test -z "$PYTHON_CMD"; then + AC_PATH_PROGS(PYTHON_CMD, python2 python, "none") + if test "$PYTHON_CMD" != "none"; then BUILD_PYTHON=1; BUILD_CLIB=1; fi + else + BUILD_PYTHON=1 + fi fi AC_SUBST(BUILD_PYTHON) -#AC_SUBST(PYTHON_CMD) - # # Matlab Interface # BUILD_MATLAB=0 -if test -z "$MATLAB_CMD"; then - AC_PATH_PROG(MATLAB_CMD, matlab, "none") - if test "$MATLAB_CMD" != "none"; then BUILD_MATLAB=1; BUILD_CLIB=1; fi +if test "$BUILD_MATLAB_TOOLBOX" != "n"; then + if test -z "$MATLAB_CMD"; then + AC_PATH_PROG(MATLAB_CMD, matlab, "none") + if test "$MATLAB_CMD" != "none"; then BUILD_MATLAB=1; BUILD_CLIB=1; fi + fi fi -#if test "$BUILD_MATLAB_TOOLBOX" = "y"; then BUILD_MATLAB=1; BUILD_CLIB=1; fi -AC_SUBST(BUILD_MATLAB) -# AC_SUBST(MATLAB_CMD) +AC_SUBST(BUILD_MATLAB) AC_SUBST(BUILD_CLIB) #------------------------------------------------- @@ -397,6 +397,7 @@ AC_OUTPUT(../Cantera/Makefile \ ../tools/src/Makefile \ ../tools/src/sample.mak \ ../tools/templates/f77/demo.mak \ + ../tools/templates/cxx/demo.mak \ ../tools/testtools/Makefile \ ../data/inputs/Makefile \ ../test_problems/Makefile \ diff --git a/configure b/configure index 31b69c763..4bb07a77b 100755 --- a/configure +++ b/configure @@ -50,9 +50,7 @@ BUILD_FORTRAN_90_INTERFACE='n' # Fortran is temporarily not working # to the path to the Python interpreter to use, if there is more than # one on your system. -#BUILD_PYTHON_INTERFACE=${BUILD_PYTHON_INTERFACE:="y"} -#PYTHON_CMD=${PYTHON_CMD:=python} - +BUILD_PYTHON_INTERFACE=${BUILD_PYTHON_INTERFACE:="y"} #----------- Matlab -------------------------------------------------- @@ -60,10 +58,7 @@ BUILD_FORTRAN_90_INTERFACE='n' # Fortran is temporarily not working # be installed on your system first, since the build process runs a # Matlab script. -# BUILD_MATLAB_TOOLBOX=${BUILD_MATLAB_TOOLBOX:="y"} -# MATLAB_CMD=${MATLAB_CMD:=matlab} - - +BUILD_MATLAB_TOOLBOX=${BUILD_MATLAB_TOOLBOX:="y"} #---------------------------------------------------------------------- # Kernel Configuration diff --git a/examples/cxx/Makefile.in b/examples/cxx/Makefile.in index d049f9446..0946fd80a 100755 --- a/examples/cxx/Makefile.in +++ b/examples/cxx/Makefile.in @@ -40,7 +40,7 @@ CXX = @CXX@ CXX_FLAGS = @CXXFLAGS@ # external libraries -EXT_LIBS = @LOCAL_LIBS@ +EXT_LIBS = @LOCAL_LIBS@ -lctcxx # Ending C++ linking libraries LCXX_END_LIBS = @LCXX_END_LIBS@ diff --git a/include/onedim.h b/include/onedim.h index 242c3013f..b12f8f06b 100755 --- a/include/onedim.h +++ b/include/onedim.h @@ -1,6 +1,7 @@ #ifndef CT_INCL_ONEDIM_H #define CT_INCL_ONEDIM_H +#include "kernel/oneD/Sim1D.h" #include "kernel/oneD/OneDim.h" #include "kernel/oneD/Resid1D.h" #include "kernel/oneD/Inlet1D.h" diff --git a/test_problems/silane_equil/Makefile.in b/test_problems/silane_equil/Makefile.in index 4ddd5cbe4..74d9593d9 100644 --- a/test_problems/silane_equil/Makefile.in +++ b/test_problems/silane_equil/Makefile.in @@ -30,7 +30,7 @@ CXX = @CXX@ CXX_FLAGS = @CXXFLAGS@ # external libraries -EXT_LIBS = @LOCAL_LIBS@ +EXT_LIBS = @LOCAL_LIBS@ -lctcxx # Ending C++ linking libraries LCXX_END_LIBS = @LCXX_END_LIBS@ diff --git a/tools/Makefile.in b/tools/Makefile.in index 4ffddb4eb..67a5587c0 100644 --- a/tools/Makefile.in +++ b/tools/Makefile.in @@ -33,4 +33,7 @@ templates-install: @INSTALL@ -m 644 templates/f77/*.cpp @prefix@/cantera/templates/f77 @INSTALL@ -m 644 templates/f77/*.f @prefix@/cantera/templates/f77 @INSTALL@ -m 644 templates/f77/*.mak @prefix@/cantera/templates/f77 + @INSTALL@ -d @prefix@/cantera/templates/cxx + @INSTALL@ -m 644 templates/cxx/*.cpp @prefix@/cantera/templates/cxx + @INSTALL@ -m 644 templates/cxx/*.mak @prefix@/cantera/templates/cxx chown -R @username@ @prefix@/cantera/templates \ No newline at end of file diff --git a/tools/src/Makefile.in b/tools/src/Makefile.in index dd8936778..71ab8008a 100755 --- a/tools/src/Makefile.in +++ b/tools/src/Makefile.in @@ -6,7 +6,7 @@ BINDIR = @ctroot@/bin build_ck = @BUILD_CK@ LCXX_FLAGS = -L$(LIBDIR) @CXXFLAGS@ -LOCAL_LIBS = -lcantera @math_libs@ @LAPACK_LIBRARY@ @BLAS_LIBRARY@ +LOCAL_LIBS = -lcantera @math_libs@ @LAPACK_LIBRARY@ @BLAS_LIBRARY@ -lctcxx LCXX_END_LIBS = @LCXX_END_LIBS@ OBJS = ck2ctml.o diff --git a/tools/templates/cxx/demo.mak.in b/tools/templates/cxx/demo.mak.in index 6c5838719..6ffc44ebc 100644 --- a/tools/templates/cxx/demo.mak.in +++ b/tools/templates/cxx/demo.mak.in @@ -33,7 +33,7 @@ CXX = @CXX@ CXX_FLAGS = @CXXFLAGS@ # external libraries -EXT_LIBS = @LOCAL_LIBS@ +EXT_LIBS = @LOCAL_LIBS@ -lctcxx # the directory where the Cantera libraries are located CANTERA_LIBDIR=@CANTERA_LIBDIR@ diff --git a/tools/templates/f77/demo.mak.in b/tools/templates/f77/demo.mak.in index d9e39dc18..6d030b2ac 100644 --- a/tools/templates/f77/demo.mak.in +++ b/tools/templates/f77/demo.mak.in @@ -47,7 +47,7 @@ CXX = @CXX@ CXX_FLAGS = @CXXFLAGS@ # external libraries -EXT_LIBS = @LOCAL_LIBS@ +EXT_LIBS = @LOCAL_LIBS@ -lctcxx # the directory where the Cantera libraries are located CANTERA_LIBDIR=@CANTERA_LIBDIR@ diff --git a/tools/testtools/Makefile.in b/tools/testtools/Makefile.in index 44988b8d0..7aa57a200 100644 --- a/tools/testtools/Makefile.in +++ b/tools/testtools/Makefile.in @@ -6,7 +6,7 @@ BINDIR = @ctroot@/bin build_ck = @BUILD_CK@ LCXX_FLAGS = -L$(LIBDIR) @CXXFLAGS@ -LOCAL_LIBS = -lcantera -lckreader @math_libs@ @LAPACK_LIBRARY@ @BLAS_LIBRARY@ +LOCAL_LIBS = -lcantera -lckreader @math_libs@ @LAPACK_LIBRARY@ @BLAS_LIBRARY@ -lctcxx LCXX_END_LIBS = @LCXX_END_LIBS@