[1D] Introduction of IonFlow flame class

tested successfully with gri30

data/inputs/gri30_ion.cti

@@ -0,0 +1,88 @@
+units(length='cm', time='s', quantity='mol', act_energy='cal/mol')
+
+ideal_gas(name='gas',
+          elements=' O  H  C  N  Ar  E',
+          species=['''gri30: H2  H  O  O2  OH  H2O  HO2  H2O2  C  CH 
+                   CH2  CH2(S)  CH3  CH4  CO  CO2  HCO  CH2O  CH2OH  CH3O 
+                   CH3OH  C2H  C2H2  C2H3  C2H4  C2H5  C2H6  HCCO  CH2CO  HCCOH 
+                   N  NH  NH2  NH3  NNH  NO  NO2  N2O  HNO  CN 
+                   HCN  H2CN  HCNN  HCNO  HOCN  HNCO  NCO  N2  AR  C3H7 
+                   C3H8  CH2CHO  CH3CHO''',
+                   'HCO+    H3O+    E'],
+          reactions=['gri30: all', 'all'],
+          transport='Mix',
+          options=['skip_undeclared_species', 'skip_undeclared_third_bodies'],
+          initial_state=state(temperature=300.0, pressure=OneAtm))
+
+#-------------------------------------------------------------------------------
+# Species data
+#-------------------------------------------------------------------------------
+
+species(name = 'HCO+',
+    atoms = ' H:1  C:1  O:1  E:-1 ',
+    thermo = (
+       NASA( [  300.00,  1000.00], [  2.473973600E+00,   8.671559000E-03, 
+               -1.003150000E-05,   6.717052700E-09,  -1.787267400E-12,
+                9.914660800E+04,   8.175711870E+00] ),
+       NASA( [ 1000.00,  5000.00], [  3.741188000E+00,   3.344151700E-03, 
+               -1.239712100E-06,   2.118938800E-10,  -1.370415000E-14,
+                9.888407800E+04,   2.078613570E+00] )
+             ),
+        transport=gas_transport(geom='linear',
+                                diam=3.59,
+                                well_depth=498.0,
+                                polar=2.5,
+                                rot_relax=0.0),
+    note = 'J12/70')
+
+species(name = 'H3O+',
+    atoms = ' H:3  O:1  E:-1 ',
+    thermo = (
+       NASA( [  298.15,  1000.00], [  3.792952700E+00,  -9.108540000E-04, 
+                1.163635490E-05,  -1.213648870E-08,   4.261596630E-12,
+                7.075124010E+04,   1.471568560E+00] ),
+       NASA( [ 1000.00,  6000.00], [  2.496477160E+00,   5.728449200E-03, 
+               -1.839532810E-06,   2.735774390E-10,  -1.540939850E-14,
+                7.097291130E+04,   7.458507790E+00] )
+             ),
+        transport=gas_transport(geom='nonlinear',
+                                diam=2.605,
+                                well_depth=572.4,
+                                dipole=1.844,
+                                polar=1.5,
+                                rot_relax=2.1),
+     note = 'TPIS89')
+
+species(name = 'E',
+    atoms = ' E:1 ',
+    thermo = (
+       NASA( [  200.00,  1000.00], [  2.500000000E+00,   0.000000000E+00, 
+                0.000000000E+00,   0.000000000E+00,   0.000000000E+00,
+               -7.453750000E+02,  -1.172469020E+01] ),
+       NASA( [ 1000.00,  6000.00], [  2.500000000E+00,   0.000000000E+00, 
+                0.000000000E+00,   0.000000000E+00,   0.000000000E+00,
+               -7.453750000E+02,  -1.172469020E+01] )
+             ),
+        transport=gas_transport(geom='atom',
+                                diam=2.05,
+                                well_depth=145.0,
+                                polar=0.667,
+                                rot_relax=0.0),
+     note = 'gas L10/92')
+
+#-------------------------------------------------------------------------------
+# Reaction data
+#-------------------------------------------------------------------------------
+
+reaction('CH + O => HCO+ + E', [2.51E+11, 0.0, 1700])
+
+reaction('HCO+ + H2O => H3O+ + CO', [1.51E+15, 0.0, 0.0])
+
+reaction('H3O+ + E => H2O + H', [2.29E+18, -0.5, 0.0])
+
+reaction('H3O+ + E => OH + H + H', [7.95E+21, -1.4, 0.0])
+
+reaction('H3O+ + E => H2 + OH', [1.25E+19, -0.5, 0.0])
+
+reaction('H3O+ + E =>  O + H2 + H', [6.0E+17, -0.3, 0.0])
+

