Ignis
Loading...
Searching...
No Matches
ignis.cc
Go to the documentation of this file.
1
2#include "ignis.h"
3#include "cantera/base/ct_defs.h"
4#include "solver_kinsol.h"
5#include "integrator_cvode.h"
6
7#include <iostream>
8#include <algorithm> // max
9#include <iomanip>
10#include <fstream>
11#include <sstream>
12#include <numeric>
13#include <cmath>
14#include <highfive/highfive.hpp>
15
16
17using namespace std;
18using soot::sootModel, soot::state, soot::gasSp, soot::gasSpMapIS, soot::gasSpMapES;
19using HighFive::File, HighFive::Group, HighFive::DataSet;
20
22
23int Func_kinsol(N_Vector varsKS, N_Vector fvec, void *user_data);
24int rhsf_cvode(realtype t, N_Vector varsCV, N_Vector dvarsdtCV, void *user_data);
25
44
45ignis::ignis(const bool _isPremixed,
46 const bool _doEnergyEqn,
47 const bool _isFlamelet,
48 const bool _doSoot,
49 const size_t _ngrd, const double _L, double _P,
50 shared_ptr<Cantera::Solution> csol,
51 string _radType,
52 const vector<double> &_yLbc, const vector<double> &_yRbc,
53 const double _TLbc, const double _TRbc,
54 shared_ptr<sootModel> _SM,
55 shared_ptr<state> _SMstate) :
56 isPremixed(_isPremixed),
57 isFlamelet(_isFlamelet),
58 doEnergyEqn(_doEnergyEqn),
59 doSoot(_doSoot),
60 ngrd(_ngrd),
61 L(_L),
62 P(_P),
63 yLbc(_yLbc),
64 yRbc(_yRbc),
65 TLbc(_TLbc),
66 TRbc(_TRbc),
67 SM(_SM),
68 SMstate(_SMstate),
69 radType(_radType) {
70
71 //----------
72
73 if(_isPremixed && _isFlamelet)
74 throw runtime_error("Cannot set both premixed and flamelet");
75 if(_isFlamelet && _L!=1.0)
76 throw runtime_error("flamelet should have L=1");
77
78 //----------
79
80 gas = csol->thermo();
81 kin = csol->kinetics();
82 trn = csol->transport();
83
84 nsp = gas->nSpecies();
85 nsoot = (doSoot ? SM->nsoot : 0);
86 nvar = nsp + nsoot + 1;
87 nvarA = nvar*ngrd;
88
89 y = vector<vector<double> >(ngrd, vector<double>(nsp, 0.0));
90 if(doSoot) sootvars = vector<vector<double> >(ngrd, vector<double>(nsoot, 0.0));
91 T.resize(ngrd, 0.0);
92
93 //---------- set grid
94
95 setGrid(L);
96
97 //----------
98
99 gas->setState_TPY(TLbc, P, &yLbc[0]);
100 hLbc = gas->enthalpy_mass();
101 cpLbc = gas->cp_mass();
102 hspLbc.resize(nsp);
103 gas->getEnthalpy_RT(&hspLbc[0]);
104 for(size_t k=0; k<nsp; k++) // --> hsp = J/kg species i
105 hspLbc[k] *= TLbc*Cantera::GasConstant/gas->molecularWeight(k);
106
107 gas->setState_TPY(TRbc, P, &yRbc[0]);
108 hRbc = gas->enthalpy_mass();
109 cpRbc = gas->cp_mass();
110 hspRbc.resize(nsp);
111 gas->getEnthalpy_RT(&hspRbc[0]);
112 for(size_t k=0; k<nsp; k++) // --> hsp = J/kg species i
113 hspRbc[k] *= TRbc*Cantera::GasConstant/gas->molecularWeight(k);
114
115 if(!isPremixed) strm = make_shared<streams>(csol, P, hLbc, hRbc, yLbc, yRbc);
116
117 hscale = max(abs(hLbc), abs(hRbc));
118 Tscale = 2500;
119
120 //---------- set sootScales, and species map
121 // Mk = integral (m^k * n(m) *dm) --> Mk ~ M0<m>^k, where <m> in an average soot size
122 // <m> = M1/M0 = rho*Ys/M0 = fv*rhos/M0
123 // Mk ~ M0(rhos*fv/M0)^k; take M0=1E10 #/m3, rhos=1850 kg/m3, fv=1E-6
124
125 if(doSoot) {
126 sootScales = vector<double>(nsoot, 1E20);
127 for(int i=1; i<nsoot; i++)
128 sootScales[i] = sootScales[i-1]*(1850*1E-6/sootScales[0]);
129 sootScales = {1E20, 1.8E-3, 3.4E-26, 6.3E-49, 1.2E-71, 2.2E-94, 4E-117, 7.4E-140}; // (same as loop above, rounded)
130 }
131 //----------
132
133 flux_y = vector<vector<double> >(ngrd+1, vector<double>(nsp, 0.0));
134 if(doSoot) flux_soot = vector<vector<double> >(ngrd+1, vector<double>(nsoot, 0.0));
135 flux_h.resize(ngrd+1);
136
137 //---------- radiation object
138
139 doRadiation = false;
140
141 if(radType == "planckmean") {
142 radProps = make_shared<rad_planck_mean>();
143 }
144 else if(radType == "wsgg") {
145 radProps = make_shared<rad_wsgg>();
146 }
147 else if(radType == "rcslw") {
148 double fvsoot = 0.0;
149 int isp;
150 isp = gas->speciesIndex("H2O");
151 double xH2O = yLbc[isp];
152 isp = gas->speciesIndex("CO2");
153 double xCO2 = yLbc[isp];
154 isp = gas->speciesIndex("CO");
155 double xCO = yLbc[isp];
156 radProps = make_shared<rad_rcslw>(4, TLbc, P, fvsoot, xH2O, xCO2, xCO);
157 }
158 else
159 radProps = make_shared<rad_planck_mean>();
160
161
162 //----------- Get the absorption & weights for the surroundings
163
164 vector<double> kabs_sur, awts_sur; // surroundings kabs and awts
165 double fvsoot = 0.0;
166 double xH2O, xCO2, xCO, xCH4;
167 int isp;
168
169 isp = gas->speciesIndex("H2O");
170 xH2O = yLbc[isp];
171 isp = gas->speciesIndex("CO2");
172 xCO2 = yLbc[isp];
173 isp = gas->speciesIndex("CO");
174 xCO = yLbc[isp];
175 isp = gas->speciesIndex("CH4");
176 xCH4 = yLbc[isp];
177
178 radProps->get_k_a(kabs_sur, awts_sur, TLbc, P, fvsoot, xH2O, xCO2, xCO, xCH4);
179
180
181
182 // set even if doRadiation is false, since we switch it on/off for some cases
183
184 Ttarget = 0.0;
185
186 //---------- hdf5 file
187
188 fh5 = make_shared<File>("data.h5", File::Truncate);
189
190}
191
201
202void ignis::setGrid(double _L) {
203
204 L = _L;
205 dx = vector<double>(ngrd, L/ngrd);
206 x.resize(ngrd);
207
208 if(isPremixed) { //----------- uniform grid
209 x[0] = dx[0]/2;
210 for (size_t i=1; i<ngrd; i++)
211 x[i] = x[i-1] + (dx[i-1]+dx[i])/2;
212 }
213 else { //--------- segmented grid
214
215 double Lfrac = 0.5; // first Lfrac fraction of the domain length
216 double Gfrac = 0.5; // gets this Gfrac fraction of the grid points
217 int n1 = ngrd*Gfrac;
218 double dx1 = L*Lfrac/n1;
219 for(int i=0; i<n1; i++)
220 dx[i] = dx1;
221
222 int n2 = ngrd-n1;
223 double dx2 = L*(1.0-Lfrac)/n2;
224 for(int i=n1; i<ngrd; i++)
225 dx[i] = dx2;
226
227 x[0] = dx[0]/2;
228 for (size_t i=1; i<ngrd; i++)
229 x[i] = x[i-1] + (dx[i-1]+dx[i])/2;
230 }
231
232 //--------- set fx
233 // interpolation factors to interior faces
234 // 0 1 2 3 4
235 // | * | * | * | * |
236 // 0 1 2 3
237 // vf[1] = v[0]*(dx[1]/(dx[0]+dx[1])) + v[1]*(dx[0]/(dx[0]+dx[1]))
238 // = v[0]*( fl[0] ) + v[1]*( 1.0 - fl[0] )
239 // = v[0]*( fl[0] ) + v[1]*( fr[0] )
240 //
241 // vf[3] = v[2]*(dx[3]/(dx[2]+dx[3])) + v[3]*(dx[2]/(dx[2]+dx[3]))
242 // = v[2]*( fl[2] ) + v[3]*( 1.0 - fl[2] )
243 // = v[2]*( fl[2] ) + v[3]*( fr[2] )
244 //
245 // vf[i] = v[i-1]*fl[i-1] + v[i]*fr[i-1]
246
247 fl.resize(ngrd-1);
248 fr.resize(ngrd-1);
249 for(size_t i=0; i<ngrd-1; i++) {
250 fl[i] = dx[i+1]/(dx[i]+dx[i+1]);
251 fr[i] = 1.0-fl[i];
252 }
253
254}
263
264void ignis::writeFileHdf5(const string gname, const string timeType) {
265
266 vector<double> mixf;
267 double mixfLbc;
268 double mixfRbc;
269 if(!isPremixed) {
270 mixf.resize(ngrd);
271 for(int i=0; i<ngrd; i++)
272 mixf[i] = strm->getMixtureFraction(&y[i][0]);
273 mixfLbc = strm->getMixtureFraction(&yLbc[0]);
274 mixfRbc = strm->getMixtureFraction(&yRbc[0]);
275 }
276
277 vector<double> rho(ngrd);
278 vector<double> h(ngrd);
279 for(int i=0; i<ngrd; i++) {
280 gas->setState_TPY(T[i], P, &y[i][0]);
281 rho[i] = gas->density();
282 h[i] = gas->enthalpy_mass();
283 }
284
285 gas->setState_TPY(TLbc, P, &yLbc[0]);
286 double rhoLbc = gas->density();
287 gas->setState_TPY(TRbc, P, &yRbc[0]);
288 double rhoRbc = gas->density();
289
290 //---------- x, mixf, T, h, rho, y, soot
291
292 vector<double> field;
293 vector<vector<double>> field2;
294
295 Group grp = fh5->createGroup(gname);
296
297 grp.createAttribute<string>("caseType", isPremixed ? "premixed" : (isFlamelet ? "flamelet" : "diffusion"));
298 grp.createAttribute<double>(isFlamelet ? "Chi0" : "L", isFlamelet ? chi0 : L);
299 grp.createAttribute<string>("timeType", timeType);
300
301 field = x; field.insert(field.begin(), 0.0); field.push_back(L);
302 DataSet dset = grp.createDataSet("x", field);
303 dset.createAttribute<string>("units", "m");
304 dset.createAttribute<string>("note", "grid cell positions");
305
306 if(!isPremixed) {
307 field = mixf; field.insert(field.begin(), mixfLbc); field.push_back(mixfRbc);
308 dset = grp.createDataSet("mixf", field);
309 dset.createAttribute<string>("units", "--");
310 }
311
312 field = T; field.insert(field.begin(), TLbc); field.push_back(isPremixed ? T.back() : TRbc);
313 dset = grp.createDataSet("T", field);
314 dset.createAttribute<string>("units", "K");
315
316 field = h; field.insert(field.begin(), hLbc); field.push_back(isPremixed ? h.back() : hRbc);
317 dset = grp.createDataSet("h", field);
318 dset.createAttribute<string>("units", "J/kg");
319
320 field = rho; field.insert(field.begin(), rhoLbc); field.push_back(isPremixed ? rho.back() : rhoRbc);
321 dset = grp.createDataSet("rho", field);
322 dset.createAttribute<string>("units", "kg/m3");
323
324 field2 = y; field2.insert(field2.begin(), yLbc); field2.push_back(isPremixed ? y.back() : yRbc);
325 dset = grp.createDataSet("y", field2);
326 dset.createAttribute<string>("units", "--");
327 dset.createAttribute<string>("note", "mass fractions");
328 dset.createAttribute<vector<string>>("species names", gas->speciesNames());
329
330 if(doSoot) {
331 vector<double> sootBC(nsoot, 0.0);
332 field2 = sootvars; field2.insert(field2.begin(), sootBC); field2.push_back(isPremixed ? sootvars.back() : sootBC);
333 dset = grp.createDataSet("soot", field2);
334 dset.createAttribute<string>("units", "mass moment: kg^k/m3");
335 }
336}
337
344
345void ignis::writeFile(const string fname) {
346
347 //-------------- compute auxiliary quantities
348
349 vector<double> mixf;
350 double mixfLbc;
351 double mixfRbc;
352 if(!isPremixed) {
353 mixf.resize(ngrd);
354 for(int i=0; i<ngrd; i++)
355 mixf[i] = strm->getMixtureFraction(&y[i][0]);
356 mixfLbc = strm->getMixtureFraction(&yLbc[0]);
357 mixfRbc = strm->getMixtureFraction(&yRbc[0]);
358 }
359
360 vector<double> rho(ngrd);
361 vector<double> h(ngrd);
362 for(int i=0; i<ngrd; i++) {
363 gas->setState_TPY(T[i], P, &y[i][0]);
364 rho[i] = gas->density();
365 h[i] = gas->enthalpy_mass();
366 }
367 gas->setState_TPY(TLbc, P, &yLbc[0]);
368 double rhoLbc = gas->density();
369 gas->setState_TPY(TRbc, P, &yRbc[0]);
370 double rhoRbc = gas->density();
371
372 //--------------
373
374 ofstream ofile(fname.c_str());
375 if(!ofile) {
376 cerr << "\n\n***************** ERROR OPENING FILE " << fname << endl << endl;
377 exit(0);
378 }
379
380 ofile << "# L (m) = " << L << endl;
381 ofile << "# P (Pa) = " << P << endl;
382
383 ofile << "#";
384 int j=1;
385 ofile << setw(15) << "00" << j++ << "_x";
386 if(!isPremixed) ofile << setw(13) << "00" << j++ << "_mixf";
387 ofile << setw(16) << "00" << j++ << "_T";
388 ofile << setw(16) << "00" << j++ << "_h";
389 ofile << setw(10) << "00" << j++ << "_density";
390 for(int k=0; k<nsp; k++) {
391 stringstream ss; ss << setfill('0') << setw(3) << j++ << "_" << gas->speciesName(k);
392 ofile << setw(19) << ss.str();
393 }
394 if(doSoot) { // todo: make sure this is calling the correct variables // jansenpb
395 for(int u=0; u<nsoot; u++) {
396 stringstream yy; yy << setfill('0') << setw(3) << j++ << "_" << "M"<<u;
397 ofile << setw(19) << yy.str();
398 }}
399
400 ofile << scientific;
401 ofile << setprecision(10);
402
403 ofile << endl;
404 ofile << setw(19) << 0;
405 if(!isPremixed) ofile << setw(19) << mixfLbc;
406 ofile << setw(19) << TLbc;
407 ofile << setw(19) << hLbc;
408 ofile << setw(19) << rhoLbc;
409 for(int k=0; k<nsp; k++)
410 ofile << setw(19) << yLbc[k];
411 if(doSoot) {
412 for(int k=0; k<nsoot; k++)
413 ofile << setw(19) << 0;
414 }
415
416 for(int i=0; i<ngrd; i++) {
417 ofile << endl;
418 ofile << setw(19) << x[i];
419 if(!isPremixed) ofile << setw(19) << mixf[i];
420 ofile << setw(19) << T[i];
421 ofile << setw(19) << h[i];
422 ofile << setw(19) << rho[i];
423 for(int k=0; k<nsp; k++)
424 ofile << setw(19) << y[i][k];
425 if(doSoot) {
426 for(int k=0; k<nsoot; k++)
427 ofile << setw(19) << sootvars[i][k];
428 }
429 }
430
431 if(isPremixed) {
432 ofile << endl;
433 ofile << setw(19) << L;
434 ofile << setw(19) << T.back();
435 ofile << setw(19) << h.back();
436 ofile << setw(19) << rho.back();
437 for(int k=0; k<nsp; k++)
438 ofile << setw(19) << y.back()[k];
439 if(doSoot) {
440 for(int k=0; k<nsoot; k++)
441 ofile << setw(19) << sootvars.back()[k];
442 }
443 }
444 else {
445 ofile << endl;
446 ofile << setw(19) << L;
447 ofile << setw(19) << mixfRbc;
448 ofile << setw(19) << TRbc;
449 ofile << setw(19) << hRbc;
450 ofile << setw(19) << rhoRbc;
451 for(int k=0; k<nsp; k++)
452 ofile << setw(19) << yRbc[k];
453 if(doSoot) {
454 for(int k=0; k<nsoot; k++)
455 ofile << setw(19) << 0.0;
456 }
457 }
458
459 ofile.close();
460}
461
468
469void ignis::storeState() {
470
471 Pstore = P;
472 ystore = y;
473 Tstore = T;
474 sootstore = sootvars;
475}
476
484
485void ignis::setIC(const std::string icType, string fname) {
486
487 if (icType == "linear") {
488 gas->setState_TPY(TLbc, P, &yLbc[0]);
489 double hLbc = gas->enthalpy_mass();
490 gas->setState_TPY(TRbc, P, &yRbc[0]);
491 double hRbc = gas->enthalpy_mass();
492 for(int i=0; i<ngrd; i++) {
493 for(int k=0; k<nsp; k++)
494 y[i][k] = yLbc[k] + x[i]/L*(yRbc[k]-yLbc[k]);
495 double h = hLbc + x[i]/L*(hRbc-hLbc);
496 gas->setMassFractions(&y[i][0]);
497 gas->setState_HP(h,P);
498 T[i] = doEnergyEqn ? gas->temperature() : LI->interp(x[i]);
499 }
500 }
501
502 //-------------------
503
504 else if (icType == "equilibrium") {
505 setIC("linear");
506 for(int i=0; i<ngrd; i++) {
507 gas->setState_TPY(T[i], P, &y[i][0]);
508 gas->equilibrate("HP");
509 gas->getMassFractions(&y[i][0]);
510 T[i] = doEnergyEqn ? gas->temperature() : LI->interp(x[i]); // redundant
511 }
512 }
513
514 //-------------------
515
516 else if (icType == "stored") {
517 P = Pstore;
518 y = ystore;
519 T = Tstore;
520 sootvars = sootstore;
521 }
522
523 //-------------------
524
525 else if (icType == "premixed") {
526 if(doEnergyEqn) {
527 gas->setMassFractions(&yLbc[0]);
528 gas->setState_HP(hLbc,P);
529 gas->equilibrate("HP");
530 for(int i=0; i<ngrd; i++) {
531 gas->getMassFractions(&y[i][0]);
532 T[i] = gas->temperature();
533 }
534 }
535 else {
536 for(int i=0; i<ngrd; i++) {
537 gas->setState_TPY(LI->interp(x[i]), P, &yLbc[0]);
538 gas->equilibrate("TP");
539 gas->getMassFractions(&y[i][0]);
540 T[i] = gas->temperature();
541 }
542 }
543 }
544
545 //-------------------
546
547 // else if (icType == "file") {{{{
548 // ifstream ifile(fname.c_str());
549 // if(!ifile) {
550 // cerr << "\n\n***************** ERROR OPENING FILE " << fname << endl << endl;
551 // exit(0);
552 // }
553 // string s1;
554 // getline(ifile, s1); // get header lines
555 // getline(ifile, s1);
556 // getline(ifile, s1);
557
558 // vector<double> XX;
559 // vector<double> HH;
560 // vector<vector<double> > YY;
561 // double d;
562 // while(!ifile.eof()) {
563 // ifile >> d; XX.push_back(d); // x
564 // ifile >> d; // T
565 // ifile >> d; HH.push_back(d); // h
566 // YY.push_back(vector<double>());
567 // for(int k=0; k<nsp; k++) {
568 // ifile >> d; YY.back().push_back(d);
569 // }
570 // }
571
572 // }}}}
573
574 //-------------------
575
576 // Soot is inialized to zero in the constructor
577
578 storeState();
579
580}
581
589
590void ignis::setFluxesUnity() {
591
592 //---------- cell center density and diffusivity
593
594 vector<double> D(ngrd);
595 vector<double> density(ngrd);
596 vector<double> nu; if(doSoot) nu.resize(ngrd);
597 vector<double> h(ngrd);
598 for(int i=0; i<ngrd; i++) {
599 gas->setMassFractions_NoNorm(&y[i][0]);
600 gas->setState_TP(T[i], P);
601 density[i] = gas->density();
602 D[i] = trn->thermalConductivity()/(density[i]*gas->cp_mass());
603 h[i] = gas->enthalpy_mass();
604 if(doSoot) nu[i] = trn->viscosity()/density[i];
605 }
606
607 //---------- interpolate to face center density and diffusivity
608
609 vector<double> D_f(ngrd+1);
610 vector<double> density_f(ngrd+1);
611 vector<double> nu_f; if(doSoot) nu_f.resize(ngrd+1);
612 vector<double> T_f; if(doSoot) T_f.resize(ngrd+1);
613
614 gas->setState_TPY(TLbc, P, &yLbc[0]); // this is only approximate for composition for premixed burner
615 density_f[0] = gas->density();
616 D_f[0] = trn->thermalConductivity()/(density_f[0]*gas->cp_mass());
617 if(doSoot){
618 nu_f[0] = trn->viscosity()/density_f[0];
619 T_f[0] = TLbc;
620 }
621
622 if(!isPremixed) {
623 gas->setState_TPY(TRbc, P, &yRbc[0]);
624 density_f.back() = gas->density();
625 D_f.back() = trn->thermalConductivity()/(density_f.back()*gas->cp_mass());
626 if(doSoot) {
627 nu_f.back() = trn->viscosity()/density_f.back();
628 T_f.back() = TRbc;
629 }
630 }
631 else{ // extrapolate half a cell
632 density_f.back() = density.back() +(density[ngrd-1] - density[ngrd-2])/(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
633 D_f.back() = D.back() +(D[ngrd-1] - D[ngrd-2])/(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
634 if(doSoot){
635 nu_f.back() = nu.back() +(nu[ngrd-1] - nu[ngrd-2])/(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
636 T_f.back() = T.back() +(T[ngrd-1] - T[ngrd-2])/(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
637 }
638 }
639
640 for(int i=1, im=0; i<ngrd; i++, im++) { // interpolate
641 density_f[i] = density[im]*fl[im] + density[i]*fr[im];
642 D_f[i] = D[im] *fl[im] + D[i] *fr[im];
643 if(doSoot) {
644 nu_f[i] = nu[im]*fl[im] + nu[i]*fr[im];
645 T_f[i] = T[im] *fl[im] + T[i] *fr[im];
646 }
647 }
648
649 //---------- fluxes y
650
651 for(int k=0; k<nsp; k++) {
652 if(isPremixed) {
653 flux_y[0][k] = mflux * yLbc[k];
654 flux_y[ngrd][k] = mflux * y[ngrd-1][k];
655 }
656 else {
657 flux_y[0][k] = -density_f[0] *D_f[0] *(y[0][k]-yLbc[k]) *2/dx[0];
658 flux_y[ngrd][k] = -density_f.back()*D_f.back()*(yRbc[k]-y[ngrd-1][k])*2/dx.back();
659 }
660 for(int i=1; i<ngrd; i++) {
661 flux_y[i][k] = -density_f[i]*D_f[i]*(y[i][k]-y[i-1][k])*2/(dx[i-1]+dx[i]);
662 if(isPremixed) flux_y[i][k] += mflux*y[i-1][k];
663 }
664 }
665
666 //---------- fluxes h
667
668 for(int i=1; i<ngrd; i++) {
669 flux_h[i] = -density_f[i]*D_f[i]*(h[i]-h[i-1])*2/(dx[i-1]+dx[i]);
670 if(isPremixed) flux_h[i] += mflux*h[i-1];
671 }
672 if(!isPremixed) {
673 flux_h.back() = -density_f.back()*D_f.back()*(hRbc-h.back())*2/dx.back();
674 flux_h[0] = -density_f[0]*D_f[0]*(h[0]-hLbc)*2/dx[0];
675 }
676 else {
677 flux_h.back() = mflux*h[ngrd-1];
678
679 flux_h[0] = 0.0;
680 vector<double> hsp(nsp); // J/kg species i
681 gas->getEnthalpy_RT(&hsp[0]); // (hhat/RT) where hhat = J/kmol
682 for(size_t k=0; k<nsp; k++){ // --> hsp = J/kg species i
683 hsp[k] *= TLbc*Cantera::GasConstant/gas->molecularWeight(k);
684 flux_h[0] += hsp[k]*flux_y[0][k];
685 }
686 gas->setState_TPY(TLbc, P, &yLbc[0]);
687 flux_h[0] -= trn->thermalConductivity() * (T[1]-TLbc)/(dx[0]*0.5);
688 }
689
690 //---------- fluxes soot
691
692 if(doSoot) { // thermophoretic
693 for(int k=0; k<nsoot; k++) {
694 if(isPremixed) {
695 flux_soot[0][k] = 0.0;
696 flux_soot[ngrd][k] = mflux/density_f.back()*sootvars[ngrd-1][k] -
697 0.556*sootvars[ngrd-1][k]*nu_f.back()/T_f.back()*
698 (T_f.back()-T.back())/dx.back()*2;
699 }
700 else {
701 flux_soot[0][k] = -0.556*sootvars[0][k]*nu_f[0]/TLbc*(T[0]-TLbc)/dx[0]*2;
702 flux_soot[ngrd][k] = -0.556*sootvars[ngrd-1][k]*nu_f.back()/TRbc*(TRbc-T.back())/dx.back()*2;
703 }
704 for(int i=1; i<ngrd; i++) {
705 flux_soot[i][k] = 0.556*nu_f[i]*(T[i]-T[i-1])/T_f[i]*2/(dx[i-1]+dx[i]);
706 flux_soot[i][k] *= (flux_soot[i][k] > 0 ? -sootvars[i][k] : -sootvars[i-1][k]); // upwind
707 if(isPremixed) flux_soot[i][k] += mflux/density_f[i] * sootvars[i-1][k]; // upwind
708 }
709 }
710 }
711
712 /*if(doSoot) { // unity Le like species
713 for(int k=0; k<nsoot; k++) {
714 if(isPremixed) {
715 flux_soot[0][k] = 0.0;
716 flux_soot[ngrd][k] = mflux/density_f.back() * sootvars[ngrd-1][k];
717 }
718 else {
719 flux_soot[0][k] = -density_f[0] *D_f[0] *(sootvars[0][k]/density[0]-0.0) *2/dx[0];
720 flux_soot[ngrd][k] = -density_f.back()*D_f.back()*(0.0-sootvars[ngrd-1][k]/density[ngrd-1])*2/dx.back();
721 }
722 for(int i=1; i<ngrd; i++) {
723 flux_soot[i][k] = -density_f[i]*D_f[i]*(sootvars[i][k]/density[i]-sootvars[i-1][k]/density[i-1])*2/(dx[i-1]+dx[i]);
724 if(isPremixed) flux_soot[i][k] += mflux/density_f[i] * sootvars[i-1][k];
725 }
726 }
727 }*/
728}
729
737
738void ignis::setFluxes() {
739
740 //---------- cell center density and diffusivity
741
742 vector<vector<double>> D(ngrd, vector<double> (nsp,0.0));
743 vector<double> density(ngrd);
744 vector<double> M(ngrd);
745 vector<double> tcond(ngrd);
746 vector<double> nu; if(doSoot) nu.resize(ngrd);
747 for(int i=0; i<ngrd; i++) {
748 gas->setMassFractions_NoNorm(&y[i][0]);
749 gas->setState_TP(T[i], P);
750 density[i] = gas->density();
751 M[i] = gas->meanMolecularWeight();
752 trn->getMixDiffCoeffs(&D[i][0]);
753 tcond[i] = trn->thermalConductivity();
754 if(doSoot)
755 nu[i] = trn->viscosity()/density[i];
756 }
757
758 //---------- interpolate to face center density and diffusivity
759
760 vector<vector<double>> D_f(ngrd+1, vector<double> (nsp,0.0));
761 vector<vector<double>> y_f(ngrd+1, vector<double> (nsp,0.0));
762 vector<double> density_f(ngrd+1);
763 vector<double> M_f(ngrd+1);
764 vector<double> T_f(ngrd+1);
765 vector<double> tcond_f(ngrd+1);
766 vector<double> nu_f; if(doSoot) nu_f.resize(ngrd+1);
767
768 gas->setState_TPY(TLbc, P, &yLbc[0]); // this is only approximate for composition for premixed burner
769 density_f[0] = gas->density();
770 M_f[0] = gas->meanMolecularWeight();
771 trn->getMixDiffCoeffs(&D_f[0][0]);
772 tcond_f[0] = trn->thermalConductivity();
773 if(doSoot) nu_f[0] = trn->viscosity()/density_f[0];
774 T_f[0] = TLbc;
775 y_f[0] = yLbc;
776
777 if(!isPremixed) {
778 gas->setState_TPY(TRbc, P, &yRbc[0]);
779 density_f.back() = gas->density();
780 M_f.back() = gas->meanMolecularWeight();
781 trn->getMixDiffCoeffs(&(D_f.back()[0]));
782 tcond_f.back() = trn->thermalConductivity();
783 T_f.back() = TRbc;
784 y_f.back() = yRbc;
785 if(doSoot)
786 nu_f.back() = trn->viscosity()/density_f.back();
787 }
788 else {
789
790 density_f.back() = density.back() +(density[ngrd-1] - density[ngrd-2])/(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
791 tcond_f.back() = tcond.back() +(tcond[ngrd-1] - tcond[ngrd-2]) /(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
792 T_f.back() = T.back() +(T[ngrd-1] - T[ngrd-2]) /(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
793 M_f.back() = M.back() +(M[ngrd-1] - M[ngrd-2]) /(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
794 for(int k=0; k<nsp; k++) {
795 D_f.back()[k] = D.back()[k] +(D[ngrd-1][k] - D[ngrd-2][k]) /(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
796 y_f.back()[k] = y.back()[k] +(y[ngrd-1][k] - y[ngrd-2][k]) /(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
797 }
798 if(doSoot)
799 nu_f.back() = nu.back() +(nu[ngrd-1] - nu[ngrd-2]) /(dx[ngrd-1]+dx[ngrd-2])*2.0*(L-x[ngrd-1]);
800 }
801
802 for (int i=1, im=0; i<ngrd; i++, im++) {
803 density_f[i] = density[im]*fl[im] + density[i]*fr[im];
804 T_f[i] = T[im] *fl[im] + T[i] *fr[im];
805 tcond_f[i] = tcond[im] *fl[im] + tcond[i] *fr[im];
806 M_f[i] = M[im] *fl[im] + M[i] *fr[im];
807 if(doSoot)
808 nu_f[i] = nu[im] *fl[im] + nu[i] *fr[im];
809 for(int k=0; k<nsp; k++) {
810 D_f[i][k] = D[im][k] *fl[im] + D[i][k] *fr[im];
811 y_f[i][k] = y[im][k] *fl[im] + y[i][k] *fr[im];
812 }
813 }
814
815 //---------- fluxes y
816
817 double jstar; // correction flux so that all fluxes sum to zero. This is equal to using a correction velocity
818 // j_i_corrected = j_i - Yi*jstar; jstar = sum(j_i).
819
820 //--------- Boundary faces
821
822 // Left boundary
823
824 if(isPremixed)
825 for(int k=0; k<nsp; k++)
826 flux_y[0][k] = mflux * yLbc[k];
827 else {
828 jstar = 0.0;
829 for(int k=0; k<nsp; k++) {
830 flux_y[0][k] = -density_f[0]*D_f[0][k]*(y[0][k]-yLbc[k])*2/dx[0]
831 -density_f[0]*D_f[0][k]*y_f[0][k]*(M[0]-M_f[0])*2/dx[0]/M_f[0];
832 jstar += flux_y[0][k];
833 }
834 for(int k=0; k<nsp; k++) {
835 flux_y[0][k] -= y_f[0][k]*jstar;
836 }
837 }
838
839 // Right boundary
840
841 if(isPremixed)
842 for(int k=0; k<nsp; k++)
843 flux_y[ngrd][k] = mflux * y[ngrd-1][k];
844 else {
845 jstar = 0.0;
846 for(int k=0; k<nsp; k++) {
847 flux_y[ngrd][k] = -density_f.back()*D_f.back()[k]*(yRbc[k]-y[ngrd-1][k])*2/dx.back()
848 -density_f.back()*D_f.back()[k]*y_f.back()[k]*(M_f.back()-M[ngrd-1])*2/dx.back()/M_f.back();
849 jstar += flux_y[ngrd][k];
850 }
851 for(int k=0; k<nsp; k++) {
852 flux_y[ngrd][k] -= y_f[ngrd][k]*jstar;
853 }
854 }
855
856 // Interior
857
858 for (int i=1; i<ngrd; i++) {
859 jstar = 0.0;
860 for(int k=0; k<nsp; k++) {
861 flux_y[i][k] = -density_f[i]*D_f[i][k]*(y[i][k]-y[i-1][k])*2/(dx[i-1]+dx[i])
862 -density_f[i]*D_f[i][k]*y_f[i][k]*(M[i]-M[i-1])*2/(dx[i-1]+dx[i])/M_f[i];
863 jstar += flux_y[i][k];
864 }
865 for(int k=0; k<nsp; k++) {
866 flux_y[i][k] -= y_f[i][k]*jstar;
867 if(isPremixed) flux_y[i][k] += mflux*y[i-1][k];
868 }
869 }
870
871 //---------- fluxes h
872
873 // thermal conductivity portion
874
875 flux_h[0] = -tcond_f[0]*(T[0]-T_f[0])*2/dx[0];
876 flux_h.back() = -tcond_f.back()*(T_f.back()-T.back())*2/dx.back();
877 for(int i=1; i<ngrd; i++)
878 flux_h[i] = -tcond_f[i]*(T[i]-T[i-1])*2/(dx[i-1]+dx[i]);
879
880 // species portion (builds in advective flux if isPremixed)
881
882 vector<double> hsp(nsp);
883
884 gas->setState_TPY(TLbc, P, &yLbc[0]);
885 gas->getEnthalpy_RT(&hsp[0]);
886 for(size_t k=0; k<nsp; k++)
887 flux_h[0] += flux_y[0][k]*hsp[k]*TLbc*Cantera::GasConstant/gas->molecularWeight(k);
888
889 gas->setState_TPY(T_f.back(), P, &(y_f.back()[0]));
890 gas->getEnthalpy_RT(&hsp[0]);
891 for(size_t k=0; k<nsp; k++)
892 flux_h.back() += flux_y.back()[k]*hsp[k]*T_f.back()*Cantera::GasConstant/gas->molecularWeight(k);
893
894 for(size_t i=1; i<ngrd; i++) {
895 gas->setState_TP(T_f[i], P); // hsp depends on T but not on y
896 gas->getEnthalpy_RT(&hsp[0]);
897 for(size_t k=0; k<nsp; k++)
898 flux_h[i] += flux_y[i][k]*hsp[k]*T_f[i]*Cantera::GasConstant/gas->molecularWeight(k);
899 }
900
901 //---------- fluxes soot
902
903 if(doSoot) { // thermophoretic
904 for(int k=0; k<nsoot; k++) {
905 if(isPremixed) {
906 flux_soot[0][k] = 0.0;
907 flux_soot[ngrd][k] = mflux/density_f.back()*sootvars[ngrd-1][k] -
908 0.556*sootvars[ngrd-1][k]*nu_f.back()/T_f.back()*
909 (T_f.back()-T.back())/dx.back()*2;
910 for(int i=1; i<ngrd; i++) {
911 flux_soot[i][k] = -0.556*nu_f[i]*(T[i]-T[i-1])/T_f[i]*2/(dx[i-1]+dx[i]) + mflux/density_f[i]; // soot velocity
912 flux_soot[i][k] *= (flux_soot[i][k] > 0 ? sootvars[i-1][k] : sootvars[i][k]); // upwind flux
913 }
914 }
915 else {
916 flux_soot[0][k] = -0.556*sootvars[0][k]*nu_f[0]/TLbc*(T[0]-TLbc)/dx[0]*2;
917 flux_soot[ngrd][k] = -0.556*sootvars[ngrd-1][k]*nu_f.back()/TRbc*(TRbc-T.back())/dx.back()*2;
918 for(int i=1; i<ngrd; i++) {
919 flux_soot[i][k] = -0.556*nu_f[i]*(T[i]-T[i-1])/T_f[i]*2/(dx[i-1]+dx[i]); // soot velocity
920 flux_soot[i][k] *= (flux_soot[i][k] > 0 ? sootvars[i-1][k] : sootvars[i][k]); // upwind
921 }
922
923 //flux_soot[0][k] += -nu_f[0]*(sootvars[0][k]-0.0)/dx[0]*2;
924 //flux_soot[ngrd][k] += -nu_f.back()*(0.0-sootvars[ngrd-1][k])/dx.back()*2;
925 //for(int i=1; i<ngrd; i++)
926 // flux_soot[i][k] += -nu_f[i]*(sootvars[i][k]-sootvars[i-1][k])/(dx[i-1]+dx[i])*2;
927 }
928 }
929 }
930
931}
932
941
942void ignis::solveSS() {
943
944 //---------- transfer variables into single array
945
946 vector<double> vars(nvarA);
947
948 for(size_t i=0; i<ngrd; i++) {
949 for(size_t k=0; k<nsp; k++)
950 vars[Ia(i,k)] = y[i][k];
951 vars[Ia(i,nvar-1)] = T[i]/Tscale; // dolh comment to remove h
952 }
953
954 //-------------- set vars0, F0 for homotopy
955
956 vars0 = vars;
957 F0.resize(nvarA);
958 Func(&vars0[0], &F0[0]);
959
960 //---------- setup solver
961
962 vector<double> scales_v(nvarA, 1.0);
963 vector<double> scales_f(nvarA, 1.0);
964 vector<double> constraints(nvarA, 0.0);
965
966 for(int i=0; i<nvarA; i++) {
967 scales_v[i] = vars[i] < 1E-3 ? 1E-3 : vars[i];
968 scales_f[i] = vars[i] < 1E-3 ? 1E-3 : vars[i];
969 }
970
971
972 double ftol = 5E-6; //5E-4; // default: 5E-6
973 double stol = 2E-11; // default: 2E-11
974 int mu = nvar*2-1;
975 int ml = nvar*2-1;
976 solver_kinsol s_kin(this, nvarA, scales_v, scales_f, constraints, mu, ml, ftol, stol);
977
978 //---------- solve
979
980 // s = 0.0;
981 // s_kin.solve(Func_kinsol, vars);
982 // s = 0.5;
983 // s_kin.solve(Func_kinsol, vars);
984
985 for(s=0.0; s<=1.0; s+=0.01) {
986 cout << endl << "s = " << s; cout.flush();
987 s_kin.solve(Func_kinsol, vars);
988 }
989
990 //---------- transfer variables back
991
992 for(size_t i=0; i<ngrd; i++) {
993 for(size_t k=0; k<nsp; k++)
994 y[i][k] = vars[Ia(i,k)];
995 T[i] = vars[Ia(i,nvar-1)]*Tscale; // dolh comment to remove h
996 }
997}
998
1007
1008int ignis::Func(const double *vars, double *F) {
1009
1010 //------------ transfer variables
1011
1012 for(size_t i=0; i<ngrd; i++) {
1013 for(size_t k=0; k<nsp; k++)
1014 y[i][k] = vars[Ia(i,k)];
1015 T[i] = vars[Ia(i,nvar-1)]*Tscale; // dolh comment to remove h
1016 }
1017
1018 //------------ set function values
1019
1020 if(doLe1)
1021 setFluxesUnity();
1022 else
1023 setFluxes();
1024
1025 vector<double> Q(ngrd);
1026 if(doRadiation) setQrad(Q);
1027
1028 vector<double> rr(nsp);
1029
1030 for(size_t i=0; i<ngrd; i++) {
1031 gas->setMassFractions_NoNorm(&y[i][0]);
1032 gas->setState_TP(T[i], P);
1033 kin->getNetProductionRates(&rr[0]);
1034 for(size_t k=0; k<nsp; k++)
1035 F[Ia(i,k)] = flux_y[i+1][k] - flux_y[i][k]
1036 - rr[k]*gas->molecularWeight(k)*dx[i];
1037 if(doEnergyEqn) {
1038 F[Ia(i,nvar-1)] = (flux_h[i+1] - flux_h[i])/hscale; // dolh comment to remove h
1039 if(doRadiation) F[Ia(i,nvar-1)] -= Q[i]/hscale;
1040 }
1041 else
1042 F[Ia(i,nvar-1)] = 0.0; // dolh comment to remove h
1043 }
1044
1045 //------------ augment for homotopy
1046
1047 for(size_t i=0; i<nvarA; i++)
1048 F[i] = F[i] + (s-1.0)*F0[i];
1049
1050 return 0;
1051}
1052
1062
1063int Func_kinsol(N_Vector varsKS, N_Vector fvec, void *user_data) {
1064
1065 ignis *flm = static_cast<ignis *>(user_data);
1066
1067 double *vars = N_VGetArrayPointer(varsKS);
1068 double *F = N_VGetArrayPointer(fvec);
1069
1070 int rv = flm->Func(vars, F);
1071
1072 return rv;
1073}
1074
1081
1082void ignis::setQrad(vector<double> &Q) {
1083
1084 vector<double> kabs, awts;
1085 double fvsoot = 0.0;
1086 double xH2O, xCO2, xCO, xCH4;
1087 int isp;
1088
1089 for(int i=0; i<ngrd; i++) {
1090 isp = gas->speciesIndex("H2O");
1091 xH2O = y[i][isp]/gas->molecularWeight(isp)*gas->meanMolecularWeight();
1092 isp = gas->speciesIndex("CO2");
1093 xCO2 = y[i][isp]/gas->molecularWeight(isp)*gas->meanMolecularWeight();
1094 isp = gas->speciesIndex("CO");
1095 xCO = y[i][isp]/gas->molecularWeight(isp)*gas->meanMolecularWeight();
1096 isp = gas->speciesIndex("CH4");
1097 xCH4 = y[i][isp]/gas->molecularWeight(isp)*gas->meanMolecularWeight();
1098 radProps->get_k_a(kabs, awts, T[i], P, fvsoot, xH2O, xCO2, xCO, xCH4);
1099
1100 int nGGa = radProps->get_nGGa();
1101 for(int j=0; j<nGGa; j++)
1102 Q[i] += -4.0*rad::sigma*kabs[j]*(awts[j]*pow(T[i],4.0)
1103 -awts_sur[0]*pow(TLbc,4.0));
1104 }
1105}
1106
1121
1122void ignis::solveUnsteady(const double nTauRun, const int nSteps, const bool doWriteTime,
1123 const double Tmin, const double Tmax) {
1124
1125 //---------- transfer variables into single array
1126 vector<double> vars(nvarA);
1127
1128 for(size_t i=0; i<ngrd; i++) {
1129 for(size_t k=0; k<nsp; k++)
1130 vars[Ia(i,k)] = y[i][k];
1131 if(doSoot) {
1132 for(size_t k=nsp; k<nsp+nsoot; k++) {
1133 vars[Ia(i,k)] = sootvars[i][k-nsp]/sootScales[k-nsp]; // dont't forget sootscales
1134 }
1135 }
1136 vars[Ia(i,nvar-1)] = T[i]/Tscale; // dolh comment to remove h
1137 }
1138
1139 //---------- setup solver
1140
1141 double rtol = 1E-4; // 1E-4 --> error at 0.01%, cvode documentation recomends < 1E-3
1142 vector<double> atol(nvarA, 1.0); // noise level of the given variable
1143
1144 for(int i=0; i<nvarA; i++)
1145 atol[i] = 1E-12;
1146 //atol[i] = vars[i] < 1E-3 ? 1E-3 : vars[i];
1147
1148
1149 int mu = nvar*2-1;
1150 int ml = nvar*2-1;
1151 integrator_cvode integ(rhsf_cvode, this, nvarA, rtol, atol, mu, ml, vars);
1152
1153 //---------- solve
1154
1155 double D = 0.00035; // avg thermal diffusivity
1156 double t = 0.0;
1157 if(isPremixed) gas->setState_TPY(TLbc, P, &yLbc[0]); // needed fro density in tend
1158
1159 double tau;
1160 if(isPremixed)
1161 tau = L/(mflux/gas->density());
1162 else if(isFlamelet)
1163 tau = 1.0/chi0;
1164 else
1165 tau = L*L/D;
1166 double tend = nTauRun * tau;
1167
1168 double dt = tend/nSteps;
1169 dT = Tmax==Tmin ? 0.0 : (Tmax-(Tmin+0.1))/nSteps;
1170 Ttarget = Tmax - dT;
1171 isave = 1;
1172
1173 for(int istep=1; istep<=nSteps; istep++, t+=dt) {
1174 integ.integrate(vars, dt);
1175 if(doWriteTime && dT <= 0.0) { // write in time; (write in Temp is in rhsf)
1176 stringstream ss;
1177 if(isFlamelet)
1178 ss << "X_" << chi0 << "U_" << setfill('0') << setw(3) << isave++;
1179 else
1180 ss << "L_" << L << "U_" << setfill('0') << setw(3) << isave++;
1181 string fname = ss.str();
1182 writeFile(fname+".dat");
1183 writeFileHdf5(fname, "unsteady");
1184 }
1185 }
1186
1187 //---------- transfer variables back
1188
1189 for(size_t i=0; i<ngrd; i++) {
1190 for(size_t k=0; k<nsp; k++)
1191 y[i][k] = vars[Ia(i,k)];
1192 if(doSoot) {
1193 for(size_t k=nsp; k<nsp+nsoot; k++) {
1194 sootvars[i][k-nsp] = vars[Ia(i,k)]*sootScales[k-nsp];
1195 if (isFlamelet) {
1196 gas->setState_TP(T[i], P);
1197 sootvars[i][k-nsp] = vars[Ia(i,k)] * gas->density();
1198 }
1199 }
1200 }
1201 T[i] = vars[Ia(i,nvar-1)]*Tscale; // dolh comment to remove h
1202 }
1203
1204 Ttarget = 0.0;
1205}
1206
1214
1215int ignis::rhsf(const double *vars, double *dvarsdt) {
1216
1217 //------------ transfer variables
1218
1219 for(size_t i=0; i<ngrd; i++) {
1220 for(size_t k=0; k<nsp; k++)
1221 y[i][k] = vars[Ia(i,k)];
1222 if(doSoot) {
1223 for(size_t k=nsp; k<nsp+nsoot; k++) {
1224 sootvars[i][k-nsp] = vars[Ia(i,k)]*sootScales[k-nsp];
1225 }
1226 }
1227 T[i] = vars[Ia(i,nvar-1)]*Tscale; // dolh comment to remove h
1228 }
1229
1230 //------------ set rates dvarsdt (dvars/dt)
1231
1232 if(doLe1)
1233 setFluxesUnity();
1234 else
1235 setFluxes();
1236
1237 vector<double> rr(nsp);
1238 vector<double> Q(ngrd);
1239 if(doRadiation) setQrad(Q);
1240
1241 vector<double> yPAH; if(doSoot) yPAH.resize(6,0.0);
1242
1243 vector<double> yGasForSM; if(doSoot) yGasForSM.resize( (size_t)gasSp::size );
1244
1245 for(size_t i=0; i<ngrd; i++) {
1246 gas->setMassFractions_NoNorm(&y[i][0]);
1247 gas->setState_TP(T[i], P);
1248 kin->getNetProductionRates(&rr[0]); // kmol/m3*s
1249 double rho = gas->density();
1250 double mu = trn->viscosity();
1251 if(doSoot) {
1252 for(size_t kSootGases=0; kSootGases<(size_t)gasSp::size; kSootGases++) {
1253 size_t kgas = gas->speciesIndex(gasSpMapIS[kSootGases]);
1254 yGasForSM[kSootGases] = (kgas != Cantera::npos) ? y[i][kgas] : 0.0;
1255 }
1256 SMstate->setState(T[i], P, rho, mu, yGasForSM, yPAH, sootvars[i], nsoot);
1257 SM->setSourceTerms(*SMstate);
1258 }
1259 for(size_t k=0; k<nsp; k++)
1260 dvarsdt[Ia(i,k)] = -(flux_y[i+1][k] - flux_y[i][k])/(rho*dx[i]) +
1261 rr[k]*gas->molecularWeight(k)/rho;
1262 if(doSoot) {
1263 for(size_t k=nsp; k<nsp+nsoot; k++) {
1264 dvarsdt[Ia(i,k)] = -(flux_soot[i+1][k-nsp] - flux_soot[i][k-nsp])/(dx[i]) + SM->sources.sootSources[k-nsp];
1265 dvarsdt[Ia(i,k)] /= sootScales[k-nsp]; //to match Tscale below
1266 }
1267 // loop over the gas species in the soot model and compare with Cantera
1268 // update the gas source terms from the soot model
1269
1270 for(size_t kSootGases=0; kSootGases<(size_t)gasSp::size; kSootGases++) {
1271 size_t kgas = gas->speciesIndex(gasSpMapIS[kSootGases]);
1272 if(kgas != Cantera::npos)
1273 dvarsdt[Ia(i,kgas)] += SM->sources.gasSources[kSootGases];
1274 }
1275 }
1276
1277 if(doEnergyEqn) {
1278 double cp = gas->cp_mass();
1279 vector<double> hsp(nsp); // J/kg species i
1280 gas->getEnthalpy_RT(&hsp[0]); // (hhat/RT) where hhat = J/kmol
1281 for(size_t k=0; k<nsp; k++) // --> hsp = J/kg species i
1282 hsp[k] *= T[i]*Cantera::GasConstant/gas->molecularWeight(k);
1283 double sum_hkdykdt = 0.0;
1284 for(size_t k=0; k<nsp; k++)
1285 sum_hkdykdt += hsp[k]*dvarsdt[Ia(i,k)];
1286 dvarsdt[Ia(i,nvar-1)] = -(flux_h[i+1] - flux_h[i])/(rho*cp*dx[i]) - sum_hkdykdt/cp;
1287 if(doRadiation) dvarsdt[Ia(i,nvar-1)] += Q[i]/(rho*cp);
1288 dvarsdt[Ia(i,nvar-1)] /= Tscale;
1289 }
1290 else
1291 dvarsdt[Ia(i,nvar-1)] = 0.0;
1292 }
1293 //-------------
1294 double TmaxLocal = *max_element(T.begin(), T.end());
1295 if(TmaxLocal <= Ttarget) {
1296 cout << endl << isave << " " << TmaxLocal << " " << Ttarget << " ";
1297 stringstream ss; ss << "L_" << L << "U_" << setfill('0') << setw(3) << isave++;
1298 string fname = ss.str();
1299 writeFile(fname + ".dat");
1300 writeFileHdf5(fname, "unsteady");
1301 Ttarget -= dT;
1302 }
1303
1304 //-------------
1305
1306 return 0;
1307}
1308
1316
1317int ignis::rhsf_flamelet(const double *vars, double *dvarsdt) {
1318
1319 //------------ transfer variables
1320
1321 for(size_t i=0; i<ngrd; i++) {
1322 for(size_t k=0; k<nsp; k++)
1323 y[i][k] = vars[Ia(i,k)];
1324 if(doSoot)
1325 for(size_t k=nsp; k<nsp+nsoot; k++)
1326 sootvars[i][k-nsp] = vars[Ia(i,k)]*sootScales[k-nsp];
1327 T[i] = vars[Ia(i,nvar-1)]*Tscale; // dolh comment to remove h
1328 }
1329
1330 //------ uncomment to solve for h instead (1 of 2), but only if doRadiation is false
1331 // for(size_t i=0; i<ngrd; i++) {
1332 // gas->setMassFractions(&y[i][0]);
1333 // gas->setState_HP(hLbc*(1.0-x[i]) + hRbc*x[i], P);
1334 // T[i] = gas->temperature();
1335 // }
1336
1337 //------------ set variables
1338
1339 vector<vector<double> > rr(ngrd, vector<double>(nsp));
1340 vector<double> rho(ngrd);
1341 vector<double> mu; if(doSoot) mu.resize(ngrd);
1342 vector<double> D; if(doSoot) D.resize(ngrd); // D mixf = thermal diffusivity
1343
1344 vector<double> yPAH; if(doSoot) yPAH.resize(6,0.0);
1345 vector<double> yGasForSM; if(doSoot) yGasForSM.resize( (size_t)gasSp::size );
1346
1347 vector<double> cp(ngrd);
1348 vector<vector<double> > hsp(ngrd, vector<double>(nsp));
1349 vector<double> hsprrSum(ngrd, 0.0);
1350
1351 for(size_t i=0; i<ngrd; i++) {
1352
1353 //------- set gas
1354
1355 gas->setMassFractions_NoNorm(&y[i][0]);
1356 gas->setState_TP(T[i], P);
1357
1358 //------- reaction rates
1359
1360 kin->getNetProductionRates(&rr[i][0]); // kmol/m3*s
1361 for(size_t k=0; k<nsp; k++)
1362 rr[i][k] *= gas->molecularWeight(k); // kg/m3*2
1363
1364 //------- species enthalpies
1365
1366 gas->getEnthalpy_RT(&hsp[i][0]); // (h/RT) where h = J/kmol
1367 for(size_t k=0; k<nsp; k++) { // --> hsp = J/kg species i
1368 hsp[i][k] *= T[i]*Cantera::GasConstant/gas->molecularWeight(k);
1369 hsprrSum[i] += hsp[i][k] * rr[i][k];
1370 }
1371
1372 //--------
1373
1374 rho[i] = gas->density(); // kg/m3
1375 cp[i] = gas->cp_mass();
1376 if(doSoot) {
1377 mu[i] = trn->viscosity();
1378 D[i] = trn->thermalConductivity()/(rho[i]*cp[i]);
1379 for (size_t k=0; k<nsoot; k++)
1380 sootvars[i][k] *= rho[i];
1381 }
1382 }
1383
1384 //------------ species
1385
1386 vector<double> d2ydz2(ngrd);
1387 vector<double> yy(ngrd); // intermediate transfer array
1388 for(size_t k=0; k<nsp; k++) {
1389 for(size_t i=0; i<ngrd; i++)
1390 yy[i] = y[i][k];
1391 setDerivative2(yLbc[k], yRbc[k], yy, d2ydz2);
1392 for(size_t i=0; i<ngrd; i++)
1393 dvarsdt[Ia(i,k)] = 0.5*chi[i]*d2ydz2[i] + rr[i][k]/rho[i];
1394 }
1395
1396 //------------ energy (temperature)
1397
1398 vector<double> d2Tdz2(ngrd);
1399 setDerivative2(TLbc, TRbc, T, d2Tdz2);
1400
1401 vector<double> dTdz(ngrd);
1402 setDerivative(TLbc, TRbc, T, dTdz);
1403
1404 vector<double> dcpdz(ngrd);
1405 setDerivative(cpLbc, cpRbc, cp, dcpdz);
1406
1407 vector<double> dydzdhdzSum(ngrd, 0.0); //todo: fill this in
1408 vector<double> dykdz(ngrd);
1409 vector<double> dhkdz(ngrd);
1410 vector<double> hh(ngrd); // intermediate transfer array
1411 for(size_t k=0; k<nsp; k++) {
1412 for(size_t i=0; i<ngrd; i++) {
1413 yy[i] = y[i][k];
1414 hh[i] = hsp[i][k];
1415 }
1416 setDerivative(yLbc[k], yRbc[k], yy, dykdz);
1417 setDerivative(hspLbc[k], hspRbc[k], hh, dhkdz);
1418 for(size_t i=0; i<ngrd; i++)
1419 dydzdhdzSum[i] += dykdz[i]*dhkdz[i];
1420 }
1421
1422 vector<double> Q(ngrd);
1423 if(doRadiation) setQrad(Q);
1424
1425 for(size_t i=0; i<ngrd; i++) {
1426 dvarsdt[Ia(i,nvar-1)] = -hsprrSum[i]/(cp[i]*rho[i]) + 0.5*chi[i]*
1427 (d2Tdz2[i] + (dTdz[i]*dcpdz[i] + dydzdhdzSum[i])/cp[i]);
1428 if(doRadiation) dvarsdt[Ia(i,nvar-1)] += Q[i]/(rho[i]*cp[i]);
1429 dvarsdt[Ia(i,nvar-1)] /= Tscale;
1430 }
1431
1432 //------ uncomment to solve for h instead (2 of 2), but only if doRadiation is false
1433 // for(size_t i=0; i<ngrd; i++)
1434 // dvarsdt[Ia(i,nvar-1)] = 0.0;
1435
1436 //---------- soot
1437
1438 if(doSoot) {
1439
1440 double LeS = 1000.0; // soot Lewis number
1441
1442 vector<double> d2sdz2(ngrd);
1443 vector<double> ss(ngrd); // intermediate transfer array
1444 for(size_t k=0; k<nsoot; k++) {
1445 for(size_t i=0; i<ngrd; i++)
1446 ss[i] = sootvars[i][k];
1447 setDerivative2(0.0, 0.0, ss, d2sdz2);
1448 for(size_t i=0; i<ngrd; i++)
1449 dvarsdt[Ia(i,k+nsp)] = 0.5*chi[i]*d2sdz2[i]/LeS; // account for diffusion
1450 }
1451
1452
1453 vector<double> C(ngrd); // velocity (convective) term
1454 vector<double> B(ngrd); // dzdy or sqrt(chi/(2*D))
1455 vector<double> rhoBD(ngrd); // rho[i]*B[i]*D
1456 vector<double> muB_T(ngrd); // 0.556*mu[i]*B[i]/T[i]
1457 double S; // source term
1458
1459 for (int i=0; i<ngrd; i++) {
1460 B[i] = sqrt(chi[i]/(2*D[i])); // fill beta
1461 rhoBD[i] = rho[i]*B[i]*D[i]; // for the setDerivative function drhoDBdz
1462 muB_T[i] = mu[i]*B[i]/T[i]; // for the setDerivative function dmuB_Tdz
1463 }
1464
1465 vector<double> drhoDBdz(ngrd);
1466 setDerivative(0.0, 0.0, rhoBD, drhoDBdz);
1467
1468 vector<double> dmuB_Tdz(ngrd);
1469 setDerivative(0.0, 0.0, muB_T, dmuB_Tdz);
1470
1471 for (int i=0; i<ngrd; i++) {
1472 C[i] = 0.556*mu[i]/(rho[i]*T[i])* B[i]*B[i]*dTdz[i]-
1473 B[i]/rho[i]*drhoDBdz[i];
1474 for(size_t kSootGases=0; kSootGases<(size_t)gasSp::size; kSootGases++) {
1475 size_t kgas = gas->speciesIndex(gasSpMapIS[kSootGases]);
1476 yGasForSM[kSootGases] = (kgas != Cantera::npos) ? y[i][kgas] : 0.0;
1477 }
1478 SMstate->setState(T[i], P, rho[i], mu[i], yGasForSM, yPAH, sootvars[i], nsoot);
1479 SM->setSourceTerms(*SMstate);
1480 for (int k=0; k<nsoot; k++) {
1481 S = 0.556*(sootvars[i][k]/rho[i])/rho[i]* (B[i]*dTdz[i]*dmuB_Tdz[i]+ B[i]*B[i]*mu[i]/T[i]*d2Tdz2[i])+
1482 SM->sources.sootSources[k]/rho[i];
1483 if (C[i] < 0) {
1484 if (i==0)
1485 dvarsdt[Ia(i,nsp+k)] += C[i]*(sootvars[i][k]/rho[i]-0)/(dx[i]/2);
1486 else
1487 dvarsdt[Ia(i,nsp+k)] += C[i]*(sootvars[i][k]/rho[i]-sootvars[i-1][k]/rho[i-1])/dx[i];
1488 }
1489 else {
1490 if (i==ngrd-1)
1491 dvarsdt[Ia(i,nsp+k)] += C[i]*(0-sootvars[i][k]/rho[i])/(dx[i]/2);
1492 else
1493 dvarsdt[Ia(i,nsp+k)] += C[i]*(sootvars[i+1][k]/rho[i+1]-sootvars[i][k]/rho[i])/dx[i];
1494 }
1495 dvarsdt[Ia(i,nsp+k)] += S;
1496 dvarsdt[Ia(i,nsp+k)] /= sootScales[k];
1497 }
1498 for(size_t kSootGases=0; kSootGases<(size_t)gasSp::size; kSootGases++) {
1499 size_t kgas = gas->speciesIndex(gasSpMapIS[kSootGases]);
1500 if(kgas != Cantera::npos)
1501 dvarsdt[Ia(i,kgas)] += SM->sources.gasSources[kSootGases];
1502 }
1503 }
1504 }
1505
1506
1507
1508 //-------------
1509
1510 double TmaxLocal = *max_element(T.begin(), T.end());
1511 if(TmaxLocal <= Ttarget) {
1512 cout << endl << isave << " " << TmaxLocal << " " << Ttarget << " ";
1513 stringstream ss; ss << "X_" << chi0 << "U_" << setfill('0') << setw(3) << isave++;
1514 string fname = ss.str();
1515 writeFile(fname + ".dat");
1516 writeFileHdf5(fname, "unsteady");
1517 Ttarget -= dT;
1518 }
1519
1520 //-------------
1521
1522 return 0;
1523}
1524
1537
1538void ignis::setDerivative(const double vL, const double vR,
1539 const vector<double> &v, vector<double> &dvdx) {
1540
1541 double vfL = vL; // first cell (left side)
1542 double vfR = v[0]*fl[0]+v[1]*fr[0];
1543 dvdx[0] = (vfR-vfL)/dx[0];
1544
1545 for(size_t i=1; i<ngrd-1; i++) { // interior cells
1546 vfL = vfR;
1547 vfR = v[i]*fl[i]+v[i+1]*fr[i];
1548 dvdx[i] = (vfR-vfL)/dx[i];
1549 }
1550 vfL = vfR; // last cell (right side)
1551 vfR = vR;
1552 dvdx[ngrd-1] = (vfR-vfL)/dx[ngrd-1];
1553}
1554
1567
1568void ignis::setDerivative2(const double vL, const double vR,
1569 const vector<double> &v, vector<double> &d2vdx2) {
1570
1571 double dvdxL = (v[0]-vL) /dx[0]*2; // first cell (left side)
1572 double dvdxR = (v[1]-v[0])/(dx[0]+dx[1])*2;
1573 d2vdx2[0] = (dvdxR - dvdxL)/dx[0];
1574
1575 for(size_t i=1; i<ngrd-1; i++) { // interior cells
1576 dvdxL = dvdxR;
1577 dvdxR = (v[i+1]-v[i])/(dx[i]+dx[i+1])*2;
1578 d2vdx2[i] = (dvdxR - dvdxL)/dx[i];
1579 }
1580 dvdxL = dvdxR; // last cell (right side)
1581 dvdxR = (vR-v[ngrd-1])/dx[ngrd-1]*2;
1582 d2vdx2[ngrd-1] = (dvdxR - dvdxL)/dx[ngrd-1];
1583}
1584
1591
1592// nimig18: https://stackoverflow.com/questions/27229371/inverse-error-function-in-c
1593
1594float myErfInv2(float x){
1595 float tt1, tt2, lnx, sgn;
1596 sgn = (x < 0) ? -1.0f : 1.0f;
1597
1598 x = (1 - x)*(1 + x); // x = 1 - x*x;
1599 lnx = logf(x);
1600
1601 tt1 = 2/(M_PI*0.147) + 0.5f * lnx;
1602 tt2 = 1/(0.147) * lnx;
1603
1604 return(sgn*sqrtf(-tt1 + sqrtf(tt1*tt1 - tt2)));
1605}
1606
1607void ignis::setChi(const double _chi0) {
1608
1609 chi0 = _chi0;
1610 chi.resize(ngrd);
1611
1612 // tanh
1613 // double d, e;
1614 // for(size_t i=0; i<ngrd; i++) {
1615 // d = 2*x[i]-1;
1616 // e = 1.0-d*d;
1617 // chi[i] = chi0*e*e;
1618 // }
1619
1620 // erf
1621 double d;
1622 for(size_t i=0; i<ngrd; i++) {
1623 d = myErfInv2(2*x[i]-1);
1624 chi[i] = chi0*exp(-2.0*d*d);
1625 }
1626}
1627
1637
1638int rhsf_cvode(realtype t, N_Vector varsCV, N_Vector dvarsdtCV, void *user_data) {
1639 ignis *flm = static_cast<ignis *>(user_data);
1640
1641 double *vars = N_VGetArrayPointer(varsCV);
1642 double *dvarsdt = N_VGetArrayPointer(dvarsdtCV);
1643
1644 int rv;
1645 if(flm->isFlamelet)
1646 rv = flm->rhsf_flamelet(vars, dvarsdt);
1647 else
1648 rv = flm->rhsf(vars, dvarsdt);
1649
1650 return rv;
1651}
ignis
Definition ignis.h:29
ignis::dx
std::vector< double > dx
grid spacing (m), nonuniform is fine
Definition ignis.h:58
ignis::setGrid
void setGrid(double _L)
Definition ignis.cc:202
ignis::T
std::vector< double > T
temperature (K)
Definition ignis.h:41
ignis::nsp
size_t nsp
number of gas species
Definition ignis.h:34
ignis::F0
std::vector< double > F0
for homotopy approaches
Definition ignis.h:67
ignis::doRadiation
bool doRadiation
radiation flag
Definition ignis.h:77
ignis::flux_soot
std::vector< std::vector< double > > flux_soot
species fluxes: [I(igrid, ksoot)]
Definition ignis.h:88
ignis::nvarA
size_t nvarA
nvar * ngrd
Definition ignis.h:37
ignis::hRbc
double hRbc
h boundary values: left and right (as needed)
Definition ignis.h:51
ignis::ngrd
size_t ngrd
number of interior grid points
Definition ignis.h:33
ignis::TRbc
double TRbc
T boundary values: left and right (as needed)
Definition ignis.h:50
ignis::L
double L
domain size (m)
Definition ignis.h:56
ignis::hscale
double hscale
scaling value for enthalpy (for solvers)
Definition ignis.h:63
ignis::hspRbc
std::vector< double > hspRbc
species enthalpies on right boundary
Definition ignis.h:54
ignis::hLbc
double hLbc
Definition ignis.h:51
ignis::Ia
size_t Ia(size_t i, size_t k)
Definition ignis.h:143
ignis::awts_sur
std::vector< double > awts_sur
awts for surroundings
Definition ignis.h:79
ignis::kin
std::shared_ptr< Cantera::Kinetics > kin
Cantera kinetics object.
Definition ignis.h:71
ignis::chi0
double chi0
Definition ignis.h:93
ignis::fl
std::vector< double > fl
fractions for interpolation
Definition ignis.h:59
ignis::s
double s
homotopy variable
Definition ignis.h:68
ignis::gas
std::shared_ptr< Cantera::ThermoPhase > gas
Cantera thermo object.
Definition ignis.h:70
ignis::writeFile
void writeFile(const std::string fname)
Definition ignis.cc:345
ignis::doLe1
bool doLe1
true if doing unity Lewis numbers (default false)
Definition ignis.h:81
ignis::yRbc
std::vector< double > yRbc
y boundary values: left and right (as needed)
Definition ignis.h:49
ignis::cpRbc
double cpRbc
cp boundary values: left and right (as needed)
Definition ignis.h:52
ignis::cpLbc
double cpLbc
Definition ignis.h:52
ignis::isFlamelet
bool isFlamelet
true for laminar flamelet (mixture fraction coordinate)
Definition ignis.h:91
ignis::sootScales
std::vector< double > sootScales
scaling value for soot variables (for solvers)
Definition ignis.h:64
ignis::strm
std::shared_ptr< streams > strm
streams object (mixture fraction, etc.)
Definition ignis.h:74
ignis::fr
std::vector< double > fr
fractions for interpolation
Definition ignis.h:60
ignis::doEnergyEqn
bool doEnergyEqn
for premixed flames: can solve energy equation or set T profile
Definition ignis.h:98
ignis::Tscale
double Tscale
scaling value for temperature (for solvers)
Definition ignis.h:62
ignis::nvar
size_t nvar
number of transported variables at each grid point
Definition ignis.h:36
ignis::setFluxesUnity
void setFluxesUnity()
Definition ignis.cc:590
ignis::setDerivative2
void setDerivative2(const double vL, const double vR, const std::vector< double > &v, std::vector< double > &d2vdx2)
Definition ignis.cc:1568
ignis::rhsf_flamelet
int rhsf_flamelet(const double *vars, double *dvarsdt)
Definition ignis.cc:1317
ignis::sootstore
std::vector< std::vector< double > > sootstore
stored soot variables
Definition ignis.h:47
ignis::dT
double dT
delta T increment for unsteady cases
Definition ignis.h:84
ignis::writeFileHdf5
void writeFileHdf5(const std::string gname, const std::string timeType)
Definition ignis.cc:264
ignis::x
std::vector< double > x
grid position values (m)
Definition ignis.h:57
ignis::fh5
std::shared_ptr< HighFive::File > fh5
hdf5 file pointer
Definition ignis.h:111
ignis::flux_y
std::vector< std::vector< double > > flux_y
species fluxes: [I(igrid, ksp)] I(igrid,ksp) maps 2D onto 1D
Definition ignis.h:87
ignis::Func
int Func(const double *vars, double *F)
Definition ignis.cc:1008
ignis::setFluxes
void setFluxes()
Definition ignis.cc:738
ignis::trn
std::shared_ptr< Cantera::Transport > trn
Cantera transport object.
Definition ignis.h:72
ignis::SM
std::shared_ptr< soot::sootModel > SM
soot model
Definition ignis.h:106
ignis::hspLbc
std::vector< double > hspLbc
species enthalpies on left boundary
Definition ignis.h:53
ignis::yLbc
std::vector< double > yLbc
Definition ignis.h:49
ignis::TLbc
double TLbc
Definition ignis.h:50
ignis::P
double P
system pressure, uniform (Pa)
Definition ignis.h:39
ignis::LI
std::shared_ptr< linearInterp > LI
interpolator for specified temperature profiles
Definition ignis.h:99
ignis::ignis
ignis(const bool _isPremixed, const bool _doEnergyEqn, const bool _isFlamelet, const bool _doSoot, const size_t _ngrd, const double _L, const double _P, std::shared_ptr< Cantera::Solution > csol, std::string _radType, const std::vector< double > &_yLbc, const std::vector< double > &_yRbc, const double _TLbc, const double _TRbc, std::shared_ptr< soot::sootModel > _SM, std::shared_ptr< soot::state > _SMstate)
Definition ignis.cc:45
ignis::y
std::vector< std::vector< double > > y
mass fractions: y[igrid][isp]
Definition ignis.h:40
ignis::ystore
std::vector< std::vector< double > > ystore
stored mass fractions
Definition ignis.h:45
ignis::storeState
void storeState()
Definition ignis.cc:469
ignis::setDerivative
void setDerivative(const double vL, const double vR, const std::vector< double > &v, std::vector< double > &dvdx)
Definition ignis.cc:1538
ignis::kabs_sur
std::vector< double > kabs_sur
kabs for surroundings
Definition ignis.h:78
ignis::radProps
std::shared_ptr< rad > radProps
radiation object
Definition ignis.h:75
ignis::Ttarget
double Ttarget
for unsteady cases, run until this max T instead of for a given time
Definition ignis.h:83
ignis::flux_h
std::vector< double > flux_h
species fluxes: [igrid]
Definition ignis.h:89
ignis::setChi
void setChi(const double _chi0)
Definition ignis.cc:1607
ignis::mflux
double mflux
premixed flame mass flux (kg/m2*s)
Definition ignis.h:96
ignis::SMstate
std::shared_ptr< soot::state > SMstate
holds state variables (gas and soot) for soot model
Definition ignis.h:107
ignis::sootvars
std::vector< std::vector< double > > sootvars
soot moments or sections; sootvars[igrid][isoot]
Definition ignis.h:42
ignis::isave
int isave
file counter for save during unsteady cases
Definition ignis.h:85
ignis::setIC
void setIC(const std::string icType, const std::string fname="")
Definition ignis.cc:485
ignis::setQrad
void setQrad(std::vector< double > &Q)
Definition ignis.cc:1082
ignis::vars0
std::vector< double > vars0
for homotopy approaches
Definition ignis.h:66
ignis::doSoot
bool doSoot
soot flag
Definition ignis.h:105
ignis::chi
std::vector< double > chi
dissipation rate profile
Definition ignis.h:92
ignis::Tstore
std::vector< double > Tstore
stored temperature
Definition ignis.h:46
ignis::rhsf
int rhsf(const double *vars, double *dvarsdt)
Definition ignis.cc:1215
ignis::solveUnsteady
void solveUnsteady(const double nTauRun, const int nsteps, const bool doWriteTime=true, const double Tmin=0, const double Tmax=0)
Definition ignis.cc:1122
ignis::isPremixed
bool isPremixed
true of the case is a premixed flame, (only left boundary condition, constant mass flux through domai...
Definition ignis.h:95
ignis::nsoot
size_t nsoot
number of soot variables
Definition ignis.h:35
ignis::radType
std::string radType
radiation model name
Definition ignis.h:76
ignis::solveSS
void solveSS()
Definition ignis.cc:942
ignis::Pstore
double Pstore
stored system pressure (for initializing from stored state)
Definition ignis.h:44
integrator_cvode
Definition integrator_cvode.h:19
integrator_cvode::integrate
int integrate(std::vector< double > &y, const realtype dt)
Definition integrator_cvode.h:91
solver_kinsol
Definition solver_kinsol.h:21
solver_kinsol::solve
int solve(int(*Func)(N_Vector, N_Vector, void *), std::vector< double > &y)
Definition solver_kinsol.h:109
rhsf_cvode
int rhsf_cvode(realtype t, N_Vector varsCV, N_Vector dvarsdtCV, void *user_data)
Definition ignis.cc:1638
myErfInv2
float myErfInv2(float x)
Definition ignis.cc:1594
Func_kinsol
int Func_kinsol(N_Vector varsKS, N_Vector fvec, void *user_data)
Definition ignis.cc:1063
integrator_cvode.h
solver_kinsol.h