H2-air flame model#

The examples include a self-contained, pure-Python model of a two-dimensional premixed H2-air flame. It provides a moderately challenging benchmark for romtools workflows while requiring only NumPy and SciPy.

The core solver and romtools QoI wrapper live in examples/models. A quick wrapper smoke test can be run from the repository root with

python examples/models/h2_air_flame_model.py

Problem definition#

The state contains hydrogen, oxygen, and water mass fractions together with a nondimensional temperature,

\[\boldsymbol{u} = [Y_{H_2}, Y_{O_2}, Y_{H_2O}, \theta]^T, \qquad \theta = T / T_{\mathrm{ref}},\]

on a rectangular domain of size 1.8 cm by 0.9 cm. The semi-discrete model combines advection, diffusion, and an Arrhenius reaction source. A hot premixed inlet occupies the middle third of the left boundary, with ambient Dirichlet data above and below it and homogeneous-Neumann conditions on the remaining boundaries.

The four model parameters are

kappa

Diffusivity.

scaled_activation_energy

Activation energy divided by 1000. A representative value is 8.0, corresponding to an activation energy of 8000.

beta_x and beta_y

Components of the advection velocity in cm/s. beta_x is restricted to nonnegative values because the benchmark inlet and x-upwind stencil assume left-to-right flow.

A representative parameter vector is

\[(\kappa, E/1000, \beta_x, \beta_y) = (2, 8, 40, 7).\]

The solver uses Crank-Nicolson time integration, central differences for diffusion, a second-order backward upwind stencil in the x direction, and an upwind stencil in the y direction. Newton systems are assembled as SciPy sparse matrices. By default, each Newton system uses the exact assembled Jacobian and is solved with GMRES, using one LU factorization of the fixed linear Crank-Nicolson operator as a preconditioner for the whole transient solve. The linear tolerance is tighter than the Newton tolerance; a GMRES solve that does not meet its true residual tolerance automatically falls back to a direct full-Jacobian solve. Set linear_solver="direct" on H2AirFlame to use direct solves throughout.

Direct use#

The core solver is independent of romtools workflow interfaces. From examples/models:

from h2_air_flame import H2AirFlame

model = H2AirFlame(
    nx=25,
    ny=13,
    dt=1.0e-4,
    t_end=5.0e-3,
    snapshot_stride=10,
)
states, times = model.solve(
    kappa=2.0,
    scaled_activation_energy=8.0,
    beta_x=40.0,
    beta_y=7.0,
)

states has shape (num_snapshots, 4, nx, ny). The temperature field is stored as theta = T / 300 K so that the four state components remain on a more useful scale for reduced-order modeling. Dimensional temperature is therefore 300 * states[:, 3].

Default solution plot#

The figure below shows the final temperature field for the default lightweight configuration and the representative parameter vector

\[(\kappa, E/1000, \beta_x, \beta_y) = (2, 8, 40, 7).\]

This corresponds to nx = 25, ny = 13, dt = 1e-4, and t_end = 5e-3.

Final temperature field for the default H2-air flame configuration

Final temperature field at t = 0.005 s for the default configuration.#

QoI model for romtools workflows#

H2AirFlameQoiModel provides the standard romtools model protocol. It uses regularly strided temperature sensors as the QoI and writes both the QoI and full state history to each run directory. Saving the state history allows the same FOM evaluations to be reused later by state-based ROM builders.

from h2_air_flame_model import H2AirFlameQoiModel

qoi_model = H2AirFlameQoiModel(
    nx=25,
    ny=13,
    dt=1.0e-4,
    t_end=5.0e-3,
    snapshot_stride=10,
    spatial_sensor_stride=4,
)

sample = {
    "kappa": 2.0,
    "scaled_activation_energy": 8.0,
    "beta_x": 40.0,
    "beta_y": 7.0,
}

qoi_model.populate_run_directory("run_0", sample)
qoi_model.run_model("run_0", sample)
qoi = qoi_model.compute_qoi("run_0", sample)

This wrapper can be passed directly to romtools workflows that consume a QoI model.

Computational cost#

The default constructor values are intended to be lightweight enough for interactive experimentation. Increasing the spatial resolution or time horizon, or reducing the time step, increases the cost of the implicit nonlinear solve. Benchmark studies should therefore report model-evaluation counts together with wall-clock cost when comparing algorithms.

Implementation#

