research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[1D] delete unnecessary constrain and improve the code structure

Browse source
This commit is contained in:
bangshiuh 2017-08-08 16:23:59 -04:00 • committed by Ray Speth
parent e78aac7b70
commit 1be8374342
6 changed files with 42 additions and 68 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
include/cantera/oneD/IonFlow.h
View file

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

4
include/cantera/oneD/StFlow.h
View file

@ -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.

6
interfaces/cython/cantera/examples/onedim/ion_flame.py
View file

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

9
interfaces/cython/cantera/onedim.py
View file

@ -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:

76
src/oneD/IonFlow.cpp
View file

@ -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<size_t>(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;

7
src/oneD/StFlow.cpp
View file

@ -244,8 +244,13 @@ void StFlow::eval(size_t jg, doublereal* xg,
size_t j0 = std::max<size_t>(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