fix mesh tally

ilhamv · Dec 4, 2024 · 3bc8f4e · 3bc8f4e
1 parent 9865bb7
commit 3bc8f4e
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 23 deletions.
diff --git a/mcdc/kernel.py b/mcdc/kernel.py
@@ -1781,12 +1781,13 @@ def mesh_get_energy_index(P_arr, mesh, mode_MG):
     # Check if outside grid
     outside = False
 
+    Ng = mesh["Ng"]
+
     if mode_MG:
-        return binary_search(P["g"], mesh["g"]), outside
+        return binary_search_with_length(P["g"], mesh["g"], Ng), outside
 
     else:
         E = P["E"]
-        Ng = mesh["Ng"]
         if E < mesh["g"][0] or E > mesh["g"][Ng]:
             outside = True
             return 0, outside
@@ -1921,7 +1922,7 @@ def score_mesh_tally(P_arr, distance, tally, data, mcdc):
         if mesh_crossed == MESH_X:
             if ux > 0.0:
                 ix += 1
-                if ix == len(mesh["x"]) - 1:
+                if ix == mesh["Nx"]:
                     break
                 idx += stride["x"]
             else:
@@ -1932,7 +1933,7 @@ def score_mesh_tally(P_arr, distance, tally, data, mcdc):
         elif mesh_crossed == MESH_Y:
             if uy > 0.0:
                 iy += 1
-                if iy == len(mesh["y"]) - 1:
+                if iy == mesh["Ny"]:
                     break
                 idx += stride["y"]
             else:
@@ -1943,7 +1944,7 @@ def score_mesh_tally(P_arr, distance, tally, data, mcdc):
         elif mesh_crossed == MESH_Z:
             if uz > 0.0:
                 iz += 1
-                if iz == len(mesh["z"]) - 1:
+                if iz == mesh["Nz"]:
                     break
                 idx += stride["z"]
             else:
@@ -1953,7 +1954,7 @@ def score_mesh_tally(P_arr, distance, tally, data, mcdc):
                 idx -= stride["z"]
         elif mesh_crossed == MESH_T:
             it += 1
-            if it == len(mesh["t"]) - 1:
+            if it == mesh["Nt"]:
                 break
             idx += stride["t"]
 

diff --git a/mcdc/main.py b/mcdc/main.py
@@ -1460,25 +1460,38 @@ def generate_hdf5(data, mcdc):
                     break
 
                 mesh = tally["filter"]
-                f.create_dataset("tallies/mesh_tally_%i/grid/t" % ID, data=mesh["t"])
-                f.create_dataset("tallies/mesh_tally_%i/grid/x" % ID, data=mesh["x"])
-                f.create_dataset("tallies/mesh_tally_%i/grid/y" % ID, data=mesh["y"])
-                f.create_dataset("tallies/mesh_tally_%i/grid/z" % ID, data=mesh["z"])
-                f.create_dataset("tallies/mesh_tally_%i/grid/mu" % ID, data=mesh["mu"])
-                f.create_dataset(
-                    "tallies/mesh_tally_%i/grid/azi" % ID, data=mesh["azi"]
-                )
-                f.create_dataset("tallies/mesh_tally_%i/grid/g" % ID, data=mesh["g"])
 
-                # Shape
-                Nmu = tally["filter"]["Nmu"]
-                N_azi = tally["filter"]["N_azi"]
-                Ng = tally["filter"]["Ng"]
-                Nx = tally["filter"]["Nx"]
-                Ny = tally["filter"]["Ny"]
-                Nz = tally["filter"]["Nz"]
-                Nt = tally["filter"]["Nt"]
+                # Filter
+                Nmu = mesh["Nmu"]
+                N_azi = mesh["N_azi"]
+                Ng = mesh["Ng"]
+                Nx = mesh["Nx"]
+                Ny = mesh["Ny"]
+                Nz = mesh["Nz"]
+                Nt = mesh["Nt"]
                 N_score = tally["N_score"]
+                #
+                f.create_dataset(
+                    "tallies/mesh_tally_%i/grid/t" % ID, data=mesh["t"][: Nt + 1]
+                )
+                f.create_dataset(
+                    "tallies/mesh_tally_%i/grid/x" % ID, data=mesh["x"][: Nx + 1]
+                )
+                f.create_dataset(
+                    "tallies/mesh_tally_%i/grid/y" % ID, data=mesh["y"][:Ny] + 1
+                )
+                f.create_dataset(
+                    "tallies/mesh_tally_%i/grid/z" % ID, data=mesh["z"][: Nz + 1]
+                )
+                f.create_dataset(
+                    "tallies/mesh_tally_%i/grid/mu" % ID, data=mesh["mu"][: Nmu + 1]
+                )
+                f.create_dataset(
+                    "tallies/mesh_tally_%i/grid/azi" % ID, data=mesh["azi"][: N_azi + 1]
+                )
+                f.create_dataset(
+                    "tallies/mesh_tally_%i/grid/g" % ID, data=mesh["g"][: Ng + 1]
+                )
 
                 if mcdc["technique"]["domain_decomposition"]:
                     Nx *= input_deck.technique["dd_mesh"]["x"].size - 1