18 IonFlow::IonFlow(IdealGasPhase* ph,
size_t nsp,
size_t points) :
19 StFlow(ph, nsp, points),
20 m_import_electron_transport(false),
25 for (
size_t k = 0; k < m_nsp; k++) {
26 m_speciesCharge.push_back(
m_thermo->charge(k));
30 for (
size_t k = 0; k < m_nsp; k++){
31 if (m_speciesCharge[k] != 0){
32 m_kCharge.push_back(k);
34 m_kNeutral.push_back(k);
40 m_kElectron =
m_thermo->speciesIndex(
"E");
44 setBounds(c_offset_E, -1.0e20, 1.0e20);
49 for (
size_t k : m_kCharge) {
50 setBounds(c_offset_Y + k, -1e-14, 1.0);
52 setBounds(c_offset_Y + m_kElectron, -1e-18, 1.0);
54 m_refiner->setActive(c_offset_E,
false);
55 m_mobility.resize(m_nsp*m_points);
56 m_do_electric_field.resize(m_points,
false);
62 m_do_species.resize(m_nsp,
true);
69 for (
size_t j = j0; j < j1; j++) {
76 m_diff[k+m_nsp*j] =
poly5(tlog, m_diff_e_fix.data());
93 for (
size_t j = j0; j < j1; j++) {
94 double wtm = m_wtm[j];
95 double rho = density(j);
96 double dz = z(j+1) - z(j);
99 m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
100 m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
106 m_flux(k,j) += sum*Y(x,k,j);
120 for (
size_t j = j0; j < j1; j++) {
121 double wtm = m_wtm[j];
122 double rho = density(j);
123 double dz = z(j+1) - z(j);
127 for (
size_t k = 0; k < m_nsp; k++) {
128 m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
129 m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
134 double E_ambi =
E(x,j);
136 double Yav = 0.5 * (Y(x,k,j) + Y(x,k,j+1));
137 double drift = rho * Yav * E_ambi
139 m_flux(k,j) += drift;
143 double sum_flux = 0.0;
144 for (
size_t k = 0; k < m_nsp; k++) {
145 sum_flux -= m_flux(k,j);
147 double sum_ion = 0.0;
153 m_flux(k,j) += Y(x,k,j) / (1-sum_ion) * sum_flux;
160 if (stage == 1 || stage == 2) {
164 "solution stage must be set to: " 165 "1) frozenIonMethod, " 166 "2) electricFieldEqnMethod");
171 double rdt,
size_t jmin,
size_t jmax)
178 for (
size_t j = jmin; j <= jmax; j++) {
184 rsd[index(c_offset_Y + k, 0)] = Y(x,k,0) - Y(x,k,1);
186 rsd[index(c_offset_E, j)] =
E(x,0);
187 diag[index(c_offset_E, j)] = 0;
188 }
else if (j == m_points - 1) {
190 diag[index(c_offset_E, j)] = 0;
200 diag[index(c_offset_E, j)] = 0;
207 bool changed =
false;
209 for (
size_t i = 0; i < m_points; i++) {
221 m_refiner->setActive(c_offset_U,
true);
222 m_refiner->setActive(c_offset_V,
true);
223 m_refiner->setActive(c_offset_T,
true);
224 m_refiner->setActive(c_offset_E,
true);
232 bool changed =
false;
234 for (
size_t i = 0; i < m_points; i++) {
246 m_refiner->setActive(c_offset_U,
false);
247 m_refiner->setActive(c_offset_V,
false);
248 m_refiner->setActive(c_offset_T,
false);
249 m_refiner->setActive(c_offset_E,
false);
260 size_t n = tfix.size();
262 for (
size_t i = 0; i < n; i++) {
263 tlog.push_back(log(tfix[i]));
266 m_diff_e_fix.resize(degree + 1);
268 polyfit(n, degree, tlog.data(), diff_e.data(), w.data(), m_diff_e_fix.data());
void solveElectricField(size_t j=npos)
set to solve electric field on a point
virtual void resize(size_t components, size_t points)
Change the grid size. Called after grid refinement.
vector_fp m_mobility
mobility
std::vector< thermo_t * > m_thermo
m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator...
virtual void updateTransport(double *x, size_t j0, size_t j1)
Update the transport properties at grid points in the range from j0 to j1, based on solution x...
bool m_import_electron_transport
flag for importing transport of electron
size_t m_kElectron
index of electron
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
doublereal temperature() const
Temperature (K).
const size_t npos
index returned by functions to indicate "no position"
R poly5(D x, R *c)
Templated evaluation of a polynomial of order 5.
virtual void _finalize(const doublereal *x)
In some cases, a domain may need to set parameters that depend on the initial solution estimate...
Headers for the Transport object, which is the virtual base class for all transport property evaluato...
double E(const double *x, size_t j) const
electric field
std::vector< bool > m_do_electric_field
flag for solving electric field or not
virtual void evalResidual(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the residual function.
double rho_e(double *x, size_t j) const
total charge density
std::vector< size_t > m_kNeutral
index of neutral species
virtual void evalResidual(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
vector_int m_speciesCharge
electrical properties
virtual void _finalize(const double *x)
In some cases, a domain may need to set parameters that depend on the initial solution estimate...
virtual void updateDiffFluxes(const double *x, size_t j0, size_t j1)
Update the diffusive mass fluxes.
virtual void setSolvingStage(const size_t phase)
set the solving stage
Base class for exceptions thrown by Cantera classes.
virtual void getMobilities(doublereal *const mobil_e)
Get the Electrical mobilities (m^2/V/s).
virtual void electricFieldMethod(const double *x, size_t j0, size_t j1)
Solving phase two: the electric field equation is added coupled by the electrical drift...
virtual void updateTransport(doublereal *x, size_t j0, size_t j1)
Update the transport properties at grid points in the range from j0 to j1, based on solution x...
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
void setGasAtMidpoint(const doublereal *x, size_t j)
Set the gas state to be consistent with the solution at the midpoint between j and j + 1...
std::vector< size_t > m_kCharge
index of species with charges
const doublereal epsilon_0
Permittivity of free space in F/m.
void fixElectricField(size_t j=npos)
set to fix voltage on a point
vector_fp m_mobi_e_fix
coefficients of polynomial fitting of fixed electron transport profile
virtual void resize(size_t components, size_t points)
Change the grid size. Called after grid refinement.
Namespace for the Cantera kernel.
double polyfit(size_t n, size_t deg, const double *xp, const double *yp, const double *wp, double *pp)
Fits a polynomial function to a set of data points.
Header for a file containing miscellaneous numerical functions.
void setElectronTransport(vector_fp &tfix, vector_fp &diff_e, vector_fp &mobi_e)
Sometimes it is desired to carry out the simulation using a specified electron transport profile...
virtual void frozenIonMethod(const double *x, size_t j0, size_t j1)
Solving phase one: the fluxes of charged species are turned off.