diff --git a/include/cantera/numerics/SteadyStateSystem.h b/include/cantera/numerics/SteadyStateSystem.h index ad07183b506..e3e346b4b89 100644 --- a/include/cantera/numerics/SteadyStateSystem.h +++ b/include/cantera/numerics/SteadyStateSystem.h @@ -6,6 +6,8 @@ #ifndef CT_STEADYSTATESYSTEM_H #define CT_STEADYSTATESYSTEM_H +#include + #include "cantera/base/ct_defs.h" #include "SystemJacobian.h" @@ -15,6 +17,20 @@ namespace Cantera class Func1; class MultiNewton; +//! Error thrown when time stepping cannot proceed and the steady-state solver +//! should be given a chance to recover by other means. +class TimeStepError : public CanteraError +{ +public: + template + TimeStepError(const string& func, const string& msg, const Args&... args) : + CanteraError(func, msg, args...) {} + + string getClass() const override { + return "TimeStepError"; + } +}; + //! Base class for representing a system of differential-algebraic equations and solving //! for its steady-state response. //! @@ -73,6 +89,11 @@ class SteadyStateSystem //! x. On return, array r contains the steady-state residual values. double ssnorm(span x, span r); + //! Transient max norm (infinity norm) of the residual evaluated using solution + //! x and the current timestep (rdt). On return, array r contains the + //! transient residual values. + double tsnorm(span x, span r); + //! Total solution vector length; size_t size() const { return m_size; @@ -194,6 +215,49 @@ class SteadyStateSystem m_tfactor = tfactor; } + //! Set the growth factor used after successful timesteps when the Jacobian is + //! re-used. + //! + //! This factor is used directly by the `fixed-growth` strategy, and as the + //! accepted growth factor or upper bound for the other named strategies. + //! + //! The default value is 1.5, matching historical behavior. + //! @param tfactor Finite growth factor applied to successful timesteps; + //! must be >= 1.0. + void setTimeStepGrowthFactor(double tfactor) { + if (!std::isfinite(tfactor) || tfactor < 1.0) { + throw CanteraError("SteadyStateSystem::setTimeStepGrowthFactor", + "Time step growth factor must be finite and >= 1.0. Got {}.", + tfactor); + } + m_tstep_growth = tfactor; + } + + //! Get the successful-step time step growth factor. + double timeStepGrowthFactor() const { + return m_tstep_growth; + } + + //! Set the strategy used to grow the timestep after a successful step that + //! reuses the current Jacobian. + //! + //! Available options are: + //! - `fixed-growth`: Always apply timeStepGrowthFactor(). + //! - `steady-norm`: Apply timeStepGrowthFactor() only if the steady-state + //! residual norm decreases. + //! - `transient-residual`: Apply timeStepGrowthFactor() only if the transient + //! residual norm decreases. + //! - `residual-ratio`: Scale the growth factor based on transient residual + //! improvement, capped by timeStepGrowthFactor(). + //! - `newton-iterations`: Apply timeStepGrowthFactor() only if the most recent + //! Newton solve used at most 3 iterations. + //! + //! The default strategy is `fixed-growth`, which matches historical behavior. + void setTimeStepGrowthStrategy(const string& strategy); + + //! Get the configured timestep growth strategy. + string timeStepGrowthStrategy() const; + //! Set the maximum number of timeteps allowed before successful steady-state solve void setMaxTimeStepCount(int nmax) { m_nsteps_max = nmax; @@ -251,10 +315,27 @@ class SteadyStateSystem virtual void clearDebugFile() {} protected: + enum class TimeStepGrowthStrategy { + fixed, + steadyNorm, + transientResidual, + residualRatio, + newtonIterations + }; + //! Evaluate the steady-state Jacobian, accessible via linearSolver() //! @param[in] x Current state vector, length size() void evalSSJacobian(span x); + static TimeStepGrowthStrategy parseTimeStepGrowthStrategy(const string& strategy); + static string timeStepGrowthStrategyName(TimeStepGrowthStrategy strategy); + + //! Determine the timestep growth factor after a successful step. + //! + //! Called only when a successful step reuses the current Jacobian. + double calculateTimeStepGrowthFactor(span x_before, + span x_after); + //! Array of number of steps to take after each unsuccessful steady-state solve //! before re-attempting the steady-state solution. For subsequent time stepping //! calls, the final value is reused. See setTimeStep(). @@ -267,6 +348,15 @@ class SteadyStateSystem //! Factor time step is multiplied by if time stepping fails ( < 1 ) double m_tfactor = 0.5; + //! Growth factor for successful steps that reuse the Jacobian. + //! + //! Used directly for `fixed-growth`, and as the base / cap value for the + //! other named growth strategies. + double m_tstep_growth = 1.5; + + //! Selected strategy for successful-step growth. + TimeStepGrowthStrategy m_tstep_growth_strategy = TimeStepGrowthStrategy::fixed; + shared_ptr> m_state; //!< Solution vector //! Work array used to hold the residual or the new solution @@ -314,6 +404,7 @@ class SteadyStateSystem double m_jacobianRelPerturb = 1e-5; //! Absolute perturbation of each component in finite difference Jacobian double m_jacobianAbsPerturb = 1e-10; + }; } diff --git a/include/cantera/oneD/MultiNewton.h b/include/cantera/oneD/MultiNewton.h index 4230e19848e..29fae5604e5 100644 --- a/include/cantera/oneD/MultiNewton.h +++ b/include/cantera/oneD/MultiNewton.h @@ -35,6 +35,11 @@ class MultiNewton return m_n; } + //! Number of Newton iterations taken in the most recent solve() call. + int lastIterations() const { + return m_lastIterations; + } + //! Compute the undamped Newton step. The residual function is evaluated //! at `x`, but the Jacobian is not recomputed. //! @since Starting in %Cantera 3.2, the Jacobian is accessed via the OneDim object. @@ -174,6 +179,9 @@ class MultiNewton //! Elapsed CPU time spent computing the Jacobian. double m_elapsed = 0.0; + + //! Number of Newton iterations taken in the last solve() call. + int m_lastIterations = 0; }; } diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index 7edeff66909..b99c507cc1d 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/include/cantera/oneD/Sim1D.h @@ -337,6 +337,25 @@ class Sim1D : public OneDim m_steady_callback = callback; } + //! Set the maximum number of regrid attempts after a timestep failure. + //! + //! This fallback is used during solve(loglevel, refine_grid=true). Set to `0` + //! to disable regrid-on-timestep-failure retries. + //! + //! @param nmax Maximum retry attempts; must be >= 0. + void setTimeStepRegridMax(int nmax) { + if (nmax < 0) { + throw CanteraError("Sim1D::setTimeStepRegridMax", + "Time step regrid retry count must be >= 0. Got {}.", nmax); + } + m_ts_regrid_max = nmax; + } + + //! Get the maximum number of regrid attempts after a timestep failure. + int timeStepRegridMax() const { + return m_ts_regrid_max; + } + protected: //! the solution vector after the last successful steady-state solve (stored //! before grid refinement) @@ -349,6 +368,9 @@ class Sim1D : public OneDim //! User-supplied function called after a successful steady-state solve. Func1* m_steady_callback; + //! 0 disables regrid-on-timestep-failure retries + int m_ts_regrid_max = 3; + private: //! Calls method _finalize in each domain. void finalize(); diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 48e1b2586c3..c5145e3e990 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -149,8 +149,14 @@ cdef extern from "cantera/oneD/Sim1D.h": void getResidual(double, span[double]) except +translate_exception void setJacAge(int, int) void setTimeStepFactor(double) + void setTimeStepGrowthFactor(double) except +translate_exception + double timeStepGrowthFactor() + void setTimeStepGrowthStrategy(const string&) except +translate_exception + string timeStepGrowthStrategy() void setMinTimeStep(double) void setMaxTimeStep(double) + void setTimeStepRegridMax(int) except +translate_exception + int timeStepRegridMax() void setMaxGridPoints(int, size_t) except +translate_exception size_t maxGridPoints(size_t) except +translate_exception void setGridMin(int, double) except +translate_exception diff --git a/interfaces/cython/cantera/_onedim.pyi b/interfaces/cython/cantera/_onedim.pyi index aea93f5b8da..f1999ac7c7b 100644 --- a/interfaces/cython/cantera/_onedim.pyi +++ b/interfaces/cython/cantera/_onedim.pyi @@ -327,6 +327,18 @@ class Sim1D: ) -> None: ... def set_max_jac_age(self, ss_age: int, ts_age: int) -> None: ... def set_time_step_factor(self, tfactor: float) -> None: ... + @property + def time_step_growth_factor(self) -> float: ... + @time_step_growth_factor.setter + def time_step_growth_factor(self, tfactor: float) -> None: ... + @property + def time_step_growth_strategy(self) -> str: ... + @time_step_growth_strategy.setter + def time_step_growth_strategy(self, strategy: str) -> None: ... + @property + def time_step_regrid(self) -> int: ... + @time_step_regrid.setter + def time_step_regrid(self, max_tries: int) -> None: ... def set_min_time_step(self, tsmin: float) -> None: ... def set_max_time_step(self, tsmax: float) -> None: ... @property diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 78e209d5c83..ea6e8e3b624 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -1425,11 +1425,73 @@ cdef class Sim1D: def set_time_step_factor(self, tfactor): """ - Set the factor by which the time step will be increased after a - successful step, or decreased after an unsuccessful one. + Set the factor by which the time step will be decreased after an + unsuccessful step. + + :param tfactor: + Multiplicative reduction factor applied after failed steps. """ self.sim.setTimeStepFactor(tfactor) + property time_step_growth_factor: + """ + Get/Set the factor by which the time step will be increased after a + successful step when the Jacobian is reused. + + This value is used directly by the ``"fixed-growth"`` strategy, and as + the accepted growth factor or cap value for the other growth + strategies. + + :param tfactor: + Finite growth factor >= 1.0. The default value is 1.5. + """ + def __get__(self): + return self.sim.timeStepGrowthFactor() + def __set__(self, tfactor): + self.sim.setTimeStepGrowthFactor(tfactor) + + property time_step_growth_strategy: + """ + Get/Set the strategy used to grow the time step after a successful + step that reuses the Jacobian. + + Available options are: + + ``"fixed-growth"``: + Always apply ``time_step_growth_factor``. + ``"steady-norm"``: + Apply growth only when the steady-state residual norm decreases. + ``"transient-residual"``: + Apply growth only when the transient residual norm decreases. + ``"residual-ratio"``: + Scale the growth factor based on transient residual improvement, + capped by ``time_step_growth_factor``. + ``"newton-iterations"``: + Apply growth only if the most recent Newton solve used at most + three iterations. + """ + def __get__(self): + return pystr(self.sim.timeStepGrowthStrategy()) + def __set__(self, strategy): + self.sim.setTimeStepGrowthStrategy(stringify(strategy)) + + property time_step_regrid: + """ + Get/Set the maximum number of regrid attempts after a time step + failure. + + This fallback is used by :meth:`Sim1D.solve` when ``refine_grid=True``. + Set to zero to disable regrid-on-failure retries. The default value is + 3. + + :param max_tries: + Maximum retry attempts. Values less than zero are invalid. + """ + def __get__(self): + return self.sim.timeStepRegridMax() + def __set__(self, max_tries): + self.sim.setTimeStepRegridMax(max_tries) + def set_min_time_step(self, tsmin): """ Set the minimum time step. """ self.sim.setMinTimeStep(tsmin) diff --git a/samples/python/onedim/diffusion_flame_timestep_regridding.py b/samples/python/onedim/diffusion_flame_timestep_regridding.py new file mode 100644 index 00000000000..9079b7032ee --- /dev/null +++ b/samples/python/onedim/diffusion_flame_timestep_regridding.py @@ -0,0 +1,110 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +""" +Timestep Regridding for a High-Pressure Diffusion Flame +======================================================= + +This example compares two attempts to solve the same high-pressure hydrogen / +oxygen counterflow diffusion flame on a deliberately coarse initial grid. With +regrid-on-timestep-failure disabled, the solver gives up after the first failed +timestepping sequence. Allowing a few regrid retries lets the solver add points +and recover a converged solution. + +Requires: cantera >= 4.0; matplotlib is optional for plotting + +.. tags:: Python, combustion, 1D flow, diffusion flame, strained flame, plotting +""" + +import time + +import numpy as np + +import cantera as ct + +try: + import matplotlib.pyplot as plt +except ImportError: + plt = None + + +WIDTH = 30e-3 +INITIAL_POINTS = 20 +PRESSURE = 7e6 +FUEL_MDOT = 0.3 +OXIDIZER_MDOT_FACTOR = 3.0 + + +def make_flame(regrid_max: int) -> ct.CounterflowDiffusionFlame: + gas = ct.Solution("h2o2.yaml") + flame = ct.CounterflowDiffusionFlame( + gas, grid=np.linspace(0.0, WIDTH, INITIAL_POINTS)) + flame.max_time_step_count = 200 + flame.set_refine_criteria(ratio=2.0, slope=0.06, curve=0.08, prune=0.02) + flame.time_step_regrid = regrid_max + + flame.P = PRESSURE + flame.fuel_inlet.X = "H2:1.0" + flame.fuel_inlet.T = 800.0 + flame.oxidizer_inlet.X = "O2:1.0" + flame.oxidizer_inlet.T = 711.0 + + rho_f = flame.fuel_inlet.phase.density + rho_o = flame.oxidizer_inlet.phase.density + flame.fuel_inlet.mdot = FUEL_MDOT + flame.oxidizer_inlet.mdot = (FUEL_MDOT / rho_f) * rho_o * OXIDIZER_MDOT_FACTOR + return flame + + +def solve_case(regrid_max: int) -> tuple[ct.CounterflowDiffusionFlame, bool, float, str]: + flame = make_flame(regrid_max) + flame.clear_stats() + t0 = time.perf_counter() + try: + flame.solve(loglevel=0, auto=False) + success = True + message = "" + except ct.CanteraError as err: + success = False + message = str(err).splitlines()[0] + elapsed = time.perf_counter() - t0 + return flame, success, elapsed, message + + +def summarize(flame: ct.CounterflowDiffusionFlame, label: str, success: bool, + elapsed: float, message: str) -> None: + print(f"\n{label}") + print(f" success: {success}") + print(f" wall time [s]: {elapsed:.2f}") + print(f" final grid points: {len(flame.grid)}") + print(f" timesteps: {sum(flame.time_step_stats)}") + print(f" Jacobian evaluations: {sum(flame.jacobian_count_stats)}") + if success: + print(f" peak temperature [K]: {np.max(flame.T):.1f}") + else: + print(f" error: {message}") + + +baseline, baseline_success, baseline_time, baseline_message = solve_case(regrid_max=0) +regridded, regridded_success, regridded_time, regridded_message = solve_case(regrid_max=3) + +summarize(baseline, "Regrid fallback disabled", baseline_success, + baseline_time, baseline_message) +summarize(regridded, "Regrid fallback enabled", regridded_success, + regridded_time, regridded_message) + +if regridded_success and plt is not None: + fig, ax = plt.subplots() + ax.plot(regridded.grid * 1e3, regridded.T, "o-", label="Recovered solution") + ax.set_xlabel("z [mm]") + ax.set_ylabel("Temperature [K]") + ax.set_title( + "Counterflow diffusion flame at " + f"{PRESSURE / 1e6:.1f} MPa with timestep regridding" + ) + ax.grid(True, alpha=0.3) + ax.legend() + fig.tight_layout() + plt.show() +elif regridded_success: + print("\nInstall matplotlib to plot the recovered temperature profile.") diff --git a/samples/python/onedim/flame_solver_timestep_controls.py b/samples/python/onedim/flame_solver_timestep_controls.py new file mode 100644 index 00000000000..bc354dda8a9 --- /dev/null +++ b/samples/python/onedim/flame_solver_timestep_controls.py @@ -0,0 +1,475 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +""" +Timestep Controls for 1D Flame Solvers +====================================== + +This example prints three compact comparisons for one-dimensional flame solvers: + +1. Named strategies for growing the timestep after a successful transient step. +2. Regridding retries after timestepping is exhausted on a difficult problem. +3. A harder 15 MPa case that is sensitive to both timestep budget and + timestep-growth strategy. + +The first table varies the timestep-growth strategy on a fixed-temperature free +flame. The second varies the regrid retry count on a coarse high-pressure +counterflow diffusion flame. The third shows a numerically difficult +variant where both `max_time_step_count` and `time_step_growth_strategy` +matter. + +Requires: cantera >= 4.0; matplotlib is optional for plotting + +.. tags:: Python, combustion, 1D flow, premixed flame, diffusion flame, strained flame, + plotting +""" + +from __future__ import annotations + +import time + +import numpy as np + +import cantera as ct + +try: + import matplotlib.pyplot as plt +except ImportError: + plt = None + + +GROWTH_CASE = { + "Tin": 250.0, + "pressure": ct.one_atm, + "reactants": "H2:3, O2:1, AR:10", + "width": 20e-3, +} +GROWTH_FACTOR = 2.0 +GROWTH_STRATEGIES = ( + "fixed-growth", + "steady-norm", + "transient-residual", + "residual-ratio", + "newton-iterations", +) + +REGRID_WIDTH = 30e-3 +REGRID_INITIAL_POINTS = 20 +REGRID_PRESSURE = 7e6 +REGRID_REFINE = (2.0, 0.06, 0.08, 0.02) +REGRID_CASES = ( + { + "label": "time_step_regrid = 0", + "plot_label": "0", + "regrid_max": 0, + "refine_criteria": REGRID_REFINE, + }, + { + "label": "time_step_regrid = 1", + "plot_label": "1", + "regrid_max": 1, + "refine_criteria": REGRID_REFINE, + }, + { + "label": "time_step_regrid = 3", + "plot_label": "3", + "regrid_max": 3, + "refine_criteria": REGRID_REFINE, + }, +) + +DIFFICULT_PRESSURE = 15e6 +DIFFICULT_BUDGET_CASES = ( + { + "label": "max_time_step_count = 200", + "plot_label": "200", + "max_time_step_count": 200, + "regrid_max": 3, + "growth_strategy": "fixed-growth", + "growth_factor": 1.5, + }, + { + "label": "max_time_step_count = 1250", + "plot_label": "1250", + "max_time_step_count": 1250, + "regrid_max": 3, + "growth_strategy": "fixed-growth", + "growth_factor": 1.5, + }, +) + +DIFFICULT_STRATEGY_CASES = ( + { + "label": "residual-ratio", + "plot_label": "ratio", + "max_time_step_count": 1250, + "regrid_max": 3, + "growth_strategy": "residual-ratio", + "growth_factor": 1.5, + }, + { + "label": "newton-iterations", + "plot_label": "newton", + "max_time_step_count": 1250, + "regrid_max": 3, + "growth_strategy": "newton-iterations", + "growth_factor": 1.5, + }, +) + + +def make_growth_flame() -> ct.FreeFlame: + gas = ct.Solution("h2o2.yaml") + gas.TPX = GROWTH_CASE["Tin"], GROWTH_CASE["pressure"], GROWTH_CASE["reactants"] + + flame = ct.FreeFlame(gas, width=GROWTH_CASE["width"]) + flame.flame.set_steady_tolerances(default=[1.0e-5, 1.0e-14]) + flame.flame.set_transient_tolerances(default=[1.0e-4, 1.0e-11]) + flame.inlet.T = GROWTH_CASE["Tin"] + flame.inlet.X = GROWTH_CASE["reactants"] + return flame + + +def run_growth_case(strategy: str) -> dict[str, object]: + flame = make_growth_flame() + dts = [] + flame.set_time_step_callback(lambda dt: dts.append(float(dt)) or 0) + flame.time_step_growth_factor = GROWTH_FACTOR + flame.time_step_growth_strategy = strategy + flame.clear_stats() + + flame.energy_enabled = False + t0 = time.perf_counter() + flame.solve(loglevel=0, refine_grid=False) + elapsed = time.perf_counter() - t0 + + return { + "strategy": strategy, + "accepted_steps": len(dts), + "dt_sum": float(np.sum(dts)), + "dt_min": float(np.min(dts)), + "dt_max": float(np.max(dts)), + "jacobians": int(sum(flame.jacobian_count_stats)), + "velocity": float(flame.velocity[0]), + "wall_time": elapsed, + "dts": np.asarray(dts, dtype=float), + } + + +def make_regrid_flame(regrid_max: int, + refine_criteria: tuple[float, float, float, float], + growth_strategy: str = "fixed-growth", + growth_factor: float = 1.5, + pressure: float = REGRID_PRESSURE, + max_time_step_count: int = 200 + ) -> ct.CounterflowDiffusionFlame: + gas = ct.Solution("h2o2.yaml") + flame = ct.CounterflowDiffusionFlame( + gas, grid=np.linspace(0.0, REGRID_WIDTH, REGRID_INITIAL_POINTS)) + flame.max_time_step_count = max_time_step_count + flame.set_refine_criteria( + ratio=refine_criteria[0], slope=refine_criteria[1], + curve=refine_criteria[2], prune=refine_criteria[3] + ) + flame.time_step_regrid = regrid_max + flame.time_step_growth_strategy = growth_strategy + flame.time_step_growth_factor = growth_factor + + flame.P = pressure + flame.fuel_inlet.X = "H2:1.0" + flame.fuel_inlet.T = 800.0 + flame.oxidizer_inlet.X = "O2:1.0" + flame.oxidizer_inlet.T = 711.0 + + rho_f = flame.fuel_inlet.phase.density + rho_o = flame.oxidizer_inlet.phase.density + flame.fuel_inlet.mdot = 0.3 + flame.oxidizer_inlet.mdot = (0.3 / rho_f) * rho_o * 3.0 + return flame + + +def run_regrid_case(case: dict[str, object]) -> dict[str, object]: + flame = make_regrid_flame(case["regrid_max"], case["refine_criteria"]) + flame.clear_stats() + + t0 = time.perf_counter() + try: + flame.solve(loglevel=0, auto=False) + success = True + except ct.CanteraError: + success = False + elapsed = time.perf_counter() - t0 + + return { + "label": case["label"], + "plot_label": case["plot_label"], + "success": success, + "grid_points": len(flame.grid), + "timesteps": int(sum(flame.time_step_stats)), + "jacobians": int(sum(flame.jacobian_count_stats)), + "peak_temperature": float(np.max(flame.T)) if success else None, + "wall_time": elapsed, + } + + +def run_difficult_case(case: dict[str, object]) -> dict[str, object]: + flame = make_regrid_flame( + case["regrid_max"], REGRID_REFINE, case["growth_strategy"], + case["growth_factor"], DIFFICULT_PRESSURE, case["max_time_step_count"]) + flame.clear_stats() + + t0 = time.perf_counter() + try: + flame.solve(loglevel=0, auto=False) + success = True + except ct.CanteraError: + success = False + elapsed = time.perf_counter() - t0 + + return { + "label": case["label"], + "plot_label": case["plot_label"], + "growth_strategy": case["growth_strategy"], + "growth_factor": case["growth_factor"], + "max_time_step_count": case["max_time_step_count"], + "regrid_max": case["regrid_max"], + "success": success, + "grid_points": len(flame.grid), + "timesteps": int(sum(flame.time_step_stats)), + "jacobians": int(sum(flame.jacobian_count_stats)), + "peak_temperature": float(np.max(flame.T)) if success else None, + "wall_time": elapsed, + } + + +def print_growth_summary(results: list[dict[str, object]]) -> None: + print("\nNamed timestep-growth strategies") + print( + " Case: fixed-temperature H2/O2/Ar free flame at " + f"{GROWTH_CASE['Tin']:.0f} K, {GROWTH_CASE['pressure'] / ct.one_atm:.1f} atm, " + f"width = {1e3 * GROWTH_CASE['width']:.0f} mm" + ) + print(" Varied option: `time_step_growth_strategy`") + print(f" Configured time_step_growth_factor = {GROWTH_FACTOR:.1f}") + + header = ( + f"{'Strategy':<20} {'Accepted':>8} {'sum(dt) [s]':>12} " + f"{'max(dt) [s]':>12} {'Jac':>5} {'u0 [m/s]':>12}" + ) + print(f"\n {header}") + print(f" {'-' * len(header)}") + for row in results: + print( + " " + f"{row['strategy']:<20} {row['accepted_steps']:>8d} " + f"{row['dt_sum']:>12.4e} {row['dt_max']:>12.4e} " + f"{row['jacobians']:>5d} {row['velocity']:>12.10f}" + ) + + print("\n Option reminder:") + print(" `fixed-growth`: always apply `time_step_growth_factor` after a successful step.") + print(" `steady-norm`: grow only if the steady-state residual norm decreases.") + print(" `transient-residual`: grow only if the transient residual norm decreases.") + print(" `residual-ratio`: scale growth based on transient residual improvement.") + print(" `newton-iterations`: grow only after small-iteration Newton solves.") + print(" Set `time_step_growth_factor = 1.0` to disable successful-step growth.") + + +def print_regrid_summary(results: list[dict[str, object]]) -> None: + print("\nRegrid retries after timestep failure") + print( + " Case: H2/O2 counterflow diffusion flame at " + f"{REGRID_PRESSURE / 1e6:.1f} MPa on an initial {REGRID_INITIAL_POINTS}-point grid" + ) + print(" Varied option: `time_step_regrid`") + print(" Growth-factor and growth-strategy settings are left at their defaults.") + + header = ( + f"{'Mode':<43} {'Success':>7} {'Grid':>6} " + f"{'Steps':>7} {'Jac':>5} {'Tmax [K]':>10} {'Wall [s]':>9}" + ) + print(f"\n {header}") + print(f" {'-' * len(header)}") + for row in results: + peak_temperature = ( + f"{row['peak_temperature']:>10.1f}" + if row["peak_temperature"] is not None else f"{'--':>10}" + ) + print( + " " + f"{row['label']:<43} {str(row['success']):>7} " + f"{row['grid_points']:>6d} {row['timesteps']:>7d} " + f"{row['jacobians']:>5d} {peak_temperature} {row['wall_time']:>9.2f}" + ) + + print("\n Option reminder:") + print(" `time_step_regrid = 0`: disable retry-after-regrid.") + print(" `time_step_regrid = n > 0`: allow up to `n` regrid retries after timestep failure.") + print(" If refinement criteria do not change the grid, the retry path aborts.") + + +def print_difficult_summary(budget_results: list[dict[str, object]], + strategy_results: list[dict[str, object]]) -> None: + print("\nA numerically difficult 15 MPa case") + print( + " Case: the same coarse H2/O2 counterflow diffusion flame at " + f"{DIFFICULT_PRESSURE / 1e6:.1f} MPa" + ) + print(" This case converges only after allowing more transient steps,") + print(" and then shows meaningful sensitivity to the growth heuristic.") + + header = ( + f"{'Mode':<25} {'Success':>7} {'Grid':>6} " + f"{'Steps':>7} {'Jac':>5} {'Tmax [K]':>10} {'Wall [s]':>9}" + ) + print("\n Recovery with a larger timestep budget") + print(" Varied option: `max_time_step_count` with `time_step_regrid = 3`") + print(" Growth settings are held at `fixed-growth` with factor 1.5.") + print(f"\n {header}") + print(f" {'-' * len(header)}") + for row in budget_results: + peak_temperature = ( + f"{row['peak_temperature']:>10.1f}" + if row["peak_temperature"] is not None else f"{'--':>10}" + ) + print( + " " + f"{row['label']:<25} {str(row['success']):>7} " + f"{row['grid_points']:>6d} {row['timesteps']:>7d} " + f"{row['jacobians']:>5d} {peak_temperature} {row['wall_time']:>9.2f}" + ) + + strategy_header = ( + f"{'Strategy':<18} {'Success':>7} {'Grid':>6} " + f"{'Steps':>7} {'Jac':>5} {'Tmax [K]':>10} {'Wall [s]':>9}" + ) + print("\n Growth-strategy sensitivity after recovery") + print( + " Varied option: `time_step_growth_strategy` with " + "`max_time_step_count = 1250`, `time_step_regrid = 3`, and " + "`time_step_growth_factor = 1.5`." + ) + print(f"\n {strategy_header}") + print(f" {'-' * len(strategy_header)}") + for row in strategy_results: + peak_temperature = ( + f"{row['peak_temperature']:>10.1f}" + if row["peak_temperature"] is not None else f"{'--':>10}" + ) + print( + " " + f"{row['growth_strategy']:<18} {str(row['success']):>7} " + f"{row['grid_points']:>6d} {row['timesteps']:>7d} " + f"{row['jacobians']:>5d} {peak_temperature} {row['wall_time']:>9.2f}" + ) + + print("\n Option reminder:") + print(" `max_time_step_count` limits accepted transient steps before a timestep attempt gives up.") + print(" Larger values can recover harder cases, but they increase wall time substantially.") + print(" On this 15 MPa case, adaptive growth strategies can also cut wall time substantially once the timestep budget is large enough.") + + +def plot_results(growth_results: list[dict[str, object]], + regrid_results: list[dict[str, object]], + difficult_strategy_results: list[dict[str, object]]) -> None: + fig, axes = plt.subplots(1, 3, figsize=(15, 4.8)) + + growth_colors = { + "fixed-growth": "#005f73", + "steady-norm": "#bb3e03", + "transient-residual": "#0a9396", + "residual-ratio": "#ca6702", + "newton-iterations": "#ae2012", + } + for row in growth_results: + steps = np.arange(1, row["accepted_steps"] + 1) + axes[0].semilogy( + steps, row["dts"], "o-", linewidth=1.6, markersize=4, + color=growth_colors[row["strategy"]], label=row["strategy"] + ) + axes[0].set_title("Accepted dt histories") + axes[0].set_xlabel("Accepted timestep number") + axes[0].set_ylabel("dt [s]") + axes[0].grid(True, which="both", alpha=0.3) + axes[0].legend(fontsize=8) + + colors = ["#ae2012" if not row["success"] else "#0a9396" for row in regrid_results] + bars = axes[1].bar( + [row["plot_label"] for row in regrid_results], + [row["timesteps"] for row in regrid_results], + color=colors, + ) + axes[1].set_title("Regrid retry modes") + axes[1].set_xlabel("time_step_regrid") + axes[1].set_ylabel("Accepted timesteps before exit") + axes[1].grid(True, axis="y", alpha=0.3) + for bar, row in zip(bars, regrid_results): + tag = "ok" if row["success"] else "fail" + axes[1].text( + bar.get_x() + bar.get_width() / 2, + bar.get_height() + 10, + f"grid={row['grid_points']}\n{tag}", + ha="center", + va="bottom", + fontsize=8, + ) + + colors = [ + growth_colors.get(row["growth_strategy"], "#0a9396") + if row["success"] else "#ae2012" + for row in difficult_strategy_results + ] + bars = axes[2].bar( + [row["plot_label"] for row in difficult_strategy_results], + [row["wall_time"] for row in difficult_strategy_results], + color=colors, + ) + axes[2].set_title("15 MPa heuristic sensitivity") + axes[2].set_xlabel("time_step_growth_strategy") + axes[2].set_ylabel("Wall time [s]") + axes[2].grid(True, axis="y", alpha=0.3) + for bar, row in zip(bars, difficult_strategy_results): + tag = "ok" if row["success"] else "fail" + axes[2].text( + bar.get_x() + bar.get_width() / 2, + bar.get_height() + 0.3, + f"jac={row['jacobians']}\n{tag}", + ha="center", + va="bottom", + fontsize=8, + ) + + fig.tight_layout() + if "agg" not in plt.get_backend().lower(): + plt.show() + + +def main() -> None: + growth_results = [run_growth_case(strategy) for strategy in GROWTH_STRATEGIES] + regrid_results = [run_regrid_case(case) for case in REGRID_CASES] + difficult_budget_results = [run_difficult_case(case) for case in DIFFICULT_BUDGET_CASES] + recovered_baseline = next( + row for row in difficult_budget_results if row["max_time_step_count"] == 1250 + ) + difficult_strategy_results = [ + { + **recovered_baseline, + "label": "fixed-growth", + "plot_label": "fixed", + }, + *[run_difficult_case(case) for case in DIFFICULT_STRATEGY_CASES], + ] + + print_growth_summary(growth_results) + print_regrid_summary(regrid_results) + print_difficult_summary(difficult_budget_results, difficult_strategy_results) + + if plt is not None: + plot_results( + growth_results, regrid_results, difficult_strategy_results) + else: + print("\nInstall matplotlib to plot the timestep histories and comparison tables.") + + +if __name__ == "__main__": + main() diff --git a/src/numerics/SteadyStateSystem.cpp b/src/numerics/SteadyStateSystem.cpp index b8dd85ceb5b..e4ffd0b41d3 100644 --- a/src/numerics/SteadyStateSystem.cpp +++ b/src/numerics/SteadyStateSystem.cpp @@ -108,6 +108,97 @@ void SteadyStateSystem::solve(int loglevel) } } +SteadyStateSystem::TimeStepGrowthStrategy +SteadyStateSystem::parseTimeStepGrowthStrategy(const string& strategy) +{ + if (strategy == "fixed-growth") { + return TimeStepGrowthStrategy::fixed; + } else if (strategy == "steady-norm") { + return TimeStepGrowthStrategy::steadyNorm; + } else if (strategy == "transient-residual") { + return TimeStepGrowthStrategy::transientResidual; + } else if (strategy == "residual-ratio") { + return TimeStepGrowthStrategy::residualRatio; + } else if (strategy == "newton-iterations") { + return TimeStepGrowthStrategy::newtonIterations; + } + throw CanteraError("SteadyStateSystem::setTimeStepGrowthStrategy", + "Unknown time step growth strategy '{}'; must be one of " + "'fixed-growth', 'steady-norm', 'transient-residual', " + "'residual-ratio', or 'newton-iterations'.", strategy); +} + +string SteadyStateSystem::timeStepGrowthStrategyName(TimeStepGrowthStrategy strategy) +{ + switch (strategy) { + case TimeStepGrowthStrategy::fixed: + return "fixed-growth"; + case TimeStepGrowthStrategy::steadyNorm: + return "steady-norm"; + case TimeStepGrowthStrategy::transientResidual: + return "transient-residual"; + case TimeStepGrowthStrategy::residualRatio: + return "residual-ratio"; + case TimeStepGrowthStrategy::newtonIterations: + return "newton-iterations"; + } + throw CanteraError("SteadyStateSystem::timeStepGrowthStrategyName", + "Unknown time step growth strategy."); +} + +void SteadyStateSystem::setTimeStepGrowthStrategy(const string& strategy) +{ + m_tstep_growth_strategy = parseTimeStepGrowthStrategy(strategy); +} + +string SteadyStateSystem::timeStepGrowthStrategy() const +{ + return timeStepGrowthStrategyName(m_tstep_growth_strategy); +} + +double SteadyStateSystem::calculateTimeStepGrowthFactor(span x_before, + span x_after) +{ + if (m_tstep_growth_strategy == TimeStepGrowthStrategy::fixed) { + return m_tstep_growth; + } + + m_work1.resize(m_size); + const double growth = m_tstep_growth; + + switch (m_tstep_growth_strategy) { + case TimeStepGrowthStrategy::fixed: + return growth; + case TimeStepGrowthStrategy::steadyNorm: { + double ss_before = ssnorm(x_before, m_work1); + double ss_after = ssnorm(x_after, m_work1); + return (ss_after < ss_before) ? growth : 1.0; + } + case TimeStepGrowthStrategy::transientResidual: { + double ts_before = tsnorm(x_before, m_work1); + double ts_after = tsnorm(x_after, m_work1); + return (ts_after < ts_before) ? growth : 1.0; + } + case TimeStepGrowthStrategy::residualRatio: { + double ts_before = tsnorm(x_before, m_work1); + double ts_after = tsnorm(x_after, m_work1); + if (!(ts_after > 0.0) || !(ts_before > ts_after)) { + return 1.0; + } + const double exponent = 0.2; + double ratio = ts_before / ts_after; + double factor = std::pow(ratio, exponent); + return std::min(growth, std::max(1.0, factor)); + } + case TimeStepGrowthStrategy::newtonIterations: { + const int max_iters_for_growth = 3; + return (newton().lastIterations() <= max_iters_for_growth) ? growth : 1.0; + } + } + throw CanteraError("SteadyStateSystem::calculateTimeStepGrowthFactor", + "Unknown time step growth strategy '{}'.", timeStepGrowthStrategy()); +} + double SteadyStateSystem::timeStep(int nsteps, double dt, span x, span r, int loglevel) { @@ -149,17 +240,16 @@ double SteadyStateSystem::timeStep(int nsteps, double dt, span x, successiveFailures = 0; m_nsteps++; n += 1; - copy(r.begin(), r.end(), x.begin()); - // No Jacobian evaluations were performed, so a larger timestep can be used if (m_jac->nEvals() == j0) { - dt *= 1.5; + dt *= calculateTimeStepGrowthFactor(x, r); } + copy(r.begin(), r.end(), x.begin()); if (m_time_step_callback) { m_time_step_callback->eval(dt); } dt = std::min(dt, m_tmax); if (m_nsteps >= m_nsteps_max) { - throw CanteraError("OneDim::timeStep", + throw TimeStepError("SteadyStateSystem::timeStep", "Took maximum number of timesteps allowed ({}) without " "reaching steady-state solution.", m_nsteps_max); } @@ -180,9 +270,8 @@ double SteadyStateSystem::timeStep(int nsteps, double dt, span x, debuglog("--> Reducing timestep", loglevel); dt *= m_tfactor; if (dt < m_tmin) { - string err_msg = fmt::format( + throw TimeStepError("SteadyStateSystem::timeStep", "Time integration failed. Minimum timestep ({}) reached.", m_tmin); - throw CanteraError("OneDim::timeStep", err_msg); } } } @@ -213,6 +302,16 @@ double SteadyStateSystem::ssnorm(span x, span r) return ss; } +double SteadyStateSystem::tsnorm(span x, span r) +{ + eval(x, r, m_rdt, 0); + double ts = 0.0; + for (size_t i = 0; i < m_size; i++) { + ts = std::max(fabs(r[i]), ts); + } + return ts; +} + void SteadyStateSystem::setTimeStep(double stepsize, span tsteps) { m_tstep = stepsize; @@ -271,7 +370,9 @@ void SteadyStateSystem::setSteadyMode() } m_rdt = 0.0; - m_jac->updateTransient(m_rdt, m_mask); + if (m_jac_ok) { + m_jac->updateTransient(m_rdt, m_mask); + } } void SteadyStateSystem::setJacAge(int ss_age, int ts_age) diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index 46b65f28099..80b71e3f1e8 100644 --- a/src/oneD/MultiNewton.cpp +++ b/src/oneD/MultiNewton.cpp @@ -208,11 +208,13 @@ int MultiNewton::solve(span x0, span x1, double s1=1.e30; copy(x0.begin(), x0.end(), m_x.begin()); + m_lastIterations = 0; double rdt = r.rdt(); int nJacReeval = 0; auto jac = r.linearSolver(); while (true) { + m_lastIterations++; // Check whether the Jacobian should be re-evaluated. if (jac->age() > m_maxAge) { if (loglevel > 1) { diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index 281e0f877cb..80402b2926c 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -348,8 +348,31 @@ void Sim1D::solve(int loglevel, bool refine_grid) clearDebugFile(); } + int retries = 0; while (new_points > 0) { - SteadyStateSystem::solve(loglevel); + if (refine_grid && retries < m_ts_regrid_max) { + try { + SteadyStateSystem::solve(loglevel); + } catch (TimeStepError&) { + if (loglevel > 0) { + writelog("\nTime stepping failed; attempting to refine the grid and retry " + "({}/{})...\n", retries+1, m_ts_regrid_max); + } + int regrid_result = refine(loglevel); + if (regrid_result != 0) { + retries++; + new_points = 1; + continue; + } + if (loglevel > 0) { + writelog("Regrid retry aborted: grid was unchanged.\n"); + } + throw; + } + } else { + SteadyStateSystem::solve(loglevel); + } + if (loglevel > 0) { writelog("\nNewton steady-state solve succeeded.\n\n"); writelog("Problem solved on ["); diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index 07fbadb1b3e..ed340a5906b 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -139,6 +139,28 @@ def test_tolerances(self): assert rtol_ss == set((5e-3, 3e-4, 7e-7)) assert rtol_ts == set((6e-3, 4e-4, 2e-7)) + def test_time_step_control_option_validation(self): + gas = ct.Solution("h2o2.yaml") + left = ct.Inlet1D(gas) + flame = ct.FreeFlow(gas) + right = ct.Outlet1D(gas) + sim = ct.Sim1D((left, flame, right)) + + with pytest.raises(ct.CanteraError, match=">= 1.0"): + sim.time_step_growth_factor = 0.9 + + with pytest.raises(ct.CanteraError, match="finite and >= 1.0"): + sim.time_step_growth_factor = float("inf") + + with pytest.raises(ct.CanteraError, match="finite and >= 1.0"): + sim.time_step_growth_factor = float("nan") + + with pytest.raises(ct.CanteraError, match="must be one of"): + sim.time_step_growth_strategy = "not-a-strategy" + + with pytest.raises(ct.CanteraError, match=">= 0"): + sim.time_step_regrid = -1 + def test_switch_transport(self): gas = ct.Solution('h2o2.yaml') gas.set_equivalence_ratio(0.9, 'H2', 'O2:0.21, N2:0.79') @@ -209,6 +231,15 @@ def solve_fixed_T(self, loglevel=0): assert not self.sim.energy_enabled + def solve_fixed_T_with_growth(self, growth_factor=1.5, strategy="fixed-growth"): + dts = [] + self.sim.set_time_step_callback(lambda dt: dts.append(float(dt)) or 0) + self.sim.time_step_growth_factor = growth_factor + self.sim.time_step_growth_strategy = strategy + self.sim.clear_stats() + self.solve_fixed_T() + return dts + def solve_mix(self, ratio=3.0, slope=0.3, curve=0.2, prune=0.0, refine=True, auto=False): # Solve with the energy equation enabled @@ -502,6 +533,59 @@ def test_timestep_jacobian_limits(self, capsys): assert "Attempt 3 timesteps" in out assert "Attempt 10 timesteps" in out + def test_time_step_growth_behavior(self): + reactants = "H2:3, O2:1, AR:10" + case = dict(p=ct.one_atm, Tin=250, reactants=reactants, width=0.02) + + self.create_sim(**case) + assert self.sim.time_step_growth_factor == approx(1.5) + assert self.sim.time_step_growth_strategy == "fixed-growth" + self.sim.time_step_growth_factor = 1.2 + assert self.sim.time_step_growth_factor == approx(1.2) + self.sim.time_step_growth_strategy = "newton-iterations" + assert self.sim.time_step_growth_strategy == "newton-iterations" + no_growth_dts = self.solve_fixed_T_with_growth(growth_factor=1.0) + reference_velocity = self.sim.velocity[0] + + self.create_sim(**case) + fixed_growth_dts = self.solve_fixed_T_with_growth(growth_factor=2.0) + assert self.sim.velocity[0] == approx(reference_velocity, rel=1e-8, abs=1e-12) + + self.create_sim(**case) + steady_norm_dts = self.solve_fixed_T_with_growth( + growth_factor=2.0, strategy="steady-norm") + assert self.sim.velocity[0] == approx(reference_velocity, rel=1e-8, abs=1e-12) + + self.create_sim(**case) + transient_residual_dts = self.solve_fixed_T_with_growth( + growth_factor=2.0, strategy="transient-residual") + assert self.sim.velocity[0] == approx(reference_velocity, rel=1e-8, abs=1e-12) + + self.create_sim(**case) + residual_ratio_dts = self.solve_fixed_T_with_growth( + growth_factor=2.0, strategy="residual-ratio") + assert self.sim.velocity[0] == approx(reference_velocity, rel=1e-8, abs=1e-12) + + self.create_sim(**case) + newton_iterations_dts = self.solve_fixed_T_with_growth( + growth_factor=2.0, strategy="newton-iterations") + assert self.sim.velocity[0] == approx(reference_velocity, rel=1e-8, abs=1e-12) + + assert len(no_growth_dts) > 0 + assert no_growth_dts == approx([no_growth_dts[0]] * len(no_growth_dts)) + assert max(fixed_growth_dts) > 100 * min(fixed_growth_dts) + assert transient_residual_dts == fixed_growth_dts + assert steady_norm_dts != fixed_growth_dts + assert sum(steady_norm_dts) < sum(fixed_growth_dts) + assert max(steady_norm_dts) < max(fixed_growth_dts) + assert residual_ratio_dts != fixed_growth_dts + assert sum(steady_norm_dts) < sum(residual_ratio_dts) < sum(fixed_growth_dts) + assert max(residual_ratio_dts) < max(fixed_growth_dts) + assert newton_iterations_dts != fixed_growth_dts + assert sum(newton_iterations_dts) < sum(steady_norm_dts) + assert min(newton_iterations_dts) < min(fixed_growth_dts) + assert max(newton_iterations_dts) < max(steady_norm_dts) + def test_mixture_averaged_case1(self): self.run_mix(phi=0.65, T=300, width=0.03, p=1.0, refine=True) @@ -1068,6 +1152,28 @@ def solve_mix(self, ratio=3.0, slope=0.1, curve=0.12, prune=0.0): assert arr.radial_pressure_gradient[0] == radial_lambda[0] assert getattr(arr, "radial-pressure-gradient")[0] == radial_lambda[0] + def make_high_pressure_regrid_case(self, pressure, regrid_max, + refine_criteria=(2.0, 0.06, 0.08, 0.02)): + gas = ct.Solution("h2o2.yaml") + flame = ct.CounterflowDiffusionFlame(gas, grid=np.linspace(0.0, 30e-3, 20)) + flame.max_time_step_count = 200 + flame.set_refine_criteria( + ratio=refine_criteria[0], slope=refine_criteria[1], + curve=refine_criteria[2], prune=refine_criteria[3] + ) + flame.time_step_regrid = regrid_max + flame.P = pressure + flame.fuel_inlet.X = "H2:1.0" + flame.fuel_inlet.T = 800.0 + flame.oxidizer_inlet.X = "O2:1.0" + flame.oxidizer_inlet.T = 711.0 + + rho_f = flame.fuel_inlet.phase.density + rho_o = flame.oxidizer_inlet.phase.density + flame.fuel_inlet.mdot = 0.3 + flame.oxidizer_inlet.mdot = (0.3 / rho_f) * rho_o * 3.0 + return flame + @pytest.mark.slow_test def test_mixture_averaged(self, request, test_data_path): referenceFile = "DiffusionFlameTest-h2-mix.csv" @@ -1136,6 +1242,43 @@ def time_step_func(dt): rtol=1e-2, atol=1e-8, xtol=1e-2) assert not bad, bad + @pytest.mark.slow_test + def test_time_step_regrid_high_pressure(self, capsys): + flame = self.make_high_pressure_regrid_case(7e6, 0) + with pytest.raises(ct.CanteraError): + flame.solve(loglevel=0, auto=False) + assert sum(flame.time_step_stats) == 200 + + flame = self.make_high_pressure_regrid_case(7e6, 1) + assert flame.time_step_regrid == 1 + with pytest.raises(ct.CanteraError): + flame.solve(loglevel=0, auto=False) + assert sum(flame.time_step_stats) == 400 + assert len(flame.grid) > 20 + + flame = self.make_high_pressure_regrid_case(7e6, 3) + flame.solve(loglevel=1, auto=False) + out = capsys.readouterr().out + + assert "Time stepping failed; attempting to refine the grid and retry" in out + assert out.count("Time stepping failed; attempting to refine the grid and retry") == 3 + assert sum(flame.time_step_stats) > 200 + assert len(flame.grid) > 20 + assert np.max(flame.T) > 3000 + + @pytest.mark.slow_test + def test_time_step_regrid_unchanged_grid_aborts(self, capsys): + flame = self.make_high_pressure_regrid_case( + 7e6, 3, refine_criteria=(1000.0, 1.0, 1.0, 0.0)) + with pytest.raises(ct.CanteraError): + flame.solve(loglevel=1, auto=False) + out = capsys.readouterr().out + + assert "Time stepping failed; attempting to refine the grid and retry" in out + assert "Regrid retry aborted: grid was unchanged." in out + assert sum(flame.time_step_stats) == 200 + assert len(flame.grid) == 20 + def run_extinction(self, mdot_fuel, mdot_ox, T_ox, width, P): self.create_sim(fuel='H2:1.0', oxidizer='O2:1.0', p=ct.one_atm*P, mdot_fuel=mdot_fuel, mdot_ox=mdot_ox, T_ox=T_ox, width=width)