From 1be837434240a8610f164ee14355d4f0eee27057 Mon Sep 17 00:00:00 2001 From: bangshiuh Date: Tue, 8 Aug 2017 16:23:59 -0400 Subject: [PATCH] [1D] delete unnecessary constrain and improve the code structure --- include/cantera/oneD/IonFlow.h | 8 +- include/cantera/oneD/StFlow.h | 4 + .../cantera/examples/onedim/ion_flame.py | 6 +- interfaces/cython/cantera/onedim.py | 9 ++- src/oneD/IonFlow.cpp | 76 +++++-------------- src/oneD/StFlow.cpp | 7 +- 6 files changed, 42 insertions(+), 68 deletions(-) diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h index af7a546ce..38eae10e0 100644 --- a/include/cantera/oneD/IonFlow.h +++ b/include/cantera/oneD/IonFlow.h @@ -46,9 +46,6 @@ public: //! set electric voltage at inlet and outlet virtual void setElectricPotential(const double v1, const double v2); - virtual void eval(size_t jg, double* xg, - double* rg, integer* diagg, double rdt); - virtual void resize(size_t components, size_t points); virtual void _finalize(const double* x); @@ -82,10 +79,11 @@ public: vector_fp& mobi_e_fixed); protected: + virtual void updateProperties(size_t jg, double* x, double* rsd, + int* diag, double rdt, size_t j0, + size_t j1, size_t jmin, size_t jmax); virtual void updateTransport(double* x, size_t j0, size_t j1); virtual void updateDiffFluxes(const double* x, size_t j0, size_t j1); - //! evaluate the residual for Poisson's equation - virtual void evalPoisson(size_t j, double* x, double* r, integer* diag, double rdt); //! 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 two: the Prager's ambipolar-diffusion model is used diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h index 0b0d76b7f..600bd3e86 100644 --- a/include/cantera/oneD/StFlow.h +++ b/include/cantera/oneD/StFlow.h @@ -242,6 +242,10 @@ protected: m_kin->getNetProductionRates(&m_wdot(0,j)); } + virtual void updateProperties(size_t jg, double* x, double* rsd, + int* diag, double rdt, size_t j0, + size_t j1, size_t jmin, size_t jmax); + /** * Update the thermodynamic properties from point j0 to point j1 * (inclusive), based on solution x. diff --git a/interfaces/cython/cantera/examples/onedim/ion_flame.py b/interfaces/cython/cantera/examples/onedim/ion_flame.py index 44e463c07..c0fb170f6 100644 --- a/interfaces/cython/cantera/examples/onedim/ion_flame.py +++ b/interfaces/cython/cantera/examples/onedim/ion_flame.py @@ -22,14 +22,14 @@ f = ct.IonFlame(gas, width=width) f.set_refine_criteria(ratio=3, slope=0.06, curve=0.12) f.show_solution() -# phase one +# stage one f.solve(loglevel=loglevel, auto=True) -# phase two +# stage two f.solve(loglevel=loglevel, stage=2, enable_energy=False) f.solve(loglevel=loglevel, stage=2, enable_energy=True) -# phase three +# stage three f.solve(loglevel=loglevel, stage=3, enable_energy=True) f.save('CH4_adiabatic.xml', 'mix', 'solution with mixture-averaged transport') diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index c8a135434..2b4bc7028 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -502,7 +502,10 @@ class IonFlame(FreeFlame): __slots__ = ('inlet', 'outlet', 'flame') def __init__(self, gas, grid=None, width=None): - self.flame = IonFlow(gas, name='flame') + if not hasattr(self, 'flame'): + # Create flame domain if not already instantiated by a child class + self.flame = IonFlow(gas, name='flame') + super(IonFlame, self).__init__(gas, grid, width) def solve(self, loglevel=1, refine_grid=True, auto=False, stage=1, enable_energy=True): @@ -544,10 +547,10 @@ class IonFlame(FreeFlame): 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) + 'phi (V)', 'E (V/m)', 'rho (kmol/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], phi[n], E[n], self.gas.density_mole] + list(getattr(self.gas, species))) csvfile.close() if not quiet: diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp index 812478fea..d7328ce9f 100644 --- a/src/oneD/IonFlow.cpp +++ b/src/oneD/IonFlow.cpp @@ -40,24 +40,8 @@ IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) : // Find the index of electron if (m_thermo->speciesIndex("E") != npos ) { m_kElectron = m_thermo->speciesIndex("E"); - setTransientTolerances(1.0e-5, 1.0e-18, c_offset_Y + m_kElectron); - setSteadyTolerances(1.0e-5, 1.0e-16, c_offset_Y + m_kElectron); - } - if (m_thermo->speciesIndex("HCO+") != npos ) { - size_t k = m_thermo->speciesIndex("HCO+"); - setTransientTolerances(1.0e-5, 1.0e-18, c_offset_Y + k); - setSteadyTolerances(1.0e-5, 1.0e-16, c_offset_Y + k); - } - if (m_thermo->speciesIndex("H3O+") != npos ) { - size_t k = m_thermo->speciesIndex("H3O+"); - setTransientTolerances(1.0e-5, 1.0e-15, c_offset_Y + k); - setSteadyTolerances(1.0e-5, 1.0e-13, c_offset_Y + k); } - // mass fraction bounds (strict bound for ions) - for (size_t k : m_kCharge) { - setBounds(c_offset_Y+k, -1.0e-20, 1.0e5); - } // no bound for electric potential setBounds(c_offset_P, -1.0e20, 1.0e20); @@ -246,43 +230,40 @@ void IonFlow::setElectricPotential(const double v1, const double v2) m_outletVoltage = v2; } -void IonFlow::eval(size_t jg, double* xg, - double* rg, integer* diagg, double rdt) +void IonFlow::updateProperties(size_t jg, double* x, double* rsd, + int* diag, double rdt, size_t j0, + size_t j1, size_t jmin, size_t jmax) { - StFlow::eval(jg, xg, rg, diagg, rdt); + StFlow::updateProperties(jg, x, rsd, diag, rdt, j0, j1, jmin, jmax); if (m_stage != 3) { return; } - // start of local part of global arrays - double* x = xg + loc(); - double* 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(jpt, 1) - 1; - jmax = std::min(jpt+1,m_points-1); - } for (size_t j = jmin; j <= jmax; j++) { if (j == 0) { rsd[index(c_offset_P, j)] = m_inletVoltage - phi(x,j); diag[index(c_offset_P, j)] = 0; - // set ions boundary for better convergence - for (size_t k : m_kCharge) { - rsd[index(c_offset_Y + k, j)] = Y(x,k,j+1) - Y(x,k,j); - } } else if (j == m_points - 1) { rsd[index(c_offset_P, j)] = m_outletVoltage - phi(x,j); diag[index(c_offset_P, j)] = 0; } else { - evalPoisson(j,x,rsd,diag,rdt); + //----------------------------------------------- + // Poisson's equation + // + // dE/dz = e/eps_0 * sum(q_k*n_k) + // + // E = -dV/dz + //----------------------------------------------- + double 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; + + // This method is used when you disable energy equation + // but still maintain the velocity profile if (!m_do_velocity[j]) { - // This method is used when you disable energy equation - // but still maintain the velocity profile rsd[index(c_offset_U, j)] = u(x,j) - u_fixed(j); diag[index(c_offset_U, j)] = 0; } @@ -290,23 +271,6 @@ void IonFlow::eval(size_t jg, double* xg, } } -void IonFlow::evalPoisson(size_t j, double* x, double* rsd, integer* diag, double rdt) -{ - //----------------------------------------------- - // Poisson's equation - // - // dE/dz = e/eps_0 * sum(q_k*n_k) - // - // E = -dV/dz - //----------------------------------------------- - double 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::solvePoissonEqn(size_t j) { bool changed = false; diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index d667822c7..660fe8075 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -244,8 +244,13 @@ void StFlow::eval(size_t jg, doublereal* xg, size_t j0 = std::max(jmin, 1) - 1; size_t j1 = std::min(jmax+1,m_points-1); - // ------------ update properties ------------ + updateProperties(jg, x, rsd, diag, rdt, j0, j1, jmin, jmax); +} +void StFlow::updateProperties(size_t jg, double* x, double* rsd, + int* diag, double rdt, size_t j0, + size_t j1, size_t jmin, size_t jmax) +{ updateThermo(x, j0, j1); if (jg == npos || m_force_full_update) { // update transport properties only if a Jacobian is not being