Core flame solver:

  1"""Pure-Python two-dimensional premixed H2-air flame model.
  2
  3Solves a transient advection-diffusion-reaction system for H2, O2, H2O mass
  4fractions and a nondimensional temperature using Crank-Nicolson time
  5integration and Newton iterations.
  6"""
  7
  8from typing import Tuple
  9
 10import numpy as np
 11import scipy.sparse
 12import scipy.sparse.linalg
 13
 14
 15class H2AirFlame:
 16    """Two-dimensional premixed H2-air advection-diffusion-reaction model.
 17
 18    The state ordering is ``[Y_H2, Y_O2, Y_H2O, theta]``, where
 19    ``theta = T / temperature_scale``. The four model parameters are the
 20    diffusivity ``kappa``, scaled activation energy ``E / 1000``, and the two
 21    advection-velocity components ``beta_x`` and ``beta_y``.
 22
 23    Parameters
 24    ----------
 25    nx, ny:
 26        Number of nodal points in the x and y directions, including boundary
 27        points.
 28    dt:
 29        Time step in seconds.
 30    t_end:
 31        Final time in seconds. ``t_end / dt`` must be an integer.
 32    snapshot_stride:
 33        Number of time steps between stored snapshots. The initial and final
 34        states are always stored.
 35    newton_relative_tolerance, newton_absolute_tolerance:
 36        Newton convergence tolerances for each implicit time step.
 37    newton_max_iterations:
 38        Maximum Newton iterations per time step.
 39    linear_solver:
 40        ``"gmres"`` uses the exact Newton Jacobian with restarted GMRES and a
 41        single LU factorization of the fixed linear Crank--Nicolson operator as
 42        a preconditioner. ``"direct"`` factors the full Newton Jacobian at
 43        every Newton iteration.
 44    gmres_relative_tolerance, gmres_absolute_tolerance:
 45        Linear-solve tolerances used by the GMRES path. They should be tighter
 46        than the Newton tolerance so that inexact linear solves do not limit
 47        nonlinear convergence.
 48    """
 49
 50    field_names = ("Y_H2", "Y_O2", "Y_H2O", "theta")
 51    parameter_names = ("kappa", "scaled_activation_energy", "beta_x", "beta_y")
 52
 53    def __init__(
 54        self,
 55        nx: int = 25,
 56        ny: int = 13,
 57        dt: float = 1.0e-4,
 58        t_end: float = 5.0e-3,
 59        snapshot_stride: int = 10,
 60        newton_relative_tolerance: float = 1.0e-6,
 61        newton_absolute_tolerance: float = 1.0e-10,
 62        newton_max_iterations: int = 80,
 63        linear_solver: str = "gmres",
 64        gmres_relative_tolerance: float = 1.0e-10,
 65        gmres_absolute_tolerance: float = 0.0,
 66        gmres_restart: int = 30,
 67        gmres_max_iterations: int = 100,
 68    ) -> None:
 69        if nx < 4:
 70            raise ValueError("nx must be at least 4 for the x-upwind stencil")
 71        if ny < 3:
 72            raise ValueError("ny must be at least 3")
 73        if dt <= 0.0:
 74            raise ValueError("dt must be positive")
 75        if t_end <= 0.0:
 76            raise ValueError("t_end must be positive")
 77        if snapshot_stride < 1:
 78            raise ValueError("snapshot_stride must be at least 1")
 79        if newton_relative_tolerance <= 0.0 or newton_absolute_tolerance < 0.0:
 80            raise ValueError(
 81                "Newton tolerances must be nonnegative and rtol must be positive"
 82            )
 83        if newton_max_iterations < 1:
 84            raise ValueError("newton_max_iterations must be at least 1")
 85        if linear_solver not in ("gmres", "direct"):
 86            raise ValueError("linear_solver must be 'gmres' or 'direct'")
 87        if (
 88            gmres_relative_tolerance <= 0.0
 89            or gmres_absolute_tolerance < 0.0
 90        ):
 91            raise ValueError("GMRES tolerances must be nonnegative and rtol positive")
 92        if gmres_restart < 1 or gmres_max_iterations < 1:
 93            raise ValueError("GMRES iteration limits must be at least 1")
 94
 95        n_steps_float = t_end / dt
 96        n_steps = int(round(n_steps_float))
 97        if not np.isclose(n_steps_float, n_steps, rtol=0.0, atol=1.0e-12):
 98            raise ValueError("t_end must be an integer multiple of dt")
 99
