SootLib
Loading...
Searching...
No Matches
sootModel_QMOM.cc
Go to the documentation of this file.
1#include "sootModel_QMOM.h"
2#include "binomial.h"
3#include <algorithm> // fill
4#include <cmath> // isfinite (not nan or inf)
5
6using namespace std;
7using namespace soot;
8
9//------ LApack function
10extern "C" void dstev_(char *JOBZ, int *N, double *D, double *E, double *Z, int *LDZ, double *WORK, int *INFO);
11
24
26 nucleationModel *nucl_,
27 growthModel *grow_,
28 oxidationModel *oxid_,
29 coagulationModel *coag_) :
30 sootModel(nsoot_, nucl_, grow_, oxid_, coag_) {
31
32 if (nsoot_%2 == 1 || nsoot_ < 1)
33 throw runtime_error("Invalid number of soot moments requested: must be even number");
34
35 if (nsoot_ > 6)
36 cerr << "Warning: QMOM inversion algorithm may behave unpredictably with "
37 "8+ soot moments. Proceed with caution." << endl;
38
40
41}
42
55
57 nucleationMech Nmech,
58 growthMech Gmech,
59 oxidationMech Omech,
60 coagulationMech Cmech) :
61 sootModel(nsoot_, Nmech, Gmech, Omech, Cmech) {
62
63 if (nsoot_%2 == 1 || nsoot_ < 1)
64 throw runtime_error("Invalid number of soot moments requested: must be even number");
65
66 if (nsoot_ > 6)
67 cerr << "Warning: QMOM inversion algorithm may behave unpredictably with "
68 "8+ soot moments. Proceed with caution." << endl;
69
71}
72
84
86
87 //---------- get weights and abscissas
88
89 getWtsAbs(state.sootVar, state.wts, state.absc); // moment downselection applied
90
91 //---------- get chemical rates
92
93 double jNuc = nucl->getNucleationSootRate(state);
94 double kGrw = grow->getGrowthSootRate(state);
95 double kOxi = oxid->getOxidationSootRate(state);
96
97 //---------- nucleation terms
98
99 vector<double> nucSrcM(nsoot, 0.);
100
102 const double mNuc = state.cMin*gasSpMW[(int)gasSp::C]/Na; // mass of a nucleated particle
103 for (size_t i=0; i<nsoot; i++)
104 nucSrcM[i] = pow(mNuc, i)*jNuc; // Nuc_rate = m_nuc^r * jNuc
105 }
106
107 //---------- PAH condensation terms
108
109 vector<double> cndSrcM(nsoot, 0);
111 for (size_t k=1; k<nsoot; k++) { // loop moments k, skip moment 0
112 for (size_t i=0; i<state.absc.size(); i++) // loop abscissas
113 cndSrcM[k] += coag->getCoagulationSootRate(state, nucl->DIMER.mDimer, state.absc[i]) *
114 state.wts[i] * ((k==1) ? 1.0 : (k==2) ? state.absc[i] : pow(state.absc[i], k-1));
115 cndSrcM[k] *= k * nucl->DIMER.nDimer * nucl->DIMER.mDimer;
116 }
117 }
118
119 //---------- growth terms
120
121 vector<double> grwSrcM(nsoot, 0);
122
123 const double Acoef = M_PI*pow(6. /(M_PI*rhoSoot), twothird); // Acoef (=) kmol^2/3 / kg^2/3
124
126 for (size_t k=1; k<nsoot; k++) // Rate_M0 = 0.0 for growth by definition
127 grwSrcM[k] = kGrw*Acoef*k*Mr(k-onethird, state.wts, state.absc); // kg^k/m3*s
128
129 //---------- oxidation terms
130
131 vector<double> oxiSrcM(nsoot, 0);
132
134 for (size_t k=1; k<nsoot; k++) // Rate_M0 = 0.0 for oxidation by definition
135 oxiSrcM[k] = -kOxi*Acoef*k*Mr(k-onethird, state.wts, state.absc); // kg^k/m3*s
136
137 //---------- coagulation terms
138
139 vector<double> coaSrcM(nsoot, 0.0); // init to 0 to handle M1
140
142
143 double s;
144
145 for(int r=0; r<nsoot; r++) { // loop all moments r
146 if (r==1) continue; // rate is zero for r==1
147
148 //------ off-diagonal terms (looping half of them, with *2 incorporated)
149
150 for(int i=1; i<state.absc.size(); i++)
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++) // not entered if r==0, then s=-1
154 s += binomial_coefficient(r,k) * pow(state.absc[i],k) * pow(state.absc[j],r-k);
155 coaSrcM[r] += coag->getCoagulationSootRate(state,state.absc[i],state.absc[j]) * state.wts[i]*state.wts[j] * s;
156 }
157
158 //------ diagonal terms
159
160 s = (r==0) ? -1.0 : 0.0;
161 for(int k=1; k<=r-1; k++) // not entered if r==0, then s=-1
162 s += binomial_coefficient(r,k);
163 for(int i=0; i<state.absc.size(); i++) {
164 coaSrcM[r] += 0.5*coag->getCoagulationSootRate(state,state.absc[i],state.absc[i]) *state.wts[i]*state.wts[i] * s *
165 ((r==0) ? 1.0 : pow(state.absc[i], r));
166 }
167 }
168 }
169
170 //---------- combine to make soot source terms
171
172 for (size_t i = 0; i < nsoot; i++)
173 sources.sootSources[i] = nucSrcM[i] + cndSrcM[i] + grwSrcM[i] + oxiSrcM[i] + coaSrcM[i]; // kg^k/m3*s
174
175 //---------- set gas source terms
176
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);
180
181 nucl->getNucleationGasRates(nucSrcM[1], nucl_gasSources);
182 grow->getGrowthGasRates( grwSrcM[1], grow_gasSources);
183 oxid->getOxidationGasRates( oxiSrcM[1], oxid_gasSources);
184
185 for (size_t sp=0; sp<(size_t)gasSp::size; sp++)
186 sources.gasSources[sp] = nucl_gasSources[sp] + grow_gasSources[sp] + oxid_gasSources[sp];
187
188 //---------- set PAH source terms
189
191 sources.pahSources = nucl->nucleationPahRxnRates; // includes both nucleation and condensation
192}
193
207
208void sootModel_QMOM::getWtsAbs(const vector<double>& M, vector<double>& weights, vector<double>& abscissas) {
209
210 fill(weights.begin(), weights.end(), 0.0); // initialize weights and abscissas
211 fill(abscissas.begin(), abscissas.end(), 0.0);
212
213 //----------
214
215 size_t Nenv = M.size()/2; // local nsoonsoot; may change with moment downselection
216
217 vector<double> w_temp(Nenv, 0);
218 vector<double> a_temp(Nenv, 0);
219
220 bool bad_values; // downselection flag
221 do { // downselection loop
222
223 bad_values = false; // reset flag
224
225 //---- 1 env (2 moment) case, return MONO output
226 if (Nenv == 1) {
227 //cout << endl << "downselected to mono" << endl;
228 w_temp[0] = M[0];
229 a_temp[0] = M[1] / M[0];
230 break;
231 }
232
233 //---- wheeler algorithm for moment inversion
234 wheeler(M, Nenv, w_temp, a_temp);
235
236 //---- check for bad values
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) )
240 bad_values = true;
241 }
242
243 //---- if we find bad values, downselect to 2 fewer moments (1 fewer env) and try again
244 if (bad_values) {
245 Nenv--;
246 w_temp.assign(Nenv, 0.0); // resizes and assigns 0.0
247 a_temp.assign(Nenv, 0.0);
248 }
249
250 }
251 while (bad_values);
252
253 //---------- assign temporary variables to output
254 for (size_t i=0; i<Nenv; i++) {
255 weights[i] = w_temp[i] >= 0 ? w_temp[i] : 0.0; // these checks should be redundant
256 abscissas[i] = a_temp[i] >= 0 ? a_temp[i] : 0.0;
257 }
258}
259
275
276void sootModel_QMOM::wheeler(const vector<double>& m, size_t N, vector<double>& w, vector<double>& x) {
277
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);
285
286 for (size_t i = 0; i <= N * 2 - 1; i++)
287 sigma[1][i] = m[i];
288
289 a[0] = m[1] / m[0];
290
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];
296 }
297
298 j_diag = a;
299 for (size_t i = 1; i < N; i++)
300 j_ldiag[i] = -sqrt(abs(b[i]));
301
302 for (size_t i = 0; i < N; i++)
303 evec[i + N * i] = 1;
304
305 //----- compute eigenvalues/eigenvectors (lapack)
306
307 char VorN = 'V';
308 vector<double> work(2*N-2);
309 int info;
310 int NN = int(N);
311 dstev_( &VorN, &NN, &j_diag[0], &j_ldiag[1], &evec[0], &NN, &work[0], &info);
312
313 //-----
314
315 x = j_diag; // j_diag are now the vector of eigenvalues.
316
317 for (size_t i = 0; i < N; i++)
318 w[i] = pow(evec[0 + i * N], 2) * m[0];
319}
320
331
332double sootModel_QMOM::Mr(double r, const vector<double>& wts, const vector<double>& absc) {
333 double Mr = 0;
334
335 for (size_t i=0; i<wts.size(); i++) {
336 if (wts[i] <= 0 || absc[i] <= 0)
337 return 0;
338 else
339 Mr += wts[i] * pow(absc[i], r);
340 }
341 return Mr;
342}
double binomial_coefficient(unsigned r, unsigned k)
Definition: binomial.h:24
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
Definition: growthModel.cc:16
virtual double getGrowthSootRate(const state &state) const =0
growthMech mechType
identity of the type of growth (child)
Definition: growthModel.h:24
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
Definition: sootModel.h:33
psdMech psdMechType
one of MONO, LOGN, QMOM, MOMIC, SECT, etc.
Definition: sootModel.h:37
sourceTerms sources
struct containing soot, gas, and pah source terms vectors
Definition: sootModel.h:42
nucleationModel * nucl
pointer to nucleation mechanism
Definition: sootModel.h:30
growthModel * grow
pointer to growth mechanism
Definition: sootModel.h:31
size_t nsoot
# of soot variables: moments or sections
Definition: sootModel.h:28
oxidationModel * oxid
pointer to oxidation mechanism
Definition: sootModel.h:32
std::vector< double > wts
moment weights
Definition: state.h:34
std::vector< double > sootVar
soot variables (moments or # in sections>
Definition: state.h:32
double cMin
soot min num carbon atoms (dynamic for PAH nucleation)
Definition: state.h:36
std::vector< double > absc
moment abscissas
Definition: state.h:33
Definition: sootDefs.h:11
oxidationMech
Definition: sootDefs.h:33
const double rhoSoot
soot particle density
Definition: sootDefs.h:21
const double onethird
Definition: sootDefs.h:24
coagulationMech
Definition: sootDefs.h:34
growthMech
Definition: sootDefs.h:32
nucleationMech
Definition: sootDefs.h:31
const double twothird
Definition: sootDefs.h:25
const double Na
Avogadro's constant: #/kmol.
Definition: sootDefs.h:15
const std::vector< double > gasSpMW
(kg/kmol); make sure the order corresponds to the gasSp enum
Definition: sootDefs.h:74
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
Definition: sootDefs.h:143
std::vector< double > gasSources
kg/m3*s
Definition: sootDefs.h:142
std::vector< double > sootSources
kg^r/m3*s (moments), or #/m3*s (sections)
Definition: sootDefs.h:141