diff --git a/src/QSM.jl b/src/QSM.jl
index 7695178..1dde493 100644
--- a/src/QSM.jl
+++ b/src/QSM.jl
@@ -80,8 +80,7 @@ const FFTW_NTHREADS = Ref{Int}(known(num_cores()))
     end
 
     function fftw_set_threading(lib::Symbol = :FFTW)
-        lib ∈ (:FFTW, :Polyester, :Threads) ||
-            throw(ArgumentError("lib must be one of :FFTW, :Polyester or :Threads, got :$(lib)"))
+        checkopts(lib, (:FFTW, :Polyester, :Threads), :lib)
 
         if lib ∈ (:Polyester, :Threads) && nthreads() < 2
             @warn "Cannot use $lib with FFTW. Defaulting to FFTW multi-threading" Threads.nthreads()
diff --git a/src/bgremove/ismv.jl b/src/bgremove/ismv.jl
index 48a8ec5..7bc4f8d 100644
--- a/src/bgremove/ismv.jl
+++ b/src/bgremove/ismv.jl
@@ -88,24 +88,24 @@ function _ismv!(
     S = _smv_kernel!(S, F̂, s, vsz, r, P)
 
     # constants
-    _δ = one(eltype(s)) - sqrt(eps(eltype(s)))
+    δ = one(eltype(s)) - sqrt(eps(eltype(s)))
 
     # erode mask
     s = tcopyto!(s, m) # in-place type conversion, reuse smv var
     m0 = tcopyto!(m0, m)
 
     F̂ = mul!(F̂, P, s)
-    @inbounds @batch for I in eachindex(F̂)
+    @batch for I in eachindex(F̂)
         F̂[I] *= S[I]
     end
 
     s = mul!(s, iP, F̂)
-    @inbounds @batch for I in eachindex(m)
-        m[I] = s[I] > _δ
+    @batch for I in eachindex(m)
+        m[I] = s[I] > δ
     end
 
     # iSMV
-    @inbounds for t in axes(f, 4)
+    for t in axes(f, 4)
         if verbose && size(f, 4) > 1
             @printf("Echo: %d/%d\n", t, size(f, 4))
         end
@@ -135,7 +135,7 @@ end
 
 
 function __ismv!(f::AbstractArray{T}, f0, S, bc, m, iP, F̂, P, tol, maxit, verbose) where {T}
-    @inbounds @batch for I in eachindex(f0)
+    @batch for I in eachindex(f0)
         f0[I] = m[I]*f[I]
     end
 
@@ -146,7 +146,7 @@ function __ismv!(f::AbstractArray{T}, f0, S, bc, m, iP, F̂, P, tol, maxit, verb
         @printf("iter%6s\tresnorm\n", "")
     end
 
-    @inbounds for i in 1:maxit
+    for i in 1:maxit
         if nr ≤ ϵ
             break
         end
diff --git a/src/bgremove/lbv.jl b/src/bgremove/lbv.jl
index bbde24e..554163f 100644
--- a/src/bgremove/lbv.jl
+++ b/src/bgremove/lbv.jl
@@ -118,7 +118,7 @@ function _lbv!(
     end
 
     # set boundaries
-    @inbounds for t in axes(f, 4)
+    for t in axes(f, 4)
         fct = @view(f[Rc,t])
         @batch for I in eachindex(fct, mc)
             fct[I] *= mc[I]
diff --git a/src/bgremove/pdf.jl b/src/bgremove/pdf.jl
index cf2a03d..f97e8db 100644
--- a/src/bgremove/pdf.jl
+++ b/src/bgremove/pdf.jl
@@ -85,8 +85,7 @@ function _pdf!(
         checkshape(W, f, (:W, :f))
     end
 
-    Dkernel ∈ (:k, :kspace, :i, :ispace) ||
-        throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)"))
+    checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel)
 
     # pad to fast fft size
     fp = padfastfft(@view(f[:,:,:,1]), pad, rfft=true)
@@ -126,7 +125,7 @@ function _pdf!(
     D = _dipole_kernel!(D, F̂, b, sz0, vsz, bdir, P, Dkernel, :rfft)
 
     # background mask
-    @inbounds @batch for I in eachindex(m)
+    @batch for I in eachindex(m)
         M̃[I] = !m[I]
     end
 
@@ -138,7 +137,7 @@ function _pdf!(
     elseif ndims(W) == 3
         # same weights for all echoes
         b = padarray!(b, W)
-        @inbounds @batch for I in eachindex(MW)
+        @batch for I in eachindex(MW)
             MW[I] = m[I] * b[I]
         end
 
@@ -147,7 +146,7 @@ function _pdf!(
         # VOID
     end
 
-    @inbounds for t in axes(f, 4)
+    for t in axes(f, 4)
         if verbose && size(f, 4) > 1
             @printf("Echo: %d/%d\n", t, size(f, 4))
         end
@@ -205,17 +204,17 @@ function _A_pdf!(Av, v, W, iP, D, F̂, P, M̃)
     v = reshape(v, size(W))
     x = reshape(Av, size(W))
 
-    @inbounds @batch for I in eachindex(x)
+    @batch for I in eachindex(x)
         x[I] = M̃[I] * v[I]
     end
 
     F̂ = mul!(F̂, P, x)
-    @inbounds @batch for I in eachindex(F̂)
+    @batch for I in eachindex(F̂)
         F̂[I] *= D[I]
     end
 
     x = mul!(x, iP, F̂)
-    @inbounds @batch for I in eachindex(x)
+    @batch for I in eachindex(x)
         x[I] *= W[I]
     end
 
@@ -227,17 +226,17 @@ function _At_pdf!(Atu, u, M̃, iP, D, F̂, P, W)
     u = reshape(u, size(W))
     x = reshape(Atu, size(W))
 
-    @inbounds @batch for I in eachindex(x)
+    @batch for I in eachindex(x)
         x[I] = W[I] * u[I]
     end
 
     F̂ = mul!(F̂, P, x)
-    @inbounds @batch for I in eachindex(F̂)
+    @batch for I in eachindex(F̂)
         F̂[I] *= D[I] # conj(D), D is real
     end
 
     x = mul!(x, iP, F̂)
-    @inbounds @batch for I in eachindex(x)
+    @batch for I in eachindex(x)
         x[I] *= M̃[I]
     end
 
diff --git a/src/bgremove/sharp.jl b/src/bgremove/sharp.jl
index 95e7ff9..be96890 100644
--- a/src/bgremove/sharp.jl
+++ b/src/bgremove/sharp.jl
@@ -80,27 +80,27 @@ function _sharp!(
     S = _smv_kernel!(S, F̂, s, vsz, r, P)
 
     # constants
-    _thr  = convert(eltype(S), thr)
-    _one  = one(eltype(S))
-    _zero = zero(eltype(F̂))
-    _δ    = one(eltype(s)) - sqrt(eps(eltype(s)))
+    thrT  = convert(eltype(S), thr)
+    oneT  = one(eltype(S))
+    zeroT = zero(eltype(F̂))
+    δ     = one(eltype(s)) - sqrt(eps(eltype(s)))
 
     # erode mask
     s = tcopyto!(s, m) # in-place type conversion, reuse smv var
 
     F̂ = mul!(F̂, P, s)
-    @inbounds @batch for I in eachindex(F̂)
+    @batch for I in eachindex(F̂)
         F̂[I] *= S[I]
-        S[I] = _one - S[I]
+        S[I] = oneT - S[I]
     end
 
     s = mul!(s, iP, F̂)
-    @inbounds @batch for I in eachindex(m)
-        m[I] = s[I] > _δ
+    @batch for I in eachindex(m)
+        m[I] = s[I] > δ
     end
 
     # SHARP
-    @inbounds for t in axes(f, 4)
+    for t in axes(f, 4)
         if t > 1
             fp = padarray!(fp, @view(f[Rc,t]))
         end
@@ -117,7 +117,7 @@ function _sharp!(
 
         F̂ = mul!(F̂, P, fp)
         @batch for I in eachindex(F̂)
-            F̂[I] = ifelse(abs(S[I]) < _thr, _zero, F̂[I]*inv(S[I]))
+            F̂[I] = ifelse(abs(S[I]) < thrT, zeroT, F̂[I]*inv(S[I]))
         end
 
         fp = mul!(fp, iP, F̂)
@@ -226,34 +226,34 @@ function _sharp!(
     fp = tfill!(fp, 0)
 
     # constants
-    _thr  = convert(eltype(S), thr)
-    _one  = one(eltype(S))
-    _zero = zero(eltype(F̂))
-    _δ    = one(eltype(s)) - sqrt(eps(eltype(s)))
+    thrT  = convert(eltype(S), thr)
+    oneT  = one(eltype(S))
+    zeroT = zero(eltype(F̂))
+    δ     = one(eltype(s)) - sqrt(eps(eltype(s)))
 
     # fft of original mask
     s = tcopyto!(s, mr) # in-place type conversion
     M̂ = mul!(M̂, P, s)
 
-    @inbounds for (i, r) in enumerate(rs)
+    for (i, r) in enumerate(rs)
         # get smv kernel
         S = _smv_kernel!(S, F̂, s, vsz, r, P)
 
         # erode mask
         @batch for I in eachindex(F̂)
             F̂[I] = S[I] * M̂[I]
-            S[I] = _one - S[I]
+            S[I] = oneT - S[I]
         end
 
         s = mul!(s, iP, F̂)
         @batch for I in eachindex(mr)
-            mr[I] = s[I] > _δ
+            mr[I] = s[I] > δ
         end
 
         # high-pass filter first (largest) kernel for deconvolution
         if i == 1
             @batch for I in eachindex(iS)
-                iS[I] = ifelse(abs(S[I]) < _thr, _zero, inv(S[I]))
+                iS[I] = ifelse(abs(S[I]) < thrT, zeroT, inv(S[I]))
             end
         end
 
@@ -301,7 +301,7 @@ function _sharp!(
     checkshape(fl, f, (:fl, :f))
     checkshape(smask, mask, (:smask, :mask))
     checkshape(axes(mask), axes(f)[1:3], (:mask, :f))
-    length(r) > 0 || throw(ArgumentError("r must not be empty"))
+    !isempty(r) || throw(ArgumentError("r must not be empty"))
 
     rs = sort!(unique(collect(r)), rev=true)
 
@@ -338,34 +338,34 @@ function _sharp!(
     fp = tfill!(fp, 0)
 
     # constants
-    _thr  = convert(eltype(S), thr)
-    _one  = one(eltype(S))
-    _zero = zero(eltype(F̂))
-    _δ    = one(eltype(s)) - sqrt(eps(eltype(s)))
+    thrT  = convert(eltype(S), thr)
+    oneT  = one(eltype(S))
+    zeroT = zero(eltype(F̂))
+    δ     = one(eltype(s)) - sqrt(eps(eltype(s)))
 
     # fft of original mask
     s = tcopyto!(s, mr) # in-place type conversion
     M̂ = mul!(M̂, P, s)
 
-    @inbounds for (i, r) in enumerate(rs)
+    for (i, r) in enumerate(rs)
         # get smv kernel
         S = _smv_kernel!(S, F̂, s, vsz, r, P)
 
         # erode mask
         @batch for I in eachindex(F̂)
             F̂[I] = S[I] * M̂[I]
-            S[I] = _one - S[I]
+            S[I] = oneT - S[I]
         end
 
         s = mul!(s, iP, F̂)
         @batch for I in eachindex(mr)
-            mr[I] = s[I] > _δ
+            mr[I] = s[I] > δ
         end
 
         # high-pass filter first (largest) kernel for deconvolution
         if i == 1
             @batch for I in eachindex(iS)
-                iS[I] = ifelse(abs(S[I]) < _thr, _zero, inv(S[I]))
+                iS[I] = ifelse(abs(S[I]) < thrT, zeroT, inv(S[I]))
             end
         end
 
@@ -391,7 +391,7 @@ function _sharp!(
 
     # deconvolution + high-pass filter
     F̂p = mul!(F̂p, P4, fp)
-    @inbounds for t in axes(F̂p, 4)
+    for t in axes(F̂p, 4)
         F̂t = @view(F̂p[:,:,:,t])
         @batch for I in eachindex(F̂t)
             F̂t[I] *= iS[I]
@@ -399,7 +399,7 @@ function _sharp!(
     end
 
     fp = mul!(fp, iP4, F̂p)
-    @inbounds for t in axes(fp, 4)
+    for t in axes(fp, 4)
         ft = @view(fp[:,:,:,t])
         @batch for I in eachindex(ft)
             ft[I] *= m[I]
diff --git a/src/inversion/direct.jl b/src/inversion/direct.jl
index d008bf1..aac5ccc 100644
--- a/src/inversion/direct.jl
+++ b/src/inversion/direct.jl
@@ -175,16 +175,9 @@ function _kdiv!(
     checkshape(x, f, (:x, :f))
     checkshape(axes(mask), axes(f)[1:3])
 
-    Dkernel ∈ (:k, :kspace, :i, :ispace) ||
-        throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)"))
-
-    method ∈ (:tkd, :tsvd, :tikh) ||
-        throw(ArgumentError("method must be one of :tkd, :tsvd, got :$(method)"))
-
-    if method == :tikh
-        reg ∈ (:identity, :gradient, :laplacian) ||
-            throw(ArgumentError("reg must be one of :identity, :gradient, :laplacian, got :$(reg)"))
-    end
+    checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel)
+    checkopts(method, (:tkd, :tsvd, :tikh), :method)
+    method == :tikh && checkopts(reg, (:identity, :gradient, :laplacian), :reg)
 
     # pad to fast fft size
     fp = padfastfft(@view(f[:,:,:,1]), pad, rfft=true)
@@ -208,7 +201,7 @@ function _kdiv!(
     # inverse k-space kernel
     iD = _kdiv_ikernel!(D, F̂, fp, vsz, P, lambda, reg, method)
 
-    @inbounds for t in axes(f, 4)
+    for t in axes(f, 4)
         fp = padarray!(fp, @view(f[:,:,:,t]))
 
         F̂ = mul!(F̂, P, fp)
@@ -243,22 +236,22 @@ function _kdiv_ikernel!(
         δ = convert(T, lambda)
         iδ = inv(δ)
 
-        @inbounds @batch for I in eachindex(D)
+        @batch for I in eachindex(D)
             D[I] = ifelse(abs(D[I]) ≤ δ, copysign(iδ, D[I]), inv(D[I]))
         end
 
     elseif method == :tsvd
         δ = convert(T, lambda)
-        _zero = zero(T)
+        zeroT = zero(T)
 
-        @inbounds @batch for I in eachindex(D)
-            D[I] = ifelse(abs(D[I]) ≤ δ, _zero, inv(D[I]))
+        @batch for I in eachindex(D)
+            D[I] = ifelse(abs(D[I]) ≤ δ, zeroT, inv(D[I]))
         end
 
 
     elseif method == :tikh
         λ = convert(T, lambda)
-        _zero = zero(T)
+        zeroT = zero(T)
 
         Γ = similar(D)
 
@@ -273,8 +266,8 @@ function _kdiv_ikernel!(
             Γ = tmap!(abs2, Γ)
         end
 
-        @inbounds @batch for I in eachindex(D)
-            D[I] = ifelse(iszero(D[I]) && iszero(Γ[I]), _zero, inv(D[I]*D[I] + λ*Γ[I])*D[I])
+        @batch for I in eachindex(D)
+            D[I] = ifelse(iszero(D[I]) && iszero(Γ[I]), zeroT, inv(D[I]*D[I] + λ*Γ[I])*D[I])
         end
 
     end
diff --git a/src/inversion/nltv.jl b/src/inversion/nltv.jl
index d192d47..0120fba 100644
--- a/src/inversion/nltv.jl
+++ b/src/inversion/nltv.jl
@@ -121,12 +121,11 @@ function _nltv!(
         end
     end
 
-    Dkernel ∈ (:k, :kspace, :i, :ispace) ||
-        throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)"))
+    checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel)
 
     # convert scalars
-    _zero = zero(T)
-    _one = one(T)
+    zeroT = zero(T)
+    oneT = one(T)
 
     ρ = convert(T, rho)
     μ = convert(T, mu)
@@ -201,7 +200,7 @@ function _nltv!(
     D = _dipole_kernel!(D, F̂, xp, sz0, vsz, bdir, P, Dkernel, :rfft)
     L = _laplace_kernel!(L, F̂, xp, vsz, P, negative=true)
 
-    @inbounds for t in axes(f, 4)
+    for t in axes(f, 4)
         if verbose && size(f, 4) > 1
             @printf("Echo: %d/%d\n", t, size(f, 4))
         end
@@ -211,7 +210,7 @@ function _nltv!(
         # iA = (μ*D^H*D - ρ*Δ)^-1
         @batch for I in eachindex(iA)
             a = μ*conj(D[I])*D[I] + ρ*L[I]
-            iA[I] = ifelse(iszero(a), _zero, inv(a))
+            iA[I] = ifelse(iszero(a), zeroT, inv(a))
         end
 
         F̂ = tfill!(F̂, 0)
@@ -234,9 +233,9 @@ function _nltv!(
             end
         end
 
-        ux = tfill!(ux, _zero)
-        uy = tfill!(uy, _zero)
-        uz = tfill!(uz, _zero)
+        ux = tfill!(ux, zeroT)
+        uy = tfill!(uy, zeroT)
+        uz = tfill!(uz, zeroT)
         dx, dy, dz = _gradfp!(dx, dy, dz, fp, vsz)
 
         if verbose
@@ -297,7 +296,7 @@ function _nltv!(
                 @batch for I in eachindex(ux)
                     λiρ = λiρWx[I]
                     u = ux[I] + dx[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     ux[I] = u - z
                     dx[I] = z - u + z
                 end
@@ -305,7 +304,7 @@ function _nltv!(
                 @batch for I in eachindex(uy)
                     λiρ = λiρWy[I]
                     u = uy[I] + dy[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     uy[I] = u - z
                     dy[I] = z - u + z
                 end
@@ -313,7 +312,7 @@ function _nltv!(
                 @batch for I in eachindex(uz)
                     λiρ = λiρWz[I]
                     u = uz[I] + dz[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     uz[I] = u - z
                     dz[I] = z - u + z
                 end
@@ -321,21 +320,21 @@ function _nltv!(
             else
                 @batch for I in eachindex(ux)
                     u = ux[I] + dx[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     ux[I] = u - z
                     dx[I] = z - u + z
                 end
 
                 @batch for I in eachindex(uy)
                     u = uy[I] + dy[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     uy[I] = u - z
                     dy[I] = z - u + z
                 end
 
                 @batch for I in eachindex(uz)
                     u = uz[I] + dz[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     uz[I] = u - z
                     dz[I] = z - u + z
                 end
@@ -368,7 +367,7 @@ function _nltv!(
                         φ = fp[I]
 
                         s, c = sincos_fast(z - φ)
-                        δ = (muladd(w, s, z) - u) / muladd(w, c, _one)
+                        δ = (muladd(w, s, z) - u) / muladd(w, c, oneT)
 
                         ze[I] = z - δ
                         threadlocal[1] = muladd(δ, δ, threadlocal[1])
diff --git a/src/inversion/rts.jl b/src/inversion/rts.jl
index ffc39a7..7bb8eb8 100644
--- a/src/inversion/rts.jl
+++ b/src/inversion/rts.jl
@@ -91,12 +91,10 @@ function _rts!(
 
     checkshape(x, f, (:x, :f))
     checkshape(axes(mask), axes(f)[1:3])
-
-    Dkernel ∈ (:k, :kspace, :i, :ispace) ||
-        throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)"))
+    checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel)
 
     # convert scalars
-    _zero = zero(T)
+    zeroT = zero(T)
 
     ρ = convert(T, rho)
     δ = convert(T, delta)
@@ -159,11 +157,11 @@ function _rts!(
     L = _laplace_kernel!(L, F̂, xp, vsz, P, negative=true)
 
     # mask of well-conditioned frequencies
-    @inbounds @batch for I in eachindex(M)
-        M[I] = ifelse(abs(D[I]) > δ, μ, _zero)
+    @batch for I in eachindex(M)
+        M[I] = ifelse(abs(D[I]) > δ, μ, zeroT)
     end
 
-    @inbounds for t in axes(f, 4)
+    for t in axes(f, 4)
         if verbose && size(f, 4) > 1
             @printf("Echo: %d/%d\n", t, size(f, 4))
         end
@@ -174,7 +172,7 @@ function _rts!(
         # Step 1: Well-conditioned
         ######################################################################
         B̂ = mul!(B̂, P, xp)
-        lsmr!(WS; lambda=_zero, atol=_zero, btol=_zero, maxit=lstol)
+        lsmr!(WS; lambda=zeroT, atol=zeroT, btol=zeroT, maxit=lstol)
 
         xp = mul!(xp, iP, X̂)
 
@@ -190,8 +188,8 @@ function _rts!(
         @batch for I in eachindex(F̂)
             a = muladd(ρ, L[I], M[I])
             if iszero(a)
-                F̂[I] = _zero
-                iA[I] = _zero
+                F̂[I] = zeroT
+                iA[I] = zeroT
             else
                 ia = inv(a)
                 F̂[I] = ia * M[I] * X̂[I]
@@ -200,9 +198,9 @@ function _rts!(
         end
 
         nr = typemax(T)
-        px = tfill!(px, _zero)
-        py = tfill!(py, _zero)
-        pz = tfill!(pz, _zero)
+        px = tfill!(px, zeroT)
+        py = tfill!(py, zeroT)
+        pz = tfill!(pz, zeroT)
         vx, vy, vz = _gradfp!(vx, vy, vz, xp, vsz)
 
         if verbose
@@ -258,7 +256,7 @@ function _rts!(
             @batch for I in eachindex(px)
                 d = dx[I]
                 p = px[I] + d
-                y = ifelse(abs(p) ≤ iρ, _zero, copysign(abs(p)-iρ, p))
+                y = ifelse(abs(p) ≤ iρ, zeroT, copysign(abs(p)-iρ, p))
                 px[I] = p - y
                 vx[I] = y - p + y
                 dx[I] = d - y
@@ -267,7 +265,7 @@ function _rts!(
             @batch for I in eachindex(py)
                 d = dy[I]
                 p = py[I] + d
-                y = ifelse(abs(p) ≤ iρ, _zero, copysign(abs(p)-iρ, p))
+                y = ifelse(abs(p) ≤ iρ, zeroT, copysign(abs(p)-iρ, p))
                 py[I] = p - y
                 vy[I] = y - p + y
                 dy[I] = d - y
@@ -276,7 +274,7 @@ function _rts!(
             @batch for I in eachindex(pz)
                 d = dz[I]
                 p = pz[I] + d
-                y = ifelse(abs(p) ≤ iρ, _zero, copysign(abs(p)-iρ, p))
+                y = ifelse(abs(p) ≤ iρ, zeroT, copysign(abs(p)-iρ, p))
                 pz[I] = p - y
                 vz[I] = y - p + y
                 dz[I] = d - y
@@ -336,7 +334,7 @@ end
 
 
 function _A_rts!(Dv, v, D)
-    @inbounds @batch for I in eachindex(Dv)
+    @batch for I in eachindex(Dv)
         Dv[I] = D[I]*v[I]
     end
     return Dv
diff --git a/src/inversion/tv.jl b/src/inversion/tv.jl
index 057e91a..95003ad 100644
--- a/src/inversion/tv.jl
+++ b/src/inversion/tv.jl
@@ -114,11 +114,10 @@ function _tv!(
         end
     end
 
-    Dkernel ∈ (:k, :kspace, :i, :ispace) ||
-        throw(ArgumentError("Dkernel must be one of :k, :kspace, :i, :ispace, got :$(Dkernel)"))
+    checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel)
 
     # convert scalars
-    _zero = zero(T)
+    zeroT = zero(T)
 
     ρ = convert(T, rho)
     μ = convert(T, mu)
@@ -189,7 +188,7 @@ function _tv!(
     D = _dipole_kernel!(D, F̂, xp, sz0, vsz, bdir, P, Dkernel, :rfft)
     L = _laplace_kernel!(L, F̂, xp, vsz, P, negative=true)
 
-    @inbounds for t in axes(f, 4)
+    for t in axes(f, 4)
         if verbose && size(f, 4) > 1
             @printf("Echo: %d/%d\n", t, size(f, 4))
         end
@@ -205,7 +204,7 @@ function _tv!(
 
             @batch for I in eachindex(iA)
                 a = μ*conj(D[I])*D[I] + ρ*L[I]
-                iA[I] = ifelse(iszero(a), _zero, inv(a))
+                iA[I] = ifelse(iszero(a), zeroT, inv(a))
             end
 
             fw = padarray!(fw, @view(W[:,:,:,min(t, end)]))
@@ -225,8 +224,8 @@ function _tv!(
             @batch for I in eachindex(F̂)
                 a = conj(D[I])*D[I] + ρ*L[I]
                 if iszero(a)
-                    F̂[I] = _zero
-                    iA[I] = _zero
+                    F̂[I] = zeroT
+                    iA[I] = zeroT
                 else
                     ia = inv(a)
                     F̂[I] = ia * conj(D[I]) * X̂[I]
@@ -246,9 +245,9 @@ function _tv!(
             end
         end
 
-        ux = tfill!(ux, _zero)
-        uy = tfill!(uy, _zero)
-        uz = tfill!(uz, _zero)
+        ux = tfill!(ux, zeroT)
+        uy = tfill!(uy, zeroT)
+        uz = tfill!(uz, zeroT)
         dx, dy, dz = _gradfp!(dx, dy, dz, xp, vsz)
 
         if verbose
@@ -309,7 +308,7 @@ function _tv!(
                 @batch for I in eachindex(ux)
                     λiρ = λiρWx[I]
                     u = ux[I] + dx[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     ux[I] = u - z
                     dx[I] = z - u + z
                 end
@@ -317,7 +316,7 @@ function _tv!(
                 @batch for I in eachindex(uy)
                     λiρ = λiρWy[I]
                     u = uy[I] + dy[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     uy[I] = u - z
                     dy[I] = z - u + z
                 end
@@ -325,7 +324,7 @@ function _tv!(
                 @batch for I in eachindex(uz)
                     λiρ = λiρWz[I]
                     u = uz[I] + dz[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     uz[I] = u - z
                     dz[I] = z - u + z
                 end
@@ -333,21 +332,21 @@ function _tv!(
             else
                 @batch for I in eachindex(ux)
                     u = ux[I] + dx[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     ux[I] = u - z
                     dx[I] = z - u + z
                 end
 
                 @batch for I in eachindex(uy)
                     u = uy[I] + dy[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     uy[I] = u - z
                     dy[I] = z - u + z
                 end
 
                 @batch for I in eachindex(uz)
                     u = uz[I] + dz[I]
-                    z = ifelse(abs(u) ≤ λiρ, _zero, copysign(abs(u)-λiρ, u))
+                    z = ifelse(abs(u) ≤ λiρ, zeroT, copysign(abs(u)-λiρ, u))
                     uz[I] = u - z
                     dz[I] = z - u + z
                 end
diff --git a/src/unwrap/laplacian.jl b/src/unwrap/laplacian.jl
index 123056d..865bf15 100644
--- a/src/unwrap/laplacian.jl
+++ b/src/unwrap/laplacian.jl
@@ -48,10 +48,9 @@ function unwrap_laplacian(
     solver::Symbol = :mgpcg
 ) where {T<:AbstractFloat, N}
     N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
-    checkshape(axes(mask), axes(phas)[1:3], (:mask, :phas))
 
-    solver ∈ (:dct, :fft, :mgpcg) ||
-        throw(ArgumentError("solver must be one of :dct, :fft, :mgpcg"))
+    checkshape(axes(mask), axes(phas)[1:3], (:mask, :phas))
+    checkopts(solver, (:dct, :fft, :mgpcg), :solver)
 
     d2uphas = wrapped_laplacian(phas, vsz)
 
@@ -78,7 +77,7 @@ function unwrap_laplacian(
         )
 
         # set boundaries
-        @inbounds for t in axes(d2uphas, 4)
+        for t in axes(d2uphas, 4)
             d2ut = @view(d2uphas[:,:,:,t])
             @batch for I in eachindex(d2ut, mask)
                 d2ut[I] *= mask[I]
@@ -113,7 +112,7 @@ function wrapped_laplacian!(
     tsz = padded_tilesize(T, (2, 2, 2), 1)
     R = vec(collect(TileIterator((2:nx-1, 2:ny-1, 2:nz-1), tsz)))
 
-    @inbounds @batch per=thread for (I, J, K) in R
+    @batch per=thread for (I, J, K) in R
         for k in K
             for j in J
                 for i in I
@@ -145,20 +144,20 @@ function wrapped_laplacian!(
     tsz = padded_tilesize(T, (2, 2, 2), 1)
     R = vec(collect(TileIterator((2:nx-1, 2:ny-1, 2:nz-1), tsz)))
 
-    @inbounds for t in axes(u, 4)
-        _u = @view(u[:,:,:,t])
-        _d2u = @view(d2u[:,:,:,t])
+    for t in axes(u, 4)
+        ut = @view(u[:,:,:,t])
+        d2ut = @view(d2u[:,:,:,t])
         @batch per=thread for (I, J, K) in R
             for k in K
                 for j in J
                     for i in I
-                        Δ =        -dz2 *_wrap(_u[i,j,k] - _u[i,j,k-1])
-                        Δ = muladd(-dy2, _wrap(_u[i,j,k] - _u[i,j-1,k]), Δ)
-                        Δ = muladd(-dx2, _wrap(_u[i,j,k] - _u[i-1,j,k]), Δ)
-                        Δ = muladd( dx2, _wrap(_u[i+1,j,k] - _u[i,j,k]), Δ)
-                        Δ = muladd( dy2, _wrap(_u[i,j+1,k] - _u[i,j,k]), Δ)
-                        Δ = muladd( dz2, _wrap(_u[i,j,k+1] - _u[i,j,k]), Δ)
-                        _d2u[i,j,k] = Δ
+                        Δ =        -dz2 *_wrap(ut[i,j,k] - ut[i,j,k-1])
+                        Δ = muladd(-dy2, _wrap(ut[i,j,k] - ut[i,j-1,k]), Δ)
+                        Δ = muladd(-dx2, _wrap(ut[i,j,k] - ut[i-1,j,k]), Δ)
+                        Δ = muladd( dx2, _wrap(ut[i+1,j,k] - ut[i,j,k]), Δ)
+                        Δ = muladd( dy2, _wrap(ut[i,j+1,k] - ut[i,j,k]), Δ)
+                        Δ = muladd( dz2, _wrap(ut[i,j,k+1] - ut[i,j,k]), Δ)
+                        d2ut[i,j,k] = Δ
                     end
                 end
             end
@@ -179,6 +178,7 @@ function wrapped_laplacian_boundary_neumann!(
     dx::NTuple{3, Real}
 ) where {T<:AbstractFloat, N}
     N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
+
     checkshape(d2u, u, (:d2u, :u))
 
     zeroT = zero(T)
@@ -187,37 +187,37 @@ function wrapped_laplacian_boundary_neumann!(
 
     R = edge_indices(axes(u)[1:3])
 
-    @inbounds for t in axes(u, 4)
-        _u = @view(u[:,:,:,t])
-        _d2u = @view(d2u[:,:,:,t])
+    for t in axes(u, 4)
+        ut = @view(u[:,:,:,t])
+        d2ut = @view(d2u[:,:,:,t])
         @batch per=thread for (i, j, k) in R
             Δ = zeroT
 
             if k > 1
-                Δ = muladd(-dz2, _wrap(_u[i,j,k] - _u[i,j,k-1]), Δ)
+                Δ = muladd(-dz2, _wrap(ut[i,j,k] - ut[i,j,k-1]), Δ)
             end
 
             if j > 1
-                Δ = muladd(-dy2, _wrap(_u[i,j,k] - _u[i,j-1,k]), Δ)
+                Δ = muladd(-dy2, _wrap(ut[i,j,k] - ut[i,j-1,k]), Δ)
             end
 
             if i > 1
-                Δ = muladd(-dx2, _wrap(_u[i,j,k] - _u[i-1,j,k]), Δ)
+                Δ = muladd(-dx2, _wrap(ut[i,j,k] - ut[i-1,j,k]), Δ)
             end
 
             if i < nx
-                Δ = muladd( dx2, _wrap(_u[i+1,j,k] - _u[i,j,k]), Δ)
+                Δ = muladd( dx2, _wrap(ut[i+1,j,k] - ut[i,j,k]), Δ)
             end
 
             if j < ny
-                Δ = muladd( dy2, _wrap(_u[i,j+1,k] - _u[i,j,k]), Δ)
+                Δ = muladd( dy2, _wrap(ut[i,j+1,k] - ut[i,j,k]), Δ)
             end
 
             if k < nz
-                Δ = muladd( dz2, _wrap(_u[i,j,k+1] - _u[i,j,k]), Δ)
+                Δ = muladd( dz2, _wrap(ut[i,j,k+1] - ut[i,j,k]), Δ)
             end
 
-            _d2u[i,j,k] = Δ
+            d2ut[i,j,k] = Δ
         end
     end
 
@@ -231,6 +231,7 @@ function wrapped_laplacian_boundary_periodic!(
     dx::NTuple{3, Real}
 ) where {T<:AbstractFloat, N}
     N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
+
     checkshape(d2u, u, (:d2u, :u))
 
     nx, ny, nz = size(u)[1:3]
@@ -238,29 +239,29 @@ function wrapped_laplacian_boundary_periodic!(
 
     R = edge_indices(axes(u)[1:3])
 
-    @inbounds for t in axes(u, 4)
-        _u = @view(u[:,:,:,t])
-        _d2u = @view(d2u[:,:,:,t])
+    for t in axes(u, 4)
+        ut = @view(u[:,:,:,t])
+        d2ut = @view(d2u[:,:,:,t])
         @batch per=thread for (i, j, k) in R
-            v = k == 1 ? _u[i,j,end] : _u[i,j,k-1]
-            Δ = -dz2 * _wrap(_u[i,j,k] - v)
+            v = k == 1 ? ut[i,j,end] : ut[i,j,k-1]
+            Δ = -dz2 * _wrap(ut[i,j,k] - v)
 
-            v = j == 1 ? _u[i,end,k] : _u[i,j-1,k]
-            Δ = muladd(-dy2, _wrap(_u[i,j,k] - v), Δ)
+            v = j == 1 ? ut[i,end,k] : ut[i,j-1,k]
+            Δ = muladd(-dy2, _wrap(ut[i,j,k] - v), Δ)
 
-            v = i == 1 ? _u[end,j,k] : _u[i-1,j,k]
-            Δ = muladd(-dx2, _wrap(_u[i,j,k] - v), Δ)
+            v = i == 1 ? ut[end,j,k] : ut[i-1,j,k]
+            Δ = muladd(-dx2, _wrap(ut[i,j,k] - v), Δ)
 
-            v = i == nx ? _u[1,j,k] : _u[i+1,j,k]
-            Δ = muladd(dx2, _wrap(v - _u[i,j,k]), Δ)
+            v = i == nx ? ut[1,j,k] : ut[i+1,j,k]
+            Δ = muladd(dx2, _wrap(v - ut[i,j,k]), Δ)
 
-            v = j == ny ? _u[i,1,k] : _u[i,j+1,k]
-            Δ = muladd(dy2, _wrap(v - _u[i,j,k]), Δ)
+            v = j == ny ? ut[i,1,k] : ut[i,j+1,k]
+            Δ = muladd(dy2, _wrap(v - ut[i,j,k]), Δ)
 
-            v = k == nz ? _u[i,j,1] : _u[i,j,k+1]
-            Δ = muladd(dz2, _wrap(v - _u[i,j,k]), Δ)
+            v = k == nz ? ut[i,j,1] : ut[i,j,k+1]
+            Δ = muladd(dz2, _wrap(v - ut[i,j,k]), Δ)
 
-            _d2u[i,j,k] = Δ
+            d2ut[i,j,k] = Δ
         end
     end
 
diff --git a/src/utils/fd.jl b/src/utils/fd.jl
index 4ef2649..81aceec 100644
--- a/src/utils/fd.jl
+++ b/src/utils/fd.jl
@@ -64,7 +64,7 @@ function _gradfp!(
     hy = convert(Tdy, inv(h[2]))
     hz = convert(Tdz, inv(h[3]))
 
-    @inbounds @batch for k in 1:nz-1
+    @batch for k in 1:nz-1
         for j in 1:ny-1
             for i in 1:nx-1
                 dx[i,j,k] = hx*(u[i+1,j,k] - u[i,j,k])
@@ -88,7 +88,7 @@ function _gradfp!(
         dz[nx,ny,k] = hz*(u[nx,ny,k+1] - u[nx,ny,k])
     end
 
-    @inbounds @batch for j in 1:ny-1
+    @batch for j in 1:ny-1
         for i in 1:nx-1
             dz[i,j,nz] = hz*(u[i,j,1] - u[i,j,nz])
             dx[i,j,nz] = hx*(u[i+1,j,nz] - u[i,j,nz])
@@ -100,17 +100,15 @@ function _gradfp!(
         dy[nx,j,nz] = hy*(u[nx,j+1,nz] - u[nx,j,nz])
     end
 
-    @inbounds @batch for i in 1:nx-1
+    @batch for i in 1:nx-1
         dz[i,ny,nz] = hz*(u[i,ny,1] - u[i,ny,nz])
         dy[i,ny,nz] = hy*(u[i,1,nz] - u[i,ny,nz])
         dx[i,ny,nz] = hx*(u[i+1,ny,nz] - u[i,ny,nz])
     end
 
-    @inbounds begin
-        dz[nx,ny,nz] = hz*(u[nx,ny,1] - u[nx,ny,nz])
-        dy[nx,ny,nz] = hy*(u[nx,1,nz] - u[nx,ny,nz])
-        dx[nx,ny,nz] = hx*(u[1,ny,nz] - u[nx,ny,nz])
-    end
+    dz[nx,ny,nz] = hz*(u[nx,ny,1] - u[nx,ny,nz])
+    dy[nx,ny,nz] = hy*(u[nx,1,nz] - u[nx,ny,nz])
+    dx[nx,ny,nz] = hx*(u[1,ny,nz] - u[nx,ny,nz])
 
     return dx, dy, dz
 end
@@ -185,14 +183,12 @@ function _gradfp_adj!(
     hy = convert(Tdy, -inv(h[2]))
     hz = convert(Tdz, -inv(h[3]))
 
-    @inbounds begin
-        d2u[1,1,1] =
-            hx*(dx[1,1,1] - dx[nx,1,1]) +
-            hy*(dy[1,1,1] - dy[1,ny,1]) +
-            hz*(dz[1,1,1] - dz[1,1,nz])
-    end
+    d2u[1,1,1] =
+        hx*(dx[1,1,1] - dx[nx,1,1]) +
+        hy*(dy[1,1,1] - dy[1,ny,1]) +
+        hz*(dz[1,1,1] - dz[1,1,nz])
 
-    @inbounds @batch for i in 2:nx
+    @batch for i in 2:nx
         ux = dx[i,1,1] - dx[i-1,1,1]
         uy = dy[i,1,1] - dy[i,ny,1]
         uz = dz[i,1,1] - dz[i,1,nz]
@@ -203,7 +199,7 @@ function _gradfp_adj!(
         d2u[i,1,1] = u
     end
 
-    @inbounds @batch for j in 2:ny
+    @batch for j in 2:ny
         uy = dy[1,j,1] - dy[1,j-1,1]
         ux = dx[1,j,1] - dx[nx,j,1]
         uz = dz[1,j,1] - dz[1,j,nz]
@@ -225,7 +221,7 @@ function _gradfp_adj!(
         end
     end
 
-    @inbounds @batch for k in 2:nz
+    @batch for k in 2:nz
         uz = dz[1,1,k] - dz[1,1,k-1]
         ux = dx[1,1,k] - dx[nx,1,k]
         uy = dy[1,1,k] - dy[1,ny,k]
@@ -322,17 +318,14 @@ function _lap!(
     u::AbstractArray{T, 3},
     h::NTuple{3, Real},
 ) where {T<:AbstractFloat}
-    idx2 = convert(T, inv(h[1]*h[1]))
-    idy2 = convert(T, inv(h[2]*h[2]))
-    idz2 = convert(T, inv(h[3]*h[3]))
+    idx2, idy2, idz2 = convert.(T, inv.(h.*h))
     D = -2*(idx2 + idy2 + idz2)
 
-    ax = map(a -> first(a)+1:last(a)-1, axes(u))
-
     tsz = padded_tilesize(T, (2, 2, 2), 1)
+    ax = map(a -> first(a)+1:last(a)-1, axes(u))
     R = vec(collect(TileIterator(ax, tsz)))
 
-    @inbounds @batch for (I, J, K) in R
+    @batch for (I, J, K) in R
         for k in K
             for j in J
                 for i in I
@@ -367,17 +360,14 @@ function _lap!(
     u::AbstractArray{T, 4},
     h::NTuple{3, Real},
 ) where {T<:AbstractFloat}
-    idx2 = convert(T, inv(h[1]*h[1]))
-    idy2 = convert(T, inv(h[2]*h[2]))
-    idz2 = convert(T, inv(h[3]*h[3]))
+    idx2, idy2, idz2 = convert.(T, inv.(h.*h))
     D = -2*(idx2 + idy2 + idz2)
 
-    ax = map(a -> first(a)+1:last(a)-1, axes(u)[1:3])
-
     tsz = padded_tilesize(T, (2, 2, 2), 1)
+    ax = map(a -> first(a)+1:last(a)-1, axes(u)[1:3])
     R = vec(collect(TileIterator(ax, tsz)))
 
-    @inbounds for t in axes(u, 4)
+    for t in axes(u, 4)
         ut = @view(u[:,:,:,t])
         d2ut = @view(d2u[:,:,:,t])
         @batch for (I, J, K) in R
diff --git a/src/utils/generate.jl b/src/utils/generate.jl
deleted file mode 100644
index aa933fe..0000000
--- a/src/utils/generate.jl
+++ /dev/null
@@ -1,84 +0,0 @@
-using FFTW
-using MAT
-using Random
-using Printf
-using StaticArrays
-
-include("models.jl")
-
-
-"""
-Generate magnetic susceptibility maps and the corresponding fields due to
-randomly placed spheres and randomly oriented infinite cylinders subject to a
-static B-field inside a 3d rectangular volume.
-
-Method
-------
-    1. Pick array size, voxel size for grid (N/n grid)
-    2. for n iterations
-        i.   pick number of tries, ntry, to check for valid config
-        ii.  randomly place sources inside grid until no valid config can be
-             found ntry times in a row
-        iii. valid config means no overlapping structures. concentric allowed.
-"""
-function generate(
-    sz::NTuple{N, Integer},
-    vsz::NTuple{3, Real};
-    outpath::Union{Nothing, AbstractString} = nothing,
-    nsources::Integer = typemax(Int)-1,
-    density::Real = 0.51,
-    bdir::NTuple{3, Real} = (0, 0, 1),
-    χμ::Float64 = 1e-6,
-    χσ::Float64 = 0.07,
-    seed::Int = 12345,
-) where {N}
-    @assert nsources > 0
-    @assert density > 0 && density < 1
-    Random.seed!(seed)
-
-    if !(norm(bdir) ≈ 1)
-        @warn "B direction vector is not normalized. Normalizing..."
-        bdir = bdir ./ norm(bdir)
-    end
-
-    x = zeros(sz)
-    f = zeros(sz)
-    m = Array{Bool}(falses(sz))
-
-    rmax = cld(minimum(sz), 5)
-    rmax = 7
-    Cs = ntuple(i -> rmax+1:sz[i]-(rmax+1)-1, Val(3))
-    Cs = ntuple(i -> 10*rmax+1:sz[i]-(10*rmax+1)-1, Val(3))
-
-    ns = 0
-    while true
-        model = rand((:sphere, :cylinder))
-        model = :sphere
-
-        r = rand(3:rmax)
-        c = rand.(Cs)
-        χ = χμ + χσ*randn()
-
-        if model == :sphere
-            x, f, m = sphere!(x, f, m, c, r, χ, bdir, vsz)
-            ns += 1
-
-        elseif model == :cylinder
-            θ = π*rand()
-            ϕ = 2π*rand()
-
-            x, f, m = cylinder!(x, f, m, c, r, θ, ϕ, χ, bdir, vsz)
-            ns += 1
-        end
-
-        if sum(m)/length(m) > density || ns >= nsources
-            break
-        end
-    end
-
-    if outpath !== nothing
-        matwrite(outpath, Dict(:x => x, :f => f, :m => m), compress=true)
-    end
-
-    return x, f, m
-end
diff --git a/src/utils/kernels.jl b/src/utils/kernels.jl
index 6f189f9..cdbe0ea 100644
--- a/src/utils/kernels.jl
+++ b/src/utils/kernels.jl
@@ -59,16 +59,13 @@ function dipole_kernel(
     transform::Union{Nothing, Symbol} = nothing,
     shift::Bool = false
 ) where {T<:AbstractFloat}
-    !any(<=(0), sz)     || throw(ArgumentError("invalid array size"))
-    !any(<=(0), vsz)    || throw(DomainError(vsz, "vsz must be > 0"))
-    !iszero(norm(bdir)) || throw(DomainError(bdir, "norm of bdir must not be zero"))
-    !any(<(0), dsz)     || throw(DomainError(dsz, "dsz must be ≥ 0"))
+    all(>(0), sz)  || throw(ArgumentError("invalid array size"))
+    all(>(0), vsz) || throw(DomainError(vsz, "vsz must be > 0"))
+    all(>(0), dsz) || throw(DomainError(dsz, "dsz must be > 0"))
+    norm(bdir) > 0 || throw(DomainError(bdir, "bdir must not be zero"))
 
-    method ∈ (:k, :kspace, :i, :ispace) ||
-        throw(ArgumentError("method must be one of :k, :kspace, :i, :ispace, got :$(method)"))
-
-    transform === nothing || transform ∈ (:fft, :rfft) ||
-        throw(ArgumentError("transform must be one of nothing, :rfft, :fft, got :$(transform)"))
+    checkopts(Dkernel, (:k, :kspace, :i, :ispace), :Dkernel)
+    transform !== nothing && checkopts(transform, (:fft, :rfft), :transform)
 
     if method == :k || method == :kspace
         D = Array{T, 3}(undef, transform == :rfft ? (sz[1]>>1 + 1, sz[2], sz[3]) : sz)
@@ -178,7 +175,7 @@ function _dipole_kernel!(
         X, Y, Z = _grid(sz, dx, fft=true, shift=shift)
         a = 1 / 3
 
-        @inbounds @batch for k in 1:nz
+        @batch for k in 1:nz
             for j in 1:ny
                 for i in 1:nx
                     k̂ = SVector{3, Float64}(X[i], Y[j], Z[k])
@@ -201,10 +198,10 @@ function _dipole_kernel!(
         X, Y, Z = _grid(size(D), dx, shift=shift)
         a = prod(dx) * inv(4*pi)
 
-        @inbounds @batch for k in 1:nz
+        @batch for k in 1:nz
             z = Z[k]
             zout = z < Z1 || z > Z2
-            @fastpow for j in 1:ny
+            for j in 1:ny
                 y = Y[j]
                 zyout = zout || y < Y1 || y > Y2
                 for i in 1:nx
@@ -217,7 +214,7 @@ function _dipole_kernel!(
                         r = SVector{3, Float64}(x, y, z)
                         rz = r⋅B
                         r2 = r⋅r
-                        D[i,j,k] = a * (3*rz^2-r2)/sqrt(r2^5)
+                        @fastpow D[i,j,k] = a * (3*rz^2-r2)/sqrt(r2^5)
                     end
                 end
             end
@@ -277,12 +274,9 @@ function smv_kernel(
     transform::Union{Nothing, Symbol} = nothing,
     shift::Bool = false
 ) where {T<:AbstractFloat}
-    !any(<=(0), sz)  || throw(ArgumentError("invalid array size"))
-    !any(<=(0), vsz) || throw(DomainError(vsz, "vsz must be > 0"))
-    r >= 0           || throw(DomainError(r, "r must be ≥ 0"))
-
-    transform === nothing || transform ∈ (:fft, :rfft) ||
-        throw(ArgumentError("transform must be one of nothing, :rfft, :fft, got :$(transform)"))
+    all(>(0), sz)  || throw(ArgumentError("invalid array size"))
+    all(>(0), vsz) || throw(DomainError(vsz, "vsz must be > 0"))
+    transform !== nothing && checkopts(transform, (:fft, :rfft), :transform)
 
     if transform === nothing
         s = Array{T, 3}(undef, sz)
@@ -379,14 +373,14 @@ function _sphere!(
 
     X, Y, Z = _grid(sz, dx, shift=shift)
 
-    _zero = zero(T)
-    _one = one(T)
+    zeroT = zero(T)
+    oneT = one(T)
 
     @inbounds for k in 1:nz
         for j in 1:ny
             for i in 1:nx
                 R = SVector{3, Float64}(X[i], Y[j], Z[k])
-                s[i,j,k] = ifelse(R⋅R <= r2, _one, _zero)
+                s[i,j,k] = ifelse(R⋅R <= r2, oneT, zeroT)
             end
         end
     end
@@ -454,11 +448,9 @@ function laplace_kernel(
     transform::Union{Nothing, Symbol} = nothing,
     shift::Bool = false
 ) where {T<:AbstractFloat}
-    !any(<=(0), sz)  || throw(ArgumentError("invalid array size"))
-    !any(<=(0), vsz) || throw(DomainError(vsz, "vsz must be > 0"))
-
-    transform === nothing || transform ∈ (:fft, :rfft) ||
-        throw(ArgumentError("transform must be one of nothing, :rfft, :fft, got :$(transform)"))
+    all(>(0), sz)  || throw(ArgumentError("invalid array size"))
+    all(>(0), vsz) || throw(DomainError(vsz, "vsz must be > 0"))
+    transform !== nothing && checkopts(transform, (:fft, :rfft), :transform)
 
     if transform === nothing
         Δ = Array{T, 3}(undef, sz)
diff --git a/src/utils/lsmr.jl b/src/utils/lsmr.jl
index ed60498..d6384b1 100644
--- a/src/utils/lsmr.jl
+++ b/src/utils/lsmr.jl
@@ -365,7 +365,7 @@ end
 @inline function _xpby_norm!(X::Array, b::Real, Y::Array{T}) where {T}
     Tr = real(T)
     bT = convert(Tr, b)
-    @inbounds @batch threadlocal=zero(T)::T for I in eachindex(Y)
+    @batch threadlocal=zero(T)::T for I in eachindex(Y)
         Y[I] = muladd(bT, Y[I], X[I])
         threadlocal = muladd(conj(Y[I]), Y[I], threadlocal)
     end
@@ -381,7 +381,7 @@ end
 
 @inline function _scal!(a, X::Array{T}) where {T}
     aT = convert(real(T), a)
-    @inbounds @batch for I in eachindex(X)
+    @batch for I in eachindex(X)
         X[I] *= aT
     end
     return X
@@ -417,7 +417,7 @@ end
     σ = ζ / (ρ * ρbar)
     ν = -θnew / ρ
 
-    @inbounds @batch threadlocal=zero(T)::T for I in eachindex(x)
+    @batch threadlocal=zero(T)::T for I in eachindex(x)
         h̄[I] = muladd(δ, h̄[I], h[I])
         x[I] = muladd(σ, h̄[I], x[I])
         h[I] = muladd(ν, h[I], v[I])
diff --git a/src/utils/models.jl b/src/utils/models.jl
deleted file mode 100644
index bf5bdbd..0000000
--- a/src/utils/models.jl
+++ /dev/null
@@ -1,199 +0,0 @@
-using FastPow
-using LinearAlgebra
-using Polyester
-using StaticArrays
-
-
-"""
-    sphere!(x, f, m, c, r, χ, bdir, vsz)
-
-Add a sphere to grid.
-
-### Arguments
-- `x::AbstractArray{T1, 3}`: magnetic susceptibilities
-- `f::AbstractArray{T2, 3}`: local field generated by magnetic susceptibilities
-- `m::AbstractArray{Bool, 3}`: binary mask of magnetic susceptibilities
-- `c::Union{NTuple{3}, SVector{3}}`: center of source to add; c[i] ∈ (1:size(x,i))
-- `r::Real`: radius of source to add
-- `χ::Real`: magnetic susceptibility of new source in ppm
-- `bdir::NTuple{3, Real}`: unit vector of B-field direction (x, y, z)
-- `vsz::NTuple{3, Real}`: voxel size
-
-### Returns
-- `x`
-- `f`
-- `m`
-
-### References
-    [1] Cheng, Y. C. N., Neelavalli, J., & Haacke, E. M. (2009). Limitations of
-    calculating field distributions and magnetic susceptibilities in MRI using
-    a Fourier based method. Physics in Medicine & Biology, 54(5), 1169.
-"""
-function sphere!(
-    x::AbstractArray{T1, 3},
-    f::AbstractArray{T2, 3},
-    m::AbstractArray{Bool, 3},
-    c::Union{NTuple{3}, SVector{3}},
-    r::Real,
-    χ::Real,
-    bdir::NTuple{3, Real},
-    vsz::NTuple{3, Real},
-) where {T1, T2}
-    size(x) == size(f) || throw(DimensionMismatch())
-    size(f) == size(m) || throw(DimensionMismatch())
-    checkbounds(Bool, x, CartesianIndex(c...)) || throw(DimensionMismatch())
-
-    sz = size(x)
-    nx, ny, nz = sz
-
-    C = SVector{3, Float64}(c)
-    B = SVector{3, Float64}(bdir)
-    ns = SVector{3, Float64}(sz)
-    dx = SVector{3, Float64}(vsz)
-
-    if !(norm(B) ≈ 1)
-        @warn "B direction vector is not normalized. Normalizing..."
-        B = B ./ norm(B)
-    end
-
-    X, Y, Z = _grid(C, ns, dx, sz)
-    C = C .+ first.((X, Y, Z))
-
-    chi = Float64(χ)
-    r2 = Float64(r^2)
-    a = Float64(1//3*χ*r^3)
-
-    @batch for k in 1:nz
-        @inbounds for j in 1:ny
-            @fastpow for i in 1:nx
-                R = SVector{3, Float64}(X[i], Y[j], Z[k]) - C
-                R2 = R⋅R
-
-                if R2 <= r2 # inside sphere, field = 0
-                    x[i,j,k] += chi
-                    m[i,j,k] = true
-                else        # outside sphere, compute field
-                    Rz = R⋅B
-                    _f = (3*Rz^2-R2)/sqrt(R2^5)
-                    f[i,j,k] = muladd(a, _f, f[i,j,k])
-                end
-            end
-        end
-    end
-
-    return x, f, m
-end
-
-
-"""
-    cylinder!(x, f, m, c, r, θ, ϕ, χ, bdir, vsz)
-
-Add an infinite cylinder to grid.
-
-### Arguments
-- `x::AbstractArray{T1, 3}`: magnetic susceptibilities
-- `f::AbstractArray{T2, 3}`: local field generated by magnetic susceptibilities
-- `m::AbstractArray{Bool, 3}`: binary mask of magnetic susceptibilities
-- `c::Union{NTuple{3}, SVector{3}}`: center of source to add; c[i] ∈ (1:size(x,i))
-- `r::Real`: radius of source to add
-- `θ::Real`: polar angle wrt B-field of cylinder
-- `ϕ::Real`: azimuthal angle wrt B-field of cylinder
-- `χ::Real`: magnetic susceptibility of new source in ppm
-- `bdir::NTuple{3, Real}`: unit vector of B-field direction (x, y, z)
-- `vsz::NTuple{3, Real}`: voxel size
-
-### Returns
-- `x`
-- `f`
-- `m`
-
-### References
-    [1] Cheng, Y. C. N., Neelavalli, J., & Haacke, E. M. (2009). Limitations of
-    calculating field distributions and magnetic susceptibilities in MRI using
-    a Fourier based method. Physics in Medicine & Biology, 54(5), 1169.
-"""
-function cylinder!(
-    x::AbstractArray{T1, 3},
-    f::AbstractArray{T2, 3},
-    m::AbstractArray{Bool, 3},
-    c::Union{NTuple{3}, SVector{3}},
-    r::Real,
-    θ::Real,
-    ϕ::Real,
-    χ::Real,
-    bdir::NTuple{3, Real},
-    vsz::NTuple{3, Real},
-) where {T1, T2}
-    size(x) == size(f) || throw(DimensionMismatch())
-    size(f) == size(m) || throw(DimensionMismatch())
-    checkbounds(Bool, x, CartesianIndex(c...)) || throw(DimensionMismatch())
-
-    sz = size(x)
-    nx, ny, nz = sz
-
-    C = SVector{3, Float64}(c)
-    ns = SVector{3, Float64}(sz)
-    dx = SVector{3, Float64}(vsz)
-
-    if !(norm(bdir) ≈ 1)
-        @warn "B direction vector is not normalized. Normalizing..."
-        bdir = bdir ./ norm(bdir)
-    end
-
-    X, Y, Z = _grid(C, ns, dx, sz)
-
-    bz = Float64(bdir[3])
-    theta = Float64(θ)
-    phi = Float64(ϕ)
-    chi = Float64(χ)
-    ρ = Float64(r)
-
-    B = cos(acos(bz) - theta)
-    F = 1/6 * chi * (3*B^2 - 1) # constant field inside cylinder
-    a = 0.5*chi * ρ^2 * sin(acos(B))^2
-
-    e = SVector{3, Float64}(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta))
-
-    @batch for k in 1:nz
-        @inbounds for j in 1:ny
-            @fastpow for i in 1:nx
-                R = SVector{3, Float64}(X[i], Y[j], Z[k])
-
-                # mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
-                if norm(R×(R-e)) <= ρ # inside cylinder, constant field
-                    f[I] += F
-                    x[I] += chi
-                    m[I] = true
-                else                     # outside cylinder, compute field
-                    if iszero(theta)
-                        j = SVector{3, Float64}(0, 1, 0)
-                        i = SVector{3, Float64}(1, 0, 0)
-                    else
-                        # j = k×z, z = (0, 0, 1), k = e
-                        j = SVector{3, Float64}(e[2], -e[1], 0)
-                        j = j .* inv(norm(j))
-                        i = j×e
-                    end
-                    Rx2 = (R⋅i)^2
-                    Ry2 = (R⋅j)^2
-                    _f = (Rx2 - Ry2)/(Rx2 + Ry2)^2
-
-                    f[i,j,k] = muladd(a, _f, f[i,j,k])
-                end
-            end
-        end
-    end
-
-    return x, f, m
-end
-
-
-function _grid(
-    c::Union{NTuple{N}, SVector{N}},
-    sz::Union{NTuple{N}, SVector{N}},
-    dx::Union{NTuple{N}, SVector{N}},
-    ns::Union{NTuple{N}, SVector{N}}
-) where {N}
-    x = dx .* fld.(sz, 2)
-    return ntuple(i -> range(-x[i], step=dx[i], length=ns[i]), Val(N))
-end
diff --git a/src/utils/multi_echo.jl b/src/utils/multi_echo.jl
index da8774e..b587cc5 100644
--- a/src/utils/multi_echo.jl
+++ b/src/utils/multi_echo.jl
@@ -35,10 +35,11 @@ function fit_echo_linear(
     NT > 1 || throw(ArgumentError("data must be multi-echo"))
 
     p = similar(phas, size(phas)[1:N-1])
-    if !phase_offset
-        return fit_echo_linear!(p, phas, W, TEs)
+
+    return if !phase_offset
+        fit_echo_linear!(p, phas, W, TEs)
     else
-        return fit_echo_linear!(p, similar(p), phas, W, TEs)
+        fit_echo_linear!(p, similar(p), phas, W, TEs)
     end
 end
 
@@ -70,8 +71,8 @@ function fit_echo_linear!(
     M > 1 || throw(ArgumentError("array must contain echoes in last dimension"))
     NT > 1 || throw(ArgumentError("data must be multi-echo"))
 
-    size(phas, M) == NT || throw(DimensionMismatch())
-    length(p) == length(phas) ÷ NT || throw(DimensionMismatch())
+    checkshape(axes(p), axes(phas)[1:M-1], (:p, :phas))
+    checkshape((NT,), (size(phas, M),), (:TEs, :phas))
     checkshape(W, phas, (:W, :phas))
 
     vphas = reshape(phas, :, NT)
@@ -81,9 +82,9 @@ function fit_echo_linear!(
     T = promote_type(Tp, Tphas, TW)
     tes = convert.(T, TEs)
 
-    _zeroTp = zero(Tp)
+    zeroTp = zero(Tp)
 
-    @inbounds @batch for I in eachindex(vp)
+    @batch for I in eachindex(vp)
         w = vW[I,1] * vW[I,1] * tes[1]
         den = w * tes[1]
         num = w * vphas[I,1]
@@ -92,7 +93,7 @@ function fit_echo_linear!(
             den = muladd(w, tes[t], den)
             num = muladd(w, vphas[I,t], num)
         end
-        vp[I] = iszero(den) ? _zeroTp : num * inv(den)
+        vp[I] = iszero(den) ? zeroTp : num * inv(den)
     end
 
     return p
@@ -130,9 +131,9 @@ function fit_echo_linear!(
     M > 1 || throw(ArgumentError("array must contain echoes in last dimension"))
     NT > 1 || throw(ArgumentError("data must be multi-echo"))
 
-    size(phas, M) == NT || throw(DimensionMismatch())
-    length(p) == length(phas) ÷ NT || throw(DimensionMismatch())
-    length(p0) == length(phas) ÷ NT || throw(DimensionMismatch())
+    checkshape(axes(p), axes(phas)[1:M-1], (:p, :phas))
+    checkshape(axes(p0), axes(phas)[1:M-1], (:p0, :phas))
+    checkshape((NT,), (size(phas, M),), (:TEs, :phas))
     checkshape(W, phas, (:W, :phas))
 
     vphas = reshape(phas, :, NT)
@@ -143,10 +144,10 @@ function fit_echo_linear!(
     T = promote_type(Tp, Tp0, Tphas, TW)
     tes = convert.(T, TEs)
 
-    _zeroTp = zero(Tp)
-    _zeroTp0 = zero(Tp0)
+    zeroTp = zero(Tp)
+    zeroTp0 = zero(Tp0)
 
-    @inbounds @batch for I in eachindex(vp)
+    @batch for I in eachindex(vp)
         w = vW[I,1]
         x = w * tes[1]
         y = w * vphas[I,1]
@@ -157,8 +158,8 @@ function fit_echo_linear!(
         end
 
         if iszero(w)
-            vp[I] = _zeroTp
-            vp0[I] = _zeroTp0
+            vp[I] = zeroTp
+            vp0[I] = zeroTp0
             continue
         end
 
@@ -183,8 +184,8 @@ function fit_echo_linear!(
             vp[I] = num * inv(den)
             vp0[I] = y - vp[I] * x
         else
-            vp[I] = _zeroTp
-            vp0[I] = _zeroTp0
+            vp[I] = zeroTp
+            vp0[I] = zeroTp0
         end
     end
 
diff --git a/src/utils/poisson_solver/mgpcg.jl b/src/utils/poisson_solver/mgpcg.jl
index ad6e59f..1c0b3bb 100644
--- a/src/utils/poisson_solver/mgpcg.jl
+++ b/src/utils/poisson_solver/mgpcg.jl
@@ -99,13 +99,13 @@ function _mgpcg!(
     maxit::Integer,
     verbose::Bool
 ) where {T<:AbstractFloat}
-    _zero = zero(T)
+    zeroT = zero(T)
     ρ = one(T)
 
     A = M.grids[1].A
     q = M.workspace.x[1]
     r = M.workspace.b[1]
-    p = tfill!(p, _zero)
+    p = tfill!(p, zeroT)
 
     # q = A * x
     # r -= q
@@ -127,7 +127,7 @@ function _mgpcg!(
 
     for i in 1:maxit
         # q = M\r
-        tfill!(q, _zero, A)
+        tfill!(q, zeroT, A)
         for _ in 1:ncycles
             cycle!(M, cycle)
         end
@@ -165,7 +165,7 @@ end
 
 function _compute_ρ(q::AbstractArray{T}, r, A::Poisson2) where {T}
     # ρ = q⋅r
-    @inbounds @batch threadlocal=zero(T)::T for I in A.RI
+    @batch threadlocal=zero(T)::T for I in A.RI
         threadlocal = muladd(q[I], r[I], threadlocal)
     end
     return sum(threadlocal::Vector{T})
@@ -174,7 +174,7 @@ end
 
 function _compute_p!(p, q, β, A::Poisson2)
     # p = q + β * p
-    @inbounds @batch for I in A.RI
+    @batch for I in A.RI
         p[I] *= β
         p[I] += q[I]
     end
@@ -190,7 +190,7 @@ function _compute_q!(q::AbstractArray{T}, A::Poisson2, p) where {T}
     idy2 = convert(T, A.idx2[2])
     idz2 = convert(T, A.idx2[3])
 
-    @inbounds @batch per=thread threadlocal=zero(T)::T for (I, J, K) in A.R8
+    @batch per=thread threadlocal=zero(T)::T for (I, J, K) in A.R8
         for k in K, j in J, i in I
             z1 = p[i,j,k-1]
             y1 = p[i,j-1,k]
@@ -224,7 +224,7 @@ function _compute_x_r!(x::AbstractArray{T}, r, p, q, α, A::Poisson2) where {T}
     # x = x + α * p
     # r = r - α * q
     # res = norm(r)
-    @inbounds @batch threadlocal=zero(T)::T for I in A.RI
+    @batch threadlocal=zero(T)::T for I in A.RI
         x[I] = muladd( α, p[I], x[I])
         r[I] = muladd(-α, q[I], r[I])
         threadlocal = muladd(r[I], r[I], threadlocal)
@@ -238,7 +238,7 @@ function _init_r!(r::AbstractArray{T}, q, A, x) where {T}
     # r -= q
     # res = norm(r)
     q = mul!(q, A, x)
-    @inbounds @batch threadlocal=zero(T)::T for I in eachindex(r)
+    @batch threadlocal=zero(T)::T for I in eachindex(r)
         r[I] -= q[I]
         threadlocal = muladd(r[I], r[I], threadlocal)
     end
diff --git a/src/utils/poisson_solver/multigrid/poisson.jl b/src/utils/poisson_solver/multigrid/poisson.jl
index 35113c6..3bb3648 100644
--- a/src/utils/poisson_solver/multigrid/poisson.jl
+++ b/src/utils/poisson_solver/multigrid/poisson.jl
@@ -68,14 +68,12 @@ function LinearAlgebra.mul!(
     x::AbstractArray{T, 3},
 ) where {T<:AbstractFloat}
     D = convert(T, A.D)
-    idx2 = convert(T, A.idx2[1])
-    idy2 = convert(T, A.idx2[2])
-    idz2 = convert(T, A.idx2[3])
+    idx2, idy2, idz2 = convert.(T, A.idx2)
 
     # @threads ~10% faster here
     # Polyester needs A.interior inside or the loop slows down by 100%???
     # ie. slow down when m = A.interior ... if m[i,j,k] ...
-    @inbounds @batch for (I, J, K) in A.R8
+    @batch for (I, J, K) in A.R8
         for k in K
             for j in J
                 for i in I
@@ -113,14 +111,12 @@ function LinearAlgebra.mul!(
     x::AbstractArray{T, 4},
 ) where {T<:AbstractFloat}
     D = convert(T, A.D)
-    idx2 = convert(T, A.idx2[1])
-    idy2 = convert(T, A.idx2[2])
-    idz2 = convert(T, A.idx2[3])
+    idx2, idy2, idz2 = convert.(T, A.idx2)
 
     # @threads ~10% faster here
     # Polyester needs A.interior inside or the loop slows down by 100%???
     # ie. slow down when m = A.interior ... if m[i,j,k] ...
-    @inbounds for t in axes(x, 4)
+    for t in axes(x, 4)
         xt = @view(x[:,:,:,t])
         d2xt = @view(d2x[:,:,:,t])
         @batch for (I, J, K) in A.R8
@@ -171,7 +167,7 @@ function residual!(
     # @threads ~10% faster here
     # Polyester needs A.interior inside or the loop slows down by 100%???
     # ie. slow down when m = A.interior ... if m[i,j,k] ...
-    @inbounds @batch for (I, J, K) in A.R8
+    @batch for (I, J, K) in A.R8
         for k in K
             for j in J
                 for i in I
@@ -216,7 +212,7 @@ function __norm_residual(
     idy2 = convert(T, -A.idx2[2])
     idz2 = convert(T, -A.idx2[3])
 
-    @inbounds @batch threadlocal=zero(T)::T for (I, J, K) in A.R8
+    @batch threadlocal=zero(T)::T for (I, J, K) in A.R8
         for k in K
             for j in J
                 for i in I
@@ -268,7 +264,7 @@ function boundary_mask!(
     end
 
     # interior
-    @inbounds @batch minbatch=8 for kc in 2:nzc-1
+    @batch minbatch=8 for kc in 2:nzc-1
         k = (kc << 1) - 1
         for jc in 2:nyc-1
             j = (jc << 1) - 1
@@ -297,7 +293,7 @@ function boundary_mask!(
         end
     end
 
-    @inbounds @batch minbatch=8 for k in 2:nz-1
+    @batch minbatch=8 for k in 2:nz-1
         for j in 2:ny-1
             for i in 2:nx-1
                 if m[i,j,k]
@@ -320,7 +316,7 @@ end
 
 function tfill!(x::AbstractArray{T, N}, v, A::Poisson2{<:AbstractFloat, N}) where {T, N}
     vT = convert(T, v)
-    @inbounds @batch minbatch=1024 for I in A.RI
+    @batch minbatch=1024 for I in A.RI
         x[I] = vT
     end
     return x
diff --git a/src/utils/poisson_solver/multigrid/smoothers.jl b/src/utils/poisson_solver/multigrid/smoothers.jl
index 1479085..6470fb6 100644
--- a/src/utils/poisson_solver/multigrid/smoothers.jl
+++ b/src/utils/poisson_solver/multigrid/smoothers.jl
@@ -29,7 +29,7 @@ struct Smoothers{N, T} <: AbstractSmoother
 end
 
 function (s::Smoothers{N})(x, A, b) where {N}
-    @inbounds for i in 1:N
+    for i in 1:N
         x = s.s[i](x, A, b)
     end
     return x
@@ -203,7 +203,7 @@ function rbgs!(
     idz2 = -iD * A.idx2[3]
 
     for s in C
-        @inbounds @batch minbatch=8 for (I, J, K) in R
+        @batch minbatch=8 for (I, J, K) in R
             for k in K
                 ks = Bool(k & 1)
                 for j in J
diff --git a/src/utils/poisson_solver/multigrid/transfer.jl b/src/utils/poisson_solver/multigrid/transfer.jl
index bf801d0..b3911ac 100644
--- a/src/utils/poisson_solver/multigrid/transfer.jl
+++ b/src/utils/poisson_solver/multigrid/transfer.jl
@@ -17,7 +17,7 @@ end
 function restrict!(mc::AbstractArray{Bool, 3}, m::AbstractArray{Bool, 3})
     nxc, nyc, nzc = size(mc)
 
-    @inbounds @batch minbatch=8 for kc in 2:nzc-1
+    @batch minbatch=8 for kc in 2:nzc-1
         k = (kc << 1) - 1
         for jc in 2:nyc-1
             j = (jc << 1) - 1
@@ -51,7 +51,7 @@ function restrict!(
     a2 = convert(T, 2//64)
     a3 = convert(T, 1//64)
 
-    @inbounds @batch minbatch=8 for (Ic, Jc, Kc) in Ac.R27
+    @batch minbatch=8 for (Ic, Jc, Kc) in Ac.R27
         for kc in Kc
             k = (kc << 1) - 1
             for jc in Jc
@@ -134,7 +134,7 @@ function prolong!(
     nxc, nyc, nzc = size(xc)
     m = A.interior
 
-    @inbounds @batch for kc in 1:nzc-1
+    @batch for kc in 1:nzc-1
         k = (kc << 1) - 1
         for jc in 1:nyc-1
             j = (jc << 1) - 1
@@ -209,7 +209,7 @@ function correct_prolong!(
     nxc, nyc, nzc = size(xc)
     m = A.interior
 
-    @inbounds @batch for kc in 1:nzc-1
+    @batch for kc in 1:nzc-1
         k = (kc << 1) - 1
         for jc in 1:nyc-1
             j = (jc << 1) - 1
diff --git a/src/utils/poisson_solver/poisson_solver.jl b/src/utils/poisson_solver/poisson_solver.jl
index 25538e9..f4a36fe 100644
--- a/src/utils/poisson_solver/poisson_solver.jl
+++ b/src/utils/poisson_solver/poisson_solver.jl
@@ -24,6 +24,7 @@ function solve_poisson_mgpcg!(
     kwargs...
 ) where {T, N}
     N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
+
     checkshape(u, d2u, (:u, :d2u))
     checkshape(axes(mask), axes(u)[1:3], (:mask, :u))
 
@@ -90,12 +91,11 @@ function solve_poisson_dct!(
     dx::NTuple{3, Real}
 ) where {N}
     N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
+
     checkshape(u, d2u, (:u, :d2u))
 
     nx, ny, nz = size(d2u)[1:3]
-    idx2 = inv(dx[1].*dx[1])
-    idy2 = inv(dx[2].*dx[2])
-    idz2 = inv(dx[3].*dx[3])
+    idx2, idy2, idz2 = inv.(dx.*dx)
 
     # extreme slowdown for certain sizes with lots of threads
     # even worse for in-place, ie dct!
@@ -109,7 +109,7 @@ function solve_poisson_dct!(
     Y = [2*(cospi(j)-1)*idy2 for j in range(0, step=1/ny, length=ny)]
     Z = [2*(cospi(k)-1)*idz2 for k in range(0, step=1/nz, length=nz)]
 
-    @inbounds for t in axes(d2û, 4)
+    for t in axes(d2û, 4)
         d2ût = @view(d2û[:,:,:,t])
         @batch for k in 1:nz
             for j in 1:ny
@@ -143,12 +143,11 @@ function solve_poisson_fft!(
     dx::NTuple{3, Real}
 ) where {N}
     N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
+
     checkshape(u, d2u, (:u, :d2u))
 
     nx, ny, nz = size(d2u)[1:3]
-    idx2 = inv(dx[1].*dx[1])
-    idy2 = inv(dx[2].*dx[2])
-    idz2 = inv(dx[3].*dx[3])
+    idx2, idy2, idz2 = inv.(dx.*dx)
 
     _rfft = iseven(nx)
     FFTW.set_num_threads(FFTW_NTHREADS[])
@@ -170,7 +169,7 @@ function solve_poisson_fft!(
     Y = [2*(cospi(j)-1)*idy2 for j in fftfreq(size(u, 2), 2)]
     Z = [2*(cospi(k)-1)*idz2 for k in fftfreq(size(u, 3), 2)]
 
-    @inbounds for t in axes(d2û, 4)
+    for t in axes(d2û, 4)
         d2ût = @view(d2û[:,:,:,t])
         @batch for k in 1:nz
             for j in 1:ny
diff --git a/src/utils/r2star.jl b/src/utils/r2star.jl
index a524bc6..47fd439 100644
--- a/src/utils/r2star.jl
+++ b/src/utils/r2star.jl
@@ -38,7 +38,7 @@ function r2star_ll(
         b = tmap(log, transpose(vmag))
         x̂ = ldiv!(A, b)
         @batch for I in eachindex(vr2s)
-            @inbounds vr2s[I] = x̂[1,I]
+            vr2s[I] = x̂[1,I]
         end
 
     else
@@ -55,14 +55,14 @@ function r2star_ll(
             bt = @view(b[t,:])
             mt = @view(vmag[:,t])
             @batch for I in eachindex(R)
-                @inbounds bt[I] = log(mt[R[I]])
+                bt[I] = log(mt[R[I]])
             end
         end
 
         x̂ = ldiv!(A, b)
 
         @batch for I in eachindex(R)
-            @inbounds vr2s[R[I]] = x̂[1,I]
+            vr2s[R[I]] = x̂[1,I]
         end
     end
 
@@ -120,7 +120,7 @@ function r2star_arlo(
     α = convert(T, 3 / (TEs[2]-TEs[1]))
 
     @batch for I in eachindex(vr2s)
-        @inbounds if mask === nothing || mask[I]
+        if mask === nothing || mask[I]
             m0 = vmag[I,1]
             m1 = vmag[I,2]
             m2 = vmag[I,3]
@@ -227,7 +227,7 @@ function r2star_crsi(
     @batch for I in eachindex(vr2s)
         if mask === nothing || mask[I]
             den = β
-            @inbounds for t in 1:NT-1
+            for t in 1:NT-1
                 p0 = vP[I,t]
                 p1 = vP[I,t+1]
                 p = pow(p0, γ0[1]) * pow(p1, γ1[1])
@@ -321,7 +321,7 @@ function r2star_numart2s(
     α = convert(T, 2*(NT - 1) / (TEs[end] - TEs[1]))
 
     @batch for I in eachindex(vr2s)
-        @inbounds if mask === nothing || mask[I]
+        if mask === nothing || mask[I]
             den = vmag[I,1]
             for t in 2:NT-1
                 den += vmag[I,t] + vmag[I,t]
diff --git a/src/utils/utils.jl b/src/utils/utils.jl
index b8fa9bc..3a71e65 100644
--- a/src/utils/utils.jl
+++ b/src/utils/utils.jl
@@ -163,15 +163,7 @@ function padarray!(
         pad == :replicate ? getindex_replicate :
         pad == :symmetric ? getindex_symmetric :
         pad == :reflect   ? getindex_reflect :
-        throw(ArgumentError(
-            "pad must be one of " *
-            ":fill, "      *
-            ":circular, "  *
-            ":replicate, " *
-            ":symmetric, " *
-            ":reflect, "   *
-            "got :$(pad)"
-        ))
+        checkopts(pad, (:fill, :circular, :replicate, :symmetric, :reflect), :pad)
 
     return _padarray_kernel!(xp, x, getindex_pad)
 end
@@ -186,9 +178,9 @@ function _padarray_kernel!(xp::AbstractArray, x::AbstractArray, getindex_pad)
     @batch for Ip in CartesianIndices(xp)
         I = Ip - ΔI
         if any(map(∉, I.I, ax))
-            @inbounds xp[Ip] = getindex_pad(x, I, lo, hi)
+            xp[Ip] = getindex_pad(x, I, lo, hi)
         else
-            @inbounds xp[Ip] = x[I]
+            xp[Ip] = x[I]
         end
     end
 
@@ -399,7 +391,7 @@ function erode_mask!(
         @batch for k in 1+t:nz-t
             for j in 1+t:ny-t
                 for i in 1+t:nx-t
-                    @inbounds m1[i,j,k] = __erode_kernel(m0, i, j, k)
+                    m1[i,j,k] = __erode_kernel(m0, i, j, k)
                 end
             end
         end
@@ -469,7 +461,7 @@ function psf2otf(
     else
         _kp = tfill!(similar(k, sz), zero(T))
         @batch minbatch=1024 for I in CartesianIndices(k)
-            @inbounds _kp[I] = k[I]
+            _kp[I] = k[I]
         end
     end
 
@@ -728,3 +720,23 @@ function checkshape(
     end
     return nothing
 end
+
+
+checkopts(::Type{Bool}, o::T, opts::NTuple{N, T}) where {N, T} = o ∈ opts
+
+function checkopts(
+    o::T,
+    opts::NTuple{N, T},
+    var::Union{Symbol, AbstractString} = :x
+) where {N, T}
+    if !checkopts(Bool, o, opts)
+        pp = map((o, opts...)) do x
+            x isa Symbol ? ":$x" :
+            x isa AbstractString ? "\"$x\"" : "$x"
+        end
+        s1 = first(pp)
+        s2 = join(Base.tail(pp), ", ", " ")
+        throw(ArgumentError("$var must be one of $s2, got $s1"))
+    end
+    return nothing
+end