Cantera 2.6.0
StFlow.cpp
Go to the documentation of this file.
1//! @file StFlow.cpp
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
7#include "cantera/oneD/refine.h"
8#include "cantera/base/ctml.h"
11#include "cantera/base/global.h"
12
13#include <set>
14
15using namespace std;
16
17namespace Cantera
18{
19
20StFlow::StFlow(ThermoPhase* ph, size_t nsp, size_t points) :
21 Domain1D(nsp+c_offset_Y, points),
22 m_press(-1.0),
23 m_nsp(nsp),
24 m_thermo(0),
25 m_kin(0),
26 m_trans(0),
27 m_epsilon_left(0.0),
28 m_epsilon_right(0.0),
29 m_do_soret(false),
30 m_do_multicomponent(false),
31 m_do_radiation(false),
32 m_kExcessLeft(0),
33 m_kExcessRight(0),
34 m_zfixed(Undef),
35 m_tfixed(-1.)
36{
37 if (ph->type() == "IdealGas") {
38 m_thermo = static_cast<IdealGasPhase*>(ph);
39 } else {
40 throw CanteraError("StFlow::StFlow",
41 "Unsupported phase type: need 'IdealGasPhase'");
42 }
43 m_type = cFlowType;
44 m_points = points;
45
46 if (ph == 0) {
47 return; // used to create a dummy object
48 }
49
50 size_t nsp2 = m_thermo->nSpecies();
51 if (nsp2 != m_nsp) {
52 m_nsp = nsp2;
53 Domain1D::resize(m_nsp+c_offset_Y, points);
54 }
55
56 // make a local copy of the species molecular weight vector
57 m_wt = m_thermo->molecularWeights();
58
59 // the species mass fractions are the last components in the solution
60 // vector, so the total number of components is the number of species
61 // plus the offset of the first mass fraction.
62 m_nv = c_offset_Y + m_nsp;
63
64 // enable all species equations by default
65 m_do_species.resize(m_nsp, true);
66
67 // but turn off the energy equation at all points
68 m_do_energy.resize(m_points,false);
69
70 m_diff.resize(m_nsp*m_points);
71 m_multidiff.resize(m_nsp*m_nsp*m_points);
72 m_flux.resize(m_nsp,m_points);
73 m_wdot.resize(m_nsp,m_points, 0.0);
74 m_ybar.resize(m_nsp);
75 m_qdotRadiation.resize(m_points, 0.0);
76
77 //-------------- default solution bounds --------------------
78 setBounds(0, -1e20, 1e20); // no bounds on u
79 setBounds(1, -1e20, 1e20); // V
80 setBounds(2, 200.0, 2*m_thermo->maxTemp()); // temperature bounds
81 setBounds(3, -1e20, 1e20); // lambda should be negative
82 setBounds(c_offset_E, -1e20, 1e20); // no bounds for inactive component
83
84 // mass fraction bounds
85 for (size_t k = 0; k < m_nsp; k++) {
86 setBounds(c_offset_Y+k, -1.0e-7, 1.0e5);
87 }
88
89 //-------------------- grid refinement -------------------------
90 m_refiner->setActive(c_offset_U, false);
91 m_refiner->setActive(c_offset_V, false);
92 m_refiner->setActive(c_offset_T, false);
93 m_refiner->setActive(c_offset_L, false);
94
95 vector_fp gr;
96 for (size_t ng = 0; ng < m_points; ng++) {
97 gr.push_back(1.0*ng/m_points);
98 }
99 setupGrid(m_points, gr.data());
100
101 // Find indices for radiating species
102 m_kRadiating.resize(2, npos);
103 m_kRadiating[0] = m_thermo->speciesIndex("CO2");
104 m_kRadiating[1] = m_thermo->speciesIndex("H2O");
105}
106
107void StFlow::resize(size_t ncomponents, size_t points)
108{
109 Domain1D::resize(ncomponents, points);
110 m_rho.resize(m_points, 0.0);
111 m_wtm.resize(m_points, 0.0);
112 m_cp.resize(m_points, 0.0);
113 m_visc.resize(m_points, 0.0);
114 m_tcon.resize(m_points, 0.0);
115
116 m_diff.resize(m_nsp*m_points);
117 if (m_do_multicomponent) {
118 m_multidiff.resize(m_nsp*m_nsp*m_points);
119 m_dthermal.resize(m_nsp, m_points, 0.0);
120 }
121 m_flux.resize(m_nsp,m_points);
122 m_wdot.resize(m_nsp,m_points, 0.0);
123 m_do_energy.resize(m_points,false);
124 m_qdotRadiation.resize(m_points, 0.0);
125 m_fixedtemp.resize(m_points);
126
127 m_dz.resize(m_points-1);
128 m_z.resize(m_points);
129}
130
131void StFlow::setupGrid(size_t n, const doublereal* z)
132{
133 resize(m_nv, n);
134
135 m_z[0] = z[0];
136 for (size_t j = 1; j < m_points; j++) {
137 if (z[j] <= z[j-1]) {
138 throw CanteraError("StFlow::setupGrid",
139 "grid points must be monotonically increasing");
140 }
141 m_z[j] = z[j];
142 m_dz[j-1] = m_z[j] - m_z[j-1];
143 }
144}
145
147{
148 double* x = xg + loc();
149 for (size_t j = 0; j < m_points; j++) {
150 double* Y = x + m_nv*j + c_offset_Y;
151 m_thermo->setMassFractions(Y);
152 m_thermo->getMassFractions(Y);
153 }
154}
155
157{
158 m_trans = &trans;
159 m_do_multicomponent = (m_trans->transportType() == "Multi" || m_trans->transportType() == "CK_Multi");
160
161 m_diff.resize(m_nsp*m_points);
162 if (m_do_multicomponent) {
163 m_multidiff.resize(m_nsp*m_nsp*m_points);
164 m_dthermal.resize(m_nsp, m_points, 0.0);
165 }
166}
167
169{
170 for (size_t j = 0; j < m_points; j++) {
171 T(x,j) = m_thermo->temperature();
172 m_thermo->getMassFractions(&Y(x, 0, j));
173 }
174}
175
176void StFlow::setGas(const doublereal* x, size_t j)
177{
178 m_thermo->setTemperature(T(x,j));
179 const doublereal* yy = x + m_nv*j + c_offset_Y;
180 m_thermo->setMassFractions_NoNorm(yy);
181 m_thermo->setPressure(m_press);
182}
183
184void StFlow::setGasAtMidpoint(const doublereal* x, size_t j)
185{
186 m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1)));
187 const doublereal* yyj = x + m_nv*j + c_offset_Y;
188 const doublereal* yyjp = x + m_nv*(j+1) + c_offset_Y;
189 for (size_t k = 0; k < m_nsp; k++) {
190 m_ybar[k] = 0.5*(yyj[k] + yyjp[k]);
191 }
192 m_thermo->setMassFractions_NoNorm(m_ybar.data());
193 m_thermo->setPressure(m_press);
194}
195
196void StFlow::_finalize(const doublereal* x)
197{
198 if (!m_do_multicomponent && m_do_soret) {
199 throw CanteraError("StFlow::_finalize",
200 "Thermal diffusion (the Soret effect) is enabled, and requires "
201 "using a multicomponent transport model.");
202 }
203
204 size_t nz = m_zfix.size();
205 bool e = m_do_energy[0];
206 for (size_t j = 0; j < m_points; j++) {
207 if (e || nz == 0) {
208 m_fixedtemp[j] = T(x, j);
209 } else {
210 double zz = (z(j) - z(0))/(z(m_points - 1) - z(0));
211 double tt = linearInterp(zz, m_zfix, m_tfix);
212 m_fixedtemp[j] = tt;
213 }
214 }
215 if (e) {
216 solveEnergyEqn();
217 }
218
219 if (domainType() == cFreeFlow) {
220 // If the domain contains the temperature fixed point, make sure that it
221 // is correctly set. This may be necessary when the grid has been modified
222 // externally.
223 if (m_tfixed != Undef) {
224 for (size_t j = 0; j < m_points; j++) {
225 if (z(j) == m_zfixed) {
226 return; // fixed point is already set correctly
227 }
228 }
229
230 for (size_t j = 0; j < m_points - 1; j++) {
231 // Find where the temperature profile crosses the current
232 // fixed temperature.
233 if ((T(x, j) - m_tfixed) * (T(x, j+1) - m_tfixed) <= 0.0) {
234 m_tfixed = T(x, j+1);
235 m_zfixed = z(j+1);
236 return;
237 }
238 }
239 }
240 }
241}
242
243void StFlow::eval(size_t jg, doublereal* xg,
244 doublereal* rg, integer* diagg, doublereal rdt)
245{
246 // if evaluating a Jacobian, and the global point is outside the domain of
247 // influence for this domain, then skip evaluating the residual
248 if (jg != npos && (jg + 1 < firstPoint() || jg > lastPoint() + 1)) {
249 return;
250 }
251
252 // start of local part of global arrays
253 doublereal* x = xg + loc();
254 doublereal* rsd = rg + loc();
255 integer* diag = diagg + loc();
256
257 size_t jmin, jmax;
258 if (jg == npos) { // evaluate all points
259 jmin = 0;
260 jmax = m_points - 1;
261 } else { // evaluate points for Jacobian
262 size_t jpt = (jg == 0) ? 0 : jg - firstPoint();
263 jmin = std::max<size_t>(jpt, 1) - 1;
264 jmax = std::min(jpt+1,m_points-1);
265 }
266
267 updateProperties(jg, x, jmin, jmax);
268 evalResidual(x, rsd, diag, rdt, jmin, jmax);
269}
270
271void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
272{
273 // properties are computed for grid points from j0 to j1
274 size_t j0 = std::max<size_t>(jmin, 1) - 1;
275 size_t j1 = std::min(jmax+1,m_points-1);
276
277 updateThermo(x, j0, j1);
278 if (jg == npos || m_force_full_update) {
279 // update transport properties only if a Jacobian is not being
280 // evaluated, or if specifically requested
281 updateTransport(x, j0, j1);
282 }
283 if (jg == npos) {
284 double* Yleft = x + index(c_offset_Y, jmin);
285 m_kExcessLeft = distance(Yleft, max_element(Yleft, Yleft + m_nsp));
286 double* Yright = x + index(c_offset_Y, jmax);
287 m_kExcessRight = distance(Yright, max_element(Yright, Yright + m_nsp));
288 }
289
290 // update the species diffusive mass fluxes whether or not a
291 // Jacobian is being evaluated
292 updateDiffFluxes(x, j0, j1);
293}
294
295void StFlow::evalResidual(double* x, double* rsd, int* diag,
296 double rdt, size_t jmin, size_t jmax)
297{
298 //----------------------------------------------------
299 // evaluate the residual equations at all required
300 // grid points
301 //----------------------------------------------------
302
303 // calculation of qdotRadiation
304
305 // The simple radiation model used was established by Y. Liu and B. Rogg [Y.
306 // Liu and B. Rogg, Modelling of thermally radiating diffusion flames with
307 // detailed chemistry and transport, EUROTHERM Seminars, 17:114-127, 1991].
308 // This model uses the optically thin limit and the gray-gas approximation
309 // to simply calculate a volume specified heat flux out of the Planck
310 // absorption coefficients, the boundary emissivities and the temperature.
311 // The model considers only CO2 and H2O as radiating species. Polynomial
312 // lines calculate the species Planck coefficients for H2O and CO2. The data
313 // for the lines is taken from the RADCAL program [Grosshandler, W. L.,
314 // RADCAL: A Narrow-Band Model for Radiation Calculations in a Combustion
315 // Environment, NIST technical note 1402, 1993]. The coefficients for the
316 // polynomials are taken from [http://www.sandia.gov/TNF/radiation.html].
317
318 if (m_do_radiation) {
319 // variable definitions for the Planck absorption coefficient and the
320 // radiation calculation:
321 doublereal k_P_ref = 1.0*OneAtm;
322
323 // polynomial coefficients:
324 const doublereal c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880,
325 0.51382, -1.86840e-5};
326 const doublereal c_CO2[6] = {18.741, -121.310, 273.500, -194.050,
327 56.310, -5.8169};
328
329 // calculation of the two boundary values
330 double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4);
331 double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4);
332
333 // loop over all grid points
334 for (size_t j = jmin; j < jmax; j++) {
335 // helping variable for the calculation
336 double radiative_heat_loss = 0;
337
338 // calculation of the mean Planck absorption coefficient
339 double k_P = 0;
340 // absorption coefficient for H2O
341 if (m_kRadiating[1] != npos) {
342 double k_P_H2O = 0;
343 for (size_t n = 0; n <= 5; n++) {
344 k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n);
345 }
346 k_P_H2O /= k_P_ref;
347 k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O;
348 }
349 // absorption coefficient for CO2
350 if (m_kRadiating[0] != npos) {
351 double k_P_CO2 = 0;
352 for (size_t n = 0; n <= 5; n++) {
353 k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n);
354 }
355 k_P_CO2 /= k_P_ref;
356 k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2;
357 }
358
359 // calculation of the radiative heat loss term
360 radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4)
361 - boundary_Rad_left - boundary_Rad_right);
362
363 // set the radiative heat loss vector
364 m_qdotRadiation[j] = radiative_heat_loss;
365 }
366 }
367
368 for (size_t j = jmin; j <= jmax; j++) {
369 //----------------------------------------------
370 // left boundary
371 //----------------------------------------------
372
373 if (j == 0) {
374 // these may be modified by a boundary object
375
376 // Continuity. This propagates information right-to-left, since
377 // rho_u at point 0 is dependent on rho_u at point 1, but not on
378 // mdot from the inlet.
379 rsd[index(c_offset_U,0)] =
380 -(rho_u(x,1) - rho_u(x,0))/m_dz[0]
381 -(density(1)*V(x,1) + density(0)*V(x,0));
382
383 // the inlet (or other) object connected to this one will modify
384 // these equations by subtracting its values for V, T, and mdot. As
385 // a result, these residual equations will force the solution
386 // variables to the values for the boundary object
387 rsd[index(c_offset_V,0)] = V(x,0);
388 if (doEnergy(0)) {
389 rsd[index(c_offset_T,0)] = T(x,0);
390 } else {
391 rsd[index(c_offset_T,0)] = T(x,0) - T_fixed(0);
392 }
393 rsd[index(c_offset_L,0)] = -rho_u(x,0);
394
395 // The default boundary condition for species is zero flux. However,
396 // the boundary object may modify this.
397 double sum = 0.0;
398 for (size_t k = 0; k < m_nsp; k++) {
399 sum += Y(x,k,0);
400 rsd[index(c_offset_Y + k, 0)] =
401 -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
402 }
403 rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum;
404
405 // set residual of poisson's equ to zero
406 rsd[index(c_offset_E, 0)] = x[index(c_offset_E, j)];
407 } else if (j == m_points - 1) {
408 evalRightBoundary(x, rsd, diag, rdt);
409 // set residual of poisson's equ to zero
410 rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
411 } else { // interior points
412 evalContinuity(j, x, rsd, diag, rdt);
413 // set residual of poisson's equ to zero
414 rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
415
416 //------------------------------------------------
417 // Radial momentum equation
418 //
419 // \rho dV/dt + \rho u dV/dz + \rho V^2
420 // = d(\mu dV/dz)/dz - lambda
421 //-------------------------------------------------
422 rsd[index(c_offset_V,j)]
423 = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j)
424 - m_rho[j]*V(x,j)*V(x,j))/m_rho[j]
425 - rdt*(V(x,j) - V_prev(j));
426 diag[index(c_offset_V, j)] = 1;
427
428 //-------------------------------------------------
429 // Species equations
430 //
431 // \rho dY_k/dt + \rho u dY_k/dz + dJ_k/dz
432 // = M_k\omega_k
433 //-------------------------------------------------
434 getWdot(x,j);
435 for (size_t k = 0; k < m_nsp; k++) {
436 double convec = rho_u(x,j)*dYdz(x,k,j);
437 double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))
438 / (z(j+1) - z(j-1));
439 rsd[index(c_offset_Y + k, j)]
440 = (m_wt[k]*(wdot(k,j))
441 - convec - diffus)/m_rho[j]
442 - rdt*(Y(x,k,j) - Y_prev(k,j));
443 diag[index(c_offset_Y + k, j)] = 1;
444 }
445
446 //-----------------------------------------------
447 // energy equation
448 //
449 // \rho c_p dT/dt + \rho c_p u dT/dz
450 // = d(k dT/dz)/dz
451 // - sum_k(\omega_k h_k_ref)
452 // - sum_k(J_k c_p_k / M_k) dT/dz
453 //-----------------------------------------------
454 if (m_do_energy[j]) {
455 setGas(x,j);
456
457 // heat release term
458 const vector_fp& h_RT = m_thermo->enthalpy_RT_ref();
459 const vector_fp& cp_R = m_thermo->cp_R_ref();
460 double sum = 0.0;
461 double sum2 = 0.0;
462 for (size_t k = 0; k < m_nsp; k++) {
463 double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
464 sum += wdot(k,j)*h_RT[k];
465 sum2 += flxk*cp_R[k]/m_wt[k];
466 }
467 sum *= GasConstant * T(x,j);
468 double dtdzj = dTdz(x,j);
469 sum2 *= GasConstant * dtdzj;
470
471 rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj
472 - divHeatFlux(x,j) - sum - sum2;
473 rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
474 rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
475 rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j]));
476 diag[index(c_offset_T, j)] = 1;
477 } else {
478 // residual equations if the energy equation is disabled
479 rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
480 diag[index(c_offset_T, j)] = 0;
481 }
482
483 rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
484 diag[index(c_offset_L, j)] = 0;
485 }
486 }
487}
488
489void StFlow::updateTransport(doublereal* x, size_t j0, size_t j1)
490{
491 if (m_do_multicomponent) {
492 for (size_t j = j0; j < j1; j++) {
493 setGasAtMidpoint(x,j);
494 doublereal wtm = m_thermo->meanMolecularWeight();
495 doublereal rho = m_thermo->density();
496 m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
497 m_trans->getMultiDiffCoeffs(m_nsp, &m_multidiff[mindex(0,0,j)]);
498
499 // Use m_diff as storage for the factor outside the summation
500 for (size_t k = 0; k < m_nsp; k++) {
501 m_diff[k+j*m_nsp] = m_wt[k] * rho / (wtm*wtm);
502 }
503
504 m_tcon[j] = m_trans->thermalConductivity();
505 if (m_do_soret) {
506 m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + j*m_nsp);
507 }
508 }
509 } else { // mixture averaged transport
510 for (size_t j = j0; j < j1; j++) {
511 setGasAtMidpoint(x,j);
512 m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
513 m_trans->getMixDiffCoeffs(&m_diff[j*m_nsp]);
514 m_tcon[j] = m_trans->thermalConductivity();
515 }
516 }
517}
518
519void StFlow::showSolution(const doublereal* x)
520{
521 writelog(" Pressure: {:10.4g} Pa\n", m_press);
522
524
525 if (m_do_radiation) {
526 writeline('-', 79, false, true);
527 writelog("\n z radiative heat loss");
528 writeline('-', 79, false, true);
529 for (size_t j = 0; j < m_points; j++) {
530 writelog("\n {:10.4g} {:10.4g}", m_z[j], m_qdotRadiation[j]);
531 }
532 writelog("\n");
533 }
534}
535
536void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
537{
538 if (m_do_multicomponent) {
539 for (size_t j = j0; j < j1; j++) {
540 double dz = z(j+1) - z(j);
541 for (size_t k = 0; k < m_nsp; k++) {
542 doublereal sum = 0.0;
543 for (size_t m = 0; m < m_nsp; m++) {
544 sum += m_wt[m] * m_multidiff[mindex(k,m,j)] * (X(x,m,j+1)-X(x,m,j));
545 }
546 m_flux(k,j) = sum * m_diff[k+j*m_nsp] / dz;
547 }
548 }
549 } else {
550 for (size_t j = j0; j < j1; j++) {
551 double sum = 0.0;
552 double wtm = m_wtm[j];
553 double rho = density(j);
554 double dz = z(j+1) - z(j);
555 for (size_t k = 0; k < m_nsp; k++) {
556 m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
557 m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
558 sum -= m_flux(k,j);
559 }
560 // correction flux to insure that \sum_k Y_k V_k = 0.
561 for (size_t k = 0; k < m_nsp; k++) {
562 m_flux(k,j) += sum*Y(x,k,j);
563 }
564 }
565 }
566
567 if (m_do_soret) {
568 for (size_t m = j0; m < j1; m++) {
569 double gradlogT = 2.0 * (T(x,m+1) - T(x,m)) /
570 ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m)));
571 for (size_t k = 0; k < m_nsp; k++) {
572 m_flux(k,m) -= m_dthermal(k,m)*gradlogT;
573 }
574 }
575 }
576}
577
578string StFlow::componentName(size_t n) const
579{
580 switch (n) {
581 case 0:
582 return "velocity";
583 case 1:
584 return "spread_rate";
585 case 2:
586 return "T";
587 case 3:
588 return "lambda";
589 case 4:
590 return "eField";
591 default:
592 if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) {
593 return m_thermo->speciesName(n - c_offset_Y);
594 } else {
595 return "<unknown>";
596 }
597 }
598}
599
600size_t StFlow::componentIndex(const std::string& name) const
601{
602 if (name=="velocity") {
603 return 0;
604 } else if (name=="spread_rate") {
605 return 1;
606 } else if (name=="T") {
607 return 2;
608 } else if (name=="lambda") {
609 return 3;
610 } else if (name == "eField") {
611 return 4;
612 } else {
613 for (size_t n=c_offset_Y; n<m_nsp+c_offset_Y; n++) {
614 if (componentName(n)==name) {
615 return n;
616 }
617 }
618 throw CanteraError("StFlow1D::componentIndex",
619 "no component named " + name);
620 }
621}
622
623bool StFlow::componentActive(size_t n) const
624{
625 switch (n) {
626 case c_offset_V: // spread_rate
627 return m_type != cFreeFlow;
628 case c_offset_L: // lambda
629 return m_type != cFreeFlow;
630 case c_offset_E: // eField
631 return false;
632 default:
633 return true;
634 }
635}
636
637void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel)
638{
639 Domain1D::restore(dom, soln, loglevel);
640 vector<string> ignored;
641 size_t nsp = m_thermo->nSpecies();
642 vector_int did_species(nsp, 0);
643
644 vector<XML_Node*> str = dom.getChildren("string");
645 for (size_t istr = 0; istr < str.size(); istr++) {
646 const XML_Node& nd = *str[istr];
647 writelog(nd["title"]+": "+nd.value()+"\n");
648 }
649
650 double pp = getFloat(dom, "pressure", "pressure");
651 setPressure(pp);
652 vector<XML_Node*> d = dom.child("grid_data").getChildren("floatArray");
653 vector_fp x;
654 size_t np = 0;
655 bool readgrid = false, wrote_header = false;
656 for (size_t n = 0; n < d.size(); n++) {
657 const XML_Node& fa = *d[n];
658 string nm = fa["title"];
659 if (nm == "z") {
660 getFloatArray(fa,x,false);
661 np = x.size();
662 if (loglevel >= 2) {
663 writelog("Grid contains {} points.\n", np);
664 }
665 readgrid = true;
666 setupGrid(np, x.data());
667 }
668 }
669 if (!readgrid) {
670 throw CanteraError("StFlow::restore",
671 "domain contains no grid points.");
672 }
673
674 debuglog("Importing datasets:\n", loglevel >= 2);
675 for (size_t n = 0; n < d.size(); n++) {
676 const XML_Node& fa = *d[n];
677 string nm = fa["title"];
678 getFloatArray(fa,x,false);
679 if (nm == "u") {
680 debuglog("axial velocity ", loglevel >= 2);
681 if (x.size() != np) {
682 throw CanteraError("StFlow::restore",
683 "axial velocity array size error");
684 }
685 for (size_t j = 0; j < np; j++) {
686 soln[index(c_offset_U,j)] = x[j];
687 }
688 } else if (nm == "z") {
689 ; // already read grid
690 } else if (nm == "V") {
691 debuglog("radial velocity ", loglevel >= 2);
692 if (x.size() != np) {
693 throw CanteraError("StFlow::restore",
694 "radial velocity array size error");
695 }
696 for (size_t j = 0; j < np; j++) {
697 soln[index(c_offset_V,j)] = x[j];
698 }
699 } else if (nm == "T") {
700 debuglog("temperature ", loglevel >= 2);
701 if (x.size() != np) {
702 throw CanteraError("StFlow::restore",
703 "temperature array size error");
704 }
705 for (size_t j = 0; j < np; j++) {
706 soln[index(c_offset_T,j)] = x[j];
707 }
708
709 // For fixed-temperature simulations, use the imported temperature
710 // profile by default. If this is not desired, call
711 // setFixedTempProfile *after* restoring the solution.
712 vector_fp zz(np);
713 for (size_t jj = 0; jj < np; jj++) {
714 zz[jj] = (grid(jj) - zmin())/(zmax() - zmin());
715 }
716 setFixedTempProfile(zz, x);
717 } else if (nm == "L") {
718 debuglog("lambda ", loglevel >= 2);
719 if (x.size() != np) {
720 throw CanteraError("StFlow::restore",
721 "lambda array size error");
722 }
723 for (size_t j = 0; j < np; j++) {
724 soln[index(c_offset_L,j)] = x[j];
725 }
726 } else if (m_thermo->speciesIndex(nm) != npos) {
727 debuglog(nm+" ", loglevel >= 2);
728 if (x.size() == np) {
729 size_t k = m_thermo->speciesIndex(nm);
730 did_species[k] = 1;
731 for (size_t j = 0; j < np; j++) {
732 soln[index(k+c_offset_Y,j)] = x[j];
733 }
734 }
735 } else {
736 ignored.push_back(nm);
737 }
738 }
739
740 if (loglevel >=2 && !ignored.empty()) {
741 writelog("\n\n");
742 writelog("Ignoring datasets:\n");
743 size_t nn = ignored.size();
744 for (size_t n = 0; n < nn; n++) {
745 writelog(ignored[n]+" ");
746 }
747 }
748
749 if (loglevel >= 1) {
750 for (size_t ks = 0; ks < nsp; ks++) {
751 if (did_species[ks] == 0) {
752 if (!wrote_header) {
753 writelog("Missing data for species:\n");
754 wrote_header = true;
755 }
756 writelog(m_thermo->speciesName(ks)+" ");
757 }
758 }
759 }
760
761 if (dom.hasChild("energy_enabled")) {
762 getFloatArray(dom, x, false, "", "energy_enabled");
763 if (x.size() == nPoints()) {
764 for (size_t i = 0; i < x.size(); i++) {
765 m_do_energy[i] = (x[i] != 0);
766 }
767 } else if (!x.empty()) {
768 throw CanteraError("StFlow::restore", "energy_enabled is length {}"
769 "but should be length {}", x.size(), nPoints());
770 }
771 }
772
773 if (dom.hasChild("species_enabled")) {
774 getFloatArray(dom, x, false, "", "species_enabled");
775 if (x.size() == m_nsp) {
776 for (size_t i = 0; i < x.size(); i++) {
777 m_do_species[i] = (x[i] != 0);
778 }
779 } else if (!x.empty()) {
780 // This may occur when restoring from a mechanism with a different
781 // number of species.
782 if (loglevel > 0) {
783 warn_user("StFlow::restore", "species_enabled is "
784 "length {} but should be length {}. Enabling all species "
785 "equations by default.", x.size(), m_nsp);
786 }
787 m_do_species.assign(m_nsp, true);
788 }
789 }
790
791 if (dom.hasChild("refine_criteria")) {
792 XML_Node& ref = dom.child("refine_criteria");
793 refiner().setCriteria(getFloat(ref, "ratio"), getFloat(ref, "slope"),
794 getFloat(ref, "curve"), getFloat(ref, "prune"));
795 refiner().setGridMin(getFloat(ref, "grid_min"));
796 }
797
798 if (domainType() == cFreeFlow) {
799 getOptionalFloat(dom, "t_fixed", m_tfixed);
800 getOptionalFloat(dom, "z_fixed", m_zfixed);
801 }
802}
803
804XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol)
805{
806 Array2D soln(m_nv, m_points, sol + loc());
807 XML_Node& flow = Domain1D::save(o, sol);
808 flow.addAttribute("type",flowType());
809
810 XML_Node& gv = flow.addChild("grid_data");
811 addFloat(flow, "pressure", m_press, "Pa", "pressure");
812
813 addFloatArray(gv,"z",m_z.size(), m_z.data(),
814 "m","length");
815 vector_fp x(soln.nColumns());
816
817 soln.getRow(c_offset_U, x.data());
818 addFloatArray(gv,"u",x.size(),x.data(),"m/s","velocity");
819
820 soln.getRow(c_offset_V, x.data());
821 addFloatArray(gv,"V",x.size(),x.data(),"1/s","rate");
822
823 soln.getRow(c_offset_T, x.data());
824 addFloatArray(gv,"T",x.size(),x.data(),"K","temperature");
825
826 soln.getRow(c_offset_L, x.data());
827 addFloatArray(gv,"L",x.size(),x.data(),"N/m^4");
828
829 for (size_t k = 0; k < m_nsp; k++) {
830 soln.getRow(c_offset_Y+k, x.data());
831 addFloatArray(gv,m_thermo->speciesName(k),
832 x.size(),x.data(),"","massFraction");
833 }
834 if (m_do_radiation) {
835 addFloatArray(gv, "radiative_heat_loss", m_z.size(),
836 m_qdotRadiation.data(), "W/m^3", "specificPower");
837 }
838 vector_fp values(nPoints());
839 for (size_t i = 0; i < nPoints(); i++) {
840 values[i] = m_do_energy[i];
841 }
842 addNamedFloatArray(flow, "energy_enabled", nPoints(), &values[0]);
843
844 values.resize(m_nsp);
845 for (size_t i = 0; i < m_nsp; i++) {
846 values[i] = m_do_species[i];
847 }
848 addNamedFloatArray(flow, "species_enabled", m_nsp, &values[0]);
849
850 XML_Node& ref = flow.addChild("refine_criteria");
851 addFloat(ref, "ratio", refiner().maxRatio());
852 addFloat(ref, "slope", refiner().maxDelta());
853 addFloat(ref, "curve", refiner().maxSlope());
854 addFloat(ref, "prune", refiner().prune());
855 addFloat(ref, "grid_min", refiner().gridMin());
856 if (m_zfixed != Undef) {
857 addFloat(flow, "z_fixed", m_zfixed, "m");
858 addFloat(flow, "t_fixed", m_tfixed, "K");
859 }
860 return flow;
861}
862
863AnyMap StFlow::serialize(const double* soln) const
864{
865 AnyMap state = Domain1D::serialize(soln);
866 state["type"] = flowType();
867 state["pressure"] = m_press;
868
869 state["phase"]["name"] = m_thermo->name();
870 AnyValue source = m_thermo->input().getMetadata("filename");
871 state["phase"]["source"] = source.empty() ? "<unknown>" : source.asString();
872
873 state["radiation-enabled"] = m_do_radiation;
874 if (m_do_radiation) {
875 state["radiative-heat-loss"] = m_qdotRadiation;
876 state["emissivity-left"] = m_epsilon_left;
877 state["emissivity-right"] = m_epsilon_right;
878 }
879
880 std::set<bool> energy_flags(m_do_energy.begin(), m_do_energy.end());
881 if (energy_flags.size() == 1) {
882 state["energy-enabled"] = m_do_energy[0];
883 } else {
884 state["energy-enabled"] = m_do_energy;
885 }
886
887 state["Soret-enabled"] = m_do_soret;
888
889 std::set<bool> species_flags(m_do_species.begin(), m_do_species.end());
890 if (species_flags.size() == 1) {
891 state["species-enabled"] = m_do_species[0];
892 } else {
893 for (size_t k = 0; k < m_nsp; k++) {
894 state["species-enabled"][m_thermo->speciesName(k)] = m_do_species[k];
895 }
896 }
897
898 state["refine-criteria"]["ratio"] = m_refiner->maxRatio();
899 state["refine-criteria"]["slope"] = m_refiner->maxDelta();
900 state["refine-criteria"]["curve"] = m_refiner->maxSlope();
901 state["refine-criteria"]["prune"] = m_refiner->prune();
902 state["refine-criteria"]["grid-min"] = m_refiner->gridMin();
903 state["refine-criteria"]["max-points"] =
904 static_cast<long int>(m_refiner->maxPoints());
905
906 if (m_zfixed != Undef) {
907 state["fixed-point"]["location"] = m_zfixed;
908 state["fixed-point"]["temperature"] = m_tfixed;
909 }
910
911 state["grid"] = m_z;
912 vector_fp data(nPoints());
913 for (size_t i = 0; i < nComponents(); i++) {
914 if (componentActive(i)) {
915 for (size_t j = 0; j < nPoints(); j++) {
916 data[j] = soln[index(i,j)];
917 }
918 state[componentName(i)] = data;
919 }
920 }
921
922 return state;
923}
924
925void StFlow::restore(const AnyMap& state, double* soln, int loglevel)
926{
927 Domain1D::restore(state, soln, loglevel);
928 m_press = state["pressure"].asDouble();
929 setupGrid(nPoints(), state["grid"].asVector<double>(nPoints()).data());
930
931 for (size_t i = 0; i < nComponents(); i++) {
932 if (!componentActive(i)) {
933 continue;
934 }
935 std::string name = componentName(i);
936 if (state.hasKey(name)) {
937 const vector_fp& data = state[name].asVector<double>(nPoints());
938 for (size_t j = 0; j < nPoints(); j++) {
939 soln[index(i,j)] = data[j];
940 }
941 } else if (loglevel) {
942 warn_user("StFlow::restore", "Saved state does not contain values for "
943 "component '{}' in domain '{}'.", name, id());
944 }
945 }
946
947 if (state.hasKey("energy-enabled")) {
948 const AnyValue& ee = state["energy-enabled"];
949 if (ee.isScalar()) {
950 m_do_energy.assign(nPoints(), ee.asBool());
951 } else {
952 m_do_energy = ee.asVector<bool>(nPoints());
953 }
954 }
955
956 if (state.hasKey("Soret-enabled")) {
957 m_do_soret = state["Soret-enabled"].asBool();
958 }
959
960 if (state.hasKey("species-enabled")) {
961 const AnyValue& se = state["species-enabled"];
962 if (se.isScalar()) {
963 m_do_species.assign(m_thermo->nSpecies(), se.asBool());
964 } else {
965 m_do_species = se.asVector<bool>(m_thermo->nSpecies());
966 }
967 }
968
969 if (state.hasKey("radiation-enabled")) {
970 m_do_radiation = state["radiation-enabled"].asBool();
971 if (m_do_radiation) {
972 m_epsilon_left = state["emissivity-left"].asDouble();
973 m_epsilon_right = state["emissivity-right"].asDouble();
974 }
975 }
976
977 if (state.hasKey("refine-criteria")) {
978 const AnyMap& criteria = state["refine-criteria"].as<AnyMap>();
979 double ratio = criteria.getDouble("ratio", m_refiner->maxRatio());
980 double slope = criteria.getDouble("slope", m_refiner->maxDelta());
981 double curve = criteria.getDouble("curve", m_refiner->maxSlope());
982 double prune = criteria.getDouble("prune", m_refiner->prune());
983 m_refiner->setCriteria(ratio, slope, curve, prune);
984
985 if (criteria.hasKey("grid-min")) {
986 m_refiner->setGridMin(criteria["grid-min"].asDouble());
987 }
988 if (criteria.hasKey("max-points")) {
989 m_refiner->setMaxPoints(criteria["max-points"].asInt());
990 }
991 }
992
993 if (state.hasKey("fixed-point")) {
994 m_zfixed = state["fixed-point"]["location"].asDouble();
995 m_tfixed = state["fixed-point"]["temperature"].asDouble();
996 }
997}
998
999void StFlow::solveEnergyEqn(size_t j)
1000{
1001 bool changed = false;
1002 if (j == npos) {
1003 for (size_t i = 0; i < m_points; i++) {
1004 if (!m_do_energy[i]) {
1005 changed = true;
1006 }
1007 m_do_energy[i] = true;
1008 }
1009 } else {
1010 if (!m_do_energy[j]) {
1011 changed = true;
1012 }
1013 m_do_energy[j] = true;
1014 }
1015 m_refiner->setActive(c_offset_U, true);
1016 m_refiner->setActive(c_offset_V, true);
1017 m_refiner->setActive(c_offset_T, true);
1018 if (changed) {
1019 needJacUpdate();
1020 }
1021}
1022
1023void StFlow::setBoundaryEmissivities(doublereal e_left, doublereal e_right)
1024{
1025 if (e_left < 0 || e_left > 1) {
1026 throw CanteraError("StFlow::setBoundaryEmissivities",
1027 "The left boundary emissivity must be between 0.0 and 1.0!");
1028 } else if (e_right < 0 || e_right > 1) {
1029 throw CanteraError("StFlow::setBoundaryEmissivities",
1030 "The right boundary emissivity must be between 0.0 and 1.0!");
1031 } else {
1032 m_epsilon_left = e_left;
1033 m_epsilon_right = e_right;
1034 }
1035}
1036
1037void StFlow::fixTemperature(size_t j)
1038{
1039 bool changed = false;
1040 if (j == npos) {
1041 for (size_t i = 0; i < m_points; i++) {
1042 if (m_do_energy[i]) {
1043 changed = true;
1044 }
1045 m_do_energy[i] = false;
1046 }
1047 } else {
1048 if (m_do_energy[j]) {
1049 changed = true;
1050 }
1051 m_do_energy[j] = false;
1052 }
1053 m_refiner->setActive(c_offset_U, false);
1054 m_refiner->setActive(c_offset_V, false);
1055 m_refiner->setActive(c_offset_T, false);
1056 if (changed) {
1057 needJacUpdate();
1058 }
1059}
1060
1061void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt)
1062{
1063 size_t j = m_points - 1;
1064
1065 // the boundary object connected to the right of this one may modify or
1066 // replace these equations. The default boundary conditions are zero u, V,
1067 // and T, and zero diffusive flux for all species.
1068
1069 rsd[index(c_offset_V,j)] = V(x,j);
1070 doublereal sum = 0.0;
1071 rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
1072 diag[index(c_offset_L, j)] = 0;
1073 for (size_t k = 0; k < m_nsp; k++) {
1074 sum += Y(x,k,j);
1075 rsd[index(k+c_offset_Y,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
1076 }
1077 rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum;
1078 diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0;
1079 if (domainType() == cAxisymmetricStagnationFlow) {
1080 rsd[index(c_offset_U,j)] = rho_u(x,j);
1081 if (m_do_energy[j]) {
1082 rsd[index(c_offset_T,j)] = T(x,j);
1083 } else {
1084 rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
1085 }
1086 } else if (domainType() == cFreeFlow) {
1087 rsd[index(c_offset_U,j)] = rho_u(x,j) - rho_u(x,j-1);
1088 rsd[index(c_offset_T,j)] = T(x,j) - T(x,j-1);
1089 }
1090}
1091
1092void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt)
1093{
1094 //algebraic constraint
1095 diag[index(c_offset_U, j)] = 0;
1096 //----------------------------------------------
1097 // Continuity equation
1098 //
1099 // d(\rho u)/dz + 2\rho V = 0
1100 //----------------------------------------------
1101 if (domainType() == cAxisymmetricStagnationFlow) {
1102 // Note that this propagates the mass flow rate information to the left
1103 // (j+1 -> j) from the value specified at the right boundary. The
1104 // lambda information propagates in the opposite direction.
1105 rsd[index(c_offset_U,j)] =
1106 -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
1107 -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
1108 } else if (domainType() == cFreeFlow) {
1109 if (grid(j) > m_zfixed) {
1110 rsd[index(c_offset_U,j)] =
1111 - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]
1112 - (density(j-1)*V(x,j-1) + density(j)*V(x,j));
1113 } else if (grid(j) == m_zfixed) {
1114 if (m_do_energy[j]) {
1115 rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
1116 } else {
1117 rsd[index(c_offset_U,j)] = (rho_u(x,j)
1118 - m_rho[0]*0.3);
1119 }
1120 } else if (grid(j) < m_zfixed) {
1121 rsd[index(c_offset_U,j)] =
1122 - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
1123 - (density(j+1)*V(x,j+1) + density(j)*V(x,j));
1124 }
1125 }
1126}
1127
1128} // namespace
Headers for the Transport object, which is the virtual base class for all transport property evaluato...
const AnyValue & getMetadata(const std::string &key) const
Get a value from the metadata applicable to the AnyMap tree containing this node.
Definition: AnyMap.cpp:577
A map of string keys to values whose type can vary at runtime.
Definition: AnyMap.h:399
double getDouble(const std::string &key, double default_) const
If key exists, return it as a double, otherwise return default_.
Definition: AnyMap.cpp:1492
bool hasKey(const std::string &key) const
Returns true if the map contains an item named key.
Definition: AnyMap.cpp:1406
A wrapper for a variable whose type is determined at runtime.
Definition: AnyMap.h:84
const std::string & asString() const
Return the held value, if it is a string.
Definition: AnyMap.cpp:716
bool & asBool()
Return the held value, if it is a bool.
Definition: AnyMap.cpp:848
bool empty() const
Return boolean indicating whether AnyValue is empty.
Definition: AnyMap.cpp:684
const std::vector< T > & asVector(size_t nMin=npos, size_t nMax=npos) const
Return the held value, if it is a vector of type T.
Definition: AnyMap.inl.h:72
bool isScalar() const
Returns true if the held value is a scalar type (such as double, long int, string,...
Definition: AnyMap.cpp:688
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:30
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:233
size_t nColumns() const
Number of columns.
Definition: Array.h:193
void getRow(size_t n, double *const rw)
Get the nth row and return it in a vector.
Definition: Array.cpp:84
void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition: Array.cpp:52
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
Base class for one-dimensional domains.
Definition: Domain1D.h:38
size_t lastPoint() const
The index of the last (that is, right-most) grid point belonging to this domain.
Definition: Domain1D.h:389
size_t nComponents() const
Number of components at each grid point.
Definition: Domain1D.h:133
virtual AnyMap serialize(const double *soln) const
Save the state of this domain as an AnyMap.
Definition: Domain1D.cpp:172
virtual void showSolution(const doublereal *x)
Print the solution.
Definition: Domain1D.cpp:261
size_t nPoints() const
Number of grid points in this domain.
Definition: Domain1D.h:155
virtual XML_Node & save(XML_Node &o, const doublereal *const sol)
Save the current solution for this domain into an XML_Node.
Definition: Domain1D.cpp:118
virtual void resize(size_t nv, size_t np)
Definition: Domain1D.cpp:41
Refiner & refiner()
Return a reference to the grid refiner.
Definition: Domain1D.h:128
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
Definition: Domain1D.cpp:131
int domainType()
Domain type flag.
Definition: Domain1D.h:53
size_t firstPoint() const
The index of the first (that is, left-most) grid point belonging to this domain.
Definition: Domain1D.h:381
void needJacUpdate()
Definition: Domain1D.cpp:110
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector,.
Definition: Domain1D.h:373
Class IdealGasPhase represents low-density gases that obey the ideal gas equation of state.
const vector_fp & cp_R_ref() const
Returns a reference to the dimensionless reference state Heat Capacity vector.
const vector_fp & enthalpy_RT_ref() const
Returns a reference to the dimensionless reference state enthalpy vector.
virtual void setPressure(doublereal p)
Set the pressure at constant temperature and composition.
std::string name() const
Return the name of the phase.
Definition: Phase.cpp:70
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:273
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
Definition: Phase.cpp:405
const vector_fp & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:509
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:200
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:751
virtual double density() const
Density (kg/m^3).
Definition: Phase.h:679
doublereal temperature() const
Temperature (K).
Definition: Phase.h:654
virtual void setTemperature(double temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:719
virtual void setMassFractions(const double *const y)
Set the mass fractions to the specified values and normalize them.
Definition: Phase.cpp:391
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:187
void getMassFractions(double *const y) const
Get the species mass fractions.
Definition: Phase.cpp:585
void setCriteria(doublereal ratio=10.0, doublereal slope=0.8, doublereal curve=0.8, doublereal prune=-0.1)
Set grid refinement criteria.
Definition: refine.cpp:24
void setGridMin(double gridmin)
Set the minimum allowable spacing between adjacent grid points [m].
Definition: refine.h:62
virtual void evalResidual(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the residual function.
Definition: StFlow.cpp:295
virtual void evalRightBoundary(double *x, double *res, int *diag, double rdt)
Evaluate all residual components at the right boundary.
Definition: StFlow.cpp:1061
virtual void evalContinuity(size_t j, double *x, double *r, int *diag, double rdt)
Evaluate the residual corresponding to the continuity equation at all interior grid points.
Definition: StFlow.cpp:1092
size_t m_kExcessLeft
Index of species with a large mass fraction at each boundary, for which the mass fraction may be calc...
Definition: StFlow.h:472
virtual AnyMap serialize(const double *soln) const
Save the state of this domain as an AnyMap.
Definition: StFlow.cpp:863
virtual void showSolution(const doublereal *x)
Print the solution.
Definition: StFlow.cpp:519
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.
Definition: StFlow.cpp:184
size_t rightExcessSpecies() const
Index of the species on the right boundary with the largest mass fraction.
Definition: StFlow.h:281
doublereal T_fixed(size_t j) const
The fixed temperature value at point j.
Definition: StFlow.h:130
StFlow(ThermoPhase *ph=0, size_t nsp=1, size_t points=1)
Create a new flow domain.
Definition: StFlow.cpp:20
void setBoundaryEmissivities(double e_left, double e_right)
Set the emissivities for the boundary values.
Definition: StFlow.cpp:1023
std::vector< size_t > m_kRadiating
Indices within the ThermoPhase of the radiating species.
Definition: StFlow.h:450
virtual void eval(size_t j, doublereal *x, doublereal *r, integer *mask, doublereal rdt)
Definition: StFlow.cpp:243
virtual XML_Node & save(XML_Node &o, const doublereal *const sol)
Save the current solution for this domain into an XML_Node.
Definition: StFlow.cpp:804
double m_tfixed
Temperature at the point used to fix the flame location.
Definition: StFlow.h:486
virtual bool componentActive(size_t n) const
Returns true if the specified component is an active part of the solver state.
Definition: StFlow.cpp:623
void setPressure(doublereal p)
Set the pressure.
Definition: StFlow.h:98
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
Definition: StFlow.cpp:637
void getWdot(doublereal *x, size_t j)
Write the net production rates at point j into array m_wdot
Definition: StFlow.h:291
virtual void _finalize(const doublereal *x)
In some cases, a domain may need to set parameters that depend on the initial solution estimate.
Definition: StFlow.cpp:196
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition: StFlow.cpp:578
void setGas(const doublereal *x, size_t j)
Set the gas object state to be consistent with the solution at point j.
Definition: StFlow.cpp:176
virtual void updateDiffFluxes(const doublereal *x, size_t j0, size_t j1)
Update the diffusive mass fluxes.
Definition: StFlow.cpp:536
void setFixedTempProfile(vector_fp &zfixed, vector_fp &tfixed)
Sometimes it is desired to carry out the simulation using a specified temperature profile,...
Definition: StFlow.h:115
virtual void resize(size_t components, size_t points)
Change the grid size. Called after grid refinement.
Definition: StFlow.cpp:107
size_t leftExcessSpecies() const
Index of the species on the left boundary with the largest mass fraction.
Definition: StFlow.h:276
double m_zfixed
Location of the point where temperature is fixed.
Definition: StFlow.h:483
void updateThermo(const doublereal *x, size_t j0, size_t j1)
Update the thermodynamic properties from point j0 to point j1 (inclusive), based on solution x.
Definition: StFlow.h:310
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.
Definition: StFlow.cpp:489
vector_fp m_qdotRadiation
radiative heat loss vector
Definition: StFlow.h:462
virtual void resetBadValues(double *xg)
Definition: StFlow.cpp:146
virtual void _getInitialSoln(double *x)
Write the initial solution estimate into array x.
Definition: StFlow.cpp:168
virtual void updateProperties(size_t jg, double *x, size_t jmin, size_t jmax)
Update the properties (thermo, transport, and diffusion flux).
Definition: StFlow.cpp:271
virtual std::string flowType() const
Return the type of flow domain being represented, either "Free Flame" or "Axisymmetric Stagnation".
Definition: StFlow.h:182
virtual void setupGrid(size_t n, const doublereal *z)
called to set up initial grid, and after grid refinement
Definition: StFlow.cpp:131
void setTransport(Transport &trans)
set the transport manager
Definition: StFlow.cpp:156
virtual size_t componentIndex(const std::string &name) const
index of component with name name.
Definition: StFlow.cpp:600
bool m_do_radiation
flag for the radiative heat loss
Definition: StFlow.h:459
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:102
virtual doublereal maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
Definition: ThermoPhase.h:218
virtual std::string type() const
String indicating the thermodynamic model implemented.
Definition: ThermoPhase.h:113
const AnyMap & input() const
Access input data associated with the phase description.
Base class for transport property managers.
virtual void getThermalDiffCoeffs(doublereal *const dt)
Return a vector of Thermal diffusion coefficients [kg/m/sec].
virtual doublereal viscosity()
virtual std::string transportType() const
Identifies the Transport object type.
virtual void getMultiDiffCoeffs(const size_t ld, doublereal *const d)
Return the Multicomponent diffusion coefficients. Units: [m^2/s].
virtual doublereal thermalConductivity()
Returns the mixture thermal conductivity in W/m/K.
virtual void getMixDiffCoeffs(doublereal *const d)
Returns a vector of mixture averaged diffusion coefficients.
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:103
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:529
void addAttribute(const std::string &attrib, const std::string &value)
Add or modify an attribute of the current node.
Definition: xml.cpp:467
std::vector< XML_Node * > getChildren(const std::string &name) const
Get a vector of pointers to XML_Node containing all of the children of the current node which match t...
Definition: xml.cpp:712
std::string value() const
Return the value of an XML node as a string.
Definition: xml.cpp:442
XML_Node & child(const size_t n) const
Return a changeable reference to the n'th child of the current node.
Definition: xml.cpp:547
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data.
Header for a file containing miscellaneous numerical functions.
This file contains definitions for utility functions and text for modules, inputfiles,...
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:164
Namespace for the Cantera kernel.
Definition: AnyMap.h:29
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:192
doublereal linearInterp(doublereal x, const vector_fp &xpts, const vector_fp &fpts)
Linearly interpolate a function defined on a discrete grid.
Definition: funcs.cpp:15
const double Undef
Fairly random number to be used to initialize variables against to see if they are subsequently defin...
Definition: ct_defs.h:161
doublereal getFloat(const XML_Node &parent, const std::string &name, const std::string &type="")
Get a floating-point value from a child element.
Definition: ctml.cpp:166
bool getOptionalFloat(const XML_Node &parent, const std::string &name, doublereal &fltRtn, const std::string &type="")
Get an optional floating-point value from a child element.
Definition: ctml.cpp:214
const double OneAtm
One atmosphere [Pa].
Definition: ct_defs.h:81
void addFloat(XML_Node &node, const std::string &titleString, const doublereal value, const std::string &unitsString="", const std::string &typeString="", const doublereal minval=Undef, const doublereal maxval=Undef)
This function adds a child node with the name, "float", with a value consisting of a single floating ...
Definition: ctml.cpp:21
void debuglog(const std::string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition: global.h:146
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:186
const double StefanBoltz
Stefan-Boltzmann constant [W/m2/K4].
Definition: ct_defs.h:125
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:184
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:113
size_t getFloatArray(const XML_Node &node, vector_fp &v, const bool convert=true, const std::string &unitsString="", const std::string &nodeName="floatArray")
This function reads the current node or a child node of the current node with the default name,...
Definition: ctml.cpp:258
void warn_user(const std::string &method, const std::string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition: global.h:245
void addFloatArray(XML_Node &node, const std::string &titleString, const size_t n, const doublereal *const values, const std::string &unitsString="", const std::string &typeString="", const doublereal minval=Undef, const doublereal maxval=Undef)
This function adds a child node with the name, "floatArray", with a value consisting of a comma separ...
Definition: ctml.cpp:42
void addNamedFloatArray(XML_Node &parentNode, const std::string &name, const size_t n, const doublereal *const vals, const std::string units="", const std::string type="", const doublereal minval=Undef, const doublereal maxval=Undef)
This function adds a child node with the name given by the first parameter with a value consisting of...
Definition: ctml.cpp:75