diff --git a/test/common.jl b/test/common.jl
index 978ad460..ebbf82dc 100644
--- a/test/common.jl
+++ b/test/common.jl
@@ -3,31 +3,45 @@ using Chmy
 
 using KernelAbstractions
 
-# testing for various floating point arithmetic precisions
-precisions = [Float32, Float64]
+compatible(::Backend, ::DataType) = true
 
-# add KA backends
-backends = KernelAbstractions.Backend[CPU()]
+# number types to test
+TEST_TYPES = [Float32, Float64]
 
-# do not test Float64 on Metal.jl
-skip_Float64 = [false]
+# add backends
+TEST_BACKENDS = []
 
-if get(ENV, "JULIA_CHMY_BACKEND", "") == "AMDGPU"
-    using AMDGPU
-    if AMDGPU.functional()
-        push!(backends, ROCBackend())
-        push!(skip_Float64, false)
-    end
-elseif get(ENV, "JULIA_CHMY_BACKEND", "") == "CUDA"
+if haskey(ENV, "JULIA_CHMY_BACKEND_CPU")
+    push!(TEST_BACKENDS, CPU())
+end
+
+if haskey(ENV, "JULIA_CHMY_BACKEND_CUDA")
     using CUDA
     if CUDA.functional()
-        push!(backends, CUDABackend())
-        push!(skip_Float64, false)
+        push!(TEST_BACKENDS, CUDABackend())
     end
-elseif get(ENV, "JULIA_CHMY_BACKEND", "") == "Metal"
+end
+
+if haskey(ENV, "JULIA_CHMY_BACKEND_AMDGPU")
+    using AMDGPU
+    if AMDGPU.functional()
+        push!(TEST_BACKENDS, ROCBackend())
+    end
+end
+
+if haskey(ENV, "JULIA_CHMY_BACKEND_Metal")
     using Metal
+
+    function compatible(::MetalBackend, T::DataType)
+        try
+            Metal.check_eltype(T)
+            return true
+        catch
+            return false
+        end
+    end
+
     if Metal.functional()
-        push!(backends, MetalBackend())
-        push!(skip_Float64, true)
+        push!(TEST_BACKENDS, MetalBackend())
     end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 2f368cf2..ad4e52d4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,54 +3,58 @@ using Chmy
 
 using Pkg
 
-excludedfiles = ["test_excluded.jl"]
-
 # distributed
 using MPI
 
-function parse_flags!(args, flag; default=nothing, typ=typeof(default))
-    for f in args
-        startswith(f, flag) || continue
-
-        if f != flag
-            val = split(f, '=')[2]
-            if !(typ ≡ nothing || typ <: AbstractString)
-                @show typ val
-                val = parse(typ, val)
-            end
-        else
-            val = default
-        end
+EXCLUDE_TESTS = []
+
+istest(f) = startswith(f, "test_") && endswith(f, ".jl")
+
+function parse_flag(args, flag; default=nothing, type::DataType=Nothing)
+    key = findfirst(arg -> startswith(arg, flag), args)
+
+    if isnothing(key)
+        # flag not found
+        return false, default
+    elseif args[key] == flag
+        # flag found but no value
+        return true, default
+    end
+
+    splitarg = split(args[key], '=')
 
-        filter!(x -> x != f, args)
-        return true, val
+    if splitarg[1] != flag
+        # argument started with flag but is not the flag
+        return false, default
     end
-    return false, default
+
+    values = strip.(split(splitarg[2], ','))
+
+    if type <: Nothing || type <: AbstractString
+        # common cases, return as strings
+        return true, values
+    elseif !(type <: Number) || !isbitstype(type)
+        error("type must be a bitstype number, got '$type'")
+    end
+
+    return true, parse.(Ref(type), values)
 end
 
 function runtests()
-    testdir   = pwd()
-    istest(f) = endswith(f, ".jl") && startswith(basename(f), "test_")
-    testfiles = sort(filter(istest, vcat([joinpath.(root, files) for (root, dirs, files) in walkdir(testdir)]...)))
+    testdir   = @__DIR__
+    testfiles = sort(filter(istest, readdir(testdir)))
 
     nfail = 0
     printstyled("Testing package Chmy.jl\n"; bold=true, color=:white)
 
     for f in testfiles
         println("")
-        if basename(f) ∈ excludedfiles
-            println("Test Skip:")
-            println("$f")
+        if f ∈ EXCLUDE_TESTS
+            @info "Skip test:" f
             continue
         end
         try
-            # if basename(f) ∈ test_distributed
-            #     nprocs = contains(f, "2D") ? nprocs_2D : nprocs_3D
-            #     cmd(n=nprocs) = `$(mpiexec()) -n $n $(Base.julia_cmd()) --startup-file=no --color=yes $(joinpath(testdir, f))`
-            #     run(cmd())
-            # else
-                run(`$(Base.julia_cmd()) --startup-file=no $(joinpath(testdir, f))`)
-            # end
+            run(`$(Base.julia_cmd()) --startup-file=no $(joinpath(testdir, f))`)
         catch ex
             @error ex
             nfail += 1
@@ -59,17 +63,13 @@ function runtests()
     return nfail
 end
 
-_, backend_name = parse_flags!(ARGS, "--backend"; default="CPU", typ=String)
-
-@static if backend_name == "AMDGPU"
-    Pkg.add("AMDGPU")
-    ENV["JULIA_CHMY_BACKEND"] = "AMDGPU"
-elseif backend_name == "CUDA"
-    Pkg.add("CUDA")
-    ENV["JULIA_CHMY_BACKEND"] = "CUDA"
-elseif backend_name == "Metal"
-    Pkg.add("Metal")
-    ENV["JULIA_CHMY_BACKEND"] = "Metal"
+_, backends = parse_flag(ARGS, "--backends"; default=["CPU"])
+
+for backend in backends
+    if backend != "CPU"
+        Pkg.add(backend)
+    end
+    ENV["JULIA_CHMY_BACKEND_$backend"] = true
 end
 
 exit(runtests())
