Cantera  4.0.0a1
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, false) != 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, false);
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, false);
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(
822 span<const double>(m_data->data() + m_loc * m_stride, nState));
823 }
824}
825
826void SolutionArray::updateState(int loc)
827{
828 setLoc(loc, false);
829 size_t nState = m_sol->thermo()->stateSize();
830 m_sol->thermo()->saveState(
831 span<double>(m_data->data() + m_loc * m_stride, nState));
832}
833
834vector<double> SolutionArray::getState(int loc)
835{
836 setLoc(loc);
837 size_t nState = m_sol->thermo()->stateSize();
838 vector<double> out(nState);
839 m_sol->thermo()->saveState(out); // thermo contains current state
840 return out;
841}
842
843void SolutionArray::setState(int loc, const vector<double>& state)
844{
845 size_t nState = m_sol->thermo()->stateSize();
846 if (state.size() != nState) {
847 throw CanteraError("SolutionArray::setState",
848 "Expected array to have length {}, but received an array of length {}.",
849 nState, state.size());
850 }
851 setLoc(loc, false);
852 m_sol->thermo()->restoreState(state);
853 m_sol->thermo()->saveState(
854 span<double>(m_data->data() + m_loc * m_stride, nState));
855}
856
857void SolutionArray::normalize() {
858 auto phase = m_sol->thermo();
859 auto nativeState = phase->nativeState();
860 if (nativeState.size() < 3) {
861 return;
862 }
863 size_t nState = phase->stateSize();
864 vector<double> out(nState);
865 if (nativeState.count("Y")) {
866 size_t offset = nativeState["Y"];
867 for (int loc = 0; loc < static_cast<int>(m_size); loc++) {
868 setLoc(loc, true); // set location and restore state
869 auto* data = m_data->data() + m_loc * m_stride + offset;
870 phase->setMassFractions(span<const double>(data, phase->nSpecies()));
871 m_sol->thermo()->saveState(out);
872 setState(loc, out);
873 }
874 } else if (nativeState.count("X")) {
875 size_t offset = nativeState["X"];
876 for (int loc = 0; loc < static_cast<int>(m_size); loc++) {
877 setLoc(loc, true); // set location and restore state
878 auto* data = m_data->data() + m_loc * m_stride + offset;
879 phase->setMoleFractions(span<const double>(data, phase->nSpecies()));
880 m_sol->thermo()->saveState(out);
881 setState(loc, out);
882 }
883 } else {
884 throw NotImplementedError("SolutionArray::normalize",
885 "Not implemented for mode '{}'.", phase->nativeMode());
886 }
887}
888
889AnyMap SolutionArray::getAuxiliary(int loc)
890{
891 setLoc(loc);
892 AnyMap out;
893 for (const auto& [key, extra] : *m_extra) {
894 if (extra.is<void>()) {
895 out[key] = extra;
896 } else if (extra.isVector<long int>()) {
897 out[key] = extra.asVector<long int>()[m_loc];
898 } else if (extra.isVector<double>()) {
899 out[key] = extra.asVector<double>()[m_loc];
900 } else if (extra.isVector<string>()) {
901 out[key] = extra.asVector<string>()[m_loc];
902 } else if (extra.isVector<vector<long int>>()) {
903 out[key] = extra.asVector<vector<long int>>()[m_loc];
904 } else if (extra.isVector<vector<double>>()) {
905 out[key] = extra.asVector<vector<double>>()[m_loc];
906 } else if (extra.isVector<vector<string>>()) {
907 out[key] = extra.asVector<vector<string>>()[m_loc];
908 } else {
909 throw NotImplementedError("SolutionArray::getAuxiliary",
910 "Unable to retrieve data for component '{}' with type '{}'.",
911 key, extra.type_str());
912 }
913 }
914 return out;
915}
916
917void SolutionArray::setAuxiliary(int loc, const AnyMap& data)
918{
919 setLoc(loc, false);
920 for (const auto& [name, value] : data) {
921 if (!m_extra->count(name)) {
922 throw CanteraError("SolutionArray::setAuxiliary",
923 "Unknown auxiliary component '{}'.", name);
924 }
925 auto& extra = m_extra->at(name);
926 if (extra.is<void>()) {
927 if (m_dataSize > 1) {
928 throw CanteraError("SolutionArray::setAuxiliary",
929 "Unable to set location for type '{}': "
930 "component is not initialized.", name);
931 }
932 _initExtra(name, value);
933 _resizeExtra(name);
934 }
935 try {
936 if (extra.isVector<long int>()) {
937 setAuxiliarySingle<long int>(m_loc, extra, value);
938 } else if (extra.isVector<double>()) {
939 setAuxiliarySingle<double>(m_loc, extra, value);
940 } else if (extra.isVector<string>()) {
941 setAuxiliarySingle<string>(m_loc, extra, value);
942 } else if (extra.isVector<vector<long int>>()) {
943 setAuxiliaryMulti<long int>(m_loc, extra, value);
944 } else if (extra.isVector<vector<double>>()) {
945 setAuxiliaryMulti<double>(m_loc, extra, value);
946 } else if (extra.isVector<vector<string>>()) {
947 setAuxiliaryMulti<string>(m_loc, extra, value);
948 } else {
949 throw CanteraError("SolutionArray::setAuxiliary",
950 "Unable to set entry for type '{}'.", extra.type_str());
951 }
952 } catch (CanteraError& err) {
953 // make failed type conversions traceable
954 throw CanteraError("SolutionArray::setAuxiliary",
955 "Encountered incompatible value for component '{}':\n{}",
956 name, err.getMessage());
957 }
958 }
959}
960
961AnyMap preamble(const string& desc)
962{
963 AnyMap data;
964 if (desc.size()) {
965 data["description"] = desc;
966 }
967 data["generator"] = "Cantera SolutionArray";
968 data["cantera-version"] = CANTERA_VERSION;
969 data["git-commit"] = gitCommit();
970
971 // Add a timestamp indicating the current time
972 time_t aclock;
973 ::time(&aclock); // Get time in seconds
974 struct tm* newtime = localtime(&aclock); // Convert time to struct tm form
975 data["date"] = stripnonprint(asctime(newtime));
976
977 return data;
978}
979
980AnyMap& openField(AnyMap& root, const string& name)
981{
982 if (!name.size()) {
983 return root;
984 }
985
986 // locate field based on 'name'
987 vector<string> tokens = tokenizePath(name);
988 AnyMap* ptr = &root; // use raw pointer to avoid copying
989 string path = "";
990 for (auto& field : tokens) {
991 path += "/" + field;
992 AnyMap& sub = *ptr;
993 if (sub.hasKey(field) && !sub[field].is<AnyMap>()) {
994 throw CanteraError("openField",
995 "Encountered invalid existing field '{}'.", path);
996 } else if (!sub.hasKey(field)) {
997 sub[field] = AnyMap();
998 }
999 ptr = &sub[field].as<AnyMap>();
1000 }
1001 return *ptr;
1002}
1003
1004void SolutionArray::writeHeader(const string& fname, const string& name,
1005 const string& desc, bool overwrite)
1006{
1007 Storage file(fname, true);
1008 if (file.checkGroup(name, true)) {
1009 if (!overwrite) {
1010 throw CanteraError("SolutionArray::writeHeader",
1011 "Group name '{}' exists; use 'overwrite' argument to overwrite.", name);
1012 }
1013 file.deleteGroup(name);
1014 file.checkGroup(name, true);
1015 }
1016 file.writeAttributes(name, preamble(desc));
1017}
1018
1019void SolutionArray::writeHeader(AnyMap& root, const string& name,
1020 const string& desc, bool overwrite)
1021{
1022 AnyMap& data = openField(root, name);
1023 if (!data.empty()) {
1024 if (!overwrite) {
1025 throw CanteraError("SolutionArray::writeHeader",
1026 "Field name '{}' exists; use 'overwrite' argument to overwrite.", name);
1027 }
1028 data.clear();
1029 }
1030 data.update(preamble(desc));
1031}
1032
1033void SolutionArray::writeEntry(const string& fname, bool overwrite, const string& basis)
1034{
1035 if (apiNdim() != 1) {
1036 throw CanteraError("SolutionArray::writeEntry",
1037 "Tabular output of CSV data only works for 1D SolutionArray objects.");
1038 }
1039 set<string> speciesNames;
1040 for (const auto& species : m_sol->thermo()->speciesNames()) {
1041 speciesNames.insert(species);
1042 }
1043 bool mole;
1044 if (basis == "") {
1045 const auto& nativeState = m_sol->thermo()->nativeState();
1046 mole = nativeState.find("X") != nativeState.end();
1047 } else if (basis == "X" || basis == "mole") {
1048 mole = true;
1049 } else if (basis == "Y" || basis == "mass") {
1050 mole = false;
1051 } else {
1052 throw CanteraError("SolutionArray::writeEntry",
1053 "Invalid species basis '{}'.", basis);
1054 }
1055
1056 auto names = componentNames();
1057 size_t last = names.size() - 1;
1058 vector<AnyValue> components;
1059 vector<bool> isSpecies;
1060 fmt::memory_buffer header;
1061 for (const auto& key : names) {
1062 string label = key;
1063 size_t col;
1064 if (speciesNames.find(key) == speciesNames.end()) {
1065 // Pre-read component vectors
1066 isSpecies.push_back(false);
1067 components.emplace_back(getComponent(key));
1068 col = components.size() - 1;
1069 if (!components[col].isVector<double>() &&
1070 !components[col].isVector<long int>() &&
1071 !components[col].isVector<string>())
1072 {
1073 throw CanteraError("SolutionArray::writeEntry",
1074 "Multi-dimensional column '{}' is not supported for CSV output.",
1075 key);
1076 }
1077 } else {
1078 // Delay reading species data as base can be either mole or mass
1079 isSpecies.push_back(true);
1080 components.emplace_back(AnyValue());
1081 col = components.size() - 1;
1082 if (mole) {
1083 label = "X_" + label;
1084 } else {
1085 label = "Y_" + label;
1086 }
1087 }
1088 if (label.find("\"") != string::npos || label.find("\n") != string::npos) {
1089 throw NotImplementedError("SolutionArray::writeEntry",
1090 "Detected column name containing double quotes or line feeds: '{}'.",
1091 label);
1092 }
1093 string sep = (col == last) ? "" : ",";
1094 if (label.find(",") != string::npos) {
1095 fmt_append(header, "\"{}\"{}", label, sep);
1096 } else {
1097 fmt_append(header, "{}{}", label, sep);
1098 }
1099 }
1100
1101 // (Most) potential exceptions have been thrown; start writing data to file
1102 if (std::ifstream(fname).good()) {
1103 if (!overwrite) {
1104 throw CanteraError("SolutionArray::writeEntry",
1105 "File '{}' already exists; use option 'overwrite' to replace CSV file.",
1106 fname);
1107 }
1108 std::remove(fname.c_str());
1109 }
1110 std::ofstream output(fname);
1111 output << to_string(header) << std::endl;
1112
1113 size_t maxLen = npos;
1114 vector<double> buf(speciesNames.size(), 0.);
1115 for (int row = 0; row < static_cast<int>(m_size); row++) {
1116 fmt::memory_buffer line;
1117 if (maxLen != npos) {
1118 line.reserve(maxLen);
1119 }
1120 setLoc(row);
1121 if (mole) {
1122 m_sol->thermo()->getMoleFractions(buf);
1123 } else {
1124 m_sol->thermo()->getMassFractions(buf);
1125 }
1126
1127 size_t idx = 0;
1128 for (size_t col = 0; col < components.size(); col++) {
1129 string sep = (col == last) ? "" : ",";
1130 if (isSpecies[col]) {
1131 fmt_append(line, "{:.9g}{}", buf[idx++], sep);
1132 } else {
1133 auto& data = components[col];
1134 if (data.isVector<double>()) {
1135 fmt_append(line, "{:.9g}{}", data.asVector<double>()[row], sep);
1136 } else if (data.isVector<long int>()) {
1137 fmt_append(line, "{}{}", data.asVector<long int>()[row], sep);
1138 } else if (data.isVector<bool>()) {
1139 fmt_append(line, "{}{}",
1140 static_cast<bool>(data.asVector<bool>()[row]), sep);
1141 } else {
1142 auto value = data.asVector<string>()[row];
1143 if (value.find("\"") != string::npos ||
1144 value.find("\n") != string::npos)
1145 {
1146 throw NotImplementedError("SolutionArray::writeEntry",
1147 "Detected value containing double quotes or line feeds: "
1148 "'{}'", value);
1149 }
1150 if (value.find(",") != string::npos) {
1151 fmt_append(line, "\"{}\"{}", value, sep);
1152 } else {
1153 fmt_append(line, "{}{}", value, sep);
1154 }
1155 }
1156 }
1157 }
1158 output << to_string(line) << std::endl;
1159 maxLen = std::max(maxLen, line.size());
1160 }
1161}
1162
1163void SolutionArray::writeEntry(const string& fname, const string& name,
1164 const string& sub, bool overwrite, int compression)
1165{
1166 if (name == "") {
1167 throw CanteraError("SolutionArray::writeEntry",
1168 "Group name specifying root location must not be empty.");
1169 }
1170 if (m_size < m_dataSize) {
1171 throw NotImplementedError("SolutionArray::writeEntry",
1172 "Unable to save sliced data.");
1173 }
1174 Storage file(fname, true);
1175 if (compression) {
1176 file.setCompressionLevel(compression);
1177 }
1178 string path = name;
1179 if (sub != "") {
1180 path += "/" + sub;
1181 } else {
1182 path += "/data";
1183 }
1184 if (file.checkGroup(path, true)) {
1185 if (!overwrite) {
1186 throw CanteraError("SolutionArray::writeEntry",
1187 "Group name '{}' exists; use 'overwrite' argument to overwrite.", name);
1188 }
1189 file.deleteGroup(path);
1190 file.checkGroup(path, true);
1191 }
1192 file.writeAttributes(path, m_meta);
1193 AnyMap more;
1194 if (apiNdim() == 1) {
1195 more["size"] = int(m_dataSize);
1196 } else {
1197 more["api-shape"] = m_apiShape;
1198 }
1199 if (!m_meta.hasKey("transport-model") && m_sol->transport()) {
1200 more["transport-model"] = m_sol->transportModel();
1201 }
1202 more["components"] = componentNames();
1203 file.writeAttributes(path, more);
1204 if (!m_dataSize) {
1205 return;
1206 }
1207
1208 const auto& nativeState = m_sol->thermo()->nativeState();
1209 size_t nSpecies = m_sol->thermo()->nSpecies();
1210 for (auto& [key, offset] : nativeState) {
1211 if (key == "X" || key == "Y") {
1212 vector<vector<double>> prop;
1213 for (size_t i = 0; i < m_size; i++) {
1214 size_t first = offset + i * m_stride;
1215 prop.emplace_back(m_data->begin() + first,
1216 m_data->begin() + first + nSpecies);
1217 }
1218 AnyValue data;
1219 data = prop;
1220 file.writeData(path, key, data);
1221 } else {
1222 auto data = getComponent(key);
1223 file.writeData(path, key, data);
1224 }
1225 }
1226
1227 for (const auto& [key, value] : *m_extra) {
1228 if (isSimpleVector(value)) {
1229 file.writeData(path, key, value);
1230 } else if (value.is<void>()) {
1231 // skip unintialized component
1232 } else {
1233 throw NotImplementedError("SolutionArray::writeEntry",
1234 "Unable to save component '{}' with data type {}.",
1235 key, value.type_str());
1236 }
1237 }
1238}
1239
1240void SolutionArray::writeEntry(AnyMap& root, const string& name, const string& sub,
1241 bool overwrite)
1242{
1243 if (name == "") {
1244 throw CanteraError("SolutionArray::writeEntry",
1245 "Field name specifying root location must not be empty.");
1246 }
1247 if (m_size < m_dataSize) {
1248 throw NotImplementedError("SolutionArray::writeEntry",
1249 "Unable to save sliced data.");
1250 }
1251 string path = name;
1252 if (sub != "") {
1253 path += "/" + sub;
1254 } else {
1255 path += "/data";
1256 }
1257 AnyMap& data = openField(root, path);
1258 bool preexisting = !data.empty();
1259 if (preexisting && !overwrite) {
1260 throw CanteraError("SolutionArray::writeEntry",
1261 "Field name '{}' exists; use 'overwrite' argument to overwrite.", name);
1262 }
1263 if (apiNdim() == 1) {
1264 data["size"] = int(m_dataSize);
1265 } else {
1266 data["api-shape"] = m_apiShape;
1267 }
1268 if (m_sol->transport() && m_sol->transportModel() != "none") {
1269 data["transport-model"] = m_sol->transportModel();
1270 }
1271 data.update(m_meta);
1272
1273 if (m_size > 1) {
1274 data["components"] = componentNames();
1275 }
1276
1277 for (auto& [_, key] : *m_order) {
1278 data[key] = m_extra->at(key);
1279 }
1280
1281 auto phase = m_sol->thermo();
1282 if (m_size == 1) {
1283 setLoc(0);
1284 data["temperature"] = phase->temperature();
1285 data["pressure"] = phase->pressure();
1286 auto surf = std::dynamic_pointer_cast<SurfPhase>(phase);
1287 auto nSpecies = phase->nSpecies();
1288 vector<double> values(nSpecies);
1289 if (surf) {
1290 surf->getCoverages(values);
1291 } else {
1292 phase->getMassFractions(values);
1293 }
1294 AnyMap items;
1295 for (size_t k = 0; k < nSpecies; k++) {
1296 if (values[k] != 0.0) {
1297 items[phase->speciesName(k)] = values[k];
1298 }
1299 }
1300 if (surf) {
1301 data["coverages"] = std::move(items);
1302 } else {
1303 data["mass-fractions"] = std::move(items);
1304 }
1305 } else if (m_size > 1) {
1306 for (auto& code : phase->nativeMode()) {
1307 string name(1, code);
1308 if (name == "X" || name == "Y") {
1309 for (auto& spc : phase->speciesNames()) {
1310 data[spc] = getComponent(spc);
1311 }
1312 data["basis"] = name == "X" ? "mole" : "mass";
1313 } else {
1314 data[name] = getComponent(name);
1315 }
1316 }
1317 }
1318
1319 static bool reg = AnyMap::addOrderingRules("SolutionArray",
1320 {{"head", "type"}, {"head", "size"}, {"head", "basis"}});
1321 if (reg) {
1322 data["__type__"] = "SolutionArray";
1323 }
1324
1325 // If this is not replacing an existing solution, put it at the end
1326 if (!preexisting) {
1327 data.setLoc(INT_MAX, 0);
1328 }
1329}
1330
1331void SolutionArray::append(const vector<double>& state, const AnyMap& extra)
1332{
1333 if (apiNdim() > 1) {
1334 throw NotImplementedError("SolutionArray::append",
1335 "Unable to append multi-dimensional arrays.");
1336 }
1337
1338 int pos = size();
1339 resize(pos + 1);
1340 try {
1341 setState(pos, state);
1342 setAuxiliary(pos, extra);
1343 } catch (CanteraError& err) {
1344 // undo resize and rethrow error
1345 resize(pos);
1346 throw CanteraError("SolutionArray::append", err.getMessage());
1347 }
1348}
1349
1350void SolutionArray::save(const string& fname, const string& name, const string& sub,
1351 const string& desc, bool overwrite, int compression,
1352 const string& basis)
1353{
1354 if (m_size < m_dataSize) {
1355 throw NotImplementedError("SolutionArray::save",
1356 "Unable to save sliced data.");
1357 }
1358 size_t dot = fname.find_last_of(".");
1359 string extension = (dot != npos) ? toLowerCopy(fname.substr(dot + 1)) : "";
1360 if (extension == "csv") {
1361 if (name != "") {
1362 warn_user("SolutionArray::save",
1363 "Parameter 'name' not used for CSV output.");
1364 }
1365 writeEntry(fname, overwrite, basis);
1366 return;
1367 }
1368 if (basis != "") {
1369 warn_user("SolutionArray::save",
1370 "Argument 'basis' is not used for HDF or YAML output.", basis);
1371 }
1372 if (extension == "h5" || extension == "hdf" || extension == "hdf5") {
1373 writeHeader(fname, name, desc, overwrite);
1374 writeEntry(fname, name, sub, true, compression);
1375 return;
1376 }
1377 if (extension == "yaml" || extension == "yml") {
1378 // Check for an existing file and load it if present
1379 AnyMap data;
1380 std::ifstream file(fname);
1381 if (file.good() && file.peek() != std::ifstream::traits_type::eof()) {
1382 data = AnyMap::fromYamlFile(fname);
1383 }
1384 writeHeader(data, name, desc, overwrite);
1385 writeEntry(data, name, sub, true);
1386
1387 // Write the output file and remove the now-outdated cached file
1388 std::ofstream out(fname);
1389 out << data.toYamlString();
1390 AnyMap::clearCachedFile(fname);
1391 return;
1392 }
1393 throw CanteraError("SolutionArray::save",
1394 "Unknown file extension '{}'.", extension);
1395}
1396
1397AnyMap SolutionArray::readHeader(const string& fname, const string& name)
1398{
1399 Storage file(fname, false);
1400 file.checkGroup(name);
1401 return file.readAttributes(name, false);
1402}
1403
1404const AnyMap& locateField(const AnyMap& root, const string& name)
1405{
1406 if (!name.size()) {
1407 return root;
1408 }
1409
1410 // locate field based on 'name'
1411 vector<string> tokens = tokenizePath(name);
1412 const AnyMap* ptr = &root; // use raw pointer to avoid copying
1413 string path = "";
1414 for (auto& field : tokens) {
1415 path += "/" + field;
1416 const AnyMap& sub = *ptr;
1417 if (!sub.hasKey(field) || !sub[field].is<AnyMap>()) {
1418 throw CanteraError("SolutionArray::locateField",
1419 "No field or solution with name '{}'.", path);
1420 }
1421 ptr = &sub[field].as<AnyMap>();
1422 }
1423 return *ptr;
1424}
1425
1426AnyMap SolutionArray::readHeader(const AnyMap& root, const string& name)
1427{
1428 auto sub = locateField(root, name);
1429 AnyMap header;
1430 for (const auto& [key, value] : sub) {
1431 if (!sub[key].is<AnyMap>()) {
1432 header[key] = value;
1433 }
1434 }
1435 return header;
1436}
1437
1438AnyMap SolutionArray::restore(const string& fname,
1439 const string& name, const string& sub)
1440{
1441 size_t dot = fname.find_last_of(".");
1442 string extension = (dot != npos) ? toLowerCopy(fname.substr(dot + 1)) : "";
1443 AnyMap header;
1444 if (extension == "csv") {
1445 throw NotImplementedError("SolutionArray::restore",
1446 "CSV import not implemented; if using Python, data can be imported via "
1447 "'read_csv' instead.");
1448 }
1449 if (extension == "h5" || extension == "hdf" || extension == "hdf5") {
1450 readEntry(fname, name, sub);
1451 header = readHeader(fname, name);
1452 } else if (extension == "yaml" || extension == "yml") {
1453 const AnyMap& root = AnyMap::fromYamlFile(fname);
1454 readEntry(root, name, sub);
1455 header = readHeader(root, name);
1456 } else {
1457 throw CanteraError("SolutionArray::restore",
1458 "Unknown file extension '{}'; supported extensions include "
1459 "'h5'/'hdf'/'hdf5' and 'yml'/'yaml'.", extension);
1460 }
1461 return header;
1462}
1463
1464void SolutionArray::_initExtra(const string& name, const AnyValue& value)
1465{
1466 if (!m_extra->count(name)) {
1467 throw CanteraError("SolutionArray::_initExtra",
1468 "Component '{}' does not exist.", name);
1469 }
1470 auto& extra = (*m_extra)[name];
1471 if (!extra.is<void>()) {
1472 throw CanteraError("SolutionArray::_initExtra",
1473 "Component '{}' is already initialized.", name);
1474 }
1475 try {
1476 if (value.is<long int>()) {
1477 extra = vector<long int>(m_dataSize, value.as<long int>());
1478 } else if (value.is<double>()) {
1479 extra = vector<double>(m_dataSize, value.as<double>());
1480 } else if (value.is<string>()) {
1481 extra = vector<string>(m_dataSize, value.as<string>());
1482 } else if (value.isVector<long int>()) {
1483 extra = vector<vector<long int>>(m_dataSize, value.asVector<long int>());
1484 } else if (value.isVector<double>()) {
1485 extra = vector<vector<double>>(m_dataSize, value.asVector<double>());
1486 } else if (value.isVector<string>()) {
1487 extra = vector<vector<string>>(m_dataSize, value.asVector<string>());
1488 } else if (value.is<void>()) {
1489 // placeholder for unknown type; settable in setComponent
1490 extra = AnyValue();
1491 } else {
1492 throw NotImplementedError("SolutionArray::_initExtra",
1493 "Unable to initialize component '{}' with type '{}'.",
1494 name, value.type_str());
1495 }
1496 } catch (CanteraError& err) {
1497 // make failed type conversions traceable
1498 throw CanteraError("SolutionArray::_initExtra",
1499 "Encountered incompatible value for initializing component '{}':\n{}",
1500 name, err.getMessage());
1501 }
1502}
1503
1504void SolutionArray::_resizeExtra(const string& name, const AnyValue& value)
1505{
1506 if (!m_extra->count(name)) {
1507 throw CanteraError("SolutionArray::_resizeExtra",
1508 "Component '{}' does not exist.", name);
1509 }
1510 auto& extra = (*m_extra)[name];
1511 if (extra.is<void>()) {
1512 // cannot resize placeholder (uninitialized component)
1513 return;
1514 }
1515
1516 try {
1517 if (extra.isVector<long int>()) {
1518 resizeSingle<long int>(extra, m_dataSize, value);
1519 } else if (extra.isVector<double>()) {
1520 resizeSingle<double>(extra, m_dataSize, value);
1521 } else if (extra.isVector<string>()) {
1522 resizeSingle<string>(extra, m_dataSize, value);
1523 } else if (extra.isVector<vector<double>>()) {
1524 resizeMulti<double>(extra, m_dataSize, value);
1525 } else if (extra.isVector<vector<long int>>()) {
1526 resizeMulti<long int>(extra, m_dataSize, value);
1527 } else if (extra.isVector<vector<string>>()) {
1528 resizeMulti<string>(extra, m_dataSize, value);
1529 } else {
1530 throw NotImplementedError("SolutionArray::_resizeExtra",
1531 "Unable to resize using type '{}'.", extra.type_str());
1532 }
1533 } catch (CanteraError& err) {
1534 // make failed type conversions traceable
1535 throw CanteraError("SolutionArray::_resizeExtra",
1536 "Encountered incompatible value for resizing component '{}':\n{}",
1537 name, err.getMessage());
1538 }
1539}
1540
1541void SolutionArray::_setExtra(const string& name, const AnyValue& data)
1542{
1543 if (!m_extra->count(name)) {
1544 throw CanteraError("SolutionArray::_setExtra",
1545 "Extra component '{}' does not exist.", name);
1546 }
1547
1548 auto& extra = m_extra->at(name);
1549 if (data.is<void>() && m_size == m_dataSize) {
1550 // reset placeholder
1551 extra = AnyValue();
1552 return;
1553 }
1554
1555 // uninitialized component
1556 if (extra.is<void>()) {
1557 if (m_size != m_dataSize) {
1558 throw CanteraError("SolutionArray::_setExtra",
1559 "Unable to replace '{}' for sliced data.", name);
1560 }
1561 if (data.vectorSize() == m_dataSize || data.matrixShape().first == m_dataSize) {
1562 // assign entire component
1563 extra = data;
1564 return;
1565 }
1566 // data contains default value
1567 if (m_dataSize == 0 && (data.isScalar() || data.vectorSize() > 0)) {
1568 throw CanteraError("SolutionArray::_setExtra",
1569 "Unable to initialize '{}' with non-empty values when SolutionArray is "
1570 "empty.", name);
1571 }
1572 if (m_dataSize && data.vectorSize() == 0) {
1573 throw CanteraError("SolutionArray::_setExtra",
1574 "Unable to initialize '{}' with empty array when SolutionArray is not "
1575 "empty." , name);
1576 }
1577 _initExtra(name, data);
1578 _resizeExtra(name, data);
1579 return;
1580 }
1581
1582 if (data.is<long int>()) {
1583 setScalar<long int>(extra, data, m_active);
1584 } else if (data.is<double>()) {
1585 setScalar<double>(extra, data, m_active);
1586 } else if (data.is<string>()) {
1587 setScalar<string>(extra, data, m_active);
1588 } else if (data.isVector<long int>()) {
1589 setSingle<long int>(extra, data, m_active);
1590 } else if (data.isVector<double>()) {
1591 setSingle<double>(extra, data, m_active);
1592 } else if (data.isVector<string>()) {
1593 setSingle<string>(extra, data, m_active);
1594 } else if (data.isVector<vector<long int>>()) {
1595 setMulti<long int>(extra, data, m_active);
1596 } else if (data.isVector<vector<double>>()) {
1597 setMulti<double>(extra, data, m_active);
1598 } else if (data.isVector<vector<string>>()) {
1599 setMulti<string>(extra, data, m_active);
1600 } else {
1601 throw NotImplementedError("SolutionArray::_setExtra",
1602 "Unable to set sliced data for component '{}' with type '{}'.",
1603 name, data.type_str());
1604 }
1605}
1606
1607string SolutionArray::_detectMode(const set<string>& names, bool native)
1608{
1609 // check set of available names against state acronyms defined by Phase::fullStates
1610 string mode = "";
1611 const auto& nativeState = m_sol->thermo()->nativeState();
1612 bool usesNativeState = false;
1613 auto surf = std::dynamic_pointer_cast<SurfPhase>(m_sol->thermo());
1614 bool found;
1615 for (const auto& item : m_sol->thermo()->fullStates()) {
1616 found = true;
1617 string name;
1618 usesNativeState = true;
1619 for (size_t i = 0; i < item.size(); i++) {
1620 // pick i-th letter from "full" state acronym
1621 name = string(1, item[i]);
1622 if (surf && (name == "X" || name == "Y")) {
1623 // override native state to enable detection of surface phases
1624 name = "C";
1625 usesNativeState = false;
1626 break;
1627 }
1628 if (names.count(name)) {
1629 // property is stored using letter acronym
1630 usesNativeState &= nativeState.count(name) > 0;
1631 } else if (aliasMap.count(name) && names.count(aliasMap.at(name))) {
1632 // property is stored using property name
1633 usesNativeState &= nativeState.count(name) > 0;
1634 } else {
1635 found = false;
1636 break;
1637 }
1638 }
1639 if (found) {
1640 mode = (name == "C") ? item.substr(0, 2) + "C" : item;
1641 break;
1642 }
1643 }
1644 if (!found) {
1645 if (surf && names.count("T") && names.count("X") && names.count("density")) {
1646 // Legacy HDF format erroneously uses density (Cantera < 3.0)
1647 return "legacySurf";
1648 }
1649 if (names.count("mass-flux") && names.count("mass-fractions")) {
1650 // Legacy YAML format stores incomplete state information (Cantera < 3.0)
1651 return "legacyInlet";
1652 }
1653 throw CanteraError("SolutionArray::_detectMode",
1654 "Detected incomplete thermodynamic information. Full states for a '{}' "
1655 "phase\nmay be defined by the following modes:\n'{}'\n"
1656 "Available data are: '{}'", m_sol->thermo()->type(),
1657 ba::join(m_sol->thermo()->fullStates(), "', '"), ba::join(names, "', '"));
1658 }
1659 if (usesNativeState && native) {
1660 return "native";
1661 }
1662 return mode;
1663}
1664
1665set<string> SolutionArray::_stateProperties(
1666 const string& mode, bool alias)
1667{
1668 set<string> states;
1669 if (mode == "native") {
1670 for (const auto& [name, offset] : m_sol->thermo()->nativeState()) {
1671 states.insert(alias ? aliasMap.at(name) : name);
1672 }
1673 } else {
1674 for (const auto& m : mode) {
1675 const string name = string(1, m);
1676 states.insert(alias ? aliasMap.at(name) : name);
1677 }
1678 }
1679
1680 return states;
1681}
1682
1683string getName(const set<string>& names, const string& name)
1684{
1685 if (names.count(name)) {
1686 return name;
1687 }
1688 const auto& alias = aliasMap.at(name);
1689 if (names.count(alias)) {
1690 return alias;
1691 }
1692 return name; // let exception be thrown elsewhere
1693}
1694
1695void SolutionArray::readEntry(const string& fname, const string& name,
1696 const string& sub)
1697{
1698 Storage file(fname, false);
1699 if (name == "") {
1700 throw CanteraError("SolutionArray::readEntry",
1701 "Group name specifying root location must not be empty.");
1702 }
1703 string path = name;
1704 if (sub != "" && file.checkGroup(name + "/" + sub, true)) {
1705 path += "/" + sub;
1706 } else if (sub == "" && file.checkGroup(name + "/data", true)) {
1707 // default data location
1708 path += "/data";
1709 }
1710 if (!file.checkGroup(path)) {
1711 throw CanteraError("SolutionArray::readEntry",
1712 "Group name specifying data entry is empty.");
1713 }
1714 m_extra->clear();
1715 auto [size, names] = file.contents(path);
1716 m_meta = file.readAttributes(path, true);
1717 if (m_meta.hasKey("size")) {
1718 // one-dimensional array
1719 resize(m_meta["size"].as<long int>());
1720 m_meta.erase("size");
1721 } else if (m_meta.hasKey("api-shape")) {
1722 // API uses multiple dimensions to interpret C++ SolutionArray
1723 setApiShape(m_meta["api-shape"].asVector<long int>());
1724 m_meta.erase("api-shape");
1725 } else {
1726 // legacy format; size needs to be detected
1727 resize(static_cast<int>(size));
1728 }
1729
1730 if (m_size == 0) {
1731 return;
1732 }
1733
1734 // determine storage mode of state data
1735 string mode = _detectMode(names);
1736 set<string> states = _stateProperties(mode);
1737 if (states.count("C")) {
1738 if (names.count("X")) {
1739 states.erase("C");
1740 states.insert("X");
1741 mode = mode.substr(0, 2) + "X";
1742 } else if (names.count("Y")) {
1743 states.erase("C");
1744 states.insert("Y");
1745 mode = mode.substr(0, 2) + "Y";
1746 }
1747 }
1748
1749 // restore state data
1750 size_t nSpecies = m_sol->thermo()->nSpecies();
1751 size_t nState = m_sol->thermo()->stateSize();
1752 const auto& nativeStates = m_sol->thermo()->nativeState();
1753 if (mode == "native") {
1754 // native state can be written directly into data storage
1755 for (const auto& [name, offset] : nativeStates) {
1756 if (name == "X" || name == "Y") {
1757 AnyValue data;
1758 data = file.readData(path, name, m_size, nSpecies);
1759 auto prop = data.asVector<vector<double>>();
1760 for (size_t i = 0; i < m_dataSize; i++) {
1761 std::copy(prop[i].begin(), prop[i].end(),
1762 m_data->data() + offset + i * m_stride);
1763 }
1764 } else {
1765 AnyValue data;
1766 data = file.readData(path, getName(names, name), m_dataSize, 0);
1767 setComponent(name, data);
1768 }
1769 }
1770 } else if (mode == "TPX") {
1771 AnyValue data;
1772 data = file.readData(path, getName(names, "T"), m_dataSize, 0);
1773 vector<double> T = std::move(data.asVector<double>());
1774 data = file.readData(path, getName(names, "P"), m_dataSize, 0);
1775 vector<double> P = std::move(data.asVector<double>());
1776 data = file.readData(path, "X", m_dataSize, nSpecies);
1777 vector<vector<double>> X = std::move(data.asVector<vector<double>>());
1778 for (size_t i = 0; i < m_dataSize; i++) {
1779 m_sol->thermo()->setMoleFractions_NoNorm(X[i]);
1780 m_sol->thermo()->setState_TP(T[i], P[i]);
1781 m_sol->thermo()->saveState(
1782 span<double>(m_data->data() + i * m_stride, nState));
1783 }
1784 } else if (mode == "TDX") {
1785 AnyValue data;
1786 data = file.readData(path, getName(names, "T"), m_dataSize, 0);
1787 vector<double> T = std::move(data.asVector<double>());
1788 data = file.readData(path, getName(names, "D"), m_dataSize, 0);
1789 vector<double> D = std::move(data.asVector<double>());
1790 data = file.readData(path, "X", m_dataSize, nSpecies);
1791 vector<vector<double>> X = std::move(data.asVector<vector<double>>());
1792 for (size_t i = 0; i < m_dataSize; i++) {
1793 m_sol->thermo()->setMoleFractions_NoNorm(X[i]);
1794 m_sol->thermo()->setState_TD(T[i], D[i]);
1795 m_sol->thermo()->saveState(
1796 span<double>(m_data->data() + i * m_stride, nState));
1797 }
1798 } else if (mode == "TPY") {
1799 AnyValue data;
1800 data = file.readData(path, getName(names, "T"), m_dataSize, 0);
1801 vector<double> T = std::move(data.asVector<double>());
1802 data = file.readData(path, getName(names, "P"), m_dataSize, 0);
1803 vector<double> P = std::move(data.asVector<double>());
1804 data = file.readData(path, "Y", m_dataSize, nSpecies);
1805 vector<vector<double>> Y = std::move(data.asVector<vector<double>>());
1806 for (size_t i = 0; i < m_dataSize; i++) {
1807 m_sol->thermo()->setMassFractions_NoNorm(Y[i]);
1808 m_sol->thermo()->setState_TP(T[i], P[i]);
1809 m_sol->thermo()->saveState(
1810 span<double>(m_data->data() + i * m_stride, nState));
1811 }
1812 } else if (mode == "legacySurf") {
1813 // erroneous TDX mode (should be TPX or TPY) - Sim1D (Cantera 2.5)
1814 AnyValue data;
1815 data = file.readData(path, getName(names, "T"), m_dataSize, 0);
1816 vector<double> T = std::move(data.asVector<double>());
1817 data = file.readData(path, "X", m_dataSize, nSpecies);
1818 vector<vector<double>> X = std::move(data.asVector<vector<double>>());
1819 for (size_t i = 0; i < m_dataSize; i++) {
1820 m_sol->thermo()->setMoleFractions_NoNorm(X[i]);
1821 m_sol->thermo()->setTemperature(T[i]);
1822 m_sol->thermo()->saveState(
1823 span<double>(m_data->data() + i * m_stride, nState));
1824 }
1825 warn_user("SolutionArray::readEntry",
1826 "Detected legacy HDF format with incomplete state information\nfor name "
1827 "'{}' (pressure missing).", path);
1828 } else if (mode == "") {
1829 throw CanteraError("SolutionArray::readEntry",
1830 "Data are not consistent with full state modes.");
1831 } else {
1832 throw NotImplementedError("SolutionArray::readEntry",
1833 "Import of '{}' data is not supported.", mode);
1834 }
1835
1836 // restore remaining data
1837 if (m_meta.hasKey("components")) {
1838 const auto& components = m_meta["components"].asVector<string>();
1839 bool back = false;
1840 for (const auto& name : components) {
1841 if (hasComponent(name, false) || name == "X" || name == "Y") {
1842 back = true;
1843 } else {
1844 auto _name = name;
1845 if (reverseAliasMap.count(name)) {
1846 _name = reverseAliasMap.at(name);
1847 }
1848 addExtra(_name, back);
1849 AnyValue data;
1850 data = file.readData(path, name, m_dataSize);
1851 setComponent(_name, data);
1852 }
1853 }
1854 m_meta.erase("components");
1855 } else {
1856 // data format used by Python h5py export (Cantera 2.5)
1857 warn_user("SolutionArray::readEntry", "Detected legacy HDF format.");
1858 for (const auto& name : names) {
1859 if (!hasComponent(name, false) && name != "X" && name != "Y") {
1860 auto _name = name;
1861 if (reverseAliasMap.count(name)) {
1862 _name = reverseAliasMap.at(name);
1863 }
1864 addExtra(_name);
1865 AnyValue data;
1866 data = file.readData(path, name, m_dataSize);
1867 setComponent(_name, data);
1868 }
1869 }
1870 }
1871
1872 if (m_meta.hasKey("transport-model")) {
1873 m_sol->setTransportModel(m_meta["transport-model"].asString());
1874 }
1875}
1876
1877void SolutionArray::readEntry(const AnyMap& root, const string& name, const string& sub)
1878{
1879 if (name == "") {
1880 throw CanteraError("SolutionArray::readEntry",
1881 "Field name specifying root location must not be empty.");
1882 }
1883 auto path = locateField(root, name);
1884 if (path.hasKey("generator") && sub != "") {
1885 // read entry from subfolder (since Cantera 3.0)
1886 path = locateField(root, name + "/" + sub);
1887 } else if (sub == "" && path.hasKey("data")) {
1888 // default data location
1889 path = locateField(root, name + "/data");
1890 }
1891
1892 // set size and initialize
1893 long size = 0;
1894 if (path.hasKey("size")) {
1895 // one-dimensional array
1896 resize(path["size"].asInt());
1897 } else if (path.hasKey("api-shape")) {
1898 // API uses multiple dimensions to interpret C++ SolutionArray
1899 auto& shape = path["api-shape"].asVector<long int>();
1900 setApiShape(shape);
1901 } else {
1902 // legacy format (Cantera 2.6)
1903 size = path.getInt("points", 0);
1904 if (!path.hasKey("T") && !path.hasKey("temperature")) {
1905 // overwrite size - Sim1D erroneously assigns '1'
1906 size = (long)0;
1907 }
1908 resize(static_cast<int>(size));
1909 }
1910 m_extra->clear();
1911
1912 // restore data
1913 set<string> exclude = {"size", "api-shape", "points", "X", "Y"};
1914 set<string> names = path.keys();
1915 size_t nState = m_sol->thermo()->stateSize();
1916 if (m_dataSize == 0) {
1917 // no data points
1918 } else if (m_dataSize == 1) {
1919 // single data point
1920 string mode = _detectMode(names, false);
1921 if (mode == "TPY") {
1922 double T = path[getName(names, "T")].asDouble();
1923 double P = path[getName(names, "P")].asDouble();
1924 auto Y = path["mass-fractions"].asMap<double>();
1925 m_sol->thermo()->setState_TPY(T, P, Y);
1926 } else if (mode == "TPC") {
1927 auto surf = std::dynamic_pointer_cast<SurfPhase>(m_sol->thermo());
1928 double T = path[getName(names, "T")].asDouble();
1929 double P = path["pressure"].asDouble();
1930 m_sol->thermo()->setState_TP(T, P);
1931 auto cov = path["coverages"].asMap<double>();
1932 surf->setCoveragesByName(cov);
1933 } else if (mode == "legacyInlet") {
1934 // missing property - Sim1D (Cantera 2.6)
1935 mode = "TPY";
1936 double T = path[getName(names, "T")].asDouble();
1937 auto Y = path["mass-fractions"].asMap<double>();
1938 m_sol->thermo()->setState_TPY(T, m_sol->thermo()->pressure(), Y);
1939 warn_user("SolutionArray::readEntry",
1940 "Detected legacy YAML format with incomplete state information\n"
1941 "for name '{}' (pressure missing).", name + "/" + sub);
1942 } else if (mode == "") {
1943 throw CanteraError("SolutionArray::readEntry",
1944 "Data are not consistent with full state modes.");
1945 } else {
1946 throw NotImplementedError("SolutionArray::readEntry",
1947 "Import of '{}' data is not supported.", mode);
1948 }
1949 m_sol->thermo()->saveState(span<double>(m_data->data(), nState));
1950 auto props = _stateProperties(mode, true);
1951 exclude.insert(props.begin(), props.end());
1952 } else {
1953 // multiple data points
1954 if (path.hasKey("components")) {
1955 const auto& components = path["components"].asVector<string>();
1956 bool back = false;
1957 for (const auto& name : components) {
1958 auto _name = name;
1959 if (hasComponent(name, false)) {
1960 back = true;
1961 } else {
1962 if (reverseAliasMap.count(name)) {
1963 _name = reverseAliasMap.at(name);
1964 }
1965 addExtra(_name, back);
1966 }
1967 setComponent(_name, path[name]);
1968 exclude.insert(name);
1969 }
1970 } else {
1971 // legacy YAML format does not provide for list of components
1972 for (const auto& [name, value] : path) {
1973 if (value.isVector<double>()) {
1974 const vector<double>& data = value.asVector<double>();
1975 if (data.size() == m_dataSize) {
1976 auto _name = name;
1977 if (!hasComponent(name, false)) {
1978 if (reverseAliasMap.count(name)) {
1979 _name = reverseAliasMap.at(name);
1980 }
1981 addExtra(_name);
1982 }
1983 setComponent(_name, value);
1984 exclude.insert(name);
1985 }
1986 }
1987 }
1988 }
1989
1990 // check that state data are complete
1991 const auto& nativeState = m_sol->thermo()->nativeState();
1992 set<string> props;
1993 set<string> missingProps;
1994 for (const auto& [name, offset] : nativeState) {
1995 if (exclude.count(name)) {
1996 props.insert(name);
1997 } else {
1998 missingProps.insert(name);
1999 }
2000 }
2001
2002 set<string> TY = {"T", "Y"};
2003 if (props == TY && missingProps.count("D") && path.hasKey("pressure")) {
2004 // missing "D" - Sim1D (Cantera 2.6)
2005 double P = path["pressure"].asDouble();
2006 const size_t offset_T = nativeState.find("T")->second;
2007 const size_t offset_D = nativeState.find("D")->second;
2008 const size_t offset_Y = nativeState.find("Y")->second;
2009 for (size_t i = 0; i < m_dataSize; i++) {
2010 double T = (*m_data)[offset_T + i * m_stride];
2011 span<double> Y(m_data->data() + offset_Y + i * m_stride,
2012 m_sol->thermo()->nSpecies());
2013 m_sol->thermo()->setState_TPY(T, P, Y);
2014 (*m_data)[offset_D + i * m_stride] = m_sol->thermo()->density();
2015 }
2016 } else if (missingProps.size()) {
2017 throw CanteraError("SolutionArray::readEntry",
2018 "Incomplete state information: missing '{}'.",
2019 ba::join(missingProps, "', '"));
2020 }
2021 }
2022
2023 // add meta data
2024 for (const auto& [name, value] : path) {
2025 if (!exclude.count(name)) {
2026 m_meta[name] = value;
2027 }
2028 }
2029
2030 if (m_meta.hasKey("transport-model")) {
2031 m_sol->setTransportModel(m_meta["transport-model"].asString());
2032 }
2033}
2034
2035namespace { // restrict scope of helper functions to local translation unit
2036
2037template<class T>
2038AnyValue getSingle(const AnyValue& extra, const vector<int>& slice)
2039{
2040 vector<T> data(slice.size());
2041 const auto& vec = extra.asVector<T>();
2042 for (size_t k = 0; k < slice.size(); ++k) {
2043 data[k] = vec[slice[k]];
2044 }
2045 AnyValue out;
2046 out = data;
2047 return out;
2048}
2049
2050template<class T>
2051AnyValue getMulti(const AnyValue& extra, const vector<int>& slice)
2052{
2053 vector<vector<T>> data(slice.size());
2054 const auto& vec = extra.asVector<vector<T>>();
2055 for (size_t k = 0; k < slice.size(); ++k) {
2056 data[k] = vec[slice[k]];
2057 }
2058 AnyValue out;
2059 out = data;
2060 return out;
2061}
2062
2063template<class T>
2064void setScalar(AnyValue& extra, const AnyValue& data, const vector<int>& slice)
2065{
2066 T value = data.as<T>();
2067 if (extra.isVector<T>()) {
2068 auto& vec = extra.asVector<T>();
2069 for (size_t k = 0; k < slice.size(); ++k) {
2070 vec[slice[k]] = value;
2071 }
2072 } else {
2073 throw CanteraError("SolutionArray::setScalar",
2074 "Incompatible input data: unable to assign '{}' data to '{}'.",
2075 data.type_str(), extra.type_str());
2076 }
2077}
2078
2079template<class T>
2080void setSingle(AnyValue& extra, const AnyValue& data, const vector<int>& slice)
2081{
2082 size_t size = slice.size();
2083 if (extra.vectorSize() == size && data.vectorSize() == size) {
2084 extra = data; // no slicing necessary; type can change
2085 return;
2086 }
2087 if (extra.matrixShape().first == size && data.vectorSize() == size) {
2088 extra = data; // no slicing necessary; shape and type can change
2089 return;
2090 }
2091 if (extra.type_str() != data.type_str()) {
2092 // do not allow changing of data type when slicing
2093 throw CanteraError("SolutionArray::setSingle",
2094 "Incompatible types: expected '{}' but received '{}'.",
2095 extra.type_str(), data.type_str());
2096 }
2097 const auto& vData = data.asVector<T>();
2098 if (vData.size() != size) {
2099 throw CanteraError("SolutionArray::setSingle",
2100 "Invalid input data size: expected {} entries but received {}.",
2101 size, vData.size());
2102 }
2103 auto& vec = extra.asVector<T>();
2104 for (size_t k = 0; k < size; ++k) {
2105 vec[slice[k]] = vData[k];
2106 }
2107}
2108
2109template<class T>
2110void setMulti(AnyValue& extra, const AnyValue& data, const vector<int>& slice)
2111{
2112 if (!data.isMatrix<T>()) {
2113 throw CanteraError("SolutionArray::setMulti",
2114 "Invalid input data shape: inconsistent number of columns.");
2115 }
2116 size_t size = slice.size();
2117 auto [rows, cols] = data.matrixShape();
2118 if (extra.matrixShape().first == size && rows == size) {
2119 extra = data; // no slicing necessary; type can change
2120 return;
2121 }
2122 if (extra.vectorSize() == size && rows == size) {
2123 extra = data; // no slicing necessary; shape and type can change
2124 return;
2125 }
2126 if (extra.type_str() != data.type_str()) {
2127 // do not allow changing of data type when slicing
2128 throw CanteraError("SolutionArray::setMulti",
2129 "Incompatible types: expected '{}' but received '{}'.",
2130 extra.type_str(), data.type_str());
2131 }
2132 if (rows != size) {
2133 throw CanteraError("SolutionArray::setMulti",
2134 "Invalid input data shape: expected {} rows but received {}.",
2135 size, rows);
2136 }
2137 if (extra.matrixShape().second != cols) {
2138 throw CanteraError("SolutionArray::setMulti",
2139 "Invalid input data shape: expected {} columns but received {}.",
2140 extra.matrixShape().second, cols);
2141 }
2142 const auto& vData = data.asVector<vector<T>>();
2143 auto& vec = extra.asVector<vector<T>>();
2144 for (size_t k = 0; k < slice.size(); ++k) {
2145 vec[slice[k]] = vData[k];
2146 }
2147}
2148
2149template<class T>
2150void resizeSingle(AnyValue& extra, size_t size, const AnyValue& value)
2151{
2152 T defaultValue;
2153 if (value.is<void>()) {
2154 defaultValue = vector<T>(1)[0];
2155 } else {
2156 defaultValue = value.as<T>();
2157 }
2158 extra.asVector<T>().resize(size, defaultValue);
2159}
2160
2161template<class T>
2162void resizeMulti(AnyValue& extra, size_t size, const AnyValue& value)
2163{
2164 vector<T> defaultValue;
2165 if (value.is<void>()) {
2166 defaultValue = vector<T>(extra.matrixShape().second);
2167 } else {
2168 defaultValue = value.as<vector<T>>();
2169 }
2170 extra.asVector<vector<T>>().resize(size, defaultValue);
2171}
2172
2173template<class T>
2174void resetSingle(AnyValue& extra, const vector<int>& slice)
2175{
2176 T defaultValue = vector<T>(1)[0];
2177 vector<T>& data = extra.asVector<T>();
2178 for (size_t k = 0; k < slice.size(); ++k) {
2179 data[slice[k]] = defaultValue;
2180 }
2181}
2182
2183template<class T>
2184void resetMulti(AnyValue& extra, const vector<int>& slice)
2185{
2186 vector<T> defaultValue = vector<T>(extra.matrixShape().second);
2187 vector<vector<T>>& data = extra.asVector<vector<T>>();
2188 for (size_t k = 0; k < slice.size(); ++k) {
2189 data[slice[k]] = defaultValue;
2190 }
2191}
2192
2193template<class T>
2194void setAuxiliarySingle(size_t loc, AnyValue& extra, const AnyValue& value)
2195{
2196 extra.asVector<T>()[loc] = value.as<T>();
2197}
2198
2199template<class T>
2200void setAuxiliaryMulti(size_t loc, AnyValue& extra, const AnyValue& data)
2201{
2202 const auto& value = data.asVector<T>();
2203 auto& vec = extra.asVector<vector<T>>();
2204 if (value.size() != vec[loc].size()) {
2205 throw CanteraError("SolutionArray::setAuxiliaryMulti",
2206 "New element size {} does not match existing column size {}.",
2207 value.size(), vec[loc].size());
2208 }
2209 vec[loc] = value;
2210}
2211
2212} // end unnamed namespace
2213
2214}
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:180
Cantera::NotImplementedError
An error indicating that an unimplemented function has been called.
Definition ctexceptions.h:209
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:584
Cantera::Storage::writeAttributes
void writeAttributes(const string &id, const AnyMap &meta)
Write attributes to a specified location.
Definition Storage.cpp:602
Cantera::Storage::deleteGroup
void deleteGroup(const string &id)
Delete group.
Definition Storage.cpp:578
Cantera::Storage::readAttributes
AnyMap readAttributes(const string &id, bool recursive) const
Read attributes from a specified location.
Definition Storage.cpp:596
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:608
Cantera::Storage::checkGroup
bool checkGroup(const string &id, bool permissive=false)
Check whether path location exists.
Definition Storage.cpp:572
Cantera::Storage::setCompressionLevel
void setCompressionLevel(int level)
Set compression level (0..9)
Definition Storage.cpp:560
Cantera::Storage::writeData
void writeData(const string &id, const string &name, const AnyValue &data)
Write dataset to a specified location.
Definition Storage.cpp:615
Cantera::tokenizePath
vector< string > tokenizePath(const string &in_val)
This function separates a string up into tokens according to the location of path separators.
Definition stringUtils.cpp:198
Cantera::toLowerCopy
string toLowerCopy(const string &input)
Convert to lower case.
Definition stringUtils.cpp:222
Cantera::len
U len(const T &container)
Get the size of a container, cast to a signed integer type.
Definition utilities.h:231
Cantera::dot
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition utilities.h:96
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:263
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:183
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...