Merge branch 'main' into subcell-limiting

bennibolm · Aug 28, 2024 · 6c3b009 · 6c3b009
2 parents 0ac57af + 90eebb0
commit 6c3b009
Show file tree

Hide file tree

Showing 113 changed files with 2,885 additions and 1,292 deletions.
diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.22.9
+        uses: crate-ci/typos@v1.23.6
diff --git a/NEWS.md b/NEWS.md
@@ -4,6 +4,14 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes in the v0.8 lifecycle
+
+#### Changed
+
+- The AMR routines for `P4estMesh` and `T8codeMesh` were changed to work on the product
+  of the Jacobian and the conserved variables instead of the conserved variables only
+  to make AMR fully conservative ([#2028]). This may change AMR results slightly.
+
 ## Changes when updating to v0.8 from v0.7.x
 
 #### Added

diff --git a/Project.toml b/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.8.5-DEV"
+version = "0.8.9-DEV"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

diff --git a/docs/src/performance.md b/docs/src/performance.md
@@ -42,7 +42,8 @@ For example, the following steps were used to benchmark the changes introduced i
    ```julia
    julia> using BenchmarkTools, Revise; using Trixi
 
-   julia> trixi_include("examples/2d/elixir_euler_sedov_blast_wave.jl")
+   julia> # nowadays "examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl"
+          trixi_include("examples/2d/elixir_euler_sedov_blast_wave.jl")
 
    julia> du_test = copy(sol.u[end]); u_test = copy(sol.u[end]);
 
@@ -65,7 +66,8 @@ For example, the following steps were used to benchmark the changes introduced i
 
    shell> git checkout 222241ff54f8a4ca9876cc1fc25ae262416a4ea0
 
-   julia> trixi_include("examples/2d/elixir_euler_sedov_blast_wave.jl")
+   julia> # nowadays "examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl"
+          trixi_include("examples/2d/elixir_euler_sedov_blast_wave.jl")
 
    julia> @benchmark Trixi.rhs!(
              $(du_test),
@@ -85,6 +87,10 @@ For example, the following steps were used to benchmark the changes introduced i
     evals/sample:     1
    ```
    Run the `@benchmark ...` commands multiple times to see whether there are any significant fluctuations.
+   Note that the elixir name has changed since
+   [PR #256](https://github.com/trixi-framework/Trixi.jl/pull/256).
+   Nowadays, the relevant elixir is
+   [`examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl).
 
 Follow these steps for both commits you want to compare. The relevant benchmark results you should typically be looking at
 are the median and mean values of the runtime and the memory/allocs estimate. In this example, the differences

diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl
@@ -1,4 +1,4 @@
-using Downloads: download
+
 using OrdinaryDiffEq
 using Trixi
 

diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl
@@ -1,7 +1,6 @@
 
 using Trixi
 using OrdinaryDiffEq
-using Downloads: download
 
 ###############################################################################
 # semidiscretization of the compressible Euler equations
@@ -62,24 +61,24 @@ volume_integral = VolumeIntegralShockCapturingHG(shock_indicator;
                                                  volume_flux_dg = volume_flux,
                                                  volume_flux_fv = surface_flux)
 
-# DG Solver                                                 
+# DG Solver
 solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
                volume_integral = volume_integral)
 
 # Mesh generated from the following gmsh geometry input file:
 # https://gist.githubusercontent.com/DanielDoehring/5ade6d93629f0d8c23a598812dbee2a9/raw/d2bc904fe92146eae1a36156e7f5c535dc1a80f1/NACA6412.geo
 mesh_file = joinpath(@__DIR__, "mesh_NACA6412.inp")
 isfile(mesh_file) ||
-    download("https://gist.githubusercontent.com/DanielDoehring/e2a389f04f1e37b33819b9637e8ee4c3/raw/4bf7607a2ce4432fdb5cb87d5e264949b11bd5d7/mesh_NACA6412.inp",
-             mesh_file)
+    Trixi.download("https://gist.githubusercontent.com/DanielDoehring/e2a389f04f1e37b33819b9637e8ee4c3/raw/4bf7607a2ce4432fdb5cb87d5e264949b11bd5d7/mesh_NACA6412.inp",
+                   mesh_file)
 
 boundary_symbols = [:PhysicalLine1, :PhysicalLine2, :PhysicalLine3, :PhysicalLine4]
 
 mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols)
 
 boundary_conditions = Dict(:PhysicalLine1 => boundary_condition_supersonic_inflow, # Left boundary
                            :PhysicalLine2 => boundary_condition_supersonic_outflow, # Right boundary
-                           :PhysicalLine3 => boundary_condition_supersonic_outflow, # Top and bottom boundary 
+                           :PhysicalLine3 => boundary_condition_supersonic_outflow, # Top and bottom boundary
                            :PhysicalLine4 => boundary_condition_slip_wall) # Airfoil
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,

diff --git a/examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec_float32.jl b/examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec_float32.jl
@@ -0,0 +1,66 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4f0)
+
+initial_condition = initial_condition_weak_blast_wave
+
+surface_flux = flux_ranocha
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(Float32, polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 1.0f0,
+                                         alpha_min = 0.001f0,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral, RealT = Float32)
+
+###############################################################################
+
+coordinates_min = (-1.0f0, -1.0f0)
+coordinates_max = (+1.0f0, +1.0f0)
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 4, initial_refinement_level = 2,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = true, RealT = Float32)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0f0, 2.0f0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0f0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        stepsize_callback)
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0f0, # 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
diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl
@@ -1,4 +1,4 @@
-using Downloads: download
+
 using OrdinaryDiffEq
 using Trixi
 

diff --git a/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -0,0 +1,117 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+
+    r0 = 0.2
+    E = 1
+    p0_inner = 3
+    p0_outer = 1
+
+    # Calculate primitive variables
+    rho = 1.1
+    v1 = 0.0
+    v2 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+# Get the DG approximation space
+
+# Activate the shock capturing + flux differencing
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+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)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+###############################################################################
+
+# Affine type mapping to take the [-1,1]^2 domain
+# and warp it as described in https://arxiv.org/abs/2012.12040
+# Warping with the coefficient 0.2 is even more extreme.
+function mapping_twist(xi, eta)
+    y = eta + 0.125 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta)
+    x = xi + 0.125 * cos(0.5 * pi * xi) * cos(2.0 * pi * y)
+    return SVector(x, y)
+end
+
+# The mesh below can be made periodic
+# Create P4estMesh with 8 x 8 trees
+trees_per_dimension = (8, 8)
+mesh = P4estMesh(trees_per_dimension, polydeg = 4,
+                 mapping = mapping_twist,
+                 initial_refinement_level = 0,
+                 periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 400
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 0.2,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 0.05,
+                                      max_level = 2, max_threshold = 0.1)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);#, maxiters=4);
+summary_callback() # print the timer summary
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -1,4 +1,4 @@
-using Downloads: download
+
 using OrdinaryDiffEq
 using Trixi
 
@@ -15,7 +15,7 @@ using Trixi
 #   Structured and Unstructured Grid Methods (2016)
 #   [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623)
 # - Deep Ray, Praveen Chandrashekar (2017)
-#   An entropy stable finite volume scheme for the 
+#   An entropy stable finite volume scheme for the
 #   two dimensional Navier–Stokes equations on triangular grids
 #   [DOI:10.1016/j.amc.2017.07.020](https://doi.org/10.1016/j.amc.2017.07.020)
 

diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl
@@ -1,5 +1,4 @@
 
-using Downloads: download
 using OrdinaryDiffEq
 using Trixi
 
@@ -18,8 +17,8 @@ solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs)
 
 default_mesh_file = joinpath(@__DIR__, "mesh_cube_with_boundaries.inp")
 isfile(default_mesh_file) ||
-    download("https://gist.githubusercontent.com/DanielDoehring/710eab379fe3042dc08af6f2d1076e49/raw/38e9803bc0dab9b32a61d9542feac5343c3e6f4b/mesh_cube_with_boundaries.inp",
-             default_mesh_file)
+    Trixi.download("https://gist.githubusercontent.com/DanielDoehring/710eab379fe3042dc08af6f2d1076e49/raw/38e9803bc0dab9b32a61d9542feac5343c3e6f4b/mesh_cube_with_boundaries.inp",
+                   default_mesh_file)
 mesh_file = default_mesh_file
 
 boundary_symbols = [:PhysicalSurface1, :PhysicalSurface2]

diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries_float32.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries_float32.jl
@@ -0,0 +1,62 @@
+# Similar to p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl
+# but using Float32 instead of the default Float64
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4f0)
+
+initial_condition = initial_condition_constant
+
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs, RealT = Float32)
+
+###############################################################################
+# Get the uncurved mesh from a file (downloads the file if not available locally)
+
+default_mesh_file = joinpath(@__DIR__, "mesh_cube_with_boundaries.inp")
+isfile(default_mesh_file) ||
+    Trixi.download("https://gist.githubusercontent.com/DanielDoehring/710eab379fe3042dc08af6f2d1076e49/raw/38e9803bc0dab9b32a61d9542feac5343c3e6f4b/mesh_cube_with_boundaries.inp",
+                   default_mesh_file)
+mesh_file = default_mesh_file
+
+boundary_symbols = [:PhysicalSurface1, :PhysicalSurface2]
+
+mesh = P4estMesh{3}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols,
+                    RealT = Float32)
+
+boundary_conditions = Dict(:PhysicalSurface1 => BoundaryConditionDirichlet(initial_condition),
+                           :PhysicalSurface2 => BoundaryConditionDirichlet(initial_condition))
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0f0, 1.0f0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.5f0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0f0, # 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