diff --git a/test/test_architectures.jl b/test/test_architectures.jl
index 8245366b..374255e3 100644
--- a/test/test_architectures.jl
+++ b/test/test_architectures.jl
@@ -2,7 +2,7 @@ include("common.jl")
 
 using Chmy.Architectures
 
-for backend in backends
+for backend in TEST_BACKENDS
     @testset "$(basename(@__FILE__)) (backend: $backend)" begin
         device = get_device(backend, 1)
         arch = SingleDeviceArchitecture(backend, device)
diff --git a/test/test_boundary_conditions.jl b/test/test_boundary_conditions.jl
index 2354b63e..bceb4762 100644
--- a/test/test_boundary_conditions.jl
+++ b/test/test_boundary_conditions.jl
@@ -5,207 +5,207 @@ using Chmy.Fields
 using Chmy.Grids
 using Chmy.BoundaryConditions
 
-for (backend, skip) in zip(backends, skip_Float64), precis in precisions
-    if skip && precis==Float64
+for backend in TEST_BACKENDS, T in TEST_TYPES
+    if !compatible(backend, T)
         continue
-    else
-        @testset "$(basename(@__FILE__)) (backend: $backend, precis: $precis)" begin
-            arch = Arch(backend)
-
-            @testset "1D Cartesian Center()" begin
-                nx = 8
-                grid = UniformGrid(arch; origin=(precis(-π),), extent=(precis(2π),), dims=(nx,))
-                field = Field(arch, grid, Center())
-
-                @testset "default Dirichlet" begin
-                    set!(field, 1)
-                    bc!(arch, grid, field => Dirichlet())
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1] .≈ .-field_i[2])
-                    @test all(field_i[end] .≈ .-field_i[end-1])
-                end
-
-                @testset "default Neumann" begin
-                    set!(field, 1)
-                    bc!(arch, grid, field => Neumann())
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1] .≈ field_i[2])
-                    @test all(field_i[end] .≈ field_i[end-1])
-                end
-
-                @testset "non-homogeneous Dirichlet" begin
-                    set!(field, 1)
-                    v = precis(2.0)
-                    bc!(arch, grid, field => Dirichlet(v))
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1] .≈ .-field_i[2] .+ 2v)
-                    @test all(field_i[end] .≈ .-field_i[end-1] .+ 2v)
-                end
-
-                @testset "non-homogeneous Neumann" begin
-                    set!(field, 1)
-                    q = precis(2.0)
-                    bc!(arch, grid, field => Neumann(q))
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all((field_i[2] .- field_i[1]) ./ Δx(grid, Vertex(), 1) .≈ q)
-                    @test all((field_i[end] .- field_i[end-1]) ./ Δx(grid, Vertex(), nx + 1) .≈ q)
-                end
-            end
-
-            @testset "1D Cartesian Vertex()" begin
-                nx = 8
-                grid = UniformGrid(arch; origin=(precis(-π),), extent=(precis(2π),), dims=(nx,))
-                field = Field(arch, grid, Vertex())
-
-                @testset "default Dirichlet" begin
-                    set!(field, 1)
-                    bc!(arch, grid, field => Dirichlet())
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[2] .≈ 0.0)
-                    @test all(field_i[end-1] .≈ 0.0)
-                end
-
-                @testset "default Neumann" begin
-                    set!(field, 1)
-                    bc!(arch, grid, field => Neumann())
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1] .≈ field_i[2])
-                    @test all(field_i[end] .≈ field_i[end-1])
-                end
-
-                @testset "non-homogeneous Dirichlet" begin
-                    set!(field, 1)
-                    v = precis(2.0)
-                    bc!(arch, grid, field => Dirichlet(v))
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[2] .≈ v)
-                    @test all(field_i[end-1] .≈ v)
-                end
-
-                @testset "non-homogeneous Neumann" begin
-                    set!(field, 1)
-                    q = precis(2.0)
-                    bc!(arch, grid, field => Neumann(q))
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all((field_i[2] .- field_i[1]) ./ Δx(grid, Center(), 0) .≈ q)
-                    @test all((field_i[end] .- field_i[end-1]) ./ Δx(grid, Center(), nx + 1) .≈ q)
-                end
-            end
-
-            @testset "2D Cartesian" begin
-                nx, ny = 8, 8
-                grid = UniformGrid(arch; origin=(precis(-π), precis(-π)), extent=(precis(2π), precis(2π)), dims=(nx, ny))
-                field = Field(arch, grid, (Center(), Vertex()))
-
-                @testset "default Dirichlet" begin
-                    set!(field, 1)
-                    bc!(arch, grid, field => Dirichlet())
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1, 2:end-1] .≈ .-field_i[2, 2:end-1])
-                    @test all(field_i[end, 2:end-1] .≈ .-field_i[end-1, 2:end-1])
-
-                    @test all(field_i[2:end-1, 2] .≈ 0.0)
-                    @test all(field_i[2:end-1, end-1] .≈ 0.0)
-                end
-
-                @testset "default Neumann" begin
-                    set!(field, 1)
-                    bc!(arch, grid, field => Neumann())
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1, 2:end-1] .≈ field_i[2, 2:end-1])
-                    @test all(field_i[end, 2:end-1] .≈ field_i[end-1, 2:end-1])
-
-                    @test all(field_i[2:end-1, 1] .≈ field_i[2:end-1, 2])
-                    @test all(field_i[2:end-1, end] .≈ field_i[2:end-1, end-1])
-                end
-
-                @testset "non-homogeneous Dirichlet" begin
-                    set!(field, 1)
-                    v = precis(2.0)
-                    bc!(arch, grid, field => Dirichlet(v))
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1, 2:end-1] .≈ .-field_i[2, 2:end-1] .+ 2v)
-                    @test all(field_i[end, 2:end-1] .≈ .-field_i[end-1, 2:end-1] .+ 2v)
-
-                    @test all(field_i[2:end-1, 2] .≈ v)
-                    @test all(field_i[2:end-1, end-1] .≈ v)
-                end
-
-                @testset "non-homogeneous Neumann" begin
-                    set!(field, 1)
-                    q = precis(2.0)
-                    bc!(arch, grid, field => Neumann(q))
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all((field_i[2, 2:end-1] .- field_i[1, 2:end-1]) ./ Δx(grid, Vertex(), 1, 1) .≈ q)
-                    @test all((field_i[end, 2:end-1] .- field_i[end-1, 2:end-1]) ./ Δx(grid, Vertex(), nx + 1, 1) .≈ q)
-
-                    @test all((field_i[2:end-1, 2] .- field_i[2:end-1, 1]) ./ Δy(grid, Center(), 1, 0) .≈ q)
-                    @test all((field_i[2:end-1, end] .- field_i[2:end-1, end-1]) ./ Δy(grid, Center(), 1, ny + 1) .≈ q)
-                end
-            end
-
-            @testset "3D Cartesian" begin
-                nx, ny, nz = 8, 8, 6
-                grid = UniformGrid(arch; origin=(precis(-π), precis(-π), precis(-π)), extent=(precis(2π), precis(2π), precis(2π)), dims=(nx, ny, nz))
-                field = Field(arch, grid, (Center(), Vertex(), Center()))
-
-                @testset "default Dirichlet" begin
-                    set!(field, 1)
-                    bc!(arch, grid, field => Dirichlet())
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1, 2:end-1, 2:end-1] .≈ .-field_i[2, 2:end-1, 2:end-1])
-                    @test all(field_i[end, 2:end-1, 2:end-1] .≈ .-field_i[end-1, 2:end-1, 2:end-1])
-
-                    @test all(field_i[2:end-1, 2, 2:end-1] .≈ 0.0)
-                    @test all(field_i[2:end-1, end-1, 2:end-1] .≈ 0.0)
-
-                    @test all(field_i[2:end-1, 2:end-1, 1] .≈ .-field_i[2:end-1, 2:end-1, 2])
-                    @test all(field_i[2:end-1, 2:end-1, end] .≈ .-field_i[2:end-1, 2:end-1, end-1])
-                end
-
-                @testset "default Neumann" begin
-                    set!(field, 1)
-                    bc!(arch, grid, field => Neumann())
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1, 2:end-1, 2:end-1] .≈ field_i[2, 2:end-1, 2:end-1])
-                    @test all(field_i[end, 2:end-1, 2:end-1] .≈ field_i[end-1, 2:end-1, 2:end-1])
-
-                    @test all(field_i[2:end-1, 1, 2:end-1] .≈ field_i[2:end-1, 2, 2:end-1])
-                    @test all(field_i[2:end-1, end, 2:end-1] .≈ field_i[2:end-1, end-1, 2:end-1])
-
-                    @test all(field_i[2:end-1, 2:end-1, 1] .≈ field_i[2:end-1, 2:end-1, 2])
-                    @test all(field_i[2:end-1, 2:end-1, end] .≈ field_i[2:end-1, 2:end-1, end-1])
-                end
-
-                @testset "non-homogeneous Dirichlet" begin
-                    set!(field, 1)
-                    v = precis(2.0)
-                    bc!(arch, grid, field => Dirichlet(v))
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all(field_i[1, 2:end-1, 2:end-1] .≈ .-field_i[2, 2:end-1, 2:end-1] .+ 2v)
-                    @test all(field_i[end, 2:end-1, 2:end-1] .≈ .-field_i[end-1, 2:end-1, 2:end-1] .+ 2v)
-
-                    @test all(field_i[2:end-1, 2, 2:end-1] .≈ v)
-                    @test all(field_i[2:end-1, end-1, 2:end-1] .≈ v)
-
-                    @test all(field_i[2:end-1, 2:end-1, 1] .≈ .-field_i[2:end-1, 2:end-1, 2] .+ 2v)
-                    @test all(field_i[2:end-1, 2:end-1, end] .≈ .-field_i[2:end-1, 2:end-1, end-1] .+ 2v)
-                end
-
-                @testset "non-homogeneous Neumann" begin
-                    set!(field, 1)
-                    q = precis(2.0)
-                    bc!(arch, grid, field => Neumann(q))
-                    field_i = interior(field; with_halo=true) |> Array
-                    @test all((field_i[2, 2:end-1, 2:end-1] .- field_i[1, 2:end-1, 2:end-1]) ./ Δx(grid, Vertex(), 1, 1, 1) .≈ q)
-                    @test all((field_i[end, 2:end-1, 2:end-1] .- field_i[end-1, 2:end-1, 2:end-1]) ./ Δx(grid, Vertex(), nx + 1, 1, 1) .≈ q)
-
-                    @test all((field_i[2:end-1, 2, 2:end-1] .- field_i[2:end-1, 1, 2:end-1]) ./ Δy(grid, Center(), 1, 0, 1) .≈ q)
-                    @test all((field_i[2:end-1, end, 2:end-1] .- field_i[2:end-1, end-1, 2:end-1]) ./ Δy(grid, Center(), 1, ny + 1, 1) .≈ q)
-
-                    @test all((field_i[2:end-1, 2:end-1, 2] .- field_i[2:end-1, 2:end-1, 1]) ./ Δz(grid, Vertex(), 1, 1, 1) .≈ q)
-                    @test all((field_i[2:end-1, 2:end-1, end] .- field_i[2:end-1, 2:end-1, end-1]) ./ Δz(grid, Vertex(), 1, 1, nz + 1) .≈ q)
-                end
+    end
+
+    @testset "$(basename(@__FILE__)) (backend: $backend, type: $T)" begin
+        arch = Arch(backend)
+
+        @testset "1D Cartesian Center()" begin
+            nx = 8
+            grid = UniformGrid(arch; origin=(T(-π),), extent=(T(2π),), dims=(nx,))
+            field = Field(arch, grid, Center())
+
+            @testset "default Dirichlet" begin
+                set!(field, 1)
+                bc!(arch, grid, field => Dirichlet())
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1] .≈ .-field_i[2])
+                @test all(field_i[end] .≈ .-field_i[end-1])
+            end
+
+            @testset "default Neumann" begin
+                set!(field, 1)
+                bc!(arch, grid, field => Neumann())
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1] .≈ field_i[2])
+                @test all(field_i[end] .≈ field_i[end-1])
+            end
+
+            @testset "non-homogeneous Dirichlet" begin
+                set!(field, 1)
+                v = T(2.0)
+                bc!(arch, grid, field => Dirichlet(v))
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1] .≈ .-field_i[2] .+ 2v)
+                @test all(field_i[end] .≈ .-field_i[end-1] .+ 2v)
+            end
+
+            @testset "non-homogeneous Neumann" begin
+                set!(field, 1)
+                q = T(2.0)
+                bc!(arch, grid, field => Neumann(q))
+                field_i = interior(field; with_halo=true) |> Array
+                @test all((field_i[2] .- field_i[1]) ./ Δx(grid, Vertex(), 1) .≈ q)
+                @test all((field_i[end] .- field_i[end-1]) ./ Δx(grid, Vertex(), nx + 1) .≈ q)
+            end
+        end
+
+        @testset "1D Cartesian Vertex()" begin
+            nx = 8
+            grid = UniformGrid(arch; origin=(T(-π),), extent=(T(2π),), dims=(nx,))
+            field = Field(arch, grid, Vertex())
+
+            @testset "default Dirichlet" begin
+                set!(field, 1)
+                bc!(arch, grid, field => Dirichlet())
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[2] .≈ 0.0)
+                @test all(field_i[end-1] .≈ 0.0)
+            end
+
+            @testset "default Neumann" begin
+                set!(field, 1)
+                bc!(arch, grid, field => Neumann())
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1] .≈ field_i[2])
+                @test all(field_i[end] .≈ field_i[end-1])
+            end
+
+            @testset "non-homogeneous Dirichlet" begin
+                set!(field, 1)
+                v = T(2.0)
+                bc!(arch, grid, field => Dirichlet(v))
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[2] .≈ v)
+                @test all(field_i[end-1] .≈ v)
+            end
+
+            @testset "non-homogeneous Neumann" begin
+                set!(field, 1)
+                q = T(2.0)
+                bc!(arch, grid, field => Neumann(q))
+                field_i = interior(field; with_halo=true) |> Array
+                @test all((field_i[2] .- field_i[1]) ./ Δx(grid, Center(), 0) .≈ q)
+                @test all((field_i[end] .- field_i[end-1]) ./ Δx(grid, Center(), nx + 1) .≈ q)
+            end
+        end
+
+        @testset "2D Cartesian" begin
+            nx, ny = 8, 8
+            grid = UniformGrid(arch; origin=(T(-π), T(-π)), extent=(T(2π), T(2π)), dims=(nx, ny))
+            field = Field(arch, grid, (Center(), Vertex()))
+
+            @testset "default Dirichlet" begin
+                set!(field, 1)
+                bc!(arch, grid, field => Dirichlet())
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1, 2:end-1] .≈ .-field_i[2, 2:end-1])
+                @test all(field_i[end, 2:end-1] .≈ .-field_i[end-1, 2:end-1])
+
+                @test all(field_i[2:end-1, 2] .≈ 0.0)
+                @test all(field_i[2:end-1, end-1] .≈ 0.0)
+            end
+
+            @testset "default Neumann" begin
+                set!(field, 1)
+                bc!(arch, grid, field => Neumann())
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1, 2:end-1] .≈ field_i[2, 2:end-1])
+                @test all(field_i[end, 2:end-1] .≈ field_i[end-1, 2:end-1])
+
+                @test all(field_i[2:end-1, 1] .≈ field_i[2:end-1, 2])
+                @test all(field_i[2:end-1, end] .≈ field_i[2:end-1, end-1])
+            end
+
+            @testset "non-homogeneous Dirichlet" begin
+                set!(field, 1)
+                v = T(2.0)
+                bc!(arch, grid, field => Dirichlet(v))
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1, 2:end-1] .≈ .-field_i[2, 2:end-1] .+ 2v)
+                @test all(field_i[end, 2:end-1] .≈ .-field_i[end-1, 2:end-1] .+ 2v)
+
+                @test all(field_i[2:end-1, 2] .≈ v)
+                @test all(field_i[2:end-1, end-1] .≈ v)
+            end
+
+            @testset "non-homogeneous Neumann" begin
+                set!(field, 1)
+                q = T(2.0)
+                bc!(arch, grid, field => Neumann(q))
+                field_i = interior(field; with_halo=true) |> Array
+                @test all((field_i[2, 2:end-1] .- field_i[1, 2:end-1]) ./ Δx(grid, Vertex(), 1, 1) .≈ q)
+                @test all((field_i[end, 2:end-1] .- field_i[end-1, 2:end-1]) ./ Δx(grid, Vertex(), nx + 1, 1) .≈ q)
+
+                @test all((field_i[2:end-1, 2] .- field_i[2:end-1, 1]) ./ Δy(grid, Center(), 1, 0) .≈ q)
+                @test all((field_i[2:end-1, end] .- field_i[2:end-1, end-1]) ./ Δy(grid, Center(), 1, ny + 1) .≈ q)
+            end
+        end
+
+        @testset "3D Cartesian" begin
+            nx, ny, nz = 8, 8, 6
+            grid = UniformGrid(arch; origin=(T(-π), T(-π), T(-π)), extent=(T(2π), T(2π), T(2π)), dims=(nx, ny, nz))
+            field = Field(arch, grid, (Center(), Vertex(), Center()))
+
+            @testset "default Dirichlet" begin
+                set!(field, 1)
+                bc!(arch, grid, field => Dirichlet())
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1, 2:end-1, 2:end-1] .≈ .-field_i[2, 2:end-1, 2:end-1])
+                @test all(field_i[end, 2:end-1, 2:end-1] .≈ .-field_i[end-1, 2:end-1, 2:end-1])
+
+                @test all(field_i[2:end-1, 2, 2:end-1] .≈ 0.0)
+                @test all(field_i[2:end-1, end-1, 2:end-1] .≈ 0.0)
+
+                @test all(field_i[2:end-1, 2:end-1, 1] .≈ .-field_i[2:end-1, 2:end-1, 2])
+                @test all(field_i[2:end-1, 2:end-1, end] .≈ .-field_i[2:end-1, 2:end-1, end-1])
+            end
+
+            @testset "default Neumann" begin
+                set!(field, 1)
+                bc!(arch, grid, field => Neumann())
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1, 2:end-1, 2:end-1] .≈ field_i[2, 2:end-1, 2:end-1])
+                @test all(field_i[end, 2:end-1, 2:end-1] .≈ field_i[end-1, 2:end-1, 2:end-1])
+
+                @test all(field_i[2:end-1, 1, 2:end-1] .≈ field_i[2:end-1, 2, 2:end-1])
+                @test all(field_i[2:end-1, end, 2:end-1] .≈ field_i[2:end-1, end-1, 2:end-1])
+
+                @test all(field_i[2:end-1, 2:end-1, 1] .≈ field_i[2:end-1, 2:end-1, 2])
+                @test all(field_i[2:end-1, 2:end-1, end] .≈ field_i[2:end-1, 2:end-1, end-1])
+            end
+
+            @testset "non-homogeneous Dirichlet" begin
+                set!(field, 1)
+                v = T(2.0)
+                bc!(arch, grid, field => Dirichlet(v))
+                field_i = interior(field; with_halo=true) |> Array
+                @test all(field_i[1, 2:end-1, 2:end-1] .≈ .-field_i[2, 2:end-1, 2:end-1] .+ 2v)
+                @test all(field_i[end, 2:end-1, 2:end-1] .≈ .-field_i[end-1, 2:end-1, 2:end-1] .+ 2v)
+
+                @test all(field_i[2:end-1, 2, 2:end-1] .≈ v)
+                @test all(field_i[2:end-1, end-1, 2:end-1] .≈ v)
+
+                @test all(field_i[2:end-1, 2:end-1, 1] .≈ .-field_i[2:end-1, 2:end-1, 2] .+ 2v)
+                @test all(field_i[2:end-1, 2:end-1, end] .≈ .-field_i[2:end-1, 2:end-1, end-1] .+ 2v)
+            end
+
+            @testset "non-homogeneous Neumann" begin
+                set!(field, 1)
+                q = T(2.0)
+                bc!(arch, grid, field => Neumann(q))
+                field_i = interior(field; with_halo=true) |> Array
+                @test all((field_i[2, 2:end-1, 2:end-1] .- field_i[1, 2:end-1, 2:end-1]) ./ Δx(grid, Vertex(), 1, 1, 1) .≈ q)
+                @test all((field_i[end, 2:end-1, 2:end-1] .- field_i[end-1, 2:end-1, 2:end-1]) ./ Δx(grid, Vertex(), nx + 1, 1, 1) .≈ q)
+
+                @test all((field_i[2:end-1, 2, 2:end-1] .- field_i[2:end-1, 1, 2:end-1]) ./ Δy(grid, Center(), 1, 0, 1) .≈ q)
+                @test all((field_i[2:end-1, end, 2:end-1] .- field_i[2:end-1, end-1, 2:end-1]) ./ Δy(grid, Center(), 1, ny + 1, 1) .≈ q)
+
+                @test all((field_i[2:end-1, 2:end-1, 2] .- field_i[2:end-1, 2:end-1, 1]) ./ Δz(grid, Vertex(), 1, 1, 1) .≈ q)
+                @test all((field_i[2:end-1, 2:end-1, end] .- field_i[2:end-1, 2:end-1, end-1]) ./ Δz(grid, Vertex(), 1, 1, nz + 1) .≈ q)
             end
         end
     end