100        self.nx = int(nx)
101        self.ny = int(ny)
102        self.dt = float(dt)
103        self.t_end = float(t_end)
104        self.snapshot_stride = int(snapshot_stride)
105        self.newton_relative_tolerance = float(newton_relative_tolerance)
106        self.newton_absolute_tolerance = float(newton_absolute_tolerance)
107        self.newton_max_iterations = int(newton_max_iterations)
108        self.linear_solver = linear_solver
109        self.gmres_relative_tolerance = float(gmres_relative_tolerance)
110        self.gmres_absolute_tolerance = float(gmres_absolute_tolerance)
111        self.gmres_restart = int(gmres_restart)
112        self.gmres_max_iterations = int(gmres_max_iterations)
113        self.linear_solver_stats = {}
114        self.num_steps = n_steps
115
116        self.length_x = 1.8  # cm
117        self.length_y = 0.9  # cm
118        self.x = np.linspace(0.0, self.length_x, self.nx)
119        self.y = np.linspace(0.0, self.length_y, self.ny)
120        self.dx = self.x[1] - self.x[0]
121        self.dy = self.y[1] - self.y[0]
122
123        self.temperature_scale = 300.0  # K
124        self.reaction_prefactor = 4.464956898694298e12
125        self.molecular_weights = np.array([2.016, 31.9, 18.0])  # g/mol
126        self.density = 1.39e-3  # g/cm^3
127        self.heat_of_reaction_scaled = 9800.0 / self.temperature_scale
128        self.gas_constant_times_temperature_scale = 8.314472 * self.temperature_scale
129        self.stoichiometric_coefficients = np.array([2.0, 1.0, -2.0])
130
131        self.ambient_temperature = 300.0
132        self.inlet_temperature = 950.0
133        self.inlet_h2_mass_fraction = 0.0282
134        self.inlet_o2_mass_fraction = 0.2259
135        self.inlet_h2o_mass_fraction = 0.0
136
137        self._gamma_1 = self.y >= 0.6
138        self._gamma_2 = np.logical_and(self.y > 0.3, self.y < 0.6)
139        self._gamma_3 = self.y <= 0.3
140        self._reaction_mask = np.ones((self.nx, self.ny))
141        self._reaction_mask[0, :] = 0.0
142        self._reaction_mask[-1, :] = 0.0
143        self._reaction_mask[:, 0] = 0.0
144        self._reaction_mask[:, -1] = 0.0
145
146    def initial_state(self) -> np.ndarray:
147        """Return the initial state with boundary conditions applied."""
148        state = np.zeros((4, self.nx, self.ny))
149        state[3, :, :] = self.ambient_temperature / self.temperature_scale
150        return self.apply_boundary_conditions(state)
151
152    def apply_boundary_conditions(self, state: np.ndarray) -> np.ndarray:
153        """Apply Dirichlet inlet and homogeneous-Neumann outlet/wall conditions."""
154        if state.shape != (4, self.nx, self.ny):
155            raise ValueError(
156                "state must have shape (4, nx, ny); got {}".format(state.shape)
157            )
158
159        state[:, :, 0] = state[:, :, 1]
160        state[:, :, -1] = state[:, :, -2]
161        state[:, -1, :] = state[:, -2, :]
162
163        state[0, 0, self._gamma_1] = 0.0
164        state[1, 0, self._gamma_1] = 0.0
165        state[2, 0, self._gamma_1] = 0.0
166        state[3, 0, self._gamma_1] = self.ambient_temperature / self.temperature_scale
167
168        state[0, 0, self._gamma_2] = self.inlet_h2_mass_fraction
169        state[1, 0, self._gamma_2] = self.inlet_o2_mass_fraction
170        state[2, 0, self._gamma_2] = self.inlet_h2o_mass_fraction
171        state[3, 0, self._gamma_2] = self.inlet_temperature / self.temperature_scale
172
173        state[0, 0, self._gamma_3] = 0.0
174        state[1, 0, self._gamma_3] = 0.0
175        state[2, 0, self._gamma_3] = 0.0
176        state[3, 0, self._gamma_3] = self.ambient_temperature / self.temperature_scale
177        return state
178
179    def _validate_parameters(
180        self,
181        kappa: float,
182        scaled_activation_energy: float,
183        beta_x: float,
184        beta_y: float,
185    ) -> Tuple[float, float, float, float]:
186        values = np.asarray(
187            [kappa, scaled_activation_energy, beta_x, beta_y], dtype=float
188        )
189        if not np.all(np.isfinite(values)):
190            raise ValueError("model parameters must be finite")
191        if kappa <= 0.0:
192            raise ValueError("kappa must be positive")
193        if scaled_activation_energy <= 0.0:
194            raise ValueError("scaled_activation_energy must be positive")
195        if beta_x < 0.0:
196            raise ValueError(
197                "beta_x must be nonnegative because the inlet and x-upwind stencil "
198                "assume left-to-right flow"
199            )
200        return tuple(float(value) for value in values)
201
202    @staticmethod
203    def _laplacian(state: np.ndarray, dx: float, dy: float) -> np.ndarray:
204        result = np.zeros_like(state)
205        result[:, 1:-1, :] += (
206            state[:, 2:, :] - 2.0 * state[:, 1:-1, :] + state[:, :-2, :]
207        ) / dx**2
208        result[:, :, 1:-1] += (
209            state[:, :, 2:] - 2.0 * state[:, :, 1:-1] + state[:, :, :-2]
210        ) / dy**2
211        return result
212
213    @staticmethod
214    def _x_upwind_derivative(state: np.ndarray, dx: float) -> np.ndarray:
215        result = np.zeros_like(state)
216        result[:, 1:, :] = (state[:, 1:, :] - state[:, :-1, :]) / dx
217        result[:, 2:, :] = (
218            1.5 * state[:, 2:, :]
219            - 2.0 * state[:, 1:-1, :]
220            + 0.5 * state[:, :-2, :]
221        ) / dx
222        return result
223
224    @staticmethod
225    def _y_upwind_derivative(
226        state: np.ndarray, dy: float, beta_y: float
227    ) -> np.ndarray:
228        result = np.zeros_like(state)
229        if beta_y >= 0.0:
230            result[:, :, 1:] = (state[:, :, 1:] - state[:, :, :-1]) / dy
231        else:
232            result[:, :, :-1] = (state[:, :, 1:] - state[:, :, :-1]) / dy
233        return result
234
235    def _linear_rhs(
236        self, state: np.ndarray, kappa: float, beta_x: float, beta_y: float
237    ) -> np.ndarray:
238        return (
239            kappa * self._laplacian(state, self.dx, self.dy)
240            - beta_x * self._x_upwind_derivative(state, self.dx)
241            - beta_y * self._y_upwind_derivative(state, self.dy, beta_y)
242        )
243
244    def _reaction_rhs(
245        self, state: np.ndarray, scaled_activation_energy: float
246    ) -> np.ndarray:
247        theta = state[3]
248        if np.any(theta <= 0.0):
249            raise FloatingPointError("temperature became nonpositive during the solve")
250
251        rho = self.density
252        weights = self.molecular_weights
253        activation_energy = scaled_activation_energy * 1.0e3
254        exponential = np.exp(
255            -activation_energy / (self.gas_constant_times_temperature_scale * theta)
256        )
257        h2_concentration = rho * state[0] / weights[0]
258        o2_concentration = rho * state[1] / weights[1]
259        rate = (
260            h2_concentration**2
261            * o2_concentration
262            * self.reaction_prefactor
263            * exponential
264        )
265
266        reaction = np.zeros_like(state)
267        for species in range(3):
268            reaction[species] = (
269                -self.stoichiometric_coefficients[species]
270                * (weights[species] / rho)
271                * rate
272            )
273        reaction[3] = reaction[2] * self.heat_of_reaction_scaled
274        return reaction
275
276    def rhs(
277        self,
278        state: np.ndarray,
279        kappa: float,
280        scaled_activation_energy: float,
281        beta_x: float,
282        beta_y: float,
283    ) -> np.ndarray:
284        """Evaluate the semi-discrete right-hand side."""
285        kappa, scaled_activation_energy, beta_x, beta_y = self._validate_parameters(
286            kappa, scaled_activation_energy, beta_x, beta_y
287        )
288        return self._linear_rhs(state, kappa, beta_x, beta_y) + self._reaction_rhs(
289            state, scaled_activation_energy
290        )
291
292    def _apply_boundary_conditions_to_residual(
293        self, state: np.ndarray, residual: np.ndarray
294    ) -> np.ndarray:
295        state_view = state.reshape((4, self.nx, self.ny))
296        residual_view = residual.reshape((4, self.nx, self.ny))
297
298        residual_view[:, :, 0] = state_view[:, :, 1] - state_view[:, :, 0]
299        residual_view[:, :, -1] = state_view[:, :, -1] - state_view[:, :, -2]
300        residual_view[:, -1, :] = state_view[:, -1, :] - state_view[:, -2, :]
301
302        residual_view[0, 0, self._gamma_1] = state_view[0, 0, self._gamma_1]
303        residual_view[1, 0, self._gamma_1] = state_view[1, 0, self._gamma_1]
304        residual_view[2, 0, self._gamma_1] = state_view[2, 0, self._gamma_1]
305        residual_view[3, 0, self._gamma_1] = (
306            state_view[3, 0, self._gamma_1]
307            - self.ambient_temperature / self.temperature_scale
308        )
309
310        residual_view[0, 0, self._gamma_2] = (
311            state_view[0, 0, self._gamma_2] - self.inlet_h2_mass_fraction
312        )
313        residual_view[1, 0, self._gamma_2] = (
314            state_view[1, 0, self._gamma_2] - self.inlet_o2_mass_fraction
315        )
316        residual_view[2, 0, self._gamma_2] = (
317            state_view[2, 0, self._gamma_2] - self.inlet_h2o_mass_fraction
318        )
319        residual_view[3, 0, self._gamma_2] = (
320            state_view[3, 0, self._gamma_2]
321            - self.inlet_temperature / self.temperature_scale
322        )
323
324        residual_view[0, 0, self._gamma_3] = state_view[0, 0, self._gamma_3]
325        residual_view[1, 0, self._gamma_3] = state_view[1, 0, self._gamma_3]
326        residual_view[2, 0, self._gamma_3] = state_view[2, 0, self._gamma_3]
327        residual_view[3, 0, self._gamma_3] = (
328            state_view[3, 0, self._gamma_3]
329            - self.ambient_temperature / self.temperature_scale
330        )
331        return residual_view.reshape(-1)
332
333    def _residual(
334        self,
335        next_state_flat: np.ndarray,
336        previous_state: np.ndarray,
337        previous_rhs: np.ndarray,
338        kappa: float,
339        scaled_activation_energy: float,
340        beta_x: float,
341        beta_y: float,
342    ) -> np.ndarray:
343        next_state = next_state_flat.reshape((4, self.nx, self.ny))
344        next_rhs = self._linear_rhs(
345            next_state, kappa, beta_x, beta_y
346        ) + self._reaction_rhs(next_state, scaled_activation_energy)
347        residual = (
348            next_state.reshape(-1)
349            - previous_state.reshape(-1)
350            - 0.5
351            * self.dt
352            * (next_rhs.reshape(-1) + previous_rhs.reshape(-1))
353        )
354        return self._apply_boundary_conditions_to_residual(next_state_flat, residual)
355
356    def _build_linear_residual_jacobian(
357        self, kappa: float, beta_x: float, beta_y: float
358    ) -> scipy.sparse.csr_matrix:
359        block_size = self.nx * self.ny
360        block = scipy.sparse.lil_matrix((block_size, block_size), dtype=float)
361        dx2_inv = 1.0 / self.dx**2
362        dy2_inv = 1.0 / self.dy**2
363        dx_inv = 1.0 / self.dx
364        dy_inv = 1.0 / self.dy
365
366        def index(i: int, j: int) -> int:
367            return i * self.ny + j
368
369        for i in range(self.nx):
370            for j in range(self.ny):
371                row = index(i, j)
372                if i == 0:
373                    block[row, row] = 1.0
374                    continue
375                if i == self.nx - 1:
376                    block[row, row] = 1.0
377                    block[row, index(i - 1, j)] = -1.0
378                    continue
379                if j == 0:
380                    block[row, row] = -1.0
381                    block[row, index(i, j + 1)] = 1.0
382                    continue
383                if j == self.ny - 1:
384                    block[row, row] = 1.0
385                    block[row, index(i, j - 1)] = -1.0
386                    continue
387
388                diagonal = 1.0 + self.dt * kappa * (dx2_inv + dy2_inv)
389                x_im1 = -0.5 * self.dt * kappa * dx2_inv
390                x_im2 = 0.0
391                if i == 1:
392                    diagonal += 0.5 * self.dt * beta_x * dx_inv
393                    x_im1 -= 0.5 * self.dt * beta_x * dx_inv
394                else:
395                    diagonal += 0.75 * self.dt * beta_x * dx_inv
396                    x_im1 -= self.dt * beta_x * dx_inv
397                    x_im2 = 0.25 * self.dt * beta_x * dx_inv
398
399                if beta_y >= 0.0:
400                    diagonal += 0.5 * self.dt * beta_y * dy_inv
401                    y_jm1 = (
402                        -0.5 * self.dt * kappa * dy2_inv
403                        - 0.5 * self.dt * beta_y * dy_inv
404                    )
405                    y_jp1 = -0.5 * self.dt * kappa * dy2_inv
406                else:
407                    diagonal -= 0.5 * self.dt * beta_y * dy_inv
408                    y_jm1 = -0.5 * self.dt * kappa * dy2_inv
409                    y_jp1 = (
410                        -0.5 * self.dt * kappa * dy2_inv
411                        + 0.5 * self.dt * beta_y * dy_inv
412                    )
413
414                block[row, row] = diagonal
415                block[row, index(i + 1, j)] = -0.5 * self.dt * kappa * dx2_inv
416                block[row, index(i - 1, j)] = x_im1
417                block[row, index(i, j - 1)] = y_jm1
418                block[row, index(i, j + 1)] = y_jp1
419                if i >= 2:
420                    block[row, index(i - 2, j)] = x_im2
421
422        return scipy.sparse.block_diag(
423            [block.tocsr() for _ in range(4)], format="csr"
424        )
425
426    def _reaction_jacobian(
427        self, state_flat: np.ndarray, scaled_activation_energy: float
428    ) -> scipy.sparse.csr_matrix:
429        state = state_flat.reshape((4, self.nx, self.ny))
430        theta = state[3]
431        if np.any(theta <= 0.0):
432            raise FloatingPointError("temperature became nonpositive during the solve")
433
434        rho = self.density
435        weights = self.molecular_weights
436        activation_energy = scaled_activation_energy * 1.0e3
437        exponential = np.exp(
438            -activation_energy / (self.gas_constant_times_temperature_scale * theta)
439        )
440        h2_concentration = rho * state[0] / weights[0]
441        o2_concentration = rho * state[1] / weights[1]
442
443        prefactor = self.reaction_prefactor * exponential
444        rate = h2_concentration**2 * o2_concentration * prefactor
445        d_rate_d_h2 = (
446            2.0
447            * (rho / weights[0])
448            * h2_concentration
449            * o2_concentration
450            * prefactor
451        )
452        d_rate_d_o2 = (rho / weights[1]) * h2_concentration**2 * prefactor
453        d_rate_d_theta = rate * activation_energy / (
454            self.gas_constant_times_temperature_scale * theta**2
455        )
456
457        d_rate_d_h2 *= self._reaction_mask
458        d_rate_d_o2 *= self._reaction_mask
459        d_rate_d_theta *= self._reaction_mask
460
461        block_size = self.nx * self.ny
462        local = np.arange(block_size, dtype=int)
463        rows = []
464        cols = []
465        data = []
466
467        for species in range(4):
468            if species < 3:
469                coefficient = (
470                    -self.stoichiometric_coefficients[species]
471                    * weights[species]
472                    / rho
473                )
474            else:
475                coefficient = (
476                    -self.stoichiometric_coefficients[2]
477                    * weights[2]
478                    / rho
479                    * self.heat_of_reaction_scaled
480                )
481
482            row = species * block_size + local
483            rows.extend([row, row, row])
484            cols.extend(
485                [
486                    0 * block_size + local,
487                    1 * block_size + local,
488                    3 * block_size + local,
489                ]
490            )
491            data.extend(
492                [
493                    coefficient * d_rate_d_h2.reshape(-1),
494                    coefficient * d_rate_d_o2.reshape(-1),
495                    coefficient * d_rate_d_theta.reshape(-1),
496                ]
497            )
498
499        return scipy.sparse.csr_matrix(
500            (np.concatenate(data), (np.concatenate(rows), np.concatenate(cols))),
501            shape=(4 * block_size, 4 * block_size),
502        )
503
504    def _newton_solve(
505        self,
506        previous_state: np.ndarray,
507        previous_rhs: np.ndarray,
508        linear_jacobian: scipy.sparse.csr_matrix,
509        kappa: float,
510        scaled_activation_energy: float,
511        beta_x: float,
512        beta_y: float,
513        linear_preconditioner=None,
514    ) -> np.ndarray:
515        iterate = previous_state.reshape(-1).copy()
516        residual = self._residual(
517            iterate,
518            previous_state,
519            previous_rhs,
520            kappa,
521            scaled_activation_energy,
522            beta_x,
523            beta_y,
524        )
525        initial_norm = np.linalg.norm(residual)
526        target = (
527            self.newton_absolute_tolerance
528            + self.newton_relative_tolerance * initial_norm
529        )
530        if initial_norm <= target:
531            return iterate.reshape(previous_state.shape)
532
533        for iteration in range(1, self.newton_max_iterations + 1):
534            reaction_jacobian = self._reaction_jacobian(
535                iterate, scaled_activation_energy
536            )
537            # Crank-Nicolson contributes one half of the new-state RHS Jacobian.
538            jacobian = linear_jacobian - 0.5 * self.dt * reaction_jacobian
539            update = self._solve_newton_linear_system(
540                jacobian, -residual, linear_preconditioner
541            )
542            if not np.all(np.isfinite(update)):
543                raise RuntimeError(
544                    "Newton linear solve produced non-finite values at iteration {}".format(
545                        iteration
546                    )
547                )
548            iterate += update
549            residual = self._residual(
550                iterate,
551                previous_state,
552                previous_rhs,
553                kappa,
554                scaled_activation_energy,
555                beta_x,
556                beta_y,
557            )
558            residual_norm = np.linalg.norm(residual)
559            if not np.isfinite(residual_norm):
560                raise RuntimeError(
561                    "Newton residual became non-finite at iteration {}".format(iteration)
562                )
563            if residual_norm <= target:
564                return iterate.reshape(previous_state.shape)
565
566        final_relative = np.linalg.norm(residual) / max(
567            initial_norm, np.finfo(float).tiny
568        )
569        raise RuntimeError(
570            "Newton solver failed to converge in {} iterations; "
571            "final relative residual={:.3e}".format(
572                self.newton_max_iterations, final_relative
573            )
574        )
575
576    def _solve_newton_linear_system(
577        self,
578        jacobian: scipy.sparse.csr_matrix,
579        right_hand_side: np.ndarray,
580        linear_preconditioner,
581    ) -> np.ndarray:
582        """Solve one exact Newton system, falling back to direct LU if needed."""
583        if linear_preconditioner is not None:
584            preconditioner = scipy.sparse.linalg.LinearOperator(
585                jacobian.shape,
586                matvec=linear_preconditioner.solve,
587                dtype=float,
588            )
589            callback_iterations = [0]
590            update, info = scipy.sparse.linalg.gmres(
591                jacobian,
592                right_hand_side,
593                M=preconditioner,
594                rtol=self.gmres_relative_tolerance,
595                atol=self.gmres_absolute_tolerance,
596                restart=self.gmres_restart,
597                maxiter=self.gmres_max_iterations,
598                callback=lambda _residual: callback_iterations.__setitem__(
599                    0, callback_iterations[0] + 1
600                ),
601                callback_type="pr_norm",
602            )
603            self.linear_solver_stats["gmres_iterations"] += callback_iterations[0]
604            linear_residual = np.linalg.norm(jacobian @ update - right_hand_side)
605            linear_target = max(
606                self.gmres_absolute_tolerance,
607                self.gmres_relative_tolerance * np.linalg.norm(right_hand_side),
608            )
609            if (
610                info == 0
611                and np.all(np.isfinite(update))
612                and linear_residual <= linear_target
613            ):
614                return update
615            self.linear_solver_stats["gmres_direct_fallbacks"] += 1
616
617        self.linear_solver_stats["direct_solves"] += 1
618        return scipy.sparse.linalg.spsolve(
619            jacobian,
620            right_hand_side,
621            permc_spec="MMD_AT_PLUS_A",
622            use_umfpack=False,
623        )
624
625    def solve(
626        self,
627        kappa: float,
628        scaled_activation_energy: float,
629        beta_x: float,
630        beta_y: float,
631    ) -> Tuple[np.ndarray, np.ndarray]:
632        """Integrate the flame model and return ``(states, times)``.
633
634        ``states`` has shape ``(num_snapshots, 4, nx, ny)``. The temperature
635        component is nondimensional; multiply ``states[:, 3]`` by
636        :attr:`temperature_scale` to obtain kelvin.
637        """
638        kappa, scaled_activation_energy, beta_x, beta_y = self._validate_parameters(
639            kappa, scaled_activation_energy, beta_x, beta_y
640        )
641        linear_jacobian = self._build_linear_residual_jacobian(
642            kappa, beta_x, beta_y
643        )
644        self.linear_solver_stats = {
645            "gmres_iterations": 0,
646            "gmres_direct_fallbacks": 0,
647            "direct_solves": 0,
648        }
649        linear_preconditioner = None
650        if self.linear_solver == "gmres":
651            linear_preconditioner = scipy.sparse.linalg.splu(
652                linear_jacobian.tocsc(), permc_spec="MMD_AT_PLUS_A"
653            )
654
655        state = self.initial_state()
656        states = [state.copy()]
657        times = [0.0]
658
659        for step in range(1, self.num_steps + 1):
660            previous_state = state.copy()
661            previous_rhs = self._linear_rhs(
662                previous_state, kappa, beta_x, beta_y
663            ) + self._reaction_rhs(previous_state, scaled_activation_energy)
664            state = self._newton_solve(
665                previous_state,
666                previous_rhs,
667                linear_jacobian,
668                kappa,
669                scaled_activation_energy,
670                beta_x,
671                beta_y,
672                linear_preconditioner,
673            )
674            if step % self.snapshot_stride == 0 or step == self.num_steps:
675                states.append(state.copy())
676                times.append(step * self.dt)
677
678        return np.stack(states, axis=0), np.asarray(times)
679
680
681def extract_temperature_sensors(
682    states: np.ndarray,
683    spatial_stride: int = 4,
684    temporal_stride: int = 1,
685    temperature_scale: float = 1.0,
686) -> np.ndarray:
687    """Extract a flattened temperature-sensor QoI from a state history.
688
689    Parameters
690    ----------
691    states:
692        State history with shape ``(num_snapshots, 4, nx, ny)``.
693    spatial_stride:
694        Sample every ``spatial_stride`` grid points in each spatial direction.
695    temporal_stride:
696        Sample every ``temporal_stride`` stored snapshots.
697    temperature_scale:
698        Optional multiplicative scale. Use ``300.0`` to convert the default
699        nondimensional temperature to kelvin.
700    """
701    states = np.asarray(states)
702    if states.ndim != 4 or states.shape[1] != 4:
703        raise ValueError("states must have shape (num_snapshots, 4, nx, ny)")
704    if spatial_stride < 1 or temporal_stride < 1:
705        raise ValueError("sensor strides must be at least 1")
706    if not np.isfinite(temperature_scale) or temperature_scale <= 0.0:
707        raise ValueError("temperature_scale must be finite and positive")
708
709    temperature = states[
710        ::temporal_stride, 3, ::spatial_stride, ::spatial_stride
711    ]
712    return (temperature_scale * temperature).reshape(-1)

