diff --git a/examples/hybrid/plane/topo_dc_2dinvariant_rhoe.jl b/examples/hybrid/plane/topo_dc_2dinvariant_rhoe.jl
index 1ae633ace0..06cbe77aa0 100644
--- a/examples/hybrid/plane/topo_dc_2dinvariant_rhoe.jl
+++ b/examples/hybrid/plane/topo_dc_2dinvariant_rhoe.jl
@@ -144,7 +144,12 @@ Y = Fields.FieldVector(Yc = Yc, uₕ = uₕ, w = w)
 energy_0 = sum(Y.Yc.ρe)
 mass_0 = sum(Y.Yc.ρ)
 
-function rhs_invariant!(dY, Y, _, t)
+function rhs_invariant!(dY, Y, p, t)
+
+    (; fρ, ᶠ∇ᵥuₕ, ᶜ∇ᵥw, ᶠ∇ᵥh_tot, ᶜ∇ₕuₕ, ᶠ∇ₕw, ᶜ∇ₕh_tot) = p
+    (; hκ₂∇²uₕ, vκ₂∇²uₕ, hκ₂∇²w, vκ₂∇²w, hκ₂∇²h_tot) = p
+    (; vκ₂∇²h_tot, cE, fu¹, fu³, fu, cω², fω², h_tot) = p
+    (; ce, cI, cT, cp, cuw, cw, u₃_bc, fuₕ, coords) = p
 
     cρ = Y.Yc.ρ # scalar on centers
     fw = Y.w # Covariant3Vector on faces
@@ -172,10 +177,10 @@ function rhs_invariant!(dY, Y, _, t)
         top = Operators.Extrapolate(),
     )
 
-    dρ .= 0 .* cρ
+    @. dρ = 0 * cρ
 
-    cw = If2c.(fw)
-    fuₕ = Ic2f.(cuₕ)
+    @. cw = If2c(fw)
+    @. fuₕ = Ic2f(cuₕ)
     # ==========
     u₁_bc = Fields.level(fuₕ, ClimaCore.Utilities.half)
     gⁱʲ =
@@ -186,25 +191,26 @@ function rhs_invariant!(dY, Y, _, t)
     g13 = gⁱʲ.components.data.:3
     g11 = gⁱʲ.components.data.:1
     g33 = gⁱʲ.components.data.:4
-    gratio = g13 ./ g33
-    u₃_bc = Geometry.Covariant3Vector.(-1 .* gratio .* u₁_bc.components.data.:1)
+    u₁_bc₁ = u₁_bc.components.data.:1
+    @. u₃_bc = Geometry.Covariant3Vector(-1 * g13 / g33 * u₁_bc₁)
     apply_boundary_w =
         Operators.SetBoundaryOperator(bottom = Operators.SetValue(u₃_bc))
     @. fw = apply_boundary_w(fw)
     Spaces.weighted_dss!(fw)
     # ==========
-    cw = If2c.(fw)
-    cuw = Geometry.Covariant13Vector.(cuₕ) .+ Geometry.Covariant13Vector.(cw)
+    @. cw = If2c(fw)
+    @. cuw = Geometry.Covariant13Vector(cuₕ) + Geometry.Covariant13Vector(cw)
 
-    ce = @. cρe / cρ
-    cI = @. ce - Φ(z) - (norm(cuw)^2) / 2
-    cT = @. cI / C_v + T_0
-    cp = @. cρ * R_d * cT
+    @. ce = cρe / cρ
+    @. cI = ce - Φ(z) - (norm(cuw)^2) / 2
+    @. cT = cI / C_v + T_0
+    @. cp = cρ * R_d * cT
 
-    h_tot = @. ce + cp / cρ # Total enthalpy at cell centers
+    @. h_tot = ce + cp / cρ # Total enthalpy at cell centers
 
     ### HYPERVISCOSITY
     # 1) compute hyperviscosity coefficients
+    # Move to cache if necessary
     χe = @. dρe = hwdiv(hgrad(h_tot)) # we store χe in dρe
     χuₕ = @. duₕ = hwgrad(hdiv(cuₕ))
 
