SolutionArray/HDF interface for 1D FlameBase

include/cantera/oneD/Domain1D.h

@@ -182,7 +182,7 @@ class Domain1D
     }
 
     //! index of component with name \a name.
-    size_t componentIndex(const std::string& name) const;
+    virtual size_t componentIndex(const std::string& name) const;
 
     void setBounds(size_t n, doublereal lower, doublereal upper) {
         m_min[n] = lower;

include/cantera/oneD/Sim1D.h

@@ -129,15 +129,22 @@ class Sim1D : public OneDim
     int refine(int loglevel=0);
 
     //! Add node for fixed temperature point of freely propagating flame
-    int setFixedTemperature(doublereal t);
+    int setFixedTemperature(double t);
+
+    //! Return temperature at the point used to fix the flame location
+    double fixedTemperature();
+
+    //! Return location of the point where temperature is fixed
+    double fixedTemperatureLocation();
 
     /**
      * Set grid refinement criteria. If dom >= 0, then the settings
      * apply only to the specified domain.  If dom < 0, the settings
      * are applied to each domain.  @see Refiner::setCriteria.
      */
-    void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
-                           doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
+    void setRefineCriteria(int dom = -1, double ratio = 10.0,
+                           double slope = 0.8, double curve = 0.8,
+                           double prune = -0.1);
 
     /**
      * Get the grid refinement criteria. dom must be greater than

include/cantera/oneD/StFlow.h

@@ -198,13 +198,24 @@ class StFlow : public Domain1D
         return m_do_radiation;
     }
 
+    //! Return radiative heat loss at grid point j
+    double radiativeHeatLoss(size_t j) const {
+        return m_qdotRadiation[j];
+    }
+
     //! Set the emissivities for the boundary values
     /*!
      * Reads the emissivities for the left and right boundary values in the
      * radiative term and writes them into the variables, which are used for the
      * calculation.
      */
-    void setBoundaryEmissivities(doublereal e_left, doublereal e_right);
+    void setBoundaryEmissivities(double e_left, double e_right);
+
+    //! Return emissivitiy at left boundary
+    double leftEmissivity() const { return m_epsilon_left; }
+
+    //! Return emissivitiy at right boundary
+    double rightEmissivity() const { return m_epsilon_right; }
 
     void fixTemperature(size_t j=npos);
 
@@ -215,7 +226,14 @@ class StFlow : public Domain1D
     //! Change the grid size. Called after grid refinement.
     virtual void resize(size_t components, size_t points);
 
-    virtual void setFixedPoint(int j0, doublereal t0) {}
+    /*!
+     * @deprecated To be removed after Cantera 2.5.
+     */
+    virtual void setFixedPoint(int j0, doublereal t0) {
+        // this does nothing and does not appear to be overloaded
+        warn_deprecated("StFlow::setFixedPoint",
+                        "To be removed after Cantera 2.5.");
+    }
 
     //! Set the gas object state to be consistent with the solution at point j.
     void setGas(const doublereal* x, size_t j);

interfaces/cython/cantera/_cantera.pxd

@@ -750,16 +750,20 @@ cdef extern from "cantera/oneD/StFlow.h":
         void setPressure(double)
         void enableRadiation(cbool)
         cbool radiationEnabled()
+        double radiativeHeatLoss(size_t)
         double pressure()
         void setFixedTempProfile(vector[double]&, vector[double]&)
         void setBoundaryEmissivities(double, double)
+        double leftEmissivity()
+        double rightEmissivity()
         void solveEnergyEqn()
         void fixTemperature()
         cbool doEnergy(size_t)
         void enableSoret(cbool) except +translate_exception
         cbool withSoret()
         void setFreeFlow()
         void setAxisymmetricFlow()
+        string flowType()
 
 
 cdef extern from "cantera/oneD/IonFlow.h":
@@ -815,6 +819,8 @@ cdef extern from "cantera/oneD/Sim1D.h":
         size_t maxGridPoints(size_t) except +translate_exception
         void setGridMin(int, double) except +translate_exception
         void setFixedTemperature(double) except +translate_exception
+        double fixedTemperature()
+        double fixedTemperatureLocation()
         void setInterrupt(CxxFunc1*) except +translate_exception
         void setTimeStepCallback(CxxFunc1*)
         void setSteadyCallback(CxxFunc1*)

interfaces/cython/cantera/composite.py

@@ -731,12 +731,16 @@ def restore_data(self, data, labels):
                 if s in valid_species:
                     state_data[-1][:, i] = data[:, valid_species[s]]
 
-        # labels may include calculated properties that must not be restored
-        calculated = self._scalar + self._n_species + self._n_reactions
-        exclude = [l for l in labels
-                   if any([v in l for v in calculated])]
-        extra = {l: list(data[:, i]) for i, l in enumerate(labels)
-                 if l not in exclude}
+        # labels may include calculated properties that must not be restored:
+        # compare column labels to names that are reserved for SolutionArray
+        # attributes (see `SolutionArray.collect_data`), i.e. scalar values,
+        # arrays with number of species, and arrays with number of reactions.
+        exclude = [lab for lab in labels
+                   if any([lab in self._scalar,
+                           '_'.join(lab.split('_')[:-1]) in self._n_species,
+                           lab.split(' ')[0] in self._n_reactions])]
+        extra = {lab: list(data[:, i]) for i, lab in enumerate(labels)
+                 if lab not in exclude}
         if len(self._extra):
             extra_lists = {k: extra[k] for k in self._extra}
         else:
@@ -861,7 +865,7 @@ def write_csv(self, filename, cols=None, *args, **kwargs):
         only with 1D `SolutionArray` objects.
         """
         data, labels = self.collect_data(cols=cols, *args, **kwargs)
-        with open(filename, 'w') as outfile:
+        with open(filename, 'w', newline='') as outfile:
             writer = _csv.writer(outfile)
             writer.writerow(labels)
             for row in data:

interfaces/cython/cantera/examples/onedim/adiabatic_flame.py

@@ -6,6 +6,11 @@
 """
 
 import cantera as ct
+try:
+    import pandas as pd
+    import tables
+except:
+    pd = None
 
 # Simulation parameters
 p = ct.one_atm  # pressure [Pa]
@@ -29,16 +34,25 @@
 f.solve(loglevel=loglevel, auto=True)
 
 # Solve with the energy equation enabled
-f.save('h2_adiabatic.xml', 'mix', 'solution with mixture-averaged transport')
+if pd:
+    f.write_hdf('h2_adiabatic.h5', key='mix', mode='w')
+else:
+    f.save('h2_adiabatic.xml', 'mix',
+           'solution with mixture-averaged transport')
+
 f.show_solution()
-print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.u[0]))
+print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.velocity[0]))
 
 # Solve with multi-component transport properties
 f.transport_model = 'Multi'
 f.solve(loglevel)  # don't use 'auto' on subsequent solves
 f.show_solution()
