Avoid adjoint sparse matrix times Diagonal (#179)

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Author: JoshuaLampert

src/equations/serre_green_naghdi_1d.jl

@@ -268,7 +268,7 @@ function assemble_system_matrix!(cache, h, D1, D1mat,
     # is not necessarily perfectly symmetric but only up to
     # round-off errors. We wrap it here to avoid issues with the
     # factorization.
-    return Symmetric(Diagonal(M_h) + D1mat' * Diagonal(M_h3_3) * D1mat)
+    return Symmetric(Diagonal(M_h) + D1mat' * (Diagonal(M_h3_3) * D1mat))
 end
 
 # variable bathymetry

src/equations/svaerd_kalisch_1d.jl

@@ -264,15 +264,15 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D,
        D1 isa FourierDerivativeOperator
         D1_central = D1
         D1mat = sparse(D1_central)
-        minus_MD1betaD1 = D1mat' * Diagonal(M_beta) * D1mat
+        minus_MD1betaD1 = D1mat' * (Diagonal(M_beta) * D1mat)
         system_matrix = let cache = (; D1, M_h, minus_MD1betaD1)
             assemble_system_matrix!(cache, h,
                                     equations, boundary_conditions)
         end
     elseif D1 isa PeriodicUpwindOperators
         D1_central = D1.central
         D1mat_minus = sparse(D1.minus)
-        minus_MD1betaD1 = D1mat_minus' * Diagonal(M_beta) * D1mat_minus
+        minus_MD1betaD1 = D1mat_minus' * (Diagonal(M_beta) * D1mat_minus)
         system_matrix = let cache = (; D1, M_h, minus_MD1betaD1)
             assemble_system_matrix!(cache, h,
                                     equations, boundary_conditions)
@@ -326,16 +326,17 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D,
        D1 isa UniformCoupledOperator
         D1_central = D1
         D1mat = sparse(D1_central)
-        minus_MD1betaD1 = sparse(D1mat' * Diagonal(M_beta) *
-                                 D1mat * Pd)[(begin + 1):(end - 1), :]
+        minus_MD1betaD1 = sparse(D1mat' * (Diagonal(M_beta) *
+                                           D1mat * Pd))[(begin + 1):(end - 1), :]
         system_matrix = let cache = (; D1, M_h, minus_MD1betaD1)
             assemble_system_matrix!(cache, h, equations, boundary_conditions)
         end
     elseif D1 isa UpwindOperators
         D1_central = D1.central
         D1mat_minus = sparse(D1.minus)
-        minus_MD1betaD1 = sparse(D1mat_minus' * Diagonal(M_beta) *
-                                 D1mat_minus * Pd)[(begin + 1):(end - 1), :]
+        minus_MD1betaD1 = sparse(D1mat_minus' * (Diagonal(M_beta) *
+                                                 D1mat_minus * Pd))[(begin + 1):(end - 1),
+                                                                    :]
         system_matrix = let cache = (; D1, M_h, minus_MD1betaD1)
             assemble_system_matrix!(cache, h, equations, boundary_conditions)
         end