Cantera  3.1.0a1
Loading...
Searching...
No Matches
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
8#include "cantera/oneD/refine.h"
12#include "cantera/base/global.h"
13
14using namespace std;
15
16namespace Cantera
17{
18
19StFlow::StFlow(ThermoPhase* ph, size_t nsp, size_t points) :
20 Domain1D(nsp+c_offset_Y, points),
21 m_nsp(nsp)
22{
23 m_points = points;
24
25 if (ph == 0) {
26 return; // used to create a dummy object
27 }
28 m_thermo = ph;
29
30 size_t nsp2 = m_thermo->nSpecies();
31 if (nsp2 != m_nsp) {
32 m_nsp = nsp2;
34 }
35
36 // make a local copy of the species molecular weight vector
37 m_wt = m_thermo->molecularWeights();
38
39 // set pressure based on associated thermo object
40 setPressure(m_thermo->pressure());
41
42 // the species mass fractions are the last components in the solution
43 // vector, so the total number of components is the number of species
44 // plus the offset of the first mass fraction.
45 m_nv = c_offset_Y + m_nsp;
46
47 // enable all species equations by default
48 m_do_species.resize(m_nsp, true);
49
50 // but turn off the energy equation at all points
51 m_do_energy.resize(m_points,false);
52
53 m_diff.resize(m_nsp*m_points);
54 m_multidiff.resize(m_nsp*m_nsp*m_points);
55 m_flux.resize(m_nsp,m_points);
56 m_wdot.resize(m_nsp,m_points, 0.0);
58 m_dhk_dz.resize(m_nsp, m_points - 1, 0.0);
59 m_ybar.resize(m_nsp);
60 m_qdotRadiation.resize(m_points, 0.0);
61
62 //-------------- default solution bounds --------------------
63 setBounds(c_offset_U, -1e20, 1e20); // no bounds on u
64 setBounds(c_offset_V, -1e20, 1e20); // V
65 setBounds(c_offset_T, 200.0, 2*m_thermo->maxTemp()); // temperature bounds
66 setBounds(c_offset_L, -1e20, 1e20); // lambda should be negative
67 setBounds(c_offset_E, -1e20, 1e20); // no bounds for inactive component
68
69 // mass fraction bounds
70 for (size_t k = 0; k < m_nsp; k++) {
71 setBounds(c_offset_Y+k, -1.0e-7, 1.0e5);
72 }
73
74 //-------------------- grid refinement -------------------------
75 m_refiner->setActive(c_offset_U, false);
76 m_refiner->setActive(c_offset_V, false);
77 m_refiner->setActive(c_offset_T, false);
78 m_refiner->setActive(c_offset_L, false);
79
80 vector<double> gr;
81 for (size_t ng = 0; ng < m_points; ng++) {
82 gr.push_back(1.0*ng/m_points);
83 }
84 setupGrid(m_points, gr.data());
85
86 // Find indices for radiating species
87 m_kRadiating.resize(2, npos);
88 m_kRadiating[0] = m_thermo->speciesIndex("CO2");
89 m_kRadiating[1] = m_thermo->speciesIndex("H2O");
90}
91
92StFlow::StFlow(shared_ptr<ThermoPhase> th, size_t nsp, size_t points)
93 : StFlow(th.get(), nsp, points)
94{
96 m_solution->setThermo(th);
97}
98
99StFlow::StFlow(shared_ptr<Solution> sol, const string& id, size_t points)
100 : StFlow(sol->thermo().get(), sol->thermo()->nSpecies(), points)
101{
102 m_solution = sol;
103 m_id = id;
104 m_kin = m_solution->kinetics().get();
105 m_trans = m_solution->transport().get();
106 if (m_trans->transportModel() == "none") {
107 throw CanteraError("StFlow::StFlow",
108 "An appropriate transport model\nshould be set when instantiating the "
109 "Solution ('gas') object.");
110 }
111 m_solution->registerChangedCallback(this, [this]() {
112 setKinetics(m_solution->kinetics());
113 setTransport(m_solution->transport());
114 });
115}
116
117StFlow::~StFlow()
118{
119 if (m_solution) {
120 m_solution->removeChangedCallback(this);
121 }
122}
123
124string StFlow::domainType() const {
125 if (m_isFree) {
126 return "free-flow";
127 }
128 if (m_usesLambda) {
129 return "axisymmetric-flow";
130 }
131 return "unstrained-flow";
132}
133
134void StFlow::setKinetics(shared_ptr<Kinetics> kin)
135{
136 m_kin = kin.get();
137 m_solution->setKinetics(kin);
138}
139
140void StFlow::setTransport(shared_ptr<Transport> trans)
141{
142 if (!trans) {
143 throw CanteraError("StFlow::setTransport", "Unable to set empty transport.");
144 }
145 m_trans = trans.get();
146 if (m_trans->transportModel() == "none") {
147 throw CanteraError("StFlow::setTransport", "Invalid Transport model 'none'.");
148 }
149 m_do_multicomponent = (m_trans->transportModel() == "multicomponent" ||
150 m_trans->transportModel() == "multicomponent-CK");
151
152 m_diff.resize(m_nsp * m_points);
153 if (m_do_multicomponent) {
154 m_multidiff.resize(m_nsp * m_nsp*m_points);
155 m_dthermal.resize(m_nsp, m_points, 0.0);
156 }
157 m_solution->setTransport(trans);
158}
159
160void StFlow::resize(size_t ncomponents, size_t points)
161{
162 Domain1D::resize(ncomponents, points);
163 m_rho.resize(m_points, 0.0);
164 m_wtm.resize(m_points, 0.0);
165 m_cp.resize(m_points, 0.0);
166 m_visc.resize(m_points, 0.0);
167 m_tcon.resize(m_points, 0.0);
168
169 m_diff.resize(m_nsp*m_points);
170 if (m_do_multicomponent) {
171 m_multidiff.resize(m_nsp*m_nsp*m_points);
172 m_dthermal.resize(m_nsp, m_points, 0.0);
173 }
174 m_flux.resize(m_nsp,m_points);
175 m_wdot.resize(m_nsp,m_points, 0.0);
176 m_hk.resize(m_nsp, m_points, 0.0);
177 m_dhk_dz.resize(m_nsp, m_points - 1, 0.0);
178 m_do_energy.resize(m_points,false);
179 m_qdotRadiation.resize(m_points, 0.0);
180 m_fixedtemp.resize(m_points);
181
182 m_dz.resize(m_points-1);
183 m_z.resize(m_points);
184}
185
186void StFlow::setupGrid(size_t n, const double* z)
187{
188 resize(m_nv, n);
189
190 m_z[0] = z[0];
191 for (size_t j = 1; j < m_points; j++) {
192 if (z[j] <= z[j-1]) {
193 throw CanteraError("StFlow::setupGrid",
194 "grid points must be monotonically increasing");
195 }
196 m_z[j] = z[j];
197 m_dz[j-1] = m_z[j] - m_z[j-1];
198 }
199}
200
202{
203 double* x = xg + loc();
204 for (size_t j = 0; j < m_points; j++) {
205 double* Y = x + m_nv*j + c_offset_Y;
206 m_thermo->setMassFractions(Y);
207 m_thermo->getMassFractions(Y);
208 }
209}
210
211void StFlow::setTransportModel(const string& trans)
212{
213 m_solution->setTransportModel(trans);
214}
215
217 return m_trans->transportModel();
218}
219
221{
222 for (size_t j = 0; j < m_points; j++) {
223 T(x,j) = m_thermo->temperature();
224 m_thermo->getMassFractions(&Y(x, 0, j));
225 m_rho[j] = m_thermo->density();
226 }
227}
228
229void StFlow::setGas(const double* x, size_t j)
230{
231 m_thermo->setTemperature(T(x,j));
232 const double* yy = x + m_nv*j + c_offset_Y;
233 m_thermo->setMassFractions_NoNorm(yy);
234 m_thermo->setPressure(m_press);
235}
236
237void StFlow::setGasAtMidpoint(const double* x, size_t j)
238{
239 m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1)));
240 const double* yyj = x + m_nv*j + c_offset_Y;
241 const double* yyjp = x + m_nv*(j+1) + c_offset_Y;
242 for (size_t k = 0; k < m_nsp; k++) {
243 m_ybar[k] = 0.5*(yyj[k] + yyjp[k]);
244 }
245 m_thermo->setMassFractions_NoNorm(m_ybar.data());
246 m_thermo->setPressure(m_press);
247}
248
249void StFlow::_finalize(const double* x)
250{
251 if (!m_do_multicomponent && m_do_soret) {
252 throw CanteraError("StFlow::_finalize",
253 "Thermal diffusion (the Soret effect) is enabled, and requires "
254 "using a multicomponent transport model.");
255 }
256
257 size_t nz = m_zfix.size();
258 bool e = m_do_energy[0];
259 for (size_t j = 0; j < m_points; j++) {
260 if (e || nz == 0) {
261 m_fixedtemp[j] = T(x, j);
262 } else {
263 double zz = (z(j) - z(0))/(z(m_points - 1) - z(0));
264 double tt = linearInterp(zz, m_zfix, m_tfix);
265 m_fixedtemp[j] = tt;
266 }
267 }
268 if (e) {
269 solveEnergyEqn();
270 }
271
272 if (m_isFree) {
273 // If the domain contains the temperature fixed point, make sure that it
274 // is correctly set. This may be necessary when the grid has been modified
275 // externally.
276 if (m_tfixed != Undef) {
277 for (size_t j = 0; j < m_points; j++) {
278 if (z(j) == m_zfixed) {
279 return; // fixed point is already set correctly
280 }
281 }
282
283 for (size_t j = 0; j < m_points - 1; j++) {
284 // Find where the temperature profile crosses the current
285 // fixed temperature.
286 if ((T(x, j) - m_tfixed) * (T(x, j+1) - m_tfixed) <= 0.0) {
287 m_tfixed = T(x, j+1);
288 m_zfixed = z(j+1);
289 return;
290 }
291 }
292 }
293 }
294}
295
296void StFlow::eval(size_t jGlobal, double* xGlobal, double* rsdGlobal,
297 integer* diagGlobal, double rdt)
298{
299 // If evaluating a Jacobian, and the global point is outside the domain of
300 // influence for this domain, then skip evaluating the residual
301 if (jGlobal != npos && (jGlobal + 1 < firstPoint() || jGlobal > lastPoint() + 1)) {
302 return;
303 }
304
305 // start of local part of global arrays
306 double* x = xGlobal + loc();
307 double* rsd = rsdGlobal + loc();
308 integer* diag = diagGlobal + loc();
309
310 size_t jmin, jmax;
311 if (jGlobal == npos) { // evaluate all points
312 jmin = 0;
313 jmax = m_points - 1;
314 } else { // evaluate points for Jacobian
315 size_t jpt = (jGlobal == 0) ? 0 : jGlobal - firstPoint();
316 jmin = std::max<size_t>(jpt, 1) - 1;
317 jmax = std::min(jpt+1,m_points-1);
318 }
319
320 updateProperties(jGlobal, x, jmin, jmax);
321
322 if (m_do_radiation) { // Calculation of qdotRadiation
323 computeRadiation(x, jmin, jmax);
324 }
325
326 evalContinuity(x, rsd, diag, rdt, jmin, jmax);
327 evalMomentum(x, rsd, diag, rdt, jmin, jmax);
328 evalEnergy(x, rsd, diag, rdt, jmin, jmax);
329 evalLambda(x, rsd, diag, rdt, jmin, jmax);
330 evalElectricField(x, rsd, diag, rdt, jmin, jmax);
331 evalSpecies(x, rsd, diag, rdt, jmin, jmax);
332}
333
334void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
335{
336 // properties are computed for grid points from j0 to j1
337 size_t j0 = std::max<size_t>(jmin, 1) - 1;
338 size_t j1 = std::min(jmax+1,m_points-1);
339
340 updateThermo(x, j0, j1);
341 if (jg == npos || m_force_full_update) {
342 // update transport properties only if a Jacobian is not being
343 // evaluated, or if specifically requested
344 updateTransport(x, j0, j1);
345 }
346 if (jg == npos) {
347 double* Yleft = x + index(c_offset_Y, jmin);
348 m_kExcessLeft = distance(Yleft, max_element(Yleft, Yleft + m_nsp));
349 double* Yright = x + index(c_offset_Y, jmax);
350 m_kExcessRight = distance(Yright, max_element(Yright, Yright + m_nsp));
351 }
352
353 // update the species diffusive mass fluxes whether or not a
354 // Jacobian is being evaluated
355 updateDiffFluxes(x, j0, j1);
356}
357
358void StFlow::computeRadiation(double* x, size_t jmin, size_t jmax)
359{
360 // Variable definitions for the Planck absorption coefficient and the
361 // radiation calculation:
362 double k_P_ref = 1.0*OneAtm;
363
364 // Polynomial coefficients:
365 const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880,
366 0.51382, -1.86840e-5};
367 const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050,
368 56.310, -5.8169};
369
370 // Calculation of the two boundary values
371 double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4);
372 double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4);
373
374 for (size_t j = jmin; j < jmax; j++) {
375 // calculation of the mean Planck absorption coefficient
376 double k_P = 0;
377 // Absorption coefficient for H2O
378 if (m_kRadiating[1] != npos) {
379 double k_P_H2O = 0;
380 for (size_t n = 0; n <= 5; n++) {
381 k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n);
382 }
383 k_P_H2O /= k_P_ref;
384 k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O;
385 }
386 // Absorption coefficient for CO2
387 if (m_kRadiating[0] != npos) {
388 double k_P_CO2 = 0;
389 for (size_t n = 0; n <= 5; n++) {
390 k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n);
391 }
392 k_P_CO2 /= k_P_ref;
393 k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2;
394 }
395
396 // Calculation of the radiative heat loss term
397 double radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4)
398 - boundary_Rad_left - boundary_Rad_right);
399
400 // set the radiative heat loss vector
401 m_qdotRadiation[j] = radiative_heat_loss;
402 }
403}
404
405void StFlow::evalContinuity(double* x, double* rsd, int* diag,
406 double rdt, size_t jmin, size_t jmax)
407{
408 // The left boundary has the same form for all cases.
409 if (jmin == 0) { // left boundary
410 rsd[index(c_offset_U,jmin)] = -(rho_u(x,jmin + 1) - rho_u(x,jmin))/m_dz[jmin]
411 -(density(jmin + 1)*V(x,jmin + 1)
412 + density(jmin)*V(x,jmin));
413 diag[index(c_offset_U,jmin)] = 0; // Algebraic constraint
414 }
415
416 if (jmax == m_points - 1) { // right boundary
417 if (m_usesLambda) { // axisymmetric flow
418 rsd[index(c_offset_U, jmax)] = rho_u(x, jmax);
419 } else { // right boundary (same for unstrained/free-flow)
420 rsd[index(c_offset_U, jmax)] = rho_u(x, jmax) - rho_u(x, jmax - 1);
421 }
422 diag[index(c_offset_U, jmax)] = 0; // Algebraic constraint
423 }
424
425 // j0 and j1 are constrained to only interior points
426 size_t j0 = std::max<size_t>(jmin, 1);
427 size_t j1 = std::min(jmax, m_points - 2);
428 if (m_usesLambda) { // "axisymmetric-flow"
429 for (size_t j = j0; j <= j1; j++) { // interior points
430 // For "axisymmetric-flow", the continuity equation propagates the
431 // mass flow rate information to the left (j+1 -> j) from the value
432 // specified at the right boundary. The lambda information propagates
433 // in the opposite direction.
434 rsd[index(c_offset_U,j)] = -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
435 -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
436 diag[index(c_offset_U, j)] = 0; // Algebraic constraint
437 }
438 } else if (m_isFree) { // "free-flow"
439 for (size_t j = j0; j <= j1; j++) {
440 // terms involving V are zero as V=0 by definition
441 if (grid(j) > m_zfixed) {
442 rsd[index(c_offset_U,j)] = -(rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1];
443 } else if (grid(j) == m_zfixed) {
444 if (m_do_energy[j]) {
445 rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
446 } else {
447 rsd[index(c_offset_U,j)] = (rho_u(x,j) - m_rho[0]*0.3); // why 0.3?
448 }
449 } else { // grid(j < m_zfixed
450 rsd[index(c_offset_U,j)] = -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j];
451 }
452 diag[index(c_offset_U, j)] = 0; // Algebraic constraint
453 }
454 } else { // "unstrained-flow" (fixed mass flow rate)
455 for (size_t j = j0; j <= j1; j++) {
456 rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
457 diag[index(c_offset_U, j)] = 0; // Algebraic constraint
458 }
459 }
460}
461
462void StFlow::evalMomentum(double* x, double* rsd, int* diag,
463 double rdt, size_t jmin, size_t jmax)
464{
465 if (!m_usesLambda) { //disable this equation
466 for (size_t j = jmin; j <= jmax; j++) {
467 rsd[index(c_offset_V, j)] = V(x, j);
468 diag[index(c_offset_V, j)] = 0;
469 }
470 return;
471 }
472
473 if (jmin == 0) { // left boundary
474 rsd[index(c_offset_V,jmin)] = V(x,jmin);
475 }
476
477 if (jmax == m_points - 1) { // right boundary
478 rsd[index(c_offset_V,jmax)] = V(x,jmax);
479 diag[index(c_offset_V,jmax)] = 0;
480 }
481
482 // j0 and j1 are constrained to only interior points
483 size_t j0 = std::max<size_t>(jmin, 1);
484 size_t j1 = std::min(jmax, m_points - 2);
485 for (size_t j = j0; j <= j1; j++) { // interior points
486 rsd[index(c_offset_V,j)] = (shear(x, j) - lambda(x, j)
487 - rho_u(x, j) * dVdz(x, j)
488 - m_rho[j] * V(x, j) * V(x, j)) / m_rho[j]
489 - rdt * (V(x, j) - V_prev(j));
490 diag[index(c_offset_V, j)] = 1;
491 }
492}
493
494void StFlow::evalLambda(double* x, double* rsd, int* diag,
495 double rdt, size_t jmin, size_t jmax)
496{
497 if (!m_usesLambda) { // disable this equation
498 for (size_t j = jmin; j <= jmax; j++) {
499 rsd[index(c_offset_L, j)] = lambda(x, j);
500 diag[index(c_offset_L, j)] = 0;
501 }
502 return;
503 }
504
505 if (jmin == 0) { // left boundary
506 rsd[index(c_offset_L, jmin)] = -rho_u(x, jmin);
507 }
508
509 if (jmax == m_points - 1) { // right boundary
510 rsd[index(c_offset_L, jmax)] = lambda(x, jmax) - lambda(x, jmax-1);
511 diag[index(c_offset_L, jmax)] = 0;
512 }
513
514 // j0 and j1 are constrained to only interior points
515 size_t j0 = std::max<size_t>(jmin, 1);
516 size_t j1 = std::min(jmax, m_points - 2);
517 for (size_t j = j0; j <= j1; j++) { // interior points
518 rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j - 1);
519 }
520}
521
522void StFlow::evalEnergy(double* x, double* rsd, int* diag,
523 double rdt, size_t jmin, size_t jmax)
524{
525 if (jmin == 0) { // left boundary
526 rsd[index(c_offset_T,jmin)] = T(x,jmin);
527 }
528
529 if (jmax == m_points - 1) { // right boundary
530 rsd[index(c_offset_T, jmax)] = T(x,jmax);
531 }
532
533 // j0 and j1 are constrained to only interior points
534 size_t j0 = std::max<size_t>(jmin, 1);
535 size_t j1 = std::min(jmax, m_points - 2);
536 for (size_t j = j0; j <= j1; j++) {
537 if (m_do_energy[j]) {
538 grad_hk(x, j);
539 double sum = 0.0;
540 for (size_t k = 0; k < m_nsp; k++) {
541 double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
542 sum += wdot(k,j)*m_hk(k,j);
543 sum += flxk * m_dhk_dz(k,j) / m_wt[k];
544 }
545
546 rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dTdz(x,j)
547 - divHeatFlux(x,j) - sum;
548 rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
549 rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
550 rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j]));
551 diag[index(c_offset_T, j)] = 1;
552 } else {
553 // residual equations if the energy equation is disabled
554 rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
555 diag[index(c_offset_T, j)] = 0;
556 }
557 }
558}
559
560void StFlow::evalSpecies(double* x, double* rsd, int* diag,
561 double rdt, size_t jmin, size_t jmax)
562{
563 if (jmin == 0) { // left boundary
564 double sum = 0.0;
565 for (size_t k = 0; k < m_nsp; k++) {
566 sum += Y(x,k,jmin);
567 rsd[index(c_offset_Y + k, jmin)] = -(m_flux(k,jmin) +
568 rho_u(x,jmin) * Y(x,k,jmin));
569 }
570 rsd[index(c_offset_Y + leftExcessSpecies(), jmin)] = 1.0 - sum;
571 }
572
573 if (jmax == m_points - 1) { // right boundary
574 double sum = 0.0;
575 for (size_t k = 0; k < m_nsp; k++) {
576 sum += Y(x,k,jmax);
577 rsd[index(k+c_offset_Y,jmax)] = m_flux(k,jmax-1) +
578 rho_u(x,jmax)*Y(x,k,jmax);
579 }
580 rsd[index(c_offset_Y + rightExcessSpecies(), jmax)] = 1.0 - sum;
581 diag[index(c_offset_Y + rightExcessSpecies(), jmax)] = 0;
582 }
583
584 // j0 and j1 are constrained to only interior points
585 size_t j0 = std::max<size_t>(jmin, 1);
586 size_t j1 = std::min(jmax, m_points - 2);
587 for (size_t j = j0; j <= j1; j++) {
588 for (size_t k = 0; k < m_nsp; k++) {
589 double convec = rho_u(x,j)*dYdz(x,k,j);
590 double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) / (z(j+1) - z(j-1));
591 rsd[index(c_offset_Y + k, j)] = (m_wt[k]*(wdot(k,j))
592 - convec - diffus)/m_rho[j]
593 - rdt*(Y(x,k,j) - Y_prev(k,j));
594 diag[index(c_offset_Y + k, j)] = 1;
595 }
596 }
597}
598
599void StFlow::evalElectricField(double* x, double* rsd, int* diag,
600 double rdt, size_t jmin, size_t jmax)
601{
602 for (size_t j = jmin; j <= jmax; j++) {
603 // The same value is used for left/right/interior points
604 rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
605 }
606}
607
608void StFlow::updateTransport(double* x, size_t j0, size_t j1)
609{
610 if (m_do_multicomponent) {
611 for (size_t j = j0; j < j1; j++) {
612 setGasAtMidpoint(x,j);
613 double wtm = m_thermo->meanMolecularWeight();
614 double rho = m_thermo->density();
615 m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
616 m_trans->getMultiDiffCoeffs(m_nsp, &m_multidiff[mindex(0,0,j)]);
617
618 // Use m_diff as storage for the factor outside the summation
619 for (size_t k = 0; k < m_nsp; k++) {
620 m_diff[k+j*m_nsp] = m_wt[k] * rho / (wtm*wtm);
621 }
622
623 m_tcon[j] = m_trans->thermalConductivity();
624 if (m_do_soret) {
625 m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + j*m_nsp);
626 }
627 }
628 } else { // mixture averaged transport
629 for (size_t j = j0; j < j1; j++) {
630 setGasAtMidpoint(x,j);
631 m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
632 m_trans->getMixDiffCoeffs(&m_diff[j*m_nsp]);
633 double rho = m_thermo->density();
634 double wtm = m_thermo->meanMolecularWeight();
635 for (size_t k=0; k < m_nsp; k++) {
636 m_diff[k+j*m_nsp] *= m_wt[k] * rho / wtm;
637 }
638 m_tcon[j] = m_trans->thermalConductivity();
639 }
640 }
641}
642
643void StFlow::show(const double* x)
644{
645 writelog(" Pressure: {:10.4g} Pa\n", m_press);
646
648
649 if (m_do_radiation) {
650 writeline('-', 79, false, true);
651 writelog("\n z radiative heat loss");
652 writeline('-', 79, false, true);
653 for (size_t j = 0; j < m_points; j++) {
654 writelog("\n {:10.4g} {:10.4g}", m_z[j], m_qdotRadiation[j]);
655 }
656 writelog("\n");
657 }
658}
659
660void StFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
661{
662 if (m_do_multicomponent) {
663 for (size_t j = j0; j < j1; j++) {
664 double dz = z(j+1) - z(j);
665 for (size_t k = 0; k < m_nsp; k++) {
666 double sum = 0.0;
667 for (size_t m = 0; m < m_nsp; m++) {
668 sum += m_wt[m] * m_multidiff[mindex(k,m,j)] * (X(x,m,j+1)-X(x,m,j));
669 }
670 m_flux(k,j) = sum * m_diff[k+j*m_nsp] / dz;
671 }
672 }
673 } else {
674 for (size_t j = j0; j < j1; j++) {
675 double sum = 0.0;
676 double dz = z(j+1) - z(j);
677 for (size_t k = 0; k < m_nsp; k++) {
678 m_flux(k,j) = m_diff[k+m_nsp*j];
679 m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
680 sum -= m_flux(k,j);
681 }
682 // correction flux to insure that \sum_k Y_k V_k = 0.
683 for (size_t k = 0; k < m_nsp; k++) {
684 m_flux(k,j) += sum*Y(x,k,j);
685 }
686 }
687 }
688
689 if (m_do_soret) {
690 for (size_t m = j0; m < j1; m++) {
691 double gradlogT = 2.0 * (T(x,m+1) - T(x,m)) /
692 ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m)));
693 for (size_t k = 0; k < m_nsp; k++) {
694 m_flux(k,m) -= m_dthermal(k,m)*gradlogT;
695 }
696 }
697 }
698}
699
700string StFlow::componentName(size_t n) const
701{
702 switch (n) {
703 case c_offset_U:
704 return "velocity";
705 case c_offset_V:
706 return "spread_rate";
707 case c_offset_T:
708 return "T";
709 case c_offset_L:
710 return "lambda";
711 case c_offset_E:
712 return "eField";
713 default:
714 if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) {
715 return m_thermo->speciesName(n - c_offset_Y);
716 } else {
717 return "<unknown>";
718 }
719 }
720}
721
722size_t StFlow::componentIndex(const string& name) const
723{
724 if (name=="velocity") {
725 return c_offset_U;
726 } else if (name=="spread_rate") {
727 return c_offset_V;
728 } else if (name=="T") {
729 return c_offset_T;
730 } else if (name=="lambda") {
731 return c_offset_L;
732 } else if (name == "eField") {
733 return c_offset_E;
734 } else {
735 for (size_t n=c_offset_Y; n<m_nsp+c_offset_Y; n++) {
736 if (componentName(n)==name) {
737 return n;
738 }
739 }
740 throw CanteraError("StFlow1D::componentIndex",
741 "no component named " + name);
742 }
743}
744
745bool StFlow::componentActive(size_t n) const
746{
747 switch (n) {
748 case c_offset_V: // spread_rate
749 return m_usesLambda;
750 case c_offset_L: // lambda
751 return m_usesLambda;
752 case c_offset_E: // eField
753 return false;
754 default:
755 return true;
756 }
757}
758
760{
761 AnyMap state = Domain1D::getMeta();
762 state["transport-model"] = m_trans->transportModel();
763
764 state["phase"]["name"] = m_thermo->name();
765 AnyValue source = m_thermo->input().getMetadata("filename");
766 state["phase"]["source"] = source.empty() ? "<unknown>" : source.asString();
767
768 state["radiation-enabled"] = m_do_radiation;
769 if (m_do_radiation) {
770 state["emissivity-left"] = m_epsilon_left;
771 state["emissivity-right"] = m_epsilon_right;
772 }
773
774 set<bool> energy_flags(m_do_energy.begin(), m_do_energy.end());
775 if (energy_flags.size() == 1) {
776 state["energy-enabled"] = m_do_energy[0];
777 } else {
778 state["energy-enabled"] = m_do_energy;
779 }
780
781 state["Soret-enabled"] = m_do_soret;
782
783 set<bool> species_flags(m_do_species.begin(), m_do_species.end());
784 if (species_flags.size() == 1) {
785 state["species-enabled"] = m_do_species[0];
786 } else {
787 for (size_t k = 0; k < m_nsp; k++) {
788 state["species-enabled"][m_thermo->speciesName(k)] = m_do_species[k];
789 }
790 }
791
792 state["refine-criteria"]["ratio"] = m_refiner->maxRatio();
793 state["refine-criteria"]["slope"] = m_refiner->maxDelta();
794 state["refine-criteria"]["curve"] = m_refiner->maxSlope();
795 state["refine-criteria"]["prune"] = m_refiner->prune();
796 state["refine-criteria"]["grid-min"] = m_refiner->gridMin();
797 state["refine-criteria"]["max-points"] =
798 static_cast<long int>(m_refiner->maxPoints());
799
800 if (m_zfixed != Undef) {
801 state["fixed-point"]["location"] = m_zfixed;
802 state["fixed-point"]["temperature"] = m_tfixed;
803 }
804
805 return state;
806}
807
808shared_ptr<SolutionArray> StFlow::asArray(const double* soln) const
809{
810 auto arr = SolutionArray::create(
811 m_solution, static_cast<int>(nPoints()), getMeta());
812 arr->addExtra("grid", false); // leading entry
813 AnyValue value;
814 value = m_z;
815 arr->setComponent("grid", value);
816 vector<double> data(nPoints());
817 for (size_t i = 0; i < nComponents(); i++) {
818 if (componentActive(i)) {
819 auto name = componentName(i);
820 for (size_t j = 0; j < nPoints(); j++) {
821 data[j] = soln[index(i, j)];
822 }
823 if (!arr->hasComponent(name)) {
824 arr->addExtra(name, componentIndex(name) > c_offset_Y);
825 }
826 value = data;
827 arr->setComponent(name, value);
828 }
829 }
830 value = m_rho;
831 arr->setComponent("D", value); // use density rather than pressure
832
833 if (m_do_radiation) {
834 arr->addExtra("radiative-heat-loss", true); // add at end
835 value = m_qdotRadiation;
836 arr->setComponent("radiative-heat-loss", value);
837 }
838
839 return arr;
840}
841
842void StFlow::fromArray(SolutionArray& arr, double* soln)
843{
844 Domain1D::setMeta(arr.meta());
845 arr.setLoc(0);
846 auto phase = arr.thermo();
847 m_press = phase->pressure();
848
849 const auto grid = arr.getComponent("grid").as<vector<double>>();
850 setupGrid(nPoints(), &grid[0]);
851
852 for (size_t i = 0; i < nComponents(); i++) {
853 if (!componentActive(i)) {
854 continue;
855 }
856 string name = componentName(i);
857 if (arr.hasComponent(name)) {
858 const vector<double> data = arr.getComponent(name).as<vector<double>>();
859 for (size_t j = 0; j < nPoints(); j++) {
860 soln[index(i,j)] = data[j];
861 }
862 } else {
863 warn_user("StFlow::fromArray", "Saved state does not contain values for "
864 "component '{}' in domain '{}'.", name, id());
865 }
866 }
867
868 updateProperties(npos, soln + loc(), 0, m_points - 1);
869 setMeta(arr.meta());
870}
871
872void StFlow::setMeta(const AnyMap& state)
873{
874 if (state.hasKey("energy-enabled")) {
875 const AnyValue& ee = state["energy-enabled"];
876 if (ee.isScalar()) {
877 m_do_energy.assign(nPoints(), ee.asBool());
878 } else {
879 m_do_energy = ee.asVector<bool>(nPoints());
880 }
881 }
882
883 setTransportModel(state.getString("transport-model", "mixture-averaged"));
884
885 if (state.hasKey("Soret-enabled")) {
886 m_do_soret = state["Soret-enabled"].asBool();
887 }
888
889 if (state.hasKey("species-enabled")) {
890 const AnyValue& se = state["species-enabled"];
891 if (se.isScalar()) {
892 m_do_species.assign(m_thermo->nSpecies(), se.asBool());
893 } else {
894 m_do_species = se.asVector<bool>(m_thermo->nSpecies());
895 }
896 }
897
898 if (state.hasKey("radiation-enabled")) {
899 m_do_radiation = state["radiation-enabled"].asBool();
900 if (m_do_radiation) {
901 m_epsilon_left = state["emissivity-left"].asDouble();
902 m_epsilon_right = state["emissivity-right"].asDouble();
903 }
904 }
905
906 if (state.hasKey("refine-criteria")) {
907 const AnyMap& criteria = state["refine-criteria"].as<AnyMap>();
908 double ratio = criteria.getDouble("ratio", m_refiner->maxRatio());
909 double slope = criteria.getDouble("slope", m_refiner->maxDelta());
910 double curve = criteria.getDouble("curve", m_refiner->maxSlope());
911 double prune = criteria.getDouble("prune", m_refiner->prune());
912 m_refiner->setCriteria(ratio, slope, curve, prune);
913
914 if (criteria.hasKey("grid-min")) {
915 m_refiner->setGridMin(criteria["grid-min"].asDouble());
916 }
917 if (criteria.hasKey("max-points")) {
918 m_refiner->setMaxPoints(criteria["max-points"].asInt());
919 }
920 }
921
922 if (state.hasKey("fixed-point")) {
923 m_zfixed = state["fixed-point"]["location"].asDouble();
924 m_tfixed = state["fixed-point"]["temperature"].asDouble();
925 }
926}
927
928void StFlow::solveEnergyEqn(size_t j)
929{
930 bool changed = false;
931 if (j == npos) {
932 for (size_t i = 0; i < m_points; i++) {
933 if (!m_do_energy[i]) {
934 changed = true;
935 }
936 m_do_energy[i] = true;
937 }
938 } else {
939 if (!m_do_energy[j]) {
940 changed = true;
941 }
942 m_do_energy[j] = true;
943 }
944 m_refiner->setActive(c_offset_U, true);
945 m_refiner->setActive(c_offset_V, true);
946 m_refiner->setActive(c_offset_T, true);
947 if (changed) {
949 }
950}
951
953{
954 throw NotImplementedError("StFlow::getSolvingStage",
955 "Not used by '{}' objects.", type());
956}
957
958void StFlow::setSolvingStage(const size_t stage)
959{
960 throw NotImplementedError("StFlow::setSolvingStage",
961 "Not used by '{}' objects.", type());
962}
963
965{
966 throw NotImplementedError("StFlow::solveElectricField",
967 "Not used by '{}' objects.", type());
968}
969
971{
972 throw NotImplementedError("StFlow::fixElectricField",
973 "Not used by '{}' objects.", type());
974}
975
976bool StFlow::doElectricField(size_t j) const
977{
978 throw NotImplementedError("StFlow::doElectricField",
979 "Not used by '{}' objects.", type());
980}
981
982void StFlow::setBoundaryEmissivities(double e_left, double e_right)
983{
984 if (e_left < 0 || e_left > 1) {
985 throw CanteraError("StFlow::setBoundaryEmissivities",
986 "The left boundary emissivity must be between 0.0 and 1.0!");
987 } else if (e_right < 0 || e_right > 1) {
988 throw CanteraError("StFlow::setBoundaryEmissivities",
989 "The right boundary emissivity must be between 0.0 and 1.0!");
990 } else {
991 m_epsilon_left = e_left;
992 m_epsilon_right = e_right;
993 }
994}
995
996void StFlow::fixTemperature(size_t j)
997{
998 bool changed = false;
999 if (j == npos) {
1000 for (size_t i = 0; i < m_points; i++) {
1001 if (m_do_energy[i]) {
1002 changed = true;
1003 }
1004 m_do_energy[i] = false;
1005 }
1006 } else {
1007 if (m_do_energy[j]) {
1008 changed = true;
1009 }
1010 m_do_energy[j] = false;
1011 }
1012 m_refiner->setActive(c_offset_U, false);
1013 m_refiner->setActive(c_offset_V, false);
1014 m_refiner->setActive(c_offset_T, false);
1015 if (changed) {
1016 needJacUpdate();
1017 }
1018}
1019
1020void StFlow::grad_hk(const double* x, size_t j)
1021{
1022 for(size_t k = 0; k < m_nsp; k++) {
1023 if (u(x, j) > 0.0) {
1024 m_dhk_dz(k,j) = (m_hk(k,j) - m_hk(k,j-1)) / m_dz[j - 1];
1025 }
1026 else {
1027 m_dhk_dz(k,j) = (m_hk(k,j+1) - m_hk(k,j)) / m_dz[j];
1028 }
1029 }
1030}
1031
1032} // namespace
Header file defining class TransportFactory (see TransportFactory)
Headers for the Transport object, which is the virtual base class for all transport property evaluato...
const AnyValue & getMetadata(const string &key) const
Get a value from the metadata applicable to the AnyMap tree containing this node.
Definition AnyMap.cpp:580
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
double getDouble(const string &key, double default_) const
If key exists, return it as a double, otherwise return default_.
Definition AnyMap.cpp:1520
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
const string & getString(const string &key, const string &default_) const
If key exists, return it as a string, otherwise return default_.
Definition AnyMap.cpp:1530
A wrapper for a variable whose type is determined at runtime.
Definition AnyMap.h:86
const string & asString() const
Return the held value, if it is a string.
Definition AnyMap.cpp:739
bool & asBool()
Return the held value, if it is a bool.
Definition AnyMap.cpp:871
bool empty() const
Return boolean indicating whether AnyValue is empty.
Definition AnyMap.cpp:647
bool isScalar() const
Returns true if the held value is a scalar type (such as double, long int, string,...
Definition AnyMap.cpp:651
const 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:109
const T & as() const
Get the value of this key as the specified type.
Definition AnyMap.inl.h:16
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition Array.h:203
virtual 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:47
Base class for exceptions thrown by Cantera classes.
Base class for one-dimensional domains.
Definition Domain1D.h:28
size_t lastPoint() const
The index of the last (that is, right-most) grid point belonging to this domain.
Definition Domain1D.h:400
shared_ptr< Solution > m_solution
Composite thermo/kinetics/transport handler.
Definition Domain1D.h:567
size_t nComponents() const
Number of components at each grid point.
Definition Domain1D.h:145
virtual void setMeta(const AnyMap &meta)
Retrieve meta data.
Definition Domain1D.cpp:155
size_t nPoints() const
Number of grid points in this domain.
Definition Domain1D.h:167
virtual void resize(size_t nv, size_t np)
Resize the domain to have nv components and np grid points.
Definition Domain1D.cpp:26
size_t m_points
Number of grid points.
Definition Domain1D.h:537
string m_id
Identity tag for the domain.
Definition Domain1D.h:560
string type() const
String indicating the domain implemented.
Definition Domain1D.h:49
size_t firstPoint() const
The index of the first (that is, left-most) grid point belonging to this domain.
Definition Domain1D.h:392
void needJacUpdate()
Set this if something has changed in the governing equations (for example, the value of a constant ha...
Definition Domain1D.cpp:95
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:384
virtual AnyMap getMeta() const
Retrieve meta data.
Definition Domain1D.cpp:103
virtual void show(std::ostream &s, const double *x)
Print the solution.
Definition Domain1D.h:459
An error indicating that an unimplemented function has been called.
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:231
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
Definition Phase.cpp:355
double temperature() const
Temperature (K).
Definition Phase.h:562
virtual void setPressure(double p)
Set the internally stored pressure (Pa) at constant temperature and composition.
Definition Phase.h:616
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:655
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:395
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:129
virtual double density() const
Density (kg/m^3).
Definition Phase.h:587
virtual void setTemperature(double temp)
Set the internally stored temperature of the phase (K).
Definition Phase.h:623
virtual void setMassFractions(const double *const y)
Set the mass fractions to the specified values and normalize them.
Definition Phase.cpp:341
void getMassFractions(double *const y) const
Get the species mass fractions.
Definition Phase.cpp:471
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:580
string name() const
Return the name of the phase.
Definition Phase.cpp:20
A container class holding arrays of state information.
void setLoc(int loc, bool restore=true)
Update the buffered location used to access SolutionArray entries.
AnyValue getComponent(const string &name) const
Retrieve a component of the SolutionArray by name.
bool hasComponent(const string &name) const
Check whether SolutionArray contains a component.
AnyMap & meta()
SolutionArray meta data.
shared_ptr< ThermoPhase > thermo()
Retrieve associated ThermoPhase object.
static shared_ptr< SolutionArray > create(const shared_ptr< Solution > &sol, int size=0, const AnyMap &meta={})
Instantiate a new SolutionArray reference.
static shared_ptr< Solution > create()
Create an empty Solution object.
Definition Solution.h:54
This class represents 1D flow domains that satisfy the one-dimensional similarity solution for chemic...
Definition StFlow.h:45
void eval(size_t jGlobal, double *xGlobal, double *rsdGlobal, integer *diagGlobal, double rdt) override
Evaluate the residual functions for axisymmetric stagnation flow.
Definition StFlow.cpp:296
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:659
void setMeta(const AnyMap &state) override
Retrieve meta data.
Definition StFlow.cpp:872
void setTransportModel(const string &trans)
Set the transport model.
Definition StFlow.cpp:211
void setTransport(shared_ptr< Transport > trans) override
Set transport model to existing instance.
Definition StFlow.cpp:140
void setKinetics(shared_ptr< Kinetics > kin) override
Set the kinetics manager.
Definition StFlow.cpp:134
void resetBadValues(double *xg) override
When called, this function should reset "bad" values in the state vector such as negative species con...
Definition StFlow.cpp:201
size_t rightExcessSpecies() const
Index of the species on the right boundary with the largest mass fraction.
Definition StFlow.h:311
vector< double > m_qdotRadiation
radiative heat loss vector
Definition StFlow.h:649
virtual void evalMomentum(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the momentum equation residual.
Definition StFlow.cpp:462
void updateThermo(const double *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:483
void resize(size_t components, size_t points) override
Change the grid size. Called after grid refinement.
Definition StFlow.cpp:160
StFlow(ThermoPhase *ph=0, size_t nsp=1, size_t points=1)
Create a new flow domain.
Definition StFlow.cpp:19
virtual void evalContinuity(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the continuity equation residual.
Definition StFlow.cpp:405
void setBoundaryEmissivities(double e_left, double e_right)
Set the emissivities for the boundary values.
Definition StFlow.cpp:982
virtual void evalEnergy(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the energy equation residual.
Definition StFlow.cpp:522
vector< double > m_diff
Array of size m_nsp by m_points for saving density times diffusion coefficient times species molar ma...
Definition StFlow.h:611
shared_ptr< SolutionArray > asArray(const double *soln) const override
Save the state of this domain as a SolutionArray.
Definition StFlow.cpp:808
size_t componentIndex(const string &name) const override
index of component with name name.
Definition StFlow.cpp:722
void setGas(const double *x, size_t j)
Set the gas object state to be consistent with the solution at point j.
Definition StFlow.cpp:229
double m_tfixed
Temperature at the point used to fix the flame location.
Definition StFlow.h:675
virtual bool componentActive(size_t n) const
Returns true if the specified component is an active part of the solver state.
Definition StFlow.cpp:745
Array2D m_hk
Array of size m_nsp by m_points for saving molar enthalpies.
Definition StFlow.h:617
virtual bool doElectricField(size_t j) const
Retrieve flag indicating whether electric field is solved or not (used by IonFlow specialization)
Definition StFlow.cpp:976
virtual void evalSpecies(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the species equations' residuals.
Definition StFlow.cpp:560
void setupGrid(size_t n, const double *z) override
called to set up initial grid, and after grid refinement
Definition StFlow.cpp:186
size_t leftExcessSpecies() const
Index of the species on the left boundary with the largest mass fraction.
Definition StFlow.h:306
Array2D m_dhk_dz
Array of size m_nsp by m_points-1 for saving enthalpy fluxes.
Definition StFlow.h:620
virtual void evalElectricField(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the electric field equation residual to be zero everywhere.
Definition StFlow.cpp:599
double m_zfixed
Location of the point where temperature is fixed.
Definition StFlow.h:672
void _finalize(const double *x) override
In some cases, a domain may need to set parameters that depend on the initial solution estimate.
Definition StFlow.cpp:249
virtual size_t getSolvingStage() const
Get the solving stage (used by IonFlow specialization)
Definition StFlow.cpp:952
size_t m_nsp
Number of species in the mechanism.
Definition StFlow.h:625
virtual void evalLambda(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the lambda equation residual.
Definition StFlow.cpp:494
void fromArray(SolutionArray &arr, double *soln) override
Restore the solution for this domain from a SolutionArray.
Definition StFlow.cpp:842
AnyMap getMeta() const override
Retrieve meta data.
Definition StFlow.cpp:759
virtual void updateDiffFluxes(const double *x, size_t j0, size_t j1)
Update the diffusive mass fluxes.
Definition StFlow.cpp:660
string componentName(size_t n) const override
Name of the nth component. May be overloaded.
Definition StFlow.cpp:700
void setGasAtMidpoint(const double *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:237
virtual void grad_hk(const double *x, size_t j)
Get the gradient of species specific molar enthalpies.
Definition StFlow.cpp:1020
string transportModel() const
Retrieve transport model.
Definition StFlow.cpp:216
void computeRadiation(double *x, size_t jmin, size_t jmax)
Computes the radiative heat loss vector over points jmin to jmax and stores the data in the qdotRadia...
Definition StFlow.cpp:358
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:334
string domainType() const override
Domain type flag.
Definition StFlow.cpp:124
void show(const double *x) override
Print the solution.
Definition StFlow.cpp:643
virtual void setSolvingStage(const size_t stage)
Solving stage mode for handling ionized species (used by IonFlow specialization)
Definition StFlow.cpp:958
void setPressure(double p)
Set the pressure.
Definition StFlow.h:111
virtual void fixElectricField(size_t j=npos)
Set to fix voltage in a point (used by IonFlow specialization)
Definition StFlow.cpp:970
virtual void updateTransport(double *x, size_t j0, size_t j1)
Update the transport properties at grid points in the range from j0 to j1, based on solution x.
Definition StFlow.cpp:608
virtual void solveElectricField(size_t j=npos)
Set to solve electric field in a point (used by IonFlow specialization)
Definition StFlow.cpp:964
void _getInitialSoln(double *x) override
Write the initial solution estimate into array x.
Definition StFlow.cpp:220
vector< size_t > m_kRadiating
Indices within the ThermoPhase of the radiating species.
Definition StFlow.h:637
double T_fixed(size_t j) const
The fixed temperature value at point j.
Definition StFlow.h:143
bool m_do_radiation
flag for the radiative heat loss
Definition StFlow.h:646
Base class for a phase with thermodynamic properties.
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
const AnyMap & input() const
Access input data associated with the phase description.
virtual void getThermalDiffCoeffs(double *const dt)
Return a vector of Thermal diffusion coefficients [kg/m/sec].
Definition Transport.h:271
virtual string transportModel() const
Identifies the model represented by this Transport object.
Definition Transport.h:93
virtual void getMixDiffCoeffs(double *const d)
Returns a vector of mixture averaged diffusion coefficients.
Definition Transport.h:315
virtual double thermalConductivity()
Returns the mixture thermal conductivity in W/m/K.
Definition Transport.h:155
virtual double viscosity()
The viscosity in Pa-s.
Definition Transport.h:122
virtual void getMultiDiffCoeffs(const size_t ld, double *const d)
Return the Multicomponent diffusion coefficients. Units: [m^2/s].
Definition Transport.h:300
Header for a file containing miscellaneous numerical functions.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
double linearInterp(double x, const vector< double > &xpts, const vector< double > &fpts)
Linearly interpolate a function defined on a discrete grid.
Definition funcs.cpp:13
const double OneAtm
One atmosphere [Pa].
Definition ct_defs.h:96
const double StefanBoltz
Stefan-Boltzmann constant [W/m2/K4].
Definition ct_defs.h:128
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition global.h:267
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
const double Undef
Fairly random number to be used to initialize variables against to see if they are subsequently defin...
Definition ct_defs.h:164
@ c_offset_U
axial velocity
Definition StFlow.h:25
@ c_offset_L
(1/r)dP/dr
Definition StFlow.h:28
@ c_offset_V
strain rate
Definition StFlow.h:26
@ c_offset_E
electric field equation
Definition StFlow.h:29
@ c_offset_Y
mass fractions
Definition StFlow.h:30
@ c_offset_T
temperature
Definition StFlow.h:27