-print('multicomponent flamespeed = {0:7f} m/s'.format(f.u[0]))
-f.save('h2_adiabatic.xml', 'multi', 'solution with multicomponent transport')
+print('multicomponent flamespeed = {0:7f} m/s'.format(f.velocity[0]))
+if pd:
+    f.write_hdf('h2_adiabatic.h5', key='multi')
+else:
+    f.save('h2_adiabatic.xml', 'multi',
+           'solution with multicomponent transport')
 
 # write the velocity, temperature, density, and mole fractions to a CSV file
 f.write_csv('h2_adiabatic.csv', quiet=False)

interfaces/cython/cantera/examples/onedim/diffusion_flame_batch.py

@@ -111,10 +111,10 @@ def interrupt_extinction(t):
     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('u', normalized_grid,
-                  f.u * rel_pressure_increase ** exp_u_p)
-    f.set_profile('V', normalized_grid,
-                  f.V * rel_pressure_increase ** exp_V_p)
+    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)
     # Update pressure curvature
     f.set_profile('lambda', normalized_grid,
                   f.L * rel_pressure_increase ** exp_lam_p)
@@ -167,8 +167,10 @@ def interrupt_extinction(t):
     f.fuel_inlet.mdot *= strain_factor ** exp_mdot_a
     f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a
     # Update velocities
-    f.set_profile('u', normalized_grid, f.u * strain_factor ** exp_u_a)
-    f.set_profile('V', normalized_grid, f.V * strain_factor ** exp_V_a)
+    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)
     # Update pressure curvature
     f.set_profile('lambda', normalized_grid, f.L * strain_factor ** exp_lam_a)
     try:
@@ -203,7 +205,7 @@ def interrupt_extinction(t):
 
     # Plot the axial velocity profiles (normalized by the fuel inlet velocity)
     # for selected pressures
