3d curl in vector hyperdiffusion

CliMA · Aug 30, 2023 · 4cc117f · 4cc117f
1 parent 60a9f63
commit 4cc117f
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 36 deletions.
diff --git a/src/cache/temporary_quantities.jl b/src/cache/temporary_quantities.jl
@@ -49,7 +49,6 @@ function temporary_quantities(atmos, center_space, face_space)
         ᶠtemp_CT3 = Fields.Field(CT3{FT}, face_space), # ᶠuₕ³
         ᶠtemp_CT12 = Fields.Field(CT12{FT}, face_space), # ᶠω¹²
         ᶠtemp_CT12ʲs = Fields.Field(NTuple{n, CT12{FT}}, face_space), # ᶠω¹²ʲs
-        ᶜtemp_C123 = Fields.Field(C123{FT}, center_space), # χ₁₂₃
         ᶠtemp_C123 = Fields.Field(C123{FT}, face_space), # χ₁₂₃
         # TODO: Remove this hack
         sfc_temp_C3 = Fields.Field(C3{FT}, Spaces.level(face_space, half)), # ρ_flux_χ

diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl
@@ -80,13 +80,8 @@ function hyperdiffusion_tendency!(Yₜ, Y, p, t)
     end
 
     # Grid scale hyperdiffusion
-    @. ᶜ∇²u = C123(wgradₕ(divₕ(p.ᶜu)))
-    # Without the C123(), the right-hand side would be a C1 or C2 in 2D space.
-    if point_type <: Geometry.Abstract3DPoint
-        @. p.ᶜtemp_C123 = C123(curlₕ(C12(p.ᶜu))) + C123(curlₕ(C3(p.ᶜu)))
-        @. ᶜ∇²u -=
-            C123(wcurlₕ(C12(p.ᶜtemp_C123))) + C123(wcurlₕ(C3(p.ᶜtemp_C123)))
-    end
+    @. ᶜ∇²u = C123(wgradₕ(divₕ(u))) - C123(wcurlₕ(C123(curlₕ(u))))
+
     if :θ in propertynames(ᶜspecific)
         @. ᶜ∇²specific_energy = wdivₕ(gradₕ(ᶜspecific.θ))
     elseif :e_tot in propertynames(ᶜspecific)
@@ -102,16 +97,9 @@ function hyperdiffusion_tendency!(Yₜ, Y, p, t)
     # Sub-grid scale hyperdiffusion
     if turbconv_model isa EDMFX
         for j in 1:n
-            @. ᶜ∇²uʲs.:($$j) = C123(wgradₕ(divₕ(p.ᶜuʲs.:($$j))))
-            # Without the C123(), the right-hand side would be a C1 or C2 in 2D space.
-            if point_type <: Geometry.Abstract3DPoint
-                @. p.ᶜtemp_C123 =
-                    C123(curlₕ(C12(p.ᶜuʲs.:($$j)))) +
-                    C123(curlₕ(C3(p.ᶜuʲs.:($$j))))
-                @. ᶜ∇²uʲs.:($$j) -=
-                    C123(wcurlₕ(C12(p.ᶜtemp_C123))) +
-                    C123(wcurlₕ(C3(p.ᶜtemp_C123)))
-            end
+            @. ᶜ∇²uʲs.:($$j) =
+                C123(wgradₕ(divₕ(p.ᶜuʲs.:($$j)))) -
+                C123(wcurlₕ(C123(curlₕ(p.ᶜuʲs.:($$j)))))
 
             if :θ in propertynames(ᶜspecificʲs.:($j))
                 @. ᶜ∇²specific_energyʲs.:($$j) =
@@ -164,18 +152,12 @@ function hyperdiffusion_tendency!(Yₜ, Y, p, t)
         end
     end
 
-    # Grid scale hyperdiffusion continued
-    @. Yₜ.c.uₕ -= κ₄ * divergence_damping_factor * C12(wgradₕ(divₕ(ᶜ∇²u)))
-    # Without the C123(), the right-hand side would be a C1 or C2 in 2D space.
-    if point_type <: Geometry.Abstract3DPoint
-        @. p.ᶜtemp_C123 = C123(curlₕ(C12(ᶜ∇²u))) + C123(curlₕ(C3(ᶜ∇²u)))
-        @. Yₜ.c.uₕ +=
-            κ₄ *
-            (C12(wcurlₕ(C12(p.ᶜtemp_C123))) + C12(wcurlₕ(C3(p.ᶜtemp_C123))))
-        # Reuse the buffer variable here so we don't need an extra one
-        @. ᶜ∇²uᵥ = C3(wcurlₕ(C12(p.ᶜtemp_C123))) + C3(wcurlₕ(C3(p.ᶜtemp_C123)))
-        @. Yₜ.f.u₃ += κ₄ * ᶠwinterp(ᶜJ * Y.c.ρ, ᶜ∇²uᵥ)
-    end
+    # re-use to store the curl-curl part
+    @. ᶜ∇²u = C123(wcurlₕ(C123(curlₕ(ᶜ∇²u))))
+    @. Yₜ.c.uₕ -=
+        κ₄ * (divergence_damping_factor * C12(wgradₕ(divₕ(ᶜ∇²u)) - C12(ᶜ∇²u)))
+    @. Yₜ.f.u₃ += κ₄ * ᶠwinterp(ᶜJ * Y.c.ρ, C3(ᶜ∇²u))
+
     ᶜρ_energyₜ = :θ in propertynames(ᶜspecific) ? Yₜ.c.ρθ : Yₜ.c.ρe_tot
     @. ᶜρ_energyₜ -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²specific_energy))
 
@@ -186,12 +168,8 @@ function hyperdiffusion_tendency!(Yₜ, Y, p, t)
     if turbconv_model isa EDMFX
         for j in 1:n
             if point_type <: Geometry.Abstract3DPoint
-                @. p.ᶜtemp_C123 =
-                    C123(curlₕ(C12(ᶜ∇²uʲs.:($$j)))) +
-                    C123(curlₕ(C3(ᶜ∇²uʲs.:($$j))))
-                # Reuse the buffer variable here so we don't need an extra one
-                @. ᶜ∇²uᵥʲs.:($$j) =
-                    C3(wcurlₕ(C12(p.ᶜtemp_C123))) + C3(wcurlₕ(C3(p.ᶜtemp_C123)))
+                # only need curl-curl part
+                @. ᶜ∇²uᵥʲs.:($$j) = C123(wcurlₕ(C123(curlₕ(ᶜ∇²uʲs.:($$j)))))
                 @. Yₜ.f.sgsʲs.:($$j).u₃ +=
                     κ₄ * ᶠwinterp(ᶜJ * Y.c.ρ, ᶜ∇²uᵥʲs.:($$j))
             end