diff --git a/include/cantera/base/SolutionArray.h b/include/cantera/base/SolutionArray.h index 8e6a6b8514..51fc725fa9 100644 --- a/include/cantera/base/SolutionArray.h +++ b/include/cantera/base/SolutionArray.h @@ -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& state, const AnyMap& extra); @@ -418,6 +424,7 @@ class SolutionArray bool m_shared = false; //!< `true` if data are shared from another object vector 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. diff --git a/include/cantera/oneD/Boundary1D.h b/include/cantera/oneD/Boundary1D.h index 6d7543d77f..3224fb998d 100644 --- a/include/cantera/oneD/Boundary1D.h +++ b/include/cantera/oneD/Boundary1D.h @@ -176,7 +176,8 @@ class Inlet1D : public Boundary1D void init() override; void eval(size_t jg, span xg, span rg, span diagg, double rdt) override; - shared_ptr toArray(bool normalize=false) override; + shared_ptr toArray(bool normalize=false, + bool skipThermo=false) override; void fromArray(const shared_ptr& arr) override; protected: @@ -222,7 +223,8 @@ class Empty1D : public Boundary1D void eval(size_t jg, span xg, span rg, span diagg, double rdt) override; - shared_ptr toArray(bool normalize=false) override; + shared_ptr toArray(bool normalize=false, + bool skipThermo=false) override; }; @@ -254,7 +256,8 @@ class Symm1D : public Boundary1D void eval(size_t jg, span xg, span rg, span diagg, double rdt) override; - shared_ptr toArray(bool normalize=false) override; + shared_ptr toArray(bool normalize=false, + bool skipThermo=false) override; }; @@ -286,7 +289,8 @@ class Outlet1D : public Boundary1D void eval(size_t jg, span xg, span rg, span diagg, double rdt) override; - shared_ptr toArray(bool normalize=false) override; + shared_ptr toArray(bool normalize=false, + bool skipThermo=false) override; }; @@ -324,7 +328,8 @@ class OutletRes1D : public Boundary1D void init() override; void eval(size_t jg, span xg, span rg, span diagg, double rdt) override; - shared_ptr toArray(bool normalize=false) override; + shared_ptr toArray(bool normalize=false, + bool skipThermo=false) override; void fromArray(const shared_ptr& arr) override; protected: @@ -362,7 +367,8 @@ class Surf1D : public Boundary1D void init() override; void eval(size_t jg, span xg, span rg, span diagg, double rdt) override; - shared_ptr toArray(bool normalize=false) override; + shared_ptr toArray(bool normalize=false, + bool skipThermo=false) override; void fromArray(const shared_ptr& arr) override; void show(span x) override; }; @@ -407,7 +413,8 @@ class ReactingSurf1D : public Boundary1D double rdt) override; double value(const string& component) const override; - shared_ptr toArray(bool normalize=false) override; + shared_ptr toArray(bool normalize=false, + bool skipThermo=false) override; void fromArray(const shared_ptr& arr) override; void _getInitialSoln(span x) override { diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index 0ec95c74f7..b85aa9d76c 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -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 toArray(bool normalize=false) { + virtual shared_ptr toArray(bool normalize=false, + bool skipThermo=false) { throw NotImplementedError("Domain1D::toArray", "Needs to be overloaded."); } diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index 4cfb48ca06..4bdee634e5 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -184,7 +184,8 @@ class Flow1D : public Domain1D span values) override; void setFlatProfile(const string& component, double value) override; - shared_ptr toArray(bool normalize=false) override; + shared_ptr toArray(bool normalize=false, + bool skipThermo=false) override; void fromArray(const shared_ptr& arr) override; //! Set flow configuration for freely-propagating flames, using an internal point diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index b99c507cc1..233b007d2d 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/include/cantera/oneD/Sim1D.h @@ -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(); }; diff --git a/src/base/SolutionArray.cpp b/src/base/SolutionArray.cpp index fc6baf7e5e..ba10dda0b8 100644 --- a/src/base/SolutionArray.cpp +++ b/src/base/SolutionArray.cpp @@ -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(phase); auto nSpecies = phase->nSpecies(); - vector 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 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(); + items[name] = val_vec[0]; + } } } if (surf) { diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp index 07b674b7ab..3f01f5a2e9 100644 --- a/src/oneD/Boundary1D.cpp +++ b/src/oneD/Boundary1D.cpp @@ -285,7 +285,7 @@ void Inlet1D::eval(size_t jg, span xg, span rg, } } -shared_ptr Inlet1D::toArray(bool normalize) +shared_ptr Inlet1D::toArray(bool normalize, bool skipThermo) { AnyMap meta = Boundary1D::getMeta(); meta["mass-flux"] = m_mdot; @@ -334,7 +334,7 @@ void Empty1D::eval(size_t jg, span xg, span rg, { } -shared_ptr Empty1D::toArray(bool normalize) +shared_ptr Empty1D::toArray(bool normalize, bool skipThermo) { AnyMap meta = Boundary1D::getMeta(); return SolutionArray::create(m_solution, 0, meta); @@ -382,7 +382,7 @@ void Symm1D::eval(size_t jg, span xg, span rg, } } -shared_ptr Symm1D::toArray(bool normalize) +shared_ptr Symm1D::toArray(bool normalize, bool skipThermo) { AnyMap meta = Boundary1D::getMeta(); return SolutionArray::create(m_solution, 0, meta); @@ -446,7 +446,7 @@ void Outlet1D::eval(size_t jg, span xg, span rg, } } -shared_ptr Outlet1D::toArray(bool normalize) +shared_ptr Outlet1D::toArray(bool normalize, bool skipThermo) { AnyMap meta = Boundary1D::getMeta(); return SolutionArray::create(m_solution, 0, meta); @@ -525,7 +525,7 @@ void OutletRes1D::eval(size_t jg, span xg, span rg, } } -shared_ptr OutletRes1D::toArray(bool normalize) +shared_ptr OutletRes1D::toArray(bool normalize, bool skipThermo) { AnyMap meta = Boundary1D::getMeta(); meta["temperature"] = m_temp; @@ -583,7 +583,7 @@ void Surf1D::eval(size_t jg, span xg, span rg, } } -shared_ptr Surf1D::toArray(bool normalize) +shared_ptr Surf1D::toArray(bool normalize, bool skipThermo) { AnyMap meta = Boundary1D::getMeta(); meta["temperature"] = m_temp; @@ -760,7 +760,7 @@ double ReactingSurf1D::value(const string& component) const return soln[index(i, 0)]; } -shared_ptr ReactingSurf1D::toArray(bool normalize) +shared_ptr ReactingSurf1D::toArray(bool normalize, bool skipThermo) { if (!m_state) { throw CanteraError("ReactingSurf1D::toArray", @@ -773,16 +773,29 @@ shared_ptr ReactingSurf1D::toArray(bool normalize) AnyValue source = m_sphase->input().getMetadata("filename"); meta["phase"]["source"] = source.empty() ? "" : source.asString(); - // set state of surface phase - m_sphase->setState_TP(m_temp, m_sphase->pressure()); - m_sphase->setCoverages(span(soln, m_nsp)); - vector 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(soln, m_nsp)); + vector data(m_sphase->stateSize()); + m_sphase->saveState(data); + arr->setState(0, data); + if (normalize) { + arr->normalize(); + } + } else { + vector 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; } diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 9b68b9b4ab..0d5e977574 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -1069,7 +1069,7 @@ void Flow1D::setFlatProfile(const string& component, double value) } } -shared_ptr Flow1D::toArray(bool normalize) +shared_ptr Flow1D::toArray(bool normalize, bool skipThermo) { if (!m_state) { throw CanteraError("Flow1D::toArray", @@ -1096,9 +1096,11 @@ shared_ptr Flow1D::toArray(bool normalize) arr->setComponent(name, value); } } - updateThermo(span(soln, size()), 0, m_points-1); - value = m_rho; - arr->setComponent("D", value); // use density rather than pressure + if (!skipThermo) { + updateThermo(span(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 @@ -1106,7 +1108,7 @@ shared_ptr Flow1D::toArray(bool normalize) arr->setComponent("radiative-heat-loss", value); } - if (normalize) { + if (normalize && !skipThermo) { arr->normalize(); } return arr; diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index 80402b2926..7b77703582 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -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; @@ -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; @@ -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); } @@ -114,10 +123,11 @@ void Sim1D::saveResidual(const string& fname, const string& name, vector 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 backup(*m_state); *m_state = res; - save(fname, name, desc, overwrite, compression); + save(fname, name, desc, overwrite, compression, "", true); *m_state = backup; } diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index 9c97fefe91..a4562a0852 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -1004,6 +1004,16 @@ def test_save_restore_hdf(self): def test_save_restore_yaml(self): self.run_save_restore("yaml") + @pytest.mark.parametrize("loglevel", [7, 8]) + def test_solve_loglevel_residual(self, monkeypatch, loglevel): + # Test high loglevels (7=solution, 8=residual) to prevent regressions (#2073) + monkeypatch.chdir(self.test_work_path) + self.run_mix(phi=1.0, T=300, width=0.05, p=1.0, refine=False) + self.sim.solve(loglevel=loglevel, refine_grid=False) + debug_file = self.test_work_path / "debug_sim1d.yaml" + assert debug_file.exists() + debug_file.unlink() + def run_save_restore(self, mode): filename = self.test_work_path / f"freeflame.{mode}" filename.unlink(missing_ok=True) @@ -2136,6 +2146,20 @@ def run_reacting_surface(self, xch4, tsurf, mdot, width): def test_reacting_surface_case1(self): self.run_reacting_surface(xch4=0.095, tsurf=900.0, mdot=0.06, width=0.1) + @pytest.mark.parametrize("loglevel", [7, 8]) + def test_reacting_surface_loglevel_residual(self, monkeypatch, loglevel): + monkeypatch.chdir(self.test_work_path) + comp = {'CH4': 0.095, 'O2': 0.21, 'N2': 0.79} + sim = self.create_reacting_surface(comp, tsurf=900.0, tinlet=300.0, width=0.1) + sim.inlet.mdot = 0.06 + sim.inlet.T = 300.0 + sim.inlet.X = comp + sim.surface.T = 900.0 + sim.solve(loglevel=loglevel, refine_grid=False) + debug_file = self.test_work_path / "debug_sim1d.yaml" + assert debug_file.exists() + debug_file.unlink() + @pytest.mark.slow_test def test_reacting_surface_case2(self): self.run_reacting_surface(xch4=0.07, tsurf=1200.0, mdot=0.2, width=0.05)