ignition_delay.py · Scientific Agent Skills

1#!/usr/bin/env python3

2"""Run closed adiabatic ideal-gas ignition with numerical and conservation checks."""

3from __future__ import annotations

4

5import argparse

6import csv

7import hashlib

8import json

9from pathlib import Path

10

11import cantera as ct

12import numpy as np

13

14

15def number(value, name):

16 if isinstance(value, bool) or not isinstance(value, (int, float)) or not np.isfinite(value) or value <= 0:

17 raise ValueError(f"{name} must be a finite positive number")

18 return float(value)

19

20

21def validate(config):

22 keys = {"mechanism", "phase", "reactor", "temperature_k", "pressure_pa", "mole_amounts",

23 "end_time_s", "samples", "rtol", "atol", "max_time_step_s", "minimum_temperature_rise_k",

24 "tracked_species", "delay_relative_tolerance"}

25 if set(config) != keys:

26 raise ValueError(f"Missing settings: {sorted(keys-set(config))}; unknown: {sorted(set(config)-keys)}")

27 for name in keys - {"mechanism", "phase", "reactor", "mole_amounts", "samples", "tracked_species"}:

28 number(config[name], name)

29 n = config["samples"]

30 if isinstance(n, bool) or not isinstance(n, int) or not 11 <= n <= 100000:

31 raise ValueError("samples must be an integer from 11 through 100000")

32 if config["reactor"] not in ("constant-volume", "constant-pressure"):

33 raise ValueError("reactor must be constant-volume or constant-pressure")

34 if not isinstance(config["mole_amounts"], dict) or not config["mole_amounts"]:

35 raise ValueError("mole_amounts must be a nonempty species-to-amount mapping")

36 amounts = config["mole_amounts"]

37 if any(isinstance(v, bool) or not isinstance(v, (int, float)) or not np.isfinite(v) or v < 0 for v in amounts.values()) or sum(amounts.values()) <= 0:

38 raise ValueError("mole_amounts must be finite, nonnegative, and have a positive sum")

39 if not isinstance(config["tracked_species"], list) or len(set(config["tracked_species"])) != len(config["tracked_species"]):

40 raise ValueError("tracked_species must be a list of unique names")

41

42

43def load_gas(config):

44 gas = ct.Solution(config["mechanism"], config["phase"])

45 if gas.thermo_model != "ideal-gas":

46 raise ValueError("This helper supports only ideal-gas phases")

47 unknown = (set(config["mole_amounts"]) | set(config["tracked_species"])) - set(gas.species_names)

48 if unknown:

49 raise ValueError(f"Unknown species: {sorted(unknown)}")

50 gas.TPX = config["temperature_k"], config["pressure_pa"], config["mole_amounts"]

51 return gas

52

53

54def delay_from_history(time, temperature, minimum_rise):

55 rise = float(np.max(temperature) - temperature[0])

56 rate = np.gradient(temperature, time, edge_order=2)

57 index = int(np.argmax(rate))

58 status = "resolved"

59 if rise < minimum_rise:

60 status = "temperature-rise-not-reached"

61 elif index < 2 or index >= len(time) - 2:

62 status = "derivative-maximum-at-boundary"

63 return {"status": status, "delay_s": float(time[index]) if status == "resolved" else None,

64 "maximum_temperature_rise_k": rise, "peak_heating_rate_k_per_s": float(rate[index]),

65 "output_spacing_s": float(time[1] - time[0])}

66

67

68def simulate(config):

69 validate(config)

70 gas = load_gas(config)

71 reactor_type = ct.IdealGasReactor if config["reactor"] == "constant-volume" else ct.IdealGasConstPressureReactor

72 reactor = reactor_type(gas, energy="on", volume=1.0, clone=True)

73 network = ct.ReactorNet([reactor])

74 network.rtol, network.atol = config["rtol"], config["atol"]

75 network.max_time_step = config["max_time_step_s"]

76 network.max_steps = 100000

77 times = np.linspace(0, config["end_time_s"], config["samples"])

78 rows, elements = [], []

79 species_indices = [gas.species_index(name) for name in config["tracked_species"]]

80 minimum_y = 0.0

81 normalization_error = 0.0

82 for time in times:

83 if time:

84 network.advance(float(time))

85 phase = reactor.phase

86 rows.append([time, reactor.T, phase.P, reactor.volume, reactor.mass,

87 reactor.mass * phase.int_energy_mass, reactor.mass * phase.enthalpy_mass,

88 *phase.X[species_indices]])

89 elements.append([phase.elemental_mass_fraction(name) for name in phase.element_names])

90 minimum_y = min(minimum_y, float(phase.Y.min()))

91 normalization_error = max(normalization_error, abs(float(phase.Y.sum()) - 1))

92 history = np.asarray(rows)

93 if not np.isfinite(history).all():

94 raise ValueError("Integrator returned a nonfinite state")

95 energy_col = 5 if config["reactor"] == "constant-volume" else 6

96 conserved = history[:, energy_col]

97 conservation = {

98 "mass_relative_drift": float(np.max(np.abs(history[:, 4] / history[0, 4] - 1))),

99 "element_mass_fraction_absolute_drift": float(np.max(np.abs(np.asarray(elements) - elements[0]))),

100 "conserved_quantity": "total_internal_energy_j" if energy_col == 5 else "total_enthalpy_j",

101 "energy_scaled_drift": float(np.max(np.abs(conserved - conserved[0])) / max(abs(conserved[0]), 1.0)),

102 "energy_denominator_j": max(abs(float(conserved[0])), 1.0),

103 "minimum_species_mass_fraction": minimum_y,

104 "mass_fraction_sum_error": normalization_error,

105 }

