CSV writer implementation

examples/hyshot_ii_jic/hyshot_ii_jic.py

@@ -11,6 +11,7 @@
 from stanshock.models.jicf import JICModel
 from stanshock.numerics.boundary_conditions import BCInput, SpecifiedFace
 from stanshock.physics.flamelet import FPVTable
+from stanshock.processing.csv_writer import CSVWriter
 from stanshock.processing.initialize import InitializeConstant
 from stanshock.processing.plot import XTDiagram
 from stanshock.system.geometry import Box
@@ -248,127 +249,6 @@ def fuel_props_from_phi(phi_gl):
     physics=fpv_table,
 )
 
-###################################################################
-
-# i_m = 1
-# z_plot = jic.z_inj[1]
-# i_z = np.argmin(np.abs(jic.z_3D_data - z_plot))
-
-# fig, ax = plt.subplots()
-# c = ax.contourf(jic.x_3D_data*scale,
-#                 jic.y_3D_data*scale,
-#                 jic.Z_3D_data[i_m,:,:,i_z].T,
-#                 levels=np.linspace(0, 0.3, 100))
-# y_cl_max = jic.y_cl(L - x_inj)[i_m]
-# y_cl_arr = np.linspace(0, y_cl_max, 1000)
-# x_cl_arr = jic.x_cl_from_y_cl(y_cl_arr[:,np.newaxis])[:,i_m]
-# ax.plot((x_cl_arr + x_inj)*scale, y_cl_arr*scale, 'r')
-# ax.set_xlabel(r'$x$ [mm]')
-# ax.set_ylabel(r'$y$ [mm]')
-# ax.set_xlim((50.0, 100.0))
-# cbar = plt.colorbar(c, ticks=[0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
-# cbar.set_label(r'$Z$')
-# ax.set_aspect('equal')
-# plt.savefig(os.path.join(figdir, "injector_xy.png"), bbox_inches='tight', dpi=300)
-
-
-# x_plot = x_inj + np.array([0.0e-3,
-#                            25.0e-3,
-#                            50.0e-3])
-# i_plot = np.array([np.argmin(np.abs(x - x_i)) for x_i in x_plot])
-
-# fig, axs = plt.subplots(len(x_plot), 1, sharex=True, figsize=(6, 8))
-# for i in range(len(x_plot)):
-#     i_x = i_plot[i]
-#     c = axs[i].contourf(jic.z_3D_data*scale,
-#                         jic.y_3D_data*scale,
-#                         jic.Z_3D_data[i_m,i_x,:,:],
-#                         levels=np.linspace(0, 0.3, 100))
-#     for i_inj in range(N_f):
-#         axs[i].scatter(jic.z_inj[i_inj]*scale,
-#                        jic.y_cl(jic.x_3D_data[i_x] - x_inj)[i_m]*scale,
-#                        color='r', marker='x')
-#     axs[i].set_ylabel(r'$y$ [mm]')
-#     axs[i].set_aspect('equal')
-#     # plt.colorbar(c, ax=axs[i_x])
-# axs[-1].set_xlabel(r'$z$ [mm]')
-# plt.tight_layout()
-# plt.savefig(os.path.join(figdir, "injector_zy.png"), bbox_inches='tight', dpi=300)
-
-# breakpoint()
-
-###################################################################
-
-# # Evaluate the mixture fraction PDFs along the axial direction
-# print("Evaluating the mixture fraction PDFs along the axial direction...")
-# Z_pdf = []
-# Z_pdf_samples = []
-# Z_pdf_weights = []
-# Z_avg_simple = np.zeros_like(x)
-# Z_avg = np.zeros_like(x)
-# Z_var = np.zeros_like(x)
-# i_last_const = np.where(x <= L_const)[0][-1]
-
-# Z_plot = np.linspace(0, 1, 1000)
-# fig, ax = plt.subplots()
-# ax.set_xlabel(r'$Z$')
-# ax.set_ylabel(r'$p(Z)$')
-
-# for i_x, x_i in enumerate(tqdm(x[:i_last_const])):
-#     func = lambda z, y : jic.Z_3D(x_i, y, z)
-#     grad_func = lambda z, y : jic.grad_Z_3D(x_i, y, z)
-#     ranges = ((-w/2, w/2), (0, h[i_x]))
-#     grid_dims = (int(np.ceil((ranges[0][1] - ranges[0][0]) / dx_Z_pdf)),
-#                  int(np.ceil((ranges[1][1] - ranges[1][0]) / dx_Z_pdf)))
-#     sampler = MonteCarloSampler2D(
-#         ranges, func, grad_func, grid_dims,
-#         alpha=0.0, num_samples_per_iter=5000, max_iters=10)
-#     Z_pdf_i, z_samples, y_samples, Z_pdf_samples_i, Z_pdf_weights_i = sampler.compute_scalar_pdf()
-#     Z_pdf.append(Z_pdf_i)
-#     Z_pdf_samples.append(Z_pdf_samples_i)
-#     Z_pdf_weights.append(Z_pdf_weights_i)
-
-#     # Plot, mapping color to x coordinate
-#     ax.plot(Z_plot, Z_pdf_i(Z_plot), c=plt.cm.viridis(i_x / (len(x) - 1)))
-#     fig.savefig(os.path.join(figdir, "Z_pdf.png"), bbox_inches='tight', dpi=300)
-
-# for i_x in range(i_last_const, len(x)):
-#     Z_pdf.append(Z_pdf[i_last_const-1])
-#     Z_pdf_samples.append(Z_pdf_samples[i_last_const-1])
-#     Z_pdf_weights.append(Z_pdf_weights[i_last_const-1])
-
-#     # Plot, mapping color to x coordinate
-#     ax.plot(Z_plot, Z_pdf[i_last_const-1](Z_plot), c=plt.cm.viridis(i_x / (len(x) - 1)))
-#     fig.savefig(os.path.join(figdir, "Z_pdf.png"), bbox_inches='tight', dpi=300)
-
-# print("Computing the mean and variance of the PDFs...")
-# for i_x in range(len(x)):
-#     Z_avg[i_x] = integrate.quad(lambda Z : Z * Z_pdf[i_x](Z), 0, 1)[0]
-#     Z_var[i_x] = integrate.quad(lambda Z : (Z - Z_avg[i_x])**2 * Z_pdf[i_x](Z), 0, 1)[0]
-
-# np.save(os.path.join(datadir, "Z_avg.npy"), Z_avg)
-# np.save(os.path.join(datadir, "Z_var.npy"), Z_var)
-# with open(os.path.join(datadir, "Z_pdf.pkl"), 'wb') as f:
-#     pickle.dump((Z_pdf_samples, Z_pdf_weights), f)
-
-###################################################################
-
-# fig, ax = plt.subplots(figsize=(4, 3.2))
-# ax.plot(x*scale, jic.Z_avg_profile, 'b')
-# ax.axhline(Z_gl, color='b', linestyle='--')
-# ax.set_ymargin(0.1)
-# ax.set_xlabel(r'$x$ [mm]')
-# ax.set_ylabel(r'$\langle Z \rangle$', color='b')
-# ax.tick_params(axis='y', labelcolor='b')
-# ax1 = ax.twinx()
-# ax1.semilogy(x*scale, jic.Z_var_profile, 'r')
-# ax1.set_ymargin(0.1)
-# ax1.set_ylabel(r"$\langle Z''^2 \rangle$", color='r')
-# ax1.tick_params(axis='y', labelcolor='r')
-# plt.savefig(os.path.join(figdir, "Z_avg_var.png"), bbox_inches='tight', dpi=300)
-
-###################################################################
-
 # Initialize and run the simulation
 ss = Combustor(
     geometry=geometry,
@@ -387,6 +267,15 @@ def fuel_props_from_phi(phi_gl):
     use_double_flux=False,
 )
 
