Merge branch 'main' into SecondOrderFiniteVolume1D

trixi-framework · Sep 24, 2024 · 0bcb5fc · 0bcb5fc
2 parents a64bb60 + 827a738
commit 0bcb5fc
Show file tree

Hide file tree

Showing 198 changed files with 4,833 additions and 2,357 deletions.
diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml
@@ -29,7 +29,7 @@ jobs:
         # TODO: Change the call below to
         #       format(".")
         run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))'
           julia  -e 'using JuliaFormatter; format(["benchmark", "examples", "ext", "src", "test", "utils"])'
       - name: Format check
         run: |

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.23.6
+        uses: crate-ci/typos@v1.24.3
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,9 +1,10 @@
 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.6-DEV"
+version = "0.8.11-DEV"
 
 [deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -29,6 +30,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -50,15 +52,16 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [weakdeps]
-Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 
 [extensions]
-TrixiMakieExt = "Makie"
 TrixiConvexECOSExt = ["Convex", "ECOS"]
+TrixiMakieExt = "Makie"
 
 [compat]
+Accessors = "0.1.12"
 CodeTracking = "1.0.5"
 ConstructionBase = "1.3"
 Convex = "0.16"
@@ -87,12 +90,13 @@ PrecompileTools = "1.1"
 Preferences = "1.3"
 Printf = "1"
 RecipesBase = "1.1"
+RecursiveArrayTools = "2.38.10"
 Reexport = "1.0"
 Requires = "1.1"
 SciMLBase = "1.90, 2"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
-StartUpDG = "0.17.7"
+StartUpDG = "0.17.7, 1.1.5"
 Static = "0.8.7"
 StaticArrayInterface = "1.4"
 StaticArrays = "1.5"
@@ -109,6 +113,6 @@ UUIDs = "1.6"
 julia = "1.8"
 
 [extras]
-Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
diff --git a/README.md b/README.md
@@ -19,28 +19,6 @@
   <img width="300px" src="https://trixi-framework.github.io/assets/logo.png">
 </p>
 
-***
-**Trixi.jl at JuliaCon 2024**<br/>
-At this year's JuliaCon in Eindhoven, Netherlands, we will be present with several contributions
-from the Trixi Framework ecosystem:
-
-* [**Julia for Particle-Based Multiphysics with TrixiParticles.jl**](https://pretalx.com/juliacon2024/talk/TPFF8L/),<br/>
-  [*Erik Faulhaber*](https://github.com/efaulhaber/), [*Niklas Neher*](https://github.com/lasnikas/),
-  10th July 2024, 11:00–11:30, Function (4.1)
-* [**Towards Aerodynamic Simulations in Julia with Trixi.jl**](https://pretalx.com/juliacon2024/talk/XH8KBG/),<br/>
-  [*Daniel Doehring*](https://github.com/danieldoehring/),
-  10th July 2024, 15:30–16:00, While Loop (4.2)
-* [**libtrixi: serving legacy codes in earth system modeling with fresh Julia CFD**](https://pretalx.com/juliacon2024/talk/SXC7LA/),<br/>
-  [*Benedict Geihe*](https://github.com/benegee/),
-  12th July 2024, 14:00–17:00, Function (4.1)
-
-The last talk is part of the 
-[Julia for High-Performance Computing](https://juliacon.org/2024/minisymposia/hpc/)
-minisymposium, which this year is hosted by our own [*Hendrik Ranocha*](https://github.com/ranocha/).
-
-We are looking forward to seeing you there ♥️
-***
-
 **Trixi.jl** is a numerical simulation framework for 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

diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl
@@ -51,7 +51,11 @@
 # The Newton-bisection algorithm is an iterative method and requires some parameters.
 # It uses a fixed maximum number of iteration steps (`max_iterations_newton = 10`) and
 # relative/absolute tolerances (`newton_tolerances = (1.0e-12, 1.0e-14)`). The given values are
-# sufficient in most cases and therefore used as default. Additionally, there is the parameter
+# sufficient in most cases and therefore used as default. If the implemented bounds checking
+# functionality indicates problems with the limiting (see [below](@ref subcell_bounds_check))
+# the Newton method with the chosen parameters might not manage to converge. If so, adapting
+# the mentioned parameters helps fix that.
+# Additionally, there is the parameter
 # `gamma_constant_newton`, which can be used to scale the antidiffusive flux for the computation
 # of the blending coefficients of nonlinear variables. The default value is `2 * ndims(equations)`,
 # as it was shown by [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) [Section 4.2.2.]
@@ -244,7 +248,7 @@ plot(sol)
 # ![blast_wave_paraview_reinterpolate=false](https://github.com/trixi-framework/Trixi.jl/assets/74359358/39274f18-0064-469c-b4da-bac4b843e116)
 
 
-# ## Bounds checking
+# ## [Bounds checking](@id subcell_bounds_check)
 # Subcell limiting is based on the fulfillment of target bounds - either global or local.
 # Although the implementation works and has been thoroughly tested, there are some cases where
 # these bounds are not met.

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/docs/src/styleguide.md b/docs/src/styleguide.md
@@ -51,7 +51,7 @@ PRs that verify that running JuliaFormatter.jl again will not change the source
 To format your contributions before created a PR (or, at least, before requesting a review
 of your PR), you need to install JuliaFormatter.jl first by running
 ```shell
-julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))'
 ```
 You can then recursively format the core Julia files in the Trixi.jl repo by executing
 ```shell

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_sedov_blast_wave_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl
@@ -0,0 +1,102 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+"""
+    initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+
+The Sedov blast wave setup based on Flash
+- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
+"""
+function initial_condition_sedov_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)
+
+    # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
+    r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+    E = 1.0
+    p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
+    p0_outer = 1.0e-5 # = true Sedov setup
+
+    # Calculate primitive variables
+    rho = 1.0
+    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_sedov_blast_wave
+
+# Get the DG approximation space
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                local_twosided_variables_cons = ["rho"],
+                                local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal,
+                                                                       min)],
+                                max_iterations_newton = 40, # Default parameters are not sufficient to fulfill bounds properly.
+                                newton_tolerances = (1.0e-14, 1.0e-15))
+
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = polydeg, initial_refinement_level = 2,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 300
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 300,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback())
+
+sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
+                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  callback = callbacks);
+summary_callback() # print the timer summary
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