106 conservation["passed"] = (conservation["mass_relative_drift"] < 1e-8 and

107 conservation["element_mass_fraction_absolute_drift"] < 1e-8 and

108 conservation["energy_scaled_drift"] < 1e-6 and

109 minimum_y > -1e-10 and normalization_error < 1e-8)

110 result = delay_from_history(history[:, 0], history[:, 1], config["minimum_temperature_rise_k"])

111 result.update(conservation=conservation, final_temperature_k=float(history[-1, 1]),

112 samples=config["samples"], rtol=config["rtol"], atol=config["atol"],

113 max_time_step_s=config["max_time_step_s"], end_time_s=config["end_time_s"],

114 thermo_temperature_range_k=[gas.min_temp, gas.max_temp],

115 within_thermo_temperature_range=bool(history[:, 1].min() >= gas.min_temp and history[:, 1].max() <= gas.max_temp))

116 return history, result

117

118

119def run(config_path, output_dir):

120 config_path, output_dir = Path(config_path), Path(output_dir)

121 if output_dir.exists():

122 raise ValueError("Output directory exists; choose a new directory")

123 config = json.loads(config_path.read_text(encoding="utf-8"))

124 validate(config)

125 if config["samples"] > 50000:

126 raise ValueError("CLI samples must be at most 50000 to permit refinement runs")

127 # Local relative mechanism paths resolve against the configuration, before built-in lookup.

128 local = config_path.parent / config["mechanism"]

129 effective = dict(config)

130 if local.is_file():

131 effective["mechanism"] = str(local.resolve())

132 gas = load_gas(effective)

133 # Hash the exact bytes written, including on platforms with CRLF text translation.

134 snapshot = gas.write_yaml(precision=17).encode("utf-8")

135 source = next((Path(d) / effective["mechanism"] for d in ct.get_data_directories()

136 if (Path(d) / effective["mechanism"]).is_file()), None)

137 variants = {

138 "baseline": effective,

139 "finer_output": {**effective, "samples": 2 * effective["samples"] - 1},

140 "tighter_solver": {**effective, "rtol": effective["rtol"] / 10,

141 "atol": effective["atol"] / 10,

142 "max_time_step_s": effective["max_time_step_s"] / 2},

143 "longer_horizon": {**effective, "end_time_s": 2 * effective["end_time_s"],

144 "samples": 2 * effective["samples"] - 1},

145 }

146 results, histories = {}, {}

147 for name, conditions in variants.items():

148 history, summary = simulate(conditions)

149 histories[name], results[name] = history, summary

150 baseline_delay = results["baseline"]["delay_s"]

151 sensitivity = {}

152 for name, result in results.items():

153 if name == "baseline":

154 continue

155 delay = result["delay_s"]

156 sensitivity[name] = (abs(delay - baseline_delay) / baseline_delay

157 if baseline_delay is not None and delay is not None else None)

158 converged = all(value is not None and value <= config["delay_relative_tolerance"] for value in sensitivity.values())

159 report = {

160 "versions": {"cantera": ct.__version__, "numpy": np.__version__}, "settings": config,

161 "config_sha256": hashlib.sha256(config_path.read_bytes()).hexdigest(),

162 "mechanism": {"requested": config["mechanism"], "phase": gas.name,

163 "source_sha256": hashlib.sha256(source.read_bytes()).hexdigest() if source else None,

164 "snapshot_sha256": hashlib.sha256(snapshot).hexdigest(),

165 "species_count": gas.n_species, "reaction_count": gas.n_reactions,

166 "initial_normalized_mole_fractions": {name: float(x) for name, x in zip(gas.species_names, gas.X) if x > 0}},

167 "definition": "Time of global maximum finite-difference dT/dt on a uniform output grid; requires the specified temperature rise and an interior maximum",

168 "runs": results, "delay_relative_changes": sensitivity,

169 "numerically_resolved": converged and all(r["conservation"]["passed"] for r in results.values()),

170 "mechanism_experimental_validity": "not_assessed",

171 }

172 output_dir.mkdir(parents=True)

173 (output_dir / "mechanism.yaml").write_bytes(snapshot)

174 (output_dir / "report.json").write_text(json.dumps(report, indent=2, allow_nan=False) + "\n", encoding="utf-8")

175 for name, history in histories.items():

176 with (output_dir / f"{name}.csv").open("w", newline="", encoding="utf-8") as handle:

177 writer = csv.writer(handle)

178 writer.writerow(["time_s", "temperature_k", "pressure_pa", "volume_m3", "mass_kg",

179 "internal_energy_j", "enthalpy_j", *[f"X_{s}" for s in config["tracked_species"]]])

180 writer.writerows(history)

181 return report

182

183

184def main():

185 parser = argparse.ArgumentParser(description=__doc__)

186 parser.add_argument("config", type=Path)

187 parser.add_argument("output", type=Path, help="New output directory")

188 args = parser.parse_args()

189 try:

190 report = run(args.config, args.output)

191 except (ValueError, OSError, KeyError, TypeError, ct.CanteraError) as error:

192 parser.error(str(error))

193 print(json.dumps({"output": str(args.output), "numerically_resolved": report["numerically_resolved"],

194 "delay_s": report["runs"]["baseline"]["delay_s"]}))

195

196

197if __name__ == "__main__":

198 main()