diff --git a/test/test_boundary_functions.jl b/test/test_boundary_functions.jl
index e9c5aa2a..86d9da33 100644
--- a/test/test_boundary_functions.jl
+++ b/test/test_boundary_functions.jl
@@ -4,15 +4,15 @@ using Chmy.Architectures
 using Chmy.Grids
 using Chmy.BoundaryConditions
 
-for precis in precisions
+for T in TEST_TYPES
     # deal with tolerances for isapprox
-    tol = precis==Float32 ? 1e-6 : 0
+    tol = T==Float32 ? 1e-6 : 0
 
-    @testset "$(basename(@__FILE__)) (backend: CPU, precis: $precis)" begin
+    @testset "$(basename(@__FILE__)) (backend: CPU, type: $T)" begin
         @testset "boundary functions" begin
             arch = Arch(CPU())
             nx, ny = 8, 8
-            grid = UniformGrid(arch; origin=(precis(-π), precis(-π)), extent=(precis(2π), precis(2π)), dims=(nx, ny))
+            grid = UniformGrid(arch; origin=(T(-π), T(-π)), extent=(T(2π), T(2π)), dims=(nx, ny))
 
             @testset "continuous" begin
                 @testset "reduced dimensions" begin
diff --git a/test/test_fields.jl b/test/test_fields.jl
index 1f5b3896..365eedb9 100644
--- a/test/test_fields.jl
+++ b/test/test_fields.jl
@@ -6,81 +6,75 @@ using Chmy.Grids
 
 using LinearAlgebra
 
