Unify analysis routines for div(B) to make them equation-agnostic (#2176

) * Unify analysis routines for div(B) to make them equation-agnostic * Apply suggestions from code review Co-authored-by: Andrew Winters <[email protected]> --------- Co-authored-by: Andrew Winters <[email protected]>
trixi-framework · Nov 23, 2024 · 62ff764 · 62ff764
1 parent b5711d4
commit 62ff764
Show file tree

Hide file tree

Showing 5 changed files with 67 additions and 70 deletions.
diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl
@@ -279,30 +279,17 @@ end
 
 function analyze(::Val{:l2_divb}, du, u, t,
                  mesh::TreeMesh{2},
-                 equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
+                 equations, dg::DGSEM, cache)
     integrate_via_indices(u, mesh, equations, dg, cache, cache,
                           dg.basis.derivative_matrix) do u, i, j, element, equations,
                                                          dg, cache, derivative_matrix
         divb = zero(eltype(u))
         for k in eachnode(dg)
-            divb += (derivative_matrix[i, k] * u[6, k, j, element] +
-                     derivative_matrix[j, k] * u[7, i, k, element])
-        end
-        divb *= cache.elements.inverse_jacobian[element]
-        divb^2
-    end |> sqrt
-end
+            B1_kj, _, _ = magnetic_field(u[:, k, j, element], equations)
+            _, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)
 
-function analyze(::Val{:l2_divb}, du, u, t,
-                 mesh::TreeMesh{2}, equations::IdealGlmMhdMulticomponentEquations2D,
-                 dg::DG, cache)
-    integrate_via_indices(u, mesh, equations, dg, cache, cache,
-                          dg.basis.derivative_matrix) do u, i, j, element, equations,
-                                                         dg, cache, derivative_matrix
-        divb = zero(eltype(u))
-        for k in eachnode(dg)
-            divb += (derivative_matrix[i, k] * u[5, k, j, element] +
-                     derivative_matrix[j, k] * u[6, i, k, element])
+            divb += (derivative_matrix[i, k] * B1_kj +
+                     derivative_matrix[j, k] * B2_ik)
         end
         divb *= cache.elements.inverse_jacobian[element]
         divb^2
@@ -312,7 +299,7 @@ end
 function analyze(::Val{:l2_divb}, du, u, t,
                  mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
                              T8codeMesh{2}},
-                 equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
+                 equations, dg::DGSEM, cache)
     @unpack contravariant_vectors = cache.elements
     integrate_via_indices(u, mesh, equations, dg, cache, cache,
                           dg.basis.derivative_matrix) do u, i, j, element, equations,
@@ -323,10 +310,13 @@ function analyze(::Val{:l2_divb}, du, u, t,
         Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
         # Compute the transformed divergence
         for k in eachnode(dg)
+            B1_kj, B2_kj, _ = magnetic_field(u[:, k, j, element], equations)
+            B1_ik, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)
+
             divb += (derivative_matrix[i, k] *
-                     (Ja11 * u[6, k, j, element] + Ja12 * u[7, k, j, element]) +
+                     (Ja11 * B1_kj + Ja12 * B2_kj) +
                      derivative_matrix[j, k] *
-                     (Ja21 * u[6, i, k, element] + Ja22 * u[7, i, k, element]))
+                     (Ja21 * B1_ik + Ja22 * B2_ik))
         end
         divb *= cache.elements.inverse_jacobian[i, j, element]
         divb^2
@@ -335,7 +325,7 @@ end
 
 function analyze(::Val{:linf_divb}, du, u, t,
                  mesh::TreeMesh{2},
-                 equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
+                 equations, dg::DGSEM, cache)
     @unpack derivative_matrix, weights = dg.basis
 
     # integrate over all elements to get the divergence-free condition errors
@@ -344,30 +334,11 @@ function analyze(::Val{:linf_divb}, du, u, t,
         for j in eachnode(dg), i in eachnode(dg)
             divb = zero(eltype(u))
             for k in eachnode(dg)
-                divb += (derivative_matrix[i, k] * u[6, k, j, element] +
-                         derivative_matrix[j, k] * u[7, i, k, element])
-            end
-            divb *= cache.elements.inverse_jacobian[element]
-            linf_divb = max(linf_divb, abs(divb))
-        end
-    end
-
-    return linf_divb
-end
-
-function analyze(::Val{:linf_divb}, du, u, t,
-                 mesh::TreeMesh{2}, equations::IdealGlmMhdMulticomponentEquations2D,
-                 dg::DG, cache)
-    @unpack derivative_matrix, weights = dg.basis
+                B1_kj, _, _ = magnetic_field(u[:, k, j, element], equations)
+                _, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)
 
