diff --git a/examples/hybrid/implicit_equation_jacobian.jl b/examples/hybrid/implicit_equation_jacobian.jl
index 4dd1174d7b..9855dab745 100644
--- 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
index 21fc6ba9d5..00c00fb5b1 100644
--- 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,7 +397,7 @@ 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(
@@ -405,7 +405,7 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
                                 ),
                             ) / 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,12 +489,12 @@ 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 + ᶜΦ)
@@ -502,7 +502,7 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
         # ∂(ᶠ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,16 +515,16 @@ 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 + ᶜΦ)
@@ -532,7 +532,7 @@ function implicit_equation_jacobian!(j, Y, p, δtγ, t)
         # ∂(ᶠ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
index 3d4725ca48..7dd68fbfbd 100644
--- 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
index 5df6eb9a83..b2e81d1b30 100644
--- 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
index 2098f6ea77..0eff7b649c 100644
--- 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
index 17ed08fdf3..d81d81428d 100644
--- 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
index c6b6c4f686..910d4838ba 100644
--- 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
index fa0492f4a9..68c677d8b9 100644
--- 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
index 1557f94944..64aff058fa 100644
--- 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
index 94779eaf3f..7d3766ab45 100644
--- 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
index e6f94e1464..a57b49ab2c 100644
--- 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
index 0af2d92171..f3188d485b 100644
--- 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)
diff --git a/test/MatrixFields/gpu_compat_bidiag_matrix_row.jl b/test/MatrixFields/gpu_compat_bidiag_matrix_row.jl
index bbdfa5247a..98f683fcac 100644
--- a/test/MatrixFields/gpu_compat_bidiag_matrix_row.jl
+++ b/test/MatrixFields/gpu_compat_bidiag_matrix_row.jl
@@ -4,29 +4,22 @@ using Revise; include(joinpath("test", "MatrixFields", "gpu_compat_bidiag_matrix
 import ClimaCore
 import ClimaComms
 ClimaComms.@import_required_backends
-if !(@isdefined(TU))
-    include(
-        joinpath(
-            pkgdir(ClimaCore),
-            "test",
-            "TestUtilities",
-            "TestUtilities.jl",
-        ),
-    )
-end
+@isdefined(TU) || include(
+    joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
+);
 import .TestUtilities as TU
 
 import ClimaCore: Spaces, Geometry, Operators, Fields, MatrixFields
 using LinearAlgebra: Adjoint
 import StaticArrays: SArray
 import ClimaCore.Geometry: AxisTensor, CovariantAxis, ContravariantAxis
+import ClimaCore.MatrixFields: ⋆
 using ClimaCore.MatrixFields:
     BandMatrixRow,
     DiagonalMatrixRow,
     BidiagonalMatrixRow,
-    TridiagonalMatrixRow,
-    MultiplyColumnwiseBandMatrixField,
-    ⋅
+    TridiagonalMatrixRow
+
 const C3 = Geometry.Covariant3Vector
 const CT3 = Geometry.Contravariant3Vector
 GFT = Float64
@@ -89,14 +82,14 @@ function foo(c, f)
     dtγ = FT(1)
 
     @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
-        dtγ * ᶠtridiagonal_matrix_c3 ⋅ DiagonalMatrixRow(adjoint(CT3(ᶠu₃))) -
+        dtγ * ᶠtridiagonal_matrix_c3 ⋆ DiagonalMatrixRow(adjoint(CT3(ᶠu₃))) -
         (I_u₃,)
 
-    @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = dtγ * ᶠtridiagonal_matrix_c3 ⋅ adj_u₃ - (I_u₃,)
+    @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = dtγ * ᶠtridiagonal_matrix_c3 ⋆ adj_u₃ - (I_u₃,)
 
     # Fails on gpu
     @. ᶠtridiagonal_matrix_c3 =
-        -(ᶠgradᵥ_matrix()) ⋅ ifelse(
+        -(ᶠgradᵥ_matrix()) ⋆ ifelse(
             ᶜu₃ʲ.components.data.:1 > 0,
             convert(BidiagonalMatrixRow{FT}, ᶜleft_bias_matrix()),
             convert(BidiagonalMatrixRow{FT}, ᶜright_bias_matrix()),
@@ -107,7 +100,7 @@ function foo(c, f)
     @. bdmr_l = convert(BidiagonalMatrixRow{FT}, ᶜleft_bias_matrix())
     @. bdmr_r = convert(BidiagonalMatrixRow{FT}, ᶜright_bias_matrix())
     @. bdmr = ifelse(ᶜu₃ʲ.components.data.:1 > 0, bdmr_l, bdmr_r)
-    @. ᶠtridiagonal_matrix_c3 = -(ᶠgradᵥ_matrix()) ⋅ bdmr
+    @. ᶠtridiagonal_matrix_c3 = -(ᶠgradᵥ_matrix()) ⋆ bdmr
 
     return nothing
 end
diff --git a/test/MatrixFields/matrix_field_test_utils.jl b/test/MatrixFields/matrix_field_test_utils.jl
index 39a6f883e6..24e766c2e8 100644
--- a/test/MatrixFields/matrix_field_test_utils.jl
+++ b/test/MatrixFields/matrix_field_test_utils.jl
@@ -21,6 +21,7 @@ import ClimaCore:
     Operators,
     Quadratures
 using ClimaCore.MatrixFields
+import ClimaCore.MatrixFields: ⋆
 
 # Test that an expression is true and that it is also type-stable.
 macro test_all(expression)
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_1.jl b/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_1.jl
index 465f9ef2da..a525cdccdc 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_1.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_1.jl
@@ -3,6 +3,7 @@ julia --project
 using Revise; include(joinpath("test", "MatrixFields", "matrix_fields_broadcasting", "test_non_scalar_1.jl"))
 =#
 import ClimaCore
+import ClimaCore.MatrixFields: ⋆
 #! format: off
 if !(@isdefined(unit_test_field_broadcast))
     include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcasting","test_non_scalar_utils.jl"))
@@ -11,14 +12,14 @@ end
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "matrix of covectors times matrix of vectors" begin
 
-    bc = @lazy @. ᶜᶠmat_AC1 ⋅ ᶠᶜmat_C12
+    bc = @lazy @. ᶜᶠmat_AC1 ⋆ ᶠᶜmat_C12
     result = materialize(bc)
 
     ref_set_result! =
         result -> (@. result =
-            ᶜᶠmat ⋅ (
-                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:1) ⋅ ᶠᶜmat2 +
-                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:2) ⋅ ᶠᶜmat3
+            ᶜᶠmat ⋆ (
+                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:1) ⋆ ᶠᶜmat2 +
+                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:2) ⋆ ᶠᶜmat3
             ))
 
     unit_test_field_broadcast(
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_2.jl b/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_2.jl
index ead7c70c07..4e9fec84fe 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_2.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_2.jl
@@ -13,17 +13,17 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
                  of numbers times matrix of covectors times matrix of \
                  vectors" begin
 
-    bc = @lazy @. ᶜᶠmat_AC1 ⋅ ᶠᶜmat_C12 ⋅ ᶜᶠmat ⋅ ᶠᶜmat_AC1 ⋅ ᶜᶠmat_C12
+    bc = @lazy @. ᶜᶠmat_AC1 ⋆ ᶠᶜmat_C12 ⋆ ᶜᶠmat ⋆ ᶠᶜmat_AC1 ⋆ ᶜᶠmat_C12
     result = materialize(bc)
 
     ref_set_result! =
         result -> (@. result =
-            ᶜᶠmat ⋅ (
-                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:1) ⋅ ᶠᶜmat2 +
-                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:2) ⋅ ᶠᶜmat3
-            ) ⋅ ᶜᶠmat ⋅ ᶠᶜmat ⋅ (
-                DiagonalMatrixRow(ᶜlg.gⁱʲ.components.data.:1) ⋅ ᶜᶠmat2 +
-                DiagonalMatrixRow(ᶜlg.gⁱʲ.components.data.:2) ⋅ ᶜᶠmat3
+            ᶜᶠmat ⋆ (
+                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:1) ⋆ ᶠᶜmat2 +
+                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:2) ⋆ ᶠᶜmat3
+            ) ⋆ ᶜᶠmat ⋆ ᶠᶜmat ⋆ (
+                DiagonalMatrixRow(ᶜlg.gⁱʲ.components.data.:1) ⋆ ᶜᶠmat2 +
+                DiagonalMatrixRow(ᶜlg.gⁱʲ.components.data.:2) ⋆ ᶜᶠmat3
             ))
 
     unit_test_field_broadcast(
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_3.jl b/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_3.jl
index 711634180a..0a2e2508d7 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_3.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_3.jl
@@ -13,19 +13,19 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
                  and covectors times matrix of numbers and vectors times \
                  vector of numbers" begin
 
-    bc = @lazy @. ᶜᶠmat_AC1_num ⋅ ᶠᶜmat_C12_AC1 ⋅ ᶜᶠmat_num_C12 ⋅ ᶠvec
+    bc = @lazy @. ᶜᶠmat_AC1_num ⋆ ᶠᶜmat_C12_AC1 ⋆ ᶜᶠmat_num_C12 ⋆ ᶠvec
     result = materialize(bc)
 
     ref_set_result! =
         result -> (@. result = tuple(
-            ᶜᶠmat ⋅ (
-                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:1) ⋅ ᶠᶜmat2 +
-                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:2) ⋅ ᶠᶜmat3
-            ) ⋅ ᶜᶠmat ⋅ ᶠvec,
-            ᶜᶠmat ⋅ ᶠᶜmat ⋅ (
-                DiagonalMatrixRow(ᶜlg.gⁱʲ.components.data.:1) ⋅ ᶜᶠmat2 +
-                DiagonalMatrixRow(ᶜlg.gⁱʲ.components.data.:2) ⋅ ᶜᶠmat3
-            ) ⋅ ᶠvec,
+            ᶜᶠmat ⋆ (
+                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:1) ⋆ ᶠᶜmat2 +
+                DiagonalMatrixRow(ᶠlg.gⁱʲ.components.data.:2) ⋆ ᶠᶜmat3
+            ) ⋆ ᶜᶠmat ⋆ ᶠvec,
+            ᶜᶠmat ⋆ ᶠᶜmat ⋆ (
+                DiagonalMatrixRow(ᶜlg.gⁱʲ.components.data.:1) ⋆ ᶜᶠmat2 +
+                DiagonalMatrixRow(ᶜlg.gⁱʲ.components.data.:2) ⋆ ᶜᶠmat3
+            ) ⋆ ᶠvec,
         ))
 
     unit_test_field_broadcast(
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_4.jl b/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_4.jl
index 5211ca51c4..9d98d0ccbd 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_4.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_non_scalar_4.jl
@@ -13,14 +13,14 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
                  times matrix of numbers times matrix of numbers times \
                  vector of nested values" begin
 
-    bc = @lazy @. ᶜᶠmat_NT ⋅ ᶠᶜmat ⋅ ᶜᶠmat ⋅ ᶠᶜmat_NT ⋅ ᶜvec_NT
+    bc = @lazy @. ᶜᶠmat_NT ⋆ ᶠᶜmat ⋆ ᶜᶠmat ⋆ ᶠᶜmat_NT ⋆ ᶜvec_NT
     result = materialize(bc)
 
     ref_set_result! =
         result -> (@. result = nested_type(
-            ᶜᶠmat ⋅ ᶠᶜmat ⋅ ᶜᶠmat ⋅ ᶠᶜmat ⋅ ᶜvec,
-            ᶜᶠmat2 ⋅ ᶠᶜmat ⋅ ᶜᶠmat ⋅ ᶠᶜmat2 ⋅ ᶜvec,
-            ᶜᶠmat3 ⋅ ᶠᶜmat ⋅ ᶜᶠmat ⋅ ᶠᶜmat3 ⋅ ᶜvec,
+            ᶜᶠmat ⋆ ᶠᶜmat ⋆ ᶜᶠmat ⋆ ᶠᶜmat ⋆ ᶜvec,
+            ᶜᶠmat2 ⋆ ᶠᶜmat ⋆ ᶜᶠmat ⋆ ᶠᶜmat2 ⋆ ᶜvec,
+            ᶜᶠmat3 ⋆ ᶠᶜmat ⋆ ᶜᶠmat ⋆ ᶠᶜmat3 ⋆ ᶜvec,
         ))
 
     unit_test_field_broadcast(
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_1.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_1.jl
index fd22239f61..9736fa3dec 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_1.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_1.jl
@@ -9,7 +9,7 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 #! format: on
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "diagonal matrix times vector" begin
-    bc = @lazy @. ᶜᶜmat ⋅ ᶜvec
+    bc = @lazy @. ᶜᶜmat ⋆ ᶜvec
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜvec)
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_10.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_10.jl
index 9d2de3d963..baf76d5faf 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_10.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_10.jl
@@ -10,7 +10,7 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "diagonal matrix times bi-diagonal matrix times \
                  tri-diagonal matrix times quad-diagonal matrix times \
                  vector, but with forced right-associativity" begin
-    bc = @lazy @. ᶜᶜmat ⋅ (ᶜᶠmat ⋅ (ᶠᶠmat ⋅ (ᶠᶜmat ⋅ ᶜvec)))
+    bc = @lazy @. ᶜᶜmat ⋆ (ᶜᶠmat ⋆ (ᶠᶠmat ⋆ (ᶠᶜmat ⋆ ᶜvec)))
     if using_cuda
         @test_throws invalid_ir_error materialize(bc)
         @warn "cuda is broken for this test, exiting."
@@ -20,9 +20,9 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
 
     input_fields = (ᶜᶜmat, ᶜᶠmat, ᶠᶠmat, ᶠᶜmat, ᶜvec)
     temp_value_fields = (
-        (@. ᶠᶜmat ⋅ ᶜvec),
-        (@. ᶠᶠmat ⋅ (ᶠᶜmat ⋅ ᶜvec)),
-        (@. ᶜᶠmat ⋅ (ᶠᶠmat ⋅ (ᶠᶜmat ⋅ ᶜvec))),
+        (@. ᶠᶜmat ⋆ ᶜvec),
+        (@. ᶠᶠmat ⋆ (ᶠᶜmat ⋆ ᶜvec)),
+        (@. ᶜᶠmat ⋆ (ᶠᶠmat ⋆ (ᶠᶜmat ⋆ ᶜvec))),
     )
     ref_set_result! =
         (
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_11.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_11.jl
index c6f624fc34..75b52b09f4 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_11.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_11.jl
@@ -8,15 +8,15 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 #! format: on
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "linear combination of matrix products and LinearAlgebra.I" begin
-    bc = @lazy @. 2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,)
+    bc = @lazy @. 2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,)
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜᶠmat, ᶠᶠmat, ᶠᶜmat)
     temp_value_fields = (
         (@. 2 * ᶠᶜmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat),
-        (@. ᶠᶠmat ⋅ ᶠᶠmat),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat),
+        (@. ᶠᶠmat ⋆ ᶠᶠmat),
     )
     ref_set_result! =
         (
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_12.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_12.jl
index 08b0ddd186..e90001d1e2 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_12.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_12.jl
@@ -9,15 +9,15 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "another linear combination of matrix products and \
                  LinearAlgebra.I" begin
-    bc = @lazy @. ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋅ ᶠᶠmat + (4I,)
+    bc = @lazy @. ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋆ ᶠᶠmat + (4I,)
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜᶠmat, ᶠᶠmat, ᶠᶜmat)
     temp_value_fields = (
-        (@. ᶠᶜmat ⋅ ᶜᶜmat),
-        (@. ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat),
+        (@. ᶠᶜmat ⋆ ᶜᶜmat),
+        (@. ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat),
         (@. ᶠᶠmat / 3),
-        (@. (ᶠᶠmat / 3) ⋅ ᶠᶠmat),
+        (@. (ᶠᶠmat / 3) ⋆ ᶠᶠmat),
     )
     ref_set_result! =
         (
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_13.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_13.jl
index 187bf7e173..688778b759 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_13.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_13.jl
@@ -9,16 +9,16 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "matrix times linear combination" begin
     bc =
-        @lazy @. ᶜᶠmat ⋅ (2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,))
+        @lazy @. ᶜᶠmat ⋆ (2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,))
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜᶠmat, ᶠᶠmat, ᶠᶜmat)
     temp_value_fields = (
         (@. 2 * ᶠᶜmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat),
-        (@. ᶠᶠmat ⋅ ᶠᶠmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,)),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat),
+        (@. ᶠᶠmat ⋆ ᶠᶠmat),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,)),
     )
     ref_set_result! =
         (
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_14.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_14.jl
index e4e6a10826..8b2212866b 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_14.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_14.jl
@@ -8,8 +8,8 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 #! format: on
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "linear combination times another linear combination" begin
-    bc = @lazy @. (2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,)) ⋅
-             (ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋅ ᶠᶠmat + (4I,))
+    bc = @lazy @. (2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,)) ⋆
+             (ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋆ ᶠᶠmat + (4I,))
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜᶠmat, ᶠᶠmat, ᶠᶜmat)
@@ -49,15 +49,15 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
 
     temp_value_fields = (
         (@. 2 * ᶠᶜmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat),
-        (@. ᶠᶠmat ⋅ ᶠᶠmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,)),
-        (@. ᶠᶜmat ⋅ ᶜᶜmat),
-        (@. ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat),
+        (@. ᶠᶠmat ⋆ ᶠᶠmat),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,)),
+        (@. ᶠᶜmat ⋆ ᶜᶜmat),
+        (@. ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat),
         (@. ᶠᶠmat / 3),
