Skip to content

Commit 122d2e0

Browse files
ilhamvilhamv
authored
Merge pull request #94 from spasmann/main
Automated iQMC Material Index
2 parents a5629d7 + 9c87821 commit 122d2e0

File tree

15 files changed

+346
-68
lines changed

15 files changed

+346
-68
lines changed

‎examples/eigenvalue/2d_c5g7_iqmc/input.py‎

Lines changed: 209 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,209 @@
1+
import h5py
2+
import numpy as np
3+
4+
import mcdc
5+
6+
# =============================================================================
7+
# Materials
8+
# =============================================================================
9+
10+
# Load material data
11+
lib = h5py.File("c5g7.h5", "r")
12+
13+
14+
# Materials
15+
def set_mat(mat):
16+
return mcdc.material(
17+
capture=mat["capture"][:],
18+
scatter=mat["scatter"][:],
19+
fission=mat["fission"][:],
20+
nu_p=mat["nu_p"][:],
21+
nu_d=mat["nu_d"][:],
22+
chi_p=mat["chi_p"][:],
23+
chi_d=mat["chi_d"][:],
24+
speed=mat["speed"],
25+
decay=mat["decay"],
26+
)
27+
28+
29+
mat_uo2 = set_mat(lib["uo2"])
30+
mat_mox43 = set_mat(lib["mox43"])
31+
mat_mox7 = set_mat(lib["mox7"])
32+
mat_mox87 = set_mat(lib["mox87"])
33+
mat_gt = set_mat(lib["gt"])
34+
mat_fc = set_mat(lib["fc"])
35+
mat_cr = set_mat(lib["cr"])
36+
mat_mod = set_mat(lib["mod"])
37+
38+
# =============================================================================
39+
# Pin cells
40+
# =============================================================================
41+
42+
pitch = 1.26
43+
radius = 0.54
44+
45+
# Surfaces
46+
cy = mcdc.surface("cylinder-z", center=[0.0, 0.0], radius=radius)
47+
48+
# Cells
49+
uo2 = mcdc.cell([-cy], mat_uo2)
50+
mox4 = mcdc.cell([-cy], mat_mox43)
51+
mox7 = mcdc.cell([-cy], mat_mox7)
52+
mox8 = mcdc.cell([-cy], mat_mox87)
53+
gt = mcdc.cell([-cy], mat_gt)
54+
fc = mcdc.cell([-cy], mat_fc)
55+
cr = mcdc.cell([-cy], mat_cr)
56+
mod = mcdc.cell([+cy], mat_mod)
57+
modi = mcdc.cell([-cy], mat_mod) # For all-water lattice
58+
59+
# Universes
60+
u = mcdc.universe([uo2, mod])["ID"]
61+
l = mcdc.universe([mox4, mod])["ID"]
62+
m = mcdc.universe([mox7, mod])["ID"]
63+
n = mcdc.universe([mox8, mod])["ID"]
64+
g = mcdc.universe([gt, mod])["ID"]
65+
f = mcdc.universe([fc, mod])["ID"]
66+
c = mcdc.universe([cr, mod])["ID"]
67+
w = mcdc.universe([modi, mod])["ID"]
68+
69+
# =============================================================================
70+
# Assemblies
71+
# =============================================================================
72+
73+
# Lattices
74+
lattice_uo2 = mcdc.lattice(
75+
x=[-pitch * 17 / 2, pitch, 17],
76+
y=[-pitch * 17 / 2, pitch, 17],
77+
universes=[
78+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
79+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
80+
[u, u, u, u, u, g, u, u, g, u, u, g, u, u, u, u, u],
81+
[u, u, u, g, u, u, u, u, u, u, u, u, u, g, u, u, u],
82+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
83+
[u, u, g, u, u, g, u, u, g, u, u, g, u, u, g, u, u],
84+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
85+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
86+
[u, u, g, u, u, g, u, u, f, u, u, g, u, u, g, u, u],
87+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
88+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
89+
[u, u, g, u, u, g, u, u, g, u, u, g, u, u, g, u, u],
90+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
91+
[u, u, u, g, u, u, u, u, u, u, u, u, u, g, u, u, u],
92+
[u, u, u, u, u, g, u, u, g, u, u, g, u, u, u, u, u],
93+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
94+
[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],
95+
],
96+
)
97+
98+
lattice_mox = mcdc.lattice(
99+
x=[-pitch * 17 / 2, pitch, 17],
100+
y=[-pitch * 17 / 2, pitch, 17],
101+
universes=[
102+
[l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
103+
[l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
104+
[l, m, m, m, m, g, m, m, g, m, m, g, m, m, m, m, l],
105+
[l, m, m, g, m, n, n, n, n, n, n, n, m, g, m, m, l],
106+
[l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
107+
[l, m, g, n, n, g, n, n, g, n, n, g, n, n, g, m, l],
108+
[l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
109+
[l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
110+
[l, m, g, n, n, g, n, n, f, n, n, g, n, n, g, m, l],
111+
[l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
112+
[l, m, m, n, n, n, n, n, n, n, n, n, n, n, m, m, l],
113+
[l, m, g, n, n, g, n, n, g, n, n, g, n, n, g, m, l],
114+
[l, m, m, m, n, n, n, n, n, n, n, n, n, m, m, m, l],
115+
[l, m, m, g, m, n, n, n, n, n, n, n, m, g, m, m, l],
116+
[l, m, m, m, m, g, m, m, g, m, m, g, m, m, m, m, l],
117+
[l, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, l],
118+
[l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l, l],
119+
],
120+
)
121+
122+
lattice_mod = mcdc.lattice(
123+
x=[-pitch * 17 / 2, pitch * 17, 1],
124+
y=[-pitch * 17 / 2, pitch * 17, 1],
125+
universes=[[w]],
126+
)
127+
128+
# Assembly cells
129+
# Surfaces
130+
x0 = mcdc.surface("plane-x", x=-pitch * 17 / 2)
131+
x1 = mcdc.surface("plane-x", x=pitch * 17 / 2)
132+
y0 = mcdc.surface("plane-y", y=-pitch * 17 / 2)
133+
y1 = mcdc.surface("plane-y", y=pitch * 17 / 2)
134+
# Cells
135+
assembly_uo2 = mcdc.cell([+x0, -x1, +y0, -y1], lattice_uo2)
136+
assembly_mox = mcdc.cell([+x0, -x1, +y0, -y1], lattice_mox)
137+
assembly_mod = mcdc.cell([+x0, -x1, +y0, -y1], lattice_mod)
138+
139+
# Set assemblies in their respective universes
140+
u_ = mcdc.universe([assembly_uo2])["ID"]
141+
m_ = mcdc.universe([assembly_mox])["ID"]
142+
w_ = mcdc.universe([assembly_mod])["ID"]
143+
144+
# =============================================================================
145+
# Root universe: core
146+
# =============================================================================
147+
148+
# Lattice
149+
lattice_core = mcdc.lattice(
150+
x=[-pitch * 17 * 3 / 2, pitch * 17, 3],
151+
y=[-pitch * 17 * 3 / 2, pitch * 17, 3],
152+
universes=[[u_, m_, w_], [m_, u_, w_], [w_, w_, w_]],
153+
)
154+
155+
# Core cell
156+
# Surfaces
157+
x0_ = mcdc.surface("plane-x", x=0.0, bc="reflective")
158+
x1_ = mcdc.surface("plane-x", x=pitch * 17 * 3, bc="vacuum")
159+
y0_ = mcdc.surface("plane-y", y=-pitch * 17 * 3, bc="vacuum")
160+
y1_ = mcdc.surface("plane-y", y=0.0, bc="reflective")
161+
# Cell
162+
core = mcdc.cell(
163+
[+x0_, -x1_, +y0_, -y1_],
164+
lattice_core,
165+
lattice_center=[pitch * 17 * 3 / 2, -pitch * 17 * 3 / 2, 0.0],
166+
)
167+
168+
# Root universe
169+
mcdc.universe([core], root=True)
170+
171+
# =============================================================================
172+
# iQMC Parameters
173+
# =============================================================================
174+
N = 500
175+
maxit = 1
176+
tol = 1e-3
177+
x_grid = np.linspace(0.0, pitch * 17 * 3, 17 * 3 + 1)
178+
y_grid = np.linspace(-pitch * 17 * 3, 0.0, 17 * 3 + 1)
179+
# x_grid = np.linspace(0.0, pitch * 17 * 3, 20)
180+
# y_grid = np.linspace(-pitch * 17 * 3, 0.0, 20)
181+
generator = "halton"
182+
solver = "power_iteration"
183+
184+
phi0 = np.ones((x_grid.size - 1, y_grid.size - 1))
185+
fixed_source = np.zeros_like(phi0)
186+
187+
mcdc.iQMC(
188+
x=x_grid,
189+
y=y_grid,
190+
phi0=phi0,
191+
fixed_source=fixed_source,
192+
eigenmode_solver=solver,
193+
# preconditioner_sweeps=4,
194+
maxitt=maxit,
195+
tol=tol,
196+
generator=generator,
197+
)
198+
199+
# =============================================================================
200+
# run mcdc
201+
# =============================================================================
202+
203+
204+
# Setting
205+
mcdc.setting(N_particle=N, output="davidson_output")
206+
mcdc.eigenmode()
207+
208+
# Run
209+
mcdc.run()

‎examples/eigenvalue/2d_c5g7_iqmc/process.py‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
import numpy as np
2+
import matplotlib.pyplot as plt
3+
import h5py
4+
5+
6+
# =============================================================================
7+
# Plot results
8+
# =============================================================================
9+
10+
# Load iqmc result
11+
with h5py.File("davidson_output.h5", "r") as f:
12+
x = f["iqmc/grid/x"][:]
13+
y = f["iqmc/grid/y"][:]
14+
phi_avg = f["tally/iqmc_flux"][:]
15+
f.close()
16+
17+
18+
dx = x[1] - x[0]
19+
x_mid = 0.5 * (x[1:] + x[:-1])
20+
y_mid = 0.5 * (y[1:] + y[:-1])
21+
Y, X = np.meshgrid(x_mid, y_mid)
22+
23+
norm = np.sum(phi_avg)
24+
phi_avg = phi_avg.sum(axis=0) / norm
25+
26+
plt.figure(dpi=300)
27+
plt.pcolormesh(X, Y, phi_avg, shading="nearest")
28+
plt.colorbar()
29+
ax = plt.gca()
30+
ax.set_aspect("equal")
31+
plt.xlabel(r"$x$ [cm]")
32+
plt.ylabel(r"$y$ [cm]")
33+
plt.title(r"Neutron Flux")
34+
plt.show()
35+
36+
37+
Z = np.log10(np.abs(phi_avg / phi_avg.min()))
38+
fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10))
39+
ax.plot_surface(Y, X, Z, edgecolor="b", color="white", linewidth=0.5)
40+
41+
ax.set_xlabel("x")
42+
ax.set_ylabel("y")
43+
ax.set_zlabel(r"log($\phi$)")
44+
45+
ax.view_init(elev=15, azim=20)

‎examples/eigenvalue/slab_2gu_iqmc/input.py‎

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -39,14 +39,13 @@
3939
# =============================================================================
4040
# iQMC Parameters
4141
# =============================================================================
42-
Nx = 10
42+
Nx = 5
4343
N = 1e3
4444
maxit = 10
4545
tol = 1e-3
4646
x = np.linspace(0.0, 6.01275, num=Nx + 1)
4747
generator = "halton"
4848
fixed_source = np.zeros(Nx)
49-
material_idx = np.zeros(Nx, dtype=int)
5049
phi0 = np.ones((Nx))
5150

5251
# =============================================================================
@@ -57,7 +56,6 @@
5756
x=x,
5857
phi0=phi0,
5958
fixed_source=fixed_source,
60-
material_idx=material_idx,
6159
maxitt=maxit,
6260
tol=tol,
6361
generator=generator,

‎examples/eigenvalue/slab_kornreich_iqmc/input.py‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -42,8 +42,6 @@
4242
Nx = len(x) - 1
4343
generator = "halton"
4444
fixed_source = np.zeros(Nx)
45-
material_idx = np.zeros(Nx, dtype=int)
46-
material_idx[15:] = 1
4745
phi0 = np.ones((Nx))
4846

4947
# =============================================================================
@@ -54,7 +52,6 @@
5452
x=x,
5553
fixed_source=fixed_source,
5654
phi0=phi0,
57-
material_idx=material_idx,
5855
maxitt=maxit,
5956
tol=tol,
6057
generator=generator,

‎examples/fixed_source/cooper2_iqmc/input.py‎

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -44,10 +44,6 @@
4444
# fixed source in lower left corner
4545
fixed_source = np.zeros((Nx, Ny))
4646
fixed_source[0 : int(0.25 * Nx), 0 : int(0.25 * Nx)] = 1
47-
# m_room
48-
material_idx = np.ones((Nx, Ny), dtype=int)
49-
# m_barrier
50-
material_idx[int(0.5 * Nx) : int(0.6 * Nx), 0 : int(0.5 * Nx)] = 0
5147

5248
phi0 = np.ones((Nx, Ny))
5349

@@ -56,7 +52,6 @@
5652
y=y,
5753
fixed_source=fixed_source,
5854
phi0=phi0,
59-
material_idx=material_idx,
6055
maxitt=maxit,
6156
tol=tol,
6257
generator=generator,

‎examples/fixed_source/inf_hdpe_iqmc/input.py‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -41,15 +41,14 @@
4141

4242
Nx = 5
4343
fixed_source = np.ones((G, Nx))
44-
material_idx = np.zeros(Nx, dtype=int)
44+
# material_idx = np.zeros(Nx, dtype=int)
4545
phi0 = np.ones((G, Nx))
4646

4747
mcdc.iQMC(
4848
g=np.ones((0, G)),
4949
x=np.linspace(LB, RB, num=Nx + 1),
5050
fixed_source=fixed_source,
5151
phi0=phi0,
52-
material_idx=material_idx,
5352
maxitt=25,
5453
tol=1e-3,
5554
generator="halton",

‎examples/fixed_source/slab_reed_iqmc/input.py‎

Lines changed: 2 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -37,8 +37,8 @@
3737
# =============================================================================
3838
# iQMC Parameters
3939
# =============================================================================
40-
N = 1e3
41-
Nx = 32
40+
N = 1e2
41+
Nx = 16
4242
maxit = 10
4343
tol = 1e-3
4444
x = np.linspace(-8, 8, num=Nx + 1)
@@ -49,21 +49,12 @@
4949
fixed_source[int(0.125 * Nx) : int(0.1875 * Nx)] = 1.0
5050
fixed_source[int(0.8125 * Nx) : int(0.875 * Nx)] = 1.0
5151

52-
material_idx = np.zeros(Nx, dtype=int)
53-
material_idx[: int(0.1875 * Nx)] = 3
54-
material_idx[int(0.1875 * Nx) : int(0.3125 * Nx)] = 2
55-
material_idx[int(0.3125 * Nx) : int(0.375 * Nx)] = 1
56-
material_idx[int(0.625 * Nx) : int(0.6875 * Nx)] = 1
57-
material_idx[int(0.6875 * Nx) : int(0.8125 * Nx)] = 2
58-
material_idx[int(0.8125 * Nx) :] = 3
59-
6052
phi0 = np.ones((Nx))
6153

6254
mcdc.iQMC(
6355
x=x,
6456
fixed_source=fixed_source,
6557
phi0=phi0,
66-
material_idx=material_idx,
6758
maxitt=maxit,
6859
tol=tol,
6960
generator=generator,

‎mcdc/input_.py‎

Lines changed: 0 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1166,7 +1166,6 @@ def iQMC(
11661166
z=None,
11671167
phi0=None,
11681168
fixed_source=None,
1169-
material_idx=None,
11701169
scramble=False,
11711170
maxitt=25,
11721171
tol=1e-6,
@@ -1212,21 +1211,8 @@ def iQMC(
12121211
fixed_source = np.expand_dims(fixed_source, axis=ax)
12131212
phi0 = np.expand_dims(phi0, axis=ax)
12141213

1215-
ax_expand = []
1216-
if t is None:
1217-
ax_expand.append(0)
1218-
if x is None:
1219-
ax_expand.append(1)
1220-
if y is None:
1221-
ax_expand.append(2)
1222-
if z is None:
1223-
ax_expand.append(3)
1224-
for ax in ax_expand:
1225-
material_idx = np.expand_dims(material_idx, axis=ax)
1226-
12271214
card["iqmc_flux"] = phi0
12281215
card["iqmc_fixed_source"] = fixed_source
1229-
card["iqmc_material_idx"] = material_idx
12301216
card["fixed_source_solver"] = fixed_source_solver
12311217
card["eigenmode_solver"] = eigenmode_solver
12321218

0 commit comments

Comments
 (0)