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