Cantera  2.0
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 
22 //------------------- importSolution ------------------------
23 
24 /**
25  * Import a previous solution to use as an initial estimate. The
26  * previous solution may have been computed using a different
27  * reaction mechanism. Species in the old and new mechanisms are
28  * matched by name, and any species in the new mechanism that were
29  * not in the old one are set to zero. The new solution is created
30  * with the same number of grid points as in the old solution.
31  */
32 void importSolution(size_t points,
33  doublereal* oldSoln, IdealGasPhase& oldmech,
34  size_t size_new, doublereal* newSoln, IdealGasPhase& newmech)
35 {
36 
37  // Number of components in old and new solutions
38  size_t nv_old = oldmech.nSpecies() + 4;
39  size_t nv_new = newmech.nSpecies() + 4;
40 
41  if (size_new < nv_new*points) {
42  throw CanteraError("importSolution",
43  "new solution array must have length "+
44  int2str(nv_new*points));
45  }
46 
47  size_t n, j, knew;
48  string nm;
49 
50  // copy u,V,T,lambda
51  for (j = 0; j < points; j++)
52  for (n = 0; n < 4; n++) {
53  newSoln[nv_new*j + n] = oldSoln[nv_old*j + n];
54  }
55 
56  // copy mass fractions
57  size_t nsp0 = oldmech.nSpecies();
58  //int nsp1 = newmech.nSpecies();
59 
60  // loop over the species in the old mechanism
61  for (size_t k = 0; k < nsp0; k++) {
62  nm = oldmech.speciesName(k); // name
63 
64  // location of this species in the new mechanism.
65  // If < 0, then the species is not in the new mechanism.
66  knew = newmech.speciesIndex(nm);
67 
68  // copy this species from the old to the new solution vectors
69  if (knew != npos) {
70  for (j = 0; j < points; j++) {
71  newSoln[nv_new*j + 4 + knew] = oldSoln[nv_old*j + 4 + k];
72  }
73  }
74  }
75 
76 
77  // normalize mass fractions
78  for (j = 0; j < points; j++) {
79  newmech.setMassFractions(&newSoln[nv_new*j + 4]);
80  newmech.getMassFractions(&newSoln[nv_new*j + 4]);
81  }
82 }
83 
84 
85 static void st_drawline()
86 {
87  writelog("\n-------------------------------------"
88  "------------------------------------------");
89 }
90 
91 StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
92  Domain1D(nsp+4, points),
93  m_inlet_u(0.0),
94  m_inlet_V(0.0),
95  m_inlet_T(-1.0),
96  m_surface_T(-1.0),
97  m_press(-1.0),
98  m_nsp(nsp),
99  m_thermo(0),
100  m_kin(0),
101  m_trans(0),
102  m_jac(0),
103  m_ok(false),
104  m_do_soret(false),
105  m_transport_option(-1)
106 {
107  m_type = cFlowType;
108 
109  m_points = points;
110  m_thermo = ph;
111 
112  if (ph == 0) {
113  return; // used to create a dummy object
114  }
115 
116  size_t nsp2 = m_thermo->nSpecies();
117  if (nsp2 != m_nsp) {
118  m_nsp = nsp2;
119  Domain1D::resize(m_nsp+4, points);
120  }
121 
122 
123  // make a local copy of the species molecular weight vector
124  m_wt = m_thermo->molecularWeights();
125 
126  // the species mass fractions are the last components in the solution
127  // vector, so the total number of components is the number of species
128  // plus the offset of the first mass fraction.
129  m_nv = c_offset_Y + m_nsp;
130 
131  // enable all species equations by default
132  m_do_species.resize(m_nsp, true);
133 
134  // but turn off the energy equation at all points
135  m_do_energy.resize(m_points,false);
136 
137  m_diff.resize(m_nsp*m_points);
138  m_multidiff.resize(m_nsp*m_nsp*m_points);
139  m_flux.resize(m_nsp,m_points);
140  m_wdot.resize(m_nsp,m_points, 0.0);
141  m_surfdot.resize(m_nsp, 0.0);
142  m_ybar.resize(m_nsp);
143 
144 
145  //-------------- default solution bounds --------------------
146 
147  vector_fp vmin(m_nv), vmax(m_nv);
148 
149  // no bounds on u
150  vmin[0] = -1.e20;
151  vmax[0] = 1.e20;
152 
153  // V
154  vmin[1] = -1.e20;
155  vmax[1] = 1.e20;
156 
157  // temperature bounds
158  vmin[2] = 200.0;
159  vmax[2]= 1.e9;
160 
161  // lamda should be negative
162  vmin[3] = -1.e20;
163  vmax[3] = 1.e20;
164 
165  // mass fraction bounds
166  for (size_t k = 0; k < m_nsp; k++) {
167  vmin[4+k] = -1.0e-5;
168  vmax[4+k] = 1.0e5;
169  }
170  setBounds(vmin.size(), DATA_PTR(vmin), vmax.size(), DATA_PTR(vmax));
171 
172 
173  //-------------------- default error tolerances ----------------
174  vector_fp rtol(m_nv, 1.0e-8);
175  vector_fp atol(m_nv, 1.0e-15);
176  setTolerances(rtol.size(), DATA_PTR(rtol), atol.size(), DATA_PTR(atol),false);
177  setTolerances(rtol.size(), DATA_PTR(rtol), atol.size(), DATA_PTR(atol),true);
178 
179  //-------------------- grid refinement -------------------------
180  m_refiner->setActive(0, false);
181  m_refiner->setActive(1, false);
182  m_refiner->setActive(2, false);
183  m_refiner->setActive(3, false);
184 
185  vector_fp gr;
186  for (size_t ng = 0; ng < m_points; ng++) {
187  gr.push_back(1.0*ng/m_points);
188  }
189  setupGrid(m_points, DATA_PTR(gr));
190  setID("stagnation flow");
191 }
192 
193 
194 /**
195  * Change the grid size. Called after grid refinement.
196  */
197 void StFlow::resize(size_t ncomponents, size_t points)
198 {
199  Domain1D::resize(ncomponents, points);
200  m_rho.resize(m_points, 0.0);
201  m_wtm.resize(m_points, 0.0);
202  m_cp.resize(m_points, 0.0);
203  m_enth.resize(m_points, 0.0);
204  m_visc.resize(m_points, 0.0);
205  m_tcon.resize(m_points, 0.0);
206 
207  if (m_transport_option == c_Mixav_Transport) {
208  m_diff.resize(m_nsp*m_points);
209  } else {
210  m_multidiff.resize(m_nsp*m_nsp*m_points);
211  m_diff.resize(m_nsp*m_points);
212  m_dthermal.resize(m_nsp, m_points, 0.0);
213  }
214  m_flux.resize(m_nsp,m_points);
215  m_wdot.resize(m_nsp,m_points, 0.0);
216  m_do_energy.resize(m_points,false);
217 
218  m_fixedy.resize(m_nsp, m_points);
219  m_fixedtemp.resize(m_points);
220 
221  m_dz.resize(m_points-1);
222  m_z.resize(m_points);
223 }
224 
225 
226 void StFlow::setupGrid(size_t n, const doublereal* z)
227 {
228  resize(m_nv, n);
229  size_t j;
230 
231  m_z[0] = z[0];
232  for (j = 1; j < m_points; j++) {
233  m_z[j] = z[j];
234  m_dz[j-1] = m_z[j] - m_z[j-1];
235  }
236 }
237 
238 
239 /**
240  * Install a transport manager.
241  */
242 void StFlow::setTransport(Transport& trans, bool withSoret)
243 {
244  m_trans = &trans;
245  m_do_soret = withSoret;
246 
247  if (m_trans->model() == cMulticomponent) {
248  m_transport_option = c_Multi_Transport;
249  m_multidiff.resize(m_nsp*m_nsp*m_points);
250  m_diff.resize(m_nsp*m_points);
251  m_dthermal.resize(m_nsp, m_points, 0.0);
252  } else if (m_trans->model() == cMixtureAveraged) {
253  m_transport_option = c_Mixav_Transport;
254  m_diff.resize(m_nsp*m_points);
255  if (withSoret)
256  throw CanteraError("setTransport",
257  "Thermal diffusion (the Soret effect) "
258  "requires using a multicomponent transport model.");
259  } else {
260  throw CanteraError("setTransport","unknown transport model.");
261  }
262 }
263 
264 void StFlow::enableSoret(bool withSoret)
265 {
266  if (m_transport_option == c_Multi_Transport) {
267  m_do_soret = withSoret;
268  } else {
269  throw CanteraError("setTransport",
270  "Thermal diffusion (the Soret effect) "
271  "requires using a multicomponent transport model.");
272  }
273 }
274 
275 
276 /**
277  * Set the gas object state to be consistent with the solution at
278  * point j.
279  */
280 void StFlow::setGas(const doublereal* x, size_t j)
281 {
282  m_thermo->setTemperature(T(x,j));
283  const doublereal* yy = x + m_nv*j + c_offset_Y;
284  m_thermo->setMassFractions_NoNorm(yy);
285  m_thermo->setPressure(m_press);
286 }
287 
288 
289 /**
290  * Set the gas state to be consistent with the solution at the
291  * midpoint between j and j + 1.
292  */
293 void StFlow::setGasAtMidpoint(const doublereal* x, size_t j)
294 {
295  m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1)));
296  const doublereal* yyj = x + m_nv*j + c_offset_Y;
297  const doublereal* yyjp = x + m_nv*(j+1) + c_offset_Y;
298  for (size_t k = 0; k < m_nsp; k++) {
299  m_ybar[k] = 0.5*(yyj[k] + yyjp[k]);
300  }
301  m_thermo->setMassFractions_NoNorm(DATA_PTR(m_ybar));
302  m_thermo->setPressure(m_press);
303 }
304 
305 
306 void StFlow::_finalize(const doublereal* x)
307 {
308  size_t k, j;
309  doublereal zz, tt;
310  size_t nz = m_zfix.size();
311  bool e = m_do_energy[0];
312  for (j = 0; j < m_points; j++) {
313  if (e || nz == 0) {
314  setTemperature(j, T(x, j));
315  } else {
316  zz = (z(j) - z(0))/(z(m_points - 1) - z(0));
317  tt = linearInterp(zz, m_zfix, m_tfix);
318  setTemperature(j, tt);
319  }
320  for (k = 0; k < m_nsp; k++) {
321  setMassFraction(j, k, Y(x, k, j));
322  }
323  }
324  if (e) {
325  solveEnergyEqn();
326  }
327 }
328 
329 
330 //------------------------------------------------------
331 
332 /**
333  * Evaluate the residual function for axisymmetric stagnation
334  * flow. If jpt is less than zero, the residual function is
335  * evaluated at all grid points. If jpt >= 0, then the residual
336  * function is only evaluated at grid points jpt-1, jpt, and
337  * jpt+1. This option is used to efficiently evaluate the
338  * Jacobian numerically.
339  *
340  */
341 
342 void AxiStagnFlow::eval(size_t jg, doublereal* xg,
343  doublereal* rg, integer* diagg, doublereal rdt)
344 {
345 
346  // if evaluating a Jacobian, and the global point is outside
347  // the domain of influence for this domain, then skip
348  // evaluating the residual
349  if (jg != npos && (jg + 1 < firstPoint() || jg > lastPoint() + 1)) {
350  return;
351  }
352 
353  // if evaluating a Jacobian, compute the steady-state residual
354  if (jg != npos) {
355  rdt = 0.0;
356  }
357 
358  // start of local part of global arrays
359  doublereal* x = xg + loc();
360  doublereal* rsd = rg + loc();
361  integer* diag = diagg + loc();
362 
363  size_t jmin, jmax;
364 
365 
366  if (jg == npos) { // evaluate all points
367  jmin = 0;
368  jmax = m_points - 1;
369  } else { // evaluate points for Jacobian
370  size_t jpt = (jg == 0) ? 0 : jg - firstPoint();
371  jmin = std::max<size_t>(jpt, 1) - 1;
372  jmax = std::min(jpt+1,m_points-1);
373  }
374 
375  // properties are computed for grid points from j0 to j1
376  size_t j0 = std::max<size_t>(jmin, 1) - 1;
377  size_t j1 = std::min(jmax+1,m_points-1);
378 
379  size_t j, k;
380 
381  //-----------------------------------------------------
382  // update properties
383  //-----------------------------------------------------
384 
385  // update thermodynamic properties only if a Jacobian is not
386  // being evaluated
387  if (jg == npos) {
388  updateThermo(x, j0, j1);
389 
390  // update transport properties only if a Jacobian is not being
391  // evaluated
392  updateTransport(x, j0, j1);
393  }
394 
395  // update the species diffusive mass fluxes whether or not a
396  // Jacobian is being evaluated
397  updateDiffFluxes(x, j0, j1);
398 
399 
400  //----------------------------------------------------
401  // evaluate the residual equations at all required
402  // grid points
403  //----------------------------------------------------
404 
405  doublereal sum, sum2, dtdzj;
406 
407  for (j = jmin; j <= jmax; j++) {
408 
409 
410  //----------------------------------------------
411  // left boundary
412  //----------------------------------------------
413 
414  if (j == 0) {
415 
416  // these may be modified by a boundary object
417 
418 
419  // Continuity. This propagates information right-to-left,
420  // since rho_u at point 0 is dependent on rho_u at point 1,
421  // but not on mdot from the inlet.
422  rsd[index(c_offset_U,0)] =
423  -(rho_u(x,1) - rho_u(x,0))/m_dz[0]
424  -(density(1)*V(x,1) + density(0)*V(x,0));
425 
426  // the inlet (or other) object connected to this one
427  // will modify these equations by subtracting its values
428  // for V, T, and mdot. As a result, these residual equations
429  // will force the solution variables to the values for
430  // the boundary object
431  rsd[index(c_offset_V,0)] = V(x,0);
432  rsd[index(c_offset_T,0)] = T(x,0);
433  rsd[index(c_offset_L,0)] = -rho_u(x,0);
434 
435  // The default boundary condition for species is zero
436  // flux. However, the boundary object may modify
437  // this.
438  sum = 0.0;
439  for (k = 0; k < m_nsp; k++) {
440  sum += Y(x,k,0);
441  rsd[index(c_offset_Y + k, 0)] =
442  -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
443  }
444  rsd[index(c_offset_Y, 0)] = 1.0 - sum;
445  }
446 
447 
448  //----------------------------------------------
449  //
450  // right boundary
451  //
452  //----------------------------------------------
453 
454  else if (j == m_points - 1) {
455 
456  // the boundary object connected to the right of this
457  // one may modify or replace these equations. The
458  // default boundary conditions are zero u, V, and T,
459  // and zero diffusive flux for all species.
460 
461  rsd[index(0,j)] = rho_u(x,j);
462  rsd[index(1,j)] = V(x,j);
463  rsd[index(2,j)] = T(x,j);
464  rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
465  diag[index(c_offset_L, j)] = 0;
466  doublereal sum = 0.0;
467  for (k = 0; k < m_nsp; k++) {
468  sum += Y(x,k,j);
469  rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
470  }
471  rsd[index(4,j)] = 1.0 - sum;
472  diag[index(4,j)] = 0;
473 
474  }
475 
476 
477  //------------------------------------------
478  // interior points
479  //------------------------------------------
480 
481  else {
482 
483  //----------------------------------------------
484  // Continuity equation
485  //
486  // Note that this propagates the mass flow rate
487  // information to the left (j+1 -> j) from the
488  // value specified at the right boundary. The
489  // lambda information propagates in the opposite
490  // direction.
491  //
492  // d(\rho u)/dz + 2\rho V = 0
493  //
494  //------------------------------------------------
495 
496  rsd[index(c_offset_U,j)] =
497  -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
498  -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
499 
500  //algebraic constraint
501  diag[index(c_offset_U, j)] = 0;
502 
503 
504  //------------------------------------------------
505  // Radial momentum equation
506  //
507  // \rho u dV/dz + \rho V^2 = d(\mu dV/dz)/dz - lambda
508  //
509  //-------------------------------------------------
510  rsd[index(c_offset_V,j)]
511  = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j)
512  - m_rho[j]*V(x,j)*V(x,j))/m_rho[j]
513  - rdt*(V(x,j) - V_prev(j));
514  diag[index(c_offset_V, j)] = 1;
515 
516 
517  //-------------------------------------------------
518  // Species equations
519  //
520  // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k
521  //
522  //-------------------------------------------------
523  getWdot(x,j);
524 
525  doublereal convec, diffus;
526  for (k = 0; k < m_nsp; k++) {
527  convec = rho_u(x,j)*dYdz(x,k,j);
528  diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))
529  /(z(j+1) - z(j-1));
530  rsd[index(c_offset_Y + k, j)]
531  = (m_wt[k]*(wdot(k,j))
532  - convec - diffus)/m_rho[j]
533  - rdt*(Y(x,k,j) - Y_prev(k,j));
534  diag[index(c_offset_Y + k, j)] = 1;
535  }
536 
537 
538  //-----------------------------------------------
539  // energy equation
540  //-----------------------------------------------
541 
542  if (m_do_energy[j]) {
543 
544  setGas(x,j);
545 
546  // heat release term
547  const vector_fp& h_RT = m_thermo->enthalpy_RT_ref();
548  const vector_fp& cp_R = m_thermo->cp_R_ref();
549 
550  sum = 0.0;
551  sum2 = 0.0;
552  doublereal flxk;
553  for (k = 0; k < m_nsp; k++) {
554  flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
555  sum += wdot(k,j)*h_RT[k];
556  sum2 += flxk*cp_R[k]/m_wt[k];
557  }
558  sum *= GasConstant * T(x,j);
559  dtdzj = dTdz(x,j);
560  sum2 *= GasConstant * dtdzj;
561 
562  rsd[index(c_offset_T, j)] =
563  - m_cp[j]*rho_u(x,j)*dtdzj
564  - divHeatFlux(x,j) - sum - sum2;
565  rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
566 
567  rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
568  diag[index(c_offset_T, j)] = 1;
569  }
570 
571  // residual equations if the energy equation is disabled
572 
573  if (!m_do_energy[j]) {
574  rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
575  diag[index(c_offset_T, j)] = 0;
576  }
577 
578  rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
579  diag[index(c_offset_L, j)] = 0;
580  }
581  }
582 }
583 
584 
585 
586 /**
587  * Update the transport properties at grid points in the range
588  * from j0 to j1, based on solution x.
589  */
590 void StFlow::updateTransport(doublereal* x, size_t j0, size_t j1)
591 {
592  if (m_transport_option == c_Mixav_Transport) {
593  for (size_t j = j0; j < j1; j++) {
594  setGasAtMidpoint(x,j);
595  m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0);
596  m_trans->getMixDiffCoeffs(DATA_PTR(m_diff) + j*m_nsp);
597  m_tcon[j] = m_trans->thermalConductivity();
598  }
599  } else if (m_transport_option == c_Multi_Transport) {
600  doublereal sum, sumx, wtm, dz;
601  doublereal eps = 1.0e-12;
602  for (size_t m = j0; m < j1; m++) {
603  setGasAtMidpoint(x,m);
604  dz = m_z[m+1] - m_z[m];
605  wtm = m_thermo->meanMolecularWeight();
606 
607  m_visc[m] = (m_dovisc ? m_trans->viscosity() : 0.0);
608 
609  m_trans->getMultiDiffCoeffs(m_nsp,
610  DATA_PTR(m_multidiff) + mindex(0,0,m));
611 
612  for (size_t k = 0; k < m_nsp; k++) {
613  sum = 0.0;
614  sumx = 0.0;
615  for (size_t j = 0; j < m_nsp; j++) {
616  if (j != k) {
617  sum += m_wt[j]*m_multidiff[mindex(k,j,m)]*
618  ((X(x,j,m+1) - X(x,j,m))/dz + eps);
619  sumx += (X(x,j,m+1) - X(x,j,m))/dz;
620  }
621  }
622  m_diff[k + m*m_nsp] = sum/(wtm*(sumx+eps));
623  }
624 
625  m_tcon[m] = m_trans->thermalConductivity();
626  if (m_do_soret) {
627  m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + m*m_nsp);
628  }
629  }
630  }
631 }
632 
633 
634 //------------------------------------------------------
635 
636 /**
637  * Evaluate the residual function for axisymmetric stagnation
638  * flow. If jpt is less than zero, the residual function is
639  * evaluated at all grid points. If jpt >= 0, then the residual
640  * function is only evaluated at grid points jpt-1, jpt, and
641  * jpt+1. This option is used to efficiently evaluate the
642  * Jacobian numerically.
643  *
644  */
645 
646 void FreeFlame::eval(size_t jg, doublereal* xg,
647  doublereal* rg, integer* diagg, doublereal rdt)
648 {
649 
650  // if evaluating a Jacobian, and the global point is outside
651  // the domain of influence for this domain, then skip
652  // evaluating the residual
653  if (jg != npos && (jg + 1 < firstPoint() || jg > lastPoint() + 1)) {
654  return;
655  }
656 
657  // if evaluating a Jacobian, compute the steady-state residual
658  if (jg != npos) {
659  rdt = 0.0;
660  }
661 
662  // start of local part of global arrays
663  doublereal* x = xg + loc();
664  doublereal* rsd = rg + loc();
665  integer* diag = diagg + loc();
666 
667  size_t jmin, jmax;
668 
669  if (jg == npos) { // evaluate all points
670  jmin = 0;
671  jmax = m_points - 1;
672  } else { // evaluate points for Jacobian
673  size_t jpt = (jg == 0) ? 0 : jg - firstPoint();
674  jmin = std::max<size_t>(jpt, 1) - 1;
675  jmax = std::min(jpt+1,m_points-1);
676  }
677 
678  // properties are computed for grid points from j0 to j1
679  size_t j0 = std::max<size_t>(jmin, 1) - 1;
680  size_t j1 = std::min(jmax+1,m_points-1);
681 
682  size_t j, k;
683 
684  //-----------------------------------------------------
685  // update properties
686  //-----------------------------------------------------
687 
688  // update thermodynamic properties only if a Jacobian is not
689  // being evaluated
690  if (jg == npos) {
691  updateThermo(x, j0, j1);
692  updateTransport(x, j0, j1);
693  }
694 
695  // update the species diffusive mass fluxes whether or not a
696  // Jacobian is being evaluated
697  updateDiffFluxes(x, j0, j1);
698 
699 
700  //----------------------------------------------------
701  // evaluate the residual equations at all required
702  // grid points
703  //----------------------------------------------------
704 
705  doublereal sum, sum2, dtdzj;
706 
707  for (j = jmin; j <= jmax; j++) {
708 
709 
710  //----------------------------------------------
711  // left boundary
712  //----------------------------------------------
713 
714  if (j == 0) {
715 
716  // these may be modified by a boundary object
717 
718  // Continuity. This propagates information right-to-left,
719  // since rho_u at point 0 is dependent on rho_u at point 1,
720  // but not on mdot from the inlet.
721  rsd[index(c_offset_U,0)] =
722  -(rho_u(x,1) - rho_u(x,0))/m_dz[0]
723  -(density(1)*V(x,1) + density(0)*V(x,0));
724 
725  // the inlet (or other) object connected to this one
726  // will modify these equations by subtracting its values
727  // for V, T, and mdot. As a result, these residual equations
728  // will force the solution variables to the values for
729  // the boundary object
730  rsd[index(c_offset_V,0)] = V(x,0);
731  rsd[index(c_offset_T,0)] = T(x,0);
732  rsd[index(c_offset_L,0)] = -rho_u(x,0);
733 
734  // The default boundary condition for species is zero
735  // flux
736  sum = 0.0;
737  for (k = 0; k < m_nsp; k++) {
738  sum += Y(x,k,0);
739  rsd[index(c_offset_Y + k, 0)] =
740  -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
741  }
742  rsd[index(c_offset_Y, 0)] = 1.0 - sum;
743  }
744 
745 
746  //----------------------------------------------
747  //
748  // right boundary
749  //
750  //----------------------------------------------
751 
752  else if (j == m_points - 1) {
753 
754  // the boundary object connected to the right of this
755  // one may modify or replace these equations. The
756  // default boundary conditions are zero u, V, and T,
757  // and zero diffusive flux for all species.
758 
759  // zero gradient
760  rsd[index(0,j)] = rho_u(x,j) - rho_u(x,j-1);
761  rsd[index(1,j)] = V(x,j);
762  rsd[index(2,j)] = T(x,j) - T(x,j-1);
763  doublereal sum = 0.0;
764  rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
765  diag[index(c_offset_L, j)] = 0;
766  for (k = 0; k < m_nsp; k++) {
767  sum += Y(x,k,j);
768  rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
769  }
770  rsd[index(4,j)] = 1.0 - sum;
771  diag[index(4,j)] = 0;
772  }
773 
774  //------------------------------------------
775  // interior points
776  //------------------------------------------
777 
778  else {
779 
780  //----------------------------------------------
781  // Continuity equation
782  //----------------------------------------------
783 
784  if (grid(j) > m_zfixed) {
785  rsd[index(c_offset_U,j)] =
786  - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]
787  - (density(j-1)*V(x,j-1) + density(j)*V(x,j));
788  }
789 
790  else if (grid(j) == m_zfixed) {
791  if (m_do_energy[j]) {
792  rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
793  } else {
794  rsd[index(c_offset_U,j)] = (rho_u(x,j)
795  - m_rho[0]*0.3);
796  }
797  } else if (grid(j) < m_zfixed) {
798  rsd[index(c_offset_U,j)] =
799  - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
800  - (density(j+1)*V(x,j+1) + density(j)*V(x,j));
801  }
802  //algebraic constraint
803  diag[index(c_offset_U, j)] = 0;
804 
805  //------------------------------------------------
806  // Radial momentum equation
807  //
808  // \rho u dV/dz + \rho V^2 = d(\mu dV/dz)/dz - lambda
809  //
810  //-------------------------------------------------
811  rsd[index(c_offset_V,j)]
812  = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j)
813  - m_rho[j]*V(x,j)*V(x,j))/m_rho[j]
814  - rdt*(V(x,j) - V_prev(j));
815  diag[index(c_offset_V, j)] = 1;
816 
817 
818  //-------------------------------------------------
819  // Species equations
820  //
821  // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k
822  //
823  //-------------------------------------------------
824  getWdot(x,j);
825 
826  doublereal convec, diffus;
827  for (k = 0; k < m_nsp; k++) {
828  convec = rho_u(x,j)*dYdz(x,k,j);
829  diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))
830  /(z(j+1) - z(j-1));
831  rsd[index(c_offset_Y + k, j)]
832  = (m_wt[k]*(wdot(k,j))
833  - convec - diffus)/m_rho[j]
834  - rdt*(Y(x,k,j) - Y_prev(k,j));
835  diag[index(c_offset_Y + k, j)] = 1;
836  }
837 
838 
839  //-----------------------------------------------
840  // energy equation
841  //-----------------------------------------------
842 
843  if (m_do_energy[j]) {
844 
845  setGas(x,j);
846 
847  // heat release term
848  const vector_fp& h_RT = m_thermo->enthalpy_RT_ref();
849  const vector_fp& cp_R = m_thermo->cp_R_ref();
850 
851  sum = 0.0;
852  sum2 = 0.0;
853  doublereal flxk;
854  for (k = 0; k < m_nsp; k++) {
855  flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
856  sum += wdot(k,j)*h_RT[k];
857  sum2 += flxk*cp_R[k]/m_wt[k];
858  }
859  sum *= GasConstant * T(x,j);
860  dtdzj = dTdz(x,j);
861  sum2 *= GasConstant * dtdzj;
862 
863  rsd[index(c_offset_T, j)] =
864  - m_cp[j]*rho_u(x,j)*dtdzj
865  - divHeatFlux(x,j) - sum - sum2;
866  rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
867 
868  rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
869  diag[index(c_offset_T, j)] = 1;
870  }
871  // residual equations if the energy equation is disabled
872  else {
873  rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
874  diag[index(c_offset_T, j)] = 0;
875  }
876 
877  rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
878  diag[index(c_offset_L, j)] = 0;
879  }
880  }
881 }
882 
883 
884 /**
885  * Print the solution.
886  */
887 void StFlow::showSolution(const doublereal* x)
888 {
889  size_t nn = m_nv/5;
890  size_t i, j, n;
891  //char* buf = new char[100];
892  char buf[100];
893 
894  // The mean molecular weight is needed to convert
895  updateThermo(x, 0, m_points-1);
896 
897  sprintf(buf, " Pressure: %10.4g Pa \n", m_press);
898  writelog(buf);
899  for (i = 0; i < nn; i++) {
900  st_drawline();
901  sprintf(buf, "\n z ");
902  writelog(buf);
903  for (n = 0; n < 5; n++) {
904  sprintf(buf, " %10s ",componentName(i*5 + n).c_str());
905  writelog(buf);
906  }
907  st_drawline();
908  for (j = 0; j < m_points; j++) {
909  sprintf(buf, "\n %10.4g ",m_z[j]);
910  writelog(buf);
911  for (n = 0; n < 5; n++) {
912  sprintf(buf, " %10.4g ",component(x, i*5+n,j));
913  writelog(buf);
914  }
915  }
916  writelog("\n");
917  }
918  size_t nrem = m_nv - 5*nn;
919  st_drawline();
920  sprintf(buf, "\n z ");
921  writelog(buf);
922  for (n = 0; n < nrem; n++) {
923  sprintf(buf, " %10s ", componentName(nn*5 + n).c_str());
924  writelog(buf);
925  }
926  st_drawline();
927  for (j = 0; j < m_points; j++) {
928  sprintf(buf, "\n %10.4g ",m_z[j]);
929  writelog(buf);
930  for (n = 0; n < nrem; n++) {
931  sprintf(buf, " %10.4g ",component(x, nn*5+n,j));
932  writelog(buf);
933  }
934  }
935  writelog("\n");
936 }
937 
938 
939 /**
940  * Update the diffusive mass fluxes.
941  */
942 void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1)
943 {
944  size_t j, k, m;
945  doublereal sum, wtm, rho, dz, gradlogT;
946 
947  switch (m_transport_option) {
948 
949  case c_Mixav_Transport:
950  case c_Multi_Transport:
951  for (j = j0; j < j1; j++) {
952  sum = 0.0;
953  wtm = m_wtm[j];
954  rho = density(j);
955  dz = z(j+1) - z(j);
956 
957  for (k = 0; k < m_nsp; k++) {
958  m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
959  m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
960  sum -= m_flux(k,j);
961  }
962  // correction flux to insure that \sum_k Y_k V_k = 0.
963  for (k = 0; k < m_nsp; k++) {
964  m_flux(k,j) += sum*Y(x,k,j);
965  }
966  }
967  break;
968 
969  default:
970  throw CanteraError("updateDiffFluxes","unknown transport model");
971  }
972 
973  if (m_do_soret) {
974  for (m = j0; m < j1; m++) {
975  gradlogT = 2.0 * (T(x,m+1) - T(x,m)) /
976  ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m)));
977  for (k = 0; k < m_nsp; k++) {
978  m_flux(k,m) -= m_dthermal(k,m)*gradlogT;
979  }
980  }
981  }
982 }
983 
984 
985 string StFlow::componentName(size_t n) const
986 {
987  switch (n) {
988  case 0:
989  return "u";
990  case 1:
991  return "V";
992  case 2:
993  return "T";
994  case 3:
995  return "lambda";
996  default:
997  if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) {
998  return m_thermo->speciesName(n - c_offset_Y);
999  } else {
1000  return "<unknown>";
1001  }
1002  }
1003 }
1004 
1005 
1006 size_t StFlow::componentIndex(string name) const
1007 {
1008 
1009 
1010  if (name=="u") {
1011  return 0;
1012  } else if (name=="V") {
1013  return 1;
1014  } else if (name=="T") {
1015  return 2;
1016  } else if (name=="lambda") {
1017  return 3;
1018  } else {
1019  for (size_t n=4; n<m_nsp+4; n++) {
1020  if (componentName(n)==name) {
1021  return n;
1022  }
1023  }
1024  }
1025 
1026  return npos;
1027 }
1028 
1029 
1030 void StFlow::restore(const XML_Node& dom, doublereal* soln)
1031 {
1032 
1033  vector<string> ignored;
1034  size_t nsp = m_thermo->nSpecies();
1035  vector_int did_species(nsp, 0);
1036 
1037  vector<XML_Node*> str;
1038  dom.getChildren("string",str);
1039  int nstr = static_cast<int>(str.size());
1040  for (int istr = 0; istr < nstr; istr++) {
1041  const XML_Node& nd = *str[istr];
1042  writelog(nd["title"]+": "+nd.value()+"\n");
1043  }
1044 
1045  //map<string, double> params;
1046  double pp = -1.0;
1047  pp = getFloat(dom, "pressure", "pressure");
1048  setPressure(pp);
1049 
1050 
1051  vector<XML_Node*> d;
1052  dom.child("grid_data").getChildren("floatArray",d);
1053  size_t nd = d.size();
1054 
1055  vector_fp x;
1056  size_t n, np = 0, j, ks, k;
1057  string nm;
1058  bool readgrid = false, wrote_header = false;
1059  for (n = 0; n < nd; n++) {
1060  const XML_Node& fa = *d[n];
1061  nm = fa["title"];
1062  if (nm == "z") {
1063  getFloatArray(fa,x,false);
1064  np = x.size();
1065  writelog("Grid contains "+int2str(np)+
1066  " points.\n");
1067  readgrid = true;
1068  setupGrid(np, DATA_PTR(x));
1069  }
1070  }
1071  if (!readgrid) {
1072  throw CanteraError("StFlow::restore",
1073  "domain contains no grid points.");
1074  }
1075 
1076  writelog("Importing datasets:\n");
1077  for (n = 0; n < nd; n++) {
1078  const XML_Node& fa = *d[n];
1079  nm = fa["title"];
1080  getFloatArray(fa,x,false);
1081  if (nm == "u") {
1082  writelog("axial velocity ");
1083  if (x.size() == np) {
1084  for (j = 0; j < np; j++) {
1085  soln[index(0,j)] = x[j];
1086  }
1087  } else {
1088  goto error;
1089  }
1090  } else if (nm == "z") {
1091  ; // already read grid
1092  } else if (nm == "V") {
1093  writelog("radial velocity ");
1094  if (x.size() == np) {
1095  for (j = 0; j < np; j++) {
1096  soln[index(1,j)] = x[j];
1097  }
1098  } else {
1099  goto error;
1100  }
1101  } else if (nm == "T") {
1102  writelog("temperature ");
1103  if (x.size() == np) {
1104  for (j = 0; j < np; j++) {
1105  soln[index(2,j)] = x[j];
1106  }
1107 
1108  // For fixed-temperature simulations, use the
1109  // imported temperature profile by default. If
1110  // this is not desired, call setFixedTempProfile
1111  // *after* restoring the solution.
1112 
1113  vector_fp zz(np);
1114  for (size_t jj = 0; jj < np; jj++) {
1115  zz[jj] = (grid(jj) - zmin())/(zmax() - zmin());
1116  }
1117  setFixedTempProfile(zz, x);
1118  } else {
1119  goto error;
1120  }
1121  } else if (nm == "L") {
1122  writelog("lambda ");
1123  if (x.size() == np) {
1124  for (j = 0; j < np; j++) {
1125  soln[index(3,j)] = x[j];
1126  }
1127  } else {
1128  goto error;
1129  }
1130  } else if (m_thermo->speciesIndex(nm) != npos) {
1131  writelog(nm+" ");
1132  if (x.size() == np) {
1133  k = m_thermo->speciesIndex(nm);
1134  did_species[k] = 1;
1135  for (j = 0; j < np; j++) {
1136  soln[index(k+4,j)] = x[j];
1137  }
1138  }
1139  } else {
1140  ignored.push_back(nm);
1141  }
1142  }
1143 
1144  if (ignored.size() != 0) {
1145  writelog("\n\n");
1146  writelog("Ignoring datasets:\n");
1147  size_t nn = ignored.size();
1148  for (size_t n = 0; n < nn; n++) {
1149  writelog(ignored[n]+" ");
1150  }
1151  }
1152 
1153  for (ks = 0; ks < nsp; ks++) {
1154  if (did_species[ks] == 0) {
1155  if (!wrote_header) {
1156  writelog("Missing data for species:\n");
1157  wrote_header = true;
1158  }
1159  writelog(m_thermo->speciesName(ks)+" ");
1160  }
1161  }
1162 
1163  return;
1164 error:
1165  throw CanteraError("StFlow::restore","Data size error");
1166 }
1167 
1168 
1169 
1170 void StFlow::save(XML_Node& o, const doublereal* const sol)
1171 {
1172  size_t k;
1173 
1174  Array2D soln(m_nv, m_points, sol + loc());
1175 
1176  XML_Node& flow = (XML_Node&)o.addChild("domain");
1177  flow.addAttribute("type",flowType());
1178  flow.addAttribute("id",m_id);
1179  flow.addAttribute("points", double(m_points));
1180  flow.addAttribute("components", double(m_nv));
1181 
1182  if (m_desc != "") {
1183  addString(flow,"description",m_desc);
1184  }
1185  XML_Node& gv = flow.addChild("grid_data");
1186  addFloat(flow, "pressure", m_press, "Pa", "pressure");
1187  addFloatArray(gv,"z",m_z.size(),DATA_PTR(m_z),
1188  "m","length");
1189  vector_fp x(soln.nColumns());
1190 
1191  soln.getRow(0,DATA_PTR(x));
1192  addFloatArray(gv,"u",x.size(),DATA_PTR(x),"m/s","velocity");
1193 
1194  soln.getRow(1,DATA_PTR(x));
1195  addFloatArray(gv,"V",
1196  x.size(),DATA_PTR(x),"1/s","rate");
1197 
1198  soln.getRow(2,DATA_PTR(x));
1199  addFloatArray(gv,"T",x.size(),DATA_PTR(x),"K","temperature",0.0);
1200 
1201  soln.getRow(3,DATA_PTR(x));
1202  addFloatArray(gv,"L",x.size(),DATA_PTR(x),"N/m^4");
1203 
1204  for (k = 0; k < m_nsp; k++) {
1205  soln.getRow(4+k,DATA_PTR(x));
1206  addFloatArray(gv,m_thermo->speciesName(k),
1207  x.size(),DATA_PTR(x),"","massFraction",0.0,1.0);
1208  }
1209 }
1210 
1211 
1212 void StFlow::setJac(MultiJac* jac)
1213 {
1214  m_jac = jac;
1215 }
1216 
1217 
1218 } // namespace