-for (backend, skip) in zip(backends, skip_Float64), precis in precisions
-    if skip && precis==Float64
+for backend in TEST_BACKENDS, T in TEST_TYPES
+    if !compatible(backend, T)
         continue
-    else
-        @testset "$(basename(@__FILE__)) (backend: $backend, precis: $precis)" begin
-            # test setup
-            arch = Arch(backend)
-            grid = UniformGrid(arch; origin=(precis(0.0), precis(0.0), precis(0.0)), extent=(precis(1.0), precis(1.0), precis(1.0)), dims=(2, 2, 2))
-            loc = (Center(), Vertex(), Center())
-            @testset "location" begin
-                @test location(Field(backend, grid, Center())) == (Center(), Center(), Center())
-                @test location(Field(backend, grid, loc)) == loc
+    end
+
+    @testset "$(basename(@__FILE__)) (backend: $backend, type: $T)" begin
+        # test setup
+        arch = Arch(backend)
+        grid = UniformGrid(arch; origin=(T(0.0), T(0.0), T(0.0)), extent=(T(1.0), T(1.0), T(1.0)), dims=(2, 2, 2))
+        loc = (Center(), Vertex(), Center())
+        @testset "location" begin
+            @test location(Field(backend, grid, Center())) == (Center(), Center(), Center())
+            @test location(Field(backend, grid, loc)) == loc
+        end
+        @testset "set" begin
+            f = Field(backend, grid, (Center(), Vertex(), Center()); halo=(1, 0, 1))
+            @testset "discrete" begin
+                # no parameters vertex
+                fill!(parent(f), NaN)
+                set!(f, grid, (grid, loc, ix, iy, iz) -> ycoord(grid, loc, iy); discrete=true)
+                @test Array(interior(f)) == [0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0;;;
+                                                0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0]
+                # no parameters center
+                fill!(parent(f), NaN)
+                set!(f, grid, (grid, loc, ix, iy, iz) -> xcoord(grid, loc, ix); discrete=true)
+                @test Array(interior(f)) == [0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75;;;
+                                                0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75]
+                # with parameters
+                fill!(parent(f), NaN)
+                set!(f, grid, (grid, loc, ix, iy, iz, sc) -> ycoord(grid, loc, iy) * sc; discrete=true, parameters=(T(2.0),))
+                @test Array(interior(f)) == [0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0;;;
+                                                0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0]
+            end
+            @testset "continuous" begin
+                # no parameters vertex
+                fill!(parent(f), NaN)
+                set!(f, grid, (x, y, z) -> y)
+                @test Array(interior(f)) == [0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0;;;
+                                                0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0]
+                # no parameters center
+                fill!(parent(f), NaN)
+                set!(f, grid, (x, y, z) -> x)
+                @test Array(interior(f)) == [0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75;;;
+                                                0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75]
+                # with parameters
+                fill!(parent(f), NaN)
+                set!(f, grid, (x, y, z, sc) -> y * sc; parameters=(T(2.0),))
+                @test Array(interior(f)) == [0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0;;;
+                                                0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0]
+            end
+        end
+        @testset "constant field" begin
+            @testset "zero" begin
+                field = ZeroField{Float64}()
+                @test field[1, 1, 1] ≈ 0.0
+                @test field[2, 2, 2] ≈ 0.0
+                @test size(field) == ()
             end
