Monte Carlo UQ for Convection-Diffusion-Reaction#

This example compares standard Monte Carlo (MC) and two-level multifidelity Monte Carlo (MFMC) for a steady convection-diffusion-reaction problem. The scientific question is simple: for a fixed computational budget, how accurately can we estimate the expected integrated right-boundary flux?

Problem definition#

On the unit square, the full-order model solves

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

with homogeneous Dirichlet boundary conditions. The advection vector is written

\[\begin{split}\mathbf{b} = b \begin{bmatrix} \cos\theta \\ \sin\theta \end{bmatrix}.\end{split}\]

The four uncertain inputs are independent uniform random variables:

Parameter

Meaning

Distribution

\(b\)

advection magnitude

\(\mathcal{U}(0.5,1.5)\)

\(\theta\)

advection direction

\(\mathcal{U}(\pi/6,\pi/3)\)

\(\nu\)

diffusion coefficient

\(\mathcal{U}(0.02,0.08)\)

\(\sigma\)

reaction coefficient

\(\mathcal{U}(0.1,0.6)\)

The scalar quantity of interest is the integral of the one-sided estimate of \(\partial u/\partial x\) along the right boundary.

Fidelity hierarchy and pilot#

The high-fidelity model uses a 21-by-21 discretization and the low-fidelity model uses a 9-by-9 discretization. The low-fidelity cost is taken to be 5% of one high-fidelity evaluation.

MFMC begins with 30 paired pilot samples. The pilot therefore costs 31.5 high-fidelity-equivalent evaluations and is included in every reported MFMC budget. For the documented seed, the pilot Pearson correlation between the high- and low-fidelity QoIs is

\[\rho_{\mathrm{pilot}} = 0.850.\]

This is the correlation used by the allocator. The example intentionally reports only the pilot correlation, because it is the diagnostic that drives the sample-allocation decision.

High- versus low-fidelity QoIs for the MFMC pilot samples

The 30 paired pilot QoIs. Their correlation is used to determine the high- and low-fidelity sample allocation.#

Reference result#

The reference mean was generated with 20,000 high-fidelity Monte Carlo samples on the 21-by-21 grid using random seed 314159:

\[\mu_{\mathrm{ref}} = -5.18032649,\]

with an estimated reference standard error of \(1.10\times10^{-2}\). Reference generation is kept separate from normal documentation and CI builds; it can be reproduced with

python examples/uq_cdr_demo/generate_reference.py

The committed reference_stats.json records the sample count, seed, mean, variance, and standard error.

Equal-cost comparison#

Because the 30-sample pilot already costs 31.5 high-fidelity-equivalent evaluations, the benchmark compares MC and MFMC at budgets of 40, 80, 160, and 320. MFMC uses the same 30-sample pilot and 5% low/high cost ratio at every budget. The pilot evaluations are included in the reported \(N_H\) and \(N_L\) totals and in the equivalent-cost constraint.

Budget

Method

\(N_H\)

\(N_L\)

Mean

Relative error

Est. standard error

40

MC

40

0

-5.2920

2.16%

0.245

40

MFMC

30

200

-5.0328

2.85%

0.109

80

MC

80

0

-5.2267

0.89%

0.168

80

MFMC

59

420

-4.9908

3.66%

0.086

160

MC

160

0

-5.3046

2.40%

0.134

160

MFMC

118

840

-5.1292

0.99%

0.074

320

MC

320

0

-5.2669

1.67%

0.086

320

MFMC

235

1700

-5.1879

0.15%

0.054

A single realized relative-error curve is not expected to be monotone because both estimators are stochastic. The estimated standard errors provide the cleaner variance-reduction diagnostic: MFMC is lower at every reported budget.

Relative error of MC and MFMC versus equivalent high-fidelity cost

Error in the estimated mean relative to the 20,000-sample high-fidelity reference calculation.#

Run the examples#

The lightweight example used by the smoke test is

python examples/uq_cdr_demo/example.py

To regenerate the full equal-cost benchmark and its figures, run

python examples/uq_cdr_demo/benchmark.py

The benchmark writes all generated results below examples/uq_cdr_demo/benchmark_output. The expensive reference calculation is not part of normal CI or documentation generation.

