diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl index 3ca76ad6aa0..c98cd55e6b5 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl @@ -1,5 +1,6 @@ using OrdinaryDiffEqSDIRK, ADTypes using Trixi +import LinearSolve: LUFactorization ############################################################################### # semidiscretization of the linear advection diffusion equation @@ -87,6 +88,8 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks time_int_tol = 1.0e-10 time_abs_tol = 1.0e-10 -sol = solve(ode, KenCarp4(autodiff = AutoFiniteDiff()); # This is an IMEX SDIRK method +# This is an IMEX SDIRK method +ode_alg = KenCarp4(autodiff = AutoFiniteDiff(), linsolve = LUFactorization()) +sol = solve(ode, ode_alg; abstol = time_abs_tol, reltol = time_int_tol, ode_default_options()..., callback = callbacks) diff --git a/src/solvers/dgsem/calc_volume_integral.jl b/src/solvers/dgsem/calc_volume_integral.jl new file mode 100644 index 00000000000..8fa83612a25 --- /dev/null +++ b/src/solvers/dgsem/calc_volume_integral.jl @@ -0,0 +1,95 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# Dimension and meshtype agnostic, i.e., valid for all 1D, 2D, and 3D meshes +function create_cache(mesh, equations, + volume_integral::VolumeIntegralFluxDifferencing, dg::DG, uEltype) + return NamedTuple() +end + +# The following `calc_volume_integral!` functions are +# dimension and meshtype agnostic, i.e., valid for all 1D, 2D, and 3D meshes. + +function calc_volume_integral!(du, u, mesh, + nonconservative_terms, equations, + volume_integral::VolumeIntegralWeakForm, + dg::DGSEM, cache) + @threaded for element in eachelement(dg, cache) + weak_form_kernel!(du, u, element, mesh, + nonconservative_terms, equations, + dg, cache) + end + + return nothing +end + +function calc_volume_integral!(du, u, mesh, + nonconservative_terms, equations, + volume_integral::VolumeIntegralFluxDifferencing, + dg::DGSEM, cache) + @threaded for element in eachelement(dg, cache) + flux_differencing_kernel!(du, u, element, mesh, nonconservative_terms, + equations, + volume_integral.volume_flux, dg, cache) + end + + return nothing +end + +function calc_volume_integral!(du, u, mesh, + nonconservative_terms, equations, + volume_integral::VolumeIntegralShockCapturingHG, + dg::DGSEM, cache) + @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral + + # Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α + alpha = @trixi_timeit timer() "blending factors" indicator(u, mesh, equations, dg, + cache) + + # For `Float64`, this gives 1.8189894035458565e-12 + # 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) + alpha_element = alpha[element] + # Clip blending factor for values close to zero (-> pure DG) + dg_only = isapprox(alpha_element, 0, atol = atol) + + if dg_only + flux_differencing_kernel!(du, u, element, mesh, + nonconservative_terms, equations, + volume_flux_dg, dg, cache) + else + # Calculate DG volume integral contribution + flux_differencing_kernel!(du, u, element, mesh, + nonconservative_terms, equations, + volume_flux_dg, dg, cache, 1 - alpha_element) + + # Calculate FV volume integral contribution + fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv, + dg, cache, element, alpha_element) + end + end + + return nothing +end + +function calc_volume_integral!(du, u, mesh, + nonconservative_terms, equations, + volume_integral::VolumeIntegralPureLGLFiniteVolume, + dg::DGSEM, cache) + @unpack volume_flux_fv = volume_integral + + # Calculate LGL FV volume integral + @threaded for element in eachelement(dg, cache) + fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv, + dg, cache, element, true) + end + + return nothing +end +end # @muladd diff --git a/src/solvers/dgsem/dgsem.jl b/src/solvers/dgsem/dgsem.jl index 7a0bece6cd7..979f2d0577d 100644 --- a/src/solvers/dgsem/dgsem.jl +++ b/src/solvers/dgsem/dgsem.jl @@ -131,4 +131,6 @@ end end return u_mean / total_volume # normalize with the total volume end + +include("calc_volume_integral.jl") end # @muladd diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 3380a2b5c66..557b5c3364f 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -28,6 +28,13 @@ end Val(ndims(contravariant_vectors) - 3))) end +# Dimension agnostic, i.e., valid for all 1D, 2D, and 3D `StructuredMesh`es. +function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, + mesh::StructuredMesh, equations, surface_integral, + dg::DG) + @assert isperiodic(mesh) +end + @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionPeriodic, diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index ee2832e66a8..c8fb266d88e 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -70,14 +70,6 @@ function calc_interface_flux!(cache, u, mesh::StructuredMesh{1}, return nothing end -# TODO: Taal dimension agnostic -function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, - mesh::StructuredMesh{1}, equations, surface_integral, - dg::DG) - @assert isperiodic(mesh) - return nothing -end - function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, mesh::StructuredMesh{1}, equations, surface_integral, dg::DG) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index b7afe36d4f3..9498f1112d5 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -128,6 +128,8 @@ end # pull the contravariant vectors and compute the average Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, element) + # average mapping terms in first coordinate direction, + # used as normal vector in the flux computation Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector @@ -144,6 +146,8 @@ end # pull the contravariant vectors and compute the average Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, element) + # average mapping terms in second coordinate direction, + # used as normal vector in the flux computation Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector @@ -195,6 +199,8 @@ end # pull the contravariant vectors and compute the average Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, element) + # average mapping terms in first coordinate direction, + # used as normal vector in the flux computation Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # Compute the contravariant nonconservative flux. fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_avg, @@ -209,6 +215,8 @@ end # pull the contravariant vectors and compute the average Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, element) + # average mapping terms in second coordinate direction, + # used as normal vector in the flux computation Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector @@ -554,14 +562,6 @@ end return nothing end -# TODO: Taal dimension agnostic -function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, - mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, - equations, surface_integral, dg::DG) - @assert isperiodic(mesh) - return nothing -end - function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, equations, surface_integral, diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl index 70f4d17db10..86fc0862444 100644 --- a/src/solvers/dgsem_structured/dg_3d.jl +++ b/src/solvers/dgsem_structured/dg_3d.jl @@ -145,6 +145,8 @@ end # pull the contravariant vectors and compute the average Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, k, element) + # average mapping terms in first coordinate direction, + # used as normal vector in the flux computation Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector @@ -161,6 +163,8 @@ end # pull the contravariant vectors and compute the average Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, k, element) + # average mapping terms in second coordinate direction, + # used as normal vector in the flux computation Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector @@ -177,6 +181,8 @@ end # pull the contravariant vectors and compute the average Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors, i, j, kk, element) + # average mapping terms in third coordinate direction, + # used as normal vector in the flux computation Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector @@ -227,6 +233,8 @@ end # pull the contravariant vectors and compute the average Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, k, element) + # average mapping terms in first coordinate direction, + # used as normal vector in the flux computation Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector @@ -242,6 +250,8 @@ end # pull the contravariant vectors and compute the average Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, k, element) + # average mapping terms in second coordinate direction, + # used as normal vector in the flux computation Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector @@ -257,6 +267,8 @@ end # pull the contravariant vectors and compute the average Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors, i, j, kk, element) + # average mapping terms in third coordinate direction, + # used as normal vector in the flux computation Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector @@ -690,14 +702,6 @@ end return nothing end -# TODO: Taal dimension agnostic -function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, - mesh::StructuredMesh{3}, equations, surface_integral, - dg::DG) - @assert isperiodic(mesh) - return nothing -end - function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, mesh::StructuredMesh{3}, equations, surface_integral, dg::DG) diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index fb7f7f9dbf1..125773c1fd5 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -5,16 +5,6 @@ @muladd begin #! format: noindent -# 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 nothing -end - function volume_jacobian(element, mesh::TreeMesh, cache) return inv(cache.elements.inverse_jacobian[element])^ndims(mesh) end @@ -25,6 +15,14 @@ end return inverse_jacobian[element] end +# Dimension agnostic, i.e., valid for all 1D, 2D, and 3D `TreeMesh`es. +function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, + mesh::TreeMesh, equations, surface_integral, dg::DG) + @assert isempty(eachboundary(dg, cache)) + + return nothing +end + # Indicators used for shock-capturing and AMR include("indicators.jl") include("indicators_1d.jl") diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index d8345c70064..b0d3a2e3db0 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -33,10 +33,6 @@ end # The methods below are specialized on the volume integral type # and called from the basic `create_cache` method at the top. -function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations, - volume_integral::VolumeIntegralFluxDifferencing, dg::DG, uEltype) - NamedTuple() -end function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) @@ -122,20 +118,6 @@ function rhs!(du, u, t, return nothing end -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralWeakForm, - dg::DGSEM, cache) - @threaded for element in eachelement(dg, cache) - weak_form_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - dg, cache) - end - - return nothing -end - #= `weak_form_kernel!` is only implemented for conserved terms as non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, @@ -164,20 +146,6 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 return nothing end -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralFluxDifferencing, - dg::DGSEM, cache) - @threaded for element in eachelement(dg, cache) - flux_differencing_kernel!(du, u, element, mesh, have_nonconservative_terms, - equations, - volume_integral.volume_flux, dg, cache) - end - - return nothing -end - @inline function flux_differencing_kernel!(du, u, element, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, @@ -244,64 +212,6 @@ end end end -# TODO: Taal dimension agnostic -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralShockCapturingHG, - dg::DGSEM, cache) - @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral - - # Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α - alpha = @trixi_timeit timer() "blending factors" indicator(u, mesh, equations, dg, - cache) - - # For `Float64`, this gives 1.8189894035458565e-12 - # 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) - alpha_element = alpha[element] - # Clip blending factor for values close to zero (-> pure DG) - dg_only = isapprox(alpha_element, 0, atol = atol) - - if dg_only - flux_differencing_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - volume_flux_dg, dg, cache) - else - # Calculate DG volume integral contribution - flux_differencing_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - volume_flux_dg, dg, cache, 1 - alpha_element) - - # Calculate FV volume integral contribution - fv_kernel!(du, u, mesh, have_nonconservative_terms, equations, - volume_flux_fv, - dg, cache, element, alpha_element) - end - end - - return nothing -end - -# TODO: Taal dimension agnostic -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralPureLGLFiniteVolume, - dg::DGSEM, cache) - @unpack volume_flux_fv = volume_integral - - # Calculate LGL FV volume integral - @threaded for element in eachelement(dg, cache) - fv_kernel!(du, u, mesh, have_nonconservative_terms, equations, volume_flux_fv, - dg, cache, element, true) - end - - return nothing -end - @inline function fv_kernel!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, have_nonconservative_terms, equations, @@ -497,14 +407,6 @@ function prolong2boundaries!(cache, u, return nothing end -# TODO: Taal dimension agnostic -function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::TreeMesh{1}, equations, surface_integral, dg::DG) - @assert isempty(eachboundary(dg, cache)) - - return nothing -end - function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple, mesh::TreeMesh{1}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @@ -639,7 +541,7 @@ function apply_jacobian!(du, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, return nothing end -# TODO: Taal dimension agnostic +# Need dimension specific version to avoid error at dispatching function calc_sources!(du, u, t, source_terms::Nothing, equations::AbstractEquations{1}, dg::DG, cache) return nothing diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index ab3ab780ad6..0b188f9fbda 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -11,7 +11,7 @@ # This method is called when a SemidiscretizationHyperbolic is constructed. # It constructs the basic `cache` used throughout the simulation to compute # the RHS etc. -function create_cache(mesh::TreeMesh{2}, equations, +function create_cache(mesh::Union{TreeMesh{2}, TreeMesh{3}}, equations, dg::DG, RealT, uEltype) # Get cells for which an element needs to be created (i.e. all leaf cells) leaf_cell_ids = local_leaf_cells(mesh.tree) @@ -34,16 +34,6 @@ function create_cache(mesh::TreeMesh{2}, equations, return cache end -# The methods below are specialized on the volume integral type -# and called from the basic `create_cache` method at the top. -function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, - equations, volume_integral::VolumeIntegralFluxDifferencing, - dg::DG, uEltype) - NamedTuple() -end - function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) @@ -181,23 +171,6 @@ function rhs!(du, u, t, return nothing end -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}, - StructuredMeshView{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, - T8codeMesh{2}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralWeakForm, - dg::DGSEM, cache) - @threaded for element in eachelement(dg, cache) - weak_form_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - dg, cache) - end - - return nothing -end - #= `weak_form_kernel!` is only implemented for conserved terms as non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, @@ -233,24 +206,6 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 return nothing end -# flux differencing volume integral. For curved meshes averaging of the -# mapping terms, stored in `cache.elements.contravariant_vectors`, is peeled apart -# from the evaluation of the physical fluxes in each Cartesian direction -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}, - StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralFluxDifferencing, - dg::DGSEM, cache) - @threaded for element in eachelement(dg, cache) - flux_differencing_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - volume_integral.volume_flux, dg, cache) - end -end - @inline function flux_differencing_kernel!(du, u, element, mesh::TreeMesh{2}, have_nonconservative_terms::False, equations, @@ -333,67 +288,6 @@ end end end -# TODO: Taal dimension agnostic -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralShockCapturingHG, - dg::DGSEM, cache) - @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral - - # Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α - alpha = @trixi_timeit timer() "blending factors" indicator(u, mesh, equations, dg, - cache) - - # For `Float64`, this gives 1.8189894035458565e-12 - # 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) - alpha_element = alpha[element] - # Clip blending factor for values close to zero (-> pure DG) - dg_only = isapprox(alpha_element, 0, atol = atol) - - if dg_only - flux_differencing_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - volume_flux_dg, dg, cache) - else - # Calculate DG volume integral contribution - flux_differencing_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - volume_flux_dg, dg, cache, 1 - alpha_element) - - # Calculate FV volume integral contribution - fv_kernel!(du, u, mesh, have_nonconservative_terms, equations, - volume_flux_fv, dg, cache, element, alpha_element) - end - end - - return nothing -end - -# TODO: Taal dimension agnostic -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralPureLGLFiniteVolume, - dg::DGSEM, cache) - @unpack volume_flux_fv = volume_integral - - # Calculate LGL FV volume integral - @threaded for element in eachelement(dg, cache) - fv_kernel!(du, u, mesh, have_nonconservative_terms, equations, volume_flux_fv, - dg, cache, element, true) - end - - return nothing -end - @inline function fv_kernel!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, @@ -690,12 +584,6 @@ function prolong2boundaries!(cache, u, return nothing end -# TODO: Taal dimension agnostic -function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::TreeMesh{2}, equations, surface_integral, dg::DG) - @assert isempty(eachboundary(dg, cache)) -end - function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple, mesh::TreeMesh{2}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @@ -1193,7 +1081,7 @@ function apply_jacobian!(du, mesh::TreeMesh{2}, return nothing end -# TODO: Taal dimension agnostic +# Need dimension specific version to avoid error at dispatching function calc_sources!(du, u, t, source_terms::Nothing, equations::AbstractEquations{2}, dg::DG, cache) return nothing diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index cdf7b75b6d7..f87fcbdcd32 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -59,6 +59,7 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}}, flux_temp_threaded, fhat_temp_threaded) end +# Subcell limiting currently only implemented for certain mesh types function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}}, diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index d97f4e8c96c..875238c5438 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -8,41 +8,6 @@ # everything related to a DG semidiscretization in 3D, # currently limited to Lobatto-Legendre nodes -# This method is called when a SemidiscretizationHyperbolic is constructed. -# It constructs the basic `cache` used throughout the simulation to compute -# the RHS etc. -function create_cache(mesh::TreeMesh{3}, equations, - dg::DG, RealT, uEltype) - # Get cells for which an element needs to be created (i.e. all leaf cells) - leaf_cell_ids = local_leaf_cells(mesh.tree) - - elements = init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype) - - interfaces = init_interfaces(leaf_cell_ids, mesh, elements) - - boundaries = init_boundaries(leaf_cell_ids, mesh, elements) - - mortars = init_mortars(leaf_cell_ids, mesh, elements, dg.mortar) - - cache = (; elements, interfaces, boundaries, mortars) - - # Add specialized parts of the cache required to compute the volume integral etc. - cache = (; cache..., - create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) - cache = (; cache..., create_cache(mesh, equations, dg.mortar, uEltype)...) - - return cache -end - -# The methods below are specialized on the volume integral type -# and called from the basic `create_cache` method at the top. -function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, - T8codeMesh{3}}, - equations, volume_integral::VolumeIntegralFluxDifferencing, - dg::DG, uEltype) - NamedTuple() -end - function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, @@ -227,21 +192,6 @@ function rhs!(du, u, t, return nothing end -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, - T8codeMesh{3}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralWeakForm, - dg::DGSEM, cache) - @threaded for element in eachelement(dg, cache) - weak_form_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - dg, cache) - end - - return nothing -end - #= `weak_form_kernel!` is only implemented for conserved terms as non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, @@ -282,21 +232,6 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 return nothing end -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, - T8codeMesh{3}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralFluxDifferencing, - dg::DGSEM, cache) - @threaded for element in eachelement(dg, cache) - flux_differencing_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - volume_integral.volume_flux, dg, cache) - end - - return nothing -end - @inline function flux_differencing_kernel!(du, u, element, mesh::TreeMesh{3}, have_nonconservative_terms::False, equations, @@ -401,65 +336,6 @@ end return nothing end -# TODO: Taal dimension agnostic -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, - T8codeMesh{3}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralShockCapturingHG, - dg::DGSEM, cache) - @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral - - # Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α - alpha = @trixi_timeit timer() "blending factors" indicator(u, mesh, equations, dg, - cache) - - # For `Float64`, this gives 1.8189894035458565e-12 - # 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) - alpha_element = alpha[element] - # Clip blending factor for values close to zero (-> pure DG) - dg_only = isapprox(alpha_element, 0, atol = atol) - - if dg_only - flux_differencing_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - volume_flux_dg, dg, cache) - else - # Calculate DG volume integral contribution - flux_differencing_kernel!(du, u, element, mesh, - have_nonconservative_terms, equations, - volume_flux_dg, dg, cache, 1 - alpha_element) - - # Calculate FV volume integral contribution - fv_kernel!(du, u, mesh, have_nonconservative_terms, equations, - volume_flux_fv, dg, cache, element, alpha_element) - end - end - - return nothing -end - -# TODO: Taal dimension agnostic -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, - T8codeMesh{3}}, - have_nonconservative_terms, equations, - volume_integral::VolumeIntegralPureLGLFiniteVolume, - dg::DGSEM, cache) - @unpack volume_flux_fv = volume_integral - - # Calculate LGL FV volume integral - @threaded for element in eachelement(dg, cache) - fv_kernel!(du, u, mesh, have_nonconservative_terms, equations, volume_flux_fv, - dg, cache, element, true) - end - - return nothing -end - @inline function fv_kernel!(du, u, mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, @@ -799,13 +675,6 @@ function prolong2boundaries!(cache, u, return nothing end -# TODO: Taal dimension agnostic -function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::TreeMesh{3}, equations, surface_integral, dg::DG) - @assert isempty(eachboundary(dg, cache)) - return nothing -end - function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple, mesh::TreeMesh{3}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @@ -1517,7 +1386,7 @@ function apply_jacobian!(du, mesh::TreeMesh{3}, return nothing end -# TODO: Taal dimension agnostic +# Need dimension specific version to avoid error at dispatching function calc_sources!(du, u, t, source_terms::Nothing, equations::AbstractEquations{3}, dg::DG, cache) return nothing diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 8dedca2b112..6e5c92bd14e 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -20,7 +20,7 @@ end # this method is used when the indicator is constructed as for AMR function create_cache(typ::Type{IndicatorHennemannGassner}, mesh, - equations::AbstractEquations{1}, dg::DGSEM, cache) + equations::AbstractEquations, dg::DGSEM, cache) create_cache(typ, equations, dg.basis) end @@ -122,7 +122,7 @@ function create_cache(::Type{IndicatorLöhner}, equations::AbstractEquations{1}, end # this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{IndicatorLöhner}, mesh, equations::AbstractEquations{1}, +function create_cache(typ::Type{IndicatorLöhner}, mesh, equations::AbstractEquations, dg::DGSEM, cache) create_cache(typ, equations, dg.basis) end @@ -172,7 +172,7 @@ function create_cache(::Type{IndicatorMax}, equations::AbstractEquations{1}, end # this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{IndicatorMax}, mesh, equations::AbstractEquations{1}, +function create_cache(typ::Type{IndicatorMax}, mesh, equations::AbstractEquations, dg::DGSEM, cache) cache = create_cache(typ, equations, dg.basis) end diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index deaeff557c2..3bdc71b0ff6 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -22,12 +22,6 @@ function create_cache(::Type{IndicatorHennemannGassner}, return (; alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded) end -# this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{IndicatorHennemannGassner}, mesh, - equations::AbstractEquations{2}, dg::DGSEM, cache) - create_cache(typ, equations, dg.basis) -end - # Use this function barrier and unpack inside to avoid passing closures to Polyester.jl # with @batch (@threaded). # Otherwise, @threaded does not work here with Julia ARM on macOS. @@ -144,12 +138,6 @@ function create_cache(::Type{IndicatorLöhner}, equations::AbstractEquations{2}, return (; alpha, indicator_threaded) end -# this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{IndicatorLöhner}, mesh, equations::AbstractEquations{2}, - dg::DGSEM, cache) - create_cache(typ, equations, dg.basis) -end - function (löhner::IndicatorLöhner)(u::AbstractArray{<:Any, 4}, mesh, equations, dg::DGSEM, cache; kwargs...) @@ -203,12 +191,6 @@ function create_cache(::Type{IndicatorMax}, equations::AbstractEquations{2}, return (; alpha, indicator_threaded) end -# this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{IndicatorMax}, mesh, equations::AbstractEquations{2}, - dg::DGSEM, cache) - cache = create_cache(typ, equations, dg.basis) -end - function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 4}, mesh, equations, dg::DGSEM, cache; kwargs...) diff --git a/src/solvers/dgsem_tree/indicators_3d.jl b/src/solvers/dgsem_tree/indicators_3d.jl index 3fd228f9883..d384c3b34e4 100644 --- a/src/solvers/dgsem_tree/indicators_3d.jl +++ b/src/solvers/dgsem_tree/indicators_3d.jl @@ -25,12 +25,6 @@ function create_cache(::Type{IndicatorHennemannGassner}, modal_tmp2_threaded) end -# this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{IndicatorHennemannGassner}, mesh, - equations::AbstractEquations{3}, dg::DGSEM, cache) - create_cache(typ, equations, dg.basis) -end - # Use this function barrier and unpack inside to avoid passing closures to Polyester.jl # with @batch (@threaded). # Otherwise, @threaded does not work here with Julia ARM on macOS. @@ -164,12 +158,6 @@ function create_cache(::Type{IndicatorLöhner}, equations::AbstractEquations{3}, return (; alpha, indicator_threaded) end -# this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{IndicatorLöhner}, mesh, equations::AbstractEquations{3}, - dg::DGSEM, cache) - create_cache(typ, equations, dg.basis) -end - function (löhner::IndicatorLöhner)(u::AbstractArray{<:Any, 5}, mesh, equations, dg::DGSEM, cache; kwargs...) @@ -231,12 +219,6 @@ function create_cache(::Type{IndicatorMax}, equations::AbstractEquations{3}, return (; alpha, indicator_threaded) end -# this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{IndicatorMax}, mesh, equations::AbstractEquations{3}, - dg::DGSEM, cache) - cache = create_cache(typ, equations, dg.basis) -end - function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 5}, mesh, equations, dg::DGSEM, cache; kwargs...) diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index d69de9e685d..acea6f53615 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -308,7 +308,6 @@ function prolong2boundaries!(cache, u, return nothing end -# TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, mesh::Union{UnstructuredMesh2D, P4estMesh, P4estMeshView, T8codeMesh}, diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index a39f7cb1751..9eb7bdfa0b2 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -5,6 +5,17 @@ @muladd begin #! format: noindent +# Used by both `dg::DGSEM` and `dg::FDSBP` +function reset_du!(du, 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) + du[.., element] .= zero(eltype(du)) + end + + return nothing +end + # define types for parabolic solvers include("solvers_parabolic.jl") diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index d2cc65230c4..dbebd2a7007 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -84,7 +84,7 @@ end # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, amr_callback) - sol = solve(ode, KenCarp4(autodiff = AutoFiniteDiff()); + sol = solve(ode, ode_alg; abstol = time_abs_tol, reltol = time_int_tol, ode_default_options()..., callback = callbacks) l2_error, linf_error = analysis_callback(sol) diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl index 9a52a7ec2ea..1b1ecc7f098 100644 --- a/test/test_special_elixirs.jl +++ b/test/test_special_elixirs.jl @@ -188,7 +188,7 @@ end J_parabolic = jacobian_ad_forward_parabolic(semi) λ_parabolic = eigvals(J_parabolic) # Parabolic spectrum is real and negative - @test maximum(real, λ_parabolic) < 10^(-14) + @test maximum(real, λ_parabolic) < 2 * 10^(-14) @test maximum(imag, λ_parabolic) < 10^(-14) end @@ -273,7 +273,7 @@ end J_parabolic = jacobian_ad_forward_parabolic(semi) λ_parabolic = eigvals(J_parabolic) # Parabolic spectrum is real and negative - @test maximum(real, λ_parabolic) < 10^(-16) + @test maximum(real, λ_parabolic) < eps(Float64) @test maximum(imag, λ_parabolic) < 10^(-15) end