From ee633f3e16314210ecdd7edbfae8b519f274aac1 Mon Sep 17 00:00:00 2001 From: bangshiuh Date: Tue, 7 Aug 2018 22:44:17 -0400 Subject: [PATCH] [1D] Change Poisson equation to first-order difference --- include/cantera/oneD/IonFlow.h | 44 ++++-------- include/cantera/oneD/StFlow.h | 2 +- interfaces/cython/cantera/_cantera.pxd | 7 +- interfaces/cython/cantera/onedim.py | 33 +++------ interfaces/cython/cantera/onedim.pyx | 11 ++- interfaces/cython/cantera/test/test_onedim.py | 4 +- src/oneD/IonFlow.cpp | 71 ++++++++----------- src/oneD/StFlow.cpp | 10 +-- 8 files changed, 64 insertions(+), 118 deletions(-) diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h index d4005caff..85319fdb2 100644 --- a/include/cantera/oneD/IonFlow.h +++ b/include/cantera/oneD/IonFlow.h @@ -34,18 +34,16 @@ public: IonFlow(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1); //! set the solving stage virtual void setSolvingStage(const size_t phase); - //! set electric voltage at inlet and outlet - virtual void setElectricPotential(const double v1, const double v2); virtual void resize(size_t components, size_t points); virtual void _finalize(const double* x); - //! set to solve Poisson's equation on a point - void solvePoissonEqn(size_t j=npos); + //! set to solve electric field on a point + void solveElectricField(size_t j=npos); //! set to fix voltage on a point - void fixElectricPotential(size_t j=npos); - bool doPoisson(size_t j) { - return m_do_poisson[j]; + void fixElectricField(size_t j=npos); + bool doElectricField(size_t j) { + return m_do_electric_field[j]; } /** @@ -66,7 +64,7 @@ public: protected: /*! * This function overloads the original function. The residual function - * of Poisson's equation is added. + * of electric field is added. */ virtual void evalResidual(double* x, double* rsd, int* diag, double rdt, size_t jmin, size_t jmax); @@ -74,10 +72,11 @@ protected: virtual void updateDiffFluxes(const double* x, size_t j0, size_t j1); //! Solving phase one: the fluxes of charged species are turned off virtual void frozenIonMethod(const double* x, size_t j0, size_t j1); - //! Solving phase three: the Poisson's equation is added coupled by the electrical drift - virtual void poissonEqnMethod(const double* x, size_t j0, size_t j1); - //! flag for solving poisson's equation or not - std::vector m_do_poisson; + //! Solving phase two: the electric field equation is added coupled + //! by the electrical drift + virtual void electricFieldMethod(const double* x, size_t j0, size_t j1); + //! flag for solving electric field or not + std::vector m_do_electric_field; //! flag for importing transport of electron bool m_import_electron_transport; @@ -101,33 +100,16 @@ protected: //! solving stage int m_stage; - //! The voltage - double m_inletVoltage; - double m_outletVoltage; - //! index of electron size_t m_kElectron; - //! fixed electric potential value - vector_fp m_fixedElecPoten; - - //! The fixed electric potential value at point j - double phi_fixed(size_t j) const { - return m_fixedElecPoten[j]; - } - - //! electric potential - double phi(const double* x, size_t j) const { - return x[index(c_offset_P, j)]; - } - //! electric field double E(const double* x, size_t j) const { - return -(phi(x,j+1)-phi(x,j))/(z(j+1)-z(j)); + return x[index(c_offset_E, j)]; } double dEdz(const double* x, size_t j) const { - return 2*(E(x,j)-E(x,j-1))/(z(j+1)-z(j-1)); + return (E(x,j)-E(x,j-1))/(z(j)-z(j-1)); } //! number density diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h index 7b1d6534d..e0b08d61f 100644 --- a/include/cantera/oneD/StFlow.h +++ b/include/cantera/oneD/StFlow.h @@ -23,7 +23,7 @@ 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_P = 4; // electric poisson's equation +const size_t c_offset_E = 4; // electric poisson's equation const size_t c_offset_Y = 5; // mass fractions class Transport; diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 1f4138bca..0c23bede9 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -701,10 +701,9 @@ cdef extern from "cantera/oneD/IonFlow.h": cdef cppclass CxxIonFlow "Cantera::IonFlow": CxxIonFlow(CxxIdealGasPhase*, int, int) void setSolvingStage(int) - void setElectricPotential(const double, const double) - void solvePoissonEqn() - void fixElectricPotential() - cbool doPoisson(size_t) + void solveElectricField() + void fixElectricField() + cbool doElectricField(size_t) cdef extern from "cantera/oneD/Sim1D.h": diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index b8225b947..dccefd4ac 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -578,52 +578,35 @@ class IonFlameBase(FlameBase): T = self.T u = self.u V = self.V - phi = self.phi E = self.E csvfile = open(filename, 'w') writer = _csv.writer(csvfile) writer.writerow(['z (m)', 'u (m/s)', 'V (1/s)', 'T (K)', - 'phi (V)', 'E (V/m)', 'rho (kg/m3)'] + self.gas.species_names) + 'E (V/m)', 'rho (kg/m3)'] + self.gas.species_names) for n in range(self.flame.n_points): self.set_gas_state(n) - writer.writerow([z[n], u[n], V[n], T[n], phi[n], E[n], self.gas.density] + + writer.writerow([z[n], u[n], V[n], T[n], E[n], self.gas.density] + list(getattr(self.gas, species))) csvfile.close() if not quiet: print("Solution saved to '{0}'.".format(filename)) @property - def poisson_enabled(self): + def electric_field_enabled(self): """ Get/Set whether or not to solve the Poisson's equation.""" - return self.flame.poisson_enabled + return self.flame.electric_field_enabled - @poisson_enabled.setter - def poisson_enabled(self, enable): - self.flame.poisson_enabled = enable - - @property - def phi(self): - """ - Array containing the electric potential at each point. - """ - return self.profile(self.flame, 'ePotential') + @electric_field_enabled.setter + def electric_field_enabled(self, enable): + self.flame.electric_field_enabled = enable @property def E(self): """ Array containing the electric field strength at each point. """ - z = self.grid - phi = self.phi - np = self.flame.n_points - Efield = [] - Efield.append((phi[0] - phi[1]) / (z[1] - z[0])) - # calculate E field strength - for n in range(1,np-1): - Efield.append((phi[n-1] - phi[n+1]) / (z[n+1] - z[n-1])) - Efield.append((phi[np-2] - phi[np-1]) / (z[np-1] - z[np-2])) - return Efield + return self.profile(self.flame, 'eField') def solve(self, loglevel=1, refine_grid=True, auto=False, stage=1, enable_energy=True): self.flame.set_solvingStage(stage) diff --git a/interfaces/cython/cantera/onedim.pyx b/interfaces/cython/cantera/onedim.pyx index 00986155e..157a4ff19 100644 --- a/interfaces/cython/cantera/onedim.pyx +++ b/interfaces/cython/cantera/onedim.pyx @@ -549,18 +549,15 @@ cdef class IonFlow(_FlowBase): def set_solvingStage(self, stage): (self.flow).setSolvingStage(stage) - def set_electricPotential(self, v_inlet, v_outlet): - (self.flow).setElectricPotential(v_inlet, v_outlet) - - property poisson_enabled: + property electric_field_enabled: """ Determines whether or not to solve the energy equation.""" def __get__(self): - return (self.flow).doPoisson(0) + return (self.flow).doElectricField(0) def __set__(self, enable): if enable: - (self.flow).solvePoissonEqn() + (self.flow).solveElectricField() else: - (self.flow).fixElectricPotential() + (self.flow).fixElectricField() cdef class Sim1D: diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py index 0f9945dbd..d9258bb09 100644 --- a/interfaces/cython/cantera/test/test_onedim.py +++ b/interfaces/cython/cantera/test/test_onedim.py @@ -961,7 +961,7 @@ class TestIonFreeFlame(utilities.CanteraTest): self.sim.solve(loglevel=0, stage=2, enable_energy=True) # Regression test - self.assertNear(max(self.sim.E), 132.1922, 1e-3) + self.assertNear(max(self.sim.E), 136.1234, 1e-3) class TestIonBurnerFlame(utilities.CanteraTest): @@ -986,4 +986,4 @@ class TestIonBurnerFlame(utilities.CanteraTest): self.sim.solve(loglevel=0, stage=2, enable_energy=True) # Regression test - self.assertNear(max(self.sim.E), 469.7287, 1e-3) + self.assertNear(max(self.sim.E), 452.9473, 1e-3) diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp index 35837c15b..860b2cb09 100644 --- a/src/oneD/IonFlow.cpp +++ b/src/oneD/IonFlow.cpp @@ -19,8 +19,6 @@ IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) : StFlow(ph, nsp, points), m_import_electron_transport(false), m_stage(1), - m_inletVoltage(0.0), - m_outletVoltage(0.0), m_kElectron(npos) { // make a local copy of species charge @@ -43,23 +41,22 @@ IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) : } // no bound for electric potential - setBounds(c_offset_P, -1.0e20, 1.0e20); + 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); - m_refiner->setActive(c_offset_P, false); + m_refiner->setActive(c_offset_E, false); m_mobility.resize(m_nsp*m_points); - m_do_poisson.resize(m_points,false); + m_do_electric_field.resize(m_points,false); } void IonFlow::resize(size_t components, size_t points){ StFlow::resize(components, points); m_mobility.resize(m_nsp*m_points); m_do_species.resize(m_nsp,true); - m_do_poisson.resize(m_points,false); - m_fixedElecPoten.resize(m_points,0.0); + m_do_electric_field.resize(m_points,false); } void IonFlow::updateTransport(double* x, size_t j0, size_t j1) @@ -83,7 +80,7 @@ void IonFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1) frozenIonMethod(x,j0,j1); } if (m_stage == 2) { - poissonEqnMethod(x,j0,j1); + electricFieldMethod(x,j0,j1); } } @@ -114,7 +111,7 @@ void IonFlow::frozenIonMethod(const double* x, size_t j0, size_t j1) } } -void IonFlow::poissonEqnMethod(const double* x, size_t j0, size_t j1) +void IonFlow::electricFieldMethod(const double* x, size_t j0, size_t j1) { for (size_t j = j0; j < j1; j++) { double wtm = m_wtm[j]; @@ -162,17 +159,10 @@ void IonFlow::setSolvingStage(const size_t stage) throw CanteraError("IonFlow::updateDiffFluxes", "solution stage must be set to: " "1) frozenIonMethod, " - "2) poissonEqnMethod"); + "2) electricFieldEqnMethod"); } } -void IonFlow::setElectricPotential(const double v1, const double v2) -{ - // This method can be used when you want to add external voltage - m_inletVoltage = v1; - m_outletVoltage = v2; -} - void IonFlow::evalResidual(double* x, double* rsd, int* diag, double rdt, size_t jmin, size_t jmax) { @@ -189,70 +179,70 @@ void IonFlow::evalResidual(double* x, double* rsd, int* diag, for (size_t k : m_kCharge) { rsd[index(c_offset_Y + k, 0)] = Y(x,k,0) - Y(x,k,1); } - rsd[index(c_offset_P, j)] = m_inletVoltage - phi(x,j); - diag[index(c_offset_P, j)] = 0; + rsd[index(c_offset_E, j)] = E(x,0); + diag[index(c_offset_E, j)] = 0; } else if (j == m_points - 1) { - rsd[index(c_offset_P, j)] = m_outletVoltage - phi(x,j); - diag[index(c_offset_P, j)] = 0; + rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0; + diag[index(c_offset_E, j)] = 0; } else { //----------------------------------------------- - // Poisson's equation + // Electric field by Gauss's law // // dE/dz = e/eps_0 * sum(q_k*n_k) // // E = -dV/dz //----------------------------------------------- - rsd[index(c_offset_P, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0; - diag[index(c_offset_P, j)] = 0; + rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0; + diag[index(c_offset_E, j)] = 0; } } } -void IonFlow::solvePoissonEqn(size_t j) +void IonFlow::solveElectricField(size_t j) { bool changed = false; if (j == npos) { for (size_t i = 0; i < m_points; i++) { - if (!m_do_poisson[i]) { + if (!m_do_electric_field[i]) { changed = true; } - m_do_poisson[i] = true; + m_do_electric_field[i] = true; } } else { - if (!m_do_poisson[j]) { + if (!m_do_electric_field[j]) { changed = true; } - m_do_poisson[j] = true; + m_do_electric_field[j] = true; } m_refiner->setActive(c_offset_U, true); m_refiner->setActive(c_offset_V, true); m_refiner->setActive(c_offset_T, true); - m_refiner->setActive(c_offset_P, true); + m_refiner->setActive(c_offset_E, true); if (changed) { needJacUpdate(); } } -void IonFlow::fixElectricPotential(size_t j) +void IonFlow::fixElectricField(size_t j) { bool changed = false; if (j == npos) { for (size_t i = 0; i < m_points; i++) { - if (m_do_poisson[i]) { + if (m_do_electric_field[i]) { changed = true; } - m_do_poisson[i] = false; + m_do_electric_field[i] = false; } } else { - if (m_do_poisson[j]) { + if (m_do_electric_field[j]) { changed = true; } - m_do_poisson[j] = false; + m_do_electric_field[j] = false; } m_refiner->setActive(c_offset_U, false); m_refiner->setActive(c_offset_V, false); m_refiner->setActive(c_offset_T, false); - m_refiner->setActive(c_offset_P, false); + m_refiner->setActive(c_offset_E, false); if (changed) { needJacUpdate(); } @@ -279,14 +269,9 @@ void IonFlow::_finalize(const double* x) { StFlow::_finalize(x); - bool p = m_do_poisson[0]; - for (size_t j = 0; j < m_points; j++) { - if (!p) { - m_fixedElecPoten[j] = phi(x, j); - } - } + bool p = m_do_electric_field[0]; if (p) { - solvePoissonEqn(); + solveElectricField(); } } diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index fa4496146..551d29669 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -398,15 +398,15 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, 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)]; + rsd[index(c_offset_E, 0)] = x[index(c_offset_E, 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)]; + rsd[index(c_offset_E, j)] = x[index(c_offset_E, 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)]; + rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)]; //------------------------------------------------ // Radial momentum equation @@ -582,7 +582,7 @@ string StFlow::componentName(size_t n) const case 3: return "lambda"; case 4: - return "ePotential"; + return "eField"; default: if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) { return m_thermo->speciesName(n - c_offset_Y); @@ -602,7 +602,7 @@ size_t StFlow::componentIndex(const std::string& name) const return 2; } else if (name=="lambda") { return 3; - } else if (name == "ePotential") { + } else if (name == "eField") { return 4; } else { for (size_t n=c_offset_Y; n