split example visualizaiton pipeline

tests/test_masked_poisson_capfft.py

@@ -1,7 +1,9 @@
 from pathlib import Path
 
 import numpy as np
+import pytest
 
+from orisys.backend_availability import CUPY_AVAILABLE
 from orisys.hhd import PoissonSolver, nHHD
 
 
@@ -124,6 +126,46 @@ def test_capfft_hhd_forces_track_sparse_reference():
     assert np.allclose(vshear[~mask], 0.0, atol=1e-5)
 
 
+def test_sparse_hhd_normal_sign_tracks_capfft_on_masked_uv_flow():
+    shape = (36, 44)
+    mask = _circle_mask(shape)
+    flow = _radial_flow(shape, scale=0.03)
+    flow[~mask] = 0.0
+
+    sparse_hhd = nHHD(
+        cuda=False,
+        grid=shape,
+        spacings=(1.0, 1.0),
+        poisson_backend="cpu",
+        masked_poisson_solver="sparse",
+    )
+    capfft_hhd = nHHD(
+        cuda=False,
+        grid=shape,
+        spacings=(1.0, 1.0),
+        poisson_backend="cpu",
+        masked_poisson_solver="capfft",
+        masked_poisson_max_boundary=40,
+        masked_poisson_max_iters=3,
+    )
+
+    sparse_hhd.get_force(flow, valid_mask=mask)
+    capfft_hhd.get_force(flow, valid_mask=mask)
+
+    sparse_vnormal = sparse_hhd.get_normal()
+    capfft_vnormal = capfft_hhd.get_normal()
+
+    # The masked sparse solver should use the same sign convention as capfft/FFT.
+    sparse_flat = sparse_vnormal[mask].reshape(-1, 2)
+    capfft_flat = capfft_vnormal[mask].reshape(-1, 2)
+    dot = float(np.sum(sparse_flat * capfft_flat))
+    assert dot > 0.0
+
+    denom = max(1.0, float(np.linalg.norm(capfft_flat)))
+    rel = float(np.linalg.norm(sparse_flat - capfft_flat) / denom)
+    assert rel < 0.35
+
+
 def test_capfft_uv_mask_fixture_improves_over_raw_fft():
     uv_path = REPO_ROOT / "unroll_map_uv.npz"
     if not uv_path.exists():
@@ -160,3 +202,91 @@ def test_capfft_uv_mask_fixture_improves_over_raw_fft():
     err_cap = float(np.linalg.norm((capfft - sparse)[mask]))
     err_raw = float(np.linalg.norm((raw_fft - sparse)[mask]))
     assert err_cap < err_raw
+
+
+@pytest.mark.skipif(not CUPY_AVAILABLE, reason="CuPy not available")
+def test_capfft_solver_cupy_backend_uses_gpu_bulk_solve_and_tracks_cpu():
+    shape = (32, 40)
+    mask = _circle_mask(shape)
+    rhs = _smooth_rhs(shape)
+    rhs_masked = np.where(mask, rhs, 0).astype(np.float32)
+
+    cpu_solver = PoissonSolver(
+        grid=shape,
+        spacings=(1.0, 1.0),
+        poisson_backend="cpu",
+        masked_poisson_solver="capfft",
+        masked_poisson_max_boundary=40,
+        masked_poisson_max_iters=3,
+    )
+    cupy_solver = PoissonSolver(
+        grid=shape,
+        spacings=(1.0, 1.0),
+        poisson_backend="cupy",
+        masked_poisson_solver="capfft",
+        masked_poisson_max_boundary=40,
+        masked_poisson_max_iters=3,
+    )
+
+    cpu_capfft = cpu_solver.solve(rhs_masked, mask=mask)
+    raw_fft = cpu_solver.solve(rhs_masked, mask=None)
+
+    original_gpu_solve = cupy_solver._rect_solver._solve_gpu_internal
+    calls = {"count": 0}
+
+    def counted_gpu_solve(f_gpu):
+        calls["count"] += 1
+        return original_gpu_solve(f_gpu)
+
+    cupy_solver._rect_solver._solve_gpu_internal = counted_gpu_solve
+    cupy_capfft = cupy_solver.solve(rhs_masked, mask=mask)
+
+    assert calls["count"] > 0
+    assert np.allclose(cupy_capfft[~mask], 0.0, atol=1e-6)
+
+    err_gpu = float(np.linalg.norm((cupy_capfft - cpu_capfft)[mask]))
+    err_raw = float(np.linalg.norm((raw_fft - cpu_capfft)[mask]))
+    assert err_gpu < err_raw
+
+    denom = max(1.0, float(np.linalg.norm(cpu_capfft[mask])))
+    rel = err_gpu / denom
+    assert rel < 1e-2
+
+
+@pytest.mark.skipif(not CUPY_AVAILABLE, reason="CuPy not available")
+def test_capfft_hhd_cupy_backend_tracks_cpu_capfft_reference():
+    shape = (36, 44)
+    mask = _circle_mask(shape)
+    flow = _radial_flow(shape, scale=0.03)
+    flow[~mask] = 0.0
+
+    cpu_hhd = nHHD(
+        cuda=False,
+        grid=shape,
+        spacings=(1.0, 1.0),
+        poisson_backend="cpu",
+        masked_poisson_solver="capfft",
+        masked_poisson_max_boundary=40,
+        masked_poisson_max_iters=3,
+    )
+    cupy_hhd = nHHD(
+        cuda=False,
+        grid=shape,
+        spacings=(1.0, 1.0),
+        poisson_backend="cupy",
+        masked_poisson_solver="capfft",
+        masked_poisson_max_boundary=40,
+        masked_poisson_max_iters=3,
+    )
+
+    _, cpu_forces = cpu_hhd.get_force(flow, valid_mask=mask)
+    _, cupy_forces = cupy_hhd.get_force(flow, valid_mask=mask)
+
+    denom = max(1.0, float(np.linalg.norm(cpu_forces)))
+    rel = float(np.linalg.norm(cupy_forces - cpu_forces) / denom)
+    assert rel < 1e-2
+
+    vnormal = cupy_hhd.get_normal()
+    vshear = cupy_hhd.get_shear()
+    assert np.allclose(vnormal[~mask], 0.0, atol=1e-5)
+    assert np.allclose(vshear[~mask], 0.0, atol=1e-5)