-    ax2.plot(f.grid / f.grid[-1], f.u / f.u[0],
+    ax2.plot(f.grid / f.grid[-1], f.velocity / f.velocity[0],
              label='{0:05.1f} bar'.format(p))
 
 ax1.legend(loc=0)
@@ -231,7 +233,7 @@ def interrupt_extinction(t):
 
     # Plot the axial velocity profiles (normalized by the fuel inlet velocity)
     # for the strain rate loop (selected)
-    ax4.plot(f.grid / f.grid[-1], f.u / f.u[0],
+    ax4.plot(f.grid / f.grid[-1], f.velocity / f.velocity[0],
              label=format(a_max, '.2e') + ' 1/s')
 
 ax3.legend(loc=0)

interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py

@@ -88,7 +88,7 @@
 # List of peak temperatures
 T_max = [np.max(f.T)]
 # List of maximum axial velocity gradients
-a_max = [np.max(np.abs(np.gradient(f.u) / np.gradient(f.grid)))]
+a_max = [np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid)))]
 
 # Simulate counterflow flames at increasing strain rates until the flame is
 # extinguished. To achieve a fast simulation, an initial coarse strain rate
@@ -109,8 +109,10 @@
     f.fuel_inlet.mdot *= strain_factor ** exp_mdot_a
     f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a
     # Update velocities
-    f.set_profile('u', normalized_grid, f.u * strain_factor ** exp_u_a)
-    f.set_profile('V', normalized_grid, f.V * strain_factor ** exp_V_a)
+    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)
     # Update pressure curvature
     f.set_profile('lambda', normalized_grid, f.L * strain_factor ** exp_lam_a)
     try:
@@ -125,7 +127,7 @@
                description='Cantera version ' + ct.__version__ +
                ', reaction mechanism ' + reaction_mechanism)
         T_max.append(np.max(f.T))
-        a_max.append(np.max(np.abs(np.gradient(f.u) / np.gradient(f.grid))))
+        a_max.append(np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid))))
         # If the temperature difference is too small and the minimum relative
         # strain rate increase is reached, abort
         if ((T_max[-2] - T_max[-1] < delta_T_min) &

interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py

@@ -24,7 +24,7 @@
 f.set_refine_criteria(ratio=3, slope=0.07, curve=0.14)
 
 f.solve(loglevel=1, auto=True)
-print('\nmixture-averaged flamespeed = {:7f} m/s\n'.format(f.u[0]))
+print('\nmixture-averaged flamespeed = {:7f} m/s\n'.format(f.velocity[0]))
 
 # Use the adjoint method to calculate sensitivities
 sens = f.get_flame_speed_reaction_sensitivities()

interfaces/cython/cantera/examples/onedim/ion_free_flame.py

@@ -31,7 +31,7 @@
 
 f.save('CH4_adiabatic.xml', 'ion', 'solution with ionized gas transport')
 f.show_solution()
-print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.u[0]))
+print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.velocity[0]))
 
 # write the velocity, temperature, density, and mole fractions to a CSV file
 f.write_csv('CH4_adiabatic.csv', quiet=False)

interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py

@@ -29,12 +29,12 @@ def derivative(x, y):
 
 def computeStrainRates(oppFlame):
     # Compute the derivative of axial velocity to obtain normal strain rate
-    strainRates = derivative(oppFlame.grid, oppFlame.u)
+    strainRates = derivative(oppFlame.grid, oppFlame.velocity)
 
     # Obtain the location of the max. strain rate upstream of the pre-heat zone.
     # This is the characteristic strain rate
     maxStrLocation = abs(strainRates).argmax()
-    minVelocityPoint = oppFlame.u[:maxStrLocation].argmin()
+    minVelocityPoint = oppFlame.velocity[:maxStrLocation].argmin()
 
     # Characteristic Strain Rate = K
     strainRatePoint = abs(strainRates[:minVelocityPoint]).argmax()
@@ -130,15 +130,17 @@ def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1,
 
     # Axial Velocity Plot
     plt.subplot(1, 2, 1)
-    plt.plot(oppFlame.grid, oppFlame.u, 'r', lw=2)
+    plt.plot(oppFlame.grid, oppFlame.velocity, 'r', lw=2)
     plt.xlim(oppFlame.grid[0], oppFlame.grid[-1])
     plt.xlabel('Distance (m)')
     plt.ylabel('Axial Velocity (m/s)')
 
     # Identify the point where the strain rate is calculated
-    plt.plot(oppFlame.grid[strainRatePoint], oppFlame.u[strainRatePoint], 'gs')
+    plt.plot(oppFlame.grid[strainRatePoint],
+             oppFlame.velocity[strainRatePoint], 'gs')
     plt.annotate('Strain-Rate point',
-                 xy=(oppFlame.grid[strainRatePoint], oppFlame.u[strainRatePoint]),
+                 xy=(oppFlame.grid[strainRatePoint],
+                     oppFlame.velocity[strainRatePoint]),
                  xytext=(0.001, 0.1),
                  arrowprops={'arrowstyle': '->'})