Add subcell limiting for nonconservative systems to the big subcell-limiting branch

.github/workflows/CacheNotebooks.yml

@@ -11,7 +11,7 @@ jobs:
     steps:
       # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
       - name: Checkout caching branch
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
         with:
           ref: tutorial_notebooks
 

.github/workflows/DocPreviewCleanup.yml

@@ -12,7 +12,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout gh-pages branch
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
         with:
           ref: gh-pages
 

.github/workflows/Documenter.yml

@@ -33,7 +33,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: '1.9'

.github/workflows/FormatCheck.yml

@@ -20,7 +20,7 @@ jobs:
         with:
           version: ${{ matrix.julia-version }}
 
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Install JuliaFormatter and format
         # This will use the latest version by default but you can set the version like so:
         #

.github/workflows/Invalidations.yml

@@ -19,12 +19,12 @@ jobs:
     - uses: julia-actions/setup-julia@v1
       with:
         version: '1'
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - uses: julia-actions/julia-buildpkg@v1
     - uses: julia-actions/julia-invalidations@v1
       id: invs_pr
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         ref: ${{ github.event.repository.default_branch }}
     - uses: julia-actions/julia-buildpkg@v1

.github/workflows/ReviewChecklist.yml

@@ -12,7 +12,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Check out repository
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
       - name: Add review checklist
         uses: trixi-framework/add-pr-review-checklist@v1
         with:

.github/workflows/SpellCheck.yml

@@ -8,6 +8,6 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout Actions Repository
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/[email protected].9
+        uses: crate-ci/[email protected].15

.github/workflows/benchmark.yml

@@ -15,7 +15,7 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
         with:
           fetch-depth: 0
       - run: |

.github/workflows/ci.yml

@@ -101,7 +101,7 @@ jobs:
             arch: x64
             trixi_test: threaded
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
@@ -175,7 +175,7 @@ jobs:
       # Instead, we use the more tedious approach described above.
       # At first, we check out the repository and download all artifacts
       # (and list files for debugging).
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: actions/download-artifact@v3
       - run: ls -R
       # Next, we merge the individual coverage files and upload

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.45-pre"
+version = "0.5.46-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

docs/src/callbacks.md

@@ -30,6 +30,11 @@ An example elixir using AMR can be found at [`examples/tree_2d_dgsem/elixir_adve
 The [`AnalysisCallback`](@ref) can be used to analyze the numerical solution, e.g. calculate
 errors or user-specified integrals, and print the results to the screen. The results can also be
 saved in a file. An example can be found at [`examples/tree_2d_dgsem/elixir_euler_vortex.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_vortex.jl).
+Note that the errors (e.g. `L2 error` or `Linf error`) are computed with respect to the initial condition.
+The percentage of the simulation time refers to the ratio of the current time and the final time, i.e. it does
+not consider the maximal number of iterations. So the simulation could finish before 100% are reached.
+Note that, e.g., due to AMR or smaller time step sizes, the simulation can actually take longer than
+the percentage indicates.
 In [Performance metrics of the `AnalysisCallback`](@ref performance-metrics) you can find a detailed
 description of the different performance metrics the `AnalysisCallback` computes.
 

docs/src/overview.md

@@ -55,14 +55,15 @@ different features on different mesh types.
 | Element type                                                 | line, square, cube |     line, quadᵃ, hexᵃ    |             quadᵃ            |     quadᵃ, hexᵃ     |    simplex, quadᵃ, hexᵃ    |
 | Adaptive mesh refinement                                     |          ✅         |             ❌            |               ❌              |          ✅          |               ❌            | [`AMRCallback`](@ref)
 | Solver type                                                  |   [`DGSEM`](@ref)  |      [`DGSEM`](@ref)     |        [`DGSEM`](@ref)       |   [`DGSEM`](@ref)   |       [`DGMulti`](@ref)    |