+# Update CSV writer initialization to match plot_state_interval
+csv_writer = CSVWriter(
+    combustor=ss,
+    filename=figdir / "data.csv",  # Will become test_00000.csv, test_00001.csv, etc.
+    interval=100,  # Same as plot_state_interval=100
+    wall_temperature=300.0,
+)
+ss.csv_writers = [csv_writer]
+
 plot_variables = [
     "density",
     "velocity",

examples/umdci/constant_area_duct.py

@@ -14,7 +14,9 @@
 )
 from stanshock.numerics.boundary_conditions import BCInput, SpecifiedFace
 from stanshock.physics.thermotable import ThermoTable
+from stanshock.processing.csv_writer import CSVWriter
 from stanshock.processing.initialize import InitializeConstant
+from stanshock.processing.plot import get_variable_info_map
 from stanshock.system.geometry import Box
 
 """
@@ -53,6 +55,7 @@
 tFinal = 0.01
 
 physics_model = ThermoTable(gas1)
+variable_info_map = get_variable_info_map(physics_model)
 BC_inlet = SpecifiedFace(
     location="left", reference_state=(gas1.density, u1, gas1.P, (1.0,))
 )
@@ -72,8 +75,20 @@
     include_diffusion=False,
     output_every=100,
     plot_state_interval=100,
+    plot_state_variables=["mach", "p", "T"],
+    plot_state_variable_info_map=variable_info_map,
     use_double_flux=False,
 )
+ss.csv_writers = [
+    CSVWriter(
+        combustor=ss,
+        filename=figdir / "csv" / "state.csv",
+        interval=ss.plot_state_interval,
+        variables=["x", "mach", "p", "T"],
+        variable_info_map=variable_info_map,
+    )
+]
+
 try:
     t0 = time.perf_counter()
     ss.advance_simulation(tFinal)

src/stanshock/components/combustor.py

@@ -36,8 +36,9 @@
 from stanshock.numerics.viscous_flux import ViscousFlux
 from stanshock.physics.chemistry_source import ChemistrySource, ConstantVolumeChemistry
 from stanshock.physics.fluid_base import FluidPhysics
+from stanshock.processing.csv_writer import CSVWriter
 from stanshock.processing.initialize import Initialization, InitializeRestart
-from stanshock.processing.plot import XTDiagram, plot_state
+from stanshock.processing.plot import VariableInfo, XTDiagram, plot_state
 from stanshock.processing.probe import Probe
 from stanshock.system.backend import Array
 from stanshock.system.base import RightHandSide
@@ -77,6 +78,8 @@ def __init__(
         include_diffusion: bool = False,  # exclude diffusion
         thickening: None = None,  # thickening function
         plot_state_interval: int = -1,  # plot the state every n iterations
+        plot_state_variables: list[str | list[str]] | None = None,
+        plot_state_variable_info_map: dict[str, VariableInfo] | None = None,
         iteration: int = 0,  # Iteration to start from
         n_restart_interval: int = -1,  # If >0, saves the fluid state to a file every n_restart_interval iterations
     ) -> None:
@@ -97,8 +100,11 @@ def __init__(
         self.include_diffusion = include_diffusion
         self.thickening = thickening
         self.plot_state_interval = plot_state_interval
+        self.plot_state_variables = plot_state_variables
+        self.plot_state_variable_info_map = plot_state_variable_info_map
         self.iteration = iteration
         self.n_restart_interval = n_restart_interval
+        self.reacting = reacting
 
         # Initialize values which are passed in as None
         self.probes: list[Probe] = [] if probes is None else probes
@@ -132,6 +138,9 @@ def __init__(
             physics=self.physics,
         )
 
+        # Initialize CSV writers
+        self.csv_writers: list[CSVWriter] = []
+
         # Set up time integrators
         integrators: list[TimeIntegrator] = []
         advection = SSPRK3(self.inviscid_flux)
@@ -183,6 +192,7 @@ def __init__(
                 geometry=self.geometry,
                 physics=self.physics,
             )
+
             integrators += [ForwardEuler(self.boundary_layer)]
 
         if source_terms is not None:
@@ -313,6 +323,7 @@ def advance_simulation(self, tFinal: float, res_p_target: float = -1.0) -> None:
             # perform other updates
             self.update_probes(iters)
             self.update_XT_diagrams(iters)
+            self.update_csv_writers(iters)
             iters += 1
             res_p = float(np.linalg.norm(p_new - p_old))
             if self.verbose and iters % self.output_every == 0:
@@ -326,7 +337,9 @@ def advance_simulation(self, tFinal: float, res_p_target: float = -1.0) -> None:
             ):
                 plot_state(
                     self,
-                    f"figures/anim/test_{iters // self.plot_state_interval:05d}.png",
+                    f"./figures/anim/test_{iters // self.plot_state_interval:05d}.png",
+                    variable_info_map=self.plot_state_variable_info_map,
+                    plot_variables=self.plot_state_variables,
                 )
 
             # Periodically save the fluid state
@@ -337,3 +350,8 @@ def advance_simulation(self, tFinal: float, res_p_target: float = -1.0) -> None:
         self.iteration = iters
         if self.n_restart_interval > 0:
             self.state.save(groupname="stop")
+
+    def update_csv_writers(self, iters: int) -> None:
+        """Update all CSV writers to the current value."""
+        for csv_writer in self.csv_writers:
+            csv_writer.update(iters)

src/stanshock/models/boundary_layer.py

@@ -2,7 +2,11 @@
 
 import numpy as np
 
-from stanshock.models.wall_models import HeatFlux, SkinFriction, get_wall_state
+from stanshock.models.wall_models import (
+    HeatFlux,
+    SkinFriction,
+    get_wall_state,
+)
 from stanshock.physics.fluid_base import FluidState
 from stanshock.system.backend import Array, Index, Unpack
 from stanshock.system.base import PrecomputeStepName, PrecomputeSteps, RightHandSide

src/stanshock/processing/csv_writer.py

@@ -0,0 +1,83 @@
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+
+from stanshock.processing.plot import VariableInfo, get_variable_info_map
+
+
+class CSVWriter:
+    """
+    Writes simulation data to CSV file at specified intervals.
+    Creates a new numbered file for each output.
+    """
+
+    def __init__(
+        self,
+        combustor,
+        filename: str | Path,
+        interval: int = 100,
+        variables: list[str] | None = None,
+        variable_info_map: dict[str, VariableInfo] | None = None,
+    ) -> None:
+        """
+        Initialize CSV writer.
+
+        Args:
+            combustor: Combustor instance
+            filename: Base output CSV file path (will be numbered: filename_00000.csv, etc.)
+            interval: Write every 'interval' iterations (0 = never)
+            wall_temperature: Wall temperature for heat flux calculation
+        """
+        self.combustor = combustor
+        self.base_filename = Path(filename)
+
+        self.interval = interval
+        self.output_counter = 0
+        self.idx = combustor.geometry.idx_cells
+        self.x = combustor.geometry.xc[combustor.geometry.idx_cells]
+        if variables is None:
+            variables = ["x", "rho", "u", "p", "a", "T"]
+        self.headers = variables
+        if variable_info_map is None:
+            variable_info_map = get_variable_info_map(combustor.physics)
+        self.variable_info_map = variable_info_map
+        self.parent = self.base_filename.parent
+        self.stem = self.base_filename.stem
+        self.suffix = self.base_filename.suffix
+        self.fmt = ["%.4e"] + [variable_info_map[x].fmt for x in self.headers[1:]]
+
+    def update(self, iteration: int) -> None:
+        """Update CSV with current state if iteration matches interval."""
+        # Same condition as plot_state in combustor.py
+        if self.interval <= 0:
+            return
+        if iteration % self.interval != 0:
+            return
+
+        self.write_current_state()
+
+    def write_current_state(self) -> None:
+        """Write current state to a new numbered CSV file."""
+
+        combustor = self.combustor
+        state = combustor.state[self.idx]
+
+        state_matrix = np.column_stack(
+            (self.x, *[self.variable_info_map[v].fun(state) for v in self.headers[1:]])
+        )
+
+        filename = self.parent / f"{self.stem}_{self.output_counter:05d}{self.suffix}"
+        filename.parent.mkdir(parents=True, exist_ok=True)
+
+        np.savetxt(
+            filename,
+            state_matrix,
+            delimiter=",",
+            header=",".join(self.headers),
+            fmt=self.fmt,
+            comments="",
+        )
+
+        self.output_counter += 1

src/stanshock/processing/plot.py

@@ -27,6 +27,7 @@ class VariableInfo:
     plot_label: str
     fun: Callable[[FluidState], Array]
     scale: float = 1.0
+    fmt: str = "%.4e"
 
 
 VariableInfoMap: TypeAlias = dict[str, VariableInfo]
@@ -56,6 +57,11 @@ def get_variable_info_map(physics: FluidPhysics) -> VariableInfoMap:
         "gamma": VariableInfo(
             short_name="g", plot_label=r"$\gamma~[\mathrm{-}]$", fun=physics.get_gamma
         ),
+        "sound speed": VariableInfo(
+            short_name="a",
+            plot_label=r"$a~[\mathrm{m/s}]$",
+            fun=physics.get_sound_speed,
+        ),
         "mach": VariableInfo(
             short_name="m",
             plot_label=r"$Ma~[\mathrm{-}]$",
@@ -407,6 +413,9 @@ def plot_state(
 
     fig.suptitle(rf"$t = {domain.t * 1.0e3:.4f}$ ms")
 
+    output_path = Path(filename)
+    output_path.parent.mkdir(parents=True, exist_ok=True)
+
     fig.tight_layout()
-    fig.savefig(filename, bbox_inches="tight", dpi=300)
+    fig.savefig(output_path, bbox_inches="tight", dpi=300)
     plt.close()