From cca0674886a922971db9ae06bb019a0e219eae7f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 29 Oct 2024 10:52:14 +0100 Subject: [PATCH] fix P4estMesh and T8codeMesh plus tests --- src/solvers/dgsem_p4est/dg_2d.jl | 51 ++++++++++++++-------- src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 23 +++++----- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 6 +-- test/test_p4est_2d.jl | 24 +++++----- test/test_parabolic_2d.jl | 4 +- test/test_t8code_2d.jl | 13 +++--- 6 files changed, 67 insertions(+), 54 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 3c868289181..acbfdbb491d 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -13,11 +13,15 @@ function create_cache(mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, uEltype, 2, nvariables(equations) * nnodes(mortar_l2)} - fstar_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] - fstar_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_primary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_primary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_secondary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_secondary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] u_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] - (; fstar_upper_threaded, fstar_lower_threaded, u_threaded) + (; fstar_primary_upper_threaded, fstar_primary_lower_threaded, + fstar_secondary_upper_threaded, fstar_secondary_lower_threaded, + u_threaded) end # index_to_start_step_2d(index::Symbol, index_range) @@ -451,13 +455,16 @@ function calc_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.mortars @unpack contravariant_vectors = cache.elements - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded, fstar_secondary_upper_threaded, fstar_secondary_lower_threaded = cache index_range = eachnode(dg) @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar = (fstar_lower_threaded[Threads.threadid()], - fstar_upper_threaded[Threads.threadid()]) + fstar_primary = (fstar_primary_lower_threaded[Threads.threadid()], + fstar_primary_upper_threaded[Threads.threadid()]) + + fstar_secondary = (fstar_secondary_lower_threaded[Threads.threadid()], + fstar_secondary_upper_threaded[Threads.threadid()]) # Get index information on the small elements small_indices = node_indices[1, mortar] @@ -480,7 +487,8 @@ function calc_mortar_flux!(surface_flux_values, contravariant_vectors, i_small, j_small, element) - calc_mortar_flux!(fstar, mesh, nonconservative_terms, equations, + calc_mortar_flux!(fstar_primary, fstar_secondary, + mesh, nonconservative_terms, equations, surface_integral, dg, cache, mortar, position, normal_direction, node) @@ -501,14 +509,15 @@ function calc_mortar_flux!(surface_flux_values, # "mortar_fluxes_to_elements!" instead. mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, - mortar, fstar, u_buffer) + mortar, fstar_primary, fstar_secondary, + u_buffer) end return nothing end # Inlined version of the mortar flux computation on small elements for conservation laws -@inline function calc_mortar_flux!(fstar, +@inline function calc_mortar_flux!(fstar_primary, fstar_secondary, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, @@ -523,12 +532,13 @@ end flux = surface_flux(u_ll, u_rr, normal_direction, equations) # Copy flux to buffer - set_node_vars!(fstar[position_index], flux, equations, dg, node_index) + set_node_vars!(fstar_primary[position_index], flux, equations, dg, node_index) + set_node_vars!(fstar_secondary[position_index], flux, equations, dg, node_index) end # Inlined version of the mortar flux computation on small elements for equations with conservative and # nonconservative terms -@inline function calc_mortar_flux!(fstar, +@inline function calc_mortar_flux!(fstar_primary, fstar_secondary, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, @@ -547,19 +557,23 @@ end # The nonconservative flux is scaled by a factor of 0.5 based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - noncons = nonconservative_flux(u_ll, u_rr, normal_direction, equations) + noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction, equations) + noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, equations) - flux_plus_noncons = flux + 0.5f0 * noncons + flux_plus_noncons_primary = flux + 0.5f0 * noncons_primary + flux_plus_noncons_secondary = flux + 0.5f0 * noncons_secondary # Copy to buffer - set_node_vars!(fstar[position_index], flux_plus_noncons, equations, dg, node_index) + set_node_vars!(fstar_primary[position_index], flux_plus_noncons_primary, equations, dg, node_index) + set_node_vars!(fstar_secondary[position_index], flux_plus_noncons_secondary, equations, dg, node_index) end @inline function mortar_fluxes_to_elements!(surface_flux_values, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, mortar_l2::LobattoLegendreMortarL2, - dg::DGSEM, cache, mortar, fstar, u_buffer) + dg::DGSEM, cache, mortar, fstar_primary, + fstar_secondary, u_buffer) @unpack neighbor_ids, node_indices = cache.mortars # Copy solution small to small @@ -570,16 +584,15 @@ end element = neighbor_ids[position, mortar] for i in eachnode(dg) for v in eachvariable(equations) - surface_flux_values[v, i, small_direction, element] = fstar[position][v, - i] + surface_flux_values[v, i, small_direction, element] = fstar_primary[position][v, i] end end end # Project small fluxes to large element. multiply_dimensionwise!(u_buffer, - mortar_l2.reverse_upper, fstar[2], - mortar_l2.reverse_lower, fstar[1]) + mortar_l2.reverse_upper, fstar_secondary[2], + mortar_l2.reverse_lower, fstar_secondary[1]) # The flux is calculated in the outward direction of the small elements, # so the sign must be switched to get the flux in outward direction diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index 23968a2bfac..88b712d1b9b 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -398,7 +398,8 @@ end mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations::AbstractEquationsParabolic, mortar_l2::LobattoLegendreMortarL2, - dg::DGSEM, cache, mortar, fstar, u_buffer) + dg::DGSEM, cache, mortar, fstar_primary, + fstar_secondary, u_buffer) @unpack neighbor_ids, node_indices = cache.mortars # Copy solution small to small small_indices = node_indices[1, mortar] @@ -408,16 +409,15 @@ end element = neighbor_ids[position, mortar] for i in eachnode(dg) for v in eachvariable(equations) - surface_flux_values[v, i, small_direction, element] = fstar[position][v, - i] + surface_flux_values[v, i, small_direction, element] = fstar_primary[position][v, i] end end end # Project small fluxes to large element. multiply_dimensionwise!(u_buffer, - mortar_l2.reverse_upper, fstar[2], - mortar_l2.reverse_lower, fstar[1]) + mortar_l2.reverse_upper, fstar_secondary[2], + mortar_l2.reverse_lower, fstar_secondary[1]) # Copy interpolated flux values from buffer to large element face in the # correct orientation. @@ -772,13 +772,13 @@ function calc_mortar_flux_divergence!(surface_flux_values, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.mortars @unpack contravariant_vectors = cache.elements - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded = cache index_range = eachnode(dg) @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar = (fstar_lower_threaded[Threads.threadid()], - fstar_upper_threaded[Threads.threadid()]) + fstar = (fstar_primary_lower_threaded[Threads.threadid()], + fstar_primary_upper_threaded[Threads.threadid()]) for position in 1:2 for node in eachnode(dg) @@ -802,7 +802,7 @@ function calc_mortar_flux_divergence!(surface_flux_values, # this reuses the hyperbolic version of `mortar_fluxes_to_elements!` mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, - mortar, fstar, u_buffer) + mortar, fstar, fstar, u_buffer) end return nothing @@ -813,7 +813,7 @@ end # The reasoning is that parabolic fluxes are treated like conservative # terms (e.g., we compute a viscous conservative "flux") and thus no # non-conservative terms are present. -@inline function calc_mortar_flux!(fstar, +@inline function calc_mortar_flux!(fstar_primary, fstar_secondary, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations::AbstractEquationsParabolic, @@ -830,7 +830,8 @@ end flux_ = 0.5f0 * (u_ll + u_rr) # Copy flux to buffer - set_node_vars!(fstar[position_index], flux_, equations, dg, node_index) + set_node_vars!(fstar_primary[position_index], flux_, equations, dg, node_index) + set_node_vars!(fstar_secondary[position_index], flux_, equations, dg, node_index) end # TODO: parabolic, finish implementing `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!` diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 3bf0cd0cab5..633c6f1d067 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -209,13 +209,13 @@ function calc_mpi_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) @unpack local_neighbor_ids, local_neighbor_positions, node_indices = cache.mpi_mortars @unpack contravariant_vectors = cache.elements - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded = cache index_range = eachnode(dg) @threaded for mortar in eachmpimortar(dg, cache) # Choose thread-specific pre-allocated container - fstar = (fstar_lower_threaded[Threads.threadid()], - fstar_upper_threaded[Threads.threadid()]) + fstar = (fstar_primary_lower_threaded[Threads.threadid()], + fstar_primary_upper_threaded[Threads.threadid()]) # Get index information on the small elements small_indices = node_indices[1, mortar] diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 46dfa542f34..7a1a3f3ac09 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -296,8 +296,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 @@ -514,8 +514,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 @@ -617,18 +617,18 @@ end @trixi_testset "elixir_mhd_rotor.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.45520942119628094, 0.8918052934746296, 0.8324715234698744, + l2=[0.4551839744017604, 0.8917986079085971, 0.832474072904728, 0.0, - 0.9801268321975156, 0.1047572273917542, - 0.15551326369065485, + 0.98015167453983, 0.10475978783943254, + 0.15551175906375883, 0.0, - 2.0604823850230503e-5], - linf=[10.194219681102396, 18.254409373691253, - 10.031954811617876, + 2.026208477271868e-5], + linf=[10.19496728149964, 18.23726813972206, + 10.04367783820621, 0.0, - 19.646870952133547, 1.39386796733875, 1.8725058402095736, + 19.63022306543678, 1.3952679820406384, 1.8716515525771589, 0.0, - 0.001619787048808356], + 0.0017266639582675424], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index ceefb65e99b..2585c350dbc 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -579,8 +579,8 @@ end @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_nonperiodic_amr.jl"), tspan=(0.0, 0.01), - l2=[0.007934195641974433], - linf=[0.11030265194954081]) + l2=[0.00793370348206051], + linf=[0.11029310107224988]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index ed9edbee3df..4419fecc183 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -305,17 +305,16 @@ end # This test is identical to the one in `test_p4est_2d.jl` besides minor # deviations in the expected error norms. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.4420732463420727, 0.8804644301158163, 0.8262542320734158, + l2=[0.4419337424073218, 0.8804938551016932, 0.8258185723720365, 0.0, - 0.9615023124248694, 0.10386709616933161, - 0.15403081916109138, + 0.961220188718187, 0.10397273631386837, 0.15408979488125943, 0.0, - 2.835066224683485e-5], - linf=[10.045486750338348, 17.998696851793447, 9.57580213608948, + 2.66769410449947e-5], + linf=[10.053140536236942, 18.17070117006211, 9.549208389448738, 0.0, - 19.431290734386764, 1.3821685025605288, 1.8186235976086789, + 19.676151923191583, 1.3896544044814965, 1.8153256887969416, 0.0, - 0.0023118793481168537], + 0.0022030404596184786], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities)