fix mesh for most of the test problems

mcdc/constant.py

@@ -67,11 +67,12 @@
 PREC = 1.0 + 1e-5  # Precision factor to determine if a distance is smaller
 BANKMAX = 100  # Default maximum active bank
 
-# Domain Decomp mesh crossing flags
+# Mesh crossing flags
 MESH_X = 0
 MESH_Y = 1
 MESH_Z = 2
 MESH_T = 3
+MESH_NONE = 0
 
 # RNG LCG parameters
 RNG_G = nb.uint64(2806196910506780709)

mcdc/iqmc/iqmc_kernel.py

@@ -300,11 +300,18 @@ def iqmc_prepare_particles(mcdc):
         P_new["ux"], P_new["uy"], P_new["uz"] = iqmc_sample_isotropic_direction(
             lds[n, 1], lds[n, 5]
         )
-        x, y, z, t, outside = mesh.get_indices(P_new, mesh)
-        q = Q[:, t, x, y, z].copy()
-        dV = iqmc_cell_volume(x, y, z, mesh)
+        x = P_new['x']
+        y = P_new['y']
+        z = P_new['z']
+        t = P_new['t']
+        ux = P_new['ux']
+        uy = P_new['uy']
+        uz = P_new['uz']
+        ix, iy, iz, it, outside = mesh_.get_indices(x, y, z, t, ux, uy, uz, mesh)
+        q = Q[:, it, ix, iy, iz].copy()
+        dV = iqmc_cell_volume(ix, iy, iz, mesh)
         # Source tilt
-        iqmc_tilt_source(t, x, y, z, P_new, q, mcdc)
+        iqmc_tilt_source(it, ix, iy, iz, P_new, q, mcdc)
         # set particle weight
         P_new["iqmc"]["w"] = q * dV * N_total / N_particle
         P_new["w"] = P_new["iqmc"]["w"].sum()
@@ -340,7 +347,15 @@ def iqmc_score_tallies(P, distance, mcdc):
     SigmaT = material["total"]
     mat_id = P["material_ID"]
 
-    x, y, z, t, outside = mesh.get_indices(P, mesh)
+    x_ = P['x']
+    y_ = P['y']
+    z_ = P['z']
+    t_ = P['t']
+    ux_ = P['ux']
+    uy_ = P['uy']
+    uz_ = P['uz']
+
+    x, y, z, t, outside = mesh_.get_indices(x_, y_, z_, t_, ux_, uy_, uz_, mesh)
     if outside:
         return
 

mcdc/iqmc/iqmc_loop.py

@@ -67,7 +67,6 @@ def iqmc_step_particle(P, prog):
 
     # Lattice or mesh crossing (skipped if surface crossing)
     elif event & EVENT_LATTICE or event & EVENT_MESH:
-        kernel.shift_particle(P, SHIFT)
         if event & EVENT_DOMAIN:
             kernel.domain_crossing(P, mcdc)
 

mcdc/kernel.py

@@ -40,18 +40,32 @@ def domain_crossing(P, mcdc):
     if mcdc["technique"]["domain_decomposition"]:
         mesh = mcdc["technique"]["dd_mesh"]
         # Determine which dimension is crossed
