diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl
index 31d3d6edc1..e01b2ff0f5 100644
--- a/src/Operators/finitedifference.jl
+++ b/src/Operators/finitedifference.jl
@@ -367,7 +367,8 @@ Base.@propagate_inbounds function stencil_interior(
 )
     a⁺ = getidx(space, arg, loc, idx + half, hidx)
     a⁻ = getidx(space, arg, loc, idx - half, hidx)
-    RecursiveApply.rdiv(a⁺ ⊞ a⁻, 2)
+    FT = Spaces.undertype(space)
+    RecursiveApply.rmul(a⁺ ⊞ a⁻, FT(0.5))
 end
 
 boundary_width(::InterpolateF2C, ::AbstractBoundaryCondition) = 0
@@ -423,7 +424,8 @@ Base.@propagate_inbounds function stencil_interior(
 )
     a⁺ = getidx(space, arg, loc, idx + half, hidx)
     a⁻ = getidx(space, arg, loc, idx - half, hidx)
-    RecursiveApply.rdiv(a⁺ ⊞ a⁻, 2)
+    FT = Spaces.undertype(space)
+    RecursiveApply.rmul(a⁺ ⊞ a⁻, FT(0.5))
 end
 boundary_width(::InterpolateC2F, ::AbstractBoundaryCondition) = 1
 
