From 5d35dc3c5cc9f7568d5b721ffc8e0157201e2c6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 7 Nov 2023 12:31:40 +0100 Subject: [PATCH 1/4] Unified divb analyze routines for equations with `magnetic_field` --- src/callbacks_step/analysis_dg2d.jl | 72 ++++++------------- src/equations/ideal_glm_mhd_2d.jl | 4 ++ .../ideal_glm_mhd_multicomponent_2d.jl | 4 ++ 3 files changed, 29 insertions(+), 51 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index a9e0cf87b0..ef92cba5ed 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -246,30 +246,16 @@ 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 - -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]) + B1_kj, _, _ = magnetic_field(view(u,:, k, j, element), equations) + _, B2_ik, _ = magnetic_field(view(u,:, i, k, element), equations) + divb += (derivative_matrix[i, k] * B1_kj + + derivative_matrix[j, k] * B2_ik) end divb *= cache.elements.inverse_jacobian[element] divb^2 @@ -279,7 +265,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, @@ -290,10 +276,12 @@ 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(view(u,:, k, j, element), equations) + B1_ik, B2_ik, _ = magnetic_field(view(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 @@ -302,29 +290,7 @@ end function analyze(::Val{:linf_divb}, du, u, t, mesh::TreeMesh{2}, - equations::IdealGlmMhdEquations2D, dg::DGSEM, cache) - @unpack derivative_matrix, weights = dg.basis - - # 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[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) + equations, dg::DGSEM, cache) @unpack derivative_matrix, weights = dg.basis # integrate over all elements to get the divergence-free condition errors @@ -333,8 +299,10 @@ 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[5, k, j, element] + - derivative_matrix[j, k] * u[6, i, k, element]) + B1_kj, _, _ = magnetic_field(view(u,:, k, j, element), equations) + _, B2_ik, _ = magnetic_field(view(u,:, i, k, element), equations) + 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)) @@ -347,7 +315,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 @@ -363,10 +331,12 @@ function analyze(::Val{:linf_divb}, du, u, t, element) # Compute the transformed divergence for k in eachnode(dg) + B1_kj, B2_kj, _ = magnetic_field(view(u,:, k, j, element), equations) + B1_ik, B2_ik, _ = magnetic_field(view(u,:, i, k, element), equations) divb += (derivative_matrix[i, k] * - (Ja11 * u[6, k, j, element] + Ja12 * u[7, k, j, element]) + - derivative_matrix[j, k] * - (Ja21 * u[6, i, k, element] + Ja22 * u[7, i, k, element])) + (Ja11 * B1_kj + Ja12 * B2_kj) + + derivative_matrix[j, k] * + (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/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index e8de0cedde..ab29252dc9 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -1338,4 +1338,8 @@ end @inline function cross_helicity(cons, ::IdealGlmMhdEquations2D) return (cons[2] * cons[6] + cons[3] * cons[7] + cons[4] * cons[8]) / cons[1] end + +# Return the magnetic field +magnetic_field(u, equations::IdealGlmMhdEquations2D) = SVector(u[6], u[7], u[8]) + end # @muladd diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl index 9b0eeb411e..2d3db16913 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_2d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_2d.jl @@ -711,4 +711,8 @@ end return SVector{ncomponents(equations), real(equations)}(u[i + 8] * v for i in eachcomponent(equations)) end + +# Return the magnetic field +magnetic_field(u, equations::IdealGlmMhdMulticomponentEquations2D) = SVector(u[5], u[6], u[7]) + end # @muladd From d6ad2f0746fccbf2e626079014b855552e698412 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 7 Nov 2023 12:38:30 +0100 Subject: [PATCH 2/4] format --- src/equations/ideal_glm_mhd_2d.jl | 1 - src/equations/ideal_glm_mhd_multicomponent_2d.jl | 5 +++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index ab29252dc9..6bb9684cf6 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -1341,5 +1341,4 @@ end # Return the magnetic field magnetic_field(u, equations::IdealGlmMhdEquations2D) = SVector(u[6], u[7], u[8]) - end # @muladd diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl index 2d3db16913..18320a2de7 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_2d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_2d.jl @@ -713,6 +713,7 @@ end end # Return the magnetic field -magnetic_field(u, equations::IdealGlmMhdMulticomponentEquations2D) = SVector(u[5], u[6], u[7]) - +function magnetic_field(u, equations::IdealGlmMhdMulticomponentEquations2D) + SVector(u[5], u[6], u[7]) +end end # @muladd From 9e8dd983b6e637ea74f943cf19650986082e03fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 7 Nov 2023 12:41:27 +0100 Subject: [PATCH 3/4] format --- src/callbacks_step/analysis_dg2d.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index ef92cba5ed..707c73b5ef 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -252,8 +252,8 @@ function analyze(::Val{:l2_divb}, du, u, t, dg, cache, derivative_matrix divb = zero(eltype(u)) for k in eachnode(dg) - B1_kj, _, _ = magnetic_field(view(u,:, k, j, element), equations) - _, B2_ik, _ = magnetic_field(view(u,:, i, k, element), equations) + B1_kj, _, _ = magnetic_field(view(u, :, k, j, element), equations) + _, B2_ik, _ = magnetic_field(view(u, :, i, k, element), equations) divb += (derivative_matrix[i, k] * B1_kj + derivative_matrix[j, k] * B2_ik) end @@ -276,8 +276,8 @@ 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(view(u,:, k, j, element), equations) - B1_ik, B2_ik, _ = magnetic_field(view(u,:, i, k, element), equations) + B1_kj, B2_kj, _ = magnetic_field(view(u, :, k, j, element), equations) + B1_ik, B2_ik, _ = magnetic_field(view(u, :, i, k, element), equations) divb += (derivative_matrix[i, k] * (Ja11 * B1_kj + Ja12 * B2_kj) + derivative_matrix[j, k] * @@ -299,8 +299,8 @@ 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) - B1_kj, _, _ = magnetic_field(view(u,:, k, j, element), equations) - _, B2_ik, _ = magnetic_field(view(u,:, i, k, element), equations) + B1_kj, _, _ = magnetic_field(view(u, :, k, j, element), equations) + _, B2_ik, _ = magnetic_field(view(u, :, i, k, element), equations) divb += (derivative_matrix[i, k] * B1_kj + derivative_matrix[j, k] * B2_ik) end @@ -331,12 +331,12 @@ function analyze(::Val{:linf_divb}, du, u, t, element) # Compute the transformed divergence for k in eachnode(dg) - B1_kj, B2_kj, _ = magnetic_field(view(u,:, k, j, element), equations) - B1_ik, B2_ik, _ = magnetic_field(view(u,:, i, k, element), equations) + B1_kj, B2_kj, _ = magnetic_field(view(u, :, k, j, element), equations) + B1_ik, B2_ik, _ = magnetic_field(view(u, :, i, k, element), equations) divb += (derivative_matrix[i, k] * - (Ja11 * B1_kj + Ja12 * B2_kj) + - derivative_matrix[j, k] * - (Ja21 * B1_ik + Ja22 * B2_ik)) + (Ja11 * B1_kj + Ja12 * B2_kj) + + derivative_matrix[j, k] * + (Ja21 * B1_ik + Ja22 * B2_ik)) end divb *= cache.elements.inverse_jacobian[i, j, element] linf_divb = max(linf_divb, abs(divb)) From 34dc448659fab472882d363ba936d128d27ea7d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 7 Nov 2023 17:56:25 +0100 Subject: [PATCH 4/4] Changed format of `magnetic_field` for consistency --- src/equations/ideal_glm_mhd_2d.jl | 4 +++- src/equations/ideal_glm_mhd_multicomponent_2d.jl | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 6bb9684cf6..4092c7bc95 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -1340,5 +1340,7 @@ end end # Return the magnetic field -magnetic_field(u, equations::IdealGlmMhdEquations2D) = SVector(u[6], u[7], u[8]) +function magnetic_field(u, equations::IdealGlmMhdEquations2D) + return SVector(u[6], u[7], u[8]) +end end # @muladd diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl index 18320a2de7..74a5873f75 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_2d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_2d.jl @@ -714,6 +714,6 @@ end # Return the magnetic field function magnetic_field(u, equations::IdealGlmMhdMulticomponentEquations2D) - SVector(u[5], u[6], u[7]) + return SVector(u[5], u[6], u[7]) end end # @muladd