-        x, y, z, t, directions = mesh_crossing_evaluate(P, mesh)
-        if len(directions) == 0:
-            return
-        elif len(directions) > 1:
-            for direction in directions[1:]:
-                if direction == MESH_X:
-                    P["x"] -= SHIFT * P["ux"] / np.abs(P["ux"])
-                if direction == MESH_Y:
-                    P["y"] -= SHIFT * P["uy"] / np.abs(P["uy"])
-                if direction == MESH_Z:
-                    P["z"] -= SHIFT * P["uz"] / np.abs(P["uz"])
-        flag = directions[0]
+        x = P['x']
+        y = P['y']
+        z = P['z']
+        t = P['t']
+        ux = P['ux']
+        uy = P['uy']
+        uz = P['uz']
+        ix, iy, iz, it, outside = mesh_.get_indices(x, y, z, t, ux, uy, uz, mesh)
+
+        d_idx = mcdc['dd_idx']
+        d_Nx = mcdc["technique"]["dd_mesh"]["x"].size - 1
+        d_Ny = mcdc["technique"]["dd_mesh"]["y"].size - 1
+        d_Nz = mcdc["technique"]["dd_mesh"]["z"].size - 1
+
+        d_iz = int(d_idx / (d_Nx * d_Ny))
+        d_iy = int((d_idx - d_Nx * d_Ny * d_iz) / d_Nx)
+        d_ix = int(d_idx - d_Nx * d_Ny * d_iz - d_Nx * d_iy)
+
+        flag = MESH_NONE
+        if d_ix != ix:
+            flag = MESH_X
+        elif d_iy != iy:
+            flag = MESH_Y
+        elif d_iz != iz:
+            flag = MESH_Z
+
         # Score on tally
         if flag == MESH_X and P["ux"] > 0:
             add_particle(P, mcdc["domain_decomp"]["bank_xp"])
@@ -408,9 +422,16 @@ def particle_in_domain(P, mcdc):
     d_iy = int((d_idx - d_Nx * d_Ny * d_iz) / d_Nx)
     d_ix = int(d_idx - d_Nx * d_Ny * d_iz - d_Nx * d_iy)
 
-    x_cell = binary_search(P["x"], mcdc["technique"]["dd_mesh"]["x"])
-    y_cell = binary_search(P["y"], mcdc["technique"]["dd_mesh"]["y"])
-    z_cell = binary_search(P["z"], mcdc["technique"]["dd_mesh"]["z"])
+    mesh = mcdc["technique"]["dd_mesh"]
+    x = P['x']
+    y = P['y']
+    z = P['z']
+    t = P['t']
+    ux = P['ux']
+    uy = P['uy']
+    uz = P['uz']
+
+    x_cell, y_cell, z_cell, t_cell, outside = mesh_.get_indices(x, y, z, t, ux, uy, uz, mesh)
 
     if d_ix == x_cell:
         if d_iy == y_cell:
@@ -1654,32 +1675,6 @@ def mesh_get_energy_index(P, mesh, mcdc):
         return binary_search(P["E"], mesh["g"]), outside
 
 
-@njit
-def mesh_crossing_evaluate(P, mesh):
-    # Shift backward
-    shift_particle(P, -2 * SHIFT)
-    x1, y1, z1, t1, outside1 = mesh_.get_indices(P, mesh)
-
-    # Double shift forward
-    shift_particle(P, 4 * SHIFT)
-    x2, y2, z2, t2, outside2 = mesh_.get_indices(P, mesh)
-
-    # Return particle to initial position
-    shift_particle(P, -2 * SHIFT)
-
-    # Determine dimension crossed
-    directions = []
-
-    if x1 != x2:
-        directions.append(MESH_X)
-    if y1 != y2:
-        directions.append(MESH_Y)
-    if z1 != z2:
-        directions.append(MESH_Z)
-
-    return x1, y1, z1, t1, directions
-
-
 # =============================================================================
 # Tally operations
 # =============================================================================

test/unit/test_mesh.py

