Merge branch 'main' into dev

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.4.39-pre"
+version = "0.4.44-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
@@ -50,7 +50,7 @@ GeometryBasics = "0.3, 0.4"
 HDF5 = "0.14, 0.15, 0.16"
 IfElse = "0.1"
 LinearMaps = "2.7, 3.0"
-LoopVectorization = "0.12.117"
+LoopVectorization = "0.12.118"
 MPI = "0.19"
 MuladdMacro = "0.2.2"
 Octavian = "0.3.5"
@@ -63,8 +63,8 @@ Requires = "1.1"
 SciMLBase = "1.21"
 Setfield = "0.8, 1"
 StartUpDG = "0.13.1"
-Static = "0.3, 0.4, 0.5, 0.6"
-StaticArrays = "1.0 - 1.4.4"
+Static = "0.3, 0.4, 0.5, 0.6, 0.7"
+StaticArrays = "1"
 StrideArrays = "0.1.18"
 StructArrays = "0.6"
 SummationByPartsOperators = "0.5.10"

README.md

@@ -17,6 +17,18 @@
   <img width="300px" src="https://trixi-framework.github.io/assets/logo.png">
 </p>
 
+***
+**Trixi.jl at JuliaCon 2022**<br/>
+At this year's JuliaCon, we will be present with two contributions that involve Trixi.jl:
+
+* [Running Julia code in parallel with MPI: Lessons learned](https://live.juliacon.org/talk/LUWYRJ), 
+  26th July 2022
+* [From Mesh Generation to Adaptive Simulation: A Journey in Julia](https://live.juliacon.org/talk/YSLKZJ), 
+  27th July 2022
+
+We are looking forward to seeing you there ♥️
+***
+
 **Trixi.jl** is a numerical simulation framework for hyperbolic conservation
 laws written in [Julia](https://julialang.org). A key objective for the
 framework is to be useful to both scientists and students. Therefore, next to
@@ -30,7 +42,7 @@ installation and postprocessing procedures. Its features include:
   * Structured and unstructured meshes
   * Hierarchical quadtree/octree grid with adaptive mesh refinement
   * Forests of quadtrees/octrees with [p4est](https://github.com/cburstedde/p4est) via [P4est.jl](https://github.com/trixi-framework/P4est.jl)
-* High-order accuracy in space in time
+* High-order accuracy in space and time
 * Discontinuous Galerkin methods
   * Kinetic energy-preserving and entropy-stable methods based on flux differencing
   * Entropy-stable shock capturing
@@ -135,7 +147,7 @@ julia> plot(sol) # No trailing semicolon, otherwise no plot is shown
 ```
 This will open a new window with a 2D visualization of the final solution:
 <p align="center">
-  <img width="300px" src="https://user-images.githubusercontent.com/72009492/130952732-633159ff-c167-4d36-ba36-f2a2eac0a8d6.PNG">
+  <img width="300px" src="https://user-images.githubusercontent.com/26361975/177492363-74cee347-7abe-4522-8b2d-0dfadc317f7e.png">
 </p>
 
 The method `trixi_include(...)` expects a single string argument with the path to a

docs/Project.toml

@@ -15,7 +15,7 @@ Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"
 CairoMakie = "0.6, 0.7, 0.8"
 Documenter = "0.27"
 ForwardDiff = "0.10"
-HOHQMesh = "0.1"
+HOHQMesh = "0.1, 0.2"
 LaTeXStrings = "1.2"
 Literate = "2.9"
 Measurements = "2.5"

docs/literate/src/files/DGSEM_FluxDiff.jl

@@ -82,7 +82,13 @@
 # When using the diagonal SBP property it is possible to rewrite the application of the derivative
 # operator $D$ in the calculation of the volume integral into a subcell based finite volume type
 # differencing formulation ([Fisher, Carpenter (2013)](https://doi.org/10.1016/j.jcp.2013.06.014)).
-# We replace $D \underline{f}$ in the strong form by $2D \underline{f}_{vol}(u^-, u^+)$ with
+# Generalizing
+# ```math
+# (D \underline{f})_i = \sum_j D_{i,j} \underline{f}_j
+# = 2\sum_j \frac{1}{2} D_{i,j} (\underline{f}_j + \underline{f}_i)
+# \eqqcolon 2\sum_j  D_{i,j} f_\text{central}(u_i, u_j),
+# ```
+# we replace $D \underline{f}$ in the strong form by $2D \underline{f}_{vol}(u^-, u^+)$ with
 # the consistent two-point volume flux $f_{vol}$ and receive the DGSEM formulation with flux differencing
 # (split form DGSEM) ([Gassner, Winters, Kopriva (2016)](https://doi.org/10.1016/j.jcp.2016.09.013)).
 

docs/literate/src/files/time_stepping.jl

@@ -46,6 +46,8 @@
 # ```math
 # \Delta t_n = \text{CFL} * \min_i \frac{\Delta x_i}{\lambda_{\max}(u_i^n)}
 # ```
+# We compute $\Delta x_i$ by scaling the element size by a factor of $1/(N+1)$, cf. 
+# [Gassner and Kopriva (2011)](https://doi.org/10.1137/100807211), Section 5.
 
 # Trixi provides such a CFL-based step size control. It is implemented as the callback
 # [`StepsizeCallback`](@ref).

docs/src/index.md

@@ -24,7 +24,7 @@ installation and postprocessing procedures. Its features include:
   * Structured and unstructured meshes
   * Hierarchical quadtree/octree grid with adaptive mesh refinement
   * Forests of quadtrees/octrees with [p4est](https://github.com/cburstedde/p4est) via [P4est.jl](https://github.com/trixi-framework/P4est.jl)
-* High-order accuracy in space in time
+* High-order accuracy in space and time
 * Discontinuous Galerkin methods
   * Kinetic energy-preserving and entropy-stable methods based on flux differencing
   * Entropy-stable shock capturing
@@ -174,7 +174,7 @@ julia> plot(sol) # No trailing semicolon, otherwise no plot is shown
 ```
 This will open a new window with a 2D visualization of the final solution:
 
-![image](https://user-images.githubusercontent.com/72009492/130952732-633159ff-c167-4d36-ba36-f2a2eac0a8d6.PNG)
+![image](https://user-images.githubusercontent.com/26361975/177492363-74cee347-7abe-4522-8b2d-0dfadc317f7e.png)
 
 The method `trixi_include(...)` expects a single string argument with the path to a
 Trixi elixir, i.e., a text file containing Julia code necessary to set up and run a

docs/src/visualization.md

@@ -221,8 +221,8 @@ This creates the following plot:
 ![ScalarPlotData2D_example](https://user-images.githubusercontent.com/1156048/133856590-a9f0be02-8200-483b-af96-eab4a69bf2c7.png)
 
 ### Plotting a 3D solution as a 2D plot
-It is possible to plot 2D slices from 3D simulation data with the same commands
-as above:
+It is possible to plot 2D slices from 3D simulation data using the [`TreeMesh`](@ref)
+with the same commands as above:
 ```julia
 julia> plot(sol) # `sol` is from a 3D simulation
 ```

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl

@@ -0,0 +1,125 @@
+using OrdinaryDiffEq
+using Trixi
+
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+"""
+    initial_condition_kelvin_helmholtz_instability_fjordholm_etal(x, t, equations::CompressibleEulerEquations2D)
+
+A version of the classical Kelvin-Helmholtz instability based on
+- Ulrik S. Fjordholm, Roger Käppeli, Siddhartha Mishra, Eitan Tadmor (2014)
+  Construction of approximate entropy measure valued
+  solutions for hyperbolic systems of conservation laws
+  [arXiv: 1402.0909](https://arxiv.org/abs/1402.0909)
+"""
+function initial_condition_kelvin_helmholtz_instability_fjordholm_etal(x, t, equations::CompressibleEulerEquations2D)
+  # typical resolution 128^2, 256^2
+  # domain size is [0,+1]^2
+  # interface is sharp, but randomly perturbed
+  # The random numbers used in the initial conditions have been generated as follows:
+  #
+  # using StableRNGs
+  #
+  # rng = StableRNG(100)
+  #
+  # a1 = rand(rng, m)
+  # a2 = rand(rng, m)
+  # a1 .= a1 / sum(a1)
+  # a2 .= a2 / sum(a2)
+  # b1 = (rand(rng, m) .- 0.5) .* pi
+  # b2 = (rand(rng, m) .- 0.5) .* pi
+
+  m = 10
+  a1 = [0.04457096674422902, 0.03891512410182607, 0.0030191053979293433, 0.0993913172320319,
+        0.1622302137588842, 0.1831383653456182, 0.11758003014101702, 0.07964318348142958,
+        0.0863245324711805, 0.18518716132585408]
+  a2 = [0.061688440856337096, 0.23000237877135882, 0.04453793881833177, 0.19251530387370916,
+        0.11107917357941084, 0.05898041974649702, 0.09949312336096268, 0.07022276346006465,
+        0.10670366489014596, 0.02477679264318211]
+  b1 = [0.06582340543754152, 0.9857886297001535, 0.8450452205037154, -1.279648120993805,
+        0.45454198915209526, -0.13359370986823993, 0.07062615913363897, -1.0097986278512623,
+        1.0810669017430343, -0.14207309803877177]
+  b2 = [-1.1376882185131414, -1.4798197129947765, 0.6139290513283818, -0.3319087388365522,
+        0.14633328999192285, -0.06373231463100072, -0.6270101051216724, 0.13941252226261905,
+        -1.0337526453303645, 1.0441408867083155]
+  Y1 = 0.0
+  Y2 = 0.0
+  for n = 1:m
+    Y1 += a1[n] * cos(b1[n] + 2 * n * pi * x[1])
+    Y2 += a2[n] * cos(b2[n] + 2 * n * pi * x[1])
+  end
+
+  J1 = 0.25
+  J2 = 0.75
+  epsilon = 0.01
+  I1 = J1 + epsilon * Y1
+  I2 = J2 + epsilon * Y2
+
+  if (x[2] > I1) && (x[2] < I2)
+    rho = 2
+    v1 = -0.5
+  else
+    rho = 1
+    v1 = 0.5
+  end
+  v2 = 0
+  p = 2.5
+
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_kelvin_helmholtz_instability_fjordholm_etal
+
+surface_flux = flux_hllc
+volume_flux  = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max=0.001,
+                                         alpha_min=0.0001,
+                                         alpha_smooth=true,
+                                         variable=density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=6,
+                n_cells_max=100_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 400
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=400,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution)
+
+###############################################################################
+# run the simulation
+sol = solve(ode, SSPRK43(),
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

src/Trixi.jl

@@ -250,6 +250,7 @@ function __init__()
 
   # FIXME upstream. This is a hacky workaround for
   #       https://github.com/trixi-framework/Trixi.jl/issues/628
+  #       https://github.com/trixi-framework/Trixi.jl/issues/1185
   # The related upstream issues appear to be
   #       https://github.com/JuliaLang/julia/issues/35800
   #       https://github.com/JuliaLang/julia/issues/32552
@@ -258,9 +259,12 @@ function __init__()
   let
     for T in (Float32, Float64)
       u_mortars_2d = zeros(T, 2, 2, 2, 2, 2)
-      view(u_mortars_2d, 1, :, 1, :, 1)
+      u_view_2d = view(u_mortars_2d, 1, :, 1, :, 1)
+      LoopVectorization.axes(u_view_2d)
+
       u_mortars_3d = zeros(T, 2, 2, 2, 2, 2, 2)
-      view(u_mortars_3d, 1, :, 1, :, :, 1)
+      u_view_3d = view(u_mortars_3d, 1, :, 1, :, :, 1)
+      LoopVectorization.axes(u_view_3d)
     end
   end
 end

src/equations/compressible_euler_2d.jl

@@ -335,13 +335,14 @@ Should be used together with [`TreeMesh`](@ref).
                                               equations::CompressibleEulerEquations2D)
   # get the appropriate normal vector from the orientation
   if orientation == 1
-    normal = SVector(1, 0)
+    normal_direction = SVector(1, 0)
   else # orientation == 2
-    normal = SVector(0, 1)
+    normal_direction = SVector(0, 1)
   end
 
   # compute and return the flux using `boundary_condition_slip_wall` routine above
-  return boundary_condition_slip_wall(u_inner, normal, x, t, surface_flux_function, equations)
+  return boundary_condition_slip_wall(u_inner, normal_direction, direction,
+                                      x, t, surface_flux_function, equations)
 end
 
 """

test/Project.toml

@@ -1,7 +1,6 @@
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
-Cassette = "7057c7e9-c182-5462-911a-8362d720325c"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -15,7 +14,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 BSON = "0.3.3"
 CairoMakie = "0.6, 0.7, 0.8"
-Cassette = "0.3.5"
 Flux = "0.13"
 ForwardDiff = "0.10"
 MPI = "0.19"

test/test_tree_2d_euler.jl

@@ -164,6 +164,13 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2
       coverage_override = (maxiters=10^5,))
   end
 
+  @trixi_testset "elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl"),
+      l2   = [0.1057230211245312, 0.10621112311257341, 0.07260957505339989, 0.11178239111065721],
+      linf = [2.998719417992662, 2.1400285015556166, 1.1569648700415078, 1.8922492268110913],
+      tspan = (0.0, 0.1))
+  end
+
   @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
       l2   = [0.055691508271624536, 0.032986009333751655, 0.05224390923711999, 0.08009536362771563],