From f18a811cf2392ff00f441671f5b8cbcb6adec667 Mon Sep 17 00:00:00 2001
From: Dave Goodwin <dgg@caltech.edu>
Date: Fri, 12 Nov 2004 23:37:32 +0000
Subject: [PATCH] changes by Karl Meredith to implement adiabatic,
 freely-propagating flames

---
 Cantera/src/oneD/Domain1D.h       |  36 +++-
 Cantera/src/oneD/Sim1D.cpp        | 136 ++++++++++++++++
 Cantera/src/oneD/Sim1D.h          | 262 +++++++++++++++---------------
 Cantera/src/oneD/StFlow.cpp       |  62 ++++++-
 Cantera/src/oneD/StFlow.h         |  21 ++-
 Cantera/src/oneD/boundaries1D.cpp |  39 ++++-
 6 files changed, 416 insertions(+), 140 deletions(-)

diff --git a/Cantera/src/oneD/Domain1D.h b/Cantera/src/oneD/Domain1D.h
index d7f75db8f..5953cc241 100644
--- a/Cantera/src/oneD/Domain1D.h
+++ b/Cantera/src/oneD/Domain1D.h
@@ -90,7 +90,7 @@ namespace Cantera {
          * Initialize. Base class method does nothing, but may be
          * overloaded. 
          */
-        virtual void init(){}
+        virtual void init(){  }
 
         virtual void setInitialState(doublereal* xlocal = 0){}
         virtual void setState(int point, const doublereal* state, doublereal* x) {}
@@ -186,6 +186,34 @@ namespace Cantera {
             }
         }
 
+        //added by Karl Meredith
+        void setTolerances(doublereal rtol, doublereal atol,int ts=0) {
+            for (int n=0;n<m_nv;n++){
+                if(ts>=0) {
+                    m_rtol_ss[n] = rtol;
+                    m_atol_ss[n] = atol;
+                }
+                if (ts <= 0) {
+                    m_rtol_ts[n] = rtol;
+                    m_atol_ts[n] = atol;
+                }
+            }
+        }
+        //added by Karl Meredith
+        void setTolerancesTS(doublereal rtol, doublereal atol) {
+            for (int n=0;n<m_nv;n++){
+                m_rtol_ts[n] = rtol;
+                m_atol_ts[n] = atol;
+            }
+        }
+        //added by Karl Meredith
+        void setTolerancesSS(doublereal rtol, doublereal atol) {
+            for (int n=0;n<m_nv;n++){
+                m_rtol_ss[n] = rtol;
+                m_atol_ss[n] = atol;
+            }
+        }
+				
         /// Relative tolerance of the nth component.
         doublereal rtol(int n) { return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]); }
 
