*** empty log message ***

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)
 

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;
     }
 }

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<Sim1D>*   Cabinet<Sim1D>::__storage = 0;
+Cabinet<Domain1D>*   Cabinet<Domain1D>::__storage = 0;
+
+
+inline Sim1D* _sim1D(int i) {
+    return Cabinet<Sim1D>::cabinet()->item(i);
+}
+
+inline Domain1D* _domain(int i) {
+    return Cabinet<Domain1D>::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<Domain1D>::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<Domain1D>::cabinet()->add(i);
+        }
+        catch (CanteraError) { return -1; }
+    }
+
+    int DLL_EXPORT surf_new() {
+        try {
+            Surf1D* i = new Surf1D();
+            return Cabinet<Domain1D>::cabinet()->add(i);
+        }
+        catch (CanteraError) { return -1; }
+    }
+
+    int DLL_EXPORT symm_new() {
+        try {
+            Symm1D* i = new Symm1D();
+            return Cabinet<Domain1D>::cabinet()->add(i);
+        }
+        catch (CanteraError) { return -1; }
+    }
+
+    int DLL_EXPORT outlet_new() {
+        try {
+            Outlet1D* i = new Outlet1D();
+            return Cabinet<Domain1D>::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<Domain1D>::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<Domain1D*> d;
+        try {
+            for (int n = 0; n < nd; n++) {
+                d.push_back(_domain(domains[n]));
+            }
+            Sim1D* s = new Sim1D(d);
+            return Cabinet<Sim1D>::cabinet()->add(s);
+        }
+        catch (CanteraError) { return -1; }
+    }
+
+    int DLL_EXPORT sim1D_del(int i) {
+        Cabinet<Sim1D>::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; }
+    }
+
+}

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

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<Resid1D*> doms;
+        vector<Domain1D*> doms;
         for (i = 0; i < nd; i++) {
             switch (types[i]) {
             case 0:

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) 

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 '-'

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');

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
+

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

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);

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;
+

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);
+

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;
+

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);

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
+
+

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
+

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
+

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
+

Cantera/matlab/cantera/1D/@Domain1D/massFlux.m

@@ -0,0 +1,5 @@
+function mdot = massFlux(d)
+% MASSFLUX - 
+%   
+mdot = domain_methods(d.dom_id, 17);
+

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
+

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);

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);
+

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

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
+

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);

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
+

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);

Cantera/matlab/cantera/1D/@Domain1D/setMdot.m

@@ -0,0 +1,4 @@
+function d = setMdot(d, mdot)
+% SETMDOT - 
+%   
+domain_methods(d.dom_id, 60, mdot);

Cantera/matlab/cantera/1D/@Domain1D/setMoleFractions.m

@@ -0,0 +1,4 @@
+function d = setMoleFractions(d, x)
+% SETMOLEFRACTIONS - 
+%   
+domain_methods(d.dom_id, 62, x);

Cantera/matlab/cantera/1D/@Domain1D/setPressure.m

@@ -0,0 +1,5 @@
+function d = setPressure(d, p)
+% SETPRESSURE - 
+%   
+domain_methods(d.dom_id, 63, p);
+

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);
+

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);

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);

Cantera/matlab/cantera/1D/@Domain1D/setupGrid.m

@@ -0,0 +1,5 @@
+function d = setupGrid(d, grid)
+% SETUPGRID - 
+%   
+domain_methods(d.dom_id, 53, grid);
+

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);
+

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
+
+

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');
+

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);
+

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

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);

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);
+

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

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

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);

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);

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);

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);
+

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
+
+

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);
+

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);
+

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

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);
+

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

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);
+
+

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);

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

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
+
+

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
+
+

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

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
+
+

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");
         }

Cantera/matlab/cantera/private/onedimmethods.cpp

@@ -0,0 +1,368 @@
+
+#include "mex.h"
+#include "ctmatutils.h"
+#include "../../../clib/src/ctonedim.h"
+
+#include <iostream>
+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;
+    }
+}

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()

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:

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()
 
 

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)

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)
 

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
+

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    

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)

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)

Cantera/src/Flow1D.h

@@ -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
-
-

Cantera/src/Jac2.h

@@ -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
-
-

Cantera/src/Resid1D.h

@@ -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
-
-

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] + 

Cantera/src/importCTML.cpp

@@ -45,14 +45,6 @@ namespace Cantera {
     typedef vector<XML_Node*> 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);

Cantera/src/misc.cpp

@@ -18,6 +18,8 @@
 #include <fstream>
 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<string> inputDirs;
         vector<string> errorMessage;
@@ -40,6 +43,7 @@ namespace Cantera {
         size_t linelen;
         bool stop_on_error;
         bool write_log_to_cout;
+        bool matlab;
         map<string, string>     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();

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
+
+

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:
 
     };

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;

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;
             }
         }

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;

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<Resid1D*> domains) :
+    OneDim::OneDim(vector<Domain1D*> 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<Resid1D*>::iterator d; 
+        vector<Domain1D*>::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");
     }
 }

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<Resid1D*> domains);
+        OneDim(vector<Domain1D*> 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<Resid1D*> m_dom, m_connect, m_bulk;
+        vector<Domain1D*> 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:
 
     };
 

Cantera/src/oneD/Sim1D.cpp

@@ -7,7 +7,15 @@
 namespace Cantera {
 
 
-    Sim1D::Sim1D(vector<Resid1D*>& domains) : OneDim(domains) {
+    static void drawline() {
+        string s(78,'.');
+        s += '\n';
+        writelog(s.c_str());
+    }
+
+    Sim1D::Sim1D() : OneDim() {}
+
+    Sim1D::Sim1D(vector<Domain1D*>& 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<XML_Node*> 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;
         }

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<Resid1D*>& domains);
+        Sim1D(vector<Domain1D*>& 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
 
 

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<string> 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<XML_Node*> 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<string, double> params;
+        getFloats(dom, params);
+        setPressure(params["pressure"]);
+
+
+        vector<XML_Node*> 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"); 

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 "<none>"; }
@@ -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

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();

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;
     }
 }

Cantera/src/oneD/refine.cpp

@@ -1,7 +1,7 @@
 
 #include <map>
 #include <algorithm>
-#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<int, int>::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<string, int>::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)
-                               
-    
-                
-            
-    
+    }    
 }

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<int, int>       m_loc;
+        map<int, int>       m_keep;
         map<string, int>    m_c;
         vector<bool>        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;
 
     };

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

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
 

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()

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 \

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