From f4c731001c8cada228f473f8f98df00d087f3450 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Thu, 26 Mar 2026 16:01:17 +0100 Subject: [PATCH 01/17] enabled MPI for 2D parabolic system on P4est mesh (from present main) --- .../analysis_surface_integral_2d.jl | 17 +- .../dgsem_p4est/dg_2d_parabolic_parallel.jl | 501 ++++++++++++++++++ src/solvers/dgsem_p4est/dg_parallel.jl | 1 + test/test_mpi.jl | 1 + test/test_mpi_p4est_parabolic_2d.jl | 126 +++++ 5 files changed, 644 insertions(+), 2 deletions(-) create mode 100644 src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl create mode 100644 test/test_mpi_p4est_parabolic_2d.jl diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index fef4b9872d1..4d99c5325c3 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -250,7 +250,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, j_node += j_node_step end end - return surface_integral + return distribute_surface_integral(surface_integral, mesh) end # 2D version of the `analyze` function for `AnalysisSurfaceIntegral` of viscous, i.e., @@ -320,6 +320,19 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, du, u, t, j_node += j_node_step end end - return surface_integral + return distribute_surface_integral(surface_integral, mesh) +end + +# Serial/default: do nothing +distribute_surface_integral(val, mesh) = val + +# Parallel: sum over all ranks +function distribute_surface_integral(val, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}) + comm = MPI.COMM_WORLD + buf = [val] + MPI.Allreduce!(buf, MPI.SUM, comm) + return buf[1] end end # muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl new file mode 100644 index 00000000000..6006df085ae --- /dev/null +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -0,0 +1,501 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin + function rhs_parabolic!(du, u, t, + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + equations_parabolic::AbstractEquationsParabolic, + boundary_conditions_parabolic, source_terms_parabolic, + dg::DG, parabolic_scheme, cache, cache_parabolic) + @unpack viscous_container = cache_parabolic + @unpack u_transformed, gradients, flux_viscous = viscous_container + + # Stage 0: local variable transform + # + @trixi_timeit timer() "transform variables" begin + transform_variables!(u_transformed, u, mesh, equations_parabolic, + dg, cache) + end + + # + # Stage 1: gradient computation + # + + # Start gradient-stage MPI receive + @trixi_timeit timer() "start MPI receive gradient" begin + start_mpi_receive!(cache.mpi_cache) + end + + # Prolong transformed variables to MPI mortars + # @trixi_timeit timer() "prolong2mpimortars gradient" begin + # prolong2mpimortars!(cache, u_transformed, mesh, equations_parabolic, + # dg.mortar, dg) + # end + + # Prolong transformed variables to MPI interfaces + @trixi_timeit timer() "prolong2mpiinterfaces gradient" begin + prolong2mpiinterfaces!(cache, u_transformed, mesh, equations_parabolic, + dg.surface_integral, dg) + end + + # Start gradient-stage MPI send + @trixi_timeit timer() "start MPI send gradient" begin + start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + end + + # Local gradient computation + @trixi_timeit timer() "calculate gradient local" begin + calc_gradient_local!(gradients, u_transformed, t, mesh, + equations_parabolic, boundary_conditions_parabolic, + dg, parabolic_scheme, cache) + end + + # Finish gradient-stage MPI receive + @trixi_timeit timer() "finish MPI receive gradient" begin + finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + end + # Finish gradient-stage MPI send + @trixi_timeit timer() "finish MPI send gradient" begin + finish_mpi_send!(cache.mpi_cache) + end + # MPI interface fluxes for gradient stage + @trixi_timeit timer() "MPI interface flux gradient" begin + calc_mpi_interface_flux_gradient!(cache.elements.surface_flux_values, + mesh, equations_parabolic, + dg, parabolic_scheme, cache) + end + + # MPI mortar fluxes for gradient stage + # @trixi_timeit timer() "MPI mortar flux gradient" begin + # calc_mpi_mortar_flux_gradient!(cache.elements.surface_flux_values, + # mesh, equations_parabolic, dg.mortar, + # dg, parabolic_scheme, cache) + # end + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" begin + calc_surface_integral_gradient!(gradients, mesh, equations_parabolic, + dg, cache) + end + + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" begin + apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg, + cache) + end + # Stage 2: local viscous flux construction + # + @trixi_timeit timer() "calculate viscous fluxes" begin + calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh, + equations_parabolic, dg, cache) + end + + # + # Stage 3: divergence of viscous fluxes + # + + # Start divergence-stage MPI receive + @trixi_timeit timer() "start MPI receive divergence" begin + start_mpi_receive!(cache.mpi_cache) + end + + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" begin + set_zero!(du, dg, cache) + end + + # Local volume integral + @trixi_timeit timer() "volume integral" begin + calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache) + end + + # Prolong viscous fluxes to MPI mortars + # @trixi_timeit timer() "prolong2mpimortars divergence" begin + # prolong2mpimortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, + # dg.mortar, dg) + # end + + # Prolong viscous fluxes to MPI interfaces + @trixi_timeit timer() "prolong2mpiinterfaces divergence" begin + prolong2mpiinterfaces!(cache, flux_viscous, mesh, equations_parabolic, dg) + end + ########################## Divergence ################################# + # Start divergence-stage MPI send + @trixi_timeit timer() "start MPI send divergence" begin + start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + end + + # Local interface fluxes + @trixi_timeit timer() "prolong2interfaces" begin + prolong2interfaces!(cache, flux_viscous, mesh, equations_parabolic, dg) + end + + @trixi_timeit timer() "interface flux" begin + calc_interface_flux!(cache.elements.surface_flux_values, mesh, + equations_parabolic, dg, parabolic_scheme, cache) + end + + # Local boundary fluxes + @trixi_timeit timer() "prolong2boundaries" begin + prolong2boundaries!(cache, flux_viscous, mesh, equations_parabolic, dg) + end + + @trixi_timeit timer() "boundary flux" begin + calc_boundary_flux_divergence!(cache, t, + boundary_conditions_parabolic, mesh, + equations_parabolic, + dg.surface_integral, dg) + end + + # Local mortar fluxes + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, + dg.mortar, dg) + end + + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux_divergence!(cache.elements.surface_flux_values, + mesh, equations_parabolic, dg.mortar, + dg, parabolic_scheme, cache) + end + + # Finish divergence-stage MPI receive + @trixi_timeit timer() "finish MPI receive divergence" begin + finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + end + + # MPI interface fluxes for divergence stage + @trixi_timeit timer() "MPI interface flux divergence" begin + calc_mpi_interface_flux_divergence!(cache.elements.surface_flux_values, + mesh, equations_parabolic, + dg, parabolic_scheme, cache) + end + + # MPI mortar fluxes for divergence stage + # @trixi_timeit timer() "MPI mortar flux divergence" begin + # calc_mpi_mortar_flux_divergence!(cache.elements.surface_flux_values, + # mesh, equations_parabolic, dg.mortar, + # dg, parabolic_scheme, cache) + # end + + # Finish divergence-stage MPI send + @trixi_timeit timer() "finish MPI send divergence" begin + finish_mpi_send!(cache.mpi_cache) + end + + # + # Stage 4: final assembly + # + @trixi_timeit timer() "surface integral" begin + calc_surface_integral!(du, u, mesh, equations_parabolic, + dg.surface_integral, dg, cache) + end + + @trixi_timeit timer() "Jacobian" begin + apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache) + end + + @trixi_timeit timer() "source terms parabolic" begin + calc_sources_parabolic!(du, u, gradients, t, source_terms_parabolic, + equations_parabolic, dg, cache) + end + + return nothing + end + + function calc_gradient_local!(gradients, u_transformed, t, + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + equations_parabolic, boundary_conditions_parabolic, + dg::DG, parabolic_scheme, cache) + # Reset gradients + @trixi_timeit timer() "reset gradients" begin + reset_gradients!(gradients, dg, cache) + end + + # Calculate volume integral + @trixi_timeit timer() "volume integral" begin + calc_volume_integral_gradient!(gradients, u_transformed, + mesh, equations_parabolic, dg, cache) + end + + # Prolong solution to interfaces. + # This reuses `prolong2interfaces` for the purely hyperbolic case. + @trixi_timeit timer() "prolong2interfaces" begin + prolong2interfaces!(cache, u_transformed, mesh, + equations_parabolic, dg) + end + + # Calculate interface fluxes for the gradients + @trixi_timeit timer() "interface flux" begin + @unpack surface_flux_values = cache.elements + calc_interface_flux_gradient!(surface_flux_values, mesh, equations_parabolic, + dg, parabolic_scheme, cache) + end + + # Prolong solution to boundaries. + # This reuses `prolong2boundaries` for the purely hyperbolic case. + @trixi_timeit timer() "prolong2boundaries" begin + prolong2boundaries!(cache, u_transformed, mesh, + equations_parabolic, dg) + end + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" begin + calc_boundary_flux_gradient!(cache, t, boundary_conditions_parabolic, + mesh, equations_parabolic, dg.surface_integral, + dg) + end + + # Prolong solution to mortars. + # This reuses `prolong2mortars` for the purely hyperbolic case. + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars!(cache, u_transformed, mesh, equations_parabolic, + dg.mortar, dg) + end + + # Calculate mortar fluxes + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux_gradient!(cache.elements.surface_flux_values, + mesh, equations_parabolic, dg.mortar, + dg, parabolic_scheme, cache) + end + + return nothing + end + + function prolong2mpiinterfaces!(cache, flux_viscous::Tuple, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, dg::DG) + @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + + flux_viscous_x, flux_viscous_y = flux_viscous + + @threaded for interface in eachmpiinterface(dg, cache) + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + local_direction = indices2direction(local_indices) + local_side = local_sides[interface] + orientationFactor = local_side == 1 ? 1 : -1 + + i_start, i_step = index_to_start_step_2d(local_indices[1], index_range) + j_start, j_step = index_to_start_step_2d(local_indices[2], index_range) + + i_elem = i_start + j_elem = j_start + + for i in eachnode(dg) + normal_direction = get_normal_direction(local_direction, + contravariant_vectors, + i_elem, j_elem, + local_element) + + for v in eachvariable(equations_parabolic) + flux_visc = SVector(flux_viscous_x[v, i_elem, j_elem, local_element], + flux_viscous_y[v, i_elem, j_elem, local_element]) + + cache.mpi_interfaces.u[local_side, v, i, interface] = orientationFactor * + dot(flux_visc, + normal_direction) + end + + i_elem += i_step + j_elem += j_step + end + end + + return nothing + end + + function calc_mpi_interface_flux_gradient!(surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, + dg::DG, parabolic_scheme, cache) + @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + index_end = last(index_range) + + @threaded for interface in eachmpiinterface(dg, cache) + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + local_direction = indices2direction(local_indices) + local_side = local_sides[interface] + + # Create the local i,j indexing on the local element used to pull normal direction information + i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], + index_range) + j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], + index_range) + + i_element = i_element_start + j_element = j_element_start + + # Initiate the node index to be used in the surface for loop, + # the surface flux storage must be indexed in alignment with the local element indexing + if :i_backward in local_indices + surface_node = index_end + surface_node_step = -1 + else + surface_node = 1 + surface_node_step = 1 + end + + for i in eachnode(dg) + normal_direction = get_normal_direction(local_direction, + contravariant_vectors, + i_element, j_element, + local_element) + + calc_mpi_interface_flux_gradient!(surface_flux_values, mesh, + equations_parabolic, + dg, parabolic_scheme, cache, + interface, normal_direction, + i, local_side, + surface_node, + local_direction, local_element) + + # Increment local element indices to pull the normal direction + i_element += i_element_step + j_element += j_element_step + + # Increment the surface node index along the local element + surface_node += surface_node_step + end + end + + return nothing + end + + @inline function calc_mpi_interface_flux_gradient!(surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, + dg::DG, parabolic_scheme, cache, + interface_index, normal_direction, + interface_node_index, local_side, + surface_node_index, + local_direction_index, + local_element_index) + @unpack u = cache.mpi_interfaces + + u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, + interface_node_index, + interface_index) + + flux_ = flux_parabolic(u_ll, u_rr, normal_direction, Gradient(), + equations_parabolic, parabolic_scheme) + + for v in eachvariable(equations_parabolic) + surface_flux_values[v, surface_node_index, + local_direction_index, local_element_index] = flux_[v] + end + + return nothing + end + + function calc_mpi_interface_flux_divergence!(surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, + dg::DG, parabolic_scheme, cache) + @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + index_end = last(index_range) + + @threaded for interface in eachmpiinterface(dg, cache) + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + local_direction = indices2direction(local_indices) + local_side = local_sides[interface] + + i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], + index_range) + j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], + index_range) + + i_element = i_element_start + j_element = j_element_start + + # Initiate the node index to be used in the surface for loop, + # the surface flux storage must be indexed in alignment with the local element indexing + if :i_backward in local_indices + surface_node = index_end + surface_node_step = -1 + else + surface_node = 1 + surface_node_step = 1 + end + + for i in eachnode(dg) + normal_direction = get_normal_direction(local_direction, + contravariant_vectors, + i_element, j_element, + local_element) + + calc_mpi_interface_flux_divergence!(surface_flux_values, mesh, + equations_parabolic, + dg, parabolic_scheme, cache, + interface, normal_direction, + i, local_side, + surface_node, + local_direction, local_element) + + i_element += i_element_step + j_element += j_element_step + + surface_node += surface_node_step + end + end + + return nothing + end + + @inline function calc_mpi_interface_flux_divergence!(surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, + dg::DG, parabolic_scheme, cache, + interface_index, normal_direction, + interface_node_index, local_side, + surface_node_index, + local_direction_index, + local_element_index) + @unpack u = cache.mpi_interfaces + + viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(u, + equations_parabolic, + dg, + interface_node_index, + interface_index) + + flux_ = flux_parabolic(viscous_flux_normal_ll, viscous_flux_normal_rr, + normal_direction, Divergence(), + equations_parabolic, parabolic_scheme) + + if local_side == 1 + flux_ = flux_parabolic(viscous_flux_normal_ll, viscous_flux_normal_rr, + normal_direction, Divergence(), + equations_parabolic, parabolic_scheme) + sign_ = 1 + else + # side 2 must also use primary-oriented normal + flux_ = flux_parabolic(viscous_flux_normal_ll, viscous_flux_normal_rr, + -normal_direction, Divergence(), + equations_parabolic, parabolic_scheme) + sign_ = -1 + end + + for v in eachvariable(equations_parabolic) + surface_flux_values[v, surface_node_index, + local_direction_index, local_element_index] = sign_ * flux_[v] + end + + return nothing + end +end #muladd diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index e7f97d96933..aea544026ef 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -608,4 +608,5 @@ end include("dg_2d_parallel.jl") include("dg_3d_parallel.jl") +include("dg_2d_parabolic_parallel.jl") end # muladd diff --git a/test/test_mpi.jl b/test/test_mpi.jl index 887eb9e01f5..3640ef8637b 100644 --- a/test/test_mpi.jl +++ b/test/test_mpi.jl @@ -26,6 +26,7 @@ CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows() if !CI_ON_WINDOWS # see comment on `CI_ON_WINDOWS` above include("test_mpi_p4est_3d.jl") include("test_mpi_t8code_3d.jl") + include("test_mpi_p4est_parabolic_2d.jl") end end # MPI diff --git a/test/test_mpi_p4est_parabolic_2d.jl b/test/test_mpi_p4est_parabolic_2d.jl new file mode 100644 index 00000000000..5d61c44a981 --- /dev/null +++ b/test/test_mpi_p4est_parabolic_2d.jl @@ -0,0 +1,126 @@ +module TestExamplesMPIP4estMesh2DParabolic + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = joinpath(examples_dir(), "p4est_2d_dgsem") + +@testset "P4estMesh MPI 2D Parabolic" begin + @trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_navierstokes_lid_driven_cavity.jl"), + initial_refinement_level=2, tspan=(0.0, 0.5), + l2=[ + 0.00028716166408816073, + 0.08101204560401647, + 0.02099595625377768, + 0.05008149754143295 + ], + linf=[ + 0.014804500261322406, + 0.9513271652357098, + 0.7223919625994717, + 1.4846907331004786 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + end + + @trixi_testset "P4estMesh2D: elixir_navierstokes_convergence_nonperiodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_navierstokes_convergence_nonperiodic.jl"), + initial_refinement_level=1, tspan=(0.0, 0.2), + l2=[ + 0.0004036496258545996, + 0.0005869762480189079, + 0.0009148853742181908, + 0.0011984191532764543 + ], + linf=[ + 0.0024993634989209923, + 0.009487866203496731, + 0.004505829506103787, + 0.011634902753554499 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + end + + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_diffusion_nonperiodic_curved.jl"), + trees_per_dimension=(1, 1), initial_refinement_level=2, + tspan=(0.0, 0.5), + l2=[0.00919917034843865], + linf=[0.14186297438393505]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + end + + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_diffusion_periodic.jl"), + trees_per_dimension=(1, 1), initial_refinement_level=2, + tspan=(0.0, 0.5), + l2=[0.0023754695605828443], + linf=[0.008154128363741964]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + end + + @trixi_testset "elixir_navierstokes_NACA0012airfoil_mach08.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_navierstokes_NACA0012airfoil_mach08.jl"), + l2=[0.000186486564226516, + 0.0005076712323400374, + 0.00038074588984354107, + 0.002128177239782089], + linf=[0.5153387072802718, + 1.199362305026636, + 0.9077214424040279, + 5.666071182328691], tspan=(0.0, 0.001), + initial_refinement_level=0) + + u_ode = copy(sol.u[end]) + du_ode = zero(u_ode) # Just a placeholder in this case + + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + + drag_p = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache, semi) + lift_p = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache, semi) + + drag_f = Trixi.analyze(drag_coefficient_shear_force, du, u, tspan[2], mesh, + equations, equations_parabolic, solver, + semi.cache, semi, semi.cache_parabolic) + lift_f = Trixi.analyze(lift_coefficient_shear_force, du, u, tspan[2], mesh, + equations, equations_parabolic, solver, + semi.cache, semi, semi.cache_parabolic) + + @test isapprox(drag_p, 0.17963843913309516, atol = 1e-13) + @test isapprox(lift_p, 0.26462588007949367, atol = 1e-13) + + @test isapprox(drag_f, 1.5427441885921553, atol = 1e-13) + @test isapprox(lift_f, 0.005621910087395724, atol = 1e-13) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + # We move these tests here to avoid modifying values used + # to compute the drag/lift coefficients above. + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + end +end +end # module From 5f779170340df15cfd33123bfe7256eae1d8c485 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Thu, 26 Mar 2026 22:18:52 +0100 Subject: [PATCH 02/17] adapted naming: viscous -> parabolic --- .../dgsem_p4est/dg_2d_parabolic_parallel.jl | 58 +++++++++---------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index 6006df085ae..012e57299a3 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -8,8 +8,8 @@ equations_parabolic::AbstractEquationsParabolic, boundary_conditions_parabolic, source_terms_parabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) - @unpack viscous_container = cache_parabolic - @unpack u_transformed, gradients, flux_viscous = viscous_container + @unpack parabolic_container = cache_parabolic + @unpack u_transformed, gradients, flux_parabolic = parabolic_container # Stage 0: local variable transform # @@ -84,15 +84,15 @@ apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg, cache) end - # Stage 2: local viscous flux construction + # Stage 2: local parabolic flux construction # - @trixi_timeit timer() "calculate viscous fluxes" begin - calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh, - equations_parabolic, dg, cache) + @trixi_timeit timer() "calculate parabolic fluxes" begin + calc_parabolic_fluxes!(flux_parabolic, gradients, u_transformed, mesh, + equations_parabolic, dg, cache) end # - # Stage 3: divergence of viscous fluxes + # Stage 3: divergence of parabolic fluxes # # Start divergence-stage MPI receive @@ -107,18 +107,18 @@ # Local volume integral @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache) + calc_volume_integral!(du, flux_parabolic, mesh, equations_parabolic, dg, cache) end - # Prolong viscous fluxes to MPI mortars + # Prolong parabolic fluxes to MPI mortars # @trixi_timeit timer() "prolong2mpimortars divergence" begin - # prolong2mpimortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, + # prolong2mpimortars_divergence!(cache, flux_parabolic, mesh, equations_parabolic, # dg.mortar, dg) # end - # Prolong viscous fluxes to MPI interfaces + # Prolong parabolic fluxes to MPI interfaces @trixi_timeit timer() "prolong2mpiinterfaces divergence" begin - prolong2mpiinterfaces!(cache, flux_viscous, mesh, equations_parabolic, dg) + prolong2mpiinterfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) end ########################## Divergence ################################# # Start divergence-stage MPI send @@ -128,7 +128,7 @@ # Local interface fluxes @trixi_timeit timer() "prolong2interfaces" begin - prolong2interfaces!(cache, flux_viscous, mesh, equations_parabolic, dg) + prolong2interfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) end @trixi_timeit timer() "interface flux" begin @@ -138,7 +138,7 @@ # Local boundary fluxes @trixi_timeit timer() "prolong2boundaries" begin - prolong2boundaries!(cache, flux_viscous, mesh, equations_parabolic, dg) + prolong2boundaries!(cache, flux_parabolic, mesh, equations_parabolic, dg) end @trixi_timeit timer() "boundary flux" begin @@ -150,7 +150,7 @@ # Local mortar fluxes @trixi_timeit timer() "prolong2mortars" begin - prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, + prolong2mortars_divergence!(cache, flux_parabolic, mesh, equations_parabolic, dg.mortar, dg) end @@ -188,7 +188,7 @@ # Stage 4: final assembly # @trixi_timeit timer() "surface integral" begin - calc_surface_integral!(du, u, mesh, equations_parabolic, + calc_surface_integral!(nothing, du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache) end @@ -222,7 +222,7 @@ # Prolong solution to interfaces. # This reuses `prolong2interfaces` for the purely hyperbolic case. @trixi_timeit timer() "prolong2interfaces" begin - prolong2interfaces!(cache, u_transformed, mesh, + prolong2interfaces!(nothing, cache, u_transformed, mesh, equations_parabolic, dg) end @@ -264,7 +264,7 @@ return nothing end - function prolong2mpiinterfaces!(cache, flux_viscous::Tuple, + function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, equations_parabolic, dg::DG) @@ -272,7 +272,7 @@ @unpack contravariant_vectors = cache.elements index_range = eachnode(dg) - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for interface in eachmpiinterface(dg, cache) local_element = local_neighbor_ids[interface] @@ -294,8 +294,8 @@ local_element) for v in eachvariable(equations_parabolic) - flux_visc = SVector(flux_viscous_x[v, i_elem, j_elem, local_element], - flux_viscous_y[v, i_elem, j_elem, local_element]) + flux_visc = SVector(flux_parabolic_x[v, i_elem, j_elem, local_element], + flux_parabolic_y[v, i_elem, j_elem, local_element]) cache.mpi_interfaces.u[local_side, v, i, interface] = orientationFactor * dot(flux_visc, @@ -468,24 +468,24 @@ local_element_index) @unpack u = cache.mpi_interfaces - viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(u, - equations_parabolic, - dg, - interface_node_index, - interface_index) + parabolic_flux_normal_ll, parabolic_flux_normal_rr = get_surface_node_vars(u, + equations_parabolic, + dg, + interface_node_index, + interface_index) - flux_ = flux_parabolic(viscous_flux_normal_ll, viscous_flux_normal_rr, + flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, normal_direction, Divergence(), equations_parabolic, parabolic_scheme) if local_side == 1 - flux_ = flux_parabolic(viscous_flux_normal_ll, viscous_flux_normal_rr, + flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, normal_direction, Divergence(), equations_parabolic, parabolic_scheme) sign_ = 1 else # side 2 must also use primary-oriented normal - flux_ = flux_parabolic(viscous_flux_normal_ll, viscous_flux_normal_rr, + flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, -normal_direction, Divergence(), equations_parabolic, parabolic_scheme) sign_ = -1 From e8d34a86ed862a6aa5ead644eacd1aa0f429911c Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Thu, 2 Apr 2026 15:50:42 +0200 Subject: [PATCH 03/17] added assertion to ensure no mortars are present --- .../dgsem_p4est/dg_2d_parabolic_parallel.jl | 866 +++++++++--------- 1 file changed, 434 insertions(+), 432 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index 012e57299a3..ac432fd31ce 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -3,499 +3,501 @@ # we need to opt-in explicitly. # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin - function rhs_parabolic!(du, u, t, - mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - equations_parabolic::AbstractEquationsParabolic, - boundary_conditions_parabolic, source_terms_parabolic, - dg::DG, parabolic_scheme, cache, cache_parabolic) - @unpack parabolic_container = cache_parabolic - @unpack u_transformed, gradients, flux_parabolic = parabolic_container - - # Stage 0: local variable transform - # - @trixi_timeit timer() "transform variables" begin - transform_variables!(u_transformed, u, mesh, equations_parabolic, - dg, cache) - end - - # - # Stage 1: gradient computation - # - - # Start gradient-stage MPI receive - @trixi_timeit timer() "start MPI receive gradient" begin - start_mpi_receive!(cache.mpi_cache) - end - - # Prolong transformed variables to MPI mortars - # @trixi_timeit timer() "prolong2mpimortars gradient" begin - # prolong2mpimortars!(cache, u_transformed, mesh, equations_parabolic, - # dg.mortar, dg) - # end - - # Prolong transformed variables to MPI interfaces - @trixi_timeit timer() "prolong2mpiinterfaces gradient" begin - prolong2mpiinterfaces!(cache, u_transformed, mesh, equations_parabolic, - dg.surface_integral, dg) - end - - # Start gradient-stage MPI send - @trixi_timeit timer() "start MPI send gradient" begin - start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) - end +#! format: noindent +function rhs_parabolic!(du, u, t, + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + equations_parabolic::AbstractEquationsParabolic, + boundary_conditions_parabolic, source_terms_parabolic, + dg::DG, parabolic_scheme, cache, cache_parabolic) + @unpack parabolic_container = cache_parabolic + @unpack u_transformed, gradients, flux_parabolic = parabolic_container + + # Stage 0: local variable transform + # + @trixi_timeit timer() "transform variables" begin + transform_variables!(u_transformed, u, mesh, equations_parabolic, + dg, cache) + end - # Local gradient computation - @trixi_timeit timer() "calculate gradient local" begin - calc_gradient_local!(gradients, u_transformed, t, mesh, - equations_parabolic, boundary_conditions_parabolic, - dg, parabolic_scheme, cache) - end + # + # Stage 1: gradient computation + # - # Finish gradient-stage MPI receive - @trixi_timeit timer() "finish MPI receive gradient" begin - finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) - end - # Finish gradient-stage MPI send - @trixi_timeit timer() "finish MPI send gradient" begin - finish_mpi_send!(cache.mpi_cache) - end - # MPI interface fluxes for gradient stage - @trixi_timeit timer() "MPI interface flux gradient" begin - calc_mpi_interface_flux_gradient!(cache.elements.surface_flux_values, - mesh, equations_parabolic, - dg, parabolic_scheme, cache) - end + # Start gradient-stage MPI receive + @trixi_timeit timer() "start MPI receive gradient" begin + start_mpi_receive!(cache.mpi_cache) + end - # MPI mortar fluxes for gradient stage - # @trixi_timeit timer() "MPI mortar flux gradient" begin - # calc_mpi_mortar_flux_gradient!(cache.elements.surface_flux_values, - # mesh, equations_parabolic, dg.mortar, - # dg, parabolic_scheme, cache) - # end - - # Calculate surface integrals - @trixi_timeit timer() "surface integral" begin - calc_surface_integral_gradient!(gradients, mesh, equations_parabolic, - dg, cache) - end + @assert isempty(eachmortar(dg, cache)) "Nonconforming meshes are not yet supported on MPI parallel P4estMesh." + # Prolong transformed variables to MPI mortars + # @trixi_timeit timer() "prolong2mpimortars gradient" begin + # prolong2mpimortars!(cache, u_transformed, mesh, equations_parabolic, + # dg.mortar, dg) + # end + + # Prolong transformed variables to MPI interfaces + @trixi_timeit timer() "prolong2mpiinterfaces gradient" begin + prolong2mpiinterfaces!(cache, u_transformed, mesh, equations_parabolic, + dg.surface_integral, dg) + end - # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" begin - apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg, - cache) - end - # Stage 2: local parabolic flux construction - # - @trixi_timeit timer() "calculate parabolic fluxes" begin - calc_parabolic_fluxes!(flux_parabolic, gradients, u_transformed, mesh, - equations_parabolic, dg, cache) - end + # Start gradient-stage MPI send + @trixi_timeit timer() "start MPI send gradient" begin + start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + end - # - # Stage 3: divergence of parabolic fluxes - # + # Local gradient computation + @trixi_timeit timer() "calculate gradient local" begin + calc_gradient_local!(gradients, u_transformed, t, mesh, + equations_parabolic, boundary_conditions_parabolic, + dg, parabolic_scheme, cache) + end - # Start divergence-stage MPI receive - @trixi_timeit timer() "start MPI receive divergence" begin - start_mpi_receive!(cache.mpi_cache) - end + # Finish gradient-stage MPI receive + @trixi_timeit timer() "finish MPI receive gradient" begin + finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + end + # Finish gradient-stage MPI send + @trixi_timeit timer() "finish MPI send gradient" begin + finish_mpi_send!(cache.mpi_cache) + end + # MPI interface fluxes for gradient stage + @trixi_timeit timer() "MPI interface flux gradient" begin + calc_mpi_interface_flux_gradient!(cache.elements.surface_flux_values, + mesh, equations_parabolic, + dg, parabolic_scheme, cache) + end - # Reset du - @trixi_timeit timer() "reset ∂u/∂t" begin - set_zero!(du, dg, cache) - end + # MPI mortar fluxes for gradient stage + # @trixi_timeit timer() "MPI mortar flux gradient" begin + # calc_mpi_mortar_flux_gradient!(cache.elements.surface_flux_values, + # mesh, equations_parabolic, dg.mortar, + # dg, parabolic_scheme, cache) + # end + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" begin + calc_surface_integral_gradient!(gradients, mesh, equations_parabolic, + dg, cache) + end - # Local volume integral - @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, flux_parabolic, mesh, equations_parabolic, dg, cache) - end + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" begin + apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg, + cache) + end + # Stage 2: local parabolic flux construction + # + @trixi_timeit timer() "calculate parabolic fluxes" begin + calc_parabolic_fluxes!(flux_parabolic, gradients, u_transformed, mesh, + equations_parabolic, dg, cache) + end - # Prolong parabolic fluxes to MPI mortars - # @trixi_timeit timer() "prolong2mpimortars divergence" begin - # prolong2mpimortars_divergence!(cache, flux_parabolic, mesh, equations_parabolic, - # dg.mortar, dg) - # end + # + # Stage 3: divergence of parabolic fluxes + # - # Prolong parabolic fluxes to MPI interfaces - @trixi_timeit timer() "prolong2mpiinterfaces divergence" begin - prolong2mpiinterfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) - end - ########################## Divergence ################################# - # Start divergence-stage MPI send - @trixi_timeit timer() "start MPI send divergence" begin - start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) - end + # Start divergence-stage MPI receive + @trixi_timeit timer() "start MPI receive divergence" begin + start_mpi_receive!(cache.mpi_cache) + end - # Local interface fluxes - @trixi_timeit timer() "prolong2interfaces" begin - prolong2interfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) - end + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" begin + set_zero!(du, dg, cache) + end - @trixi_timeit timer() "interface flux" begin - calc_interface_flux!(cache.elements.surface_flux_values, mesh, - equations_parabolic, dg, parabolic_scheme, cache) - end + # Local volume integral + @trixi_timeit timer() "volume integral" begin + calc_volume_integral!(du, flux_parabolic, mesh, equations_parabolic, dg, cache) + end - # Local boundary fluxes - @trixi_timeit timer() "prolong2boundaries" begin - prolong2boundaries!(cache, flux_parabolic, mesh, equations_parabolic, dg) - end + # Prolong parabolic fluxes to MPI mortars + # @trixi_timeit timer() "prolong2mpimortars divergence" begin + # prolong2mpimortars_divergence!(cache, flux_parabolic, mesh, equations_parabolic, + # dg.mortar, dg) + # end - @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux_divergence!(cache, t, - boundary_conditions_parabolic, mesh, - equations_parabolic, - dg.surface_integral, dg) - end + # Prolong parabolic fluxes to MPI interfaces + @trixi_timeit timer() "prolong2mpiinterfaces divergence" begin + prolong2mpiinterfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) + end + ########################## Divergence ################################# + # Start divergence-stage MPI send + @trixi_timeit timer() "start MPI send divergence" begin + start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + end - # Local mortar fluxes - @trixi_timeit timer() "prolong2mortars" begin - prolong2mortars_divergence!(cache, flux_parabolic, mesh, equations_parabolic, - dg.mortar, dg) - end + # Local interface fluxes + @trixi_timeit timer() "prolong2interfaces" begin + prolong2interfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) + end - @trixi_timeit timer() "mortar flux" begin - calc_mortar_flux_divergence!(cache.elements.surface_flux_values, - mesh, equations_parabolic, dg.mortar, - dg, parabolic_scheme, cache) - end + @trixi_timeit timer() "interface flux" begin + calc_interface_flux!(cache.elements.surface_flux_values, mesh, + equations_parabolic, dg, parabolic_scheme, cache) + end - # Finish divergence-stage MPI receive - @trixi_timeit timer() "finish MPI receive divergence" begin - finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) - end + # Local boundary fluxes + @trixi_timeit timer() "prolong2boundaries" begin + prolong2boundaries!(cache, flux_parabolic, mesh, equations_parabolic, dg) + end - # MPI interface fluxes for divergence stage - @trixi_timeit timer() "MPI interface flux divergence" begin - calc_mpi_interface_flux_divergence!(cache.elements.surface_flux_values, - mesh, equations_parabolic, - dg, parabolic_scheme, cache) - end + @trixi_timeit timer() "boundary flux" begin + calc_boundary_flux_divergence!(cache, t, + boundary_conditions_parabolic, mesh, + equations_parabolic, + dg.surface_integral, dg) + end - # MPI mortar fluxes for divergence stage - # @trixi_timeit timer() "MPI mortar flux divergence" begin - # calc_mpi_mortar_flux_divergence!(cache.elements.surface_flux_values, - # mesh, equations_parabolic, dg.mortar, - # dg, parabolic_scheme, cache) - # end + # Local mortar fluxes + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars_divergence!(cache, flux_parabolic, mesh, equations_parabolic, + dg.mortar, dg) + end - # Finish divergence-stage MPI send - @trixi_timeit timer() "finish MPI send divergence" begin - finish_mpi_send!(cache.mpi_cache) - end + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux_divergence!(cache.elements.surface_flux_values, + mesh, equations_parabolic, dg.mortar, + dg, parabolic_scheme, cache) + end - # - # Stage 4: final assembly - # - @trixi_timeit timer() "surface integral" begin - calc_surface_integral!(nothing, du, u, mesh, equations_parabolic, - dg.surface_integral, dg, cache) - end + # Finish divergence-stage MPI receive + @trixi_timeit timer() "finish MPI receive divergence" begin + finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + end - @trixi_timeit timer() "Jacobian" begin - apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache) - end + # MPI interface fluxes for divergence stage + @trixi_timeit timer() "MPI interface flux divergence" begin + calc_mpi_interface_flux_divergence!(cache.elements.surface_flux_values, + mesh, equations_parabolic, + dg, parabolic_scheme, cache) + end - @trixi_timeit timer() "source terms parabolic" begin - calc_sources_parabolic!(du, u, gradients, t, source_terms_parabolic, - equations_parabolic, dg, cache) - end + # MPI mortar fluxes for divergence stage + # @trixi_timeit timer() "MPI mortar flux divergence" begin + # calc_mpi_mortar_flux_divergence!(cache.elements.surface_flux_values, + # mesh, equations_parabolic, dg.mortar, + # dg, parabolic_scheme, cache) + # end - return nothing + # Finish divergence-stage MPI send + @trixi_timeit timer() "finish MPI send divergence" begin + finish_mpi_send!(cache.mpi_cache) end - function calc_gradient_local!(gradients, u_transformed, t, - mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - equations_parabolic, boundary_conditions_parabolic, - dg::DG, parabolic_scheme, cache) - # Reset gradients - @trixi_timeit timer() "reset gradients" begin - reset_gradients!(gradients, dg, cache) - end + # + # Stage 4: final assembly + # + @trixi_timeit timer() "surface integral" begin + calc_surface_integral!(nothing, du, u, mesh, equations_parabolic, + dg.surface_integral, dg, cache) + end - # Calculate volume integral - @trixi_timeit timer() "volume integral" begin - calc_volume_integral_gradient!(gradients, u_transformed, - mesh, equations_parabolic, dg, cache) - end + @trixi_timeit timer() "Jacobian" begin + apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache) + end - # Prolong solution to interfaces. - # This reuses `prolong2interfaces` for the purely hyperbolic case. - @trixi_timeit timer() "prolong2interfaces" begin - prolong2interfaces!(nothing, cache, u_transformed, mesh, - equations_parabolic, dg) - end + @trixi_timeit timer() "source terms parabolic" begin + calc_sources_parabolic!(du, u, gradients, t, source_terms_parabolic, + equations_parabolic, dg, cache) + end - # Calculate interface fluxes for the gradients - @trixi_timeit timer() "interface flux" begin - @unpack surface_flux_values = cache.elements - calc_interface_flux_gradient!(surface_flux_values, mesh, equations_parabolic, - dg, parabolic_scheme, cache) - end + return nothing +end - # Prolong solution to boundaries. - # This reuses `prolong2boundaries` for the purely hyperbolic case. - @trixi_timeit timer() "prolong2boundaries" begin - prolong2boundaries!(cache, u_transformed, mesh, - equations_parabolic, dg) - end +function calc_gradient_local!(gradients, u_transformed, t, + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + equations_parabolic, boundary_conditions_parabolic, + dg::DG, parabolic_scheme, cache) + # Reset gradients + @trixi_timeit timer() "reset gradients" begin + reset_gradients!(gradients, dg, cache) + end - # Calculate boundary fluxes - @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux_gradient!(cache, t, boundary_conditions_parabolic, - mesh, equations_parabolic, dg.surface_integral, - dg) - end + # Calculate volume integral + @trixi_timeit timer() "volume integral" begin + calc_volume_integral_gradient!(gradients, u_transformed, + mesh, equations_parabolic, dg, cache) + end - # Prolong solution to mortars. - # This reuses `prolong2mortars` for the purely hyperbolic case. - @trixi_timeit timer() "prolong2mortars" begin - prolong2mortars!(cache, u_transformed, mesh, equations_parabolic, - dg.mortar, dg) - end + # Prolong solution to interfaces. + # This reuses `prolong2interfaces` for the purely hyperbolic case. + @trixi_timeit timer() "prolong2interfaces" begin + prolong2interfaces!(nothing, cache, u_transformed, mesh, + equations_parabolic, dg) + end - # Calculate mortar fluxes - @trixi_timeit timer() "mortar flux" begin - calc_mortar_flux_gradient!(cache.elements.surface_flux_values, - mesh, equations_parabolic, dg.mortar, - dg, parabolic_scheme, cache) - end + # Calculate interface fluxes for the gradients + @trixi_timeit timer() "interface flux" begin + @unpack surface_flux_values = cache.elements + calc_interface_flux_gradient!(surface_flux_values, mesh, equations_parabolic, + dg, parabolic_scheme, cache) + end - return nothing + # Prolong solution to boundaries. + # This reuses `prolong2boundaries` for the purely hyperbolic case. + @trixi_timeit timer() "prolong2boundaries" begin + prolong2boundaries!(cache, u_transformed, mesh, + equations_parabolic, dg) end - function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, - equations_parabolic, dg::DG) - @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces - @unpack contravariant_vectors = cache.elements - index_range = eachnode(dg) + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" begin + calc_boundary_flux_gradient!(cache, t, boundary_conditions_parabolic, + mesh, equations_parabolic, dg.surface_integral, + dg) + end - flux_parabolic_x, flux_parabolic_y = flux_parabolic + # Prolong solution to mortars. + # This reuses `prolong2mortars` for the purely hyperbolic case. + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars!(cache, u_transformed, mesh, equations_parabolic, + dg.mortar, dg) + end - @threaded for interface in eachmpiinterface(dg, cache) - local_element = local_neighbor_ids[interface] - local_indices = node_indices[interface] - local_direction = indices2direction(local_indices) - local_side = local_sides[interface] - orientationFactor = local_side == 1 ? 1 : -1 + # Calculate mortar fluxes + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux_gradient!(cache.elements.surface_flux_values, + mesh, equations_parabolic, dg.mortar, + dg, parabolic_scheme, cache) + end - i_start, i_step = index_to_start_step_2d(local_indices[1], index_range) - j_start, j_step = index_to_start_step_2d(local_indices[2], index_range) + return nothing +end + +function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, dg::DG) + @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + + flux_parabolic_x, flux_parabolic_y = flux_parabolic + + @threaded for interface in eachmpiinterface(dg, cache) + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + local_direction = indices2direction(local_indices) + local_side = local_sides[interface] + orientationFactor = local_side == 1 ? 1 : -1 + + i_start, i_step = index_to_start_step_2d(local_indices[1], index_range) + j_start, j_step = index_to_start_step_2d(local_indices[2], index_range) + + i_elem = i_start + j_elem = j_start + + for i in eachnode(dg) + normal_direction = get_normal_direction(local_direction, + contravariant_vectors, + i_elem, j_elem, + local_element) + + for v in eachvariable(equations_parabolic) + flux_visc = SVector(flux_parabolic_x[v, i_elem, j_elem, local_element], + flux_parabolic_y[v, i_elem, j_elem, local_element]) + + cache.mpi_interfaces.u[local_side, v, i, interface] = orientationFactor * + dot(flux_visc, + normal_direction) + end - i_elem = i_start - j_elem = j_start + i_elem += i_step + j_elem += j_step + end + end - for i in eachnode(dg) - normal_direction = get_normal_direction(local_direction, - contravariant_vectors, - i_elem, j_elem, - local_element) + return nothing +end - for v in eachvariable(equations_parabolic) - flux_visc = SVector(flux_parabolic_x[v, i_elem, j_elem, local_element], - flux_parabolic_y[v, i_elem, j_elem, local_element]) +function calc_mpi_interface_flux_gradient!(surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, + dg::DG, parabolic_scheme, cache) + @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + index_end = last(index_range) + + @threaded for interface in eachmpiinterface(dg, cache) + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + local_direction = indices2direction(local_indices) + local_side = local_sides[interface] + + # Create the local i,j indexing on the local element used to pull normal direction information + i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], + index_range) + j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], + index_range) + + i_element = i_element_start + j_element = j_element_start + + # Initiate the node index to be used in the surface for loop, + # the surface flux storage must be indexed in alignment with the local element indexing + if :i_backward in local_indices + surface_node = index_end + surface_node_step = -1 + else + surface_node = 1 + surface_node_step = 1 + end - cache.mpi_interfaces.u[local_side, v, i, interface] = orientationFactor * - dot(flux_visc, - normal_direction) - end + for i in eachnode(dg) + normal_direction = get_normal_direction(local_direction, + contravariant_vectors, + i_element, j_element, + local_element) - i_elem += i_step - j_elem += j_step - end - end + calc_mpi_interface_flux_gradient!(surface_flux_values, mesh, + equations_parabolic, + dg, parabolic_scheme, cache, + interface, normal_direction, + i, local_side, + surface_node, + local_direction, local_element) - return nothing - end - - function calc_mpi_interface_flux_gradient!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, - equations_parabolic, - dg::DG, parabolic_scheme, cache) - @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces - @unpack contravariant_vectors = cache.elements - index_range = eachnode(dg) - index_end = last(index_range) - - @threaded for interface in eachmpiinterface(dg, cache) - local_element = local_neighbor_ids[interface] - local_indices = node_indices[interface] - local_direction = indices2direction(local_indices) - local_side = local_sides[interface] - - # Create the local i,j indexing on the local element used to pull normal direction information - i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], - index_range) - j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], - index_range) - - i_element = i_element_start - j_element = j_element_start - - # Initiate the node index to be used in the surface for loop, - # the surface flux storage must be indexed in alignment with the local element indexing - if :i_backward in local_indices - surface_node = index_end - surface_node_step = -1 - else - surface_node = 1 - surface_node_step = 1 - end + # Increment local element indices to pull the normal direction + i_element += i_element_step + j_element += j_element_step - for i in eachnode(dg) - normal_direction = get_normal_direction(local_direction, - contravariant_vectors, - i_element, j_element, - local_element) - - calc_mpi_interface_flux_gradient!(surface_flux_values, mesh, - equations_parabolic, - dg, parabolic_scheme, cache, - interface, normal_direction, - i, local_side, - surface_node, - local_direction, local_element) - - # Increment local element indices to pull the normal direction - i_element += i_element_step - j_element += j_element_step - - # Increment the surface node index along the local element - surface_node += surface_node_step - end + # Increment the surface node index along the local element + surface_node += surface_node_step end - - return nothing end - @inline function calc_mpi_interface_flux_gradient!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, - equations_parabolic, - dg::DG, parabolic_scheme, cache, - interface_index, normal_direction, - interface_node_index, local_side, - surface_node_index, - local_direction_index, - local_element_index) - @unpack u = cache.mpi_interfaces - - u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, - interface_node_index, - interface_index) - - flux_ = flux_parabolic(u_ll, u_rr, normal_direction, Gradient(), - equations_parabolic, parabolic_scheme) + return nothing +end + +@inline function calc_mpi_interface_flux_gradient!(surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, + dg::DG, parabolic_scheme, cache, + interface_index, normal_direction, + interface_node_index, local_side, + surface_node_index, + local_direction_index, + local_element_index) + @unpack u = cache.mpi_interfaces + + u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, + interface_node_index, + interface_index) + + flux_ = flux_parabolic(u_ll, u_rr, normal_direction, Gradient(), + equations_parabolic, parabolic_scheme) + + for v in eachvariable(equations_parabolic) + surface_flux_values[v, surface_node_index, + local_direction_index, local_element_index] = flux_[v] + end - for v in eachvariable(equations_parabolic) - surface_flux_values[v, surface_node_index, - local_direction_index, local_element_index] = flux_[v] + return nothing +end + +function calc_mpi_interface_flux_divergence!(surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, + dg::DG, parabolic_scheme, cache) + @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + index_end = last(index_range) + + @threaded for interface in eachmpiinterface(dg, cache) + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + local_direction = indices2direction(local_indices) + local_side = local_sides[interface] + + i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], + index_range) + j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], + index_range) + + i_element = i_element_start + j_element = j_element_start + + # Initiate the node index to be used in the surface for loop, + # the surface flux storage must be indexed in alignment with the local element indexing + if :i_backward in local_indices + surface_node = index_end + surface_node_step = -1 + else + surface_node = 1 + surface_node_step = 1 end - return nothing - end - - function calc_mpi_interface_flux_divergence!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, - equations_parabolic, - dg::DG, parabolic_scheme, cache) - @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces - @unpack contravariant_vectors = cache.elements - index_range = eachnode(dg) - index_end = last(index_range) - - @threaded for interface in eachmpiinterface(dg, cache) - local_element = local_neighbor_ids[interface] - local_indices = node_indices[interface] - local_direction = indices2direction(local_indices) - local_side = local_sides[interface] - - i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], - index_range) - j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], - index_range) - - i_element = i_element_start - j_element = j_element_start - - # Initiate the node index to be used in the surface for loop, - # the surface flux storage must be indexed in alignment with the local element indexing - if :i_backward in local_indices - surface_node = index_end - surface_node_step = -1 - else - surface_node = 1 - surface_node_step = 1 - end - - for i in eachnode(dg) - normal_direction = get_normal_direction(local_direction, - contravariant_vectors, - i_element, j_element, - local_element) + for i in eachnode(dg) + normal_direction = get_normal_direction(local_direction, + contravariant_vectors, + i_element, j_element, + local_element) - calc_mpi_interface_flux_divergence!(surface_flux_values, mesh, - equations_parabolic, - dg, parabolic_scheme, cache, - interface, normal_direction, - i, local_side, - surface_node, - local_direction, local_element) + calc_mpi_interface_flux_divergence!(surface_flux_values, mesh, + equations_parabolic, + dg, parabolic_scheme, cache, + interface, normal_direction, + i, local_side, + surface_node, + local_direction, local_element) - i_element += i_element_step - j_element += j_element_step + i_element += i_element_step + j_element += j_element_step - surface_node += surface_node_step - end + surface_node += surface_node_step end - - return nothing end - @inline function calc_mpi_interface_flux_divergence!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, - equations_parabolic, - dg::DG, parabolic_scheme, cache, - interface_index, normal_direction, - interface_node_index, local_side, - surface_node_index, - local_direction_index, - local_element_index) - @unpack u = cache.mpi_interfaces - - parabolic_flux_normal_ll, parabolic_flux_normal_rr = get_surface_node_vars(u, - equations_parabolic, - dg, - interface_node_index, - interface_index) - + return nothing +end + +@inline function calc_mpi_interface_flux_divergence!(surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations_parabolic, + dg::DG, parabolic_scheme, cache, + interface_index, normal_direction, + interface_node_index, local_side, + surface_node_index, + local_direction_index, + local_element_index) + @unpack u = cache.mpi_interfaces + + parabolic_flux_normal_ll, parabolic_flux_normal_rr = get_surface_node_vars(u, + equations_parabolic, + dg, + interface_node_index, + interface_index) + + flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, + normal_direction, Divergence(), + equations_parabolic, parabolic_scheme) + + if local_side == 1 flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, normal_direction, Divergence(), equations_parabolic, parabolic_scheme) + sign_ = 1 + else + # side 2 must also use primary-oriented normal + flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, + -normal_direction, Divergence(), + equations_parabolic, parabolic_scheme) + sign_ = -1 + end - if local_side == 1 - flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, - normal_direction, Divergence(), - equations_parabolic, parabolic_scheme) - sign_ = 1 - else - # side 2 must also use primary-oriented normal - flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, - -normal_direction, Divergence(), - equations_parabolic, parabolic_scheme) - sign_ = -1 - end - - for v in eachvariable(equations_parabolic) - surface_flux_values[v, surface_node_index, - local_direction_index, local_element_index] = sign_ * flux_[v] - end - - return nothing + for v in eachvariable(equations_parabolic) + surface_flux_values[v, surface_node_index, + local_direction_index, local_element_index] = sign_ * flux_[v] end + + return nothing +end end #muladd From 206eea318376960e5a4fe7dccf7689a6a636e206 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Thu, 2 Apr 2026 16:33:42 +0200 Subject: [PATCH 04/17] cleanup style --- src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index ac432fd31ce..82649acf1e5 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -281,7 +281,7 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, local_indices = node_indices[interface] local_direction = indices2direction(local_indices) local_side = local_sides[interface] - orientationFactor = local_side == 1 ? 1 : -1 + orientation_factor = local_side == 1 ? 1 : -1 i_start, i_step = index_to_start_step_2d(local_indices[1], index_range) j_start, j_step = index_to_start_step_2d(local_indices[2], index_range) @@ -299,7 +299,7 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, flux_visc = SVector(flux_parabolic_x[v, i_elem, j_elem, local_element], flux_parabolic_y[v, i_elem, j_elem, local_element]) - cache.mpi_interfaces.u[local_side, v, i, interface] = orientationFactor * + cache.mpi_interfaces.u[local_side, v, i, interface] = orientation_factor * dot(flux_visc, normal_direction) end @@ -476,10 +476,6 @@ end interface_node_index, interface_index) - flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, - normal_direction, Divergence(), - equations_parabolic, parabolic_scheme) - if local_side == 1 flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, normal_direction, Divergence(), From 606174f317020f41f2ab66fa8ca6543716795f86 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Fri, 3 Apr 2026 22:33:11 +0200 Subject: [PATCH 05/17] increased allocation limit (checked similar allocations in MPI as hyperbolic communication, now just 3 x communication) --- test/test_mpi_p4est_parabolic_2d.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/test_mpi_p4est_parabolic_2d.jl b/test/test_mpi_p4est_parabolic_2d.jl index 5d61c44a981..3c97c9b41c5 100644 --- a/test/test_mpi_p4est_parabolic_2d.jl +++ b/test/test_mpi_p4est_parabolic_2d.jl @@ -26,8 +26,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "p4est_2d_dgsem") ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 1000) - @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + @test_allocations(Trixi.rhs!, semi, sol, 1500) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1500) end @trixi_testset "P4estMesh2D: elixir_navierstokes_convergence_nonperiodic.jl" begin @@ -48,8 +48,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "p4est_2d_dgsem") ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 1000) - @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + @test_allocations(Trixi.rhs!, semi, sol, 1500) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1500) end @trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_curved.jl" begin @@ -61,8 +61,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "p4est_2d_dgsem") linf=[0.14186297438393505]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 1000) - @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + @test_allocations(Trixi.rhs!, semi, sol, 1500) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1500) end @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin @@ -74,8 +74,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "p4est_2d_dgsem") linf=[0.008154128363741964]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 1000) - @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + @test_allocations(Trixi.rhs!, semi, sol, 1500) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1500) end @trixi_testset "elixir_navierstokes_NACA0012airfoil_mach08.jl" begin @@ -119,8 +119,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "p4est_2d_dgsem") # (e.g., from type instabilities) # We move these tests here to avoid modifying values used # to compute the drag/lift coefficients above. - @test_allocations(Trixi.rhs!, semi, sol, 1000) - @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + @test_allocations(Trixi.rhs!, semi, sol, 1500) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1500) end end end # module From 4a02f5432dec8b2d8b3e20bda25a0148d27b54cb Mon Sep 17 00:00:00 2001 From: "T. Jeremy P. Karpowski" <49682045+TJP-Karpowski@users.noreply.github.com> Date: Wed, 8 Apr 2026 10:57:37 +0200 Subject: [PATCH 06/17] Apply suggestions from code review Removed dispatch to unsupoorted T8codeMeshParallel Co-authored-by: Daniel Doehring --- .../dgsem_p4est/dg_2d_parabolic_parallel.jl | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index 82649acf1e5..08bce319ec0 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -5,7 +5,7 @@ @muladd begin #! format: noindent function rhs_parabolic!(du, u, t, - mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + mesh::P4estMeshParallel{2}, equations_parabolic::AbstractEquationsParabolic, boundary_conditions_parabolic, source_terms_parabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) @@ -207,7 +207,7 @@ function rhs_parabolic!(du, u, t, end function calc_gradient_local!(gradients, u_transformed, t, - mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + mesh::P4estMeshParallel{2}, equations_parabolic, boundary_conditions_parabolic, dg::DG, parabolic_scheme, cache) # Reset gradients @@ -267,8 +267,7 @@ function calc_gradient_local!(gradients, u_transformed, t, end function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @unpack contravariant_vectors = cache.elements @@ -313,8 +312,7 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, end function calc_mpi_interface_flux_gradient!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG, parabolic_scheme, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @@ -374,8 +372,7 @@ function calc_mpi_interface_flux_gradient!(surface_flux_values, end @inline function calc_mpi_interface_flux_gradient!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG, parabolic_scheme, cache, interface_index, normal_direction, @@ -401,8 +398,7 @@ end end function calc_mpi_interface_flux_divergence!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG, parabolic_scheme, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @@ -459,8 +455,7 @@ function calc_mpi_interface_flux_divergence!(surface_flux_values, end @inline function calc_mpi_interface_flux_divergence!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG, parabolic_scheme, cache, interface_index, normal_direction, From 7b36fdccdb840bbf2550faf0297ddd7421c0a9fa Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Wed, 8 Apr 2026 11:43:50 +0200 Subject: [PATCH 07/17] added comments to clarify sign-flip --- .../dgsem_p4est/dg_2d_parabolic_parallel.jl | 31 ++++++++++++------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index 08bce319ec0..e160b19c4c1 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -5,7 +5,7 @@ @muladd begin #! format: noindent function rhs_parabolic!(du, u, t, - mesh::P4estMeshParallel{2}, + mesh::P4estMeshParallel{2}, equations_parabolic::AbstractEquationsParabolic, boundary_conditions_parabolic, source_terms_parabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) @@ -207,7 +207,7 @@ function rhs_parabolic!(du, u, t, end function calc_gradient_local!(gradients, u_transformed, t, - mesh::P4estMeshParallel{2}, + mesh::P4estMeshParallel{2}, equations_parabolic, boundary_conditions_parabolic, dg::DG, parabolic_scheme, cache) # Reset gradients @@ -267,7 +267,7 @@ function calc_gradient_local!(gradients, u_transformed, t, end function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, - mesh::P4estMeshParallel{2}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @unpack contravariant_vectors = cache.elements @@ -276,10 +276,17 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for interface in eachmpiinterface(dg, cache) + # Copy solution data from the local element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. local_element = local_neighbor_ids[interface] local_indices = node_indices[interface] local_direction = indices2direction(local_indices) local_side = local_sides[interface] + # store the normal flux with respect to the primary normal direction, + # which is the negative of the secondary normal direction orientation_factor = local_side == 1 ? 1 : -1 i_start, i_step = index_to_start_step_2d(local_indices[1], index_range) @@ -297,7 +304,7 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, for v in eachvariable(equations_parabolic) flux_visc = SVector(flux_parabolic_x[v, i_elem, j_elem, local_element], flux_parabolic_y[v, i_elem, j_elem, local_element]) - + #writes data to the mpi cache exchanged with the other ranks cache.mpi_interfaces.u[local_side, v, i, interface] = orientation_factor * dot(flux_visc, normal_direction) @@ -312,7 +319,7 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, end function calc_mpi_interface_flux_gradient!(surface_flux_values, - mesh::P4estMeshParallel{2}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG, parabolic_scheme, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @@ -372,7 +379,7 @@ function calc_mpi_interface_flux_gradient!(surface_flux_values, end @inline function calc_mpi_interface_flux_gradient!(surface_flux_values, - mesh::P4estMeshParallel{2}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG, parabolic_scheme, cache, interface_index, normal_direction, @@ -398,7 +405,7 @@ end end function calc_mpi_interface_flux_divergence!(surface_flux_values, - mesh::P4estMeshParallel{2}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG, parabolic_scheme, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @@ -455,7 +462,7 @@ function calc_mpi_interface_flux_divergence!(surface_flux_values, end @inline function calc_mpi_interface_flux_divergence!(surface_flux_values, - mesh::P4estMeshParallel{2}, + mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG, parabolic_scheme, cache, interface_index, normal_direction, @@ -475,18 +482,20 @@ end flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, normal_direction, Divergence(), equations_parabolic, parabolic_scheme) - sign_ = 1 + orientation_factor = 1 else # side 2 must also use primary-oriented normal flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, -normal_direction, Divergence(), equations_parabolic, parabolic_scheme) - sign_ = -1 + # Sign flip required for divergence calculation since the divergence interface flux + # involves the normal direction. + orientation_factor = -1 end for v in eachvariable(equations_parabolic) surface_flux_values[v, surface_node_index, - local_direction_index, local_element_index] = sign_ * flux_[v] + local_direction_index, local_element_index] = orientation_factor * flux_[v] end return nothing From cb9f5191e8c68d77d99ba96012b5c2e37dd5f864 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Wed, 8 Apr 2026 15:52:06 +0200 Subject: [PATCH 08/17] extended comment on orientation_factor --- src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index e160b19c4c1..56fc7175700 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -305,6 +305,8 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, flux_visc = SVector(flux_parabolic_x[v, i_elem, j_elem, local_element], flux_parabolic_y[v, i_elem, j_elem, local_element]) #writes data to the mpi cache exchanged with the other ranks + #Side 1 and 2 must be consistent, thus the orientation_factor changes the orientation + # So flux entering side 1 leaves side 2. cache.mpi_interfaces.u[local_side, v, i, interface] = orientation_factor * dot(flux_visc, normal_direction) @@ -489,7 +491,7 @@ end -normal_direction, Divergence(), equations_parabolic, parabolic_scheme) # Sign flip required for divergence calculation since the divergence interface flux - # involves the normal direction. + # involves the normal direction. local_side=2 must be inverted to stay consistent. orientation_factor = -1 end From 0da652eb16e9f74907a87d99d42900dd7041ddf3 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Wed, 8 Apr 2026 18:08:53 +0200 Subject: [PATCH 09/17] renamed surface integral MPI accumulate --- src/callbacks_step/analysis_surface_integral_2d.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index 2efb827e572..03f42ac0110 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -250,7 +250,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, j_node += j_node_step end end - return distribute_surface_integral(surface_integral, mesh) + return accumulate(surface_integral, mesh) end # 2D version of the `analyze` function for `AnalysisSurfaceIntegral` of viscous, i.e., @@ -320,16 +320,14 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, du, u, t, j_node += j_node_step end end - return distribute_surface_integral(surface_integral, mesh) + return accumulate(surface_integral, mesh) end # Serial/default: do nothing -distribute_surface_integral(val, mesh) = val +accumulate(val, mesh) = val # Parallel: sum over all ranks -function distribute_surface_integral(val, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}) +function accumulate(val, mesh::Union{P4estMeshParallel{2}}) comm = MPI.COMM_WORLD buf = [val] MPI.Allreduce!(buf, MPI.SUM, comm) From 68208b0596068be5c89e13026869398bd45b54a0 Mon Sep 17 00:00:00 2001 From: "T. Jeremy P. Karpowski" <49682045+TJP-Karpowski@users.noreply.github.com> Date: Wed, 8 Apr 2026 19:56:06 +0200 Subject: [PATCH 10/17] Apply suggestions from code review adhere to mpi internal functions Co-authored-by: Hendrik Ranocha --- src/callbacks_step/analysis_surface_integral_2d.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index 03f42ac0110..f93d0148c63 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -320,7 +320,10 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, du, u, t, j_node += j_node_step end end - return accumulate(surface_integral, mesh) + if mpi_isparallel() + surface_integral = MPI.Allreduce!(Ref(surface_integral), +, mpi_comm())[] + end + return surface_integral end # Serial/default: do nothing @@ -328,9 +331,9 @@ accumulate(val, mesh) = val # Parallel: sum over all ranks function accumulate(val, mesh::Union{P4estMeshParallel{2}}) - comm = MPI.COMM_WORLD - buf = [val] - MPI.Allreduce!(buf, MPI.SUM, comm) - return buf[1] + comm = mpi_comm() + buf = Ref(val) + global_val = MPI.Allreduce!(buf, +, comm)[] + return global_val end end # muladd From 71dd53dae9e99650df6c6511da9d1630ef8cfbc9 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Wed, 8 Apr 2026 20:22:19 +0200 Subject: [PATCH 11/17] removed unnecessary function --- .../analysis_surface_integral_2d.jl | 24 ++-- .../dgsem_p4est/dg_2d_parabolic_parallel.jl | 111 +++++------------- 2 files changed, 39 insertions(+), 96 deletions(-) diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index f93d0148c63..003c279f8bd 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -250,7 +250,10 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, j_node += j_node_step end end - return accumulate(surface_integral, mesh) + if mpi_isparallel() + surface_integral = MPI.Allreduce!(Ref(surface_integral), +, mpi_comm())[] + end + return surface_integral end # 2D version of the `analyze` function for `AnalysisSurfaceIntegral` of viscous, i.e., @@ -320,20 +323,9 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, du, u, t, j_node += j_node_step end end - if mpi_isparallel() - surface_integral = MPI.Allreduce!(Ref(surface_integral), +, mpi_comm())[] - end - return surface_integral -end - -# Serial/default: do nothing -accumulate(val, mesh) = val - -# Parallel: sum over all ranks -function accumulate(val, mesh::Union{P4estMeshParallel{2}}) - comm = mpi_comm() - buf = Ref(val) - global_val = MPI.Allreduce!(buf, +, comm)[] - return global_val + if mpi_isparallel() + surface_integral = MPI.Allreduce!(Ref(surface_integral), +, mpi_comm())[] + end + return surface_integral end end # muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index 56fc7175700..9c20ea7b497 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -326,6 +326,7 @@ function calc_mpi_interface_flux_gradient!(surface_flux_values, dg::DG, parabolic_scheme, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @unpack contravariant_vectors = cache.elements + @unpack u = cache.mpi_interfaces index_range = eachnode(dg) index_end = last(index_range) @@ -360,13 +361,17 @@ function calc_mpi_interface_flux_gradient!(surface_flux_values, i_element, j_element, local_element) - calc_mpi_interface_flux_gradient!(surface_flux_values, mesh, - equations_parabolic, - dg, parabolic_scheme, cache, - interface, normal_direction, - i, local_side, - surface_node, - local_direction, local_element) + u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, + i, + interface) + + flux_ = flux_parabolic(u_ll, u_rr, normal_direction, Gradient(), + equations_parabolic, parabolic_scheme) + + for v in eachvariable(equations_parabolic) + surface_flux_values[v, surface_node, + local_direction, local_element] = flux_[v] + end # Increment local element indices to pull the normal direction i_element += i_element_step @@ -380,38 +385,13 @@ function calc_mpi_interface_flux_gradient!(surface_flux_values, return nothing end -@inline function calc_mpi_interface_flux_gradient!(surface_flux_values, - mesh::P4estMeshParallel{2}, - equations_parabolic, - dg::DG, parabolic_scheme, cache, - interface_index, normal_direction, - interface_node_index, local_side, - surface_node_index, - local_direction_index, - local_element_index) - @unpack u = cache.mpi_interfaces - - u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, - interface_node_index, - interface_index) - - flux_ = flux_parabolic(u_ll, u_rr, normal_direction, Gradient(), - equations_parabolic, parabolic_scheme) - - for v in eachvariable(equations_parabolic) - surface_flux_values[v, surface_node_index, - local_direction_index, local_element_index] = flux_[v] - end - - return nothing -end - function calc_mpi_interface_flux_divergence!(surface_flux_values, mesh::P4estMeshParallel{2}, equations_parabolic, dg::DG, parabolic_scheme, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @unpack contravariant_vectors = cache.elements + @unpack u = cache.mpi_interfaces index_range = eachnode(dg) index_end = last(index_range) @@ -445,13 +425,24 @@ function calc_mpi_interface_flux_divergence!(surface_flux_values, i_element, j_element, local_element) - calc_mpi_interface_flux_divergence!(surface_flux_values, mesh, - equations_parabolic, - dg, parabolic_scheme, cache, - interface, normal_direction, - i, local_side, - surface_node, - local_direction, local_element) + parabolic_flux_normal_ll, parabolic_flux_normal_rr = get_surface_node_vars(u, + equations_parabolic, + dg, + i, + interface) + + # side 2 must also use primary-oriented normal + # Sign flip required for divergence calculation since the divergence interface flux + # involves the normal direction. local_side=2 must be inverted to stay consistent. + orientation_factor = local_side == 1 ? 1 : -1 + flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, + orientation_factor * normal_direction, Divergence(), + equations_parabolic, parabolic_scheme) + + for v in eachvariable(equations_parabolic) + surface_flux_values[v, surface_node, + local_direction, local_element] = orientation_factor * flux_[v] + end i_element += i_element_step j_element += j_element_step @@ -462,44 +453,4 @@ function calc_mpi_interface_flux_divergence!(surface_flux_values, return nothing end - -@inline function calc_mpi_interface_flux_divergence!(surface_flux_values, - mesh::P4estMeshParallel{2}, - equations_parabolic, - dg::DG, parabolic_scheme, cache, - interface_index, normal_direction, - interface_node_index, local_side, - surface_node_index, - local_direction_index, - local_element_index) - @unpack u = cache.mpi_interfaces - - parabolic_flux_normal_ll, parabolic_flux_normal_rr = get_surface_node_vars(u, - equations_parabolic, - dg, - interface_node_index, - interface_index) - - if local_side == 1 - flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, - normal_direction, Divergence(), - equations_parabolic, parabolic_scheme) - orientation_factor = 1 - else - # side 2 must also use primary-oriented normal - flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, - -normal_direction, Divergence(), - equations_parabolic, parabolic_scheme) - # Sign flip required for divergence calculation since the divergence interface flux - # involves the normal direction. local_side=2 must be inverted to stay consistent. - orientation_factor = -1 - end - - for v in eachvariable(equations_parabolic) - surface_flux_values[v, surface_node_index, - local_direction_index, local_element_index] = orientation_factor * flux_[v] - end - - return nothing -end end #muladd From 91772d635577210194b5df4b13b4b387033c24c9 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Wed, 8 Apr 2026 21:15:18 +0200 Subject: [PATCH 12/17] moved mpi send to end --- src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index 9c20ea7b497..1009f41d908 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -57,10 +57,6 @@ function rhs_parabolic!(du, u, t, @trixi_timeit timer() "finish MPI receive gradient" begin finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end - # Finish gradient-stage MPI send - @trixi_timeit timer() "finish MPI send gradient" begin - finish_mpi_send!(cache.mpi_cache) - end # MPI interface fluxes for gradient stage @trixi_timeit timer() "MPI interface flux gradient" begin calc_mpi_interface_flux_gradient!(cache.elements.surface_flux_values, @@ -86,6 +82,12 @@ function rhs_parabolic!(du, u, t, apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg, cache) end + + # Finish gradient-stage MPI send + @trixi_timeit timer() "finish MPI send gradient" begin + finish_mpi_send!(cache.mpi_cache) + end + # Stage 2: local parabolic flux construction # @trixi_timeit timer() "calculate parabolic fluxes" begin From 7c1638fe5147b02f8b01bef485074d5c20c3b964 Mon Sep 17 00:00:00 2001 From: "T. Jeremy P. Karpowski" <49682045+TJP-Karpowski@users.noreply.github.com> Date: Thu, 9 Apr 2026 10:41:59 +0200 Subject: [PATCH 13/17] Apply suggestions from code review applied suggested comment improvments Co-authored-by: Daniel Doehring --- .../dgsem_p4est/dg_2d_parabolic_parallel.jl | 56 ++++++++----------- 1 file changed, 22 insertions(+), 34 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index 1009f41d908..b44b75b51e2 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -12,18 +12,14 @@ function rhs_parabolic!(du, u, t, @unpack parabolic_container = cache_parabolic @unpack u_transformed, gradients, flux_parabolic = parabolic_container - # Stage 0: local variable transform - # @trixi_timeit timer() "transform variables" begin transform_variables!(u_transformed, u, mesh, equations_parabolic, dg, cache) end - # - # Stage 1: gradient computation - # + ### Gradient computation ### - # Start gradient-stage MPI receive + # Start gradient MPI receive @trixi_timeit timer() "start MPI receive gradient" begin start_mpi_receive!(cache.mpi_cache) end @@ -41,7 +37,7 @@ function rhs_parabolic!(du, u, t, dg.surface_integral, dg) end - # Start gradient-stage MPI send + # Start gradient MPI send @trixi_timeit timer() "start MPI send gradient" begin start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end @@ -53,18 +49,19 @@ function rhs_parabolic!(du, u, t, dg, parabolic_scheme, cache) end - # Finish gradient-stage MPI receive + # Finish gradient MPI receive @trixi_timeit timer() "finish MPI receive gradient" begin finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end - # MPI interface fluxes for gradient stage - @trixi_timeit timer() "MPI interface flux gradient" begin - calc_mpi_interface_flux_gradient!(cache.elements.surface_flux_values, - mesh, equations_parabolic, - dg, parabolic_scheme, cache) - end + + # MPI interface fluxes for gradients + @trixi_timeit timer() "MPI interface flux gradient" begin + calc_mpi_interface_flux_gradient!(cache.elements.surface_flux_values, + mesh, equations_parabolic, + dg, parabolic_scheme, cache) + end - # MPI mortar fluxes for gradient stage + # MPI mortar fluxes for gradients # @trixi_timeit timer() "MPI mortar flux gradient" begin # calc_mpi_mortar_flux_gradient!(cache.elements.surface_flux_values, # mesh, equations_parabolic, dg.mortar, @@ -95,11 +92,9 @@ function rhs_parabolic!(du, u, t, equations_parabolic, dg, cache) end - # - # Stage 3: divergence of parabolic fluxes - # + ### Divergence of parabolic/gradient fluxes ### - # Start divergence-stage MPI receive + # Start divergence MPI receive @trixi_timeit timer() "start MPI receive divergence" begin start_mpi_receive!(cache.mpi_cache) end @@ -124,8 +119,7 @@ function rhs_parabolic!(du, u, t, @trixi_timeit timer() "prolong2mpiinterfaces divergence" begin prolong2mpiinterfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) end - ########################## Divergence ################################# - # Start divergence-stage MPI send + # Start divergence MPI send @trixi_timeit timer() "start MPI send divergence" begin start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end @@ -169,28 +163,25 @@ function rhs_parabolic!(du, u, t, finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end - # MPI interface fluxes for divergence stage + # MPI interface fluxes for divergence @trixi_timeit timer() "MPI interface flux divergence" begin calc_mpi_interface_flux_divergence!(cache.elements.surface_flux_values, mesh, equations_parabolic, dg, parabolic_scheme, cache) end - # MPI mortar fluxes for divergence stage + # MPI mortar fluxes for divergence # @trixi_timeit timer() "MPI mortar flux divergence" begin # calc_mpi_mortar_flux_divergence!(cache.elements.surface_flux_values, # mesh, equations_parabolic, dg.mortar, # dg, parabolic_scheme, cache) # end - # Finish divergence-stage MPI send + # Finish divergence MPI send @trixi_timeit timer() "finish MPI send divergence" begin finish_mpi_send!(cache.mpi_cache) end - # - # Stage 4: final assembly - # @trixi_timeit timer() "surface integral" begin calc_surface_integral!(nothing, du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache) @@ -306,8 +297,7 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, for v in eachvariable(equations_parabolic) flux_visc = SVector(flux_parabolic_x[v, i_elem, j_elem, local_element], flux_parabolic_y[v, i_elem, j_elem, local_element]) - #writes data to the mpi cache exchanged with the other ranks - #Side 1 and 2 must be consistent, thus the orientation_factor changes the orientation + # Side 1 and 2 must be consistent, thus the orientation_factor changes the orientation # So flux entering side 1 leaves side 2. cache.mpi_interfaces.u[local_side, v, i, interface] = orientation_factor * dot(flux_visc, @@ -364,8 +354,7 @@ function calc_mpi_interface_flux_gradient!(surface_flux_values, local_element) u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, - i, - interface) + i, interface) flux_ = flux_parabolic(u_ll, u_rr, normal_direction, Gradient(), equations_parabolic, parabolic_scheme) @@ -433,9 +422,8 @@ function calc_mpi_interface_flux_divergence!(surface_flux_values, i, interface) - # side 2 must also use primary-oriented normal - # Sign flip required for divergence calculation since the divergence interface flux - # involves the normal direction. local_side=2 must be inverted to stay consistent. + # Sign flip for `local_side = 2` required for divergence calculation since the divergence interface flux + # involves the normal direction. `local_side=2` is thus flipped (opposite of primary side) orientation_factor = local_side == 1 ? 1 : -1 flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, orientation_factor * normal_direction, Divergence(), From 444fc42dd95c14cefc8ad4de9fcc08e6a7d74a08 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Thu, 9 Apr 2026 11:08:43 +0200 Subject: [PATCH 14/17] formatting applied --- .../dgsem_p4est/dg_2d_parabolic_parallel.jl | 44 +++++++++---------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index b44b75b51e2..26dc9fdfbb8 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -17,9 +17,9 @@ function rhs_parabolic!(du, u, t, dg, cache) end - ### Gradient computation ### + ### Gradient computation ### - # Start gradient MPI receive + # Start gradient MPI receive @trixi_timeit timer() "start MPI receive gradient" begin start_mpi_receive!(cache.mpi_cache) end @@ -37,7 +37,7 @@ function rhs_parabolic!(du, u, t, dg.surface_integral, dg) end - # Start gradient MPI send + # Start gradient MPI send @trixi_timeit timer() "start MPI send gradient" begin start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end @@ -49,19 +49,19 @@ function rhs_parabolic!(du, u, t, dg, parabolic_scheme, cache) end - # Finish gradient MPI receive + # Finish gradient MPI receive @trixi_timeit timer() "finish MPI receive gradient" begin finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end - - # MPI interface fluxes for gradients - @trixi_timeit timer() "MPI interface flux gradient" begin - calc_mpi_interface_flux_gradient!(cache.elements.surface_flux_values, - mesh, equations_parabolic, - dg, parabolic_scheme, cache) - end - # MPI mortar fluxes for gradients + # MPI interface fluxes for gradients + @trixi_timeit timer() "MPI interface flux gradient" begin + calc_mpi_interface_flux_gradient!(cache.elements.surface_flux_values, + mesh, equations_parabolic, + dg, parabolic_scheme, cache) + end + + # MPI mortar fluxes for gradients # @trixi_timeit timer() "MPI mortar flux gradient" begin # calc_mpi_mortar_flux_gradient!(cache.elements.surface_flux_values, # mesh, equations_parabolic, dg.mortar, @@ -92,9 +92,9 @@ function rhs_parabolic!(du, u, t, equations_parabolic, dg, cache) end - ### Divergence of parabolic/gradient fluxes ### + ### Divergence of parabolic/gradient fluxes ### - # Start divergence MPI receive + # Start divergence MPI receive @trixi_timeit timer() "start MPI receive divergence" begin start_mpi_receive!(cache.mpi_cache) end @@ -119,7 +119,7 @@ function rhs_parabolic!(du, u, t, @trixi_timeit timer() "prolong2mpiinterfaces divergence" begin prolong2mpiinterfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) end - # Start divergence MPI send + # Start divergence MPI send @trixi_timeit timer() "start MPI send divergence" begin start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end @@ -163,21 +163,21 @@ function rhs_parabolic!(du, u, t, finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end - # MPI interface fluxes for divergence + # MPI interface fluxes for divergence @trixi_timeit timer() "MPI interface flux divergence" begin calc_mpi_interface_flux_divergence!(cache.elements.surface_flux_values, mesh, equations_parabolic, dg, parabolic_scheme, cache) end - # MPI mortar fluxes for divergence + # MPI mortar fluxes for divergence # @trixi_timeit timer() "MPI mortar flux divergence" begin # calc_mpi_mortar_flux_divergence!(cache.elements.surface_flux_values, # mesh, equations_parabolic, dg.mortar, # dg, parabolic_scheme, cache) # end - # Finish divergence MPI send + # Finish divergence MPI send @trixi_timeit timer() "finish MPI send divergence" begin finish_mpi_send!(cache.mpi_cache) end @@ -297,7 +297,7 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, for v in eachvariable(equations_parabolic) flux_visc = SVector(flux_parabolic_x[v, i_elem, j_elem, local_element], flux_parabolic_y[v, i_elem, j_elem, local_element]) - # Side 1 and 2 must be consistent, thus the orientation_factor changes the orientation + # Side 1 and 2 must be consistent, thus the orientation_factor changes the orientation # So flux entering side 1 leaves side 2. cache.mpi_interfaces.u[local_side, v, i, interface] = orientation_factor * dot(flux_visc, @@ -354,7 +354,7 @@ function calc_mpi_interface_flux_gradient!(surface_flux_values, local_element) u_ll, u_rr = get_surface_node_vars(u, equations_parabolic, dg, - i, interface) + i, interface) flux_ = flux_parabolic(u_ll, u_rr, normal_direction, Gradient(), equations_parabolic, parabolic_scheme) @@ -422,8 +422,8 @@ function calc_mpi_interface_flux_divergence!(surface_flux_values, i, interface) - # Sign flip for `local_side = 2` required for divergence calculation since the divergence interface flux - # involves the normal direction. `local_side=2` is thus flipped (opposite of primary side) + # Sign flip for `local_side = 2` required for divergence calculation since the divergence interface flux + # involves the normal direction. `local_side=2` is thus flipped (opposite of primary side) orientation_factor = local_side == 1 ? 1 : -1 flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, orientation_factor * normal_direction, Divergence(), From b4df801387baccbc5877955454a57f9a3bce49d3 Mon Sep 17 00:00:00 2001 From: TJP-Karpowski Date: Sat, 18 Apr 2026 10:54:57 +0200 Subject: [PATCH 15/17] added name to Authors --- AUTHORS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/AUTHORS.md b/AUTHORS.md index b54b8c45f41..bf659a68874 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -37,6 +37,7 @@ are listed in alphabetical order: * Gregor Gassner * Lucas Gemein * Sven Goldberg +* T. Jeremy P. Karpowski * Joshua Lampert * Tristan Montoya * Julia Odenthal From cc6d70f2243e55a12518e6f7a167dfd05616350b Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 20 Apr 2026 09:08:30 +0200 Subject: [PATCH 16/17] Apply suggestions from code review --- .../dgsem_p4est/dg_2d_parabolic_parallel.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index 26dc9fdfbb8..8679bfefbf1 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -80,13 +80,12 @@ function rhs_parabolic!(du, u, t, cache) end - # Finish gradient-stage MPI send + # Finish gradient MPI send @trixi_timeit timer() "finish MPI send gradient" begin finish_mpi_send!(cache.mpi_cache) end - # Stage 2: local parabolic flux construction - # + # Local parabolic flux construction @trixi_timeit timer() "calculate parabolic fluxes" begin calc_parabolic_fluxes!(flux_parabolic, gradients, u_transformed, mesh, equations_parabolic, dg, cache) @@ -119,6 +118,7 @@ function rhs_parabolic!(du, u, t, @trixi_timeit timer() "prolong2mpiinterfaces divergence" begin prolong2mpiinterfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) end + # Start divergence MPI send @trixi_timeit timer() "start MPI send divergence" begin start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) @@ -158,7 +158,7 @@ function rhs_parabolic!(du, u, t, dg, parabolic_scheme, cache) end - # Finish divergence-stage MPI receive + # Finish divergence MPI receive @trixi_timeit timer() "finish MPI receive divergence" begin finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) end @@ -365,6 +365,7 @@ function calc_mpi_interface_flux_gradient!(surface_flux_values, end # Increment local element indices to pull the normal direction + # from the element data i_element += i_element_step j_element += j_element_step @@ -422,8 +423,9 @@ function calc_mpi_interface_flux_divergence!(surface_flux_values, i, interface) - # Sign flip for `local_side = 2` required for divergence calculation since the divergence interface flux - # involves the normal direction. `local_side=2` is thus flipped (opposite of primary side) + # Sign flip for `local_side = 2` required for divergence calculation since + # the divergence interface flux involves the normal direction. + # `local_side=2` is thus flipped (opposite of primary side) orientation_factor = local_side == 1 ? 1 : -1 flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, orientation_factor * normal_direction, Divergence(), From 4bff172752e0c48cc61888ddfd956ec0a847439c Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 20 Apr 2026 09:13:36 +0200 Subject: [PATCH 17/17] Update src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl --- src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index 8679bfefbf1..f7f95710c63 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -297,8 +297,9 @@ function prolong2mpiinterfaces!(cache, flux_parabolic::Tuple, for v in eachvariable(equations_parabolic) flux_visc = SVector(flux_parabolic_x[v, i_elem, j_elem, local_element], flux_parabolic_y[v, i_elem, j_elem, local_element]) - # Side 1 and 2 must be consistent, thus the orientation_factor changes the orientation - # So flux entering side 1 leaves side 2. + # Side 1 and 2 must be consistent, i.e., with their outward-pointing normals. + # Thus, the `orientation_factor` changes the logic such that the + # flux which enters side 1 leaves side 2. cache.mpi_interfaces.u[local_side, v, i, interface] = orientation_factor * dot(flux_visc, normal_direction)