@@ -426,7 +454,11 @@ namespace Cantera {
         //      "base class method _finalize called!");
         //}
 
+		//added by Karl Meredith
+		doublereal m_zfixed;
+		doublereal m_tfixed;
 
+		bool m_adiabatic;
     protected:
 
         doublereal m_rdt;
@@ -448,6 +480,8 @@ namespace Cantera {
         string m_id, m_desc;
         Refiner* m_refiner;
 
+
+
     private:
 
     };
diff --git a/Cantera/src/oneD/Sim1D.cpp b/Cantera/src/oneD/Sim1D.cpp
index 8c4d617be..b1f8619c1 100644
--- a/Cantera/src/oneD/Sim1D.cpp
+++ b/Cantera/src/oneD/Sim1D.cpp
@@ -33,6 +33,7 @@ namespace Cantera {
         m_xnew.resize(size(), 0.0);
         for (int n = 0; n < m_nd; n++) {
             domain(n)._getInitialSoln(m_x.begin() + start(n));
+			domain(n).m_adiabatic=false;
         }
 
         // set some defaults
@@ -45,6 +46,21 @@ namespace Cantera {
 
     }
 
+	
+    // added by Karl Meredith
+    void Sim1D::setInitialGuess(string component, vector_fp& locs, vector_fp& vals){
+        
+        for (int dom=0;dom<m_nd;dom++){
+            Domain1D& d = domain(dom);
+            int ncomp=d.nComponents();
+            for (int comp=0;comp<ncomp;comp++){
+                if(d.componentName(comp)==component){
+                    setProfile(dom,comp,locs,vals);
+                }
+            }
+        }
+    }
+
 
     /**
      * Set a single value in the solution vector.
@@ -396,6 +412,126 @@ namespace Cantera {
     }
 
 
+    /**
+     * Add node for fixed temperature point of freely propagating flame
+     */
+	//added by Karl Meredith
+    int Sim1D::setFixedTemperature(doublereal t) {
+        int np = 0;
+        vector_fp znew, xnew;
+        doublereal xmid;
+        doublereal zfixed,interp_factor;
+        doublereal z1,z2,t1,t2;
+        int strt, n, m, i;
+        int m1,m2;
+        vector_int dsize;
+        bool addnewpt=false;
+
+
+        for (n = 0; n < m_nd; n++) {
+            strt = znew.size();
+            Domain1D& d = domain(n);
+            
+            int comp = d.nComponents();
+            
+            // loop over points in the current grid to determine where new point is needed.
+            int npnow = d.nPoints();
+            int nstart = znew.size();
+            for (m = 0; m < npnow-1; m++) {
+                cout << "T["<<m<<"]="<<value(n,2,m)<<endl;
+                if(value(n,2,m)==t){
+                    zfixed=d.grid(m);
+                    //set d.zfixed, d.ztemp
+                    d.m_zfixed=zfixed;
+                    d.m_tfixed=t;
+                    cout << "T already fixed at "<<d.grid(m)<<endl;
+                    addnewpt=false;
+                    break;
+                }
+                else if((value(n,2,m)<t) && (value(n,2,m+1)>t)){
+                    cout << "T in between "<<value(n,2,m)<<" and "<<value(n,2,m+1)<<endl;
+                    z1=d.grid(m);
+                    m1=m;
+                    m2=m+1;
+                    z2=d.grid(m+1);
+                    t1=value(n,2,m);
+                    t2=value(n,2,m+1);
+                    
+                    zfixed=(z1-z2)/(t1-t2)*(t-t2)+z2;
+                    cout << zfixed<<endl;
+                    //set d.zfixed, d.ztemp;
+                    d.m_zfixed=zfixed;
+                    d.m_tfixed=t;
+                    addnewpt=true;
+                    break;
+                    //copy solution domain and push back values
+                }
+            }
+            
+            
+            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));
+                }
+                if(m==m1&&addnewpt){
+                    //add new point at zfixed
+                    znew.push_back(zfixed);
+                    np++;
+                    interp_factor=(zfixed-z2)/(z1-z2);
+                    // for each component, linearly interpolate
+                    // the solution to this point
+                    for (i = 0; i < comp; i++) {
+                        xmid = interp_factor*(value(n, i, m) - value(n, i, m+1))+value(n,i,m+1);
+                        xnew.push_back(xmid);
+                    }
+                }
+		
+		
+            }	
+            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.
+        
+        int gridstart = 0, gridsize;
+        for (n = 0; n < m_nd; n++) {
+            Domain1D& d = domain(n);
+            //            Refiner& r = d.refiner();
+            gridsize = dsize[n]; // d.nPoints() + r.nNewPoints();
+            d.setupGrid(gridsize, znew.begin() + gridstart);
+            gridstart += gridsize;
+        }
+        
+        // Replace the current solution vector with the new one
+        m_x.resize(xnew.size());
+        copy(xnew.begin(), xnew.end(), m_x.begin());
+
+        // resize the work array
+        m_xnew.resize(xnew.size());
+
+        copy(xnew.begin(), xnew.end(), m_xnew.begin());
+
+        resize();
+        finalize();
+        return np;
+    }
+
+    //added by Karl Meredith
+    void Sim1D::setAdiabaticFlame(void){
+        int n;
+        for (n = 0; n < m_nd; n++) {
+            Domain1D& d = domain(n);
+            d.m_adiabatic=true;
+        }
+    }
+
     /**
      * Set grid refinement criteria. If dom >= 0, then the settings
      * apply only to the specified domain.  If dom < 0, the settings
diff --git a/Cantera/src/oneD/Sim1D.h b/Cantera/src/oneD/Sim1D.h
index bdff0552e..c3a00e0a8 100644
--- a/Cantera/src/oneD/Sim1D.h
+++ b/Cantera/src/oneD/Sim1D.h
@@ -1,127 +1,135 @@
-/**
- * @file Sim1D.h
- */
-
-#ifndef CT_SIM1D_H
-#define CT_SIM1D_H
-
-#include "OneDim.h"
-#include "../funcs.h"
-
-namespace Cantera {
-
-    /**
-     * 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 is provided to make
-         * the class default-constructible, but is not meant to be
-         * used in most applications.  Use the next constructor
-         * instead.
-         */
-        Sim1D();
-
-
-        /**
-         * Standard constructor. 
-         * @param domains A vector of pointers to the domains to be linked together.
-         * 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<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);
-
-        /// 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);
-        void showSolution();
-
-        const doublereal* solution() { return m_x.begin(); }
-
-        void setTimeStep(doublereal stepsize, int n, integer* tsteps);
-
-        //void setMaxTimeStep(doublereal tmax) { m_maxtimestep = tmax; }
-
-        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, doublereal prune = -0.1);
-        void setMaxGridPoints(int dom = -1, int npoints = 300);
-
-        void restore(string fname, string id);
-        void getInitialSoln();
-
-        void setSolution(const doublereal* soln) {
-            copy(soln, soln + m_x.size(), m_x.begin());
-        }
-
-        const doublereal* solution() const { return m_x.begin(); }
-
-    protected:
-
-        vector_fp m_x;          // the solution vector
-        vector_fp m_xnew;       // a work array used to hold the residual 
-                                //      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
-
-
-    private:
-
-        /// Calls method _finalize in each domain.
-        void finalize();
-
-        void newtonSolve(int loglevel);
-
-
-    };
-
-}
-#endif
-
-
+/**
+ * @file Sim1D.h
+ */
+
+#ifndef CT_SIM1D_H
+#define CT_SIM1D_H
+
+#include "OneDim.h"
+#include "../funcs.h"
+
+namespace Cantera {
+
+    /**
+     * 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 is provided to make
+         * the class default-constructible, but is not meant to be
+         * used in most applications.  Use the next constructor
+         * instead.
+         */
+        Sim1D();
+
+
+        /**
+         * Standard constructor. 
+         * @param domains A vector of pointers to the domains to be linked together.
+         * 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<Domain1D*>& domains);
+
+        /// Destructor. Does nothing.
+        virtual ~Sim1D() {}
+
+        /**
+         * @name Setting initial values
+         *
+         * These methods are used to set the initial values of 
+         * solution components.
+         */
+        //@{
+
+        /// Set initial guess based on equilibrium
+        //added by Karl Meredith
+        void setInitialGuess(string component, vector_fp& locs, vector_fp& vals);
+
+        /// 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);
+
+        /// 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);
+        void showSolution();
+
+        const doublereal* solution() { return m_x.begin(); }
+
+        void setTimeStep(doublereal stepsize, int n, integer* tsteps);
+
+        //void setMaxTimeStep(doublereal tmax) { m_maxtimestep = tmax; }
+
+        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);
+
+        //added by Karl Meredith
+        int setFixedTemperature(doublereal t);
+        //added by Karl Meredith
+        void setAdiabaticFlame(void);
+
+        /// Set the criteria for grid refinement.
+        void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
+            doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
+        void setMaxGridPoints(int dom = -1, int npoints = 300);
+
+        void restore(string fname, string id);
+        void getInitialSoln();
+
+        void setSolution(const doublereal* soln) {
+            copy(soln, soln + m_x.size(), m_x.begin());
+        }
+
+        const doublereal* solution() const { return m_x.begin(); }
+
+    protected:
+
+        vector_fp m_x;          // the solution vector
+        vector_fp m_xnew;       // a work array used to hold the residual 
+                                //      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
+
+
+    private:
+
+        /// Calls method _finalize in each domain.
+        void finalize();
+
+        void newtonSolve(int loglevel);
+
+
+    };
+
+}
+#endif
+
+
diff --git a/Cantera/src/oneD/StFlow.cpp b/Cantera/src/oneD/StFlow.cpp
index fc5c61e28..77c36bfeb 100644
--- a/Cantera/src/oneD/StFlow.cpp
+++ b/Cantera/src/oneD/StFlow.cpp
@@ -241,7 +241,15 @@ namespace Cantera {
     }        
         
 
