diff --git a/examples/eigenvalue/2d_c5g7/2d_c5g7_xs.h5 b/examples/eigenvalue/2d_c5g7/2d_c5g7_xs.h5 new file mode 100644 index 00000000..90946512 Binary files /dev/null and b/examples/eigenvalue/2d_c5g7/2d_c5g7_xs.h5 differ diff --git a/examples/eigenvalue/2d_c5g7/c5g7.h5 b/examples/eigenvalue/2d_c5g7/c5g7.h5 deleted file mode 100644 index 218b300f..00000000 Binary files a/examples/eigenvalue/2d_c5g7/c5g7.h5 and /dev/null differ diff --git a/examples/eigenvalue/2d_c5g7/input.py b/examples/eigenvalue/2d_c5g7/input.py index d2bdc838..078abab9 100644 --- a/examples/eigenvalue/2d_c5g7/input.py +++ b/examples/eigenvalue/2d_c5g7/input.py @@ -8,7 +8,7 @@ # ============================================================================= # Load material data -lib = h5py.File("c5g7.h5", "r") +lib = h5py.File("2d_c5g7_xs.h5", "r") # Materials @@ -17,22 +17,17 @@ def set_mat(mat): capture=mat["capture"][:], scatter=mat["scatter"][:], fission=mat["fission"][:], - nu_p=mat["nu_p"][:], - nu_d=mat["nu_d"][:], - chi_p=mat["chi_p"][:], - chi_d=mat["chi_d"][:], - speed=mat["speed"], - decay=mat["decay"], + nu_p=mat["nu"][:], + chi_p=mat["chi"][:], ) mat_uo2 = set_mat(lib["uo2"]) mat_mox43 = set_mat(lib["mox43"]) -mat_mox7 = set_mat(lib["mox7"]) +mat_mox7 = set_mat(lib["mox70"]) mat_mox87 = set_mat(lib["mox87"]) mat_gt = set_mat(lib["gt"]) mat_fc = set_mat(lib["fc"]) -mat_cr = set_mat(lib["cr"]) mat_mod = set_mat(lib["mod"]) # ============================================================================= @@ -52,7 +47,6 @@ def set_mat(mat): mox8 = mcdc.cell([-cy], mat_mox87) gt = mcdc.cell([-cy], mat_gt) fc = mcdc.cell([-cy], mat_fc) -cr = mcdc.cell([-cy], mat_cr) mod = mcdc.cell([+cy], mat_mod) modi = mcdc.cell([-cy], mat_mod) # For all-water lattice @@ -63,7 +57,6 @@ def set_mat(mat): n = mcdc.universe([mox8, mod])["ID"] g = mcdc.universe([gt, mod])["ID"] f = mcdc.universe([fc, mod])["ID"] -c = mcdc.universe([cr, mod])["ID"] w = mcdc.universe([modi, mod])["ID"] # ============================================================================= diff --git a/examples/eigenvalue/2d_c5g7_iqmc/2d_c5g7_xs.h5 b/examples/eigenvalue/2d_c5g7_iqmc/2d_c5g7_xs.h5 new file mode 100644 index 00000000..90946512 Binary files /dev/null and b/examples/eigenvalue/2d_c5g7_iqmc/2d_c5g7_xs.h5 differ diff --git a/examples/eigenvalue/2d_c5g7_iqmc/input.py b/examples/eigenvalue/2d_c5g7_iqmc/input.py index 1d7c4f13..4e5cfe37 100644 --- a/examples/eigenvalue/2d_c5g7_iqmc/input.py +++ b/examples/eigenvalue/2d_c5g7_iqmc/input.py @@ -164,7 +164,7 @@ def set_mat(mat): # ============================================================================= # iQMC Parameters # ============================================================================= -N = 1e3 +N = 1e4 maxit = 10 tol = 1e-3 Nx = 17 * 3 * 2 diff --git a/examples/eigenvalue/2d_c5g7_iqmc/process.py b/examples/eigenvalue/2d_c5g7_iqmc/process.py index 398e7e69..440793b1 100644 --- a/examples/eigenvalue/2d_c5g7_iqmc/process.py +++ b/examples/eigenvalue/2d_c5g7_iqmc/process.py @@ -8,10 +8,10 @@ # ============================================================================= # Load iqmc result -with h5py.File("PI_output.h5", "r") as f: +with h5py.File("output.h5", "r") as f: x = f["iqmc/grid/x"][:] y = f["iqmc/grid/y"][:] - phi_avg = f["iqmc/flux"][:] + phi_avg = f["iqmc/tally/flux"][:] f.close() diff --git a/examples/eigenvalue/slab_2gu_iqmc/input.py b/examples/eigenvalue/slab_2gu_iqmc/input.py index d42109d7..b948221a 100644 --- a/examples/eigenvalue/slab_2gu_iqmc/input.py +++ b/examples/eigenvalue/slab_2gu_iqmc/input.py @@ -9,17 +9,6 @@ # "Analytical Benchmark Test Set For Criticality Code Verification" # Set materials - -# 1G-PU Slab data -# m1 = mcdc.material( -# capture=np.array([0.019584]), -# scatter=np.array([[0.225216]]), -# fission=np.array([0.081600]), -# nu_p=np.array([2.84]), -# chi_p=np.array([1.0]), -# ) -# R = [0.0, 4.513502] - # 2G-U Slab data m1 = mcdc.material( capture=np.array([0.01344, 0.00384]), diff --git a/examples/eigenvalue/slab_kornreich_iqmc/input.py b/examples/eigenvalue/slab_kornreich_iqmc/input.py index e5ea8460..5b8a20ed 100644 --- a/examples/eigenvalue/slab_kornreich_iqmc/input.py +++ b/examples/eigenvalue/slab_kornreich_iqmc/input.py @@ -35,14 +35,13 @@ # ============================================================================= # iQMC Parameters # ============================================================================= -N = 1000 +N = 5000 maxit = 25 tol = 1e-3 x = np.arange(0.0, 2.6, 0.1) Nx = len(x) - 1 generator = "halton" solver = "power_iteration" -inner_solver = "source_iteration" fixed_source = np.zeros(Nx) phi0 = np.ones((Nx)) @@ -58,8 +57,7 @@ tol=tol, generator=generator, eigenmode_solver=solver, - fixed_source_solver=inner_solver, - # score=['tilt-x'] + score=["tilt-x"], ) # Setting mcdc.setting(N_particle=N) diff --git a/examples/eigenvalue/slab_kornreich_iqmc/process.py b/examples/eigenvalue/slab_kornreich_iqmc/process.py index 90cbe50e..708587eb 100644 --- a/examples/eigenvalue/slab_kornreich_iqmc/process.py +++ b/examples/eigenvalue/slab_kornreich_iqmc/process.py @@ -3,6 +3,26 @@ import h5py import sys +# ============================================================================= +# Import data +# ============================================================================= + +with h5py.File("output.h5", "r") as f: + # Note the spatial (dx) and source strength (100+1) normalization + keff = f["k_eff"][()] + phi_avg = f["iqmc/tally/flux"][:] + sweeps = f["iqmc/sweep_count"][()] + x = f["iqmc/grid/x"][:] + dx = x[1] - x[0] + x_mid = 0.5 * (x[:-1] + x[1:]) + f.close() + +tmp = 0.5 * (phi_avg[1:] + phi_avg[:-1]) +norm = np.sum(tmp * dx) +phi_avg /= norm +print("Keff = ", keff) +print("Number of QMC Transport Sweeps = ", sweeps) + # ============================================================================= # Reference solution # ============================================================================= @@ -69,27 +89,7 @@ plt.figure(dpi=300, figsize=(8, 5)) plt.plot(x_exact, phi_exact, label="sol") -with h5py.File("output.h5", "r") as f: - # Note the spatial (dx) and source strength (100+1) normalization - keff = f["k_eff"][()] - phi_avg = f["iqmc/flux"][:] - sweeps = f["iqmc/sweep_count"][()] - x = f["iqmc/grid/x"][:] - dx = x[1] - x[0] - x_mid = 0.5 * (x[:-1] + x[1:]) - f.close() - -tmp = 0.5 * (phi_avg[1:] + phi_avg[:-1]) -norm = np.sum(tmp * dx) -phi_avg /= norm -print("Keff = ", keff) -print(sweeps) -plt.plot(x_mid, phi_avg, label="PI") - - -# ============================================================================= -# Finish Plot -# ============================================================================= +plt.plot(x_mid, phi_avg) plt.title("Kornreich et al. Slab") plt.ylabel(r"$\phi(x)$") plt.xlabel(r"$x$") diff --git a/examples/fixed_source/cooper2_iqmc/input.py b/examples/fixed_source/cooper2_iqmc/input.py index 3bd5b5a3..56c0d0b6 100644 --- a/examples/fixed_source/cooper2_iqmc/input.py +++ b/examples/fixed_source/cooper2_iqmc/input.py @@ -35,13 +35,12 @@ # ============================================================================= N = 1e4 Nx = Ny = 40 -maxit = 15 -tol = 1e-9 +maxit = 30 +tol = 1e-6 x = np.linspace(0, 4, num=Nx + 1) y = np.linspace(0, 4, num=Ny + 1) generator = "halton" solver = "gmres" -tilt = 1 # fixed source in lower left corner fixed_source = np.zeros((Nx, Ny)) @@ -57,8 +56,8 @@ maxitt=maxit, tol=tol, generator=generator, - source_tilt=tilt, fixed_source_solver=solver, + score=["tilt-x", "tilt-y"], ) # ============================================================================= diff --git a/examples/fixed_source/cooper2_iqmc/process.py b/examples/fixed_source/cooper2_iqmc/process.py index 1c58da05..918b15b9 100644 --- a/examples/fixed_source/cooper2_iqmc/process.py +++ b/examples/fixed_source/cooper2_iqmc/process.py @@ -1,36 +1,14 @@ import matplotlib.pyplot as plt - -plt.rcParams.update({"font.size": 14}) - -import matplotlib as mpl import h5py import numpy as np -# from mpl_toolkits.mplot3d import axes3d, Axes3D - - -# Load iqmc result -# with h5py.File("output.h5", "r") as f: -# x = f["iqmc/grid/x"][:] -# dx = [x[1:] - x[:-1]][-1] -# x_mid = 0.5 * (x[:-1] + x[1:]) - -# phi = f["iqmc/flux"][:] -# f.close() - # Load iqmc result -with h5py.File("output2.h5", "r") as f: +with h5py.File("output.h5", "r") as f: meshx = f["iqmc/grid/x"][:] meshy = f["iqmc/grid/y"][:] dx = [meshx[1:] - meshx[:-1]][-1] x_mid = 0.5 * (meshx[:-1] + meshx[1:]) - lowX = meshx[:-1] - highX = meshx[1:] - - phi = f["iqmc/flux"][:] - q = f["iqmc/source"][:][0, 0, :, :, 0] - qdotx = f["iqmc/Q11/x"][:][0, 0, :, :, 0] - qdoty = f["iqmc/Q11/x"][:][0, 0, :, :, 0] + phi = f["iqmc/tally/flux"][:] f.close() @@ -46,115 +24,4 @@ ax.set_zlabel(r"log($\phi$)", rotation=180) ax.view_init(elev=15, azim=20) -# iqmc_label = mpl.lines.Line2D([0], [0], linestyle="none", c="b", marker="o") -# mc_label = mpl.lines.Line2D([0], [0], linestyle="none", c="r", marker="o") -# ax.legend([iqmc_label, mc_label], ["iQMC", "MC/DC"], numpoints=1) -# plt.tight_layout() -plt.show() - -# ============================================================================= -# Source Plots -# ============================================================================= -X, Y = np.meshgrid(x_mid, x_mid) -Z = np.log10(q) -fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10)) -ax.plot_surface(Y, X, Z, edgecolor="b", linewidth=0.5, cmap="viridis") -ax.set_xlabel("x") -ax.set_ylabel("y") -ax.set_zlabel(r"Source", rotation=180) -ax.view_init(elev=15, azim=20) - - -num_steps = 800 -y_mid = x_mid.copy() -Nx = num_steps -Ny = num_steps -x = np.linspace(lowX[0], highX[-1], num=num_steps) -y = np.linspace(lowX[0], highX[-1], num=num_steps) -# Nx = meshx.size -# Ny = meshy.size -# x = meshx -# y = meshy -z1 = np.zeros(shape=(Nx, Ny)) -z2 = np.zeros(shape=(Nx, Ny)) -z3 = np.zeros(shape=(Nx, Ny)) -z4 = np.zeros(shape=(Nx, Ny)) -points = np.zeros(shape=(int(Nx * Ny), 2)) -count = 0 -for i in range(Nx): - zonex = np.argmax((x[i] > lowX) * (x[i] <= highX)) - midx = x_mid[zonex] - for j in range(Ny): - zoney = np.argmax((y[j] > lowX) * (y[j] <= highX)) - midy = y_mid[zoney] - points[count, :] = [x[i], y[j]] - count += 1 - z1[i, j] = q[zonex, zoney] - z2[i, j] = q[zonex, zoney] + qdotx[zonex, zoney] * (x[i] - midx) - z3[i, j] = q[zonex, zoney] + qdoty[zonex, zoney] * (y[j] - midy) - z4[i, j] = ( - q[zonex, zoney] - + qdotx[zonex, zoney] * (x[i] - midx) - + qdoty[zonex, zoney] * (y[j] - midy) - ) - -X, Y = np.meshgrid(x, y) -Z1 = np.log10(z1) -fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10)) -ax.plot_surface(Y, X, z1, edgecolor="b", linewidth=0.5, cmap="viridis") -ax.set_xlabel("x") -ax.set_ylabel("y") -ax.set_zlabel(r"Source", rotation=180) -ax.view_init(elev=15, azim=20) -plt.show() - -Z2 = np.log10(z2) -fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10)) -ax.plot_surface(Y, X, z2, edgecolor="b", linewidth=0.5, cmap="viridis") -ax.set_xlabel("x") -ax.set_ylabel("y") -ax.set_zlabel(r"Source - Tilted X", rotation=180) -ax.view_init(elev=15, azim=20) -plt.show() - -Z3 = np.log10(z3) -fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10)) -ax.plot_surface(Y, X, z3, edgecolor="b", linewidth=0.5, cmap="viridis") -ax.set_xlabel("x") -ax.set_ylabel("y") -ax.set_zlabel(r"Source - Tilted Y", rotation=180) -ax.view_init(elev=15, azim=20) -plt.show() - -Z4 = np.log10(z4) -fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10)) -ax.plot_surface(Y, X, z4, edgecolor="b", linewidth=0.5, cmap="viridis") -ax.set_xlabel("x") -ax.set_ylabel("y") -ax.set_zlabel(r"Source - Tilted X,Y", rotation=180) -ax.view_init(elev=15, azim=20) -plt.show() - -# ============================================================================= -# Scipy Interpolate plots -# ============================================================================= - -from scipy.interpolate import RegularGridInterpolator - -X, Y = np.meshgrid(x_mid, x_mid) -interp = RegularGridInterpolator( - (x_mid, x_mid), q, method="linear", bounds_error=False, fill_value=None -) -result = interp(points) -result = np.reshape(result, (Nx, Ny)) -X, Y = np.meshgrid(x, y) -fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, dpi=300, figsize=(12, 10)) -ax.plot_surface(Y, X, np.log10(result), edgecolor="b", linewidth=0.5, cmap="viridis") -ax.set_xlabel("x") -ax.set_ylabel("y") -ax.set_zlabel(r"Source") -ax.view_init(elev=15, azim=20) -# plt.xlim((0,10)) -# plt.ylim((0,10)) -plt.title("Source - Scipy Interp. X, Y, XY") plt.show() diff --git a/examples/fixed_source/inf_hdpe_iqmc/HDPE/D_12G_HDPE.csv b/examples/fixed_source/inf_hdpe_iqmc/HDPE/D_12G_HDPE.csv deleted file mode 100644 index 9be3ca93..00000000 --- a/examples/fixed_source/inf_hdpe_iqmc/HDPE/D_12G_HDPE.csv +++ /dev/null @@ -1,12 +0,0 @@ -6.269442407290791053e+00 -5.393556856621185780e+00 -4.274584869582430002e+00 -2.870042958497598384e+00 -1.883455832731082191e+00 -1.203933532770413928e+00 -7.292351459099583044e-01 -4.799050830986973382e-01 -3.635901034315170977e-01 -3.155700105004432543e-01 -2.970505307233495818e-01 -2.893465974506593263e-01 diff --git a/examples/fixed_source/inf_hdpe_iqmc/HDPE/Scat_12G_HDPE.csv b/examples/fixed_source/inf_hdpe_iqmc/HDPE/Scat_12G_HDPE.csv deleted file mode 100644 index a7a29a64..00000000 --- a/examples/fixed_source/inf_hdpe_iqmc/HDPE/Scat_12G_HDPE.csv +++ /dev/null @@ -1,12 +0,0 @@ -1.334151993189741300e-02,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 -1.110650719937275599e-02,1.484916919295316273e-02,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 -1.053837413948423458e-02,1.870389755333026760e-02,2.253221526634388572e-02,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 -8.567210411278577617e-03,1.445275107954792360e-02,2.989333155960674851e-02,4.662811179135290668e-02,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 -4.100305239438383452e-03,6.277755132324696653e-03,1.388445482413599476e-02,4.050851568681469728e-02,6.830348806947485196e-02,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 -2.699131538537213101e-03,3.876106808595280399e-03,7.234055136811592089e-03,2.055544038772569651e-02,7.140274085797508608e-02,1.282851746225865674e-01,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 -1.232602951783819687e-03,1.611284864892969680e-03,2.591424755391343385e-03,7.597068544289146096e-03,2.499072058264244223e-02,9.925423388777832034e-02,2.081172197778242416e-01,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 -4.124100622592566595e-04,5.610792256655693739e-04,9.537006738845994151e-04,2.777942594939304243e-03,9.212126316212066249e-03,3.499212930004284849e-02,1.645866901017251638e-01,3.098376823745143671e-01,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 -1.269631712627834524e-04,1.874394874881495847e-04,3.497909936655281657e-04,1.016901202879264780e-03,3.388132579926314732e-03,1.287102888213230087e-02,5.850494914972402538e-02,2.499335598841757045e-01,4.061470151844644993e-01,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 -4.146669257787518264e-05,6.482893694346368802e-05,1.279564254308332546e-04,3.716503739851435340e-04,1.241373087200701563e-03,4.715963703931650049e-03,2.143766568623491769e-02,8.932269185453069194e-02,3.269727852597164031e-01,4.677505835076816165e-01,-0.000000000000000000e+00,-0.000000000000000000e+00 -1.446866336083170517e-05,2.325226208930476768e-05,4.705339633955184063e-05,1.366266833613154674e-04,4.568220819077066486e-04,1.735484664356683410e-03,7.889279017900367730e-03,3.287225379229724576e-02,1.179691845796226185e-01,3.751320745062756701e-01,4.981553769961465927e-01,-0.000000000000000000e+00 -7.090191146320424270e-06,1.153951502062724407e-05,2.360856783335970704e-05,6.854340796791712307e-05,2.292836838304767482e-04,8.710696050484713175e-04,3.959788896768070936e-03,1.649937630428657820e-02,5.921179016651618304e-02,1.850030256356753888e-01,5.403779829123985801e-01,8.081681095670927295e-01 diff --git a/examples/fixed_source/inf_hdpe_iqmc/HDPE/Siga_12G_HDPE.csv b/examples/fixed_source/inf_hdpe_iqmc/HDPE/Siga_12G_HDPE.csv deleted file mode 100644 index 8ce56620..00000000 --- a/examples/fixed_source/inf_hdpe_iqmc/HDPE/Siga_12G_HDPE.csv +++ /dev/null @@ -1,12 +0,0 @@ -1.039773737506125310e-03 -1.298598835897585291e-03 -1.251628915459136954e-03 -2.232399903718903303e-06 -2.354251726366774652e-06 -2.220566003781338333e-06 -2.810203855975818708e-06 -5.465927093183564288e-06 -1.465115894524601705e-04 --3.686367604710159840e-03 --1.413002836367041766e-04 -7.135437977273252999e-05 diff --git a/examples/fixed_source/inf_hdpe_iqmc/HDPE/group_centers_12G_HDPE.csv b/examples/fixed_source/inf_hdpe_iqmc/HDPE/group_centers_12G_HDPE.csv deleted file mode 100644 index a3042ae2..00000000 --- a/examples/fixed_source/inf_hdpe_iqmc/HDPE/group_centers_12G_HDPE.csv +++ /dev/null @@ -1,12 +0,0 @@ -1.525000000000000000e+01 -1.175000000000000000e+01 -8.035000000000000142e+00 -4.467499999999999361e+00 -2.108999999999999986e+00 -9.264999999999999902e-01 -3.419999999999999707e-01 -1.257999999999999952e-01 -4.619999999999999801e-02 -1.695999999999999938e-02 -6.234999999999999223e-03 -1.901999999999999828e-03 diff --git a/examples/fixed_source/inf_hdpe_iqmc/HDPE/group_edges_12G_HDPE.csv b/examples/fixed_source/inf_hdpe_iqmc/HDPE/group_edges_12G_HDPE.csv deleted file mode 100644 index 822d2c32..00000000 --- a/examples/fixed_source/inf_hdpe_iqmc/HDPE/group_edges_12G_HDPE.csv +++ /dev/null @@ -1,13 +0,0 @@ -1.700000000000000000e+01 -1.350000000000000000e+01 -1.000000000000000000e+01 -6.069999999999999396e+00 -2.864999999999999769e+00 -1.352999999999999980e+00 -5.000000000000000000e-01 -1.839999999999999969e-01 -6.759999999999999343e-02 -2.479999999999999913e-02 -9.119999999999999635e-03 -3.349999999999999679e-03 -4.539999999999999765e-04 diff --git a/examples/fixed_source/inf_hdpe_iqmc/HDPE/nuSigf_12G_HDPE.csv b/examples/fixed_source/inf_hdpe_iqmc/HDPE/nuSigf_12G_HDPE.csv deleted file mode 100644 index ec17ce77..00000000 --- a/examples/fixed_source/inf_hdpe_iqmc/HDPE/nuSigf_12G_HDPE.csv +++ /dev/null @@ -1,12 +0,0 @@ --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 --0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00,-0.000000000000000000e+00 diff --git a/examples/fixed_source/inf_hdpe_iqmc/input.py b/examples/fixed_source/inf_hdpe_iqmc/input.py deleted file mode 100644 index de9b4e30..00000000 --- a/examples/fixed_source/inf_hdpe_iqmc/input.py +++ /dev/null @@ -1,69 +0,0 @@ -import numpy as np -import os -import mcdc - -# Infinite medium of high-density polyethylene (HDPE) - -# ============================================================================= -# Import Cross Section Data -# ============================================================================= -G = 12 # G must be 12, 70, or 618 - -script_dir = os.path.dirname(__file__) # <-- absolute dir the script is in -rel_path = "./HDPE/" -abs_file_path = os.path.join(script_dir, rel_path) -D = np.genfromtxt(abs_file_path + "D_{}G_HDPE.csv".format(G), delimiter=",") -SigmaA = np.genfromtxt(abs_file_path + "Siga_{}G_HDPE.csv".format(G), delimiter=",") -SigmaS = np.genfromtxt(abs_file_path + "Scat_{}G_HDPE.csv".format(G), delimiter=",") -SigmaS = np.flip(SigmaS, 1) - -# ============================================================================= -# Set Model -# ============================================================================= - -# x-bounds -LB = 0.0 -RB = 5.0 -# Set material -m1 = mcdc.material(capture=SigmaA, scatter=SigmaS) - -# Set surfaces -s1 = mcdc.surface("plane-x", x=LB, bc="reflective") -s2 = mcdc.surface("plane-x", x=RB, bc="reflective") - -# Set cells -mcdc.cell([+s1, -s2], m1) - - -# ============================================================================= -# iQMC Parameters -# ============================================================================= - -Nx = 5 -fixed_source = np.ones((G, Nx)) -# material_idx = np.zeros(Nx, dtype=int) -phi0 = np.ones((G, Nx)) - -mcdc.iQMC( - g=np.ones((0, G)), - x=np.linspace(LB, RB, num=Nx + 1), - fixed_source=fixed_source, - phi0=phi0, - maxitt=25, - tol=1e-3, - generator="halton", -) - -# ============================================================================= -# Set tally, setting, and run mcdc -# ============================================================================= - -# weight roulette -chance = 0.9 -threshold = 1e-3 -mcdc.weight_roulette(chance, threshold) -# Setting -mcdc.setting(N_particle=1e3) - -# Run -mcdc.run() diff --git a/examples/fixed_source/inf_hdpe_iqmc/process.py b/examples/fixed_source/inf_hdpe_iqmc/process.py deleted file mode 100644 index fb71aee3..00000000 --- a/examples/fixed_source/inf_hdpe_iqmc/process.py +++ /dev/null @@ -1,27 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import h5py - - -with h5py.File("output.h5", "r") as f: - phi = f["tally/iqmc_flux"][:] - x = f["iqmc/grid/x"][:] - dx = x[1] - x[0] - x_mid = 0.5 * (x[:-1] + x[1:]) - f.close() - - -# ============================================================================= -# Plot -# ============================================================================= - -Nx = phi.shape[1] -G = phi.shape[0] -# Flux - spatial average -plt.figure(dpi=300, figsize=(8, 5)) -# plt.plot(x_mid,phi_ref,label='Sol') -for i in range(G): - plt.plot(x_mid, phi[i, :], label="iQMC") -plt.ylabel(r"$\phi(x)$") -plt.xlabel(r"$x$") -plt.grid() diff --git a/examples/fixed_source/slab_reed_iqmc/input.py b/examples/fixed_source/slab_reed_iqmc/input.py index e35146f9..e470fb6d 100644 --- a/examples/fixed_source/slab_reed_iqmc/input.py +++ b/examples/fixed_source/slab_reed_iqmc/input.py @@ -37,8 +37,8 @@ # ============================================================================= # iQMC Parameters # ============================================================================= -N = 1000 -Nx = 32 +N = 10000 +Nx = 64 maxit = 20 tol = 1e-4 x = np.linspace(-8, 8, num=Nx + 1) diff --git a/examples/fixed_source/slab_reed_iqmc/process.py b/examples/fixed_source/slab_reed_iqmc/process.py index d1955653..f393b0c1 100644 --- a/examples/fixed_source/slab_reed_iqmc/process.py +++ b/examples/fixed_source/slab_reed_iqmc/process.py @@ -3,15 +3,27 @@ import h5py from scipy.integrate import quad +# ============================================================================= +# Import data +# ============================================================================= + +with h5py.File("output.h5", "r") as f: + phi = f["iqmc/tally/flux"][:] + x = f["iqmc/grid/x"][:] + mesh = f["iqmc/grid/x"][:] + dx = x[1] - x[0] + x_mid = 0.5 * (x[:-1] + x[1:]) + f.close() + # ============================================================================= # Reference solution (not accurate enough for N_hist > 1E7) # ============================================================================= def reeds_sol(Nx=80, LB=-8.0, RB=8.0): - # ============================================================================= + # ========================================================================= # Reference solution - # ============================================================================= + # ========================================================================= def phi1(x): return ( 1.0 @@ -105,20 +117,6 @@ def f_phi(x1, x2): return phi_ref -with h5py.File("output.h5", "r") as f: - phi = f["iqmc/tally/flux"][:] - # q = f["iqmc/tally/source/constant"][:][0, 0, :, 0, 0] - # qdot = f["iqmc/tally/source_x"][:][0, 0, :, 0, 0] - x = f["iqmc/grid/x"][:] - mesh = f["iqmc/grid/x"][:] - dx = x[1] - x[0] - x_mid = 0.5 * (x[:-1] + x[1:]) - lowX = x[:-1] - highX = x[1:] - Nx = x_mid.size - f.close() - -Nx = 27 phi_ref = reeds_sol(Nx=x_mid.size, LB=-8.0, RB=8.0) @@ -128,33 +126,9 @@ def f_phi(x1, x2): # Flux - spatial average plt.figure(dpi=300, figsize=(8, 5)) -# plt.plot(x_mid,phi_ref,label='Sol') +plt.plot(x_mid, phi_ref, label="Sol") plt.plot(x_mid, phi, label="iQMC") plt.ylabel(r"$\phi(x)$") plt.xlabel(r"$x$") plt.grid() - -# ============================================================================= -# Plot piecewise source -# ============================================================================= -source_tilt = False - -# plt.figure(figsize=(6, 4), dpi=300) -# plt.title("MCDC") -# x = np.linspace(4, 8, num=1000) -# n = len(x) -# conditions = [(mesh[i] <= x) & (x <= mesh[i + 1]) for i in range(48, Nx)] -# y1 = np.piecewise(x, conditions, q[48:Nx]) -# plt.plot(x, y1, label=r"$a_j$") -# if source_tilt: -# y2 = np.zeros_like(x) -# for i in range(n): -# zone = np.argmax((x[i] > lowX) * (x[i] <= highX)) -# mid = x_mid[zone] -# y2[i] = q[zone] + qdot[zone] * (x[i] - mid) -# plt.plot(x, y2, label=r"$a_j + b_j(x)$") -# for i in range(48, Nx): -# plt.axvline(mesh[i], linestyle="--", color="black") -# plt.legend() -# plt.tight_layout() -# plt.ylim((-0.1, 3.0)) +plt.legend()