From d1ca10ff38d6c01afcb79eb6891c4ba85e1c36fc Mon Sep 17 00:00:00 2001 From: Nadexterbrown Date: Thu, 28 Aug 2025 12:42:39 -0400 Subject: [PATCH 1/6] Added problem_setup.py for accurate mixture parameter extraction using cantera --- problem_setup.py | 268 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 268 insertions(+) create mode 100644 problem_setup.py diff --git a/problem_setup.py b/problem_setup.py new file mode 100644 index 0000000..14f5bf3 --- /dev/null +++ b/problem_setup.py @@ -0,0 +1,268 @@ +# problem_setup.py +# Minimal, robust problem setup for CLTORC-style runs. +# - Validates inputs +# - Computes stoichiometric O2 from fuel +# - Builds air with N2/O2 = 3.76 unless you override nitrogenAmount +# - Supports two φ conventions: scale fuel (default) or scale oxidizer +# - Provides a ready-to-use ct.Solution() with TPX set + +from __future__ import annotations +from dataclasses import dataclass, asdict +from typing import Dict, List, Optional, Tuple +import math + +try: + import cantera as ct +except Exception as e: + raise RuntimeError( + "Cantera is required. Install with `pip install cantera` " + "and ensure the mechanism file path is valid." + ) from e + + +# --- Fuel database: fuel -> required O2 moles per 1 mole fuel at stoich --- +# Using standard complete combustion to CO2 + H2O (no dissociation): +# CxHy + (x + y/4) O2 -> x CO2 + (y/2) H2O +# Feel free to extend as needed. +_FUEL_O2_STOICH = { + "H2": lambda: 0.5, # H2 + 0.5 O2 -> H2O + "CH4": lambda: 2.0, # CH4 + 2 O2 -> CO2 + 2 H2O + "C2H6": lambda: 3.5, # C2H6 + 3.5 O2 -> 2 CO2 + 3 H2O + "C3H8": lambda: 5.0, + "C4H10": lambda: 6.5, + "CO": lambda: 0.5, # CO + 0.5 O2 -> CO2 +} + +_AIR_N2_PER_O2 = 3.76 # dry air default + + +def _stoich_o2_for_fuel(fuel: str) -> float: + if fuel not in _FUEL_O2_STOICH: + raise ValueError( + f"Unknown/unsupported fuel '{fuel}'. " + f"Add it to _FUEL_O2_STOICH or provide your own logic." + ) + return float(_FUEL_O2_STOICH[fuel]()) + + +def _to_si_temperature(T: float, unit: str) -> float: + unit = unit.lower() + if unit in ("k", "kelvin"): + return float(T) + if unit in ("c", "celsius"): + return float(T) + 273.15 + if unit in ("f", "fahrenheit"): + return (float(T) - 32.0) * (5.0 / 9.0) + 273.15 + raise ValueError(f"Unsupported temperature unit '{unit}'") + + +def _to_si_pressure(P: float, unit: str) -> float: + unit = unit.lower() + if unit in ("pa", "pascal", "pascals"): + return float(P) + if unit in ("kpa",): + return float(P) * 1.0e3 + if unit in ("mpa",): + return float(P) * 1.0e6 + if unit in ("bar",): + return float(P) * 1.0e5 + if unit in ("atm", "atmosphere", "atmospheres"): + return float(P) * ct.one_atm + if unit in ("psi",): + return float(P) * 6894.757293168 + raise ValueError(f"Unsupported pressure unit '{unit}'") + + +@dataclass +class MyClass: + T: Optional[float] = None # [K] + P: Optional[float] = None # [Pa] + Phi: Optional[float] = None # equivalence ratio φ + Fuel: Optional[str] = None # e.g., "H2", "CH4", "C2H6", "C4H10", ... + mech: Optional[str] = None # Cantera mechanism file + species: Optional[List[str]] = None + oxygenAmount: Optional[float] = None # moles of O2 per 1 fuel at stoich + nitrogenAmount: Optional[float] = None # moles of N2 per stoich O2, unless overridden + X: Dict[str, float] = None + + def __post_init__(self): + if self.X is None: + self.X = {} + + # ---- composition/update helpers ---- + def update_composition( + self, + *, + phi_convention: str = "scale_fuel", + # "scale_fuel": n_fuel = φ, n_O2 = ν_O2(stoich), n_N2 = 3.76*ν_O2 or override + # "scale_oxidizer": n_fuel = 1.0, n_O2 = ν_O2 / φ, n_N2 scaled with O2 + use_dry_air: bool = True + ) -> None: + if self.Fuel is None: + raise ValueError("Fuel not set.") + if self.Phi is None or self.Phi <= 0.0: + raise ValueError("Phi must be > 0.") + # compute stoich O2 for 1 mole fuel + self.oxygenAmount = _stoich_o2_for_fuel(self.Fuel) + + # Set N2 baseline if not explicitly provided: + # Default: dry air with 3.76 N2 per O2 at stoich (unless nitrogenAmount was set) + n2_per_o2 = _AIR_N2_PER_O2 if use_dry_air else 0.0 + if self.nitrogenAmount is None or self.nitrogenAmount < 0: + self.nitrogenAmount = n2_per_o2 * self.oxygenAmount + + # Build mixture by φ convention + phi = float(self.Phi) + if phi_convention not in ("scale_fuel", "scale_oxidizer"): + raise ValueError("phi_convention must be 'scale_fuel' or 'scale_oxidizer'") + + if phi_convention == "scale_fuel": + n_fuel = phi + n_O2 = self.oxygenAmount + n_N2 = self.nitrogenAmount + else: # scale_oxidizer + n_fuel = 1.0 + n_O2 = self.oxygenAmount / phi + # keep N2 tied to O2 in air proportion + n_N2 = (self.nitrogenAmount / self.oxygenAmount) * n_O2 if self.oxygenAmount > 0 else 0.0 + + # Build the X dict; Cantera will normalize internally. + comp = {self.Fuel: n_fuel} + if n_O2 > 0.0: + comp["O2"] = n_O2 + if n_N2 > 0.0: + comp["N2"] = n_N2 + + self.X = comp + + def load_mechanism_species(self) -> None: + if not self.mech: + raise ValueError("Mechanism file not specified.") + gas = ct.Solution(self.mech) + self.species = list(gas.species_names) + # We intentionally don't keep the gas object here + del gas + + # ---- convenient builders ---- + def build_gas(self) -> ct.Solution: + if self.any_none(self.T, self.P, self.Phi, self.Fuel, self.mech, self.X): + raise RuntimeError("Missing required parameters before building gas.") + gas = ct.Solution(self.mech) + gas.TPX = self.T, self.P, self.X + return gas + + @staticmethod + def any_none(*args) -> bool: + return any(a is None for a in args) + + def summary(self) -> str: + lines = [ + "---- Problem Setup Summary ----", + f"T [K] : {self.T}", + f"P [Pa] : {self.P}", + f"Phi : {self.Phi}", + f"Fuel : {self.Fuel}", + f"mech : {self.mech}", + f"O2_stoich : {self.oxygenAmount}", + f"N2 amount : {self.nitrogenAmount}", + f"X (raw) : {self.X}", + "--------------------------------", + ] + return "\n".join(lines) + + +# Global params object, like you had +input_params = MyClass() + + +def initialize_parameters( + T: float, + P: float, + Phi: float, + Fuel: str, + mech: str, + *, + nitrogenAmount: Optional[float] = None, + temperature_unit: str = "K", + pressure_unit: str = "Pa", + phi_convention: str = "scale_fuel", + use_dry_air: bool = False, +) -> None: + """ + Initialize the global input_params with validated, SI-consistent values + and a chemically consistent composition for the requested φ. + + Parameters + ---------- + T : float + Temperature (units via temperature_unit) + P : float + Pressure (units via pressure_unit) + Phi : float + Equivalence ratio (>0) + Fuel : str + Fuel key present in _FUEL_O2_STOICH (extend as needed) + mech : str + Cantera mechanism file path/name (e.g., 'gri30.yaml') + nitrogenAmount : float or None + If provided, overrides default air-based N2 moles. Interpreted as + moles of N2 associated with stoich O2 for 1 mole fuel (before φ convention). + temperature_unit : {'K','C','F'} + pressure_unit : {'Pa','kPa','MPa','bar','atm','psi'} + phi_convention : {'scale_fuel','scale_oxidizer'} + - 'scale_fuel' : n_fuel = φ, n_O2 = ν_O2(stoich), N2 tied to stoich O2 + - 'scale_oxidizer' : n_fuel = 1, n_O2 = ν_O2/φ, N2 tied to adjusted O2 + use_dry_air : bool + If True, use N2/O2=3.76 unless nitrogenAmount explicitly provided. + """ + # Convert to SI + T_SI = _to_si_temperature(T, temperature_unit) + P_SI = _to_si_pressure(P, pressure_unit) + + # Quick sanity checks + if not (0.0 < T_SI < 1.0e4): + raise ValueError(f"Unreasonable temperature after conversion: {T_SI} K") + if not (1.0 <= P_SI < 1.0e9): + raise ValueError(f"Unreasonable pressure after conversion: {P_SI} Pa") + if Phi <= 0.0: + raise ValueError("Phi must be > 0.") + + # Populate + input_params.T = float(T_SI) + input_params.P = float(P_SI) + input_params.Phi = float(Phi) + input_params.Fuel = str(Fuel) + input_params.mech = str(mech) + input_params.nitrogenAmount = float(nitrogenAmount) if nitrogenAmount is not None else None + + # Chemistry pieces + input_params.update_composition(phi_convention=phi_convention, use_dry_air=use_dry_air) + input_params.load_mechanism_species() + + +# ---- Convenience getters (optional) ---- +def get_params_dict() -> Dict: + return asdict(input_params) + + +def build_gas() -> ct.Solution: + """Shortcut: build a Cantera gas from the global input_params.""" + return input_params.build_gas() + + +# ---- Example usage (remove or keep as a reference) ---- +if __name__ == "__main__": + # Example: φ=1.2 H2/air at 300 K and 1 atm using gri30 + initialize_parameters( + T=300.0, + P=1.0, + Phi=1.0, + Fuel="H2", + mech="./Li-Dryer-H2-mechanism.yaml", + temperature_unit="K", + pressure_unit="atm", # convenience; converted to Pa internally + phi_convention="scale_fuel", # or "scale_oxidizer" + ) + print(input_params.summary()) + gas = build_gas() + print(f"Built gas with {gas.T:.1f} K, {gas.P/ct.one_atm:.3f} atm, X={gas.X}") From 2373c53ee69b1dba82f8ee9c45f7d2bcd54fadfe Mon Sep 17 00:00:00 2001 From: Nadexterbrown Date: Thu, 28 Aug 2025 12:45:23 -0400 Subject: [PATCH 2/6] Added piston_setup.py for boundary condition selection between predetermined (constant, oscillating, etc.) and external file driven piston behavior. --- piston_setup.py | 331 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 331 insertions(+) create mode 100644 piston_setup.py diff --git a/piston_setup.py b/piston_setup.py new file mode 100644 index 0000000..820e532 --- /dev/null +++ b/piston_setup.py @@ -0,0 +1,331 @@ +import numpy as np +import pandas as pd +from scipy.interpolate import interp1d +from pathlib import Path +import re +from typing import Optional, Union, Dict, Any + +__all__ = ['simple_piston', 'datafile_piston', 'initialize_piston_mode', 'piston_velocity'] + +# Global variables for piston mode and data +_piston_mode = 'simple' +_piston_data_loader = None +_piston_params = None + + +class PistonDataLoader: + """ + Loads and interpolates piston velocity data from file. + """ + + def __init__(self, filepath: Union[str, Path], + time_column: str = None, + velocity_column: str = None, + interpolation_kind: str = 'linear'): + """ + Initialize the piston data loader. + + Parameters: + ----------- + filepath : str or Path + Path to the data file + time_column : str, optional + Name of the time column (auto-detected if None) + velocity_column : str, optional + Name of the velocity column (auto-detected if None) + interpolation_kind : str + Type of interpolation ('linear', 'cubic', etc.) + """ + self.filepath = Path(filepath) + self.interpolation_kind = interpolation_kind + self.interpolator = None + self.time_data = None + self.velocity_data = None + + # Load and process data + self._load_data(time_column, velocity_column) + + def _load_data_file(self, file_path): + """Load data file with custom header parsing.""" + # Read the file, extracting the second row as column names + with open(file_path, 'r') as f: + lines = f.readlines() + + # Extract the header row (second line) and clean column names + raw_headers = lines[1].strip() + + # Use regex to split headers by multiple spaces while keeping words together + headers = re.split(r'\s{2,}', raw_headers) # Split on 2+ spaces + headers[0] = headers[0].lstrip('#') # Remove the '#' from the first column name + + # Read the data, skipping the first two rows to exclude the initial comment and headers + return pd.read_csv(file_path, delim_whitespace=True, skiprows=2, names=headers) + + def _find_closest_columns(self, df, search_string, num_matches=1): + """ + Finds the closest matching column names in a DataFrame to a given search string. + + Parameters: + - df: pandas DataFrame. + - search_string: string to match against column headers. + - num_matches: number of closest matches to return. + + Returns: + - List of column names that best match the search string. + """ + import difflib + + columns = df.columns + matches = difflib.get_close_matches(search_string, columns, n=num_matches, cutoff=0.2) + return matches + + def _load_data(self, time_column=None, velocity_column=None): + """Load and prepare data for interpolation.""" + try: + # Load data using custom parser + data = self._load_data_file(self.filepath) + + # Auto-detect columns if not provided + if time_column is None: + time_matches = self._find_closest_columns(data, 'time', 1) + if not time_matches: + # Try other common time column names + for name in ['t', 'Time', 'TIME', 'time_s', 'time[s]']: + if name in data.columns: + time_column = name + break + if time_column is None: + raise ValueError("Could not auto-detect time column") + else: + time_column = time_matches[0] + + if velocity_column is None: + velocity_matches = self._find_closest_columns(data, 'velocity', 1) + if not velocity_matches: + # Try other common velocity column names + for name in ['flame_velocity', 'vel', 'Velocity', 'VELOCITY', 'u', 'speed']: + if name in data.columns: + velocity_column = name + break + if velocity_column is None: + raise ValueError("Could not auto-detect velocity column") + else: + velocity_column = velocity_matches[0] + + # Validate columns exist + if time_column not in data.columns: + raise ValueError(f"Time column '{time_column}' not found. Available: {list(data.columns)}") + if velocity_column not in data.columns: + raise ValueError(f"Velocity column '{velocity_column}' not found. Available: {list(data.columns)}") + + # Extract data + self.time_data = data[time_column].values + self.velocity_data = data[velocity_column].values + + # Clean data + mask = ~(np.isnan(self.time_data) | np.isnan(self.velocity_data)) + self.time_data = self.time_data[mask] + self.velocity_data = self.velocity_data[mask] + + # Sort by time + sort_idx = np.argsort(self.time_data) + self.time_data = self.time_data[sort_idx] + self.velocity_data = self.velocity_data[sort_idx] + + # Create interpolator + self.interpolator = interp1d( + self.time_data, + self.velocity_data, + kind=self.interpolation_kind, + bounds_error=False, + fill_value='extrapolate' + ) + + print(f"Piston data loaded successfully:") + print(f" File: {self.filepath}") + print(f" Points: {len(self.time_data)}") + print(f" Time range: [{self.time_data.min():.6f}, {self.time_data.max():.6f}]") + print(f" Velocity range: [{self.velocity_data.min():.6f}, {self.velocity_data.max():.6f}]") + print(f" Time column: '{time_column}'") + print(f" Velocity column: '{velocity_column}'") + + except Exception as e: + raise RuntimeError(f"Failed to load piston data from {self.filepath}: {str(e)}") + + def get_velocity(self, time: float) -> float: + """Get interpolated velocity at given time.""" + if self.interpolator is None: + raise RuntimeError("Data not loaded properly") + return float(self.interpolator(time)) + + +def simple_piston(t: float, piston_params: dict) -> float: + """ + Time-dependent piston velocity u_p(t) with an internal C^2 ramp. + + Parameters + ---------- + t : float + Current time [s]. + piston_params : dict + Dictionary containing: + kind : {"constant","sine"} + Type of motion. + U : float + Constant velocity [m/s] for kind="constant". + A : float + Amplitude [m/s] for kind="sine". + f : float + Frequency [Hz] for kind="sine". + ramp_time : float + Ramp-in duration [s] for smooth start (0 disables ramp). + U0 : float + Mean velocity [m/s] for kind="sine". + + Returns + ------- + float + Wall speed u_p(t) [m/s]. + """ + + kind = piston_params.get("kind", "constant") + U = piston_params.get("U", 1.0) + A = piston_params.get("A", 0.0) + f = piston_params.get("f", 0.0) + ramp_time = piston_params.get("ramp_time", 0.0) + + def ramp_smooth(t, ramp_time): + """ + Sigmoid/tanh ramp from 0 to 1 over [0, ramp_time]. + Centered at 0.5*ramp_time for symmetric rise. + """ + if ramp_time <= 0.0: + return 1.0 + arg = 6.0 * (t - 0.5 * ramp_time) / ramp_time # 95% rise within ramp_time + return 0.5 * (1.0 + np.tanh(arg)) + + if kind == "sine": + vl = U + A * np.sin(2.0 * np.pi * f * t) + else: # constant + vl = U + + if ramp_time == 0.0: + return vl + else: + return vl * ramp_smooth(t, ramp_time) + + +def datafile_piston(t: float, piston_params: dict = None) -> float: + """ + Get piston velocity from interpolated data file. + + Parameters + ---------- + t : float + Current time [s] + piston_params : dict, optional + Not used for datafile mode, kept for interface compatibility + + Returns + ------- + float + Interpolated wall speed u_p(t) [m/s] + """ + global _piston_data_loader + + if _piston_data_loader is None: + raise RuntimeError("Piston data not initialized. Call initialize_piston_mode() with mode='datafile' first.") + + return _piston_data_loader.get_velocity(t) + + +def initialize_piston_mode(mode: str, + piston_params: dict = None, + datafile_path: str = None, + time_column: str = None, + velocity_column: str = None, + interpolation_kind: str = 'linear') -> None: + """ + Initialize piston mode at simulation start. + + Parameters + ---------- + mode : {'simple', 'datafile'} + Piston velocity mode + piston_params : dict, optional + Parameters for simple mode + datafile_path : str, optional + Path to data file (required for datafile mode) + time_column : str, optional + Time column name (auto-detected if None) + velocity_column : str, optional + Velocity column name (auto-detected if None) + interpolation_kind : str, optional + Interpolation method for datafile mode + """ + global _piston_mode, _piston_data_loader, _piston_params + + if mode not in ['simple', 'datafile']: + raise ValueError("Mode must be 'simple' or 'datafile'") + + _piston_mode = mode + _piston_params = piston_params or {} + + if mode == 'datafile': + if datafile_path is None: + raise ValueError("datafile_path is required for datafile mode") + + print(f"Initializing piston in datafile mode...") + _piston_data_loader = PistonDataLoader( + filepath=datafile_path, + time_column=time_column, + velocity_column=velocity_column, + interpolation_kind=interpolation_kind + ) + else: + print(f"Initializing piston in simple mode with parameters: {_piston_params}") + _piston_data_loader = None + + +def piston_velocity(t: float, piston_params: dict = None) -> float: + """ + Main piston velocity function that switches between modes. + + This is the function you call from your solver. + + Parameters + ---------- + t : float + Current time [s] + piston_params : dict, optional + Parameters (used only in simple mode) + + Returns + ------- + float + Piston velocity [m/s] + """ + global _piston_mode, _piston_params + + if _piston_mode == 'simple': + # Use provided params or global params + params = piston_params if piston_params is not None else _piston_params + return simple_piston(t, params) + + elif _piston_mode == 'datafile': + return datafile_piston(t, piston_params) + + else: + raise RuntimeError(f"Unknown piston mode: {_piston_mode}. Call initialize_piston_mode() first.") + + +# Example usage +if __name__ == "__main__": + print("Piston velocity functions initialized.") + print("\nUsage examples:") + print("1. Simple mode:") + print(" initialize_piston_mode('simple', {'kind': 'sine', 'U': 10, 'A': 5, 'f': 2})") + print(" velocity = piston_velocity(0.5)") + print("\n2. Datafile mode:") + print(" initialize_piston_mode('datafile', datafile_path='piston_data.dat')") + print(" velocity = piston_velocity(0.5)") \ No newline at end of file From be0d8e94a527a3e55cce5c58ea140b21a61db37f Mon Sep 17 00:00:00 2001 From: Nadexterbrown Date: Thu, 28 Aug 2025 12:45:40 -0400 Subject: [PATCH 3/6] Constant velocity piston example --- constant_velocity_piston.py | 384 ++++++++++++++++++++++++++++++++++++ 1 file changed, 384 insertions(+) create mode 100644 constant_velocity_piston.py diff --git a/constant_velocity_piston.py b/constant_velocity_piston.py new file mode 100644 index 0000000..b4affc7 --- /dev/null +++ b/constant_velocity_piston.py @@ -0,0 +1,384 @@ +import numpy as np +import cantera as ct +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm +from typing import Tuple, List, Optional + +from problem_setup import initialize_parameters, input_params # your dataclass + +from lagrangian import ( + PlanarGeometry, + LagrangianMesh, + LagrangianSolver, + SimulationParameters +) + + +def run_piston_problem( + n_cells: int = 100, + target_velocity: float = 1.0, + ramp_rate: Optional[float] = None, + t_end: float = 0.5, + cfl: float = 0.3, + output_frequency: int = 10, + *, + # ---- pass ONE of these sources of thermodynamic state ---- + gas: "ct.Solution | None" = None, + params: "MyClass | None" = None, + uniform_state: "tuple[float, float] | None" = None, # (rho0 [kg/m^3], p0 [Pa]) + override_gamma: Optional[float] = None, + initial: str = "uniform" # or "riemann_cltorc" +) -> Tuple[np.ndarray, List[np.ndarray]]: + """ + Run a piston-in-gas simulation with thermodynamics supplied by caller. + No internal initialization of input_params happens here. + + Provide exactly one of: + - gas: a Cantera Solution with TPX set + - params: your MyClass object; we call params.build_gas() + - uniform_state: (rho0, p0) in SI + + If override_gamma is None, gamma = cp/cv from the gas (when available). + """ + + # --------------------------- + # Resolve thermodynamic state + # --------------------------- + provided = sum(x is not None for x in (gas, params, uniform_state)) + if provided != 1: + raise ValueError("Provide exactly ONE of gas=, params=, or uniform_state=.") + + if params is not None: + if not hasattr(params, "build_gas"): + raise TypeError("params must be a MyClass with .build_gas()") + if ct is None: + raise RuntimeError("Cantera required when using params= or gas=.") + gas = params.build_gas() + + if gas is not None: + p0 = float(gas.P) # Pa + rho0 = float(gas.density) # kg/m^3 + gamma = float(override_gamma) if override_gamma is not None else float(gas.cp_mass / gas.cv_mass) + else: + # uniform_state provided + rho0, p0 = map(float, uniform_state) + if override_gamma is None: + raise ValueError("uniform_state provided but override_gamma is None. Specify gamma.") + gamma = float(override_gamma) + + if not (gamma > 1.0): + raise ValueError(f"gamma must be > 1. Got {gamma}.") + + # specific internal energy for a gamma-law closure + e0 = p0 / ((gamma - 1.0) * rho0) + + # --------------------------- + # Mesh / geometry + # --------------------------- + geometry = PlanarGeometry() + x_max = 1.0 + r = np.linspace(0.0, x_max, n_cells + 1) # interface positions + + # --------------------------- + # Initial conditions + # --------------------------- + rho = np.empty(n_cells) + u = np.empty(n_cells + 1) + e = np.empty(n_cells) + + rho[:] = rho0 + u[:] = 0.0 + e[:] = e0 + + # Create mesh and solver + mesh = LagrangianMesh(r, rho, u, e, geometry, gamma) + + params = SimulationParameters( + gamma=gamma, + cfl=cfl, + t_end=t_end, + output_frequency=output_frequency, + use_artificial_viscosity=True + ) + + # Create a custom solver that overrides the boundary conditions + class PistonSolver(LagrangianSolver): + def apply_boundary_conditions(self) -> None: + """Apply piston boundary condition on the left and fixed wall on the right.""" + # Left boundary (piston) + if ramp_rate is None: + # Constant velocity (no ramp) + self.mesh.u[0] = target_velocity + else: + # Ramp up velocity based on current time + current_velocity = min(self.time * ramp_rate, target_velocity) + self.mesh.u[0] = current_velocity + + # Right boundary (fixed wall) + self.mesh.u[-1] = 0.0 + + solver = PistonSolver(mesh, params) + + # Solve + return solver.solve() + + +def plot_xt_diagram( + times: np.ndarray, + solutions: List[np.ndarray], + variable: str = 'density', + output_file: Optional[str] = None +) -> None: + """ + Plot an x-t diagram of the solution. + + Args: + times: Array of simulation times + solutions: List of solution snapshots (r, u, rho, e, p) + variable: Which variable to plot ('density', 'velocity', 'pressure', 'energy') + output_file: Optional file to save plot to + """ + # Extract data for all time steps + n_times = len(times) + + # Get variable index + var_idx = {'density': 2, 'velocity': 1, 'pressure': 4, 'energy': 3}[variable] + var_name = {'density': 'Density', 'velocity': 'Velocity', + 'pressure': 'Pressure', 'energy': 'Internal Energy'}[variable] + + # Create arrays for x, t, and the variable + x_values = [] + t_values = [] + var_values = [] + + for i, time in enumerate(times): + r, u, rho, e, p = solutions[i] + + # For cell-centered quantities (density, pressure, energy) + if var_idx in [2, 3, 4]: + x = 0.5 * (r[:-1] + r[1:]) # Cell centers + var = solutions[i][var_idx] + # For interface quantities (velocity) + else: + x = r # Interface positions + var = solutions[i][var_idx] + + # Add to arrays + for j in range(len(x)): + x_values.append(x[j]) + t_values.append(time) + var_values.append(var[j]) + + # Convert to numpy arrays + x_values = np.array(x_values) + t_values = np.array(t_values) + var_values = np.array(var_values) + + # Create figure + plt.figure(figsize=(10, 8)) + + # Create scatter plot with color representing the variable + if variable == 'density': + # Use log scale for density to better visualize compression + scatter = plt.scatter(x_values, t_values, c=var_values, cmap='viridis', + s=10, alpha=0.8, norm=LogNorm()) + else: + scatter = plt.scatter(x_values, t_values, c=var_values, cmap='viridis', + s=10, alpha=0.8) + + # Add colorbar + cbar = plt.colorbar(scatter) + cbar.set_label(var_name) + + # Set labels and title + plt.xlabel('Position (x)') + plt.ylabel('Time (t)') + plt.title(f'X-T Diagram: {var_name} Evolution for Piston Problem') + + # Save or show + if output_file: + plt.savefig(output_file, dpi=300, bbox_inches='tight') + print(f"Plot saved to {output_file}") + + plt.tight_layout() + plt.show() + + +def plot_final_state( + times: np.ndarray, + solutions: List[np.ndarray], + output_file: Optional[str] = None +) -> None: + """ + Plot the final state of the simulation. + + Args: + times: Array of simulation times + solutions: List of solution snapshots (r, u, rho, e, p) + output_file: Optional file to save plot to + """ + # Extract final solution + r_final, u_final, rho_final, e_final, p_final = solutions[-1] + + # Create figure with subplots + fig, axs = plt.subplots(2, 2, figsize=(12, 10)) + + # Plot density + axs[0, 0].plot(0.5 * (r_final[:-1] + r_final[1:]), rho_final, 'b-', linewidth=2) + axs[0, 0].set_title('Density') + axs[0, 0].set_ylabel('Density') + axs[0, 0].grid(True) + + # Plot velocity + axs[0, 1].plot(r_final, u_final, 'r-', linewidth=2) + axs[0, 1].set_title('Velocity') + axs[0, 1].set_ylabel('Velocity') + axs[0, 1].grid(True) + + # Plot pressure + axs[1, 0].plot(0.5 * (r_final[:-1] + r_final[1:]), p_final, 'g-', linewidth=2) + axs[1, 0].set_title('Pressure') + axs[1, 0].set_xlabel('Position') + axs[1, 0].set_ylabel('Pressure') + axs[1, 0].grid(True) + + # Plot internal energy + axs[1, 1].plot(0.5 * (r_final[:-1] + r_final[1:]), e_final, 'm-', linewidth=2) + axs[1, 1].set_title('Specific Internal Energy') + axs[1, 1].set_xlabel('Position') + axs[1, 1].set_ylabel('Energy') + axs[1, 1].grid(True) + + plt.tight_layout() + + # Add main title + plt.suptitle(f'Piston Problem - Final State at t = {times[-1]:.3f}', fontsize=16) + plt.subplots_adjust(top=0.93) + + # Save or show + if output_file: + plt.savefig(output_file, dpi=300, bbox_inches='tight') + print(f"Plot saved to {output_file}") + + plt.show() + + +def create_animation_frames( + times: np.ndarray, + solutions: List[np.ndarray], + output_dir: str = 'frames', + variable: str = 'density' +) -> None: + """ + Create animation frames of the solution evolution. + + Args: + times: Array of simulation times + solutions: List of solution snapshots (r, u, rho, e, p) + output_dir: Directory to save frames + variable: Which variable to plot ('density', 'velocity', 'pressure', 'energy') + """ + import os + + # Create output directory if it doesn't exist + os.makedirs(output_dir, exist_ok=True) + + # Get variable index + var_idx = {'density': 2, 'velocity': 1, 'pressure': 4, 'energy': 3}[variable] + var_name = {'density': 'Density', 'velocity': 'Velocity', + 'pressure': 'Pressure', 'energy': 'Internal Energy'}[variable] + + # Initialize y-axis limits + y_min = float('inf') + y_max = float('-inf') + + # Create frames + for i, time in enumerate(times): + r, u, rho, e, p = solutions[i] + + # For cell-centered quantities (density, pressure, energy) + if var_idx in [2, 3, 4]: + x = 0.5 * (r[:-1] + r[1:]) # Cell centers + var = solutions[i][var_idx] + # For interface quantities (velocity) + else: + x = r # Interface positions + var = solutions[i][var_idx] + + # Create figure + plt.figure(figsize=(10, 6)) + plt.plot(x, var, 'b-', linewidth=2) + plt.title(f'{var_name} at t = {time:.3f}') + plt.xlabel('Position') + plt.ylabel(var_name) + plt.grid(True) + + # Set consistent y-axis limits + if i == 0: + y_min = min(var) + y_max = max(var) + else: + y_min = min(y_min, min(var)) + y_max = max(y_max, max(var)) + + plt.ylim(y_min, y_max * 1.1) + + # Save frame + plt.savefig(f'{output_dir}/frame_{i:04d}.png', dpi=200, bbox_inches='tight') + plt.close() + + print(f"Created {len(times)} animation frames in {output_dir}") + + +if __name__ == "__main__": + # Gas Parameters + initialize_parameters( + T=293.0, + P=1.0, + Phi=1.0, + Fuel="H2", + # nitrogenAmount=0.5 * 3.76, + mech="./chemical_mechanism_files/Li-Dryer-H2-mechanism.yaml", + temperature_unit="K", + pressure_unit="atm", # convenience; converted to Pa internally + phi_convention="scale_fuel", # or "scale_oxidizer" + use_dry_air=True # Use dry air for nitrogen/oxygen ratio + ) + print(input_params.summary()) + gas = input_params.build_gas() + print(f"Built gas with {gas.T:.1f} K, {gas.P / ct.one_atm:.3f} atm, X={gas.X}") + + # Simulation parameters + n_cells = 1000 + target_velocity = 1500.0 # Positive means moving into the gas + ramp_rate = 0.0 # Velocity units per time unit (reaches target_velocity at t=0.5) + gamma = 1.4 + t_end = 0.0005 + cfl = 0.1 + output_frequency = 5 # Save more frequently for better visualization + + print(f"Running piston problem simulation with target velocity = {target_velocity}") + print(f"Velocity ramp rate = {ramp_rate} (velocity units per time unit)") + + # Run simulation + times, solutions = run_piston_problem( + n_cells=n_cells, + target_velocity=target_velocity, + t_end=t_end, + cfl=cfl, + output_frequency=output_frequency, + gas=gas + ) + + print(f"Simulation completed.") + print(f"Number of time steps: {len(times)}") + print(f"Final time: {times[-1]:.6f}") + + # Plot x-t diagram for density + plot_xt_diagram(times, solutions, variable='density') + + # Plot final state + plot_final_state(times, solutions) + + # Optionally create animation frames + # create_animation_frames(times, solutions, output_dir='frames', variable='density') \ No newline at end of file From ae9c16631039520bc11f9a32153519b5b913c969 Mon Sep 17 00:00:00 2001 From: Nadexterbrown Date: Thu, 28 Aug 2025 12:46:14 -0400 Subject: [PATCH 4/6] Preset oscillating (sine) velocity piston example. --- oscillating_velocity_piston.py | 620 +++++++++++++++++++++++++++++++++ 1 file changed, 620 insertions(+) create mode 100644 oscillating_velocity_piston.py diff --git a/oscillating_velocity_piston.py b/oscillating_velocity_piston.py new file mode 100644 index 0000000..23fbc24 --- /dev/null +++ b/oscillating_velocity_piston.py @@ -0,0 +1,620 @@ +import numpy as np +import cantera as ct +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm +from typing import Tuple, List, Optional, Union, Dict + +# Updated import to use the new piston functions file +from piston_setup import simple_piston, datafile_piston, initialize_piston_mode, piston_velocity +from problem_setup import initialize_parameters, input_params # your dataclass + +from lagrangian import ( + PlanarGeometry, + LagrangianMesh, + LagrangianSolver, + SimulationParameters +) + + +######################################################################################################################## +# Piston Problem Simulation +######################################################################################################################## + +def run_piston_problem( + n_cells: int = 100, + t_end: float = 0.5, + cfl: float = 0.3, + output_frequency: int = 10, + *, + # ---- pass ONE of these sources of thermodynamic state ---- + gas: "ct.Solution | None" = None, + params: "MyClass | None" = None, + uniform_state: "tuple[float, float] | None" = None, # (rho0 [kg/m^3], p0 [Pa]) + override_gamma: Optional[float] = None, + initial: str = "uniform", # or "riemann_cltorc" + piston_enabled: bool = True # New parameter to enable/disable piston +) -> Tuple[np.ndarray, List[np.ndarray]]: + """ + Run a piston-in-gas simulation with thermodynamics supplied by caller. + No internal initialization of input_params happens here. + + Provide exactly one of: + - gas: a Cantera Solution with TPX set + - params: your MyClass object; we call params.build_gas() + - uniform_state: (rho0, p0) in SI + + If override_gamma is None, gamma = cp/cv from the gas (when available). + """ + + # --------------------------- + # Resolve thermodynamic state + # --------------------------- + provided = sum(x is not None for x in (gas, params, uniform_state)) + if provided != 1: + raise ValueError("Provide exactly ONE of gas=, params=, or uniform_state=.") + + if params is not None: + if not hasattr(params, "build_gas"): + raise TypeError("params must be a MyClass with .build_gas()") + if ct is None: + raise RuntimeError("Cantera required when using params= or gas=.") + gas = params.build_gas() + + if gas is not None: + p0 = float(gas.P) # Pa + rho0 = float(gas.density) # kg/m^3 + gamma = float(override_gamma) if override_gamma is not None else float(gas.cp_mass / gas.cv_mass) + else: + # uniform_state provided + rho0, p0 = map(float, uniform_state) + if override_gamma is None: + raise ValueError("uniform_state provided but override_gamma is None. Specify gamma.") + gamma = float(override_gamma) + + if not (gamma > 1.0): + raise ValueError(f"gamma must be > 1. Got {gamma}.") + + # specific internal energy for a gamma-law closure + e0 = p0 / ((gamma - 1.0) * rho0) + + # --------------------------- + # Mesh / geometry + # --------------------------- + geometry = PlanarGeometry() + x_max = 2.5 # domain length + r = np.linspace(0.0, x_max, n_cells + 1) # interface positions + + # --------------------------- + # Initial conditions + # --------------------------- + rho = np.empty(n_cells) + u = np.empty(n_cells + 1) + e = np.empty(n_cells) + + rho[:] = rho0 + u[:] = 0.0 + e[:] = e0 + + # Create mesh and solver + mesh = LagrangianMesh(r, rho, u, e, geometry, gamma) + + sim_params = SimulationParameters( + gamma=gamma, + cfl=cfl, + t_end=t_end, + output_frequency=output_frequency, + use_artificial_viscosity=True + ) + + # Create a custom solver that overrides the boundary conditions + class PistonSolver(LagrangianSolver): + def __init__(self, mesh, params, piston_enabled=True): + super().__init__(mesh, params) + self.piston_enabled = piston_enabled + + def apply_boundary_conditions(self) -> None: + """Apply piston boundary condition on the left and fixed wall on the right.""" + # Left boundary (piston) + if not self.piston_enabled: + self.mesh.u[0] = 0.0 # Fixed wall if piston disabled + else: + # Use the global piston_velocity function + # No need to pass piston_params since it's handled globally + self.mesh.u[0] = piston_velocity(self.time) + + # Right boundary (fixed wall) + self.mesh.u[-1] = 0.0 + + solver = PistonSolver(mesh, sim_params, piston_enabled=piston_enabled) + + # Solve + return solver.solve() + + +######################################################################################################################## +# Visualization Functions +######################################################################################################################## + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm +from typing import List, Optional, Union, Dict +import cantera as ct + + +def plot_xt_diagram( + times: np.ndarray, + solutions: List[np.ndarray], # (r, u, rho, e, p) + variable: str = 'density', + output_file: Optional[str] = None, + *, + gas: Optional[ct.Solution] = None, + mech: Optional[str] = None, + X: Optional[Union[str, Dict[str, float]]] = None, + log_for: Optional[set] = None, + var_limit: Optional[Tuple[float, float]] = None, + coord: str = 'x', # 'x' or 'phi' + return_phi: bool = False +) -> Optional[List[tuple]]: + """ + Plot an x–t diagram. If coord='phi', use the mass-weighted Lagrangian + coordinate phi(x,t)=∫ rho dx on the horizontal axis. + + solutions[i] is (r, u, rho, e, p) with shapes: + r: (Ni+1,) interfaces + u: (Ni+1,) interface velocity + rho: (Ni,) cell-centered density + e: (Ni,) cell-centered internal energy + p: (Ni,) cell-centered pressure + """ + + # ----------------------- + # Setup gas / composition (only if needed) + # ----------------------- + native = {'density': 2, 'velocity': 1, 'pressure': 4, 'energy': 3} + derived_requested = variable not in native + + if gas is None and derived_requested: + if mech is None: + try: + mech = input_params.mech + except Exception: + raise ValueError("Need `mech` (or a ready `gas`) to reconstruct state with Cantera.") + if X is None: + try: + X = input_params.X + except Exception: + raise ValueError("Need mixture composition `X` (or set input_params.X).") + gas = ct.Solution(mech) + X_str = ",".join([f"{k}:{v}" for k, v in X.items()]) if isinstance(X, dict) else X + else: + X_str = gas.X.copy() + + # Pretty names + pretty = { + 'density': 'Density', 'velocity': 'Velocity', 'pressure': 'Pressure', 'energy': 'Internal Energy', + 'temperature': 'Temperature', 'cp': 'cp', 'cv': 'cv', 'gamma': 'Gamma', 'R': 'Gas constant (mix)', + 'sound_speed': 'Sound speed', 'mu': 'Viscosity', 'lambda': 'Thermal conductivity', + 'enthalpy': 'Specific enthalpy', 'entropy': 'Specific entropy', 'Mach': 'Mach number' + } + var_name = pretty.get(variable, variable) + + # Defaults for log scaling + default_log_set = {'density', 'pressure', 'cp', 'cv', 'mu', 'lambda'} + if log_for is None: + log_for = default_log_set + + # Species request parsing + def _parse_species_token(vname: str): + if vname.startswith('Y(') and vname.endswith(')'): return ('Y', vname[2:-1]) + if vname.startswith('X(') and vname.endswith(')'): return ('X', vname[2:-1]) + if vname.startswith('Y_'): return ('Y', vname[2:]) + if vname.startswith('X_'): return ('X', vname[2:]) + return (None, None) + + req_kind, req_species = _parse_species_token(variable) + + # ----------------------- + # φ helper + # ----------------------- + def _compute_phi( + r_iface: np.ndarray, + rho_cell: np.ndarray, + x_piston: float | None = None, + *, x_units: str = "m", rho_units: str = "kg/m^3" + ): + xu = {"m": 1.0, "cm": 1e-2, "mm": 1e-3}[x_units] + ru = {"kg/m^3": 1.0, "g/cm^3": 1000.0}[rho_units] + + r = r_iface.astype(float) * xu + rho = np.maximum(rho_cell, 0.0).astype(float) * ru + + dx = r[1:] - r[:-1] # (Ni,) + m = rho * dx # mass/area per cell (kg/m^2) + + if x_piston is None: + x_piston = r[0] + else: + x_piston = float(x_piston) * xu + + idx = np.searchsorted(r, x_piston, side="right") - 1 + idx = max(0, min(idx, len(dx) - 1)) + + m0 = rho[idx] * max(r[idx + 1] - x_piston, 0.0) + + phi_iface = np.zeros_like(r) + if idx + 1 < len(r): + cs = np.cumsum(m[idx + 1:]) # length Ni - idx - 1 + phi_iface[idx + 1:] = m0 + np.concatenate(([0.0], cs)) # <-- fixed + + phi_cell = 0.5 * (phi_iface[:-1] + phi_iface[1:]) + return phi_iface, phi_cell + + # ----------------------- + # Accumulators + # ----------------------- + x_values = [] + t_values = [] + var_values = [] + phi_list = [] # Will store (phi_iface, phi_cell) per time step if requested + + # ----------------------- + # Main loop + # ----------------------- + for i, t in enumerate(times): + r, ufaces, rho, e, p = solutions[i] + Ni = rho.shape[0] + + # Compute φ if needed (or if caller wants it returned) + if coord == 'phi' or return_phi: + phi_iface, phi_cell = _compute_phi(r, rho) + if return_phi: + phi_list.append((phi_iface.copy(), phi_cell.copy())) + + # Determine alignment and base x-grid + if variable in native: + if variable == 'velocity': + # interface quantity + x_plot = phi_iface if coord == 'phi' else r + vals = ufaces + for j in range(len(x_plot)): + x_values.append(x_plot[j]) + t_values.append(t) + var_values.append(vals[j]) + else: + # cell-centered + x_centers = phi_cell if coord == 'phi' else 0.5 * (r[:-1] + r[1:]) + if variable == 'density': + vals = rho + elif variable == 'pressure': + vals = p + elif variable == 'energy': + vals = e + for j in range(Ni): + x_values.append(x_centers[j]) + t_values.append(t) + var_values.append(vals[j]) + + else: + # Derived thermodynamic variable → cell-centered + x_centers = phi_cell if coord == 'phi' else 0.5 * (r[:-1] + r[1:]) + # For Mach we need a cell-centered u + u_cell = 0.5 * (ufaces[:-1] + ufaces[1:]) if variable == 'Mach' else None + + vals = np.empty(Ni, dtype=float) + for j in range(Ni): + gas.DPX = rho[j], p[j], X_str + if req_kind == 'Y': + k = gas.species_index(req_species); + vals[j] = gas.Y[k] + elif req_kind == 'X': + k = gas.species_index(req_species); + vals[j] = gas.X[k] + elif variable == 'temperature': + vals[j] = gas.T + elif variable == 'cp': + vals[j] = gas.cp_mass + elif variable == 'cv': + vals[j] = gas.cv_mass + elif variable == 'gamma': + vals[j] = gas.cp_mass / gas.cv_mass + elif variable == 'R': + vals[j] = gas.gas_constant / gas.mean_molecular_weight + elif variable == 'sound_speed': + vals[j] = gas.sound_speed + elif variable == 'mu': + vals[j] = gas.viscosity + elif variable == 'lambda': + vals[j] = gas.thermal_conductivity + elif variable == 'enthalpy': + vals[j] = gas.h_mass + elif variable == 'entropy': + vals[j] = gas.s_mass + elif variable == 'Mach': + a = gas.sound_speed + vals[j] = abs(u_cell[j]) / max(a, 1e-16) + else: + raise ValueError( + f"Unknown derived variable '{variable}'. " + f"Supported: temperature, cp, cv, gamma, R, sound_speed, mu, lambda, " + f"enthalpy, entropy, Mach, and species via Y(NAME)/X(NAME)/Y_NAME/X_NAME." + ) + + for j in range(Ni): + x_values.append(x_centers[j]) + t_values.append(t) + var_values.append(vals[j]) + + # ----------------------- + # Plot + # ----------------------- + x_values = np.asarray(x_values) + t_values = np.asarray(t_values) + var_values = np.asarray(var_values) + + plt.figure(figsize=(10, 8)) + use_log = variable in log_for and np.all(var_values > 0.0) + if use_log: + sc = plt.scatter(x_values, t_values, c=var_values, cmap='rainbow', + s=10, alpha=0.8, norm=LogNorm()) + else: + if var_limit is None: + var_limit = (np.min(var_values), np.max(var_values)) + sc = plt.scatter(x_values, t_values, c=var_values, cmap='rainbow', + s=10, alpha=0.8, vmin=var_limit[0], vmax=var_limit[1]) + + cbar = plt.colorbar(sc) + cbar.set_label(var_name) + + if coord == 'phi': + plt.xlabel(r'$\phi$ (mass per unit area)') + title_axis = 'φ' + else: + plt.xlabel('Position (x)') + title_axis = 'x' + + plt.ylabel('Time (t)') + plt.title(f'{title_axis}–T Diagram: {var_name}') + + if output_file: + plt.savefig(output_file, dpi=300, bbox_inches='tight') + print(f"Plot saved to {output_file}") + + plt.tight_layout() + plt.show() + + if return_phi: + return phi_list + return None + + +def plot_final_state( + times: np.ndarray, + solutions: List[np.ndarray], + output_file: Optional[str] = None +) -> None: + """ + Plot the final state of the simulation. + + Args: + times: Array of simulation times + solutions: List of solution snapshots (r, u, rho, e, p) + output_file: Optional file to save plot to + """ + # Extract final solution + r_final, u_final, rho_final, e_final, p_final = solutions[-1] + + # Create figure with subplots + fig, axs = plt.subplots(2, 2, figsize=(12, 10)) + + # Plot density + axs[0, 0].plot(0.5 * (r_final[:-1] + r_final[1:]), rho_final, 'b-', linewidth=2) + axs[0, 0].set_title('Density') + axs[0, 0].set_ylabel('Density') + axs[0, 0].grid(True) + + # Plot velocity + axs[0, 1].plot(r_final, u_final, 'r-', linewidth=2) + axs[0, 1].set_title('Velocity') + axs[0, 1].set_ylabel('Velocity') + axs[0, 1].grid(True) + + # Plot pressure + axs[1, 0].plot(0.5 * (r_final[:-1] + r_final[1:]), p_final, 'g-', linewidth=2) + axs[1, 0].set_title('Pressure') + axs[1, 0].set_xlabel('Position') + axs[1, 0].set_ylabel('Pressure') + axs[1, 0].grid(True) + + # Plot internal energy + axs[1, 1].plot(0.5 * (r_final[:-1] + r_final[1:]), e_final, 'm-', linewidth=2) + axs[1, 1].set_title('Specific Internal Energy') + axs[1, 1].set_xlabel('Position') + axs[1, 1].set_ylabel('Energy') + axs[1, 1].grid(True) + + plt.tight_layout() + + # Add main title + plt.suptitle(f'Piston Problem - Final State at t = {times[-1]:.3f}', fontsize=16) + plt.subplots_adjust(top=0.93) + + # Save or show + if output_file: + plt.savefig(output_file, dpi=300, bbox_inches='tight') + print(f"Plot saved to {output_file}") + + plt.show() + + +def create_animation_frames( + times: np.ndarray, + solutions: List[np.ndarray], + output_dir: str = 'frames', + variable: str = 'density' +) -> None: + """ + Create animation frames of the solution evolution. + + Args: + times: Array of simulation times + solutions: List of solution snapshots (r, u, rho, e, p) + output_dir: Directory to save frames + variable: Which variable to plot ('density', 'velocity', 'pressure', 'energy') + """ + import os + + # Create output directory if it doesn't exist + os.makedirs(output_dir, exist_ok=True) + + # Get variable index + var_idx = {'density': 2, 'velocity': 1, 'pressure': 4, 'energy': 3}[variable] + var_name = {'density': 'Density', 'velocity': 'Velocity', + 'pressure': 'Pressure', 'energy': 'Internal Energy'}[variable] + + # Initialize y-axis limits + y_min = float('inf') + y_max = float('-inf') + + # Create frames + for i, time in enumerate(times): + r, u, rho, e, p = solutions[i] + + # For cell-centered quantities (density, pressure, energy) + if var_idx in [2, 3, 4]: + x = 0.5 * (r[:-1] + r[1:]) # Cell centers + var = solutions[i][var_idx] + # For interface quantities (velocity) + else: + x = r # Interface positions + var = solutions[i][var_idx] + + # Create figure + plt.figure(figsize=(10, 6)) + plt.plot(x, var, 'b-', linewidth=2) + plt.title(f'{var_name} at t = {time:.3f}') + plt.xlabel('Position') + plt.ylabel(var_name) + plt.grid(True) + + # Set consistent y-axis limits + if i == 0: + y_min = min(var) + y_max = max(var) + else: + y_min = min(y_min, min(var)) + y_max = max(y_max, max(var)) + + plt.ylim(y_min, y_max * 1.1) + + # Save frame + plt.savefig(f'{output_dir}/frame_{i:04d}.png', dpi=200, bbox_inches='tight') + plt.close() + + print(f"Created {len(times)} animation frames in {output_dir}") + + +def plot_piston_velocity_vs_time(times: np.ndarray, output_file: Optional[str] = None) -> None: + """ + Plot the piston velocity as a function of time to visualize the boundary condition. + + Args: + times: Array of simulation times + output_file: Optional file to save plot to + """ + # Calculate piston velocity at each time step + piston_velocities = [piston_velocity(t) for t in times] + + plt.figure(figsize=(10, 6)) + plt.plot(times, piston_velocities, 'b-', linewidth=2) + plt.title('Piston Velocity vs Time') + plt.xlabel('Time [s]') + plt.ylabel('Piston Velocity [m/s]') + plt.grid(True) + + if output_file: + plt.savefig(output_file, dpi=300, bbox_inches='tight') + print(f"Piston velocity plot saved to {output_file}") + + plt.show() + + +######################################################################################################################## +# Main Execution +######################################################################################################################## + +if __name__ == "__main__": + # Gas Parameters + initialize_parameters( + T=503.0, + P=10.0e5, + Phi=1.0, + Fuel="H2", + # nitrogenAmount=0.5 * 3.76, + mech="./chemical_mechanism_files/Li-Dryer-H2-mechanism.yaml", + temperature_unit="K", + pressure_unit="Pa", # convenience; converted to Pa internally + phi_convention="scale_fuel", # or "scale_oxidizer" + use_dry_air=False # Use dry air for nitrogen/oxygen ratio + ) + print(input_params.summary()) + gas = input_params.build_gas() + print(f"Built gas with {gas.T:.1f} K, {gas.P / ct.one_atm:.3f} atm, X={gas.X}") + + # Simulation parameters + n_cells = 1000 + gamma = 1.4 + t_end = 0.0025537811000000002 + cfl = 0.3 + output_frequency = 5 # Save more frequently for better visualization + + # ======================================== + # Choose your piston mode here! + # ======================================== + + # Option 1: Simple piston (analytical functions) + print("Initializing piston in SIMPLE mode...") + initialize_piston_mode( + mode='simple', + piston_params={ + 'kind': 'sine', # or 'constant' + 'U': 1626.35, # mean velocity [m/s] + 'A': 0.2 * 1626.35, # amplitude [m/s] + 'f': 45.4e3, # frequency [Hz] - high frequency for short simulation time + 'ramp_time': 0.0 # smooth startup [s] - short ramp for short simulation + } + ) + + # Run simulation + print("Starting simulation...") + times, solutions = run_piston_problem( + n_cells=n_cells, + t_end=t_end, + cfl=cfl, + output_frequency=output_frequency, + gas=gas, + piston_enabled=True # Set to False to disable piston (fixed wall instead) + ) + + print(f"Simulation completed.") + print(f"Number of time steps: {len(times)}") + print(f"Final time: {times[-1]:.6f} s") + + # Plot piston velocity vs time to verify boundary condition + plot_piston_velocity_vs_time(times) + + # Plot x-t diagram for temperature + plot_xt_diagram(times, solutions, variable='temperature', gas=gas) + plot_xt_diagram(times, solutions, coord='phi', variable='temperature', gas=gas) + + plot_xt_diagram(times, solutions, variable='pressure', gas=gas) + plot_xt_diagram(times, solutions, coord='phi', variable='pressure', gas=gas) + + plot_xt_diagram(times, solutions, variable='velocity', gas=gas) + plot_xt_diagram(times, solutions, coord='phi', variable='velocity', gas=gas) + + # Plot final state + plot_final_state(times, solutions) + + # Optionally create animation frames + # create_animation_frames(times, solutions, output_dir='frames', variable='density') \ No newline at end of file From 80667e567d9610510dbc160de6883f670a64a617 Mon Sep 17 00:00:00 2001 From: Nadexterbrown Date: Thu, 28 Aug 2025 12:46:42 -0400 Subject: [PATCH 5/6] Added chemical mechanism directory and files --- .../Li-Dryer-H2-mechanism.yaml | 365 ++++++++++++++++++ 1 file changed, 365 insertions(+) create mode 100644 chemical_mechanisms_files/Li-Dryer-H2-mechanism.yaml diff --git a/chemical_mechanisms_files/Li-Dryer-H2-mechanism.yaml b/chemical_mechanisms_files/Li-Dryer-H2-mechanism.yaml new file mode 100644 index 0000000..a519927 --- /dev/null +++ b/chemical_mechanisms_files/Li-Dryer-H2-mechanism.yaml @@ -0,0 +1,365 @@ +description: |- + <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + H2/O2 oxidation reaction mechanism -- + (c) Li, Zhao, Kazakov, and Dryer, Princeton University, 2003. + + !!!!!!!!!!!!!!! IMPORTANT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + HOW TO USE THIS MECHANISM: + + Due to + (1) limitations of CHEMKIN-II format (specifically, an inability to implement + temperature-dependent collision efficiencies in falloff reactions) + and + (2) lack of fundamental understanding of the mixing rules for the falloff + reactions with the bath gases that have different broadening factors, + + the present implementation represents a compromise (approximate) formulation. + + As a consequence, PRIOR TO ITS USE IN THE CALCULATIONS, THIS FILE HAS TO BE + MODIFIED. DEPENDING ON WHAT BATH GAS (DILUTANT) IS MOST ABUNDANT IN YOUR SYSTEM + (THE PRESENT CHOICES ARE N2, AR, OR HE), YOU SHOULD UNCOMMENT THE CORRESPONDING + BLOCK FOR THE REACTION H+O2(+M)=HO2(+M), AND COMMENT THE BLOCK FOR OTHER DILUTANT(S). + AS GIVEN, THE MAIN DILUTANT IS SET TO BE N2. + + + HOW TO REFERENCE THIS MECHANISM: + + Li, J., Zhao, Z., Kazakov, A., and Dryer, F.L. "An Updated Comprehensive Kinetic Model + of Hydrogen Combustion", Int. J. Chem. Kinet. 2004 (in press). + + + HOW TO CONTACT THE AUTHORS: + + Prof. Frederick L. Dryer + D-329-D Engineering Quadrangle + Mechanical and Aerospace Engineering + Princeton University + Princeton, NJ 08544-5263 + Phone: 609-258-5206 + Lab: 609-258-0316 + FAX: 609-258-1939 + Email: fldryer@Princeton.EDU + + ********************************************************************************************** + Development notes: + + The following H2/O2 mechanism is based on Mueller et al's (Int.J.Chem.Kinet.1999,31:113) + Changes: + + 1.update the standard heat of formation of OH at 0K to 8.85kcal/mol (Ruscic et al, + J. Phys. Chem. A, 2002, 106:2727) + + 2.update the rate constant of H+O2=O+OH as proposed by Hessler (J. Phys. Chem. A, 1998, + 102:4517) + + 3.update the low-pressure-limit rate constant of H+O2(+M)=HO2(+M) with bath gases: H2, + O2, N2, AR, HE, H2O as proposed by Michael et al (J. Phys. Chem. A, 2002,106:5297). + The third-body efficiency of H2, O2, and H2O are taken as the average value over + the temperature range of 300-3000K. + The Fc in Troe's form with N2 and AR/HE as bath gas are different, so the fall-off + kinetics is expressed in two sets, for N2 and AR/HE, respectively. + + 4.for all other recombination reactions, assume the third-body efficiency of HE is + the same as AR. + + 5.modify the A factor of the rate constant of H+OH+M=H2O+M to 3.8E+22. + + END OF NOTES + ********************************************************************************************** + +generator: ck2yaml +input-files: [mechanism.inp, therm.dat, tran.dat] +cantera-version: 2.6.0 +date: Wed, 11 May 2022 17:40:07 -0700 + +units: {length: cm, time: s, quantity: mol, activation-energy: cal/mol} + +phases: +- name: gas + thermo: ideal-gas + elements: [H, O, N] + species: [H2, O2, H2O, H, O, OH, HO2, H2O2, N2] + kinetics: gas + transport: mixture-averaged + state: {T: 300.0, P: 1 atm} + +species: +- name: H2 + composition: {H: 2} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [3.29812431, 8.24944174e-04, -8.14301529e-07, -9.47543433e-11, 4.13487224e-13, + -1012.52087, -3.29409409] + - [2.99142337, 7.00064411e-04, -5.63382869e-08, -9.23157818e-12, 1.58275179e-15, + -835.033997, -1.35511017] + note: '121286' + transport: + model: gas + geometry: linear + well-depth: 38.0 + diameter: 2.92 + polarizability: 0.79 + rotational-relaxation: 280.0 +- name: O2 + composition: {O: 2} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [3.2129364, 1.12748635e-03, -5.75615047e-07, 1.31387723e-09, -8.76855392e-13, + -1005.24902, 6.03473759] + - [3.69757819, 6.13519689e-04, -1.25884199e-07, 1.77528148e-11, -1.13643531e-15, + -1233.93018, 3.18916559] + note: '121386' + transport: + model: gas + geometry: linear + well-depth: 107.4 + diameter: 3.458 + polarizability: 1.6 + rotational-relaxation: 3.8 +- name: H2O + composition: {H: 2, O: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [3.38684249, 3.47498246e-03, -6.35469633e-06, 6.96858127e-09, -2.50658847e-12, + -3.02081133e+04, 2.59023285] + - [2.67214561, 3.05629289e-03, -8.73026011e-07, 1.20099639e-10, -6.39161787e-15, + -2.9899209e+04, 6.86281681] + note: '20387' + transport: + model: gas + geometry: nonlinear + well-depth: 572.4 + diameter: 2.605 + dipole: 1.844 + rotational-relaxation: 4.0 +- name: H + composition: {H: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [2.5, 0.0, 0.0, 0.0, 0.0, 2.5471627e+04, -0.460117608] + - [2.5, 0.0, 0.0, 0.0, 0.0, 2.5471627e+04, -0.460117638] + note: '120186' + transport: + model: gas + geometry: atom + well-depth: 145.0 + diameter: 2.05 +- name: O + composition: {O: 1} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [2.94642878, -1.63816649e-03, 2.4210317e-06, -1.60284319e-09, 3.89069636e-13, + 2.91476445e+04, 2.96399498] + - [2.54205966, -2.75506191e-05, -3.10280335e-09, 4.55106742e-12, -4.3680515e-16, + 2.92308027e+04, 4.92030811] + note: '120186' + transport: + model: gas + geometry: atom + well-depth: 80.0 + diameter: 2.75 +- name: OH + composition: {O: 1, H: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 6000.0] + data: + - [4.12530561, -3.22544939e-03, 6.52764691e-06, -5.79853643e-09, 2.06237379e-12, + 3346.30913, -0.69043296] + - [2.86472886, 1.05650448e-03, -2.59082758e-07, 3.05218674e-11, -1.33195876e-15, + 3683.62875, 5.70164073] + note: S 9/01 + transport: + model: gas + geometry: linear + well-depth: 80.0 + diameter: 2.75 +- name: HO2 + composition: {H: 1, O: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [4.30179801, -4.74912051e-03, 2.11582891e-05, -2.42763894e-08, 9.29225124e-12, + 294.80804, 3.71666245] + - [4.0172109, 2.23982013e-03, -6.3365815e-07, 1.1424637e-10, -1.07908535e-14, + 111.856713, 3.78510215] + note: L 5/89 + transport: + model: gas + geometry: nonlinear + well-depth: 107.4 + diameter: 3.458 + rotational-relaxation: 1.0 +- name: H2O2 + composition: {H: 2, O: 2} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [3.38875365, 6.56922581e-03, -1.48501258e-07, -4.62580552e-09, 2.47151475e-12, + -1.76631465e+04, 6.7853632] + - [4.57316685, 4.33613639e-03, -1.47468882e-06, 2.34890357e-10, -1.43165356e-14, + -1.80069609e+04, 0.501136959] + note: '120186' + transport: + model: gas + geometry: nonlinear + well-depth: 107.4 + diameter: 3.458 + rotational-relaxation: 3.8 +- name: N2 + composition: {N: 2} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [3.298677, 1.40824e-03, -3.963222e-06, 5.641515e-09, -2.444855e-12, + -1020.9, 3.950372] + - [2.92664, 1.487977e-03, -5.684761e-07, 1.009704e-10, -6.753351e-15, + -922.7977, 5.980528] + note: '121286' + transport: + model: gas + geometry: linear + well-depth: 97.53 + diameter: 3.621 + polarizability: 1.76 + rotational-relaxation: 4.0 + +reactions: +- equation: H + O2 <=> O + OH # Reaction 1 + rate-constant: {A: 3.547e+15, b: -0.406, Ea: 1.6599e+04} + note: |- + H2-O2 Chain Reactions + Hessler, J. Phys. Chem. A, 102:4517 (1998) +- equation: O + H2 <=> H + OH # Reaction 2 + rate-constant: {A: 5.08e+04, b: 2.67, Ea: 6290.0} + note: Sutherland et al., 21st Symposium, p. 929 (1986) +- equation: H2 + OH <=> H2O + H # Reaction 3 + rate-constant: {A: 2.16e+08, b: 1.51, Ea: 3430.0} + note: Michael and Sutherland, J. Phys. Chem. 92:3853 (1988) +- equation: O + H2O <=> OH + OH # Reaction 4 + rate-constant: {A: 2.97e+06, b: 2.02, Ea: 1.34e+04} + note: Sutherland et al., 23rd Symposium, p. 51 (1990) +- equation: H2 + M <=> H + H + M # Reaction 5 + type: three-body + rate-constant: {A: 4.577e+19, b: -1.4, Ea: 1.0438e+05} + efficiencies: {H2: 2.5, H2O: 12.0} + note: |- + H2-O2 Dissociation Reactions + Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) +- equation: O + O + M <=> O2 + M # Reaction 6 + type: three-body + rate-constant: {A: 6.165e+15, b: -0.5, Ea: 0.0} + efficiencies: {H2: 2.5, H2O: 12.0} + note: |2- + CO/1.9/ CO2/3.8/ + AR/0.0/ HE/0.0/ + Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) + H2+AR=H+H+AR 5.84e18 -1.1 1.0438E+05 + H2+HE=H+H+HE 5.84e18 -1.1 1.0438E+05 + Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) +- equation: O + H + M <=> OH + M # Reaction 7 + type: three-body + rate-constant: {A: 4.714e+18, b: -1.0, Ea: 0.0} + efficiencies: {H2: 2.5, H2O: 12.0} + note: |2- + AR/0.0/ HE/0.0/ + CO/1.9/ CO2/3.8/ + Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) + O+O+AR=O2+AR 1.886E+13 0.00 -1.788E+03 + O+O+HE=O2+HE 1.886E+13 0.00 -1.788E+03 + Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) +- equation: H + OH + M <=> H2O + M # Reaction 8 + type: three-body + rate-constant: {A: 3.8e+22, b: -2.0, Ea: 0.0} + efficiencies: {H2: 2.5, H2O: 12.0} + note: |2- + AR/0.75/ HE/0.75/ + CO/1.9/ CO2/3.8/ + Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) + H+OH+M=H2O+M 2.212E+22 -2.00 0.000E+00 +- equation: H + O2 (+M) <=> HO2 (+M) # Reaction 9 + type: falloff + low-P-rate-constant: {A: 6.366e+20, b: -1.72, Ea: 524.8} + high-P-rate-constant: {A: 1.475e+12, b: 0.6, Ea: 0.0} + Troe: {A: 0.8, T3: 1.0e-30, T1: 1.0e+30} + efficiencies: {H2: 2.0, H2O: 11.0, O2: 0.78} + note: |2- + AR/0.38/ HE/0.38/ + CO/1.9/ CO2/3.8/ + Formation and Consumption of HO2 + Cobos et al., J. Phys. Chem. 89:342 (1985) for kinf + Michael, et al., J. Phys. Chem. A, 106:5297 (2002) for k0 + ****************************************************************************** + MAIN BATH GAS IS N2 (comment this reaction otherwise) +- equation: HO2 + H <=> H2 + O2 # Reaction 10 + rate-constant: {A: 1.66e+13, b: 0.0, Ea: 823.0} + note: |- + CO/1.9/ CO2/3.8/ + ****************************************************************************** + MAIN BATH GAS IS AR OR HE (comment this reaction otherwise) + H+O2(+M)=HO2(+M) 1.475E+12 0.60 0.00E+00 + LOW/9.042E+19 -1.50 4.922E+02/ + TROE/0.5 1E-30 1E+30/ + H2/3.0/ H2O/16/ O2/1.1/ CO/2.7/ CO2/5.4/ HE/1.2/ + Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) [modified] +- equation: HO2 + H <=> OH + OH # Reaction 11 + rate-constant: {A: 7.079e+13, b: 0.0, Ea: 295.0} + note: Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) [modified] +- equation: HO2 + O <=> O2 + OH # Reaction 12 + rate-constant: {A: 3.25e+13, b: 0.0, Ea: 0.0} + note: Baulch et al., J. Phys. Chem. Ref Data, 21:411 (1992) +- equation: HO2 + OH <=> H2O + O2 # Reaction 13 + rate-constant: {A: 2.89e+13, b: 0.0, Ea: -497.0} + note: Keyser, J. Phys. Chem. 92:1193 (1988) +- equation: HO2 + HO2 <=> H2O2 + O2 # Reaction 14 + duplicate: true + rate-constant: {A: 4.2e+14, b: 0.0, Ea: 1.1982e+04} + note: |- + Formation and Consumption of H2O2 + Hippler et al., J. Chem. Phys. 93:1755 (1990) +- equation: HO2 + HO2 <=> H2O2 + O2 # Reaction 15 + duplicate: true + rate-constant: {A: 1.3e+11, b: 0.0, Ea: -1629.3} +- equation: H2O2 (+M) <=> OH + OH (+M) # Reaction 16 + type: falloff + low-P-rate-constant: {A: 1.202e+17, b: 0.0, Ea: 4.55e+04} + high-P-rate-constant: {A: 2.951e+14, b: 0.0, Ea: 4.843e+04} + Troe: {A: 0.5, T3: 1.0e-30, T1: 1.0e+30} + efficiencies: {H2: 2.5, H2O: 12.0} + note: |- + Brouwer et al., J. Chem. Phys. 86:6171 (1987) for kinf + Warnatz, J. in Combustion chemistry (1984) for k0 +- equation: H2O2 + H <=> H2O + OH # Reaction 17 + rate-constant: {A: 2.41e+13, b: 0.0, Ea: 3970.0} + note: |2- + CO/1.9/ CO2/3.8/ + AR/0.64/ HE/0.64/ + Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) +- equation: H2O2 + H <=> HO2 + H2 # Reaction 18 + rate-constant: {A: 4.82e+13, b: 0.0, Ea: 7950.0} + note: Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) +- equation: H2O2 + O <=> OH + HO2 # Reaction 19 + rate-constant: {A: 9.55e+06, b: 2.0, Ea: 3970.0} + note: Tsang and Hampson, J. Phys. Chem. Ref. Data, 15:1087 (1986) +- equation: H2O2 + OH <=> HO2 + H2O # Reaction 20 + duplicate: true + rate-constant: {A: 1.0e+12, b: 0.0, Ea: 0.0} + note: Hippler and Troe, J. Chem. Phys. Lett. 192:333 (1992) +- equation: H2O2 + OH <=> HO2 + H2O # Reaction 21 + duplicate: true + rate-constant: {A: 5.8e+14, b: 0.0, Ea: 9557.0} From 0e970a38ceac0fe3ab2d961ccfa0e6be6401f202 Mon Sep 17 00:00:00 2001 From: Nadexterbrown Date: Thu, 28 Aug 2025 12:47:19 -0400 Subject: [PATCH 6/6] Data processing from external files using piston boundary conditions. --- piston_data_processing.py | 615 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 615 insertions(+) create mode 100644 piston_data_processing.py diff --git a/piston_data_processing.py b/piston_data_processing.py new file mode 100644 index 0000000..499be8a --- /dev/null +++ b/piston_data_processing.py @@ -0,0 +1,615 @@ +import numpy as np +import cantera as ct +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm +from typing import Tuple, List, Optional, Union, Dict + +# Updated import to use the new piston functions file +from piston_setup import simple_piston, datafile_piston, initialize_piston_mode, piston_velocity +from problem_setup import initialize_parameters, input_params # your dataclass + +from lagrangian import ( + PlanarGeometry, + LagrangianMesh, + LagrangianSolver, + SimulationParameters +) + + +######################################################################################################################## +# Piston Problem Simulation +######################################################################################################################## + +def run_piston_problem( + n_cells: int = 100, + t_end: float = 0.5, + cfl: float = 0.3, + output_frequency: int = 10, + *, + # ---- pass ONE of these sources of thermodynamic state ---- + gas: "ct.Solution | None" = None, + params: "MyClass | None" = None, + uniform_state: "tuple[float, float] | None" = None, # (rho0 [kg/m^3], p0 [Pa]) + override_gamma: Optional[float] = None, + initial: str = "uniform", # or "riemann_cltorc" + piston_enabled: bool = True # New parameter to enable/disable piston +) -> Tuple[np.ndarray, List[np.ndarray]]: + """ + Run a piston-in-gas simulation with thermodynamics supplied by caller. + No internal initialization of input_params happens here. + + Provide exactly one of: + - gas: a Cantera Solution with TPX set + - params: your MyClass object; we call params.build_gas() + - uniform_state: (rho0, p0) in SI + + If override_gamma is None, gamma = cp/cv from the gas (when available). + """ + + # --------------------------- + # Resolve thermodynamic state + # --------------------------- + provided = sum(x is not None for x in (gas, params, uniform_state)) + if provided != 1: + raise ValueError("Provide exactly ONE of gas=, params=, or uniform_state=.") + + if params is not None: + if not hasattr(params, "build_gas"): + raise TypeError("params must be a MyClass with .build_gas()") + if ct is None: + raise RuntimeError("Cantera required when using params= or gas=.") + gas = params.build_gas() + + if gas is not None: + p0 = float(gas.P) # Pa + rho0 = float(gas.density) # kg/m^3 + gamma = float(override_gamma) if override_gamma is not None else float(gas.cp_mass / gas.cv_mass) + else: + # uniform_state provided + rho0, p0 = map(float, uniform_state) + if override_gamma is None: + raise ValueError("uniform_state provided but override_gamma is None. Specify gamma.") + gamma = float(override_gamma) + + if not (gamma > 1.0): + raise ValueError(f"gamma must be > 1. Got {gamma}.") + + # specific internal energy for a gamma-law closure + e0 = p0 / ((gamma - 1.0) * rho0) + + # --------------------------- + # Mesh / geometry + # --------------------------- + geometry = PlanarGeometry() + x_max = 2.5 # domain length + r = np.linspace(0.0, x_max, n_cells + 1) # interface positions + + # --------------------------- + # Initial conditions + # --------------------------- + rho = np.empty(n_cells) + u = np.empty(n_cells + 1) + e = np.empty(n_cells) + + rho[:] = rho0 + u[:] = 0.0 + e[:] = e0 + + # Create mesh and solver + mesh = LagrangianMesh(r, rho, u, e, geometry, gamma) + + sim_params = SimulationParameters( + gamma=gamma, + cfl=cfl, + t_end=t_end, + output_frequency=output_frequency, + use_artificial_viscosity=True + ) + + # Create a custom solver that overrides the boundary conditions + class PistonSolver(LagrangianSolver): + def __init__(self, mesh, params, piston_enabled=True): + super().__init__(mesh, params) + self.piston_enabled = piston_enabled + + def apply_boundary_conditions(self) -> None: + """Apply piston boundary condition on the left and fixed wall on the right.""" + # Left boundary (piston) + if not self.piston_enabled: + self.mesh.u[0] = 0.0 # Fixed wall if piston disabled + else: + # Use the global piston_velocity function + # No need to pass piston_params since it's handled globally + self.mesh.u[0] = piston_velocity(self.time) + + # Right boundary (fixed wall) + self.mesh.u[-1] = 0.0 + + solver = PistonSolver(mesh, sim_params, piston_enabled=piston_enabled) + + # Solve + return solver.solve() + + +######################################################################################################################## +# Visualization Functions +######################################################################################################################## + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm +from typing import List, Optional, Union, Dict +import cantera as ct + + +def plot_xt_diagram( + times: np.ndarray, + solutions: List[np.ndarray], # (r, u, rho, e, p) + variable: str = 'density', + output_file: Optional[str] = None, + *, + gas: Optional[ct.Solution] = None, + mech: Optional[str] = None, + X: Optional[Union[str, Dict[str, float]]] = None, + log_for: Optional[set] = None, + var_limit: Optional[Tuple[float, float]] = None, + coord: str = 'x', # 'x' or 'phi' + return_phi: bool = False +) -> Optional[List[tuple]]: + """ + Plot an x–t diagram. If coord='phi', use the mass-weighted Lagrangian + coordinate phi(x,t)=∫ rho dx on the horizontal axis. + + solutions[i] is (r, u, rho, e, p) with shapes: + r: (Ni+1,) interfaces + u: (Ni+1,) interface velocity + rho: (Ni,) cell-centered density + e: (Ni,) cell-centered internal energy + p: (Ni,) cell-centered pressure + """ + + # ----------------------- + # Setup gas / composition (only if needed) + # ----------------------- + native = {'density': 2, 'velocity': 1, 'pressure': 4, 'energy': 3} + derived_requested = variable not in native + + if gas is None and derived_requested: + if mech is None: + try: + mech = input_params.mech + except Exception: + raise ValueError("Need `mech` (or a ready `gas`) to reconstruct state with Cantera.") + if X is None: + try: + X = input_params.X + except Exception: + raise ValueError("Need mixture composition `X` (or set input_params.X).") + gas = ct.Solution(mech) + X_str = ",".join([f"{k}:{v}" for k, v in X.items()]) if isinstance(X, dict) else X + else: + X_str = gas.X.copy() + + # Pretty names + pretty = { + 'density': 'Density', 'velocity': 'Velocity', 'pressure': 'Pressure', 'energy': 'Internal Energy', + 'temperature': 'Temperature', 'cp': 'cp', 'cv': 'cv', 'gamma': 'Gamma', 'R': 'Gas constant (mix)', + 'sound_speed': 'Sound speed', 'mu': 'Viscosity', 'lambda': 'Thermal conductivity', + 'enthalpy': 'Specific enthalpy', 'entropy': 'Specific entropy', 'Mach': 'Mach number' + } + var_name = pretty.get(variable, variable) + + # Defaults for log scaling + default_log_set = {'density', 'pressure', 'cp', 'cv', 'mu', 'lambda'} + if log_for is None: + log_for = default_log_set + + # Species request parsing + def _parse_species_token(vname: str): + if vname.startswith('Y(') and vname.endswith(')'): return ('Y', vname[2:-1]) + if vname.startswith('X(') and vname.endswith(')'): return ('X', vname[2:-1]) + if vname.startswith('Y_'): return ('Y', vname[2:]) + if vname.startswith('X_'): return ('X', vname[2:]) + return (None, None) + + req_kind, req_species = _parse_species_token(variable) + + # ----------------------- + # φ helper + # ----------------------- + def _compute_phi( + r_iface: np.ndarray, + rho_cell: np.ndarray, + x_piston: float | None = None, + *, x_units: str = "m", rho_units: str = "kg/m^3" + ): + xu = {"m": 1.0, "cm": 1e-2, "mm": 1e-3}[x_units] + ru = {"kg/m^3": 1.0, "g/cm^3": 1000.0}[rho_units] + + r = r_iface.astype(float) * xu + rho = np.maximum(rho_cell, 0.0).astype(float) * ru + + dx = r[1:] - r[:-1] # (Ni,) + m = rho * dx # mass/area per cell (kg/m^2) + + if x_piston is None: + x_piston = r[0] + else: + x_piston = float(x_piston) * xu + + idx = np.searchsorted(r, x_piston, side="right") - 1 + idx = max(0, min(idx, len(dx) - 1)) + + m0 = rho[idx] * max(r[idx + 1] - x_piston, 0.0) + + phi_iface = np.zeros_like(r) + if idx + 1 < len(r): + cs = np.cumsum(m[idx + 1:]) # length Ni - idx - 1 + phi_iface[idx + 1:] = m0 + np.concatenate(([0.0], cs)) # <-- fixed + + phi_cell = 0.5 * (phi_iface[:-1] + phi_iface[1:]) + return phi_iface, phi_cell + + # ----------------------- + # Accumulators + # ----------------------- + x_values = [] + t_values = [] + var_values = [] + phi_list = [] # Will store (phi_iface, phi_cell) per time step if requested + + # ----------------------- + # Main loop + # ----------------------- + for i, t in enumerate(times): + r, ufaces, rho, e, p = solutions[i] + Ni = rho.shape[0] + + # Compute φ if needed (or if caller wants it returned) + if coord == 'phi' or return_phi: + phi_iface, phi_cell = _compute_phi(r, rho) + if return_phi: + phi_list.append((phi_iface.copy(), phi_cell.copy())) + + # Determine alignment and base x-grid + if variable in native: + if variable == 'velocity': + # interface quantity + x_plot = phi_iface if coord == 'phi' else r + vals = ufaces + for j in range(len(x_plot)): + x_values.append(x_plot[j]) + t_values.append(t) + var_values.append(vals[j]) + else: + # cell-centered + x_centers = phi_cell if coord == 'phi' else 0.5 * (r[:-1] + r[1:]) + if variable == 'density': + vals = rho + elif variable == 'pressure': + vals = p + elif variable == 'energy': + vals = e + for j in range(Ni): + x_values.append(x_centers[j]) + t_values.append(t) + var_values.append(vals[j]) + + else: + # Derived thermodynamic variable → cell-centered + x_centers = phi_cell if coord == 'phi' else 0.5 * (r[:-1] + r[1:]) + # For Mach we need a cell-centered u + u_cell = 0.5 * (ufaces[:-1] + ufaces[1:]) if variable == 'Mach' else None + + vals = np.empty(Ni, dtype=float) + for j in range(Ni): + gas.DPX = rho[j], p[j], X_str + if req_kind == 'Y': + k = gas.species_index(req_species); + vals[j] = gas.Y[k] + elif req_kind == 'X': + k = gas.species_index(req_species); + vals[j] = gas.X[k] + elif variable == 'temperature': + vals[j] = gas.T + elif variable == 'cp': + vals[j] = gas.cp_mass + elif variable == 'cv': + vals[j] = gas.cv_mass + elif variable == 'gamma': + vals[j] = gas.cp_mass / gas.cv_mass + elif variable == 'R': + vals[j] = gas.gas_constant / gas.mean_molecular_weight + elif variable == 'sound_speed': + vals[j] = gas.sound_speed + elif variable == 'mu': + vals[j] = gas.viscosity + elif variable == 'lambda': + vals[j] = gas.thermal_conductivity + elif variable == 'enthalpy': + vals[j] = gas.h_mass + elif variable == 'entropy': + vals[j] = gas.s_mass + elif variable == 'Mach': + a = gas.sound_speed + vals[j] = abs(u_cell[j]) / max(a, 1e-16) + else: + raise ValueError( + f"Unknown derived variable '{variable}'. " + f"Supported: temperature, cp, cv, gamma, R, sound_speed, mu, lambda, " + f"enthalpy, entropy, Mach, and species via Y(NAME)/X(NAME)/Y_NAME/X_NAME." + ) + + for j in range(Ni): + x_values.append(x_centers[j]) + t_values.append(t) + var_values.append(vals[j]) + + # ----------------------- + # Plot + # ----------------------- + x_values = np.asarray(x_values) + t_values = np.asarray(t_values) + var_values = np.asarray(var_values) + + plt.figure(figsize=(10, 8)) + use_log = variable in log_for and np.all(var_values > 0.0) + if use_log: + sc = plt.scatter(x_values, t_values, c=var_values, cmap='rainbow', + s=10, alpha=0.8, norm=LogNorm()) + else: + if var_limit is None: + var_limit = (np.min(var_values), np.max(var_values)) + sc = plt.scatter(x_values, t_values, c=var_values, cmap='rainbow', + s=10, alpha=0.8, vmin=var_limit[0], vmax=var_limit[1]) + + cbar = plt.colorbar(sc) + cbar.set_label(var_name) + + if coord == 'phi': + plt.xlabel(r'$\phi$ (mass per unit area)') + title_axis = 'φ' + else: + plt.xlabel('Position (x)') + title_axis = 'x' + + plt.ylabel('Time (t)') + plt.title(f'{title_axis}–T Diagram: {var_name}') + + if output_file: + plt.savefig(output_file, dpi=300, bbox_inches='tight') + print(f"Plot saved to {output_file}") + + plt.tight_layout() + plt.show() + + if return_phi: + return phi_list + return None + + +def plot_final_state( + times: np.ndarray, + solutions: List[np.ndarray], + output_file: Optional[str] = None +) -> None: + """ + Plot the final state of the simulation. + + Args: + times: Array of simulation times + solutions: List of solution snapshots (r, u, rho, e, p) + output_file: Optional file to save plot to + """ + # Extract final solution + r_final, u_final, rho_final, e_final, p_final = solutions[-1] + + # Create figure with subplots + fig, axs = plt.subplots(2, 2, figsize=(12, 10)) + + # Plot density + axs[0, 0].plot(0.5 * (r_final[:-1] + r_final[1:]), rho_final, 'b-', linewidth=2) + axs[0, 0].set_title('Density') + axs[0, 0].set_ylabel('Density') + axs[0, 0].grid(True) + + # Plot velocity + axs[0, 1].plot(r_final, u_final, 'r-', linewidth=2) + axs[0, 1].set_title('Velocity') + axs[0, 1].set_ylabel('Velocity') + axs[0, 1].grid(True) + + # Plot pressure + axs[1, 0].plot(0.5 * (r_final[:-1] + r_final[1:]), p_final, 'g-', linewidth=2) + axs[1, 0].set_title('Pressure') + axs[1, 0].set_xlabel('Position') + axs[1, 0].set_ylabel('Pressure') + axs[1, 0].grid(True) + + # Plot internal energy + axs[1, 1].plot(0.5 * (r_final[:-1] + r_final[1:]), e_final, 'm-', linewidth=2) + axs[1, 1].set_title('Specific Internal Energy') + axs[1, 1].set_xlabel('Position') + axs[1, 1].set_ylabel('Energy') + axs[1, 1].grid(True) + + plt.tight_layout() + + # Add main title + plt.suptitle(f'Piston Problem - Final State at t = {times[-1]:.3f}', fontsize=16) + plt.subplots_adjust(top=0.93) + + # Save or show + if output_file: + plt.savefig(output_file, dpi=300, bbox_inches='tight') + print(f"Plot saved to {output_file}") + + plt.show() + + +def create_animation_frames( + times: np.ndarray, + solutions: List[np.ndarray], + output_dir: str = 'frames', + variable: str = 'density' +) -> None: + """ + Create animation frames of the solution evolution. + + Args: + times: Array of simulation times + solutions: List of solution snapshots (r, u, rho, e, p) + output_dir: Directory to save frames + variable: Which variable to plot ('density', 'velocity', 'pressure', 'energy') + """ + import os + + # Create output directory if it doesn't exist + os.makedirs(output_dir, exist_ok=True) + + # Get variable index + var_idx = {'density': 2, 'velocity': 1, 'pressure': 4, 'energy': 3}[variable] + var_name = {'density': 'Density', 'velocity': 'Velocity', + 'pressure': 'Pressure', 'energy': 'Internal Energy'}[variable] + + # Initialize y-axis limits + y_min = float('inf') + y_max = float('-inf') + + # Create frames + for i, time in enumerate(times): + r, u, rho, e, p = solutions[i] + + # For cell-centered quantities (density, pressure, energy) + if var_idx in [2, 3, 4]: + x = 0.5 * (r[:-1] + r[1:]) # Cell centers + var = solutions[i][var_idx] + # For interface quantities (velocity) + else: + x = r # Interface positions + var = solutions[i][var_idx] + + # Create figure + plt.figure(figsize=(10, 6)) + plt.plot(x, var, 'b-', linewidth=2) + plt.title(f'{var_name} at t = {time:.3f}') + plt.xlabel('Position') + plt.ylabel(var_name) + plt.grid(True) + + # Set consistent y-axis limits + if i == 0: + y_min = min(var) + y_max = max(var) + else: + y_min = min(y_min, min(var)) + y_max = max(y_max, max(var)) + + plt.ylim(y_min, y_max * 1.1) + + # Save frame + plt.savefig(f'{output_dir}/frame_{i:04d}.png', dpi=200, bbox_inches='tight') + plt.close() + + print(f"Created {len(times)} animation frames in {output_dir}") + + +def plot_piston_velocity_vs_time(times: np.ndarray, output_file: Optional[str] = None) -> None: + """ + Plot the piston velocity as a function of time to visualize the boundary condition. + + Args: + times: Array of simulation times + output_file: Optional file to save plot to + """ + # Calculate piston velocity at each time step + piston_velocities = [piston_velocity(t) for t in times] + + plt.figure(figsize=(10, 6)) + plt.plot(times, piston_velocities, 'b-', linewidth=2) + plt.title('Piston Velocity vs Time') + plt.xlabel('Time [s]') + plt.ylabel('Piston Velocity [m/s]') + plt.grid(True) + + if output_file: + plt.savefig(output_file, dpi=300, bbox_inches='tight') + print(f"Piston velocity plot saved to {output_file}") + + plt.show() + + +######################################################################################################################## +# Main Execution +######################################################################################################################## + +if __name__ == "__main__": + # Gas Parameters + initialize_parameters( + T=503.0, + P=10.0e5, + Phi=1.0, + Fuel="H2", + # nitrogenAmount=0.5 * 3.76, + mech="./chemical_mechanism_files/Li-Dryer-H2-mechanism.yaml", + temperature_unit="K", + pressure_unit="Pa", # convenience; converted to Pa internally + phi_convention="scale_fuel", # or "scale_oxidizer" + use_dry_air=False # Use dry air for nitrogen/oxygen ratio + ) + print(input_params.summary()) + gas = input_params.build_gas() + print(f"Built gas with {gas.T:.1f} K, {gas.P / ct.one_atm:.3f} atm, X={gas.X}") + + # Simulation parameters + n_cells = 1000 + gamma = 1.4 + t_end = 0.0025537811000000002 + cfl = 0.3 + output_frequency = 5 # Save more frequently for better visualization + + # Datafile piston (uncomment to use instead of simple mode) + print("Initializing piston in DATAFILE mode...") + initialize_piston_mode( + mode='datafile', + datafile_path='Pele-Data/H2-O2-Phi-1.0-P-10-bar-T-503-K/Level-6/Wave-Tracking-Concatenated-Results-V34.txt', + #datafile_path='' + # Optional: specify column names (auto-detected if not provided) + time_column='Time [s]', # will auto-find if None + velocity_column='Flame Velocity [m/s]', # will auto-find if None + interpolation_kind='linear' # 'cubic', 'quadratic', etc. + ) + + # Run simulation + print("Starting simulation...") + times, solutions = run_piston_problem( + n_cells=n_cells, + t_end=t_end, + cfl=cfl, + output_frequency=output_frequency, + gas=gas, + piston_enabled=True # Set to False to disable piston (fixed wall instead) + ) + + print(f"Simulation completed.") + print(f"Number of time steps: {len(times)}") + print(f"Final time: {times[-1]:.6f} s") + + # Plot piston velocity vs time to verify boundary condition + plot_piston_velocity_vs_time(times) + + # Plot x-t diagram for temperature + plot_xt_diagram(times, solutions, variable='temperature', gas=gas) + plot_xt_diagram(times, solutions, coord='phi', variable='temperature', gas=gas) + + plot_xt_diagram(times, solutions, variable='pressure', gas=gas) + plot_xt_diagram(times, solutions, coord='phi', variable='pressure', gas=gas) + + plot_xt_diagram(times, solutions, variable='velocity', gas=gas) + plot_xt_diagram(times, solutions, coord='phi', variable='velocity', gas=gas) + + # Plot final state + plot_final_state(times, solutions) + + # Optionally create animation frames + # create_animation_frames(times, solutions, output_dir='frames', variable='density') \ No newline at end of file