include/cantera/oneD/IonFlow.h

@@ -0,0 +1,130 @@
+//! @file IonFlow.h
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at http://www.cantera.org/license.txt for license and copyright information.
+
+#include "Domain1D.h"
+#include "cantera/base/Array.h"
+#include "cantera/thermo/IdealGasPhase.h"
+#include "cantera/kinetics/Kinetics.h"
+#include "cantera/oneD/StFlow.h"
+#include "cantera/oneD/Sim1D.h"
+#include "cantera/IdealGasMix.h"
+
+namespace Cantera
+{
+/**
+ * A class for ion flow.
+ * @ingroup onedim
+ */
+class IonFlow : public FreeFlame
+{
+public:
+    IonFlow(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1);
+
+    //! Turn electric field effect on/off
+    virtual void enableElectric(bool withElectric);
+    bool withElectric() const {
+        return m_do_electric;
+    }
+
+    virtual void setSolvingPhase(const size_t phase);
+
+    std::vector<size_t> chargeList() const {
+        return m_kCharge;
+    }
+
+    virtual void eval(size_t jg, doublereal* xg,
+              doublereal* rg, integer* diagg, doublereal rdt);
+
+    virtual void resize(size_t components, size_t points);
+
+    virtual void _finalize(const doublereal* x);
+
+    void solveSpeciesEqn(size_t k=npos);
+    void fixSpeciesMassFrac(size_t k=npos);
+    void solvePoissonEqn(size_t j=npos);
+    void fixElectricPotential(size_t j=npos);
+    void solveVelocity(size_t j=npos);
+    void fixVelocity(size_t j=npos);
+
+protected:
+    virtual void updateTransport(doublereal* x, size_t j0, size_t j1);
+    virtual void updateDiffFluxes(const doublereal* x, size_t j0, size_t j1);
+    virtual void evalPoisson(size_t j, doublereal* x, doublereal* r, integer* diag, doublereal rdt);
+    virtual void phaseOneDiffFluxes(const doublereal* x, size_t j0, size_t j1);       
+    virtual void phaseTwoDiffFluxes(const doublereal* x, size_t j0, size_t j1);
+    virtual void phaseThreeDiffFluxes(const doublereal* x, size_t j0, size_t j1);
+                   
+    bool m_do_electric;
+    std::vector<bool> m_do_velocity;
+    std::vector<bool> m_do_poisson;
+
+    // !electrical properties
+    vector_int m_speciesCharge;
+
+    // !index of species with charges
+    std::vector<size_t> m_kCharge;
+
+    // !index of neutral species
+    std::vector<size_t> m_kNeutral;
+
+    // mobility
+    vector_fp m_mobi;
+
+    // mass fraction of ion by equlibrium
+    Array2D m_yCharge;
+
+    // IonFlow solving phase
+    int m_solnPhase;
+
+    // !index of electron
+    size_t m_kElectron;
+
+    // fixed mass fraction value
+    vector_fp m_fixedMassFrac;
+
+    // fixed electric potential value
+    vector_fp m_fixedElecPoten;
+
+    // fixed velocity value
+    vector_fp m_fixedVelocity;
+
+    //! The fixed electric potential value at point j
+    doublereal phi_fixed(size_t j) const {
+        return m_fixedElecPoten[j];
+    }
+
+    //! The fixed mass fraction value at point j.
+    doublereal Y_fixed(size_t k, size_t j) const {
+        return m_fixedMassFrac[m_points*k+j];
+    }
+
+    //! The fixed velocity value at point j
+    doublereal u_fixed(size_t j) const {
+        return m_fixedVelocity[j];
+    }    
+
+    // electric potential
+    doublereal phi(const doublereal* x, size_t j) const {
+        return x[index(c_offset_P, j)];
+    }  
+
+    //electric field
+    doublereal E(const doublereal* x, size_t j) const {
+        return -(phi(x,j+1)-phi(x,j))/(z(j+1)-z(j));
+    }
+
+    doublereal dEdz(const doublereal* x, size_t j) const {
+        return 2*(E(x,j)-E(x,j-1))/(z(j+1)-z(j-1));
+    }
+
+    // number density
+    doublereal ND(const doublereal* x, size_t k, size_t j) const {
+        return Avogadro * m_rho[j] * Y(x,k,j) / m_wt[k];
+    }
+};
+
+}
+
+

