10extern "C" void dstev_(
char *JOBZ,
int *N,
double *D,
double *E,
double *Z,
int *LDZ,
double *WORK,
int *INFO);
30 sootModel(nsoot_, nucl_, grow_, oxid_, coag_) {
32 if (nsoot_%2 == 1 || nsoot_ < 1)
33 throw runtime_error(
"Invalid number of soot moments requested: must be even number");
36 cerr <<
"Warning: QMOM inversion algorithm may behave unpredictably with "
37 "8+ soot moments. Proceed with caution." << endl;
61 sootModel(nsoot_, Nmech, Gmech, Omech, Cmech) {
63 if (nsoot_%2 == 1 || nsoot_ < 1)
64 throw runtime_error(
"Invalid number of soot moments requested: must be even number");
67 cerr <<
"Warning: QMOM inversion algorithm may behave unpredictably with "
68 "8+ soot moments. Proceed with caution." << endl;
99 vector<double> nucSrcM(
nsoot, 0.);
103 for (
size_t i=0; i<
nsoot; i++)
104 nucSrcM[i] = pow(mNuc, i)*jNuc;
109 vector<double> cndSrcM(
nsoot, 0);
111 for (
size_t k=1; k<
nsoot; k++) {
112 for (
size_t i=0; i<
state.
absc.size(); i++)
121 vector<double> grwSrcM(
nsoot, 0);
126 for (
size_t k=1; k<
nsoot; k++)
131 vector<double> oxiSrcM(
nsoot, 0);
134 for (
size_t k=1; k<
nsoot; k++)
139 vector<double> coaSrcM(
nsoot, 0.0);
145 for(
int r=0; r<
nsoot; r++) {
151 for(
int j=0; j<i; j++) {
152 s = (r==0) ? -1.0 : 0.0;
153 for(
int k=1; k<=r-1; k++)
160 s = (r==0) ? -1.0 : 0.0;
161 for(
int k=1; k<=r-1; k++)
172 for (
size_t i = 0; i <
nsoot; i++)
173 sources.
sootSources[i] = nucSrcM[i] + cndSrcM[i] + grwSrcM[i] + oxiSrcM[i] + coaSrcM[i];
177 vector<double> nucl_gasSources((
size_t)
gasSp::size, 0.0);
178 vector<double> grow_gasSources((
size_t)
gasSp::size, 0.0);
179 vector<double> oxid_gasSources((
size_t)
gasSp::size, 0.0);
186 sources.
gasSources[sp] = nucl_gasSources[sp] + grow_gasSources[sp] + oxid_gasSources[sp];
210 fill(weights.begin(), weights.end(), 0.0);
211 fill(abscissas.begin(), abscissas.end(), 0.0);
215 size_t Nenv = M.size()/2;
217 vector<double> w_temp(Nenv, 0);
218 vector<double> a_temp(Nenv, 0);
229 a_temp[0] = M[1] / M[0];
234 wheeler(M, Nenv, w_temp, a_temp);
237 for (
size_t i=0; i<Nenv; i++) {
238 if ( (!isfinite(w_temp[i]) || w_temp[i]<0) ||
239 (!isfinite(a_temp[i]) || a_temp[i]<0) )
246 w_temp.assign(Nenv, 0.0);
247 a_temp.assign(Nenv, 0.0);
254 for (
size_t i=0; i<Nenv; i++) {
255 weights[i] = w_temp[i] >= 0 ? w_temp[i] : 0.0;
256 abscissas[i] = a_temp[i] >= 0 ? a_temp[i] : 0.0;
278 vector<vector<double>> sigma(N + 1, vector<double>(N * 2, 0));
279 vector<double> a(N, 0);
280 vector<double> b(N, 0);
281 vector<double> eval(N);
282 vector<double> evec(N * N);
283 vector<double> j_diag(N);
284 vector<double> j_ldiag(N);
286 for (
size_t i = 0; i <= N * 2 - 1; i++)
291 for (
size_t i = 1; i < N; i++) {
292 for (
size_t j = i; j < N * 2 - i; j++)
293 sigma[i + 1][j] = sigma[i][j + 1] - a[i - 1] * sigma[i][j] - b[i - 1] * sigma[i - 1][j];
294 a[i] = -sigma[i][i] / sigma[i][i - 1] + sigma[i + 1][i + 1] / sigma[i + 1][i];
295 b[i] = sigma[i + 1][i] / sigma[i][i - 1];
299 for (
size_t i = 1; i < N; i++)
300 j_ldiag[i] = -sqrt(abs(b[i]));
302 for (
size_t i = 0; i < N; i++)
308 vector<double> work(2*N-2);
311 dstev_( &VorN, &NN, &j_diag[0], &j_ldiag[1], &evec[0], &NN, &work[0], &info);
317 for (
size_t i = 0; i < N; i++)
318 w[i] = pow(evec[0 + i * N], 2) * m[0];
335 for (
size_t i=0; i<wts.size(); i++) {
336 if (wts[i] <= 0 || absc[i] <= 0)
339 Mr += wts[i] * pow(absc[i], r);
double binomial_coefficient(unsigned r, unsigned k)
coagulationMech mechType
identity of the type of coagulation (child)
virtual double getCoagulationSootRate(const state &state, double m1, double m2) const =0
void getGrowthGasRates(const double &msootDotGrow, std::vector< double > &gasSourcesGrow) const
virtual double getGrowthSootRate(const state &state) const =0
growthMech mechType
identity of the type of growth (child)
std::vector< double > nucleationPahRxnRates
mole ratios for PAH gas species rate coupling
void getNucleationGasRates(const double &msootDotNucl, std::vector< double > &gasSourcesNucl) const
dimerStruct DIMER
used for PAH nucleation only
nucleationMech mechType
identity of the type of nucleation (child)
virtual double getNucleationSootRate(state &state)=0
void getOxidationGasRates(const double &msootDotOxid, std::vector< double > &gasSourcesOxid) const
oxidationMech mechType
identity of the type of oxidation (child)
virtual double getOxidationSootRate(const state &state) const =0
static void wheeler(const std::vector< double > &m, size_t N, std::vector< double > &w, std::vector< double > &x)
static double Mr(double r, const std::vector< double > &wts, const std::vector< double > &absc)
virtual void setSourceTerms(state &state)
static void getWtsAbs(const std::vector< double > &M, std::vector< double > &weights, std::vector< double > &abscissas)
sootModel_QMOM(size_t nsoot_, nucleationModel *nucl_, growthModel *grow_, oxidationModel *oxid_, coagulationModel *coag_)
coagulationModel * coag
pointer to coagulation mechanism
psdMech psdMechType
one of MONO, LOGN, QMOM, MOMIC, SECT, etc.
sourceTerms sources
struct containing soot, gas, and pah source terms vectors
nucleationModel * nucl
pointer to nucleation mechanism
growthModel * grow
pointer to growth mechanism
size_t nsoot
# of soot variables: moments or sections
oxidationModel * oxid
pointer to oxidation mechanism
std::vector< double > wts
moment weights
std::vector< double > sootVar
soot variables (moments or # in sections>
double cMin
soot min num carbon atoms (dynamic for PAH nucleation)
std::vector< double > absc
moment abscissas
const double rhoSoot
soot particle density
const double Na
Avogadro's constant: #/kmol.
const std::vector< double > gasSpMW
(kg/kmol); make sure the order corresponds to the gasSp enum
void dstev_(char *JOBZ, int *N, double *D, double *E, double *Z, int *LDZ, double *WORK, int *INFO)
std::vector< double > pahSources
kg/m3*s
std::vector< double > gasSources
kg/m3*s
std::vector< double > sootSources
kg^r/m3*s (moments), or #/m3*s (sections)