diff --git a/include/cantera/base/SolutionArray.h b/include/cantera/base/SolutionArray.h index e0e5229aada..87e36a21134 100644 --- a/include/cantera/base/SolutionArray.h +++ b/include/cantera/base/SolutionArray.h @@ -129,8 +129,10 @@ 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) const; + bool hasComponent(const string& name, bool checkAlias=true) const; /** * Retrieve a component of the SolutionArray by name. @@ -404,6 +406,9 @@ class SolutionArray vector m_active; //!< Vector of locations referencing active entries }; +//! Return mapping of component alias names to standardized component names. +const map& _componentAliasMap(); + } #endif 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 f7d53003d0c..05bd0620e7e 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -196,8 +196,23 @@ class Domain1D m_name[n] = name; } - //! index of component with name `name`. - virtual size_t componentIndex(const string& name) const; + /** + * Index of component with name `name`. + * @param name name of component + * @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. + * @param name name of component + * @param checkAlias if `true` (default), check alias mapping + * + * @since New in %Cantera 3.2. + */ + virtual bool hasComponent(const string& name, bool checkAlias=true) const; /** * Set the upper and lower bounds for a solution component, n. @@ -348,6 +363,60 @@ class Domain1D return x[index(n,j)]; } + /** + * Retrieve component 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::getValues", "Needs to be overloaded."); + } + + /** + * Specify component 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::setValues", "Needs to be overloaded."); + } + + /** + * Specify a profile for a 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 + * 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. + */ + 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. + * @param component Name of the component. + * @param value Constant value. + * + * @since New in %Cantera 3.2. + */ + virtual void setFlatProfile(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 @@ -355,9 +424,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. @@ -368,7 +440,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 /*! @@ -376,9 +450,16 @@ 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."); + // create a shared_ptr that skips deletion of the held pointer + shared_ptr arr_ptr(&arr, [](SolutionArray*){}); + fromArray(arr_ptr); } //! Restore the solution for this domain from a SolutionArray. @@ -387,7 +468,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. @@ -499,6 +582,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 5eaeb444683..060ade16d61 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" @@ -175,15 +176,23 @@ 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; void show(const double* x) override; - shared_ptr asArray(const double* soln) const override; - void fromArray(SolutionArray& arr, double* soln) 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 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. @@ -548,15 +557,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 - * 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. @@ -571,7 +580,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 +683,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/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..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); @@ -92,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 3019e6f326e..b412ba22295 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -23,7 +23,14 @@ 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 + 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 + 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) @@ -37,8 +44,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 bfb3ff05e13..b7aea240569 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -86,6 +86,95 @@ 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)) + + @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. + + :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_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. @@ -252,20 +341,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): @@ -456,12 +531,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 named + ``spreadRate`` in the C++ code and 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: """ @@ -935,7 +1087,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 `Domain1D.get_values`. """ + 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 @@ -958,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) @@ -980,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/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/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index b5467776c3f..17e253a8528 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) @@ -257,6 +261,9 @@ def boundary_emissivities(self, epsilon): @property def grid(self): """ Array of grid point positions along the flame. """ + warnings.warn("To be removed in a future version of Cantera. " + "Replaceable by 'Domain1D.grid'.", + PendingDeprecationWarning) return self.flame.grid @property @@ -271,14 +278,20 @@ def P(self, P): @property def T(self): """ Array containing the temperature [K] at each grid point. """ - return self.profile(self.flame, 'T') + warnings.warn("To be removed in a future version of Cantera. " + "Replaceable by 'Domain1D.T'.", + PendingDeprecationWarning) + 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.profile(self.flame, 'velocity') + warnings.warn("To be removed in a future version of Cantera. " + "Replaceable by 'Domain1D.velocity'.", + PendingDeprecationWarning) + return self.flame.get_values("velocity") @property def spread_rate(self): @@ -286,33 +299,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') + warnings.warn("To be removed in a future version of Cantera. " + "Replaceable by 'Domain1D.spread_rate'.", + PendingDeprecationWarning) + return self.flame.get_values("spreadRate") @property 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') + warnings.warn("To be removed in a future version of Cantera. " + "Replaceable by 'Domain1D.radial_pressure_gradient'.", + PendingDeprecationWarning) + return self.flame.get_values("Lambda") @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``. """ - 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') + warnings.warn("To be removed in a future version of Cantera. " + "Replaceable by 'Domain1D.electric_field'.", + PendingDeprecationWarning) + return self.flame.get_values("eField") @property def Uo(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. """ - return self.profile(self.flame, 'Uo') + warnings.warn("To be removed in a future version of Cantera. " + "Replaceable by 'Domain1D.oxidizer_velocity'.", + PendingDeprecationWarning) + return self.flame.get_values("Uo") @property def left_control_point_temperature(self): @@ -664,8 +690,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 @@ -680,7 +706,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): @@ -859,10 +885,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): @@ -924,9 +950,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) @@ -1045,12 +1071,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('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): - 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 @@ -1199,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') @@ -1338,20 +1364,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("spreadRate", locs, [0.0, 0.0]) class CounterflowPremixedFlame(FlameBase): @@ -1439,10 +1464,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 @@ -1453,9 +1478,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("spreadRate", [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): @@ -1526,10 +1551,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 @@ -1537,6 +1562,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("spreadRate", [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 0a38fda0728..238692c984d 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctdomain_auto.yaml @@ -19,6 +19,13 @@ recipes: - name: nPoints - name: componentName - name: componentIndex + wraps: componentIndex(const string&) +- 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 c86588e69e3..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 # New in Cantera 3.2 +- 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/samples/python/onedim/diffusion_flame_batch.py b/samples/python/onedim/diffusion_flame_batch.py index 58dda388b1c..986dcd6e2a5 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( + "spreadRate", f.flame.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.flame.set_values( + "Lambda", f.flame.radial_pressure_gradient * rel_pressure_increase ** exp_lam_p) try: # Try solving the flame @@ -186,17 +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("spreadRate", f.flame.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature - f.set_profile('lambda', normalized_grid, f.L * 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 43d09603b00..ac377dee306 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,17 +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("spreadRate", f.flame.spread_rate * strain_factor ** exp_V_a) # Update pressure curvature - f.set_profile('lambda', normalized_grid, f.L * 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: diff --git a/src/base/SolutionArray.cpp b/src/base/SolutionArray.cpp index d501d0b66ec..5b88dd716c3 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,8 +40,28 @@ const std::map aliasMap = { {"Q", "vapor-fraction"}, }; -namespace Cantera +const map reverseAliasMap = { + // reserved names used by states + {"temperature", "T"}, + {"pressure", "P"}, + {"density", "D"}, + // reserved names used for 1-D objects + {"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 + {"E", "eField"}, + {"electric-field", "eField"}, + {"oxidizer-velocity", "Uo"}, +}; + +} // end unnamed namespace + +const map& _componentAliasMap() { + return reverseAliasMap; +} SolutionArray::SolutionArray(const shared_ptr& sol, int size, const AnyMap& meta) @@ -649,7 +674,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 @@ -663,21 +688,28 @@ bool SolutionArray::hasComponent(const string& name) const // reserved names return false; } + if (checkAlias && 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 (!hasComponent(name, false) && 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 +741,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(); @@ -1799,13 +1831,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"); @@ -1813,11 +1849,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); } } } @@ -1908,12 +1948,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 { @@ -1922,10 +1966,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); } } 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); } } diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp index b4a41b3976c..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 ---------------- @@ -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 != "") { @@ -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 d3a77b8f114..9dee0a695ba 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) @@ -146,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 asArray."); - } - 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) @@ -285,6 +272,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 a9943682a4c..97301e7cf0c 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -16,6 +16,19 @@ using namespace std; namespace Cantera { +namespace { // restrict scope of auxiliary variables to local translation unit + +const map componentMap = { + {"velocity", c_offset_U}, // axial velocity [m/s] + {"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 + {"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) @@ -60,7 +73,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 +133,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 +565,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 +619,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 +636,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 +644,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 +663,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; } @@ -809,11 +828,11 @@ 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: - return "lambda"; + return "Lambda"; case c_offset_E: return "eField"; case c_offset_Uo: @@ -827,37 +846,48 @@ 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 { - for (size_t n=c_offset_Y; n Flow1D::asArray(const double* soln) 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()); + } + auto i = componentIndex(component); + if (!componentActive(i)) { + warn_user( + "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(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::_etValues", values.size(), nPoints()); + } + auto i = componentIndex(component); + if (!componentActive(i)) { + throw CanteraError( + "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(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 " + "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* soln = m_state->data() + m_iloc; + 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(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 '{}'.", + component, domainType()); + } + double* soln = m_state->data() + m_iloc; + for (size_t j = 0; j < nPoints(); j++) { + soln[index(i,j)] = v; + } +} + +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 @@ -950,27 +1078,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")) { + // edge case: 'lambda' is renamed to 'Lambda' in Cantera 3.2 + const auto data = arr->getComponent("lambda").as>(); for (size_t j = 0; j < nPoints(); j++) { soln[index(i,j)] = data[j]; } @@ -980,7 +1125,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/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..5fe6728c049 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& 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 " @@ -98,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; @@ -113,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; @@ -127,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); } @@ -271,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{}", @@ -300,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{}", @@ -323,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); 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/data/flame_lambda.h5 b/test/data/flame_lambda.h5 new file mode 100644 index 00000000000..c42e5f6bb6d Binary files /dev/null and b/test/data/flame_lambda.h5 differ 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/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 3bb891423dd..2d545e57835 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -113,7 +113,7 @@ def test_tolerances(self): # Some things don't work until the domains have been added to a Sim1D sim = ct.Sim1D((left, flame, right)) - with pytest.raises(ct.CanteraError, match='no component'): + with pytest.raises(ct.CanteraError, match='No component'): flame.set_steady_tolerances(foobar=(3e-4, 3e-6)) flame.set_steady_tolerances(default=(5e-3, 5e-5), @@ -252,6 +252,30 @@ def test_fixed_temperature(self): assert self.sim.fixed_temperature == approx(tfixed) assert self.sim.fixed_temperature_location == approx(zfixed) + # test some alias names + 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.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'"): + self.sim.flame.set_values("L", zeros) + with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): + self.sim.flame.set_values("Lambda", zeros) + with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): + self.sim.flame.set_values("electric-field", zeros) + with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): + self.sim.flame.set_values("E", zeros) + with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): + self.sim.flame.set_values("eField", zeros) + with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): + self.sim.flame.set_values("oxidizer-velocity", zeros) + with pytest.raises(ct.CanteraError, match="is not used by 'free-flow'"): + self.sim.flame.set_values("Uo", zeros) + @pytest.mark.slow_test def test_auto_width(self): Tin = 300 @@ -701,7 +725,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") @@ -729,24 +752,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 @@ -1036,6 +1055,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.radial_pressure_gradient[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 +1290,45 @@ def test_bad_boundary_conditions(self): sim.fuel_inlet.X = "H2: 1.0" sim.set_initial_guess() + def run_restore_diffusionflame(self, fname): + gas = ct.Solution("h2o2.yaml") + sim = ct.CounterflowDiffusionFlame(gas) + sim.restore(fname) + 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.flame.set_values("electric-field", zeros) + with pytest.raises(ct.CanteraError, match="is not used by 'axisymmetric-flow'"): + sim.flame.set_values("oxidizer-velocity", zeros) + + def run_restore_solutionarray(self, fname): + gas = ct.Solution("h2o2.yaml") + arr = ct.SolutionArray(gas) + arr.restore(fname, "solution/flame") + + radial_lambda = np.unique(arr.radial_pressure_gradient) + assert len(radial_lambda) == 1 + 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] + assert getattr(arr, "radial-pressure-gradient")[0] == radial_lambda[0] + + def test_restore_lambda_yaml(self): + self.run_restore_diffusionflame("flame_lambda.yaml") + self.run_restore_solutionarray("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_diffusionflame("flame_lambda.h5") + self.run_restore_solutionarray("flame_lambda.h5") + def test_two_point_control(self): """ Computes a solution using a standard counterflow diffusion flame as a @@ -1331,7 +1397,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.flame.oxidizer_velocity == approx(sim.velocity[-1]) temperature_4 = sim.T # Check the difference between the un-perturbed two-point solution and the @@ -1376,7 +1442,7 @@ def test_two_point_control(self): assert (sim.right_control_point_temperature == approx(original_settings['right-temperature'], 1e-4)) - assert sim.Uo == approx(sim.velocity[-1]) + assert sim.flame.oxidizer_velocity == approx(sim.velocity[-1]) # Test - Check error conditions sim = ct.CounterflowDiffusionFlame(gas, width=width) @@ -1493,7 +1559,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.L, sim.L[0]) + radial_lambda = sim.flame.radial_pressure_gradient + assert np.allclose(radial_lambda, radial_lambda[0]) return sim return _run_case @@ -1610,7 +1677,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.L, sim.L[0]) + radial_lambda = sim.flame.radial_pressure_gradient + assert np.allclose(radial_lambda, radial_lambda[0]) return sim return _run_case @@ -1657,7 +1725,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.L, 0) return _run_case @pytest.mark.parametrize( @@ -1750,7 +1817,8 @@ 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]) + radial_lambda = sim.flame.radial_pressure_gradient + assert np.allclose(radial_lambda, radial_lambda[0]) self.sim = sim def test_stagnation_case1(self): @@ -1937,7 +2005,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.L, sim.L[0]) + radial_lambda = sim.flame.radial_pressure_gradient + assert np.allclose(radial_lambda, radial_lambda[0]) return sim def test_restart(self): @@ -2009,14 +2078,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.E)) == 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.E) == approx(149.63179056676853, rel=1e-3) + assert max(self.sim.flame.electric_field) == approx(149.6317905667685, rel=1e-3) + with pytest.warns(PendingDeprecationWarning, match="'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.flame.set_values("oxidizer-velocity", zeros) class TestIonBurnerFlame: @@ -2038,5 +2112,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.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) diff --git a/test_problems/cxx_samples/cxx_flamespeed_blessed.txt b/test_problems/cxx_samples/cxx_flamespeed_blessed.txt index cf454a77899..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