WW

mcdc/global_.py

@@ -82,6 +82,7 @@ def reset(self):
             "ww": np.ones([1, 1, 1, 1]),
             "ww_width": 2.5,
             "ww_mesh": make_card_mesh(),
+            "ww_auto": False,
             "domain_decomposition": False,
             "dd_idx": 0,
             "dd_mesh": make_card_mesh(),

mcdc/input_.py

@@ -1217,7 +1217,7 @@ def time_census(t):
     """
 
     # Remove census beyond the final tally time grid point
-    while True:
+    while False:
         if t[-1] >= global_.input_deck.tally["mesh"]["t"][-1]:
             t = t[:-1]
         else:
@@ -1274,6 +1274,19 @@ def weight_window(x=None, y=None, z=None, t=None, window=None, width=None):
         card["ww_mesh"]["t"] = t
 
     # Set window
+    if window is None:
+        card["ww_auto"] = True
+        window_ax = []
+        if t is not None:
+            window_ax.append(len(t) - 1)
+        if x is not None:
+            window_ax.append(len(x) - 1)
+        if y is not None:
+            window_ax.append(len(y) - 1)
+        if z is not None:
+            window_ax.append(len(z) - 1)
+        window = np.ones(window_ax)
+
     ax_expand = []
     if t is None:
         ax_expand.append(0)

mcdc/kernel.py

@@ -3135,7 +3135,7 @@ def weight_window(P_arr, prog):
         P["w"] = w_target
 
         # Splitting (keep the original particle)
-        n_split = math.floor(p)
+        n_split = min(math.floor(p), 5)
         for i in range(n_split - 1):
             split_as_record(P_new_arr, P_arr)
             adapt.add_active(P_new_arr, prog)
@@ -3157,6 +3157,84 @@ def weight_window(P_arr, prog):
             P["w"] = w_target
 
 
+# =============================================================================
+# Weight widow
+# =============================================================================
+
+
+def get_flux(idx, data, mcdc):
+    # Determine the correct tally list based on tally_type
+    tallies = mcdc["mesh_tallies"]
+
+    # Loop over the tallies to find the requested score
+    for ID, tally in enumerate(tallies):
+        if mcdc["technique"]["iQMC"]:
+            break
+
+        mesh = tally["filter"]
+
+        # Shape based on technique and tally type
+        N_sensitivity = 0
+        Ns = 1 + N_sensitivity
+
+        Nmu = len(mesh["mu"]) - 1
+        N_azi = len(mesh["azi"]) - 1
+        Ng = len(mesh["g"]) - 1
+        Nt = len(mesh["t"]) - 1
+
+        Nx = len(mesh["x"]) - 1
+        Ny = len(mesh["y"]) - 1
+        Nz = len(mesh["z"]) - 1
+
+        N_score = tally["N_score"]
+
+        if not mcdc["technique"]["uq"]:
+            shape = (3, Ns, Nmu, N_azi, Ng, Nt, Nx, Ny, Nz, N_score)
+        else:
+            shape = (5, Ns, Nmu, N_azi, Ng, Nt, Nx, Ny, Nz, N_score)
+
+        # Reshape tally
+        N_bin = tally["N_bin"]
+        start = tally["stride"]["tally"]
+        tally_bin = np.ascontiguousarray(data[TALLY][:, start : start + N_bin])
+        tally_bin = tally_bin.reshape(shape)
+
+        # Roll tally so that score is in the front
+        tally_bin = np.rollaxis(tally_bin, 9, 0)
+
+        # Iterate over scores to find the requested score
+        for i in range(N_score):
+            score_type = tally["scores"][i]
+            if score_type == SCORE_FLUX:
+                score_tally_bin = np.squeeze(tally_bin[i])
+                mean = score_tally_bin[TALLY_SUM]
+                sdev = score_tally_bin[TALLY_SUM_SQ]
+                old_flux = mean[idx][:]
+                if Nx == 1 and old_flux.ndim < 3:
+                    old_flux = np.expand_dims(old_flux, axis=0)
+                if Ny == 1 and old_flux.ndim < 3:
+                    old_flux = np.expand_dims(old_flux, axis=1)
+                if Nz == 1 and old_flux.ndim < 3:
+                    old_flux = np.expand_dims(old_flux, axis=2)
+                return old_flux
+
+
+@njit
+def update_weight_window(census_timestep, data, mcdc):
+    # Target weight
+    if census_timestep > 0:
+        window_centers = mcdc["technique"]["ww"]
+        with objmode(old_flux="float64[:,:,:]"):
+            old_flux = get_flux(census_timestep, data, mcdc)
+            if np.max(old_flux) > 0:
+                min_flux = np.min(old_flux[old_flux != 0])
+                old_flux[old_flux == 0] = min_flux / 2
+            else:
+                old_flux = np.ones_like(old_flux)
+        window_centers[census_timestep] = old_flux
+        window_centers[census_timestep] /= np.max(window_centers[census_timestep])
+
+
 # =============================================================================
 # Weight Roulette
 # =============================================================================

mcdc/loop.py

@@ -89,6 +89,8 @@ def loop_fixed_source(data_arr, mcdc_arr):
 
         # Loop over time censuses
         for idx_census in range(mcdc["setting"]["N_census"]):
+            if mcdc["technique"]["ww_auto"]:
+                kernel.update_weight_window(idx_census, data, mcdc)
             mcdc["idx_census"] = idx_census
             seed_census = kernel.split_seed(seed_batch, SEED_SPLIT_CENSUS)
 

mcdc/main.py

@@ -945,6 +945,7 @@ def prepare():
 
     # WW windows
     mcdc["technique"]["ww"] = input_deck.technique["ww"]
+    mcdc["technique"]["ww_auto"] = input_deck.technique["ww_auto"]
     mcdc["technique"]["ww_width"] = input_deck.technique["ww_width"]
 
     # =========================================================================
@@ -1339,6 +1340,15 @@ def generate_hdf5(data, mcdc):
                     input_deck.technique, input_group.create_group("technique")
                 )
 
+            if mcdc["technique"]["weight_window"]:
+                T = mcdc["technique"]
+                f.create_dataset("ww/grid/t", data=T["ww_mesh"]["t"])
+                f.create_dataset("ww/grid/x", data=T["ww_mesh"]["x"])
+                f.create_dataset("ww/grid/y", data=T["ww_mesh"]["y"])
+                f.create_dataset("ww/grid/z", data=T["ww_mesh"]["z"])
+                f.create_dataset("ww/center", data=np.squeeze(T["ww"]))
+                f.create_dataset("ww/width", data=T["ww_width"])
+
             # Mesh tallies
             for ID, tally in enumerate(mcdc["mesh_tallies"]):
                 if mcdc["technique"]["iQMC"]:

mcdc/type_.py

@@ -1023,6 +1023,7 @@ def make_type_technique(input_deck):
     mesh, Nx, Ny, Nz, Nt, Nmu, N_azi, Ng = make_type_mesh(card["ww_mesh"])
     struct += [("ww_mesh", mesh)]
     struct += [("ww_width", float64)]
+    struct += [("ww_auto", bool_)]
 
     # Window
     struct += [("ww", float64, (Nt, Nx, Ny, Nz))]