Merge pull request #470 from steps-re/fix/neutron-physics-sampling

fix: correct several neutron sampling / tally bugs

Author: ilhamv

mcdc/object_/tally.py

@@ -198,7 +198,7 @@ def __init__(
             self.filter_direction = True
         if polar_reference is not None:
             polar_reference = np.array(polar_reference)
-            self.polar_reference /= polar_reference / np.linalg.norm(polar_reference)
+            self.polar_reference = polar_reference / np.linalg.norm(polar_reference)
         if energy is not None:
             if type(energy) == str and energy == "all_groups":
                 G = simulation.materials[0].G

mcdc/transport/distribution.py

@@ -154,6 +154,10 @@ def sample_direction(polar_cosine, azimuthal, polar_coordinate, rng_state):
     wy = polar_coordinate[1]
     wz = polar_coordinate[2]
     if abs(wz) >= 0.999:
+        # Axis nearly parallel to z: use a fixed transverse basis
+        ux, uy, uz = 1.0, 0.0, 0.0
+        vx, vy, vz = 0.0, 1.0, 0.0
+    else:
         inv = 1.0 / math.sqrt(wx * wx + wy * wy)
 
         ux = -wy * inv
@@ -163,10 +167,6 @@ def sample_direction(polar_cosine, azimuthal, polar_coordinate, rng_state):
         vx = -wz * wx * inv
         vy = -wz * wy * inv
         vz = math.sqrt(wx * wx + wy * wy)
-    else:
-        # Axis nearly parallel to z
-        ux, uy, uz = 1.0, 0.0, 0.0
-        vx, vy, vz = 0.0, 1.0, 0.0
 
     # Rotate into lab frame
     s = math.sqrt(max(0.0, 1.0 - mu * mu))
@@ -271,7 +271,7 @@ def sample_white_direction(nx, ny, nz, rng_state):
     sin_azi = math.sin(azi)
     Ac = (1.0 - mu**2) ** 0.5
 
-    if nz != 1.0:
+    if abs(nz) != 1.0:
         B = (1.0 - nz**2) ** 0.5
         C = Ac / B
 

mcdc/transport/physics/neutron/native.py

@@ -828,7 +828,7 @@ def fission(
         prompt = True
         delayed_group = -1
         xi = rng.lcg(particle_container_new)
-        total = nu_p
+        total = nu_p / nu
         if xi > total:
             prompt = False
             # Determine delayed group
@@ -912,7 +912,7 @@ def fission(
             ux_new, uy_new, uz_new = sample_isotropic_direction(particle_container_new)
 
             # Sample emission time
-            decay_rate = mcdc_get.nuclide.neutron_fission_delayed_fractions(
+            decay_rate = mcdc_get.nuclide.neutron_fission_delayed_decay_rates(
                 delayed_group, nuclide, data
             )
             if not prompt: