Cantera  2.4.0
boundaries1D.cpp
Go to the documentation of this file.
1 //! @file boundaries1D.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at http://www.cantera.org/license.txt for license and copyright information.
5 
6 #include "cantera/oneD/Inlet1D.h"
7 #include "cantera/oneD/OneDim.h"
8 #include "cantera/base/ctml.h"
9 #include "cantera/oneD/StFlow.h"
10 
11 using namespace std;
12 
13 namespace Cantera
14 {
15 
16 Bdry1D::Bdry1D() : Domain1D(1, 1, 0.0),
17  m_flow_left(0), m_flow_right(0),
18  m_ilr(0), m_left_nv(0), m_right_nv(0),
19  m_left_loc(0), m_right_loc(0),
20  m_left_points(0),
21  m_left_nsp(0), m_right_nsp(0),
22  m_sp_left(0), m_sp_right(0),
23  m_start_left(0), m_start_right(0),
24  m_phase_left(0), m_phase_right(0), m_temp(0.0), m_mdot(0.0)
25 {
26  m_type = cConnectorType;
27 }
28 
29 void Bdry1D::_init(size_t n)
30 {
31  if (m_index == npos) {
32  throw CanteraError("Bdry1D::_init",
33  "install in container before calling init.");
34  }
35 
36  // A boundary object contains only one grid point
37  resize(n,1);
38 
39  m_left_nsp = 0;
40  m_right_nsp = 0;
41 
42  // check for left and right flow objects
43  if (m_index > 0) {
44  Domain1D& r = container().domain(m_index-1);
45  if (!r.isConnector()) { // flow domain
46  m_flow_left = (StFlow*)&r;
47  m_left_nv = m_flow_left->nComponents();
48  m_left_points = m_flow_left->nPoints();
49  m_left_loc = container().start(m_index-1);
50  m_left_nsp = m_left_nv - c_offset_Y;
51  m_phase_left = &m_flow_left->phase();
52  } else {
53  throw CanteraError("Bdry1D::_init",
54  "Boundary domains can only be connected on the left to flow "
55  "domains, not type {} domains.", r.domainType());
56  }
57  }
58 
59  // if this is not the last domain, see what is connected on the right
60  if (m_index + 1 < container().nDomains()) {
61  Domain1D& r = container().domain(m_index+1);
62  if (!r.isConnector()) { // flow domain
63  m_flow_right = (StFlow*)&r;
64  m_right_nv = m_flow_right->nComponents();
65  m_right_loc = container().start(m_index+1);
66  m_right_nsp = m_right_nv - c_offset_Y;
67  m_phase_right = &m_flow_right->phase();
68  } else {
69  throw CanteraError("Bdry1D::_init",
70  "Boundary domains can only be connected on the right to flow "
71  "domains, not type {} domains.", r.domainType());
72  }
73  }
74 }
75 
76 // ---------------- Inlet1D methods ----------------
77 
78 Inlet1D::Inlet1D()
79  : m_V0(0.0)
80  , m_nsp(0)
81  , m_flow(0)
82 {
83  m_type = cInletType;
84  m_xstr = "";
85 }
86 
87 void Inlet1D::showSolution(const double* x)
88 {
89  writelog(" Mass Flux: {:10.4g} kg/m^2/s \n", m_mdot);
90  writelog(" Temperature: {:10.4g} K \n", m_temp);
91  if (m_flow) {
92  writelog(" Mass Fractions: \n");
93  for (size_t k = 0; k < m_flow->phase().nSpecies(); k++) {
94  if (m_yin[k] != 0.0) {
95  writelog(" {:>16s} {:10.4g} \n",
96  m_flow->phase().speciesName(k), m_yin[k]);
97  }
98  }
99  }
100  writelog("\n");
101 }
102 
103 void Inlet1D::setMoleFractions(const std::string& xin)
104 {
105  m_xstr = xin;
106  if (m_flow) {
107  m_flow->phase().setMoleFractionsByName(xin);
108  m_flow->phase().getMassFractions(m_yin.data());
109  needJacUpdate();
110  }
111 }
112 
113 void Inlet1D::setMoleFractions(const doublereal* xin)
114 {
115  if (m_flow) {
116  m_flow->phase().setMoleFractions(xin);
117  m_flow->phase().getMassFractions(m_yin.data());
118  needJacUpdate();
119  }
120 }
121 
122 void Inlet1D::init()
123 {
124  _init(0);
125 
126  // if a flow domain is present on the left, then this must be a right inlet.
127  // Note that an inlet object can only be a terminal object - it cannot have
128  // flows on both the left and right
129  if (m_flow_left) {
130  m_ilr = RightInlet;
131  m_flow = m_flow_left;
132  } else if (m_flow_right) {
133  m_ilr = LeftInlet;
134  m_flow = m_flow_right;
135  } else {
136  throw CanteraError("Inlet1D::init","no flow!");
137  }
138 
139  // components = u, V, T, lambda, + mass fractions
140  m_nsp = m_flow->nComponents() - c_offset_Y;
141  m_yin.resize(m_nsp, 0.0);
142  if (m_xstr != "") {
143  setMoleFractions(m_xstr);
144  } else {
145  m_yin[0] = 1.0;
146  }
147 }
148 
149 void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg,
150  integer* diagg, doublereal rdt)
151 {
152  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
153  return;
154  }
155 
156  if (m_ilr == LeftInlet) {
157  // Array elements corresponding to the first point of the flow domain
158  double* xb = xg + m_flow->loc();
159  double* rb = rg + m_flow->loc();
160 
161  // The first flow residual is for u. This, however, is not modified by
162  // the inlet, since this is set within the flow domain from the
163  // continuity equation.
164 
165  // spreading rate. The flow domain sets this to V(0),
166  // so for finite spreading rate subtract m_V0.
167  rb[c_offset_V] -= m_V0;
168 
169  if (m_flow->doEnergy(0)) {
170  // The third flow residual is for T, where it is set to T(0). Subtract
171  // the local temperature to hold the flow T to the inlet T.
172  rb[c_offset_T] -= m_temp;
173  }
174 
175  if (m_flow->fixed_mdot()) {
176  // The flow domain sets this to -rho*u. Add mdot to specify the mass
177  // flow rate.
178  rb[c_offset_L] += m_mdot;
179  } else {
180  // if the flow is a freely-propagating flame, mdot is not specified.
181  // Set mdot equal to rho*u, and also set lambda to zero.
182  m_mdot = m_flow->density(0)*xb[0];
183  rb[c_offset_L] = xb[c_offset_L];
184  }
185 
186  // add the convective term to the species residual equations
187  for (size_t k = 0; k < m_nsp; k++) {
188  if (k != m_flow_right->leftExcessSpecies()) {
189  rb[c_offset_Y+k] += m_mdot*m_yin[k];
190  }
191  }
192 
193  } else {
194  // right inlet
195  // Array elements corresponding to the flast point in the flow domain
196  double* rb = rg + loc() - m_flow->nComponents();
197  rb[c_offset_V] -= m_V0;
198  if (m_flow->doEnergy(m_flow->nPoints() - 1)) {
199  rb[c_offset_T] -= m_temp; // T
200  }
201  rb[c_offset_U] += m_mdot; // u
202  for (size_t k = 0; k < m_nsp; k++) {
203  if (k != m_flow_left->rightExcessSpecies()) {
204  rb[c_offset_Y+k] += m_mdot * m_yin[k];
205  }
206  }
207  }
208 }
209 
210 XML_Node& Inlet1D::save(XML_Node& o, const doublereal* const soln)
211 {
212  XML_Node& inlt = Domain1D::save(o, soln);
213  inlt.addAttribute("type","inlet");
214  addFloat(inlt, "temperature", m_temp);
215  addFloat(inlt, "mdot", m_mdot);
216  for (size_t k=0; k < m_nsp; k++) {
217  addFloat(inlt, "massFraction", m_yin[k], "",
218  m_flow->phase().speciesName(k));
219  }
220  return inlt;
221 }
222 
223 void Inlet1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
224 {
225  Domain1D::restore(dom, soln, loglevel);
226  m_mdot = getFloat(dom, "mdot");
227  m_temp = getFloat(dom, "temperature");
228 
229  m_yin.assign(m_nsp, 0.0);
230 
231  for (size_t i = 0; i < dom.nChildren(); i++) {
232  const XML_Node& node = dom.child(i);
233  if (node.name() == "massFraction") {
234  size_t k = m_flow->phase().speciesIndex(node.attrib("type"));
235  if (k != npos) {
236  m_yin[k] = node.fp_value();
237  }
238  }
239  }
240  resize(0, 1);
241 }
242 
243 // ------------- Empty1D -------------
244 
245 void Empty1D::init()
246 {
247  _init(0);
248 }
249 
250 void Empty1D::eval(size_t jg, doublereal* xg, doublereal* rg,
251  integer* diagg, doublereal rdt)
252 {
253 }
254 
255 XML_Node& Empty1D::save(XML_Node& o, const doublereal* const soln)
256 {
257  XML_Node& symm = Domain1D::save(o, soln);
258  symm.addAttribute("type","empty");
259  return symm;
260 }
261 
262 void Empty1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
263 {
264  Domain1D::restore(dom, soln, loglevel);
265  resize(0, 1);
266 }
267 
268 // -------------- Symm1D --------------
269 
270 void Symm1D::init()
271 {
272  _init(0);
273 }
274 
275 void Symm1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg,
276  doublereal rdt)
277 {
278  if (jg != npos && (jg + 2< firstPoint() || jg > lastPoint() + 2)) {
279  return;
280  }
281 
282  // start of local part of global arrays
283  doublereal* x = xg + loc();
284  doublereal* r = rg + loc();
285  integer* diag = diagg + loc();
286 
287  if (m_flow_right) {
288  size_t nc = m_flow_right->nComponents();
289  double* xb = x;
290  double* rb = r;
291  int* db = diag;
292  db[c_offset_V] = 0;
293  db[c_offset_T] = 0;
294  rb[c_offset_V] = xb[c_offset_V] - xb[c_offset_V + nc]; // zero dV/dz
295  if (m_flow_right->doEnergy(0)) {
296  rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T + nc]; // zero dT/dz
297  }
298  }
299 
300  if (m_flow_left) {
301  size_t nc = m_flow_left->nComponents();
302  double* xb = x - nc;
303  double* rb = r - nc;
304  int* db = diag - nc;
305  db[c_offset_V] = 0;
306  db[c_offset_T] = 0;
307  rb[c_offset_V] = xb[c_offset_V] - xb[c_offset_V - nc]; // zero dV/dz
308  if (m_flow_left->doEnergy(m_flow_left->nPoints() - 1)) {
309  rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T - nc]; // zero dT/dz
310  }
311  }
312 }
313 
314 XML_Node& Symm1D::save(XML_Node& o, const doublereal* const soln)
315 {
316  XML_Node& symm = Domain1D::save(o, soln);
317  symm.addAttribute("type","symmetry");
318  return symm;
319 }
320 
321 void Symm1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
322 {
323  Domain1D::restore(dom, soln, loglevel);
324  resize(0, 1);
325 }
326 
327 // -------- Outlet1D --------
328 
329 OutletRes1D::OutletRes1D()
330  : m_nsp(0)
331  , m_flow(0)
332 {
333  m_type = cOutletResType;
334  m_xstr = "";
335 }
336 
338 {
339  _init(0);
340 
341  if (m_flow_right) {
342  m_flow_right->setViscosityFlag(false);
343  }
344  if (m_flow_left) {
345  m_flow_left->setViscosityFlag(false);
346  }
347 }
348 
349 void Outlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg,
350  doublereal rdt)
351 {
352  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
353  return;
354  }
355 
356  // start of local part of global arrays
357  doublereal* x = xg + loc();
358  doublereal* r = rg + loc();
359  integer* diag = diagg + loc();
360 
361  if (m_flow_right) {
362  size_t nc = m_flow_right->nComponents();
363  double* xb = x;
364  double* rb = r;
365  rb[c_offset_U] = xb[c_offset_L];
366  if (m_flow_right->doEnergy(0)) {
367  rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T + nc];
368  }
369  for (size_t k = c_offset_Y; k < nc; k++) {
370  rb[k] = xb[k] - xb[k + nc];
371  }
372  }
373 
374  if (m_flow_left) {
375  size_t nc = m_flow_left->nComponents();
376  double* xb = x - nc;
377  double* rb = r - nc;
378  int* db = diag - nc;
379 
380  // zero Lambda
381  if (m_flow_left->fixed_mdot()) {
382  rb[c_offset_U] = xb[c_offset_L];
383  }
384 
385  if (m_flow_left->doEnergy(m_flow_left->nPoints()-1)) {
386  rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T - nc]; // zero T gradient
387  }
388  size_t kSkip = c_offset_Y + m_flow_left->rightExcessSpecies();
389  for (size_t k = c_offset_Y; k < nc; k++) {
390  if (k != kSkip) {
391  rb[k] = xb[k] - xb[k - nc]; // zero mass fraction gradient
392  db[k] = 0;
393  }
394  }
395  }
396 }
397 
398 XML_Node& Outlet1D::save(XML_Node& o, const doublereal* const soln)
399 {
400  XML_Node& outlt = Domain1D::save(o, soln);
401  outlt.addAttribute("type","outlet");
402  return outlt;
403 }
404 
405 void Outlet1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
406 {
407  Domain1D::restore(dom, soln, loglevel);
408  resize(0, 1);
409 }
410 
411 // -------- OutletRes1D --------
412 
413 void OutletRes1D::setMoleFractions(const std::string& xres)
414 {
415  m_xstr = xres;
416  if (m_flow) {
417  m_flow->phase().setMoleFractionsByName(xres);
418  m_flow->phase().getMassFractions(m_yres.data());
419  needJacUpdate();
420  }
421 }
422 
423 void OutletRes1D::setMoleFractions(const doublereal* xres)
424 {
425  if (m_flow) {
426  m_flow->phase().setMoleFractions(xres);
427  m_flow->phase().getMassFractions(m_yres.data());
428  needJacUpdate();
429  }
430 }
431 
433 {
434  _init(0);
435 
436  if (m_flow_left) {
437  m_flow = m_flow_left;
438  } else if (m_flow_right) {
439  m_flow = m_flow_right;
440  } else {
441  throw CanteraError("OutletRes1D::init","no flow!");
442  }
443 
444  m_nsp = m_flow->nComponents() - c_offset_Y;
445  m_yres.resize(m_nsp, 0.0);
446  if (m_xstr != "") {
447  setMoleFractions(m_xstr);
448  } else {
449  m_yres[0] = 1.0;
450  }
451 }
452 
453 void OutletRes1D::eval(size_t jg, doublereal* xg, doublereal* rg,
454  integer* diagg, doublereal rdt)
455 {
456  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
457  return;
458  }
459 
460  // start of local part of global arrays
461  doublereal* x = xg + loc();
462  doublereal* r = rg + loc();
463  integer* diag = diagg + loc();
464 
465  if (m_flow_right) {
466  size_t nc = m_flow_right->nComponents();
467  double* xb = x;
468  double* rb = r;
469 
470  // this seems wrong...
471  // zero Lambda
472  rb[c_offset_U] = xb[c_offset_L];
473 
474  if (m_flow_right->doEnergy(0)) {
475  // zero gradient for T
476  rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T + nc];
477  }
478 
479  // specified mass fractions
480  for (size_t k = c_offset_Y; k < nc; k++) {
481  rb[k] = xb[k] - m_yres[k-c_offset_Y];
482  }
483  }
484 
485  if (m_flow_left) {
486  size_t nc = m_flow_left->nComponents();
487  double* xb = x - nc;
488  double* rb = r - nc;
489  int* db = diag - nc;
490 
491  if (!m_flow_left->fixed_mdot()) {
492  ;
493  } else {
494  rb[c_offset_U] = xb[c_offset_L]; // zero Lambda
495  }
496  if (m_flow_left->doEnergy(m_flow_left->nPoints()-1)) {
497  rb[c_offset_T] = xb[c_offset_T] - m_temp; // zero dT/dz
498  }
499  size_t kSkip = m_flow_left->rightExcessSpecies();
500  for (size_t k = c_offset_Y; k < nc; k++) {
501  if (k != kSkip) {
502  rb[k] = xb[k] - m_yres[k-c_offset_Y]; // fixed Y
503  db[k] = 0;
504  }
505  }
506  }
507 }
508 
509 XML_Node& OutletRes1D::save(XML_Node& o, const doublereal* const soln)
510 {
511  XML_Node& outlt = Domain1D::save(o, soln);
512  outlt.addAttribute("type","outletres");
513  addFloat(outlt, "temperature", m_temp, "K");
514  for (size_t k=0; k < m_nsp; k++) {
515  addFloat(outlt, "massFraction", m_yres[k], "",
516  m_flow->phase().speciesName(k));
517  }
518  return outlt;
519 }
520 
521 void OutletRes1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
522 {
523  Domain1D::restore(dom, soln, loglevel);
524  m_temp = getFloat(dom, "temperature");
525 
526  m_yres.assign(m_nsp, 0.0);
527  for (size_t i = 0; i < dom.nChildren(); i++) {
528  const XML_Node& node = dom.child(i);
529  if (node.name() == "massFraction") {
530  size_t k = m_flow->phase().speciesIndex(node.attrib("type"));
531  if (k != npos) {
532  m_yres[k] = node.fp_value();
533  }
534  }
535  }
536 
537  resize(0, 1);
538 }
539 
540 // -------- Surf1D --------
541 
543 {
544  _init(0);
545 }
546 
547 void Surf1D::eval(size_t jg, doublereal* xg, doublereal* rg,
548  integer* diagg, doublereal rdt)
549 {
550  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
551  return;
552  }
553 
554  // start of local part of global arrays
555  doublereal* x = xg + loc();
556  doublereal* r = rg + loc();
557 
558  if (m_flow_right) {
559  double* rb = r;
560  double* xb = x;
561  rb[c_offset_T] = xb[c_offset_T] - m_temp; // specified T
562  }
563 
564  if (m_flow_left) {
565  size_t nc = m_flow_left->nComponents();
566  double* rb = r - nc;
567  double* xb = x - nc;
568  rb[c_offset_T] = xb[c_offset_T] - m_temp; // specified T
569  }
570 }
571 
572 XML_Node& Surf1D::save(XML_Node& o, const doublereal* const soln)
573 {
574  XML_Node& inlt = Domain1D::save(o, soln);
575  inlt.addAttribute("type","surface");
576  addFloat(inlt, "temperature", m_temp);
577  return inlt;
578 }
579 
580 void Surf1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
581 {
582  Domain1D::restore(dom, soln, loglevel);
583  m_temp = getFloat(dom, "temperature");
584  resize(0, 1);
585 }
586 
587 void Surf1D::showSolution_s(std::ostream& s, const double* x)
588 {
589  s << "------------------- Surface " << domainIndex() << " ------------------- " << std::endl;
590  s << " temperature: " << m_temp << " K" << std::endl;
591 }
592 
593 // -------- ReactingSurf1D --------
594 
595 ReactingSurf1D::ReactingSurf1D()
596  : m_kin(0)
597  , m_surfindex(0)
598  , m_nsp(0)
599 {
600  m_type = cSurfType;
601 }
602 
603 void ReactingSurf1D::setKineticsMgr(InterfaceKinetics* kin)
604 {
605  m_kin = kin;
606  m_surfindex = kin->surfacePhaseIndex();
607  m_sphase = (SurfPhase*)&kin->thermo(m_surfindex);
608  m_nsp = m_sphase->nSpecies();
609  m_enabled = true;
610 }
611 
612 string ReactingSurf1D::componentName(size_t n) const
613 {
614  if (n < m_nsp) {
615  return m_sphase->speciesName(n);
616  } else {
617  return "<unknown>";
618  }
619 }
620 
622 {
623  m_nv = m_nsp;
624  _init(m_nsp);
625  m_fixed_cov.resize(m_nsp, 0.0);
626  m_fixed_cov[0] = 1.0;
627  m_work.resize(m_kin->nTotalSpecies(), 0.0);
628 
629  for (size_t n = 0; n < m_nsp; n++) {
630  setBounds(n, -1.0e-5, 2.0);
631  }
632 }
633 
635  double* x = xg + loc();
636  m_sphase->setCoverages(x);
637  m_sphase->getCoverages(x);
638 }
639 
640 void ReactingSurf1D::eval(size_t jg, doublereal* xg, doublereal* rg,
641  integer* diagg, doublereal rdt)
642 {
643  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
644  return;
645  }
646 
647  // start of local part of global arrays
648  doublereal* x = xg + loc();
649  doublereal* r = rg + loc();
650  integer* diag = diagg + loc();
651 
652  // set the coverages
653  doublereal sum = 0.0;
654  for (size_t k = 0; k < m_nsp; k++) {
655  m_work[k] = x[k];
656  sum += x[k];
657  }
658  m_sphase->setTemperature(m_temp);
659  m_sphase->setCoveragesNoNorm(m_work.data());
660 
661  // set the left gas state to the adjacent point
662 
663  size_t leftloc = 0, rightloc = 0;
664  size_t pnt = 0;
665 
666  if (m_flow_left) {
667  leftloc = m_flow_left->loc();
668  pnt = m_flow_left->nPoints() - 1;
669  m_flow_left->setGas(xg + leftloc, pnt);
670  }
671 
672  if (m_flow_right) {
673  rightloc = m_flow_right->loc();
674  m_flow_right->setGas(xg + rightloc, 0);
675  }
676 
677  m_kin->getNetProductionRates(m_work.data());
678  doublereal rs0 = 1.0/m_sphase->siteDensity();
679  size_t ioffset = m_kin->kineticsSpeciesIndex(0, m_surfindex);
680 
681  if (m_enabled) {
682  doublereal maxx = -1.0;
683  for (size_t k = 0; k < m_nsp; k++) {
684  r[k] = m_work[k + ioffset] * m_sphase->size(k) * rs0;
685  r[k] -= rdt*(x[k] - prevSoln(k,0));
686  diag[k] = 1;
687  maxx = std::max(x[k], maxx);
688  }
689  r[0] = 1.0 - sum;
690  diag[0] = 0;
691  } else {
692  for (size_t k = 0; k < m_nsp; k++) {
693  r[k] = x[k] - m_fixed_cov[k];
694  diag[k] = 0;
695  }
696  }
697 
698  if (m_flow_right) {
699  double* rb = r + m_nsp;
700  double* xb = x + m_nsp;
701  rb[c_offset_T] = xb[c_offset_T] - m_temp; // specified T
702  }
703  if (m_flow_left) {
704  size_t nc = m_flow_left->nComponents();
705  const vector_fp& mwleft = m_phase_left->molecularWeights();
706  double* rb = r - nc;
707  double* xb = x - nc;
708  rb[c_offset_T] = xb[c_offset_T] - m_temp; // specified T
709  size_t nSkip = m_flow_left->rightExcessSpecies();
710  for (size_t nl = 0; nl < m_left_nsp; nl++) {
711  if (nl != nSkip) {
712  rb[c_offset_Y+nl] += m_work[nl]*mwleft[nl];
713  }
714  }
715  }
716 }
717 
718 XML_Node& ReactingSurf1D::save(XML_Node& o, const doublereal* const soln)
719 {
720  const doublereal* s = soln + loc();
721  XML_Node& dom = Domain1D::save(o, soln);
722  dom.addAttribute("type","surface");
723  addFloat(dom, "temperature", m_temp, "K");
724  for (size_t k=0; k < m_nsp; k++) {
725  addFloat(dom, "coverage", s[k], "", m_sphase->speciesName(k));
726  }
727  return dom;
728 }
729 
730 void ReactingSurf1D::restore(const XML_Node& dom, doublereal* soln,
731  int loglevel)
732 {
733  Domain1D::restore(dom, soln, loglevel);
734  m_temp = getFloat(dom, "temperature");
735 
736  m_fixed_cov.assign(m_nsp, 0.0);
737  for (size_t i = 0; i < dom.nChildren(); i++) {
738  const XML_Node& node = dom.child(i);
739  if (node.name() == "coverage") {
740  size_t k = m_sphase->speciesIndex(node.attrib("type"));
741  if (k != npos) {
742  m_fixed_cov[k] = soln[k] = node.fp_value();
743  }
744  }
745  }
746  m_sphase->setCoverages(&m_fixed_cov[0]);
747 
748  resize(m_nsp, 1);
749 }
750 
751 void ReactingSurf1D::showSolution(const double* x)
752 {
753  writelog(" Temperature: {:10.4g} K \n", m_temp);
754  writelog(" Coverages: \n");
755  for (size_t k = 0; k < m_nsp; k++) {
756  writelog(" {:>20s} {:10.4g} \n", m_sphase->speciesName(k), x[k]);
757  }
758  writelog("\n");
759 }
760 }
void setCoveragesNoNorm(const doublereal *theta)
Set the surface site fractions to a specified state.
Definition: SurfPhase.cpp:253
virtual double size(size_t k) const
Returns the number of sites occupied by one molecule of species k.
Definition: SurfPhase.h:317
const vector_fp & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:437
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
std::string name() const
Returns the name of the XML node.
Definition: xml.h:370
void getMassFractions(doublereal *const y) const
Get the species mass fractions.
Definition: Phase.cpp:508
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
size_t speciesIndex(const std::string &name) const
Returns the index of a species named &#39;name&#39; within the Phase object.
Definition: Phase.cpp:175
virtual void setMoleFractions(const std::string &xin)
Set the mole fractions by specifying a std::string.
virtual void eval(size_t jg, doublereal *xg, doublereal *rg, integer *diagg, doublereal rdt)
Evaluate the residual function at point j.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
size_t firstPoint() const
The index of the first (i.e., left-most) grid point belonging to this domain.
Definition: Domain1D.h:361
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:153
void setCoverages(const doublereal *theta)
Set the surface site fractions to a specified state.
Definition: SurfPhase.cpp:236
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:266
STL namespace.
virtual void eval(size_t jg, doublereal *xg, doublereal *rg, integer *diagg, doublereal rdt)
Evaluate the residual function at point j.
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
Definition: Domain1D.cpp:123
virtual void eval(size_t jg, doublereal *xg, doublereal *rg, integer *diagg, doublereal rdt)
Evaluate the residual function at point j.
virtual XML_Node & save(XML_Node &o, const doublereal *const soln)
Save the current solution for this domain into an XML_Node.
virtual void init()
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:110
virtual void getNetProductionRates(doublereal *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition: Kinetics.cpp:418
virtual XML_Node & save(XML_Node &o, const doublereal *const soln)
Save the current solution for this domain into an XML_Node.
size_t lastPoint() const
The index of the last (i.e., right-most) grid point belonging to this domain.
Definition: Domain1D.h:369
virtual void eval(size_t jg, doublereal *xg, doublereal *rg, integer *diagg, doublereal rdt)
Evaluate the residual function at point j.
Boundary objects for one-dimensional simulations.
virtual void resize(size_t nv, size_t np)
Definition: Domain1D.cpp:33
virtual void showSolution(const doublereal *x)
Print the solution.
void getCoverages(doublereal *theta) const
Return a vector of surface coverages.
Definition: SurfPhase.cpp:261
virtual XML_Node & save(XML_Node &o, const doublereal *const soln)
Save the current solution for this domain into an XML_Node.
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition: Kinetics.h:239
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:191
void setMoleFractionsByName(const compositionMap &xMap)
Set the species mole fractions by name.
Definition: Phase.cpp:292
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:261
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values.
Definition: Phase.cpp:251
virtual XML_Node & save(XML_Node &o, const doublereal *const soln)
Save the current solution for this domain into an XML_Node.
void addAttribute(const std::string &attrib, const std::string &value)
Add or modify an attribute of the current node.
Definition: xml.cpp:474
size_t domainIndex()
The left-to-right location of this domain.
Definition: Domain1D.h:59
double prevSoln(size_t n, size_t j) const
Value of component n at point j in the previous solution.
Definition: Domain1D.h:404
XML_Node & child(const size_t n) const
Return a changeable reference to the n&#39;th child of the current node.
Definition: xml.cpp:546
virtual void init()
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:166
void addFloat(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:19
size_t nComponents() const
Number of components at each grid point.
Definition: Domain1D.h:134
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:500
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:157
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:637
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
doublereal getFloat(const XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
Definition: ctml.cpp:164
size_t rightExcessSpecies() const
Index of the species on the right boundary with the largest mass fraction.
Definition: StFlow.h:258
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
void needJacUpdate()
Definition: Domain1D.cpp:102
size_t nPoints() const
Number of grid points in this domain.
Definition: Domain1D.h:156
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
doublereal fp_value() const
Return the value of an XML node as a single double.
Definition: xml.cpp:454
size_t nChildren(bool discardComments=false) const
Return the number of children.
Definition: xml.cpp:556
virtual void init()
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:353
doublereal siteDensity()
Returns the site density.
Definition: SurfPhase.h:312
virtual void resetBadValues(double *xg)