Merge pull request #282 from northroj/feature/northroj/time_cell_tallies

Feature/northroj/time cell tallies

mcdc/constant.py

@@ -94,7 +94,7 @@
 
 # Misc.
 TINY = 1e-10
-COINCIDENCE_TOLERANCE = TINY
+COINCIDENCE_TOLERANCE = TINY * 1e1
 COINCIDENCE_TOLERANCE_TIME = TINY * 1e-2
 INF = 1e10
 PI = math.acos(-1.0)

mcdc/kernel.py

@@ -2041,25 +2041,75 @@ def score_cell_tally(P_arr, distance, tally, data, mcdc):
     P = P_arr[0]
     tally_bin = data[TALLY]
     material = mcdc["materials"][P["material_ID"]]
+    mesh = tally["filter"]
     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
-
-        adapt.global_add(tally_bin, (TALLY_SCORE, cell_idx + i), round(score))
-        # tally_bin[TALLY_SCORE, cell_idx + i] += score
+    # Particle/track properties
+    ut = 1.0 / physics.get_speed(P_arr, mcdc)
+    t = P["t"]
+    t_final = t + ut * distance
+    Nt = mesh["Nt"]
+
+    # Get starting indices
+    g, outside_energy = mesh_get_energy_index(P_arr, mesh, mcdc["setting"]["mode_MG"])
+
+    # Outside grid?
+    if (
+        t < mesh["t"][0] - COINCIDENCE_TOLERANCE_TIME
+        or t > mesh["t"][Nt] + COINCIDENCE_TOLERANCE_TIME
+        or (abs(t - mesh["t"][Nt]) < COINCIDENCE_TOLERANCE_TIME)
+        or outside_energy
+    ):
+        return
+
+    it = mesh_.structured._grid_index(
+        t, 1.0, mesh["t"], Nt + 1, COINCIDENCE_TOLERANCE_TIME
+    )
+
+    # Tally index
+    idx = stride["tally"] + g * stride["g"] + it * stride["t"]
+
+    # Sweep through the distance
+    distance_swept = 0.0
+    while distance_swept < distance - COINCIDENCE_TOLERANCE:
+
+        # Find distance to mesh grids
+        dt = (min(mesh["t"][it + 1], t_final) - t) / ut
+
+        # Get the grid crossed
+        distance_scored = INF
+        mesh_crossed = MESH_NONE
+        if dt <= distance_scored:
+            mesh_crossed = MESH_T
+            distance_scored = dt
+
+        # Score
+        flux = distance_scored * P["w"]
+        for i in range(tally["N_score"]):
+            score_type = tally["scores"][i]
+            score = 0
+            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
+            adapt.global_add(tally_bin, (TALLY_SCORE, idx + i), round(score))
+
+        # Accumulate distance swept
+        distance_swept += distance_scored
+
+        # Move the 4D (just t for this tally) position
+        t += distance_scored * ut
+
+        # Increment index and check if out of bounds
+        if mesh_crossed == MESH_T:
+            it += 1
+            if it == Nt:
+                break
+            idx += stride["t"]
 
 
 @njit

mcdc/main.py

@@ -1037,6 +1037,8 @@ def prepare():
         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
+        mcdc["cell_tallies"][i]["filter"]["Ng"] = Ng
+        mcdc["cell_tallies"][i]["filter"]["Nt"] = Nt
 
         # Update N_bin
         mcdc["cell_tallies"][i]["N_bin"] *= N_score
@@ -1941,13 +1943,26 @@ def generate_hdf5(data, mcdc):
                 if mcdc["technique"]["iQMC"]:
                     break
 
+                mesh = tally["filter"]
+
+                # Get grid
+                Nt = mesh["Nt"]
+                Ng = mesh["Ng"]
+                #
+                grid_t = mesh["t"][: Nt + 1]
+                grid_g = mesh["g"][: Ng + 1]
+
+                # Save to dataset
+                f.create_dataset("tallies/cell_tally_%i/grid/t" % ID, data=grid_t)
+                f.create_dataset("tallies/cell_tally_%i/grid/g" % ID, data=grid_g)
+
                 # Shape
                 N_score = tally["N_score"]
 
                 if not mcdc["technique"]["uq"]:
-                    shape = (3, N_score)
+                    shape = (3, Ng, Nt, N_score)
                 else:
-                    shape = (5, N_score)
+                    shape = (5, Ng, Nt, N_score)
 
                 # Reshape tally
                 N_bin = tally["N_bin"]
