From 7407988b2c208ffac5993cd07d4d12dbb9f21e00 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Fri, 16 Jun 2023 16:54:00 -0700 Subject: [PATCH 1/2] Cache intermediate variables in topo example --- .../hybrid/plane/topo_dc_2dinvariant_rhoe.jl | 197 ++++++++++++++---- 1 file changed, 155 insertions(+), 42 deletions(-) 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(), From 9a3ca976d7c15d6443961b92235500ce8ae57720 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Fri, 16 Jun 2023 16:57:18 -0700 Subject: [PATCH 2/2] Apply propagate inbounds to some SEM operators --- src/Operators/spectralelement.jl | 98 +++++++++++++++++++++++++++----- src/Spaces/dss_transform.jl | 53 +++++++++++++---- 2 files changed, 125 insertions(+), 26 deletions(-) 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(