Merge branch 'trixi-framework:main' into main

.gitignore

@@ -28,4 +28,6 @@ coverage_report/
 .DS_Store
 
 run
-run/*
+run/*
+
+LocalPreferences.toml

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.23-pre"
+version = "0.5.25-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
@@ -53,7 +53,7 @@ HDF5 = "0.14, 0.15, 0.16"
 IfElse = "0.1"
 LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.118"
-MPI = "0.20 - 0.20.8"
+MPI = "0.20"
 MuladdMacro = "0.2.2"
 Octavian = "0.3.5"
 OffsetArrays = "1.3"

docs/literate/src/files/scalar_linear_advection_1d.jl

@@ -44,10 +44,11 @@ dx = (coordinates_max - coordinates_min) / n_elements # length of one element
 # ```
 # Here, $J$ is the Jacobian determinant of the transformation.
 
-# Using this transformation, we can transform our equation for all elements $Q_l$.
+# Using this transformation, we can transform our equation for each element $Q_l$.
 # ```math
 # \frac{dx}{2} u_t^{Q_l} + u_\xi^{Q_l} = 0 \text{, for }t\in\mathbb{R}^+,\; \xi\in[-1, 1]
 # ```
+# Here, $u_t^{Q_l}$ and $u_\xi^{Q_l}$ denote the time and spatial derivatives of the solution on the element $Q_l$.
 
 
 # ### ii. Polynomial approach

docs/src/development.md

@@ -256,17 +256,13 @@ Breakpoints can be set by adding a line with the ```@infiltrate``` macro at the
 in the code. Use [Revise](@ref interactive-use-of-julia) if you want to set and delete breakpoints
 in your package without having to restart Julia.
 
-!!! note
-    When running Julia inside a package environment, the ```@infiltrate``` macro only works if `Infiltrator`
-    has been added to the dependencies. Another work around when using Revise is to first load the
-    package and then add breakpoints with `Main.@infiltrate` to the code. If this is not
-    desired, the functional form
-    ```julia
-    if isdefined(Main, :Infiltrator)
-      Main.Infiltrator.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__)
-    end
-    ```
-    can be used to set breakpoints when working with Trixi.jl or other packages.
+!!! note "Use `@autoinfiltrate` when debugging Trixi.jl"
+    When running Julia inside a package environment, e.g., inside the source
+    code of Trixi.jl itself, the `@infiltrate` macro only works if
+    `Infiltrator` has been added to the package dependencies. To avoid this,
+    you can use the (non-exported) `@autoinfiltrate` macro
+    in Trixi.jl, which only requires Infiltrator.jl to be available in the
+    current environment stack and will auto-load it for you.
 
 Triggering the breakpoint starts a REPL session where it is possible to interact with the current
 local scope. Possible commands are:

docs/src/parallelization.md

@@ -162,3 +162,13 @@ If you use error-based step size control (see also the section on [error-based a
 together with MPI you need to pass `internalnorm=ode_norm` and you should pass
 `unstable_check=ode_unstable_check` to OrdinaryDiffEq's [`solve`](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/),
 which are both included in [`ode_default_options`](@ref).
+
+### Using parallel input and output
+Trixi.jl allows parallel I/O using MPI by leveraging parallel HDF5.jl. To enable this, you first need
+to use a system-provided MPI library, see also [here](@ref parallel_system_MPI) and you need to tell
+[HDF5.jl](https://github.com/JuliaIO/HDF5.jl) to use this library.
+To do so, set the environment variable `JULIA_HDF5_PATH` to the local path
+that contains the `libhdf5.so` shared object file and build HDF5.jl by executing `using Pkg; Pkg.build("HDF5")`.
+For more information see also the [documentation of HDF5.jl](https://juliaio.github.io/HDF5.jl/stable/mpi/).
+
+If you do not perform these steps to use parallel HDF5 or if the HDF5 is not MPI-enabled, Trixi.jl will fall back on a less efficient I/O mechanism. In that case, all disk I/O is performed only on rank zero and data is distributed to/gathered from the other ranks using regular MPI communication.

docs/src/troubleshooting.md

@@ -175,14 +175,16 @@ to execute the following Julia code.
 
 ```julia
 using Preferences, UUIDs
-set_preferences!(UUID("1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"), "PrecompileNonStiff" => true)
-set_preferences!(UUID("1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"), "PrecompileStiff" => false)
-set_preferences!(UUID("1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"), "PrecompileAutoSwitch" => false)
-set_preferences!(UUID("1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"), "PrecompileLowStorage" => true)
-set_preferences!(UUID("1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"), "PrecompileDefaultSpecialize" => true)
-set_preferences!(UUID("1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"), "PrecompileAutoSpecialize" => false)
-set_preferences!(UUID("1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"), "PrecompileFunctionWrapperSpecialize" => false)
-set_preferences!(UUID("1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"), "PrecompileNoSpecialize" => false)
+let uuid = UUID("1dea7af3-3e70-54e6-95c3-0bf5283fa5ed")
+  set_preferences!(uuid, "PrecompileAutoSpecialize" => false)
+  set_preferences!(uuid, "PrecompileAutoSwitch" => false)
+  set_preferences!(uuid, "PrecompileDefaultSpecialize" => true)
+  set_preferences!(uuid, "PrecompileFunctionWrapperSpecialize" => false)
+  set_preferences!(uuid, "PrecompileLowStorage" => true)
+  set_preferences!(uuid, "PrecompileNoSpecialize" => false)
+  set_preferences!(uuid, "PrecompileNonStiff" => true)
+  set_preferences!(uuid, "PrecompileStiff" => false)
+end
 ```
 
 This disables precompilation of all implicit methods. This should usually not affect

examples/dgmulti_2d/elixir_euler_shockcapturing.jl

@@ -0,0 +1,51 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_weak_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux  = flux_ranocha
+
+polydeg = 3
+basis = DGMultiBasis(Quad(), polydeg, approximation_type=GaussSBP())
+
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max=0.5,
+                                         alpha_min=0.001,
+                                         alpha_smooth=true,
+                                         variable=density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+dg = DGMulti(basis,
+             surface_integral = SurfaceIntegralWeakForm(surface_flux),
+             volume_integral = volume_integral)
+
+cells_per_dimension = (8, 8)
+mesh = DGMultiMesh(dg, cells_per_dimension, periodicity=true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)
+
+tspan = (0.0, 0.15)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,
+            ode_default_options()..., callback=callbacks);
+
+summary_callback() # print the timer summary
+

examples/dgmulti_2d/elixir_euler_shockcapturing_curved.jl

@@ -0,0 +1,56 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_weak_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux  = flux_ranocha
+
+polydeg = 3
+basis = DGMultiBasis(Quad(), polydeg, approximation_type=GaussSBP())
+
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max=0.5,
+                                         alpha_min=0.001,
+                                         alpha_smooth=true,
+                                         variable=density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+dg = DGMulti(basis,
+             surface_integral = SurfaceIntegralWeakForm(surface_flux),
+             volume_integral = volume_integral)
+
+function mapping(xi, eta)
+  x = xi  + 0.1 * sin(pi * xi) * sin(pi * eta)
+  y = eta + 0.1 * sin(pi * xi) * sin(pi * eta)
+  return SVector(x, y)
+end
+cells_per_dimension = (16, 16)
+mesh = DGMultiMesh(dg, cells_per_dimension, mapping)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)
+
+tspan = (0.0, 0.15)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,
+            ode_default_options()..., callback=callbacks);
+
+summary_callback() # print the timer summary
+

examples/dgmulti_2d/elixir_mhd_reflective_wall.jl

@@ -0,0 +1,105 @@
+
+using OrdinaryDiffEq
+using Trixi
+using LinearAlgebra: norm, dot # for use in the MHD boundary condition
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+equations = IdealGlmMhdEquations2D(1.4)
+
+function initial_condition_perturbation(x, t, equations::IdealGlmMhdEquations2D)
+  # pressure perturbation in a vertically magnetized field on the domain [-1, 1]^2
+
+  r2 = (x[1] + 0.25)^2 + (x[2] + 0.25)^2
+
+  rho = 1.0
+  v1  = 0.0
+  v2  = 0.0
+  v3  = 0.0
+  p   = 1 + 0.5 * exp(-100 * r2)
+
+  # the pressure and magnetic field are chosen to be strongly
+  # magnetized, such that p / ||B||^2 ≈ 0.01.
+  B1  = 0.0
+  B2  = 40.0 / sqrt(4.0 * pi)
+  B3  = 0.0
+
+  psi = 0.0
+  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
+end
+initial_condition = initial_condition_perturbation
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell)
+volume_flux  = (flux_hindenlang_gassner, flux_nonconservative_powell)
+
+solver = DGMulti(polydeg=3, element_type = Quad(), approximation_type = GaussSBP(),
+                 surface_integral = SurfaceIntegralWeakForm(surface_flux),
+                 volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+x_neg(x, tol=50*eps()) = abs(x[1] + 1) < tol
+x_pos(x, tol=50*eps()) = abs(x[1] - 1) < tol
+y_neg(x, tol=50*eps()) = abs(x[2] + 1) < tol
+y_pos(x, tol=50*eps()) = abs(x[2] - 1) < tol
+is_on_boundary = Dict(:x_neg => x_neg, :x_pos => x_pos, :y_neg => y_neg, :y_pos => y_pos)
+
+cells_per_dimension = (16, 16)
+mesh = DGMultiMesh(solver, cells_per_dimension; periodicity=(false, false), is_on_boundary)
+
+# Create a "reflective-like" boundary condition by mirroring the velocity but leaving the magnetic field alone.
+# Note that this boundary condition is probably not entropy stable.
+function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, x, t,
+                                               surface_flux_function, equations::IdealGlmMhdEquations2D)
+
+  # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later
+  norm_ = norm(normal_direction)
+  normal = normal_direction / norm_
+
+  # compute the primitive variables
+  rho, v1, v2, v3, p, B1, B2, B3, psi = cons2prim(u_inner, equations)
+
+  v_normal = dot(normal, SVector(v1, v2))
+  u_mirror = prim2cons(SVector(rho, v1 - 2 * v_normal * normal[1],
+                                    v2 - 2 * v_normal * normal[2],
+                                    v3, p, B1, B2, B3, psi), equations)
+
+  return surface_flux_function(u_inner, u_mirror, normal, equations) * norm_
+end
+
+boundary_conditions = (; x_neg=boundary_condition_velocity_slip_wall,
+                         x_pos=boundary_condition_velocity_slip_wall,
+                         y_neg=boundary_condition_do_nothing,
+                         y_pos=BoundaryConditionDirichlet(initial_condition))
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions=boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.075)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(solver))
+alive_callback = AliveCallback(alive_interval=10)
+
+cfl = 0.5
+stepsize_callback = StepsizeCallback(cfl=cfl)
+glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1e-5, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+summary_callback() # print the timer summary

examples/tree_3d_fdsbp/elixir_advection_extended.jl

@@ -0,0 +1,56 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
+                            derivative_order=1, accuracy_order=4,
+                            xmin=0.0, xmax=1.0, N=10)
+solver = FDSBP(D_SBP,
+               surface_integral=SurfaceIntegralStrongForm(flux_lax_friedrichs),
+               volume_integral=VolumeIntegralStrongForm())
+
+coordinates_min = (-1.0, -1.0, -1.0)
+coordinates_max = ( 1.0,  1.0,  1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=1,
+                n_cells_max=30_000,
+                periodicity=true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(energy_total,))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-9, reltol=1.0e-9,
+            ode_default_options()..., callback=callbacks)
+summary_callback()