From 57b38ce40fdc4dd29196ee2008e804a6ce5fdbea Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Tue, 14 Aug 2018 14:16:54 -0400
Subject: [PATCH] [1D] Use tighter tolerances and bounds for charged species

Looser tolerances can lead to instabilities, especially in cases where negative
concentrations of charged species are found at the end of the first solving
stage.
---
 interfaces/cython/cantera/onedim.pyx | 25 ++++++++++++++++++++++---
 src/oneD/IonFlow.cpp                 | 10 +++++++---
 2 files changed, 29 insertions(+), 6 deletions(-)

diff --git a/interfaces/cython/cantera/onedim.pyx b/interfaces/cython/cantera/onedim.pyx
index e3a3f48bd..505ad13f8 100644
--- a/interfaces/cython/cantera/onedim.pyx
+++ b/interfaces/cython/cantera/onedim.pyx
@@ -23,7 +23,7 @@ cdef class Domain1D:
 
         self.gas = phase
         self.gas._references[self._weakref_proxy] = True
-        self.have_user_tolerances = False
+        self.set_default_tolerances()
 
     property index:
         """
@@ -140,6 +140,14 @@ cdef class Domain1D:
         for name,(lower,upper) in kwargs.items():
             self.domain.setTransientTolerances(lower, upper, self.component_index(name))
 
+    def set_default_tolerances(self):
+        """
+        Set all tolerances to their default values
+        """
+        self.set_steady_tolerances(default=(1e-4, 1e-9))
+        self.set_transient_tolerances(default=(1e-4, 1e-11))
+        self.have_user_tolerances = False
+
     def bounds(self, component):
         """
         Return the (lower, upper) bounds for a solution component.
@@ -576,6 +584,18 @@ cdef class IonFlow(_FlowBase):
             else:
                 (<CxxIonFlow*>self.flow).fixElectricField()
 
+    def set_default_tolerances(self):
+        super().set_default_tolerances()
+        chargetol = {}
+        for S in self.gas.species():
+            if S.composition == {'E': 1.0}:
+                chargetol[S.name] = (1e-5, 1e-20)
+            elif S.charge != 0:
+                chargetol[S.name] = (1e-5, 1e-16)
+        self.set_steady_tolerances(**chargetol)
+        self.set_transient_tolerances(**chargetol)
+        self.have_user_tolerances = False
+
 
 cdef class Sim1D:
     """
@@ -898,8 +918,7 @@ cdef class Sim1D:
             rtol_ts_final = [dom.transient_reltol() for dom in self.domains]
 
         for dom in self.domains:
-            dom.set_steady_tolerances(default=(1e-4, 1e-9))
-            dom.set_transient_tolerances(default=(1e-4, 1e-11))
+            dom.set_default_tolerances()
 
         # Do initial steps without Soret diffusion
         soret_doms = [dom for dom in self.domains if getattr(dom, 'soret_enabled', False)]
diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp
index 860b2cb09..f7dd79e51 100644
--- a/src/oneD/IonFlow.cpp
+++ b/src/oneD/IonFlow.cpp
@@ -43,9 +43,13 @@ IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
     // no bound for electric potential
     setBounds(c_offset_E, -1.0e20, 1.0e20);
 
-    // Set tighter negative species limit on electron concentration to avoid
-    // instabilities
-    setBounds(c_offset_Y + m_kElectron, -1e-16, 1.0);
+    // Set tighter negative species limit on charged species to avoid
+    // instabilities. Tolerance on electrons is even tighter to account for the
+    // low "molecular" weight.
+    for (size_t k : m_kCharge) {
+        setBounds(c_offset_Y + k, -1e-14, 1.0);
+    }
+    setBounds(c_offset_Y + m_kElectron, -1e-18, 1.0);
 
     m_refiner->setActive(c_offset_E, false);
     m_mobility.resize(m_nsp*m_points);