Skip to content

Fix 1D solver crash when loglevel > 7 (#2073) #2162

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
hazalkoom wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Fix 1D solver crash when loglevel > 7 (#2073) #2162

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions include/cantera/base/SolutionArray.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,12 @@ class SolutionArray
//! Set auxiliary data for a given location.
void setAuxiliary(int loc, const AnyMap& data);

//! Get skipThermo setting
bool skipThermo() const { return m_skipThermo; }

//! Set skipThermo setting
void setSkipThermo(bool skip) { m_skipThermo = skip; }

//! Append location entry at end of SolutionArray.
void append(const vector<double>& state, const AnyMap& extra);

Expand Down Expand Up @@ -418,6 +424,7 @@ class SolutionArray

bool m_shared = false; //!< `true` if data are shared from another object
vector<int> m_active; //!< Vector of locations referencing active entries
bool m_skipThermo = false; //!< `true` if thermodynamic calculations are bypassed
};

//! Return mapping of component alias names to standardized component names.
Expand Down
21 changes: 14 additions & 7 deletions include/cantera/oneD/Boundary1D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -176,7 +176,8 @@ class Inlet1D : public Boundary1D
void init() override;
void eval(size_t jg, span<const double> xg, span<double> rg, span<int> diagg,
double rdt) override;
shared_ptr<SolutionArray> toArray(bool normalize=false) override;
shared_ptr<SolutionArray> toArray(bool normalize=false,
bool skipThermo=false) override;
void fromArray(const shared_ptr<SolutionArray>& arr) override;

protected:
Expand Down Expand Up @@ -222,7 +223,8 @@ class Empty1D : public Boundary1D
void eval(size_t jg, span<const double> xg, span<double> rg, span<int> diagg,
double rdt) override;

shared_ptr<SolutionArray> toArray(bool normalize=false) override;
shared_ptr<SolutionArray> toArray(bool normalize=false,
bool skipThermo=false) override;
};


Expand Down Expand Up @@ -254,7 +256,8 @@ class Symm1D : public Boundary1D
void eval(size_t jg, span<const double> xg, span<double> rg, span<int> diagg,
double rdt) override;

shared_ptr<SolutionArray> toArray(bool normalize=false) override;
shared_ptr<SolutionArray> toArray(bool normalize=false,
bool skipThermo=false) override;
};


Expand Down Expand Up @@ -286,7 +289,8 @@ class Outlet1D : public Boundary1D
void eval(size_t jg, span<const double> xg, span<double> rg, span<int> diagg,
double rdt) override;

shared_ptr<SolutionArray> toArray(bool normalize=false) override;
shared_ptr<SolutionArray> toArray(bool normalize=false,
bool skipThermo=false) override;
};


Expand Down Expand Up @@ -324,7 +328,8 @@ class OutletRes1D : public Boundary1D
void init() override;
void eval(size_t jg, span<const double> xg, span<double> rg, span<int> diagg,
double rdt) override;
shared_ptr<SolutionArray> toArray(bool normalize=false) override;
shared_ptr<SolutionArray> toArray(bool normalize=false,
bool skipThermo=false) override;
void fromArray(const shared_ptr<SolutionArray>& arr) override;

protected:
Expand Down Expand Up @@ -362,7 +367,8 @@ class Surf1D : public Boundary1D
void init() override;
void eval(size_t jg, span<const double> xg, span<double> rg, span<int> diagg,
double rdt) override;
shared_ptr<SolutionArray> toArray(bool normalize=false) override;
shared_ptr<SolutionArray> toArray(bool normalize=false,
bool skipThermo=false) override;
void fromArray(const shared_ptr<SolutionArray>& arr) override;
void show(span<const double> x) override;
};
Expand Down Expand Up @@ -407,7 +413,8 @@ class ReactingSurf1D : public Boundary1D
double rdt) override;