-| Domain                                                       |      hypercube     |     mapped hypercube     |           arbitrary          |      arbitrary      |       arbitraryᵇ   |
+| Domain                                                       |      hypercube     |     mapped hypercube     |           arbitrary          |      arbitrary      |       arbitrary    |
 | Weak form                                                    |          ✅         |             ✅            |               ✅              |          ✅          |               ✅            | [`VolumeIntegralWeakForm`](@ref)
 | Flux differencing                                            |          ✅         |             ✅            |               ✅              |          ✅          |               ✅            | [`VolumeIntegralFluxDifferencing`](@ref)
 | Shock capturing                                              |          ✅         |             ✅            |               ✅              |          ✅          |               ❌            | [`VolumeIntegralShockCapturingHG`](@ref)
 | Nonconservative equations                                    |          ✅         |             ✅            |               ✅              |          ✅          |               ✅            | e.g., GLM MHD or shallow water equations
+| Parabolic termsᵇ                                             |          ✅         |             ✅            |               ❌              |          ✅          |               ✅            | e.g., [`CompressibleNavierStokesDiffusion2D`](@ref)
 
 ᵃ: quad = quadrilateral, hex = hexahedron
-ᵇ: curved meshes supported for `SBP` and `GaussSBP` approximation types for `VolumeIntegralFluxDifferencing` solvers on quadrilateral and hexahedral `DGMultiMesh`es (non-conservative terms not yet supported)
+ᵇ: Parabolic terms do not currently support adaptivity. 
 
 ## Time integration methods
 

docs/src/parallelization.md

