Steady convection-diffusion-reaction model#

The examples include a self-contained steady convection-diffusion-reaction (CDR) model that is used throughout the romtools workflow demonstrations. It is intentionally small and inexpensive, making it useful for testing romtools interfaces and workflows without an external simulation code.

The implementation lives in examples/models and requires only NumPy and SciPy.

Problem definition#

The model solves a scalar steady CDR equation on the unit square,

\[\nu \nabla^2 u - \boldsymbol{b}\cdot\nabla u - \sigma u = -1,\]

with homogeneous boundary data represented by the discrete operators. The four model parameters are

bmag

Magnitude of the advection velocity.

theta

Advection direction in radians. The velocity vector is

\[\boldsymbol{b} = b_{\mathrm{mag}} [\cos(\theta),\,\sin(\theta)]^T.\]
nu

Diffusion coefficient.

sigma

Reaction coefficient.

The spatial discretization uses a second-order central stencil for diffusion and upwind finite differences for advection. The resulting sparse linear system is assembled with SciPy and solved with scipy.sparse.linalg.spsolve.

Direct use#

The PDE solver can be used independently of the romtools workflow interfaces. From examples/models:

import numpy as np
import steady_cdr as cdr

system = cdr.AdvectionDiffusionSystem(nx=21, ny=21)
bmag = 0.5
theta = np.pi / 3.0
b = bmag * np.array([np.cos(theta), np.sin(theta)])

state = cdr.solveFom(
    system,
    b=b,
    nu=1.0e-3,
    sigma=1.0,
)

state contains the solution values at the nx * ny interior grid points.

romtools QoI model#

SteadyCdrQoiModel wraps the solver using the standard romtools QoI-model protocol. It implements populate_run_directory, run_model, and compute_qoi so the model can be passed directly to romtools workflows. For each evaluation, the wrapper stores the state, model parameters, and a one-sided boundary-derivative QoI in solution.npz.

A representative parameter sample is

sample = {
    "bmag": 0.5,
    "theta": np.pi / 3.0,
    "nu": 1.0e-3,
    "sigma": 1.0,
}

The complete wrapper can be exercised from the repository root with

python examples/models/steady_cdr_model.py

The script performs one model evaluation and reports the size and norm of the resulting QoI.

Implementation#