double value(const string& component) const override;
shared_ptr<SolutionArray> toArray(bool normalize=false) override;
shared_ptr<SolutionArray> toArray(bool normalize=false,
bool skipThermo=false) override;
void fromArray(const shared_ptr<SolutionArray>& arr) override;

void _getInitialSoln(span<double> x) override {
Expand Down
6 changes: 5 additions & 1 deletion include/cantera/oneD/Domain1D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -519,10 +519,14 @@ class Domain1D
* This method serves as an external interface for high-level API's; it does not
* provide direct access to memory.
* @param normalize If true, normalize concentrations (default=false)
* @param skipThermo If true, skip evaluation of thermodynamic properties
* (default=false). Used when saving residuals, where solution data does
* not represent a valid thermodynamic state.
*
* @since New in %Cantera 3.0.
*/
virtual shared_ptr<SolutionArray> toArray(bool normalize=false) {
virtual shared_ptr<SolutionArray> toArray(bool normalize=false,
bool skipThermo=false) {
throw NotImplementedError("Domain1D::toArray", "Needs to be overloaded.");
}

Expand Down
3 changes: 2 additions & 1 deletion include/cantera/oneD/Flow1D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,8 @@ class Flow1D : public Domain1D
span<const double> values) override;
void setFlatProfile(const string& component, double value) override;

shared_ptr<SolutionArray> toArray(bool normalize=false) override;
shared_ptr<SolutionArray> toArray(bool normalize=false,
bool skipThermo=false) override;
void fromArray(const shared_ptr<SolutionArray>& arr) override;

//! Set flow configuration for freely-propagating flames, using an internal point
Expand Down
5 changes: 5 additions & 0 deletions include/cantera/oneD/Sim1D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -372,6 +372,11 @@ class Sim1D : public OneDim
int m_ts_regrid_max = 3;

private:
//! Internal save with option to skip thermodynamic evaluation.
void save(const string& fname, const string& name, const string& desc,
bool overwrite, int compression, const string& basis,
bool skipThermo);

//! Calls method _finalize in each domain.
void finalize();
};
Expand Down
42 changes: 29 additions & 13 deletions src/base/SolutionArray.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1279,27 +1279,43 @@ void SolutionArray::writeEntry(AnyMap& root, const string& name, const string& s
data["components"] = componentNames();
}

auto phase = m_sol->thermo();
for (auto& [_, key] : *m_order) {
if (m_skipThermo && phase->speciesIndex(key, false) != npos) {
continue;
}
data[key] = m_extra->at(key);
}

