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

Jump to bottom

WIP: Unify divb analyze routines for equations with magnetic field #1712

Draft
wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

WIP: Unify divb analyze routines for equations with magnetic field #1712

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5d35dc3
Unified divb analyze routines for equations with `magnetic_field`
amrueda Nov 7, 2023
d6ad2f0
format
amrueda Nov 7, 2023
9e8dd98
format
amrueda Nov 7, 2023
34dc448
Changed format of `magnetic_field` for consistency
amrueda Nov 7, 2023
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: 20 additions & 50 deletions src/callbacks_step/analysis_dg2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Comment on lines +255 to +256
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
B1_kj, _, _ = magnetic_field(view(u, :, k, j, element), equations)
_, B2_ik, _ = magnetic_field(view(u, :, i, k, element), equations)
B1_kj, _, _ = magnetic_field(get_node_vars(u, equations, dg, k, j, element), equations)
_, B2_ik, _ = magnetic_field(get_node_vars(u, equations, dg, i, k, element), equations)

To avoid the ugly view? If this works, it should also be applied to the other places in this PR

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It works for a short simulation, but it seems that it creates a lot of allocations. I tested with the following:

julia> using Trixi
julia> using BenchmarkTools
julia> include("../examples/tree_2d_dgsem/elixir_mhd_ec.jl")
julia> @benchmark Trixi.analyze($Val(:l2_divb), $ode.u0, $ode.u0, 0.0, $mesh, $equations, $solver, $semi.cache)

I got the following for the implementation with view on Rocinante:

BenchmarkTools.Trial: 4272 samples with 1 evaluation.
 Range (min … max):  661.995 μs … 9.138 ms  ┊ GC (min … max):  0.00% … 88.29%
 Time  (median):     898.936 μs             ┊ GC (median):     0.00%
 Time  (mean ± σ):     1.164 ms ± 1.153 ms  ┊ GC (mean ± σ):  20.53% ± 17.92%

  █▆▆▆▅▅▃▁                                                    ▁
  █████████▅▃▅▅▃▁▁▁▁▁▆▃▇▁▁▁▁▁▆▇▇▇▆▁▁▄▄▃▃▃▃▄▅▁▄▄▃▃▄▄▅▅▅▇▄▆▆▆▆▆ █
  662 μs       Histogram: log(frequency) by time      7.53 ms <

 Memory estimate: 2.13 MiB, allocs estimate: 28685.

But the implementation with get_node_vars crashes with the following (I removed many lines in the middle of the error message for readability):

[644396] signal (11.1): Segmentation fault
in expression starting at REPL[4]:1
getindex at ./essentials.jl:14 [inlined]
#309 at /mnt/hd1/scratch/aruedara/20231107_divB/Trixi.jl/src/solvers/dg.jl:542 [inlined]
macro expansion at ./ntuple.jl:72 [inlined]
ntuple at ./ntuple.jl:69 [inlined]
get_node_vars at /mnt/hd1/scratch/aruedara/20231107_divB/Trixi.jl/src/solvers/dg.jl:542 [inlined]
.
.
.
jl_apply at /cache/build/default-amdci5-5/julialang/julia-release-1-dot-9/src/julia.h:1880 [inlined]
true_main at /cache/build/default-amdci5-5/julialang/julia-release-1-dot-9/src/jlapi.c:573
jl_repl_entrypoint at /cache/build/default-amdci5-5/julialang/julia-release-1-dot-9/src/jlapi.c:717
main at julia (unknown line)
__libc_start_main at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
unknown function (ip: 0x4010b8)
Allocations: 83483411 (Pool: 83454153; Big: 29258); GC: 106
Segmentation fault (core dumped)

Is this test OK? should I keep the implementation with view?

Copy link
Member

Choose a reason for hiding this comment

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

No, I think there seems to be an underlying issue that I just do not understand.

Can you try to separate the get_node_vars call from the magnetic_field call? That is, to use something like

Suggested change
B1_kj, _, _ = magnetic_field(view(u, :, k, j, element), equations)
_, B2_ik, _ = magnetic_field(view(u, :, i, k, element), equations)
u_local_kj = get_node_vars(u, equations, dg, k, j, element)
u_local_ik = get_node_vars(u, equations, dg, i, k, element)
B1_kj, _, _ = magnetic_field(u_local_kj, equations)
_, B2_ik, _ = magnetic_field(u_local_ik, equations)

such that we know where exactly the issue comes from?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That works 🙂 :