@@ -186,6 +186,11 @@ julia> set_preferences!(
            "libhdf5" => "/path/to/your/libhdf5.so",
            "libhdf5_hl" => "/path/to/your/libhdf5_hl.so", force = true)
 ```
+Alternatively, with HDF5.jl v0.17.1 or higher you can use
+```julia
+julia> using HDF5
+julia> HDF5.API.set_libraries!("/path/to/your/libhdf5.so", "/path/to/your/libhdf5_hl.so")
+```
 For more information see also the
 [documentation of HDF5.jl](https://juliaio.github.io/HDF5.jl/stable/mpi/). In total, you should
 have a file called LocalPreferences.toml in the project directory that contains a section

examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl

@@ -0,0 +1,215 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+prandtl_number() = 0.72
+mu() = 0.01
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(),
+                                                          gradient_variables=GradientVariablesPrimitive())
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
+               volume_integral=VolumeIntegralWeakForm())
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=3, initial_refinement_level=2,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 periodicity=(false, false))
+
+# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D`
+#       since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`)
+#       and by the initial condition (which passes in `CompressibleEulerEquations2D`).
+# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000)
+function initial_condition_navier_stokes_convergence_test(x, t, equations)
+  # Amplitude and shift
+  A = 0.5
+  c = 2.0
+
+  # convenience values for trig. functions
+  pi_x = pi * x[1]
+  pi_y = pi * x[2]
+  pi_t = pi * t
+
+  rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+  v1  = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t)
+  v2  = v1
+  p   = rho^2
+
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations)
+  y = x[2]
+
+  # TODO: parabolic
+  # we currently need to hardcode these parameters until we fix the "combined equation" issue
+  # see also https://github.com/trixi-framework/Trixi.jl/pull/1160
+  inv_gamma_minus_one = inv(equations.gamma - 1)
+  Pr = prandtl_number()
+  mu_ = mu()
+
+  # Same settings as in `initial_condition`
+  # Amplitude and shift
+  A = 0.5
+  c = 2.0
+
+  # convenience values for trig. functions
+  pi_x = pi * x[1]
+  pi_y = pi * x[2]
+  pi_t = pi * t
+
+  # compute the manufactured solution and all necessary derivatives
+  rho    =  c  + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+  rho_t  = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t)
+  rho_x  =  pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t)
+  rho_y  = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t)
+  rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+  rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+
+  v1    =       sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
+  v1_t  = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t)
+  v1_x  =  pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
+  v1_y  =       sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
+  v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
+  v1_xy =  pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
+  v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0)
+                         - A * A * log(y + 2.0) * exp(-A * (y - 1.0))
+                         - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t))
+  v2    = v1
+  v2_t  = v1_t
+  v2_x  = v1_x
+  v2_y  = v1_y
+  v2_xx = v1_xx
+  v2_xy = v1_xy
+  v2_yy = v1_yy
+
+  p    = rho * rho
+  p_t  = 2.0 * rho * rho_t
+  p_x  = 2.0 * rho * rho_x
+  p_y  = 2.0 * rho * rho_y
+  p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x
+  p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y
+
+  # Note this simplifies slightly because the ansatz assumes that v1 = v2
+  E   = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2)
+  E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t
+  E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x
+  E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y
+
+  # Some convenience constants
+  T_const = equations.gamma * inv_gamma_minus_one / Pr
+  inv_rho_cubed = 1.0 / (rho^3)
+
+  # compute the source terms
+  # density equation
+  du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y
+
+  # x-momentum equation
+  du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2
+                                        + 2.0   * rho  * v1 * v1_x
+                                        + rho_y * v1   * v2
+                                        + rho   * v1_y * v2
+                                        + rho   * v1   * v2_y
+    # stress tensor from x-direction
+                      - 4.0 / 3.0 * v1_xx * mu_
+                      + 2.0 / 3.0 * v2_xy * mu_
+                      - v1_yy             * mu_
+                      - v2_xy             * mu_ )
+  # y-momentum equation
+  du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1    * v2
+                                        + rho   * v1_x  * v2
+                                        + rho   * v1    * v2_x
+                                        +         rho_y * v2^2
+                                        + 2.0   * rho   * v2 * v2_y
+    # stress tensor from y-direction
+                      - v1_xy             * mu_
+                      - v2_xx             * mu_
+                      - 4.0 / 3.0 * v2_yy * mu_
+                      + 2.0 / 3.0 * v1_xy * mu_ )
+  # total energy equation
+  du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x)
+              + v2_y * (E + p) + v2 * (E_y + p_y)
+    # stress tensor and temperature gradient terms from x-direction
+                                - 4.0 / 3.0 * v1_xx * v1   * mu_
+                                + 2.0 / 3.0 * v2_xy * v1   * mu_
+                                - 4.0 / 3.0 * v1_x  * v1_x * mu_
+                                + 2.0 / 3.0 * v2_y  * v1_x * mu_
+                                - v1_xy     * v2           * mu_
+                                - v2_xx     * v2           * mu_
+                                - v1_y      * v2_x         * mu_
+                                - v2_x      * v2_x         * mu_
+         - T_const * inv_rho_cubed * (        p_xx * rho   * rho
+                                      - 2.0 * p_x  * rho   * rho_x
+                                      + 2.0 * p    * rho_x * rho_x
+                                      -       p    * rho   * rho_xx ) * mu_
+    # stress tensor and temperature gradient terms from y-direction
+                                - v1_yy     * v1           * mu_
+                                - v2_xy     * v1           * mu_
+                                - v1_y      * v1_y         * mu_
+                                - v2_x      * v1_y         * mu_
+                                - 4.0 / 3.0 * v2_yy * v2   * mu_
+                                + 2.0 / 3.0 * v1_xy * v2   * mu_
+                                - 4.0 / 3.0 * v2_y  * v2_y * mu_
+                                + 2.0 / 3.0 * v1_x  * v2_y * mu_
+         - T_const * inv_rho_cubed * (        p_yy * rho   * rho
+                                      - 2.0 * p_y  * rho   * rho_y
+                                      + 2.0 * p    * rho_y * rho_y
+                                      -       p    * rho   * rho_yy ) * mu_ )
+
+  return SVector(du1, du2, du3, du4)
+end
+
+initial_condition = initial_condition_navier_stokes_convergence_test
+
+# BC types
+velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3])
+heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
+boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom)
+
+boundary_condition_left_right = BoundaryConditionDirichlet(initial_condition_navier_stokes_convergence_test)
+
+# define inviscid boundary conditions
+boundary_conditions = Dict(:x_neg => boundary_condition_left_right,
+                           :x_pos => boundary_condition_left_right,
+                           :y_neg => boundary_condition_slip_wall,
+                           :y_pos => boundary_condition_slip_wall)
+
+# define viscous boundary conditions
+boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_left_right,
+                                     :x_pos => boundary_condition_left_right,
+                                     :y_neg => boundary_condition_top_bottom,
+                                     :y_pos => boundary_condition_top_bottom)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver;
+                                             boundary_conditions=(boundary_conditions, boundary_conditions_parabolic),
+                                             source_terms=source_terms_navier_stokes_convergence_test)
+
+# ###############################################################################
+# # ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5,
+            ode_default_options()..., callback=callbacks)
+summary_callback() # print the timer summary
+