@@ -216,10 +222,10 @@ function rhs_invariant!(dY, Y, _, t)
     @. duₕ = -κ₄ * (hwgrad(hdiv(χuₕ)))
 
     # 1) Mass conservation
-    dw .= fw .* 0
+    @. dw = fw * 0
 
     # 1.a) horizontal divergence
-    dρ .-= hdiv.(cρ .* (cuw))
+    @. dρ -= hdiv(cρ * (cuw))
 
     # 1.b) vertical divergence
     vdivf2c = Operators.DivergenceF2C(
@@ -235,9 +241,9 @@ function rhs_invariant!(dY, Y, _, t)
     # negation
 
     # explicit part
-    dρ .-= vdivf2c.(Ic2f.(cρ .* cuₕ))
+    @. dρ -= vdivf2c(Ic2f(cρ * cuₕ))
     # implicit part
-    dρ .-= vdivf2c.(Ic2f.(cρ) .* fw)
+    @. dρ -= vdivf2c(Ic2f(cρ) * fw)
 
     # 2) Momentum equation
 
@@ -249,11 +255,11 @@ function rhs_invariant!(dY, Y, _, t)
         top = Operators.SetCurl(Geometry.Contravariant2Vector(0.0)),
     )
 
-    fω² = hcurl.(fw)
-    fω² .+= vcurlc2f.(cuₕ)
+    @. fω² = hcurl(fw)
+    @. fω² += vcurlc2f(cuₕ)
 
-    cω² = hcurl.(cw) # Compute new center curl
-    cω² .+= If2c.(vcurlc2f.(cuₕ)) # Compute new centerl curl
+    @. cω² = hcurl(cw) # Compute new center curl
+    @. cω² += If2c(vcurlc2f(cuₕ)) # Compute new centerl curl
 
     # Linearly interpolate the horizontal velocities multiplying the curl terms from centers to faces 
     # (i.e., u_f^i = (u_c^{i+1} + u_c^{i})/2)
@@ -262,11 +268,12 @@ function rhs_invariant!(dY, Y, _, t)
     # cross product
     # convert to contravariant
     # these will need to be modified with topography
-    fu =
-        Geometry.Covariant13Vector.(Ic2f.(cuₕ)) .+
-        Geometry.Covariant13Vector.(fw)
-    fu¹ = Geometry.project.(Ref(Geometry.Contravariant1Axis()), fu)
-    fu³ = Geometry.project.(Ref(Geometry.Contravariant3Axis()), fu)
+    @. fu =
+        Geometry.Covariant13Vector(Ic2f(cuₕ)) + Geometry.Covariant13Vector(fw)
+    contra1 = (Geometry.Contravariant1Axis(),)
+    contra3 = (Geometry.Contravariant3Axis(),)
+    @. fu¹ = Geometry.project(contra1, fu)
+    @. fu³ = Geometry.project(contra3, fu)
 
     @. duₕ -= If2c(fω² × fu³)
     @. dw -= fω² × fu¹ # Covariant3Vector on faces
@@ -278,7 +285,7 @@ function rhs_invariant!(dY, Y, _, t)
     )
     @. dw -= vgradc2f(cp) / Ic2f(cρ)
 
-    cE = @. (norm(cuw)^2) / 2 + Φ(z)
+    @. cE = (norm(cuw)^2) / 2 + Φ(z)
     @. duₕ -= hgrad(cE)
     @. dw -= vgradc2f(cE)
 
@@ -289,23 +296,23 @@ function rhs_invariant!(dY, Y, _, t)
 
     # Uniform 2nd order diffusion
     ∂c = Operators.GradientF2C()
-    fρ = @. Ic2f(cρ)
+    @. fρ = Ic2f(cρ)
     κ₂ = 0.0 # m^2/s
 
-    ᶠ∇ᵥuₕ = @. vgradc2f(cuₕ.components.data.:1)
-    ᶜ∇ᵥw = @. ∂c(fw.components.data.:1)
-    ᶠ∇ᵥh_tot = @. vgradc2f(h_tot)
+    @. ᶠ∇ᵥuₕ = vgradc2f(cuₕ.components.data.:1)
+    @. ᶜ∇ᵥw = ∂c(fw.components.data.:1)
+    @. ᶠ∇ᵥh_tot = vgradc2f(h_tot)
 
