From a424eb1684d1578f31b49f0d34562b9faa7eb07e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 21 Mar 2024 12:39:22 +0100 Subject: [PATCH 01/24] PartitionRHS_Comp_1D --- src/solvers/dgsem_structured/dg_1d.jl | 22 ++++--- src/solvers/dgsem_tree/dg.jl | 12 +++- src/solvers/dgsem_tree/dg_1d.jl | 72 ++++++++++++---------- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 74 +++++++++++++---------- src/solvers/dgsem_tree/dg_2d.jl | 3 +- src/solvers/dgsem_tree/dg_3d.jl | 3 +- src/solvers/fdsbp_tree/fdsbp_1d.jl | 22 +++---- 7 files changed, 122 insertions(+), 86 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index 3d63cc5af36..abdd0fa82cb 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -9,19 +9,23 @@ function rhs!(du, u, t, mesh::StructuredMesh{1}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} + + element_range = eachelement(dg, cache) + # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_range) # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(du, u, mesh, have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache) + dg.volume_integral, dg, cache, element_range) end # Calculate interface and boundary fluxes @trixi_timeit timer() "interface flux" begin - calc_interface_flux!(cache, u, mesh, equations, dg.surface_integral, dg) + calc_interface_flux!(cache, u, mesh, equations, dg.surface_integral, dg, + element_range) end # Calculate boundary fluxes @@ -33,25 +37,27 @@ function rhs!(du, u, t, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_surface_integral!(du, u, mesh, equations, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, element_range) end # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache) + @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, + element_range) # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range) end return nothing end function calc_interface_flux!(cache, u, mesh::StructuredMesh{1}, - equations, surface_integral, dg::DG) + equations, surface_integral, dg::DG, + element_range) @unpack surface_flux = surface_integral - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range left_element = cache.elements.left_neighbors[1, element] if left_element > 0 # left_element = 0 at boundaries diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index ef9a42b4c1a..3b22d40234b 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -15,17 +15,25 @@ function reset_du!(du, dg, cache) return du end +function reset_du!(du, dg, cache, element_range) + @threaded for element in element_range + du[.., element] .= zero(eltype(du)) + end + + return du +end + # pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache) # # Given blending factors `alpha` and the solver `dg`, fill # `element_ids_dg` with the IDs of elements using a pure DG scheme and # `element_ids_dgfv` with the IDs of elements using a blended DG-FV scheme. function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg::DG, - cache) + cache, element_range) empty!(element_ids_dg) empty!(element_ids_dgfv) - for element in eachelement(dg, cache) + for element in element_range # Clip blending factor for values close to zero (-> pure DG) dg_only = isapprox(alpha[element], 0, atol = 1e-12) if dg_only diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 4a0747d1c09..114da5211ea 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -78,33 +78,38 @@ function rhs!(du, u, t, mesh::TreeMesh{1}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} + + element_range = eachelement(dg, cache) + interface_range = eachinterface(dg, cache) + boundary_range = eachboundary(dg, cache) + # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_range) # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(du, u, mesh, have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache) + dg.volume_integral, dg, cache, element_range) end # Prolong solution to interfaces @trixi_timeit timer() "prolong2interfaces" begin prolong2interfaces!(cache, u, mesh, equations, - dg.surface_integral, dg) + dg.surface_integral, dg, interface_range) end # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, interface_range) end # Prolong solution to boundaries @trixi_timeit timer() "prolong2boundaries" begin prolong2boundaries!(cache, u, mesh, equations, - dg.surface_integral, dg) + dg.surface_integral, dg, boundary_range) end # Calculate boundary fluxes @@ -116,15 +121,16 @@ function rhs!(du, u, t, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_surface_integral!(du, u, mesh, equations, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, element_range) end # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache) + @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, + element_range) # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range) end return nothing @@ -134,8 +140,8 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, - dg::DGSEM, cache) - @threaded for element in eachelement(dg, cache) + dg::DGSEM, cache, element_range) + @threaded for element in element_range weak_form_kernel!(du, u, element, mesh, nonconservative_terms, equations, dg, cache) @@ -176,8 +182,8 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, - dg::DGSEM, cache) - @threaded for element in eachelement(dg, cache) + dg::DGSEM, cache, element_range) + @threaded for element in element_range flux_differencing_kernel!(du, u, element, mesh, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache) @@ -255,7 +261,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, - dg::DGSEM, cache) + dg::DGSEM, cache, element_range) @unpack element_ids_dg, element_ids_dgfv = cache @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral @@ -264,7 +270,8 @@ function calc_volume_integral!(du, u, cache) # Determine element ids for DG-only and blended DG-FV volume integral - pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache) + pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache, + element_range) # Loop over pure DG elements @trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg) @@ -297,11 +304,11 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, - dg::DGSEM, cache) + dg::DGSEM, cache, element_range) @unpack volume_flux_fv = volume_integral # Calculate LGL FV volume integral - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv, dg, cache, element, true) end @@ -390,12 +397,13 @@ end # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, - mesh::TreeMesh{1}, equations, surface_integral, dg::DG) + mesh::TreeMesh{1}, equations, surface_integral, dg::DG, + interface_range) @unpack interfaces = cache @unpack neighbor_ids = interfaces interfaces_u = interfaces.u - @threaded for interface in eachinterface(dg, cache) + @threaded for interface in interface_range left_element = neighbor_ids[1, interface] right_element = neighbor_ids[2, interface] @@ -412,11 +420,12 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::False, equations, - surface_integral, dg::DG, cache) + surface_integral, dg::DG, cache, + interface_range) @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces - @threaded for interface in eachinterface(dg, cache) + @threaded for interface in interface_range # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -441,11 +450,11 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::True, equations, - surface_integral, dg::DG, cache) + surface_integral, dg::DG, cache, interface_range) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u, neighbor_ids, orientations = cache.interfaces - @threaded for interface in eachinterface(dg, cache) + @threaded for interface in interface_range # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -481,11 +490,12 @@ function calc_interface_flux!(surface_flux_values, end function prolong2boundaries!(cache, u, - mesh::TreeMesh{1}, equations, surface_integral, dg::DG) + mesh::TreeMesh{1}, equations, surface_integral, dg::DG, + boundary_range) @unpack boundaries = cache @unpack neighbor_sides = boundaries - @threaded for boundary in eachboundary(dg, cache) + @threaded for boundary in boundary_range element = boundaries.neighbor_ids[boundary] # boundary in x-direction @@ -603,7 +613,7 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A end function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - equations, surface_integral, dg::DGSEM, cache) + equations, surface_integral, dg::DGSEM, cache, element_range) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements @@ -613,7 +623,7 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1 # into FMAs (see comment at the top of the file). factor_1 = boundary_interpolation[1, 1] factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range for v in eachvariable(equations) # surface at -x du[v, 1, element] = (du[v, 1, element] - @@ -629,10 +639,10 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1 end function apply_jacobian!(du, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - equations, dg::DG, cache) + equations, dg::DG, cache, element_range) @unpack inverse_jacobian = cache.elements - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range factor = -inverse_jacobian[element] for i in eachnode(dg) @@ -647,15 +657,15 @@ end # TODO: Taal dimension agnostic function calc_sources!(du, u, t, source_terms::Nothing, - equations::AbstractEquations{1}, dg::DG, cache) + equations::AbstractEquations{1}, dg::DG, cache, elements) return nothing end function calc_sources!(du, u, t, source_terms, - equations::AbstractEquations{1}, dg::DG, cache) + equations::AbstractEquations{1}, dg::DG, cache, element_range) @unpack node_coordinates = cache.elements - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range for i in eachnode(dg) u_local = get_node_vars(u, equations, dg, i, element) x_local = get_node_coords(node_coordinates, equations, dg, diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 0017f9ca88e..5c68eb94022 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -20,22 +20,29 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, @unpack viscous_container = cache_parabolic @unpack u_transformed, gradients, flux_viscous = viscous_container + element_range = eachelement(dg, cache) + interface_range = eachinterface(dg, cache) + boundary_range = eachboundary(dg, cache) + # Convert conservative variables to a form more suitable for viscous flux calculations @trixi_timeit timer() "transform variables" begin transform_variables!(u_transformed, u, mesh, equations_parabolic, - dg, parabolic_scheme, cache, cache_parabolic) + dg, parabolic_scheme, cache, cache_parabolic, + element_range) end # Compute the gradients of the transformed variables @trixi_timeit timer() "calculate gradient" begin calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, - boundary_conditions_parabolic, dg, cache, cache_parabolic) + boundary_conditions_parabolic, dg, cache, cache_parabolic, + element_range, interface_range, boundary_range) end # Compute and store the viscous fluxes @trixi_timeit timer() "calculate viscous fluxes" begin calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh, - equations_parabolic, dg, cache, cache_parabolic) + equations_parabolic, dg, cache, cache_parabolic, + element_range) end # The remainder of this function is essentially a regular rhs! for @@ -53,29 +60,29 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # TODO: parabolic; reconsider current data structure reuse strategy # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_range) # Calculate volume integral @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache) + calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache, element_range) end # Prolong solution to interfaces @trixi_timeit timer() "prolong2interfaces" begin prolong2interfaces!(cache_parabolic, flux_viscous, mesh, equations_parabolic, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, interface_range) end # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache_parabolic.elements.surface_flux_values, mesh, - equations_parabolic, dg, cache_parabolic) + equations_parabolic, dg, cache_parabolic, interface_range) end # Prolong solution to boundaries @trixi_timeit timer() "prolong2boundaries" begin prolong2boundaries!(cache_parabolic, flux_viscous, mesh, equations_parabolic, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, boundary_range) end # Calculate boundary fluxes @@ -89,12 +96,12 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_surface_integral!(du, u, mesh, equations_parabolic, - dg.surface_integral, dg, cache_parabolic) + dg.surface_integral, dg, cache_parabolic, element_range) end # Apply Jacobian from mapping to reference element @trixi_timeit timer() "Jacobian" begin - apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic) + apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic, element_range) end return nothing @@ -105,10 +112,11 @@ end # TODO: can we avoid copying data? function transform_variables!(u_transformed, u, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DG, parabolic_scheme, cache, cache_parabolic) + dg::DG, parabolic_scheme, cache, cache_parabolic, + element_range) transformation = gradient_variable_transformation(equations_parabolic) - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range # Calculate volume terms in one element for i in eachnode(dg) u_node = get_node_vars(u, equations_parabolic, dg, i, element) @@ -123,10 +131,10 @@ end function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache) + dg::DGSEM, cache, element_range) @unpack derivative_dhat = dg.basis - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range # Calculate volume terms in one element for i in eachnode(dg) flux_1_node = get_node_vars(flux_viscous, equations_parabolic, dg, i, @@ -147,12 +155,12 @@ end function prolong2interfaces!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) + surface_integral, dg::DG, cache, interface_range) @unpack interfaces = cache_parabolic @unpack neighbor_ids = interfaces interfaces_u = interfaces.u - @threaded for interface in eachinterface(dg, cache) + @threaded for interface in interface_range left_element = neighbor_ids[1, interface] right_element = neighbor_ids[2, interface] @@ -170,10 +178,10 @@ end # This is the version used when calculating the divergence of the viscous fluxes function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, equations_parabolic, - dg::DG, cache_parabolic) + dg::DG, cache_parabolic, interface_range) @unpack neighbor_ids, orientations = cache_parabolic.interfaces - @threaded for interface in eachinterface(dg, cache_parabolic) + @threaded for interface in interface_range # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -206,12 +214,12 @@ end function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) + surface_integral, dg::DG, cache, boundary_range) @unpack boundaries = cache_parabolic @unpack neighbor_sides, neighbor_ids = boundaries boundaries_u = boundaries.u - @threaded for boundary in eachboundary(dg, cache_parabolic) + @threaded for boundary in boundary_range element = neighbor_ids[boundary] if neighbor_sides[boundary] == 1 @@ -233,8 +241,8 @@ end function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache, cache_parabolic) - @threaded for element in eachelement(dg, cache) + dg::DG, cache, cache_parabolic, element_range) + @threaded for element in element_range for i in eachnode(dg) # Get solution and gradients u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, element) @@ -405,17 +413,18 @@ end # Calculate the gradient of the transformed variables function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, equations_parabolic, - boundary_conditions_parabolic, dg::DG, cache, cache_parabolic) + boundary_conditions_parabolic, dg::DG, cache, cache_parabolic, + element_range, interface_range, boundary_range) # Reset du @trixi_timeit timer() "reset gradients" begin - reset_du!(gradients, dg, cache) + reset_du!(gradients, dg, cache, element_range) end # Calculate volume integral @trixi_timeit timer() "volume integral" begin @unpack derivative_dhat = dg.basis - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range # Calculate volume terms in one element for i in eachnode(dg) @@ -436,14 +445,14 @@ function calc_gradient!(gradients, u_transformed, t, u_transformed, mesh, equations_parabolic, dg.surface_integral, - dg) + dg, interface_range) # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin @unpack surface_flux_values = cache_parabolic.elements @unpack neighbor_ids, orientations = cache_parabolic.interfaces - @threaded for interface in eachinterface(dg, cache_parabolic) + @threaded for interface in interface_range # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -471,7 +480,7 @@ function calc_gradient!(gradients, u_transformed, t, u_transformed, mesh, equations_parabolic, dg.surface_integral, - dg) + dg, boundary_range) # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!(cache_parabolic, @@ -493,7 +502,7 @@ function calc_gradient!(gradients, u_transformed, t, # into FMAs (see comment at the top of the file). factor_1 = boundary_interpolation[1, 1] factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range for v in eachvariable(equations_parabolic) # surface at -x gradients[v, 1, element] = (gradients[v, 1, element] - @@ -512,7 +521,7 @@ function calc_gradient!(gradients, u_transformed, t, # Apply Jacobian from mapping to reference element @trixi_timeit timer() "Jacobian" begin apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg, - cache_parabolic) + cache_parabolic, element_range) end return nothing @@ -549,10 +558,11 @@ end # `du/dt + df/dx = dg/dx + source(x,t)`, # where f(u) is the inviscid flux and g(u) is the viscous flux. function apply_jacobian_parabolic!(du, mesh::TreeMesh{1}, - equations::AbstractEquationsParabolic, dg::DG, cache) + equations::AbstractEquationsParabolic, dg::DG, cache, + element_range) @unpack inverse_jacobian = cache.elements - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range factor = inverse_jacobian[element] for i in eachnode(dg) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 547ed352ef3..ca210018d92 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -344,7 +344,8 @@ function calc_volume_integral!(du, u, cache) # Determine element ids for DG-only and blended DG-FV volume integral - pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache) + pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache, + eachelement(dg, cache)) # Loop over pure DG elements @trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg) diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 02ff338e912..b02be68ad91 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -394,7 +394,8 @@ function calc_volume_integral!(du, u, cache) # Determine element ids for DG-only and blended DG-FV volume integral - pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache) + pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache, + eachelement(dg, cache)) # Loop over pure DG elements @trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg) diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl index 0de0cff4851..3ecc0d9c1b1 100644 --- a/src/solvers/fdsbp_tree/fdsbp_1d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_1d.jl @@ -44,7 +44,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralStrongForm, - dg::FDSBP, cache) + dg::FDSBP, cache, element_range) D = dg.basis # SBP derivative operator @unpack f_threaded = cache @@ -67,7 +67,7 @@ function calc_volume_integral!(du, u, # Use the tensor product structure to compute the discrete derivatives of # the fluxes line-by-line and add them to `du` for each element. - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range f_element = f_threaded[Threads.threadid()] u_element = view(u_vectors, :, element) @@ -91,7 +91,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralUpwind, - dg::FDSBP, cache) + dg::FDSBP, cache, element_range) # Assume that # dg.basis isa SummationByPartsOperators.UpwindOperators D_minus = dg.basis.minus # Upwind SBP D^- derivative operator @@ -118,7 +118,7 @@ function calc_volume_integral!(du, u, # Use the tensor product structure to compute the discrete derivatives of # the fluxes line-by-line and add them to `du` for each element. - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range # f_minus_plus_element wraps the storage provided by f_minus_element and # f_plus_element such that we can use a single plain broadcasting below. # f_minus_element and f_plus_element are updated in broadcasting calls @@ -141,12 +141,12 @@ end function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralStrongForm, - dg::DG, cache) + dg::DG, cache, element_range) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range # surface at -x u_node = get_node_vars(u, equations, dg, 1, element) f_node = flux(u_node, 1, equations) @@ -182,11 +182,11 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::False, equations, surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache) + dg::FDSBP, cache, interface_range) @unpack splitting = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces - @threaded for interface in eachinterface(dg, cache) + @threaded for interface in interface_range # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -222,13 +222,13 @@ end # side of the element are computed in the upwind direction. function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache) + dg::FDSBP, cache, element_range) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements @unpack splitting = surface_integral - @threaded for element in eachelement(dg, cache) + @threaded for element in element_range # surface at -x u_node = get_node_vars(u, equations, dg, 1, element) f_node = splitting(u_node, Val{:plus}(), 1, equations) @@ -250,7 +250,7 @@ end # Periodic FDSBP operators need to use a single element without boundaries function calc_surface_integral!(du, u, mesh::TreeMesh1D, equations, surface_integral::SurfaceIntegralUpwind, - dg::PeriodicFDSBP, cache) + dg::PeriodicFDSBP, cache, element_range) @assert nelements(dg, cache) == 1 return nothing end From 20acb37f7748fe41fb42b882554e6918906b308e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 21 Mar 2024 12:45:53 +0100 Subject: [PATCH 02/24] partitioned rhs structured & tree --- src/solvers/dgsem_structured/dg_1d.jl | 50 ++++++++++++++++- src/solvers/dgsem_tree/dg_1d.jl | 67 ++++++++++++++++++++++- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 8 ++- 3 files changed, 119 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index abdd0fa82cb..615204307b8 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -9,7 +9,6 @@ function rhs!(du, u, t, mesh::StructuredMesh{1}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} - element_range = eachelement(dg, cache) # Reset du @@ -52,6 +51,55 @@ function rhs!(du, u, t, return nothing end +# RHS for PERK integrator +function rhs!(du, u, t, + mesh::StructuredMesh{1}, equations, + initial_condition, boundary_conditions, source_terms::Source, + dg::DG, cache, + element_range, + # Interfaces, boundaries, boundary_orientations, mortars not present for structured mesh + _, _, _, _) where {Source} + + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_range) + + # Calculate volume integral + @trixi_timeit timer() "volume integral" begin + calc_volume_integral!(du, u, mesh, + have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache, element_range) + end + + # Calculate interface and boundary fluxes + @trixi_timeit timer() "interface flux" begin + calc_interface_flux!(cache, u, mesh, equations, dg.surface_integral, dg, + element_range) + end + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" begin + calc_boundary_flux!(cache, u, t, boundary_conditions, mesh, equations, + dg.surface_integral, dg) + end + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" begin + calc_surface_integral!(du, u, mesh, equations, + dg.surface_integral, dg, cache, element_range) + end + + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, + element_range) + + # Calculate source terms + @trixi_timeit timer() "source terms" begin + calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range) + end + + return nothing +end + function calc_interface_flux!(cache, u, mesh::StructuredMesh{1}, equations, surface_integral, dg::DG, element_range) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 114da5211ea..d87a07beb9d 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -78,7 +78,6 @@ function rhs!(du, u, t, mesh::TreeMesh{1}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} - element_range = eachelement(dg, cache) interface_range = eachinterface(dg, cache) boundary_range = eachboundary(dg, cache) @@ -136,6 +135,69 @@ function rhs!(du, u, t, return nothing end +# RHS for PERK integrator +function rhs!(du, u, t, + mesh::TreeMesh{1}, equations, + initial_condition, boundary_conditions, source_terms::Source, + dg::DG, cache, + element_range, interface_range, boundary_range, + boundary_orientation_range, + # Mortars not appearant in 1D + _) where {Source} + + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, element_range) + + # Calculate volume integral + @trixi_timeit timer() "volume integral" begin + calc_volume_integral!(du, u, mesh, + have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache, element_range) + end + + # Prolong solution to interfaces, i.e., reconstruct interface/trace values + @trixi_timeit timer() "prolong2interfaces" begin + prolong2interfaces!(cache, u, mesh, equations, + dg.surface_integral, dg, interface_range) + end + + # Calculate interface fluxes + @trixi_timeit timer() "interface flux" begin + calc_interface_flux!(cache.elements.surface_flux_values, mesh, + have_nonconservative_terms(equations), equations, + dg.surface_integral, dg, cache, interface_range) + end + + # Prolong solution to boundaries + @trixi_timeit timer() "prolong2boundaries" begin + prolong2boundaries!(cache, u, mesh, equations, + dg.surface_integral, dg, boundary_range) + end + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" begin + calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations, + dg.surface_integral, dg, boundary_orientation_range) + end + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" begin + calc_surface_integral!(du, u, mesh, equations, + dg.surface_integral, dg, cache, element_range) + end + + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, + element_range) + + # Calculate source terms + @trixi_timeit timer() "source terms" begin + calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range) + end + + return nothing +end + function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, @@ -613,7 +675,8 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A end function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - equations, surface_integral, dg::DGSEM, cache, element_range) + equations, surface_integral, dg::DGSEM, cache, + element_range) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 5c68eb94022..74330ffc0ad 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -64,7 +64,8 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # Calculate volume integral @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache, element_range) + calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache, + element_range) end # Prolong solution to interfaces @@ -101,7 +102,8 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # Apply Jacobian from mapping to reference element @trixi_timeit timer() "Jacobian" begin - apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic, element_range) + apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic, + element_range) end return nothing @@ -112,7 +114,7 @@ end # TODO: can we avoid copying data? function transform_variables!(u_transformed, u, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DG, parabolic_scheme, cache, cache_parabolic, + dg::DG, parabolic_scheme, cache, cache_parabolic, element_range) transformation = gradient_variable_transformation(equations_parabolic) From e56855b98f51e215ac5c46b171297395b65219e7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 5 Jun 2024 09:21:05 +0200 Subject: [PATCH 03/24] remove PERK Multi rhs, make range argument optional --- src/solvers/dgsem_structured/dg_1d.jl | 67 ++------------ src/solvers/dgsem_tree/dg.jl | 12 +-- src/solvers/dgsem_tree/dg_1d.jl | 108 ++++------------------ src/solvers/dgsem_tree/dg_1d_parabolic.jl | 46 ++++----- src/solvers/fdsbp_tree/fdsbp_1d.jl | 12 +-- 5 files changed, 55 insertions(+), 190 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index 615204307b8..d115a7c6fd0 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -9,71 +9,19 @@ function rhs!(du, u, t, mesh::StructuredMesh{1}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} - element_range = eachelement(dg, cache) - - # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_range) - - # Calculate volume integral - @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, u, mesh, - have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache, element_range) - end - - # Calculate interface and boundary fluxes - @trixi_timeit timer() "interface flux" begin - calc_interface_flux!(cache, u, mesh, equations, dg.surface_integral, dg, - element_range) - end - - # Calculate boundary fluxes - @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux!(cache, u, t, boundary_conditions, mesh, equations, - dg.surface_integral, dg) - end - - # Calculate surface integrals - @trixi_timeit timer() "surface integral" begin - calc_surface_integral!(du, u, mesh, equations, - dg.surface_integral, dg, cache, element_range) - end - - # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, - element_range) - - # Calculate source terms - @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range) - end - - return nothing -end - -# RHS for PERK integrator -function rhs!(du, u, t, - mesh::StructuredMesh{1}, equations, - initial_condition, boundary_conditions, source_terms::Source, - dg::DG, cache, - element_range, - # Interfaces, boundaries, boundary_orientations, mortars not present for structured mesh - _, _, _, _) where {Source} - # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_range) + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(du, u, mesh, have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache, element_range) + dg.volume_integral, dg, cache) end # Calculate interface and boundary fluxes @trixi_timeit timer() "interface flux" begin - calc_interface_flux!(cache, u, mesh, equations, dg.surface_integral, dg, - element_range) + calc_interface_flux!(cache, u, mesh, equations, dg.surface_integral, dg) end # Calculate boundary fluxes @@ -85,16 +33,15 @@ function rhs!(du, u, t, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_surface_integral!(du, u, mesh, equations, - dg.surface_integral, dg, cache, element_range) + dg.surface_integral, dg, cache) end # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, - element_range) + @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache) # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range) + calc_sources!(du, u, t, source_terms, equations, dg, cache) end return nothing @@ -102,7 +49,7 @@ end function calc_interface_flux!(cache, u, mesh::StructuredMesh{1}, equations, surface_integral, dg::DG, - element_range) + element_range=eachelement(dg, cache)) @unpack surface_flux = surface_integral @threaded for element in element_range diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index 3b22d40234b..2ec83365ae0 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -7,15 +7,7 @@ # du .= zero(eltype(du)) doesn't scale when using multiple threads. # See https://github.com/trixi-framework/Trixi.jl/pull/924 for a performance comparison. -function reset_du!(du, dg, cache) - @threaded for element in eachelement(dg, cache) - du[.., element] .= zero(eltype(du)) - end - - return du -end - -function reset_du!(du, dg, cache, element_range) +function reset_du!(du, dg, cache, element_range=eachelement(dg, cache)) @threaded for element in element_range du[.., element] .= zero(eltype(du)) end @@ -29,7 +21,7 @@ end # `element_ids_dg` with the IDs of elements using a pure DG scheme and # `element_ids_dgfv` with the IDs of elements using a blended DG-FV scheme. function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg::DG, - cache, element_range) + cache, element_range=eachelement(dg, cache)) empty!(element_ids_dg) empty!(element_ids_dgfv) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index d87a07beb9d..c06cdd8fca3 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -78,121 +78,53 @@ function rhs!(du, u, t, mesh::TreeMesh{1}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} - element_range = eachelement(dg, cache) - interface_range = eachinterface(dg, cache) - boundary_range = eachboundary(dg, cache) - # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_range) + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(du, u, mesh, have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache, element_range) + dg.volume_integral, dg, cache) end # Prolong solution to interfaces @trixi_timeit timer() "prolong2interfaces" begin prolong2interfaces!(cache, u, mesh, equations, - dg.surface_integral, dg, interface_range) - end - - # Calculate interface fluxes - @trixi_timeit timer() "interface flux" begin - calc_interface_flux!(cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, - dg.surface_integral, dg, cache, interface_range) - end - - # Prolong solution to boundaries - @trixi_timeit timer() "prolong2boundaries" begin - prolong2boundaries!(cache, u, mesh, equations, - dg.surface_integral, dg, boundary_range) - end - - # Calculate boundary fluxes - @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) end - # Calculate surface integrals - @trixi_timeit timer() "surface integral" begin - calc_surface_integral!(du, u, mesh, equations, - dg.surface_integral, dg, cache, element_range) - end - - # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, - element_range) - - # Calculate source terms - @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range) - end - - return nothing -end - -# RHS for PERK integrator -function rhs!(du, u, t, - mesh::TreeMesh{1}, equations, - initial_condition, boundary_conditions, source_terms::Source, - dg::DG, cache, - element_range, interface_range, boundary_range, - boundary_orientation_range, - # Mortars not appearant in 1D - _) where {Source} - - # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, element_range) - - # Calculate volume integral - @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, u, mesh, - have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache, element_range) - end - - # Prolong solution to interfaces, i.e., reconstruct interface/trace values - @trixi_timeit timer() "prolong2interfaces" begin - prolong2interfaces!(cache, u, mesh, equations, - dg.surface_integral, dg, interface_range) - end - # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg.surface_integral, dg, cache, interface_range) + dg.surface_integral, dg, cache) end # Prolong solution to boundaries @trixi_timeit timer() "prolong2boundaries" begin prolong2boundaries!(cache, u, mesh, equations, - dg.surface_integral, dg, boundary_range) + dg.surface_integral, dg) end # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" begin calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations, - dg.surface_integral, dg, boundary_orientation_range) + dg.surface_integral, dg) end # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_surface_integral!(du, u, mesh, equations, - dg.surface_integral, dg, cache, element_range) + dg.surface_integral, dg, cache) end # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, - element_range) + @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache) # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range) + calc_sources!(du, u, t, source_terms, equations, dg, cache) end return nothing @@ -202,7 +134,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, - dg::DGSEM, cache, element_range) + dg::DGSEM, cache, element_range=eachelement(dg, cache)) @threaded for element in element_range weak_form_kernel!(du, u, element, mesh, nonconservative_terms, equations, @@ -244,7 +176,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, - dg::DGSEM, cache, element_range) + dg::DGSEM, cache, element_range=eachelement(dg, cache)) @threaded for element in element_range flux_differencing_kernel!(du, u, element, mesh, nonconservative_terms, equations, @@ -323,7 +255,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, - dg::DGSEM, cache, element_range) + dg::DGSEM, cache, element_range=eachelement(dg, cache)) @unpack element_ids_dg, element_ids_dgfv = cache @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral @@ -366,7 +298,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, - dg::DGSEM, cache, element_range) + dg::DGSEM, cache, element_range=eachelement(dg, cache)) @unpack volume_flux_fv = volume_integral # Calculate LGL FV volume integral @@ -460,7 +392,7 @@ end # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::TreeMesh{1}, equations, surface_integral, dg::DG, - interface_range) + interface_range=eachinterface(dg, cache)) @unpack interfaces = cache @unpack neighbor_ids = interfaces interfaces_u = interfaces.u @@ -483,7 +415,7 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, - interface_range) + interface_range=eachinterface(dg, cache)) @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @@ -512,7 +444,7 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::True, equations, - surface_integral, dg::DG, cache, interface_range) + surface_integral, dg::DG, cache, interface_range=eachinterface(dg, cache)) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u, neighbor_ids, orientations = cache.interfaces @@ -553,7 +485,7 @@ end function prolong2boundaries!(cache, u, mesh::TreeMesh{1}, equations, surface_integral, dg::DG, - boundary_range) + boundary_range=eachboundary(dg, cache)) @unpack boundaries = cache @unpack neighbor_sides = boundaries @@ -676,7 +608,7 @@ end function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations, surface_integral, dg::DGSEM, cache, - element_range) + element_range=eachelement(dg, cache)) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements @@ -702,7 +634,7 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1 end function apply_jacobian!(du, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - equations, dg::DG, cache, element_range) + equations, dg::DG, cache, element_range=eachelement(dg, cache)) @unpack inverse_jacobian = cache.elements @threaded for element in element_range @@ -720,12 +652,12 @@ end # TODO: Taal dimension agnostic function calc_sources!(du, u, t, source_terms::Nothing, - equations::AbstractEquations{1}, dg::DG, cache, elements) + equations::AbstractEquations{1}, dg::DG, cache, element_range=eachelement(dg, cache)) return nothing end function calc_sources!(du, u, t, source_terms, - equations::AbstractEquations{1}, dg::DG, cache, element_range) + equations::AbstractEquations{1}, dg::DG, cache, element_range=eachelement(dg, cache)) @unpack node_coordinates = cache.elements @threaded for element in element_range diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 74330ffc0ad..e06dffc0b10 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -20,29 +20,23 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, @unpack viscous_container = cache_parabolic @unpack u_transformed, gradients, flux_viscous = viscous_container - element_range = eachelement(dg, cache) - interface_range = eachinterface(dg, cache) - boundary_range = eachboundary(dg, cache) - # Convert conservative variables to a form more suitable for viscous flux calculations @trixi_timeit timer() "transform variables" begin transform_variables!(u_transformed, u, mesh, equations_parabolic, - dg, parabolic_scheme, cache, cache_parabolic, - element_range) + dg, parabolic_scheme, cache, cache_parabolic) end # Compute the gradients of the transformed variables @trixi_timeit timer() "calculate gradient" begin calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, boundary_conditions_parabolic, dg, cache, cache_parabolic, - element_range, interface_range, boundary_range) + element_range, interface_range) end # Compute and store the viscous fluxes @trixi_timeit timer() "calculate viscous fluxes" begin calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh, - equations_parabolic, dg, cache, cache_parabolic, - element_range) + equations_parabolic, dg, cache, cache_parabolic) end # The remainder of this function is essentially a regular rhs! for @@ -60,30 +54,29 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # TODO: parabolic; reconsider current data structure reuse strategy # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_range) + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) # Calculate volume integral @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache, - element_range) + calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache) end # Prolong solution to interfaces @trixi_timeit timer() "prolong2interfaces" begin prolong2interfaces!(cache_parabolic, flux_viscous, mesh, equations_parabolic, - dg.surface_integral, dg, cache, interface_range) + dg.surface_integral, dg, cache) end # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache_parabolic.elements.surface_flux_values, mesh, - equations_parabolic, dg, cache_parabolic, interface_range) + equations_parabolic, dg, cache_parabolic) end # Prolong solution to boundaries @trixi_timeit timer() "prolong2boundaries" begin prolong2boundaries!(cache_parabolic, flux_viscous, mesh, equations_parabolic, - dg.surface_integral, dg, cache, boundary_range) + dg.surface_integral, dg, cache) end # Calculate boundary fluxes @@ -97,13 +90,12 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_surface_integral!(du, u, mesh, equations_parabolic, - dg.surface_integral, dg, cache_parabolic, element_range) + dg.surface_integral, dg, cache_parabolic) end # Apply Jacobian from mapping to reference element @trixi_timeit timer() "Jacobian" begin - apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic, - element_range) + apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic) end return nothing @@ -115,7 +107,7 @@ end function transform_variables!(u_transformed, u, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, dg::DG, parabolic_scheme, cache, cache_parabolic, - element_range) + element_range=eachelement(dg, cache)) transformation = gradient_variable_transformation(equations_parabolic) @threaded for element in element_range @@ -133,7 +125,7 @@ end function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache, element_range) + dg::DGSEM, cache, element_range=eachelement(dg, cache)) @unpack derivative_dhat = dg.basis @threaded for element in element_range @@ -157,7 +149,7 @@ end function prolong2interfaces!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache, interface_range) + surface_integral, dg::DG, cache, interface_range=eachinterface(dg, cache)) @unpack interfaces = cache_parabolic @unpack neighbor_ids = interfaces interfaces_u = interfaces.u @@ -180,7 +172,7 @@ end # This is the version used when calculating the divergence of the viscous fluxes function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, equations_parabolic, - dg::DG, cache_parabolic, interface_range) + dg::DG, cache_parabolic, interface_range=eachinterface(dg, cache)) @unpack neighbor_ids, orientations = cache_parabolic.interfaces @threaded for interface in interface_range @@ -216,7 +208,7 @@ end function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache, boundary_range) + surface_integral, dg::DG, cache, boundary_range=eachboundary(dg, cache)) @unpack boundaries = cache_parabolic @unpack neighbor_sides, neighbor_ids = boundaries boundaries_u = boundaries.u @@ -243,7 +235,7 @@ end function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache, cache_parabolic, element_range) + dg::DG, cache, cache_parabolic, element_range=eachelement(dg, cache)) @threaded for element in element_range for i in eachnode(dg) # Get solution and gradients @@ -416,7 +408,9 @@ end function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, equations_parabolic, boundary_conditions_parabolic, dg::DG, cache, cache_parabolic, - element_range, interface_range, boundary_range) + element_range=eachelement(dg, cache) + interface_range=eachinterface(dg, cache), + boundary_range=eachboundary(dg, cache)) # Reset du @trixi_timeit timer() "reset gradients" begin @@ -561,7 +555,7 @@ end # where f(u) is the inviscid flux and g(u) is the viscous flux. function apply_jacobian_parabolic!(du, mesh::TreeMesh{1}, equations::AbstractEquationsParabolic, dg::DG, cache, - element_range) + element_range=eachelement(dg, cache)) @unpack inverse_jacobian = cache.elements @threaded for element in element_range diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl index 3ecc0d9c1b1..829e51b29ca 100644 --- a/src/solvers/fdsbp_tree/fdsbp_1d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_1d.jl @@ -44,7 +44,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralStrongForm, - dg::FDSBP, cache, element_range) + dg::FDSBP, cache, element_range=eachelement(dg, cache)) D = dg.basis # SBP derivative operator @unpack f_threaded = cache @@ -91,7 +91,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralUpwind, - dg::FDSBP, cache, element_range) + dg::FDSBP, cache, element_range=eachelement(dg, cache)) # Assume that # dg.basis isa SummationByPartsOperators.UpwindOperators D_minus = dg.basis.minus # Upwind SBP D^- derivative operator @@ -141,7 +141,7 @@ end function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralStrongForm, - dg::DG, cache, element_range) + dg::DG, cache, element_range=eachelement(dg, cache)) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements @@ -182,7 +182,7 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::False, equations, surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache, interface_range) + dg::FDSBP, cache, interface_range=eachinterface(dg, cache)) @unpack splitting = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @@ -222,7 +222,7 @@ end # side of the element are computed in the upwind direction. function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache, element_range) + dg::FDSBP, cache, element_range=eachelement(dg, cache)) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements @@ -250,7 +250,7 @@ end # Periodic FDSBP operators need to use a single element without boundaries function calc_surface_integral!(du, u, mesh::TreeMesh1D, equations, surface_integral::SurfaceIntegralUpwind, - dg::PeriodicFDSBP, cache, element_range) + dg::PeriodicFDSBP, cache, element_range=eachelement(dg, cache)) @assert nelements(dg, cache) == 1 return nothing end From ea127785c9f268d2f522144dfcb69800aea503eb Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 5 Jun 2024 09:25:15 +0200 Subject: [PATCH 04/24] remove changes for 2d 3d SC --- src/solvers/dgsem_tree/dg_2d.jl | 3 +-- src/solvers/dgsem_tree/dg_3d.jl | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index ca210018d92..547ed352ef3 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -344,8 +344,7 @@ function calc_volume_integral!(du, u, cache) # Determine element ids for DG-only and blended DG-FV volume integral - pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache, - eachelement(dg, cache)) + pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache) # Loop over pure DG elements @trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg) diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index b02be68ad91..02ff338e912 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -394,8 +394,7 @@ function calc_volume_integral!(du, u, cache) # Determine element ids for DG-only and blended DG-FV volume integral - pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache, - eachelement(dg, cache)) + pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache) # Loop over pure DG elements @trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg) From cf3813de02c1adde74a5e811cc326668b4d9208d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 5 Jun 2024 09:26:19 +0200 Subject: [PATCH 05/24] fmt --- src/solvers/dgsem_structured/dg_1d.jl | 2 +- src/solvers/dgsem_tree/dg.jl | 4 ++-- src/solvers/dgsem_tree/dg_1d.jl | 28 +++++++++++++---------- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 24 +++++++++++-------- src/solvers/fdsbp_tree/fdsbp_1d.jl | 15 +++++++----- 5 files changed, 42 insertions(+), 31 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index d115a7c6fd0..e06e218f15d 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -49,7 +49,7 @@ end function calc_interface_flux!(cache, u, mesh::StructuredMesh{1}, equations, surface_integral, dg::DG, - element_range=eachelement(dg, cache)) + element_range = eachelement(dg, cache)) @unpack surface_flux = surface_integral @threaded for element in element_range diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index 2ec83365ae0..201e425183b 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -7,7 +7,7 @@ # du .= zero(eltype(du)) doesn't scale when using multiple threads. # See https://github.com/trixi-framework/Trixi.jl/pull/924 for a performance comparison. -function reset_du!(du, dg, cache, element_range=eachelement(dg, cache)) +function reset_du!(du, dg, cache, element_range = eachelement(dg, cache)) @threaded for element in element_range du[.., element] .= zero(eltype(du)) end @@ -21,7 +21,7 @@ end # `element_ids_dg` with the IDs of elements using a pure DG scheme and # `element_ids_dgfv` with the IDs of elements using a blended DG-FV scheme. function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg::DG, - cache, element_range=eachelement(dg, cache)) + cache, element_range = eachelement(dg, cache)) empty!(element_ids_dg) empty!(element_ids_dgfv) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index c06cdd8fca3..f72e81634a3 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -134,7 +134,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, - dg::DGSEM, cache, element_range=eachelement(dg, cache)) + dg::DGSEM, cache, element_range = eachelement(dg, cache)) @threaded for element in element_range weak_form_kernel!(du, u, element, mesh, nonconservative_terms, equations, @@ -176,7 +176,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, - dg::DGSEM, cache, element_range=eachelement(dg, cache)) + dg::DGSEM, cache, element_range = eachelement(dg, cache)) @threaded for element in element_range flux_differencing_kernel!(du, u, element, mesh, nonconservative_terms, equations, @@ -255,7 +255,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, - dg::DGSEM, cache, element_range=eachelement(dg, cache)) + dg::DGSEM, cache, element_range = eachelement(dg, cache)) @unpack element_ids_dg, element_ids_dgfv = cache @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral @@ -298,7 +298,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, - dg::DGSEM, cache, element_range=eachelement(dg, cache)) + dg::DGSEM, cache, element_range = eachelement(dg, cache)) @unpack volume_flux_fv = volume_integral # Calculate LGL FV volume integral @@ -392,7 +392,7 @@ end # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::TreeMesh{1}, equations, surface_integral, dg::DG, - interface_range=eachinterface(dg, cache)) + interface_range = eachinterface(dg, cache)) @unpack interfaces = cache @unpack neighbor_ids = interfaces interfaces_u = interfaces.u @@ -415,7 +415,7 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, - interface_range=eachinterface(dg, cache)) + interface_range = eachinterface(dg, cache)) @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @@ -444,7 +444,8 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::True, equations, - surface_integral, dg::DG, cache, interface_range=eachinterface(dg, cache)) + surface_integral, dg::DG, cache, + interface_range = eachinterface(dg, cache)) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u, neighbor_ids, orientations = cache.interfaces @@ -485,7 +486,7 @@ end function prolong2boundaries!(cache, u, mesh::TreeMesh{1}, equations, surface_integral, dg::DG, - boundary_range=eachboundary(dg, cache)) + boundary_range = eachboundary(dg, cache)) @unpack boundaries = cache @unpack neighbor_sides = boundaries @@ -608,7 +609,7 @@ end function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations, surface_integral, dg::DGSEM, cache, - element_range=eachelement(dg, cache)) + element_range = eachelement(dg, cache)) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements @@ -634,7 +635,8 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1 end function apply_jacobian!(du, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - equations, dg::DG, cache, element_range=eachelement(dg, cache)) + equations, dg::DG, cache, + element_range = eachelement(dg, cache)) @unpack inverse_jacobian = cache.elements @threaded for element in element_range @@ -652,12 +654,14 @@ end # TODO: Taal dimension agnostic function calc_sources!(du, u, t, source_terms::Nothing, - equations::AbstractEquations{1}, dg::DG, cache, element_range=eachelement(dg, cache)) + equations::AbstractEquations{1}, dg::DG, cache, + element_range = eachelement(dg, cache)) return nothing end function calc_sources!(du, u, t, source_terms, - equations::AbstractEquations{1}, dg::DG, cache, element_range=eachelement(dg, cache)) + equations::AbstractEquations{1}, dg::DG, cache, + element_range = eachelement(dg, cache)) @unpack node_coordinates = cache.elements @threaded for element in element_range diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index e06dffc0b10..7ba9ae5ceb8 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -107,7 +107,7 @@ end function transform_variables!(u_transformed, u, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, dg::DG, parabolic_scheme, cache, cache_parabolic, - element_range=eachelement(dg, cache)) + element_range = eachelement(dg, cache)) transformation = gradient_variable_transformation(equations_parabolic) @threaded for element in element_range @@ -125,7 +125,7 @@ end function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache, element_range=eachelement(dg, cache)) + dg::DGSEM, cache, element_range = eachelement(dg, cache)) @unpack derivative_dhat = dg.basis @threaded for element in element_range @@ -149,7 +149,8 @@ end function prolong2interfaces!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache, interface_range=eachinterface(dg, cache)) + surface_integral, dg::DG, cache, + interface_range = eachinterface(dg, cache)) @unpack interfaces = cache_parabolic @unpack neighbor_ids = interfaces interfaces_u = interfaces.u @@ -172,7 +173,8 @@ end # This is the version used when calculating the divergence of the viscous fluxes function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, equations_parabolic, - dg::DG, cache_parabolic, interface_range=eachinterface(dg, cache)) + dg::DG, cache_parabolic, + interface_range = eachinterface(dg, cache)) @unpack neighbor_ids, orientations = cache_parabolic.interfaces @threaded for interface in interface_range @@ -208,7 +210,8 @@ end function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache, boundary_range=eachboundary(dg, cache)) + surface_integral, dg::DG, cache, + boundary_range = eachboundary(dg, cache)) @unpack boundaries = cache_parabolic @unpack neighbor_sides, neighbor_ids = boundaries boundaries_u = boundaries.u @@ -235,7 +238,8 @@ end function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache, cache_parabolic, element_range=eachelement(dg, cache)) + dg::DG, cache, cache_parabolic, + element_range = eachelement(dg, cache)) @threaded for element in element_range for i in eachnode(dg) # Get solution and gradients @@ -408,9 +412,9 @@ end function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, equations_parabolic, boundary_conditions_parabolic, dg::DG, cache, cache_parabolic, - element_range=eachelement(dg, cache) - interface_range=eachinterface(dg, cache), - boundary_range=eachboundary(dg, cache)) + element_range = eachelement(dg, cache) + interface_range = eachinterface(dg, cache), + boundary_range = eachboundary(dg, cache)) # Reset du @trixi_timeit timer() "reset gradients" begin @@ -555,7 +559,7 @@ end # where f(u) is the inviscid flux and g(u) is the viscous flux. function apply_jacobian_parabolic!(du, mesh::TreeMesh{1}, equations::AbstractEquationsParabolic, dg::DG, cache, - element_range=eachelement(dg, cache)) + element_range = eachelement(dg, cache)) @unpack inverse_jacobian = cache.elements @threaded for element in element_range diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl index 829e51b29ca..4df910fcef4 100644 --- a/src/solvers/fdsbp_tree/fdsbp_1d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_1d.jl @@ -44,7 +44,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralStrongForm, - dg::FDSBP, cache, element_range=eachelement(dg, cache)) + dg::FDSBP, cache, element_range = eachelement(dg, cache)) D = dg.basis # SBP derivative operator @unpack f_threaded = cache @@ -91,7 +91,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralUpwind, - dg::FDSBP, cache, element_range=eachelement(dg, cache)) + dg::FDSBP, cache, element_range = eachelement(dg, cache)) # Assume that # dg.basis isa SummationByPartsOperators.UpwindOperators D_minus = dg.basis.minus # Upwind SBP D^- derivative operator @@ -141,7 +141,7 @@ end function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralStrongForm, - dg::DG, cache, element_range=eachelement(dg, cache)) + dg::DG, cache, element_range = eachelement(dg, cache)) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements @@ -182,7 +182,8 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::False, equations, surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache, interface_range=eachinterface(dg, cache)) + dg::FDSBP, cache, + interface_range = eachinterface(dg, cache)) @unpack splitting = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @@ -222,7 +223,8 @@ end # side of the element are computed in the upwind direction. function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache, element_range=eachelement(dg, cache)) + dg::FDSBP, cache, + element_range = eachelement(dg, cache)) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements @@ -250,7 +252,8 @@ end # Periodic FDSBP operators need to use a single element without boundaries function calc_surface_integral!(du, u, mesh::TreeMesh1D, equations, surface_integral::SurfaceIntegralUpwind, - dg::PeriodicFDSBP, cache, element_range=eachelement(dg, cache)) + dg::PeriodicFDSBP, cache, + element_range = eachelement(dg, cache)) @assert nelements(dg, cache) == 1 return nothing end From 749c67004dd8c5f52738344200477699fc09fd08 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 5 Jun 2024 09:36:41 +0200 Subject: [PATCH 06/24] typo --- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 7ba9ae5ceb8..9165d2f4540 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -412,7 +412,7 @@ end function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, equations_parabolic, boundary_conditions_parabolic, dg::DG, cache, cache_parabolic, - element_range = eachelement(dg, cache) + element_range = eachelement(dg, cache), interface_range = eachinterface(dg, cache), boundary_range = eachboundary(dg, cache)) From 24c8c78b02e62425f8dcd57b01be72ec0b290938 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 5 Jun 2024 10:25:33 +0200 Subject: [PATCH 07/24] remove --- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 9165d2f4540..d059104419c 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -29,8 +29,7 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # Compute the gradients of the transformed variables @trixi_timeit timer() "calculate gradient" begin calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, - boundary_conditions_parabolic, dg, cache, cache_parabolic, - element_range, interface_range) + boundary_conditions_parabolic, dg, cache, cache_parabolic) end # Compute and store the viscous fluxes From 947687ecc8d4d97ea3372d6a3a34bbb1e69c75bb Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 5 Jun 2024 11:09:11 +0200 Subject: [PATCH 08/24] debug --- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index d059104419c..704f2583f16 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -444,7 +444,8 @@ function calc_gradient!(gradients, u_transformed, t, u_transformed, mesh, equations_parabolic, dg.surface_integral, - dg, interface_range) + dg, cache, + interface_range) # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin @@ -479,7 +480,8 @@ function calc_gradient!(gradients, u_transformed, t, u_transformed, mesh, equations_parabolic, dg.surface_integral, - dg, boundary_range) + dg, cache, + boundary_range) # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!(cache_parabolic, From 0e76dca1df918a6a10f4a7324222609c6c7eb7df Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 5 Jun 2024 11:12:51 +0200 Subject: [PATCH 09/24] fmt --- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 704f2583f16..35a0f8a533e 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -69,7 +69,8 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache_parabolic.elements.surface_flux_values, mesh, - equations_parabolic, dg, cache_parabolic) + equations_parabolic, dg, cache_parabolic, + eachinterface(dg, cache)) end # Prolong solution to boundaries @@ -173,7 +174,7 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, equations_parabolic, dg::DG, cache_parabolic, - interface_range = eachinterface(dg, cache)) + interface_range) @unpack neighbor_ids, orientations = cache_parabolic.interfaces @threaded for interface in interface_range From 622cc8e7f189ce29b949328fc923343b0d9d3b92 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 30 Oct 2024 11:07:14 +0100 Subject: [PATCH 10/24] "range" to "indices" --- src/solvers/dgsem_structured/dg_1d.jl | 4 +- src/solvers/dgsem_tree/dg.jl | 8 ++-- src/solvers/dgsem_tree/dg_1d.jl | 46 +++++++++++----------- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 48 +++++++++++------------ src/solvers/fdsbp_tree/fdsbp_1d.jl | 22 +++++------ 5 files changed, 64 insertions(+), 64 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index e06e218f15d..5f5573cb32a 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -49,10 +49,10 @@ end function calc_interface_flux!(cache, u, mesh::StructuredMesh{1}, equations, surface_integral, dg::DG, - element_range = eachelement(dg, cache)) + element_indices = eachelement(dg, cache)) @unpack surface_flux = surface_integral - @threaded for element in element_range + @threaded for element in element_indices left_element = cache.elements.left_neighbors[1, element] if left_element > 0 # left_element = 0 at boundaries diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index 3b2447d8ee6..2d5c2af28b4 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -7,8 +7,8 @@ # du .= zero(eltype(du)) doesn't scale when using multiple threads. # See https://github.com/trixi-framework/Trixi.jl/pull/924 for a performance comparison. -function reset_du!(du, dg, cache, element_range = eachelement(dg, cache)) - @threaded for element in element_range +function reset_du!(du, dg, cache, element_indices = eachelement(dg, cache)) + @threaded for element in element_indices du[.., element] .= zero(eltype(du)) end @@ -21,11 +21,11 @@ end # `element_ids_dg` with the IDs of elements using a pure DG scheme and # `element_ids_dgfv` with the IDs of elements using a blended DG-FV scheme. function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg::DG, - cache, element_range = eachelement(dg, cache)) + cache, element_indices = eachelement(dg, cache)) empty!(element_ids_dg) empty!(element_ids_dgfv) - for element in element_range + for element in element_indices # Clip blending factor for values close to zero (-> pure DG) dg_only = isapprox(alpha[element], 0, atol = 1e-12) if dg_only diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index f72e81634a3..c2d1d6908ba 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -134,8 +134,8 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, - dg::DGSEM, cache, element_range = eachelement(dg, cache)) - @threaded for element in element_range + dg::DGSEM, cache, element_indices = eachelement(dg, cache)) + @threaded for element in element_indices weak_form_kernel!(du, u, element, mesh, nonconservative_terms, equations, dg, cache) @@ -176,8 +176,8 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, - dg::DGSEM, cache, element_range = eachelement(dg, cache)) - @threaded for element in element_range + dg::DGSEM, cache, element_indices = eachelement(dg, cache)) + @threaded for element in element_indices flux_differencing_kernel!(du, u, element, mesh, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache) @@ -255,7 +255,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, - dg::DGSEM, cache, element_range = eachelement(dg, cache)) + dg::DGSEM, cache, element_indices = eachelement(dg, cache)) @unpack element_ids_dg, element_ids_dgfv = cache @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral @@ -265,7 +265,7 @@ function calc_volume_integral!(du, u, # Determine element ids for DG-only and blended DG-FV volume integral pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache, - element_range) + element_indices) # Loop over pure DG elements @trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg) @@ -298,11 +298,11 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, - dg::DGSEM, cache, element_range = eachelement(dg, cache)) + dg::DGSEM, cache, element_indices = eachelement(dg, cache)) @unpack volume_flux_fv = volume_integral # Calculate LGL FV volume integral - @threaded for element in element_range + @threaded for element in element_indices fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv, dg, cache, element, true) end @@ -392,12 +392,12 @@ end # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::TreeMesh{1}, equations, surface_integral, dg::DG, - interface_range = eachinterface(dg, cache)) + interface_indices = eachinterface(dg, cache)) @unpack interfaces = cache @unpack neighbor_ids = interfaces interfaces_u = interfaces.u - @threaded for interface in interface_range + @threaded for interface in interface_indices left_element = neighbor_ids[1, interface] right_element = neighbor_ids[2, interface] @@ -415,11 +415,11 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, - interface_range = eachinterface(dg, cache)) + interface_indices = eachinterface(dg, cache)) @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces - @threaded for interface in interface_range + @threaded for interface in interface_indices # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -445,11 +445,11 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, - interface_range = eachinterface(dg, cache)) + interface_indices = eachinterface(dg, cache)) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u, neighbor_ids, orientations = cache.interfaces - @threaded for interface in interface_range + @threaded for interface in interface_indices # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -486,11 +486,11 @@ end function prolong2boundaries!(cache, u, mesh::TreeMesh{1}, equations, surface_integral, dg::DG, - boundary_range = eachboundary(dg, cache)) + boundary_indices = eachboundary(dg, cache)) @unpack boundaries = cache @unpack neighbor_sides = boundaries - @threaded for boundary in boundary_range + @threaded for boundary in boundary_indices element = boundaries.neighbor_ids[boundary] # boundary in x-direction @@ -609,7 +609,7 @@ end function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations, surface_integral, dg::DGSEM, cache, - element_range = eachelement(dg, cache)) + element_indices = eachelement(dg, cache)) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements @@ -619,7 +619,7 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1 # into FMAs (see comment at the top of the file). factor_1 = boundary_interpolation[1, 1] factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in element_range + @threaded for element in element_indices for v in eachvariable(equations) # surface at -x du[v, 1, element] = (du[v, 1, element] - @@ -636,10 +636,10 @@ end function apply_jacobian!(du, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations, dg::DG, cache, - element_range = eachelement(dg, cache)) + element_indices = eachelement(dg, cache)) @unpack inverse_jacobian = cache.elements - @threaded for element in element_range + @threaded for element in element_indices factor = -inverse_jacobian[element] for i in eachnode(dg) @@ -655,16 +655,16 @@ end # TODO: Taal dimension agnostic function calc_sources!(du, u, t, source_terms::Nothing, equations::AbstractEquations{1}, dg::DG, cache, - element_range = eachelement(dg, cache)) + element_indices = eachelement(dg, cache)) return nothing end function calc_sources!(du, u, t, source_terms, equations::AbstractEquations{1}, dg::DG, cache, - element_range = eachelement(dg, cache)) + element_indices = eachelement(dg, cache)) @unpack node_coordinates = cache.elements - @threaded for element in element_range + @threaded for element in element_indices for i in eachnode(dg) u_local = get_node_vars(u, equations, dg, i, element) x_local = get_node_coords(node_coordinates, equations, dg, diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 35a0f8a533e..985f6e8b6c0 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -107,10 +107,10 @@ end function transform_variables!(u_transformed, u, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, dg::DG, parabolic_scheme, cache, cache_parabolic, - element_range = eachelement(dg, cache)) + element_indices = eachelement(dg, cache)) transformation = gradient_variable_transformation(equations_parabolic) - @threaded for element in element_range + @threaded for element in element_indices # Calculate volume terms in one element for i in eachnode(dg) u_node = get_node_vars(u, equations_parabolic, dg, i, element) @@ -125,10 +125,10 @@ end function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache, element_range = eachelement(dg, cache)) + dg::DGSEM, cache, element_indices = eachelement(dg, cache)) @unpack derivative_dhat = dg.basis - @threaded for element in element_range + @threaded for element in element_indices # Calculate volume terms in one element for i in eachnode(dg) flux_1_node = get_node_vars(flux_viscous, equations_parabolic, dg, i, @@ -150,12 +150,12 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache, - interface_range = eachinterface(dg, cache)) + interface_indices = eachinterface(dg, cache)) @unpack interfaces = cache_parabolic @unpack neighbor_ids = interfaces interfaces_u = interfaces.u - @threaded for interface in interface_range + @threaded for interface in interface_indices left_element = neighbor_ids[1, interface] right_element = neighbor_ids[2, interface] @@ -174,10 +174,10 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, equations_parabolic, dg::DG, cache_parabolic, - interface_range) + interface_indices) @unpack neighbor_ids, orientations = cache_parabolic.interfaces - @threaded for interface in interface_range + @threaded for interface in interface_indices # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -211,12 +211,12 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache, - boundary_range = eachboundary(dg, cache)) + boundary_indices = eachboundary(dg, cache)) @unpack boundaries = cache_parabolic @unpack neighbor_sides, neighbor_ids = boundaries boundaries_u = boundaries.u - @threaded for boundary in boundary_range + @threaded for boundary in boundary_indices element = neighbor_ids[boundary] if neighbor_sides[boundary] == 1 @@ -239,8 +239,8 @@ end function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache, cache_parabolic, - element_range = eachelement(dg, cache)) - @threaded for element in element_range + element_indices = eachelement(dg, cache)) + @threaded for element in element_indices for i in eachnode(dg) # Get solution and gradients u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, element) @@ -412,19 +412,19 @@ end function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, equations_parabolic, boundary_conditions_parabolic, dg::DG, cache, cache_parabolic, - element_range = eachelement(dg, cache), - interface_range = eachinterface(dg, cache), - boundary_range = eachboundary(dg, cache)) + element_indices = eachelement(dg, cache), + interface_indices = eachinterface(dg, cache), + boundary_indices = eachboundary(dg, cache)) # Reset du @trixi_timeit timer() "reset gradients" begin - reset_du!(gradients, dg, cache, element_range) + reset_du!(gradients, dg, cache, element_indices) end # Calculate volume integral @trixi_timeit timer() "volume integral" begin @unpack derivative_dhat = dg.basis - @threaded for element in element_range + @threaded for element in element_indices # Calculate volume terms in one element for i in eachnode(dg) @@ -446,14 +446,14 @@ function calc_gradient!(gradients, u_transformed, t, equations_parabolic, dg.surface_integral, dg, cache, - interface_range) + interface_indices) # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin @unpack surface_flux_values = cache_parabolic.elements @unpack neighbor_ids, orientations = cache_parabolic.interfaces - @threaded for interface in interface_range + @threaded for interface in interface_indices # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -482,7 +482,7 @@ function calc_gradient!(gradients, u_transformed, t, equations_parabolic, dg.surface_integral, dg, cache, - boundary_range) + boundary_indices) # Calculate boundary fluxes @trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!(cache_parabolic, @@ -504,7 +504,7 @@ function calc_gradient!(gradients, u_transformed, t, # into FMAs (see comment at the top of the file). factor_1 = boundary_interpolation[1, 1] factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in element_range + @threaded for element in element_indices for v in eachvariable(equations_parabolic) # surface at -x gradients[v, 1, element] = (gradients[v, 1, element] - @@ -523,7 +523,7 @@ function calc_gradient!(gradients, u_transformed, t, # Apply Jacobian from mapping to reference element @trixi_timeit timer() "Jacobian" begin apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg, - cache_parabolic, element_range) + cache_parabolic, element_indices) end return nothing @@ -561,10 +561,10 @@ end # where f(u) is the inviscid flux and g(u) is the viscous flux. function apply_jacobian_parabolic!(du, mesh::TreeMesh{1}, equations::AbstractEquationsParabolic, dg::DG, cache, - element_range = eachelement(dg, cache)) + element_indices = eachelement(dg, cache)) @unpack inverse_jacobian = cache.elements - @threaded for element in element_range + @threaded for element in element_indices factor = inverse_jacobian[element] for i in eachnode(dg) diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl index 4df910fcef4..345b2af8bbd 100644 --- a/src/solvers/fdsbp_tree/fdsbp_1d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_1d.jl @@ -44,7 +44,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralStrongForm, - dg::FDSBP, cache, element_range = eachelement(dg, cache)) + dg::FDSBP, cache, element_indices = eachelement(dg, cache)) D = dg.basis # SBP derivative operator @unpack f_threaded = cache @@ -67,7 +67,7 @@ function calc_volume_integral!(du, u, # Use the tensor product structure to compute the discrete derivatives of # the fluxes line-by-line and add them to `du` for each element. - @threaded for element in element_range + @threaded for element in element_indices f_element = f_threaded[Threads.threadid()] u_element = view(u_vectors, :, element) @@ -91,7 +91,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralUpwind, - dg::FDSBP, cache, element_range = eachelement(dg, cache)) + dg::FDSBP, cache, element_indices = eachelement(dg, cache)) # Assume that # dg.basis isa SummationByPartsOperators.UpwindOperators D_minus = dg.basis.minus # Upwind SBP D^- derivative operator @@ -118,7 +118,7 @@ function calc_volume_integral!(du, u, # Use the tensor product structure to compute the discrete derivatives of # the fluxes line-by-line and add them to `du` for each element. - @threaded for element in element_range + @threaded for element in element_indices # f_minus_plus_element wraps the storage provided by f_minus_element and # f_plus_element such that we can use a single plain broadcasting below. # f_minus_element and f_plus_element are updated in broadcasting calls @@ -141,12 +141,12 @@ end function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralStrongForm, - dg::DG, cache, element_range = eachelement(dg, cache)) + dg::DG, cache, element_indices = eachelement(dg, cache)) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements - @threaded for element in element_range + @threaded for element in element_indices # surface at -x u_node = get_node_vars(u, equations, dg, 1, element) f_node = flux(u_node, 1, equations) @@ -183,11 +183,11 @@ function calc_interface_flux!(surface_flux_values, nonconservative_terms::False, equations, surface_integral::SurfaceIntegralUpwind, dg::FDSBP, cache, - interface_range = eachinterface(dg, cache)) + interface_indices = eachinterface(dg, cache)) @unpack splitting = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces - @threaded for interface in interface_range + @threaded for interface in interface_indices # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -224,13 +224,13 @@ end function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralUpwind, dg::FDSBP, cache, - element_range = eachelement(dg, cache)) + element_indices = eachelement(dg, cache)) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements @unpack splitting = surface_integral - @threaded for element in element_range + @threaded for element in element_indices # surface at -x u_node = get_node_vars(u, equations, dg, 1, element) f_node = splitting(u_node, Val{:plus}(), 1, equations) @@ -253,7 +253,7 @@ end function calc_surface_integral!(du, u, mesh::TreeMesh1D, equations, surface_integral::SurfaceIntegralUpwind, dg::PeriodicFDSBP, cache, - element_range = eachelement(dg, cache)) + element_indices = eachelement(dg, cache)) @assert nelements(dg, cache) == 1 return nothing end From 5e17307c258c477fe9db71d2a25a5cc240290875 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 30 Oct 2024 11:29:32 +0100 Subject: [PATCH 11/24] fmt --- src/solvers/dgsem_tree/dg_1d.jl | 12 ++++++++---- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 3 ++- src/solvers/fdsbp_tree/fdsbp_1d.jl | 6 ++++-- 3 files changed, 14 insertions(+), 7 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index c2d1d6908ba..b88c905e9c4 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -134,7 +134,8 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, - dg::DGSEM, cache, element_indices = eachelement(dg, cache)) + dg::DGSEM, cache, + element_indices = eachelement(dg, cache)) @threaded for element in element_indices weak_form_kernel!(du, u, element, mesh, nonconservative_terms, equations, @@ -176,7 +177,8 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, - dg::DGSEM, cache, element_indices = eachelement(dg, cache)) + dg::DGSEM, cache, + element_indices = eachelement(dg, cache)) @threaded for element in element_indices flux_differencing_kernel!(du, u, element, mesh, nonconservative_terms, equations, @@ -255,7 +257,8 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, - dg::DGSEM, cache, element_indices = eachelement(dg, cache)) + dg::DGSEM, cache, + element_indices = eachelement(dg, cache)) @unpack element_ids_dg, element_ids_dgfv = cache @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral @@ -298,7 +301,8 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, - dg::DGSEM, cache, element_indices = eachelement(dg, cache)) + dg::DGSEM, cache, + element_indices = eachelement(dg, cache)) @unpack volume_flux_fv = volume_integral # Calculate LGL FV volume integral diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 985f6e8b6c0..03b75cdeb38 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -125,7 +125,8 @@ end function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache, element_indices = eachelement(dg, cache)) + dg::DGSEM, cache, + element_indices = eachelement(dg, cache)) @unpack derivative_dhat = dg.basis @threaded for element in element_indices diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl index 345b2af8bbd..8b006d7bcbe 100644 --- a/src/solvers/fdsbp_tree/fdsbp_1d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_1d.jl @@ -44,7 +44,8 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralStrongForm, - dg::FDSBP, cache, element_indices = eachelement(dg, cache)) + dg::FDSBP, cache, + element_indices = eachelement(dg, cache)) D = dg.basis # SBP derivative operator @unpack f_threaded = cache @@ -91,7 +92,8 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::False, equations, volume_integral::VolumeIntegralUpwind, - dg::FDSBP, cache, element_indices = eachelement(dg, cache)) + dg::FDSBP, cache, + element_indices = eachelement(dg, cache)) # Assume that # dg.basis isa SummationByPartsOperators.UpwindOperators D_minus = dg.basis.minus # Upwind SBP D^- derivative operator From f3f1cadcee54705b6c2b28e01f3694ccece62b0e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 31 Oct 2024 15:07:48 +0100 Subject: [PATCH 12/24] specialize ambiguous funcs --- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 03b75cdeb38..31a85e0fc13 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -150,7 +150,7 @@ end function prolong2interfaces!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache, + surface_integral, dg::DG, cache::NamedTuple, interface_indices = eachinterface(dg, cache)) @unpack interfaces = cache_parabolic @unpack neighbor_ids = interfaces @@ -211,7 +211,7 @@ end function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache, + surface_integral, dg::DG, cache::NamedTuple, boundary_indices = eachboundary(dg, cache)) @unpack boundaries = cache_parabolic @unpack neighbor_sides, neighbor_ids = boundaries From 347f09b0f007e73c53756b1d8be5c0e75d05c3a8 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 4 Nov 2024 09:23:23 +0100 Subject: [PATCH 13/24] structured 1D rhs with indices --- src/solvers/dgsem_structured/dg_1d.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index bc05c3c52c6..f3ac09b616a 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -8,20 +8,22 @@ function rhs!(du, u, t, mesh::StructuredMesh{1}, equations, boundary_conditions, source_terms::Source, - dg::DG, cache) where {Source} + dg::DG, cache, + element_indices = eachelement(dg, cache)) where {Source} # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_indices) # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(du, u, mesh, have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache) + dg.volume_integral, dg, cache, element_indices) end # Calculate interface and boundary fluxes @trixi_timeit timer() "interface flux" begin - calc_interface_flux!(cache, u, mesh, equations, dg.surface_integral, dg) + calc_interface_flux!(cache, u, mesh, equations, dg.surface_integral, dg, + element_indices) end # Calculate boundary fluxes @@ -33,15 +35,17 @@ function rhs!(du, u, t, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_surface_integral!(du, u, mesh, equations, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, + element_indices) end # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache) + @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, + element_indices) # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(du, u, t, source_terms, equations, dg, cache, element_indices) end return nothing From 90237472cfcd67f1efdae2060e38d4d6343b9fd2 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 4 Nov 2024 09:32:26 +0100 Subject: [PATCH 14/24] 1d tree rhs with indices --- src/solvers/dgsem_structured/dg_1d.jl | 5 ++++- src/solvers/dgsem_tree/dg_1d.jl | 23 ++++++++++++++--------- 2 files changed, 18 insertions(+), 10 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index f3ac09b616a..cbbc87cd8ae 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -9,7 +9,10 @@ function rhs!(du, u, t, mesh::StructuredMesh{1}, equations, boundary_conditions, source_terms::Source, dg::DG, cache, - element_indices = eachelement(dg, cache)) where {Source} + element_indices = eachelement(dg, cache), + interface_indices = nothing, + boundary_indices = nothing, + mortar_indices = nothing) where {Source} # Reset du @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_indices) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 53e2948dce8..ed20a829ed6 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -77,34 +77,38 @@ end function rhs!(du, u, t, mesh::TreeMesh{1}, equations, boundary_conditions, source_terms::Source, - dg::DG, cache) where {Source} + dg::DG, cache, + element_indices = eachelement(dg, cache), + interface_indices = eachinterface(dg, cache), + boundary_indices = eachboundary(dg, cache), + mortar_indices = nothing) where {Source} # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_indices) # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(du, u, mesh, have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache) + dg.volume_integral, dg, cache, element_indices) end # Prolong solution to interfaces @trixi_timeit timer() "prolong2interfaces" begin prolong2interfaces!(cache, u, mesh, equations, - dg.surface_integral, dg) + dg.surface_integral, dg, interface_indices) end # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, interface_indices) end # Prolong solution to boundaries @trixi_timeit timer() "prolong2boundaries" begin prolong2boundaries!(cache, u, mesh, equations, - dg.surface_integral, dg) + dg.surface_integral, dg, boundary_indices) end # Calculate boundary fluxes @@ -116,15 +120,16 @@ function rhs!(du, u, t, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_surface_integral!(du, u, mesh, equations, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, element_indices) end # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache) + @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache, + element_indices) # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(du, u, t, source_terms, equations, dg, cache, element_indices) end return nothing From f468f46e7b5fcaaf797133d0e812021c952a33d7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 4 Nov 2024 09:36:39 +0100 Subject: [PATCH 15/24] SC 1D indices --- src/solvers/dgsem_tree/dg_1d.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index df3c72e51b3..67fe50e72eb 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -258,6 +258,8 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, + dg::DGSEM, cache, + element_indices = eachelement(dg, cache)) @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral # Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α @@ -267,7 +269,7 @@ function calc_volume_integral!(du, u, # For `Float32`, this gives 1.1920929f-5 RealT = eltype(alpha) atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0)) - @threaded for element in eachelement(dg, cache) + @threaded for element in element_indices alpha_element = alpha[element] # Clip blending factor for values close to zero (-> pure DG) dg_only = isapprox(alpha_element, 0, atol = atol) From 40dc503a7557af49e11402bfffa0ef0d56f3e933 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 4 Nov 2024 10:05:01 +0100 Subject: [PATCH 16/24] parabolic 1D --- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 32 +++++++++++++++-------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 8815bbb947d..d6695eaa212 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -16,26 +16,33 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, boundary_conditions_parabolic, source_terms, - dg::DG, parabolic_scheme, cache, cache_parabolic) + dg::DG, parabolic_scheme, cache, cache_parabolic, + element_indices = eachelement(dg, cache), + interface_indices = eachinterface(dg, cache), + boundary_indices = eachboundary(dg, cache), + mortar_indices = nothing) @unpack viscous_container = cache_parabolic @unpack u_transformed, gradients, flux_viscous = viscous_container # Convert conservative variables to a form more suitable for viscous flux calculations @trixi_timeit timer() "transform variables" begin transform_variables!(u_transformed, u, mesh, equations_parabolic, - dg, parabolic_scheme, cache, cache_parabolic) + dg, parabolic_scheme, cache, cache_parabolic, + element_indices) end # Compute the gradients of the transformed variables @trixi_timeit timer() "calculate gradient" begin calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, - boundary_conditions_parabolic, dg, cache, cache_parabolic) + boundary_conditions_parabolic, dg, cache, cache_parabolic, + element_indices, interface_indices, boundary_indices) end # Compute and store the viscous fluxes @trixi_timeit timer() "calculate viscous fluxes" begin calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh, - equations_parabolic, dg, cache, cache_parabolic) + equations_parabolic, dg, cache, cache_parabolic, + element_indices) end # The remainder of this function is essentially a regular rhs! for @@ -53,30 +60,31 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # TODO: parabolic; reconsider current data structure reuse strategy # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_indices) # Calculate volume integral @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache) + calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache, + element_indices) end # Prolong solution to interfaces @trixi_timeit timer() "prolong2interfaces" begin prolong2interfaces!(cache_parabolic, flux_viscous, mesh, equations_parabolic, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, interface_indices) end # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache_parabolic.elements.surface_flux_values, mesh, equations_parabolic, dg, cache_parabolic, - eachinterface(dg, cache)) + interface_indices) end # Prolong solution to boundaries @trixi_timeit timer() "prolong2boundaries" begin prolong2boundaries!(cache_parabolic, flux_viscous, mesh, equations_parabolic, - dg.surface_integral, dg, cache) + dg.surface_integral, dg, cache, boundary_indices) end # Calculate boundary fluxes @@ -90,12 +98,14 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_surface_integral!(du, u, mesh, equations_parabolic, - dg.surface_integral, dg, cache_parabolic) + dg.surface_integral, dg, cache_parabolic, + element_indices) end # Apply Jacobian from mapping to reference element @trixi_timeit timer() "Jacobian" begin - apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic) + apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic, + element_indices) end return nothing From 64b5e036e0ba448a5090d8d1ef1cd69001e34a18 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 4 Nov 2024 10:18:12 +0100 Subject: [PATCH 17/24] consistency --- src/solvers/fdsbp_tree/fdsbp_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl index 8b006d7bcbe..6e1894600b0 100644 --- a/src/solvers/fdsbp_tree/fdsbp_1d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_1d.jl @@ -170,7 +170,7 @@ end # Periodic FDSBP operators need to use a single element without boundaries function calc_surface_integral!(du, u, mesh::TreeMesh1D, equations, surface_integral::SurfaceIntegralStrongForm, - dg::PeriodicFDSBP, cache) + dg::PeriodicFDSBP, cache, element_indices = nothing) @assert nelements(dg, cache) == 1 return nothing end From 0619874b3ad1373f33add59aea9b828d30b01b23 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 19 Jul 2025 15:55:21 -0400 Subject: [PATCH 18/24] local conservation --- src/solvers/dgsem_tree/dg.jl | 1 + src/solvers/dgsem_tree/dg_1d.jl | 2 ++ 2 files changed, 3 insertions(+) diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index b051f432ed3..76568a1971d 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -8,6 +8,7 @@ # du .= zero(eltype(du)) doesn't scale when using multiple threads. # See https://github.com/trixi-framework/Trixi.jl/pull/924 for a performance comparison. function reset_du!(du, dg, cache, element_indices = eachelement(dg, cache)) + # TODO: Reset surface flux values here? @threaded for element in element_indices du[.., element] .= zero(eltype(du)) end diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 67fe50e72eb..86349960ce7 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -436,6 +436,8 @@ function calc_interface_flux!(surface_flux_values, # Copy flux to left and right element storage for v in eachvariable(equations) + # TODO: Accumulate fluxes, i.e., `+=`? + # Reset would then be handled by `reset_du!` surface_flux_values[v, left_direction, left_id] = flux[v] surface_flux_values[v, right_direction, right_id] = flux[v] end From 93a0e7e710065b95fbb5ea73c5fbf2a783b77511 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 31 Oct 2025 15:54:40 +0100 Subject: [PATCH 19/24] mrege --- src/solvers/solvers.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index 9eb7bdfa0b2..8f26814446b 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -6,10 +6,10 @@ #! format: noindent # Used by both `dg::DGSEM` and `dg::FDSBP` -function reset_du!(du, dg, cache) +function reset_du!(du, dg, cache, element_indices = eachelement(dg, cache)) # du .= zero(eltype(du)) doesn't scale when using multiple threads. # See https://github.com/trixi-framework/Trixi.jl/pull/924 for a performance comparison. - @threaded for element in eachelement(dg, cache) + @threaded for element in element_indices du[.., element] .= zero(eltype(du)) end From abddfe12deb7df5ac09fa42a9192cd13f54b8480 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 31 Oct 2025 15:58:54 +0100 Subject: [PATCH 20/24] fmt --- src/solvers/dgsem_tree/dg_1d.jl | 2 +- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 37 +++++------------------ 2 files changed, 9 insertions(+), 30 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 5337209f9d2..7cd03b0c845 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -415,7 +415,7 @@ end end function prolong2interfaces!(cache, u, mesh::TreeMesh{1}, equations, dg::DG, - interface_indices = eachinterface(dg, cache)) + interface_indices = eachinterface(dg, cache)) @unpack interfaces = cache @unpack neighbor_ids = interfaces interfaces_u = interfaces.u diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 8a2e33ca665..06b2be24400 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -456,10 +456,11 @@ end function calc_gradient_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, equations_parabolic, dg::DG, parabolic_scheme, - cache, cache_parabolic) + cache, cache_parabolic, + interface_indices = eachinterface(dg, cache)) @unpack neighbor_ids, orientations = cache_parabolic.interfaces - @threaded for interface in eachinterface(dg, cache_parabolic) + @threaded for interface in interface_indices # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -489,7 +490,8 @@ end function calc_gradient_surface_integral!(gradients, mesh::TreeMesh{1}, # for dispatch only equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache, cache_parabolic) + dg::DGSEM, cache, cache_parabolic, + element_indices = eachelement(dg, cache)) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache_parabolic.elements @@ -499,7 +501,7 @@ function calc_gradient_surface_integral!(gradients, # into FMAs (see comment at the top of the file). factor_1 = boundary_interpolation[1, 1] factor_2 = boundary_interpolation[nnodes(dg), 2] - @threaded for element in eachelement(dg, cache) + @threaded for element in element_indices for v in eachvariable(equations_parabolic) # surface at -x gradients[v, 1, element] = (gradients[v, 1, element] - @@ -544,32 +546,9 @@ function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, # Calculate interface fluxes @trixi_timeit timer() "interface flux" begin @unpack surface_flux_values = cache_parabolic.elements - @unpack neighbor_ids, orientations = cache_parabolic.interfaces - - @threaded for interface in interface_indices - # Get neighboring elements - left_id = neighbor_ids[1, interface] - right_id = neighbor_ids[2, interface] - - # Determine interface direction with respect to elements: - # orientation = 1: left -> 2, right -> 1 - left_direction = 2 * orientations[interface] - right_direction = 2 * orientations[interface] - 1 - - # Call pointwise Riemann solver - u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, - equations_parabolic, dg, interface) - flux = 0.5f0 * (u_ll + u_rr) - - # Copy flux to left and right element storage - for v in eachvariable(equations_parabolic) - surface_flux_values[v, left_direction, left_id] = flux[v] - surface_flux_values[v, right_direction, right_id] = flux[v] - end - end calc_gradient_interface_flux!(surface_flux_values, mesh, equations_parabolic, dg, parabolic_scheme, - cache, cache_parabolic) + cache, cache_parabolic, interface_indices) end # Prolong solution to boundaries @@ -592,7 +571,7 @@ function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1}, # Calculate surface integrals @trixi_timeit timer() "surface integral" begin calc_gradient_surface_integral!(gradients, mesh, equations_parabolic, - dg, cache, cache_parabolic) + dg, cache, cache_parabolic, element_indices) end # Apply Jacobian from mapping to reference element From 28fb4090d798872b4964a017fad33d441aeb6e1a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 31 Oct 2025 15:59:19 +0100 Subject: [PATCH 21/24] fix --- src/solvers/dgsem_structured/dg_1d.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index 9d94cc4b575..63c1599d963 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -7,7 +7,8 @@ function calc_interface_flux!(cache, u, mesh::StructuredMesh{1}, nonconservative_terms, # can be True/False - equations, surface_integral, dg::DG) + equations, surface_integral, dg::DG, + element_indices = eachelement(dg, cache)) @unpack surface_flux = surface_integral @threaded for element in element_indices @@ -70,10 +71,11 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, end function apply_jacobian!(du, mesh::StructuredMesh{1}, - equations, dg::DG, cache) + equations, dg::DG, cache, + element_indices = eachelement(dg, cache)) @unpack inverse_jacobian = cache.elements - @threaded for element in eachelement(dg, cache) + @threaded for element in element_indices for i in eachnode(dg) factor = -inverse_jacobian[i, element] From b3400e9ffbb4db0315ec22aa331575f2f163f470 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 31 Oct 2025 16:06:17 +0100 Subject: [PATCH 22/24] each --- src/solvers/dgsem_tree/dg_1d.jl | 2 +- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 7cd03b0c845..c5c6720e143 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -255,7 +255,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1} @unpack x_interfaces, volume_flux_fv, reconstruction_mode, slope_limiter = volume_integral # Calculate LGL second-order FV volume integral - @threaded for element in eachelement(dg, cache) + @threaded for element in element_indices fvO2_kernel!(du, u, mesh, have_nonconservative_terms, equations, volume_flux_fv, dg, cache, element, diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 06b2be24400..32e8e293209 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -432,10 +432,11 @@ end function calc_gradient_volume_integral!(gradients, u_transformed, mesh::TreeMesh{1}, # for dispatch only equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache) + dg::DGSEM, cache, + element_indices = eachelement(dg, cache)) @unpack derivative_dhat = dg.basis - @threaded for element in eachelement(dg, cache) + @threaded for element in element_indices # Calculate volume terms in one element, # corresponds to `kernel` functions for the hyperbolic part of the flux for i in eachnode(dg) From a6d936a01daba536ea5b8ae6fa4b88b86b75ded1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 31 Oct 2025 16:09:05 +0100 Subject: [PATCH 23/24] Tree only --- src/solvers/dgsem_structured/dg_1d.jl | 10 ++++----- src/solvers/fdsbp_tree/fdsbp_1d.jl | 29 +++++++++++---------------- 2 files changed, 16 insertions(+), 23 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index 63c1599d963..cb98c45aed3 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -7,11 +7,10 @@ function calc_interface_flux!(cache, u, mesh::StructuredMesh{1}, nonconservative_terms, # can be True/False - equations, surface_integral, dg::DG, - element_indices = eachelement(dg, cache)) + equations, surface_integral, dg::DG) @unpack surface_flux = surface_integral - @threaded for element in element_indices + @threaded for element in eachelement(dg, cache) left_element = cache.elements.left_neighbors[1, element] # => `element` is the right element of the interface @@ -71,11 +70,10 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, end function apply_jacobian!(du, mesh::StructuredMesh{1}, - equations, dg::DG, cache, - element_indices = eachelement(dg, cache)) + equations, dg::DG, cache) @unpack inverse_jacobian = cache.elements - @threaded for element in element_indices + @threaded for element in eachelement(dg, cache) for i in eachnode(dg) factor = -inverse_jacobian[i, element] diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl index ed7f1878fb0..93336c2f70d 100644 --- a/src/solvers/fdsbp_tree/fdsbp_1d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_1d.jl @@ -44,8 +44,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, have_nonconservative_terms::False, equations, volume_integral::VolumeIntegralStrongForm, - dg::FDSBP, cache, - element_indices = eachelement(dg, cache)) + dg::FDSBP, cache) D = dg.basis # SBP derivative operator @unpack f_threaded = cache @@ -68,7 +67,7 @@ function calc_volume_integral!(du, u, # Use the tensor product structure to compute the discrete derivatives of # the fluxes line-by-line and add them to `du` for each element. - @threaded for element in element_indices + @threaded for element in eachelement(dg, cache) f_element = f_threaded[Threads.threadid()] u_element = view(u_vectors, :, element) @@ -92,8 +91,7 @@ function calc_volume_integral!(du, u, mesh::TreeMesh{1}, have_nonconservative_terms::False, equations, volume_integral::VolumeIntegralUpwind, - dg::FDSBP, cache, - element_indices = eachelement(dg, cache)) + dg::FDSBP, cache) # Assume that # dg.basis isa SummationByPartsOperators.UpwindOperators D_minus = dg.basis.minus # Upwind SBP D^- derivative operator @@ -120,7 +118,7 @@ function calc_volume_integral!(du, u, # Use the tensor product structure to compute the discrete derivatives of # the fluxes line-by-line and add them to `du` for each element. - @threaded for element in element_indices + @threaded for element in eachelement(dg, cache) # f_minus_plus_element wraps the storage provided by f_minus_element and # f_plus_element such that we can use a single plain broadcasting below. # f_minus_element and f_plus_element are updated in broadcasting calls @@ -143,12 +141,12 @@ end function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralStrongForm, - dg::DG, cache, element_indices = eachelement(dg, cache)) + dg::DG, cache) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements - @threaded for element in element_indices + @threaded for element in eachelement(dg, cache) # surface at -x u_node = get_node_vars(u, equations, dg, 1, element) f_node = flux(u_node, 1, equations) @@ -170,7 +168,7 @@ end # Periodic FDSBP operators need to use a single element without boundaries function calc_surface_integral!(du, u, mesh::TreeMesh1D, equations, surface_integral::SurfaceIntegralStrongForm, - dg::PeriodicFDSBP, cache, element_indices = nothing) + dg::PeriodicFDSBP, cache) @assert nelements(dg, cache) == 1 return nothing end @@ -184,12 +182,11 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, have_nonconservative_terms::False, equations, surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache, - interface_indices = eachinterface(dg, cache)) + dg::FDSBP, cache) @unpack splitting = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces - @threaded for interface in interface_indices + @threaded for interface in eachinterface(dg, cache) # Get neighboring elements left_id = neighbor_ids[1, interface] right_id = neighbor_ids[2, interface] @@ -225,14 +222,13 @@ end # side of the element are computed in the upwind direction. function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache, - element_indices = eachelement(dg, cache)) + dg::FDSBP, cache) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack surface_flux_values = cache.elements @unpack splitting = surface_integral - @threaded for element in element_indices + @threaded for element in eachelement(dg, cache) # surface at -x u_node = get_node_vars(u, equations, dg, 1, element) f_node = splitting(u_node, Val{:plus}(), 1, equations) @@ -254,8 +250,7 @@ end # Periodic FDSBP operators need to use a single element without boundaries function calc_surface_integral!(du, u, mesh::TreeMesh1D, equations, surface_integral::SurfaceIntegralUpwind, - dg::PeriodicFDSBP, cache, - element_indices = eachelement(dg, cache)) + dg::PeriodicFDSBP, cache) @assert nelements(dg, cache) == 1 return nothing end From f0fb7be7db27bd7f2b71abbee6a32bc30bbab348 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 31 Oct 2025 16:18:38 +0100 Subject: [PATCH 24/24] semi rhs! --- .../semidiscretization_hyperbolic.jl | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 2a563c02229..b51b10feab7 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -409,4 +409,25 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t) return nothing end + +# `rhs!` for partitioned Runge-Kutta methods, such as the Paired Explicit Runge-Kutta (PERK) methods +function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t, + element_indices, interface_indices, boundary_indices, mortar_indices) + @unpack mesh, equations, boundary_conditions, source_terms, solver, cache = semi + + u = wrap_array(u_ode, mesh, equations, solver, cache) + du = wrap_array(du_ode, mesh, equations, solver, cache) + + # TODO: Taal decide, do we need to pass the mesh? + time_start = time_ns() + @trixi_timeit timer() "rhs! (part.)" rhs!(du, u, t, mesh, equations, + boundary_conditions, source_terms, + solver, cache, + element_indices, interface_indices, + boundary_indices, mortar_indices) + runtime = time_ns() - time_start + put!(semi.performance_counter, runtime) + + return nothing +end end # @muladd