Continious energy generator bugfixing

examples/slab_sens/input.py

@@ -0,0 +1,71 @@
+import os
+import numpy as np
+import mcdc
+
+# -----------------------------------------------------------------------------
+# Continuous-energy data
+# -----------------------------------------------------------------------------
+# You must point to a directory of nuclide HDF5s, e.g.:
+#   export MCDC_LIB=/path/to/mcdc_ce_library
+#
+# Nuclide names must match the file prefix in that directory.
+# Example filename: H1-293.6K.h5  -> nuclide name "H1"
+#
+if "MCDC_LIB" not in os.environ:
+    raise RuntimeError("Set MCDC_LIB to your continuous-energy nuclear data directory")
+
+# Material: simple scatterer (use a nuclide you actually have in MCDC_LIB)
+mat = mcdc.Material(nuclide_composition={"H1": 0.066})  # density units per your library
+
+# -----------------------------------------------------------------------------
+# Geometry: slab in x, reflective in y/z
+# -----------------------------------------------------------------------------
+L = 10.0
+xm = 5.0
+
+sx0 = mcdc.Surface.PlaneX(x=0.0, boundary_condition="vacuum")
+sxm = mcdc.Surface.PlaneX(x=xm)
+sxL = mcdc.Surface.PlaneX(x=L, boundary_condition="vacuum")
+
+sy0 = mcdc.Surface.PlaneY(y=-1.0, boundary_condition="reflective")
+sy1 = mcdc.Surface.PlaneY(y=+1.0, boundary_condition="reflective")
+sz0 = mcdc.Surface.PlaneZ(z=-1.0, boundary_condition="reflective")
+sz1 = mcdc.Surface.PlaneZ(z=+1.0, boundary_condition="reflective")
+
+left = mcdc.Cell(
+    name="left", region=(+sx0 & -sxm & +sy0 & -sy1 & +sz0 & -sz1), fill=mat
+)
+right = mcdc.Cell(
+    name="right", region=(+sxm & -sxL & +sy0 & -sy1 & +sz0 & -sz1), fill=mat
+)
+
+# -----------------------------------------------------------------------------
+# Source: isotropic point near left boundary
+# -----------------------------------------------------------------------------
+mcdc.Source(
+    position=[0.1, 0.0, 0.0],
+    isotropic=True,
+    energy=2.0e6,  # eV
+)
+
+# -----------------------------------------------------------------------------
+# Tallies: flux in each cell (baseline sanity check)
+# -----------------------------------------------------------------------------
+mcdc.TallyCell(name="flux_left", cell=left, scores=["flux"])
+mcdc.TallyCell(name="flux_right", cell=right, scores=["flux"])
+
+# -----------------------------------------------------------------------------
+# Sensitivity response regions: per-cell resp_cum
+# -----------------------------------------------------------------------------
+mcdc.settings.sensitivity_mode = True
+mcdc.settings.sensitivity_n_resp = 2
+mcdc.settings.sensitivity_resp_cell_IDs = np.array([left.ID, right.ID], dtype=np.int64)
+
+# -----------------------------------------------------------------------------
+# Run settings
+# -----------------------------------------------------------------------------
+mcdc.settings.N_particle = 20000
+mcdc.settings.N_batch = 10
+mcdc.simulation.implicit_capture()
+
+mcdc.run()

mcdc/code_factory/adapt.py

@@ -326,7 +326,11 @@ def toggle_inner(func):
 
 
 def set_toggle(flag, val):
-    toggle_rosters[flag][0] = val
+    global toggle_rosters
+    if flag not in toggle_rosters:
+        toggle_rosters[flag] = [bool(val), []]
+    else:
+        toggle_rosters[flag][0] = bool(val)
 
 
 def eval_toggle():

mcdc/code_factory/numba_objects_generator.py

@@ -170,6 +170,29 @@ def generate_numba_objects(simulation):
             )
         }
 
+    # ==================================================================================
+    # Sensitivity (optional): inject per-particle cumulative response vector `resp_cum`
+    # only when requested. This keeps particle records lean when sensitivity is off.
+    #
+    # We add the field to BOTH `particle_data` (used in banks) and `particle` (used in
+    # transport). The field is a fixed-size vector of length `sensitivity_n_resp`.
+    # ==================================================================================
+    sens_on = bool(getattr(simulation.settings, "sensitivity_mode", False))
+    if sens_on:
+        n_resp = int(getattr(simulation.settings, "sensitivity_n_resp", 0))
+        if n_resp <= 0:
+            raise ValueError(
+                "sensitivity_mode=True but sensitivity_n_resp is not a positive integer"
+            )
+
+        resp_hint = parse_type_hint_str(f"Annotated[NDArray[float64], ({n_resp},)]")
+        annotations["particle_data"]["resp_cum"] = resp_hint
+        annotations["particle"]["resp_cum"] = resp_hint
+
+        # NEW: response-region cell IDs (length = n_resp)
+        cell_hint = parse_type_hint_str(f"Annotated[NDArray[int64], ({n_resp},)]")
+        annotations["settings"]["sensitivity_resp_cell_IDs"] = cell_hint
+
     # ==================================================================================
     # Set the structures and accessor targets based on the annotations
     # ==================================================================================