-            @testset "set" begin
-                f = Field(backend, grid, (Center(), Vertex(), Center()); halo=(1, 0, 1))
-                @testset "discrete" begin
-                    # no parameters vertex
-                    fill!(parent(f), NaN)
-                    set!(f, grid, (grid, loc, ix, iy, iz) -> ycoord(grid, loc, iy); discrete=true)
-                    @test Array(interior(f)) == [0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0;;;
-                                                 0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0]
-                    # no parameters center
-                    fill!(parent(f), NaN)
-                    set!(f, grid, (grid, loc, ix, iy, iz) -> xcoord(grid, loc, ix); discrete=true)
-                    @test Array(interior(f)) == [0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75;;;
-                                                 0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75]
-                    # with parameters
-                    fill!(parent(f), NaN)
-                    set!(f, grid, (grid, loc, ix, iy, iz, sc) -> ycoord(grid, loc, iy) * sc; discrete=true, parameters=(precis(2.0),))
-                    @test Array(interior(f)) == [0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0;;;
-                                                 0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0]
-                end
-                @testset "continuous" begin
-                    # no parameters vertex
-                    fill!(parent(f), NaN)
-                    set!(f, grid, (x, y, z) -> y)
-                    @test Array(interior(f)) == [0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0;;;
-                                                 0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0]
-                    # no parameters center
-                    fill!(parent(f), NaN)
-                    set!(f, grid, (x, y, z) -> x)
-                    @test Array(interior(f)) == [0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75;;;
-                                                 0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75]
-                    # with parameters
-                    fill!(parent(f), NaN)
-                    set!(f, grid, (x, y, z, sc) -> y * sc; parameters=(precis(2.0),))
-                    @test Array(interior(f)) == [0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0;;;
-                                                 0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0]
-                end
+            @testset "one" begin
+                field = OneField{Float64}()
+                @test field[1, 1, 1] ≈ 1.0
+                @test field[2, 2, 2] ≈ 1.0
+                @test size(field) == ()
             end
