diff --git a/examples/fixed_source/cooper1/input.py b/examples/fixed_source/cooper1/input.py new file mode 100644 index 00000000..082b2670 --- /dev/null +++ b/examples/fixed_source/cooper1/input.py @@ -0,0 +1,63 @@ +import numpy as np +import mcdc + + +# ============================================================================= +# Set model +# ============================================================================= +# A problem based on Problem 1 of [Coper NSE 2001] +# https://ans.tandfonline.com/action/showCitFormats?doi=10.13182/NSE00-34 + +# Set materials +SigmaT = 0.01 +c = 0.8 +m_duct = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]])) +SigmaT = 1.0 +m_room = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]])) + +# Set surfaces +sx1 = mcdc.surface("plane-x", x=0.0, bc="reflective") +sx2 = mcdc.surface("plane-x", x=2.0) +sx3 = mcdc.surface("plane-x", x=2.4) +sx4 = mcdc.surface("plane-x", x=4.0, bc="vacuum") +sy1 = mcdc.surface("plane-y", y=0.0, bc="reflective") +sy2 = mcdc.surface("plane-y", y=1.0) +sy3 = mcdc.surface("plane-y", y=2.6) +sy4 = mcdc.surface("plane-y", y=3) +sy5 = mcdc.surface("plane-y", y=4.0, bc="vacuum") + +# Set cells +# Room +mcdc.cell(+sx1 & -sx2 & +sy1 & -sy5, m_room) +mcdc.cell(+sx1 & -sx4 & +sy4 & -sy5, m_room) +mcdc.cell(+sx2 & -sx3 & +sy1 & -sy2, m_room) +mcdc.cell(+sx3 & -sx4 & +sy1 & -sy3, m_room) + +# Duct Channel +channel_1 = +sx2 & -sx3 & +sy2 & -sy3 +channel_2 = +sx2 & -sx4 & +sy3 & -sy4 +mcdc.cell(channel_1 | channel_2, m_duct) + +# ============================================================================= +# Set source +# ============================================================================= +# Uniform isotropic source throughout the domain + +mcdc.source(x=[0.0, 1.0], y=[0.0, 1.0], isotropic=True) + +# ============================================================================= +# Set tally, setting, and run mcdc +# ============================================================================= + +mcdc.tally.mesh_tally( + scores=["flux"], + x=np.linspace(0.0, 4.0, 40), + y=np.linspace(0.0, 4.0, 40), +) + +# Setting +mcdc.setting(N_particle=50) +mcdc.implicit_capture() + +# Run +mcdc.run() diff --git a/examples/fixed_source/cooper1/process.py b/examples/fixed_source/cooper1/process.py new file mode 100644 index 00000000..7bbfe3c9 --- /dev/null +++ b/examples/fixed_source/cooper1/process.py @@ -0,0 +1,38 @@ +import matplotlib.pyplot as plt +import h5py +import numpy as np + +# Load result +with h5py.File("output.h5", "r") as f: + tally = f["tallies/mesh_tally_0"] + x = tally["grid/x"][:] + dx = [x[1:] - x[:-1]][-1] + x_mid = 0.5 * (x[:-1] + x[1:]) + + phi = tally["flux/mean"][:] + phi_sd = tally["flux/sdev"][:] + + +# Plot result +X, Y = np.meshgrid(x_mid, x_mid) +Z = phi +Z = np.log10(np.abs(Z)) +fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) +ax.plot_surface(Y, X, Z, edgecolor="k", color="white") +ax.set_xlabel("x") +ax.set_ylabel("y") +ax.zaxis.set_rotate_label(False) +ax.set_zlabel("Log10 of scalar flux", rotation=90) +ax.view_init(elev=18, azim=38) +plt.show() + +Z = phi_sd / phi +Z = np.log10(np.abs(Z)) +fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) +ax.plot_surface(Y, X, Z, edgecolor="k", color="white") +ax.set_xlabel("x") +ax.set_ylabel("y") +ax.zaxis.set_rotate_label(False) +ax.set_zlabel("Log10 of scalar flux rel. stdev.", rotation=90) +ax.view_init(elev=18, azim=38) +plt.show() diff --git a/examples/fixed_source/cooper2/input.py b/examples/fixed_source/cooper2/input.py index 1c2f1b5a..b2363580 100644 --- a/examples/fixed_source/cooper2/input.py +++ b/examples/fixed_source/cooper2/input.py @@ -41,13 +41,14 @@ # Set tally, setting, and run mcdc # ============================================================================= -# Tally: cell-average and cell-edge angular fluxes and currents mcdc.tally.mesh_tally( - scores=["flux"], x=np.linspace(0.0, 4.0, 40), y=np.linspace(0.0, 4.0, 40) + scores=["flux"], + x=np.linspace(0.0, 4.0, 40), + y=np.linspace(0.0, 4.0, 40), ) # Setting -mcdc.setting(N_particle=1e2) +mcdc.setting(N_particle=50) mcdc.implicit_capture() # Run diff --git a/examples/fixed_source/cooper_combo/input.py b/examples/fixed_source/cooper_combo/input.py new file mode 100644 index 00000000..952a6083 --- /dev/null +++ b/examples/fixed_source/cooper_combo/input.py @@ -0,0 +1,68 @@ +import numpy as np +import mcdc + + +# ============================================================================= +# Set model +# ============================================================================= +# A problem based on a combination of Problems 1 & 2 of [Coper NSE 2001] +# https://ans.tandfonline.com/action/showCitFormats?doi=10.13182/NSE00-34 + +# Set materials +SigmaT = 5.0 +c = 0.8 +m_barrier = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]])) +SigmaT = 1.0 +m_room = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]])) +SigmaT = 0.01 +m_duct = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]])) + + +# Set surfaces +sx1 = mcdc.surface("plane-x", x=0.0, bc="reflective") +sx2 = mcdc.surface("plane-x", x=2.0) +sx3 = mcdc.surface("plane-x", x=2.4) +sx4 = mcdc.surface("plane-x", x=4.0, bc="vacuum") +sy1 = mcdc.surface("plane-y", y=0.0, bc="reflective") +sy2 = mcdc.surface("plane-y", y=1.0) +sy3 = mcdc.surface("plane-y", y=2.6) +sy4 = mcdc.surface("plane-y", y=3) +sy5 = mcdc.surface("plane-y", y=4.0, bc="vacuum") + +# Set cells +# Room +mcdc.cell(+sx1 & -sx2 & +sy1 & -sy5, m_room) +mcdc.cell(+sx1 & -sx4 & +sy4 & -sy5, m_room) +mcdc.cell(+sx3 & -sx4 & +sy1 & -sy3, m_room) + +# Barrier +mcdc.cell(+sx2 & -sx3 & +sy1 & -sy2, m_barrier) + +# Duct Channel +channel_1 = +sx2 & -sx3 & +sy2 & -sy3 +channel_2 = +sx2 & -sx4 & +sy3 & -sy4 +mcdc.cell(channel_1 | channel_2, m_duct) + +# ============================================================================= +# Set source +# ============================================================================= +# Uniform isotropic source throughout the domain + +mcdc.source(x=[0.0, 1.0], y=[0.0, 1.0], isotropic=True) + +# ============================================================================= +# Set tally, setting, and run mcdc +# ============================================================================= + +mcdc.tally.mesh_tally( + scores=["flux"], + x=np.linspace(0.0, 4.0, 40), + y=np.linspace(0.0, 4.0, 40), +) + +# Setting +mcdc.setting(N_particle=50) +mcdc.implicit_capture() + +# Run +mcdc.run() diff --git a/examples/fixed_source/cooper_combo/process.py b/examples/fixed_source/cooper_combo/process.py new file mode 100644 index 00000000..8ab5fad0 --- /dev/null +++ b/examples/fixed_source/cooper_combo/process.py @@ -0,0 +1,39 @@ +import matplotlib.pyplot as plt +import h5py +import numpy as np + + +# Load result +with h5py.File("output.h5", "r") as f: + tally = f["tallies/mesh_tally_0"] + x = tally["grid/x"][:] + dx = [x[1:] - x[:-1]][-1] + x_mid = 0.5 * (x[:-1] + x[1:]) + + phi = tally["flux/mean"][:] + phi_sd = tally["flux/sdev"][:] + + +# Plot result +X, Y = np.meshgrid(x_mid, x_mid) +Z = phi +Z = np.log10(np.abs(Z)) +fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) +ax.plot_surface(Y, X, Z, edgecolor="k", color="white") +ax.set_xlabel("x") +ax.set_ylabel("y") +ax.zaxis.set_rotate_label(False) +ax.set_zlabel("Log10 of scalar flux", rotation=90) +ax.view_init(elev=18, azim=38) +plt.show() + +Z = phi_sd / phi +Z = np.log10(np.abs(Z)) +fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) +ax.plot_surface(Y, X, Z, edgecolor="k", color="white") +ax.set_xlabel("x") +ax.set_ylabel("y") +ax.zaxis.set_rotate_label(False) +ax.set_zlabel("Log10 of scalar flux rel. stdev.", rotation=90) +ax.view_init(elev=18, azim=38) +plt.show() diff --git a/examples/fixed_source/kobayashi3-TD/input.py b/examples/fixed_source/kobayashi3-TD/input.py index 360e0620..4bdea735 100644 --- a/examples/fixed_source/kobayashi3-TD/input.py +++ b/examples/fixed_source/kobayashi3-TD/input.py @@ -30,17 +30,17 @@ # Set cells # Source -mcdc.cell(+sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2, m) +source_cell = mcdc.cell(+sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2, m) # Voids channel_1 = +sx1 & -sx2 & +sy2 & -sy3 & +sz1 & -sz2 channel_2 = +sx1 & -sx3 & +sy3 & -sy4 & +sz1 & -sz2 channel_3 = +sx3 & -sx4 & +sy3 & -sy4 & +sz1 & -sz3 channel_4 = +sx3 & -sx4 & +sy3 & -sy5 & +sz3 & -sz4 void_channel = channel_1 | channel_2 | channel_3 | channel_4 -mcdc.cell(void_channel, m_void) +void_cell = mcdc.cell(void_channel, m_void) # Shield box = +sx1 & -sx5 & +sy1 & -sy5 & +sz1 & -sz5 -mcdc.cell(box & ~void_channel, m) +shield_cell = mcdc.cell(box & ~void_channel, m) # ============================================================================= # Set source @@ -58,13 +58,17 @@ # Tally: z-integrated flux (X-Y section view) mcdc.tally.mesh_tally( scores=["flux"], - x=np.linspace(0.0, 60.0, 61), - y=np.linspace(0.0, 100.0, 101), + x=np.linspace(0.0, 60.0, 31), + y=np.linspace(0.0, 100.0, 51), t=np.linspace(0.0, 200.0, 21), ) +mcdc.tally.cell_tally(source_cell, scores=["flux"]) +mcdc.tally.cell_tally(void_cell, scores=["flux"]) +mcdc.tally.cell_tally(shield_cell, scores=["flux"]) + # Setting -mcdc.setting(N_particle=1e2) +mcdc.setting(N_particle=80) # Run mcdc.run() diff --git a/examples/fixed_source/kobayashi3-TD/process.py b/examples/fixed_source/kobayashi3-TD/process.py index 45334989..6e37fc1c 100644 --- a/examples/fixed_source/kobayashi3-TD/process.py +++ b/examples/fixed_source/kobayashi3-TD/process.py @@ -25,6 +25,12 @@ phi = flux["mean"][:] phi_sd = flux["sdev"][:] + for i in range(len(f["input_deck"]["cell_tallies"])): + flux_score = f[f"tallies/cell_tally_{i}/flux"] + print( + f'cell {i+1} mean = {flux_score["mean"][()]}, sdev = {flux_score["sdev"][()]}' + ) + fig, ax = plt.subplots() cax = ax.pcolormesh(X, Y, phi[0], vmin=phi[0].min(), vmax=phi[0].max()) text = ax.text(0.02, 1.02, "", transform=ax.transAxes) diff --git a/examples/fixed_source/sphere_in_cube/input.py b/examples/fixed_source/sphere_in_cube/input.py new file mode 100644 index 00000000..61157306 --- /dev/null +++ b/examples/fixed_source/sphere_in_cube/input.py @@ -0,0 +1,57 @@ +import numpy as np +import mcdc + +# ============================================================================= +# Set model +# ============================================================================= +# Homogeneous pure-fission sphere inside a pure-scattering cube + +# Set materials +pure_f = mcdc.material(fission=np.array([1.0]), nu_p=np.array([1.2])) +pure_s = mcdc.material(scatter=np.array([[1.0]])) + +# Set surfaces +sx1 = mcdc.surface("plane-x", x=0.0, bc="vacuum") +sx2 = mcdc.surface("plane-x", x=4.0, bc="vacuum") +sy1 = mcdc.surface("plane-y", y=0.0, bc="vacuum") +sy2 = mcdc.surface("plane-y", y=4.0, bc="vacuum") +sz1 = mcdc.surface("plane-z", z=0.0, bc="vacuum") +sz2 = mcdc.surface("plane-z", z=4.0, bc="vacuum") +sphere = mcdc.surface("sphere", center=[2.0, 2.0, 2.0], radius=1.5) +inside_sphere = -sphere +inside_box = +sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2 + +# Set cells +# Source +mcdc.cell(inside_box & ~inside_sphere, pure_s) + +# Sphere +sphere_cell = mcdc.cell(inside_sphere, pure_f) + +# ============================================================================= +# Set source +# ============================================================================= +# The source pulses in t=[0,5] + +mcdc.source(x=[0.0, 4.0], y=[0.0, 4.0], z=[0.0, 4.0], time=[0.0, 50.0], isotropic=True) + +# ============================================================================= +# Set tally, setting, and run mcdc +# ============================================================================= +mcdc.tally.mesh_tally( + scores=["fission"], + x=np.linspace(0.0, 4.0, 2), + y=np.linspace(0.0, 4.0, 2), + z=np.linspace(0.0, 4.0, 2), + # t=np.linspace(0.0, 200.0, 2), +) + + +mcdc.tally.cell_tally(sphere_cell, scores=["fission"]) + +# Setting +mcdc.setting(N_particle=100) +mcdc.implicit_capture() + +# Run +mcdc.run() diff --git a/examples/fixed_source/sphere_in_cube/process.py b/examples/fixed_source/sphere_in_cube/process.py new file mode 100644 index 00000000..3cc7373f --- /dev/null +++ b/examples/fixed_source/sphere_in_cube/process.py @@ -0,0 +1,23 @@ +import h5py + +# Load results +with h5py.File("output.h5", "r") as f: + # print(f["tallies"].keys()) + print(f["input_deck"]["cell_tallies"].keys()) + + for i in range(len(f["input_deck"]["cell_tallies"])): + fission_score = f[f"tallies/cell_tally_{i}/fission"] + + print( + f'for sphere {i+1}, mean = {fission_score["mean"][()]}, sdev = {fission_score["sdev"][()]}' + ) + + # print(fission_score["mean"][()]) + # print(fission_score["sdev"][()]) + + # print(f"fission_score mean = {fission_score["mean"][()]}") + # print(f"fission_score mean = {fission_score["sdev"][()]}") + + # cell = f["tallies/cell_tally_0/fission"] + # print(f'sphere1 mean = {cell["mean"][()]}') + # print(f'sphere2 sdev = {cell["sdev"][()]}') diff --git a/mcdc/card.py b/mcdc/card.py index 6920faf3..0b300c16 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -211,6 +211,8 @@ def __init__(self): self.fill_ID = None self.translation = np.array([0.0, 0.0, 0.0]) self.surface_IDs = np.zeros(0, dtype=int) + self.N_tally = 0 + self.tally_IDs = [] self._region_RPN = [] # Reverse Polish Notation def set_region_RPN(self): @@ -383,3 +385,12 @@ def __init__(self, surface_ID): # Set card data self.surface_ID = surface_ID self.N_bin = 1 + + +class CellTallyCard(TallyCard): + def __init__(self, cell_ID): + TallyCard.__init__(self, "Cell tally") + + # Set card data + self.cell_ID = cell_ID + self.N_bin = 1 diff --git a/mcdc/constant.py b/mcdc/constant.py index a584fe5e..0e7d2fe5 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -113,9 +113,9 @@ NEUTRON_MASS = 1.67492749804e-27 # kg EV_TO_J = 1.6022e-19 SQRT_E_TO_SPEED = math.sqrt(2.0 * EV_TO_J / NEUTRON_MASS) * 100 -BOLTZMAN_K = 8.61733326e-5 # eV/K +BOLTZMANN_K = 8.61733326e-5 # eV/K T_ROOM = 294 # K -E_THERMAL_THRESHOLD = 400 * BOLTZMAN_K * T_ROOM +E_THERMAL_THRESHOLD = 400 * BOLTZMANN_K * T_ROOM # Cross Section Type XS_TOTAL = 0 diff --git a/mcdc/global_.py b/mcdc/global_.py index 2c1efe1e..3e3d8f58 100644 --- a/mcdc/global_.py +++ b/mcdc/global_.py @@ -35,6 +35,7 @@ def reset(self): self.sources = [] self.mesh_tallies = [] self.surface_tallies = [] + self.cell_tallies = [] self.setting = { "tag": "Setting", diff --git a/mcdc/kernel.py b/mcdc/kernel.py index d3921224..032214d2 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -1980,6 +1980,31 @@ def score_surface_tally(P_arr, surface, tally, data, mcdc): adapt.global_add(tally_bin, (TALLY_SCORE, idx + i), round(score)) +@njit +def score_cell_tally(P_arr, distance, tally, data, mcdc): + P = P_arr[0] + tally_bin = data[TALLY] + material = mcdc["materials"][P["material_ID"]] + stride = tally["stride"] + cell_idx = stride["tally"] + score = 0 + + # Score + flux = distance * P["w"] + for i in range(tally["N_score"]): + score_type = tally["scores"][i] + if score_type == SCORE_FLUX: + score = flux + elif score_type == SCORE_TOTAL: + SigmaT = get_MacroXS(XS_TOTAL, material, P_arr, mcdc) + score = flux * SigmaT + elif score_type == SCORE_FISSION: + SigmaF = get_MacroXS(XS_FISSION, material, P_arr, mcdc) + score = flux * SigmaF + + tally_bin[TALLY_SCORE, cell_idx] += score + + @njit def tally_reduce(data, mcdc): tally_bin = data[TALLY] @@ -2000,12 +2025,20 @@ def tally_reduce(data, mcdc): @njit def tally_accumulate(data, mcdc): + + # print(f'fixed_source_data = {data[TALLY][TALLY_SCORE][-1]}') + tally_bin = data[TALLY] + # print(data[0][0, 30000], data[0][0, 30001]) N_bin = tally_bin.shape[1] + # print(f'N_bin = {N_bin}') + for i in range(N_bin): # Accumulate score and square of score into sum and sum_sq score = tally_bin[TALLY_SCORE, i] + # if (score != 0 and i > 29999): + # print(f'score = {score}, i = {i}') tally_bin[TALLY_SUM, i] += score tally_bin[TALLY_SUM_SQ, i] += score * score @@ -2362,8 +2395,16 @@ def move_to_event(P_arr, data, mcdc): # Score tracklength tallies if mcdc["cycle_active"]: + # Mesh tallies for tally in mcdc["mesh_tallies"]: score_mesh_tally(P_arr, distance, tally, data, mcdc) + + # Cell tallies + cell = mcdc["cells"][P["cell_ID"]] + for i in range(cell["N_tally"]): + ID = cell["tally_IDs"][i] + tally = mcdc["cell_tallies"][ID] + score_cell_tally(P_arr, distance, tally, data, mcdc) if mcdc["setting"]["mode_eigenvalue"]: eigenvalue_tally(P_arr, distance, mcdc) @@ -2403,6 +2444,7 @@ def surface_crossing(P_arr, data, prog): geometry.surface_bc(P_arr, surface) # Score tally + # N_tally is an int64, tally_IDs is a numpy.ndarray for i in range(surface["N_tally"]): ID = surface["tally_IDs"][i] tally = mcdc["surface_tallies"][ID] diff --git a/mcdc/main.py b/mcdc/main.py index 5f587d5d..414a9ad5 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -352,11 +352,13 @@ def prepare(): type_.make_type_nuclide(input_deck) type_.make_type_material(input_deck) type_.make_type_surface(input_deck) + type_.make_type_cell(input_deck) type_.make_type_universe(input_deck) type_.make_type_lattice(input_deck) type_.make_type_source(input_deck) type_.make_type_mesh_tally(input_deck) type_.make_type_surface_tally(input_deck) + type_.make_type_cell_tally(input_deck) type_.make_type_setting(input_deck) type_.make_type_uq(input_deck) type_.make_type_domain_decomp(input_deck) @@ -506,7 +508,7 @@ def prepare(): surface_data_idx = 0 region_data_idx = 0 for i in range(N_cell): - for name in ["ID", "fill_ID", "translation"]: + for name in ["ID", "fill_ID", "translation", "N_tally"]: copy_field(mcdc["cells"][i], input_deck.cells[i], name) # Fill type @@ -539,6 +541,11 @@ def prepare(): ) region_data_idx += N_RPN + 1 + # Variables with possible different sizes + for name in ["tally_IDs"]: + N = len(getattr(input_deck.cells[i], name)) + mcdc["cells"][i][name][:N] = getattr(input_deck.cells[i], name) + # ========================================================================= # Universes # ========================================================================= @@ -599,6 +606,7 @@ def prepare(): N_mesh_tally = len(input_deck.mesh_tallies) N_surface_tally = len(input_deck.surface_tallies) + N_cell_tally = len(input_deck.cell_tallies) tally_size = 0 # Mesh tallies @@ -746,6 +754,70 @@ def prepare(): mcdc["surface_tallies"][i]["stride"]["tally"] = tally_size tally_size += mcdc["surface_tallies"][i]["N_bin"] + # Cell tallies + for i in range(N_cell_tally): + copy_field(mcdc["cell_tallies"][i], input_deck.cell_tallies[i], "N_bin") + + # Filters (variables with possible different sizes) + for name in ["t", "mu", "azi", "g"]: + N = len(getattr(input_deck.cell_tallies[i], name)) + mcdc["cell_tallies"][i]["filter"][name][:N] = getattr( + input_deck.cell_tallies[i], name + ) + + # Differentiating the tallies by cell_ID + mcdc["cell_tallies"][i]["filter"]["cell_ID"] = getattr( + input_deck.cell_tallies[i], "cell_ID" + ) + + # Set tally scores and their strides + N_score = len(input_deck.cell_tallies[i].scores) + mcdc["cell_tallies"][i]["N_score"] = N_score + for j in range(N_score): + score_name = input_deck.cell_tallies[i].scores[j] + score_type = None + if score_name == "flux": + score_type = SCORE_FLUX + elif score_name == "fission": + score_type = SCORE_FISSION + elif score_name == "net-current": + score_type = SCORE_NET_CURRENT + mcdc["cell_tallies"][i]["scores"][j] = score_type + + # Filter grid sizes + Nmu = len(input_deck.cell_tallies[i].mu) - 1 + N_azi = len(input_deck.cell_tallies[i].azi) - 1 + Ng = len(input_deck.cell_tallies[i].g) - 1 + Nt = len(input_deck.cell_tallies[i].t) - 1 + + # Update N_bin + mcdc["cell_tallies"][i]["N_bin"] *= N_score + + # Filter strides + stride = N_score + if Nt > 1: + mcdc["mesh_tallies"][i]["stride"]["t"] = stride + stride *= Nt + if Ng > 1: + mcdc["mesh_tallies"][i]["stride"]["g"] = stride + stride *= Ng + if N_azi > 1: + mcdc["mesh_tallies"][i]["stride"]["azi"] = stride + stride *= N_azi + if Nmu > 1: + mcdc["mesh_tallies"][i]["stride"]["mu"] = stride + stride *= Nmu + + # Set tally stride and accumulate total tally size + mcdc["cell_tallies"][i]["stride"]["tally"] = tally_size + tally_size += mcdc["cell_tallies"][i]["N_bin"] + + # Set tally data + if not input_deck.technique["uq"]: + tally = np.zeros((3, tally_size), dtype=type_.float64) + else: + tally = np.zeros((5, tally_size), dtype=type_.float64) + # ========================================================================= # Establish Data Type from Tally Info and Construct Tallies # ========================================================================= @@ -786,11 +858,18 @@ def prepare(): t_limit = max( [ tally["filter"]["t"][-1] - for tally in list(mcdc["mesh_tallies"]) + list(mcdc["surface_tallies"]) + for tally in list(mcdc["mesh_tallies"]) + + list(mcdc["surface_tallies"]) + + list(mcdc["cell_tallies"]) ] ) - if len(input_deck.mesh_tallies) + len(input_deck.surface_tallies) == 0: + if ( + len(input_deck.mesh_tallies) + + len(input_deck.surface_tallies) + + len(input_deck.cell_tallies) + == 0 + ): t_limit = INF # Check if time boundary is above the final tally mesh time grid @@ -1246,10 +1325,15 @@ def generate_hdf5(data, mcdc): cardlist_to_h5group(input_deck.universes, input_group, "universe") cardlist_to_h5group(input_deck.lattices, input_group, "lattice") cardlist_to_h5group(input_deck.sources, input_group, "source") - cardlist_to_h5group(input_deck.mesh_tallies, input_group, "mesh_tallie") + cardlist_to_h5group( + input_deck.mesh_tallies, input_group, "mesh_tallies" + ) cardlist_to_h5group( input_deck.surface_tallies, input_group, "surface_tally" ) + cardlist_to_h5group( + input_deck.cell_tallies, input_group, "cell_tallies" + ) dict_to_h5group(input_deck.setting, input_group.create_group("setting")) dict_to_h5group( input_deck.technique, input_group.create_group("technique") @@ -1373,6 +1457,61 @@ def generate_hdf5(data, mcdc): uq_var = tot_var - mc_var f.create_dataset(group_name + "uq_var", data=uq_var) + # Cell tallies + for ID, tally in enumerate(mcdc["cell_tallies"]): + if mcdc["technique"]["iQMC"]: + break + + # Shape + N_score = tally["N_score"] + + if not mcdc["technique"]["uq"]: + shape = (3, N_score) + else: + shape = (5, N_score) + + # Reshape tally + N_bin = tally["N_bin"] + start = tally["stride"]["tally"] + tally_bin = data[TALLY][:, start : start + N_bin] + + # print(f'data = {data[TALLY][2]}') + # if data[TALLY][1].all() == data[TALLY][2].all: + # print('hell yeah') + # print(f'data keys = {data[TALLY].dtype.names}') + + # print(f'tally_bin = {tally_bin}') + + tally_bin = tally_bin.reshape(shape) + + # Roll tally so that score is in the front + tally_bin = np.rollaxis(tally_bin, 1, 0) + + # Iterate over scores + for i in range(N_score): + score_type = tally["scores"][i] + score_tally_bin = np.squeeze(tally_bin[i]) + # print(f'score_tally_bin = {score_tally_bin}') + if score_type == SCORE_FLUX: + score_name = "flux" + elif score_type == SCORE_NET_CURRENT: + score_name = "net-current" + elif score_type == SCORE_FISSION: + score_name = "fission" + group_name = "tallies/cell_tally_%i/%s/" % (ID, score_name) + + mean = score_tally_bin[TALLY_SUM] + # print(f'ID = {ID}, score_name = {score_name}, mean = {mean}') + sdev = score_tally_bin[TALLY_SUM_SQ] + + f.create_dataset(group_name + "mean", data=mean) + f.create_dataset(group_name + "sdev", data=sdev) + if mcdc["technique"]["uq"]: + mc_var = score_tally_bin[TALLY_UQ_BATCH_VAR] + tot_var = score_tally_bin[TALLY_UQ_BATCH] + uq_var = tot_var - mc_var + f.create_dataset(group_name + "uq_var", data=uq_var) + # Eigenvalues if mcdc["setting"]["mode_eigenvalue"]: if mcdc["technique"]["iQMC"]: diff --git a/mcdc/tally.py b/mcdc/tally.py index 0630e71b..f798af84 100644 --- a/mcdc/tally.py +++ b/mcdc/tally.py @@ -8,6 +8,7 @@ from mcdc.card import ( MeshTallyCard, SurfaceTallyCard, + CellTallyCard, ) from mcdc.constant import ( INF, @@ -157,3 +158,34 @@ def surface_tally( global_.input_deck.surface_tallies.append(card) return card + + +def cell_tally(cell, scores=["flux"]): + + # Make tally card + card = CellTallyCard(cell.ID) + + # Set ID + card.ID = len(global_.input_deck.cell_tallies) + + # Set cell + card.cell_ID = cell.ID + cell.tally_IDs.append(card.ID) + cell.N_tally += 1 + + # Calculate total number bins + card.N_bin = 1 + + # Scores + for s in scores: + score_checked = check_support( + "score type", + s, + ["flux", "net-current", "fission"], + ) + card.scores.append(score_checked) + + # Add to deck + global_.input_deck.cell_tallies.append(card) + + return card diff --git a/mcdc/type_.py b/mcdc/type_.py index 27189baf..20606f23 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -44,6 +44,7 @@ setting = None mesh_tally = None surface_tally = None +cell_tally = None technique = None global_ = None @@ -551,20 +552,30 @@ def make_type_surface(input_deck): # ============================================================================== -cell = into_dtype( - [ - ("ID", int64), - # Fill status - ("fill_type", int64), - ("fill_ID", int64), - ("fill_translated", bool), - # Local coordinate modifier - ("translation", float64, (3,)), - # Data indices - ("surface_data_idx", int64), - ("region_data_idx", int64), - ] -) +def make_type_cell(input_deck): + global cell + + # Maximum number tallies + Nmax_tally = 0 + for cell in input_deck.cells: + Nmax_tally = max(Nmax_tally, len(cell.tally_IDs)) + + cell = into_dtype( + [ + ("ID", int64), + # Fill status + ("fill_type", int64), + ("fill_ID", int64), + ("fill_translated", bool), + # Local coordinate modifier + ("translation", float64, (3,)), + # Data indices + ("surface_data_idx", int64), + ("region_data_idx", int64), + ("N_tally", int64), + ("tally_IDs", int64, (Nmax_tally,)), + ] + ) # ============================================================================== @@ -794,7 +805,9 @@ def make_type_surface_tally(input_deck): Nmax_azi = 2 Nmax_g = 2 Nmax_score = 1 - for card in input_deck.mesh_tallies: + + # IDK if this is right, but I changed this to input_deck.surface_tallies + for card in input_deck.surface_tallies: Nmax_t = max(Nmax_t, len(card.t)) Nmax_mu = max(Nmax_mu, len(card.mu)) Nmax_azi = max(Nmax_azi, len(card.azi)) @@ -832,6 +845,55 @@ def make_type_surface_tally(input_deck): surface_tally = into_dtype(struct) +def make_type_cell_tally(input_deck): + global cell_tally + struct = [] + + # Maximum number of grid for each mesh coordinate and filter + Nmax_t = 2 + Nmax_mu = 2 + Nmax_azi = 2 + Nmax_g = 2 + Nmax_score = 1 + + for card in input_deck.cell_tallies: + Nmax_t = max(Nmax_t, len(card.t)) + Nmax_mu = max(Nmax_mu, len(card.mu)) + Nmax_azi = max(Nmax_azi, len(card.azi)) + Nmax_g = max(Nmax_g, len(card.g)) + Nmax_score = max(Nmax_score, len(card.scores)) + + # Set the filter + filter_ = [ + ("cell_ID", int64), + ("t", float64, (Nmax_t,)), + ("mu", float64, (Nmax_mu,)), + ("azi", float64, (Nmax_azi,)), + ("g", float64, (Nmax_g,)), + ] + struct = [("filter", filter_)] + + # Tally strides + stride = [ + ("tally", int64), + ("sensitivity", int64), + ("mu", int64), + ("azi", int64), + ("g", int64), + ("t", int64), + ] + struct += [("stride", stride)] + + # Total number of bins + struct += [("N_bin", int64)] + + # Scores + struct += [("N_score", int64), ("scores", int64, (Nmax_score,))] + + # Make tally structure + cell_tally = into_dtype(struct) + + # ============================================================================== # Setting # ============================================================================== @@ -1255,6 +1317,7 @@ def make_type_global(input_deck): N_lattice = len(input_deck.lattices) N_mesh_tally = len(input_deck.mesh_tallies) N_surface_tally = len(input_deck.surface_tallies) + N_cell_tally = len(input_deck.cell_tallies) # Cell data sizes N_cell_surface = 0 @@ -1325,6 +1388,7 @@ def make_type_global(input_deck): ("sources", source, (N_source,)), ("mesh_tallies", mesh_tally, (N_mesh_tally,)), ("surface_tallies", surface_tally, (N_surface_tally,)), + ("cell_tallies", cell_tally, (N_cell_tally,)), ("setting", setting), ("technique", technique), ("domain_decomp", domain_decomp), diff --git a/test/regression/kobayashi3-TD/answer.h5 b/test/regression/kobayashi3-TD/answer.h5 index ad289a0f..bd56f20c 100644 Binary files a/test/regression/kobayashi3-TD/answer.h5 and b/test/regression/kobayashi3-TD/answer.h5 differ diff --git a/test/regression/kobayashi3-TD/input.py b/test/regression/kobayashi3-TD/input.py index 607f7630..4bdea735 100644 --- a/test/regression/kobayashi3-TD/input.py +++ b/test/regression/kobayashi3-TD/input.py @@ -29,22 +29,18 @@ sz5 = mcdc.surface("plane-z", z=60.0, bc="vacuum") # Set cells -# Soruce -mcdc.cell(+sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2, m) +# Source +source_cell = mcdc.cell(+sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2, m) # Voids -mcdc.cell(+sx1 & -sx2 & +sy2 & -sy3 & +sz1 & -sz2, m_void) -mcdc.cell(+sx1 & -sx3 & +sy3 & -sy4 & +sz1 & -sz2, m_void) -mcdc.cell(+sx3 & -sx4 & +sy3 & -sy4 & +sz1 & -sz3, m_void) -mcdc.cell(+sx3 & -sx4 & +sy3 & -sy5 & +sz3 & -sz4, m_void) +channel_1 = +sx1 & -sx2 & +sy2 & -sy3 & +sz1 & -sz2 +channel_2 = +sx1 & -sx3 & +sy3 & -sy4 & +sz1 & -sz2 +channel_3 = +sx3 & -sx4 & +sy3 & -sy4 & +sz1 & -sz3 +channel_4 = +sx3 & -sx4 & +sy3 & -sy5 & +sz3 & -sz4 +void_channel = channel_1 | channel_2 | channel_3 | channel_4 +void_cell = mcdc.cell(void_channel, m_void) # Shield -mcdc.cell(+sx1 & -sx3 & +sy1 & -sy5 & +sz2 & -sz5, m) -mcdc.cell(+sx2 & -sx5 & +sy1 & -sy3 & +sz1 & -sz2, m) -mcdc.cell(+sx3 & -sx5 & +sy1 & -sy3 & +sz2 & -sz5, m) -mcdc.cell(+sx3 & -sx5 & +sy4 & -sy5 & +sz1 & -sz3, m) -mcdc.cell(+sx4 & -sx5 & +sy4 & -sy5 & +sz3 & -sz5, m) -mcdc.cell(+sx4 & -sx5 & +sy3 & -sy4 & +sz1 & -sz5, m) -mcdc.cell(+sx3 & -sx4 & +sy3 & -sy5 & +sz4 & -sz5, m) -mcdc.cell(+sx1 & -sx3 & +sy4 & -sy5 & +sz1 & -sz2, m) +box = +sx1 & -sx5 & +sy1 & -sy5 & +sz1 & -sz5 +shield_cell = mcdc.cell(box & ~void_channel, m) # ============================================================================= # Set source @@ -62,14 +58,17 @@ # Tally: z-integrated flux (X-Y section view) mcdc.tally.mesh_tally( scores=["flux"], - x=np.linspace(0.0, 60.0, 61), - y=np.linspace(0.0, 100.0, 101), + x=np.linspace(0.0, 60.0, 31), + y=np.linspace(0.0, 100.0, 51), t=np.linspace(0.0, 200.0, 21), ) +mcdc.tally.cell_tally(source_cell, scores=["flux"]) +mcdc.tally.cell_tally(void_cell, scores=["flux"]) +mcdc.tally.cell_tally(shield_cell, scores=["flux"]) + # Setting mcdc.setting(N_particle=80) -mcdc.implicit_capture() # Run mcdc.run() diff --git a/test/regression/sphere_in_cube/answer.h5 b/test/regression/sphere_in_cube/answer.h5 new file mode 100644 index 00000000..b48e38da Binary files /dev/null and b/test/regression/sphere_in_cube/answer.h5 differ diff --git a/test/regression/sphere_in_cube/input.py b/test/regression/sphere_in_cube/input.py new file mode 100644 index 00000000..61157306 --- /dev/null +++ b/test/regression/sphere_in_cube/input.py @@ -0,0 +1,57 @@ +import numpy as np +import mcdc + +# ============================================================================= +# Set model +# ============================================================================= +# Homogeneous pure-fission sphere inside a pure-scattering cube + +# Set materials +pure_f = mcdc.material(fission=np.array([1.0]), nu_p=np.array([1.2])) +pure_s = mcdc.material(scatter=np.array([[1.0]])) + +# Set surfaces +sx1 = mcdc.surface("plane-x", x=0.0, bc="vacuum") +sx2 = mcdc.surface("plane-x", x=4.0, bc="vacuum") +sy1 = mcdc.surface("plane-y", y=0.0, bc="vacuum") +sy2 = mcdc.surface("plane-y", y=4.0, bc="vacuum") +sz1 = mcdc.surface("plane-z", z=0.0, bc="vacuum") +sz2 = mcdc.surface("plane-z", z=4.0, bc="vacuum") +sphere = mcdc.surface("sphere", center=[2.0, 2.0, 2.0], radius=1.5) +inside_sphere = -sphere +inside_box = +sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2 + +# Set cells +# Source +mcdc.cell(inside_box & ~inside_sphere, pure_s) + +# Sphere +sphere_cell = mcdc.cell(inside_sphere, pure_f) + +# ============================================================================= +# Set source +# ============================================================================= +# The source pulses in t=[0,5] + +mcdc.source(x=[0.0, 4.0], y=[0.0, 4.0], z=[0.0, 4.0], time=[0.0, 50.0], isotropic=True) + +# ============================================================================= +# Set tally, setting, and run mcdc +# ============================================================================= +mcdc.tally.mesh_tally( + scores=["fission"], + x=np.linspace(0.0, 4.0, 2), + y=np.linspace(0.0, 4.0, 2), + z=np.linspace(0.0, 4.0, 2), + # t=np.linspace(0.0, 200.0, 2), +) + + +mcdc.tally.cell_tally(sphere_cell, scores=["fission"]) + +# Setting +mcdc.setting(N_particle=100) +mcdc.implicit_capture() + +# Run +mcdc.run()