-        (@. (ᶠᶠmat / 3) ⋅ ᶠᶠmat),
-        (@. ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋅ ᶠᶠmat + (4I,)),
+        (@. (ᶠᶠmat / 3) ⋆ ᶠᶠmat),
+        (@. ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋆ ᶠᶠmat + (4I,)),
     )
 
     unit_test_field_broadcast_vs_array_reference(
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_15.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_15.jl
index 5570f6e934..879920f257 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_15.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_15.jl
@@ -9,10 +9,10 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "matrix times matrix times linear combination times matrix \
                  times another linear combination times matrix" begin
-    bc = @lazy @. ᶠᶜmat ⋅ ᶜᶠmat ⋅
-             (2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,)) ⋅
-             ᶠᶠmat ⋅
-             (ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋅ ᶠᶠmat + (4I,)) ⋅
+    bc = @lazy @. ᶠᶜmat ⋆ ᶜᶠmat ⋆
+             (2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,)) ⋆
+             ᶠᶠmat ⋆
+             (ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋆ ᶠᶠmat + (4I,)) ⋆
              ᶠᶠmat
     result = materialize(bc)
 
@@ -60,26 +60,26 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
         end
 
     temp_value_fields = (
-        (@. ᶠᶜmat ⋅ ᶜᶠmat),
+        (@. ᶠᶜmat ⋆ ᶜᶠmat),
         (@. 2 * ᶠᶜmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat),
-        (@. ᶠᶠmat ⋅ ᶠᶠmat),
-        (@. 2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,)),
-        (@. ᶠᶜmat ⋅ ᶜᶠmat ⋅
-            (2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,))),
-        (@. ᶠᶜmat ⋅ ᶜᶠmat ⋅
-            (2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,)) ⋅
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat),
+        (@. ᶠᶠmat ⋆ ᶠᶠmat),
+        (@. 2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,)),
+        (@. ᶠᶜmat ⋆ ᶜᶠmat ⋆
+            (2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,))),
+        (@. ᶠᶜmat ⋆ ᶜᶠmat ⋆
+            (2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,)) ⋆
             ᶠᶠmat),
