Iterative Source Reweighting

This one-group fixed-source example demonstrates partial model updates across several runs of one mcdc.Simulation. A homogeneous slab contains symmetric left and right sources. Each iteration changes only their relative probabilities, compiles a complete new snapshot, and writes a separate output file.

The three source mixtures are 20/80, 50/50, and 80/20. Because the geometry and material are symmetric, the 20/80 and 80/20 flux profiles should approximately mirror one another, while the 50/50 profile should be approximately symmetric about the slab midpoint.

Full Input

import numpy as np

import mcdc

simulation = mcdc.Simulation("Iterative source reweighting")

# Homogeneous one-group slab
material = mcdc.MaterialMG(
    capture=np.array([0.2]),
    scatter=np.array([[0.8]]),
)
left_boundary = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum")
right_boundary = mcdc.Surface.PlaneZ(z=10.0, boundary_condition="vacuum")
slab = mcdc.Cell(
    region=+left_boundary & -right_boundary,
    fill=material,
)
simulation.set_model([slab])

# Symmetric source regions
source_left = mcdc.Source(
    name="Left source",
    z=[1.0, 2.0],
    isotropic=True,
    energy_group=0,
)
source_right = mcdc.Source(
    name="Right source",
    z=[8.0, 9.0],
    isotropic=True,
    energy_group=0,
)
simulation.set_sources([source_left, source_right])

# One tally is reused by every iteration
mesh = mcdc.MeshStructured(z=np.linspace(0.0, 10.0, 51))
flux_tally = mcdc.Tally(
    name="source_mix_flux",
    mesh=mesh,
    scores=["flux"],
)
simulation.set_tallies([flux_tally])

simulation.settings.N_particle = 20_000
simulation.settings.N_batch = 10

# Partially update the same owned model and compile a fresh snapshot each time.
source_mixes = (
    ("left_20", 0.2),
    ("left_50", 0.5),
    ("left_80", 0.8),
)

for case_name, left_fraction in source_mixes:
    source_left.probability = left_fraction
    source_right.probability = 1.0 - left_fraction
    simulation.settings.output_name = f"source_mix_{case_name}"

    simulation.compile()
    simulation.run()

Post-processing

from pathlib import Path

import h5py
import matplotlib.pyplot as plt
import numpy as np

cases = (
    ("left_20", "20% left / 80% right"),
    ("left_50", "50% left / 50% right"),
    ("left_80", "80% left / 20% right"),
)

profiles = {}
uncertainties = {}
z = None

for case_name, _ in cases:
    output_path = Path(f"source_mix_{case_name}.h5")
    with h5py.File(output_path, "r") as output:
        tally = output["tallies/source_mix_flux"]
        case_z = tally["grid/z"][:]
        flux = tally["flux/mean"][:]
        flux_sdev = tally["flux/sdev"][:]

    if z is None:
        z = case_z
    elif not np.array_equal(z, case_z):
        raise ValueError(f"Inconsistent tally grid in {output_path}")

    dz = case_z[1:] - case_z[:-1]
    profiles[case_name] = flux / dz
    uncertainties[case_name] = flux_sdev / dz

z_mid = 0.5 * (z[:-1] + z[1:])
dz = z[1:] - z[:-1]
left_half = z_mid < 5.0
right_half = ~left_half

print("Integrated flux comparison")
print("--------------------------")
for case_name, label in cases:
    profile = profiles[case_name]
    left_flux = np.sum(profile[left_half] * dz[left_half])
    right_flux = np.sum(profile[right_half] * dz[right_half])
    print(
        f"{label:24s}  left={left_flux:.6e}  right={right_flux:.6e}  "
        f"left/right={left_flux / right_flux:.4f}"
    )

mirror_difference = np.linalg.norm(profiles["left_20"] - profiles["left_80"][::-1])
mirror_scale = np.linalg.norm(0.5 * (profiles["left_20"] + profiles["left_80"][::-1]))
balanced_difference = np.linalg.norm(profiles["left_50"] - profiles["left_50"][::-1])
balanced_scale = np.linalg.norm(profiles["left_50"])

print()
print(
    "20/80 versus mirrored 80/20 relative RMS difference: "
    f"{mirror_difference / mirror_scale:.4e}"
)
print(
    "50/50 profile relative left-right asymmetry: "
    f"{balanced_difference / balanced_scale:.4e}"
)

figure, axis = plt.subplots()
for case_name, label in cases:
    profile = profiles[case_name]
    uncertainty = uncertainties[case_name]
    axis.plot(z_mid, profile, label=label)
    axis.fill_between(
        z_mid,
        profile - uncertainty,
        profile + uncertainty,
        alpha=0.15,
    )

axis.axvline(5.0, color="black", linestyle="--", linewidth=1.0)
axis.set_xlabel("z [cm]")
axis.set_ylabel("Flux")
axis.grid()
axis.legend()
figure.tight_layout()
figure.savefig("iterative_source_comparison.png", dpi=150)
plt.close(figure)

The post-processing script overlays the three spatial flux profiles, prints the integrated flux in each half of the slab, and reports two symmetry comparisons.

How to Run

From inside examples/iterative_source_reweighting:

python input.py
python process-output.py

The calculation writes source_mix_left_20.h5, source_mix_left_50.h5, and source_mix_left_80.h5. Post-processing writes iterative_source_comparison.png.