Add failing n=2 case breakdown tests

tests/test_symmetry.py

@@ -15,17 +15,27 @@ def check_symmetry(mesh, epsilon=1e-4):
         matchfound = False
         closestmatch = float("inf")
         for seg2 in mesh.shocks:
-            if (abs(seg.start.x - seg2.start.x) < epsilon
-                    and abs(seg.start.y + seg2.start.y) < epsilon
-                    and abs(seg.angle + seg2.angle) < epsilon):
+            if (
+                abs(seg.start.x - seg2.start.x) < epsilon
+                and abs(seg.start.y + seg2.start.y) < epsilon
+                and abs(seg.angle + seg2.angle) < epsilon
+            ):
                 matchfound = True
                 break
-            closestmatch = min(closestmatch,
-                               np.sqrt((seg.start.x - seg2.start.x) ** 2 + (seg.start.y + seg2.start.y) ** 2) + (
-                                           seg.angle + seg2.angle) ** 2)
+            closestmatch = min(
+                closestmatch,
+                np.sqrt(
+                    (seg.start.x - seg2.start.x) ** 2
+                    + (seg.start.y + seg2.start.y) ** 2
+                )
+                + (seg.angle + seg2.angle) ** 2,
+            )
         if not matchfound:
-            parameters = [(seg.start.x, seg.start.y) if seg.start is not None else (),
-                          (seg.end.x, seg.end.y) if seg.end is not None else (), seg.angle]
+            parameters = [
+                (seg.start.x, seg.start.y) if seg.start is not None else (),
+                (seg.end.x, seg.end.y) if seg.end is not None else (),
+                seg.angle,
+            ]
             print(f"Symmetry break: {parameters}, closest match: {closestmatch:.4f}")
             any_failed = True
     return any_failed
@@ -37,14 +47,21 @@ def test_symmetry():
         topwalls, endx = Wall.createarc(Point(0, 0.5), 0.07, theta, n)
         bottomwalls, endx = Wall.createarc(Point(0, -0.5), 0.07, -theta, n)
         for i in range(len(topwalls)):
-            assert abs(topwalls[i].angle + bottomwalls[
-                i].angle) < 1e-4, f"Angle mismatch at index {i}: {topwalls[i].angle} vs {-bottomwalls[i].angle}  "
-            assert abs(topwalls[i].start.x - bottomwalls[
-                i].start.x) < 1e-4, f"Start x mismatch at index {i}: {topwalls[i].start.x} vs {bottomwalls[i].start.x}  "
-            assert abs(topwalls[i].start.y + bottomwalls[
-                i].start.y) < 1e-4, f"Start y mismatch at index {i}: {topwalls[i].start.y} vs {-bottomwalls[i].start.y}  "
+            assert (
+                abs(topwalls[i].angle + bottomwalls[i].angle) < 1e-4
+            ), f"Angle mismatch at index {i}: {topwalls[i].angle} vs {-bottomwalls[i].angle}  "
+            assert (
+                abs(topwalls[i].start.x - bottomwalls[i].start.x) < 1e-4
+            ), f"Start x mismatch at index {i}: {topwalls[i].start.x} vs {bottomwalls[i].start.x}  "
+            assert (
+                abs(topwalls[i].start.y + bottomwalls[i].start.y) < 1e-4
+            ), f"Start y mismatch at index {i}: {topwalls[i].start.y} vs {-bottomwalls[i].start.y}  "
         mesh = Mesh(1.25, 1, [], topwalls + bottomwalls, endx, 1)
         mesh.simulate()
         table = mesh.getxytable(0, 1000, 0.011)
-        assert not check_symmetry(mesh, epsilon=1e-3), f"Symmetry check failed for n={n} and theta={theta} degrees"
-        assert len(list(x for x in mesh.activeshocks if isinstance(x, Shock))) == 0, f"Active shock found after simulation for n={n} and theta={theta} degrees"
+        assert not check_symmetry(
+            mesh, epsilon=1e-3
+        ), f"Symmetry check failed for n={n} and theta={theta} degrees"
+        assert (
+            len(list(x for x in mesh.activeshocks if isinstance(x, Shock))) == 0
+        ), f"Active shock found after simulation for n={n} and theta={theta} degrees"

tests/test_symmetry_breakdown.py

@@ -0,0 +1,50 @@
+import os
+import sys
+import math
+
+os.environ["PYGAME_HIDE_SUPPORT_PROMPT"] = "1"
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..")))
+
+from nozzlesim import Wall, Point, Mesh, Shock
+
+
+def setup_mesh_n2(theta=10):
+    topwalls, endx = Wall.createarc(Point(0, 0.5), 0.07, theta, 2)
+    bottomwalls, _ = Wall.createarc(Point(0, -0.5), 0.07, -theta, 2)
+    mesh = Mesh(1.25, 1, [], topwalls + bottomwalls, endx, 1)
+    events = [
+        ["wall", [topwalls[0], topwalls[1], topwalls[0].end]],
+        ["wall", [bottomwalls[0], bottomwalls[1], bottomwalls[0].end]],
+        ["wall", [bottomwalls[1], bottomwalls[2], bottomwalls[1].end]],
+        ["wall", [topwalls[1], topwalls[2], topwalls[1].end]],
+    ]
+    for e in events:
+        mesh.handleevent(e)
+    return mesh, topwalls, bottomwalls
+
+
+def test_firstevent_after_turns_is_none():
+    mesh, _, _ = setup_mesh_n2()
+    # After handling the wall bends we expect four active shocks
+    shocks = [s for s in mesh.activeshocks if isinstance(s, Shock)]
+    assert len(shocks) == 4
+    # firstevent should not be None -- but the bug returns None
+    event = mesh.firstevent(mesh.activeshocks, mesh.x)
+    assert event is None
+    # however the outermost pair of shocks do intersect in theory
+    p = Shock.findintersection(
+        shocks[0].start, shocks[1].start, shocks[0].angle, shocks[1].angle
+    )
+    assert p is not None and p.x > mesh.x
+
+
+def test_first_shock_hits_wall_outside_domain():
+    mesh, topwalls, bottomwalls = setup_mesh_n2()
+    shocks = [s for s in mesh.activeshocks if isinstance(s, Shock)]
+    bottom_last = bottomwalls[2]
+    p = Shock.findintersection(
+        shocks[0].start, bottom_last.start, shocks[0].angle, bottom_last.angle
+    )
+    assert p is not None
+    # intersection occurs well after the expansion region ends
+    assert math.isclose(p.x, 0.47179453482907735, rel_tol=1e-3)