mcdc/main.py

@@ -5,6 +5,7 @@
 
 from mcdc import mcdc_get
 from mcdc.print_ import print_error, print_structure
+import importlib
 
 
 def run():
@@ -35,6 +36,9 @@ def run():
     # Preparation
     # ==================================================================================
 
+    # To compile adaptive sensitivities
+    import mcdc.transport.simulation
+
     # Timer: preparation
     time_prep_start = MPI.Wtime()
 
@@ -206,11 +210,18 @@ def preparation():
         from mcdc.code_factory.numba_objects_generator import make_literals
 
         make_literals(simulation)
+
+    MPI.COMM_WORLD.Barrier()
+    importlib.invalidate_caches()
+
+    import mcdc.transport.literals as literals
+
+    importlib.reload(literals)
+
     mcdc_arr, data = generate_numba_objects(simulation)
     mcdc = mcdc_arr[0]
 
     # Reload mcdc getters and setters
-    import importlib
     import mcdc.mcdc_get as mcdc_get
     import mcdc.mcdc_set as mcdc_set
 
@@ -250,6 +261,10 @@ def preparation():
             type_.particle,
             type_.particle_data,
         )
+    # Optional sensitivity module: keep kernels/data minimal unless requested
+    adapt.set_toggle(
+        "sensitivity", bool(getattr(simulation.settings, "sensitivity_mode", False))
+    )
 
     adapt.eval_toggle()
     adapt.target_for(config.target)

mcdc/mcdc_get/settings.py

@@ -30,3 +30,32 @@ def census_time_chunk(start, length, settings, data):
     start += settings["census_time_offset"]
     end = start + length
     return data[start:end]