romtools QoI wrapper:

 1"""romtools QoI-model wrapper for the pure-Python H2-air flame benchmark."""
 2
 3import os
 4from typing import Dict
 5
 6import numpy as np
 7
 8from h2_air_flame import H2AirFlame, extract_temperature_sensors
 9
10
11class H2AirFlameQoiModel:
12    """Sensor-based ``QoiModel`` wrapper around :class:`H2AirFlame`.
13
14    The QoI consists of nondimensional temperature samples on regular spatial
15    and temporal strides. Full state histories are saved to each run directory
16    so the same FOM evaluations can also be reused for ROM construction.
17    """
18
19    parameter_names = H2AirFlame.parameter_names
20
21    def __init__(
22        self,
23        nx: int = 25,
24        ny: int = 13,
25        dt: float = 1.0e-4,
26        t_end: float = 5.0e-3,
27        snapshot_stride: int = 10,
28        spatial_sensor_stride: int = 4,
29        temporal_sensor_stride: int = 1,
30    ) -> None:
31        self._model_kwargs = {
32            "nx": nx,
33            "ny": ny,
34            "dt": dt,
35            "t_end": t_end,
36            "snapshot_stride": snapshot_stride,
37        }
38        self.spatial_sensor_stride = int(spatial_sensor_stride)
39        self.temporal_sensor_stride = int(temporal_sensor_stride)
40        if self.spatial_sensor_stride < 1 or self.temporal_sensor_stride < 1:
41            raise ValueError("sensor strides must be at least 1")
42
43    def populate_run_directory(self, run_directory: str, parameter_sample: Dict) -> None:
44        """Create the run directory required by the workflow model protocol."""
45        del parameter_sample
46        os.makedirs(run_directory, exist_ok=True)
47
48    def _parameters(self, parameter_sample: Dict):
49        missing = [name for name in self.parameter_names if name not in parameter_sample]
50        if missing:
51            raise KeyError(
52                "parameter_sample is missing required keys: {}".format(missing)
53            )
54        return [float(parameter_sample[name]) for name in self.parameter_names]
55
56    def run_model(self, run_directory: str, parameter_sample: Dict) -> int:
57        """Run the flame model and save states, times, parameters, and QoI."""
58        model = H2AirFlame(**self._model_kwargs)
59        parameters = self._parameters(parameter_sample)
60        states, times = model.solve(*parameters)
61        qoi = extract_temperature_sensors(
62            states,
63            spatial_stride=self.spatial_sensor_stride,
64            temporal_stride=self.temporal_sensor_stride,
65        )
66        np.savez(
67            os.path.join(run_directory, "solution.npz"),
68            states=states,
69            times=times,
70            qoi=qoi,
71            parameter_names=np.asarray(self.parameter_names),
72            parameters=np.asarray(parameters),
73        )
74        return 0
75
76    def compute_qoi(self, run_directory: str, parameter_sample: Dict) -> np.ndarray:
77        """Read and return the flattened temperature-sensor QoI."""
78        del parameter_sample
79        with np.load(os.path.join(run_directory, "solution.npz")) as data:
80            return np.asarray(data["qoi"], dtype=float)
81
82
83if __name__ == "__main__":
84    import tempfile
85
86    sample = {
87        "kappa": 2.0,
88        "scaled_activation_energy": 8.0,
89        "beta_x": 40.0,
90        "beta_y": 7.0,
91    }
92    with tempfile.TemporaryDirectory() as directory:
93        model = H2AirFlameQoiModel(
94            nx=12, ny=7, t_end=2.0e-4, snapshot_stride=1
95        )
96        model.populate_run_directory(directory, sample)
97        model.run_model(directory, sample)
98        qoi = model.compute_qoi(directory, sample)
99        print("H2-air flame QoI size:", qoi.size)