[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

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

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