-        (@. ᶠᶜmat ⋅ ᶜᶜmat),
-        (@. ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat),
+        (@. ᶠᶜmat ⋆ ᶜᶜmat),
+        (@. ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat),
         (@. ᶠᶠmat / 3),
-        (@. (ᶠᶠmat / 3) ⋅ ᶠᶠmat),
-        (@. ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋅ ᶠᶠmat + (4I,)),
-        (@. ᶠᶜmat ⋅ ᶜᶠmat ⋅
-            (2 * ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat + ᶠᶠmat ⋅ ᶠᶠmat / 3 - (4I,)) ⋅
-            ᶠᶠmat ⋅
-            (ᶠᶜmat ⋅ ᶜᶜmat ⋅ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋅ ᶠᶠmat + (4I,))),
+        (@. (ᶠᶠmat / 3) ⋆ ᶠᶠmat),
+        (@. ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋆ ᶠᶠmat + (4I,)),
+        (@. ᶠᶜmat ⋆ ᶜᶠmat ⋆
+            (2 * ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat + ᶠᶠmat ⋆ ᶠᶠmat / 3 - (4I,)) ⋆
+            ᶠᶠmat ⋆
+            (ᶠᶜmat ⋆ ᶜᶜmat ⋆ ᶜᶠmat * 2 - (ᶠᶠmat / 3) ⋆ ᶠᶠmat + (4I,))),
     )
 
     unit_test_field_broadcast_vs_array_reference(
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_16.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_16.jl
index 9a46a2f0f4..762de2ca17 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_16.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_16.jl
@@ -8,9 +8,9 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 #! format: on
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "matrix constructions and multiplications" begin
-    bc = @lazy @. BidiagonalMatrixRow(ᶜᶠmat ⋅ ᶠvec, ᶜᶜmat ⋅ ᶜvec) ⋅
-             TridiagonalMatrixRow(ᶠvec, ᶠᶜmat ⋅ ᶜvec, 1) ⋅ ᶠᶠmat ⋅
-             DiagonalMatrixRow(DiagonalMatrixRow(ᶠvec) ⋅ ᶠvec)
+    bc = @lazy @. BidiagonalMatrixRow(ᶜᶠmat ⋆ ᶠvec, ᶜᶜmat ⋆ ᶜvec) ⋆
+             TridiagonalMatrixRow(ᶠvec, ᶠᶜmat ⋆ ᶜvec, 1) ⋆ ᶠᶠmat ⋆
+             DiagonalMatrixRow(DiagonalMatrixRow(ᶠvec) ⋆ ᶠvec)
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜᶠmat, ᶠᶠmat, ᶠᶜmat, ᶜvec, ᶠvec)
@@ -44,14 +44,14 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
         end
 
     temp_value_fields = (
-        (@. BidiagonalMatrixRow(ᶜᶠmat ⋅ ᶠvec, ᶜᶜmat ⋅ ᶜvec)),
-        (@. TridiagonalMatrixRow(ᶠvec, ᶠᶜmat ⋅ ᶜvec, 1)),
-        (@. BidiagonalMatrixRow(ᶜᶠmat ⋅ ᶠvec, ᶜᶜmat ⋅ ᶜvec) ⋅
-            TridiagonalMatrixRow(ᶠvec, ᶠᶜmat ⋅ ᶜvec, 1)),
-        (@. BidiagonalMatrixRow(ᶜᶠmat ⋅ ᶠvec, ᶜᶜmat ⋅ ᶜvec) ⋅
-            TridiagonalMatrixRow(ᶠvec, ᶠᶜmat ⋅ ᶜvec, 1) ⋅ ᶠᶠmat),
+        (@. BidiagonalMatrixRow(ᶜᶠmat ⋆ ᶠvec, ᶜᶜmat ⋆ ᶜvec)),
+        (@. TridiagonalMatrixRow(ᶠvec, ᶠᶜmat ⋆ ᶜvec, 1)),
+        (@. BidiagonalMatrixRow(ᶜᶠmat ⋆ ᶠvec, ᶜᶜmat ⋆ ᶜvec) ⋆
+            TridiagonalMatrixRow(ᶠvec, ᶠᶜmat ⋆ ᶜvec, 1)),
+        (@. BidiagonalMatrixRow(ᶜᶠmat ⋆ ᶠvec, ᶜᶜmat ⋆ ᶜvec) ⋆
+            TridiagonalMatrixRow(ᶠvec, ᶠᶜmat ⋆ ᶜvec, 1) ⋆ ᶠᶠmat),
         (@. DiagonalMatrixRow(ᶠvec)),
-        (@. DiagonalMatrixRow(DiagonalMatrixRow(ᶠvec) ⋅ ᶠvec)),
+        (@. DiagonalMatrixRow(DiagonalMatrixRow(ᶠvec) ⋆ ᶠvec)),
     )
     unit_test_field_broadcast_vs_array_reference(
         result,
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_2.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_2.jl
index ef579d3727..3d23fa84e8 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_2.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_2.jl
@@ -8,7 +8,7 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 #! format: on
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "tri-diagonal matrix times vector" begin
-    bc = @lazy @. ᶠᶠmat ⋅ ᶠvec
+    bc = @lazy @. ᶠᶠmat ⋆ ᶠvec
     result = materialize(bc)
 
     input_fields = (ᶠᶠmat, ᶠvec)
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_3.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_3.jl
index 2dc6056fbe..bda0742f80 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_3.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_3.jl
@@ -8,7 +8,7 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 #! format: on
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "quad-diagonal matrix times vector" begin
-    bc = @lazy @. ᶠᶜmat ⋅ ᶜvec
+    bc = @lazy @. ᶠᶜmat ⋆ ᶜvec
     result = materialize(bc)
 
     input_fields = (ᶠᶜmat, ᶜvec)
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_4.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_4.jl
index 93d53fc6e2..1daa08cac9 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_4.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_4.jl
@@ -8,7 +8,7 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 #! format: on
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "diagonal matrix times bi-diagonal matrix" begin
-    bc = @lazy @. ᶜᶜmat ⋅ ᶜᶠmat
+    bc = @lazy @. ᶜᶜmat ⋆ ᶜᶠmat
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜᶠmat)
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_5.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_5.jl
index fc47f30bd4..71c8925b06 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_5.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_5.jl
@@ -8,7 +8,7 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 #! format: on
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "tri-diagonal matrix times tri-diagonal matrix" begin
-    bc = @lazy @. ᶠᶠmat ⋅ ᶠᶠmat
+    bc = @lazy @. ᶠᶠmat ⋆ ᶠᶠmat
     result = materialize(bc)
 
     input_fields = (ᶠᶠmat,)
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_6.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_6.jl
index 4f43539453..469846855c 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_6.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_6.jl
@@ -8,7 +8,7 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 #! format: on
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "quad-diagonal matrix times diagonal matrix" begin
-    bc = @lazy @. ᶠᶜmat ⋅ ᶜᶜmat
+    bc = @lazy @. ᶠᶜmat ⋆ ᶜᶜmat
     result = materialize(bc)
 
     input_fields = (ᶠᶜmat, ᶜᶜmat)
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_7.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_7.jl
index b59b8beb54..4bf39f1fc7 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_7.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_7.jl
@@ -9,11 +9,11 @@ include(joinpath(pkgdir(ClimaCore),"test","MatrixFields","matrix_fields_broadcas
 test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "diagonal matrix times bi-diagonal matrix times \
                  tri-diagonal matrix times quad-diagonal matrix" begin
-    bc = @lazy @. ᶜᶜmat ⋅ ᶜᶠmat ⋅ ᶠᶠmat ⋅ ᶠᶜmat
+    bc = @lazy @. ᶜᶜmat ⋆ ᶜᶠmat ⋆ ᶠᶠmat ⋆ ᶠᶜmat
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜᶠmat, ᶠᶠmat, ᶠᶜmat)
-    temp_value_fields = ((@. ᶜᶜmat ⋅ ᶜᶠmat), (@. ᶜᶜmat ⋅ ᶜᶠmat ⋅ ᶠᶠmat))
+    temp_value_fields = ((@. ᶜᶜmat ⋆ ᶜᶠmat), (@. ᶜᶜmat ⋆ ᶜᶠmat ⋆ ᶠᶠmat))
     ref_set_result! =
         (_result, _ᶜᶜmat, _ᶜᶠmat, _ᶠᶠmat, _ᶠᶜmat, _temp1, _temp2) -> begin
             mul!(_temp1, _ᶜᶜmat, _ᶜᶠmat)
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_8.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_8.jl
index 8608f34d2b..d12b7d878e 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_8.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_8.jl
@@ -10,7 +10,7 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "diagonal matrix times bi-diagonal matrix times \
                  tri-diagonal matrix times quad-diagonal matrix, but with \
                  forced right-associativity" begin
-    bc = @lazy @. ᶜᶜmat ⋅ (ᶜᶠmat ⋅ (ᶠᶠmat ⋅ ᶠᶜmat))
+    bc = @lazy @. ᶜᶜmat ⋆ (ᶜᶠmat ⋆ (ᶠᶠmat ⋆ ᶠᶜmat))
     if using_cuda
         @test_throws invalid_ir_error materialize(bc)
         @warn "cuda is broken for this test, exiting."
@@ -19,7 +19,7 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜᶠmat, ᶠᶠmat, ᶠᶜmat)
-    temp_value_fields = ((@. ᶠᶠmat ⋅ ᶠᶜmat), (@. ᶜᶠmat ⋅ (ᶠᶠmat ⋅ ᶠᶜmat)))
+    temp_value_fields = ((@. ᶠᶠmat ⋆ ᶠᶜmat), (@. ᶜᶠmat ⋆ (ᶠᶠmat ⋆ ᶠᶜmat)))
     ref_set_result! =
         (_result, _ᶜᶜmat, _ᶜᶠmat, _ᶠᶠmat, _ᶠᶜmat, _temp1, _temp2) -> begin
             mul!(_temp1, _ᶠᶠmat, _ᶠᶜmat)
