Merge branch 'main' into sc/converters_coupling

Project.toml

@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.5.46-pre"
+version = "0.5.47-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

docs/src/parallelization.md

@@ -70,7 +70,32 @@ the same for P4est.jl and T8code.jl. This could e.g. be `libp4est.so` that usual
 in `lib/` or `local/lib/` in the installation directory of `t8code`.
 In total, in your active Julia project you should have a LocalPreferences.toml file with sections
 `[MPIPreferences]`, `[T8code]` and `[P4est]` as well as an entry `MPIPreferences` in your
-Project.toml to use a custom MPI installation.
+Project.toml to use a custom MPI installation. A `LocalPreferences.toml` file 
+created as described above might look something like the following:
+```toml
+[HDF5]
+libhdf5 = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so"
+libhdf5_hl = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5_hl.so"
+
+[MPIPreferences]
+__clear__ = ["preloads_env_switch"]
+_format = "1.0"
+abi = "OpenMPI"
+binary = "system"
+cclibs = []
+libmpi = "/lib/x86_64-linux-gnu/libmpi.so"
+mpiexec = "mpiexec"
+preloads = []
+
+[P4est]
+libp4est = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libp4est.so"
+libsc = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libsc.so"
+
+[T8code]
+libp4est = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libp4est.so"
+libsc = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libsc.so"
+libt8 = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libt8.so"
+```
 
 
 ### [Usage](@id parallel_usage)

examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl

@@ -27,15 +27,15 @@ See Section IV c on the paper below for details.
     phi = pi / 6
     sin_phi, cos_phi = sincos(phi)
 
-    rho =  8
+    rho =  8.0
     v1  =  8.25 * cos_phi
     v2  = -8.25 * sin_phi
     p   =  116.5
   else
     rho = 1.4
-    v1  = 0
-    v2  = 0
-    p   = 1
+    v1  = 0.0
+    v2  = 0.0
+    p   = 1.0
   end
 
   prim = SVector(rho, v1, v2, p)

src/callbacks_step/save_solution.jl

@@ -211,11 +211,11 @@ end
         get_element_variables!(element_variables, u_ode, semi)
         callbacks = integrator.opts.callback
         if callbacks isa CallbackSet
-            for cb in callbacks.continuous_callbacks
+            foreach(callbacks.continuous_callbacks) do cb
                 get_element_variables!(element_variables, u_ode, semi, cb;
                                        t = integrator.t, iter = iter)
             end
-            for cb in callbacks.discrete_callbacks
+            foreach(callbacks.discrete_callbacks) do cb
                 get_element_variables!(element_variables, u_ode, semi, cb;
                                        t = integrator.t, iter = iter)
             end

src/callbacks_step/summary.jl

@@ -168,16 +168,17 @@ function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator;
 
     callbacks = integrator.opts.callback
     if callbacks isa CallbackSet
-        for cb in callbacks.continuous_callbacks
+        foreach(callbacks.continuous_callbacks) do cb
             show(io_context, MIME"text/plain"(), cb)
             println(io, "\n")
         end
-        for cb in callbacks.discrete_callbacks
+        foreach(callbacks.discrete_callbacks) do cb
             # Do not show ourselves
-            cb.affect! === summary_callback && continue
+            cb.affect! === summary_callback && return nothing
 
             show(io_context, MIME"text/plain"(), cb)
             println(io, "\n")
+            return nothing
         end
     else
         show(io_context, MIME"text/plain"(), callbacks)

src/solvers/dg.jl

