Cantera  3.2.0a4
Loading...
Searching...
No Matches
SolutionArray.cpp
Go to the documentation of this file.
1/**
2 * @file SolutionArray.cpp
3 * Definition file for class SolutionArray.
4 */
5
6// This file is part of Cantera. See License.txt in the top-level directory or
7// at https://cantera.org/license.txt for license and copyright information.
8
9#include "cantera/base/SolutionArray.h"
10#include "cantera/base/Solution.h"
11#include "cantera/base/Storage.h"
12#include "cantera/base/stringUtils.h"
13#include "cantera/thermo/ThermoPhase.h"
14#include "cantera/thermo/SurfPhase.h"
15#include "cantera/base/utilities.h"
16#include <boost/algorithm/string.hpp>
17#include <boost/range/adaptor/reversed.hpp>
18#include <fstream>
19#include <sstream>
20
21
22namespace ba = boost::algorithm;
23
24namespace Cantera
25{
26
27namespace { // restrict scope of auxiliary variables to local translation unit
28
29const map<string, string> aliasMap = {
30 {"T", "temperature"},
31 {"P", "pressure"},
32 {"D", "density"},
33 {"Y", "mass-fractions"},
34 {"X", "mole-fractions"},
35 {"C", "coverages"},
36 {"U", "specific-internal-energy"},
37 {"V", "specific-volume"},
38 {"H", "specific-enthalpy"},
39 {"S", "specific-entropy"},
40 {"Q", "vapor-fraction"},
41};
42
43const map<string, string> reverseAliasMap = {
44 // reserved names used by states
45 {"temperature", "T"},
46 {"pressure", "P"},
47 {"density", "D"},
48 // reserved names used for 1-D objects
49 {"spread-rate", "spreadRate"},
50 {"spread_rate", "spreadRate"}, // snake_case version used prior to Cantera 3.2
51 {"L", "Lambda"},
52 {"radial-pressure-gradient", "Lambda"},
53 {"lambda", "Lambda"}, // lower-case version used prior to Cantera 3.2
54 {"E", "eField"},
55 {"electric-field", "eField"},
56 {"oxidizer-velocity", "Uo"},
57};
58
59} // end unnamed namespace
60
61const map<string, string>& _componentAliasMap()
62{
63 return reverseAliasMap;
64}
65
66SolutionArray::SolutionArray(const shared_ptr<Solution>& sol,
67 int size, const AnyMap& meta)
68 : m_sol(sol)
69 , m_size(size)
70 , m_dataSize(size)
71 , m_meta(meta)
72{
73 if (!m_sol || !m_sol->thermo()) {
74 throw CanteraError("SolutionArray::SolutionArray",
75 "Unable to create SolutionArray from invalid Solution object.");
76 }
77 m_stride = m_sol->thermo()->stateSize();
78 m_sol->thermo()->addSpeciesLock();
79 m_data = make_shared<vector<double>>(m_dataSize * m_stride, 0.);
80 m_extra = make_shared<map<string, AnyValue>>();
81 m_order = make_shared<map<int, string>>();
82 for (size_t i = 0; i < m_dataSize; ++i) {
83 m_active.push_back(static_cast<int>(i));
84 }
85 reset();
86 m_apiShape.resize(1);
87 m_apiShape[0] = static_cast<long>(m_dataSize);
88}
89
90SolutionArray::SolutionArray(const SolutionArray& other,
91 const vector<int>& selected)
92 : m_sol(other.m_sol)
93 , m_size(selected.size())
94 , m_dataSize(other.m_data->size())
95 , m_stride(other.m_stride)
96 , m_data(other.m_data)
97 , m_extra(other.m_extra)
98 , m_order(other.m_order)
99 , m_shared(true)
100{
101 m_sol->thermo()->addSpeciesLock();
102 if (!other.m_shared) {
103 // direct slicing is possible
104 m_active = selected;
105 } else {
106 // slicing a previously sliced SolutionArray
107 m_active.clear();
108 m_active.reserve(selected.size());
109 for (auto loc : selected) {
110 m_active.push_back(other.m_active.at(loc));
111 }
112 }
113 for (auto loc : m_active) {
114 if (loc < 0 || loc >= (int)m_dataSize) {
115 IndexError("SolutionArray::SolutionArray", "indices", loc, m_dataSize);
116 }
117 }
118 set<int> unique(selected.begin(), selected.end());
119 if (unique.size() < selected.size()) {
120 throw CanteraError("SolutionArray::SolutionArray", "Indices must be unique.");
121 }
122}
123
124SolutionArray::~SolutionArray()
125{
126 m_sol->thermo()->removeSpeciesLock();
127}
128
129namespace { // restrict scope of helper functions to local translation unit
130
131template<class T>
132void resetSingle(AnyValue& extra, const vector<int>& slice);
133
134template<class T>
135AnyValue getSingle(const AnyValue& extra, const vector<int>& slice);
136
137template<class T>
138void setSingle(AnyValue& extra, const AnyValue& data, const vector<int>& slice);
139
140template<class T>
141void resizeSingle(AnyValue& extra, size_t size, const AnyValue& value);
142
143template<class T>
144void resetMulti(AnyValue& extra, const vector<int>& slice);
145
146template<class T>
147AnyValue getMulti(const AnyValue& extra, const vector<int>& slice);
148
149template<class T>
150void setMulti(AnyValue& extra, const AnyValue& data, const vector<int>& slice);
151
152template<class T>
153void resizeMulti(AnyValue& extra, size_t size, const AnyValue& value);
154
155template<class T>
156void setAuxiliarySingle(size_t loc, AnyValue& extra, const AnyValue& value);
157
158template<class T>
159void setAuxiliaryMulti(size_t loc, AnyValue& extra, const AnyValue& data);
160
161template<class T>
162void setScalar(AnyValue& extra, const AnyValue& data, const vector<int>& slice);
163
164} // end unnamed namespace
165
166void SolutionArray::reset()
167{
168 size_t nState = m_sol->thermo()->stateSize();
169 vector<double> state(nState);
170 m_sol->thermo()->saveState(state); // thermo contains current state
171 for (size_t k = 0; k < m_size; ++k) {
172 std::copy(state.begin(), state.end(), m_data->data() + m_active[k] * m_stride);
173 }
174 for (auto& [key, extra] : *m_extra) {
175 if (extra.is<void>()) {
176 // cannot reset placeholder (uninitialized component)
177 } else if (extra.isVector<double>()) {
178 resetSingle<double>(extra, m_active);
179 } else if (extra.isVector<long int>()) {
180 resetSingle<long int>(extra, m_active);
181 } else if (extra.isVector<string>()) {
182 resetSingle<string>(extra, m_active);
183 } else if (extra.isVector<vector<double>>()) {
184 resetMulti<double>(extra, m_active);
185 } else if (extra.isVector<vector<long int>>()) {
186 resetMulti<long int>(extra, m_active);
187 } else if (extra.isVector<vector<string>>()) {
188 resetMulti<string>(extra, m_active);
189 } else {
190 throw NotImplementedError("SolutionArray::reset",
191 "Unable to reset component '{}' with type '{}'.",
192 key, extra.type_str());
193 }
194 }
195}
196
197void SolutionArray::resize(int size)
198{
199 if (apiNdim() > 1) {
200 throw CanteraError("SolutionArray::resize",
201 "Resize is ambiguous for multi-dimensional arrays; use setApiShape "
202 "instead.");
203 }
204 if (m_data.use_count() > 1) {
205 throw CanteraError("SolutionArray::resize",
206 "Unable to resize as data are shared by multiple objects.");
207 }
208 _resize(static_cast<size_t>(size));
209 m_apiShape[0] = static_cast<long>(size);
210}
211
212void SolutionArray::setApiShape(const vector<long int>& shape)
213{
214 size_t size = 1;
215 for (auto dim : shape) {
216 size *= dim;
217 }
218 if (m_shared && size != m_size) {
219 throw CanteraError("SolutionArray::setApiShape",
220 "Unable to set shape of shared data as sizes are inconsistent:\n"
221 "active size is {} but shape implies {}.", m_size, size);
222 }
223 if (!m_shared && size != m_dataSize) {
224 if (m_data.use_count() > 1) {
225 throw CanteraError("SolutionArray::setApiShape",
226 "Unable to set shape as data are shared by multiple objects.");
227 }
228 _resize(size);
229 }
230 m_apiShape = shape;
231}
232
233void SolutionArray::_resize(size_t size)
234{
235 m_size = size;
236 m_dataSize = size;
237 m_data->resize(m_dataSize * m_stride, 0.);
238 for (auto& [key, data] : *m_extra) {
239 _resizeExtra(key);
240 }
241 m_active.clear();
242 for (size_t i = 0; i < m_dataSize; ++i) {
243 m_active.push_back(static_cast<int>(i));
244 }
245}
246
247namespace { // restrict scope of helper functions to local translation unit
248
249vector<string> doubleColumn(string name, const vector<double>& comp,
250 int rows, int width)
251{
252 // extract data for processing
253 vector<double> data;
254 vector<string> raw;
255 string notation = fmt::format("{{:{}.{}g}}", width, (width - 1) / 2);
256 int csize = static_cast<int>(comp.size());
257 int dots = csize + 1;
258 if (csize <= rows) {
259 for (const auto& val : comp) {
260 data.push_back(val);
261 raw.push_back(boost::trim_copy(fmt::format(fmt::runtime(notation), val)));
262 }
263 } else {
264 dots = (rows + 1) / 2;
265 for (int row = 0; row < dots; row++) {
266 data.push_back(comp[row]);
267 raw.push_back(boost::trim_copy(
268 fmt::format(fmt::runtime(notation), comp[row])));
269 }
270 for (int row = csize - rows / 2; row < csize; row++) {
271 data.push_back(comp[row]);
272 raw.push_back(boost::trim_copy(
273 fmt::format(fmt::runtime(notation), comp[row])));
274 }
275 }
276
277 // determine notation; all entries use identical formatting
278 bool isFloat = false;
279 bool isScientific = false;
280 size_t head = 0;
281 size_t tail = 0;
282 size_t exp = 0;
283 for (const auto& repr : raw) {
284 string name = repr;
285 if (name[0] == '-') {
286 // leading negative sign is not considered
287 name = name.substr(1);
288 }
289 if (name.find('e') != string::npos) {
290 // scientific notation
291 if (!isScientific) {
292 head = 1;
293 tail = name.find('e') - name.find('.');
294 exp = 4; // size of exponent
295 } else {
296 tail = std::max(tail, name.find('e') - name.find('.'));
297 }
298 isFloat = true;
299 isScientific = true;
300 } else if (name.find('.') != string::npos) {
301 // floating point notation
302 isFloat = true;
303 if (!isScientific) {
304 head = std::max(head, name.find('.'));
305 tail = std::max(tail, name.size() - name.find('.'));
306 }
307 } else {
308 head = std::max(head, name.size());
309 }
310 }
311 size_t maxLen = std::max((size_t)4, head + tail + exp + isFloat + 1);
312 size_t over = std::max(0, (int)name.size() - (int)maxLen);
313 if (isScientific) {
314 // at least one entry has scientific notation
315 notation = fmt::format(" {{:>{}.{}e}}", over + maxLen, tail);
316 } else if (isFloat) {
317 // at least one entry is a floating point
318 notation = fmt::format(" {{:>{}.{}f}}", over + maxLen, tail);
319 } else {
320 // all entries are integers
321 notation = fmt::format(" {{:>{}.0f}}", over + maxLen);
322 }
323 maxLen = fmt::format(fmt::runtime(notation), 0.).size();
324
325 // assemble output
326 string section = fmt::format("{{:>{}}}", maxLen);
327 vector<string> col = {fmt::format(fmt::runtime(section), name)};
328 int count = 0;
329 for (const auto& val : data) {
330 col.push_back(fmt::format(fmt::runtime(notation), val));
331 count++;
332 if (count == dots) {
333 col.push_back(fmt::format(fmt::runtime(section), "..."));
334 }
335 }
336 return col;
337}
338
339vector<string> integerColumn(string name, const vector<long int>& comp,
340 int rows, int width)
341{
342 // extract data for processing
343 vector<long int> data;
344 string notation = fmt::format("{{:{}}}", width);
345 size_t maxLen = 2; // minimum column width is 2
346 int csize = static_cast<int>(comp.size());
347 int dots = csize + 1;
348 if (csize <= rows) {
349 for (const auto& val : comp) {
350 data.push_back(val);
351 string formatted = boost::trim_copy(
352 fmt::format(fmt::runtime(notation), val));
353 if (formatted[0] == '-') {
354 formatted = formatted.substr(1);
355 }
356 maxLen = std::max(maxLen, formatted.size());
357 }
358 } else {
359 dots = (rows + 1) / 2;
360 for (int row = 0; row < dots; row++) {
361 data.push_back(comp[row]);
362 string formatted = boost::trim_copy(
363 fmt::format(fmt::runtime(notation), comp[row]));
364 if (formatted[0] == '-') {
365 formatted = formatted.substr(1);
366 }
367 maxLen = std::max(maxLen, formatted.size());
368 }
369 for (int row = csize - rows / 2; row < csize; row++) {
370 data.push_back(comp[row]);
371 string formatted = boost::trim_copy(
372 fmt::format(fmt::runtime(notation), comp[row]));
373 if (formatted[0] == '-') {
374 formatted = formatted.substr(1);
375 }
376 maxLen = std::max(maxLen, formatted.size());
377 }
378 }
379
380 if (name == "") {
381 // index column
382 notation = fmt::format("{{:<{}}}", maxLen);
383 } else {
384 // regular column
385 maxLen = std::max(maxLen, name.size());
386 notation = fmt::format(" {{:>{}}}", maxLen + 1);
387 }
388
389 // assemble output
390 vector<string> col = {fmt::format(fmt::runtime(notation), name)};
391 int count = 0;
392 for (const auto& val : data) {
393 col.push_back(fmt::format(fmt::runtime(notation), val));
394 count++;
395 if (count == dots) {
396 col.push_back(fmt::format(fmt::runtime(notation), ".."));
397 }
398 }
399 return col;
400}
401
402vector<string> stringColumn(string name, const vector<string>& comp,
403 int rows, int width)
404{
405 // extract data for processing
406 vector<string> data;
407 string notation = fmt::format("{{:{}}}", width);
408 size_t maxLen = 3; // minimum column width is 3
409 int csize = static_cast<int>(comp.size());
410 int dots = csize + 1;
411 if (csize <= rows) {
412 for (const auto& val : comp) {
413 data.push_back(val);
414 maxLen = std::max(maxLen,
415 boost::trim_copy(fmt::format(fmt::runtime(notation), val)).size());
416 }
417 } else {
418 dots = (rows + 1) / 2;
419 for (int row = 0; row < dots; row++) {
420 data.push_back(comp[row]);
421 maxLen = std::max(maxLen,
422 boost::trim_copy(
423 fmt::format(fmt::runtime(notation), comp[row])).size());
424 }
425 for (int row = csize - rows / 2; row < csize; row++) {
426 data.push_back(comp[row]);
427 maxLen = std::max(maxLen,
428 boost::trim_copy(
429 fmt::format(fmt::runtime(notation), comp[row])).size());
430 }
431 }
432
433 // assemble output
434 notation = fmt::format(" {{:>{}}}", maxLen);
435 vector<string> col = {fmt::format(fmt::runtime(notation), name)};
436 int count = 0;
437 for (const auto& val : data) {
438 col.push_back(fmt::format(fmt::runtime(notation), val));
439 count++;
440 if (count == dots) {
441 col.push_back(fmt::format(fmt::runtime(notation), "..."));
442 }
443 }
444 return col;
445}
446
447vector<string> formatColumn(string name, const AnyValue& comp, int rows, int width)
448{
449 if (comp.isVector<double>()) {
450 return doubleColumn(name, comp.asVector<double>(), rows, width);
451 }
452 if (comp.isVector<long int>()) {
453 return integerColumn(name, comp.asVector<long int>(), rows, width);
454 }
455 if (comp.isVector<string>()) {
456 return stringColumn(name, comp.asVector<string>(), rows, width);
457 }
458
459 // create alternative representation
460 string repr;
461 int size;
462 if (comp.isVector<vector<double>>()) {
463 repr = "[ <double> ]";
464 size = len(comp.asVector<vector<double>>());
465 } else if (comp.isVector<vector<long int>>()) {
466 repr = "[ <long int> ]";
467 size = len(comp.asVector<vector<long int>>());
468 } else if (comp.isVector<vector<string>>()) {
469 repr = "[ <string> ]";
470 size = len(comp.asVector<vector<string>>());
471 } else {
472 throw CanteraError(
473 "formatColumn", "Encountered invalid data for component '{}'.", name);
474 }
475 size_t maxLen = std::max(repr.size(), name.size());
476
477 // assemble output
478 string notation = fmt::format(" {{:>{}}}", maxLen);
479 repr = fmt::format(fmt::runtime(notation), repr);
480 vector<string> col = {fmt::format(fmt::runtime(notation), name)};
481 if (size <= rows) {
482 for (int row = 0; row < size; row++) {
483 col.push_back(repr);
484 }
485 } else {
486 int dots = (rows + 1) / 2;
487 for (int row = 0; row < dots; row++) {
488 col.push_back(repr);
489 }
490 col.push_back(fmt::format(fmt::runtime(notation), "..."));
491 for (int row = size - rows / 2; row < size; row++) {
492 col.push_back(repr);
493 }
494 }
495 return col;
496}
497
498} // end unnamed namespace
499
500string SolutionArray::info(const vector<string>& keys, int rows, int width)
501{
502 fmt::memory_buffer b;
503 int col_width = 12; // targeted maximum column width
504 vector<string> components;
505 if (keys.size()) {
506 components = keys;
507 } else {
508 components = componentNames();
509 }
510 try {
511 // build columns
512 vector<long int> index;
513 for (const auto ix : m_active) {
514 index.push_back(ix);
515 }
516 vector<vector<string>> cols = {integerColumn("", index, rows, col_width)};
517 vector<vector<string>> tail; // trailing columns in reverse order
518 size_t size = cols.back().size();
519
520 // assemble columns fitting within a maximum width; if this width is exceeded,
521 // a "..." separator is inserted close to the center. Accordingly, the matrix
522 // needs to be assembled from two halves.
523 int front = 0;
524 int back = len(components) - 1;
525 int fLen = len(cols.back()[0]);
526 int bLen = 0;
527 int sep = 5; // separator width
528 bool done = false;
529 while (!done && front <= back) {
530 string key;
531 while (bLen + sep <= fLen && front <= back) {
532 // add trailing columns
533 key = components[back];
534 auto col = formatColumn(key, getComponent(key), rows, col_width);
535 if (len(col[0]) + fLen + bLen + sep > width) {
536 done = true;
537 break;
538 }
539 tail.push_back(col);
540 bLen += len(tail.back()[0]);
541 back--;
542 }
543 if (done || front > back) {
544 break;
545 }
546 while (fLen <= bLen + sep && front <= back) {
547 // add leading columns
548 key = components[front];
549 auto col = formatColumn(key, getComponent(key), rows, col_width);
550 if (len(col[0]) + fLen + bLen + sep > width) {
551 done = true;
552 break;
553 }
554 cols.push_back(col);
555 fLen += len(cols.back()[0]);
556 front++;
557 }
558 }
559 if (cols.size() + tail.size() < components.size() + 1) {
560 // add separator
561 cols.push_back(vector<string>(size + 1, " ..."));
562 }
563 // copy trailing columns
564 cols.insert(cols.end(), tail.rbegin(), tail.rend());
565
566 // assemble formatted output
567 for (size_t row = 0; row < size; row++) {
568 for (const auto& col : cols) {
569 fmt_append(b, col[row]);
570 }
571 fmt_append(b, "\n");
572 }
573
574 // add size information
575 fmt_append(b, "\n[{} rows x {} components; state='{}']",
576 m_size, components.size(), m_sol->thermo()->nativeMode());
577
578 } catch (CanteraError& err) {
579 return to_string(b) + err.what();
580 }
581 return to_string(b);
582}
583
584shared_ptr<ThermoPhase> SolutionArray::thermo()
585{
586 return m_sol->thermo();
587}
588
589string SolutionArray::transportModel()
590{
591 return m_sol->transportModel();
592}
593
594vector<string> SolutionArray::componentNames() const
595{
596 vector<string> components;
597 // leading auxiliary components
598 int pos = 0;
599 while (m_order->count(pos)) {
600 components.push_back(m_order->at(pos));
601 pos++;
602 }
603
604 // state information
605 auto phase = m_sol->thermo();
606 for (auto code : phase->nativeMode()) {
607 string name = string(1, code);
608 if (name == "X" || name == "Y") {
609 for (auto& spc : phase->speciesNames()) {
610 components.push_back(spc);
611 }
612 } else {
613 components.push_back(name);
614 }
615 }
616
617 // trailing auxiliary components
618 pos = -1;
619 while (m_order->count(pos)) {
620 components.push_back(m_order->at(pos));
621 pos--;
622 }
623
624 return components;
625}
626
627void SolutionArray::addExtra(const string& name, bool back)
628{
629 if (m_extra->count(name)) {
630 throw CanteraError("SolutionArray::addExtra",
631 "Component '{}' already exists.", name);
632 }
633 (*m_extra)[name] = AnyValue();
634 if (back) {
635 if (m_order->size()) {
636 // add name at end of back components
637 m_order->emplace(m_order->begin()->first - 1, name);
638 } else {
639 // first name after state components
640 m_order->emplace(-1, name);
641 }
642 } else {
643 if (m_order->size()) {
644 // insert name at end of front components
645 m_order->emplace(m_order->rbegin()->first + 1, name);
646 } else {
647 // name in leading position
648 m_order->emplace(0, name);
649 }
650 }
651}
652
653vector<string> SolutionArray::listExtra(bool all) const
654{
655 vector<string> names;
656 int pos = 0;
657 while (m_order->count(pos)) {
658 const auto& name = m_order->at(pos);
659 if (all || !m_extra->at(name).is<void>()) {
660 names.push_back(name);
661 }
662 pos++;
663 }
664
665 // trailing auxiliary components
666 pos = -1;
667 while (m_order->count(pos)) {
668 const auto& name = m_order->at(pos);
669 if (all || !m_extra->at(name).is<void>()) {
670 names.push_back(name);
671 }
672 pos--;
673 }
674 return names;
675}
676
677bool SolutionArray::hasComponent(const string& name, bool checkAlias) const
678{
679 if (m_extra->count(name)) {
680 // auxiliary data
681 return true;
682 }
683 if (m_sol->thermo()->speciesIndex(name) != npos) {
684 // species
685 return true;
686 }
687 if (name == "X" || name == "Y") {
688 // reserved names
689 return false;
690 }
691 if (checkAlias && reverseAliasMap.count(name)) {
692 // registered alias
693 return true;
694 }
695 // native state
696 return (m_sol->thermo()->nativeState().count(name));
697}
698
699AnyValue SolutionArray::getComponent(const string& name) const
700{
701 string _name = name;
702 if (!hasComponent(name, false) && reverseAliasMap.count(name)) {
703 _name = reverseAliasMap.at(name);
704 } else if (!hasComponent(name)) {
705 throw CanteraError("SolutionArray::getComponent",
706 "Unknown component '{}'.", name);
707 }
708
709 AnyValue out;
710 if (m_extra->count(_name)) {
711 // extra component
712 const auto& extra = m_extra->at(_name);
713 if (extra.is<void>()) {
714 return AnyValue();
715 }
716 if (m_size == m_dataSize) {
717 return extra; // slicing not necessary
718 }
719 if (extra.isVector<long int>()) {
720 return getSingle<long int>(extra, m_active);
721 }
722 if (extra.isVector<double>()) {
723 return getSingle<double>(extra, m_active);
724 }
725 if (extra.isVector<string>()) {
726 return getSingle<string>(extra, m_active);
727 }
728 if (extra.isVector<vector<double>>()) {
729 return getMulti<double>(extra, m_active);
730 }
731 if (extra.isVector<vector<long int>>()) {
732 return getMulti<long int>(extra, m_active);
733 }
734 if (extra.isVector<vector<string>>()) {
735 return getMulti<string>(extra, m_active);
736 }
737 throw NotImplementedError("SolutionArray::getComponent",
738 "Unable to get sliced data for component '{}' with type '{}'.",
739 name, extra.type_str());
740 }
741
742 // component is part of state information
743 vector<double> data(m_size);
744 size_t ix = m_sol->thermo()->speciesIndex(_name);
745 if (ix == npos) {
746 // state other than species
747 ix = m_sol->thermo()->nativeState()[_name];
748 } else {
749 // species information
750 ix += m_stride - m_sol->thermo()->nSpecies();
751 }
752 for (size_t k = 0; k < m_size; ++k) {
753 data[k] = (*m_data)[m_active[k] * m_stride + ix];
754 }
755 out = data;
756 return out;
757}
758
759bool isSimpleVector(const AnyValue& any) {
760 return any.isVector<double>() || any.isVector<long int>() ||
761 any.isVector<string>() || any.isVector<bool>() ||
762 any.isVector<vector<double>>() || any.isVector<vector<long int>>() ||
763 any.isVector<vector<string>>() || any.isVector<vector<bool>>();
764}
765
766void SolutionArray::setComponent(const string& name, const AnyValue& data)
767{
768 if (!hasComponent(name)) {
769 throw CanteraError("SolutionArray::setComponent",
770 "Unknown component '{}'.", name);
771 }
772 if (m_extra->count(name)) {
773 _setExtra(name, data);
774 return;
775 }
776
777 size_t size = data.vectorSize();
778 if (size == npos) {
779 throw CanteraError("SolutionArray::setComponent",
780 "Invalid type of component '{}': expected simple array type, "
781 "but received '{}'.", name, data.type_str());
782 }
783 if (size != m_size) {
784 throw CanteraError("SolutionArray::setComponent",
785 "Invalid size of component '{}': expected size {} but received {}.",
786 name, m_size, size);
787 }
788
789 auto& vec = data.asVector<double>();
790 size_t ix = m_sol->thermo()->speciesIndex(name);
791 if (ix == npos) {
792 ix = m_sol->thermo()->nativeState()[name];
793 } else {
794 ix += m_stride - m_sol->thermo()->nSpecies();
795 }
796 for (size_t k = 0; k < m_size; ++k) {
797 (*m_data)[m_active[k] * m_stride + ix] = vec[k];
798 }
799}
800
801void SolutionArray::setLoc(int loc, bool restore)
802{
803 size_t loc_ = static_cast<size_t>(loc);
804 if (m_size == 0) {
805 throw CanteraError("SolutionArray::setLoc",
806 "Unable to set location in empty SolutionArray.");
807 } else if (loc < 0) {
808 if (m_loc == npos) {
809 throw CanteraError("SolutionArray::setLoc",
810 "Both current and buffered indices are invalid.");
811 }
812 return;
813 } else if (static_cast<size_t>(m_active[loc_]) == m_loc) {
814 return;
815 } else if (loc_ >= m_size) {
816 throw IndexError("SolutionArray::setLoc", "indices", loc_, m_size);
817 }
818 m_loc = static_cast<size_t>(m_active[loc_]);
819 if (restore) {
820 size_t nState = m_sol->thermo()->stateSize();
821 m_sol->thermo()->restoreState(nState, m_data->data() + m_loc * m_stride);
822 }
823}
824
825void SolutionArray::updateState(int loc)
826{
827 setLoc(loc, false);
828 size_t nState = m_sol->thermo()->stateSize();
829 m_sol->thermo()->saveState(nState, m_data->data() + m_loc * m_stride);
830}
831
832vector<double> SolutionArray::getState(int loc)
833{
834 setLoc(loc);
835 size_t nState = m_sol->thermo()->stateSize();
836 vector<double> out(nState);
837 m_sol->thermo()->saveState(out); // thermo contains current state
838 return out;
839}
840
841void SolutionArray::setState(int loc, const vector<double>& state)
842{
843 size_t nState = m_sol->thermo()->stateSize();
844 if (state.size() != nState) {
845 throw CanteraError("SolutionArray::setState",
846 "Expected array to have length {}, but received an array of length {}.",
847 nState, state.size());
848 }
849 setLoc(loc, false);
850 m_sol->thermo()->restoreState(state);
851 m_sol->thermo()->saveState(nState, m_data->data() + m_loc * m_stride);
852}
853
854void SolutionArray::normalize() {
855 auto phase = m_sol->thermo();
856 auto nativeState = phase->nativeState();
857 if (nativeState.size() < 3) {
858 return;
859 }
860 size_t nState = phase->stateSize();
861 vector<double> out(nState);
862 if (nativeState.count("Y")) {
863 size_t offset = nativeState["Y"];
864 for (int loc = 0; loc < static_cast<int>(m_size); loc++) {
865 setLoc(loc, true); // set location and restore state
866 phase->setMassFractions(m_data->data() + m_loc * m_stride + offset);
867 m_sol->thermo()->saveState(out);
868 setState(loc, out);
869 }
870 } else if (nativeState.count("X")) {
871 size_t offset = nativeState["X"];
872 for (int loc = 0; loc < static_cast<int>(m_size); loc++) {
873 setLoc(loc, true); // set location and restore state
874 phase->setMoleFractions(m_data->data() + m_loc * m_stride + offset);
875 m_sol->thermo()->saveState(out);
876 setState(loc, out);
877 }
878 } else {
879 throw NotImplementedError("SolutionArray::normalize",
880 "Not implemented for mode '{}'.", phase->nativeMode());
881 }
882}
883
884AnyMap SolutionArray::getAuxiliary(int loc)
885{
886 setLoc(loc);
887 AnyMap out;
888 for (const auto& [key, extra] : *m_extra) {
889 if (extra.is<void>()) {
890 out[key] = extra;
891 } else if (extra.isVector<long int>()) {
892 out[key] = extra.asVector<long int>()[m_loc];
893 } else if (extra.isVector<double>()) {
894 out[key] = extra.asVector<double>()[m_loc];
895 } else if (extra.isVector<string>()) {
896 out[key] = extra.asVector<string>()[m_loc];
897 } else if (extra.isVector<vector<long int>>()) {
898 out[key] = extra.asVector<vector<long int>>()[m_loc];
899 } else if (extra.isVector<vector<double>>()) {
900 out[key] = extra.asVector<vector<double>>()[m_loc];
901 } else if (extra.isVector<vector<string>>()) {
902 out[key] = extra.asVector<vector<string>>()[m_loc];
903 } else {
904 throw NotImplementedError("SolutionArray::getAuxiliary",
905 "Unable to retrieve data for component '{}' with type '{}'.",
906 key, extra.type_str());
907 }
908 }
909 return out;
910}
911
912void SolutionArray::setAuxiliary(int loc, const AnyMap& data)
913{
914 setLoc(loc, false);
915 for (const auto& [name, value] : data) {
916 if (!m_extra->count(name)) {
917 throw CanteraError("SolutionArray::setAuxiliary",
918 "Unknown auxiliary component '{}'.", name);
919 }
920 auto& extra = m_extra->at(name);
921 if (extra.is<void>()) {
922 if (m_dataSize > 1) {
923 throw CanteraError("SolutionArray::setAuxiliary",
924 "Unable to set location for type '{}': "
925 "component is not initialized.", name);
926 }
927 _initExtra(name, value);
928 _resizeExtra(name);
929 }
930 try {
931 if (extra.isVector<long int>()) {
932 setAuxiliarySingle<long int>(m_loc, extra, value);
933 } else if (extra.isVector<double>()) {
934 setAuxiliarySingle<double>(m_loc, extra, value);
935 } else if (extra.isVector<string>()) {
936 setAuxiliarySingle<string>(m_loc, extra, value);
937 } else if (extra.isVector<vector<long int>>()) {
938 setAuxiliaryMulti<long int>(m_loc, extra, value);
939 } else if (extra.isVector<vector<double>>()) {
940 setAuxiliaryMulti<double>(m_loc, extra, value);
941 } else if (extra.isVector<vector<string>>()) {
942 setAuxiliaryMulti<string>(m_loc, extra, value);
943 } else {
944 throw CanteraError("SolutionArray::setAuxiliary",
945 "Unable to set entry for type '{}'.", extra.type_str());
946 }
947 } catch (CanteraError& err) {
948 // make failed type conversions traceable
949 throw CanteraError("SolutionArray::setAuxiliary",
950 "Encountered incompatible value for component '{}':\n{}",
951 name, err.getMessage());
952 }
953 }
954}
955
956AnyMap preamble(const string& desc)
957{
958 AnyMap data;
959 if (desc.size()) {
960 data["description"] = desc;
961 }
962 data["generator"] = "Cantera SolutionArray";
963 data["cantera-version"] = CANTERA_VERSION;
964 data["git-commit"] = gitCommit();
965
966 // Add a timestamp indicating the current time
967 time_t aclock;
968 ::time(&aclock); // Get time in seconds
969 struct tm* newtime = localtime(&aclock); // Convert time to struct tm form
970 data["date"] = stripnonprint(asctime(newtime));
971
972 return data;
973}
974
975AnyMap& openField(AnyMap& root, const string& name)
976{
977 if (!name.size()) {
978 return root;
979 }
980
981 // locate field based on 'name'
982 vector<string> tokens;
983 tokenizePath(name, tokens);
984 AnyMap* ptr = &root; // use raw pointer to avoid copying
985 string path = "";
986 for (auto& field : tokens) {
987 path += "/" + field;
988 AnyMap& sub = *ptr;
989 if (sub.hasKey(field) && !sub[field].is<AnyMap>()) {
990 throw CanteraError("openField",
991 "Encountered invalid existing field '{}'.", path);
992 } else if (!sub.hasKey(field)) {
993 sub[field] = AnyMap();
994 }
995 ptr = &sub[field].as<AnyMap>();
996 }
997 return *ptr;
998}
999
1000void SolutionArray::writeHeader(const string& fname, const string& name,
1001 const string& desc, bool overwrite)
1002{
1003 Storage file(fname, true);
1004 if (file.checkGroup(name, true)) {
1005 if (!overwrite) {
1006 throw CanteraError("SolutionArray::writeHeader",
1007 "Group name '{}' exists; use 'overwrite' argument to overwrite.", name);
1008 }
1009 file.deleteGroup(name);
1010 file.checkGroup(name, true);
1011 }
1012 file.writeAttributes(name, preamble(desc));
1013}
1014
1015void SolutionArray::writeHeader(AnyMap& root, const string& name,
1016 const string& desc, bool overwrite)
1017{
1018 AnyMap& data = openField(root, name);
1019 if (!data.empty()) {
1020 if (!overwrite) {
1021 throw CanteraError("SolutionArray::writeHeader",
1022 "Field name '{}' exists; use 'overwrite' argument to overwrite.", name);
1023 }
1024 data.clear();
1025 }
1026 data.update(preamble(desc));
1027}
1028
1029void SolutionArray::writeEntry(const string& fname, bool overwrite, const string& basis)
1030{
1031 if (apiNdim() != 1) {
1032 throw CanteraError("SolutionArray::writeEntry",
1033 "Tabular output of CSV data only works for 1D SolutionArray objects.");
1034 }
1035 set<string> speciesNames;
1036 for (const auto& species : m_sol->thermo()->speciesNames()) {
1037 speciesNames.insert(species);
1038 }
1039 bool mole;
1040 if (basis == "") {
1041 const auto& nativeState = m_sol->thermo()->nativeState();
1042 mole = nativeState.find("X") != nativeState.end();
1043 } else if (basis == "X" || basis == "mole") {
1044 mole = true;
1045 } else if (basis == "Y" || basis == "mass") {
1046 mole = false;
1047 } else {
1048 throw CanteraError("SolutionArray::writeEntry",
1049 "Invalid species basis '{}'.", basis);
1050 }
1051
1052 auto names = componentNames();
1053 size_t last = names.size() - 1;
1054 vector<AnyValue> components;
1055 vector<bool> isSpecies;
1056 fmt::memory_buffer header;
1057 for (const auto& key : names) {
1058 string label = key;
1059 size_t col;
1060 if (speciesNames.find(key) == speciesNames.end()) {
1061 // Pre-read component vectors
1062 isSpecies.push_back(false);
1063 components.emplace_back(getComponent(key));
1064 col = components.size() - 1;
1065 if (!components[col].isVector<double>() &&
1066 !components[col].isVector<long int>() &&
1067 !components[col].isVector<string>())
1068 {
1069 throw CanteraError("SolutionArray::writeEntry",
1070 "Multi-dimensional column '{}' is not supported for CSV output.",
1071 key);
1072 }
1073 } else {
1074 // Delay reading species data as base can be either mole or mass
1075 isSpecies.push_back(true);
1076 components.emplace_back(AnyValue());
1077 col = components.size() - 1;
1078 if (mole) {
1079 label = "X_" + label;
1080 } else {
1081 label = "Y_" + label;
1082 }
1083 }
1084 if (label.find("\"") != string::npos || label.find("\n") != string::npos) {
1085 throw NotImplementedError("SolutionArray::writeEntry",
1086 "Detected column name containing double quotes or line feeds: '{}'.",
1087 label);
1088 }
1089 string sep = (col == last) ? "" : ",";
1090 if (label.find(",") != string::npos) {
1091 fmt_append(header, "\"{}\"{}", label, sep);
1092 } else {
1093 fmt_append(header, "{}{}", label, sep);
1094 }
1095 }
1096
1097 // (Most) potential exceptions have been thrown; start writing data to file
1098 if (std::ifstream(fname).good()) {
1099 if (!overwrite) {
1100 throw CanteraError("SolutionArray::writeEntry",
1101 "File '{}' already exists; use option 'overwrite' to replace CSV file.",
1102 fname);
1103 }
1104 std::remove(fname.c_str());
1105 }
1106 std::ofstream output(fname);
1107 output << to_string(header) << std::endl;
1108
1109 size_t maxLen = npos;
1110 vector<double> buf(speciesNames.size(), 0.);
1111 for (int row = 0; row < static_cast<int>(m_size); row++) {
1112 fmt::memory_buffer line;
1113 if (maxLen != npos) {
1114 line.reserve(maxLen);
1115 }
1116 setLoc(row);
1117 if (mole) {
1118 m_sol->thermo()->getMoleFractions(buf.data());
1119 } else {
1120 m_sol->thermo()->getMassFractions(buf.data());
1121 }
1122
1123 size_t idx = 0;
1124 for (size_t col = 0; col < components.size(); col++) {
1125 string sep = (col == last) ? "" : ",";
1126 if (isSpecies[col]) {
1127 fmt_append(line, "{:.9g}{}", buf[idx++], sep);
1128 } else {
1129 auto& data = components[col];
1130 if (data.isVector<double>()) {
1131 fmt_append(line, "{:.9g}{}", data.asVector<double>()[row], sep);
1132 } else if (data.isVector<long int>()) {
1133 fmt_append(line, "{}{}", data.asVector<long int>()[row], sep);
1134 } else if (data.isVector<bool>()) {
1135 fmt_append(line, "{}{}",
1136 static_cast<bool>(data.asVector<bool>()[row]), sep);
1137 } else {
1138 auto value = data.asVector<string>()[row];
1139 if (value.find("\"") != string::npos ||
1140 value.find("\n") != string::npos)
1141 {
1142 throw NotImplementedError("SolutionArray::writeEntry",
1143 "Detected value containing double quotes or line feeds: "
1144 "'{}'", value);
1145 }
1146 if (value.find(",") != string::npos) {
1147 fmt_append(line, "\"{}\"{}", value, sep);
1148 } else {
1149 fmt_append(line, "{}{}", value, sep);
1150 }
1151 }
1152 }
1153 }
1154 output << to_string(line) << std::endl;
1155 maxLen = std::max(maxLen, line.size());
1156 }
1157}
1158
1159void SolutionArray::writeEntry(const string& fname, const string& name,
1160 const string& sub, bool overwrite, int compression)
1161{
1162 if (name == "") {
1163 throw CanteraError("SolutionArray::writeEntry",
1164 "Group name specifying root location must not be empty.");
1165 }
1166 if (m_size < m_dataSize) {
1167 throw NotImplementedError("SolutionArray::writeEntry",
1168 "Unable to save sliced data.");
1169 }
1170 Storage file(fname, true);
1171 if (compression) {
1172 file.setCompressionLevel(compression);
1173 }
1174 string path = name;
1175 if (sub != "") {
1176 path += "/" + sub;
1177 } else {
1178 path += "/data";
1179 }
1180 if (file.checkGroup(path, true)) {
1181 if (!overwrite) {
1182 throw CanteraError("SolutionArray::writeEntry",
1183 "Group name '{}' exists; use 'overwrite' argument to overwrite.", name);
1184 }
1185 file.deleteGroup(path);
1186 file.checkGroup(path, true);
1187 }
1188 file.writeAttributes(path, m_meta);
1189 AnyMap more;
1190 if (apiNdim() == 1) {
1191 more["size"] = int(m_dataSize);
1192 } else {
1193 more["api-shape"] = m_apiShape;
1194 }
1195 if (!m_meta.hasKey("transport-model") && m_sol->transport()) {
1196 more["transport-model"] = m_sol->transportModel();
1197 }
1198 more["components"] = componentNames();
1199 file.writeAttributes(path, more);
1200 if (!m_dataSize) {
1201 return;
1202 }
1203
1204 const auto& nativeState = m_sol->thermo()->nativeState();
1205 size_t nSpecies = m_sol->thermo()->nSpecies();
1206 for (auto& [key, offset] : nativeState) {
1207 if (key == "X" || key == "Y") {
1208 vector<vector<double>> prop;
1209 for (size_t i = 0; i < m_size; i++) {
1210 size_t first = offset + i * m_stride;
1211 prop.emplace_back(m_data->begin() + first,
1212 m_data->begin() + first + nSpecies);
1213 }
1214 AnyValue data;
1215 data = prop;
1216 file.writeData(path, key, data);
1217 } else {
1218 auto data = getComponent(key);
1219 file.writeData(path, key, data);
1220 }
1221 }
1222
1223 for (const auto& [key, value] : *m_extra) {
1224 if (isSimpleVector(value)) {
1225 file.writeData(path, key, value);
1226 } else if (value.is<void>()) {
1227 // skip unintialized component
1228 } else {
1229 throw NotImplementedError("SolutionArray::writeEntry",
1230 "Unable to save component '{}' with data type {}.",
1231 key, value.type_str());
1232 }
1233 }
1234}
1235
1236void SolutionArray::writeEntry(AnyMap& root, const string& name, const string& sub,
1237 bool overwrite)
1238{
1239 if (name == "") {
1240 throw CanteraError("SolutionArray::writeEntry",
1241 "Field name specifying root location must not be empty.");
1242 }
1243 if (m_size < m_dataSize) {
1244 throw NotImplementedError("SolutionArray::writeEntry",
1245 "Unable to save sliced data.");
1246 }
1247 string path = name;
1248 if (sub != "") {
1249 path += "/" + sub;
1250 } else {
1251 path += "/data";
1252 }
1253 AnyMap& data = openField(root, path);
1254 bool preexisting = !data.empty();
1255 if (preexisting && !overwrite) {
1256 throw CanteraError("SolutionArray::writeEntry",
1257 "Field name '{}' exists; use 'overwrite' argument to overwrite.", name);
1258 }
1259 if (apiNdim() == 1) {
1260 data["size"] = int(m_dataSize);
1261 } else {
1262 data["api-shape"] = m_apiShape;
1263 }
1264 if (m_sol->transport() && m_sol->transportModel() != "none") {
1265 data["transport-model"] = m_sol->transportModel();
1266 }
1267 data.update(m_meta);
1268
1269 if (m_size > 1) {
1270 data["components"] = componentNames();
1271 }
1272
1273 for (auto& [_, key] : *m_order) {
1274 data[key] = m_extra->at(key);
1275 }
1276
1277 auto phase = m_sol->thermo();
1278 if (m_size == 1) {
1279 setLoc(0);
1280 data["temperature"] = phase->temperature();
1281 data["pressure"] = phase->pressure();
1282 auto surf = std::dynamic_pointer_cast<SurfPhase>(phase);
1283 auto nSpecies = phase->nSpecies();
1284 vector<double> values(nSpecies);
1285 if (surf) {
1286 surf->getCoverages(&values[0]);
1287 } else {
1288 phase->getMassFractions(&values[0]);
1289 }
1290 AnyMap items;
1291 for (size_t k = 0; k < nSpecies; k++) {
1292 if (values[k] != 0.0) {
1293 items[phase->speciesName(k)] = values[k];
1294 }
1295 }
1296 if (surf) {
1297 data["coverages"] = std::move(items);
1298 } else {
1299 data["mass-fractions"] = std::move(items);
1300 }
1301 } else if (m_size > 1) {
1302 for (auto& code : phase->nativeMode()) {
1303 string name(1, code);
1304 if (name == "X" || name == "Y") {
1305 for (auto& spc : phase->speciesNames()) {
1306 data[spc] = getComponent(spc);
1307 }
1308 data["basis"] = name == "X" ? "mole" : "mass";
1309 } else {
1310 data[name] = getComponent(name);
1311 }
1312 }
1313 }
1314
1315 static bool reg = AnyMap::addOrderingRules("SolutionArray",
1316 {{"head", "type"}, {"head", "size"}, {"head", "basis"}});
1317 if (reg) {
1318 data["__type__"] = "SolutionArray";
1319 }
1320
1321 // If this is not replacing an existing solution, put it at the end
1322 if (!preexisting) {
1323 data.setLoc(INT_MAX, 0);
1324 }
1325}
1326
1327void SolutionArray::append(const vector<double>& state, const AnyMap& extra)
1328{
1329 if (apiNdim() > 1) {
1330 throw NotImplementedError("SolutionArray::append",
1331 "Unable to append multi-dimensional arrays.");
1332 }
1333
1334 int pos = size();
1335 resize(pos + 1);
1336 try {
1337 setState(pos, state);
1338 setAuxiliary(pos, extra);
1339 } catch (CanteraError& err) {
1340 // undo resize and rethrow error
1341 resize(pos);
1342 throw CanteraError("SolutionArray::append", err.getMessage());
1343 }
1344}
1345
1346void SolutionArray::save(const string& fname, const string& name, const string& sub,
1347 const string& desc, bool overwrite, int compression,
1348 const string& basis)
1349{
1350 if (m_size < m_dataSize) {
1351 throw NotImplementedError("SolutionArray::save",
1352 "Unable to save sliced data.");
1353 }
1354 size_t dot = fname.find_last_of(".");
1355 string extension = (dot != npos) ? toLowerCopy(fname.substr(dot + 1)) : "";
1356 if (extension == "csv") {
1357 if (name != "") {
1358 warn_user("SolutionArray::save",
1359 "Parameter 'name' not used for CSV output.");
1360 }
1361 writeEntry(fname, overwrite, basis);
1362 return;
1363 }
1364 if (basis != "") {
1365 warn_user("SolutionArray::save",
1366 "Argument 'basis' is not used for HDF or YAML output.", basis);
1367 }
1368 if (extension == "h5" || extension == "hdf" || extension == "hdf5") {
1369 writeHeader(fname, name, desc, overwrite);
1370 writeEntry(fname, name, sub, true, compression);
1371 return;
1372 }
1373 if (extension == "yaml" || extension == "yml") {
1374 // Check for an existing file and load it if present
1375 AnyMap data;
1376 std::ifstream file(fname);
1377 if (file.good() && file.peek() != std::ifstream::traits_type::eof()) {
1378 data = AnyMap::fromYamlFile(fname);
1379 }
1380 writeHeader(data, name, desc, overwrite);
1381 writeEntry(data, name, sub, true);
1382
1383 // Write the output file and remove the now-outdated cached file
1384 std::ofstream out(fname);
1385 out << data.toYamlString();
1386 AnyMap::clearCachedFile(fname);
1387 return;
1388 }
1389 throw CanteraError("SolutionArray::save",
1390 "Unknown file extension '{}'.", extension);
1391}
1392
1393AnyMap SolutionArray::readHeader(const string& fname, const string& name)
1394{
1395 Storage file(fname, false);
1396 file.checkGroup(name);
1397 return file.readAttributes(name, false);
1398}
1399
1400const AnyMap& locateField(const AnyMap& root, const string& name)
1401{
1402 if (!name.size()) {
1403 return root;
1404 }
1405
1406 // locate field based on 'name'
1407 vector<string> tokens;
1408 tokenizePath(name, tokens);
1409 const AnyMap* ptr = &root; // use raw pointer to avoid copying
1410 string path = "";
1411 for (auto& field : tokens) {
1412 path += "/" + field;
1413 const AnyMap& sub = *ptr;
1414 if (!sub.hasKey(field) || !sub[field].is<AnyMap>()) {
1415 throw CanteraError("SolutionArray::locateField",
1416 "No field or solution with name '{}'.", path);
1417 }
1418 ptr = &sub[field].as<AnyMap>();
1419 }
1420 return *ptr;
1421}
1422
1423AnyMap SolutionArray::readHeader(const AnyMap& root, const string& name)
1424{
1425 auto sub = locateField(root, name);
1426 AnyMap header;
1427 for (const auto& [key, value] : sub) {
1428 if (!sub[key].is<AnyMap>()) {
1429 header[key] = value;
1430 }
1431 }
1432 return header;
1433}
1434
1435AnyMap SolutionArray::restore(const string& fname,
1436 const string& name, const string& sub)
1437{
1438 size_t dot = fname.find_last_of(".");
1439 string extension = (dot != npos) ? toLowerCopy(fname.substr(dot + 1)) : "";
1440 AnyMap header;
1441 if (extension == "csv") {
1442 throw NotImplementedError("SolutionArray::restore",
1443 "CSV import not implemented; if using Python, data can be imported via "
1444 "'read_csv' instead.");
1445 }
1446 if (extension == "h5" || extension == "hdf" || extension == "hdf5") {
1447 readEntry(fname, name, sub);
1448 header = readHeader(fname, name);
1449 } else if (extension == "yaml" || extension == "yml") {
1450 const AnyMap& root = AnyMap::fromYamlFile(fname);
1451 readEntry(root, name, sub);
1452 header = readHeader(root, name);
1453 } else {
1454 throw CanteraError("SolutionArray::restore",
1455 "Unknown file extension '{}'; supported extensions include "
1456 "'h5'/'hdf'/'hdf5' and 'yml'/'yaml'.", extension);
1457 }
1458 return header;
1459}
1460
1461void SolutionArray::_initExtra(const string& name, const AnyValue& value)
1462{
1463 if (!m_extra->count(name)) {
1464 throw CanteraError("SolutionArray::_initExtra",
1465 "Component '{}' does not exist.", name);
1466 }
1467 auto& extra = (*m_extra)[name];
1468 if (!extra.is<void>()) {
1469 throw CanteraError("SolutionArray::_initExtra",
1470 "Component '{}' is already initialized.", name);
1471 }
1472 try {
1473 if (value.is<long int>()) {
1474 extra = vector<long int>(m_dataSize, value.as<long int>());
1475 } else if (value.is<double>()) {
1476 extra = vector<double>(m_dataSize, value.as<double>());
1477 } else if (value.is<string>()) {
1478 extra = vector<string>(m_dataSize, value.as<string>());
1479 } else if (value.isVector<long int>()) {
1480 extra = vector<vector<long int>>(m_dataSize, value.asVector<long int>());
1481 } else if (value.isVector<double>()) {
1482 extra = vector<vector<double>>(m_dataSize, value.asVector<double>());
1483 } else if (value.isVector<string>()) {
1484 extra = vector<vector<string>>(m_dataSize, value.asVector<string>());
1485 } else if (value.is<void>()) {
1486 // placeholder for unknown type; settable in setComponent
1487 extra = AnyValue();
1488 } else {
1489 throw NotImplementedError("SolutionArray::_initExtra",
1490 "Unable to initialize component '{}' with type '{}'.",
1491 name, value.type_str());
1492 }
1493 } catch (CanteraError& err) {
1494 // make failed type conversions traceable
1495 throw CanteraError("SolutionArray::_initExtra",
1496 "Encountered incompatible value for initializing component '{}':\n{}",
1497 name, err.getMessage());
1498 }
1499}
1500
1501void SolutionArray::_resizeExtra(const string& name, const AnyValue& value)
1502{
1503 if (!m_extra->count(name)) {
1504 throw CanteraError("SolutionArray::_resizeExtra",
1505 "Component '{}' does not exist.", name);
1506 }
1507 auto& extra = (*m_extra)[name];
1508 if (extra.is<void>()) {
1509 // cannot resize placeholder (uninitialized component)
1510 return;
1511 }
1512
1513 try {
1514 if (extra.isVector<long int>()) {
1515 resizeSingle<long int>(extra, m_dataSize, value);
1516 } else if (extra.isVector<double>()) {
1517 resizeSingle<double>(extra, m_dataSize, value);
1518 } else if (extra.isVector<string>()) {
1519 resizeSingle<string>(extra, m_dataSize, value);
1520 } else if (extra.isVector<vector<double>>()) {
1521 resizeMulti<double>(extra, m_dataSize, value);
1522 } else if (extra.isVector<vector<long int>>()) {
1523 resizeMulti<long int>(extra, m_dataSize, value);
1524 } else if (extra.isVector<vector<string>>()) {
1525 resizeMulti<string>(extra, m_dataSize, value);
1526 } else {
1527 throw NotImplementedError("SolutionArray::_resizeExtra",
1528 "Unable to resize using type '{}'.", extra.type_str());
1529 }
1530 } catch (CanteraError& err) {
1531 // make failed type conversions traceable
1532 throw CanteraError("SolutionArray::_resizeExtra",
1533 "Encountered incompatible value for resizing component '{}':\n{}",
1534 name, err.getMessage());
1535 }
1536}
1537
1538void SolutionArray::_setExtra(const string& name, const AnyValue& data)
1539{
1540 if (!m_extra->count(name)) {
1541 throw CanteraError("SolutionArray::_setExtra",
1542 "Extra component '{}' does not exist.", name);
1543 }
1544
1545 auto& extra = m_extra->at(name);
1546 if (data.is<void>() && m_size == m_dataSize) {
1547 // reset placeholder
1548 extra = AnyValue();
1549 return;
1550 }
1551
1552 // uninitialized component
1553 if (extra.is<void>()) {
1554 if (m_size != m_dataSize) {
1555 throw CanteraError("SolutionArray::_setExtra",
1556 "Unable to replace '{}' for sliced data.", name);
1557 }
1558 if (data.vectorSize() == m_dataSize || data.matrixShape().first == m_dataSize) {
1559 // assign entire component
1560 extra = data;
1561 return;
1562 }
1563 // data contains default value
1564 if (m_dataSize == 0 && (data.isScalar() || data.vectorSize() > 0)) {
1565 throw CanteraError("SolutionArray::_setExtra",
1566 "Unable to initialize '{}' with non-empty values when SolutionArray is "
1567 "empty.", name);
1568 }
1569 if (m_dataSize && data.vectorSize() == 0) {
1570 throw CanteraError("SolutionArray::_setExtra",
1571 "Unable to initialize '{}' with empty array when SolutionArray is not "
1572 "empty." , name);
1573 }
1574 _initExtra(name, data);
1575 _resizeExtra(name, data);
1576 return;
1577 }
1578
1579 if (data.is<long int>()) {
1580 setScalar<long int>(extra, data, m_active);
1581 } else if (data.is<double>()) {
1582 setScalar<double>(extra, data, m_active);
1583 } else if (data.is<string>()) {
1584 setScalar<string>(extra, data, m_active);
1585 } else if (data.isVector<long int>()) {
1586 setSingle<long int>(extra, data, m_active);
1587 } else if (data.isVector<double>()) {
1588 setSingle<double>(extra, data, m_active);
1589 } else if (data.isVector<string>()) {
1590 setSingle<string>(extra, data, m_active);
1591 } else if (data.isVector<vector<long int>>()) {
1592 setMulti<long int>(extra, data, m_active);
1593 } else if (data.isVector<vector<double>>()) {
1594 setMulti<double>(extra, data, m_active);
1595 } else if (data.isVector<vector<string>>()) {
1596 setMulti<string>(extra, data, m_active);
1597 } else {
1598 throw NotImplementedError("SolutionArray::_setExtra",
1599 "Unable to set sliced data for component '{}' with type '{}'.",
1600 name, data.type_str());
1601 }
1602}
1603
1604string SolutionArray::_detectMode(const set<string>& names, bool native)
1605{
1606 // check set of available names against state acronyms defined by Phase::fullStates
1607 string mode = "";
1608 const auto& nativeState = m_sol->thermo()->nativeState();
1609 bool usesNativeState = false;
1610 auto surf = std::dynamic_pointer_cast<SurfPhase>(m_sol->thermo());
1611 bool found;
1612 for (const auto& item : m_sol->thermo()->fullStates()) {
1613 found = true;
1614 string name;
1615 usesNativeState = true;
1616 for (size_t i = 0; i < item.size(); i++) {
1617 // pick i-th letter from "full" state acronym
1618 name = string(1, item[i]);
1619 if (surf && (name == "X" || name == "Y")) {
1620 // override native state to enable detection of surface phases
1621 name = "C";
1622 usesNativeState = false;
1623 break;
1624 }
1625 if (names.count(name)) {
1626 // property is stored using letter acronym
1627 usesNativeState &= nativeState.count(name) > 0;
1628 } else if (aliasMap.count(name) && names.count(aliasMap.at(name))) {
1629 // property is stored using property name
1630 usesNativeState &= nativeState.count(name) > 0;
1631 } else {
1632 found = false;
1633 break;
1634 }
1635 }
1636 if (found) {
1637 mode = (name == "C") ? item.substr(0, 2) + "C" : item;
1638 break;
1639 }
1640 }
1641 if (!found) {
1642 if (surf && names.count("T") && names.count("X") && names.count("density")) {
1643 // Legacy HDF format erroneously uses density (Cantera < 3.0)
1644 return "legacySurf";
1645 }
1646 if (names.count("mass-flux") && names.count("mass-fractions")) {
1647 // Legacy YAML format stores incomplete state information (Cantera < 3.0)
1648 return "legacyInlet";
1649 }
1650 throw CanteraError("SolutionArray::_detectMode",
1651 "Detected incomplete thermodynamic information. Full states for a '{}' "
1652 "phase\nmay be defined by the following modes:\n'{}'\n"
1653 "Available data are: '{}'", m_sol->thermo()->type(),
1654 ba::join(m_sol->thermo()->fullStates(), "', '"), ba::join(names, "', '"));
1655 }
1656 if (usesNativeState && native) {
1657 return "native";
1658 }
1659 return mode;
1660}
1661
1662set<string> SolutionArray::_stateProperties(
1663 const string& mode, bool alias)
1664{
1665 set<string> states;
1666 if (mode == "native") {
1667 for (const auto& [name, offset] : m_sol->thermo()->nativeState()) {
1668 states.insert(alias ? aliasMap.at(name) : name);
1669 }
1670 } else {
1671 for (const auto& m : mode) {
1672 const string name = string(1, m);
1673 states.insert(alias ? aliasMap.at(name) : name);
1674 }
1675 }
1676
1677 return states;
1678}
1679
1680string getName(const set<string>& names, const string& name)
1681{
1682 if (names.count(name)) {
1683 return name;
1684 }
1685 const auto& alias = aliasMap.at(name);
1686 if (names.count(alias)) {
1687 return alias;
1688 }
1689 return name; // let exception be thrown elsewhere
1690}
1691
1692void SolutionArray::readEntry(const string& fname, const string& name,
1693 const string& sub)
1694{
1695 Storage file(fname, false);
1696 if (name == "") {
1697 throw CanteraError("SolutionArray::readEntry",
1698 "Group name specifying root location must not be empty.");
1699 }
1700 string path = name;
1701 if (sub != "" && file.checkGroup(name + "/" + sub, true)) {
1702 path += "/" + sub;
1703 } else if (sub == "" && file.checkGroup(name + "/data", true)) {
1704 // default data location
1705 path += "/data";
1706 }
1707 if (!file.checkGroup(path)) {
1708 throw CanteraError("SolutionArray::readEntry",
1709 "Group name specifying data entry is empty.");
1710 }
1711 m_extra->clear();
1712 auto [size, names] = file.contents(path);
1713 m_meta = file.readAttributes(path, true);
1714 if (m_meta.hasKey("size")) {
1715 // one-dimensional array
1716 resize(m_meta["size"].as<long int>());
1717 m_meta.erase("size");
1718 } else if (m_meta.hasKey("api-shape")) {
1719 // API uses multiple dimensions to interpret C++ SolutionArray
1720 setApiShape(m_meta["api-shape"].asVector<long int>());
1721 m_meta.erase("api-shape");
1722 } else {
1723 // legacy format; size needs to be detected
1724 resize(static_cast<int>(size));
1725 }
1726
1727 if (m_size == 0) {
1728 return;
1729 }
1730
1731 // determine storage mode of state data
1732 string mode = _detectMode(names);
1733 set<string> states = _stateProperties(mode);
1734 if (states.count("C")) {
1735 if (names.count("X")) {
1736 states.erase("C");
1737 states.insert("X");
1738 mode = mode.substr(0, 2) + "X";
1739 } else if (names.count("Y")) {
1740 states.erase("C");
1741 states.insert("Y");
1742 mode = mode.substr(0, 2) + "Y";
1743 }
1744 }
1745
1746 // restore state data
1747 size_t nSpecies = m_sol->thermo()->nSpecies();
1748 size_t nState = m_sol->thermo()->stateSize();
1749 const auto& nativeStates = m_sol->thermo()->nativeState();
1750 if (mode == "native") {
1751 // native state can be written directly into data storage
1752 for (const auto& [name, offset] : nativeStates) {
1753 if (name == "X" || name == "Y") {
1754 AnyValue data;
1755 data = file.readData(path, name, m_size, nSpecies);
1756 auto prop = data.asVector<vector<double>>();
1757 for (size_t i = 0; i < m_dataSize; i++) {
1758 std::copy(prop[i].begin(), prop[i].end(),
1759 m_data->data() + offset + i * m_stride);
1760 }
1761 } else {
1762 AnyValue data;
1763 data = file.readData(path, getName(names, name), m_dataSize, 0);
1764 setComponent(name, data);
1765 }
1766 }
1767 } else if (mode == "TPX") {
1768 AnyValue data;
1769 data = file.readData(path, getName(names, "T"), m_dataSize, 0);
1770 vector<double> T = std::move(data.asVector<double>());
1771 data = file.readData(path, getName(names, "P"), m_dataSize, 0);
1772 vector<double> P = std::move(data.asVector<double>());
1773 data = file.readData(path, "X", m_dataSize, nSpecies);
1774 vector<vector<double>> X = std::move(data.asVector<vector<double>>());
1775 for (size_t i = 0; i < m_dataSize; i++) {
1776 m_sol->thermo()->setMoleFractions_NoNorm(X[i].data());
1777 m_sol->thermo()->setState_TP(T[i], P[i]);
1778 m_sol->thermo()->saveState(nState, m_data->data() + i * m_stride);
1779 }
1780 } else if (mode == "TDX") {
1781 AnyValue data;
1782 data = file.readData(path, getName(names, "T"), m_dataSize, 0);
1783 vector<double> T = std::move(data.asVector<double>());
1784 data = file.readData(path, getName(names, "D"), m_dataSize, 0);
1785 vector<double> D = std::move(data.asVector<double>());
1786 data = file.readData(path, "X", m_dataSize, nSpecies);
1787 vector<vector<double>> X = std::move(data.asVector<vector<double>>());
1788 for (size_t i = 0; i < m_dataSize; i++) {
1789 m_sol->thermo()->setMoleFractions_NoNorm(X[i].data());
1790 m_sol->thermo()->setState_TD(T[i], D[i]);
1791 m_sol->thermo()->saveState(nState, m_data->data() + i * m_stride);
1792 }
1793 } else if (mode == "TPY") {
1794 AnyValue data;
1795 data = file.readData(path, getName(names, "T"), m_dataSize, 0);
1796 vector<double> T = std::move(data.asVector<double>());
1797 data = file.readData(path, getName(names, "P"), m_dataSize, 0);
1798 vector<double> P = std::move(data.asVector<double>());
1799 data = file.readData(path, "Y", m_dataSize, nSpecies);
1800 vector<vector<double>> Y = std::move(data.asVector<vector<double>>());
1801 for (size_t i = 0; i < m_dataSize; i++) {
1802 m_sol->thermo()->setMassFractions_NoNorm(Y[i].data());
1803 m_sol->thermo()->setState_TP(T[i], P[i]);
1804 m_sol->thermo()->saveState(nState, m_data->data() + i * m_stride);
1805 }
1806 } else if (mode == "legacySurf") {
1807 // erroneous TDX mode (should be TPX or TPY) - Sim1D (Cantera 2.5)
1808 AnyValue data;
1809 data = file.readData(path, getName(names, "T"), m_dataSize, 0);
1810 vector<double> T = std::move(data.asVector<double>());
1811 data = file.readData(path, "X", m_dataSize, nSpecies);
1812 vector<vector<double>> X = std::move(data.asVector<vector<double>>());
1813 for (size_t i = 0; i < m_dataSize; i++) {
1814 m_sol->thermo()->setMoleFractions_NoNorm(X[i].data());
1815 m_sol->thermo()->setTemperature(T[i]);
1816 m_sol->thermo()->saveState(nState, m_data->data() + i * m_stride);
1817 }
1818 warn_user("SolutionArray::readEntry",
1819 "Detected legacy HDF format with incomplete state information\nfor name "
1820 "'{}' (pressure missing).", path);
1821 } else if (mode == "") {
1822 throw CanteraError("SolutionArray::readEntry",
1823 "Data are not consistent with full state modes.");
1824 } else {
1825 throw NotImplementedError("SolutionArray::readEntry",
1826 "Import of '{}' data is not supported.", mode);
1827 }
1828
1829 // restore remaining data
1830 if (m_meta.hasKey("components")) {
1831 const auto& components = m_meta["components"].asVector<string>();
1832 bool back = false;
1833 for (const auto& name : components) {
1834 if (hasComponent(name, false) || name == "X" || name == "Y") {
1835 back = true;
1836 } else {
1837 auto _name = name;
1838 if (reverseAliasMap.count(name)) {
1839 _name = reverseAliasMap.at(name);
1840 }
1841 addExtra(_name, back);
1842 AnyValue data;
1843 data = file.readData(path, name, m_dataSize);
1844 setComponent(_name, data);
1845 }
1846 }
1847 m_meta.erase("components");
1848 } else {
1849 // data format used by Python h5py export (Cantera 2.5)
1850 warn_user("SolutionArray::readEntry", "Detected legacy HDF format.");
1851 for (const auto& name : names) {
1852 if (!hasComponent(name, false) && name != "X" && name != "Y") {
1853 auto _name = name;
1854 if (reverseAliasMap.count(name)) {
1855 _name = reverseAliasMap.at(name);
1856 }
1857 addExtra(_name);
1858 AnyValue data;
1859 data = file.readData(path, name, m_dataSize);
1860 setComponent(_name, data);
1861 }
1862 }
1863 }
1864
1865 if (m_meta.hasKey("transport-model")) {
1866 m_sol->setTransportModel(m_meta["transport-model"].asString());
1867 }
1868}
1869
1870void SolutionArray::readEntry(const AnyMap& root, const string& name, const string& sub)
1871{
1872 if (name == "") {
1873 throw CanteraError("SolutionArray::readEntry",
1874 "Field name specifying root location must not be empty.");
1875 }
1876 auto path = locateField(root, name);
1877 if (path.hasKey("generator") && sub != "") {
1878 // read entry from subfolder (since Cantera 3.0)
1879 path = locateField(root, name + "/" + sub);
1880 } else if (sub == "" && path.hasKey("data")) {
1881 // default data location
1882 path = locateField(root, name + "/data");
1883 }
1884
1885 // set size and initialize
1886 long size = 0;
1887 if (path.hasKey("size")) {
1888 // one-dimensional array
1889 resize(path["size"].asInt());
1890 } else if (path.hasKey("api-shape")) {
1891 // API uses multiple dimensions to interpret C++ SolutionArray
1892 auto& shape = path["api-shape"].asVector<long int>();
1893 setApiShape(shape);
1894 } else {
1895 // legacy format (Cantera 2.6)
1896 size = path.getInt("points", 0);
1897 if (!path.hasKey("T") && !path.hasKey("temperature")) {
1898 // overwrite size - Sim1D erroneously assigns '1'
1899 size = (long)0;
1900 }
1901 resize(static_cast<int>(size));
1902 }
1903 m_extra->clear();
1904
1905 // restore data
1906 set<string> exclude = {"size", "api-shape", "points", "X", "Y"};
1907 set<string> names = path.keys();
1908 size_t nState = m_sol->thermo()->stateSize();
1909 if (m_dataSize == 0) {
1910 // no data points
1911 } else if (m_dataSize == 1) {
1912 // single data point
1913 string mode = _detectMode(names, false);
1914 if (mode == "TPY") {
1915 double T = path[getName(names, "T")].asDouble();
1916 double P = path[getName(names, "P")].asDouble();
1917 auto Y = path["mass-fractions"].asMap<double>();
1918 m_sol->thermo()->setState_TPY(T, P, Y);
1919 } else if (mode == "TPC") {
1920 auto surf = std::dynamic_pointer_cast<SurfPhase>(m_sol->thermo());
1921 double T = path[getName(names, "T")].asDouble();
1922 double P = path["pressure"].asDouble();
1923 m_sol->thermo()->setState_TP(T, P);
1924 auto cov = path["coverages"].asMap<double>();
1925 surf->setCoveragesByName(cov);
1926 } else if (mode == "legacyInlet") {
1927 // missing property - Sim1D (Cantera 2.6)
1928 mode = "TPY";
1929 double T = path[getName(names, "T")].asDouble();
1930 auto Y = path["mass-fractions"].asMap<double>();
1931 m_sol->thermo()->setState_TPY(T, m_sol->thermo()->pressure(), Y);
1932 warn_user("SolutionArray::readEntry",
1933 "Detected legacy YAML format with incomplete state information\n"
1934 "for name '{}' (pressure missing).", name + "/" + sub);
1935 } else if (mode == "") {
1936 throw CanteraError("SolutionArray::readEntry",
1937 "Data are not consistent with full state modes.");
1938 } else {
1939 throw NotImplementedError("SolutionArray::readEntry",
1940 "Import of '{}' data is not supported.", mode);
1941 }
1942 m_sol->thermo()->saveState(nState, m_data->data());
1943 auto props = _stateProperties(mode, true);
1944 exclude.insert(props.begin(), props.end());
1945 } else {
1946 // multiple data points
1947 if (path.hasKey("components")) {
1948 const auto& components = path["components"].asVector<string>();
1949 bool back = false;
1950 for (const auto& name : components) {
1951 auto _name = name;
1952 if (hasComponent(name, false)) {
1953 back = true;
1954 } else {
1955 if (reverseAliasMap.count(name)) {
1956 _name = reverseAliasMap.at(name);
1957 }
1958 addExtra(_name, back);
1959 }
1960 setComponent(_name, path[name]);
1961 exclude.insert(name);
1962 }
1963 } else {
1964 // legacy YAML format does not provide for list of components
1965 for (const auto& [name, value] : path) {
1966 if (value.isVector<double>()) {
1967 const vector<double>& data = value.asVector<double>();
1968 if (data.size() == m_dataSize) {
1969 auto _name = name;
1970 if (!hasComponent(name, false)) {
1971 if (reverseAliasMap.count(name)) {
1972 _name = reverseAliasMap.at(name);
1973 }
1974 addExtra(_name);
1975 }
1976 setComponent(_name, value);
1977 exclude.insert(name);
1978 }
1979 }
1980 }
1981 }
1982
1983 // check that state data are complete
1984 const auto& nativeState = m_sol->thermo()->nativeState();
1985 set<string> props;
1986 set<string> missingProps;
1987 for (const auto& [name, offset] : nativeState) {
1988 if (exclude.count(name)) {
1989 props.insert(name);
1990 } else {
1991 missingProps.insert(name);
1992 }
1993 }
1994
1995 set<string> TY = {"T", "Y"};
1996 if (props == TY && missingProps.count("D") && path.hasKey("pressure")) {
1997 // missing "D" - Sim1D (Cantera 2.6)
1998 double P = path["pressure"].asDouble();
1999 const size_t offset_T = nativeState.find("T")->second;
2000 const size_t offset_D = nativeState.find("D")->second;
2001 const size_t offset_Y = nativeState.find("Y")->second;
2002 for (size_t i = 0; i < m_dataSize; i++) {
2003 double T = (*m_data)[offset_T + i * m_stride];
2004 m_sol->thermo()->setState_TPY(
2005 T, P, m_data->data() + offset_Y + i * m_stride);
2006 (*m_data)[offset_D + i * m_stride] = m_sol->thermo()->density();
2007 }
2008 } else if (missingProps.size()) {
2009 throw CanteraError("SolutionArray::readEntry",
2010 "Incomplete state information: missing '{}'.",
2011 ba::join(missingProps, "', '"));
2012 }
2013 }
2014
2015 // add meta data
2016 for (const auto& [name, value] : path) {
2017 if (!exclude.count(name)) {
2018 m_meta[name] = value;
2019 }
2020 }
2021
2022 if (m_meta.hasKey("transport-model")) {
2023 m_sol->setTransportModel(m_meta["transport-model"].asString());
2024 }
2025}
2026
2027namespace { // restrict scope of helper functions to local translation unit
2028
2029template<class T>
2030AnyValue getSingle(const AnyValue& extra, const vector<int>& slice)
2031{
2032 vector<T> data(slice.size());
2033 const auto& vec = extra.asVector<T>();
2034 for (size_t k = 0; k < slice.size(); ++k) {
2035 data[k] = vec[slice[k]];
2036 }
2037 AnyValue out;
2038 out = data;
2039 return out;
2040}
2041
2042template<class T>
2043AnyValue getMulti(const AnyValue& extra, const vector<int>& slice)
2044{
2045 vector<vector<T>> data(slice.size());
2046 const auto& vec = extra.asVector<vector<T>>();
2047 for (size_t k = 0; k < slice.size(); ++k) {
2048 data[k] = vec[slice[k]];
2049 }
2050 AnyValue out;
2051 out = data;
2052 return out;
2053}
2054
2055template<class T>
2056void setScalar(AnyValue& extra, const AnyValue& data, const vector<int>& slice)
2057{
2058 T value = data.as<T>();
2059 if (extra.isVector<T>()) {
2060 auto& vec = extra.asVector<T>();
2061 for (size_t k = 0; k < slice.size(); ++k) {
2062 vec[slice[k]] = value;
2063 }
2064 } else {
2065 throw CanteraError("SolutionArray::setScalar",
2066 "Incompatible input data: unable to assign '{}' data to '{}'.",
2067 data.type_str(), extra.type_str());
2068 }
2069}
2070
2071template<class T>
2072void setSingle(AnyValue& extra, const AnyValue& data, const vector<int>& slice)
2073{
2074 size_t size = slice.size();
2075 if (extra.vectorSize() == size && data.vectorSize() == size) {
2076 extra = data; // no slicing necessary; type can change
2077 return;
2078 }
2079 if (extra.matrixShape().first == size && data.vectorSize() == size) {
2080 extra = data; // no slicing necessary; shape and type can change
2081 return;
2082 }
2083 if (extra.type_str() != data.type_str()) {
2084 // do not allow changing of data type when slicing
2085 throw CanteraError("SolutionArray::setSingle",
2086 "Incompatible types: expected '{}' but received '{}'.",
2087 extra.type_str(), data.type_str());
2088 }
2089 const auto& vData = data.asVector<T>();
2090 if (vData.size() != size) {
2091 throw CanteraError("SolutionArray::setSingle",
2092 "Invalid input data size: expected {} entries but received {}.",
2093 size, vData.size());
2094 }
2095 auto& vec = extra.asVector<T>();
2096 for (size_t k = 0; k < size; ++k) {
2097 vec[slice[k]] = vData[k];
2098 }
2099}
2100
2101template<class T>
2102void setMulti(AnyValue& extra, const AnyValue& data, const vector<int>& slice)
2103{
2104 if (!data.isMatrix<T>()) {
2105 throw CanteraError("SolutionArray::setMulti",
2106 "Invalid input data shape: inconsistent number of columns.");
2107 }
2108 size_t size = slice.size();
2109 auto [rows, cols] = data.matrixShape();
2110 if (extra.matrixShape().first == size && rows == size) {
2111 extra = data; // no slicing necessary; type can change
2112 return;
2113 }
2114 if (extra.vectorSize() == size && rows == size) {
2115 extra = data; // no slicing necessary; shape and type can change
2116 return;
2117 }
2118 if (extra.type_str() != data.type_str()) {
2119 // do not allow changing of data type when slicing
2120 throw CanteraError("SolutionArray::setMulti",
2121 "Incompatible types: expected '{}' but received '{}'.",
2122 extra.type_str(), data.type_str());
2123 }
2124 if (rows != size) {
2125 throw CanteraError("SolutionArray::setMulti",
2126 "Invalid input data shape: expected {} rows but received {}.",
2127 size, rows);
2128 }
2129 if (extra.matrixShape().second != cols) {
2130 throw CanteraError("SolutionArray::setMulti",
2131 "Invalid input data shape: expected {} columns but received {}.",
2132 extra.matrixShape().second, cols);
2133 }
2134 const auto& vData = data.asVector<vector<T>>();
2135 auto& vec = extra.asVector<vector<T>>();
2136 for (size_t k = 0; k < slice.size(); ++k) {
2137 vec[slice[k]] = vData[k];
2138 }
2139}
2140
2141template<class T>
2142void resizeSingle(AnyValue& extra, size_t size, const AnyValue& value)
2143{
2144 T defaultValue;
2145 if (value.is<void>()) {
2146 defaultValue = vector<T>(1)[0];
2147 } else {
2148 defaultValue = value.as<T>();
2149 }
2150 extra.asVector<T>().resize(size, defaultValue);
2151}
2152
2153template<class T>
2154void resizeMulti(AnyValue& extra, size_t size, const AnyValue& value)
2155{
2156 vector<T> defaultValue;
2157 if (value.is<void>()) {
2158 defaultValue = vector<T>(extra.matrixShape().second);
2159 } else {
2160 defaultValue = value.as<vector<T>>();
2161 }
2162 extra.asVector<vector<T>>().resize(size, defaultValue);
2163}
2164
2165template<class T>
2166void resetSingle(AnyValue& extra, const vector<int>& slice)
2167{
2168 T defaultValue = vector<T>(1)[0];
2169 vector<T>& data = extra.asVector<T>();
2170 for (size_t k = 0; k < slice.size(); ++k) {
2171 data[slice[k]] = defaultValue;
2172 }
2173}
2174
2175template<class T>
2176void resetMulti(AnyValue& extra, const vector<int>& slice)
2177{
2178 vector<T> defaultValue = vector<T>(extra.matrixShape().second);
2179 vector<vector<T>>& data = extra.asVector<vector<T>>();
2180 for (size_t k = 0; k < slice.size(); ++k) {
2181 data[slice[k]] = defaultValue;
2182 }
2183}
2184
2185template<class T>
2186void setAuxiliarySingle(size_t loc, AnyValue& extra, const AnyValue& value)
2187{
2188 extra.asVector<T>()[loc] = value.as<T>();
2189}
2190
2191template<class T>
2192void setAuxiliaryMulti(size_t loc, AnyValue& extra, const AnyValue& data)
2193{
2194 const auto& value = data.asVector<T>();
2195 auto& vec = extra.asVector<vector<T>>();
2196 if (value.size() != vec[loc].size()) {
2197 throw CanteraError("SolutionArray::setAuxiliaryMulti",
2198 "New element size {} does not match existing column size {}.",
2199 value.size(), vec[loc].size());
2200 }
2201 vec[loc] = value;
2202}
2203
2204} // end unnamed namespace
2205
2206}
SurfPhase.h
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase,...
ThermoPhase.h
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
Cantera::AnyBase::setLoc
void setLoc(int line, int column)
For values which are derived from an input file, set the line and column of this value in that file.
Definition AnyMap.cpp:617
Cantera::AnyMap
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:431
Cantera::AnyMap::empty
bool empty() const
Return boolean indicating whether AnyMap is empty.
Definition AnyMap.cpp:1468
Cantera::AnyMap::clear
void clear()
Erase all items in the mapping.
Definition AnyMap.cpp:1488
Cantera::AnyMap::at
const AnyValue & at(const string &key) const
Get the value of the item stored in key.
Definition AnyMap.cpp:1458
Cantera::AnyMap::update
void update(const AnyMap &other, bool keepExisting=true)
Add items from other to this AnyMap.
Definition AnyMap.cpp:1493
Cantera::AnyValue
A wrapper for a variable whose type is determined at runtime.
Definition AnyMap.h:88
Cantera::AnyValue::isVector
bool isVector() const
Returns true if the held value is a vector of the specified type, such as vector<double>.
Definition AnyMap.inl.h:75
Cantera::AnyValue::matrixShape
pair< size_t, size_t > matrixShape() const
Returns rows and columns of a matrix.
Definition AnyMap.cpp:714
Cantera::AnyValue::vectorSize
size_t vectorSize() const
Returns size of the held vector.
Definition AnyMap.cpp:698
Cantera::AnyValue::isScalar
bool isScalar() const
Returns true if the held value is a scalar type (such as double, long int, string,...
Definition AnyMap.cpp:694
Cantera::AnyValue::is
bool is() const
Returns true if the held value is of the specified type.
Definition AnyMap.inl.h:68
Cantera::AnyValue::asVector
const vector< T > & asVector(size_t nMin=npos, size_t nMax=npos) const
Return the held value, if it is a vector of type T.
Definition AnyMap.inl.h:109
Cantera::AnyValue::as
const T & as() const
Get the value of this key as the specified type.
Definition AnyMap.inl.h:16
Cantera::AnyValue::type_str
string type_str() const
Returns a string specifying the type of the held value.
Definition AnyMap.cpp:686
Cantera::CanteraError
Base class for exceptions thrown by Cantera classes.
Definition ctexceptions.h:66
Cantera::CanteraError::what
const char * what() const override
Get a description of the error.
Definition ctexceptions.cpp:29
Cantera::CanteraError::getMessage
virtual string getMessage() const
Method overridden by derived classes to format the error message.
Definition ctexceptions.cpp:52
Cantera::IndexError
An array index is out of range.
Definition ctexceptions.h:172
Cantera::NotImplementedError
An error indicating that an unimplemented function has been called.
Definition ctexceptions.h:201
Cantera::Storage
A wrapper class handling storage to HDF.
Definition Storage.h:39
Cantera::Storage::contents
pair< size_t, set< string > > contents(const string &id) const
Retrieve contents of file from a specified location.
Definition Storage.cpp:611
Cantera::Storage::writeAttributes
void writeAttributes(const string &id, const AnyMap &meta)
Write attributes to a specified location.
Definition Storage.cpp:629
Cantera::Storage::deleteGroup
void deleteGroup(const string &id)
Delete group.
Definition Storage.cpp:605
Cantera::Storage::readAttributes
AnyMap readAttributes(const string &id, bool recursive) const
Read attributes from a specified location.
Definition Storage.cpp:623
Cantera::Storage::readData
AnyValue readData(const string &id, const string &name, size_t rows, size_t cols=npos) const
Read dataset from a specified location.
Definition Storage.cpp:635
Cantera::Storage::checkGroup
bool checkGroup(const string &id, bool permissive=false)
Check whether path location exists.
Definition Storage.cpp:599
Cantera::Storage::setCompressionLevel
void setCompressionLevel(int level)
Set compression level (0..9)
Definition Storage.cpp:587
Cantera::Storage::writeData
void writeData(const string &id, const string &name, const AnyValue &data)
Write dataset to a specified location.
Definition Storage.cpp:642
Cantera::tokenizePath
void tokenizePath(const string &in_val, vector< string > &v)
This function separates a string up into tokens according to the location of path separators.
Definition stringUtils.cpp:196
Cantera::toLowerCopy
string toLowerCopy(const string &input)
Convert to lower case.
Definition stringUtils.cpp:219
Cantera::len
U len(const T &container)
Get the size of a container, cast to a signed integer type.
Definition utilities.h:198
Cantera::dot
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition utilities.h:82
Cantera::warn_user
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition global.h:268
Cantera
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
Cantera::npos
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
Cantera::_componentAliasMap
const map< string, string > & _componentAliasMap()
Return mapping of component alias names to standardized component names.
Definition SolutionArray.cpp:61
Cantera::offset
offset
Offsets of solution components in the 1D solution array.
Definition Flow1D.h:25
stringUtils.h
Contains declarations for string manipulation functions within Cantera.
utilities.h
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...