diff --git a/Convective-Exit.jl b/Convective-Exit.jl
deleted file mode 100644
index 9207ff6c..00000000
--- a/Convective-Exit.jl
+++ /dev/null
@@ -1,35 +0,0 @@
-using WaterLily,SpecialFunctions,ForwardDiff,StaticArrays
-
-function lamb_dipole(N;D=N/3,U=1,exit=false)
-    β = 2.4394π/D
-    @show besselj1(β*D/2)
-    C = -2U/(β*besselj0(β*D/2))
-    function ψ(x,y)
-        r = √(x^2+y^2)
-        ifelse(r ≥ D/2, U*((D/2r)^2-1)*y, C*besselj1(β*r)*y/r)
-    end
-    center = SA[N/2,N/2]
-    function uλ(i,xy)
-        x,y = xy-center
-        ifelse(i==1,ForwardDiff.derivative(y->ψ(x,y),y)+1+U,-ForwardDiff.derivative(x->ψ(x,y),x))
-    end
-    Simulation((N, N), (1,0), D; uλ, exit)
-end
-
-include("examples/TwoD_plots.jl")
-
-sim = lamb_dipole(64,U=0.25);
-@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
-flood(sim.flow.σ,clims=(-5,5))
-sim_step!(sim,1.2,remeasure=false)
-@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
-flood(sim.flow.σ,clims=(-5,5))
-
-sim = lamb_dipole(64,U=0.25,exit=true);
-sim_step!(sim,1.2,remeasure=false)
-@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
-flood(sim.flow.σ,clims=(-5,5))
-
-U=2
-sim = lamb_dipole(64,U=U,exit=true);
-sim_gif!(sim,duration=2.5/(1+U),step=0.05,clims=(-20,20))
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 61ffe13e..725575ee 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,7 +11,6 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/runtests.jl b/test/runtests.jl
index ac4c6c9c..be7aaae5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -32,7 +32,8 @@ arrays = setup_backends()
     @test ex == :(a[I, i] = Math.add(b[I], func(I, q)))
     @test sym == [:a, :I, :i, :(p.b), :q]
 
-    for f ∈ arrays
+    # for f ∈ arrays
+    for f ∈ [Array]
         p = Float64[i+j  for i ∈ 1:4, j ∈ 1:5] |> f
         @test inside(p) == CartesianIndices((2:3,2:4))
         @test L₂(p) == 187
@@ -46,14 +47,19 @@ arrays = setup_backends()
         σ = rand(Ng...) |> f # scalar
         BC!(u, U)
         BC!(σ)
-        @allowscalar() do
-            @test all(u[1, :, 1] .== U[1]) && all(u[2, :, 1] .== U[1]) && all(u[end, :, 1] .== U[1]) &&
+        @allowscalar @test all(u[1, :, 1] .== U[1]) && all(u[2, :, 1] .== U[1]) && all(u[end, :, 1] .== U[1]) &&
                 all(u[3:end-1, 1, 1] .== u[3:end-1, 2, 1]) && all(u[3:end-1, end, 1] .== u[3:end-1, end-1, 1])
-            @test all(u[:, 1, 2] .== U[2]) && all(u[:, 2, 2] .== U[2]) && all(u[:, end, 2] .== U[2]) &&
+        @allowscalar @test all(u[:, 1, 2] .== U[2]) && all(u[:, 2, 2] .== U[2]) && all(u[:, end, 2] .== U[2]) &&
                 all(u[1, 3:end-1, 2] .== u[2, 3:end-1, 2]) && all(u[end, 3:end-1, 2] .== u[end-1, 3:end-1, 2])
-            @test all(σ[1, 2:end-1] .== σ[2, 2:end-1]) && all(σ[end, 2:end-1] .== σ[end-1, 2:end-1]) &&
+        @allowscalar @test all(σ[1, 2:end-1] .== σ[2, 2:end-1]) && all(σ[end, 2:end-1] .== σ[end-1, 2:end-1]) &&
                 all(σ[2:end-1, 1] .== σ[2:end-1, 2]) && all(σ[2:end-1, end] .== σ[2:end-1, end-1])
-        end
+
+        @allowscalar u[end,:,1] .= 3
+        BC!(u, U, true) # save exit values
+        @allowscalar @test all(u[end, :, 1] .== 3)
+
+        WaterLily.exitBC!(u,u,U,0) # conservative exit check
+        @allowscalar @test all(u[end,2:end-1, 1] .== U[1])
     end
 end
 
@@ -205,13 +211,13 @@ end
     end
 end
 
-function sphere_sim(radius = 8; mem=Array)
-    body = AutoBody((x,t)-> √sum(abs2,x .- 2radius) - radius)
-    return Simulation(radius.*(6,4),(1,0),radius; body, ν=radius/250, T=Float32, mem)
+function sphere_sim(radius = 8; mem=Array, exit=false)
+    body = AutoBody((x,t)-> √sum(abs2,x .- (2radius+1.5)) - radius)
+    return Simulation(radius.*(6,4),(1,0),radius; body, ν=radius/250, T=Float32, mem, exit)
 end
 @testset "WaterLily.jl" begin
-    for mem ∈ arrays
-        sim = sphere_sim(32;mem);
+    for mem ∈ arrays, exit ∈ (true,false)
+        sim = sphere_sim(;mem,exit);
         @test sim_time(sim) == 0
         sim_step!(sim,0.1,remeasure=false)
         @test length(sim.flow.Δt)-1 == length(sim.pois.n)÷2