Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 2 commits into from
Nov 23, 2024
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
734e3f5
Unify analysis routines for div(B) to make them equation-agnostic
amrueda Nov 22, 2024
7d88f6e
Apply suggestions from code review
amrueda Nov 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 22 additions & 48 deletions src/callbacks_step/analysis_dg2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Comment on lines +288 to +289
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This allocates a lot. Can we please use get_node_vars instead?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I implemented a version of the TreeMesh{3} routines that uses get_node_vars here (new branch).

If I run the same benchmark as above, I get the exact same number of allocations as before:

julia> @benchmark Trixi.analyze(Val{:l2_divb}(), $du, $u, 0.0, $mesh, $equations, $solver, $semi.cache)
BenchmarkTools.Trial: 4430 samples with 1 evaluation.
 Range (min … max):  930.719 μs …   3.922 ms  ┊ GC (min … max):  0.00% … 68.17%
 Time  (median):     979.486 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):     1.126 ms ± 599.382 μs  ┊ GC (mean ± σ):  12.17% ± 15.93%

  ██▄▃▁                                                     ▁▁▂ ▁
  █████▃▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅███ █
  931 μs        Histogram: log(frequency) by time       3.73 ms <

 Memory estimate: 2.75 MiB, allocs estimate: 40976.

julia> @benchmark Trixi.analyze(Val{:linf_divb}(), $du, $u, 0.0, $mesh, $equations, $solver, $semi.cache)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  51.021 μs … 101.782 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     51.941 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   52.111 μs ±   2.429 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄ ▂▂ █▃▅█▁                                                   ▂
  █▄██▅██████▇▅█▅▇▇▆▅▅▇▅▅▄▆▄▆▆▅▇█▇▅▆▅▅▆▅▅▅▅▃▄▅▄▅▄▅▆▅▆▆▆▅▃▅▅▁▃▄ █
  51 μs         Histogram: log(frequency) by time        58 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.


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
Expand All @@ -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,
Expand All @@ -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)
Comment on lines +313 to +314
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here


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
Expand All @@ -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
Expand All @@ -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)
Comment on lines +337 to +338
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same


# 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))
Expand All @@ -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

Expand All @@ -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)
Comment on lines +370 to +371
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same


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))
Expand Down
57 changes: 35 additions & 22 deletions src/callbacks_step/analysis_dg3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Comment on lines +317 to +319
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same


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
Expand All @@ -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,
Expand All @@ -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)
Comment on lines +348 to +350
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same


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

Expand All @@ -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)
Comment on lines +375 to +377
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same


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))
Expand All @@ -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
Expand All @@ -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)
Comment on lines +417 to +419
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same


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))
Expand Down
3 changes: 3 additions & 0 deletions src/equations/ideal_glm_mhd_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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])
amrueda marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be exported and become official API like velocity, density etc.?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the future, that will be useful. I'll export it in a future PR.


# Set initial conditions at physical location `x` for time `t`
"""
initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D)
Expand Down
3 changes: 3 additions & 0 deletions src/equations/ideal_glm_mhd_3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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])
amrueda marked this conversation as resolved.
Show resolved Hide resolved

# Set initial conditions at physical location `x` for time `t`
"""
initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D)
Expand Down
4 changes: 4 additions & 0 deletions src/equations/ideal_glm_mhd_multicomponent_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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],
amrueda marked this conversation as resolved.
Show resolved Hide resolved
u[7])

"""
initial_condition_convergence_test(x, t, equations::IdealGlmMhdMulticomponentEquations2D)

Expand Down
Loading