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