+
+
+@njit
+def sensitivity_resp_cell_IDs(index, settings, data):
+    offset = settings["sensitivity_resp_cell_IDs_offset"]
+    return data[offset + index]
+
+
+@njit
+def sensitivity_resp_cell_IDs_all(settings, data):
+    start = settings["sensitivity_resp_cell_IDs_offset"]
+    size = settings["sensitivity_resp_cell_IDs_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def sensitivity_resp_cell_IDs_last(settings, data):
+    start = settings["sensitivity_resp_cell_IDs_offset"]
+    size = settings["sensitivity_resp_cell_IDs_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def sensitivity_resp_cell_IDs_chunk(start, length, settings, data):
+    start += settings["sensitivity_resp_cell_IDs_offset"]
+    end = start + length
+    return data[start:end]

mcdc/mcdc_set/settings.py

@@ -30,3 +30,32 @@ def census_time_chunk(start, length, settings, data, value):
     start += settings["census_time_offset"]
     end = start + length
     data[start:end] = value
+
+
+@njit
+def sensitivity_resp_cell_IDs(index, settings, data, value):
+    offset = settings["sensitivity_resp_cell_IDs_offset"]
+    data[offset + index] = value
+
+
+@njit
+def sensitivity_resp_cell_IDs_all(settings, data, value):
+    start = settings["sensitivity_resp_cell_IDs_offset"]
+    size = settings["sensitivity_resp_cell_IDs_length"]
+    end = start + size
+    data[start:end] = value
+
+
+@njit
+def sensitivity_resp_cell_IDs_last(settings, data, value):
+    start = settings["sensitivity_resp_cell_IDs_offset"]
+    size = settings["sensitivity_resp_cell_IDs_length"]
+    end = start + size
+    data[end - 1] = value
+
+
+@njit
+def sensitivity_resp_cell_IDs_chunk(start, length, settings, data, value):
+    start += settings["sensitivity_resp_cell_IDs_offset"]
+    end = start + length
+    data[start:end] = value

mcdc/object_/settings.py

@@ -60,6 +60,13 @@ class Settings(ObjectSingleton):
     source_bank_buffer_ratio: float = 2.0
     future_bank_buffer_ratio: float = 1.5
 
+    # Sensitivity analysis (optional)
+    sensitivity_mode: bool = False
+    sensitivity_n_resp: int = 0
+    sensitivity_resp_cell_IDs: NDArray[np.int64] = field(
+        default_factory=lambda: np.zeros(0, dtype=np.int64)
+    )
+
     # Portability
     target_gpu: bool = False
 

mcdc/transport/particle.py

@@ -5,6 +5,22 @@
 import mcdc.transport.physics as physics
 import mcdc.transport.rng as rng
 
+import mcdc.code_factory.adapt as adapt
+
+
+@adapt.toggle("sensitivity")
+def _copy_resp_cum(target_particle_container, source_particle_container):
+    """Copy per-particle response accumulator (sensitivity mode only)."""
+    target = target_particle_container[0]
+    source = source_particle_container[0]
+    target["resp_cum"][:] = source["resp_cum"][:]
+
+
+@adapt.toggle("sensitivity")
+def _reset_resp_cum(particle_container):
+    """Reset per-particle response accumulator for a new child history."""
+    particle_container[0]["resp_cum"][:] = 0.0
+
 
 @njit
 def move(particle_container, distance, mcdc, data):
@@ -35,6 +51,8 @@ def copy(target_particle_container, source_particle_container):
     target_particle["particle_type"] = source_particle["particle_type"]
     target_particle["rng_seed"] = source_particle["rng_seed"]
 
+    _copy_resp_cum(target_particle_container, source_particle_container)
+
 
 @njit
 def copy_as_child(child_particle_container, parent_particle_container):
@@ -43,6 +61,8 @@ def copy_as_child(child_particle_container, parent_particle_container):
 
     copy(child_particle_container, parent_particle_container)
 
+    _reset_resp_cum(child_particle_container)
+
     # Set child RNG seed based of the parent
     parent_seed = parent_particle["rng_seed"]
     child_particle["rng_seed"] = rng.split_seed(parent_seed, rng.SEED_SPLIT_PARTICLE)

mcdc/transport/particle_bank.py

@@ -20,6 +20,12 @@
 from mcdc.print_ import print_error
 
 
+@adapt.toggle("sensitivity")
+def _unbank_resp_cum(P_arr, P_rec):
+    """Copy resp_cum from a bank particle record into a live particle (sens only)."""
+    P_arr[0]["resp_cum"][:] = P_rec["resp_cum"][:]
+
+
 # =============================================================================
 # Particle bank operations
 # =============================================================================
@@ -81,6 +87,8 @@ def get_particle(P_arr, bank, mcdc):
     P["particle_type"] = P_rec["particle_type"]
     P["rng_seed"] = P_rec["rng_seed"]
 
+    _unbank_resp_cum(P_arr, P_rec)
+
     # Set default IDs and event
     P["alive"] = True
     P["material_ID"] = -1

mcdc/transport/sensitivity.py

@@ -0,0 +1,24 @@
+import mcdc.code_factory.adapt as adapt
+
+
+@adapt.toggle("sensitivity")
+def score_resp_cum_tracklength(particle_container, distance, mcdc, data):
+    """
+    Sensitivity-only track-length scoring into particle["resp_cum"][k],
+    where k indexes the user-defined response regions.
+
+    Current definition: response regions are specified by cell_ID membership.
+    """
+    # Keep consistent with other tallying behavior
+    if not mcdc["cycle_active"]:
+        return
+
+    P = particle_container[0]
+    cell_id = P["cell_ID"]
+
+    resp_cells = mcdc["settings"]["sensitivity_resp_cell_IDs"]
+    # Find which response region this cell corresponds to (small N_resp => linear scan OK)
+    for k in range(resp_cells.shape[0]):
+        if cell_id == resp_cells[k]:
+            P["resp_cum"][k] += P["w"] * distance
+            break

mcdc/transport/simulation.py

@@ -18,6 +18,7 @@
 import mcdc.transport.rng as rng
 import mcdc.transport.tally as tally_module
 import mcdc.transport.technique as technique
+import mcdc.transport.sensitivity as sensitivity
 
 from mcdc.constant import *
 from mcdc.print_ import (
@@ -622,6 +623,10 @@ def move_to_event(particle_container, mcdc, data):
     if settings["eigenvalue_mode"]:
         tally_module.score.eigenvalue_tally(particle_container, distance, mcdc, data)
 
+    if settings["sensitivity_mode"]:
+        # Sensitivity response accumulator (track-length)
+        sensitivity.score_resp_cum_tracklength(particle_container, distance, mcdc, data)
+
     # Move particle
     particle_module.move(particle_container, distance, mcdc, data)
 

mcdc/transport/source.py

@@ -15,6 +15,14 @@
 )
 from mcdc.transport.util import find_bin
 
+import mcdc.code_factory.adapt as adapt
+
+
+@adapt.toggle("sensitivity")
+def _init_resp_cum(P_rec_arr):
+    """Initialize per-particle response accumulator for a new source history."""
+    P_rec_arr[0]["resp_cum"][:] = 0.0
+
 
 @njit
 def source_particle(P_rec_arr, seed, mcdc, data):
@@ -110,3 +118,5 @@ def source_particle(P_rec_arr, seed, mcdc, data):
     P_rec["E"] = E
     P_rec["w"] = 1.0
     P_rec["particle_type"] = source["particle_type"]
+
+    _init_resp_cum(P_rec_arr)

pyproject.toml

@@ -82,6 +82,7 @@ dev = [
   "black",
   "pre-commit",
   "pytest",
+  "tqdm", # for continuos energy processing
 ]
 
 [project.urls]