+/*	void StFlow::init() {
+	
 
+		cout << m_do_energy.begin()<< endl;
+	//	this->_getInitialSoln();
+		cout << "Initializing StFlow\n";
+	}
+
+*/
     void StFlow::setupGrid(int n, const doublereal* z) {        
         resize(n);
         int j;
@@ -466,6 +474,7 @@ namespace Cantera {
                 }
                 rsd[index(4,j)] = 1.0 - sum;
                 diag[index(4,j)] = 0;
+
             }
 
 
@@ -487,9 +496,36 @@ namespace Cantera {
                 //    d(\rho u)/dz + 2\rho V = 0
                 //
                 //------------------------------------------------
-                rsd[index(c_offset_U,j)] = 
-                    -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
-                    -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
+
+                //added by Karl Meredith
+                if(!m_adiabatic){
+                    rsd[index(c_offset_U,j)] = 
+                        -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
+                        -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
+                }
+                else{
+                    //we want mdot to propagate outward from fixed T point.
+                    if(grid(j)>m_zfixed){
+                        rsd[index(c_offset_U,j)] = 
+                            -(rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]
+                            -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
+                        //algebraic constraint    
+                        diag[index(c_offset_U, j)] = 0;
+                    }
+                    else if(grid(j)==m_zfixed){
+                        rsd[index(c_offset_U,j)] = 0.001*(T(x,j)-m_tfixed);
+                        //algebraic constraint    
+                        diag[index(c_offset_U, j)] = 0;    
+                    }
+                    else if(grid(j)<m_zfixed){
+                        rsd[index(c_offset_U,j)] = 
+                            -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
+                            -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
+                        //algebraic constraint    
+                        diag[index(c_offset_U, j)] = 0;
+                    }
+                }
+                //end of 'added by Karl Meredith'
 
 
                 //------------------------------------------------             
@@ -750,6 +786,26 @@ namespace Cantera {
     }
 
 
+	//added by Karl Meredith
+    int StFlow::componentIndex(string name) const {
+        
+
+        if(name=="u") {return 0;}
+        else if (name=="V") {return 1;}
+        else if (name=="T") {return 2;}
+        else if (name=="lambda") {return 3;}
+        else {
+            for (int n=4;n<m_nsp+4;n++){
+                if(componentName(n)==name){
+                    return n;
+                }	
+            }
+        }
+
+        return -1;
+    }
+
+
     void StFlow::restore(const XML_Node& dom, doublereal* soln) {
 
         vector<string> ignored;
diff --git a/Cantera/src/oneD/StFlow.h b/Cantera/src/oneD/StFlow.h
index 7ba66e3b8..bd576dd95 100644
--- a/Cantera/src/oneD/StFlow.h
+++ b/Cantera/src/oneD/StFlow.h
@@ -86,7 +86,10 @@ namespace Cantera {
         thermo_t& phase() { return *m_thermo; }
         kinetics_t& kinetics() { return *m_kin; }
 
-        /**
+        virtual void init(){
+		}
+		
+		/**
          * Set the thermo manager. Note that the flow equations assume
          * the ideal gas equation.
          */
@@ -169,6 +172,9 @@ namespace Cantera {
 
 
         virtual string componentName(int n) const;
+	    
+        //added by Karl Meredith
+        int componentIndex(string name) const;
 
 
         virtual void showSolution(const doublereal* x);
@@ -236,6 +242,14 @@ namespace Cantera {
         void setGas(const doublereal* x,int j);
         void setGasAtMidpoint(const doublereal* x,int j);
 
+        //Karl Meredith
+        //        doublereal density_unprotected(int j) const {
+        //    return m_rho[j];
+        // }
+        doublereal density(int j) const {
+            return m_rho[j];
+        }
+
 
     protected:
 
@@ -329,10 +343,6 @@ namespace Cantera {
             return m_wtm[j]*Y(x,k,j)/m_wt[k];
         }
 
-        doublereal density(int j) const {
-            return m_rho[j];
-        }
-
         doublereal flux(int k, int j) const {
             return m_flux(k, j);
         }
@@ -447,6 +457,7 @@ namespace Cantera {
 
         doublereal m_efctr;
 
+
     private:
 
         vector_fp m_ybar;
diff --git a/Cantera/src/oneD/boundaries1D.cpp b/Cantera/src/oneD/boundaries1D.cpp
index 71552be4f..fbac9675f 100644
--- a/Cantera/src/oneD/boundaries1D.cpp
+++ b/Cantera/src/oneD/boundaries1D.cpp
@@ -171,7 +171,22 @@ namespace Cantera {
         doublereal *xb, *rb;
 
         // residual equations for the two local variables
-        r[0] = m_mdot - x[0];
+
+        //added by Karl Meredith
+        if (m_adiabatic) 
+            // For the adiabatic case, the mass  flow rate is not known. For
+            // this case, set mdot (x[0]) to match rho*u in the flow domain
+            // dgg: I think this formulation will only work if the adiabatic
+            // inlet is on the left, i.e., the flow is left-to-right.
+            
+            r[0] = m_flow->density(0)*x[2] - x[0];
+        else
+            // Specified mass flow rate
+            r[0] = m_mdot    - x[0];
+
+        // The inlet temperature is always specified. For the adiabatic 
+        // case, this is the temperature of the gas far upstream of the 
+        // flame.
         r[1] = m_temp - x[1];
 
         // both are algebraic constraints
@@ -189,11 +204,18 @@ namespace Cantera {
             // spreading rate. Flow domain sets this to V(0),
             // so for finite spreading rate subtract m_V0.
             rb[1] -= m_V0;
-
-            rb[3] += x[0];       // lambda
+			
+			//added by Karl Meredith
+			if(m_adiabatic){
+				rb[3]=xb[3]; //zeroo lambda  (Avoids last species on last node being singular ???)
+			}
+			else{
+				rb[3] += x[0];       // lambda
+			}
             for (k = 1; k < m_nsp; k++) {
                 if (m_flow->doSpecies(k)) {
                     rb[4+k] += x[0]*m_yin[k];
+					//writelog("Left "+int2str(k)+" "+fp2str(m_yin[k])+"\n");
                 }
             }
         }
@@ -210,6 +232,7 @@ namespace Cantera {
                 if (m_flow->doSpecies(k)) {
                     //                    rb[4+k] += x[0]*(-xb[4+k] + m_yin[k]);
                     rb[4+k] += x[0]*(m_yin[k]);
+					//writelog("Right "+int2str(k)+" "+fp2str(m_yin[k])+"\n");
                 }
             }
         }                
@@ -438,7 +461,15 @@ namespace Cantera {
             db = diag - nc;
 
             // zero Lambda
-            rb[0] = xb[3];               // zero Lambda
+
+            //added by Karl Meredith
+            if (m_adiabatic) {
+                rb[0] = xb[0] - xb[0-nc];	//zero U gradient  
+                //(This makes it so that U at last node is not undefined)
+            }
+            else{
+                rb[0] = xb[3];          // zero Lambda
+            }
             rb[2] = xb[2] - xb[2 - nc];  // zero T gradient
             for (k = 5; k < nc; k++) {
                 rb[k] = xb[k] - xb[k - nc]; // zero mass fraction gradient