-    # integrate over all elements to get the divergence-free condition errors
-    linf_divb = zero(eltype(u))
-    for element in eachelement(dg, cache)
-        for j in eachnode(dg), i in eachnode(dg)
-            divb = zero(eltype(u))
-            for k in eachnode(dg)
-                divb += (derivative_matrix[i, k] * u[5, k, j, element] +
-                         derivative_matrix[j, k] * u[6, i, k, element])
+                divb += (derivative_matrix[i, k] * B1_kj +
+                         derivative_matrix[j, k] * B2_ik)
             end
             divb *= cache.elements.inverse_jacobian[element]
             linf_divb = max(linf_divb, abs(divb))
@@ -380,7 +351,7 @@ end
 function analyze(::Val{:linf_divb}, du, u, t,
                  mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
                              T8codeMesh{2}},
-                 equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
+                 equations, dg::DGSEM, cache)
     @unpack derivative_matrix, weights = dg.basis
     @unpack contravariant_vectors = cache.elements
 
@@ -396,10 +367,13 @@ function analyze(::Val{:linf_divb}, du, u, t,
                                                   element)
             # Compute the transformed divergence
             for k in eachnode(dg)
+                B1_kj, B2_kj, _ = magnetic_field(u[:, k, j, element], equations)
+                B1_ik, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)
+
                 divb += (derivative_matrix[i, k] *
-                         (Ja11 * u[6, k, j, element] + Ja12 * u[7, k, j, element]) +
+                         (Ja11 * B1_kj + Ja12 * B2_kj) +
                          derivative_matrix[j, k] *
-                         (Ja21 * u[6, i, k, element] + Ja22 * u[7, i, k, element]))
+                         (Ja21 * B1_ik + Ja22 * B2_ik))
             end
             divb *= cache.elements.inverse_jacobian[i, j, element]
             linf_divb = max(linf_divb, abs(divb))

diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl
@@ -307,16 +307,20 @@ function analyze(::typeof(entropy_timederivative), du, u, t,
 end
 
 function analyze(::Val{:l2_divb}, du, u, t,
-                 mesh::TreeMesh{3}, equations::IdealGlmMhdEquations3D,
+                 mesh::TreeMesh{3}, equations,
                  dg::DGSEM, cache)
     integrate_via_indices(u, mesh, equations, dg, cache, cache,
                           dg.basis.derivative_matrix) do u, i, j, k, element, equations,
                                                          dg, cache, derivative_matrix
         divb = zero(eltype(u))
         for l in eachnode(dg)
-            divb += (derivative_matrix[i, l] * u[6, l, j, k, element] +
-                     derivative_matrix[j, l] * u[7, i, l, k, element] +
-                     derivative_matrix[k, l] * u[8, i, j, l, element])
+            B_ljk = magnetic_field(u[:, l, j, k, element], equations)
+            B_ilk = magnetic_field(u[:, i, l, k, element], equations)
+            B_ijl = magnetic_field(u[:, i, j, l, element], equations)
+
+            divb += (derivative_matrix[i, l] * B_ljk[1] +
+                     derivative_matrix[j, l] * B_ilk[2] +
+                     derivative_matrix[k, l] * B_ijl[3])
         end
         divb *= cache.elements.inverse_jacobian[element]
         divb^2
@@ -325,7 +329,7 @@ end
 
 function analyze(::Val{:l2_divb}, du, u, t,
                  mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
-                 equations::IdealGlmMhdEquations3D,
+                 equations,
                  dg::DGSEM, cache)
     @unpack contravariant_vectors = cache.elements
     integrate_via_indices(u, mesh, equations, dg, cache, cache,
@@ -341,23 +345,24 @@ function analyze(::Val{:l2_divb}, du, u, t,
                                                     element)
         # Compute the transformed divergence
         for l in eachnode(dg)
+            B_ljk = magnetic_field(u[:, l, j, k, element], equations)
+            B_ilk = magnetic_field(u[:, i, l, k, element], equations)
+            B_ijl = magnetic_field(u[:, i, j, l, element], equations)
+
             divb += (derivative_matrix[i, l] *
-                     (Ja11 * u[6, l, j, k, element] + Ja12 * u[7, l, j, k, element] +
-                      Ja13 * u[8, l, j, k, element]) +
+                     (Ja11 * B_ljk[1] + Ja12 * B_ljk[2] + Ja13 * B_ljk[3]) +
                      derivative_matrix[j, l] *
-                     (Ja21 * u[6, i, l, k, element] + Ja22 * u[7, i, l, k, element] +
-                      Ja23 * u[8, i, l, k, element]) +
+                     (Ja21 * B_ilk[1] + Ja22 * B_ilk[2] + Ja23 * B_ilk[3]) +
                      derivative_matrix[k, l] *
-                     (Ja31 * u[6, i, j, l, element] + Ja32 * u[7, i, j, l, element] +
-                      Ja33 * u[8, i, j, l, element]))
+                     (Ja31 * B_ijl[1] + Ja32 * B_ijl[2] + Ja33 * B_ijl[3]))
         end
         divb *= cache.elements.inverse_jacobian[i, j, k, element]
         divb^2
     end |> sqrt
 end
 
 function analyze(::Val{:linf_divb}, du, u, t,
-                 mesh::TreeMesh{3}, equations::IdealGlmMhdEquations3D,
+                 mesh::TreeMesh{3}, equations,
                  dg::DGSEM, cache)
     @unpack derivative_matrix, weights = dg.basis
 
@@ -367,9 +372,13 @@ function analyze(::Val{:linf_divb}, du, u, t,
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
             divb = zero(eltype(u))
             for l in eachnode(dg)
-                divb += (derivative_matrix[i, l] * u[6, l, j, k, element] +
-                         derivative_matrix[j, l] * u[7, i, l, k, element] +
-                         derivative_matrix[k, l] * u[8, i, j, l, element])
+                B_ljk = magnetic_field(u[:, l, j, k, element], equations)
+                B_ilk = magnetic_field(u[:, i, l, k, element], equations)
+                B_ijl = magnetic_field(u[:, i, j, l, element], equations)
+
+                divb += (derivative_matrix[i, l] * B_ljk[1] +
+                         derivative_matrix[j, l] * B_ilk[2] +
+                         derivative_matrix[k, l] * B_ijl[3])
             end
             divb *= cache.elements.inverse_jacobian[element]
             linf_divb = max(linf_divb, abs(divb))
@@ -386,7 +395,7 @@ end
 
 function analyze(::Val{:linf_divb}, du, u, t,
                  mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
-                 equations::IdealGlmMhdEquations3D,
+                 equations,
                  dg::DGSEM, cache)
     @unpack derivative_matrix, weights = dg.basis
     @unpack contravariant_vectors = cache.elements
@@ -405,12 +414,16 @@ function analyze(::Val{:linf_divb}, du, u, t,
                                                         k, element)
             # Compute the transformed divergence
             for l in eachnode(dg)
-                divb += (derivative_matrix[i, l] * (Ja11 * u[6, l, j, k, element] +
-                          Ja12 * u[7, l, j, k, element] + Ja13 * u[8, l, j, k, element]) +
-                         derivative_matrix[j, l] * (Ja21 * u[6, i, l, k, element] +
-                          Ja22 * u[7, i, l, k, element] + Ja23 * u[8, i, l, k, element]) +
-                         derivative_matrix[k, l] * (Ja31 * u[6, i, j, l, element] +
-                          Ja32 * u[7, i, j, l, element] + Ja33 * u[8, i, j, l, element]))
+                B_ljk = magnetic_field(u[:, l, j, k, element], equations)
+                B_ilk = magnetic_field(u[:, i, l, k, element], equations)
+                B_ijl = magnetic_field(u[:, i, j, l, element], equations)
+
+                divb += (derivative_matrix[i, l] * (Ja11 * B_ljk[1] +
+                          Ja12 * B_ljk[2] + Ja13 * B_ljk[3]) +
+                         derivative_matrix[j, l] * (Ja21 * B_ilk[1] +
+                          Ja22 * B_ilk[2] + Ja23 * B_ilk[3]) +
+                         derivative_matrix[k, l] * (Ja31 * B_ijl[1] +
+                          Ja32 * B_ijl[2] + Ja33 * B_ijl[3]))
             end
             divb *= cache.elements.inverse_jacobian[i, j, k, element]
             linf_divb = max(linf_divb, abs(divb))

diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
@@ -46,6 +46,9 @@ function default_analysis_integrals(::IdealGlmMhdEquations2D)
     (entropy_timederivative, Val(:l2_divb), Val(:linf_divb))
 end
 
+# Helper function to extract the magnetic field vector from the conservative variables
+magnetic_field(u, equations::IdealGlmMhdEquations2D) = SVector(u[6], u[7], u[8])
+
 # Set initial conditions at physical location `x` for time `t`
 """
     initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D)

diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl
@@ -44,6 +44,9 @@ function default_analysis_integrals(::IdealGlmMhdEquations3D)
     (entropy_timederivative, Val(:l2_divb), Val(:linf_divb))
 end
 
+# Helper function to extract the magnetic field vector from the conservative variables
+magnetic_field(u, equations::IdealGlmMhdEquations3D) = SVector(u[6], u[7], u[8])
+
 # Set initial conditions at physical location `x` for time `t`
 """
     initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D)

diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl
@@ -101,6 +101,10 @@ function default_analysis_integrals(::IdealGlmMhdMulticomponentEquations2D)
     (entropy_timederivative, Val(:l2_divb), Val(:linf_divb))
 end
 
+# Helper function to extract the magnetic field vector from the conservative variables
+magnetic_field(u, equations::IdealGlmMhdMulticomponentEquations2D) = SVector(u[5], u[6],
+                                                                             u[7])
+
 """
     initial_condition_convergence_test(x, t, equations::IdealGlmMhdMulticomponentEquations2D)