-    ᶜ∇ₕuₕ = @. hgrad(cuₕ.components.data.:1)
-    ᶠ∇ₕw = @. hgrad(fw.components.data.:1)
-    ᶜ∇ₕh_tot = @. hgrad(h_tot)
+    @. ᶜ∇ₕuₕ = hgrad(cuₕ.components.data.:1)
+    @. ᶠ∇ₕw = hgrad(fw.components.data.:1)
+    @. ᶜ∇ₕh_tot = hgrad(h_tot)
 
-    hκ₂∇²uₕ = @. hwdiv(κ₂ * ᶜ∇ₕuₕ)
-    vκ₂∇²uₕ = @. vdivf2c(κ₂ * ᶠ∇ᵥuₕ)
-    hκ₂∇²w = @. hwdiv(κ₂ * ᶠ∇ₕw)
-    vκ₂∇²w = @. vdivc2f(κ₂ * ᶜ∇ᵥw)
-    hκ₂∇²h_tot = @. hwdiv(cρ * κ₂ * ᶜ∇ₕh_tot)
-    vκ₂∇²h_tot = @. vdivf2c(fρ * κ₂ * ᶠ∇ᵥh_tot)
+    @. hκ₂∇²uₕ = hwdiv(κ₂ * ᶜ∇ₕuₕ)
+    @. vκ₂∇²uₕ = vdivf2c(κ₂ * ᶠ∇ᵥuₕ)
+    @. hκ₂∇²w = hwdiv(κ₂ * ᶠ∇ₕw)
+    @. vκ₂∇²w = vdivc2f(κ₂ * ᶜ∇ᵥw)
+    @. hκ₂∇²h_tot = hwdiv(cρ * κ₂ * ᶜ∇ₕh_tot)
+    @. vκ₂∇²h_tot = vdivf2c(fρ * κ₂ * ᶠ∇ᵥh_tot)
 
     dfw = dY.w.components.data.:1
     dcu = dY.uₕ.components.data.:1
@@ -325,14 +332,120 @@ function rhs_invariant!(dY, Y, _, t)
     return dY
 end
 
