CliMA · charleskawczynski · Mar 8, 2025
diff --git a/examples/hybrid/implicit_equation_jacobian.jl b/examples/hybrid/implicit_equation_jacobian.jl
@@ -2,6 +2,7 @@ import LinearAlgebra: ldiv!
 using ClimaCore: Spaces, Fields, Operators
 using ClimaCore.Utilities: half
 using ClimaCore.MatrixFields
+import ClimaCore.MatrixFields: ⋆
 using ClimaCore.MatrixFields: @name
 
 """

diff --git a/examples/hybrid/staggered_nonhydrostatic_model.jl b/examples/hybrid/staggered_nonhydrostatic_model.jl
@@ -331,11 +331,11 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
     #     norm_sqr(C123(ᶜuₕ) + C123(ᶜinterp(ᶠw))) / 2 =
     #     ACT12(ᶜuₕ) * ᶜuₕ / 2 + ACT3(ᶜinterp(ᶠw)) * ᶜinterp(ᶠw) / 2
     # ∂(ᶜK)/∂(ᶠw) = ACT3(ᶜinterp(ᶠw)) * ᶜinterp_matrix()
-    @. ∂ᶜK∂ᶠw = DiagonalMatrixRow(adjoint(CT3(ᶜinterp(ᶠw)))) ⋅ ᶜinterp_matrix()
+    @. ∂ᶜK∂ᶠw = DiagonalMatrixRow(adjoint(CT3(ᶜinterp(ᶠw)))) ⋆ ᶜinterp_matrix()
 
     # ᶜρₜ = -ᶜdivᵥ(ᶠinterp(ᶜρ) * ᶠw)
     # ∂(ᶜρₜ)/∂(ᶠw) = -ᶜdivᵥ_matrix() * ᶠinterp(ᶜρ) * ᶠg³³
-    @. ∂ᶜρₜ∂ᶠ𝕄 = -(ᶜdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ) * g³³(ᶠgⁱʲ))
+    @. ∂ᶜρₜ∂ᶠ𝕄 = -(ᶜdivᵥ_matrix()) ⋆ DiagonalMatrixRow(ᶠinterp(ᶜρ) * g³³(ᶠgⁱʲ))
 
     if :ρθ in propertynames(Y.c)
         ᶜρθ = Y.c.ρθ
@@ -349,14 +349,14 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
             # ᶜρθₜ = -ᶜdivᵥ(ᶠinterp(ᶜρθ) * ᶠw)
             # ∂(ᶜρθₜ)/∂(ᶠw) = -ᶜdivᵥ_matrix() * ᶠinterp(ᶜρθ) * ᶠg³³
             @. ∂ᶜ𝔼ₜ∂ᶠ𝕄 =
-                -(ᶜdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρθ) * g³³(ᶠgⁱʲ))
+                -(ᶜdivᵥ_matrix()) ⋆ DiagonalMatrixRow(ᶠinterp(ᶜρθ) * g³³(ᶠgⁱʲ))
         else
             # ᶜρθₜ = -ᶜdivᵥ(ᶠinterp(ᶜρ) * ᶠupwind_product(ᶠw, ᶜρθ / ᶜρ))
             # ∂(ᶜρθₜ)/∂(ᶠw) =
             #     -ᶜdivᵥ_matrix() * ᶠinterp(ᶜρ) *
             #     ∂(ᶠupwind_product(ᶠw, ᶜρθ / ᶜρ))/∂(ᶠw)
             @. ∂ᶜ𝔼ₜ∂ᶠ𝕄 =
-                -(ᶜdivᵥ_matrix()) ⋅ DiagonalMatrixRow(
+                -(ᶜdivᵥ_matrix()) ⋆ DiagonalMatrixRow(
                     ᶠinterp(ᶜρ) *
                     vec_data(ᶠno_flux(ᶠupwind_product(ᶠw + εw, ᶜρθ / ᶜρ))) /
                     vec_data(CT3(ᶠw + εw)) * g³³(ᶠgⁱʲ),
@@ -381,10 +381,10 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
                 # ∂(ᶜp)/∂(ᶠw) = ∂(ᶜp)/∂(ᶜK) * ∂(ᶜK)/∂(ᶠw)
                 # ∂(ᶜp)/∂(ᶜK) = -ᶜρ * R_d / cv_d
                 @. ∂ᶜ𝔼ₜ∂ᶠ𝕄 =
-                    -(ᶜdivᵥ_matrix()) ⋅ (
+                    -(ᶜdivᵥ_matrix()) ⋆ (
                         DiagonalMatrixRow(ᶠinterp(ᶜρe + ᶜp) * g³³(ᶠgⁱʲ)) +
-                        DiagonalMatrixRow(CT3(ᶠw)) ⋅ ᶠinterp_matrix() ⋅
-                        DiagonalMatrixRow(-(ᶜρ * R_d / cv_d)) ⋅ ∂ᶜK∂ᶠw
+                        DiagonalMatrixRow(CT3(ᶠw)) ⋆ ᶠinterp_matrix() ⋆
+                        DiagonalMatrixRow(-(ᶜρ * R_d / cv_d)) ⋆ ∂ᶜK∂ᶠw
                     )
             else
                 # ᶜρeₜ =
@@ -397,15 +397,15 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
                 # ∂(ᶜp)/∂(ᶠw) = ∂(ᶜp)/∂(ᶜK) * ∂(ᶜK)/∂(ᶠw)
                 # ∂(ᶜp)/∂(ᶜK) = -ᶜρ * R_d / cv_d
                 @. ∂ᶜ𝔼ₜ∂ᶠ𝕄 =
-                    -(ᶜdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ)) ⋅ (
+                    -(ᶜdivᵥ_matrix()) ⋆ DiagonalMatrixRow(ᶠinterp(ᶜρ)) ⋆ (
                         DiagonalMatrixRow(
                             vec_data(
                                 ᶠno_flux(
                                     ᶠupwind_product(ᶠw + εw, (ᶜρe + ᶜp) / ᶜρ),
                                 ),
                             ) / vec_data(CT3(ᶠw + εw)) * g³³(ᶠgⁱʲ),
                         ) +
-                        ᶠno_flux_row(ᶠupwind_product_matrix(ᶠw)) ⋅
+                        ᶠno_flux_row(ᶠupwind_product_matrix(ᶠw)) ⋆
                         (-R_d / cv_d * ∂ᶜK∂ᶠw)
                     )
             end
@@ -414,11 +414,11 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
             # ∂ᶜ𝔼ₜ∂ᶠ𝕄 has 3 diagonals instead of 5
             if isnothing(ᶠupwind_product)
                 @. ∂ᶜ𝔼ₜ∂ᶠ𝕄 =
-                    -(ᶜdivᵥ_matrix()) ⋅
+                    -(ᶜdivᵥ_matrix()) ⋆
                     DiagonalMatrixRow(ᶠinterp(ᶜρe + ᶜp) * g³³(ᶠgⁱʲ))
             else
                 @. ∂ᶜ𝔼ₜ∂ᶠ𝕄 =
-                    -(ᶜdivᵥ_matrix()) ⋅ DiagonalMatrixRow(
+                    -(ᶜdivᵥ_matrix()) ⋆ DiagonalMatrixRow(
                         ᶠinterp(ᶜρ) * vec_data(
                             ᶠno_flux(ᶠupwind_product(ᶠw + εw, (ᶜρe + ᶜp) / ᶜρ)),
                         ) / vec_data(CT3(ᶠw + εw)) * g³³(ᶠgⁱʲ),
@@ -443,9 +443,9 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
             #     -ᶜdivᵥ_matrix() * ᶠinterp(ᶜρe_int + ᶜp) * ᶠg³³ +
             #     ᶜinterp_matrix() * adjoint(ᶠgradᵥ(ᶜp)) * ᶠg³³
             @. ∂ᶜ𝔼ₜ∂ᶠ𝕄 =
-                -(ᶜdivᵥ_matrix()) ⋅
+                -(ᶜdivᵥ_matrix()) ⋆
                 DiagonalMatrixRow(ᶠinterp(ᶜρe_int + ᶜp) * g³³(ᶠgⁱʲ)) +
-                ᶜinterp_matrix() ⋅
+                ᶜinterp_matrix() ⋆
                 DiagonalMatrixRow(adjoint(ᶠgradᵥ(ᶜp)) * g³³(ᶠgⁱʲ))
         else
             # ᶜρe_intₜ =
@@ -456,12 +456,12 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
             #     ∂(ᶠupwind_product(ᶠw, (ᶜρe_int + ᶜp) / ᶜρ))/∂(ᶠw) +
             #     ᶜinterp_matrix() * adjoint(ᶠgradᵥ(ᶜp)) * ᶠg³³
             @. ∂ᶜ𝔼ₜ∂ᶠ𝕄 =
-                -(ᶜdivᵥ_matrix()) ⋅ DiagonalMatrixRow(
+                -(ᶜdivᵥ_matrix()) ⋆ DiagonalMatrixRow(
                     ᶠinterp(ᶜρ) * vec_data(
                         ᶠno_flux(ᶠupwind_product(ᶠw + εw, (ᶜρe_int + ᶜp) / ᶜρ)),
                     ) / vec_data(CT3(ᶠw + εw)) * g³³(ᶠgⁱʲ),
                 ) +
-                ᶜinterp_matrix() ⋅
+                ᶜinterp_matrix() ⋆
                 DiagonalMatrixRow(adjoint(ᶠgradᵥ(ᶜp)) * g³³(ᶠgⁱʲ))
         end
     end
@@ -480,7 +480,7 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
         # ∂(ᶠgradᵥ(ᶜp))/∂(ᶜρθ) =
         #     ᶠgradᵥ_matrix() * γ * R_d * (ᶜρθ * R_d / p_0)^(γ - 1)
         @. ∂ᶠ𝕄ₜ∂ᶜ𝔼 =
-            -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix() ⋅
+            -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋆ ᶠgradᵥ_matrix() ⋆
             DiagonalMatrixRow(γ * R_d * (ᶜρθ * R_d / p_0)^(γ - 1))
 
         if flags.∂ᶠ𝕄ₜ∂ᶜρ_mode == :exact
@@ -489,20 +489,20 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
             # ∂(ᶠwₜ)/∂(ᶠinterp(ᶜρ)) = ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ)^2
             # ∂(ᶠinterp(ᶜρ))/∂(ᶜρ) = ᶠinterp_matrix()
             @. ∂ᶠ𝕄ₜ∂ᶜρ =
-                DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ)^2) ⋅ ᶠinterp_matrix()
+                DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ)^2) ⋆ ᶠinterp_matrix()
         elseif flags.∂ᶠ𝕄ₜ∂ᶜρ_mode == :hydrostatic_balance
             # same as above, but we assume that ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ) =
             # -ᶠgradᵥ(ᶜΦ)
             @. ∂ᶠ𝕄ₜ∂ᶜρ =
-                -DiagonalMatrixRow(ᶠgradᵥ(ᶜΦ) / ᶠinterp(ᶜρ)) ⋅ ᶠinterp_matrix()
+                -DiagonalMatrixRow(ᶠgradᵥ(ᶜΦ) / ᶠinterp(ᶜρ)) ⋆ ᶠinterp_matrix()
         end
     elseif :ρe in propertynames(Y.c)
         # ᶠwₜ = -ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ) - ᶠgradᵥ(ᶜK + ᶜΦ)
         # ∂(ᶠwₜ)/∂(ᶜρe) = ∂(ᶠwₜ)/∂(ᶠgradᵥ(ᶜp)) * ∂(ᶠgradᵥ(ᶜp))/∂(ᶜρe)
         # ∂(ᶠwₜ)/∂(ᶠgradᵥ(ᶜp)) = -1 / ᶠinterp(ᶜρ)
         # ∂(ᶠgradᵥ(ᶜp))/∂(ᶜρe) = ᶠgradᵥ_matrix() * R_d / cv_d
         @. ∂ᶠ𝕄ₜ∂ᶜ𝔼 =
-            -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋅ (ᶠgradᵥ_matrix() * R_d / cv_d)
+            -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋆ (ᶠgradᵥ_matrix() * R_d / cv_d)
 
         if flags.∂ᶠ𝕄ₜ∂ᶜρ_mode == :exact
             # ᶠwₜ = -ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ) - ᶠgradᵥ(ᶜK + ᶜΦ)
@@ -515,24 +515,24 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
             # ∂(ᶠwₜ)/∂(ᶠinterp(ᶜρ)) = ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ)^2
             # ∂(ᶠinterp(ᶜρ))/∂(ᶜρ) = ᶠinterp_matrix()
             @. ∂ᶠ𝕄ₜ∂ᶜρ =
-                -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix() ⋅
+                -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋆ ᶠgradᵥ_matrix() ⋆
                 DiagonalMatrixRow(R_d * (-(ᶜK + ᶜΦ) / cv_d + T_tri)) +
-                DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ)^2) ⋅ ᶠinterp_matrix()
+                DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ)^2) ⋆ ᶠinterp_matrix()
         elseif flags.∂ᶠ𝕄ₜ∂ᶜρ_mode == :hydrostatic_balance
             # same as above, but we assume that ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ) =
             # -ᶠgradᵥ(ᶜΦ) and that ᶜK is negligible compared ot ᶜΦ
             @. ∂ᶠ𝕄ₜ∂ᶜρ =
-                -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix() ⋅
+                -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋆ ᶠgradᵥ_matrix() ⋆
                 DiagonalMatrixRow(R_d * (-(ᶜΦ) / cv_d + T_tri)) -
-                DiagonalMatrixRow(ᶠgradᵥ(ᶜΦ) / ᶠinterp(ᶜρ)) ⋅ ᶠinterp_matrix()
+                DiagonalMatrixRow(ᶠgradᵥ(ᶜΦ) / ᶠinterp(ᶜρ)) ⋆ ᶠinterp_matrix()
         end
     elseif :ρe_int in propertynames(Y.c)
         # ᶠwₜ = -ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ) - ᶠgradᵥ(ᶜK + ᶜΦ)
         # ∂(ᶠwₜ)/∂(ᶜρe_int) = ∂(ᶠwₜ)/∂(ᶠgradᵥ(ᶜp)) * ∂(ᶠgradᵥ(ᶜp))/∂(ᶜρe_int)
         # ∂(ᶠwₜ)/∂(ᶠgradᵥ(ᶜp)) = -1 / ᶠinterp(ᶜρ)
         # ∂(ᶠgradᵥ(ᶜp))/∂(ᶜρe_int) = ᶠgradᵥ_matrix() * R_d / cv_d
         @. ∂ᶠ𝕄ₜ∂ᶜ𝔼 =
-            DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅ (ᶠgradᵥ_matrix() * R_d / cv_d)
+            DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋆ (ᶠgradᵥ_matrix() * R_d / cv_d)
 
         if flags.∂ᶠ𝕄ₜ∂ᶜρ_mode == :exact
             # ᶠwₜ = -ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ) - ᶠgradᵥ(ᶜK + ᶜΦ)
@@ -544,16 +544,16 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
             # ∂(ᶠwₜ)/∂(ᶠinterp(ᶜρ)) = ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ)^2
             # ∂(ᶠinterp(ᶜρ))/∂(ᶜρ) = ᶠinterp_matrix()
             @. ∂ᶠ𝕄ₜ∂ᶜρ =
-                -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋅
+                -DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋆
                 (ᶠgradᵥ_matrix() * R_d * T_tri) +
-                DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ)^2) ⋅ ᶠinterp_matrix()
+                DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ)^2) ⋆ ᶠinterp_matrix()
         elseif flags.∂ᶠ𝕄ₜ∂ᶜρ_mode == :hydrostatic_balance
             # same as above, but we assume that ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ) =
             # -ᶠgradᵥ(ᶜΦ)
             @. ∂ᶠ𝕄ₜ∂ᶜρ =
-                DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅
+                DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋆
                 (ᶠgradᵥ_matrix() * R_d * T_tri) -
-                DiagonalMatrixRow(ᶠgradᵥ(ᶜΦ) / ᶠinterp(ᶜρ)) ⋅ ᶠinterp_matrix()
+                DiagonalMatrixRow(ᶠgradᵥ(ᶜΦ) / ᶠinterp(ᶜρ)) ⋆ ᶠinterp_matrix()
         end
     end
 
@@ -571,13 +571,13 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
     # ∂(ᶠwₜ)/∂(ᶠgradᵥ(ᶜK + ᶜΦ)) = -1
     # ∂(ᶠgradᵥ(ᶜK + ᶜΦ))/∂(ᶜK) = ᶠgradᵥ_matrix()
     if :ρθ in propertynames(Y.c) || :ρe_int in propertynames(Y.c)
-        @. ∂ᶠ𝕄ₜ∂ᶠ𝕄 = -(ᶠgradᵥ_matrix()) ⋅ ∂ᶜK∂ᶠw
+        @. ∂ᶠ𝕄ₜ∂ᶠ𝕄 = -(ᶠgradᵥ_matrix()) ⋆ ∂ᶜK∂ᶠw
     elseif :ρe in propertynames(Y.c)
         @. ∂ᶠ𝕄ₜ∂ᶠ𝕄 =
             -(
-                DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix() ⋅
+                DiagonalMatrixRow(1 / ᶠinterp(ᶜρ)) ⋆ ᶠgradᵥ_matrix() ⋆
                 DiagonalMatrixRow(-(ᶜρ * R_d / cv_d)) + ᶠgradᵥ_matrix()
-            ) ⋅ ∂ᶜK∂ᶠw
+            ) ⋆ ∂ᶜK∂ᶠw
     end
 
     I = one(∂R∂Y)

diff --git a/src/MatrixFields/MatrixFields.jl b/src/MatrixFields/MatrixFields.jl
@@ -78,7 +78,7 @@ export DiagonalMatrixRow,
     QuaddiagonalMatrixRow,
     PentadiagonalMatrixRow
 export FieldVectorKeys, FieldMatrixKeys, FieldVectorView, FieldMatrix
-export FieldMatrixWithSolver, ⋅
+export FieldMatrixWithSolver
 
 include("band_matrix_row.jl")
 

diff --git a/src/MatrixFields/matrix_multiplication.jl b/src/MatrixFields/matrix_multiplication.jl
@@ -200,6 +200,7 @@ This means that we can express the bounds on the interior values of ``i`` as
 """
 struct MultiplyColumnwiseBandMatrixField <: Operators.FiniteDifferenceOperator end
 const ⋅ = MultiplyColumnwiseBandMatrixField()