diff --git a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_9.jl b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_9.jl
index 0f71848f34..5999586c8a 100644
--- a/test/MatrixFields/matrix_fields_broadcasting/test_scalar_9.jl
+++ b/test/MatrixFields/matrix_fields_broadcasting/test_scalar_9.jl
@@ -10,14 +10,14 @@ test_opt = get(ENV, "BUILDKITE", "") == "true"
 @testset "diagonal matrix times bi-diagonal matrix times \
                  tri-diagonal matrix times quad-diagonal matrix times \
                  vector" begin
-    bc = @lazy @. ᶜᶜmat ⋅ ᶜᶠmat ⋅ ᶠᶠmat ⋅ ᶠᶜmat ⋅ ᶜvec
+    bc = @lazy @. ᶜᶜmat ⋆ ᶜᶠmat ⋆ ᶠᶠmat ⋆ ᶠᶜmat ⋆ ᶜvec
     result = materialize(bc)
 
     input_fields = (ᶜᶜmat, ᶜᶠmat, ᶠᶠmat, ᶠᶜmat, ᶜvec)
     temp_value_fields = (
-        (@. ᶜᶜmat ⋅ ᶜᶠmat),
-        (@. ᶜᶜmat ⋅ ᶜᶠmat ⋅ ᶠᶠmat),
-        (@. ᶜᶜmat ⋅ ᶜᶠmat ⋅ ᶠᶠmat ⋅ ᶠᶜmat),
+        (@. ᶜᶜmat ⋆ ᶜᶠmat),
+        (@. ᶜᶜmat ⋆ ᶜᶠmat ⋆ ᶠᶠmat),
+        (@. ᶜᶜmat ⋆ ᶜᶠmat ⋆ ᶠᶠmat ⋆ ᶠᶜmat),
     )
     ref_set_result! =
         (
diff --git a/test/MatrixFields/matrix_multiplication_at_boundaries.jl b/test/MatrixFields/matrix_multiplication_at_boundaries.jl
index 0f54ffdaf1..535a2590d4 100644
--- a/test/MatrixFields/matrix_multiplication_at_boundaries.jl
+++ b/test/MatrixFields/matrix_multiplication_at_boundaries.jl
@@ -65,6 +65,6 @@ end
         (ᶠᶠmatrix_without_outside_entries, ᶠᶠmatrix_with_outside_entries),
         (ᶠᶠmatrix_without_outside_entries, ᶠᶜmatrix_with_outside_entries),
     )