include/cantera/oneD/StFlow.h

@@ -23,7 +23,8 @@ const size_t c_offset_U = 0; // axial velocity
 const size_t c_offset_V = 1; // strain rate
 const size_t c_offset_T = 2; // temperature
 const size_t c_offset_L = 3; // (1/r)dP/dr
-const size_t c_offset_Y = 4; // mass fractions
+const size_t c_offset_P = 4; // electric poisson's equation
+const size_t c_offset_Y = 5; // mass fractions
 
 class Transport;
 
@@ -179,7 +180,7 @@ public:
     }
 
     //! Change the grid size. Called after grid refinement.
-    void resize(size_t components, size_t points);
+    virtual void resize(size_t components, size_t points);
 
     virtual void setFixedPoint(int j0, doublereal t0) {}
 
@@ -343,7 +344,7 @@ protected:
     }
 
     //! Update the diffusive mass fluxes.
-    void updateDiffFluxes(const doublereal* x, size_t j0, size_t j1);
+    virtual void updateDiffFluxes(const doublereal* x, size_t j0, size_t j1);
 
     //---------------------------------------------------------
     //             member data
@@ -414,7 +415,7 @@ protected:
 
     //! Update the transport properties at grid points in the range from `j0`
     //! to `j1`, based on solution `x`.
-    void updateTransport(doublereal* x, size_t j0, size_t j1);
+    virtual void updateTransport(doublereal* x, size_t j0, size_t j1);
 
 private:
     vector_fp m_ybar;

src/oneD/IonFlow.cpp