@@ -467,7 +469,8 @@ Base.@propagate_inbounds function stencil_left_boundary(
         getidx(space, bc.val, loc, nothing, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    a⁺ ⊟ RecursiveApply.rdiv(v₃, 2)
+    FT = Spaces.undertype(space)
+    a⁺ ⊟ RecursiveApply.rmul(v₃, FT(0.5))
 end
 Base.@propagate_inbounds function stencil_right_boundary(
     ::InterpolateC2F,
@@ -484,7 +487,8 @@ Base.@propagate_inbounds function stencil_right_boundary(
         getidx(space, bc.val, loc, nothing, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    a⁻ ⊞ RecursiveApply.rdiv(v₃, 2)
+    FT = Spaces.undertype(space)
+    a⁻ ⊞ RecursiveApply.rmul(v₃, FT(0.5))
 end
 
 Base.@propagate_inbounds function stencil_left_boundary(
@@ -1180,7 +1184,8 @@ Base.@propagate_inbounds function stencil_left_boundary(
         getidx(space, bc.val, loc, nothing, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    a⁺ ⊟ RecursiveApply.rdiv(v₃, 2)
+    FT = Spaces.undertype(space)
+    a⁺ ⊟ RecursiveApply.rmul(v₃, FT(0.5))
 end
 Base.@propagate_inbounds function stencil_right_boundary(
     ::WeightedInterpolateC2F,
@@ -1198,7 +1203,8 @@ Base.@propagate_inbounds function stencil_right_boundary(
         getidx(space, bc.val, loc, nothing, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    a⁻ ⊞ RecursiveApply.rdiv(v₃, 2)
+    FT = Spaces.undertype(space)
+    a⁻ ⊞ RecursiveApply.rmul(v₃, FT(0.5))
 end
 
 Base.@propagate_inbounds function stencil_left_boundary(
@@ -1916,7 +1922,8 @@ Base.@propagate_inbounds function stencil_interior(
         getidx(space, velocity, loc, idx, hidx),
         Geometry.LocalGeometry(space, idx, hidx),
     )
-    ∂θ₃ = RecursiveApply.rdiv(θ⁺ ⊟ θ⁻, 2)
+    FT = Spaces.undertype(space)
+    ∂θ₃ = RecursiveApply.rmul(θ⁺ ⊟ θ⁻, FT(0.5))
     return w³ ⊠ ∂θ₃
 end
 boundary_width(::AdvectionF2F, ::AbstractBoundaryCondition) = 1
@@ -1986,7 +1993,8 @@ Base.@propagate_inbounds function stencil_interior(
     )
     ∂θ₃⁺ = θ⁺ ⊟ θ
     ∂θ₃⁻ = θ ⊟ θ⁻
-    return RecursiveApply.rdiv((w³⁺ ⊠ ∂θ₃⁺) ⊞ (w³⁻ ⊠ ∂θ₃⁻), 2)
+    FT = Spaces.undertype(space)
+    return RecursiveApply.rmul((w³⁺ ⊠ ∂θ₃⁺) ⊞ (w³⁻ ⊠ ∂θ₃⁻), FT(0.5))
 end
 
 boundary_width(::AdvectionC2C, ::AbstractBoundaryCondition) = 1
@@ -2014,7 +2022,8 @@ Base.@propagate_inbounds function stencil_left_boundary(
     )
     ∂θ₃⁺ = θ⁺ ⊟ θ
     ∂θ₃⁻ = 2 ⊠ (θ ⊟ θ⁻)
-    return RecursiveApply.rdiv((w³⁺ ⊠ ∂θ₃⁺) ⊞ (w³⁻ ⊠ ∂θ₃⁻), 2)
+    FT = Spaces.undertype(space)
+    return RecursiveApply.rmul((w³⁺ ⊠ ∂θ₃⁺) ⊞ (w³⁻ ⊠ ∂θ₃⁻), FT(0.5))
 end
 Base.@propagate_inbounds function stencil_right_boundary(
     ::AdvectionC2C,
@@ -2040,7 +2049,8 @@ Base.@propagate_inbounds function stencil_right_boundary(
     )
     ∂θ₃⁺ = 2 ⊠ (θ⁺ ⊟ θ)
     ∂θ₃⁻ = θ ⊟ θ⁻
-    return RecursiveApply.rdiv((w³⁺ ⊠ ∂θ₃⁺) ⊞ (w³⁻ ⊠ ∂θ₃⁻), 2)
+    FT = Spaces.undertype(space)
+    return RecursiveApply.rmul((w³⁺ ⊠ ∂θ₃⁺) ⊞ (w³⁻ ⊠ ∂θ₃⁻), FT(0.5))
 end
 
 Base.@propagate_inbounds function stencil_left_boundary(
diff --git a/src/Operators/numericalflux.jl b/src/Operators/numericalflux.jl
index f0d9c0e3ca..8818b05622 100644
--- a/src/Operators/numericalflux.jl
+++ b/src/Operators/numericalflux.jl
@@ -91,10 +91,14 @@ end
 function (fn::RusanovNumericalFlux)(normal, argvals⁻, argvals⁺)
     y⁻ = argvals⁻[1]
     y⁺ = argvals⁺[1]
-    Favg =
-        RecursiveApply.rdiv(fn.fluxfn(argvals⁻...) ⊞ fn.fluxfn(argvals⁺...), 2)
     λ = max(fn.wavespeedfn(argvals⁻...), fn.wavespeedfn(argvals⁺...))
-    return RecursiveApply.rmap(f -> f' * normal, Favg) ⊞ (λ / 2) ⊠ (y⁻ ⊟ y⁺)
+    FT = typeof(λ)
+    Favg = RecursiveApply.rmul(
+        fn.fluxfn(argvals⁻...) ⊞ fn.fluxfn(argvals⁺...),
+        FT(0.5),
+    )
+    return RecursiveApply.rmap(f -> f' * normal, Favg) ⊞
+           (FT(0.5) * λ) ⊠ (y⁻ ⊟ y⁺)
 end
 
 
diff --git a/src/Operators/operator2stencil.jl b/src/Operators/operator2stencil.jl
index 111a7e4f25..b5fa35d823 100644
--- a/src/Operators/operator2stencil.jl
+++ b/src/Operators/operator2stencil.jl
@@ -107,8 +107,11 @@ function stencil_interior(
     hidx,
     arg,
 )
-    val⁻ = RecursiveApply.rdiv(getidx(space, arg, loc, idx - half, hidx), 2)
-    val⁺ = RecursiveApply.rdiv(getidx(space, arg, loc, idx + half, hidx), 2)
+    FT = Spaces.undertype(space)
+    val⁻ =
+        RecursiveApply.rmul(getidx(space, arg, loc, idx - half, hidx), FT(0.5))
+    val⁺ =
+        RecursiveApply.rmul(getidx(space, arg, loc, idx + half, hidx), FT(0.5))
     return StencilCoefs{-half, half}((val⁻, val⁺))
 end
 function stencil_left_boundary(
@@ -231,8 +234,9 @@ function stencil_interior(
     θ⁻ = getidx(space, arg, loc, idx - 1, hidx)
     θ = getidx(space, arg, loc, idx, hidx)
     θ⁺ = getidx(space, arg, loc, idx + 1, hidx)
-    val⁻ = RecursiveApply.rdiv(w³⁻ ⊠ (θ ⊟ θ⁻), 2)
-    val⁺ = RecursiveApply.rdiv(w³⁺ ⊠ (θ⁺ ⊟ θ), 2)
+    FT = Spaces.undertype(space)
+    val⁻ = RecursiveApply.rmul(w³⁻ ⊠ (θ ⊟ θ⁻), FT(0.5))
+    val⁺ = RecursiveApply.rmul(w³⁺ ⊠ (θ⁺ ⊟ θ), FT(0.5))
     return StencilCoefs{-half, half}((val⁻, val⁺))
 end
 # TODO: find out why using Base.@propagate_inbounds blows up compilation time
@@ -258,7 +262,8 @@ function stencil_left_boundary(
     θ = getidx(space, arg, loc, idx, hidx)
     θ⁺ = getidx(space, arg, loc, idx + 1, hidx)
     val⁻ = w³⁻ ⊠ (θ ⊟ θ⁻)
-    val⁺ = RecursiveApply.rdiv(w³⁺ ⊠ (θ⁺ ⊟ θ), 2)
+    FT = Spaces.undertype(space)
+    val⁺ = RecursiveApply.rmul(w³⁺ ⊠ (θ⁺ ⊟ θ), FT(0.5))
     return StencilCoefs{-half, half}((val⁻, val⁺))
 end
 # TODO: find out why using Base.@propagate_inbounds blows up compilation time
@@ -283,7 +288,8 @@ function stencil_right_boundary(
     θ⁻ = getidx(space, arg, loc, idx - 1, hidx)
     θ = getidx(space, arg, loc, idx, hidx)
     θ⁺ = getidx(space, bc.val, loc, nothing, hidx)
-    val⁻ = RecursiveApply.rdiv(w³⁻ ⊠ (θ ⊟ θ⁻), 2)
+    FT = Spaces.undertype(space)
+    val⁻ = RecursiveApply.rmul(w³⁻ ⊠ (θ ⊟ θ⁻), FT(0.5))
     val⁺ = w³⁺ ⊠ (θ⁺ ⊟ θ)
     return StencilCoefs{-half, half}((val⁻, val⁺))
 end
diff --git a/src/Operators/spectralelement.jl b/src/Operators/spectralelement.jl
index d0e953c06c..d85743a56a 100644
--- a/src/Operators/spectralelement.jl
+++ b/src/Operators/spectralelement.jl
@@ -734,7 +734,7 @@ function apply_operator(op::Divergence{(1,)}, space, slabidx, arg)
     @inbounds for i in 1:Nq
         ij = CartesianIndex((i,))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        out[i] = RecursiveApply.rdiv(out[i], local_geometry.J)
+        out[i] = RecursiveApply.rmul(out[i], local_geometry.invJ)
     end
     return Field(SArray(out), space)
 end
@@ -777,7 +777,7 @@ Base.@propagate_inbounds function apply_operator(
     @inbounds for j in 1:Nq, i in 1:Nq
         ij = CartesianIndex((i, j))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        out[i, j] = RecursiveApply.rdiv(out[i, j], local_geometry.J)
+        out[i, j] = RecursiveApply.rmul(out[i, j], local_geometry.invJ)
     end
     return Field(SArray(out), space)
 end
@@ -846,7 +846,7 @@ Base.@propagate_inbounds function operator_evaluate(
     for k in 1:Nq
         DJv = DJv ⊞ D[j, k] ⊠ Jv²[i, k, vt]
     end
-    return RecursiveApply.rdiv(DJv, local_geometry.J)
+    return RecursiveApply.rmul(DJv, local_geometry.invJ)
 end
 
 """
@@ -1213,7 +1213,7 @@ function apply_operator(op::WeakGradient{(1,)}, space, slabidx, arg)
     @inbounds for i in 1:Nq
         ij = CartesianIndex((i,))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        W = local_geometry.WJ / local_geometry.J
+        W = local_geometry.WJ * local_geometry.invJ
         Wx = W ⊠ get_node(space, arg, ij, slabidx)
         for ii in 1:Nq
             Dᵀ₁Wf = Geometry.Covariant1Vector(D[i, ii]) ⊗ Wx
@@ -1223,7 +1223,7 @@ function apply_operator(op::WeakGradient{(1,)}, space, slabidx, arg)
     @inbounds for i in 1:Nq
         ij = CartesianIndex((i,))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        W = local_geometry.WJ / local_geometry.J
+        W = local_geometry.WJ * local_geometry.invJ
         out[i] = RecursiveApply.rdiv(out[i], W)
     end
     return Field(SArray(out), space)
@@ -1242,7 +1242,7 @@ function apply_operator(op::WeakGradient{(1, 2)}, space, slabidx, arg)
     @inbounds for j in 1:Nq, i in 1:Nq
         ij = CartesianIndex((i, j))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        W = local_geometry.WJ / local_geometry.J
+        W = local_geometry.WJ * local_geometry.invJ
         Wx = W ⊠ get_node(space, arg, ij, slabidx)
         for ii in 1:Nq
             Dᵀ₁Wf = Geometry.Covariant12Vector(D[i, ii], zero(eltype(D))) ⊗ Wx
@@ -1256,7 +1256,7 @@ function apply_operator(op::WeakGradient{(1, 2)}, space, slabidx, arg)
     @inbounds for j in 1:Nq, i in 1:Nq
         ij = CartesianIndex((i, j))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        W = local_geometry.WJ / local_geometry.J
+        W = local_geometry.WJ * local_geometry.invJ
         out[i, j] = RecursiveApply.rdiv(out[i, j], W)
     end
     return Field(SArray(out), space)
@@ -1287,7 +1287,7 @@ Base.@propagate_inbounds function operator_fill_shmem!(
 )
     vt = threadIdx().z
     local_geometry = get_local_geometry(space, ij, slabidx)
-    W = local_geometry.WJ / local_geometry.J
+    W = local_geometry.WJ * local_geometry.invJ
     i, j = ij.I
     Wf[i, j, vt] = W ⊠ arg
 end
@@ -1308,7 +1308,7 @@ Base.@propagate_inbounds function operator_evaluate(
     D = Quadratures.differentiation_matrix(FT, QS)
 
     local_geometry = get_local_geometry(space, ij, slabidx)
-    W = local_geometry.WJ / local_geometry.J
+    W = local_geometry.WJ * local_geometry.invJ
 
     Dᵀ₁Wf = D[1, i] ⊠ Wf[1, j, vt]
     Dᵀ₂Wf = D[1, j] ⊠ Wf[i, 1, vt]
@@ -1423,7 +1423,7 @@ function apply_operator(op::Curl{(1,)}, space, slabidx, arg)
     @inbounds for i in 1:Nq
         ij = CartesianIndex((i,))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        out[i] = RecursiveApply.rdiv(out[i], local_geometry.J)
+        out[i] = RecursiveApply.rmul(out[i], local_geometry.invJ)
     end
     return Field(SArray(out), space)
 end
@@ -1511,7 +1511,7 @@ Base.@propagate_inbounds function operator_evaluate(
             D₂v₁ = D₂v₁ ⊞ D[j, k] ⊠ v₁[i, k, vt]
         end
         return Geometry.Contravariant3Vector(
-            RecursiveApply.rdiv(D₁v₂ ⊟ D₂v₁, local_geometry.J),
+            RecursiveApply.rmul(D₁v₂ ⊟ D₂v₁, local_geometry.invJ),
         )
     elseif length(work) == 1
         (v₃,) = work
@@ -1522,8 +1522,8 @@ Base.@propagate_inbounds function operator_evaluate(
             D₂v₃ = D₂v₃ ⊞ D[j, k] ⊠ v₃[i, k, vt]
         end
         return Geometry.Contravariant12Vector(
-            RecursiveApply.rdiv(D₂v₃, local_geometry.J),
-            ⊟(RecursiveApply.rdiv(D₁v₃, local_geometry.J)),
+            RecursiveApply.rmul(D₂v₃, local_geometry.invJ),
+            ⊟(RecursiveApply.rmul(D₁v₃, local_geometry.invJ)),
         )
     else
         v₁, v₂, v₃ = work
@@ -1538,9 +1538,9 @@ Base.@propagate_inbounds function operator_evaluate(
             D₂v₃ = D₂v₃ ⊞ D[j, k] ⊠ v₃[i, k, vt]
         end
         return Geometry.Contravariant123Vector(
-            RecursiveApply.rdiv(D₂v₃, local_geometry.J),
-            ⊟(RecursiveApply.rdiv(D₁v₃, local_geometry.J)),
-            RecursiveApply.rdiv(D₁v₂ ⊟ D₂v₁, local_geometry.J),
+            RecursiveApply.rmul(D₂v₃, local_geometry.invJ),
+            ⊟(RecursiveApply.rmul(D₁v₃, local_geometry.invJ)),
+            RecursiveApply.rmul(D₁v₂ ⊟ D₂v₁, local_geometry.invJ),
         )
     end
 end
@@ -1621,7 +1621,7 @@ function apply_operator(op::Curl{(1, 2)}, space, slabidx, arg)
     @inbounds for j in 1:Nq, i in 1:Nq
         ij = CartesianIndex((i, j))
         local_geometry = get_local_geometry(space, ij, slabidx)
-        out[i, j] = RecursiveApply.rdiv(out[i, j], local_geometry.J)
+        out[i, j] = RecursiveApply.rmul(out[i, j], local_geometry.invJ)
     end
     return Field(SArray(out), space)
 end
@@ -1684,7 +1684,7 @@ function apply_operator(op::WeakCurl{(1,)}, space, slabidx, arg)
             ij = CartesianIndex((i,))
             local_geometry = get_local_geometry(space, ij, slabidx)
             v = get_node(space, arg, ij, slabidx)
-            W = local_geometry.WJ / local_geometry.J
+            W = local_geometry.WJ * local_geometry.invJ
             Wv₃ = W ⊠ Geometry.covariant3(v, local_geometry)
             for ii in 1:Nq
                 Dᵀ₁Wv₃ = D[i, ii] ⊠ Wv₃
@@ -1696,7 +1696,7 @@ function apply_operator(op::WeakCurl{(1,)}, space, slabidx, arg)
             ij = CartesianIndex((i,))
             local_geometry = get_local_geometry(space, ij, slabidx)
             v = get_node(space, arg, ij, slabidx)
-            W = local_geometry.WJ / local_geometry.J
+            W = local_geometry.WJ * local_geometry.invJ
             Wv₂ = W ⊠ Geometry.covariant2(v, local_geometry)
             for ii in 1:Nq
                 Dᵀ₁Wv₂ = D[i, ii] ⊠ Wv₂
@@ -1708,7 +1708,7 @@ function apply_operator(op::WeakCurl{(1,)}, space, slabidx, arg)
             ij = CartesianIndex((i,))
             local_geometry = get_local_geometry(space, ij, slabidx)
             v = get_node(space, arg, ij, slabidx)
-            W = local_geometry.WJ / local_geometry.J
+            W = local_geometry.WJ * local_geometry.invJ
             Wv₃ = W ⊠ Geometry.covariant3(v, local_geometry)
             Wv₂ = W ⊠ Geometry.covariant2(v, local_geometry)
             for ii in 1:Nq
@@ -1745,7 +1745,7 @@ function apply_operator(op::WeakCurl{(1, 2)}, space, slabidx, arg)
             ij = CartesianIndex((i, j))
             local_geometry = get_local_geometry(space, ij, slabidx)
             v = get_node(space, arg, ij, slabidx)
-            W = local_geometry.WJ / local_geometry.J
+            W = local_geometry.WJ * local_geometry.invJ
             Wv₁ = W ⊠ Geometry.covariant1(v, local_geometry)
             for jj in 1:Nq
                 Dᵀ₂Wv₁ = D[j, jj] ⊠ Wv₁
@@ -1764,7 +1764,7 @@ function apply_operator(op::WeakCurl{(1, 2)}, space, slabidx, arg)
             ij = CartesianIndex((i, j))
             local_geometry = get_local_geometry(space, ij, slabidx)
             v = get_node(space, arg, ij, slabidx)
-            W = local_geometry.WJ / local_geometry.J
+            W = local_geometry.WJ * local_geometry.invJ
             Wv₃ = W ⊠ Geometry.covariant3(v, local_geometry)
             for ii in 1:Nq
                 Dᵀ₁Wv₃ = D[i, ii] ⊠ Wv₃
@@ -1784,7 +1784,7 @@ function apply_operator(op::WeakCurl{(1, 2)}, space, slabidx, arg)
             ij = CartesianIndex((i, j))
             local_geometry = get_local_geometry(space, ij, slabidx)
             v = get_node(space, arg, ij, slabidx)
-            W = local_geometry.WJ / local_geometry.J
+            W = local_geometry.WJ * local_geometry.invJ
             Wv₁ = W ⊠ Geometry.covariant1(v, local_geometry)
             Wv₂ = W ⊠ Geometry.covariant2(v, local_geometry)
             Wv₃ = W ⊠ Geometry.covariant3(v, local_geometry)
@@ -1862,7 +1862,7 @@ Base.@propagate_inbounds function operator_fill_shmem!(
     vt = threadIdx().z
     i, j = ij.I
     local_geometry = get_local_geometry(space, ij, slabidx)
-    W = local_geometry.WJ / local_geometry.J
+    W = local_geometry.WJ * local_geometry.invJ
     RT = operator_return_eltype(op, typeof(arg))
     if RT <: Geometry.Contravariant3Vector
         Wv₁, Wv₂ = work
diff --git a/src/Remapping/interpolate_array.jl b/src/Remapping/interpolate_array.jl
index 06fde5df90..fffc9174ba 100644
--- a/src/Remapping/interpolate_array.jl
+++ b/src/Remapping/interpolate_array.jl
@@ -69,7 +69,8 @@ function interpolate_slab_level(
     end
     f_lo = interpolate_slab(field, Fields.SlabIndex(v_lo, h), Is)
     f_hi = interpolate_slab(field, Fields.SlabIndex(v_hi, h), Is)
-    return ((1 - ξ3) * f_lo + (1 + ξ3) * f_hi) / 2
+    FT = Spaces.undertype(axes(field))
+    return FT(0.5) * ((1 - ξ3) * f_lo + (1 + ξ3) * f_hi)
 end
 
 """
diff --git a/src/Spaces/finitedifference.jl b/src/Spaces/finitedifference.jl
index 38302bcba6..a8210dcf39 100644
--- a/src/Spaces/finitedifference.jl
+++ b/src/Spaces/finitedifference.jl
@@ -62,7 +62,7 @@ function FiniteDifferenceSpace{S}(
         # at the moment we use a "discrete Jacobian"
         # ideally we should use the continuous quantity via the derivative of the warp function
         # could we just define this then as deriv on the mesh element coordinates?
-        coord = (coord⁺ + coord⁻) / 2
+        coord = FT(0.5) * (coord⁺ + coord⁻)
         Δcoord = coord⁺ - coord⁻
         J = Δcoord
         WJ = Δcoord
@@ -88,22 +88,22 @@ function FiniteDifferenceSpace{S}(
                 Δcoord⁻ =
                     Geometry.component(face_coordinates[end], 1) -
                     Geometry.component(face_coordinates[end - 1], 1)
-                J = (Δcoord⁺ + Δcoord⁻) / 2
+                J = FT(0.5) * (Δcoord⁺ + Δcoord⁻)
                 WJ = J
             else
                 coord⁺ = Geometry.component(face_coordinates[2], 1)
                 J = coord⁺ - coord
-                WJ = J / 2
+                WJ = FT(0.5) * J
             end
         elseif !Topologies.isperiodic(topology) && i == nface
             # top face
             coord⁻ = Geometry.component(face_coordinates[i - 1], 1)
             J = coord - coord⁻
-            WJ = J / 2
+            WJ = FT(0.5) * J
         else
             coord⁺ = Geometry.component(face_coordinates[i + 1], 1)
             coord⁻ = Geometry.component(face_coordinates[i - 1], 1)
-            J = (coord⁺ - coord⁻) / 2
+            J = FT(0.5) * (coord⁺ - coord⁻)
             WJ = J
         end
         ∂x∂ξ = SMatrix{1, 1}(J)
diff --git a/src/Spaces/spectralelement.jl b/src/Spaces/spectralelement.jl
index 2d7212befb..65dfa52540 100644
--- a/src/Spaces/spectralelement.jl
+++ b/src/Spaces/spectralelement.jl
@@ -118,10 +118,10 @@ function SpectralElementSpace1D(
             vcoords = Topologies.vertex_coordinates(topology, elem)
             x = Geometry.linear_interpolate(vcoords, ξ)
             ∂x∂ξ =
-                (
+                FT(0.5) * (
                     Geometry.component(vcoords[2], 1) -
                     Geometry.component(vcoords[1], 1)
-                ) / 2
+                )
             J = abs(∂x∂ξ)
             WJ = J * quad_weights[i]
             local_geometry_slab[i] = Geometry.LocalGeometry(