@@ -40,6 +40,11 @@ standard textbooks.
   Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and
   Applications
   [doi: 10.1007/978-0-387-72067-8](https://doi.org/10.1007/978-0-387-72067-8)
+
+`VolumeIntegralWeakForm()` is only implemented for conserved terms as 
+non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, 
+see [`VolumeIntegralFluxDifferencing`](@ref).
+This treatment is required to achieve, e.g., entropy-stability or well-balancedness.
 """
 struct VolumeIntegralWeakForm <: AbstractVolumeIntegral end
 

src/solvers/dgsem_structured/dg_2d.jl

@@ -49,6 +49,13 @@ function rhs!(du, u, t,
     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, 
+see `flux_differencing_kernel!`.
+This treatment is required to achieve, e.g., entropy-stability or well-balancedness.
+See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064
+=#
 @inline function weak_form_kernel!(du, u,
                                    element,
                                    mesh::Union{StructuredMesh{2}, UnstructuredMesh2D,

src/solvers/dgsem_structured/dg_3d.jl

@@ -49,6 +49,13 @@ function rhs!(du, u, t,
     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, 
+see `flux_differencing_kernel!`.
+This treatment is required to achieve, e.g., entropy-stability or well-balancedness.
+See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064
+=#
 @inline function weak_form_kernel!(du, u,
                                    element,
                                    mesh::Union{StructuredMesh{3}, P4estMesh{3}},

src/solvers/dgsem_tree/dg_1d.jl

@@ -144,6 +144,13 @@ function calc_volume_integral!(du, u,
     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, 
+see `flux_differencing_kernel!`.
+This treatment is required to achieve, e.g., entropy-stability or well-balancedness.
+See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064
+=#
 @inline function weak_form_kernel!(du, u,
                                    element, mesh::Union{TreeMesh{1}, StructuredMesh{1}},
                                    nonconservative_terms::False, equations,

src/solvers/dgsem_tree/dg_2d.jl

@@ -194,6 +194,13 @@ function calc_volume_integral!(du, u,
     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, 
+see `flux_differencing_kernel!`.
+This treatment is required to achieve, e.g., entropy-stability or well-balancedness.
+See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064
+=#
 @inline function weak_form_kernel!(du, u,
                                    element, mesh::TreeMesh{2},
                                    nonconservative_terms::False, equations,

src/solvers/dgsem_tree/dg_3d.jl

@@ -223,6 +223,13 @@ function calc_volume_integral!(du, u,
     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, 
+see `flux_differencing_kernel!`.
+This treatment is required to achieve, e.g., entropy-stability or well-balancedness.
+See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064
+=#
 @inline function weak_form_kernel!(du, u,
                                    element, mesh::TreeMesh{3},
                                    nonconservative_terms::False, equations,

src/time_integration/methods_2N.jl

@@ -119,10 +119,10 @@ function solve(ode::ODEProblem, alg::T;
 
     # initialize callbacks
     if callback isa CallbackSet
-        for cb in callback.continuous_callbacks
+        foreach(callback.continuous_callbacks) do cb
             error("unsupported")
         end
-        for cb in callback.discrete_callbacks
+        foreach(callback.discrete_callbacks) do cb
             cb.initialize(cb, integrator.u, integrator.t, integrator)
         end
     elseif !isnothing(callback)
@@ -172,10 +172,11 @@ function solve!(integrator::SimpleIntegrator2N)
 
         # handle callbacks
         if callbacks isa CallbackSet
-            for cb in callbacks.discrete_callbacks
+            foreach(callbacks.discrete_callbacks) do cb
                 if cb.condition(integrator.u, integrator.t, integrator)
                     cb.affect!(integrator)
                 end
+                return nothing
             end
         end
 

src/time_integration/methods_3Sstar.jl

@@ -189,10 +189,10 @@ function solve(ode::ODEProblem, alg::T;
 
     # initialize callbacks
     if callback isa CallbackSet
-        for cb in callback.continuous_callbacks
+        foreach(callback.continuous_callbacks) do cb
             error("unsupported")
         end
-        for cb in callback.discrete_callbacks
+        foreach(callback.discrete_callbacks) do cb
             cb.initialize(cb, integrator.u, integrator.t, integrator)
         end
     elseif !isnothing(callback)
@@ -248,10 +248,11 @@ function solve!(integrator::SimpleIntegrator3Sstar)
 
         # handle callbacks
         if callbacks isa CallbackSet
-            for cb in callbacks.discrete_callbacks
+            foreach(callbacks.discrete_callbacks) do cb
                 if cb.condition(integrator.u, integrator.t, integrator)
                     cb.affect!(integrator)
                 end
+                return nothing
             end
         end
 

src/time_integration/methods_SSP.jl

@@ -124,10 +124,10 @@ function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP;
 
     # initialize callbacks
     if callback isa CallbackSet
-        for cb in callback.continuous_callbacks
+        foreach(callback.continuous_callbacks) do cb
             error("unsupported")
         end
-        for cb in callback.discrete_callbacks
+        foreach(callback.discrete_callbacks) do cb
             cb.initialize(cb, integrator.u, integrator.t, integrator)
         end
     elseif !isnothing(callback)
@@ -148,7 +148,7 @@ function solve!(integrator::SimpleIntegratorSSP)
     callbacks = integrator.opts.callback
 
     integrator.finalstep = false
-    while !integrator.finalstep
+    @trixi_timeit timer() "main loop" while !integrator.finalstep
         if isnan(integrator.dt)
             error("time step size `dt` is NaN")
         end
@@ -184,7 +184,7 @@ function solve!(integrator::SimpleIntegratorSSP)
 
         # handle callbacks
         if callbacks isa CallbackSet
-            for cb in callbacks.discrete_callbacks
+            foreach(callbacks.discrete_callbacks) do cb
                 if cb.condition(integrator.u, integrator.t, integrator)
                     cb.affect!(integrator)
                 end

test/test_dgmulti_1d.jl

@@ -19,6 +19,14 @@ isdir(outdir) && rm(outdir, recursive=true)
       l2 = [2.9953644500009865e-5],
       linf = [4.467840577382365e-5]
     )
+    # Ensure that we do not have excessive memory allocations 
+    # (e.g., from type instabilities) 
+    let 
+      t = sol.t[end] 
+      u_ode = sol.u[end] 
+      du_ode = similar(u_ode) 
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 
+    end
   end
 
   @trixi_testset "elixir_euler_flux_diff.jl " begin
@@ -28,11 +36,27 @@ isdir(outdir) && rm(outdir, recursive=true)
       l2 = [7.853842541289665e-7, 9.609905503440606e-7, 2.832322219966481e-6] ./ sqrt(2.0),
       linf = [1.5003758788711963e-6, 1.802998748523521e-6, 4.83599270806323e-6]
     )
+    # Ensure that we do not have excessive memory allocations 
+    # (e.g., from type instabilities) 
+    let 
+      t = sol.t[end] 
+      u_ode = sol.u[end] 
+      du_ode = similar(u_ode) 
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 
+    end
   end
 
   @trixi_testset "elixir_euler_flux_diff.jl (convergence)" begin
     mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_euler_flux_diff.jl"), 3)
     @test isapprox(mean_convergence[:l2], [4.1558759698638434, 3.977911306037128, 4.041421206468769], rtol=0.05)
+    # Ensure that we do not have excessive memory allocations 
+    # (e.g., from type instabilities) 
+    let 
+      t = sol.t[end] 
+      u_ode = sol.u[end] 
+      du_ode = similar(u_ode) 
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 
+    end
   end
 
   @trixi_testset "elixir_euler_flux_diff.jl (SBP) " begin
@@ -42,6 +66,14 @@ isdir(outdir) && rm(outdir, recursive=true)
       l2 = [6.437827414849647e-6, 2.1840558851820947e-6, 1.3245669629438228e-5],
       linf = [2.0715843751295537e-5, 8.519520630301258e-6, 4.2642194098885255e-5]
     )
+    # Ensure that we do not have excessive memory allocations 
+    # (e.g., from type instabilities) 
+    let 
+      t = sol.t[end] 
+      u_ode = sol.u[end] 
+      du_ode = similar(u_ode) 
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 
+    end
   end
 
   @trixi_testset "elixir_euler_flux_diff.jl (FD SBP)" begin
@@ -56,6 +88,14 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
     show(stdout, semi.solver.basis)
     show(stdout, MIME"text/plain"(), semi.solver.basis)
+    # Ensure that we do not have excessive memory allocations 
+    # (e.g., from type instabilities) 
+    let 
+      t = sol.t[end] 
+      u_ode = sol.u[end] 
+      du_ode = similar(u_ode) 
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 
+    end
   end
 
   @trixi_testset "elixir_euler_fdsbp_periodic.jl" begin
@@ -65,6 +105,14 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
     show(stdout, semi.solver.basis)
     show(stdout, MIME"text/plain"(), semi.solver.basis)
+    # Ensure that we do not have excessive memory allocations 
+    # (e.g., from type instabilities) 
+    let 
+      t = sol.t[end] 
+      u_ode = sol.u[end] 
+      du_ode = similar(u_ode) 
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 
+    end
   end
 
   @trixi_testset "DGMulti with periodic SBP unit test" begin