auto phase = m_sol->thermo();
if (m_size == 1) {
setLoc(0);
data["temperature"] = phase->temperature();
data["pressure"] = phase->pressure();
auto surf = std::dynamic_pointer_cast<SurfPhase>(phase);
auto nSpecies = phase->nSpecies();
vector<double> values(nSpecies);
if (surf) {
surf->getCoverages(values);
} else {
phase->getMassFractions(values);
}
AnyMap items;
for (size_t k = 0; k < nSpecies; k++) {
if (values[k] != 0.0) {
items[phase->speciesName(k)] = values[k];
if (!m_skipThermo) {
setLoc(0);
data["temperature"] = phase->temperature();
data["pressure"] = phase->pressure();
vector<double> values(nSpecies);
if (surf) {
surf->getCoverages(values);
} else {
phase->getMassFractions(values);
}
for (size_t k = 0; k < nSpecies; k++) {
if (values[k] != 0.0) {
items[phase->speciesName(k)] = values[k];
}
}
} else {
setLoc(0, false);
data["temperature"] = phase->temperature();
data["pressure"] = phase->pressure();
for (size_t k = 0; k < nSpecies; k++) {
string name = phase->speciesName(k);
if (m_extra->count(name)) {
auto& val_vec = m_extra->at(name).asVector<double>();
items[name] = val_vec[0];
}
}
}
if (surf) {
Expand Down
45 changes: 29 additions & 16 deletions src/oneD/Boundary1D.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,7 @@ void Inlet1D::eval(size_t jg, span<const double> xg, span<double> rg,
}
}

shared_ptr<SolutionArray> Inlet1D::toArray(bool normalize)
shared_ptr<SolutionArray> Inlet1D::toArray(bool normalize, bool skipThermo)
{
AnyMap meta = Boundary1D::getMeta();
meta["mass-flux"] = m_mdot;
Expand Down Expand Up @@ -334,7 +334,7 @@ void Empty1D::eval(size_t jg, span<const double> xg, span<double> rg,
{
}

shared_ptr<SolutionArray> Empty1D::toArray(bool normalize)
shared_ptr<SolutionArray> Empty1D::toArray(bool normalize, bool skipThermo)
{
AnyMap meta = Boundary1D::getMeta();
return SolutionArray::create(m_solution, 0, meta);
Expand Down Expand Up @@ -382,7 +382,7 @@ void Symm1D::eval(size_t jg, span<const double> xg, span<double> rg,
}
}

shared_ptr<SolutionArray> Symm1D::toArray(bool normalize)
shared_ptr<SolutionArray> Symm1D::toArray(bool normalize, bool skipThermo)
{
AnyMap meta = Boundary1D::getMeta();
return SolutionArray::create(m_solution, 0, meta);
Expand Down Expand Up @@ -446,7 +446,7 @@ void Outlet1D::eval(size_t jg, span<const double> xg, span<double> rg,
}
}

shared_ptr<SolutionArray> Outlet1D::toArray(bool normalize)
shared_ptr<SolutionArray> Outlet1D::toArray(bool normalize, bool skipThermo)
{
AnyMap meta = Boundary1D::getMeta();
return SolutionArray::create(m_solution, 0, meta);
Expand Down Expand Up @@ -525,7 +525,7 @@ void OutletRes1D::eval(size_t jg, span<const double> xg, span<double> rg,
}
}

shared_ptr<SolutionArray> OutletRes1D::toArray(bool normalize)
shared_ptr<SolutionArray> OutletRes1D::toArray(bool normalize, bool skipThermo)
{
AnyMap meta = Boundary1D::getMeta();
meta["temperature"] = m_temp;
Expand Down Expand Up @@ -583,7 +583,7 @@ void Surf1D::eval(size_t jg, span<const double> xg, span<double> rg,
}
}

shared_ptr<SolutionArray> Surf1D::toArray(bool normalize)
shared_ptr<SolutionArray> Surf1D::toArray(bool normalize, bool skipThermo)
{
AnyMap meta = Boundary1D::getMeta();
meta["temperature"] = m_temp;
Expand Down Expand Up @@ -760,7 +760,7 @@ double ReactingSurf1D::value(const string& component) const
return soln[index(i, 0)];
}

shared_ptr<SolutionArray> ReactingSurf1D::toArray(bool normalize)
shared_ptr<SolutionArray> ReactingSurf1D::toArray(bool normalize, bool skipThermo)
{
if (!m_state) {
throw CanteraError("ReactingSurf1D::toArray",
Expand All @@ -773,16 +773,29 @@ shared_ptr<SolutionArray> ReactingSurf1D::toArray(bool normalize)
AnyValue source = m_sphase->input().getMetadata("filename");
meta["phase"]["source"] = source.empty() ? "<unknown>" : source.asString();

// set state of surface phase
m_sphase->setState_TP(m_temp, m_sphase->pressure());
m_sphase->setCoverages(span<const double>(soln, m_nsp));
vector<double> data(m_sphase->stateSize());
m_sphase->saveState(data);

auto arr = SolutionArray::create(m_solution, 1, meta);
arr->setState(0, data);
if (normalize) {
arr->normalize();
if (!skipThermo) {
// set state of surface phase
m_sphase->setState_TP(m_temp, m_sphase->pressure());
m_sphase->setCoverages(span<const double>(soln, m_nsp));
vector<double> data(m_sphase->stateSize());
m_sphase->saveState(data);
arr->setState(0, data);
if (normalize) {
arr->normalize();
}
} else {
vector<double> val(1);
AnyValue value;
for (size_t k = 0; k < m_nsp; k++) {
auto name = componentName(k);
if (!arr->hasComponent(name)) {
arr->addExtra(name);
}
val[0] = soln[k];
value = val;
arr->setComponent(name, value);
}
}
return arr;
}
Expand Down
12 changes: 7 additions & 5 deletions src/oneD/Flow1D.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1069,7 +1069,7 @@ void Flow1D::setFlatProfile(const string& component, double value)
}
}

shared_ptr<SolutionArray> Flow1D::toArray(bool normalize)
shared_ptr<SolutionArray> Flow1D::toArray(bool normalize, bool skipThermo)
{
if (!m_state) {
throw CanteraError("Flow1D::toArray",
Expand All @@ -1096,17 +1096,19 @@ shared_ptr<SolutionArray> Flow1D::toArray(bool normalize)
arr->setComponent(name, value);
}
}
updateThermo(span<const double>(soln, size()), 0, m_points-1);
value = m_rho;
arr->setComponent("D", value); // use density rather than pressure
if (!skipThermo) {
updateThermo(span<const double>(soln, size()), 0, m_points-1);
value = m_rho;
arr->setComponent("D", value); // use density rather than pressure
}

if (m_do_radiation) {
arr->addExtra("radiative-heat-loss", true); // add at end
value = m_qdotRadiation;
arr->setComponent("radiative-heat-loss", value);
}

if (normalize) {
if (normalize && !skipThermo) {
arr->normalize();
}
return arr;
Expand Down
20 changes: 15 additions & 5 deletions src/oneD/Sim1D.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,12 +61,19 @@ double Sim1D::_workValue(size_t dom, size_t comp, size_t localPoint) const

void Sim1D::save(const string& fname, const string& name, const string& desc,
bool overwrite, int compression, const string& basis)
{
save(fname, name, desc, overwrite, compression, basis, false);
}

void Sim1D::save(const string& fname, const string& name, const string& desc,
bool overwrite, int compression, const string& basis,
bool skipThermo)
{
size_t dot = fname.find_last_of(".");
string extension = (dot != npos) ? toLowerCopy(fname.substr(dot+1)) : "";
if (extension == "csv") {
for (auto dom : m_dom) {
auto arr = dom->toArray();
auto arr = dom->toArray(false, skipThermo);
if (dom->size() > 1) {
arr->writeEntry(fname, overwrite, basis);
break;
Expand All @@ -81,7 +88,8 @@ void Sim1D::save(const string& fname, const string& name, const string& desc,
if (extension == "h5" || extension == "hdf" || extension == "hdf5") {
SolutionArray::writeHeader(fname, name, desc, overwrite);
for (auto dom : m_dom) {
auto arr = dom->toArray();
auto arr = dom->toArray(false, skipThermo);
arr->setSkipThermo(skipThermo);
arr->writeEntry(fname, name, dom->id(), overwrite, compression);
}
return;
Expand All @@ -95,7 +103,8 @@ void Sim1D::save(const string& fname, const string& name, const string& desc,
SolutionArray::writeHeader(data, name, desc, overwrite);

for (auto dom : m_dom) {
auto arr = dom->toArray();
auto arr = dom->toArray(false, skipThermo);
arr->setSkipThermo(skipThermo);
arr->writeEntry(data, name, dom->id(), overwrite);
}

Expand All @@ -114,10 +123,11 @@ void Sim1D::saveResidual(const string& fname, const string& name,
vector<double> res(m_state->size(), -999);
OneDim::eval(npos, *m_state, res, 0.0);
// Temporarily put the residual into m_state, since this is the vector that the
// save() function reads.
// toArray() functions read. Pass skipThermo=true since residual values
// are not valid thermodynamic states.
vector<double> backup(*m_state);
*m_state = res;
save(fname, name, desc, overwrite, compression);
save(fname, name, desc, overwrite, compression, "", true);
*m_state = backup;
}

Expand Down
Loading
Loading