From 083129bbf3c360cde0ee30ffa4b0314fb6a34160 Mon Sep 17 00:00:00 2001 From: Kamil Kania Date: Mon, 20 Jul 2026 16:28:14 +0200 Subject: [PATCH 1/3] CI/docs: Execution of Python snippets from first_mcdc.rst --- .github/workflows/docs_test.yml | 23 ++++- docs/source/user/first_mcdc.rst | 2 +- extract.py | 15 +++ test.py | 166 ++++++++++++++++++++++++++++++++ 4 files changed, 204 insertions(+), 2 deletions(-) create mode 100644 extract.py create mode 100644 test.py diff --git a/.github/workflows/docs_test.yml b/.github/workflows/docs_test.yml index 24e1ca536..fc5a8476f 100644 --- a/.github/workflows/docs_test.yml +++ b/.github/workflows/docs_test.yml @@ -39,4 +39,25 @@ jobs: - name: Build the docs working-directory: docs - run: make html \ No newline at end of file + run: make html + + - name: Run Python snippets + run: | + python -c ' + from docutils.core import publish_doctree + from docutils import nodes + from docutils.parsers.rst import roles + + def dummy_role(name, rawtext, text, lineno, inliner, options=None, content=None): + return [nodes.Text(text)], [] + + roles.register_local_role("doc", dummy_role) + + rst = open("docs/source/user/first_mcdc.rst").read() + + for node in publish_doctree(rst).findall(nodes.literal_block): + if node.get("classes") == ['code', 'python3']: + print(node.astext()) + ' > tmp.py + cat -n tmp.py + python tmp.py \ No newline at end of file diff --git a/docs/source/user/first_mcdc.rst b/docs/source/user/first_mcdc.rst index 5c1d1ad92..ed76bd50f 100644 --- a/docs/source/user/first_mcdc.rst +++ b/docs/source/user/first_mcdc.rst @@ -207,7 +207,7 @@ and GPUs (AMD and Nvidia) and supports threading with MPI (Python or compiled mo Other guides are included to execute in these modes but for the sake of this first MC/DC simulation we will simply execute in Python mode (slower, no acceleration) simply with -.. code-block:: python3 +.. code-block:: bash python input.py diff --git a/extract.py b/extract.py new file mode 100644 index 000000000..5814404bb --- /dev/null +++ b/extract.py @@ -0,0 +1,15 @@ +from docutils.core import publish_doctree +from docutils import nodes +from docutils.parsers.rst import roles + +def dummy_role(name, rawtext, text, lineno, inliner, options=None, content=None): + return [nodes.Text(text)], [] + +roles.register_local_role("doc", dummy_role) + +rst = open("docs/source/user/first_mcdc.rst").read() + +for node in publish_doctree(rst).findall(nodes.literal_block): + if node.get("classes") == ['code', 'python3']: + print(node.astext()) + diff --git a/test.py b/test.py new file mode 100644 index 000000000..2f8567d7e --- /dev/null +++ b/test.py @@ -0,0 +1,166 @@ +import numpy as np + +import mcdc +# Set materials +m1 = mcdc.MaterialMG(capture=np.array([1.0])) +m2 = mcdc.MaterialMG(capture=np.array([1.5])) +m3 = mcdc.MaterialMG(capture=np.array([2.0])) +# Set surfaces +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=2.0) +s3 = mcdc.Surface.PlaneZ(z=4.0) +s4 = mcdc.Surface.PlaneZ(z=6.0, boundary_condition="vacuum") +mcdc.Cell(region=+s1 & -s2, fill=m2) +mcdc.Cell(region=+s2 & -s3, fill=m3) +mcdc.Cell(region=+s3 & -s4, fill=m1) +mcdc.Source(z=[0.0, 6.0], isotropic=True, energy_group=0) +# Tally: cell-average fluxes and collisions +mesh = mcdc.MeshStructured(z=np.linspace(0.0, 6.0, 61)) +mcdc.Tally( + mesh=mesh, + scores=["flux", "collision"], + mu=np.linspace(-1.0, 1.0, 32 + 1), +) + +# Tally: current crossing a cell boundary +mcdc.Tally( + cell=my_cell, + scores=["current-net", "current-in", "current-out"], +) +mcdc.settings.N_particle = 1000 +mcdc.run() +import numpy as np +import mcdc + +# ============================================================================= +# Set model +# ============================================================================= +# Three slab layers with different purely-absorbing materials + +# Set materials +m1 = mcdc.MaterialMG(capture=np.array([1.0])) +m2 = mcdc.MaterialMG(capture=np.array([1.5])) +m3 = mcdc.MaterialMG(capture=np.array([2.0])) + +# Set surfaces +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=2.0) +s3 = mcdc.Surface.PlaneZ(z=4.0) +s4 = mcdc.Surface.PlaneZ(z=6.0, boundary_condition="vacuum") + +# Set cells +mcdc.Cell(region=+s1 & -s2, fill=m2) +mcdc.Cell(region=+s2 & -s3, fill=m3) +mcdc.Cell(region=+s3 & -s4, fill=m1) + +# ============================================================================= +# Set source +# ============================================================================= +# Uniform isotropic source throughout the domain + +mcdc.Source(z=[0.0, 6.0], isotropic=True, energy_group=0) + +# ============================================================================= +# Set tally, setting, and run mcdc +# ============================================================================= + +# Tally: cell-average fluxes and collisions +mesh = mcdc.MeshStructured(z=np.linspace(0.0, 6.0, 61)) +mcdc.Tally( + mesh=mesh, + scores=["flux", "collision"], + mu=np.linspace(-1.0, 1.0, 32 + 1), +) + +# Setting +mcdc.settings.N_particle = 1000 + +# Run +mcdc.run() +import h5py +import numpy as np +# Load results +with h5py.File("output.h5", "r") as f: + # The tally name matches the auto-generated name (e.g., "mesh_tally_0") + tally_name = list(f["tallies"].keys())[0] + tally = f[f"tallies/{tally_name}"] + + z = tally["grid/z"][:] + dz = z[1:] - z[:-1] + z_mid = 0.5 * (z[:-1] + z[1:]) + + mu = tally["grid/mu"][:] + dmu = mu[1:] - mu[:-1] + mu_mid = 0.5 * (mu[:-1] + mu[1:]) + + psi = tally["flux/mean"][:] + psi_sd = tally["flux/sdev"][:] +from reference import reference +import matplotlib.pyplot as plt +import numpy as np + +I = len(z) - 1 +N = len(mu) - 1 + +# Scalar flux +phi = np.zeros(I) +phi_sd = np.zeros(I) +for i in range(I): + phi[i] += np.sum(psi[i, :]) + phi_sd[i] += np.linalg.norm(psi_sd[i, :]) + +# Normalize +phi /= dz +phi_sd /= dz +J /= dz +J_sd /= dz +for n in range(N): + psi[:, n] = psi[:, n] / dz / dmu[n] + psi_sd[:, n] = psi_sd[:, n] / dz / dmu[n] + +# Reference solution +phi_ref, J_ref, psi_ref = reference(z, mu) + +# Flux - spatial average +plt.plot(z_mid, phi, "-b", label="MC") +plt.fill_between(z_mid, phi - phi_sd, phi + phi_sd, alpha=0.2, color="b") +plt.plot(z_mid, phi_ref, "--r", label="Ref.") +plt.xlabel(r"$z$, cm") +plt.ylabel("Flux") +plt.ylim([0.06, 0.16]) +plt.grid() +plt.legend() +plt.title(r"$\bar{\phi}_i$") +plt.show() + +# Current - spatial average +plt.plot(z_mid, J, "-b", label="MC") +plt.fill_between(z_mid, J - J_sd, J + J_sd, alpha=0.2, color="b") +plt.plot(z_mid, J_ref, "--r", label="Ref.") +plt.xlabel(r"$z$, cm") +plt.ylabel("Current") +plt.ylim([-0.03, 0.045]) +plt.grid() +plt.legend() +plt.title(r"$\bar{J}_i$") +plt.show() + +# Angular flux - spatial average +vmin = min(np.min(psi_ref), np.min(psi)) +vmax = max(np.max(psi_ref), np.max(psi)) +fig, ax = plt.subplots(1, 2, sharey=True) +Z, MU = np.meshgrid(z_mid, mu_mid) +im = ax[0].pcolormesh(MU.T, Z.T, psi_ref, vmin=vmin, vmax=vmax) +ax[0].set_xlabel(r"Polar cosine, $\mu$") +ax[0].set_ylabel(r"$z$") +ax[0].set_title(r"\psi") +ax[0].set_title(r"$\bar{\psi}_i(\mu)$ [Ref.]") +ax[1].pcolormesh(MU.T, Z.T, psi, vmin=vmin, vmax=vmax) +ax[1].set_xlabel(r"Polar cosine, $\mu$") +ax[1].set_ylabel(r"$z$") +ax[1].set_title(r"$\bar{\psi}_i(\mu)$ [MC]") +fig.subplots_adjust(right=0.8) +cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) +cbar = fig.colorbar(im, cax=cbar_ax) +cbar.set_label("Angular flux") +plt.show() From 6bd965653e89f3a5d9e8de491986587d2901024b Mon Sep 17 00:00:00 2001 From: Kamil Kania Date: Mon, 20 Jul 2026 16:35:06 +0200 Subject: [PATCH 2/3] Changed single quotation marks to double --- .github/workflows/docs_test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/docs_test.yml b/.github/workflows/docs_test.yml index fc5a8476f..52c66bc0c 100644 --- a/.github/workflows/docs_test.yml +++ b/.github/workflows/docs_test.yml @@ -56,7 +56,7 @@ jobs: rst = open("docs/source/user/first_mcdc.rst").read() for node in publish_doctree(rst).findall(nodes.literal_block): - if node.get("classes") == ['code', 'python3']: + if node.get("classes") == ["code", "python3"]: print(node.astext()) ' > tmp.py cat -n tmp.py From f0a3b3bc2158f880ac2c7df7f192691b690e9062 Mon Sep 17 00:00:00 2001 From: Kamil Kania Date: Mon, 20 Jul 2026 16:40:44 +0200 Subject: [PATCH 3/3] Removed unintended files --- extract.py | 15 ----- test.py | 166 ----------------------------------------------------- 2 files changed, 181 deletions(-) delete mode 100644 extract.py delete mode 100644 test.py diff --git a/extract.py b/extract.py deleted file mode 100644 index 5814404bb..000000000 --- a/extract.py +++ /dev/null @@ -1,15 +0,0 @@ -from docutils.core import publish_doctree -from docutils import nodes -from docutils.parsers.rst import roles - -def dummy_role(name, rawtext, text, lineno, inliner, options=None, content=None): - return [nodes.Text(text)], [] - -roles.register_local_role("doc", dummy_role) - -rst = open("docs/source/user/first_mcdc.rst").read() - -for node in publish_doctree(rst).findall(nodes.literal_block): - if node.get("classes") == ['code', 'python3']: - print(node.astext()) - diff --git a/test.py b/test.py deleted file mode 100644 index 2f8567d7e..000000000 --- a/test.py +++ /dev/null @@ -1,166 +0,0 @@ -import numpy as np - -import mcdc -# Set materials -m1 = mcdc.MaterialMG(capture=np.array([1.0])) -m2 = mcdc.MaterialMG(capture=np.array([1.5])) -m3 = mcdc.MaterialMG(capture=np.array([2.0])) -# Set surfaces -s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") -s2 = mcdc.Surface.PlaneZ(z=2.0) -s3 = mcdc.Surface.PlaneZ(z=4.0) -s4 = mcdc.Surface.PlaneZ(z=6.0, boundary_condition="vacuum") -mcdc.Cell(region=+s1 & -s2, fill=m2) -mcdc.Cell(region=+s2 & -s3, fill=m3) -mcdc.Cell(region=+s3 & -s4, fill=m1) -mcdc.Source(z=[0.0, 6.0], isotropic=True, energy_group=0) -# Tally: cell-average fluxes and collisions -mesh = mcdc.MeshStructured(z=np.linspace(0.0, 6.0, 61)) -mcdc.Tally( - mesh=mesh, - scores=["flux", "collision"], - mu=np.linspace(-1.0, 1.0, 32 + 1), -) - -# Tally: current crossing a cell boundary -mcdc.Tally( - cell=my_cell, - scores=["current-net", "current-in", "current-out"], -) -mcdc.settings.N_particle = 1000 -mcdc.run() -import numpy as np -import mcdc - -# ============================================================================= -# Set model -# ============================================================================= -# Three slab layers with different purely-absorbing materials - -# Set materials -m1 = mcdc.MaterialMG(capture=np.array([1.0])) -m2 = mcdc.MaterialMG(capture=np.array([1.5])) -m3 = mcdc.MaterialMG(capture=np.array([2.0])) - -# Set surfaces -s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") -s2 = mcdc.Surface.PlaneZ(z=2.0) -s3 = mcdc.Surface.PlaneZ(z=4.0) -s4 = mcdc.Surface.PlaneZ(z=6.0, boundary_condition="vacuum") - -# Set cells -mcdc.Cell(region=+s1 & -s2, fill=m2) -mcdc.Cell(region=+s2 & -s3, fill=m3) -mcdc.Cell(region=+s3 & -s4, fill=m1) - -# ============================================================================= -# Set source -# ============================================================================= -# Uniform isotropic source throughout the domain - -mcdc.Source(z=[0.0, 6.0], isotropic=True, energy_group=0) - -# ============================================================================= -# Set tally, setting, and run mcdc -# ============================================================================= - -# Tally: cell-average fluxes and collisions -mesh = mcdc.MeshStructured(z=np.linspace(0.0, 6.0, 61)) -mcdc.Tally( - mesh=mesh, - scores=["flux", "collision"], - mu=np.linspace(-1.0, 1.0, 32 + 1), -) - -# Setting -mcdc.settings.N_particle = 1000 - -# Run -mcdc.run() -import h5py -import numpy as np -# Load results -with h5py.File("output.h5", "r") as f: - # The tally name matches the auto-generated name (e.g., "mesh_tally_0") - tally_name = list(f["tallies"].keys())[0] - tally = f[f"tallies/{tally_name}"] - - z = tally["grid/z"][:] - dz = z[1:] - z[:-1] - z_mid = 0.5 * (z[:-1] + z[1:]) - - mu = tally["grid/mu"][:] - dmu = mu[1:] - mu[:-1] - mu_mid = 0.5 * (mu[:-1] + mu[1:]) - - psi = tally["flux/mean"][:] - psi_sd = tally["flux/sdev"][:] -from reference import reference -import matplotlib.pyplot as plt -import numpy as np - -I = len(z) - 1 -N = len(mu) - 1 - -# Scalar flux -phi = np.zeros(I) -phi_sd = np.zeros(I) -for i in range(I): - phi[i] += np.sum(psi[i, :]) - phi_sd[i] += np.linalg.norm(psi_sd[i, :]) - -# Normalize -phi /= dz -phi_sd /= dz -J /= dz -J_sd /= dz -for n in range(N): - psi[:, n] = psi[:, n] / dz / dmu[n] - psi_sd[:, n] = psi_sd[:, n] / dz / dmu[n] - -# Reference solution -phi_ref, J_ref, psi_ref = reference(z, mu) - -# Flux - spatial average -plt.plot(z_mid, phi, "-b", label="MC") -plt.fill_between(z_mid, phi - phi_sd, phi + phi_sd, alpha=0.2, color="b") -plt.plot(z_mid, phi_ref, "--r", label="Ref.") -plt.xlabel(r"$z$, cm") -plt.ylabel("Flux") -plt.ylim([0.06, 0.16]) -plt.grid() -plt.legend() -plt.title(r"$\bar{\phi}_i$") -plt.show() - -# Current - spatial average -plt.plot(z_mid, J, "-b", label="MC") -plt.fill_between(z_mid, J - J_sd, J + J_sd, alpha=0.2, color="b") -plt.plot(z_mid, J_ref, "--r", label="Ref.") -plt.xlabel(r"$z$, cm") -plt.ylabel("Current") -plt.ylim([-0.03, 0.045]) -plt.grid() -plt.legend() -plt.title(r"$\bar{J}_i$") -plt.show() - -# Angular flux - spatial average -vmin = min(np.min(psi_ref), np.min(psi)) -vmax = max(np.max(psi_ref), np.max(psi)) -fig, ax = plt.subplots(1, 2, sharey=True) -Z, MU = np.meshgrid(z_mid, mu_mid) -im = ax[0].pcolormesh(MU.T, Z.T, psi_ref, vmin=vmin, vmax=vmax) -ax[0].set_xlabel(r"Polar cosine, $\mu$") -ax[0].set_ylabel(r"$z$") -ax[0].set_title(r"\psi") -ax[0].set_title(r"$\bar{\psi}_i(\mu)$ [Ref.]") -ax[1].pcolormesh(MU.T, Z.T, psi, vmin=vmin, vmax=vmax) -ax[1].set_xlabel(r"Polar cosine, $\mu$") -ax[1].set_ylabel(r"$z$") -ax[1].set_title(r"$\bar{\psi}_i(\mu)$ [MC]") -fig.subplots_adjust(right=0.8) -cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) -cbar = fig.colorbar(im, cax=cbar_ax) -cbar.set_label("Angular flux") -plt.show()