+function get_cache(Y)
+    cρ = Y.Yc.ρ # scalar on centers
+    fw = Y.w # Covariant3Vector on faces
+    cuₕ = Y.uₕ # Covariant1Vector on centers
+    cρe = Y.Yc.ρe
+
+    z = coords.z
+
+    # 0) update w at the bottom
+    # fw = -g^31 cuₕ/ g^33
+
+    hwdiv = Operators.WeakDivergence()
+    hgrad = Operators.Gradient()
+    hcurl = Operators.Curl()
+
+    If2c = Operators.InterpolateF2C()
+    Ic2f = Operators.InterpolateC2F(
+        bottom = Operators.Extrapolate(),
+        top = Operators.Extrapolate(),
+    )
+
+    fuₕ = Ic2f.(cuₕ)
+    # ==========
+    u₁_bc = Fields.level(fuₕ, ClimaCore.Utilities.half)
+    gⁱʲ =
+        Fields.level(
+            Fields.local_geometry_field(hv_face_space),
+            ClimaCore.Utilities.half,
+        ).gⁱʲ
+    g13 = gⁱʲ.components.data.:3
+    g11 = gⁱʲ.components.data.:1
+    g33 = gⁱʲ.components.data.:4
+
+    u₁_bc₁ = u₁_bc.components.data.:1
+    u₃_bc = @. Geometry.Covariant3Vector(-1 * g13 / g33 * u₁_bc₁)
+    # ==========
+    cw = If2c.(fw)
+    cuw = Geometry.Covariant13Vector.(cuₕ) .+ Geometry.Covariant13Vector.(cw)
+
+    ce = @. cρe / cρ
+    cI = @. ce - Φ(z) - (norm(cuw)^2) / 2
+    cT = @. cI / C_v + T_0
+    cp = @. cρ * R_d * cT
+    h_tot = @. ce + cp / cρ # Total enthalpy at cell centers
+    # 1.b) vertical divergence
+    vdivf2c = Operators.DivergenceF2C(
+        top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
+        bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
+    )
+    vdivc2f = Operators.DivergenceC2F(
+        top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
+        bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
+    )
+
+    fω² = hcurl.(fw)
+    cω² = hcurl.(cw) # Compute new center curl
+
+    # Linearly interpolate the horizontal velocities multiplying the curl terms from centers to faces 
+    # (i.e., u_f^i = (u_c^{i+1} + u_c^{i})/2)
+    # Leave the horizontal curl terms (living on faces) untouched.
+
+    # cross product
+    # convert to contravariant
+    # these will need to be modified with topography
+    fu =
+        Geometry.Covariant13Vector.(Ic2f.(cuₕ)) .+
+        Geometry.Covariant13Vector.(fw)
+    fu¹ = Geometry.project.(Ref(Geometry.Contravariant1Axis()), fu)
+    fu³ = Geometry.project.(Ref(Geometry.Contravariant3Axis()), fu)
+
+    vgradc2f = Operators.GradientC2F(
+        bottom = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
+        top = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
+    )
+
+    cE = @. (norm(cuw)^2) / 2 + Φ(z)
+
+    # Uniform 2nd order diffusion
+    ∂c = Operators.GradientF2C()
+    fρ = @. Ic2f(cρ)
+
+    ᶠ∇ᵥuₕ = @. vgradc2f(cuₕ.components.data.:1)
+    ᶜ∇ᵥw = @. ∂c(fw.components.data.:1)
+    ᶠ∇ᵥh_tot = @. vgradc2f(h_tot)
+
+    ᶜ∇ₕuₕ = @. hgrad(cuₕ.components.data.:1)
+    ᶠ∇ₕw = @. hgrad(fw.components.data.:1)
+    ᶜ∇ₕh_tot = @. hgrad(h_tot)
+
+    hκ₂∇²uₕ = @. hwdiv(ᶜ∇ₕuₕ)
+    vκ₂∇²uₕ = @. vdivf2c(ᶠ∇ᵥuₕ)
+    hκ₂∇²w = @. hwdiv(ᶠ∇ₕw)
+    vκ₂∇²w = @. vdivc2f(ᶜ∇ᵥw)
+    hκ₂∇²h_tot = @. hwdiv(cρ * ᶜ∇ₕh_tot)
+    vκ₂∇²h_tot = @. vdivf2c(fρ * ᶠ∇ᵥh_tot)
+
+    p1 = (; fρ, ᶠ∇ᵥuₕ, ᶜ∇ᵥw, ᶠ∇ᵥh_tot, ᶜ∇ₕuₕ, ᶠ∇ₕw, ᶜ∇ₕh_tot)
+    p2 = (; hκ₂∇²uₕ, vκ₂∇²uₕ, hκ₂∇²w, vκ₂∇²w, hκ₂∇²h_tot)
+    p3 = (; vκ₂∇²h_tot, cE, fu¹, fu³, fu, cω², fω², h_tot)
+    p4 = (; ce, cI, cT, cp, cuw, cw, u₃_bc, fuₕ, coords)
+
+    return merge(p1, p2, p3, p4)
+end
+
+p = get_cache(Y);
+
 dYdt = similar(Y);
-rhs_invariant!(dYdt, Y, nothing, 0.0);
+rhs_invariant!(dYdt, Y, p, 0.0);
 
 # run!
 using OrdinaryDiffEq
 timeend = 3600.0 * 3
 Δt = 0.05