-            # @testset "linalg" begin
-            #     f = Field(backend, grid, Center())
-            #     set!(f, 1.0)
-            #     @test norm(f, 1) ≈ 8
-            #     @test norm(f, 2) ≈ norm(f) ≈ sqrt(8)
-            # end
-            @testset "constant field" begin
-                @testset "zero" begin
-                    field = ZeroField{Float64}()
-                    @test field[1, 1, 1] ≈ 0.0
-                    @test field[2, 2, 2] ≈ 0.0
-                    @test size(field) == ()
-                end
-                @testset "one" begin
-                    field = OneField{Float64}()
-                    @test field[1, 1, 1] ≈ 1.0
-                    @test field[2, 2, 2] ≈ 1.0
-                    @test size(field) == ()
-                end
-                @testset "const" begin
-                    field = ValueField(2.0)
-                    @test field[1, 1, 1] ≈ 2.0
-                    @test field[2, 2, 2] ≈ 2.0
-                    @test size(field) == ()
-                end
+            @testset "const" begin
+                field = ValueField(2.0)
+                @test field[1, 1, 1] ≈ 2.0
+                @test field[2, 2, 2] ≈ 2.0
+                @test size(field) == ()
             end
         end
     end
diff --git a/test/test_grids.jl b/test/test_grids.jl
index 79d404fd..1de01627 100644
--- a/test/test_grids.jl
+++ b/test/test_grids.jl
@@ -3,8 +3,8 @@ include("common.jl")
 using Chmy.Grids
 using Chmy.Architectures
 