@@ -0,0 +1,461 @@
+//! @file IonFlow.cpp
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at http://www.cantera.org/license.txt for license and copyright information.
+
+#include "cantera/oneD/IonFlow.h"
+#include "cantera/oneD/StFlow.h"
+#include "cantera/base/ctml.h"
+#include "cantera/transport/TransportBase.h"
+#include "cantera/numerics/funcs.h"
+#include "cantera/oneD/Domain1D.h"
+
+
+using namespace std;
+
+namespace Cantera
+{
+
+IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
+    FreeFlame(ph, nsp, points),
+    m_do_electric(false),
+    m_solnPhase(1)
+{
+    // make a local copy of species charge
+    for (size_t k = 0; k < m_nsp; k++) {
+        m_speciesCharge.push_back(m_thermo->charge(k));
+    }
+
+    // Find indices for charge of species
+    for (size_t k = 0; k < m_nsp; k++){
+        if (m_speciesCharge[k] != 0){
+            m_kCharge.push_back(k);
+        } else {
+            m_kNeutral.push_back(k);
+        }
+    }
+
+    // Find the index of electron 
+    if (m_thermo->speciesIndex("E") < m_nsp ) {
+        m_kElectron = m_thermo->speciesIndex("E");
+    }
+
+    // mass fraction bounds (strict bound for ions)
+    for (size_t k : m_kCharge) {
+        setBounds(c_offset_Y+k, -1.0e-20, 1e-5);
+    }
+    setBounds(c_offset_P, -1.0e20, 1.0e20);
+    m_refiner->setActive(c_offset_P, false);
+
+    m_mobi.resize(m_nsp*m_points);
+    m_do_poisson.resize(m_points,false);
+    m_do_velocity.resize(m_points,true);
+}
+
+void IonFlow::resize(size_t components, size_t points){
+    StFlow::resize(components, points);
+    m_mobi.resize(m_nsp*m_points);
+    m_do_species.resize(m_nsp,true);
+    m_do_poisson.resize(m_points,false);
+    m_do_velocity.resize(m_points,true);
+    m_fixedMassFrac.resize(m_points*m_nsp);
+    m_fixedElecPoten.resize(m_points,0.0);
+    m_fixedVelocity.resize(m_points);
+}
+
+void IonFlow::updateTransport(doublereal* x, size_t j0, size_t j1)
+{
+    StFlow::updateTransport(x,j0,j1);
+    for (size_t j = j0; j < j1; j++) {
+        setGasAtMidpoint(x,j);
+        m_trans->getMobilities(&m_mobi[j*m_nsp]);
+        m_mobi[m_kElectron+m_nsp*j] = 0.4;
+        m_diff[m_kElectron+m_nsp*j] = 0.4*(Boltzmann * T(x,j)) / ElectronCharge;
+    }
+}
+void IonFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
+{
+    if (m_solnPhase == 1) {
+        phaseOneDiffFluxes(x,j0,j1);
+    } else if (m_solnPhase == 2) {
+        phaseTwoDiffFluxes(x,j0,j1);
+    } else {
+        phaseThreeDiffFluxes(x,j0,j1);
+    }
+}
+
+void IonFlow::phaseOneDiffFluxes(const doublereal* x, size_t j0, size_t j1)
+{
+    for (size_t j = j0; j < j1; j++) {
+        double wtm = m_wtm[j];
+        double rho = density(j);
+        double dz = z(j+1) - z(j);
+        double sum = 0.0;
+        for (size_t k : m_kNeutral) {
+            m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
+            m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
+            sum -= m_flux(k,j);
+        }
+
+        // correction flux to insure that \sum_k Y_k V_k = 0.
+        for (size_t k : m_kNeutral) {  
+            m_flux(k,j) += sum*Y(x,k,j);
+        }
+
+        // flux for ions
+        for (size_t k : m_kCharge) {
+            m_flux(k,j) = 0;
+        }
+    }
+}
+
+void IonFlow::phaseTwoDiffFluxes(const doublereal* x, size_t j0, size_t j1)
+{
+    for (size_t j = j0; j < j1; j++) {
+        double wtm = m_wtm[j];
+        double rho = density(j);
+        double dz = z(j+1) - z(j);
+        // mixture-average diffusion
+        double sum_flux = 0.0;
+        for (size_t k = 0; k < m_nsp; k++) {       
+            m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
+            m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
+            sum_flux -= m_flux(k,j);
+        }
+
+        // ambipolar diffusion
+        double sum_chargeFlux = 0.0;
+        double sum = 0.0;
+        for (size_t k : m_kCharge) {
+            double Xav = 0.5 * (X(x,k,j+1) + X(x,k,j));
+            int q_k = m_speciesCharge[k];
+            sum_chargeFlux += m_speciesCharge[k] / m_wt[k] * m_flux(k,j);
+            sum += m_mobi[k+m_nsp*j] * Xav * q_k * q_k;
+        }
+        double drift;
+        double sum_drift = 0.0;
+        for (size_t k : m_kCharge) {
+            double Xav = 0.5 * (X(x,k,j+1) + X(x,k,j));
+            int q_k = m_speciesCharge[k];
+            drift = q_k * q_k * m_mobi[k+m_nsp*j] * Xav / sum;
+            drift *= -sum_chargeFlux * m_wt[k] / q_k;
+            m_flux(k,j) += drift;
+            sum_drift -= drift;
+        }
+
+        // correction flux
+        for (size_t k = 0; k < m_nsp; k++) {       
+            m_flux(k,j) += Y(x,k,j) * sum_flux;
+        }
+    }
+}
+
+void IonFlow::phaseThreeDiffFluxes(const doublereal* x, size_t j0, size_t j1)
+{
+    for (size_t j = j0; j < j1; j++) {
+        double wtm = m_wtm[j];
+        double rho = density(j);
+        double dz = z(j+1) - z(j);
+        // mixture-average diffusion
+        double sum = 0.0;
+        for (size_t k = 0; k < m_nsp; k++) {       
+            m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
+            m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
+            sum -= m_flux(k,j); 
+        }
+
+        // correction flux
+        for (size_t k = 0; k < m_nsp; k++) {       
+            m_flux(k,j) += Y(x,k,j) * sum;
+        }
+
+        // ambipolar diffusion
+        double drift;
+        double E_ambi = E(x,j);
+        sum = 0.0;
+        for (size_t k : m_kCharge) {
+            double Yav = 0.5 * (Y(x,k,j) + Y(x,k,j+1));
+            drift = rho * Yav * E_ambi;
+            drift *= m_speciesCharge[k] * m_mobi[k+m_nsp*j]; 
+            m_flux(k,j) += drift;
+            sum -= drift;
+        }
+
+        // correction drift
+        for (size_t k : m_kCharge) {
+            m_flux(k,j) += Y(x,k,j) * sum;
+        }
+    }   
+}
+
+void IonFlow::enableElectric(bool withElectric)
+{
+    m_do_electric = withElectric;
+}
+
+void IonFlow::setSolvingPhase(const size_t phase)
+{
+    m_solnPhase = phase;
+}
+
+void IonFlow::eval(size_t jg, doublereal* xg,
+                  doublereal* rg, integer* diagg, doublereal rdt)
+{
+    StFlow::eval(jg, xg, rg, diagg, rdt);
+    // start of local part of global arrays
+    doublereal* x = xg + loc();
+    doublereal* rsd = rg + loc();
+    integer* diag = diagg + loc();
+
+    size_t jmin, jmax;
+    if (jg == npos) { // evaluate all points
+        jmin = 0;
+        jmax = m_points - 1;
+    } else { // evaluate points for Jacobian
+        size_t jpt = (jg == 0) ? 0 : jg - firstPoint();
+        jmin = std::max<size_t>(jpt, 1) - 1;
+        jmax = std::min(jpt+1,m_points-1);
+    }
+    // the boundary points are not applied
+    for (size_t j = jmin; j <= jmax; j++) {
+        if (j == 0) {
+            rsd[index(c_offset_P, j)] = -phi(x,j);
+            diag[index(c_offset_P, j)] = 0;
+            for ( size_t k : m_kCharge) {
+                rsd[index(c_offset_Y + k, j)] = Y(x,k,j); 
+                diag[index(c_offset_Y + k, j)] = 0;
+            }
+        } else if (j == m_points - 1) {
+            rsd[index(c_offset_P, j)] = -phi(x,j);
+            diag[index(c_offset_P, j)] = 0;
+            for ( size_t k : m_kCharge) {
+                rsd[index(c_offset_Y + k, j)] = Y(x,k,j); 
+                diag[index(c_offset_Y + k, j)] = 0;
+            }
+        } else {
+            if (!m_do_velocity[j]) {
+                rsd[index(c_offset_U, j)] = u(x,j) - u_fixed(j);
+                diag[index(c_offset_U, j)] = 0;
+            }
+            for (size_t k = 0; k < m_nsp; k++) {
+                if (!m_do_species[k]) {
+                rsd[index(c_offset_Y + k, j)] = Y(x,k,j) - Y_fixed(k,j);
+                rsd[index(c_offset_Y + k, j)] -= rdt*(Y(x,k,j) - Y_prev(k,j));
+                diag[index(c_offset_Y + k, j)] = 1;            
+                }
+            }
+        }
+    }
+    
+    // convinent method due to interference
+    for (size_t j = jmin; j <= jmax; j++) {
+        if (j == 0) {
+            rsd[index(c_offset_P, j)] = -phi(x,j);
+            diag[index(c_offset_P, j)] = 0;
+        } else if (j == m_points - 1) {
+            rsd[index(c_offset_P, j)] = -phi(x,j);
+            diag[index(c_offset_P, j)] = 0;
+        } else {
+            if (m_do_poisson[j]) {
+                evalPoisson(j,x,rsd,diag,rdt);
+            } else {
+                rsd[index(c_offset_P, j)] = phi(x,j) - phi_fixed(j);
+                diag[index(c_offset_P, j)] = 0;
+            }
+        }
+    }
+}
+
+void IonFlow::evalPoisson(size_t j, doublereal* x, doublereal* rsd, integer* diag, doublereal rdt)
+{
+    //-----------------------------------------------
+    //    Poisson's equation
+    //
+    //    dE/dz = e/eps_0 * sum(q_k*n_k)
+    //
+    //    E = -dV/dz 
+    //-----------------------------------------------
+    doublereal chargeDensity = 0.0;
+    for (size_t k : m_kCharge) {
+        chargeDensity += m_speciesCharge[k] * ElectronCharge * ND(x,k,j);
+    }
+    rsd[index(c_offset_P, j)] = dEdz(x,j) - chargeDensity / epsilon_0;
+    diag[index(c_offset_P, j)] = 0;
+}
+
+void IonFlow::solveSpeciesEqn(size_t k)
+{
+    bool changed = false;
+    if (k == npos) {
+        for (size_t i = 0; i < m_nsp; i++) {
+            if (!m_do_energy[i]) {
+                changed = true;
+            }
+            m_do_species[i] = true;
+        }
+    } else {
+        if (!m_do_species[k]) {
+            changed = true;
+        }
+        m_do_species[k] = true;
+    }
+    if (changed) {
+        needJacUpdate();
+    }
+}
+
+void IonFlow::fixSpeciesMassFrac(size_t k)
+{
+    bool changed = false;
+    if (k == npos) {
+        for (size_t i = 0; i < m_nsp; i++) {
+            if (m_do_species[i]) {
+                changed = true;
+            }
+            m_do_species[i] = false;
+        }
+    } else {
+        if (m_do_species[k]) {
+            changed = true;
+        }
+        m_do_species[k] = false;
+    }
+    if (changed) {
+        needJacUpdate();
+    }
+}
+
+void IonFlow::solvePoissonEqn(size_t j)
+{
+    bool changed = false;
+    if (j == npos) {
+        for (size_t i = 0; i < m_points; i++) {
+            if (!m_do_poisson[i]) {
+                changed = true;
+            }
+            m_do_poisson[i] = true;
+        }
+    } else {
+        if (!m_do_poisson[j]) {
+            changed = true;
+        }
+        m_do_poisson[j] = true;
+    }
+    m_refiner->setActive(0, true);
+    m_refiner->setActive(1, true);
+    m_refiner->setActive(2, true);
+    if (changed) {
+        needJacUpdate();
+    }
+}
+
+void IonFlow::fixElectricPotential(size_t j)
+{
+    bool changed = false;
+    if (j == npos) {
+        for (size_t i = 0; i < m_points; i++) {
+            if (m_do_poisson[i]) {
+                changed = true;
+            }
+            m_do_poisson[i] = false;
+        }
+    } else {
+        if (m_do_poisson[j]) {
+            changed = true;
+        }
+        m_do_poisson[j] = false;
+    }
+    m_refiner->setActive(0, false);
+    m_refiner->setActive(1, false);
+    m_refiner->setActive(2, false);
+    if (changed) {
+        needJacUpdate();
+    }
+}
+
+void IonFlow::solveVelocity(size_t j)
+{
+    bool changed = false;
+    if (j == npos) {
+        for (size_t i = 0; i < m_points; i++) {
+            if (!m_do_velocity[i]) {
+                changed = true;
+            }
+            m_do_velocity[i] = true;
+        }
+    } else {
+        if (!m_do_velocity[j]) {
+            changed = true;
+        }
+        m_do_velocity[j] = true;
+    }
+    m_refiner->setActive(0, true);
+    m_refiner->setActive(1, true);
+    m_refiner->setActive(2, true);
+    if (changed) {
+        needJacUpdate();
+    }
+}
+
+void IonFlow::fixVelocity(size_t j)
+{
+    bool changed = false;
+    if (j == npos) {
+        for (size_t i = 0; i < m_points; i++) {
+            if (m_do_velocity[i]) {
+                changed = true;
+            }
+            m_do_velocity[i] = false;
+        }
+    } else {
+        if (m_do_velocity[j]) {
+            changed = true;
+        }
+        m_do_velocity[j] = false;
+    }
+    m_refiner->setActive(0, false);
+    m_refiner->setActive(1, false);
+    m_refiner->setActive(2, false);
+    if (changed) {
+        needJacUpdate();
+    }
+}
+
+void IonFlow::_finalize(const doublereal* x)
+{
+    FreeFlame::_finalize(x);
+
+    for (size_t k = 0; k < m_nsp; k++) {
+        bool y = m_do_species[k];
+        if (!y) {
+            for (size_t j = 0; j < m_points; j++) {
+                m_fixedMassFrac[m_points*k+j] = Y(x,k,j);
+            }
+        }
+    }
+
+    // This method is still not tested
+    // not sure why you want to return to original state
+    // if not doing on point zero
+    bool p = m_do_poisson[0];
+    for (size_t j = 0; j < m_points; j++) {
+        if (!p) {
+            m_fixedElecPoten[j] = phi(x, j);
+        }
+    }
+    if (p) {
+        solvePoissonEqn();
+    }
+
+    bool v = m_do_velocity[0];
+    for (size_t j = 0; j < m_points; j++) {
+        if (!v) {
+            m_fixedVelocity[j] = u(x,j);
+        }
+    }
+    if (v) {
+        solveVelocity();
+    }
+}
+
+}