Implementation#

 1"""Monte Carlo UQ for the steady convection-diffusion-reaction example."""
 2
 3import argparse
 4import os
 5from pathlib import Path
 6import sys
 7from typing import Optional
 8
 9import numpy as np
10
11from romtools.workflows.parameter_spaces import UniformParameterSpace
12from romtools.workflows.sampling_methods import MonteCarloSampler
13from romtools.workflows.uq import run_monte_carlo, run_multifidelity_monte_carlo
14
15
16EXAMPLE_DIRECTORY = Path(__file__).resolve().parent
17MODEL_DIRECTORY = EXAMPLE_DIRECTORY.parent / "models"
18if str(MODEL_DIRECTORY) not in sys.path:
19    sys.path.insert(0, str(MODEL_DIRECTORY))
20
21from steady_cdr_model import SteadyCdrQoiModel  # noqa: E402
22
23
24class IntegratedFluxCdrModel(SteadyCdrQoiModel):
25    """Return the integrated right-boundary flux as a scalar QoI."""
26
27    def compute_qoi(self, run_directory: str, parameter_sample: dict) -> np.ndarray:
28        with np.load(os.path.join(run_directory, "solution.npz")) as solution:
29            boundary_flux = solution["qoi"]
30        return np.array([np.sum(boundary_flux) * self.system.dy])
31
32
33def build_parameter_space() -> UniformParameterSpace:
34    """Construct independent uniform distributions for the CDR parameters."""
35    return UniformParameterSpace(
36        parameter_names=["bmag", "theta", "nu", "sigma"],
37        lower_bounds=[0.5, np.pi / 6.0, 0.02, 0.1],
38        upper_bounds=[1.5, np.pi / 3.0, 0.08, 0.6],
39        sampler=MonteCarloSampler,
40    )
41
42
43def main(output_directory: Optional[Path] = None) -> None:
44    """Run single- and multifidelity CDR mean-estimation examples."""
45    output_directory = (
46        EXAMPLE_DIRECTORY / "uq_cdr_output"
47        if output_directory is None
48        else Path(output_directory).resolve()
49    )
50    output_directory.mkdir(parents=True, exist_ok=True)
51    parameter_space = build_parameter_space()
52
53    mc_result = run_monte_carlo(
54        model=IntegratedFluxCdrModel(nx=21, ny=21),
55        parameter_space=parameter_space,
56        absolute_work_dir=str(output_directory / "monte_carlo"),
57        number_of_samples=40,
58        random_seed=7,
59    )
60
61    mfmc_directory = output_directory / "multifidelity"
62    mfmc_result = run_multifidelity_monte_carlo(
63        high_fidelity_model=IntegratedFluxCdrModel(nx=21, ny=21),
64        low_fidelity_model=IntegratedFluxCdrModel(nx=9, ny=9),
65        parameter_space=parameter_space,
66        absolute_work_dir=str(mfmc_directory),
67        pilot_sample_count=30,
68        high_fidelity_equivalent_budget=40.0,
69        low_to_high_fidelity_cost_ratio=0.05,
70        allocation_qoi_index=0,
71        random_seed=7,
72    )
73    with np.load(mfmc_directory / "uq_stats.npz") as statistics:
74        pilot_correlation = float(statistics["pilot_correlations"][0])
75
76    print("CDR integrated right-boundary flux")
77    print(
78        f"MC:   mean={mc_result.mean[0]:.6e}, "
79        f"standard error={mc_result.standard_error[0]:.3e}"
80    )
81    print(
82        f"MFMC: mean={mfmc_result.mean[0]:.6e}, "
83        f"standard error={mfmc_result.standard_error[0]:.3e}, "
84        f"N_H={mfmc_result.high_fidelity_sample_count}, "
85        f"N_L={mfmc_result.low_fidelity_sample_count}, "
86        f"pilot correlation={pilot_correlation:.3f}"
87    )
88    print(f"Results written to {output_directory}")
89
90
91if __name__ == "__main__":
92    parser = argparse.ArgumentParser(description=__doc__)
93    parser.add_argument(
94        "--output-directory",
95        type=Path,
96        help="Output path (default: examples/uq_cdr_demo/uq_cdr_output)",
97    )
98    arguments = parser.parse_args()
99    main(arguments.output_directory)