From 1ff297ec3b04218c513ca5e1a0d84b922de0579e Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 4 Oct 2025 09:41:14 -0500 Subject: [PATCH 01/22] [oneD] Rename 'lambda' to 'Lambda' --- include/cantera/oneD/Flow1D.h | 15 +++++--- interfaces/cython/cantera/onedim.py | 10 +++--- .../python/onedim/diffusion_flame_batch.py | 4 +-- .../onedim/diffusion_flame_extinction.py | 2 +- src/oneD/Boundary1D.cpp | 2 +- src/oneD/Flow1D.cpp | 34 ++++++++++++------- .../cxx_samples/cxx_flamespeed_blessed.txt | 2 +- 7 files changed, 42 insertions(+), 27 deletions(-) diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index 5eaeb444683..ee416547bf2 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -548,15 +548,15 @@ class Flow1D : public Domain1D double rdt, size_t jmin, size_t jmax); /** - * Evaluate the lambda equation residual. + * Evaluate the radial pressure gradient equation residual. * * @f[ * \frac{d\Lambda}{dz} = 0 * @f] * - * The lambda equation serves as an eigenvalue that allows the momentum + * The Lambda equation serves as an eigenvalue that allows the momentum * equation and continuity equations to be simultaneously satisfied in - * axisymmetric flows. The lambda equation propagates information from + * axisymmetric flows. The Lambda equation propagates information from * left-to-right. The default boundary condition is @f$ \Lambda = 0 @f$ * at the left boundary. The equation is first order and so only one * boundary condition is needed. @@ -571,7 +571,7 @@ class Flow1D : public Domain1D * * @f[ * \rho c_p u \frac{dT}{dz} = - * \frac{d}{dz}\left( \lambda \frac{dT}{dz} \right) + * \frac{d}{dz}\left( \Lambda \frac{dT}{dz} \right) * - \sum_k h_kW_k\dot{\omega}_k * - \sum_k j_k \frac{dh_k}{dz} * @f] @@ -674,7 +674,12 @@ class Flow1D : public Domain1D //! Get the radial pressure gradient [N/m⁴] at point `j` from the local state vector //! `x` - double lambda(const double* x, size_t j) const { + //! @deprecated To be removed after %Cantera 3.2. Renamed to Lambda(). + double lambda(const double* x, size_t j) const; + + //! Get the radial pressure gradient [N/m⁴] at point `j` from the local state vector + //! `x` + double Lambda(const double* x, size_t j) const { return x[index(c_offset_L, j)]; } diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index b5467776c3f..0353cb70cfd 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -292,9 +292,9 @@ def spread_rate(self): def L(self): """ Array containing the radial pressure gradient (1/r)(dP/dr) [N/m⁴] at - each point. Note: This value is named 'lambda' in the C++ code. + each point. Note: This value is named `Lambda` in the C++ code. """ - return self.profile(self.flame, 'lambda') + return self.profile(self.flame, "Lambda") @property def E(self): @@ -1047,7 +1047,7 @@ def set_initial_guess(self, data=None, group=None): self.set_profile('velocity', [0.0, 1.0], [u0f, -u0o]) self.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) - self.set_profile("lambda", [0.0, 1.0], [L, L]) + self.set_profile("Lambda", [0.0, 1.0], [L, L]) self.set_profile('T', zrel, T) for k,spec in enumerate(self.gas.species_names): self.set_profile(spec, zrel, Y[:,k]) @@ -1455,7 +1455,7 @@ def set_initial_guess(self, equilibrate=True, data=None, group=None): self.set_profile('velocity', [0.0, 1.0], [uu, -ub]) self.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) - self.set_profile("lambda", [0.0, 1.0], [L, L]) + self.set_profile("Lambda", [0.0, 1.0], [L, L]) class CounterflowTwinPremixedFlame(FlameBase): @@ -1539,4 +1539,4 @@ def set_initial_guess(self, data=None, group=None): self.set_profile('velocity', [0.0, 1.0], [uu, 0]) self.set_profile('spread_rate', [0.0, 1.0], [0.0, a]) - self.set_profile("lambda", [0.0, 1.0], [L, L]) + self.set_profile("Lambda", [0.0, 1.0], [L, L]) diff --git a/samples/python/onedim/diffusion_flame_batch.py b/samples/python/onedim/diffusion_flame_batch.py index 58dda388b1c..224ce3d098a 100644 --- a/samples/python/onedim/diffusion_flame_batch.py +++ b/samples/python/onedim/diffusion_flame_batch.py @@ -141,7 +141,7 @@ def names(test): f.set_profile('spread_rate', normalized_grid, f.spread_rate * rel_pressure_increase ** exp_V_p) # Update pressure curvature - f.set_profile('lambda', normalized_grid, + f.set_profile("Lambda", normalized_grid, f.L * rel_pressure_increase ** exp_lam_p) try: @@ -196,7 +196,7 @@ def names(test): f.set_profile('spread_rate', normalized_grid, f.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature - f.set_profile('lambda', normalized_grid, f.L * strain_factor ** exp_lam_a) + f.set_profile("Lambda", normalized_grid, f.L * strain_factor ** exp_lam_a) try: # Try solving the flame f.solve(loglevel=0) diff --git a/samples/python/onedim/diffusion_flame_extinction.py b/samples/python/onedim/diffusion_flame_extinction.py index 43d09603b00..d17b2fc204c 100644 --- a/samples/python/onedim/diffusion_flame_extinction.py +++ b/samples/python/onedim/diffusion_flame_extinction.py @@ -134,7 +134,7 @@ def names(test): f.set_profile('spread_rate', normalized_grid, f.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature - f.set_profile('lambda', normalized_grid, f.L * strain_factor ** exp_lam_a) + f.set_profile("Lambda", normalized_grid, f.L * strain_factor ** exp_lam_a) try: f.solve(loglevel=0) except ct.CanteraError as e: diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp index b4a41b3976c..7d387559084 100644 --- a/src/oneD/Boundary1D.cpp +++ b/src/oneD/Boundary1D.cpp @@ -165,7 +165,7 @@ void Inlet1D::init() throw CanteraError("Inlet1D::init", "Inlet1D is not properly connected."); } - // components = u, V, T, lambda, + mass fractions + // components = u, V, T, Lambda, + mass fractions m_nsp = m_flow->phase().nSpecies(); m_yin.resize(m_nsp, 0.0); if (m_xstr != "") { diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index a9943682a4c..59bf3aa2fe9 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -60,7 +60,7 @@ Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) : setBounds(c_offset_U, -1e20, 1e20); // no bounds on u setBounds(c_offset_V, -1e20, 1e20); // no bounds on V setBounds(c_offset_T, 200.0, 2*m_thermo->maxTemp()); // temperature bounds - setBounds(c_offset_L, -1e20, 1e20); // lambda should be negative + setBounds(c_offset_L, -1e20, 1e20); // Lambda should be negative setBounds(c_offset_E, -1e20, 1e20); // no bounds on electric field setBounds(c_offset_Uo, -1e20, 1e20); // no bounds on Uo // mass fraction bounds @@ -120,6 +120,12 @@ Flow1D::~Flow1D() } } +double Flow1D::lambda(const double* x, size_t j) const { + warn_deprecated("Flow1D::lambda", + "To be removed after Cantera 3.2. Component 'lambda' is renamed to 'Lambda'."); + return x[index(c_offset_L, j)]; +} + string Flow1D::domainType() const { if (m_isFree) { return "free-flow"; @@ -546,7 +552,7 @@ void Flow1D::evalContinuity(double* x, double* rsd, int* diag, for (size_t j = j0; j <= j1; j++) { // interior points // For "axisymmetric-flow", the continuity equation propagates the // mass flow rate information to the left (j+1 -> j) from the value - // specified at the right boundary. The lambda information propagates + // specified at the right boundary. The Lambda information propagates // in the opposite direction. rsd[index(c_offset_U, j)] = -(rho_u(x, j+1) - rho_u(x, j))/m_dz[j] -(density(j+1)*V(x, j+1) + density(j)*V(x, j)); @@ -600,7 +606,7 @@ void Flow1D::evalMomentum(double* x, double* rsd, int* diag, size_t j0 = std::max(jmin, 1); size_t j1 = std::min(jmax, m_points-2); for (size_t j = j0; j <= j1; j++) { // interior points - rsd[index(c_offset_V, j)] = (shear(x, j) - lambda(x, j) + rsd[index(c_offset_V, j)] = (shear(x, j) - Lambda(x, j) - rho_u(x, j) * dVdz(x, j) - m_rho[j] * V(x, j) * V(x, j)) / m_rho[j]; if (!m_twoPointControl) { @@ -617,7 +623,7 @@ void Flow1D::evalLambda(double* x, double* rsd, int* diag, { if (!m_usesLambda) { // disable this equation for (size_t j = jmin; j <= jmax; j++) { - rsd[index(c_offset_L, j)] = lambda(x, j); + rsd[index(c_offset_L, j)] = Lambda(x, j); diag[index(c_offset_L, j)] = 0; } return; @@ -625,14 +631,14 @@ void Flow1D::evalLambda(double* x, double* rsd, int* diag, if (jmin == 0) { // left boundary if (m_twoPointControl) { - rsd[index(c_offset_L, jmin)] = lambda(x, jmin+1) - lambda(x, jmin); + rsd[index(c_offset_L, jmin)] = Lambda(x, jmin+1) - Lambda(x, jmin); } else { rsd[index(c_offset_L, jmin)] = -rho_u(x, jmin); } } if (jmax == m_points - 1) { // right boundary - rsd[index(c_offset_L, jmax)] = lambda(x, jmax) - lambda(x, jmax-1); + rsd[index(c_offset_L, jmax)] = Lambda(x, jmax) - Lambda(x, jmax-1); diag[index(c_offset_L, jmax)] = 0; } @@ -644,12 +650,12 @@ void Flow1D::evalLambda(double* x, double* rsd, int* diag, if (z(j) == m_zLeft) { rsd[index(c_offset_L, j)] = T(x,j) - m_tLeft; } else if (z(j) > m_zLeft) { - rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1); + rsd[index(c_offset_L, j)] = Lambda(x, j) - Lambda(x, j-1); } else if (z(j) < m_zLeft) { - rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j+1); + rsd[index(c_offset_L, j)] = Lambda(x, j) - Lambda(x, j+1); } } else { - rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1); + rsd[index(c_offset_L, j)] = Lambda(x, j) - Lambda(x, j-1); } diag[index(c_offset_L, j)] = 0; } @@ -813,7 +819,7 @@ string Flow1D::componentName(size_t n) const case c_offset_T: return "T"; case c_offset_L: - return "lambda"; + return "Lambda"; case c_offset_E: return "eField"; case c_offset_Uo: @@ -835,12 +841,16 @@ size_t Flow1D::componentIndex(const string& name) const return c_offset_V; } else if (name=="T") { return c_offset_T; - } else if (name=="lambda") { + } else if (name=="Lambda") { return c_offset_L; } else if (name == "eField") { return c_offset_E; } else if (name == "Uo") { return c_offset_Uo; + } else if (name=="lambda") { + warn_deprecated("Flow1D::componentIndex", + "Component 'lambda' is renamed to 'Lambda'."); + return c_offset_L; } else { for (size_t n=c_offset_Y; n Date: Sat, 4 Oct 2025 10:19:54 -0500 Subject: [PATCH 02/22] [oneD] Edge case for importing 'lambda' --- src/oneD/Flow1D.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 59bf3aa2fe9..d2f10d95482 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -984,6 +984,12 @@ void Flow1D::fromArray(SolutionArray& arr, double* soln) for (size_t j = 0; j < nPoints(); j++) { soln[index(i,j)] = data[j]; } + } else if (name == "Lambda" && arr.hasComponent("lambda")) { + // edge case: 'lambda' is renamed to 'Lambda' in Cantera 3.2 + const vector data = arr.getComponent("lambda").as>(); + for (size_t j = 0; j < nPoints(); j++) { + soln[index(i,j)] = data[j]; + } } else { warn_user("Flow1D::fromArray", "Saved state does not contain values for " "component '{}' in domain '{}'.", name, id()); From 8582bc3713464d5e0d8b8cee0526f2f2dda782af Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 4 Oct 2025 16:18:58 -0500 Subject: [PATCH 03/22] [SolutionArray] Enable data path for HDF5 files --- src/base/Storage.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/base/Storage.cpp b/src/base/Storage.cpp index 564853480fe..9abfef9b25a 100644 --- a/src/base/Storage.cpp +++ b/src/base/Storage.cpp @@ -7,6 +7,7 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/base/AnyMap.h" +#include "cantera/base/global.h" #include "cantera/base/Storage.h" #if CT_USE_HDF5 @@ -59,7 +60,8 @@ Storage::Storage(string fname, bool write) : m_write(write) if (m_write) { m_file = make_unique(fname, h5::File::OpenOrCreate); } else { - m_file = make_unique(fname, h5::File::ReadOnly); + auto fullName = findInputFile(fname); + m_file = make_unique(fullName, h5::File::ReadOnly); } } From 513ca17e4dbfc52093a9adf9d100815b56514025 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 4 Oct 2025 13:16:21 -0500 Subject: [PATCH 04/22] [SolutionArray] Implement component alias names --- interfaces/cython/cantera/composite.py | 6 +++++ src/base/SolutionArray.cpp | 32 ++++++++++++++++++++++---- 2 files changed, 33 insertions(+), 5 deletions(-) diff --git a/interfaces/cython/cantera/composite.py b/interfaces/cython/cantera/composite.py index 86ffd1ffca2..1e62e4c4d2f 100644 --- a/interfaces/cython/cantera/composite.py +++ b/interfaces/cython/cantera/composite.py @@ -677,6 +677,9 @@ def __getitem__(self, index): return out def __getattr__(self, name): + if (not self._has_component(name) and + self._has_component(name.replace("_", "-"))): + name = name.replace("_", "-") if self._has_component(name): out = self._get_component(name) out.setflags(write=False) @@ -688,6 +691,9 @@ def __getattr__(self, name): f"{self.__class__.__name__!r} object has no attribute '{name}'") def __setattr__(self, name, value): + if (not self._has_component(name) and + self._has_component(name.replace("_", "-"))): + name = name.replace("_", "-") if self._has_extra(name): if not self.shape: # scalar diff --git a/src/base/SolutionArray.cpp b/src/base/SolutionArray.cpp index d501d0b66ec..183893f8980 100644 --- a/src/base/SolutionArray.cpp +++ b/src/base/SolutionArray.cpp @@ -35,6 +35,21 @@ const std::map aliasMap = { {"Q", "vapor-fraction"}, }; +const std::map reverseAliasMap = { + // reserved names used by states + {"temperature", "T"}, + {"pressure", "P"}, + {"density", "D"}, + // reserved names used for 1-D objects + {"spread-rate", "spread_rate"}, + {"L", "Lambda"}, + {"radial-pressure-gradient", "Lambda"}, + {"lambda", "Lambda"}, // lower-case version used prior to Cantera 3.2 + {"E", "eField"}, + {"electric-field", "eField"}, + {"oxidizer-velocity", "Uo"}, +}; + namespace Cantera { @@ -663,21 +678,28 @@ bool SolutionArray::hasComponent(const string& name) const // reserved names return false; } + if (reverseAliasMap.count(name)) { + // registered alias + return true; + } // native state return (m_sol->thermo()->nativeState().count(name)); } AnyValue SolutionArray::getComponent(const string& name) const { - if (!hasComponent(name)) { + string _name = name; + if (reverseAliasMap.count(name)) { + _name = reverseAliasMap.at(name); + } else if (!hasComponent(name)) { throw CanteraError("SolutionArray::getComponent", "Unknown component '{}'.", name); } AnyValue out; - if (m_extra->count(name)) { + if (m_extra->count(_name)) { // extra component - const auto& extra = m_extra->at(name); + const auto& extra = m_extra->at(_name); if (extra.is()) { return AnyValue(); } @@ -709,10 +731,10 @@ AnyValue SolutionArray::getComponent(const string& name) const // component is part of state information vector data(m_size); - size_t ix = m_sol->thermo()->speciesIndex(name); + size_t ix = m_sol->thermo()->speciesIndex(_name); if (ix == npos) { // state other than species - ix = m_sol->thermo()->nativeState()[name]; + ix = m_sol->thermo()->nativeState()[_name]; } else { // species information ix += m_stride - m_sol->thermo()->nSpecies(); From c698e66f13d6795d1f87cf29917d852afed0ea4c Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 4 Oct 2025 16:35:15 -0500 Subject: [PATCH 05/22] [SolutionArray] Implement legacy import --- include/cantera/base/SolutionArray.h | 2 +- src/base/SolutionArray.cpp | 44 +++++++++++++++++++--------- 2 files changed, 31 insertions(+), 15 deletions(-) diff --git a/include/cantera/base/SolutionArray.h b/include/cantera/base/SolutionArray.h index e0e5229aada..564a0b112ae 100644 --- a/include/cantera/base/SolutionArray.h +++ b/include/cantera/base/SolutionArray.h @@ -130,7 +130,7 @@ class SolutionArray * Check whether SolutionArray contains a component. * A component is a property defining state or auxiliary variable. */ - bool hasComponent(const string& name) const; + bool hasComponent(const string& name, bool checkAlias=true) const; /** * Retrieve a component of the SolutionArray by name. diff --git a/src/base/SolutionArray.cpp b/src/base/SolutionArray.cpp index 183893f8980..d9a249b28fe 100644 --- a/src/base/SolutionArray.cpp +++ b/src/base/SolutionArray.cpp @@ -664,7 +664,7 @@ vector SolutionArray::listExtra(bool all) const return names; } -bool SolutionArray::hasComponent(const string& name) const +bool SolutionArray::hasComponent(const string& name, bool checkAlias) const { if (m_extra->count(name)) { // auxiliary data @@ -678,7 +678,7 @@ bool SolutionArray::hasComponent(const string& name) const // reserved names return false; } - if (reverseAliasMap.count(name)) { + if (checkAlias && reverseAliasMap.count(name)) { // registered alias return true; } @@ -1821,13 +1821,17 @@ void SolutionArray::readEntry(const string& fname, const string& name, const auto& components = m_meta["components"].asVector(); bool back = false; for (const auto& name : components) { - if (hasComponent(name) || name == "X" || name == "Y") { + if (hasComponent(name, false) || name == "X" || name == "Y") { back = true; } else { - addExtra(name, back); + auto _name = name; + if (reverseAliasMap.count(name)) { + _name = reverseAliasMap.at(name); + } + addExtra(_name, back); AnyValue data; data = file.readData(path, name, m_dataSize); - setComponent(name, data); + setComponent(_name, data); } } m_meta.erase("components"); @@ -1835,11 +1839,15 @@ void SolutionArray::readEntry(const string& fname, const string& name, // data format used by Python h5py export (Cantera 2.5) warn_user("SolutionArray::readEntry", "Detected legacy HDF format."); for (const auto& name : names) { - if (!hasComponent(name) && name != "X" && name != "Y") { - addExtra(name); + if (!hasComponent(name, false) && name != "X" && name != "Y") { + auto _name = name; + if (reverseAliasMap.count(name)) { + _name = reverseAliasMap.at(name); + } + addExtra(_name); AnyValue data; data = file.readData(path, name, m_dataSize); - setComponent(name, data); + setComponent(_name, data); } } } @@ -1930,12 +1938,16 @@ void SolutionArray::readEntry(const AnyMap& root, const string& name, const stri const auto& components = path["components"].asVector(); bool back = false; for (const auto& name : components) { - if (hasComponent(name)) { + auto _name = name; + if (hasComponent(name, false)) { back = true; } else { - addExtra(name, back); + if (reverseAliasMap.count(name)) { + _name = reverseAliasMap.at(name); + } + addExtra(_name, back); } - setComponent(name, path[name]); + setComponent(_name, path[name]); exclude.insert(name); } } else { @@ -1944,10 +1956,14 @@ void SolutionArray::readEntry(const AnyMap& root, const string& name, const stri if (value.isVector()) { const vector& data = value.asVector(); if (data.size() == m_dataSize) { - if (!hasComponent(name)) { - addExtra(name); + auto _name = name; + if (!hasComponent(name, false)) { + if (reverseAliasMap.count(name)) { + _name = reverseAliasMap.at(name); + } + addExtra(_name); } - setComponent(name, value); + setComponent(_name, value); exclude.insert(name); } } From d72a51d71baa633f0a974bcbd664b14175dcb8f6 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 4 Oct 2025 19:09:22 -0500 Subject: [PATCH 06/22] [oneD] Implement component alias names --- include/cantera/base/SolutionArray.h | 3 ++ include/cantera/oneD/Domain1D.h | 7 +++- include/cantera/oneD/Flow1D.h | 4 +- interfaces/cython/cantera/_onedim.pxd | 1 + interfaces/cython/cantera/_onedim.pyx | 4 ++ interfaces/cython/cantera/onedim.py | 24 ++++++++++-- src/base/SolutionArray.cpp | 5 +++ src/oneD/Domain1D.cpp | 15 +++++++- src/oneD/Flow1D.cpp | 55 ++++++++++++++++++--------- 9 files changed, 92 insertions(+), 26 deletions(-) diff --git a/include/cantera/base/SolutionArray.h b/include/cantera/base/SolutionArray.h index 564a0b112ae..7db418b5dea 100644 --- a/include/cantera/base/SolutionArray.h +++ b/include/cantera/base/SolutionArray.h @@ -404,6 +404,9 @@ class SolutionArray vector m_active; //!< Vector of locations referencing active entries }; +//! Return mapping of component alias names to standardized component names. +const std::map& _componentAliasMap(); + } #endif diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index f7d53003d0c..ad0e7b1107b 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -197,7 +197,12 @@ class Domain1D } //! index of component with name `name`. - virtual size_t componentIndex(const string& name) const; + virtual size_t componentIndex(const string& name, bool checkAlias=true) const; + + /** + * Check whether the Domain contains a component. + */ + virtual bool hasComponent(const string& name, bool checkAlias=true) const; /** * Set the upper and lower bounds for a solution component, n. diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index ee416547bf2..37203ad51fb 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -175,7 +175,9 @@ class Flow1D : public Domain1D string componentName(size_t n) const override; - size_t componentIndex(const string& name) const override; + size_t componentIndex(const string& name, bool checkAlias=true) const override; + + bool hasComponent(const string& name, bool checkAlias=true) const override; //! Returns true if the specified component is an active part of the solver state virtual bool componentActive(size_t n) const; diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 3019e6f326e..5e3bcfb525e 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -24,6 +24,7 @@ cdef extern from "cantera/oneD/Domain1D.h": size_t nPoints() string componentName(size_t) except +translate_exception size_t componentIndex(string) except +translate_exception + cbool hasComponent(string) except +translate_exception void setBounds(size_t, double, double) double upperBound(size_t) double lowerBound(size_t) diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index bfb3ff05e13..11f3c77cff8 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -86,6 +86,10 @@ cdef class Domain1D: """Index of the component with name 'name'""" return self.domain.componentIndex(stringify(name)) + def _has_component(self, str name): + """Check whether `Domain1D` has component""" + return self.domain.hasComponent(stringify(name)) + def set_bounds(self, *, default=None, Y=None, **kwargs): """ Set the lower and upper bounds on the solution. diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index 0353cb70cfd..23c0de55592 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -289,7 +289,7 @@ def spread_rate(self): return self.profile(self.flame, 'spread_rate') @property - def L(self): + def radial_pressure_gradient(self): """ Array containing the radial pressure gradient (1/r)(dP/dr) [N/m⁴] at each point. Note: This value is named `Lambda` in the C++ code. @@ -297,7 +297,7 @@ def L(self): return self.profile(self.flame, "Lambda") @property - def E(self): + def electric_field(self): """ Array containing the electric field strength at each point. """ @@ -307,13 +307,31 @@ def E(self): return self.profile(self.flame, 'eField') @property - def Uo(self): + def oxidizer_velocity(self): """ Array containing the oxidizer velocity (right boundary velocity) [m/s] at each point. Note: This value is only defined when using two-point control. """ return self.profile(self.flame, 'Uo') + def __getattr__(self, name): + try: + return object.__getattribute__(self, name) + except AttributeError: + pass + + # Fallback to flame component lookup + flame = object.__getattribute__(self, 'flame') + + if (not flame._has_component(name) and + flame._has_component(name.replace("_", "-"))): + name = name.replace("_", "-") + if flame._has_component(name): + return self.profile(flame, name) + + raise AttributeError( + f"{flame.__class__.__name__!r} object has no attribute {name!r}") + @property def left_control_point_temperature(self): """ Get/Set the left control point temperature [K] """ diff --git a/src/base/SolutionArray.cpp b/src/base/SolutionArray.cpp index d9a249b28fe..c9c3a27354a 100644 --- a/src/base/SolutionArray.cpp +++ b/src/base/SolutionArray.cpp @@ -53,6 +53,11 @@ const std::map reverseAliasMap = { namespace Cantera { +const std::map& _componentAliasMap() +{ + return reverseAliasMap; +} + SolutionArray::SolutionArray(const shared_ptr& sol, int size, const AnyMap& meta) : m_sol(sol) diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index d3a77b8f114..8df69d20497 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -73,7 +73,7 @@ string Domain1D::componentName(size_t n) const } } -size_t Domain1D::componentIndex(const string& name) const +size_t Domain1D::componentIndex(const string& name, bool checkAlias) const { for (size_t n = 0; n < nComponents(); n++) { if (name == componentName(n)) { @@ -81,7 +81,18 @@ size_t Domain1D::componentIndex(const string& name) const } } throw CanteraError("Domain1D::componentIndex", - "no component named "+name); + "no component named '{}'", name); +} + +bool Domain1D::hasComponent(const string& name, bool checkAlias) const +{ + for (size_t n = 0; n < nComponents(); n++) { + if (name == componentName(n)) { + return true; + } + } + throw CanteraError("Domain1D::hasComponent", + "no component named '{}'", name); } void Domain1D::setTransientTolerances(double rtol, double atol, size_t n) diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index d2f10d95482..7007f181f12 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -16,6 +16,15 @@ using namespace std; namespace Cantera { +const std::map componentMap = { + {"velocity", c_offset_U}, // axial velocity [m/s] + {"spread_rate", c_offset_V}, // strain rate + {"T", c_offset_T}, // temperature [kelvin] + {"Lambda", c_offset_L}, // (1/r)dP/dr + {"eField", c_offset_E}, // electric field + {"Uo", c_offset_Uo}, // oxidizer axial velocity [m/s] +}; + Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) : Domain1D(nsp+c_offset_Y, points), m_nsp(nsp) @@ -833,33 +842,41 @@ string Flow1D::componentName(size_t n) const } } -size_t Flow1D::componentIndex(const string& name) const +size_t Flow1D::componentIndex(const string& name, bool checkAlias) const { - if (name=="velocity") { - return c_offset_U; - } else if (name=="spread_rate") { - return c_offset_V; - } else if (name=="T") { - return c_offset_T; - } else if (name=="Lambda") { - return c_offset_L; - } else if (name == "eField") { - return c_offset_E; - } else if (name == "Uo") { - return c_offset_Uo; - } else if (name=="lambda") { - warn_deprecated("Flow1D::componentIndex", - "Component 'lambda' is renamed to 'Lambda'."); - return c_offset_L; + if (componentMap.count(name)) { + return componentMap.at(name); } else { for (size_t n=c_offset_Y; n Date: Sun, 5 Oct 2025 08:25:08 -0500 Subject: [PATCH 07/22] [oneD] Implement set/getValues --- include/cantera/oneD/Domain1D.h | 26 +++++++++++++ include/cantera/oneD/Flow1D.h | 5 +++ include/cantera/oneD/OneDim.h | 2 +- include/cantera/oneD/Sim1D.h | 22 +++++++++++ interfaces/cython/cantera/_onedim.pxd | 2 + interfaces/cython/cantera/_onedim.pyx | 38 ++++++++++++++++++ interfaces/cython/cantera/onedim.py | 39 ++++++++++++------- .../src/sourcegen/headers/ctonedim_auto.yaml | 2 + src/oneD/Domain1D.cpp | 2 +- src/oneD/Flow1D.cpp | 34 ++++++++++++++++ src/oneD/OneDim.cpp | 2 +- src/oneD/Sim1D.cpp | 14 +++++++ 12 files changed, 172 insertions(+), 16 deletions(-) diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index ad0e7b1107b..bc52bcf9161 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -353,6 +353,32 @@ class Domain1D return x[index(n,j)]; } + /** + * Retrieve component values from the solution vector. + * @param soln local solution vector for this domain + * @param component name of component + * @param[out] values vector of values + * + * @since New in %Cantera 3.2. + */ + virtual void getValues(const double* soln, const string& component, + vector& values) const { + throw NotImplementedError("Domain1D::getValues", "Needs to be overloaded."); + } + + /** + * Specify component values in the solution vector. + * @param soln local solution vector for this domain + * @param component name of component + * @param[in] values vector of values + * + * @since New in %Cantera 3.2. + */ + virtual void setValues(double* soln, const string& component, + const vector& values) { + throw NotImplementedError("Domain1D::setValues", "Needs to be overloaded."); + } + //! Save the state of this domain as a SolutionArray. /*! * @param soln local solution vector for this domain diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index 37203ad51fb..16c7781ee27 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -184,6 +184,11 @@ class Flow1D : public Domain1D void show(const double* x) override; + void getValues(const double* soln, const string& component, + vector& values) const override; + void setValues(double* soln, const string& component, + const vector& values) override; + shared_ptr asArray(const double* soln) const override; void fromArray(SolutionArray& arr, double* soln) override; diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index 1a824c442c9..6669289c662 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -80,7 +80,7 @@ class OneDim : public SteadyStateSystem } //! Get the index of the domain named `name`. - size_t domainIndex(const string& name); + size_t domainIndex(const string& name) const; //! Check that the specified domain index is in range. //! Throws an exception if n is greater than nDomains()-1 diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index 035fb1d0a3c..5d3b9546551 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/include/cantera/oneD/Sim1D.h @@ -79,6 +79,28 @@ class Sim1D : public OneDim //! @param localPoint grid point within the domain double workValue(size_t dom, size_t comp, size_t localPoint) const; + /** + * Retrieve component values from the solution vector. + * @param dom name of domain + * @param component name of component + * @param[out] values vector of values + * + * @since New in %Cantera 3.2. + */ + void getValues(const string& dom, const string& component, + vector& values) const; + + /** + * Specify component values in the solution vector. + * @param dom name of domain + * @param component name of component + * @param[in] values vector of values + * + * @since New in %Cantera 3.2. + */ + void setValues(const string& dom, const string& component, + const vector& values); + /** * Specify a profile for one component of one domain. * @param dom domain number, beginning with 0 for the leftmost domain. diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 5e3bcfb525e..be20c26499f 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -111,6 +111,8 @@ cdef extern from "cantera/oneD/Sim1D.h": cdef cppclass CxxSim1D "Cantera::Sim1D": CxxSim1D(vector[shared_ptr[CxxDomain1D]]&) except +translate_exception void setValue(size_t, size_t, size_t, double) except +translate_exception + void getValues(string&, string&, vector[double]&) except +translate_exception + void setValues(string&, string&, vector[double]&) except +translate_exception void setProfile(size_t, size_t, vector[double]&, vector[double]&) except +translate_exception void setFlatProfile(size_t, size_t, double) except +translate_exception void show() except +translate_exception diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 11f3c77cff8..de07fa400b4 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -929,6 +929,44 @@ cdef class Sim1D: dom, comp = self._get_indices(domain, component) return self.sim.workValue(dom, comp, point) + def get_values(self, Domain1D domain, str component): + """ + Retrieve spatial profile of one component in one domain. + + :param domain: + Domain1D object + :param component: + component name + + >>> T = s.get_values(flow, 'T') + + .. versionadded:: 3.2 + """ + cdef vector[double] values + values.resize(domain.n_points) + self.sim.getValues(stringify(domain.name), stringify(component), values) + return np.asarray(values) + + def set_values(self, Domain1D domain, str component, values): + """ + Specify spatial profile of one component in one domain. + + :param domain: + Domain1D object + :param component: + component name + :param values: + array containing values + + >>> s.set_values(flow, 'T', T) + + .. versionadded:: 3.2 + """ + cdef vector[double] values_vec + for v in values: + values_vec.push_back(v) + self.sim.setValues(stringify(domain.name), stringify(component), values) + def profile(self, domain, component): """ Spatial profile of one component in one domain. diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index 23c0de55592..7340dda9770 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -271,14 +271,14 @@ def P(self, P): @property def T(self): """ Array containing the temperature [K] at each grid point. """ - return self.profile(self.flame, 'T') + return self.get_values(self.flame, "T") @property def velocity(self): """ Array containing the velocity [m/s] normal to the flame at each point. """ - return self.profile(self.flame, 'velocity') + return self.get_values(self.flame, "velocity") @property def spread_rate(self): @@ -286,33 +286,46 @@ def spread_rate(self): Array containing the tangential velocity gradient [1/s] (that is, radial velocity divided by radius) at each point. """ - return self.profile(self.flame, 'spread_rate') + return self.get_values(self.flame, "spread_rate") @property def radial_pressure_gradient(self): """ Array containing the radial pressure gradient (1/r)(dP/dr) [N/m⁴] at - each point. Note: This value is named `Lambda` in the C++ code. + each point. Note: This value is named ``Lambda`` in the C++ code. + + .. versionchanged:: 3.2 + + Renamed from ``L``, which remains accessible as an attribute. """ - return self.profile(self.flame, "Lambda") + return self.get_values(self.flame, "Lambda") @property def electric_field(self): """ Array containing the electric field strength at each point. + Note: This value is named ``eField`` in the C++ code and is only defined if + the transport model is ``ionized_gas``. + + .. versionchanged:: 3.2 + + Renamed from ``E``, which remains accessible as an attribute. """ - if self.flame.transport_model != "ionized-gas": - raise AttributeError( - "Electric field is only defined for transport model 'ionized_gas'.") - return self.profile(self.flame, 'eField') + return self.get_values(self.flame, "eField") @property def oxidizer_velocity(self): """ Array containing the oxidizer velocity (right boundary velocity) [m/s] at - each point. Note: This value is only defined when using two-point control. + each point. + Note: This value is named ``Uo`` in the C++ code and is only defined when using + two-point control. + + .. versionchanged:: 3.2 + + Renamed from ``Uo``, which remains accessible as an attribute. """ - return self.profile(self.flame, 'Uo') + return self.get_values(self.flame, "Uo") def __getattr__(self, name): try: @@ -321,13 +334,13 @@ def __getattr__(self, name): pass # Fallback to flame component lookup - flame = object.__getattribute__(self, 'flame') + flame = object.__getattribute__(self, "flame") if (not flame._has_component(name) and flame._has_component(name.replace("_", "-"))): name = name.replace("_", "-") if flame._has_component(name): - return self.profile(flame, name) + return self.get_values(flame, name) raise AttributeError( f"{flame.__class__.__name__!r} object has no attribute {name!r}") diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml index c86588e69e3..cbaf19e6bbe 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml @@ -10,6 +10,8 @@ recipes: - name: newSim1D uses: domain - name: setValue +- name: getValues +- name: setValues - name: setProfile - name: setFlatProfile - name: setInitialGuess # New in Cantera 3.2 diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index 8df69d20497..f4149e4fb9b 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -160,7 +160,7 @@ AnyMap Domain1D::getMeta() const shared_ptr Domain1D::toArray(bool normalize) const { if (!m_state) { throw CanteraError("Domain1D::toArray", - "Domain needs to be installed in a container before calling asArray."); + "Domain needs to be installed in a container before calling toArray."); } auto ret = asArray(m_state->data() + m_iloc); if (normalize) { diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 7007f181f12..d13e7484664 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -946,6 +946,40 @@ AnyMap Flow1D::getMeta() const return state; } +void Flow1D::getValues(const double* soln, const string& component, + vector& values) const +{ + if (values.size() != nPoints()) { + throw ArraySizeError("Flow1D::getValues", values.size(), nPoints()); + } + auto i = componentIndex(component); + if (!componentActive(i)) { + warn_user( + "Flow1D::getValues", "Component '{}' is not used by '{}'.", + component, domainType()); + } + for (size_t j = 0; j < nPoints(); j++) { + values[j] = soln[index(i,j)]; + } +} + +void Flow1D::setValues(double* soln, const string& component, + const vector& values) +{ + if (values.size() != nPoints()) { + throw ArraySizeError("Flow1D::setValues", values.size(), nPoints()); + } + auto i = componentIndex(component); + if (!componentActive(i)) { + throw CanteraError( + "Flow1D::setValues", "Component '{}' is not used by '{}'.", + component, domainType()); + } + for (size_t j = 0; j < nPoints(); j++) { + soln[index(i,j)] = values[j]; + } +} + shared_ptr Flow1D::asArray(const double* soln) const { auto arr = SolutionArray::create( diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 0afa0d7a017..2b8eec13440 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -26,7 +26,7 @@ OneDim::OneDim(vector>& domains) resize(); } -size_t OneDim::domainIndex(const string& name) +size_t OneDim::domainIndex(const string& name) const { for (size_t n = 0; n < m_dom.size(); n++) { if (domain(n).id() == name) { diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index 50528596061..c5f758a93c8 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -72,6 +72,20 @@ double Sim1D::workValue(size_t dom, size_t comp, size_t localPoint) const return m_xnew[iloc]; } +void Sim1D::getValues(const string& dom, const string& component, + vector& values) const +{ + size_t dIdx = domainIndex(dom); + domain(dIdx).getValues(m_state->data() + domain(dIdx).loc(), component, values); +} + +void Sim1D::setValues(const string& dom, const string& component, + const vector& values) +{ + size_t dIdx = domainIndex(dom); + domain(dIdx).setValues(m_state->data() + domain(dIdx).loc(), component, values); +} + void Sim1D::setProfile(size_t dom, size_t comp, const vector& pos, const vector& values) { From 0c8cbeead265a2e63677c079430e4378fd8f66a0 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 4 Oct 2025 10:06:35 -0500 Subject: [PATCH 08/22] [unittest] Test deprecated 'lambda' and other updates --- test/data/flame_lambda.h5 | Bin 0 -> 37112 bytes test/data/flame_lambda.yaml | 86 ++++++++++++++++++++++++++++++++++++ test/python/test_onedim.py | 39 ++++++++++++++-- 3 files changed, 122 insertions(+), 3 deletions(-) create mode 100644 test/data/flame_lambda.h5 create mode 100644 test/data/flame_lambda.yaml diff --git a/test/data/flame_lambda.h5 b/test/data/flame_lambda.h5 new file mode 100644 index 0000000000000000000000000000000000000000..c42e5f6bb6d9516da9593b695b82cfe969c1180c GIT binary patch literal 37112 zcmeHQ4|Ej88J|l`2uX|qL4r+Wc|sEUyd)+-Fv9Jn33?zVk)VQF$R@W*4(^|Omyl2_ zDApoXOrc1$qN!;uTCI6X`zZ1!!9Nd>Rw^yVXDz9sri!ALDn*pew=-Yvvhi{ug%rpP zkDHyDZ@-!O&Ft*%@B3zRudVQ!3u7;hB@kj_$OR-(Kg6GtJRECaIhe-sGJas<@k$>Q0{;j};`zh4exux++0&;JDvA0r=BJU4D0UQs86qB32-v3Qmq3T-@?PV7+QsJ? zg5mQxYy1w68`C&1t(b9&v5CZVo4<*C{-zIMidwk7aPxQ9cbK5W*#mf^fk#}h-urA5Vsuma<{DMd0tJm1FC z)rNM2P&a6^@w&$wSquj-o>+BEOeBhopeWlFKMWTcZ8GUAMIt?BqHF{ZcYXchg93T+ z#0IBP**INKmBVkY@VHzKJ`cj_8ww|LeTF3wNb ztKs+v{k@;_5x@lm^u;oMLv_>Fe!oZKkcHEy^95WVa0mf?@wz^{>Z{Nk-Y)7loIagb z;QD}T2*vXzrmzOqal#`U&l__klSh7-9L?ejg!;bU@gcSao5<8RJ0xZ2fd8 zKA*YLS%duOFiGw@e_|Y&c+$Z0vlev7v`$?mA8h*%&X5D*9m1Ox&C0fB%(Kp-Fx5C{ka1Ofwq z0P?5EimV*dd!_h(Tf*OyLrF6hhmf+IP!i`?{=o5AiRT25)f10&eZWnOmh0kAke5L| z1=IcK`I0z3CEX@iPoA&3Te(OG1Ox&Cfx(8rZ>N2pOupPC>-}!ktgQc^vh-;H!$cDF zwH}`ba6d`s)M1_DvI9FX%oLjQsm-cPjBwK&df$8o;y zx`EY$^#<&neY0nyeJh@@=ModSh@?_07Ptdd)tgnfQtkvRQ4%#o!Ltz9#o04a47nzM zR_Sz_H)r;Z5a|W@ycCF|zV$k5x9Vq*G8R!Cu-}-3yw13ohz2PFc%9AV{ioAV4^s?x z1e|x>pR&1dKO*)XdI9$mA;kOr-jz(@#mg&um$UNXhN|A>@V?S?j^9vly{hKxRT=Mg zz>j8#?l(8wA zxYTk@U7T5<_!em%k1NyVa69OB$0{|mT=A)ys=L;giLRJ-wKmh~D9`lP`>Q?foQx~8 zGkp%fYW6A>i>Rr+OE8n}NmVkdvpv}v^@@u>FQAcva z6W!qfJgcVI9SZD{3{b8538%`pvRZ(P6+9|pA>y70gyK<`cA5*4KYC~KrF1`?Rj&Z-d+_}?pPznl ztlq4Y`?>sDX)Ljv!G>lSc#f$`}7{yWqM+*iP^l%TJ7r5;W72d#DeUFQPh zGeY0j2U{=ARj``E+@a*GgKeM;w~sT_KJ%e8sxT@DRV>jzFk&wS8W9d0b~s z_qy76GplDauP>f3OBC!YhJDR0rOu3dn?W+vbYTWQ36mFe5@#JT( zUteCb!+5^W>Crc6hdck_hamhjXdNG_*AM);nOa|1FT=Ovg>7%Vd2s1l(~Rw?sL^Ws z(vLtp5J!ds{YHkv(+3<_JI7h^zGvL;PuE#}D)e{XH`*~UhP4*gU9dfXJp$8UqX4@G zreVzweqIIJ1~dlCao;$mVQVNn(_$KIE5Y*6hTt^Vnu6t_uHZD-mSCEpK4@!j8f;0y z{Q(4W0kmVriXxQOlukUa{onFP&vxvcfWEc}DVwM=f&u12rczubX4J zw&8}$ewZ-MGV{QOqN4rHQ#&5qRCf8@uDlH?^A?&9AC~vjZAyu{bB4U)k*yarZ_ScF zuDNmlL#YMw_imf;*yicSEPJnb)xSA=vwY~KrCSbsu+Q?^17ly7|IlQy*%}T#J)?PQ z;>OEIzWlR}Eu;TBds+P2W_gD5aNDHDJLPA8U9fz6JVY~9-lGq=ef zr0(7Qgtb$SYx1vqqbg3mCE|(a$Go)9a&lPO_@uG(r+ues{C%kcK2uFpO#0>{G{>T&%Y}V z{oOYyt-E67x~S(?jC^S5)SVAU%SYQ* zZOMCW&skFaKYHroLdkO3>x^}*u%*Spuq8~%7-^5UP6w!X$Ew(m~x&b(>si%SdQcC1}( zjeRZQ=boACtc{mO&bV#S53F)-#?yaW&fHXl_M*3kLcZJI6C!(p94!3Hi-$JB+R>7f=LVg0$tw;}`e;E3e8U4@D zuMF$MNx|MuF*kfXaO3xB8EehiJ$?g+hoqbfcrbtGXucbO?vq7_*bq|--rXVW<90Nh zmqxpqq)So%Y0gujoRPPhi*k$jL5K)A2MAO0S&aIO_61xWce$9{K>ZHfyed&xpU(E{ z?A!ltSL5xh^?m)Gt5n0}*egx|OU%o>cV?7h!A;tb7pNwhXYk|MTH1N65A56nN zB{&UyaIie`!6*-9z#n58j)T*{CBuI#M}8X9@T?lFANc9hVFyM1z@G=kw{QG3j&Hd8 zNI;HOLPP(bx7+Vw8t{(JQa#@{>R7@5>DTyOkVJp)4Z-Jx_SL=1;r*oZf!^grY5x!i zwciS^MF`@1t_%mfZ{hYJyuR^$5!Mn2IIoTSyjZFDo-5)YVu64_Kp-Fx5C{ka1Ofs9 zfq+0jARrJJSOi|&?^-oH$X~Gc+E3qi-g>V`lnMj{0s;YnfIvVXa3&ELI6GFyIbp}@ miATB)Kl&9S?~i_k#&rMfe4QMhl5WG@)6R#UGu_0_SmM7WwyW|0 literal 0 HcmV?d00001 diff --git a/test/data/flame_lambda.yaml b/test/data/flame_lambda.yaml new file mode 100644 index 00000000000..7ebd5328da5 --- /dev/null +++ b/test/data/flame_lambda.yaml @@ -0,0 +1,86 @@ +solution: + description: | + Test import for legacy component names 'lambda' and 'spread_rate' + (replaced by 'Lambda' and 'spreadRate' in Cantera 3.2) + generator: Cantera SolutionArray + cantera-version: 3.2.0a4 + git-commit: d7f7a0a6a + date: Tue Oct 7 21:22:10 2025 + fuel_inlet: + type: inlet + size: 1 + transport-model: mixture-averaged + points: 1 + mass-flux: 0.5 + temperature: 300.0 + pressure: 1.01325e+05 + mass-fractions: + H2: 1.0 + flame: + type: axisymmetric-flow + size: 6 + basis: mass + transport-model: mixture-averaged + points: 6 + tolerances: + transient-abstol: 1.0e-11 + steady-abstol: 1.0e-09 + transient-reltol: 1.0e-04 + steady-reltol: 1.0e-04 + phase: + name: ohmech + source: /opt/homebrew/Caskroom/miniforge/base/envs/cantera-dev/lib/python3.12/site-packages/cantera/data/h2o2.yaml + radiation-enabled: false + energy-enabled: true + Soret-enabled: false + flux-gradient-basis: 1 + refine-criteria: + ratio: 10.0 + slope: 0.8 + curve: 0.8 + prune: -1.0e-03 + grid-min: 1.0e-10 + max-points: 1000 + components: [grid, velocity, spread_rate, lambda, T, D, H2, H, O, O2, OH, H2O, HO2, + H2O2, AR, N2] + grid: [0.0, 3.6e-03, 7.2e-03, 0.0108, 0.0144, 0.018] + velocity: [6.102932415387452, 5.212628821879346, -2.947066925624583, + -2.680752244268045, -1.763938145228785, -2.308011820997561] + spread_rate: [-4.833540765078792e-18, 978.4682832024049, 1362.608989936244, + 438.7100305672880, 163.2616573240902, 0.0] + lambda: [-1.355348889404357e+05, -1.355348889404357e+05, -1.355348889404357e+05, + -1.355348889404357e+05, -1.355348889404357e+05, -1.355348889404357e+05] + T: [300.0, 468.5529213409615, 3094.093759960199, 1056.968577408735, + 305.5510768535476, 300.0] + D: [0.08192782845493371, 0.05724023330548490, 0.04940148245658487, + 0.3522103611612766, 1.275683705975624, 1.299820032422256] + H2: [0.9995132911440364, 0.9038302495109068, 0.03823126844637609, + 3.019659739704242e-04, 6.437965548973601e-06, 6.235770419963188e-08] + H: [9.832982544077039e-06, 7.783198338165095e-04, 8.591098766876597e-03, + 9.982933233000733e-07, 1.178976989033755e-12, 1.789939048465165e-14] + O: [1.606510173235227e-06, 2.5538425392075e-04, 0.02931115777838727, + 1.299441663252144e-04, 3.415964904326214e-07, 1.358647080040236e-09] + O2: [2.741298317654846e-06, 5.844355679688863e-04, 0.05150688826729902, + 0.9441585603177651, 0.9995862996501134, 0.9999988134546678] + OH: [1.451168485312993e-07, 2.314791522403334e-05, 0.1144112864470808, + 1.748691889619610e-04, 4.237601773617709e-07, 1.651940626064034e-09] + H2O: [4.723617132018648e-04, 0.09452393390090491, 0.7578920308453377, + 0.05497698284946553, 4.052338848858928e-04, 1.118028628452684e-06] + HO2: [1.889998199384880e-08, 4.030613957970947e-06, 5.227209720209965e-05, + 2.055924382742969e-04, 7.574965166245186e-07, 1.893602903100419e-09] + H2O2: [2.334896930790474e-09, 4.984029919993926e-07, 3.997351820551009e-06, + 5.108677191798846e-05, 5.056450883263629e-07, 1.254791074753632e-09] + AR: [1.375664417203775e-16, 1.338359588790799e-16, 6.533245011697561e-18, + 1.734586902703745e-19, -4.913882049989299e-22, -1.177879523773831e-24] + N2: [2.980371007988044e-27, -4.286122226344077e-25, -9.912199957010441e-26, + 1.080732196310891e-26, 5.715707190638608e-29, 1.441959686430471e-31] + oxidizer_inlet: + type: inlet + size: 1 + transport-model: mixture-averaged + points: 1 + mass-flux: 3.0 + temperature: 300.0 + pressure: 1.01325e+05 + mass-fractions: + O2: 1.0 diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index 3bb891423dd..b5525744e4d 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -1036,6 +1036,14 @@ def solve_mix(self, ratio=3.0, slope=0.1, curve=0.12, prune=0.0): assert self.sim.energy_enabled assert self.sim.transport_model == 'mixture-averaged' + arr = self.sim.to_array() + radial_lambda = np.unique(arr.Lambda) + assert len(radial_lambda) == 1 + assert getattr(arr, "lambda")[0] == radial_lambda[0] + assert arr.L[0] == radial_lambda[0] + assert arr.radial_pressure_gradient[0] == radial_lambda[0] + assert getattr(arr, "radial-pressure-gradient")[0] == radial_lambda[0] + @pytest.mark.slow_test def test_mixture_averaged(self, request, test_data_path): referenceFile = "DiffusionFlameTest-h2-mix.csv" @@ -1263,6 +1271,31 @@ def test_bad_boundary_conditions(self): sim.fuel_inlet.X = "H2: 1.0" sim.set_initial_guess() + def run_restore_lambda(self, fname): + gas = ct.Solution("h2o2.yaml") + sim = ct.CounterflowDiffusionFlame(gas) + sim.restore(fname) + assert np.isclose(np.unique(sim.L)[0], -90257.82873678) + + arr = ct.SolutionArray(gas) + arr.restore(fname, "solution/flame") + + radial_lambda = np.unique(arr.Lambda) + assert len(radial_lambda) == 1 + assert radial_lambda[0] == pytest.approx(-90257.82873678) + assert getattr(arr, "lambda")[0] == radial_lambda[0] + assert arr.L[0] == radial_lambda[0] + assert arr.radial_pressure_gradient[0] == radial_lambda[0] + assert getattr(arr, "radial-pressure-gradient")[0] == radial_lambda[0] + + def test_restore_lambda_yaml(self): + self.run_restore_lambda("flame_lambda.yaml") + + @pytest.mark.skipif("native" not in ct.hdf_support(), + reason="Cantera compiled without HDF support") + def test_restore_lambda_hdf(self): + self.run_restore_lambda("flame_lambda.h5") + def test_two_point_control(self): """ Computes a solution using a standard counterflow diffusion flame as a @@ -1331,7 +1364,7 @@ def test_two_point_control(self): sim.right_control_point_temperature += temperature_decrement sim.left_control_point_temperature += temperature_decrement sim.solve(loglevel=0, refine_grid=False) - assert sim.Uo == approx(sim.velocity[-1]) + assert sim.oxidizer_velocity == approx(sim.velocity[-1]) temperature_4 = sim.T # Check the difference between the un-perturbed two-point solution and the @@ -1493,7 +1526,7 @@ def _run_case(phi, T, width, P): sim.solve(loglevel=0, auto=True) assert all(sim.T >= T - 1e-3) assert all(sim.spread_rate >= -1e-9) - assert np.allclose(sim.L, sim.L[0]) + assert np.allclose(sim.radial_pressure_gradient, sim.L[0]) return sim return _run_case @@ -2009,7 +2042,7 @@ def test_ion_profile(self): # stage one assert self.sim.electric_field_enabled is False self.sim.solve(loglevel=0, auto=True) - assert max(np.abs(self.sim.E)) == approx(0.0, abs=1e-3) + assert max(np.abs(self.sim.electric_field)) == approx(0.0, abs=1e-3) #stage two self.sim.electric_field_enabled = True From 75dee32bb3a247a0d99e8691ba0821166b443389 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 5 Oct 2025 11:02:54 -0500 Subject: [PATCH 09/22] [onedim] Address remaining papercuts --- include/cantera/oneD/Domain1D.h | 8 +- interfaces/cython/cantera/_onedim.pxd | 4 +- interfaces/cython/cantera/_onedim.pyx | 6 ++ interfaces/cython/cantera/onedim.py | 26 ++++++- .../src/sourcegen/headers/ctdomain_auto.yaml | 1 + .../python/onedim/diffusion_flame_batch.py | 5 +- .../onedim/diffusion_flame_extinction.py | 3 +- src/base/SolutionArray.cpp | 2 +- src/oneD/Flow1D.cpp | 11 ++- test/python/test_onedim.py | 76 ++++++++++++++----- 10 files changed, 107 insertions(+), 35 deletions(-) diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index bc52bcf9161..8191f1dbd59 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -196,11 +196,17 @@ class Domain1D m_name[n] = name; } - //! index of component with name `name`. + /** + * Index of component with name `name`. + * @param name name of component + * @param checkAlias if `True` (default), check alias mapping + */ virtual size_t componentIndex(const string& name, bool checkAlias=true) const; /** * Check whether the Domain contains a component. + * @param name name of component + * @param checkAlias if `True` (default), check alias mapping */ virtual bool hasComponent(const string& name, bool checkAlias=true) const; diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index be20c26499f..57e6db0bd87 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -23,8 +23,8 @@ cdef extern from "cantera/oneD/Domain1D.h": size_t nComponents() size_t nPoints() string componentName(size_t) except +translate_exception - size_t componentIndex(string) except +translate_exception - cbool hasComponent(string) except +translate_exception + size_t componentIndex(string&) except +translate_exception + cbool hasComponent(string&) except +translate_exception void setBounds(size_t, double, double) double upperBound(size_t) double lowerBound(size_t) diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index de07fa400b4..f2484d662d4 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -977,7 +977,13 @@ cdef class Sim1D: component name or index >>> T = s.profile(flow, 'T') + + .. deprecated:: 3.2 + + To be removed after Cantera 3.2. Replaceable by `get_values`. """ + warnings.warn("To be removed after Cantera 3.2. Replaceable by 'get_values'.", + DeprecationWarning) idom, kcomp = self._get_indices(domain, component) dom = self.domains[idom] cdef int j diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index 7340dda9770..793ee841b63 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -296,7 +296,8 @@ def radial_pressure_gradient(self): .. versionchanged:: 3.2 - Renamed from ``L``, which remains accessible as an attribute. + Renamed from ``L``, which remains accessible as an attribute; note: may be + shadowed by a species name ``L``. """ return self.get_values(self.flame, "Lambda") @@ -307,10 +308,26 @@ def electric_field(self): Note: This value is named ``eField`` in the C++ code and is only defined if the transport model is ``ionized_gas``. - .. versionchanged:: 3.2 + .. versionadded:: 3.2 + + Renamed from ``E``, which remains accessible as an attribute; note: may be + shadowed by a species name ``E``. + """ + return self.get_values(self.flame, "eField") + + @property + def E(self): + """ + Array containing the electric field strength at each point. + Note: This value is named ``eField`` in the C++ code and is only defined if + the transport model is ``ionized_gas``. + + .. deprecated:: 3.2 - Renamed from ``E``, which remains accessible as an attribute. + To be removed after Cantera 3.2. Renamed to `electric_field`. """ + warnings.warn("To be removed after Cantera 3.2. Renamed to 'electric_field'.", + DeprecationWarning) return self.get_values(self.flame, "eField") @property @@ -323,7 +340,8 @@ def oxidizer_velocity(self): .. versionchanged:: 3.2 - Renamed from ``Uo``, which remains accessible as an attribute. + Renamed from ``Uo``, which remains accessible as an attribute; note: may be + shadowed by a species name ``Uo``. """ return self.get_values(self.flame, "Uo") diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml index 0a38fda0728..acb61bcc3af 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml @@ -19,6 +19,7 @@ recipes: - name: nPoints - name: componentName - name: componentIndex + wraps: componentIndex(const char*) - name: setBounds - name: lowerBound - name: upperBound diff --git a/samples/python/onedim/diffusion_flame_batch.py b/samples/python/onedim/diffusion_flame_batch.py index 224ce3d098a..69a8210871c 100644 --- a/samples/python/onedim/diffusion_flame_batch.py +++ b/samples/python/onedim/diffusion_flame_batch.py @@ -142,7 +142,7 @@ def names(test): f.spread_rate * rel_pressure_increase ** exp_V_p) # Update pressure curvature f.set_profile("Lambda", normalized_grid, - f.L * rel_pressure_increase ** exp_lam_p) + f.radial_pressure_gradient * rel_pressure_increase ** exp_lam_p) try: # Try solving the flame @@ -196,7 +196,8 @@ def names(test): f.set_profile('spread_rate', normalized_grid, f.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature - f.set_profile("Lambda", normalized_grid, f.L * strain_factor ** exp_lam_a) + f.set_profile("Lambda", normalized_grid, + f.radial_pressure_gradient * strain_factor ** exp_lam_a) try: # Try solving the flame f.solve(loglevel=0) diff --git a/samples/python/onedim/diffusion_flame_extinction.py b/samples/python/onedim/diffusion_flame_extinction.py index d17b2fc204c..074f4b5e51b 100644 --- a/samples/python/onedim/diffusion_flame_extinction.py +++ b/samples/python/onedim/diffusion_flame_extinction.py @@ -134,7 +134,8 @@ def names(test): f.set_profile('spread_rate', normalized_grid, f.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature - f.set_profile("Lambda", normalized_grid, f.L * strain_factor ** exp_lam_a) + f.set_profile("Lambda", normalized_grid, + f.radial_pressure_gradient * strain_factor ** exp_lam_a) try: f.solve(loglevel=0) except ct.CanteraError as e: diff --git a/src/base/SolutionArray.cpp b/src/base/SolutionArray.cpp index c9c3a27354a..23c50e94cec 100644 --- a/src/base/SolutionArray.cpp +++ b/src/base/SolutionArray.cpp @@ -694,7 +694,7 @@ bool SolutionArray::hasComponent(const string& name, bool checkAlias) const AnyValue SolutionArray::getComponent(const string& name) const { string _name = name; - if (reverseAliasMap.count(name)) { + if (!hasComponent(name, false) && reverseAliasMap.count(name)) { _name = reverseAliasMap.at(name); } else if (!hasComponent(name)) { throw CanteraError("SolutionArray::getComponent", diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index d13e7484664..1762b09a9b9 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -846,11 +846,10 @@ size_t Flow1D::componentIndex(const string& name, bool checkAlias) const { if (componentMap.count(name)) { return componentMap.at(name); - } else { - for (size_t n=c_offset_Y; n= T - 1e-3) assert all(sim.spread_rate >= -1e-9) - assert np.allclose(sim.radial_pressure_gradient, sim.L[0]) + assert np.allclose(sim.radial_pressure_gradient, sim.radial_pressure_gradient[0]) return sim return _run_case @@ -1643,7 +1678,7 @@ def _run_case(phi, T, width, P): sim.solve(loglevel=0, auto=True) assert all(sim.T >= T - 1e-3) assert all(sim.spread_rate >= -1e-9) - assert np.allclose(sim.L, sim.L[0]) + assert np.allclose(sim.radial_pressure_gradient, sim.radial_pressure_gradient[0]) return sim return _run_case @@ -1690,7 +1725,7 @@ def _run_case(phi, T, width, P): sim.burner.mdot = gas.density * 0.15 sim.solve(loglevel=0, auto=True) assert sim.T[1] > T - assert np.allclose(sim.L, 0) + assert np.allclose(sim.radial_pressure_gradient, 0) return _run_case @pytest.mark.parametrize( @@ -1783,7 +1818,7 @@ def run_stagnation(self, xh2, mdot, width): sim.solve(loglevel=0, auto=True) assert sim.T.max() > tburner + tsurf - assert np.allclose(sim.L, sim.L[0]) + assert np.allclose(sim.radial_pressure_gradient, sim.radial_pressure_gradient[0]) self.sim = sim def test_stagnation_case1(self): @@ -1970,7 +2005,7 @@ def solve(self, phi, T, width, P): sim.reactants.mdot = gas.density * axial_velocity sim.solve(loglevel=0, auto=True) assert sim.T[-1] > T + 100 - assert np.allclose(sim.L, sim.L[0]) + assert np.allclose(sim.radial_pressure_gradient, sim.radial_pressure_gradient[0]) return sim def test_restart(self): @@ -2049,7 +2084,12 @@ def test_ion_profile(self): self.sim.solve(loglevel=0) # Regression test - assert max(self.sim.E) == approx(149.63179056676853, rel=1e-3) + assert max(self.sim.electric_field) == approx(149.63179056676853, rel=1e-3) + with pytest.warns(DeprecationWarning, match="Renamed to 'electric_field'"): + assert np.allclose(self.sim.E, self.sim.electric_field) + zeros = self.sim.T + with pytest.raises(ct.CanteraError, match="is not used by 'free-ion-flow'"): + self.sim.set_values(self.sim.flame, "oxidizer-velocity", zeros) class TestIonBurnerFlame: @@ -2071,5 +2111,5 @@ def test_ion_profile(self): self.sim.solve(loglevel=0, auto=True) # Regression test - assert max(self.sim.E) == approx(591.76, rel=1e-2) + assert max(self.sim.electric_field) == approx(591.76, rel=1e-2) assert max(self.sim.X[self.gas.species_index('E')]) == approx(8.024e-9, rel=1e-2) From 062542125999e658d38bf54f34099caa3fb9578d Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 5 Oct 2025 13:51:07 -0500 Subject: [PATCH 10/22] [oneD] Tie set/getValues to Domain1D --- include/cantera/oneD/Domain1D.h | 26 ++++++++++++++++++-------- include/cantera/oneD/Flow1D.h | 8 ++++---- include/cantera/oneD/OneDim.h | 26 ++++++++++++++++++++++++++ include/cantera/oneD/Sim1D.h | 4 ++-- src/oneD/Domain1D.cpp | 16 ++++++++++++++++ src/oneD/Flow1D.cpp | 14 +++++++------- src/oneD/Sim1D.cpp | 4 ++-- 7 files changed, 75 insertions(+), 23 deletions(-) diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index 8191f1dbd59..17e99dab014 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -361,27 +361,37 @@ class Domain1D /** * Retrieve component values from the solution vector. - * @param soln local solution vector for this domain * @param component name of component * @param[out] values vector of values * * @since New in %Cantera 3.2. */ - virtual void getValues(const double* soln, const string& component, - vector& values) const { - throw NotImplementedError("Domain1D::getValues", "Needs to be overloaded."); - } + void getValues(const string& component, vector& values) const; /** * Specify component values in the solution vector. - * @param soln local solution vector for this domain * @param component name of component * @param[in] values vector of values * * @since New in %Cantera 3.2. */ - virtual void setValues(double* soln, const string& component, - const vector& values) { + void setValues(const string& component, const vector& values); + + /** + * Retrieve component values from the solution vector. + * @since New in %Cantera 3.2. + */ + virtual void _getValues(const double* soln, const string& component, + vector& values) const { + throw NotImplementedError("Domain1D::getValues", "Needs to be overloaded."); + } + + /** + * Specify component values in the solution vector. + * @since New in %Cantera 3.2. + */ + virtual void _setValues(double* soln, const string& component, + const vector& values) { throw NotImplementedError("Domain1D::setValues", "Needs to be overloaded."); } diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index 16c7781ee27..86eebdaf4e5 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -184,10 +184,10 @@ class Flow1D : public Domain1D void show(const double* x) override; - void getValues(const double* soln, const string& component, - vector& values) const override; - void setValues(double* soln, const string& component, - const vector& values) override; + void _getValues(const double* soln, const string& component, + vector& values) const override; + void _setValues(double* soln, const string& component, + const vector& values) override; shared_ptr asArray(const double* soln) const override; void fromArray(SolutionArray& arr, double* soln) override; diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index 6669289c662..a1880ec32a1 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -141,6 +141,32 @@ class OneDim : public SteadyStateSystem double upperBound(size_t i) const override; double lowerBound(size_t i) const override; + /** + * Retrieve component values from the solution vector. + * @param dom name of domain + * @param component name of component + * @param[out] values vector of values + * + * @since New in %Cantera 3.2. + */ + virtual void getValues(const string& dom, const string& component, + vector& values) const { + throw NotImplementedError("OneDim::getValues", "Needs to be overloaded."); + } + + /** + * Specify component values in the solution vector. + * @param dom name of domain + * @param component name of component + * @param[in] values vector of values + * + * @since New in %Cantera 3.2. + */ + virtual void setValues(const string& dom, const string& component, + const vector& values) { + throw NotImplementedError("OneDim::setValues", "Needs to be overloaded."); + } + /** * Initialize all domains. On the first call, this methods calls the init * method of each domain, proceeding from left to right. Subsequent calls diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index 5d3b9546551..4c225970700 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/include/cantera/oneD/Sim1D.h @@ -88,7 +88,7 @@ class Sim1D : public OneDim * @since New in %Cantera 3.2. */ void getValues(const string& dom, const string& component, - vector& values) const; + vector& values) const override; /** * Specify component values in the solution vector. @@ -99,7 +99,7 @@ class Sim1D : public OneDim * @since New in %Cantera 3.2. */ void setValues(const string& dom, const string& component, - const vector& values); + const vector& values) override; /** * Specify a profile for one component of one domain. diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index f4149e4fb9b..1d341483192 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -95,6 +95,22 @@ bool Domain1D::hasComponent(const string& name, bool checkAlias) const "no component named '{}'", name); } +void Domain1D::getValues(const string& component, vector& values) const +{ + if (!m_container) { + throw CanteraError("Domain1D::getValues", "Domain is not connected."); + } + m_container->getValues(id(), component, values); +} + +void Domain1D::setValues(const string& component, const vector& values) +{ + if (!m_container) { + throw CanteraError("Domain1D::setValues", "Domain is not connected."); + } + m_container->setValues(id(), component, values); +} + void Domain1D::setTransientTolerances(double rtol, double atol, size_t n) { if (n == npos) { diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 1762b09a9b9..8bcfb49183b 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -945,16 +945,16 @@ AnyMap Flow1D::getMeta() const return state; } -void Flow1D::getValues(const double* soln, const string& component, - vector& values) const +void Flow1D::_getValues(const double* soln, const string& component, + vector& values) const { if (values.size() != nPoints()) { - throw ArraySizeError("Flow1D::getValues", values.size(), nPoints()); + throw ArraySizeError("Flow1D::_getValues", values.size(), nPoints()); } auto i = componentIndex(component); if (!componentActive(i)) { warn_user( - "Flow1D::getValues", "Component '{}' is not used by '{}'.", + "Flow1D::_getValues", "Component '{}' is not used by '{}'.", component, domainType()); } for (size_t j = 0; j < nPoints(); j++) { @@ -962,16 +962,16 @@ void Flow1D::getValues(const double* soln, const string& component, } } -void Flow1D::setValues(double* soln, const string& component, +void Flow1D::_setValues(double* soln, const string& component, const vector& values) { if (values.size() != nPoints()) { - throw ArraySizeError("Flow1D::setValues", values.size(), nPoints()); + throw ArraySizeError("Flow1D::_setValues", values.size(), nPoints()); } auto i = componentIndex(component); if (!componentActive(i)) { throw CanteraError( - "Flow1D::setValues", "Component '{}' is not used by '{}'.", + "Flow1D::_setValues", "Component '{}' is not used by '{}'.", component, domainType()); } for (size_t j = 0; j < nPoints(); j++) { diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index c5f758a93c8..a3251147ed5 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -76,14 +76,14 @@ void Sim1D::getValues(const string& dom, const string& component, vector& values) const { size_t dIdx = domainIndex(dom); - domain(dIdx).getValues(m_state->data() + domain(dIdx).loc(), component, values); + domain(dIdx)._getValues(m_state->data() + domain(dIdx).loc(), component, values); } void Sim1D::setValues(const string& dom, const string& component, const vector& values) { size_t dIdx = domainIndex(dom); - domain(dIdx).setValues(m_state->data() + domain(dIdx).loc(), component, values); + domain(dIdx)._setValues(m_state->data() + domain(dIdx).loc(), component, values); } void Sim1D::setProfile(size_t dom, size_t comp, From 925338c336649ba209c254454c29bff3d309ad7c Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 5 Oct 2025 13:51:32 -0500 Subject: [PATCH 11/22] [oneD] Update APIs for set/getValues --- interfaces/cython/cantera/_onedim.pxd | 4 +- interfaces/cython/cantera/_onedim.pyx | 163 +++++++++++++----- interfaces/cython/cantera/onedim.py | 79 +++------ .../src/sourcegen/headers/ctdomain_auto.yaml | 2 + 4 files changed, 150 insertions(+), 98 deletions(-) diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 57e6db0bd87..8374031cd41 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -24,6 +24,8 @@ cdef extern from "cantera/oneD/Domain1D.h": size_t nPoints() string componentName(size_t) except +translate_exception size_t componentIndex(string&) except +translate_exception + void getValues(string&, vector[double]&) except +translate_exception + void setValues(string&, vector[double]&) except +translate_exception cbool hasComponent(string&) except +translate_exception void setBounds(size_t, double, double) double upperBound(size_t) @@ -111,8 +113,6 @@ cdef extern from "cantera/oneD/Sim1D.h": cdef cppclass CxxSim1D "Cantera::Sim1D": CxxSim1D(vector[shared_ptr[CxxDomain1D]]&) except +translate_exception void setValue(size_t, size_t, size_t, double) except +translate_exception - void getValues(string&, string&, vector[double]&) except +translate_exception - void setValues(string&, string&, vector[double]&) except +translate_exception void setProfile(size_t, size_t, vector[double]&, vector[double]&) except +translate_exception void setFlatProfile(size_t, size_t, double) except +translate_exception void show() except +translate_exception diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index f2484d662d4..cbadbeaa716 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -90,6 +90,40 @@ cdef class Domain1D: """Check whether `Domain1D` has component""" return self.domain.hasComponent(stringify(name)) + def get_values(self, str component): + """ + Retrieve spatial profile of a component. + + :param component: + component name + + >>> T = d.get_values('T') + + .. versionadded:: 3.2 + """ + cdef vector[double] values + values.resize(self.n_points) + self.domain.getValues(stringify(component), values) + return np.asarray(values) + + def set_values(self, str component, values): + """ + Specify spatial profile of a component. + + :param component: + component name + :param values: + array containing values + + >>> d.set_values('T', T) + + .. versionadded:: 3.2 + """ + cdef vector[double] values_vec + for v in values: + values_vec.push_back(v) + self.domain.setValues(stringify(component), values) + def set_bounds(self, *, default=None, Y=None, **kwargs): """ Set the lower and upper bounds on the solution. @@ -460,12 +494,89 @@ cdef class FlowBase(Domain1D): self.P = self.gas.P self.flow.solveEnergyEqn() - property P: - """ Pressure [Pa] """ - def __get__(self): - return self.flow.pressure() - def __set__(self, P): - self.flow.setPressure(P) + def __getattr__(self, name): + component_name = name + if (not self._has_component(name) and + self._has_component(name.replace("_", "-"))): + component_name = name.replace("_", "-") + + if self._has_component(component_name): + return self.get_values(component_name) + + raise AttributeError( + f"{self.__class__.__name__!r} object has no attribute {name!r}") + + @property + def P(self): + """Pressure [Pa]""" + return self.flow.pressure() + + @P.setter + def P(self, P): + self.flow.setPressure(P) + + @property + def T(self): + """ + Array containing the temperature [K] at each grid point. + + .. versionadded:: 3.2 + """ + return self.get_values("T") + + @property + def velocity(self): + """ + Array containing the velocity [m/s] normal to the flame at each point. + + .. versionadded:: 3.2 + """ + return self.get_values("velocity") + + @property + def spread_rate(self): + """ + Array containing the tangential velocity gradient [1/s] (that is, radial + velocity divided by radius) at each point. Note: This value is only + defined for axisymmetric flows. + + .. versionadded:: 3.2 + """ + return self.get_values("spread_rate") + + @property + def radial_pressure_gradient(self): + """ + Array containing the radial pressure gradient (1/r)(dP/dr) [N/m⁴] at + each point. Note: This value is named ``Lambda`` in the C++ code and is only + defined for axisymmetric flows. + + .. versionadded:: 3.2 + """ + return self.get_values("Lambda") + + @property + def electric_field(self): + """ + Array containing the electric field strength at each point. + Note: This value is named ``eField`` in the C++ code and is only defined if + the transport model is ``ionized_gas``. + + .. versionadded:: 3.2 + """ + return self.get_values("eField") + + @property + def oxidizer_velocity(self): + """ + Array containing the oxidizer velocity (right boundary velocity) [m/s] at + each point. + Note: This value is named ``Uo`` in the C++ code and is only defined when using + two-point control. + + .. versionchanged:: 3.2 + """ + return self.get_values("Uo") property transport_model: """ @@ -929,44 +1040,6 @@ cdef class Sim1D: dom, comp = self._get_indices(domain, component) return self.sim.workValue(dom, comp, point) - def get_values(self, Domain1D domain, str component): - """ - Retrieve spatial profile of one component in one domain. - - :param domain: - Domain1D object - :param component: - component name - - >>> T = s.get_values(flow, 'T') - - .. versionadded:: 3.2 - """ - cdef vector[double] values - values.resize(domain.n_points) - self.sim.getValues(stringify(domain.name), stringify(component), values) - return np.asarray(values) - - def set_values(self, Domain1D domain, str component, values): - """ - Specify spatial profile of one component in one domain. - - :param domain: - Domain1D object - :param component: - component name - :param values: - array containing values - - >>> s.set_values(flow, 'T', T) - - .. versionadded:: 3.2 - """ - cdef vector[double] values_vec - for v in values: - values_vec.push_back(v) - self.sim.setValues(stringify(domain.name), stringify(component), values) - def profile(self, domain, component): """ Spatial profile of one component in one domain. @@ -980,7 +1053,7 @@ cdef class Sim1D: .. deprecated:: 3.2 - To be removed after Cantera 3.2. Replaceable by `get_values`. + To be removed after Cantera 3.2. Replaceable by `Domain1D.get_values`. """ warnings.warn("To be removed after Cantera 3.2. Replaceable by 'get_values'.", DeprecationWarning) diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index 793ee841b63..3563ca3b23f 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -271,14 +271,14 @@ def P(self, P): @property def T(self): """ Array containing the temperature [K] at each grid point. """ - return self.get_values(self.flame, "T") + return self.flame.get_values("T") @property def velocity(self): """ Array containing the velocity [m/s] normal to the flame at each point. """ - return self.get_values(self.flame, "velocity") + return self.flame.get_values("velocity") @property def spread_rate(self): @@ -286,34 +286,24 @@ def spread_rate(self): Array containing the tangential velocity gradient [1/s] (that is, radial velocity divided by radius) at each point. """ - return self.get_values(self.flame, "spread_rate") + return self.flame.get_values("spread_rate") @property - def radial_pressure_gradient(self): + def L(self): """ Array containing the radial pressure gradient (1/r)(dP/dr) [N/m⁴] at each point. Note: This value is named ``Lambda`` in the C++ code. - .. versionchanged:: 3.2 - - Renamed from ``L``, which remains accessible as an attribute; note: may be - shadowed by a species name ``L``. - """ - return self.get_values(self.flame, "Lambda") - - @property - def electric_field(self): - """ - Array containing the electric field strength at each point. - Note: This value is named ``eField`` in the C++ code and is only defined if - the transport model is ``ionized_gas``. - - .. versionadded:: 3.2 + .. deprecated:: 3.2 - Renamed from ``E``, which remains accessible as an attribute; note: may be - shadowed by a species name ``E``. + To be removed after Cantera 3.2. Replaceable by + `Domain1D.radial_pressure_gradient`. + Note: replacement may be shadowed by a species name ``L``. """ - return self.get_values(self.flame, "eField") + warnings.warn("To be removed after Cantera 3.2. " + "Replaceable by 'Domain1D.radial_pressure_gradient'.", + DeprecationWarning) + return self.flame.get_values("Lambda") @property def E(self): @@ -324,44 +314,31 @@ def E(self): .. deprecated:: 3.2 - To be removed after Cantera 3.2. Renamed to `electric_field`. + To be removed after Cantera 3.2. Replaceable by `Domain1D.electric_field`. + Note: replacement may be shadowed by a species name ``E``. """ - warnings.warn("To be removed after Cantera 3.2. Renamed to 'electric_field'.", - DeprecationWarning) - return self.get_values(self.flame, "eField") + warnings.warn("To be removed after Cantera 3.2. " + "Replaceable by 'Domain1D.electric_field'.", DeprecationWarning) + return self.flame.get_values("eField") @property - def oxidizer_velocity(self): + def Uo(self): """ Array containing the oxidizer velocity (right boundary velocity) [m/s] at each point. Note: This value is named ``Uo`` in the C++ code and is only defined when using two-point control. - .. versionchanged:: 3.2 + .. deprecated:: 3.2 - Renamed from ``Uo``, which remains accessible as an attribute; note: may be - shadowed by a species name ``Uo``. + To be removed after Cantera 3.2. Replaceable by + `Domain1D.oxidizer_velocity`. + Note: replacement may be shadowed by a species name ``Uo``. """ - return self.get_values(self.flame, "Uo") - - def __getattr__(self, name): - try: - return object.__getattribute__(self, name) - except AttributeError: - pass - - # Fallback to flame component lookup - flame = object.__getattribute__(self, "flame") - - if (not flame._has_component(name) and - flame._has_component(name.replace("_", "-"))): - name = name.replace("_", "-") - if flame._has_component(name): - return self.get_values(flame, name) - - raise AttributeError( - f"{flame.__class__.__name__!r} object has no attribute {name!r}") + warnings.warn("To be removed after Cantera 3.2. " + "Replaceable by 'Domain1D.oxidizer_velocity'.", + DeprecationWarning) + return self.flame.get_values("Uo") @property def left_control_point_temperature(self): @@ -1248,10 +1225,10 @@ def strain_rate(self, definition, fuel=None, oxidizer='O2', stoich=None): return np.abs(np.interp(z_stoich, self.grid, d_u_d_z)) elif definition == 'potential_flow_fuel': - return np.sqrt(- self.L[0] / self.density[0]) + return np.sqrt(- self.flame.radial_pressure_gradient[0] / self.density[0]) elif definition == 'potential_flow_oxidizer': - return np.sqrt(- self.L[0] / self.density[-1]) + return np.sqrt(- self.flame.radial_pressure_gradient[0] / self.density[-1]) else: raise ValueError('Definition "' + definition + '" is not available') diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml index acb61bcc3af..40bc212e843 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml @@ -20,6 +20,8 @@ recipes: - name: componentName - name: componentIndex wraps: componentIndex(const char*) +- name: getValues +- name: setValues - name: setBounds - name: lowerBound - name: upperBound From 52c7d7f0ad3d3241b5ed29bd6f8f6a9bbaed50f2 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 5 Oct 2025 13:51:57 -0500 Subject: [PATCH 12/22] [unittest] Update oneD unit tests --- test/python/test_onedim.py | 67 ++++++++++++++++++++------------------ 1 file changed, 36 insertions(+), 31 deletions(-) diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index 7696c5155cf..6311049ef2b 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -253,26 +253,26 @@ def test_fixed_temperature(self): assert self.sim.fixed_temperature_location == approx(zfixed) # test some alias names - assert np.allclose(self.sim.T, self.sim.temperature) + assert np.allclose(self.sim.T, self.sim.flame.temperature) zeros = self.sim.T * 0 with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): - self.sim.set_values(self.sim.flame, "spread_rate", zeros) + self.sim.flame.set_values("spread_rate", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): - self.sim.set_values(self.sim.flame, "radial-pressure-gradient", zeros) + self.sim.flame.set_values("radial-pressure-gradient", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): - self.sim.set_values(self.sim.flame, "L", zeros) + self.sim.flame.set_values("L", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): - self.sim.set_values(self.sim.flame, "Lambda", zeros) + self.sim.flame.set_values("Lambda", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): - self.sim.set_values(self.sim.flame, "electric-field", zeros) + self.sim.flame.set_values("electric-field", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): - self.sim.set_values(self.sim.flame, "E", zeros) + self.sim.flame.set_values("E", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): - self.sim.set_values(self.sim.flame, "eField", zeros) + self.sim.flame.set_values("eField", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): - self.sim.set_values(self.sim.flame, "oxidizer-velocity", zeros) + self.sim.flame.set_values("oxidizer-velocity", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): - self.sim.set_values(self.sim.flame, "Uo", zeros) + self.sim.flame.set_values("Uo", zeros) @pytest.mark.slow_test def test_auto_width(self): @@ -1297,16 +1297,17 @@ def run_restore_diffusionflame(self, fname): gas = ct.Solution("h2o2.yaml") sim = ct.CounterflowDiffusionFlame(gas) sim.restore(fname) - assert np.isclose(np.unique(sim.radial_pressure_gradient)[0], -90257.82873678) - assert np.allclose(sim.T, sim.temperature) - assert np.allclose(sim.radial_pressure_gradient, sim.L) - assert np.allclose(sim.radial_pressure_gradient, sim.Lambda) - assert np.allclose(sim.radial_pressure_gradient, getattr(sim, "lambda")) + radial_lambda = sim.flame.radial_pressure_gradient + assert np.isclose(np.unique(radial_lambda)[0], -1.355348889404357e+05) + assert np.allclose(sim.T, sim.flame.temperature) + assert np.allclose(radial_lambda, sim.flame.L) + assert np.allclose(radial_lambda, sim.flame.Lambda) + assert np.allclose(radial_lambda, getattr(sim.flame, "lambda")) zeros = sim.T * 0 with pytest.raises(ct.CanteraError, match="is not used by 'axisymmetric-flow'"): - sim.set_values(sim.flame, "electric-field", zeros) + sim.flame.set_values("electric-field", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'axisymmetric-flow'"): - sim.set_values(sim.flame, "oxidizer-velocity", zeros) + sim.flame.set_values("oxidizer-velocity", zeros) def run_restore_solutionarray(self, fname): gas = ct.Solution("h2o2.yaml") @@ -1315,7 +1316,7 @@ def run_restore_solutionarray(self, fname): radial_lambda = np.unique(arr.radial_pressure_gradient) assert len(radial_lambda) == 1 - assert radial_lambda[0] == pytest.approx(-90257.82873678) + assert radial_lambda[0] == pytest.approx(-1.355348889404357e+05) assert getattr(arr, "lambda")[0] == radial_lambda[0] assert radial_lambda[0] == arr.L[0] assert radial_lambda[0] == arr.Lambda[0] @@ -1399,7 +1400,7 @@ def test_two_point_control(self): sim.right_control_point_temperature += temperature_decrement sim.left_control_point_temperature += temperature_decrement sim.solve(loglevel=0, refine_grid=False) - assert sim.oxidizer_velocity == approx(sim.velocity[-1]) + assert sim.flame.oxidizer_velocity == approx(sim.velocity[-1]) temperature_4 = sim.T # Check the difference between the un-perturbed two-point solution and the @@ -1444,7 +1445,7 @@ def test_two_point_control(self): assert (sim.right_control_point_temperature == approx(original_settings['right-temperature'], 1e-4)) - assert sim.oxidizer_velocity == approx(sim.velocity[-1]) + assert sim.flame.oxidizer_velocity == approx(sim.velocity[-1]) # Test - Check error conditions sim = ct.CounterflowDiffusionFlame(gas, width=width) @@ -1561,7 +1562,8 @@ def _run_case(phi, T, width, P): sim.solve(loglevel=0, auto=True) assert all(sim.T >= T - 1e-3) assert all(sim.spread_rate >= -1e-9) - assert np.allclose(sim.radial_pressure_gradient, sim.radial_pressure_gradient[0]) + radial_lambda = sim.flame.radial_pressure_gradient + assert np.allclose(radial_lambda, radial_lambda[0]) return sim return _run_case @@ -1678,7 +1680,8 @@ def _run_case(phi, T, width, P): sim.solve(loglevel=0, auto=True) assert all(sim.T >= T - 1e-3) assert all(sim.spread_rate >= -1e-9) - assert np.allclose(sim.radial_pressure_gradient, sim.radial_pressure_gradient[0]) + radial_lambda = sim.flame.radial_pressure_gradient + assert np.allclose(radial_lambda, radial_lambda[0]) return sim return _run_case @@ -1725,7 +1728,7 @@ def _run_case(phi, T, width, P): sim.burner.mdot = gas.density * 0.15 sim.solve(loglevel=0, auto=True) assert sim.T[1] > T - assert np.allclose(sim.radial_pressure_gradient, 0) + assert np.allclose(sim.flame.radial_pressure_gradient, 0) return _run_case @pytest.mark.parametrize( @@ -1818,7 +1821,8 @@ def run_stagnation(self, xh2, mdot, width): sim.solve(loglevel=0, auto=True) assert sim.T.max() > tburner + tsurf - assert np.allclose(sim.radial_pressure_gradient, sim.radial_pressure_gradient[0]) + radial_lambda = sim.flame.radial_pressure_gradient + assert np.allclose(radial_lambda, radial_lambda[0]) self.sim = sim def test_stagnation_case1(self): @@ -2005,7 +2009,8 @@ def solve(self, phi, T, width, P): sim.reactants.mdot = gas.density * axial_velocity sim.solve(loglevel=0, auto=True) assert sim.T[-1] > T + 100 - assert np.allclose(sim.radial_pressure_gradient, sim.radial_pressure_gradient[0]) + radial_lambda = sim.flame.radial_pressure_gradient + assert np.allclose(radial_lambda, radial_lambda[0]) return sim def test_restart(self): @@ -2077,19 +2082,19 @@ def test_ion_profile(self): # stage one assert self.sim.electric_field_enabled is False self.sim.solve(loglevel=0, auto=True) - assert max(np.abs(self.sim.electric_field)) == approx(0.0, abs=1e-3) + assert max(np.abs(self.sim.flame.electric_field)) == approx(0.0, abs=1e-3) #stage two self.sim.electric_field_enabled = True self.sim.solve(loglevel=0) # Regression test - assert max(self.sim.electric_field) == approx(149.63179056676853, rel=1e-3) - with pytest.warns(DeprecationWarning, match="Renamed to 'electric_field'"): - assert np.allclose(self.sim.E, self.sim.electric_field) + assert max(self.sim.flame.electric_field) == approx(149.6317905667685, rel=1e-3) + with pytest.warns(DeprecationWarning, match="by 'Domain1D.electric_field'"): + assert np.allclose(self.sim.E, self.sim.flame.electric_field) zeros = self.sim.T with pytest.raises(ct.CanteraError, match="is not used by 'free-ion-flow'"): - self.sim.set_values(self.sim.flame, "oxidizer-velocity", zeros) + self.sim.flame.set_values("oxidizer-velocity", zeros) class TestIonBurnerFlame: @@ -2111,5 +2116,5 @@ def test_ion_profile(self): self.sim.solve(loglevel=0, auto=True) # Regression test - assert max(self.sim.electric_field) == approx(591.76, rel=1e-2) + assert max(self.sim.flame.electric_field) == approx(591.76, rel=1e-2) assert max(self.sim.X[self.gas.species_index('E')]) == approx(8.024e-9, rel=1e-2) From bf7d64f5d866d5299821128242db3e8188bc8537 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 5 Oct 2025 14:05:02 -0500 Subject: [PATCH 13/22] [samples] Simplify profile setting --- .../python/onedim/diffusion_flame_batch.py | 25 ++++++++----------- .../onedim/diffusion_flame_extinction.py | 13 ++++------ 2 files changed, 15 insertions(+), 23 deletions(-) diff --git a/samples/python/onedim/diffusion_flame_batch.py b/samples/python/onedim/diffusion_flame_batch.py index 69a8210871c..75dfc55931f 100644 --- a/samples/python/onedim/diffusion_flame_batch.py +++ b/samples/python/onedim/diffusion_flame_batch.py @@ -16,7 +16,7 @@ This example can, for example, be used to iterate to a counterflow diffusion flame to an awkward pressure and strain rate, or to create the basis for a flamelet table. -Requires: cantera >= 3.0, matplotlib >= 2.0 +Requires: cantera >= 3.2, matplotlib >= 2.0 .. tags:: Python, combustion, 1D flow, extinction, diffusion flame, strained flame, saving output, plotting @@ -131,18 +131,16 @@ def names(test): rel_pressure_increase = p / p_previous # Update grid f.flame.grid *= rel_pressure_increase ** exp_d_p - normalized_grid = f.grid / (f.grid[-1] - f.grid[0]) # Update mass fluxes f.fuel_inlet.mdot *= rel_pressure_increase ** exp_mdot_p f.oxidizer_inlet.mdot *= rel_pressure_increase ** exp_mdot_p # Update velocities - f.set_profile('velocity', normalized_grid, - f.velocity * rel_pressure_increase ** exp_u_p) - f.set_profile('spread_rate', normalized_grid, - f.spread_rate * rel_pressure_increase ** exp_V_p) + f.flame.set_values("velocity", f.flame.velocity * rel_pressure_increase ** exp_u_p) + f.flame.set_values( + "spread_rate", f.flame.spread_rate * rel_pressure_increase ** exp_V_p) # Update pressure curvature - f.set_profile("Lambda", normalized_grid, - f.radial_pressure_gradient * rel_pressure_increase ** exp_lam_p) + f.flame.set_values( + "Lambda", f.flame.radial_pressure_gradient * rel_pressure_increase ** exp_lam_p) try: # Try solving the flame @@ -186,18 +184,15 @@ def names(test): # Create an initial guess based on the previous solution # Update grid f.flame.grid *= strain_factor ** exp_d_a - normalized_grid = f.grid / (f.grid[-1] - f.grid[0]) # Update mass fluxes f.fuel_inlet.mdot *= strain_factor ** exp_mdot_a f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a # Update velocities - f.set_profile('velocity', normalized_grid, - f.velocity * strain_factor ** exp_u_a) - f.set_profile('spread_rate', normalized_grid, - f.spread_rate * strain_factor ** exp_V_a) + f.flame.set_values("velocity", f.flame.velocity * strain_factor ** exp_u_a) + f.flame.set_values("spread_rate", f.flame.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature - f.set_profile("Lambda", normalized_grid, - f.radial_pressure_gradient * strain_factor ** exp_lam_a) + f.flame.set_values( + "Lambda", f.flame.radial_pressure_gradient * strain_factor ** exp_lam_a) try: # Try solving the flame f.solve(loglevel=0) diff --git a/samples/python/onedim/diffusion_flame_extinction.py b/samples/python/onedim/diffusion_flame_extinction.py index 074f4b5e51b..8e6df42b722 100644 --- a/samples/python/onedim/diffusion_flame_extinction.py +++ b/samples/python/onedim/diffusion_flame_extinction.py @@ -12,7 +12,7 @@ (doi:10.1155/2014/484372). Please refer to this publication for a detailed explanation. Also, please don't forget to cite it if you make use of it. -Requires: cantera >= 3.0, matplotlib >= 2.0 +Requires: cantera >= 3.2, matplotlib >= 2.0 .. tags:: Python, combustion, 1D flow, diffusion flame, strained flame, extinction, saving output, plotting @@ -124,18 +124,15 @@ def names(test): # Update grid # Note that grid scaling changes the diffusion flame width f.flame.grid *= strain_factor ** exp_d_a - normalized_grid = f.grid / (f.grid[-1] - f.grid[0]) # Update mass fluxes f.fuel_inlet.mdot *= strain_factor ** exp_mdot_a f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a # Update velocities - f.set_profile('velocity', normalized_grid, - f.velocity * strain_factor ** exp_u_a) - f.set_profile('spread_rate', normalized_grid, - f.spread_rate * strain_factor ** exp_V_a) + f.flame.set_values("velocity", f.flame.velocity * strain_factor ** exp_u_a) + f.flame.set_values('spread_rate', f.flame.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature - f.set_profile("Lambda", normalized_grid, - f.radial_pressure_gradient * strain_factor ** exp_lam_a) + f.flame.set_values( + "Lambda", f.flame.radial_pressure_gradient * strain_factor ** exp_lam_a) try: f.solve(loglevel=0) except ct.CanteraError as e: From 1df7042d1665301dc2a9931724ee6316fdc11086 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 5 Oct 2025 15:37:39 -0500 Subject: [PATCH 14/22] [Python] Use grid setters/getters --- interfaces/cython/cantera/_onedim.pxd | 6 +++--- interfaces/cython/cantera/_onedim.pyx | 27 +++++++++++++-------------- 2 files changed, 16 insertions(+), 17 deletions(-) diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 8374031cd41..f571ff920a8 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -24,9 +24,11 @@ cdef extern from "cantera/oneD/Domain1D.h": size_t nPoints() string componentName(size_t) except +translate_exception size_t componentIndex(string&) except +translate_exception + cbool hasComponent(string&) except +translate_exception + vector[double] grid() except +translate_exception + void setupGrid(vector[double]&) except +translate_exception void getValues(string&, vector[double]&) except +translate_exception void setValues(string&, vector[double]&) except +translate_exception - cbool hasComponent(string&) except +translate_exception void setBounds(size_t, double, double) double upperBound(size_t) double lowerBound(size_t) @@ -40,8 +42,6 @@ cdef extern from "cantera/oneD/Domain1D.h": double steady_atol(size_t) double transient_rtol(size_t) double transient_atol(size_t) - double z(size_t) - void setupGrid(size_t, double*) except +translate_exception void setID(string) string& id() string domainType "type"() diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index cbadbeaa716..4e211543bb6 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -90,6 +90,19 @@ cdef class Domain1D: """Check whether `Domain1D` has component""" return self.domain.hasComponent(stringify(name)) + @property + def grid(self): + """The grid for this domain.""" + cdef vector[double] grid_vec = self.domain.grid() + return np.asarray(grid_vec) + + @grid.setter + def grid(self, grid): + cdef vector[double] grid_vec + for g in grid: + grid_vec.push_back(g) + self.domain.setupGrid(grid_vec) + def get_values(self, str component): """ Retrieve spatial profile of a component. @@ -290,20 +303,6 @@ cdef class Domain1D: else: return self.domain.transient_atol(self.component_index(component)) - property grid: - """ The grid for this domain """ - def __get__(self): - cdef np.ndarray[np.double_t, ndim=1] grid = np.empty(self.n_points) - cdef int i - for i in range(self.n_points): - grid[i] = self.domain.z(i) - return grid - - def __set__(self, grid): - cdef np.ndarray[np.double_t, ndim=1] data = \ - np.ascontiguousarray(grid, dtype=np.double) - self.domain.setupGrid(len(data), &data[0]) - property name: """ The name / id of this domain """ def __get__(self): From fda66cfef40e8243c96c65639f5ff0f13533d62f Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 5 Oct 2025 16:22:37 -0500 Subject: [PATCH 15/22] [oneD] Make setProfile accessible from Domain1D --- include/cantera/oneD/Domain1D.h | 51 +++++++++++++- include/cantera/oneD/Flow1D.h | 3 + include/cantera/oneD/OneDim.h | 41 ++++++++--- include/cantera/oneD/Sim1D.h | 30 +++----- interfaces/cython/cantera/_onedim.pxd | 2 + interfaces/cython/cantera/_onedim.pyx | 46 +++++++++++++ interfaces/cython/cantera/onedim.py | 69 ++++++++++--------- .../src/sourcegen/headers/ctdomain_auto.yaml | 4 ++ .../src/sourcegen/headers/ctonedim_auto.yaml | 4 +- src/oneD/Domain1D.cpp | 20 ++++++ src/oneD/Flow1D.cpp | 43 +++++++++++- src/oneD/OneDim.cpp | 28 ++++++++ src/oneD/Sim1D.cpp | 20 ++---- 13 files changed, 278 insertions(+), 83 deletions(-) diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index 17e99dab014..f3ea1c47ae5 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -360,7 +360,7 @@ class Domain1D } /** - * Retrieve component values from the solution vector. + * Retrieve component values. * @param component name of component * @param[out] values vector of values * @@ -369,7 +369,7 @@ class Domain1D void getValues(const string& component, vector& values) const; /** - * Specify component values in the solution vector. + * Specify component values. * @param component name of component * @param[in] values vector of values * @@ -377,6 +377,33 @@ class Domain1D */ void setValues(const string& component, const vector& values); + /** + * Specify a profile for a component. + * @param component name of component + * @param[in] pos A vector of relative positions, beginning with 0.0 at the + * left of the domain, and ending with 1.0 at the right of the domain. + * @param[in] values A vector of values corresponding to the relative position + * locations. + * + * Note that the vector pos and values can have lengths different than the + * number of grid points, but their lengths must be equal. The values at + * the grid points will be linearly interpolated based on the (pos, + * values) specification. + * + * @since New in %Cantera 3.2. + */ + void setProfile(const string& component, + const vector& pos, const vector& values); + + /** + * Specify a flat profile for a component. + * @param component name of component + * @param value constant value + * + * @since New in %Cantera 3.2. + */ + void setFlatProfile(const string& component, double v); + /** * Retrieve component values from the solution vector. * @since New in %Cantera 3.2. @@ -395,6 +422,24 @@ class Domain1D throw NotImplementedError("Domain1D::setValues", "Needs to be overloaded."); } + /** + * Specify a profile for a component in the solution vector. + * @since New in %Cantera 3.2. + */ + virtual void _setProfile(double* soln, const string& component, + const vector& pos, const vector& values) { + throw NotImplementedError("Domain1D::setProfile", "Needs to be overloaded."); + } + + /** + * Specify a flat profile for a component in the solution vector. + * @since New in %Cantera 3.2. + */ + virtual void _setFlatProfile(double* soln, const string& component, double v) { + throw NotImplementedError( + "Domain1D::setFlatProfile", "Needs to be overloaded."); + } + //! Save the state of this domain as a SolutionArray. /*! * @param soln local solution vector for this domain @@ -546,6 +591,8 @@ class Domain1D //! @param name Name of the component //! @param values Array of length nPoints() containing the initial values //! @param soln Pointer to the local portion of the system state vector + //! @deprecated To be removed after %Cantera 3.2. + //! Replaceable by version using vectors arguments. void setProfile(const string& name, double* values, double* soln); //! Access the array of grid coordinates [m] diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index 86eebdaf4e5..6fbf424468f 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -188,6 +188,9 @@ class Flow1D : public Domain1D vector& values) const override; void _setValues(double* soln, const string& component, const vector& values) override; + void _setProfile(double* soln, const string& component, + const vector& pos, const vector& values) override; + void _setFlatProfile(double* soln, const string& component, double v) override; shared_ptr asArray(const double* soln) const override; void fromArray(SolutionArray& arr, double* soln) override; diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index a1880ec32a1..eec0e58be89 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -149,10 +149,8 @@ class OneDim : public SteadyStateSystem * * @since New in %Cantera 3.2. */ - virtual void getValues(const string& dom, const string& component, - vector& values) const { - throw NotImplementedError("OneDim::getValues", "Needs to be overloaded."); - } + void getValues(const string& dom, const string& component, + vector& values) const; /** * Specify component values in the solution vector. @@ -162,10 +160,37 @@ class OneDim : public SteadyStateSystem * * @since New in %Cantera 3.2. */ - virtual void setValues(const string& dom, const string& component, - const vector& values) { - throw NotImplementedError("OneDim::setValues", "Needs to be overloaded."); - } + void setValues(const string& dom, const string& component, + const vector& values); + + /** + * Specify a profile for one component of one domain. + * @param dom name of domain + * @param component name of component + * @param[in] pos A vector of relative positions, beginning with 0.0 at the + * left of the domain, and ending with 1.0 at the right of the domain. + * @param[in] values A vector of values corresponding to the relative position + * locations. + * + * Note that the vector pos and values can have lengths different than the + * number of grid points, but their lengths must be equal. The values at + * the grid points will be linearly interpolated based on the (pos, + * values) specification. + * + * @since New in %Cantera 3.2. + */ + void setProfile(const string& dom, const string& component, + const vector& pos, const vector& values); + + /** + * Specify a flat profile for one component of one domain. + * @param dom name of domain + * @param component name of component + * @param value constant value + * + * @since New in %Cantera 3.2. + */ + void setFlatProfile(const string& dom, const string& component, double v); /** * Initialize all domains. On the first call, this methods calls the init diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index 4c225970700..cacd65e255f 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/include/cantera/oneD/Sim1D.h @@ -49,6 +49,9 @@ class Sim1D : public OneDim * left of the domain, and ending with 1.0 at the right of the domain. * @param vals A vector of values corresponding to the relative position * locations. + * + * @deprecated To be removed after %Cantera 3.2. + * Replaceable by Domain1D::setProfile. */ void setInitialGuess(const string& component, vector& locs, vector& vals); @@ -79,28 +82,6 @@ class Sim1D : public OneDim //! @param localPoint grid point within the domain double workValue(size_t dom, size_t comp, size_t localPoint) const; - /** - * Retrieve component values from the solution vector. - * @param dom name of domain - * @param component name of component - * @param[out] values vector of values - * - * @since New in %Cantera 3.2. - */ - void getValues(const string& dom, const string& component, - vector& values) const override; - - /** - * Specify component values in the solution vector. - * @param dom name of domain - * @param component name of component - * @param[in] values vector of values - * - * @since New in %Cantera 3.2. - */ - void setValues(const string& dom, const string& component, - const vector& values) override; - /** * Specify a profile for one component of one domain. * @param dom domain number, beginning with 0 for the leftmost domain. @@ -114,11 +95,16 @@ class Sim1D : public OneDim * number of grid points, but their lengths must be equal. The values at * the grid points will be linearly interpolated based on the (pos, * values) specification. + * + * @deprecated To be removed after %Cantera 3.2. + * Replaceable by Domain1D::setProfile. */ void setProfile(size_t dom, size_t comp, const vector& pos, const vector& values); //! Set component 'comp' of domain 'dom' to value 'v' at all points. + //! @deprecated To be removed after %Cantera 3.2. + //! Replaceable by Domain1D::setProfile. void setFlatProfile(size_t dom, size_t comp, double v); //! @} diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index f571ff920a8..b412ba22295 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -29,6 +29,8 @@ cdef extern from "cantera/oneD/Domain1D.h": void setupGrid(vector[double]&) except +translate_exception void getValues(string&, vector[double]&) except +translate_exception void setValues(string&, vector[double]&) except +translate_exception + void setProfile(string&, vector[double]&, vector[double]&) except +translate_exception + void setFlatProfile(string&, double) except +translate_exception void setBounds(size_t, double, double) double upperBound(size_t) double lowerBound(size_t) diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 4e211543bb6..64e1efe06de 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -137,6 +137,44 @@ cdef class Domain1D: values_vec.push_back(v) self.domain.setValues(stringify(component), values) + def set_profile(self, component, positions, values): + """ + Set an initial estimate for a profile of one component in one domain. + + :param component: + component name + :param positions: + sequence of relative positions, from 0 on the left to 1 on the right + :param values: + sequence of values at the relative positions specified in ``positions`` + + >>> d.set_profile('T', [0.0, 0.2, 1.0], [400.0, 800.0, 1500.0]) + + .. versionadded:: 3.2 + """ + cdef vector[double] pos_vec, val_vec + for p in positions: + pos_vec.push_back(p) + for v in values: + val_vec.push_back(v) + + self.domain.setProfile(stringify(component), pos_vec, val_vec) + + def set_flat_profile(self, component, value): + """ + Set a flat profile for a component. + + :param component: + component name + :param v: + value + + >>> d.set_flat_profile('u', -3.0) + + .. versionadded:: 3.2 + """ + self.domain.setFlatProfile(stringify(component), value) + def set_bounds(self, *, default=None, Y=None, **kwargs): """ Set the lower and upper bounds on the solution. @@ -1078,6 +1116,10 @@ cdef class Sim1D: sequence of values at the relative positions specified in ``positions`` >>> s.set_profile(d, 'T', [0.0, 0.2, 1.0], [400.0, 800.0, 1500.0]) + + .. deprecated:: 3.2 + + To be removed after Cantera 3.2. Replaceable by Domain1D.set_profile. """ dom, comp = self._get_indices(domain, component) @@ -1100,6 +1142,10 @@ cdef class Sim1D: value >>> s.set_flat_profile(d, 'u', -3.0) + + .. deprecated:: 3.2 + + To be removed after Cantera 3.2. Replaceable by Domain1D.set_flat_profile. """ dom, comp = self._get_indices(domain, component) self.sim.setFlatProfile(dom, comp, value) diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index 3563ca3b23f..16f411b6ab9 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -169,6 +169,10 @@ def set_profile(self, component, positions, values): sequence of values at the relative positions specified in ``positions`` >>> f.set_profile('T', [0.0, 0.2, 1.0], [400.0, 800.0, 1500.0]) + + .. deprecated:: 3.2 + + To be removed after Cantera 3.2. Replaceable by Domain1D.set_profile. """ super().set_profile(self.flame, component, positions, values) @@ -690,8 +694,8 @@ def set_initial_guess(self, locs=[0.0, 0.3, 0.5, 1.0], data=None, group=None): Yeq = self.gas.Y u1 = self.inlet.mdot / self.gas.density - self.set_profile('velocity', locs, [u0, u0, u1, u1]) - self.set_profile('T', locs, [T0, T0, Teq, Teq]) + self.flame.set_profile('velocity', locs, [u0, u0, u1, u1]) + self.flame.set_profile('T', locs, [T0, T0, Teq, Teq]) # Pick the location of the fixed temperature point, using an existing # point if a reasonable choice exists @@ -706,7 +710,7 @@ def set_initial_guess(self, locs=[0.0, 0.3, 0.5, 1.0], data=None, group=None): self.fixed_temperature = Tmid for n in range(self.gas.n_species): - self.set_profile(self.gas.species_name(n), + self.flame.set_profile(self.gas.species_name(n), locs, [Y0[n], Y0[n], Yeq[n], Yeq[n]]) def solve(self, loglevel=1, refine_grid=True, auto=False, stage=None): @@ -885,10 +889,10 @@ def set_initial_guess(self, data=None, group=None): u1 = self.burner.mdot / self.gas.density locs = [0.0, 0.2, 1.0] - self.set_profile('velocity', locs, [u0, u1, u1]) - self.set_profile('T', locs, [T0, Teq, Teq]) + self.flame.set_profile('velocity', locs, [u0, u1, u1]) + self.flame.set_profile('T', locs, [T0, Teq, Teq]) for n in range(self.gas.n_species): - self.set_profile(self.gas.species_name(n), + self.flame.set_profile(self.gas.species_name(n), locs, [Y0[n], Yeq[n], Yeq[n]]) def solve(self, loglevel=1, refine_grid=True, auto=False, stage=None): @@ -950,9 +954,9 @@ def check_blowoff(t): except FlameBlowoff: # The eventual solution for a blown off flame is the non-reacting # solution, so just set the state to this now - self.set_flat_profile(self.flame, 'T', self.T[0]) + self.flame.set_flat_profile("T", self.T[0]) for k,spec in enumerate(self.gas.species_names): - self.set_flat_profile(self.flame, spec, self.burner.Y[k]) + self.flame.set_flat_profile(spec, self.burner.Y[k]) self.set_steady_callback(original_callback) super().solve(loglevel, False, False) @@ -1071,12 +1075,12 @@ def set_initial_guess(self, data=None, group=None): T[-1] = T0o zrel = (zz - zz[0])/dz - self.set_profile('velocity', [0.0, 1.0], [u0f, -u0o]) - self.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) - self.set_profile("Lambda", [0.0, 1.0], [L, L]) - self.set_profile('T', zrel, T) + self.flame.set_profile('velocity', [0.0, 1.0], [u0f, -u0o]) + self.flame.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) + self.flame.set_profile("Lambda", [0.0, 1.0], [L, L]) + self.flame.set_profile('T', zrel, T) for k,spec in enumerate(self.gas.species_names): - self.set_profile(spec, zrel, Y[:,k]) + self.flame.set_profile(spec, zrel, Y[:,k]) def extinct(self): return max(self.T) - max(self.fuel_inlet.T, self.oxidizer_inlet.T) < 10 @@ -1364,20 +1368,19 @@ def set_initial_guess(self, products='inlet', data=None, group=None): Teq = self.gas.T Yeq = self.gas.Y locs = np.array([0.0, 0.3, 0.7, 1.0]) - self.set_profile('T', locs, [T0, Teq, Teq, self.surface.T]) + self.flame.set_profile('T', locs, [T0, Teq, Teq, self.surface.T]) for k in range(self.gas.n_species): - self.set_profile(self.gas.species_name(k), locs, - [Y0[k], Yeq[k], Yeq[k], Yeq[k]]) + self.flame.set_profile(self.gas.species_name(k), locs, + [Y0[k], Yeq[k], Yeq[k], Yeq[k]]) else: locs = np.array([0.0, 1.0]) - self.set_profile('T', locs, [T0, self.surface.T]) + self.flame.set_profile('T', locs, [T0, self.surface.T]) for k in range(self.gas.n_species): - self.set_profile(self.gas.species_name(k), locs, - [Y0[k], Y0[k]]) + self.flame.set_profile(self.gas.species_name(k), locs, [Y0[k], Y0[k]]) locs = np.array([0.0, 1.0]) - self.set_profile('velocity', locs, [u0, 0.0]) - self.set_profile('spread_rate', locs, [0.0, 0.0]) + self.flame.set_profile('velocity', locs, [u0, 0.0]) + self.flame.set_profile('spread_rate', locs, [0.0, 0.0]) class CounterflowPremixedFlame(FlameBase): @@ -1465,10 +1468,10 @@ def set_initial_guess(self, equilibrate=True, data=None, group=None): "must be positive") locs = np.array([0.0, 0.4, 0.6, 1.0]) - self.set_profile('T', locs, [Tu, Tu, Teq, Tb]) + self.flame.set_profile('T', locs, [Tu, Tu, Teq, Tb]) for k in range(self.gas.n_species): - self.set_profile(self.gas.species_name(k), locs, - [Yu[k], Yu[k], Yeq[k], Yb[k]]) + self.flame.set_profile(self.gas.species_name(k), locs, + [Yu[k], Yu[k], Yeq[k], Yb[k]]) # estimate strain rate self.gas.TPY = Teq, self.flame.P, Yeq @@ -1479,9 +1482,9 @@ def set_initial_guess(self, equilibrate=True, data=None, group=None): # estimate stagnation point x0 = rhou*uu * dz / (rhou*uu + rhob*ub) - self.set_profile('velocity', [0.0, 1.0], [uu, -ub]) - self.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) - self.set_profile("Lambda", [0.0, 1.0], [L, L]) + self.flame.set_profile('velocity', [0.0, 1.0], [uu, -ub]) + self.flame.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) + self.flame.set_profile("Lambda", [0.0, 1.0], [L, L]) class CounterflowTwinPremixedFlame(FlameBase): @@ -1552,10 +1555,10 @@ def set_initial_guess(self, data=None, group=None): Yb = self.gas.Y locs = np.array([0.0, 0.4, 0.6, 1.0]) - self.set_profile('T', locs, [Tu, Tu, Tb, Tb]) + self.flame.set_profile('T', locs, [Tu, Tu, Tb, Tb]) for k in range(self.gas.n_species): - self.set_profile(self.gas.species_name(k), locs, - [Yu[k], Yu[k], Yb[k], Yb[k]]) + self.flame.set_profile(self.gas.species_name(k), locs, + [Yu[k], Yu[k], Yb[k], Yb[k]]) # estimate strain rate zz = self.flame.grid @@ -1563,6 +1566,6 @@ def set_initial_guess(self, data=None, group=None): a = 2 * uu / dz L = - rhou * a**2 - self.set_profile('velocity', [0.0, 1.0], [uu, 0]) - self.set_profile('spread_rate', [0.0, 1.0], [0.0, a]) - self.set_profile("Lambda", [0.0, 1.0], [L, L]) + self.flame.set_profile('velocity', [0.0, 1.0], [uu, 0]) + self.flame.set_profile('spread_rate', [0.0, 1.0], [0.0, a]) + self.flame.set_profile("Lambda", [0.0, 1.0], [L, L]) diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml index 40bc212e843..1b1067ef9b5 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml @@ -22,6 +22,10 @@ recipes: wraps: componentIndex(const char*) - name: getValues - name: setValues +- name: setProfile + # alternative is deprecated in Cantera 3.2 + wraps: setProfile(const string&, const vector&, const vector&) +- name: setFlatProfile - name: setBounds - name: lowerBound - name: upperBound diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml index cbaf19e6bbe..084bca26a91 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml @@ -10,11 +10,9 @@ recipes: - name: newSim1D uses: domain - name: setValue -- name: getValues -- name: setValues - name: setProfile - name: setFlatProfile -- name: setInitialGuess # New in Cantera 3.2 +- name: setInitialGuess - name: show - name: setTimeStep - name: getInitialSoln diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index 1d341483192..28e79e2deb8 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -111,6 +111,23 @@ void Domain1D::setValues(const string& component, const vector& values) m_container->setValues(id(), component, values); } +void Domain1D::setProfile(const string& component, + const vector& pos, const vector& values) +{ + if (!m_container) { + throw CanteraError("Domain1D::setProfile", "Domain is not connected."); + } + m_container->setProfile(id(), component, pos, values); +} + +void Domain1D::setFlatProfile(const string& component, double v) +{ + if (!m_container) { + throw CanteraError("Domain1D::setFlatProfile", "Domain is not connected."); + } + m_container->setFlatProfile(id(), component, v); +} + void Domain1D::setTransientTolerances(double rtol, double atol, size_t n) { if (n == npos) { @@ -312,6 +329,9 @@ void Domain1D::show(const double* x) void Domain1D::setProfile(const string& name, double* values, double* soln) { + warn_deprecated( + "Domain1D::setProfile", "To be removed after Cantera 3.2. Replaceable by " + "version using vector arguments."); for (size_t n = 0; n < m_nv; n++) { if (name == componentName(n)) { for (size_t j = 0; j < m_points; j++) { diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 8bcfb49183b..5df36c8630e 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -963,7 +963,7 @@ void Flow1D::_getValues(const double* soln, const string& component, } void Flow1D::_setValues(double* soln, const string& component, - const vector& values) + const vector& values) { if (values.size() != nPoints()) { throw ArraySizeError("Flow1D::_setValues", values.size(), nPoints()); @@ -979,6 +979,47 @@ void Flow1D::_setValues(double* soln, const string& component, } } +void Flow1D::_setProfile(double* soln, const string& component, + const vector& pos, const vector& values) +{ + if (pos.size() != values.size()) { + throw CanteraError( + "Flow1D::_setProfile", "Vectors for positions and values must have same " + "size.\nSizes are {} and {}, respectively.", pos.size(), values.size()); + } + if (pos.front() != 0.0 || pos.back() != 1.0) { + throw CanteraError("Flow1D::_setProfile", + "'pos' vector must span the range [0, 1]. Got a vector spanning " + "[{}, {}] instead.", pos.front(), pos.back()); + } + auto i = componentIndex(component); + if (!componentActive(i)) { + throw CanteraError( + "Flow1D::_setProfile", "Component '{}' is not used by '{}'.", + component, domainType()); + } + double z0 = zmin(); + double zDelta = zmax() - z0; + for (size_t j = 0; j < nPoints(); j++) { + double frac = (z(j) - z0)/zDelta; + double v = linearInterp(frac, pos, values); + soln[index(i,j)] = v; + } +} + +void Flow1D::_setFlatProfile(double* soln, const string& component, double v) +{ + auto i = componentIndex(component); + if (!componentActive(i)) { + throw CanteraError( + "Flow1D::_setFlatProfile", "Component '{}' is not used by '{}'.", + component, domainType()); + } + for (size_t j = 0; j < nPoints(); j++) { + soln[index(i,j)] = v; + } +} + shared_ptr Flow1D::asArray(const double* soln) const { auto arr = SolutionArray::create( diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 2b8eec13440..87f30228c32 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -72,6 +72,34 @@ double OneDim::lowerBound(size_t i) const return dom.lowerBound(k); } +void OneDim::getValues(const string& dom, const string& component, + vector& values) const +{ + size_t dIdx = domainIndex(dom); + domain(dIdx)._getValues(m_state->data() + domain(dIdx).loc(), component, values); +} + +void OneDim::setValues(const string& dom, const string& component, + const vector& values) +{ + size_t dIdx = domainIndex(dom); + domain(dIdx)._setValues(m_state->data() + domain(dIdx).loc(), component, values); +} + +void OneDim::setProfile(const string& dom, const string& component, + const vector& pos, const vector& values) +{ + size_t dIdx = domainIndex(dom); + domain(dIdx)._setProfile(m_state->data() + domain(dIdx).loc(), component, + pos, values); +} + +void OneDim::setFlatProfile(const string& dom, const string& component, double v) +{ + size_t dIdx = domainIndex(dom); + domain(dIdx)._setFlatProfile(m_state->data() + domain(dIdx).loc(), component, v); +} + void OneDim::addDomain(shared_ptr d) { // if 'd' is not the first domain, link it to the last domain diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index a3251147ed5..d2c20abdef4 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -37,6 +37,8 @@ Sim1D::Sim1D(vector>& domains) : void Sim1D::setInitialGuess(const string& component, vector& locs, vector& vals) { + warn_deprecated("Sim1D::setInitialGuess", "To be removed after Cantera 3.2. " + "Replaceable by Domain1D::setProfile."); for (size_t dom=0; dom& values) const -{ - size_t dIdx = domainIndex(dom); - domain(dIdx)._getValues(m_state->data() + domain(dIdx).loc(), component, values); -} - -void Sim1D::setValues(const string& dom, const string& component, - const vector& values) -{ - size_t dIdx = domainIndex(dom); - domain(dIdx)._setValues(m_state->data() + domain(dIdx).loc(), component, values); -} - void Sim1D::setProfile(size_t dom, size_t comp, const vector& pos, const vector& values) { + warn_deprecated("Sim1D::setProfile", "To be removed after Cantera 3.2. " + "Replaceable by Domain1D::setProfile."); if (pos.front() != 0.0 || pos.back() != 1.0) { throw CanteraError("Sim1D::setProfile", "`pos` vector must span the range [0, 1]. Got a vector spanning " @@ -337,6 +327,8 @@ void Sim1D::_restore(const string& fname, const string& name) void Sim1D::setFlatProfile(size_t dom, size_t comp, double v) { + warn_deprecated("Sim1D::setFlatProfile", "To be removed after Cantera 3.2. " + "Replaceable by Domain1D::setProfile."); size_t np = domain(dom).nPoints(); for (size_t n = 0; n < np; n++) { setValue(dom, comp, n, v); From 03d6f2d2de6827952a07db5b21578509480604e8 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sun, 5 Oct 2025 17:01:39 -0500 Subject: [PATCH 16/22] [unittest/samples] Use updated methods --- interfaces/cython/cantera/onedim.py | 25 ------------------- .../src/sourcegen/headers/ctonedim_auto.yaml | 6 ++--- samples/cxx/flamespeed/flamespeed.cpp | 8 +++--- test/clib/test_ctonedim.cpp | 6 ++--- test/clib_legacy/test_ctonedim.cpp | 2 ++ test/oneD/test_oneD.cpp | 6 ++--- test/python/test_onedim.py | 9 +------ 7 files changed, 16 insertions(+), 46 deletions(-) diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index 16f411b6ab9..e83f4239ecb 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -297,16 +297,7 @@ def L(self): """ Array containing the radial pressure gradient (1/r)(dP/dr) [N/m⁴] at each point. Note: This value is named ``Lambda`` in the C++ code. - - .. deprecated:: 3.2 - - To be removed after Cantera 3.2. Replaceable by - `Domain1D.radial_pressure_gradient`. - Note: replacement may be shadowed by a species name ``L``. """ - warnings.warn("To be removed after Cantera 3.2. " - "Replaceable by 'Domain1D.radial_pressure_gradient'.", - DeprecationWarning) return self.flame.get_values("Lambda") @property @@ -315,14 +306,7 @@ def E(self): Array containing the electric field strength at each point. Note: This value is named ``eField`` in the C++ code and is only defined if the transport model is ``ionized_gas``. - - .. deprecated:: 3.2 - - To be removed after Cantera 3.2. Replaceable by `Domain1D.electric_field`. - Note: replacement may be shadowed by a species name ``E``. """ - warnings.warn("To be removed after Cantera 3.2. " - "Replaceable by 'Domain1D.electric_field'.", DeprecationWarning) return self.flame.get_values("eField") @property @@ -332,16 +316,7 @@ def Uo(self): each point. Note: This value is named ``Uo`` in the C++ code and is only defined when using two-point control. - - .. deprecated:: 3.2 - - To be removed after Cantera 3.2. Replaceable by - `Domain1D.oxidizer_velocity`. - Note: replacement may be shadowed by a species name ``Uo``. """ - warnings.warn("To be removed after Cantera 3.2. " - "Replaceable by 'Domain1D.oxidizer_velocity'.", - DeprecationWarning) return self.flame.get_values("Uo") @property diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml index 084bca26a91..b81fdd62942 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml @@ -10,9 +10,9 @@ recipes: - name: newSim1D uses: domain - name: setValue -- name: setProfile -- name: setFlatProfile -- name: setInitialGuess +- name: setProfile # deprecated +- name: setFlatProfile # deprecated +- name: setInitialGuess # deprecated - name: show - name: setTimeStep - name: getInitialSoln diff --git a/samples/cxx/flamespeed/flamespeed.cpp b/samples/cxx/flamespeed/flamespeed.cpp index 8b9b6daa938..39872894d2f 100644 --- a/samples/cxx/flamespeed/flamespeed.cpp +++ b/samples/cxx/flamespeed/flamespeed.cpp @@ -8,7 +8,7 @@ * * flamespeed [equivalence_ratio] [refine_grid] [loglevel] * - * Requires: cantera >= 3.1 + * Requires: cantera >= 3.2 * * .. tags:: C++, combustion, 1D flow, premixed flame, flame speed, saving output */ @@ -96,13 +96,13 @@ int flamespeed(double phi, bool refine_grid, int loglevel) double uout = inlet->mdot()/rho_out; value = {uin, uin, uout, uout}; - flame.setInitialGuess("velocity",locs,value); + flow->setProfile("velocity", locs, value); value = {temp, temp, Tad, Tad}; - flame.setInitialGuess("T",locs,value); + flow->setProfile("T", locs, value); for (size_t i=0; ispeciesName(i),locs,value); + flow->setProfile(gas->speciesName(i), locs, value); } inlet->setMoleFractions(x.data()); diff --git a/test/clib/test_ctonedim.cpp b/test/clib/test_ctonedim.cpp index 2a4413439d4..56cbc73825d 100644 --- a/test/clib/test_ctonedim.cpp +++ b/test/clib/test_ctonedim.cpp @@ -156,9 +156,9 @@ TEST(ctonedim, freeflame) // set up initial guess vector locs{0.0, 0.3, 0.7, 1.0}; vector value{uin, uin, uout, uout}; - sim1D_setInitialGuess(flame, "velocity", 4, locs.data(), 4, value.data()); + domain_setProfile(flow, "velocity", 4, locs.data(), 4, value.data()); value = {T, T, Tad, Tad}; - sim1D_setInitialGuess(flame, "T", 4, locs.data(), 4, value.data()); + domain_setProfile(flow, "T", 4, locs.data(), 4, value.data()); for (size_t i = 0; i < nsp; i++) { value = {yin[i], yin[i], yout[i], yout[i]}; int32_t buflen = thermo_speciesName(ph, i, 0, 0); @@ -166,7 +166,7 @@ TEST(ctonedim, freeflame) thermo_speciesName(ph, i, buflen, buf.data()); string name(buf.data()); ASSERT_EQ(name, gas->speciesName(i)); - sim1D_setInitialGuess(flame, buf.data(), 4, locs.data(), 4, value.data()); + domain_setProfile(flow, buf.data(), 4, locs.data(), 4, value.data()); } // simulation settings diff --git a/test/clib_legacy/test_ctonedim.cpp b/test/clib_legacy/test_ctonedim.cpp index 4a30c4da770..63580f71848 100644 --- a/test/clib_legacy/test_ctonedim.cpp +++ b/test/clib_legacy/test_ctonedim.cpp @@ -99,6 +99,7 @@ TEST(ctonedim, catcomb) TEST(ctonedim, freeflame) { ct_resetStorage(); + suppress_deprecation_warnings(); auto gas = newThermo("h2o2.yaml", "ohmech"); int sol = soln_newSolution("h2o2.yaml", "ohmech", "default"); @@ -203,4 +204,5 @@ TEST(ctonedim, freeflame) ASSERT_GE(T, Tprev); Tprev = T; } + make_deprecation_warnings_fatal(); } diff --git a/test/oneD/test_oneD.cpp b/test/oneD/test_oneD.cpp index 2bfe3703144..79fc66636fd 100644 --- a/test/oneD/test_oneD.cpp +++ b/test/oneD/test_oneD.cpp @@ -61,12 +61,12 @@ TEST(onedim, freeflame) // set up initial guess vector locs{0.0, 0.3, 0.7, 1.0}; vector value{uin, uin, uout, uout}; - flame->setInitialGuess("velocity", locs, value); + flow->setProfile("velocity", locs, value); value = {T, T, Tad, Tad}; - flame->setInitialGuess("T", locs, value); + flow->setProfile("T", locs, value); for (size_t i = 0; i < nsp; i++) { value = {yin[i], yin[i], yout[i], yout[i]}; - flame->setInitialGuess(gas->speciesName(i), locs, value); + flow->setProfile(gas->speciesName(i), locs, value); } // simulation settings diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index 6311049ef2b..211ad9001bd 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -723,7 +723,6 @@ def run_fixed_restore(self, mode): Y1 = self.sim.Y u1 = self.sim.velocity - V1 = self.sim.spread_rate P1 = self.sim.P T1 = self.sim.T self.sim.save(filename, "test") @@ -751,24 +750,20 @@ def run_fixed_restore(self, mode): Y2 = self.sim.Y u2 = self.sim.velocity - V2 = self.sim.spread_rate T2 = self.sim.T assert Y1 == approx(Y2) assert u1 == approx(u2) - assert V1 == approx(V2) assert T1 == approx(T2) self.solve_fixed_T() Y3 = self.sim.Y u3 = self.sim.velocity - V3 = self.sim.spread_rate # TODO: These tolerances seem too loose, but the tests fail on some # systems with tighter tolerances. assert Y1 == approx(Y3, rel=3e-3) assert u1 == approx(u3, rel=1e-3) - assert V1 == approx(V3, rel=1e-3) assert not self.sim.radiation_enabled assert not self.sim.soret_enabled @@ -1728,7 +1723,6 @@ def _run_case(phi, T, width, P): sim.burner.mdot = gas.density * 0.15 sim.solve(loglevel=0, auto=True) assert sim.T[1] > T - assert np.allclose(sim.flame.radial_pressure_gradient, 0) return _run_case @pytest.mark.parametrize( @@ -2090,8 +2084,7 @@ def test_ion_profile(self): # Regression test assert max(self.sim.flame.electric_field) == approx(149.6317905667685, rel=1e-3) - with pytest.warns(DeprecationWarning, match="by 'Domain1D.electric_field'"): - assert np.allclose(self.sim.E, self.sim.flame.electric_field) + assert np.allclose(self.sim.E, self.sim.flame.electric_field) zeros = self.sim.T with pytest.raises(ct.CanteraError, match="is not used by 'free-ion-flow'"): self.sim.flame.set_values("oxidizer-velocity", zeros) From dc648bd4a7762b56f8106446236648761919679e Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Mon, 6 Oct 2025 10:00:29 -0500 Subject: [PATCH 17/22] [oneD] Streamline access to work vector --- include/cantera/oneD/Domain1D.h | 49 +++++++------------------------ include/cantera/oneD/Flow1D.h | 12 ++++---- include/cantera/oneD/OneDim.h | 51 --------------------------------- src/oneD/Domain1D.cpp | 33 --------------------- src/oneD/Flow1D.cpp | 48 +++++++++++++++++++++---------- src/oneD/OneDim.cpp | 28 ------------------ 6 files changed, 49 insertions(+), 172 deletions(-) diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index f3ea1c47ae5..50314186a9d 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -366,7 +366,9 @@ class Domain1D * * @since New in %Cantera 3.2. */ - void getValues(const string& component, vector& values) const; + virtual void getValues(const string& component, vector& values) const { + throw NotImplementedError("Domain1D::setProfile", "Needs to be overloaded."); + } /** * Specify component values. @@ -375,7 +377,9 @@ class Domain1D * * @since New in %Cantera 3.2. */ - void setValues(const string& component, const vector& values); + virtual void setValues(const string& component, const vector& values) { + throw NotImplementedError("Domain1D::setProfile", "Needs to be overloaded."); + } /** * Specify a profile for a component. @@ -392,8 +396,10 @@ class Domain1D * * @since New in %Cantera 3.2. */ - void setProfile(const string& component, - const vector& pos, const vector& values); + virtual void setProfile(const string& component, + const vector& pos, const vector& values) { + throw NotImplementedError("Domain1D::setProfile", "Needs to be overloaded."); + } /** * Specify a flat profile for a component. @@ -402,40 +408,7 @@ class Domain1D * * @since New in %Cantera 3.2. */ - void setFlatProfile(const string& component, double v); - - /** - * Retrieve component values from the solution vector. - * @since New in %Cantera 3.2. - */ - virtual void _getValues(const double* soln, const string& component, - vector& values) const { - throw NotImplementedError("Domain1D::getValues", "Needs to be overloaded."); - } - - /** - * Specify component values in the solution vector. - * @since New in %Cantera 3.2. - */ - virtual void _setValues(double* soln, const string& component, - const vector& values) { - throw NotImplementedError("Domain1D::setValues", "Needs to be overloaded."); - } - - /** - * Specify a profile for a component in the solution vector. - * @since New in %Cantera 3.2. - */ - virtual void _setProfile(double* soln, const string& component, - const vector& pos, const vector& values) { - throw NotImplementedError("Domain1D::setProfile", "Needs to be overloaded."); - } - - /** - * Specify a flat profile for a component in the solution vector. - * @since New in %Cantera 3.2. - */ - virtual void _setFlatProfile(double* soln, const string& component, double v) { + virtual void setFlatProfile(const string& component, double v) { throw NotImplementedError( "Domain1D::setFlatProfile", "Needs to be overloaded."); } diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index 6fbf424468f..d49fb2f4807 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -184,13 +184,11 @@ class Flow1D : public Domain1D void show(const double* x) override; - void _getValues(const double* soln, const string& component, - vector& values) const override; - void _setValues(double* soln, const string& component, - const vector& values) override; - void _setProfile(double* soln, const string& component, - const vector& pos, const vector& values) override; - void _setFlatProfile(double* soln, const string& component, double v) override; + void getValues(const string& component, vector& values) const override; + void setValues(const string& component, const vector& values) override; + void setProfile(const string& component, + const vector& pos, const vector& values) override; + void setFlatProfile(const string& component, double v) override; shared_ptr asArray(const double* soln) const override; void fromArray(SolutionArray& arr, double* soln) override; diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index eec0e58be89..6669289c662 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -141,57 +141,6 @@ class OneDim : public SteadyStateSystem double upperBound(size_t i) const override; double lowerBound(size_t i) const override; - /** - * Retrieve component values from the solution vector. - * @param dom name of domain - * @param component name of component - * @param[out] values vector of values - * - * @since New in %Cantera 3.2. - */ - void getValues(const string& dom, const string& component, - vector& values) const; - - /** - * Specify component values in the solution vector. - * @param dom name of domain - * @param component name of component - * @param[in] values vector of values - * - * @since New in %Cantera 3.2. - */ - void setValues(const string& dom, const string& component, - const vector& values); - - /** - * Specify a profile for one component of one domain. - * @param dom name of domain - * @param component name of component - * @param[in] pos A vector of relative positions, beginning with 0.0 at the - * left of the domain, and ending with 1.0 at the right of the domain. - * @param[in] values A vector of values corresponding to the relative position - * locations. - * - * Note that the vector pos and values can have lengths different than the - * number of grid points, but their lengths must be equal. The values at - * the grid points will be linearly interpolated based on the (pos, - * values) specification. - * - * @since New in %Cantera 3.2. - */ - void setProfile(const string& dom, const string& component, - const vector& pos, const vector& values); - - /** - * Specify a flat profile for one component of one domain. - * @param dom name of domain - * @param component name of component - * @param value constant value - * - * @since New in %Cantera 3.2. - */ - void setFlatProfile(const string& dom, const string& component, double v); - /** * Initialize all domains. On the first call, this methods calls the init * method of each domain, proceeding from left to right. Subsequent calls diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index 28e79e2deb8..a7111863df9 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -95,39 +95,6 @@ bool Domain1D::hasComponent(const string& name, bool checkAlias) const "no component named '{}'", name); } -void Domain1D::getValues(const string& component, vector& values) const -{ - if (!m_container) { - throw CanteraError("Domain1D::getValues", "Domain is not connected."); - } - m_container->getValues(id(), component, values); -} - -void Domain1D::setValues(const string& component, const vector& values) -{ - if (!m_container) { - throw CanteraError("Domain1D::setValues", "Domain is not connected."); - } - m_container->setValues(id(), component, values); -} - -void Domain1D::setProfile(const string& component, - const vector& pos, const vector& values) -{ - if (!m_container) { - throw CanteraError("Domain1D::setProfile", "Domain is not connected."); - } - m_container->setProfile(id(), component, pos, values); -} - -void Domain1D::setFlatProfile(const string& component, double v) -{ - if (!m_container) { - throw CanteraError("Domain1D::setFlatProfile", "Domain is not connected."); - } - m_container->setFlatProfile(id(), component, v); -} - void Domain1D::setTransientTolerances(double rtol, double atol, size_t n) { if (n == npos) { diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 5df36c8630e..026e0ec625b 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -945,59 +945,72 @@ AnyMap Flow1D::getMeta() const return state; } -void Flow1D::_getValues(const double* soln, const string& component, - vector& values) const +void Flow1D::getValues(const string& component, vector& values) const { + if (!m_state) { + throw CanteraError("Flow1D::getValues", + "Domain needs to be installed in a container."); + } if (values.size() != nPoints()) { - throw ArraySizeError("Flow1D::_getValues", values.size(), nPoints()); + throw ArraySizeError("Flow1D::getValues", values.size(), nPoints()); } auto i = componentIndex(component); if (!componentActive(i)) { warn_user( - "Flow1D::_getValues", "Component '{}' is not used by '{}'.", + "Flow1D::getValues", "Component '{}' is not used by '{}'.", component, domainType()); } + const double* soln = m_state->data() + m_iloc; for (size_t j = 0; j < nPoints(); j++) { values[j] = soln[index(i,j)]; } } -void Flow1D::_setValues(double* soln, const string& component, - const vector& values) +void Flow1D::setValues(const string& component, const vector& values) { + if (!m_state) { + throw CanteraError("Flow1D::setValues", + "Domain needs to be installed in a container."); + } if (values.size() != nPoints()) { - throw ArraySizeError("Flow1D::_setValues", values.size(), nPoints()); + throw ArraySizeError("Flow1D::_etValues", values.size(), nPoints()); } auto i = componentIndex(component); if (!componentActive(i)) { throw CanteraError( - "Flow1D::_setValues", "Component '{}' is not used by '{}'.", + "Flow1D::setValues", "Component '{}' is not used by '{}'.", component, domainType()); } + double* soln = m_state->data() + m_iloc; for (size_t j = 0; j < nPoints(); j++) { soln[index(i,j)] = values[j]; } } -void Flow1D::_setProfile(double* soln, const string& component, - const vector& pos, const vector& values) +void Flow1D::setProfile(const string& component, + const vector& pos, const vector& values) { + if (!m_state) { + throw CanteraError("Flow1D::setProfile", + "Domain needs to be installed in a container."); + } if (pos.size() != values.size()) { throw CanteraError( - "Flow1D::_setProfile", "Vectors for positions and values must have same " + "Flow1D::setProfile", "Vectors for positions and values must have same " "size.\nSizes are {} and {}, respectively.", pos.size(), values.size()); } if (pos.front() != 0.0 || pos.back() != 1.0) { - throw CanteraError("Flow1D::_setProfile", + throw CanteraError("Flow1D::setProfile", "'pos' vector must span the range [0, 1]. Got a vector spanning " "[{}, {}] instead.", pos.front(), pos.back()); } auto i = componentIndex(component); if (!componentActive(i)) { throw CanteraError( - "Flow1D::_setProfile", "Component '{}' is not used by '{}'.", + "Flow1D::setProfile", "Component '{}' is not used by '{}'.", component, domainType()); } + double* soln = m_state->data() + m_iloc; double z0 = zmin(); double zDelta = zmax() - z0; for (size_t j = 0; j < nPoints(); j++) { @@ -1007,14 +1020,19 @@ void Flow1D::_setProfile(double* soln, const string& component, } } -void Flow1D::_setFlatProfile(double* soln, const string& component, double v) +void Flow1D::setFlatProfile(const string& component, double v) { + if (!m_state) { + throw CanteraError("Flow1D::setFlatProfile", + "Domain needs to be installed in a container."); + } auto i = componentIndex(component); if (!componentActive(i)) { throw CanteraError( - "Flow1D::_setFlatProfile", "Component '{}' is not used by '{}'.", + "Flow1D::setFlatProfile", "Component '{}' is not used by '{}'.", component, domainType()); } + double* soln = m_state->data() + m_iloc; for (size_t j = 0; j < nPoints(); j++) { soln[index(i,j)] = v; } diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 87f30228c32..2b8eec13440 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -72,34 +72,6 @@ double OneDim::lowerBound(size_t i) const return dom.lowerBound(k); } -void OneDim::getValues(const string& dom, const string& component, - vector& values) const -{ - size_t dIdx = domainIndex(dom); - domain(dIdx)._getValues(m_state->data() + domain(dIdx).loc(), component, values); -} - -void OneDim::setValues(const string& dom, const string& component, - const vector& values) -{ - size_t dIdx = domainIndex(dom); - domain(dIdx)._setValues(m_state->data() + domain(dIdx).loc(), component, values); -} - -void OneDim::setProfile(const string& dom, const string& component, - const vector& pos, const vector& values) -{ - size_t dIdx = domainIndex(dom); - domain(dIdx)._setProfile(m_state->data() + domain(dIdx).loc(), component, - pos, values); -} - -void OneDim::setFlatProfile(const string& dom, const string& component, double v) -{ - size_t dIdx = domainIndex(dom); - domain(dIdx)._setFlatProfile(m_state->data() + domain(dIdx).loc(), component, v); -} - void OneDim::addDomain(shared_ptr d) { // if 'd' is not the first domain, link it to the last domain From ed0b37e562a490d23ad8b30e8768e2fef7332aa2 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Mon, 6 Oct 2025 10:32:07 -0500 Subject: [PATCH 18/22] [oneD] Simplify to/fromArray --- include/cantera/oneD/Boundary1D.h | 24 +++++----- include/cantera/oneD/Domain1D.h | 25 ++++++++--- include/cantera/oneD/Flow1D.h | 5 ++- src/oneD/Boundary1D.cpp | 73 ++++++++++++++++++++----------- src/oneD/Domain1D.cpp | 24 ---------- src/oneD/Flow1D.cpp | 41 ++++++++++++----- src/oneD/Sim1D.cpp | 10 ++--- 7 files changed, 116 insertions(+), 86 deletions(-) diff --git a/include/cantera/oneD/Boundary1D.h b/include/cantera/oneD/Boundary1D.h index 64cc323986e..c4948e57ef6 100644 --- a/include/cantera/oneD/Boundary1D.h +++ b/include/cantera/oneD/Boundary1D.h @@ -109,7 +109,7 @@ class Boundary1D : public Domain1D void setupGrid(size_t n, const double* z) override {} - void fromArray(SolutionArray& arr, double* soln) override; + void fromArray(const shared_ptr& arr) override; protected: //! Initialize member variables based on the adjacent domains. @@ -173,8 +173,8 @@ class Inlet1D : public Boundary1D } void init() override; void eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) override; - shared_ptr asArray(const double* soln) const override; - void fromArray(SolutionArray& arr, double* soln) override; + shared_ptr toArray(bool normalize=false) const override; + void fromArray(const shared_ptr& arr) override; protected: //! A marker that indicates whether this is a left inlet or a right inlet. @@ -216,7 +216,7 @@ class Empty1D : public Boundary1D void eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) override; - shared_ptr asArray(const double* soln) const override; + shared_ptr toArray(bool normalize=false) const override; }; /** @@ -245,7 +245,7 @@ class Symm1D : public Boundary1D void eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) override; - shared_ptr asArray(const double* soln) const override; + shared_ptr toArray(bool normalize=false) const override; }; @@ -275,7 +275,7 @@ class Outlet1D : public Boundary1D void eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) override; - shared_ptr asArray(const double* soln) const override; + shared_ptr toArray(bool normalize=false) const override; }; @@ -312,8 +312,8 @@ class OutletRes1D : public Boundary1D void init() override; void eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) override; - shared_ptr asArray(const double* soln) const override; - void fromArray(SolutionArray& arr, double* soln) override; + shared_ptr toArray(bool normalize=false) const override; + void fromArray(const shared_ptr& arr) override; protected: size_t m_nsp = 0; //!< Number of species in the adjacent flow domain @@ -347,8 +347,8 @@ class Surf1D : public Boundary1D void init() override; void eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) override; - shared_ptr asArray(const double* soln) const override; - void fromArray(SolutionArray& arr, double* soln) override; + shared_ptr toArray(bool normalize=false) const override; + void fromArray(const shared_ptr& arr) override; void show(const double* x) override; }; @@ -390,8 +390,8 @@ class ReactingSurf1D : public Boundary1D void eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) override; - shared_ptr asArray(const double* soln) const override; - void fromArray(SolutionArray& arr, double* soln) override; + shared_ptr toArray(bool normalize=false) const override; + void fromArray(const shared_ptr& arr) override; void _getInitialSoln(double* x) override { m_sphase->getCoverages(x); diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index 50314186a9d..e5d3e50bf11 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -420,9 +420,12 @@ class Domain1D * directly in future revisions, where a non-const version will be implemented. * * @since New in %Cantera 3.0. + * @deprecated To be removed after %Cantera 3.2. Replaceable by toArray(). */ - virtual shared_ptr asArray(const double* soln) const { - throw NotImplementedError("Domain1D::asArray", "Needs to be overloaded."); + shared_ptr asArray(const double* soln) const { + warn_deprecated("Domain1D::asArray", + "To be removed after Cantera 3.2. Replaceable by 'toArray'."); + return toArray(); } //! Save the state of this domain to a SolutionArray. @@ -433,7 +436,9 @@ class Domain1D * * @since New in %Cantera 3.0. */ - shared_ptr toArray(bool normalize=false) const; + virtual shared_ptr toArray(bool normalize=false) const { + throw NotImplementedError("Domain1D::toArray", "Needs to be overloaded."); + } //! Restore the solution for this domain from a SolutionArray /*! @@ -441,9 +446,15 @@ class Domain1D * @param[out] soln Value of the solution vector, local to this domain * * @since New in %Cantera 3.0. + * @deprecated To be removed after %Cantera 3.2. + * Replaceable by version that does not require solution vector. */ - virtual void fromArray(SolutionArray& arr, double* soln) { - throw NotImplementedError("Domain1D::fromArray", "Needs to be overloaded."); + void fromArray(SolutionArray& arr, double* soln) { + warn_deprecated("Domain1D::fromArray", + "To be removed after Cantera 3.2. Replaceable by version that does not " + "require solution vector."); + shared_ptr arr_ptr(&arr, [](SolutionArray*){}); + fromArray(arr_ptr); } //! Restore the solution for this domain from a SolutionArray. @@ -452,7 +463,9 @@ class Domain1D * @param arr SolutionArray defining the state of this domain * @since New in %Cantera 3.0. */ - void fromArray(const shared_ptr& arr); + virtual void fromArray(const shared_ptr& arr) { + throw NotImplementedError("Domain1D::fromArray", "Needs to be overloaded."); + } //! Return thermo/kinetics/transport manager used in the domain //! @since New in %Cantera 3.0. diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index d49fb2f4807..cb8a1526300 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -7,6 +7,7 @@ #define CT_FLOW1D_H #include "Domain1D.h" +#include "OneDim.h" #include "cantera/base/Array.h" #include "cantera/base/Solution.h" #include "cantera/thermo/ThermoPhase.h" @@ -190,8 +191,8 @@ class Flow1D : public Domain1D const vector& pos, const vector& values) override; void setFlatProfile(const string& component, double v) override; - shared_ptr asArray(const double* soln) const override; - void fromArray(SolutionArray& arr, double* soln) override; + shared_ptr toArray(bool normalize=false) const override; + void fromArray(const shared_ptr& arr) override; //! Set flow configuration for freely-propagating flames, using an internal point //! with a fixed temperature as the condition to determine the inlet mass flux. diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp index 7d387559084..9b35cbb41aa 100644 --- a/src/oneD/Boundary1D.cpp +++ b/src/oneD/Boundary1D.cpp @@ -73,9 +73,9 @@ void Boundary1D::_init(size_t n) } } -void Boundary1D::fromArray(SolutionArray& arr, double* soln) +void Boundary1D::fromArray(const shared_ptr& arr) { - setMeta(arr.meta()); + setMeta(arr->meta()); } // ---------------- Inlet1D methods ---------------- @@ -264,7 +264,7 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg, } } -shared_ptr Inlet1D::asArray(const double* soln) const +shared_ptr Inlet1D::toArray(bool normalize) const { AnyMap meta = Boundary1D::getMeta(); meta["mass-flux"] = m_mdot; @@ -278,21 +278,24 @@ shared_ptr Inlet1D::asArray(const double* soln) const phase->saveState(data); arr->setState(0, data); + if (normalize) { + arr->normalize(); + } return arr; } -void Inlet1D::fromArray(SolutionArray& arr, double* soln) +void Inlet1D::fromArray(const shared_ptr& arr) { - Boundary1D::setMeta(arr.meta()); - arr.setLoc(0); - auto phase = arr.thermo(); - auto meta = arr.meta(); + Boundary1D::setMeta(arr->meta()); + arr->setLoc(0); + auto phase = arr->thermo(); + auto meta = arr->meta(); m_temp = phase->temperature(); if (meta.hasKey("mass-flux")) { m_mdot = meta.at("mass-flux").asDouble(); } else { // convert data format used by Python h5py export (Cantera < 3.0) - auto aux = arr.getAuxiliary(0); + auto aux = arr->getAuxiliary(0); m_mdot = phase->density() * aux.at("velocity").as(); } phase->getMassFractions(m_yin.data()); @@ -310,7 +313,7 @@ void Empty1D::eval(size_t jg, double* xg, double* rg, { } -shared_ptr Empty1D::asArray(const double* soln) const +shared_ptr Empty1D::toArray(bool normalize) const { AnyMap meta = Boundary1D::getMeta(); return SolutionArray::create(m_solution, 0, meta); @@ -362,7 +365,7 @@ void Symm1D::eval(size_t jg, double* xg, double* rg, integer* diagg, } } -shared_ptr Symm1D::asArray(const double* soln) const +shared_ptr Symm1D::toArray(bool normalize) const { AnyMap meta = Boundary1D::getMeta(); return SolutionArray::create(m_solution, 0, meta); @@ -429,7 +432,7 @@ void Outlet1D::eval(size_t jg, double* xg, double* rg, integer* diagg, } } -shared_ptr Outlet1D::asArray(const double* soln) const +shared_ptr Outlet1D::toArray(bool normalize) const { AnyMap meta = Boundary1D::getMeta(); return SolutionArray::create(m_solution, 0, meta); @@ -511,7 +514,7 @@ void OutletRes1D::eval(size_t jg, double* xg, double* rg, } } -shared_ptr OutletRes1D::asArray(const double* soln) const +shared_ptr OutletRes1D::toArray(bool normalize) const { AnyMap meta = Boundary1D::getMeta(); meta["temperature"] = m_temp; @@ -525,14 +528,17 @@ shared_ptr OutletRes1D::asArray(const double* soln) const phase->saveState(data); arr->setState(0, data); + if (normalize) { + arr->normalize(); + } return arr; } -void OutletRes1D::fromArray(SolutionArray& arr, double* soln) +void OutletRes1D::fromArray(const shared_ptr& arr) { - Boundary1D::setMeta(arr.meta()); - arr.setLoc(0); - auto phase = arr.thermo(); + Boundary1D::setMeta(arr->meta()); + arr->setLoc(0); + auto phase = arr->thermo(); m_temp = phase->temperature(); auto Y = phase->massFractions(); std::copy(Y, Y + m_nsp, &m_yres[0]); @@ -570,16 +576,16 @@ void Surf1D::eval(size_t jg, double* xg, double* rg, } } -shared_ptr Surf1D::asArray(const double* soln) const +shared_ptr Surf1D::toArray(bool normalize) const { AnyMap meta = Boundary1D::getMeta(); meta["temperature"] = m_temp; return SolutionArray::create(m_solution, 0, meta); } -void Surf1D::fromArray(SolutionArray& arr, double* soln) +void Surf1D::fromArray(const shared_ptr& arr) { - auto meta = arr.meta(); + auto meta = arr->meta(); m_temp = meta["temperature"].asDouble(); meta.erase("temperature"); Boundary1D::setMeta(meta); @@ -744,8 +750,13 @@ void ReactingSurf1D::eval(size_t jg, double* xg, double* rg, } } -shared_ptr ReactingSurf1D::asArray(const double* soln) const +shared_ptr ReactingSurf1D::toArray(bool normalize) const { + if (!m_state) { + throw CanteraError("ReactingSurf1D::toArray", + "Domain needs to be installed in a container before calling toArray."); + } + double* soln = m_state->data() + m_iloc; AnyMap meta = Boundary1D::getMeta(); meta["temperature"] = m_temp; meta["phase"]["name"] = m_sphase->name(); @@ -760,20 +771,32 @@ shared_ptr ReactingSurf1D::asArray(const double* soln) const auto arr = SolutionArray::create(m_solution, 1, meta); arr->setState(0, data); + if (normalize) { + arr->normalize(); + } return arr; } -void ReactingSurf1D::fromArray(SolutionArray& arr, double* soln) +void ReactingSurf1D::fromArray(const shared_ptr& arr) { - Boundary1D::setMeta(arr.meta()); - arr.setLoc(0); - auto surf = std::dynamic_pointer_cast(arr.thermo()); + if (!m_state) { + throw CanteraError("Domain1D::fromArray", + "Domain needs to be installed in a container before calling fromArray."); + } + resize(nComponents(), arr->size()); + m_container->resize(); + double* soln = m_state->data() + m_iloc; + + Boundary1D::setMeta(arr->meta()); + arr->setLoc(0); + auto surf = std::dynamic_pointer_cast(arr->thermo()); if (!surf) { throw CanteraError("ReactingSurf1D::fromArray", "Restoring of coverages requires surface phase"); } m_temp = surf->temperature(); surf->getCoverages(soln); + _finalize(soln); } void ReactingSurf1D::show(const double* x) diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index a7111863df9..9dee0a695ba 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -157,30 +157,6 @@ AnyMap Domain1D::getMeta() const return state; } -shared_ptr Domain1D::toArray(bool normalize) const { - if (!m_state) { - throw CanteraError("Domain1D::toArray", - "Domain needs to be installed in a container before calling toArray."); - } - auto ret = asArray(m_state->data() + m_iloc); - if (normalize) { - ret->normalize(); - } - return ret; -} - -void Domain1D::fromArray(const shared_ptr& arr) -{ - if (!m_state) { - throw CanteraError("Domain1D::fromArray", - "Domain needs to be installed in a container before calling fromArray."); - } - resize(nComponents(), arr->size()); - m_container->resize(); - fromArray(*arr, m_state->data() + m_iloc); - _finalize(m_state->data() + m_iloc); -} - void Domain1D::setMeta(const AnyMap& meta) { auto set_tols = [&](const AnyValue& tols, const string& which, vector& out) diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 026e0ec625b..e71dee5321e 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -1038,8 +1038,13 @@ void Flow1D::setFlatProfile(const string& component, double v) } } -shared_ptr Flow1D::asArray(const double* soln) const +shared_ptr Flow1D::toArray(bool normalize) const { + if (!m_state) { + throw CanteraError("Flow1D::toArray", + "Domain needs to be installed in a container before calling toArray."); + } + double* soln = m_state->data() + m_iloc; auto arr = SolutionArray::create( m_solution, static_cast(nPoints()), getMeta()); arr->addExtra("grid", false); // leading entry @@ -1069,33 +1074,44 @@ shared_ptr Flow1D::asArray(const double* soln) const arr->setComponent("radiative-heat-loss", value); } + if (normalize) { + arr->normalize(); + } return arr; } -void Flow1D::fromArray(SolutionArray& arr, double* soln) +void Flow1D::fromArray(const shared_ptr& arr) { - Domain1D::setMeta(arr.meta()); - arr.setLoc(0); - auto phase = arr.thermo(); + if (!m_state) { + throw CanteraError("Domain1D::fromArray", + "Domain needs to be installed in a container before calling fromArray."); + } + resize(nComponents(), arr->size()); + m_container->resize(); + double* soln = m_state->data() + m_iloc; + + Domain1D::setMeta(arr->meta()); + arr->setLoc(0); + auto phase = arr->thermo(); m_press = phase->pressure(); - const auto grid = arr.getComponent("grid").as>(); + const auto grid = arr->getComponent("grid").as>(); setupGrid(nPoints(), &grid[0]); - setMeta(arr.meta()); // can affect which components are active + setMeta(arr->meta()); // can affect which components are active for (size_t i = 0; i < nComponents(); i++) { if (!componentActive(i)) { continue; } string name = componentName(i); - if (arr.hasComponent(name)) { - const vector data = arr.getComponent(name).as>(); + if (arr->hasComponent(name)) { + const vector data = arr->getComponent(name).as>(); for (size_t j = 0; j < nPoints(); j++) { soln[index(i,j)] = data[j]; } - } else if (name == "Lambda" && arr.hasComponent("lambda")) { + } else if (name == "Lambda" && arr->hasComponent("lambda")) { // edge case: 'lambda' is renamed to 'Lambda' in Cantera 3.2 - const vector data = arr.getComponent("lambda").as>(); + const auto data = arr->getComponent("lambda").as>(); for (size_t j = 0; j < nPoints(); j++) { soln[index(i,j)] = data[j]; } @@ -1105,7 +1121,8 @@ void Flow1D::fromArray(SolutionArray& arr, double* soln) } } - updateProperties(npos, soln + loc(), 0, m_points - 1); + updateProperties(npos, soln, 0, m_points - 1); + _finalize(soln); } void Flow1D::setMeta(const AnyMap& state) diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index d2c20abdef4..5fe6728c049 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -102,7 +102,7 @@ void Sim1D::save(const string& fname, const string& name, const string& desc, string extension = (dot != npos) ? toLowerCopy(fname.substr(dot+1)) : ""; if (extension == "csv") { for (auto dom : m_dom) { - auto arr = dom->asArray(m_state->data() + dom->loc()); + auto arr = dom->toArray(); if (dom->size() > 1) { arr->writeEntry(fname, overwrite, basis); break; @@ -117,7 +117,7 @@ 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->asArray(m_state->data() + dom->loc()); + auto arr = dom->toArray(); arr->writeEntry(fname, name, dom->id(), overwrite, compression); } return; @@ -131,7 +131,7 @@ 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->asArray(m_state->data() + dom->loc()); + auto arr = dom->toArray(); arr->writeEntry(data, name, dom->id(), overwrite); } @@ -275,7 +275,7 @@ AnyMap Sim1D::restore(const string& fname, const string& name) m_xlast_ts.clear(); for (auto dom : m_dom) { try { - dom->fromArray(*arrs[dom->id()], m_state->data() + dom->loc()); + dom->fromArray(arrs[dom->id()]); } catch (CanteraError& err) { throw CanteraError("Sim1D::restore", "Encountered exception when restoring domain '{}' from HDF:\n{}", @@ -304,7 +304,7 @@ AnyMap Sim1D::restore(const string& fname, const string& name) m_xlast_ts.clear(); for (auto dom : m_dom) { try { - dom->fromArray(*arrs[dom->id()], m_state->data() + dom->loc()); + dom->fromArray(arrs[dom->id()]); } catch (CanteraError& err) { throw CanteraError("Sim1D::restore", "Encountered exception when restoring domain '{}' from YAML:\n{}", From 4075baa74e1e4a12499864af493291f47af61cc7 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Tue, 7 Oct 2025 20:43:32 -0500 Subject: [PATCH 19/22] Address review comments --- include/cantera/base/SolutionArray.h | 4 ++- include/cantera/oneD/Domain1D.h | 29 +++++++++++-------- include/cantera/oneD/Flow1D.h | 6 ++-- interfaces/cython/cantera/_onedim.pyx | 9 ++++-- interfaces/cython/cantera/onedim.py | 2 +- .../src/sourcegen/headers/ctdomain_auto.yaml | 2 +- src/base/SolutionArray.cpp | 14 +++++---- src/oneD/Flow1D.cpp | 16 ++++++---- 8 files changed, 50 insertions(+), 32 deletions(-) diff --git a/include/cantera/base/SolutionArray.h b/include/cantera/base/SolutionArray.h index 7db418b5dea..87e36a21134 100644 --- a/include/cantera/base/SolutionArray.h +++ b/include/cantera/base/SolutionArray.h @@ -129,6 +129,8 @@ class SolutionArray /** * Check whether SolutionArray contains a component. * A component is a property defining state or auxiliary variable. + * @param name name of component + * @param checkAlias if `true` (default), check alias mapping */ bool hasComponent(const string& name, bool checkAlias=true) const; @@ -405,7 +407,7 @@ class SolutionArray }; //! Return mapping of component alias names to standardized component names. -const std::map& _componentAliasMap(); +const map& _componentAliasMap(); } diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index e5d3e50bf11..05bd0620e7e 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -199,14 +199,18 @@ class Domain1D /** * Index of component with name `name`. * @param name name of component - * @param checkAlias if `True` (default), check alias mapping + * @param checkAlias if `true` (default), check alias mapping + * + * @since New in %Cantera 3.2. */ virtual size_t componentIndex(const string& name, bool checkAlias=true) const; /** - * Check whether the Domain contains a component. + * Check whether the Domain contains a component. * @param name name of component - * @param checkAlias if `True` (default), check alias mapping + * @param checkAlias if `true` (default), check alias mapping + * + * @since New in %Cantera 3.2. */ virtual bool hasComponent(const string& name, bool checkAlias=true) const; @@ -361,29 +365,29 @@ class Domain1D /** * Retrieve component values. - * @param component name of component - * @param[out] values vector of values + * @param component Name of the component. + * @param[out] values Vector of length nPoints() containing values at grid points. * * @since New in %Cantera 3.2. */ virtual void getValues(const string& component, vector& values) const { - throw NotImplementedError("Domain1D::setProfile", "Needs to be overloaded."); + throw NotImplementedError("Domain1D::getValues", "Needs to be overloaded."); } /** * Specify component values. - * @param component name of component - * @param[in] values vector of values + * @param component Name of the component. + * @param[in] values Vector of length nPoints() containing values at grid points. * * @since New in %Cantera 3.2. */ virtual void setValues(const string& component, const vector& values) { - throw NotImplementedError("Domain1D::setProfile", "Needs to be overloaded."); + throw NotImplementedError("Domain1D::setValues", "Needs to be overloaded."); } /** * Specify a profile for a component. - * @param component name of component + * @param component Name of the component. * @param[in] pos A vector of relative positions, beginning with 0.0 at the * left of the domain, and ending with 1.0 at the right of the domain. * @param[in] values A vector of values corresponding to the relative position @@ -403,8 +407,8 @@ class Domain1D /** * Specify a flat profile for a component. - * @param component name of component - * @param value constant value + * @param component Name of the component. + * @param value Constant value. * * @since New in %Cantera 3.2. */ @@ -453,6 +457,7 @@ class Domain1D warn_deprecated("Domain1D::fromArray", "To be removed after Cantera 3.2. Replaceable by version that does not " "require solution vector."); + // create a shared_ptr that skips deletion of the held pointer shared_ptr arr_ptr(&arr, [](SolutionArray*){}); fromArray(arr_ptr); } diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index cb8a1526300..060ade16d61 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -563,9 +563,9 @@ class Flow1D : public Domain1D * \frac{d\Lambda}{dz} = 0 * @f] * - * The Lambda equation serves as an eigenvalue that allows the momentum - * equation and continuity equations to be simultaneously satisfied in - * axisymmetric flows. The Lambda equation propagates information from + * The radial pressure gradient @f$ \Lambda @f$ serves as an eigenvalue that allows + * the momentum and continuity equations to be simultaneously satisfied in + * axisymmetric flows. This equation propagates information from * left-to-right. The default boundary condition is @f$ \Lambda = 0 @f$ * at the left boundary. The equation is first order and so only one * boundary condition is needed. diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 64e1efe06de..dc7ace20d44 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -532,6 +532,9 @@ cdef class FlowBase(Domain1D): self.flow.solveEnergyEqn() def __getattr__(self, name): + # used to access fields by alias rather than canonical names; + # also provides access to species that are valid Python identifiers + # (replicates SolutionArray behavior) component_name = name if (not self._has_component(name) and self._has_component(name.replace("_", "-"))): @@ -597,7 +600,7 @@ cdef class FlowBase(Domain1D): """ Array containing the electric field strength at each point. Note: This value is named ``eField`` in the C++ code and is only defined if - the transport model is ``ionized_gas``. + the transport model is ``ionized-gas``. .. versionadded:: 3.2 """ @@ -1092,8 +1095,8 @@ cdef class Sim1D: To be removed after Cantera 3.2. Replaceable by `Domain1D.get_values`. """ - warnings.warn("To be removed after Cantera 3.2. Replaceable by 'get_values'.", - DeprecationWarning) + warnings.warn("To be removed after Cantera 3.2. Replaceable by " + "'Domain1D.get_values'.", DeprecationWarning) idom, kcomp = self._get_indices(domain, component) dom = self.domains[idom] cdef int j diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index e83f4239ecb..17700b3eb29 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -305,7 +305,7 @@ def E(self): """ Array containing the electric field strength at each point. Note: This value is named ``eField`` in the C++ code and is only defined if - the transport model is ``ionized_gas``. + the transport model is ``ionized-gas``. """ return self.flame.get_values("eField") diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml index 1b1067ef9b5..238692c984d 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml @@ -19,7 +19,7 @@ recipes: - name: nPoints - name: componentName - name: componentIndex - wraps: componentIndex(const char*) + wraps: componentIndex(const string&) - name: getValues - name: setValues - name: setProfile diff --git a/src/base/SolutionArray.cpp b/src/base/SolutionArray.cpp index 23c50e94cec..02030e364ea 100644 --- a/src/base/SolutionArray.cpp +++ b/src/base/SolutionArray.cpp @@ -21,7 +21,12 @@ namespace ba = boost::algorithm; -const std::map aliasMap = { +namespace Cantera +{ + +namespace { // restrict scope of auxiliary variables to local translation unit + +const map aliasMap = { {"T", "temperature"}, {"P", "pressure"}, {"D", "density"}, @@ -35,7 +40,7 @@ const std::map aliasMap = { {"Q", "vapor-fraction"}, }; -const std::map reverseAliasMap = { +const map reverseAliasMap = { // reserved names used by states {"temperature", "T"}, {"pressure", "P"}, @@ -50,10 +55,9 @@ const std::map reverseAliasMap = { {"oxidizer-velocity", "Uo"}, }; -namespace Cantera -{ +} // end unnamed namespace -const std::map& _componentAliasMap() +const map& _componentAliasMap() { return reverseAliasMap; } diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index e71dee5321e..809d99a579a 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -16,7 +16,9 @@ using namespace std; namespace Cantera { -const std::map componentMap = { +namespace { // restrict scope of auxiliary variables to local translation unit + +const map componentMap = { {"velocity", c_offset_U}, // axial velocity [m/s] {"spread_rate", c_offset_V}, // strain rate {"T", c_offset_T}, // temperature [kelvin] @@ -25,6 +27,8 @@ const std::map componentMap = { {"Uo", c_offset_Uo}, // oxidizer axial velocity [m/s] }; +} // end unnamed namespace + Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) : Domain1D(nsp+c_offset_Y, points), m_nsp(nsp) @@ -847,13 +851,13 @@ size_t Flow1D::componentIndex(const string& name, bool checkAlias) const if (componentMap.count(name)) { return componentMap.at(name); } - for (size_t n=c_offset_Y; n Date: Tue, 7 Oct 2025 21:36:51 -0500 Subject: [PATCH 20/22] [oneD] Rename spread_rate to spreadRate This removes inconsistencies in naming convention, where the snake_case 'spread_rate' is an outlier. Note that the other component specific to axisymmetric flows is also renamed in the same PR. --- interfaces/cython/cantera/_onedim.pyx | 4 ++-- interfaces/cython/cantera/onedim.py | 16 ++++++++-------- samples/python/onedim/diffusion_flame_batch.py | 4 ++-- .../python/onedim/diffusion_flame_extinction.py | 2 +- src/base/SolutionArray.cpp | 3 ++- src/oneD/Flow1D.cpp | 6 +++--- test/python/test_onedim.py | 2 ++ .../cxx_samples/cxx_flamespeed_blessed.txt | 2 +- 8 files changed, 21 insertions(+), 18 deletions(-) diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index dc7ace20d44..3692a8b0ea8 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -577,8 +577,8 @@ cdef class FlowBase(Domain1D): def spread_rate(self): """ Array containing the tangential velocity gradient [1/s] (that is, radial - velocity divided by radius) at each point. Note: This value is only - defined for axisymmetric flows. + velocity divided by radius) at each point. Note: This value is named + ``spreadRate`` in the C++ code and is only defined for axisymmetric flows. .. versionadded:: 3.2 """ diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index 17700b3eb29..23bb13598e1 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -290,7 +290,7 @@ def spread_rate(self): Array containing the tangential velocity gradient [1/s] (that is, radial velocity divided by radius) at each point. """ - return self.flame.get_values("spread_rate") + return self.flame.get_values("spreadRate") @property def L(self): @@ -1051,7 +1051,7 @@ def set_initial_guess(self, data=None, group=None): zrel = (zz - zz[0])/dz self.flame.set_profile('velocity', [0.0, 1.0], [u0f, -u0o]) - self.flame.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) + self.flame.set_profile('spreadRate', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) self.flame.set_profile("Lambda", [0.0, 1.0], [L, L]) self.flame.set_profile('T', zrel, T) for k,spec in enumerate(self.gas.species_names): @@ -1354,8 +1354,8 @@ def set_initial_guess(self, products='inlet', data=None, group=None): self.flame.set_profile(self.gas.species_name(k), locs, [Y0[k], Y0[k]]) locs = np.array([0.0, 1.0]) - self.flame.set_profile('velocity', locs, [u0, 0.0]) - self.flame.set_profile('spread_rate', locs, [0.0, 0.0]) + self.flame.set_profile("velocity", locs, [u0, 0.0]) + self.flame.set_profile("spreadRate", locs, [0.0, 0.0]) class CounterflowPremixedFlame(FlameBase): @@ -1457,8 +1457,8 @@ def set_initial_guess(self, equilibrate=True, data=None, group=None): # estimate stagnation point x0 = rhou*uu * dz / (rhou*uu + rhob*ub) - self.flame.set_profile('velocity', [0.0, 1.0], [uu, -ub]) - self.flame.set_profile('spread_rate', [0.0, x0/dz, 1.0], [0.0, a, 0.0]) + self.flame.set_profile("velocity", [0.0, 1.0], [uu, -ub]) + self.flame.set_profile("spreadRate", [0.0, x0/dz, 1.0], [0.0, a, 0.0]) self.flame.set_profile("Lambda", [0.0, 1.0], [L, L]) @@ -1541,6 +1541,6 @@ def set_initial_guess(self, data=None, group=None): a = 2 * uu / dz L = - rhou * a**2 - self.flame.set_profile('velocity', [0.0, 1.0], [uu, 0]) - self.flame.set_profile('spread_rate', [0.0, 1.0], [0.0, a]) + self.flame.set_profile("velocity", [0.0, 1.0], [uu, 0]) + self.flame.set_profile("spreadRate", [0.0, 1.0], [0.0, a]) self.flame.set_profile("Lambda", [0.0, 1.0], [L, L]) diff --git a/samples/python/onedim/diffusion_flame_batch.py b/samples/python/onedim/diffusion_flame_batch.py index 75dfc55931f..986dcd6e2a5 100644 --- a/samples/python/onedim/diffusion_flame_batch.py +++ b/samples/python/onedim/diffusion_flame_batch.py @@ -137,7 +137,7 @@ def names(test): # Update velocities f.flame.set_values("velocity", f.flame.velocity * rel_pressure_increase ** exp_u_p) f.flame.set_values( - "spread_rate", f.flame.spread_rate * rel_pressure_increase ** exp_V_p) + "spreadRate", f.flame.spread_rate * rel_pressure_increase ** exp_V_p) # Update pressure curvature f.flame.set_values( "Lambda", f.flame.radial_pressure_gradient * rel_pressure_increase ** exp_lam_p) @@ -189,7 +189,7 @@ def names(test): f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a # Update velocities f.flame.set_values("velocity", f.flame.velocity * strain_factor ** exp_u_a) - f.flame.set_values("spread_rate", f.flame.spread_rate * strain_factor ** exp_V_a) + f.flame.set_values("spreadRate", f.flame.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature f.flame.set_values( "Lambda", f.flame.radial_pressure_gradient * strain_factor ** exp_lam_a) diff --git a/samples/python/onedim/diffusion_flame_extinction.py b/samples/python/onedim/diffusion_flame_extinction.py index 8e6df42b722..ac377dee306 100644 --- a/samples/python/onedim/diffusion_flame_extinction.py +++ b/samples/python/onedim/diffusion_flame_extinction.py @@ -129,7 +129,7 @@ def names(test): f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a # Update velocities f.flame.set_values("velocity", f.flame.velocity * strain_factor ** exp_u_a) - f.flame.set_values('spread_rate', f.flame.spread_rate * strain_factor ** exp_V_a) + f.flame.set_values("spreadRate", f.flame.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature f.flame.set_values( "Lambda", f.flame.radial_pressure_gradient * strain_factor ** exp_lam_a) diff --git a/src/base/SolutionArray.cpp b/src/base/SolutionArray.cpp index 02030e364ea..5b88dd716c3 100644 --- a/src/base/SolutionArray.cpp +++ b/src/base/SolutionArray.cpp @@ -46,7 +46,8 @@ const map reverseAliasMap = { {"pressure", "P"}, {"density", "D"}, // reserved names used for 1-D objects - {"spread-rate", "spread_rate"}, + {"spread-rate", "spreadRate"}, + {"spread_rate", "spreadRate"}, // snake_case version used prior to Cantera 3.2 {"L", "Lambda"}, {"radial-pressure-gradient", "Lambda"}, {"lambda", "Lambda"}, // lower-case version used prior to Cantera 3.2 diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 809d99a579a..97301e7cf0c 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -20,7 +20,7 @@ namespace { // restrict scope of auxiliary variables to local translation unit const map componentMap = { {"velocity", c_offset_U}, // axial velocity [m/s] - {"spread_rate", c_offset_V}, // strain rate + {"spreadRate", c_offset_V}, // strain rate {"T", c_offset_T}, // temperature [kelvin] {"Lambda", c_offset_L}, // (1/r)dP/dr {"eField", c_offset_E}, // electric field @@ -828,7 +828,7 @@ string Flow1D::componentName(size_t n) const case c_offset_U: return "velocity"; case c_offset_V: - return "spread_rate"; + return "spreadRate"; case c_offset_T: return "T"; case c_offset_L: @@ -885,7 +885,7 @@ bool Flow1D::hasComponent(const string& name, bool checkAlias) const bool Flow1D::componentActive(size_t n) const { switch (n) { - case c_offset_V: // spread_rate + case c_offset_V: // spreadRate return m_usesLambda; case c_offset_L: // Lambda return m_usesLambda; diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index 211ad9001bd..d933eccee53 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -257,6 +257,8 @@ def test_fixed_temperature(self): zeros = self.sim.T * 0 with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): self.sim.flame.set_values("spread_rate", zeros) + with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): + self.sim.flame.set_values("spreadRate", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): self.sim.flame.set_values("radial-pressure-gradient", zeros) with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): diff --git a/test_problems/cxx_samples/cxx_flamespeed_blessed.txt b/test_problems/cxx_samples/cxx_flamespeed_blessed.txt index 83d2b312e03..2940f5cf542 100644 --- a/test_problems/cxx_samples/cxx_flamespeed_blessed.txt +++ b/test_problems/cxx_samples/cxx_flamespeed_blessed.txt @@ -17,7 +17,7 @@ phi = 0.9, Tad = 2133.7925650475245 Pressure: 1.013e+05 Pa ------------------------------------------------------------------------------- - z velocity spread_rate T Lambda eField + z velocity spreadRate T Lambda eField ------------------------------------------------------------------------------- 0 0.3 0 300 0 0 0.02 0.3 0 300 0 0 From b0cc6b2cf3e07de13c06c3356dbbef96e23b6cab Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Wed, 8 Oct 2025 06:12:17 -0500 Subject: [PATCH 21/22] [oneD] Deprecate checkers used by legacy CLib --- include/cantera/oneD/Domain1D.h | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index 05bd0620e7e..fb78631fd6a 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -151,7 +151,10 @@ class Domain1D //! Check that the specified component index is in range. //! Throws an exception if n is greater than nComponents()-1 + //! @deprecated To be removed after %Cantera 3.2. Only used by legacy CLib. void checkComponentIndex(size_t n) const { + warn_deprecated("Domain1D::checkPointIndex", + "To be removed after Cantera 3.2. Only used by legacy CLib."); if (n >= m_nv) { throw IndexError("Domain1D::checkComponentIndex", "points", n, m_nv); } @@ -160,7 +163,10 @@ class Domain1D //! Check that an array size is at least nComponents(). //! Throws an exception if nn is less than nComponents(). Used before calls //! which take an array pointer. + //! @deprecated To be removed after %Cantera 3.2. Unused. void checkComponentArraySize(size_t nn) const { + warn_deprecated("Domain1D::checkComponentArraySize", + "To be removed after Cantera 3.2. Unused."); if (m_nv > nn) { throw ArraySizeError("Domain1D::checkComponentArraySize", nn, m_nv); } @@ -173,7 +179,10 @@ class Domain1D //! Check that the specified point index is in range. //! Throws an exception if n is greater than nPoints()-1 + //! @deprecated To be removed after %Cantera 3.2. Only used by legacy CLib. void checkPointIndex(size_t n) const { + warn_deprecated("Domain1D::checkPointIndex", + "To be removed after Cantera 3.2. Only used by legacy CLib."); if (n >= m_points) { throw IndexError("Domain1D::checkPointIndex", "points", n, m_points); } @@ -182,7 +191,10 @@ class Domain1D //! Check that an array size is at least nPoints(). //! Throws an exception if nn is less than nPoints(). Used before calls //! which take an array pointer. + //! @deprecated To be removed after %Cantera 3.2. Unused. void checkPointArraySize(size_t nn) const { + warn_deprecated("Domain1D::checkPointArraySize", + "To be removed after Cantera 3.2. Unused."); if (m_points > nn) { throw ArraySizeError("Domain1D::checkPointArraySize", nn, m_points); } From 026cb536b9e3760fb1aaa2827dafc5d697822d35 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Wed, 8 Oct 2025 05:50:24 -0500 Subject: [PATCH 22/22] [oneD] Implement single value access for Domain1D --- include/cantera/oneD/Boundary1D.h | 2 + include/cantera/oneD/Domain1D.h | 64 ++++++++++++++-- include/cantera/oneD/Flow1D.h | 3 + include/cantera/oneD/OneDim.h | 5 ++ include/cantera/oneD/Sim1D.h | 55 +++++++++++++- interfaces/cython/cantera/_onedim.pxd | 3 + interfaces/cython/cantera/_onedim.pyx | 73 +++++++++++++++++-- interfaces/cython/cantera/onedim.py | 14 ++-- .../src/sourcegen/headers/ctdomain_auto.yaml | 5 ++ .../src/sourcegen/headers/ctonedim_auto.yaml | 6 +- samples/cxx/flamespeed/flamespeed.cpp | 20 ++--- src/oneD/Boundary1D.cpp | 19 +++++ src/oneD/Domain1D.cpp | 6 +- src/oneD/Flow1D.cpp | 53 ++++++++++++++ src/oneD/Sim1D.cpp | 35 ++++----- test/clib/test_ctonedim.cpp | 4 +- test/oneD/test_oneD.cpp | 5 +- test/python/test_onedim.py | 6 +- 18 files changed, 311 insertions(+), 67 deletions(-) diff --git a/include/cantera/oneD/Boundary1D.h b/include/cantera/oneD/Boundary1D.h index c4948e57ef6..0f59ff29a10 100644 --- a/include/cantera/oneD/Boundary1D.h +++ b/include/cantera/oneD/Boundary1D.h @@ -384,12 +384,14 @@ class ReactingSurf1D : public Boundary1D } string componentName(size_t n) const override; + size_t componentIndex(const string& name, bool checkAlias=true) const override; void init() override; void resetBadValues(double* xg) override; void eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) override; + double value(const string& component, size_t localPoint=0) const override; shared_ptr toArray(bool normalize=false) const override; void fromArray(const shared_ptr& arr) override; diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index fb78631fd6a..d80a9f2cdce 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -370,9 +370,43 @@ class Domain1D * @param x solution vector * @param n component index * @param j grid point index + * + * @deprecated To be removed after %Cantera 3.2. Replaceable with version accessing + * component by name. */ double value(const double* x, size_t n, size_t j) const { - return x[index(n,j)]; + warn_deprecated("Domain1D::value", + "Replaceable with version accessing component by name"); + return x[index(n,j)]; // caution: this assumes m_iloc = 0 + } + + /** + * Set a single component value in a flow domain or at a boundary. + * Use getValues() for efficient access across an entire Flow1D domain. + * @param component Name of the component. + * @param localPoint Grid point within a Flow1D domain, beginning with 0 for + * the leftmost grid point. Unused for Boundary1D domains; defaults to zero. + * + * @since New in %Cantera 3.2. + */ + virtual double value(const string& component, size_t localPoint=0) const { + throw NotImplementedError("Domain1D::value", + "Not implemented for domain type '{}'.", domainType()); + } + + /** + * Set a single component value in a flow domain or at a boundary. + * Use setValues() for efficient access across an entire Flow1D domain. + * @param component Name of the component. + * @param value Value of the component. + * @param localPoint Grid point within a Flow1D domain, beginning with 0 for + * the leftmost grid point. Unused for Boundary1D domains; defaults to zero. + * + * @since New in %Cantera 3.2. + */ + virtual void setValue(const string& component, double value, size_t localPoint=0) { + throw NotImplementedError("Domain1D::setValue", + "Not implemented for domain type '{}'.", domainType()); } /** @@ -383,7 +417,8 @@ class Domain1D * @since New in %Cantera 3.2. */ virtual void getValues(const string& component, vector& values) const { - throw NotImplementedError("Domain1D::getValues", "Needs to be overloaded."); + throw NotImplementedError("Domain1D::getValues", + "Not implemented for domain type '{}'.", domainType()); } /** @@ -394,7 +429,23 @@ class Domain1D * @since New in %Cantera 3.2. */ virtual void setValues(const string& component, const vector& values) { - throw NotImplementedError("Domain1D::setValues", "Needs to be overloaded."); + throw NotImplementedError("Domain1D::setValues", + "Not implemented for domain type '{}'.", domainType()); + } + + /** + * Retrieve internal work array values for a component. + * After calling Sim1D::eval(), this array contains the values of the residual + * function. + * @param component Name of the component. + * @param[out] values Vector of length nPoints() containing residuals at grid + * points. + * + * @since New in %Cantera 3.2. + */ + virtual void getResiduals(const string& component, vector& values) const { + throw NotImplementedError("Domain1D::getResiduals", + "Not applicable or not implemented for domain type '{}'.", domainType()); } /** @@ -414,7 +465,8 @@ class Domain1D */ virtual void setProfile(const string& component, const vector& pos, const vector& values) { - throw NotImplementedError("Domain1D::setProfile", "Needs to be overloaded."); + throw NotImplementedError("Domain1D::setProfile", + "Not implemented for domain type '{}'.", domainType()); } /** @@ -425,8 +477,8 @@ class Domain1D * @since New in %Cantera 3.2. */ virtual void setFlatProfile(const string& component, double v) { - throw NotImplementedError( - "Domain1D::setFlatProfile", "Needs to be overloaded."); + throw NotImplementedError("Domain1D::setFlatProfile", + "Not implemented for domain type '{}'.", domainType()); } //! Save the state of this domain as a SolutionArray. diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index 060ade16d61..b6876dc1e3a 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -185,8 +185,11 @@ class Flow1D : public Domain1D void show(const double* x) override; + double value(const string& component, size_t localPoint=0) const override; + void setValue(const string& component, double value, size_t localPoint=0) override; void getValues(const string& component, vector& values) const override; void setValues(const string& component, const vector& values) override; + void getResiduals(const string& component, vector& values) const override; void setProfile(const string& component, const vector& pos, const vector& values) override; void setFlatProfile(const string& component, double v) override; diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index 6669289c662..99ef28cfa7c 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -249,6 +249,11 @@ class OneDim : public SteadyStateSystem return m_timeSteps; } + //! Access internal work array. + const vector& _workVector() const { + return m_xnew; + } + protected: //! All domains comprising the system vector> m_dom; diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index cacd65e255f..8ffd5f44f20 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/include/cantera/oneD/Sim1D.h @@ -63,8 +63,13 @@ class Sim1D : public OneDim * @param localPoint grid point within the domain, beginning with 0 for * the leftmost grid point in the domain. * @param value the value. + * @deprecated To be removed after %Cantera 3.2. Replaceable by Domain1D::setValue() */ - void setValue(size_t dom, size_t comp, size_t localPoint, double value); + void setValue(size_t dom, size_t comp, size_t localPoint, double value) { + warn_deprecated("Sim1D::setValue", + "To be removed after Cantera 3.2. Access from Domain1D object instead."); + _setValue(dom, comp, localPoint, value); + } /** * Get one entry in the solution vector. @@ -72,16 +77,60 @@ class Sim1D : public OneDim * @param comp component number * @param localPoint grid point within the domain, beginning with 0 for * the leftmost grid point in the domain. + * @deprecated To be removed after %Cantera 3.2. Replaceable by Domain1D::value() */ - double value(size_t dom, size_t comp, size_t localPoint) const; + double value(size_t dom, size_t comp, size_t localPoint) const { + warn_deprecated("Sim1D::setValue", + "To be removed after Cantera 3.2. Access from Domain1D object instead."); + return _value(dom, comp, localPoint); + } //! Get an entry in the work vector, which may contain either a new system state //! or the current residual of the system. //! @param dom domain index //! @param comp component index //! @param localPoint grid point within the domain - double workValue(size_t dom, size_t comp, size_t localPoint) const; + //! @deprecated To be removed after %Cantera 3.2. Replaceable by + //! Domain1D::getResiduals() + double workValue(size_t dom, size_t comp, size_t localPoint) const { + warn_deprecated("Sim1D::setValue", + "To be removed after Cantera 3.2. Access residuals via Domain1D instead."); + return _workValue(dom, comp, localPoint); + } + +protected: + /** + * Set a single value in the solution vector. + * @param dom domain number, beginning with 0 for the leftmost domain. + * @param comp component number + * @param localPoint grid point within the domain, beginning with 0 for + * the leftmost grid point in the domain. + * @param value the value. + * @since New in %Cantera 3.2. Previously part of public interface. + */ + void _setValue(size_t dom, size_t comp, size_t localPoint, double value); + /** + * Get one entry in the solution vector. + * @param dom domain number, beginning with 0 for the leftmost domain. + * @param comp component number + * @param localPoint grid point within the domain, beginning with 0 for + * the leftmost grid point in the domain. + * @since New in %Cantera 3.2. Previously part of public interface. + */ + double _value(size_t dom, size_t comp, size_t localPoint) const; + + /** + * Get an entry in the work vector, which may contain either a new system state + * or the current residual of the system. + * @param dom domain index + * @param comp component index + * @param localPoint grid point within the domain + * @since New in %Cantera 3.2. Previously part of public interface. + */ + double _workValue(size_t dom, size_t comp, size_t localPoint) const; + +public: /** * Specify a profile for one component of one domain. * @param dom domain number, beginning with 0 for the leftmost domain. diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index b412ba22295..30b6cc9c0e0 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -27,8 +27,11 @@ cdef extern from "cantera/oneD/Domain1D.h": cbool hasComponent(string&) except +translate_exception vector[double] grid() except +translate_exception void setupGrid(vector[double]&) except +translate_exception + double value(string&, size_t) except +translate_exception + void setValue(string&, double, size_t) except +translate_exception void getValues(string&, vector[double]&) except +translate_exception void setValues(string&, vector[double]&) except +translate_exception + void getResiduals(string&, vector[double]&) except +translate_exception void setProfile(string&, vector[double]&, vector[double]&) except +translate_exception void setFlatProfile(string&, double) except +translate_exception void setBounds(size_t, double, double) diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 3692a8b0ea8..56cfa7e457c 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -103,6 +103,38 @@ cdef class Domain1D: grid_vec.push_back(g) self.domain.setupGrid(grid_vec) + def value(self, str component, point=0): + """ + Component value at one point. + + :param component: + component name + :param point: + grid point number within ``domain`` starting with 0 on the left + + >>> t = d.value("T", 6) + + .. versionadded:: 3.2 + """ + return self.domain.value(stringify(component), point) + + def set_value(self, str component, value, point=0): + """ + Set the value of one component at one point to 'value'. + + :param component: + component name + :param value: + numerical value + :param point: + grid point number within ``domain`` starting with 0 on the left + + >>> d.set("T", 500., 5) + + .. versionadded:: 3.2 + """ + return self.domain.setValue(stringify(component), value, point) + def get_values(self, str component): """ Retrieve spatial profile of a component. @@ -110,7 +142,7 @@ cdef class Domain1D: :param component: component name - >>> T = d.get_values('T') + >>> T = d.get_values("T") .. versionadded:: 3.2 """ @@ -128,7 +160,7 @@ cdef class Domain1D: :param values: array containing values - >>> d.set_values('T', T) + >>> d.set_values("T", T) .. versionadded:: 3.2 """ @@ -137,6 +169,23 @@ cdef class Domain1D: values_vec.push_back(v) self.domain.setValues(stringify(component), values) + def get_residuals(self, str component): + """ + Retrieve internal work array value at one point. After calling `Sim1D.eval`, + this array contains the values of the residual function. + + :param component: + component name + + >>> T = d.get_residuals("T") + + .. versionadded:: 3.2 + """ + cdef vector[double] values + values.resize(self.n_points) + self.domain.getResiduals(stringify(component), values) + return np.asarray(values) + def set_profile(self, component, positions, values): """ Set an initial estimate for a profile of one component in one domain. @@ -1028,6 +1077,11 @@ cdef class Sim1D: grid point number within ``domain`` starting with 0 on the left >>> t = s.value('flow', 'T', 6) + + .. deprecated:: 3.2 + + To be removed after Cantera 3.2. Replaceable by `Domain1D.value` or + `Domain1D.getValues`. """ dom, comp = self._get_indices(domain, component) return self.sim.value(dom, comp, point) @@ -1048,6 +1102,11 @@ cdef class Sim1D: >>> s.set(d, 3, 5, 6.7) >>> s.set(1, 0, 5, 6.7) >>> s.set('flow', 'T', 5, 500) + + .. deprecated:: 3.2 + + To be removed after Cantera 3.2. Replaceable by `Domain1D.setValue` or + `Domain1D.setValues`. """ dom, comp = self._get_indices(domain, component) self.sim.setValue(dom, comp, point, value) @@ -1056,7 +1115,7 @@ cdef class Sim1D: """ Evaluate the governing equations using the current solution estimate, storing the residual in the array which is accessible with the - `work_value` function. + `Domain1D.get_residuals` function. :param rdt: Reciprocal of the time-step @@ -1076,6 +1135,10 @@ cdef class Sim1D: grid point number in the domain, starting with zero at the left >>> t = s.value(flow, 'T', 6) + + .. deprecated:: 3.2 + + To be removed after Cantera 3.2. Replaceable by `Domain1D.get_residuals`. """ dom, comp = self._get_indices(domain, component) return self.sim.workValue(dom, comp, point) @@ -1122,7 +1185,7 @@ cdef class Sim1D: .. deprecated:: 3.2 - To be removed after Cantera 3.2. Replaceable by Domain1D.set_profile. + To be removed after Cantera 3.2. Replaceable by `Domain1D.set_profile`. """ dom, comp = self._get_indices(domain, component) @@ -1148,7 +1211,7 @@ cdef class Sim1D: .. deprecated:: 3.2 - To be removed after Cantera 3.2. Replaceable by Domain1D.set_flat_profile. + To be removed after Cantera 3.2. Replaceable by `Domain1D.set_flat_profile`. """ dom, comp = self._get_indices(domain, component) self.sim.setFlatProfile(dom, comp, value) diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index 23bb13598e1..5a34c9dc4d7 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -399,11 +399,9 @@ def set_gas_state(self, point): Set the state of the `Solution` object used for calculations to the temperature and composition at the point with index ``point``. """ - k0 = self.flame.component_index(self.gas.species_name(0)) - Y = [self.value(self.flame, k, point) - for k in range(k0, k0 + self.gas.n_species)] + Y = [self.flame.value(k, point) for k in self.gas.species_names] self.gas.set_unnormalized_mass_fractions(Y) - self.gas.TP = self.value(self.flame, 'T', point), self.P + self.gas.TP = self.flame.value("T", point), self.P def to_array(self, domain=None, normalize=False): """ @@ -1246,8 +1244,8 @@ def mixture_fraction(self, m): >>> f.mixture_fraction('Bilger') """ - Yf = [self.value(self.flame, k, 0) for k in self.gas.species_names] - Yo = [self.value(self.flame, k, self.flame.n_points - 1) + Yf = [self.flame.value(k, 0) for k in self.gas.species_names] + Yo = [self.flame.value(k, self.flame.n_points - 1) for k in self.gas.species_names] vals = np.empty(self.flame.n_points) @@ -1258,8 +1256,8 @@ def mixture_fraction(self, m): @property def equivalence_ratio(self): - Yf = [self.value(self.flame, k, 0) for k in self.gas.species_names] - Yo = [self.value(self.flame, k, self.flame.n_points - 1) + Yf = [self.flame.value(k, 0) for k in self.gas.species_names] + Yo = [self.flame.value(k, self.flame.n_points - 1) for k in self.gas.species_names] vals = np.empty(self.flame.n_points) diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml index 238692c984d..466fbe0cd11 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml @@ -20,8 +20,13 @@ recipes: - name: componentName - name: componentIndex wraps: componentIndex(const string&) +- name: value + # alternative is deprecated in Cantera 3.2 + wraps: value(const string&, size_t) +- name: setValue - name: getValues - name: setValues +- name: getResiduals - name: setProfile # alternative is deprecated in Cantera 3.2 wraps: setProfile(const string&, const vector&, const vector&) diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml index b81fdd62942..97f97575490 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctonedim_auto.yaml @@ -9,7 +9,7 @@ parents: [OneDim, SteadyStateSystem] # List of parent classes recipes: - name: newSim1D uses: domain -- name: setValue +- name: setValue # deprecated - name: setProfile # deprecated - name: setFlatProfile # deprecated - name: setInitialGuess # deprecated @@ -26,8 +26,8 @@ recipes: wraps: _restore # skips AnyMap return (not implemented in CLib) - name: writeStats - name: domainIndex -- name: value -- name: workValue +- name: value # deprecated +- name: workValue # deprecated - name: eval wraps: eval(double*, double*, double, int) - name: setMaxJacAge diff --git a/samples/cxx/flamespeed/flamespeed.cpp b/samples/cxx/flamespeed/flamespeed.cpp index 39872894d2f..3d6e1ffc29e 100644 --- a/samples/cxx/flamespeed/flamespeed.cpp +++ b/samples/cxx/flamespeed/flamespeed.cpp @@ -139,24 +139,21 @@ int flamespeed(double phi, bool refine_grid, int loglevel) flow->solveEnergyEqn(); flame.solve(loglevel,refine_grid); - double flameSpeed_mix = flame.value(flowdomain, - flow->componentIndex("velocity"),0); + double flameSpeed_mix = flow->value("velocity", 0); print("Flame speed with mixture-averaged transport: {} m/s\n", flameSpeed_mix); flame.save(fileName, "mix", "Solution with mixture-averaged transport", true); // now switch to multicomponent transport flow->setTransportModel("multicomponent"); flame.solve(loglevel, refine_grid); - double flameSpeed_multi = flame.value(flowdomain, - flow->componentIndex("velocity"),0); + double flameSpeed_multi = flow->value("velocity", 0); print("Flame speed with multicomponent transport: {} m/s\n", flameSpeed_multi); flame.save(fileName, "multi", "Solution with multicomponent transport", true); // now enable Soret diffusion flow->enableSoret(true); flame.solve(loglevel, refine_grid); - double flameSpeed_full = flame.value(flowdomain, - flow->componentIndex("velocity"),0); + double flameSpeed_full = flow->value("velocity"); print("Flame speed with multicomponent transport + Soret: {} m/s\n", flameSpeed_full); flame.save(fileName, "soret", @@ -167,13 +164,10 @@ int flamespeed(double phi, bool refine_grid, int loglevel) print("\n{:9s}\t{:8s}\t{:5s}\t{:7s}\n", "z (m)", "T (K)", "U (m/s)", "Y(CO)"); for (size_t n = 0; n < flow->nPoints(); n++) { - Tvec.push_back(flame.value(flowdomain,flow->componentIndex("T"),n)); - COvec.push_back(flame.value(flowdomain, - flow->componentIndex("CO"),n)); - CO2vec.push_back(flame.value(flowdomain, - flow->componentIndex("CO2"),n)); - Uvec.push_back(flame.value(flowdomain, - flow->componentIndex("velocity"),n)); + Tvec.push_back(flow->value("T", n)); + COvec.push_back(flow->value("CO", n)); + CO2vec.push_back(flow->value("CO2", n)); + Uvec.push_back(flow->value("velocity", n)); zvec.push_back(flow->z(n)); print("{:9.6f}\t{:8.3f}\t{:5.3f}\t{:7.5f}\n", flow->z(n), Tvec[n], Uvec[n], COvec[n]); diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp index 9b35cbb41aa..94c752b32e5 100644 --- a/src/oneD/Boundary1D.cpp +++ b/src/oneD/Boundary1D.cpp @@ -647,6 +647,11 @@ string ReactingSurf1D::componentName(size_t n) const } } +size_t ReactingSurf1D::componentIndex(const string& name, bool checkAlias) const +{ + return m_sphase->speciesIndex(name); +} + void ReactingSurf1D::init() { m_nv = m_nsp; @@ -750,6 +755,20 @@ void ReactingSurf1D::eval(size_t jg, double* xg, double* rg, } } +double ReactingSurf1D::value(const string& component, size_t localPoint) const +{ + if (!m_state) { + throw CanteraError("ReactingSurf1D::value", + "Domain needs to be installed in a container."); + } + if (localPoint) { + throw IndexError("ReactingSurf1D::value", "", localPoint, 1); + } + auto i = componentIndex(component); + const double* soln = m_state->data() + m_iloc; + return soln[index(i, 0)]; +} + shared_ptr ReactingSurf1D::toArray(bool normalize) const { if (!m_state) { diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index 9dee0a695ba..25b52cdaf77 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -247,8 +247,7 @@ void Domain1D::show(const double* x) for (size_t j = 0; j < m_points; j++) { writelog("\n {:10.4g} ", m_z[j]); for (size_t n = 0; n < 5; n++) { - double v = value(x, i*5+n, j); - writelog(" {:10.4g} ", v); + writelog(" {:10.4g} ", x[index(i*5 + n, j)]); } } writelog("\n"); @@ -263,8 +262,7 @@ void Domain1D::show(const double* x) for (size_t j = 0; j < m_points; j++) { writelog("\n {:10.4g} ", m_z[j]); for (size_t n = 0; n < nrem; n++) { - double v = value(x, nn*5+n, j); - writelog(" {:10.4g} ", v); + writelog(" {:10.4g} ", x[index(nn*5 + n, j)]); } } writelog("\n"); diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 97301e7cf0c..0451197bdd5 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -949,6 +949,38 @@ AnyMap Flow1D::getMeta() const return state; } +double Flow1D::value(const string& component, size_t localPoint) const +{ + if (!m_state) { + throw CanteraError("Flow1D::value", + "Domain needs to be installed in a container."); + } + auto i = componentIndex(component); + if (!componentActive(i)) { + warn_user( + "Flow1D::values", "Component '{}' is not used by '{}'.", + component, domainType()); + } + const double* soln = m_state->data() + m_iloc; + return soln[index(i, localPoint)]; +} + +void Flow1D::setValue(const string& component, double value, size_t localPoint) +{ + if (!m_state) { + throw CanteraError("Flow1D::setValue", + "Domain needs to be installed in a container."); + } + auto i = componentIndex(component); + if (!componentActive(i)) { + throw CanteraError( + "Flow1D::setValue", "Component '{}' is not used by '{}'.", + component, domainType()); + } + double* soln = m_state->data() + m_iloc; + soln[index(i, localPoint)] = value; +} + void Flow1D::getValues(const string& component, vector& values) const { if (!m_state) { @@ -991,6 +1023,27 @@ void Flow1D::setValues(const string& component, const vector& values) } } +void Flow1D::getResiduals(const string& component, vector& values) const +{ + if (!m_state) { + throw CanteraError("Flow1D::getResiduals", + "Domain needs to be installed in a container."); + } + if (values.size() != nPoints()) { + throw ArraySizeError("Flow1D::getResiduals", values.size(), nPoints()); + } + auto i = componentIndex(component); + if (!componentActive(i)) { + warn_user( + "Flow1D::getResiduals", "Component '{}' is not used by '{}'.", + component, domainType()); + } + const double* soln = m_container->_workVector().data() + m_iloc; + for (size_t j = 0; j < nPoints(); j++) { + values[j] = soln[index(i,j)]; + } +} + void Flow1D::setProfile(const string& component, const vector& pos, const vector& values) { diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index 5fe6728c049..b5308a23523 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -50,7 +50,7 @@ void Sim1D::setInitialGuess(const string& component, vector& locs, } } -void Sim1D::setValue(size_t dom, size_t comp, size_t localPoint, double value) +void Sim1D::_setValue(size_t dom, size_t comp, size_t localPoint, double value) { size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint); AssertThrowMsg(iloc < m_state->size(), "Sim1D::setValue", @@ -58,7 +58,7 @@ void Sim1D::setValue(size_t dom, size_t comp, size_t localPoint, double value) (*m_state)[iloc] = value; } -double Sim1D::value(size_t dom, size_t comp, size_t localPoint) const +double Sim1D::_value(size_t dom, size_t comp, size_t localPoint) const { size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint); AssertThrowMsg(iloc < m_state->size(), "Sim1D::value", @@ -66,7 +66,7 @@ double Sim1D::value(size_t dom, size_t comp, size_t localPoint) const return (*m_state)[iloc]; } -double Sim1D::workValue(size_t dom, size_t comp, size_t localPoint) const +double Sim1D::_workValue(size_t dom, size_t comp, size_t localPoint) const { size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint); AssertThrowMsg(iloc < m_state->size(), "Sim1D::workValue", @@ -91,7 +91,7 @@ void Sim1D::setProfile(size_t dom, size_t comp, double zpt = d.z(n); double frac = (zpt - z0)/(z1 - z0); double v = linearInterp(frac, pos, values); - setValue(dom, comp, n, v); + _setValue(dom, comp, n, v); } } @@ -331,7 +331,7 @@ void Sim1D::setFlatProfile(size_t dom, size_t comp, double v) "Replaceable by Domain1D::setProfile."); size_t np = domain(dom).nPoints(); for (size_t n = 0; n < np; n++) { - setValue(dom, comp, n, v); + _setValue(dom, comp, n, v); } } @@ -467,7 +467,7 @@ int Sim1D::refine(int loglevel) // do the same for the solution at this point for (size_t i = 0; i < comp; i++) { - xnew.push_back(value(n, i, m)); + xnew.push_back(_value(n, i, m)); } // now check whether a new point is needed in the interval to @@ -482,7 +482,7 @@ int Sim1D::refine(int loglevel) // for each component, linearly interpolate the solution to this // point for (size_t i = 0; i < comp; i++) { - double xmid = 0.5*(value(n, i, m) + value(n, i, m+1)); + double xmid = 0.5*(_value(n, i, m) + _value(n, i, m+1)); xnew.push_back(xmid); } } @@ -556,8 +556,8 @@ int Sim1D::setFixedTemperature(double t) if (d_free && d_free->isFree()) { for (size_t m = 0; m < npnow - 1; m++) { bool fixedpt = false; - double t1 = value(n, c_offset_T, m); - double t2 = value(n, c_offset_T, m + 1); + double t1 = _value(n, c_offset_T, m); + double t2 = _value(n, c_offset_T, m + 1); // threshold to avoid adding new point too close to existing point double thresh = min(1., 1.e-1 * (t2 - t1)); z1 = d.z(m); @@ -589,7 +589,7 @@ int Sim1D::setFixedTemperature(double t) // do the same for the solution at this point for (size_t i = 0; i < comp; i++) { - xnew.push_back(value(n, i, m)); + xnew.push_back(_value(n, i, m)); } if (m == mfixed) { // add new point at zfixed (mfixed is not npos) @@ -599,7 +599,8 @@ int Sim1D::setFixedTemperature(double t) // for each component, linearly interpolate // the solution to this point for (size_t i = 0; i < comp; i++) { - double xmid = interp_factor*(value(n, i, m) - value(n, i, m+1)) + value(n,i,m+1); + double xmid = interp_factor*( + _value(n, i, m) - _value(n, i, m+1)) + _value(n,i,m+1); xnew.push_back(xmid); } } @@ -674,8 +675,8 @@ void Sim1D::setLeftControlPoint(double temperature) double current_val, next_val; for (size_t m = 0; m < np-1; m++) { - current_val = value(n,c_offset_T,m); - next_val = value(n,c_offset_T,m+1); + current_val = _value(n,c_offset_T,m); + next_val = _value(n,c_offset_T,m+1); if ((current_val - temperature) * (next_val - temperature) < 0.0) { // Pick the coordinate of the point with the temperature closest // to the desired temperature @@ -687,7 +688,7 @@ void Sim1D::setLeftControlPoint(double temperature) index = m+1; } d_axis.setLeftControlPointCoordinate(d_axis.z(index)); - d_axis.setLeftControlPointTemperature(value(n,c_offset_T,index)); + d_axis.setLeftControlPointTemperature(_value(n,c_offset_T,index)); return; } } @@ -725,8 +726,8 @@ void Sim1D::setRightControlPoint(double temperature) double current_val, next_val; for (size_t m = np-1; m > 0; m--) { - current_val = value(n,c_offset_T,m); - next_val = value(n,c_offset_T,m-1); + current_val = _value(n,c_offset_T,m); + next_val = _value(n,c_offset_T,m-1); if ((current_val - temperature) * (next_val - temperature) < 0.0) { // Pick the coordinate of the point with the temperature closest // to the desired temperature @@ -738,7 +739,7 @@ void Sim1D::setRightControlPoint(double temperature) index = m-1; } d_axis.setRightControlPointCoordinate(d_axis.z(index)); - d_axis.setRightControlPointTemperature(value(n,c_offset_T,index)); + d_axis.setRightControlPointTemperature(_value(n,c_offset_T,index)); return; } } diff --git a/test/clib/test_ctonedim.cpp b/test/clib/test_ctonedim.cpp index 56cbc73825d..f74e38ad1d4 100644 --- a/test/clib/test_ctonedim.cpp +++ b/test/clib/test_ctonedim.cpp @@ -194,9 +194,9 @@ TEST(ctonedim, freeflame) ASSERT_EQ(domain_nComponents(flow), 16); int32_t comp = domain_componentIndex(flow, "T"); ASSERT_EQ(comp, 2); - double Tprev = sim1D_value(flame, dom, comp, 0); + double Tprev = domain_value(dom, "T", 0); for (int32_t n = 0; n < domain_nPoints(flow); n++) { - T = sim1D_value(flame, dom, comp, n); + T = domain_value(dom, "T", n); ASSERT_GE(T + 1e-3, Tprev); Tprev = T; } diff --git a/test/oneD/test_oneD.cpp b/test/oneD/test_oneD.cpp index 79fc66636fd..19d007df95b 100644 --- a/test/oneD/test_oneD.cpp +++ b/test/oneD/test_oneD.cpp @@ -87,10 +87,9 @@ TEST(onedim, freeflame) } ASSERT_EQ(flow->nPoints(), static_cast(nz + 1)); - size_t comp = flow->componentIndex("T"); - double Tprev = flame->value(dom, comp, 0); + double Tprev = flow->value("T", 0); for (size_t n = 0; n < flow->nPoints(); n++) { - T = flame->value(dom, comp, n); + T = flow->value("T", n); ASSERT_GE(T, Tprev); Tprev = T; } diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index d933eccee53..fc2094baf66 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -1987,9 +1987,9 @@ def check_save_restore(self, jet, tol_T=None, tol_X=None): assert settings[k] == v assert list(jet.surface.surface.X) == approx(list(self.sim.surface.surface.X)) - for i in range(self.sim.surface.n_components): - assert self.sim.value("surface", i, 0) == \ - approx(jet.value("surface", i, 0), tol_X) + for k in self.sim.surface.phase.species_names: + assert self.sim.surface.value(k) == \ + approx(jet.surface.value(k), tol_X) jet.solve(loglevel=0)