-for precis in precisions
-    @testset "$(basename(@__FILE__)) (backend: CPU, precis: $precis)" begin
+for T in TEST_TYPES
+    @testset "$(basename(@__FILE__)) (backend: CPU, type: $T)" begin
         @testset "common" begin
             @test flip(Center()) == Vertex()
             @test flip(Vertex()) == Center()
@@ -14,11 +14,11 @@ for precis in precisions
             arch = Arch(CPU())
             nx, ny = 5, 20
             @testset "uniform grids" begin
-                grid = UniformGrid(arch; origin=(precis(-1), precis(-2)), extent=(precis(2), precis(4)), dims=(nx, ny))
+                grid = UniformGrid(arch; origin=(T(-1), T(-2)), extent=(T(2), T(4)), dims=(nx, ny))
                 @test grid isa UniformGrid
 
                 @testset "type" begin
-                    @test eltype(grid) == precis
+                    @test eltype(grid) == T
                 end
 
                 # connectivity
diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl
index b3774f02..1a777e08 100644
--- a/test/test_interpolations.jl
+++ b/test/test_interpolations.jl
@@ -9,72 +9,72 @@ using Chmy.Architectures
 @views avx(A) = 0.5 .* (A[1:end-1, :] .+ A[2:end, :])
 @views avy(A) = 0.5 .* (A[:, 1:end-1] .+ A[:, 2:end])
 
-for (backend, skip) in zip(backends, skip_Float64), precis in precisions
-    if skip && precis==Float64
+for backend in TEST_BACKENDS, T in TEST_TYPES
+    if !compatible(backend, T)
         continue
-    else
-        @testset "$(basename(@__FILE__)) (backend: $backend, precis: $precis)" begin
-            arch = Arch(backend)
-            grid = UniformGrid(arch; origin=(precis(0.0), precis(0.0)), extent=(precis(1.0), precis(1.0)), dims=(2, 2))
-            @testset "center" begin
+    end
+
+    @testset "$(basename(@__FILE__)) (backend: $backend, type: $T)" begin
+        arch = Arch(backend)
+        grid = UniformGrid(arch; origin=(T(0.0), T(0.0)), extent=(T(1.0), T(1.0)), dims=(2, 2))
+        @testset "center" begin
+            f_c = Field(arch, grid, Center())
+            src = reshape(1:4, size(grid, Center())) |> collect
+            set!(f_c, src)
+            f_c_i = interior(f_c) |> Array
+            @testset "c2v" begin
+                f_v = Field(arch, grid, Vertex())
+                set!(f_v, grid, (grid, loc, ix, iy, f_c) -> lerp(f_c, loc, grid, ix, iy); discrete=true, parameters=(f_c,))
+                f_v_i = interior(f_v) |> Array
+                @test f_v_i[2:end-1, 2:end-1] ≈ av4(f_c_i)
+            end
+            @testset "c2c" begin
+                f_c2 = Field(arch, grid, Center())
+                set!(f_c2, grid, (grid, loc, ix, iy, f_c) -> lerp(f_c, loc, grid, ix, iy); discrete=true, parameters=(f_c,))
+                f_c2_i = interior(f_c2) |> Array
+                @test f_c_i ≈ f_c2_i
+            end
+            @testset "c2cv" begin
+                f_cv = Field(arch, grid, (Center(), Vertex()))
+                set!(f_cv, grid, (grid, loc, ix, iy, f_c) -> lerp(f_c, loc, grid, ix, iy); discrete=true, parameters=(f_c,))
+                f_cv_i = interior(f_cv) |> Array
+                @test f_cv_i[:, 2:end-1] ≈ avy(f_c_i)
+            end
+            @testset "c2vc" begin
+                f_vc = Field(arch, grid, (Vertex(), Center()))
+                set!(f_vc, grid, (grid, loc, ix, iy, f_c) -> lerp(f_c, loc, grid, ix, iy); discrete=true, parameters=(f_c,))
+                f_vc_i = interior(f_vc) |> Array
+                @test f_vc_i[2:end-1, :] ≈ avx(f_c_i)
+            end
+        end
+        @testset "vertex" begin
+            f_v = Field(arch, grid, Vertex())
+            src = reshape(1:9, size(grid, Vertex())) |> collect
+            set!(f_v, src)
+            f_v_i = interior(f_v) |> Array
+            @testset "v2c" begin
                 f_c = Field(arch, grid, Center())
