From 5d62cd5542490193f7dbb94fb0d5c5f0f6ec2a79 Mon Sep 17 00:00:00 2001
From: Ilham Variansyah <variansi@oregonstate.edu>
Date: Wed, 21 Aug 2024 17:32:22 +0700
Subject: [PATCH] fix quadric surface coincidence issue

---
 mcdc/geometry.py | 64 +++++++++---------------------------------------
 mcdc/loop.py     |  6 ++---
 2 files changed, 14 insertions(+), 56 deletions(-)

diff --git a/mcdc/geometry.py b/mcdc/geometry.py
index 52f633d0..5bb0c9fa 100644
--- a/mcdc/geometry.py
+++ b/mcdc/geometry.py
@@ -251,57 +251,6 @@ def distance_to_nearest_surface(particle, cell, mcdc):
     return distance, surface_ID
 
 
-@njit
-def distance_to_lattice_grid(particle, lattice):
-    # TODO: docs
-    d = INF
-    d = min(d, lattice_grid_distance(x, ux, lattice["x0"], lattice["dx"]))
-    d = min(d, lattice_grid_distance(y, uy, lattice["y0"], lattice["dy"]))
-    d = min(d, lattice_grid_distance(z, uz, lattice["z0"], lattice["dz"]))
-    return d
-
-
-# ======================================================================================
-# Lattice operations
-# ======================================================================================
-
-
-@njit
-def lattice_get_indices(x, y, z, ux, uy, uz, lattice):
-    # TODO: docs
-    ix = lattice_get_grid_index(x, ux, lattice["x0"], lattice["dx"])
-    iy = lattice_get_grid_index(y, uy, lattice["y0"], lattice["dy"])
-    iz = lattice_get_grid_index(z, uz, lattice["z0"], lattice["dz"])
-    return ix, iy, iz
-
-
-@njit
-def lattice_grid_distance(value, direction, x0, dx):
-    # TODO: docs
-    if direction == 0.0:
-        return INF
-    idx = math.floor((value - x0) / dx)
-    if direction > 0.0:
-        idx += 1
-    ref = x0 + idx * dx
-    dist = (ref - value) / direction
-    return dist
-
-
-@njit
-def lattice_get_grid_index(value, direction, x0, dx):
-    idx = int64(math.floor((x - lattice["x0"]) / lattice["dx"]))
-
-    # Coinciding cases
-    if direction > 0.0:
-        if abs(grid[idx + 1] - value) < COINCIDENCE_TOLERANCE:
-            idx += 1
-    else:
-        if abs(grid[idx] - value) < COINCIDENCE_TOLERANCE:
-            idx -= 1
-    return idx
-
-
 # =============================================================================
 # Surface operations
 # =============================================================================
@@ -360,10 +309,15 @@ def surface_distance(particle, surface, mcdc):
     uy = particle["uy"]
     uz = particle["uz"]
 
-    # Check if coincident
+    # Check if coincident and leaving the surface
     evaluation = surface_evaluate(particle, surface)
+    coincident = False
     if abs(evaluation) < COINCIDENCE_TOLERANCE:
-        return INF
+        coincident = True
+        if surface["linear"]:
+            return INF
+        elif surface_normal_component(particle, surface) > 0.0:
+            return INF
 
     # Surface coefficients
     G = surface["G"]
@@ -423,6 +377,10 @@ def surface_distance(particle, surface, mcdc):
         root_1 = (-b + sqrt) / denom
         root_2 = (-b - sqrt) / denom
 
+        # Coincident treatment
+        if coincident:
+            return max(root_1, root_2)
+
         # Negative roots, moving away from the surface
         if root_1 < 0.0:
             root_1 = INF
diff --git a/mcdc/loop.py b/mcdc/loop.py
index 19de82a3..5f47d10c 100644
--- a/mcdc/loop.py
+++ b/mcdc/loop.py
@@ -431,10 +431,7 @@ def loop_particle(P, prog):
     mcdc = adapt.device(prog)
 
     while P["alive"]:
-        print('start', P['x'], P['y'], P['z'], P['t'])
         step_particle(P, prog)
-        print('end', P['x'], P['y'], P['z'], P['t'])
-        input()
 
 
 @njit(cache=caching)
@@ -445,6 +442,9 @@ def step_particle(P, prog):
     kernel.move_to_event(P, mcdc)
     event = P["event"]
 
+    if event == EVENT_LOST:
+        return
+
     # The & operator here is a bitwise and.
     # It is used to determine if an event type is part of the particle event.