Make volume integrals, indicators, etc dimension agnostic

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)

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

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

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,

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)

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,

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)

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")