-                src = reshape(1:4, size(grid, Center())) |> collect
-                set!(f_c, src)
+                set!(f_c, grid, (grid, loc, ix, iy, f_v) -> lerp(f_v, loc, grid, ix, iy); discrete=true, parameters=(f_v,))
                 f_c_i = interior(f_c) |> Array
-                @testset "c2v" begin
-                    f_v = Field(arch, grid, Vertex())
-                    set!(f_v, grid, (grid, loc, ix, iy, f_c) -> lerp(f_c, loc, grid, ix, iy); discrete=true, parameters=(f_c,))
-                    f_v_i = interior(f_v) |> Array
-                    @test f_v_i[2:end-1, 2:end-1] ≈ av4(f_c_i)
-                end
-                @testset "c2c" begin
-                    f_c2 = Field(arch, grid, Center())
-                    set!(f_c2, grid, (grid, loc, ix, iy, f_c) -> lerp(f_c, loc, grid, ix, iy); discrete=true, parameters=(f_c,))
-                    f_c2_i = interior(f_c2) |> Array
-                    @test f_c_i ≈ f_c2_i
-                end
-                @testset "c2cv" begin
-                    f_cv = Field(arch, grid, (Center(), Vertex()))
-                    set!(f_cv, grid, (grid, loc, ix, iy, f_c) -> lerp(f_c, loc, grid, ix, iy); discrete=true, parameters=(f_c,))
-                    f_cv_i = interior(f_cv) |> Array
-                    @test f_cv_i[:, 2:end-1] ≈ avy(f_c_i)
-                end
-                @testset "c2vc" begin
-                    f_vc = Field(arch, grid, (Vertex(), Center()))
-                    set!(f_vc, grid, (grid, loc, ix, iy, f_c) -> lerp(f_c, loc, grid, ix, iy); discrete=true, parameters=(f_c,))
-                    f_vc_i = interior(f_vc) |> Array
-                    @test f_vc_i[2:end-1, :] ≈ avx(f_c_i)
-                end
+                @test f_c_i ≈ av4(f_v_i)
             end
-            @testset "vertex" begin
-                f_v = Field(arch, grid, Vertex())
-                src = reshape(1:9, size(grid, Vertex())) |> collect
-                set!(f_v, src)
-                f_v_i = interior(f_v) |> Array
-                @testset "v2c" begin
-                    f_c = Field(arch, grid, Center())
-                    set!(f_c, grid, (grid, loc, ix, iy, f_v) -> lerp(f_v, loc, grid, ix, iy); discrete=true, parameters=(f_v,))
-                    f_c_i = interior(f_c) |> Array
-                    @test f_c_i ≈ av4(f_v_i)
-                end
-                @testset "v2v" begin
-                    f_v2 = Field(arch, grid, Vertex())
-                    set!(f_v2, grid, (grid, loc, ix, iy, f_v) -> lerp(f_v, loc, grid, ix, iy); discrete=true, parameters=(f_v,))
-                    f_v2_i = interior(f_v2) |> Array
-                    @test f_v2_i ≈ f_v_i
-                end
-                @testset "v2cv" begin
-                    f_cv = Field(arch, grid, (Center(), Vertex()))
-                    set!(f_cv, grid, (grid, loc, ix, iy, f_v) -> lerp(f_v, loc, grid, ix, iy); discrete=true, parameters=(f_v,))
-                    f_cv_i = interior(f_cv) |> Array
-                    @test f_cv_i ≈ avx(f_v_i)
-                end
-                @testset "v2vc" begin
-                    f_vc = Field(arch, grid, (Vertex(), Center()))
-                    set!(f_vc, grid, (grid, loc, ix, iy, f_v) -> lerp(f_v, loc, grid, ix, iy); discrete=true, parameters=(f_v,))
-                    f_vc_i = interior(f_vc) |> Array
-                    @test f_vc_i ≈ avy(f_v_i)
-                end
+            @testset "v2v" begin
+                f_v2 = Field(arch, grid, Vertex())
+                set!(f_v2, grid, (grid, loc, ix, iy, f_v) -> lerp(f_v, loc, grid, ix, iy); discrete=true, parameters=(f_v,))
+                f_v2_i = interior(f_v2) |> Array
+                @test f_v2_i ≈ f_v_i
+            end
+            @testset "v2cv" begin
+                f_cv = Field(arch, grid, (Center(), Vertex()))
+                set!(f_cv, grid, (grid, loc, ix, iy, f_v) -> lerp(f_v, loc, grid, ix, iy); discrete=true, parameters=(f_v,))
+                f_cv_i = interior(f_cv) |> Array
+                @test f_cv_i ≈ avx(f_v_i)
+            end
+            @testset "v2vc" begin
+                f_vc = Field(arch, grid, (Vertex(), Center()))
+                set!(f_vc, grid, (grid, loc, ix, iy, f_v) -> lerp(f_v, loc, grid, ix, iy); discrete=true, parameters=(f_v,))
+                f_vc_i = interior(f_vc) |> Array
+                @test f_vc_i ≈ avy(f_v_i)
             end
         end
     end