src/oneD/StFlow.cpp

@@ -132,7 +132,8 @@ void StFlow::setupGrid(size_t n, const doublereal* z)
     }
 }
 
-void StFlow::resetBadValues(double* xg) {
+void StFlow::resetBadValues(double* xg) 
+{
     double* x = xg + loc();
     for (size_t j = 0; j < m_points; j++) {
         double* Y = x + m_nv*j + c_offset_Y;
@@ -368,10 +369,17 @@ void StFlow::eval(size_t jg, doublereal* xg,
                     -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
             }
             rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum;
+
+            // set residual of poisson's equ to zero
+            rsd[index(c_offset_P, 0)] = x[index(c_offset_P, j)];
         } else if (j == m_points - 1) {
             evalRightBoundary(x, rsd, diag, rdt);
+            // set residual of poisson's equ to zero
+            rsd[index(c_offset_P, j)] = x[index(c_offset_P, j)];
         } else { // interior points
             evalContinuity(j, x, rsd, diag, rdt);
+            // set residual of poisson's equ to zero
+            rsd[index(c_offset_P, j)] = x[index(c_offset_P, j)];
 
             //------------------------------------------------
             //    Radial momentum equation
@@ -546,6 +554,8 @@ string StFlow::componentName(size_t n) const
         return "T";
     case 3:
         return "lambda";
+    case 4:
+        return "ePotential";
     default:
         if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) {
             return m_thermo->speciesName(n - c_offset_Y);
@@ -565,6 +575,8 @@ size_t StFlow::componentIndex(const std::string& name) const
         return 2;
     } else if (name=="lambda") {
         return 3;
+    } else if (name == "ePotential") {
+        return 4;
     } else {
         for (size_t n=c_offset_Y; n<m_nsp+c_offset_Y; n++) {
             if (componentName(n)==name) {