Delta tracking

docs/source/pythonapi/index.rst

@@ -63,6 +63,7 @@ Defining techniques
    mcdc.weighted_emission
    mcdc.weight_roulette
    mcdc.weight_window
+   mcdc.delta_tracking
 
 
 

mcdc/__init__.py

@@ -20,6 +20,7 @@
     weight_roulette,
     IC_generator,
     uq,
+    delta_tracking,
     reset,
     domain_decomposition,
     make_particle_bank,

mcdc/config.py

@@ -67,6 +67,12 @@
 parser.add_argument("--caching", action="store_true")
 parser.add_argument("--no_caching", dest="caching", action="store_false")
 parser.add_argument("--runtime_output", default=False, action="store_true")
+parser.add_argument(
+    "--delta_tracking",
+    default=False,
+    action="store_true",
+    help="Set delta tracking to on",
+)
 parser.set_defaults(caching=False)
 args, unargs = parser.parse_known_args()
 
@@ -77,6 +83,7 @@
 clear_cache = args.clear_cache
 
 from mpi4py import MPI
+
 import shutil
 
 src_path = os.path.dirname(os.path.abspath(__file__))

mcdc/constant.py

@@ -83,6 +83,7 @@
 EVENT_TIME_CENSUS = 1 << 9
 EVENT_TIME_BOUNDARY = 1 << 10
 EVENT_IQMC_MESH = 1 << 11
+EVENT_PHANTOM_COLLISION = 1 << 12
 
 # Gyration raius type
 GYRATION_RADIUS_ALL = 0
@@ -148,6 +149,7 @@
 XS_NU_FISSION_PROMPT = 5
 XS_NU_FISSION_DELAYED = 6
 XS_NU_SCATTER = 7
+XS_MAJORANT = 8
 
 NU_FISSION = 0
 NU_FISSION_PROMPT = 1

mcdc/global_.py

@@ -153,6 +153,11 @@ def reset(self):
             "IC_cycle_stretch": 1.0,
             "branchless_collision": False,
             "uq": False,
+            "delta_tracking": False,
+            "collision_estimator": False,
+            "dt_cutoff_energy": 0,
+            "dt_material_exclusion": INF,
+            "plot_majorant": True,
         }
 
         self.uq_deltas = {

mcdc/input_.py

@@ -227,6 +227,7 @@ def material(
     chi_d=None,
     speed=None,
     decay=None,
+    name=None,
 ):
     """
     Create a material
@@ -258,6 +259,8 @@ def material(
         Energy group speed [cm/s].
     decay : numpy.ndarray (1D), optional
         Precursor group decay constant [/s].
+    name : str, optional
+        Name of the material (e.g. moderator)
 
     Returns
     -------
@@ -324,7 +327,7 @@ def material(
                     print_error(
                         "Continuous energy data directory not configured \n       "
                         "see https://cement-psaapgithubio.readthedocs.io/en/latest"
-                        "/install.html#configuring-continuous-energy-library \n"
+                        "./install.html#configuring-continuous-energy-library \n"
                     )
 
                 # Fissionable flag
@@ -446,8 +449,7 @@ def material(
 
 
 def surface(type_, bc="interface", **kw):
-    """
-    Create a surface to define the region of a cell.
+    """Create a surface to define the region of a cell.
 
     Parameters
     ----------
@@ -1245,6 +1247,40 @@ def branchless_collision():
     card["weighted_emission"] = False
 
 
+def delta_tracking(
+    collision_estimator=False,
+    cutoff_energy=0,
+    material_exclusion=INF,
+    plot_majorant=True,
+):
+    """
+    Activate delta tracking and set hybrid parameters
+
+    Parameters
+    ----------
+    collision_estimator : bool
+        Use a collision estimator for flux tallies
+    cutoff_energy : float
+        A cut off energy, underwich do surface tracking (**continuous energy only**) (look at the majorant plot and make an informed decision!)
+    material_exclusion : int
+        Material ID of a single region not to do delta tracking in (TODO: Make multiple regions) (**continuous energy only**)
+    plot_majorant : bool
+        Saves a .pdf plot of the majorant and nuclides when produced
+
+    Notes
+    -----
+    TODO: turn off other incompatible tracking techniques if any (k-eigenvalue, fission flux tallies, )
+    Doesn't work with k-eigenvalue simulations
+
+    """
+    card = global_.input_deck.technique
+    card["delta_tracking"] = True
+    card["collision_estimator"] = collision_estimator
+    card["dt_cutoff_energy"] = cutoff_energy
+    card["dt_material_exclusion"] = material_exclusion
+    card["plot_majorant"] = plot_majorant
+
+
 def time_census(t, tally_frequency=None):
     """
     Set time-census boundaries.
@@ -1253,7 +1289,7 @@ def time_census(t, tally_frequency=None):
     ----------
     t : array_like[float]
         The time-census boundaries.
-    tally_frecuency : integer, optional
+    tally_frequency : integer, optional
         Number of uniform tally time mesh bins in census-based tallying.
         This overrides manual tally time mesh definitions.