BenchmarkTools.Trial: 3875 samples with 1 evaluation.
 Range (min … max):  685.359 μs … 9.795 ms  ┊ GC (min … max):  0.00% … 83.22%
 Time  (median):     992.257 μs             ┊ GC (median):     0.00%
 Time  (mean ± σ):     1.284 ms ± 1.303 ms  ┊ GC (mean ± σ):  20.88% ± 17.77%

  █▅▆▆▆▅▄▂▁                                                   ▁
  █████████▆▆▃▃▁▁▁▁▁▁▁▆▅▆▁▁▁▁▁▄▆▇▅▅▃▃▃▅▃▆▅▅▅▁▁▁▃▅▃▃▅▄▇▆▇▆▆▆▄▆ █
  685 μs       Histogram: log(frequency) by time       8.2 ms <

 Memory estimate: 2.13 MiB, allocs estimate: 28685.

Is this preferred to the view implementation?

Copy link
Member

Choose a reason for hiding this comment

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

Is this preferred to the view implementation?

Yes. We introduced the get_node_vars function and friends specifically for the purpose of not having to use views (which used to allocate) and also to hide memory layout details from the application.

@ranocha out of curiosity, do you have an explanation why my original formulation caused a segfault?

Copy link
Member

Choose a reason for hiding this comment

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

You called it with u_ode as argument but the function is written to accept u? But I would have expected that both versions lead to the same problem

Copy link
Contributor Author

Choose a reason for hiding this comment

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

You called it with u_ode as argument but the function is written to accept u?

No, I called it exactly as detailed above, with ode.u0 for both u and du (du is not used in the function).

I'm getting a weird behavior... Yesterday, when I tried the first code suggestion of Michael, I got the error message twice in a row. Today, I decided to give it a go again, and the benchmark is working fine (now 12 times in a row). Maybe the problem yesterday was that someone else was running a simulation on Rocinante(?).

Copy link
Member

Choose a reason for hiding this comment

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

That's what I'm saying - ode.u0 is a vector, we accept the wrapped array

Copy link
Member

Choose a reason for hiding this comment

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

Do you run it with any --check-bounds flags?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ahh, ok! I didn't realize that... No, I ran it with --check-bounds=no. Without that flag it crashes with

BoundsError: attempt to access 36864-element Vector{Float64} at index [1, 2, 1, 1]

I corrected my testing procedure and now I'm using:

julia> using Trixi
julia> using BenchmarkTools
julia> include("../examples/tree_2d_dgsem/elixir_mhd_ec.jl")
julia> u = Trixi.wrap_array(ode.u0, mesh, equations, solver, semi.cache)
julia> @benchmark Trixi.analyze($Val(:l2_divb), $u, $u, 0.0, $mesh, $equations, $solver, $semi.cache)

This seems to work fine using the first suggestion by @sloede (checking bounds):

BenchmarkTools.Trial: 5838 samples with 1 evaluation.
 Range (min … max):  700.557 μs …   3.900 ms  ┊ GC (min … max):  0.00% … 69.72%
 Time  (median):     724.738 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   853.904 μs ± 585.713 μs  ┊ GC (mean ± σ):  14.73% ± 16.32%

  █▂                                                          ▃ ▁
  ██▃▃▁▁▁▁▆▄▅▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇█ █
  701 μs        Histogram: log(frequency) by time       3.55 ms <

 Memory estimate: 2.31 MiB, allocs estimate: 32785.

For reference, this is the original implementation in Trixi:

BenchmarkTools.Trial: 5729 samples with 1 evaluation.
 Range (min … max):  711.988 μs …   3.850 ms  ┊ GC (min … max):  0.00% … 71.50%
 Time  (median):     738.548 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   870.079 μs ± 592.862 μs  ┊ GC (mean ± σ):  14.63% ± 16.28%

  █▃                                                          ▃ ▁
  ██▄▃▁▁▆▆▄▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄█ █
  712 μs        Histogram: log(frequency) by time        3.6 ms <

 Memory estimate: 2.31 MiB, allocs estimate: 32785.

So it seems that it works fine with get_node_vars.

divb += (derivative_matrix[i, k] * B1_kj +
derivative_matrix[j, k] * B2_ik)
end
divb *= cache.elements.inverse_jacobian[element]
divb^2
Expand All @@ -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,
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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))
Expand All @@ -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

Expand All @@ -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]) +
(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
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 @@ -1338,4 +1338,7 @@ 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
5 changes: 5 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 @@ -711,4 +711,9 @@ end
return SVector{ncomponents(equations), real(equations)}(u[i + 8] * v
for i in eachcomponent(equations))
end

# Return the magnetic field
function magnetic_field(u, equations::IdealGlmMhdMulticomponentEquations2D)
SVector(u[5], u[6], u[7])
end
amrueda marked this conversation as resolved.
Show resolved Hide resolved
end # @muladd
Loading