Cantera  2.1.2
boundaries1D.cpp
Go to the documentation of this file.
1 /**
2  * @file boundaries1D.cpp
3  */
4 // Copyright 2002-3 California Institute of Technology
5 
6 #include "cantera/oneD/Inlet1D.h"
7 
8 using namespace std;
9 
10 namespace Cantera
11 {
12 
13 Bdry1D::Bdry1D() : Domain1D(1, 1, 0.0),
14  m_flow_left(0), m_flow_right(0),
15  m_ilr(0), m_left_nv(0), m_right_nv(0),
16  m_left_loc(0), m_right_loc(0),
17  m_left_points(0), m_nv(0),
18  m_left_nsp(0), m_right_nsp(0),
19  m_sp_left(0), m_sp_right(0),
20  m_start_left(0), m_start_right(0),
21  m_phase_left(0), m_phase_right(0), m_temp(0.0), m_mdot(0.0)
22 {
23  m_type = cConnectorType;
24 }
25 
26 void Bdry1D::
27 _init(size_t n)
28 {
29  if (m_index == npos) {
30  throw CanteraError("Bdry1D",
31  "install in container before calling init.");
32  }
33 
34  // A boundary object contains only one grid point
35  resize(n,1);
36 
37  m_left_nsp = 0;
38  m_right_nsp = 0;
39 
40  // check for left and right flow objects
41  if (m_index > 0) {
42  Domain1D& r = container().domain(m_index-1);
43  if (r.domainType() == cFlowType) {
44  m_flow_left = (StFlow*)&r;
45  m_left_nv = m_flow_left->nComponents();
46  m_left_points = m_flow_left->nPoints();
47  m_left_loc = container().start(m_index-1);
48  m_left_nsp = m_left_nv - 4;
49  m_phase_left = &m_flow_left->phase();
50  } else
51  throw CanteraError("Bdry1D::init",
52  "Boundary domains can only be "
53  "connected on the left to flow domains, not type "+int2str(r.domainType())
54  + " domains.");
55  }
56 
57  // if this is not the last domain, see what is connected on
58  // the right
59  if (m_index + 1 < container().nDomains()) {
60  Domain1D& r = container().domain(m_index+1);
61  if (r.domainType() == cFlowType) {
62  m_flow_right = (StFlow*)&r;
63  m_right_nv = m_flow_right->nComponents();
64  m_right_loc = container().start(m_index+1);
65  m_right_nsp = m_right_nv - 4;
66  m_phase_right = &m_flow_right->phase();
67  } else
68  throw CanteraError("Bdry1D::init",
69  "Boundary domains can only be "
70  "connected on the right to flow domains, not type "+int2str(r.domainType())
71  + " domains.");
72  }
73 }
74 
75 //----------------------------------------------------------
76 // Inlet1D methods
77 //----------------------------------------------------------
78 
79 void Inlet1D::
80 setMoleFractions(const std::string& xin)
81 {
82  m_xstr = xin;
83  if (m_flow) {
84  m_flow->phase().setMoleFractionsByName(xin);
85  m_flow->phase().getMassFractions(DATA_PTR(m_yin));
86  needJacUpdate();
87  }
88 }
89 
90 void Inlet1D::
91 setMoleFractions(doublereal* xin)
92 {
93  if (m_flow) {
94  m_flow->phase().setMoleFractions(xin);
95  m_flow->phase().getMassFractions(DATA_PTR(m_yin));
96  needJacUpdate();
97  }
98 }
99 
100 string Inlet1D::
101 componentName(size_t n) const
102 {
103  switch (n) {
104  case 0:
105  return "mdot";
106  case 1:
107  return "temperature";
108  default:
109  break;
110  }
111  return "unknown";
112 }
113 
114 void Inlet1D::
116 {
117  _init(2);
118 
119  setBounds(0, -1e5, 1e5); // mdot
120  setBounds(1, 200.0, 1e5); // T
121 
122  // set tolerances
123  setSteadyTolerances(1e-4, 1e-5);
124  setTransientTolerances(1e-4, 1e-5);
125 
126  // if a flow domain is present on the left, then this must be
127  // a right inlet. Note that an inlet object can only be a
128  // terminal object - it cannot have flows on both the left and
129  // right
130  if (m_flow_left) {
131  m_ilr = RightInlet;
132  m_flow = m_flow_left;
133  } else if (m_flow_right) {
134  m_ilr = LeftInlet;
135  m_flow = m_flow_right;
136  } else {
137  throw CanteraError("Inlet1D::init","no flow!");
138  }
139 
140  // components = u, V, T, lambda, + mass fractions
141  m_nsp = m_flow->nComponents() - 4;
142  m_yin.resize(m_nsp, 0.0);
143  if (m_xstr != "") {
144  setMoleFractions(m_xstr);
145  } else {
146  m_yin[0] = 1.0;
147  }
148 }
149 
150 void Inlet1D::
151 eval(size_t jg, doublereal* xg, doublereal* rg,
152  integer* diagg, doublereal rdt)
153 {
154  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
155  return;
156  }
157 
158  // start of local part of global arrays
159  doublereal* x = xg + loc();
160  doublereal* r = rg + loc();
161  integer* diag = diagg + loc();
162  doublereal* xb, *rb;
163 
164  // residual equations for the two local variables
165 
166  r[0] = m_mdot - x[0];
167 
168  // Temperature
169  r[1] = m_temp - x[1];
170 
171  // both are algebraic constraints
172  diag[0] = 0;
173  diag[1] = 0;
174 
175 
176  // if it is a left inlet, then the flow solution vector
177  // starts 2 to the right in the global solution vector
178  if (m_ilr == LeftInlet) {
179  xb = x + 2;
180  rb = r + 2;
181 
182  // The first flow residual is for u. This, however, is not
183  // modified by the inlet, since this is set within the flow
184  // domain from the continuity equation.
185 
186  // spreading rate. The flow domain sets this to V(0),
187  // so for finite spreading rate subtract m_V0.
188  rb[1] -= m_V0;
189 
190  // The third flow residual is for T, where it is set to
191  // T(0). Subtract the local temperature to hold the flow
192  // T to the inlet T.
193  rb[2] -= x[1];
194 
195  // The flow domain sets this to -rho*u. Add mdot to
196  // specify the mass flow rate.
197  rb[3] += x[0];
198 
199  // add the convective term to the species residual equations
200  for (size_t k = 1; k < m_nsp; k++) {
201  rb[4+k] += x[0]*m_yin[k];
202  }
203 
204  // if the flow is a freely-propagating flame, mdot is not
205  // specified. Set mdot equal to rho*u, and also set
206  // lambda to zero.
207  if (!m_flow->fixed_mdot()) {
208  m_mdot = m_flow->density(0)*xb[0];
209  r[0] = m_mdot - x[0];
210  rb[3] = xb[3];
211  }
212  }
213 
214  // right inlet.
215  else {
216  size_t boffset = m_flow->nComponents();
217  xb = x - boffset;
218  rb = r - boffset;
219  rb[1] -= m_V0;
220  rb[2] -= x[1]; // T
221  rb[0] += x[0]; // u
222  for (size_t k = 1; k < m_nsp; k++) {
223  rb[4+k] += x[0]*(m_yin[k]);
224  }
225  }
226 }
227 
229 save(XML_Node& o, const doublereal* const soln)
230 {
231  const doublereal* s = soln + loc();
232  XML_Node& inlt = Domain1D::save(o, soln);
233  inlt.addAttribute("type","inlet");
234  for (size_t k = 0; k < nComponents(); k++) {
235  ctml::addFloat(inlt, componentName(k), s[k]);
236  }
237  for (size_t k=0; k < m_nsp; k++) {
238  ctml::addFloat(inlt, "massFraction", m_yin[k], "",
239  m_flow->phase().speciesName(k));
240  }
241  return inlt;
242 }
243 
244 void Inlet1D::
245 restore(const XML_Node& dom, doublereal* soln, int loglevel)
246 {
247  Domain1D::restore(dom, soln, loglevel);
248  soln[0] = m_mdot = ctml::getFloat(dom, "mdot", "massflowrate");
249  soln[1] = m_temp = ctml::getFloat(dom, "temperature", "temperature");
250 
251  m_yin.assign(m_nsp, 0.0);
252 
253  for (size_t i = 0; i < dom.nChildren(); i++) {
254  const XML_Node& node = dom.child(i);
255  if (node.name() == "massFraction") {
256  size_t k = m_flow->phase().speciesIndex(node.attrib("type"));
257  if (k != npos) {
258  m_yin[k] = node.fp_value();
259  }
260  }
261  }
262  resize(2,1);
263 }
264 
265 //--------------------------------------------------
266 // Empty1D
267 //--------------------------------------------------
268 
269 string Empty1D::componentName(size_t n) const
270 {
271  switch (n) {
272  case 0:
273  return "dummy";
274  default:
275  break;
276  }
277  return "<unknown>";
278 }
279 
280 void Empty1D::
281 init() //_init(1);
282 {
283  setBounds(0, -1.0, 1.0);
284 
285  // set tolerances
286  setSteadyTolerances(1e-4, 1e-4);
287  setTransientTolerances(1e-4, 1e-4);
288 }
289 
290 void Empty1D::
291 eval(size_t jg, doublereal* xg, doublereal* rg,
292  integer* diagg, doublereal rdt)
293 {
294  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
295  return;
296  }
297 
298  // start of local part of global arrays
299  doublereal* x = xg + loc();
300  doublereal* r = rg + loc();
301  integer* diag = diagg + loc();
302  // integer *db;
303 
304  r[0] = x[0];
305  diag[0] = 0;
306 }
307 
309 save(XML_Node& o, const doublereal* const soln)
310 {
311  XML_Node& symm = Domain1D::save(o, soln);
312  symm.addAttribute("type","empty");
313  return symm;
314 }
315 
316 void Empty1D::
317 restore(const XML_Node& dom, doublereal* soln, int loglevel)
318 {
319  Domain1D::restore(dom, soln, loglevel);
320  resize(1,1);
321 }
322 
323 //--------------------------------------------------
324 // Symm1D
325 //--------------------------------------------------
326 
327 string Symm1D::componentName(size_t n) const
328 {
329  switch (n) {
330  case 0:
331  return "dummy";
332  default:
333  break;
334  }
335  return "<unknown>";
336 }
337 
338 void Symm1D::
340 {
341  _init(1);
342  setBounds(0, -1.0, 1.0);
343 
344  // set tolerances
345  setSteadyTolerances(1e-4, 1e-4);
346  setTransientTolerances(1e-4, 1e-4);
347 }
348 
349 void Symm1D::
350 eval(size_t jg, doublereal* xg, doublereal* rg,
351  integer* diagg, doublereal rdt)
352 {
353  if (jg != npos && (jg + 2< firstPoint() || jg > lastPoint() + 2)) {
354  return;
355  }
356 
357  // start of local part of global arrays
358  doublereal* x = xg + loc();
359  doublereal* r = rg + loc();
360  integer* diag = diagg + loc();
361  doublereal* xb, *rb;
362  integer* db;
363 
364  r[0] = x[0];
365  diag[0] = 0;
366  size_t nc;
367 
368  if (m_flow_right) {
369  nc = m_flow_right->nComponents();
370  xb = x + 1;
371  rb = r + 1;
372  db = diag + 1;
373  db[1] = 0;
374  db[2] = 0;
375  rb[1] = xb[1] - xb[1 + nc]; // zero dV/dz
376  rb[2] = xb[2] - xb[2 + nc]; // zero dT/dz
377  }
378 
379  if (m_flow_left) {
380  nc = m_flow_left->nComponents();
381  xb = x - nc;
382  rb = r - nc;
383  db = diag - nc;
384  db[1] = 0;
385  db[2] = 0;
386  rb[1] = xb[1] - xb[1 - nc]; // zero dV/dz
387  rb[2] = xb[2] - xb[2 - nc]; // zero dT/dz
388  }
389 }
390 
392 save(XML_Node& o, const doublereal* const soln)
393 {
394  XML_Node& symm = Domain1D::save(o, soln);
395  symm.addAttribute("type","symmetry");
396  return symm;
397 }
398 
399 void Symm1D::
400 restore(const XML_Node& dom, doublereal* soln, int loglevel)
401 {
402  Domain1D::restore(dom, soln, loglevel);
403  resize(1,1);
404 }
405 
406 //--------------------------------------------------
407 // Outlet1D
408 //--------------------------------------------------
409 
410 string Outlet1D::componentName(size_t n) const
411 {
412  switch (n) {
413  case 0:
414  return "outlet dummy";
415  default:
416  break;
417  }
418  return "<unknown>";
419 }
420 
421 void Outlet1D::
423 {
424  _init(1);
425  setBounds(0, -1.0, 1.0);
426 
427  // set tolerances
428  setSteadyTolerances(1e-4, 1e-4);
429  setTransientTolerances(1e-4, 1e-4);
430  if (m_flow_right) {
431  m_flow_right->setViscosityFlag(false);
432  }
433  if (m_flow_left) {
434  m_flow_left->setViscosityFlag(false);
435  }
436 }
437 
438 void Outlet1D::
439 eval(size_t jg, doublereal* xg, doublereal* rg,
440  integer* diagg, doublereal rdt)
441 {
442  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
443  return;
444  }
445 
446  // start of local part of global arrays
447  doublereal* x = xg + loc();
448  doublereal* r = rg + loc();
449  integer* diag = diagg + loc();
450  doublereal* xb, *rb;
451  integer* db;
452 
453  r[0] = x[0];
454  diag[0] = 0;
455  size_t nc, k;
456 
457  if (m_flow_right) {
458  nc = m_flow_right->nComponents();
459  xb = x + 1;
460  rb = r + 1;
461  db = diag + 1;
462  rb[0] = xb[3];
463  rb[2] = xb[2] - xb[2 + nc];
464  for (k = 4; k < nc; k++) {
465  //if (m_flow_right->doSpecies(k-4)) {
466  rb[k] = xb[k] - xb[k + nc];
467  //}
468  }
469  }
470 
471  if (m_flow_left) {
472  nc = m_flow_left->nComponents();
473  xb = x - nc;
474  rb = r - nc;
475  db = diag - nc;
476 
477  // zero Lambda
478 
479  if (!m_flow_left->fixed_mdot()) {
480  ; // rb[0] = xb[0] - xb[0-nc]; //zero U gradient
481  } else {
482  rb[0] = xb[3]; // zero Lambda
483  }
484 
485  rb[2] = xb[2] - xb[2 - nc]; // zero T gradient
486  for (k = 5; k < nc; k++) {
487  rb[k] = xb[k] - xb[k - nc]; // zero mass fraction gradient
488  db[k] = 0;
489  }
490  }
491 }
492 
494 save(XML_Node& o, const doublereal* const soln)
495 {
496  XML_Node& outlt = Domain1D::save(o, soln);
497  outlt.addAttribute("type","outlet");
498  return outlt;
499 }
500 
501 void Outlet1D::
502 restore(const XML_Node& dom, doublereal* soln, int loglevel)
503 {
504  Domain1D::restore(dom, soln, loglevel);
505  resize(1,1);
506 }
507 
508 //--------------------------------------------------
509 // OutletRes1D
510 //--------------------------------------------------
511 
512 void OutletRes1D::
513 setMoleFractions(const std::string& xres)
514 {
515  m_xstr = xres;
516  if (m_flow) {
517  m_flow->phase().setMoleFractionsByName(xres);
518  m_flow->phase().getMassFractions(DATA_PTR(m_yres));
519  needJacUpdate();
520  }
521 }
522 
523 void OutletRes1D::
524 setMoleFractions(doublereal* xres)
525 {
526  if (m_flow) {
527  m_flow->phase().setMoleFractions(xres);
528  m_flow->phase().getMassFractions(DATA_PTR(m_yres));
529  needJacUpdate();
530  }
531 }
532 
533 string OutletRes1D::componentName(size_t n) const
534 {
535  switch (n) {
536  case 0:
537  return "dummy";
538  default:
539  break;
540  }
541  return "<unknown>";
542 }
543 
544 void OutletRes1D::
546 {
547  _init(1);
548  // set bounds (dummy)
549  setBounds(0, -1.0, 1.0);
550 
551  // set tolerances
552  setSteadyTolerances(1e-4, 1e-4);
553  setTransientTolerances(1e-4, 1e-4);
554 
555  if (m_flow_left) {
556  m_flow = m_flow_left;
557  } else if (m_flow_right) {
558  m_flow = m_flow_right;
559  } else {
560  throw CanteraError("OutletRes1D::init","no flow!");
561  }
562 
563  m_nsp = m_flow->nComponents() - 4;
564  m_yres.resize(m_nsp, 0.0);
565  if (m_xstr != "") {
566  setMoleFractions(m_xstr);
567  } else {
568  m_yres[0] = 1.0;
569  }
570 }
571 
572 void OutletRes1D::
573 eval(size_t jg, doublereal* xg, doublereal* rg,
574  integer* diagg, doublereal rdt)
575 {
576 
577  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
578  return;
579  }
580 
581  // start of local part of global arrays
582  doublereal* x = xg + loc();
583  doublereal* r = rg + loc();
584  integer* diag = diagg + loc();
585  doublereal* xb, *rb;
586  integer* db;
587 
588  // drive dummy component to zero
589  r[0] = x[0];
590  diag[0] = 0;
591  size_t nc, k;
592 
593  if (m_flow_right) {
594  nc = m_flow_right->nComponents();
595  xb = x + 1;
596  rb = r + 1;
597  db = diag + 1;
598 
599  // this seems wrong...
600  // zero Lambda
601  rb[0] = xb[3];
602 
603  // zero gradient for T
604  rb[2] = xb[2] - xb[2 + nc];
605 
606  // specified mass fractions
607  for (k = 4; k < nc; k++) {
608  rb[k] = xb[k] - m_yres[k-4];
609  }
610  }
611 
612  if (m_flow_left) {
613 
614  nc = m_flow_left->nComponents();
615  xb = x - nc;
616  rb = r - nc;
617  db = diag - nc;
618 
619  if (!m_flow_left->fixed_mdot()) {
620  ;
621  } else {
622  rb[0] = xb[3]; // zero Lambda
623  }
624  rb[2] = xb[2] - m_temp; //xb[2] - xb[2 - nc]; // zero dT/dz
625  for (k = 5; k < nc; k++) {
626  rb[k] = xb[k] - m_yres[k-4]; // fixed Y
627  db[k] = 0;
628  }
629  }
630 }
631 
633 save(XML_Node& o, const doublereal* const soln)
634 {
635  XML_Node& outlt = Domain1D::save(o, soln);
636  outlt.addAttribute("type","outletres");
637  ctml::addFloat(outlt, "temperature", m_temp, "K");
638  for (size_t k=0; k < m_nsp; k++) {
639  ctml::addFloat(outlt, "massFraction", m_yres[k], "",
640  m_flow->phase().speciesName(k));
641  }
642  return outlt;
643 }
644 
645 void OutletRes1D::
646 restore(const XML_Node& dom, doublereal* soln, int loglevel)
647 {
648  Domain1D::restore(dom, soln, loglevel);
649  m_temp = ctml::getFloat(dom, "temperature");
650 
651  m_yres.assign(m_nsp, 0.0);
652  for (size_t i = 0; i < dom.nChildren(); i++) {
653  const XML_Node& node = dom.child(i);
654  if (node.name() == "massFraction") {
655  size_t k = m_flow->phase().speciesIndex(node.attrib("type"));
656  if (k != npos) {
657  m_yres[k] = node.fp_value();
658  }
659  }
660  }
661 
662  resize(1,1);
663 }
664 
665 //-----------------------------------------------------------
666 // Surf1D
667 //-----------------------------------------------------------
668 
669 string Surf1D::componentName(size_t n) const
670 {
671  switch (n) {
672  case 0:
673  return "temperature";
674  default:
675  break;
676  }
677  return "<unknown>";
678 }
679 
680 void Surf1D::
682 {
683  _init(1);
684  // set bounds (T)
685  setBounds(0, 200.0, 1e5);
686 
687  // set tolerances
688  setSteadyTolerances(1e-4, 1e-4);
689  setTransientTolerances(1e-4, 1e-4);
690 }
691 
692 void Surf1D::
693 eval(size_t jg, doublereal* xg, doublereal* rg,
694  integer* diagg, doublereal rdt)
695 {
696  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
697  return;
698  }
699 
700  // start of local part of global arrays
701  doublereal* x = xg + loc();
702  doublereal* r = rg + loc();
703  integer* diag = diagg + loc();
704  doublereal* xb, *rb;
705 
706  r[0] = x[0] - m_temp;
707  diag[0] = 0;
708  size_t nc;
709 
710  if (m_flow_right) {
711  rb = r + 1;
712  xb = x + 1;
713  rb[2] = xb[2] - x[0]; // specified T
714  }
715 
716  if (m_flow_left) {
717  nc = m_flow_left->nComponents();
718  rb = r - nc;
719  xb = x - nc;
720  rb[2] = xb[2] - x[0]; // specified T
721  }
722 }
723 
725 save(XML_Node& o, const doublereal* const soln)
726 {
727  const doublereal* s = soln + loc();
728  //XML_Node& inlt = o.addChild("inlet");
729  XML_Node& inlt = Domain1D::save(o, soln);
730  inlt.addAttribute("type","surface");
731  for (size_t k = 0; k < nComponents(); k++) {
732  ctml::addFloat(inlt, componentName(k), s[k]);
733  }
734  return inlt;
735 }
736 
737 void Surf1D::
738 restore(const XML_Node& dom, doublereal* soln, int loglevel)
739 {
740  Domain1D::restore(dom, soln, loglevel);
741  soln[0] = m_temp = ctml::getFloat(dom, "temperature", "temperature");
742  resize(1,1);
743 }
744 
745 //-----------------------------------------------------------
746 // ReactingSurf1D
747 //-----------------------------------------------------------
748 
749 string ReactingSurf1D::componentName(size_t n) const
750 {
751  if (n == 0) {
752  return "temperature";
753  } else if (n < m_nsp + 1) {
754  return m_sphase->speciesName(n-1);
755  } else {
756  return "<unknown>";
757  }
758 }
759 
760 void ReactingSurf1D::
762 {
763  m_nv = m_nsp + 1;
764  _init(m_nsp+1);
765  m_fixed_cov.resize(m_nsp, 0.0);
766  m_fixed_cov[0] = 1.0;
767  m_work.resize(m_kin->nTotalSpecies(), 0.0);
768 
769  setBounds(0, 200.0, 1e5);
770  for (size_t n = 0; n < m_nsp; n++) {
771  setBounds(n+1, -1.0e-5, 2.0);
772  }
773  setSteadyTolerances(1.0e-5, 1.0e-9);
774  setTransientTolerances(1.0e-5, 1.0e-9);
775  setSteadyTolerances(1.0e-5, 1.0e-4, 0);
776  setTransientTolerances(1.0e-5, 1.0e-4, 0);
777 }
778 
779 void ReactingSurf1D::
780 eval(size_t jg, doublereal* xg, doublereal* rg,
781  integer* diagg, doublereal rdt)
782 {
783  if (jg != npos && (jg + 2 < firstPoint() || jg > lastPoint() + 2)) {
784  return;
785  }
786 
787  // start of local part of global arrays
788  doublereal* x = xg + loc();
789  doublereal* r = rg + loc();
790  integer* diag = diagg + loc();
791  doublereal* xb, *rb;
792 
793  // specified surface temp
794  r[0] = x[0] - m_temp;
795 
796  // set the coverages
797  doublereal sum = 0.0;
798  for (size_t k = 0; k < m_nsp; k++) {
799  m_work[k] = x[k+1];
800  sum += x[k+1];
801  }
802  m_sphase->setTemperature(x[0]);
803  m_sphase->setCoverages(DATA_PTR(m_work));
804  //m_kin->advanceCoverages(1.0);
805  //m_sphase->getCoverages(m_fixed_cov.begin());
806 
807  // set the left gas state to the adjacent point
808 
809  size_t leftloc = 0, rightloc = 0;
810  size_t pnt = 0;
811 
812  if (m_flow_left) {
813  leftloc = m_flow_left->loc();
814  pnt = m_flow_left->nPoints() - 1;
815  m_flow_left->setGas(xg + leftloc, pnt);
816  }
817 
818  if (m_flow_right) {
819  rightloc = m_flow_right->loc();
820  m_flow_right->setGas(xg + rightloc, 0);
821  }
822 
823  m_kin->getNetProductionRates(DATA_PTR(m_work));
824  doublereal rs0 = 1.0/m_sphase->siteDensity();
825 
826  //scale(m_work.begin(), m_work.end(), m_work.begin(), m_mult[0]);
827 
828  // bool enabled = true;
829  size_t ioffset = m_kin->kineticsSpeciesIndex(0, m_surfindex);
830 
831  if (m_enabled) {
832  doublereal maxx = -1.0;
833  for (size_t k = 0; k < m_nsp; k++) {
834  r[k+1] = m_work[k + ioffset] * m_sphase->size(k) * rs0;
835  r[k+1] -= rdt*(x[k+1] - prevSoln(k+1,0));
836  diag[k+1] = 1;
837  if (x[k+1] > maxx) {
838  maxx = x[k+1];
839  }
840  }
841  r[1] = 1.0 - sum;
842  diag[1] = 0;
843  } else {
844  for (size_t k = 0; k < m_nsp; k++) {
845  r[k+1] = x[k+1] - m_fixed_cov[k];
846  diag[k+1] = 0;
847  }
848  }
849 
850  if (m_flow_right) {
851  rb = r + 1;
852  xb = x + 1;
853  rb[2] = xb[2] - x[0]; // specified T
854  }
855  size_t nc;
856  if (m_flow_left) {
857  nc = m_flow_left->nComponents();
858  const doublereal* mwleft = DATA_PTR(m_phase_left->molecularWeights());
859  rb =r - nc;
860  xb = x - nc;
861  rb[2] = xb[2] - x[0]; // specified T
862  for (size_t nl = 1; nl < m_left_nsp; nl++) {
863  rb[4+nl] += m_work[nl]*mwleft[nl];
864  }
865  }
866 }
867 
869 save(XML_Node& o, const doublereal* const soln)
870 {
871  const doublereal* s = soln + loc();
872  XML_Node& dom = Domain1D::save(o, soln);
873  dom.addAttribute("type","surface");
874  ctml::addFloat(dom, "temperature", s[0], "K");
875  for (size_t k=0; k < m_nsp; k++) {
876  ctml::addFloat(dom, "coverage", s[k+1], "",
877  m_sphase->speciesName(k));
878  }
879  return dom;
880 }
881 
882 void ReactingSurf1D::
883 restore(const XML_Node& dom, doublereal* soln, int loglevel)
884 {
885  Domain1D::restore(dom, soln, loglevel);
886  soln[0] = m_temp = ctml::getFloat(dom, "temperature");
887 
888  m_fixed_cov.assign(m_nsp, 0.0);
889  for (size_t i = 0; i < dom.nChildren(); i++) {
890  const XML_Node& node = dom.child(i);
891  if (node.name() == "coverage") {
892  size_t k = m_sphase->speciesIndex(node.attrib("type"));
893  if (k != npos) {
894  m_fixed_cov[k] = soln[k+1] = node.fp_value();
895  }
896  }
897  }
898  m_sphase->setCoverages(&m_fixed_cov[0]);
899 
900  resize(m_nsp+1,1);
901 }
902 }
virtual void eval(size_t jg, doublereal *xg, doublereal *rg, integer *diagg, doublereal rdt)
Evaluate the residual function at point j.
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
size_t nPoints() const
Number of grid points in this domain.
Definition: Domain1D.h:186
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
virtual void getNetProductionRates(doublereal *net)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:534
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
void getMassFractions(doublereal *const y) const
Get the species mass fractions.
Definition: Phase.cpp:561
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:173
size_t firstPoint() const
The index of the first (i.e., left-most) grid point belonging to this domain.
Definition: Domain1D.h:469
virtual XML_Node & save(XML_Node &o, const doublereal *const soln)
Save the current solution for this domain into an XML_Node.
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
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
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
void setCoverages(const doublereal *theta)
Set the surface site fractions to a specified state.
Definition: SurfPhase.cpp:283
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
doublereal size(size_t k) const
This routine returns the size of species k.
Definition: Phase.h:398
virtual void init()
size_t nTotalSpecies() const
The total number of species in all phases participating in the kinetics mechanism.
Definition: Kinetics.h:300
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 eval(size_t jg, doublereal *xg, doublereal *rg, integer *diagg, doublereal rdt)
Evaluate the residual function at point j.
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
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
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 soln)
Save the current solution for this domain into an XML_Node.
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
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.
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
size_t nComponents() const
Number of components at each grid point.
Definition: Domain1D.h:164
virtual XML_Node & save(XML_Node &o, const doublereal *const soln)
Save the current solution for this domain into an XML_Node.
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:229
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.
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
Definition: OneDim.h:78
size_t lastPoint() const
The index of the last (i.e., right-most) grid point belonging to this domain.
Definition: Domain1D.h:477
Domain1D & domain(size_t i) const
Return a reference to domain i.
Definition: OneDim.h:54
virtual void init()
virtual XML_Node & save(XML_Node &o, const doublereal *const soln)
Save the current solution for this domain into an XML_Node.
std::string name() const
Returns the name of the XML node.
Definition: xml.h:390
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
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 setMoleFractionsByName(compositionMap &xMap)
Set the species mole fractions by name.
Definition: Phase.cpp:354
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values There is no restriction on the sum of the mole fractio...
Definition: Phase.cpp:306
virtual void resize(size_t nv, size_t np)
Definition: Domain1D.h:135
virtual XML_Node & save(XML_Node &o, const doublereal *const soln)
Save the current solution for this domain into an XML_Node.
virtual void setMoleFractions(const std::string &xin)
Set the mole fractions by specifying a std::string.
void addAttribute(const std::string &attrib, const std::string &value)
Add or modify an attribute of the current node.
Definition: xml.cpp:518
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:331
virtual void eval(size_t jg, doublereal *xg, doublereal *rg, integer *diagg, doublereal rdt)
Evaluate the residual function at point j.
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:238
Domain1D(size_t nv=1, size_t points=1, doublereal time=0.0)
Constructor.
Definition: Domain1D.h:46
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)
Set tolerances for steady-state mode.
Definition: Domain1D.cpp:88
virtual void eval(size_t jg, doublereal *xg, doublereal *rg, integer *diagg, doublereal rdt)
Evaluate the residual function at point j.
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:562
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
virtual void init()
const OneDim & container() const
The container holding this domain.
Definition: Domain1D.h:83
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
void needJacUpdate()
Definition: OneDim.cpp:323
doublereal fp_value() const
Return the value of an XML node as a single double.
Definition: xml.cpp:498
virtual XML_Node & save(XML_Node &o, const doublereal *const soln)
Save the current solution for this domain into an XML_Node.
double prevSoln(size_t n, size_t j) const
Value of component n at point j in the previous solution.
Definition: Domain1D.h:512
void setTransientTolerances(doublereal rtol, doublereal atol, size_t n=npos)
Set tolerances for time-stepping mode.
Definition: Domain1D.cpp:68
virtual void init()
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:246
size_t nChildren(bool discardComments=false) const
Return the number of children.
Definition: xml.cpp:594
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
doublereal siteDensity()
Returns the site density.
Definition: SurfPhase.h:415