From b4f34045a857c9be8ce0fe06e071b65ba66551a5 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Thu, 6 Dec 2012 17:31:44 +0000
Subject: [PATCH] [1D] Added a minimum grid spacing parameter

---
 include/cantera/oneD/Sim1D.h  |  8 ++++++++
 include/cantera/oneD/refine.h | 13 +++++++++++++
 src/oneD/Sim1D.cpp            | 14 ++++++++++++++
 src/oneD/refine.cpp           |  8 +++++---
 4 files changed, 40 insertions(+), 3 deletions(-)

diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h
index c3e89622a..b8281b5c3 100644
--- a/include/cantera/oneD/Sim1D.h
+++ b/include/cantera/oneD/Sim1D.h
@@ -107,6 +107,14 @@ public:
                            doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
     void setMaxGridPoints(int dom = -1, int npoints = 300);
 
+    //! Set the minimum grid spacing in the specified domain(s).
+    /*!
+     *  @param dom Domain index. If dom == -1, the specified spacing
+                   is applied to all domains.
+        @param gridmin The minimum allowable grid spacing [m]
+    */
+    void setGridMin(int dom, double gridmin);
+
     void restore(const std::string& fname, const std::string& id);
     void getInitialSoln();
 
diff --git a/include/cantera/oneD/refine.h b/include/cantera/oneD/refine.h
index 2d6cfad95..30255df71 100644
--- a/include/cantera/oneD/refine.h
+++ b/include/cantera/oneD/refine.h
@@ -31,6 +31,18 @@ public:
     void setMaxPoints(int npmax) {
         m_npmax = npmax;
     }
+
+    //! Set the minimum allowable spacing between adjacent grid points [m].
+    void setGridMin(double gridmin) {
+        m_gridmin = gridmin;
+    }
+
+    //! Returns the the minimum allowable spacing between adjacent
+    //! grid points [m].
+    double gridMin() const {
+        return m_gridmin;
+    }
+
     int analyze(size_t 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);
@@ -69,6 +81,7 @@ protected:
     Domain1D*            m_domain;
     size_t m_nv, m_npmax;
     doublereal          m_thresh;
+    doublereal m_gridmin; //!< minimum grid spacing [m]
 
 };
 
diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp
index de05987c2..40f6a15f8 100644
--- a/src/oneD/Sim1D.cpp
+++ b/src/oneD/Sim1D.cpp
@@ -624,6 +624,20 @@ void Sim1D::setRefineCriteria(int dom, doublereal ratio,
     }
 }
 
+void Sim1D::setGridMin(int dom, double gridmin)
+{
+    if (dom >= 0) {
+        Refiner& r = domain(dom).refiner();
+        r.setGridMin(gridmin);
+    } else {
+        for (size_t n = 0; n < m_nd; n++) {
+            Refiner& r = domain(n).refiner();
+            r.setGridMin(gridmin);
+        }
+    }
+}
+
+
 void Sim1D::setMaxGridPoints(int dom, int npoints)
 {
     if (dom >= 0) {
diff --git a/src/oneD/refine.cpp b/src/oneD/refine.cpp
index 44755f263..623b10b0e 100644
--- a/src/oneD/refine.cpp
+++ b/src/oneD/refine.cpp
@@ -18,7 +18,8 @@ static void r_drawline()
 
 Refiner::Refiner(Domain1D& domain) :
     m_ratio(10.0), m_slope(0.8), m_curve(0.8), m_prune(-0.001),
-    m_min_range(0.01), m_domain(&domain), m_npmax(3000)
+    m_min_range(0.01), m_domain(&domain), m_npmax(3000),
+    m_gridmin(5e-6)
 {
     m_nv = m_domain->nComponents();
     m_active.resize(m_nv, true);
@@ -107,7 +108,7 @@ int Refiner::analyze(size_t n, const doublereal* z,
                 dmax = m_slope*(vmax - vmin) + m_thresh;
                 for (j = 0; j < n-1; j++) {
                     r = fabs(v[j+1] - v[j])/dmax;
-                    if (r > 1.0) {
+                    if (r > 1.0 && dz[j] >= 2 * m_gridmin) {
                         m_loc[j] = 1;
                         m_c[name] = 1;
                         //if (int(m_loc.size()) + n > m_npmax) goto done;
@@ -139,7 +140,8 @@ int Refiner::analyze(size_t n, const doublereal* z,
                 dmax = m_curve*(smax - smin); // + 0.5*m_curve*(smax + smin);
                 for (j = 0; j < n-2; j++) {
                     r = fabs(s[j+1] - s[j]) / (dmax + m_thresh/dz[j]);
-                    if (r > 1.0) {
+                    if (r > 1.0 && dz[j] >= 2 * m_gridmin &&
+                        dz[j+1] >= 2 * m_gridmin) {
                         m_c[name] = 1;
                         m_loc[j] = 1;
                         m_loc[j+1] = 1;