Core CDR operators and solver:

  1"""Steady convection-diffusion-reaction (CDR) operators and solve utilities."""
  2
  3import numpy as np
  4import scipy.sparse
  5import scipy.sparse.linalg
  6
  7
  8class _SparseMatrixBuilder:
  9    """Incremental sparse matrix builder."""
 10
 11    def __init__(self, n_rows: int, n_cols: int):
 12        self._n_rows = n_rows
 13        self._n_cols = n_cols
 14        self._row_indices = []
 15        self._col_indices = []
 16        self._values = []
 17
 18    def add_entry(self, row_index: int, col_index: int, value: float) -> None:
 19        self._row_indices.append(row_index)
 20        self._col_indices.append(col_index)
 21        self._values.append(value)
 22
 23    def assemble(self) -> scipy.sparse.csr_matrix:
 24        return scipy.sparse.csr_matrix(
 25            (self._values, (self._row_indices, self._col_indices)),
 26            shape=(self._n_rows, self._n_cols),
 27        )
 28
 29
 30def build_advection_matrices(nx: int, ny: int, dx: float, dy: float):
 31    """Build first-order advection operators using 2nd-order upwind stencils."""
 32    ax_builder = _SparseMatrixBuilder(nx * ny, nx * ny)
 33    ay_builder = _SparseMatrixBuilder(nx * ny, nx * ny)
 34
 35    for j in range(1, ny):
 36        for i in range(1, nx):
 37            index = i + j * nx
 38            index_im1 = index - 1
 39            index_im2 = index - 2
 40            index_jm1 = index - nx
 41            index_jm2 = index - 2 * nx
 42
 43            ax_builder.add_entry(index, index, 1.5 / dx)
 44            ax_builder.add_entry(index, index_im1, -2.0 / dx)
 45            if i > 1:
 46                ax_builder.add_entry(index, index_im2, 0.5 / dx)
 47
 48            ay_builder.add_entry(index, index, 1.5 / dy)
 49            ay_builder.add_entry(index, index_jm1, -2.0 / dy)
 50            if j > 1:
 51                ay_builder.add_entry(index, index_jm2, 0.5 / dy)
 52
 53    for i in range(1, nx):
 54        index = i
 55        index_im1 = index - 1
 56        index_im2 = index - 2
 57
 58        ax_builder.add_entry(index, index, 1.5 / dx)
 59        ax_builder.add_entry(index, index_im1, -2.0 / dx)
 60        if i > 1:
 61            ax_builder.add_entry(index, index_im2, 0.5 / dx)
 62        ay_builder.add_entry(index, index, 1.0 / dy)
 63
 64    for j in range(1, ny):
 65        index = j * nx
 66        index_jm1 = index - nx
 67        index_jm2 = index - 2 * nx
 68
 69        ax_builder.add_entry(index, index, 1.0 / dx)
 70        ay_builder.add_entry(index, index, 1.5 / dy)
 71        ay_builder.add_entry(index, index_jm1, -2.0 / dy)
 72        if j > 1:
 73            ay_builder.add_entry(index, index_jm2, 0.5 / dy)
 74
 75    ax_builder.add_entry(0, 0, 1.0 / dx)
 76    ay_builder.add_entry(0, 0, 1.0 / dy)
 77
 78    return ax_builder.assemble(), ay_builder.assemble()
 79
 80
 81def build_diffusion_matrix(nx: int, ny: int, dx: float, dy: float):
 82    """Build Laplacian operator using 2nd-order central stencils."""
 83    matrix_builder = _SparseMatrixBuilder(nx * ny, nx * ny)
 84    for j in range(ny):
 85        for i in range(nx):
 86            index = i + j * nx
 87            index_im1 = index - 1
 88            index_ip1 = index + 1
 89            index_jm1 = index - nx
 90            index_jp1 = index + nx
 91
 92            matrix_builder.add_entry(index, index, -2.0 / dx**2 - 2.0 / dy**2)
 93            if i > 0:
 94                matrix_builder.add_entry(index, index_im1, 1.0 / dx**2)
 95            if i < nx - 1:
 96                matrix_builder.add_entry(index, index_ip1, 1.0 / dx**2)
 97            if j > 0:
 98                matrix_builder.add_entry(index, index_jm1, 1.0 / dy**2)
 99            if j < ny - 1:
100                matrix_builder.add_entry(index, index_jp1, 1.0 / dy**2)
101    return matrix_builder.assemble()
102
103
104def build_qoi_vector(nx: int, ny: int, dx: float, dy: float):
105    """Build QoI vector for integral of du/dx along the right boundary."""
106    qoi_vector = np.zeros(nx * ny)
107    for j in range(ny):
108        index = nx - 1 + j * nx
109        index_im1 = index - 1
110        qoi_vector[index] += -2.0 / dx * dy
111        qoi_vector[index_im1] += 0.5 / dx * dy
112    return qoi_vector
113
114
115class AdvectionDiffusionSystem:
116    """Container for steady CDR operators on a uniform grid."""
117
118    def __init__(self, nx: int, ny: int):
119        self.Lx = 1.0
120        self.Ly = 1.0
121        self.Nx = nx
122        self.Ny = ny
123        self.N = nx * ny
124        self.dx = self.Lx / (self.Nx + 1)
125        self.dy = self.Ly / (self.Ny + 1)
126        self.g = np.ones(self.N)
127        self.x = np.linspace(self.dx, self.Lx - self.dx, self.Nx)
128        self.y = np.linspace(self.dy, self.Ly - self.dy, self.Ny)
129
130        self.A_diffusion = build_diffusion_matrix(self.Nx, self.Ny, self.dx, self.dy)
131        self.A_advection_x, self.A_advection_y = build_advection_matrices(
132            self.Nx, self.Ny, self.dx, self.dy
133        )
134        self.I = scipy.sparse.identity(self.N, format="csr")
135        self.C = build_qoi_vector(self.Nx, self.Ny, self.dx, self.dy)
136
137
138def solveFom(system: AdvectionDiffusionSystem, b: np.ndarray, nu: float, sigma: float):
139    """Solve the steady CDR system for a given parameter sample."""
140    lhs = (
141        system.A_diffusion * nu
142        - b[0] * system.A_advection_x
143        - b[1] * system.A_advection_y
144        - sigma * system.I
145    )
146    rhs = -system.g
147    return scipy.sparse.linalg.spsolve(lhs, rhs)
148