@@ -0,0 +1,193 @@
+import numpy as np
+
+from mcdc.constant import (
+    COINCIDENCE_TOLERANCE,
+    INF,
+)
+
+from mcdc.mesh import (
+    get_indices,
+    get_grid_index,
+    nearest_distance_to_grid,
+)
+
+
+def test_get_grid_index():
+    grid = np.array([-6.0, -3.0, -1.0, 0.0, 1.0, 3.0, 6.0])
+    tiny = COINCIDENCE_TOLERANCE * 0.8
+
+    # Inside bin, going right
+    assert np.isclose(get_grid_index(-3.2, 0.4, grid), 0)
+    assert np.isclose(get_grid_index(3.2, 0.4, grid), 5)
+
+    # Inside bin, going left
+    assert np.isclose(get_grid_index(-3.2, -0.4, grid), 0)
+    assert np.isclose(get_grid_index(3.2, -0.4, grid), 5)
+
+    # At internal grid, going right
+    assert np.isclose(get_grid_index(3.0, 0.4, grid), 5)
+    assert np.isclose(get_grid_index(-3.0, 0.4, grid), 1)
+
+    # At internal grid, going left
+    assert np.isclose(get_grid_index(3.0, -0.4, grid), 4)
+    assert np.isclose(get_grid_index(-3.0, -0.4, grid), 0)
+
+    # At left-most grid, going right
+    assert np.isclose(get_grid_index(-6.0, 0.4, grid), 0)
+
+    # At right-most grid, going left
+    assert np.isclose(get_grid_index(6.0, -0.4, grid), 5)
+
+    # At internal grid (within tolerance), going right
+    assert np.isclose(get_grid_index(3.0 + tiny, 0.4, grid), 5)
+    assert np.isclose(get_grid_index(3.0 - tiny, 0.4, grid), 5)
+    assert np.isclose(get_grid_index(-3.0 + tiny, 0.4, grid), 1)
+    assert np.isclose(get_grid_index(-3.0 - tiny, 0.4, grid), 1)
+
+    # At internal grid (within tolerance), going left
+    assert np.isclose(get_grid_index(3.0 + tiny, -0.4, grid), 4)
+    assert np.isclose(get_grid_index(3.0 - tiny, -0.4, grid), 4)
+    assert np.isclose(get_grid_index(-3.0 + tiny, -0.4, grid), 0)
+    assert np.isclose(get_grid_index(-3.0 - tiny, -0.4, grid), 0)
+
+    # At left-most grid (within tolerance), going right
+    assert np.isclose(get_grid_index(-6.0 + tiny, 0.4, grid), 0)
+    assert np.isclose(get_grid_index(-6.0 - tiny, 0.4, grid), 0)
+
+    # At right-most grid (within tolerance), going left
+    assert np.isclose(get_grid_index(6.0 + tiny, -0.4, grid), 5)
+    assert np.isclose(get_grid_index(6.0 - tiny, -0.4, grid), 5)
+
+
+def test_get_indices():
+    mesh = {}
+    mesh['x'] = np.array([-6.0, -3.0, -1.0, 0.0, 1.0, 3.0, 6.0])
+    mesh['y'] = np.array([-6.0, -3.0, -1.0, 0.0, 1.0, 3.0, 6.0])
+    mesh['z'] = np.array([-6.0, -3.0, -1.0, 0.0, 1.0, 3.0, 6.0])
+    mesh['t'] = np.array([-6.0, -3.0, -1.0, 0.0, 1.0, 3.0, 6.0])
+    tiny = COINCIDENCE_TOLERANCE * 0.8
+
+    # Inside bin
+    ix, iy, iz, it, outside = get_indices(-5.0, -2.0, -0.5, 4.0, 1.0, 1.0, 1.0, mesh)
+    assert ix == 0 and iy == 1 and iz == 2 and it == 5 and not outside
+    ix, iy, iz, it, outside = get_indices(-5.0, -2.0, -0.5, 4.0, -1.0, -1.0, -1.0, mesh)
+    assert ix == 0 and iy == 1 and iz == 2 and it == 5 and not outside
+
+    # Outside
+    ix, iy, iz, it, outside = get_indices(-5.0, -2.0, -0.5, 10.0, 1.0, 1.0, 1.0, mesh)
+    assert ix == -1 and iy == -1 and iz == -1 and it == -1 and outside
+    ix, iy, iz, it, outside = get_indices(-5.0, -2.0, -0.5, 10.0, -1.0, -1.0, -1.0, mesh)
+    assert ix == -1 and iy == -1 and iz == -1 and it == -1 and outside
+
+    # At internal grid
+    ix, iy, iz, it, outside = get_indices(-3.0, -1.0, 1.0, 3.0, 1.0, 1.0, 1.0, mesh)
+    assert ix == 1 and iy == 2 and iz == 4 and it == 5 and not outside
+    ix, iy, iz, it, outside = get_indices(-3.0, -1.0, 1.0, 3.0, -1.0, -1.0, -1.0, mesh)
+    assert ix == 0 and iy == 1 and iz == 3 and it == 5 and not outside
+
+    # At left-most grid
+    ix, iy, iz, it, outside = get_indices(-6.0, -6.0, -6.0, 0.0, 1.0, 1.0, 1.0, mesh)
+    assert ix == 0 and iy == 0 and iz == 0 and it == 3 and not outside
+
+    # At right-most grid
+    ix, iy, iz, it, outside = get_indices(6.0, 6.0, 6.0, 0.0, -1.0, -1.0, -1.0, mesh)
+    assert ix == 5 and iy == 5 and iz == 5 and it == 3 and not outside
+
+    # At internal grid (within tolerance)
+    ix, iy, iz, it, outside = get_indices(-3.0 + tiny, -1.0 + tiny, 1.0 + tiny, 3.0 + tiny, 1.0, 1.0, 1.0, mesh)
+    assert ix == 1 and iy == 2 and iz == 4 and it == 5 and not outside
+    ix, iy, iz, it, outside = get_indices(-3.0 - tiny, -1.0 - tiny, 1.0 - tiny, 3.0 - tiny, 1.0, 1.0, 1.0, mesh)
+    assert ix == 1 and iy == 2 and iz == 4 and it == 5 and not outside
+    ix, iy, iz, it, outside = get_indices(-3.0 + tiny, -1.0 + tiny, 1.0 + tiny, 3.0 + tiny, -1.0, -1.0, -1.0, mesh)
+    assert ix == 0 and iy == 1 and iz == 3 and it == 5 and not outside
+    ix, iy, iz, it, outside = get_indices(-3.0 - tiny, -1.0 - tiny, 1.0 - tiny, 3.0 - tiny, -1.0, -1.0, -1.0, mesh)
+    assert ix == 0 and iy == 1 and iz == 3 and it == 5 and not outside
+
+    # At left-most grid (within tolerance)
+    ix, iy, iz, it, outside = get_indices(-6.0 + tiny, -6.0 + tiny, -6.0 + tiny, 0.0 + tiny, 1.0, 1.0, 1.0, mesh)
+    assert ix == 0 and iy == 0 and iz == 0 and it == 3 and not outside
+    ix, iy, iz, it, outside = get_indices(-6.0 - tiny, -6.0 - tiny, -6.0 - tiny, 0.0 - tiny, 1.0, 1.0, 1.0, mesh)
+    assert ix == 0 and iy == 0 and iz == 0 and it == 3 and not outside
+    ix, iy, iz, it, outside = get_indices(-6.0 + tiny, -6.0 + tiny, -6.0 + tiny, 0.0 + tiny, -1.0, -1.0, -1.0, mesh)
+    assert ix == -1 and iy == -1 and iz == -1 and it == -1 and outside
+    ix, iy, iz, it, outside = get_indices(-6.0 - tiny, -6.0 - tiny, -6.0 - tiny, 0.0 - tiny, -1.0, -1.0, -1.0, mesh)
+    assert ix == -1 and iy == -1 and iz == -1 and it == -1 and outside
+
+    # At right-most grid (within tolerance)
+    ix, iy, iz, it, outside = get_indices(6.0 + tiny, 6.0 + tiny, 6.0 + tiny, 0.0 + tiny, 1.0, 1.0, 1.0, mesh)
+    assert ix == -1 and iy == -1 and iz == -1 and it == -1 and outside
+    ix, iy, iz, it, outside = get_indices(6.0 - tiny, 6.0 - tiny, 6.0 - tiny, 0.0 - tiny, 1.0, 1.0, 1.0, mesh)
+    assert ix == -1 and iy == -1 and iz == -1 and it == -1 and outside
+    ix, iy, iz, it, outside = get_indices(6.0 + tiny, 6.0 + tiny, 6.0 + tiny, 0.0 + tiny, -1.0, -1.0, -1.0, mesh)
+    assert ix == 5 and iy == 5 and iz == 5 and it == 3 and not outside
+    ix, iy, iz, it, outside = get_indices(6.0 - tiny, 6.0 - tiny, 6.0 - tiny, 0.0 - tiny, -1.0, -1.0, -1.0, mesh)
+    assert ix == 5 and iy == 5 and iz == 5 and it == 3 and not outside
+
+
+def test_nearest_distance_to_grid():
+    grid = np.array([-6.0, -3.0, -1.0, 0.0, 1.0, 3.0, 6.0])
+    tiny = COINCIDENCE_TOLERANCE * 0.8
+
+    # Inside bin, going right
+    assert np.isclose(nearest_distance_to_grid(-3.2, 0.4, grid), 0.2/0.4)
+    assert np.isclose(nearest_distance_to_grid(3.2, 0.4, grid), 2.8/0.4)
+
+    # Inside bin, going left
+    assert np.isclose(nearest_distance_to_grid(-3.2, -0.4, grid), 2.8/0.4)
+    assert np.isclose(nearest_distance_to_grid(3.2, -0.4, grid), 0.2/0.4)
+
+    # Outside, moving away
+    assert np.isclose(nearest_distance_to_grid(8.0, 0.4, grid), INF)
+    assert np.isclose(nearest_distance_to_grid(-8.0, -0.4, grid), INF)
+
+    # Outside, moving closer
+    assert np.isclose(nearest_distance_to_grid(8.0, -0.4, grid), 2.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(-8.0, 0.4, grid), 2.0/0.4)
+
+    # At internal grid, going right
+    assert np.isclose(nearest_distance_to_grid(3.0, 0.4, grid), 3.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(-3.0, 0.4, grid), 2.0/0.4)
+
+    # At internal grid, going left
+    assert np.isclose(nearest_distance_to_grid(3.0, -0.4, grid), 2.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(-3.0, -0.4, grid), 3.0/0.4)
+
+    # At left-most grid, going right
+    assert np.isclose(nearest_distance_to_grid(-6.0, 0.4, grid), 3.0/0.4)
+
+    # At left-most grid, going left
+    assert np.isclose(nearest_distance_to_grid(-6.0, -0.4, grid), INF)
+
+    # At right-most grid, going right
+    assert np.isclose(nearest_distance_to_grid(6.0, 0.4, grid), INF)
+
+    # At right-most grid, going left
+    assert np.isclose(nearest_distance_to_grid(6.0, -0.4, grid), 3.0/0.4)
+
+    # At internal grid (within tolerance), going right
+    assert np.isclose(nearest_distance_to_grid(3.0 + tiny, 0.4, grid), 3.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(3.0 - tiny, 0.4, grid), 3.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(-3.0 + tiny, 0.4, grid), 2.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(-3.0 - tiny, 0.4, grid), 2.0/0.4)
+
+    # At internal grid (within tolerance), going left
+    assert np.isclose(nearest_distance_to_grid(3.0 + tiny, -0.4, grid), 2.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(3.0 - tiny, -0.4, grid), 2.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(-3.0 + tiny, -0.4, grid), 3.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(-3.0 - tiny, -0.4, grid), 3.0/0.4)
+
+    # At left-most grid (within tolerance), going right
+    assert np.isclose(nearest_distance_to_grid(-6.0 + tiny, 0.4, grid), 3.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(-6.0 - tiny, 0.4, grid), 3.0/0.4)
+
+    # At left-most grid (within tolerance), going left
+    assert np.isclose(nearest_distance_to_grid(-6.0 + tiny, -0.4, grid), INF)
+    assert np.isclose(nearest_distance_to_grid(-6.0 - tiny, -0.4, grid), INF)
+
+    # At right-most grid (within tolerance), going right
+    assert np.isclose(nearest_distance_to_grid(6.0 + tiny, 0.4, grid), INF)
+    assert np.isclose(nearest_distance_to_grid(6.0 - tiny, 0.4, grid), INF)
+
+    # At right-most grid (within tolerance), going left
+    assert np.isclose(nearest_distance_to_grid(6.0 + tiny, -0.4, grid), 3.0/0.4)
+    assert np.isclose(nearest_distance_to_grid(6.0 - tiny, -0.4, grid), 3.0/0.4)