-prob = ODEProblem(rhs_invariant!, Y, (0.0, timeend))
+prob = ODEProblem(rhs_invariant!, Y, (0.0, timeend), p)
 integrator = OrdinaryDiffEq.init(
     prob,
     SSPRK33(),
diff --git a/src/Operators/spectralelement.jl b/src/Operators/spectralelement.jl
index 5af4673ee3..8e2b838d6f 100644
--- a/src/Operators/spectralelement.jl
+++ b/src/Operators/spectralelement.jl
@@ -507,7 +507,12 @@ Divergence{()}(space) = Divergence{operator_axes(space)}()
 operator_return_eltype(op::Divergence{I}, ::Type{S}) where {I, S} =
     RecursiveApply.rmaptype(Geometry.divergence_result_type, S)
 
-function apply_operator(op::Divergence{(1,)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::Divergence{(1,)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -671,7 +676,12 @@ WeakDivergence{()}(space) = WeakDivergence{operator_axes(space)}()
 operator_return_eltype(::WeakDivergence{I}, ::Type{S}) where {I, S} =
     RecursiveApply.rmaptype(Geometry.divergence_result_type, S)
 
-function apply_operator(op::WeakDivergence{(1,)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::WeakDivergence{(1,)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -701,7 +711,12 @@ function apply_operator(op::WeakDivergence{(1,)}, space, slabidx, arg)
     return Field(SArray(out), space)
 end
 
-function apply_operator(op::WeakDivergence{(1, 2)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::WeakDivergence{(1, 2)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -812,7 +827,12 @@ Gradient{()}(space) = Gradient{operator_axes(space)}()
 operator_return_eltype(::Gradient{I}, ::Type{S}) where {I, S} =
     RecursiveApply.rmaptype(T -> Geometry.gradient_result_type(Val(I), T), S)
 
-function apply_operator(op::Gradient{(1,)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::Gradient{(1,)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -935,7 +955,12 @@ WeakGradient{()}(space) = WeakGradient{operator_axes(space)}()
 operator_return_eltype(::WeakGradient{I}, ::Type{S}) where {I, S} =
     RecursiveApply.rmaptype(T -> Geometry.gradient_result_type(Val(I), T), S)
 
-function apply_operator(op::WeakGradient{(1,)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::WeakGradient{(1,)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -963,7 +988,12 @@ function apply_operator(op::WeakGradient{(1,)}, space, slabidx, arg)
     return Field(SArray(out), space)
 end
 
-function apply_operator(op::WeakGradient{(1, 2)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::WeakGradient{(1, 2)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -1085,7 +1115,12 @@ Curl{()}(space) = Curl{operator_axes(space)}()
 operator_return_eltype(::Curl{I}, ::Type{S}) where {I, S} =
     RecursiveApply.rmaptype(T -> Geometry.curl_result_type(Val(I), T), S)
 
-function apply_operator(op::Curl{(1,)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::Curl{(1,)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -1116,7 +1151,12 @@ function apply_operator(op::Curl{(1,)}, space, slabidx, arg)
     return Field(SArray(out), space)
 end
 
-function apply_operator(op::Curl{(1, 2)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::Curl{(1, 2)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -1280,7 +1320,12 @@ WeakCurl{()}(space) = WeakCurl{operator_axes(space)}()
 operator_return_eltype(::WeakCurl{I}, ::Type{S}) where {I, S} =
     RecursiveApply.rmaptype(T -> Geometry.curl_result_type(Val(I), T), S)
 
-function apply_operator(op::WeakCurl{(1,)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::WeakCurl{(1,)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -1313,7 +1358,12 @@ function apply_operator(op::WeakCurl{(1,)}, space, slabidx, arg)
     return Field(SArray(out), space)
 end
 
-function apply_operator(op::WeakCurl{(1, 2)}, space, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::WeakCurl{(1, 2)},
+    space,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space)
     QS = Spaces.quadrature_style(space)
     Nq = Quadratures.degrees_of_freedom(QS)
@@ -1464,7 +1514,12 @@ struct Interpolate{I, S} <: TensorOperator
 end
 Interpolate(space) = Interpolate{operator_axes(space), typeof(space)}(space)
 
-function apply_operator(op::Interpolate{(1,)}, space_out, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::Interpolate{(1,)},
+    space_out,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space_out)
     space_in = axes(arg)
     QS_in = Spaces.quadrature_style(space_in)
@@ -1491,7 +1546,12 @@ function apply_operator(op::Interpolate{(1,)}, space_out, slabidx, arg)
     return Field(SArray(out), space_out)
 end
 
-function apply_operator(op::Interpolate{(1, 2)}, space_out, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::Interpolate{(1, 2)},
+    space_out,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space_out)
     space_in = axes(arg)
     QS_in = Spaces.quadrature_style(space_in)
@@ -1553,7 +1613,12 @@ struct Restrict{I, S} <: TensorOperator
 end
 Restrict(space) = Restrict{operator_axes(space), typeof(space)}(space)
 
-function apply_operator(op::Restrict{(1,)}, space_out, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::Restrict{(1,)},
+    space_out,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space_out)
     space_in = axes(arg)
     QS_in = Spaces.quadrature_style(space_in)
@@ -1584,7 +1649,12 @@ function apply_operator(op::Restrict{(1,)}, space_out, slabidx, arg)
     return Field(SArray(out), space_out)
 end
 
-function apply_operator(op::Restrict{(1, 2)}, space_out, slabidx, arg)
+Base.@propagate_inbounds function apply_operator(
+    op::Restrict{(1, 2)},
+    space_out,
+    slabidx,
+    arg,
+)
     FT = Spaces.undertype(space_out)
     space_in = axes(arg)
     QS_in = Spaces.quadrature_style(space_in)
diff --git a/src/Spaces/dss_transform.jl b/src/Spaces/dss_transform.jl
index 8bf6c0e277..3e02dfecff 100644
--- a/src/Spaces/dss_transform.jl
+++ b/src/Spaces/dss_transform.jl
@@ -12,17 +12,37 @@ Transformations only apply to vector quantities.
 
 See [`Spaces.weighted_dss!`](@ref).
 """
-dss_transform(arg, local_geometry, weight, i, j) =
+Base.@propagate_inbounds dss_transform(arg, local_geometry, weight, i, j) =
     dss_transform(arg[i, j], local_geometry[i, j], weight[i, j])
-dss_transform(arg, local_geometry, weight::Nothing, i, j) =
-    dss_transform(arg[i, j], local_geometry[i, j], 1)
-dss_transform(arg, local_geometry::Nothing, weight::Nothing, i, j) = arg[i, j]
+Base.@propagate_inbounds dss_transform(
+    arg,
+    local_geometry,
+    weight::Nothing,
+    i,
+    j,
+) = dss_transform(arg[i, j], local_geometry[i, j], 1)
+Base.@propagate_inbounds dss_transform(
+    arg,
+    local_geometry::Nothing,
+    weight::Nothing,
+    i,
+    j,
+) = arg[i, j]
 
-dss_transform(arg, local_geometry, weight, i) =
+Base.@propagate_inbounds dss_transform(arg, local_geometry, weight, i) =
     dss_transform(arg[i], local_geometry[i], weight[i])
-dss_transform(arg, local_geometry, weight::Nothing, i) =
-    dss_transform(arg[i], local_geometry[i], 1)
-dss_transform(arg, local_geometry::Nothing, weight::Nothing, i) = arg[i]
+Base.@propagate_inbounds dss_transform(
+    arg,
+    local_geometry,
+    weight::Nothing,
+    i,
+) = dss_transform(arg[i], local_geometry[i], 1)
+Base.@propagate_inbounds dss_transform(
+    arg,
+    local_geometry::Nothing,
+    weight::Nothing,
+    i,
+) = arg[i]
 
 @inline function dss_transform(
     arg::Tuple{},
@@ -128,12 +148,21 @@ Transform `targ[I...]` back to a value of type `T` after performing direct stiff
 
 See [`Spaces.weighted_dss!`](@ref).
 """
-dss_untransform(::Type{T}, targ, local_geometry, i, j) where {T} =
-    dss_untransform(T, targ, local_geometry[i, j])
+Base.@propagate_inbounds dss_untransform(
+    ::Type{T},
+    targ,
+    local_geometry,
+    i,
+    j,
+) where {T} = dss_untransform(T, targ, local_geometry[i, j])
 dss_untransform(::Type{T}, targ, local_geometry::Nothing, i, j) where {T} =
     dss_untransform(T, targ, local_geometry)
-dss_untransform(::Type{T}, targ, local_geometry, i) where {T} =
-    dss_untransform(T, targ, local_geometry[i])
+Base.@propagate_inbounds dss_untransform(
+    ::Type{T},
+    targ,
+    local_geometry,
+    i,
+) where {T} = dss_untransform(T, targ, local_geometry[i])
 dss_untransform(::Type{T}, targ, local_geometry::Nothing, i) where {T} =
     dss_untransform(T, targ, local_geometry)
 @inline function dss_untransform(