-        @test !any(isnan, parent(@. matrix_field1 ⋅ matrix_field2))
+        @test !any(isnan, parent(@. matrix_field1 ⋆ matrix_field2))
     end
 end
diff --git a/test/MatrixFields/matrix_multiplication_recursion.jl b/test/MatrixFields/matrix_multiplication_recursion.jl
index 2e504f33d0..49fb7de026 100644
--- a/test/MatrixFields/matrix_multiplication_recursion.jl
+++ b/test/MatrixFields/matrix_multiplication_recursion.jl
@@ -1,17 +1,11 @@
 import ClimaCore
-if !(@isdefined(TU))
-    include(
-        joinpath(
-            pkgdir(ClimaCore),
-            "test",
-            "TestUtilities",
-            "TestUtilities.jl",
-        ),
-    )
-end
+@isdefined(TU) || include(
+    joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
+);
 import .TestUtilities as TU
 import ClimaCore: Spaces, Geometry, Operators, Fields, MatrixFields
-import ClimaCore.MatrixFields: ⋅
+import ClimaCore.MatrixFields: ⋆
+
 import ClimaComms
 ClimaComms.@import_required_backends
 
@@ -74,9 +68,9 @@ fill!(parent(dest_field), NaN)
 # Test field update with multiple nested operations
 function update_field!(dest_field, K, dψdϑ_res)
     @. dest_field =