+const ⋆ = MultiplyColumnwiseBandMatrixField()
 
 Operators.strip_space(op::MultiplyColumnwiseBandMatrixField, _) = op
 

diff --git a/test/Fields/benchmark_fieldvectors.jl b/test/Fields/benchmark_fieldvectors.jl
@@ -16,17 +16,10 @@ if ClimaComms.device() isa ClimaComms.CUDADevice
 else
     device_name = "CPU"
 end
-if !(@isdefined(TU))
-    include(
-        joinpath(
-            pkgdir(ClimaCore),
-            "test",
-            "TestUtilities",
-            "TestUtilities.jl",
-        ),
-    )
-    import .TestUtilities as TU
-end
+@isdefined(TU) || include(
+    joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
+);
+import .TestUtilities as TU
 
 include(joinpath(pkgdir(ClimaCore), "benchmarks/scripts/benchmark_utils.jl"))
 

diff --git a/test/Fields/convergence_field_integrals.jl b/test/Fields/convergence_field_integrals.jl
@@ -28,9 +28,9 @@ using LinearAlgebra: norm
 using Statistics: mean
 using ForwardDiff
 
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
+);
 import .TestUtilities as TU
 
 """

diff --git a/test/Fields/field_opt.jl b/test/Fields/field_opt.jl
@@ -25,9 +25,9 @@ using LinearAlgebra: norm
 using Statistics: mean
 using ForwardDiff
 
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
+);
 import .TestUtilities as TU
 
 # https://github.com/CliMA/ClimaCore.jl/issues/946

diff --git a/test/Fields/unit_field.jl b/test/Fields/unit_field.jl
@@ -29,9 +29,9 @@ using LinearAlgebra: norm
 using Statistics: mean
 using ForwardDiff
 
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
+);
 import .TestUtilities as TU
 
 function spectral_space_2D(; n1 = 1, n2 = 1, Nij = 4)

diff --git a/test/Fields/utils_field_multi_broadcast_fusion.jl b/test/Fields/utils_field_multi_broadcast_fusion.jl
@@ -31,17 +31,10 @@ using LinearAlgebra: norm
 using Statistics: mean
 using ForwardDiff
 
-if !(@isdefined(TU))
-    include(
-        joinpath(
-            pkgdir(ClimaCore),
-            "test",
-            "TestUtilities",
-            "TestUtilities.jl",
-        ),
-    )
-    import .TestUtilities as TU
-end
+@isdefined(TU) || include(
+    joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
+);
+import .TestUtilities as TU
 
 @show ClimaComms.device()
 

diff --git a/test/InputOutput/unit_read_type.jl b/test/InputOutput/unit_read_type.jl
@@ -3,10 +3,10 @@ import ClimaCore
 import ClimaCore: Fields, InputOutput
 using ClimaComms
 ClimaComms.@import_required_backends
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
-import .TestUtilities as TU
+);
+import .TestUtilities as TU;
 
 compare_read_type(x) = InputOutput.read_type(string(eltype(x))) == eltype(x)
 @testset "Read field element types" begin

diff --git a/test/MatrixFields/field_names.jl b/test/MatrixFields/field_names.jl
@@ -5,18 +5,19 @@ import ClimaCore.MatrixFields: @name, is_subset_that_covers_set
 
 include("matrix_field_test_utils.jl")
 
-struct Foo{T}
+struct FieldNameFoo{T}
     _value::T
 end
-Base.propertynames(::Foo) = (:value,)
-Base.getproperty(foo::Foo, s::Symbol) =
+Base.propertynames(::FieldNameFoo) = (:value,)
+Base.getproperty(foo::FieldNameFoo, s::Symbol) =
     s == :value ? getfield(foo, :_value) : error("Invalid property name")
-Base.convert(::Type{Foo{T}}, foo::Foo) where {T} = Foo{T}(foo.value)
-Base.zero(::Type{Foo{T}}) where {T} = Foo(zero(T))
+Base.convert(::Type{FieldNameFoo{T}}, foo::FieldNameFoo) where {T} = FieldNameFoo{T}(foo.value)
+Base.zero(::Type{FieldNameFoo{T}}) where {T} = FieldNameFoo(zero(T))
 
-const x = (; foo = Foo(0), a = (; b = 1, c = ((; d = 2), (;), (3, ()))))
+get_x() = (; foo = FieldNameFoo(0), a = (; b = 1, c = ((; d = 2), (;), (3, ()))))
 
 @testset "FieldName Unit Tests" begin
+    x = get_x()
     @test_all @name() == MatrixFields.FieldName()
     @test_all @name(a.c.:(1).d) == MatrixFields.FieldName(:a, :c, 1, :d)
     @test_all @name(a.c.:(3).:(1)) == MatrixFields.FieldName(:a, :c, 3, 1)
@@ -81,6 +82,7 @@ const x = (; foo = Foo(0), a = (; b = 1, c = ((; d = 2), (;), (3, ()))))
 end
 
 @testset "FieldNameTree Unit Tests" begin
+    x = get_x()
     name_tree = MatrixFields.FieldNameTree(x)
 
     @test_all MatrixFields.FieldNameTree(x) == name_tree
@@ -110,6 +112,7 @@ end
 end
 
 @testset "FieldNameSet Unit Tests" begin
+    x = get_x()
     name_tree = MatrixFields.FieldNameTree(x)
 
     vector_keys(names...) = MatrixFields.FieldVectorKeys(names, name_tree)

diff --git a/test/MatrixFields/flat_spaces.jl b/test/MatrixFields/flat_spaces.jl
@@ -1,11 +1,12 @@
 import ClimaCore
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
+);
 import .TestUtilities as TU
 
 include("matrix_field_test_utils.jl")
-import ClimaCore.MatrixFields: @name, ⋅
+import ClimaCore.MatrixFields: @name
+import ClimaCore.MatrixFields: ⋆
 
 @testset "Matrix Fields with Spectral Element and Point Spaces" begin
     get_j_field(space, FT) = fill(MatrixFields.DiagonalMatrixRow(FT(1)), space)