@@ -1956,7 +1971,7 @@ def generate_hdf5(data, mcdc):
                 tally_bin = tally_bin.reshape(shape)
 
                 # Roll tally so that score is in the front
-                tally_bin = np.rollaxis(tally_bin, 1, 0)
+                tally_bin = np.rollaxis(tally_bin, 3, 0)
 
                 # Iterate over scores
                 for i in range(N_score):

mcdc/tally.py

@@ -161,14 +161,28 @@ def surface_tally(
     return card
 
 
-def cell_tally(cell, scores=["flux"]):
+def cell_tally(
+    cell,
+    t=np.array([-INF, INF]),
+    g=np.array([-INF, INF]),
+    E=np.array([0.0, INF]),
+    scores=["flux"],
+):
     """
     Create a tally card to collect MC solutions.
 
     Parameters
     ----------
     cell : CellCard
         Cell to which the tally is attached to
+    t : array_like[float], optional
+        Times that demarcate tally bins.
+    g : array_like[float] or str, optional
+        Energy group halves that demarcate energy group tally bins.
+        String value "all" can be used to tally each individual group.
+    E : array_like[float], optional
+        Energies that demarcate energy tally bins. This overrides `g` in
+        continuous-energy mode.
     scores : list of str {"flux", "net-current"}
         List of physical quantities to be scored.
 
@@ -184,13 +198,25 @@ def cell_tally(cell, scores=["flux"]):
     # Set ID
     card.ID = len(global_.input_deck.cell_tallies)
 
+    card.t = t
+    # Set energy group grid
+    if type(g) == type("string") and g == "all":
+        G = global_.input_deck.materials[0].G
+        card.g = np.linspace(0, G, G + 1) - 0.5
+    else:
+        card.g = g
+    if global_.input_deck.setting["mode_CE"]:
+        card.g = E
+
     # 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
+    Nt = len(card.t) - 1
+    Ng = len(card.g) - 1
+    card.N_bin = Nt * Ng
 
     # Scores
     for s in scores:

mcdc/type_.py

@@ -889,6 +889,8 @@ def make_type_cell_tally(input_deck):
         ("mu", float64, (Nmax_mu,)),
         ("azi", float64, (Nmax_azi,)),
         ("g", float64, (Nmax_g,)),
+        ("Nt", int64),
+        ("Ng", int64),
     ]
     struct = [("filter", filter_)]
 

test/regression/inf_shem361_census/input.py

@@ -0,0 +1,65 @@
+import numpy as np
+import sys
+
+import mcdc
+
+# This regression test adds time census and time/energy binned cell tallies to the inf_shem361 test
+
+# =============================================================================
+# Set model
+# =============================================================================
+# The infinite homogenous medium is modeled with reflecting slab
+
+# Load material data
+with np.load("SHEM-361.npz") as data:
+    SigmaC = data["SigmaC"] * 5  # /cm
+    SigmaS = data["SigmaS"]
+    SigmaF = data["SigmaF"]
+    nu_p = data["nu_p"]
+    nu_d = data["nu_d"]
+    chi_p = data["chi_p"]
+    chi_d = data["chi_d"]
+    G = data["G"]
+    speed = data["v"]
+    lamd = data["lamd"]
+
+# Set material
+m = mcdc.material(
+    capture=SigmaC,
+    scatter=SigmaS,
+    fission=SigmaF,
+    nu_p=nu_p,
+    chi_p=chi_p,
+    nu_d=nu_d,
+    chi_d=chi_d,
+)
+
+# Set surfaces
+s1 = mcdc.surface("plane-x", x=-1e10, bc="reflective")
+s2 = mcdc.surface("plane-x", x=1e10, bc="reflective")
+
+# Set cells
+c = mcdc.cell(+s1 & -s2, m)
+
+# =============================================================================
+# Set initial source
+# =============================================================================
+
+energy = np.zeros(G)
+energy[-1] = 1.0
+source = mcdc.source(energy=energy)
+
+# =============================================================================
+# Set problem and tally, and then run mcdc
+# =============================================================================
+
+# Tally
+# mcdc.tally.mesh_tally(scores=["flux"], g="all")
+mcdc.tally.cell_tally(c, scores=["flux"], g="all", t=np.linspace(0.0, 20.0, 21)[1:-1])
+
+# Setting
+mcdc.setting(N_particle=1e2, active_bank_buff=1000)
+mcdc.time_census(np.linspace(0.0, 20.0, 21)[1:-1])
+
+# Run
+mcdc.run()