romtools QoI wrapper:

 1"""ROM-tools model wrapper for the steady convection-diffusion-reaction equation."""
 2
 3import os
 4import numpy as np
 5
 6try:
 7    import steady_cdr as cdr
 8except ImportError:  # pragma: no cover
 9    from . import steady_cdr as cdr
10
11
12class SteadyCdrQoiModel:
13    """QoiModel-compatible wrapper for the steady CDR equation."""
14
15    def __init__(self, nx: int = 31, ny: int = 31) -> None:
16        self.system = cdr.AdvectionDiffusionSystem(nx, ny)
17
18    def populate_run_directory(self, run_directory: str, parameter_sample: dict) -> None:
19        with open(os.path.join(run_directory, "params.txt"), "w", encoding="utf-8") as handle:
20            for key, value in parameter_sample.items():
21                handle.write(f"{key}: {value}\n")
22
23    def run_model(self, run_directory: str, parameter_sample: dict) -> int:
24        b_magnitude = float(parameter_sample["bmag"])
25        theta = float(parameter_sample["theta"])
26        nu = float(parameter_sample["nu"])
27        sigma = float(parameter_sample["sigma"])
28
29        b_vec = np.zeros(2)
30        b_vec[0] = b_magnitude * np.cos(theta)
31        b_vec[1] = b_magnitude * np.sin(theta)
32
33        state = cdr.solveFom(self.system, b_vec, nu, sigma)
34
35        state_grid = np.zeros((self.system.Nx + 2, self.system.Ny + 2))
36        state_grid[1:-1, 1:-1] = state.reshape(self.system.Nx, self.system.Ny)
37
38        # QoI: one-sided estimate of du/dx along the right boundary.
39        qoi = -state_grid[:, -2] / self.system.dx
40        params = np.array([b_magnitude, theta, nu, sigma], dtype=float)
41
42        np.savez(
43            os.path.join(run_directory, "solution.npz"),
44            state=state,
45            qoi=qoi,
46            params=params,
47        )
48        return 0
49
50    def compute_qoi(self, run_directory: str, parameter_sample: dict) -> np.ndarray:
51        solution = np.load(os.path.join(run_directory, "solution.npz"))
52        return solution["qoi"]
53
54
55CdrQoiModel = SteadyCdrQoiModel
56
57
58def _main():
59    run_directory = os.path.join(os.path.dirname(__file__), "steady_cdr_output")
60    os.makedirs(run_directory, exist_ok=True)
61
62    model = SteadyCdrQoiModel(nx=21, ny=21)
63    sample = {"bmag": 0.5, "theta": np.pi / 3.0, "nu": 1e-3, "sigma": 1.0}
64    model.populate_run_directory(run_directory, sample)
65    model.run_model(run_directory, sample)
66    qoi = model.compute_qoi(run_directory, sample)
67    print(f"Saved {os.path.join(run_directory, 'solution.npz')}")
68    print(f"QoI shape: {qoi.shape}, QoI L2 norm: {np.linalg.norm(qoi):.6e}")
69
70
71if __name__ == "__main__":
72    _main()