-        divf2c_matrix() ⋅ (
-            MatrixFields.DiagonalMatrixRow(interpc2f_op(-K)) ⋅
-            gradc2f_matrix() ⋅ MatrixFields.DiagonalMatrixRow(dψdϑ_res) +
+        divf2c_matrix() ⋆ (
+            MatrixFields.DiagonalMatrixRow(interpc2f_op(-K)) ⋆
+            gradc2f_matrix() ⋆ MatrixFields.DiagonalMatrixRow(dψdϑ_res) +
             MatrixFields.LowerDiagonalMatrixRow(
                 topBC_op(Geometry.Covariant3Vector(zero(interpc2f_op(K)))),
             )
diff --git a/test/MatrixFields/operator_matrices.jl b/test/MatrixFields/operator_matrices.jl
index f7642718cb..5516ba2622 100644
--- a/test/MatrixFields/operator_matrices.jl
+++ b/test/MatrixFields/operator_matrices.jl
@@ -7,6 +7,7 @@ import LinearAlgebra: I
 
 import ClimaCore.RecursiveApply: rzero
 import ClimaCore
+import ClimaCore.MatrixFields: ⋆
 import ClimaCore.Operators:
     SetValue,
     SetGradient,
@@ -41,13 +42,13 @@ import ClimaCore.Operators:
 
 include("matrix_field_test_utils.jl")
 
-apply_op_matrix(::Nothing, op_matrix, arg) = @lazy @. op_matrix() ⋅ arg
+apply_op_matrix(::Nothing, op_matrix, arg) = @lazy @. op_matrix() ⋆ arg
 apply_op_matrix(boundary_op, op_matrix, arg) =
-    @lazy @. boundary_op(op_matrix() ⋅ arg)
+    @lazy @. boundary_op(op_matrix() ⋆ arg)
 apply_op_matrix(::Nothing, op_matrix, arg1, arg2) =
-    @lazy @. op_matrix(arg1) ⋅ arg2
+    @lazy @. op_matrix(arg1) ⋆ arg2
 apply_op_matrix(boundary_op, op_matrix, arg1, arg2) =
-    @lazy @. boundary_op(op_matrix(arg1) ⋅ arg2)
+    @lazy @. boundary_op(op_matrix(arg1) ⋆ arg2)
 
 apply_op(::Nothing, op, args...) = @lazy @. op(args...)
 apply_op(boundary_op, op, args...) = @lazy @. boundary_op(op(args...))
@@ -225,36 +226,36 @@ end
     ᶜdiv_matrix = MatrixFields.operator_matrix(ᶜdiv)
     ᶠcurl_matrix = MatrixFields.operator_matrix(ᶠcurl)
 
-    @test_throws "does not contain any Fields" @. ᶜlbias_matrix() ⋅
+    @test_throws "does not contain any Fields" @. ᶜlbias_matrix() ⋆
                                                   ᶠinterp_matrix()
 
     ᶜ0 = @. zero(ᶜscalar)
     ᶜ1 = @. one(ᶜscalar)
     ᶠ1 = @. one(ᶠscalar)
     for get_result in (
-        @lazy(@. ᶜlbias_matrix() ⋅ ᶠinterp_matrix() + DiagonalMatrixRow(ᶜ0)),
-        @lazy(@. DiagonalMatrixRow(ᶜ0) + ᶜlbias_matrix() ⋅ ᶠinterp_matrix()),
-        @lazy(@. ᶜlbias_matrix() ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(ᶜ1)),
-        @lazy(@. ᶜlbias_matrix() ⋅ DiagonalMatrixRow(ᶠ1) ⋅ ᶠinterp_matrix()),
-        @lazy(@. DiagonalMatrixRow(ᶜ1) ⋅ ᶜlbias_matrix() ⋅ ᶠinterp_matrix()),
+        @lazy(@. ᶜlbias_matrix() ⋆ ᶠinterp_matrix() + DiagonalMatrixRow(ᶜ0)),
+        @lazy(@. DiagonalMatrixRow(ᶜ0) + ᶜlbias_matrix() ⋆ ᶠinterp_matrix()),
+        @lazy(@. ᶜlbias_matrix() ⋆ ᶠinterp_matrix() ⋆ DiagonalMatrixRow(ᶜ1)),
+        @lazy(@. ᶜlbias_matrix() ⋆ DiagonalMatrixRow(ᶠ1) ⋆ ᶠinterp_matrix()),
+        @lazy(@. DiagonalMatrixRow(ᶜ1) ⋆ ᶜlbias_matrix() ⋆ ᶠinterp_matrix()),
     )
         test_field_broadcast(;
             test_name = "product of two lazy operator matrices",
             get_result,
-            set_result = @lazy(@. ᶜlbias_matrix() ⋅ ᶠinterp_matrix()),
+            set_result = @lazy(@. ᶜlbias_matrix() ⋆ ᶠinterp_matrix()),
         )
     end
 
     test_field_broadcast(;
         test_name = "product of six operator matrices",
         get_result = @lazy(
-            @. ᶜflux_correct_matrix(ᶠuvw) ⋅ ᶜadvect_matrix(ᶠuvw) ⋅
-               ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠrbias_matrix() ⋅ ᶜlbias_matrix() ⋅
+            @. ᶜflux_correct_matrix(ᶠuvw) ⋆ ᶜadvect_matrix(ᶠuvw) ⋆
+               ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠrbias_matrix() ⋆ ᶜlbias_matrix() ⋆
                ᶠinterp_matrix()
         ),
         set_result = @lazy(
-            @. ᶜflux_correct_matrix(ᶠuvw) ⋅ ᶜadvect_matrix(ᶠuvw) ⋅
-               ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠrbias_matrix() ⋅ ᶜlbias_matrix() ⋅
+            @. ᶜflux_correct_matrix(ᶠuvw) ⋆ ᶜadvect_matrix(ᶠuvw) ⋆
+               ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠrbias_matrix() ⋆ ᶜlbias_matrix() ⋆
                ᶠinterp_matrix()
         ),
     )
@@ -263,14 +264,14 @@ end
         test_name = "applying six operators to a nested field using operator \
                      matrices",
         get_result = @lazy(
-            @. ᶜflux_correct_matrix(ᶠuvw) ⋅ ᶜadvect_matrix(ᶠuvw) ⋅
-               ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠrbias_matrix() ⋅ ᶜlbias_matrix() ⋅
-               ᶠinterp_matrix() ⋅ ᶜnested
+            @. ᶜflux_correct_matrix(ᶠuvw) ⋆ ᶜadvect_matrix(ᶠuvw) ⋆
+               ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠrbias_matrix() ⋆ ᶜlbias_matrix() ⋆
+               ᶠinterp_matrix() ⋆ ᶜnested
         ),
         set_result = @lazy(
-            @. ᶜflux_correct_matrix(ᶠuvw) ⋅ ᶜadvect_matrix(ᶠuvw) ⋅
-               ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠrbias_matrix() ⋅ ᶜlbias_matrix() ⋅
-               ᶠinterp_matrix() ⋅ ᶜnested
+            @. ᶜflux_correct_matrix(ᶠuvw) ⋆ ᶜadvect_matrix(ᶠuvw) ⋆
+               ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠrbias_matrix() ⋆ ᶜlbias_matrix() ⋆
+               ᶠinterp_matrix() ⋆ ᶜnested
         ),
         ref_set_result = @lazy(
             @. ᶜflux_correct(
@@ -287,21 +288,21 @@ end
         test_name = "applying six operators to a nested field using operator \
                      matrices, but with forced right associativity",
         get_result = @lazy(
-            @. ᶜflux_correct_matrix(ᶠuvw) ⋅ (
-                ᶜadvect_matrix(ᶠuvw) ⋅ (
-                    ᶜwinterp_matrix(ᶠscalar) ⋅ (
-                        ᶠrbias_matrix() ⋅
-                        (ᶜlbias_matrix() ⋅ (ᶠinterp_matrix() ⋅ ᶜnested))
+            @. ᶜflux_correct_matrix(ᶠuvw) ⋆ (
+                ᶜadvect_matrix(ᶠuvw) ⋆ (
+                    ᶜwinterp_matrix(ᶠscalar) ⋆ (
+                        ᶠrbias_matrix() ⋆
+                        (ᶜlbias_matrix() ⋆ (ᶠinterp_matrix() ⋆ ᶜnested))
                     )
                 )
             )
         ),
         set_result = @lazy(
-            @. ᶜflux_correct_matrix(ᶠuvw) ⋅ (
-                ᶜadvect_matrix(ᶠuvw) ⋅ (
-                    ᶜwinterp_matrix(ᶠscalar) ⋅ (
-                        ᶠrbias_matrix() ⋅
-                        (ᶜlbias_matrix() ⋅ (ᶠinterp_matrix() ⋅ ᶜnested))
+            @. ᶜflux_correct_matrix(ᶠuvw) ⋆ (
+                ᶜadvect_matrix(ᶠuvw) ⋆ (
+                    ᶜwinterp_matrix(ᶠscalar) ⋆ (
+                        ᶠrbias_matrix() ⋆
+                        (ᶜlbias_matrix() ⋆ (ᶠinterp_matrix() ⋆ ᶜnested))
                     )
                 )
             )
@@ -325,13 +326,13 @@ end
     # false positive, a compiler issue, or a sign that the code can be improved?
     for get_result in (
         @lazy(
-            @. (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠcurl_matrix() *
+            @. (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠcurl_matrix() *
                (c12_a,) +
                (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5
         ),
         @lazy(
-            @. ᶜdiv_matrix() ⋅ DiagonalMatrixRow(ᶠscalar) ⋅ ᶠgrad_matrix() ⋅ (
-                (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠcurl_matrix() *
+            @. ᶜdiv_matrix() ⋆ DiagonalMatrixRow(ᶠscalar) ⋆ ᶠgrad_matrix() ⋆ (
+                (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠcurl_matrix() *
                 (c12_a,) +
                 (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5
             )
@@ -345,20 +346,20 @@ end
         test_name = "non-trivial combination of operator matrices and other \
                      matrix fields",
         get_result = @lazy(
-            @. ᶠupwind_matrix(ᶠuvw) ⋅ (
-                ᶜdiv_matrix() ⋅ DiagonalMatrixRow(ᶠscalar) ⋅ ᶠgrad_matrix() ⋅
+            @. ᶠupwind_matrix(ᶠuvw) ⋆ (
+                ᶜdiv_matrix() ⋆ DiagonalMatrixRow(ᶠscalar) ⋆ ᶠgrad_matrix() ⋆
                 (
-                    (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠcurl_matrix() *
+                    (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠcurl_matrix() *
                     (c12_a,) +
                     (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5
                 ) - (2I,)
             )
         ),
         set_result = @lazy(
-            @. ᶠupwind_matrix(ᶠuvw) ⋅ (
-                ᶜdiv_matrix() ⋅ DiagonalMatrixRow(ᶠscalar) ⋅ ᶠgrad_matrix() ⋅
+            @. ᶠupwind_matrix(ᶠuvw) ⋆ (
+                ᶜdiv_matrix() ⋆ DiagonalMatrixRow(ᶠscalar) ⋆ ᶠgrad_matrix() ⋆
                 (
-                    (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠcurl_matrix() *
+                    (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠcurl_matrix() *
                     (c12_a,) +
                     (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5
                 ) - (2I,)
@@ -368,31 +369,31 @@ end
 
     # TODO: This case's reference function takes too long to compile on both
     # CPUs and GPUs (more than half an hour), as of Julia 1.9. This might be
-    # happening because of excessive inlining---aside from ⋅, all other finite
+    # happening because of excessive inlining---aside from ⋆, all other finite
     # difference operators use @propagate_inbounds. So, the reference function
     # is currently disabled, although the test does pass when it is enabled.
     test_field_broadcast(;
         test_name = "applying a non-trivial sequence of operations to a scalar \
                      field using operator matrices and other matrix fields",
         get_result = @lazy(
-            @. ᶠupwind_matrix(ᶠuvw) ⋅ (
-                ᶜdiv_matrix() ⋅ DiagonalMatrixRow(ᶠscalar) ⋅ ᶠgrad_matrix() ⋅
+            @. ᶠupwind_matrix(ᶠuvw) ⋆ (
+                ᶜdiv_matrix() ⋆ DiagonalMatrixRow(ᶠscalar) ⋆ ᶠgrad_matrix() ⋆
                 (
-                    (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠcurl_matrix() *
+                    (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠcurl_matrix() *
                     (c12_a,) +
                     (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5
                 ) - (2I,)
-            ) ⋅ ᶜscalar
+            ) ⋆ ᶜscalar
         ),
         set_result = @lazy(
-            @. ᶠupwind_matrix(ᶠuvw) ⋅ (
-                ᶜdiv_matrix() ⋅ DiagonalMatrixRow(ᶠscalar) ⋅ ᶠgrad_matrix() ⋅
+            @. ᶠupwind_matrix(ᶠuvw) ⋆ (
+                ᶜdiv_matrix() ⋆ DiagonalMatrixRow(ᶠscalar) ⋆ ᶠgrad_matrix() ⋆
                 (
-                    (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋅ ᶠcurl_matrix() *
+                    (c12_b',) * ᶜwinterp_matrix(ᶠscalar) ⋆ ᶠcurl_matrix() *
                     (c12_a,) +
                     (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5
                 ) - (2I,)
-            ) ⋅ ᶜscalar
+            ) ⋆ ᶜscalar
         ),
         # ref_set_result = @lazy(@. ᶠupwind(
         #     ᶠuvw,
diff --git a/test/Operators/finitedifference/benchmark_stencils_utils.jl b/test/Operators/finitedifference/benchmark_stencils_utils.jl
index b5aedab5b8..42c28a6a69 100755
--- a/test/Operators/finitedifference/benchmark_stencils_utils.jl
+++ b/test/Operators/finitedifference/benchmark_stencils_utils.jl
@@ -10,9 +10,9 @@ using ClimaCore.Geometry: ⊗
 import ClimaCore.DataLayouts
 
 import ClimaCore
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
+);
 import .TestUtilities as TU
 
 @show ClimaComms.device() isa ClimaComms.CUDADevice
diff --git a/test/Operators/integrals.jl b/test/Operators/integrals.jl
index da7f893ba7..45162fe132 100644
--- a/test/Operators/integrals.jl
+++ b/test/Operators/integrals.jl
@@ -12,9 +12,9 @@ import ClimaCore.Operators:
     column_reduce!,
     column_accumulate!
 
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
+);
 import .TestUtilities as TU
 
 are_boundschecks_forced = Base.JLOptions().check_bounds == 1
diff --git a/test/Operators/spectralelement/covar_deriv_ops.jl b/test/Operators/spectralelement/covar_deriv_ops.jl
index 20cb379f17..445eefad2c 100644
--- a/test/Operators/spectralelement/covar_deriv_ops.jl
+++ b/test/Operators/spectralelement/covar_deriv_ops.jl
@@ -25,7 +25,6 @@ zlim = (FT(0), FT(1))
 helem = 5
 velem = 5
 npoly = 5
-ndims = 3
 stretch = Meshes.Uniform()
 comms_context = ClimaComms.SingletonCommsContext(device)
 FT = eltype(xlim)
diff --git a/test/Spaces/extruded_cuda.jl b/test/Spaces/extruded_cuda.jl
index ef1766a6f0..03764d1b66 100644
--- a/test/Spaces/extruded_cuda.jl
+++ b/test/Spaces/extruded_cuda.jl
@@ -12,9 +12,9 @@ import ClimaCore:
     Domains, Topologies, Meshes, Spaces, Geometry, column, Fields, Grids
 using Test
 
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
+);
 import .TestUtilities as TU
 
 compare(cpu, gpu) = all(parent(cpu) .≈ Array(parent(gpu)))
diff --git a/test/Spaces/opt_spaces.jl b/test/Spaces/opt_spaces.jl
index d236242ea7..a652711c69 100644
--- a/test/Spaces/opt_spaces.jl
+++ b/test/Spaces/opt_spaces.jl
@@ -6,9 +6,9 @@ using Revise; include(joinpath("test", "Spaces", "opt_spaces.jl"))
 import ClimaCore
 import ClimaCore: Spaces, Grids, Topologies
 using Test
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
+);
 import .TestUtilities as TU
 import ClimaComms
 ClimaComms.@import_required_backends
diff --git a/test/Spaces/unit_dss.jl b/test/Spaces/unit_dss.jl
index c079433776..13e1daa457 100644
--- a/test/Spaces/unit_dss.jl
+++ b/test/Spaces/unit_dss.jl
@@ -20,9 +20,9 @@ import ClimaCore:
     Topologies,
     DataLayouts
 
-include(
+@isdefined(TU) || include(
     joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
-)
+);
 import .TestUtilities as TU
 
 function get_space_cs(::Type{FT}; context, R = 300.0) where {FT}
@@ -34,7 +34,7 @@ function get_space_cs(::Type{FT}; context, R = 300.0) where {FT}
     return space
 end
 
-function one_to_n(a::AbstractArray)
+function one_to_n_dss(a::AbstractArray)
     _a = Array(a)
     Random.seed!(1234)
     for i in 1:length(_a)
@@ -44,7 +44,7 @@ function one_to_n(a::AbstractArray)
 end
 
 function test_dss_count(f::Fields.Field, buff::Topologies.DSSBuffer, nc)
-    parent(f) .= one_to_n(parent(f))
+    parent(f) .= one_to_n_dss(parent(f))
     @test allunique(parent(f))
     cf = copy(f)
     Spaces.weighted_dss!(f => buff)
@@ -53,7 +53,7 @@ function test_dss_count(f::Fields.Field, buff::Topologies.DSSBuffer, nc)
     return (; n_dss_affected)
 end
 
-function get_space_and_buffers(::Type{FT}; context) where {FT}
+function get_space_and_buffers3(::Type{FT}; context) where {FT}
     init_state_covariant12(local_geometry, p) =
         Geometry.Covariant12Vector(1.0, -1.0)
     init_state_covariant123(local_geometry, p) =
@@ -102,7 +102,7 @@ end
 @testset "DSS of AxisTensors on Cubed Sphere" begin
     FT = Float64
     device = ClimaComms.device()
-    nt = get_space_and_buffers(FT; context = ClimaComms.context(device))
+    nt = get_space_and_buffers3(FT; context = ClimaComms.context(device))
 
     # test DSS for a Covariant12Vector
     # ensure physical velocity is continuous across SE boundary for initial state
diff --git a/test/runtests.jl b/test/runtests.jl
index 7d82c25493..f9dea243ba 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 #=
-julia --project
+julia --check-bounds=yes --project
+using TestEnv; TestEnv.activate();
 using Revise; include(joinpath("test", "runtests.jl"))
 =#
 using Test
@@ -140,7 +141,7 @@ fail_fast = true
 # Use prevent_leaky_tests = !isnothing(get(ENV, "CI", nothing))
 # once https://github.com/CliMA/ClimaCore.jl/issues/1826 is fixed:
 # prevent_leaky_tests = !isnothing(get(ENV, "CI", nothing)) # prevent leaky tests on CI)
-prevent_leaky_tests = true
+prevent_leaky_tests = false
 # tabulate_tests(unit_tests; include_timings = false) # uncomment for preview
 
 # If `fail_fast` is `true`, then let's error on invalid tests before starting them.
diff --git a/test/tabulated_tests.jl b/test/tabulated_tests.jl
index 391bc856fb..6aaf42f1c8 100644
--- a/test/tabulated_tests.jl
+++ b/test/tabulated_tests.jl
@@ -224,6 +224,7 @@ function run_unit_test!(
     prevent_leaky_tests::Bool = false,
 )
     @debug "--- About to test $(test.filename)"
+    # @info "--- About to test $(test.filename)"
     stats = if prevent_leaky_tests
         # This is a home-brewed `SafeTestsets`, but it allows us to (inside a module)
         @timevd eval(Meta.parse(test_expr_safe(test)))
@@ -231,6 +232,7 @@ function run_unit_test!(
         @timevd eval(Meta.parse(test_expr(test)))
         # test.elapsed = @elapsed include(test.filename)
     end
+    # @info "--- finished evaling test $(test.filename)"
     (; compile_elapsedtimes, elapsed) = stats
     compile_time = first(compile_elapsedtimes) / 1e9
     recompile_time = last(compile_elapsedtimes) / 1e9