Cantera  2.1.2
StFlow.cpp
Go to the documentation of this file.
1 /**
2  * @file StFlow.cpp
3  */
4 // Copyright 2002 California Institute of Technology
5
6 #include <stdlib.h>
7 #include <time.h>
8
9 #include "cantera/oneD/StFlow.h"
10 #include "cantera/base/ctml.h"
11 #include "cantera/oneD/MultiJac.h"
12
13 #include <cstdio>
14
15 using namespace ctml;
16 using namespace std;
17
18 namespace Cantera
19 {
20
21 void importSolution(size_t points,
22  doublereal* oldSoln, IdealGasPhase& oldmech,
23  size_t size_new, doublereal* newSoln, IdealGasPhase& newmech)
24 {
25  // Number of components in old and new solutions
26  size_t nv_old = oldmech.nSpecies() + 4;
27  size_t nv_new = newmech.nSpecies() + 4;
28
29  if (size_new < nv_new*points) {
30  throw CanteraError("importSolution",
31  "new solution array must have length "+
32  int2str(nv_new*points));
33  }
34
35  size_t n, j, knew;
36  string nm;
37
38  // copy u,V,T,lambda
39  for (j = 0; j < points; j++)
40  for (n = 0; n < 4; n++) {
41  newSoln[nv_new*j + n] = oldSoln[nv_old*j + n];
42  }
43
44  // copy mass fractions
45  size_t nsp0 = oldmech.nSpecies();
46  //int nsp1 = newmech.nSpecies();
47
48  // loop over the species in the old mechanism
49  for (size_t k = 0; k < nsp0; k++) {
50  nm = oldmech.speciesName(k); // name
51
52  // location of this species in the new mechanism.
53  // If < 0, then the species is not in the new mechanism.
54  knew = newmech.speciesIndex(nm);
55
56  // copy this species from the old to the new solution vectors
57  if (knew != npos) {
58  for (j = 0; j < points; j++) {
59  newSoln[nv_new*j + 4 + knew] = oldSoln[nv_old*j + 4 + k];
60  }
61  }
62  }
63
64
65  // normalize mass fractions
66  for (j = 0; j < points; j++) {
67  newmech.setMassFractions(&newSoln[nv_new*j + 4]);
68  newmech.getMassFractions(&newSoln[nv_new*j + 4]);
69  }
70 }
71
72 static void st_drawline()
73 {
74  writelog("\n-------------------------------------"
75  "------------------------------------------");
76 }
77
78 StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
79  Domain1D(nsp+4, points),
80  m_inlet_u(0.0),
81  m_inlet_V(0.0),
82  m_inlet_T(-1.0),
83  m_surface_T(-1.0),
84  m_press(-1.0),
85  m_nsp(nsp),
86  m_thermo(0),
87  m_kin(0),
88  m_trans(0),
89  m_jac(0),
90  m_ok(false),
91  m_do_soret(false),
92  m_transport_option(-1)
93 {
94  m_type = cFlowType;
95
96  m_points = points;
97  m_thermo = ph;
98
99  m_zfixed = Undef;
100  m_tfixed = Undef;
101
102  if (ph == 0) {
103  return; // used to create a dummy object
104  }
105
106  size_t nsp2 = m_thermo->nSpecies();
107  if (nsp2 != m_nsp) {
108  m_nsp = nsp2;
109  Domain1D::resize(m_nsp+4, points);
110  }
111
112
113  // make a local copy of the species molecular weight vector
114  m_wt = m_thermo->molecularWeights();
115
116  // the species mass fractions are the last components in the solution
117  // vector, so the total number of components is the number of species
118  // plus the offset of the first mass fraction.
119  m_nv = c_offset_Y + m_nsp;
120
121  // enable all species equations by default
122  m_do_species.resize(m_nsp, true);
123
124  // but turn off the energy equation at all points
125  m_do_energy.resize(m_points,false);
126
127  m_diff.resize(m_nsp*m_points);
128  m_multidiff.resize(m_nsp*m_nsp*m_points);
129  m_flux.resize(m_nsp,m_points);
130  m_wdot.resize(m_nsp,m_points, 0.0);
131  m_surfdot.resize(m_nsp, 0.0);
132  m_ybar.resize(m_nsp);
133
134
135  //-------------- default solution bounds --------------------
136
137  setBounds(0, -1e20, 1e20); // no bounds on u
138  setBounds(1, -1e20, 1e20); // V
139  setBounds(2, 200.0, 1e9); // temperature bounds
140  setBounds(3, -1e20, 1e20); // lamda should be negative
141
142  // mass fraction bounds
143  for (size_t k = 0; k < m_nsp; k++) {
144  setBounds(4+k, -1.0e-5, 1.0e5);
145  }
146
147  //-------------------- default error tolerances ----------------
148  setTransientTolerances(1.0e-8, 1.0e-15);
150
151  //-------------------- grid refinement -------------------------
152  m_refiner->setActive(0, false);
153  m_refiner->setActive(1, false);
154  m_refiner->setActive(2, false);
155  m_refiner->setActive(3, false);
156
157  vector_fp gr;
158  for (size_t ng = 0; ng < m_points; ng++) {
159  gr.push_back(1.0*ng/m_points);
160  }
161  setupGrid(m_points, DATA_PTR(gr));
162  setID("stagnation flow");
163 }
164
165 void StFlow::resize(size_t ncomponents, size_t points)
166 {
167  Domain1D::resize(ncomponents, points);
168  m_rho.resize(m_points, 0.0);
169  m_wtm.resize(m_points, 0.0);
170  m_cp.resize(m_points, 0.0);
171  m_enth.resize(m_points, 0.0);
172  m_visc.resize(m_points, 0.0);
173  m_tcon.resize(m_points, 0.0);
174
175  if (m_transport_option == c_Mixav_Transport) {
176  m_diff.resize(m_nsp*m_points);
177  } else {
178  m_multidiff.resize(m_nsp*m_nsp*m_points);
179  m_diff.resize(m_nsp*m_points);
180  m_dthermal.resize(m_nsp, m_points, 0.0);
181  }
182  m_flux.resize(m_nsp,m_points);
183  m_wdot.resize(m_nsp,m_points, 0.0);
184  m_do_energy.resize(m_points,false);
185
186  m_fixedy.resize(m_nsp, m_points);
187  m_fixedtemp.resize(m_points);
188
189  m_dz.resize(m_points-1);
190  m_z.resize(m_points);
191 }
192
193 void StFlow::setupGrid(size_t n, const doublereal* z)
194 {
195  resize(m_nv, n);
196  size_t j;
197
198  m_z[0] = z[0];
199  for (j = 1; j < m_points; j++) {
200  m_z[j] = z[j];
201  m_dz[j-1] = m_z[j] - m_z[j-1];
202  }
203 }
204
205 void StFlow::setTransport(Transport& trans, bool withSoret)
206 {
207  m_trans = &trans;
208  m_do_soret = withSoret;
209
210  if (m_trans->model() == cMulticomponent) {
211  m_transport_option = c_Multi_Transport;
212  m_multidiff.resize(m_nsp*m_nsp*m_points);
213  m_diff.resize(m_nsp*m_points);
214  m_dthermal.resize(m_nsp, m_points, 0.0);
215  } else if (m_trans->model() == cMixtureAveraged) {
216  m_transport_option = c_Mixav_Transport;
217  m_diff.resize(m_nsp*m_points);
218  if (withSoret)
219  throw CanteraError("setTransport",
220  "Thermal diffusion (the Soret effect) "
221  "requires using a multicomponent transport model.");
222  } else {
223  throw CanteraError("setTransport","unknown transport model.");
224  }
225 }
226
227 void StFlow::enableSoret(bool withSoret)
228 {
229  if (m_transport_option == c_Multi_Transport) {
230  m_do_soret = withSoret;
231  } else {
232  throw CanteraError("setTransport",
233  "Thermal diffusion (the Soret effect) "
234  "requires using a multicomponent transport model.");
235  }
236 }
237
238 void StFlow::setGas(const doublereal* x, size_t j)
239 {
240  m_thermo->setTemperature(T(x,j));
241  const doublereal* yy = x + m_nv*j + c_offset_Y;
242  m_thermo->setMassFractions_NoNorm(yy);
243  m_thermo->setPressure(m_press);
244 }
245
246 void StFlow::setGasAtMidpoint(const doublereal* x, size_t j)
247 {
248  m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1)));
249  const doublereal* yyj = x + m_nv*j + c_offset_Y;
250  const doublereal* yyjp = x + m_nv*(j+1) + c_offset_Y;
251  for (size_t k = 0; k < m_nsp; k++) {
252  m_ybar[k] = 0.5*(yyj[k] + yyjp[k]);
253  }
254  m_thermo->setMassFractions_NoNorm(DATA_PTR(m_ybar));
255  m_thermo->setPressure(m_press);
256 }
257
258 void StFlow::_finalize(const doublereal* x)
259 {
260  size_t k, j;
261  doublereal zz, tt;
262  size_t nz = m_zfix.size();
263  bool e = m_do_energy[0];
264  for (j = 0; j < m_points; j++) {
265  if (e || nz == 0) {
266  setTemperature(j, T(x, j));
267  } else {
268  zz = (z(j) - z(0))/(z(m_points - 1) - z(0));
269  tt = linearInterp(zz, m_zfix, m_tfix);
270  setTemperature(j, tt);
271  }
272  for (k = 0; k < m_nsp; k++) {
273  setMassFraction(j, k, Y(x, k, j));
274  }
275  }
276  if (e) {
277  solveEnergyEqn();
278  }
279
280  // If the domain contains the temperature fixed point, make sure that it
281  // is correctly set. This may be necessary when the grid has been modified
282  // externally.
283  if (m_tfixed != Undef) {
284  bool found_zfix = false;
285  for (size_t j = 0; j < m_points; j++) {
286  if (z(j) == m_zfixed) {
287  found_zfix = true;
288  break;
289  }
290  }
291  if (!found_zfix) {
292  for (size_t j = 0; j < m_points - 1; j++) {
293  // Find where the temperature profile crosses the current
294  // fixed temperature.
295  if ((T(x, j) - m_tfixed) * (T(x, j+1) - m_tfixed) <= 0.0) {
296  m_tfixed = T(x, j+1);
297  m_zfixed = z(j+1);
298  break;
299  }
300  }
301  }
302  }
303 }
304
305 void StFlow::eval(size_t jg, doublereal* xg,
306  doublereal* rg, integer* diagg, doublereal rdt)
307 {
308  // if evaluating a Jacobian, and the global point is outside
309  // the domain of influence for this domain, then skip
310  // evaluating the residual
311  if (jg != npos && (jg + 1 < firstPoint() || jg > lastPoint() + 1)) {
312  return;
313  }
314
315  // if evaluating a Jacobian, compute the steady-state residual
316  if (jg != npos) {
317  rdt = 0.0;
318  }
319
320  // start of local part of global arrays
321  doublereal* x = xg + loc();
322  doublereal* rsd = rg + loc();
323  integer* diag = diagg + loc();
324
325  size_t jmin, jmax;
326
327  if (jg == npos) { // evaluate all points
328  jmin = 0;
329  jmax = m_points - 1;
330  } else { // evaluate points for Jacobian
331  size_t jpt = (jg == 0) ? 0 : jg - firstPoint();
332  jmin = std::max<size_t>(jpt, 1) - 1;
333  jmax = std::min(jpt+1,m_points-1);
334  }
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  size_t j, k;
341
342  //-----------------------------------------------------
343  // update properties
344  //-----------------------------------------------------
345
346  // update thermodynamic and transport properties only if a Jacobian is not
347  // being evaluated
348  if (jg == npos) {
349  updateThermo(x, j0, j1);
350  updateTransport(x, j0, j1);
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
358  //----------------------------------------------------
359  // evaluate the residual equations at all required
360  // grid points
361  //----------------------------------------------------
362
363  doublereal sum, sum2, dtdzj;
364
365  for (j = jmin; j <= jmax; j++) {
366
367
368  //----------------------------------------------
369  // left boundary
370  //----------------------------------------------
371
372  if (j == 0) {
373
374  // these may be modified by a boundary object
375
376  // Continuity. This propagates information right-to-left,
377  // since rho_u at point 0 is dependent on rho_u at point 1,
378  // but not on 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
384  // will modify these equations by subtracting its values
385  // for V, T, and mdot. As a result, these residual equations
386  // will force the solution variables to the values for
387  // the boundary object
388  rsd[index(c_offset_V,0)] = V(x,0);
389  rsd[index(c_offset_T,0)] = T(x,0);
390  rsd[index(c_offset_L,0)] = -rho_u(x,0);
391
392  // The default boundary condition for species is zero
393  // flux. However, the boundary object may modify
394  // this.
395  sum = 0.0;
396  for (k = 0; k < m_nsp; k++) {
397  sum += Y(x,k,0);
398  rsd[index(c_offset_Y + k, 0)] =
399  -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
400  }
401  rsd[index(c_offset_Y, 0)] = 1.0 - sum;
402  }
403
404
405  else if (j == m_points - 1) {
406  evalRightBoundary(x, rsd, diag, rdt);
407
408  } else { // interior points
409  evalContinuity(j, x, rsd, diag, rdt);
410
411  //------------------------------------------------
413  //
414  // \rho dV/dt + \rho u dV/dz + \rho V^2
415  // = d(\mu dV/dz)/dz - lambda
416  //
417  //-------------------------------------------------
418  rsd[index(c_offset_V,j)]
419  = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j)
420  - m_rho[j]*V(x,j)*V(x,j))/m_rho[j]
421  - rdt*(V(x,j) - V_prev(j));
422  diag[index(c_offset_V, j)] = 1;
423
424  //-------------------------------------------------
425  // Species equations
426  //
427  // \rho dY_k/dt + \rho u dY_k/dz + dJ_k/dz
428  // = M_k\omega_k
429  //
430  //-------------------------------------------------
431  getWdot(x,j);
432
433  doublereal convec, diffus;
434  for (k = 0; k < m_nsp; k++) {
435  convec = rho_u(x,j)*dYdz(x,k,j);
436  diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))
437  /(z(j+1) - z(j-1));
438  rsd[index(c_offset_Y + k, j)]
439  = (m_wt[k]*(wdot(k,j))
440  - convec - diffus)/m_rho[j]
441  - rdt*(Y(x,k,j) - Y_prev(k,j));
442  diag[index(c_offset_Y + k, j)] = 1;
443  }
444
445  //-----------------------------------------------
446  // energy equation
447  //
448  // \rho c_p dT/dt + \rho c_p u dT/dz
449  // = d(k dT/dz)/dz
450  // - sum_k(\omega_k h_k_ref)
451  // - sum_k(J_k c_p_k / M_k) dT/dz
452  //-----------------------------------------------
453
454  if (m_do_energy[j]) {
455
456  setGas(x,j);
457
458  // heat release term
459  const vector_fp& h_RT = m_thermo->enthalpy_RT_ref();
460  const vector_fp& cp_R = m_thermo->cp_R_ref();
461
462  sum = 0.0;
463  sum2 = 0.0;
464  doublereal flxk;
465  for (k = 0; k < m_nsp; k++) {
466  flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
467  sum += wdot(k,j)*h_RT[k];
468  sum2 += flxk*cp_R[k]/m_wt[k];
469  }
470  sum *= GasConstant * T(x,j);
471  dtdzj = dTdz(x,j);
472  sum2 *= GasConstant * dtdzj;
473
474  rsd[index(c_offset_T, j)] =
475  - m_cp[j]*rho_u(x,j)*dtdzj
476  - divHeatFlux(x,j) - sum - sum2;
477  rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
478
479  rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
480  diag[index(c_offset_T, j)] = 1;
481  } else {
482  // residual equations if the energy equation is disabled
483  rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
484  diag[index(c_offset_T, j)] = 0;
485  }
486
487  rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
488  diag[index(c_offset_L, j)] = 0;
489  }
490  }
491 }
492
493 void StFlow::updateTransport(doublereal* x, size_t j0, size_t j1)
494 {
495  if (m_transport_option == c_Mixav_Transport) {
496  for (size_t j = j0; j < j1; j++) {
497  setGasAtMidpoint(x,j);
498  m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
499  m_trans->getMixDiffCoeffs(DATA_PTR(m_diff) + j*m_nsp);
500  m_tcon[j] = m_trans->thermalConductivity();
501  }
502  } else if (m_transport_option == c_Multi_Transport) {
503  doublereal sum, sumx, wtm, dz;
504  doublereal eps = 1.0e-12;
505  for (size_t m = j0; m < j1; m++) {
506  setGasAtMidpoint(x,m);
507  dz = m_z[m+1] - m_z[m];
508  wtm = m_thermo->meanMolecularWeight();
509
510  m_visc[m] = (m_dovisc ? m_trans->viscosity() : 0.0);
511
512  m_trans->getMultiDiffCoeffs(m_nsp,
513  DATA_PTR(m_multidiff) + mindex(0,0,m));
514
515  for (size_t k = 0; k < m_nsp; k++) {
516  sum = 0.0;
517  sumx = 0.0;
518  for (size_t j = 0; j < m_nsp; j++) {
519  if (j != k) {
520  sum += m_wt[j]*m_multidiff[mindex(k,j,m)]*
521  ((X(x,j,m+1) - X(x,j,m))/dz + eps);
522  sumx += (X(x,j,m+1) - X(x,j,m))/dz;
523  }
524  }
525  m_diff[k + m*m_nsp] = sum/(wtm*(sumx+eps));
526  }
527
528  m_tcon[m] = m_trans->thermalConductivity();
529  if (m_do_soret) {
530  m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + m*m_nsp);
531  }
532  }
533  }
534 }
535
536 void StFlow::showSolution(const doublereal* x)
537 {
538  size_t nn = m_nv/5;
539  size_t i, j, n;
540  //char* buf = new char[100];
541  char buf[100];
542
543  // The mean molecular weight is needed to convert
544  updateThermo(x, 0, m_points-1);
545
546  sprintf(buf, " Pressure: %10.4g Pa \n", m_press);
547  writelog(buf);
548  for (i = 0; i < nn; i++) {
549  st_drawline();
550  sprintf(buf, "\n z ");
551  writelog(buf);
552  for (n = 0; n < 5; n++) {
553  sprintf(buf, " %10s ",componentName(i*5 + n).c_str());
554  writelog(buf);
555  }
556  st_drawline();
557  for (j = 0; j < m_points; j++) {
558  sprintf(buf, "\n %10.4g ",m_z[j]);
559  writelog(buf);
560  for (n = 0; n < 5; n++) {
561  sprintf(buf, " %10.4g ",component(x, i*5+n,j));
562  writelog(buf);
563  }
564  }
565  writelog("\n");
566  }
567  size_t nrem = m_nv - 5*nn;
568  st_drawline();
569  sprintf(buf, "\n z ");
570  writelog(buf);
571  for (n = 0; n < nrem; n++) {
572  sprintf(buf, " %10s ", componentName(nn*5 + n).c_str());
573  writelog(buf);
574  }
575  st_drawline();
576  for (j = 0; j < m_points; j++) {
577  sprintf(buf, "\n %10.4g ",m_z[j]);
578  writelog(buf);
579  for (n = 0; n < nrem; n++) {
580  sprintf(buf, " %10.4g ",component(x, nn*5+n,j));
581  writelog(buf);
582  }
583  }
584  writelog("\n");
585 }
586
587 void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
588 {
589  size_t j, k, m;
590  doublereal sum, wtm, rho, dz, gradlogT;
591
592  switch (m_transport_option) {
593
594  case c_Mixav_Transport:
595  case c_Multi_Transport:
596  for (j = j0; j < j1; j++) {
597  sum = 0.0;
598  wtm = m_wtm[j];
599  rho = density(j);
600  dz = z(j+1) - z(j);
601
602  for (k = 0; k < m_nsp; k++) {
603  m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
604  m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
605  sum -= m_flux(k,j);
606  }
607  // correction flux to insure that \sum_k Y_k V_k = 0.
608  for (k = 0; k < m_nsp; k++) {
609  m_flux(k,j) += sum*Y(x,k,j);
610  }
611  }
612  break;
613
614  default:
615  throw CanteraError("updateDiffFluxes","unknown transport model");
616  }
617
618  if (m_do_soret) {
619  for (m = j0; m < j1; m++) {
620  gradlogT = 2.0 * (T(x,m+1) - T(x,m)) /
621  ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m)));
622  for (k = 0; k < m_nsp; k++) {
624  }
625  }
626  }
627 }
628
629 string StFlow::componentName(size_t n) const
630 {
631  switch (n) {
632  case 0:
633  return "u";
634  case 1:
635  return "V";
636  case 2:
637  return "T";
638  case 3:
639  return "lambda";
640  default:
641  if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) {
642  return m_thermo->speciesName(n - c_offset_Y);
643  } else {
644  return "<unknown>";
645  }
646  }
647 }
648
649 size_t StFlow::componentIndex(const std::string& name) const
650 {
651  if (name=="u") {
652  return 0;
653  } else if (name=="V") {
654  return 1;
655  } else if (name=="T") {
656  return 2;
657  } else if (name=="lambda") {
658  return 3;
659  } else {
660  for (size_t n=4; n<m_nsp+4; n++) {
661  if (componentName(n)==name) {
662  return n;
663  }
664  }
665  }
666
667  return npos;
668 }
669
670 void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel)
671 {
672  Domain1D::restore(dom, soln, loglevel);
673  vector<string> ignored;
674  size_t nsp = m_thermo->nSpecies();
675  vector_int did_species(nsp, 0);
676
677  vector<XML_Node*> str;
678  dom.getChildren("string",str);
679  int nstr = static_cast<int>(str.size());
680  for (int istr = 0; istr < nstr; istr++) {
681  const XML_Node& nd = *str[istr];
682  writelog(nd["title"]+": "+nd.value()+"\n");
683  }
684
685  //map<string, double> params;
686  double pp = -1.0;
687  pp = getFloat(dom, "pressure", "pressure");
688  setPressure(pp);
689
690  getOptionalFloat(dom, "t_fixed", m_tfixed);
691  getOptionalFloat(dom, "z_fixed", m_zfixed);
692
693  vector<XML_Node*> d;
694  dom.child("grid_data").getChildren("floatArray",d);
695  size_t nd = d.size();
696
697  vector_fp x;
698  size_t n, np = 0, j, ks, k;
699  string nm;
701  for (n = 0; n < nd; n++) {
702  const XML_Node& fa = *d[n];
703  nm = fa["title"];
704  if (nm == "z") {
705  getFloatArray(fa,x,false);
706  np = x.size();
707  writelog("Grid contains "+int2str(np)+" points.\n", loglevel >= 2);
709  setupGrid(np, DATA_PTR(x));
710  }
711  }
713  throw CanteraError("StFlow::restore",
714  "domain contains no grid points.");
715  }
716
717  writelog("Importing datasets:\n", loglevel >= 2);
718  for (n = 0; n < nd; n++) {
719  const XML_Node& fa = *d[n];
720  nm = fa["title"];
721  getFloatArray(fa,x,false);
722  if (nm == "u") {
723  writelog("axial velocity ", loglevel >= 2);
724  if (x.size() != np) {
725  throw CanteraError("StFlow::restore",
726  "axial velocity array size error");
727  }
728  for (j = 0; j < np; j++) {
729  soln[index(0,j)] = x[j];
730  }
731  } else if (nm == "z") {
733  } else if (nm == "V") {
734  writelog("radial velocity ", loglevel >= 2);
735  if (x.size() != np) {
736  throw CanteraError("StFlow::restore",
737  "radial velocity array size error");
738  }
739  for (j = 0; j < np; j++) {
740  soln[index(1,j)] = x[j];
741  }
742  } else if (nm == "T") {
743  writelog("temperature ", loglevel >= 2);
744  if (x.size() != np) {
745  throw CanteraError("StFlow::restore",
746  "temperature array size error");
747  }
748  for (j = 0; j < np; j++) {
749  soln[index(2,j)] = x[j];
750  }
751
752  // For fixed-temperature simulations, use the
753  // imported temperature profile by default. If
754  // this is not desired, call setFixedTempProfile
755  // *after* restoring the solution.
756
757  vector_fp zz(np);
758  for (size_t jj = 0; jj < np; jj++) {
759  zz[jj] = (grid(jj) - zmin())/(zmax() - zmin());
760  }
761  setFixedTempProfile(zz, x);
762  } else if (nm == "L") {
763  writelog("lambda ", loglevel >= 2);
764  if (x.size() != np) {
765  throw CanteraError("StFlow::restore",
766  "lambda arary size error");
767  }
768  for (j = 0; j < np; j++) {
769  soln[index(3,j)] = x[j];
770  }
771  } else if (m_thermo->speciesIndex(nm) != npos) {
772  writelog(nm+" ", loglevel >= 2);
773  if (x.size() == np) {
774  k = m_thermo->speciesIndex(nm);
775  did_species[k] = 1;
776  for (j = 0; j < np; j++) {
777  soln[index(k+4,j)] = x[j];
778  }
779  }
780  } else {
781  ignored.push_back(nm);
782  }
783  }
784
785  if (loglevel >=2 && !ignored.empty()) {
786  writelog("\n\n");
787  writelog("Ignoring datasets:\n");
788  size_t nn = ignored.size();
789  for (size_t n = 0; n < nn; n++) {
790  writelog(ignored[n]+" ");
791  }
792  }
793
794  if (loglevel >= 1) {
795  for (ks = 0; ks < nsp; ks++) {
796  if (did_species[ks] == 0) {
798  writelog("Missing data for species:\n");
800  }
801  writelog(m_thermo->speciesName(ks)+" ");
802  }
803  }
804  }
805
806  if (dom.hasChild("energy_enabled")) {
807  getFloatArray(dom, x, false, "", "energy_enabled");
808  if (x.size() == nPoints()) {
809  for (size_t i = 0; i < x.size(); i++) {
810  m_do_energy[i] = x[i];
811  }
812  } else if (!x.empty()) {
813  throw CanteraError("StFlow::restore", "energy_enabled is length" +
814  int2str(x.size()) + "but should be length" +
815  int2str(nPoints()));
816  }
817  }
818
819  if (dom.hasChild("species_enabled")) {
820  getFloatArray(dom, x, false, "", "species_enabled");
821  if (x.size() == m_nsp) {
822  for (size_t i = 0; i < x.size(); i++) {
823  m_do_species[i] = x[i];
824  }
825  } else if (!x.empty()) {
826  // This may occur when restoring from a mechanism with a different
827  // number of species.
828  if (loglevel > 0) {
829  writelog("\nWarning: StFlow::restore: species_enabled is length " +
830  int2str(x.size()) + " but should be length " +
831  int2str(m_nsp) + ". Enabling all species equations by default.");
832  }
833  m_do_species.assign(m_nsp, true);
834  }
835  }
836
837  if (dom.hasChild("refine_criteria")) {
838  XML_Node& ref = dom.child("refine_criteria");
839  refiner().setCriteria(getFloat(ref, "ratio"), getFloat(ref, "slope"),
840  getFloat(ref, "curve"), getFloat(ref, "prune"));
841  refiner().setGridMin(getFloat(ref, "grid_min"));
842  }
843 }
844
845 XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol)
846 {
847  size_t k;
848
849  Array2D soln(m_nv, m_points, sol + loc());
850
851  XML_Node& flow = Domain1D::save(o, sol);
853
854  if (m_desc != "") {
856  }
858  addFloat(flow, "pressure", m_press, "Pa", "pressure");
859
860  if (m_zfixed != Undef) {
863  }
864
866  "m","length");
867  vector_fp x(soln.nColumns());
868
869  soln.getRow(0,DATA_PTR(x));
871
872  soln.getRow(1,DATA_PTR(x));
874  x.size(),DATA_PTR(x),"1/s","rate");
875
876  soln.getRow(2,DATA_PTR(x));
878
879  soln.getRow(3,DATA_PTR(x));
881
882  for (k = 0; k < m_nsp; k++) {
883  soln.getRow(4+k,DATA_PTR(x));
885  x.size(),DATA_PTR(x),"","massFraction");
886  }
887
888  vector_fp values(nPoints());
889  for (size_t i = 0; i < nPoints(); i++) {
890  values[i] = m_do_energy[i];
891  }
893
894  values.resize(m_nsp);
895  for (size_t i = 0; i < m_nsp; i++) {
896  values[i] = m_do_species[i];
897  }
899
906
907  return flow;
908 }
909
910 void StFlow::setJac(MultiJac* jac)
911 {
912  m_jac = jac;
913 }
914
915 void AxiStagnFlow::evalRightBoundary(doublereal* x, doublereal* rsd,
916  integer* diag, doublereal rdt)
917 {
918  size_t j = m_points - 1;
919  // the boundary object connected to the right of this one may modify or
920  // replace these equations. The default boundary conditions are zero u, V,
921  // and T, and zero diffusive flux for all species.
922
923  rsd[index(0,j)] = rho_u(x,j);
924  rsd[index(1,j)] = V(x,j);
925  rsd[index(2,j)] = T(x,j);
926  rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
927  diag[index(c_offset_L, j)] = 0;
928  doublereal sum = 0.0;
929  for (size_t k = 0; k < m_nsp; k++) {
930  sum += Y(x,k,j);
931  rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
932  }
933  rsd[index(4,j)] = 1.0 - sum;
934  diag[index(4,j)] = 0;
935 }
936
937 void AxiStagnFlow::evalContinuity(size_t j, doublereal* x, doublereal* rsd,
938  integer* diag, doublereal rdt)
939 {
940  //----------------------------------------------
941  // Continuity equation
942  //
943  // Note that this propagates the mass flow rate information to the left
944  // (j+1 -> j) from the value specified at the right boundary. The
945  // lambda information propagates in the opposite direction.
946  //
947  // d(\rho u)/dz + 2\rho V = 0
948  //
949  //------------------------------------------------
950
951  rsd[index(c_offset_U,j)] =
952  -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
953  -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
954
955  //algebraic constraint
956  diag[index(c_offset_U, j)] = 0;
957 }
958
959 void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd,
960  integer* diag, doublereal rdt)
961 {
962  size_t j = m_points - 1;
963
964  // the boundary object connected to the right of this one may modify or
965  // replace these equations. The default boundary conditions are zero u, V,
966  // and T, and zero diffusive flux for all species.
967
969  rsd[index(0,j)] = rho_u(x,j) - rho_u(x,j-1);
970  rsd[index(1,j)] = V(x,j);
971  rsd[index(2,j)] = T(x,j) - T(x,j-1);
972  doublereal sum = 0.0;
973  rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
974  diag[index(c_offset_L, j)] = 0;
975  for (size_t k = 0; k < m_nsp; k++) {
976  sum += Y(x,k,j);
977  rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
978  }
979  rsd[index(4,j)] = 1.0 - sum;
980  diag[index(4,j)] = 0;
981 }
982
983 void FreeFlame::evalContinuity(size_t j, doublereal* x, doublereal* rsd,
984  integer* diag, doublereal rdt)
985 {
986  //----------------------------------------------
987  // Continuity equation
988  //
989  // d(\rho u)/dz + 2\rho V = 0
990  //
991  //----------------------------------------------
992
993  if (grid(j) > m_zfixed) {
994  rsd[index(c_offset_U,j)] =
995  - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]
996  - (density(j-1)*V(x,j-1) + density(j)*V(x,j));
997  }
998
999  else if (grid(j) == m_zfixed) {
1000  if (m_do_energy[j]) {
1001  rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
1002  } else {
1003  rsd[index(c_offset_U,j)] = (rho_u(x,j)
1004  - m_rho[0]*0.3);
1005  }
1006  } else if (grid(j) < m_zfixed) {
1007  rsd[index(c_offset_U,j)] =
1008  - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
1009  - (density(j+1)*V(x,j+1) + density(j)*V(x,j));
1010  }
1011  //algebraic constraint
1012  diag[index(c_offset_U, j)] = 0;
1013 }
1014
1015 } // namespace
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
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:315
size_t nPoints() const
Number of grid points in this domain.
Definition: Domain1D.h:186
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
virtual void evalRightBoundary(doublereal *x, doublereal *res, integer *diag, doublereal rdt)=0
Evaluate all residual components at the right boundary.
virtual int model() const
Transport model.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition: Array.h:128
void getMassFractions(doublereal *const y) const
Get the species mass fractions.
Definition: Phase.cpp:561
virtual void evalRightBoundary(doublereal *x, doublereal *res, integer *diag, doublereal rdt)
Evaluate all residual components at the right boundary.
Definition: StFlow.cpp:959
Class IdealGasPhase represents low-density gases that obey the ideal gas equation of state...
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
size_t firstPoint() const
The index of the first (i.e., left-most) grid point belonging to this domain.
Definition: Domain1D.h:469
void addNamedFloatArray(Cantera::XML_Node &node, const std::string &name, const int n, const doublereal *const vals, const std::string units, const std::string type, const doublereal minval, const doublereal maxval)
This function adds a child node with the name given by the first parameter with a value consisting of...
Definition: ctml.cpp:102
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:258
virtual void evalContinuity(size_t j, doublereal *x, doublereal *r, integer *diag, doublereal rdt)
Evaluate the residual corresponding to the continuity equation at all interior grid points...
Definition: StFlow.cpp:937
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:461
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
bool getOptionalFloat(const Cantera::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:335
Base class for transport property managers.
void setPressure(doublereal p)
Set the pressure.
Definition: StFlow.h:93
doublereal getFloat(const Cantera::XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
Definition: ctml.cpp:267
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
Definition: Domain1D.cpp:167
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:363
virtual void evalContinuity(size_t j, doublereal *x, doublereal *r, integer *diag, doublereal rdt)=0
Evaluate the residual corresponding to the continuity equation at all interior grid points...
void addString(Cantera::XML_Node &node, const std::string &titleString, const std::string &valueString, const std::string &typeString)
This function adds a child node with the name string with a string value to the current node...
Definition: ctml.cpp:141
XML_Node & child(const size_t n) const
Return a changeable reference to the n'th child of the current node.
Definition: xml.cpp:584
const doublereal Undef
Fairly random number to be used to initialize variables against to see if they are subsequently defin...
Definition: ct_defs.h:147
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:33
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:154
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
Definition: StFlow.cpp:670
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:167
void setBounds(size_t nl, const doublereal *lower, size_t nu, const doublereal *upper)
Set the lower and upper bounds for each solution component.
Definition: Domain1D.h:240
Base class for one-dimensional domains.
Definition: Domain1D.h:37
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:229
doublereal T_fixed(size_t j) const
The fixed temperature value at point j.
Definition: StFlow.h:150
void setCriteria(doublereal ratio=10.0, doublereal slope=0.8, doublereal curve=0.8, doublereal prune=-0.1)
Set grid refinement criteria.
Definition: refine.h:31
size_t lastPoint() const
The index of the last (i.e., right-most) grid point belonging to this domain.
Definition: Domain1D.h:477
void setID(const std::string &s)
Specify an identifying tag for this domain.
Definition: Domain1D.h:517
virtual void evalRightBoundary(doublereal *x, doublereal *res, integer *diag, doublereal rdt)
Evaluate all residual components at the right boundary.
Definition: StFlow.cpp:915
size_t nColumns() const
Number of columns.
Definition: Array.h:322
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:125
void addFloatArray(Cantera::XML_Node &node, const std::string &title, const size_t n, const doublereal *const vals, const std::string &units, const std::string &type, const doublereal minval, const doublereal maxval)
This function adds a child node with the name, "floatArray", with a value consisting of a comma separ...
Definition: ctml.cpp:68
void addFloat(Cantera::XML_Node &node, const std::string &title, const doublereal val, const std::string &units, const std::string &type, const doublereal minval, const doublereal maxval)
This function adds a child node with the name, "float", with a value consisting of a single floating ...
Definition: ctml.cpp:42
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
void setTemperature(size_t j, doublereal t)
Definition: StFlow.h:134
virtual doublereal thermalConductivity()
Returns the mixture thermal conductivity in W/m/K.
virtual doublereal viscosity()
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:574
virtual void resize(size_t nv, size_t np)
Definition: Domain1D.h:135
void setGridMin(double gridmin)
Set the minimum allowable spacing between adjacent grid points [m].
Definition: refine.h:50
void addAttribute(const std::string &attrib, const std::string &value)
Add or modify an attribute of the current node.
Definition: xml.cpp:518
Refiner & refiner()
Return a reference to the grid refiner.
Definition: Domain1D.h:159
void getWdot(doublereal *x, size_t j)
Write the net production rates at point j into array m_wdot
Definition: StFlow.h:306
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition: StFlow.cpp:629
virtual void setMassFractions_NoNorm(const doublereal *const y)
Set the mass fractions to the specified values without normalizing.
Definition: Phase.cpp:388
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:252
virtual void evalContinuity(size_t j, doublereal *x, doublereal *r, integer *diag, doublereal rdt)
Evaluate the residual corresponding to the continuity equation at all interior grid points...
Definition: StFlow.cpp:983
Class MultiJac evaluates the Jacobian of a system of equations defined by a residual function supplie...
Definition: MultiJac.h:26
virtual void getMultiDiffCoeffs(const size_t ld, doublereal *const d)
Return the Multicomponent diffusion coefficients. Units: [m^2/s].
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:238
std::string value() const
Return the value of an XML node as a string.
Definition: xml.cpp:488
virtual void showSolution(const doublereal *x)
Print the solution.
Definition: StFlow.cpp:536
const vector_fp & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:505
void setSteadyTolerances(doublereal rtol, doublereal atol, size_t n=npos)
Definition: Domain1D.cpp:88
const vector_fp & cp_R_ref() const
Returns a reference to the dimensionless reference state Heat Capacity vector.
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:493
virtual void getMixDiffCoeffs(doublereal *const d)
Returns a vector of mixture averaged diffusion coefficients.
virtual void setupGrid(size_t n, const doublereal *z)
called to set up initial grid, and after grid refinement
Definition: StFlow.cpp:193
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:165
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:246
void setTransport(Transport &trans, bool withSoret=false)
set the transport manager
Definition: StFlow.cpp:205
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:562
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:588
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:66
const vector_fp & enthalpy_RT_ref() const
Returns a reference to the dimensionless reference state enthalpy vector.
void updateDiffFluxes(const doublereal *x, size_t j0, size_t j1)
Update the diffusive mass fluxes.
Definition: StFlow.cpp:587
void importSolution(size_t points, doublereal *oldSoln, IdealGasPhase &oldmech, size_t size_new, doublereal *newSoln, IdealGasPhase &newmech)
Import a previous solution to use as an initial estimate.
Definition: StFlow.cpp:21
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
virtual void setMassFractions(const doublereal *const y)
Set the mass fractions to the specified values and normalize them.
Definition: Phase.cpp:374
virtual void getThermalDiffCoeffs(doublereal *const dt)
Return a vector of Thermal diffusion coefficients [kg/m/sec].
void getRow(size_t n, doublereal *const rw)
Get the nth row and return it in a vector.
Definition: Array.h:195
void writelog(const std::string &msg)
Write a message to the screen.
Definition: global.cpp:43
void resize(size_t components, size_t points)
Change the grid size. Called after grid refinement.
Definition: StFlow.cpp:165
void setMassFraction(size_t j, size_t k, doublereal y)
Definition: StFlow.h:144
void setTransientTolerances(doublereal rtol, doublereal atol, size_t n=npos)
Set tolerances for time-stepping mode.
Definition: Domain1D.cpp:68
virtual void eval(size_t j, doublereal *x, doublereal *r, integer *mask, doublereal rdt)
Definition: StFlow.cpp:305
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:246
size_t getFloatArray(const Cantera::XML_Node &node, std::vector< doublereal > &v, const bool convert, const std::string &unitsString, const std::string &nodeName)
This function reads the current node or a child node of the current node with the default name...
Definition: ctml.cpp:419
void getChildren(const std::string &name, std::vector< XML_Node * > &children) const
Get a vector of pointers to XML_Node containing all of the children of the current node which matches...
Definition: xml.cpp:916
doublereal linearInterp(doublereal x, const vector_fp &xpts, const vector_fp &fpts)
Linearly interpolate a function defined on a discrete grid.
Definition: funcs.cpp:40
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:845
virtual void setPressure(doublereal p)
Set the pressure at constant temperature and composition.