diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml
index a7c4f9caf99..b9bdca5aeac 100644
--- a/.github/workflows/Downgrade.yml
+++ b/.github/workflows/Downgrade.yml
@@ -27,9 +27,9 @@ jobs:
   downgrade_test:
     if: "!contains(github.event.head_commit.message, 'skip ci')"
     # We could also include the Julia version as in
-    # name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - Julia ${{ matrix.version }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    # name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - Julia ${{ matrix.version }} - ${{ matrix.arch }}
     # to be more specific. However, that requires us updating the required CI tests whenever we update Julia.
-    name: Downgrade ${{ matrix.trixi_test }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: Downgrade ${{ matrix.trixi_test }} - ${{ matrix.os }} - ${{ matrix.arch }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
@@ -76,11 +76,11 @@ jobs:
           projects: ., test
       - uses: julia-actions/julia-buildpkg@v1
         env:
-          PYTHON: ""
+          PYTHON: ''
       - name: Run tests without coverage
         uses: julia-actions/julia-runtest@v1
         with:
           coverage: false
         env:
-          PYTHON: ""
+          PYTHON: ''
           TRIXI_TEST: ${{ matrix.trixi_test }}
diff --git a/.github/workflows/GPUCompat.yml b/.github/workflows/GPUCompat.yml
new file mode 100644
index 00000000000..335e1c83c4c
--- /dev/null
+++ b/.github/workflows/GPUCompat.yml
@@ -0,0 +1,86 @@
+name: GPU Package Compatibility
+
+on:
+  pull_request:
+    paths-ignore:
+      - 'AUTHORS.md'
+      - 'CITATION.bib'
+      - 'CONTRIBUTING.md'
+      - 'LICENSE.md'
+      - 'NEWS.md'
+      - 'README.md'
+      - '.zenodo.json'
+      - '.github/workflows/benchmark.yml'
+      - '.github/workflows/CompatHelper.yml'
+      - '.github/workflows/TagBot.yml'
+      - 'benchmark/**'
+      - 'docs/**'
+      - 'utils/**'
+  workflow_dispatch:
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  test:
+    if: "!contains(github.event.head_commit.message, 'skip ci')"
+    name: ${{ matrix.os }} - ${{ matrix.arch }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - version: '1.10'
+            os: ubuntu-latest
+            arch: x64
+          - version: '1.10'
+            os: windows-latest
+            arch: x64
+          # CUDA.jl only supports 64-bit Linux and Windows, see https://github.com/JuliaGPU/CUDA.jl?tab=readme-ov-file#requirements
+    steps:
+      - name: Checkout repository
+        uses: actions/checkout@v4
+
+      - name: Set up Julia
+        uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+
+      - name: Display version info
+        run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)'
+
+      - name: Cache Julia packages
+        uses: julia-actions/cache@v2
+
+      - name: Build project
+        uses: julia-actions/julia-buildpkg@v1
+
+      # Only CUDA.jl is needed for GPU compatibility test now
+      - name: Add CUDA.jl to environment
+        run: |
+          julia --project=. -e '
+          using Pkg; 
+          Pkg.activate(temp=true);
+          Pkg.develop(PackageSpec(path=pwd())); 
+          Pkg.add("CUDA");
+          Pkg.update()'
+
+    #   - name: Add Metal.jl to environment
+    #     run: |
+    #       julia --project=. -e '
+    #       using Pkg;
+    #       Pkg.activate(temp=true);
+    #       Pkg.develop(PackageSpec(path=pwd()));
+    #       Pkg.add("Metal");
+    #       Pkg.update()'
+
+    #   - name: Add AMDGPU.jl to environment
+    #     run: |
+    #       julia --project=. -e '
+    #       using Pkg;
+    #       Pkg.activate(temp=true);
+    #       Pkg.develop(PackageSpec(path=pwd()));
+    #       Pkg.add("AMDGPU");
+    #       Pkg.update()'
diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
index a40e187e82a..3d557915755 100644
--- 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.27.0
+        uses: crate-ci/typos@v1.28.1
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0c636ee8b0b..8862747fefc 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -44,9 +44,9 @@ jobs:
   test:
     if: "!contains(github.event.head_commit.message, 'skip ci')"
     # We could also include the Julia version as in
-    # name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - Julia ${{ matrix.version }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    # name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - Julia ${{ matrix.version }} - ${{ matrix.arch }}
     # to be more specific. However, that requires us updating the required CI tests whenever we update Julia.
-    name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - ${{ matrix.arch }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
@@ -115,27 +115,27 @@ jobs:
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
         env:
-          PYTHON: ""
+          PYTHON: ''
       - name: Run tests without coverage
         uses: julia-actions/julia-runtest@v1
         with:
           coverage: false
         env:
-          PYTHON: ""
+          PYTHON: ''
           TRIXI_TEST: ${{ matrix.trixi_test }}
       - name: Run tests with coverage
         uses: julia-actions/julia-runtest@v1
         with:
           coverage: true
         env:
-          PYTHON: ""
+          PYTHON: ''
           TRIXI_TEST: ${{ matrix.trixi_test }}
       - uses: julia-actions/julia-processcoverage@v1
         with:
           directories: src,examples,ext
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
-          file: ./lcov.info
+          files: ./lcov.info
           flags: unittests
           name: codecov-umbrella
           fail_ci_if_error: true
diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml
index 939eb533ffc..13a01984827 100644
--- a/.github/workflows/downstream.yml
+++ b/.github/workflows/downstream.yml
@@ -46,9 +46,9 @@ jobs:
   test:
     if: "!contains(github.event.head_commit.message, 'skip ci')"
     # We could also include the Julia version as in
-    # name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - Julia ${{ matrix.version }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    # name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - Julia ${{ matrix.version }} - ${{ matrix.arch }}
     # to be more specific. However, that requires us updating the required CI tests whenever we update Julia.
-    name: ${{ matrix.package }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: ${{ matrix.package }} - ${{ matrix.os }} - ${{ matrix.arch }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
diff --git a/NEWS.md b/NEWS.md
index 2f71c45e57e..c187f229519 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,8 @@ for human readability.
 - New time integrator `PairedExplicitRK3`, implementing the third-order paired explicit Runge-Kutta
   method with [Convex.jl](https://github.com/jump-dev/Convex.jl), [ECOS.jl](https://github.com/jump-dev/ECOS.jl),
   and [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) ([#2008])
+- `LobattoLegendreBasis` and related datastructures made fully floating-type general,
+  enabling calculations with higher than double (`Float64`) precision ([#2128])
 
 ## Changes when updating to v0.9 from v0.8.x
 
diff --git a/Project.toml b/Project.toml
index ac1af2c49eb..8a3a4772613 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <m.schlottke-lakemper@hlrs.de>", "Gregor Gassner <ggassner@uni-koeln.de>", "Hendrik Ranocha <mail@ranocha.de>", "Andrew R. Winters <andrew.ross.winters@liu.se>", "Jesse Chan <jesse.chan@rice.edu>"]
-version = "0.9.8-DEV"
+version = "0.9.9-DEV"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -88,7 +88,7 @@ MuladdMacro = "0.2.2"
 NLsolve = "4.5.1"
 Octavian = "0.3.21"
 OffsetArrays = "1.12"
-P4est = "0.4.9"
+P4est = "0.4.12"
 Polyester = "0.7.10"
 PrecompileTools = "1.1"
 Preferences = "1.3"
@@ -108,7 +108,7 @@ StaticArrays = "1.5"
 StrideArrays = "0.1.26"
 StructArrays = "0.6.11"
 SummationByPartsOperators = "0.5.41"
-T8code = "0.7"
+T8code = "0.7.2"
 TimerOutputs = "0.5.7"
 Triangulate = "2.2"
 TriplotBase = "0.1"
diff --git a/README.md b/README.md
index 70c28799f31..543a0339e41 100644
--- a/README.md
+++ b/README.md
@@ -33,6 +33,7 @@ installation and postprocessing procedures. Its features include:
   * 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 and time
+  * Arbitrary floating-point precision
 * Discontinuous Galerkin methods
   * Kinetic energy-preserving and entropy-stable methods based on flux differencing
   * Entropy-stable shock capturing
@@ -260,6 +261,25 @@ To get in touch with the developers,
 [join us on Slack](https://join.slack.com/t/trixi-framework/shared_invite/zt-sgkc6ppw-6OXJqZAD5SPjBYqLd8MU~g)
 or [create an issue](https://github.com/trixi-framework/Trixi.jl/issues/new).
 
+## Participating research groups
+Participating research groups in alphabetical order:
+
+[Applied and Computational Mathematics, RWTH Aachen University :de:](https://www.acom.rwth-aachen.de)
+
+[Applied Mathematics, Department of Mathematics, University of Hamburg :de:](https://www.math.uni-hamburg.de/en/forschung/bereiche/am.html)
+
+[Division of Applied Mathematics, Department of Mathematics, Linköping University :sweden:](https://liu.se/en/employee/andwi94)
+
+[Computational and Applied Mathematics, Rice University :us:](https://jlchan.github.io/)
+
+[High-Performance Computing, Institute of Software Technology, German Aerospace Center (DLR) :de:](https://www.dlr.de/en/sc/about-us/departments/high-performance-computing)
+
+[High-Performance Scientific Computing, University of Augsburg :de:](https://hpsc.math.uni-augsburg.de)
+
+[Numerical Mathematics, Institute of Mathematics, Johannes Gutenberg University Mainz :de:](https://ranocha.de)
+
+[Numerical Simulation, Department of Mathematics and Computer Science, University of Cologne :de:](https://www.mi.uni-koeln.de/NumSim/)
+
 
 ## Acknowledgments
 <p align="center" style="font-size:0;"><!--
@@ -278,7 +298,7 @@ Forschungsgemeinschaft](https://www.dfg.de/) (DFG, German Research Foundation)
 through the following grants:
 * Excellence Strategy EXC 2044-390685587, Mathematics Münster: Dynamics-Geometry-Structure.
 * Research unit FOR 5409 "Structure-Preserving Numerical Methods for Bulk- and
-  Interface Coupling of Heterogeneous Models (SNuBIC)" (project number 463312734).
+  Interface Coupling of Heterogeneous Models ([SNuBIC](https://www.snubic.io))" (project number 463312734).
 * Individual grant no. 528753982.
 
 This project has benefited from funding from the [European Research Council](https://erc.europa.eu)
diff --git a/docs/Project.toml b/docs/Project.toml
index f3a1e66d788..d7cbd579133 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,8 @@
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Changelog = "5217a498-cd5d-4ec6-b8c2-9b85a09b6e3e"
+Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
+ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HOHQMesh = "e4f4c7b8-17cb-445a-93c5-f69190ed6c8c"
@@ -16,7 +18,9 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 [compat]
 CairoMakie = "0.6, 0.7, 0.8, 0.9, 0.10, 0.11"
 Changelog = "1.1"
+Convex = "0.16"
 Documenter = "1"
+ECOS = "1.1.2"
 ForwardDiff = "0.10"
 HOHQMesh = "0.1, 0.2"
 LaTeXStrings = "1.2"
diff --git a/docs/make.jl b/docs/make.jl
index 4eb8fb11b65..c65a20133a8 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -117,6 +117,8 @@ open(joinpath(@__DIR__, "src", "changelog.md"), "w") do io
         println(io, line)
     end
 end
+# Remove temporary file
+rm(joinpath(@__DIR__, "src", "changelog_tmp.md"))
 
 # Make documentation
 makedocs(
diff --git a/docs/src/index.md b/docs/src/index.md
index 88752cc7f11..79da54f3cf9 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -10,6 +10,7 @@
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 [![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT)
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3996439.svg)](https://doi.org/10.5281/zenodo.3996439)
+[![OpenSSF Best Practices](https://www.bestpractices.dev/projects/8695/badge)](https://www.bestpractices.dev/projects/8695)
 
 [**Trixi.jl**](https://github.com/trixi-framework/Trixi.jl)
 is a numerical simulation framework for conservation
@@ -26,6 +27,7 @@ installation and postprocessing procedures. Its features include:
   * 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 and time
+  * Arbitrary floating-point precision
 * Discontinuous Galerkin methods
   * Kinetic energy-preserving and entropy-stable methods based on flux differencing
   * Entropy-stable shock capturing
@@ -52,7 +54,7 @@ installation and postprocessing procedures. Its features include:
   * Hyperbolic diffusion equations for elliptic problems
   * Lattice-Boltzmann equations (D2Q9 and D3Q27 schemes)
   * Shallow water equations
-  * Scalar advection
+  * Several scalar conservation laws (e.g., linear advection, Burgers' equation, LWR traffic flow)
 * Multi-physics simulations
   * [Self-gravitating gas dynamics](https://github.com/trixi-framework/paper-self-gravitating-gas-dynamics)
 * Shared-memory parallelization via multithreading
@@ -98,11 +100,11 @@ If the installed version does not match the current release, please check the
 [Troubleshooting](@ref old-release) section.
 
 The commands above can also be used to update Trixi.jl. A brief list of notable
-changes to Trixi.jl is available in [`NEWS.md`](https://github.com/trixi-framework/Trixi.jl/NEWS.md).
+changes to Trixi.jl is available in the [`Changelog`](@ref).
 
 ### [For developers](@id for-developers)
-If you plan on editing Trixi.jl itself, you can download Trixi.jl to a local folder
-and use the code from the cloned directory:
+If you plan on editing Trixi.jl itself, you can download Trixi.jl locally and use the
+code from the cloned directory:
 ```bash
 git clone git@github.com:trixi-framework/Trixi.jl.git
 cd Trixi.jl
@@ -332,7 +334,7 @@ Schlottke-Lakemper](https://www.uni-augsburg.de/fakultaet/mntf/math/prof/hpsc)
 (University of Augsburg, Germany) and
 [Gregor Gassner](https://www.mi.uni-koeln.de/NumSim/gregor-gassner)
 (University of Cologne, Germany). Together with [Hendrik Ranocha](https://ranocha.de)
-(Johannes Gutenberg University Mainz, Germany) and [Andrew Winters](https://liu.se/en/employee/andwi94)
+(Johannes Gutenberg University Mainz, Germany), [Andrew Winters](https://liu.se/en/employee/andwi94)
 (Linköping University, Sweden), and [Jesse Chan](https://jlchan.github.io) (Rice University, US),
 they are the principal developers of Trixi.jl.
 The full list of contributors can be found under [Authors](@ref).
@@ -349,6 +351,26 @@ To get in touch with the developers,
 or [create an issue](https://github.com/trixi-framework/Trixi.jl/issues/new).
 
 
+## Participating research groups
+Participating research groups in alphabetical order:
+
+[Applied and Computational Mathematics, RWTH Aachen University 🇩🇪](https://www.acom.rwth-aachen.de)
+
+[Applied Mathematics, Department of Mathematics, University of Hamburg 🇩🇪](https://www.math.uni-hamburg.de/en/forschung/bereiche/am.html)
+
+[Division of Applied Mathematics, Department of Mathematics, Linköping University 🇸🇪](https://liu.se/en/employee/andwi94)
+
+[Computational and Applied Mathematics, Rice University 🇺🇸](https://jlchan.github.io/)
+
+[High-Performance Computing, Institute of Software Technology, German Aerospace Center (DLR) 🇩🇪](https://www.dlr.de/en/sc/about-us/departments/high-performance-computing)
+
+[High-Performance Scientific Computing, University of Augsburg 🇩🇪](https://hpsc.math.uni-augsburg.de)
+
+[Numerical Mathematics, Institute of Mathematics, Johannes Gutenberg University Mainz 🇩🇪](https://ranocha.de)
+
+[Numerical Simulation, Department of Mathematics and Computer Science, University of Cologne 🇩🇪](https://www.mi.uni-koeln.de/NumSim/)
+
+
 ## Acknowledgments
 
 ```@raw html
@@ -377,7 +399,7 @@ Forschungsgemeinschaft](https://www.dfg.de/) (DFG, German Research Foundation)
 through the following grants:
 * Excellence Strategy EXC 2044-390685587, Mathematics Münster: Dynamics-Geometry-Structure.
 * Research unit FOR 5409 "Structure-Preserving Numerical Methods for Bulk- and
-  Interface Coupling of Heterogeneous Models (SNuBIC)" (project number 463312734).
+  Interface Coupling of Heterogeneous Models ([SNuBIC](https://www.snubic.io))" (project number 463312734).
 * Individual grant no. 528753982.
 
 This project has benefited from funding from the [European Research Council](https://erc.europa.eu)
diff --git a/docs/src/time_integration.md b/docs/src/time_integration.md
index 600354bb81d..3cdb7fa75ae 100644
--- a/docs/src/time_integration.md
+++ b/docs/src/time_integration.md
@@ -1,11 +1,13 @@
 # [Time integration methods](@id time-integration)
 
+## Methods from OrdinaryDiffEq.jl
+
 Trixi.jl is compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/).
-In particular, explicit Runge-Kutta methods from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
+In particular, [explicit Runge-Kutta methods](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Explicit-Runge-Kutta-Methods) from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
 are tested extensively.
 Interesting classes of time integration schemes are
-- [Explicit low-storage Runge-Kutta methods](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods)
-- [Strong stability preserving methods](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws))
+- [Explicit low-storage Runge-Kutta methods](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Low-Storage-Methods)
+- [Strong stability preserving (SSP) methods](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws))
 
 Some common options for `solve` from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
 are the following. Further documentation can be found in the
@@ -21,7 +23,7 @@ are the following. Further documentation can be found in the
   out of memory or start to swap).
 - You can set the maximal number of time steps via `maxiters=...`.
 - SSP methods and many low-storage methods from OrdinaryDiffEq.jl support
-  `stage_limiter!`s and `step_limiter!`s, e.g., [`PositivityPreservingLimiterZhangShu`](@ref)
+  `stage_limiter!`s and `step_limiter!`s, e.g., [`PositivityPreservingLimiterZhangShu`](@ref) and [`EntropyBoundedLimiter`](@ref)
   from Trixi.jl.
 - If you start Julia with multiple threads and want to use them also in the time 
   integration method from OrdinaryDiffEq.jl, you need to pass the keyword argument
@@ -31,7 +33,7 @@ are the following. Further documentation can be found in the
   For more information on using thread-based parallelism in Trixi.jl, please refer to 
   [Shared-memory parallelization with threads](@ref).
 - If you use error-based step size control (see also the section on
-  [error-based adaptive step sizes](@ref adaptive_step_sizes) together with MPI, you need to
+  [error-based adaptive step sizes](@ref adaptive_step_sizes)) 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).
@@ -42,3 +44,172 @@ are the following. Further documentation can be found in the
     of stages, e.g. to allow for interpolation (dense output), root-finding for continuous callbacks,
     and error-based time step control. In general, you often should not need to worry about this if you
     use Trixi.jl.
+
+## Custom Optimized Schemes
+
+### Stabilized Explicit Runge-Kutta Methods
+
+Optimized explicit schemes aim to maximize the timestep $\Delta t$ for a given simulation setup.
+Formally, this boils down to an optimization problem of the form 
+```math
+\underset{P_{p;S} \, \in \, \mathcal{P}_{p;S}}{\max} \Delta t \text{ such that } \big \vert P_{p;S}(\Delta t \lambda_m) \big \vert \leq 1, \quad  m = 1 , \dots , M \tag{1}
+```
+where $p$ denotes the order of consistency of the scheme, $S$ is the number of stage evaluations and $M$ denotes the number of eigenvalues $\lambda_m$ of the Jacobian matrix $J \coloneqq \frac{\partial \boldsymbol F}{\partial \boldsymbol U}$ of the right-hand side of the [semidiscretized PDE](https://trixi-framework.github.io/Trixi.jl/stable/overview/#overview-semidiscretizations):
+```math
+\dot{\boldsymbol U} = \boldsymbol F(\boldsymbol U) \tag{2} \: .
+```
+In particular, for $S > p$ the Runge-Kutta method includes some free coefficients which may be used to adapt the domain of absolute stability to the problem at hand.
+Since Trixi.jl [supports exact computation of the Jacobian $J$ by means of automatic differentiation](https://trixi-framework.github.io/Trixi.jl/stable/tutorials/differentiable_programming/), we have access to the Jacobian of a given simulation setup.
+For small (say, up to roughly $10^4$ DoF) systems, the spectrum $\boldsymbol \sigma = \left \{ \lambda_m \right \}_{m=1, \dots, M}$ can be computed directly using [`LinearAlgebra.eigvals(J)`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.eigvals).
+For larger systems, we recommend the procedure outlined in section 4.1 of [Doehring et al. (2024)](https://doi.org/10.1016/j.jcp.2024.113223). This approach computes a reduced set of (estimated) eigenvalues $\widetilde{\boldsymbol \sigma}$ around the convex hull of the spectrum by means of the [Arnoldi method](https://github.com/JuliaLinearAlgebra/Arpack.jl).
+
+The optimization problem (1) can be solved using the algorithms described in [Ketcheson, Ahmadia (2012)](http://dx.doi.org/10.2140/camcos.2012.7.247) for a moderate number of stages $S$ or [Doehring, Gassner, Torrilhon (2024)](https://doi.org/10.1007/s10915-024-02478-5) for a large number of stages $S$.
+In Trixi.jl, the former approach is implemented by means of convex optimization using the [Convex.jl](https://github.com/jump-dev/Convex.jl) package.
+
+The resulting stability polynomial $P_{p;S}$ is then used to construct an optimized Runge-Kutta method.
+Trixi.jl implements the [Paired-Explicit Runge-Kutta (PERK)](https://doi.org/10.1016/j.jcp.2019.05.014) method, a low-storage, multirate-ready method with optimized domain of absolute stability.
+
+### Paired-Explicit Runge-Kutta (PERK) Schemes
+
+Paired-Explicit Runge-Kutta (PERK) or `PairedExplicitRK` schemes are an advanced class of numerical methods designed to efficiently solve ODEs.
+In the [original publication](https://doi.org/10.1016/j.jcp.2019.05.014), second-order schemes were introduced, which have been extended to [third](https://doi.org/10.1016/j.jcp.2022.111470) and [fourth](https://doi.org/10.48550/arXiv.2408.05470) order in subsequent works.
+
+By construction, PERK schemes are suited for integrating multirate systems, i.e., systems with varying characteristics speeds throughout the domain.
+Nevertheless, due to their optimized stability properties and low-storage nature, the PERK schemes are also highly efficient when applied as standalone methods. In Trixi.jl, the _standalone_ PERK integrators are implemented such that all stages of the method are active.
+
+#### Tutorial: Using `PairedExplicitRK2`
+
+In this tutorial, we will demonstrate how you can use the second-order PERK time integrator. You need the packages `Convex.jl` and `ECOS.jl`, so be sure they are added to your environment.
+
+First, you need to load the necessary packages:
+
+```@example PERK-example-1
+using Convex, ECOS
+using OrdinaryDiffEq
+using Trixi
+```
+
+Then, define the ODE problem and the semidiscretization setup. For this example, we will use a simple advection problem.
+
+```@example PERK-example-1
+# Define the mesh
+cells_per_dimension = 100
+coordinates_min = 0.0
+coordinates_max = 1.0
+mesh = StructuredMesh(cells_per_dimension, 
+                      coordinates_min, coordinates_max)
+
+# Define the equation and initial condition
+advection_velocity = 1.0
+equations = LinearScalarAdvectionEquation1D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+# Define the solver
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+# Define the semidiscretization
+semi = SemidiscretizationHyperbolic(mesh, 
+                                    equations, initial_condition, 
+                                    solver)
+```
+
+After that, we will define the necessary [callbacks](@ref callbacks-id) for the simulation. Callbacks are used to perform actions at specific points during the integration process.
+
+```@example PERK-example-1
+# Define some standard callbacks
+summary_callback  = SummaryCallback()
+alive_callback    = AliveCallback()
+analysis_callback = AnalysisCallback(semi, interval = 200)
+# For this optimized method we can use a relatively large CFL number
+stepsize_callback = StepsizeCallback(cfl = 2.5)
+
+# Create a CallbackSet to collect all callbacks
+callbacks = CallbackSet(summary_callback, 
+                        alive_callback, 
+                        analysis_callback, 
+                        stepsize_callback)
+```
+
+Now, we define the ODE problem by specifying the time span over which the ODE will be solved. 
+The `tspan` parameter is a tuple `(t_start, t_end)` that defines the start and end times for the simulation. 
+The `semidiscretize` function is used to create the ODE problem from the simulation setup.
+
+```@example PERK-example-1
+# Define the time span
+tspan = (0.0, 1.0)
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, tspan)
+```
+
+Next, we will construct the time integrator. In order to do this, you need the following components:
+
+  - Number of stages: The number of stages $S$ in the Runge-Kutta method. 
+  In this example, we use `6` stages.
+  - Time span (`tspan`): A tuple `(t_start, t_end)` that defines the time span over which the ODE will be solved. 
+  This defines the bounds for the bisection routine for the optimal timestep $\Delta t$ used in calculating the polynomial coefficients at optimization stage. 
+  This variable is already defined in step 5.
+  - Semidiscretization (`semi`): The semidiscretization setup that includes the mesh, equations, initial condition, and solver. In this example, this variable is already defined in step 3.
+  In the background, we compute from `semi` the Jacobian $J$ evaluated at the initial condition using [`jacobian_ad_forward`](https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/#Trixi.jacobian_ad_forward-Tuple{Trixi.AbstractSemidiscretization}).
+  This is then followed by the computation of the spectrum $\boldsymbol \sigma(J)$ using `LinearAlgebra.eigvals`.
+  Equipped with the spectrum, the optimal stability polynomial is computed, from which the corresponding Runge-Kutta method is derived. Other constructors (if the coefficients $\boldsymbol{\alpha}$ of the stability polynomial are already available, or if a reduced spectrum $\widetilde{\boldsymbol{\sigma}}$ should be used) are discussed below.
+
+```@example PERK-example-1
+# Construct second order-explicit Runge-Kutta method with 6 stages for given simulation setup (`semi`)
+# `tspan` provides the bounds for the bisection routine that is used to calculate the maximum timestep.
+ode_algorithm = Trixi.PairedExplicitRK2(6, tspan, semi)
+```
+
+With everything set up, you can now use `Trixi.solve` to solve the ODE problem. The `solve` function takes the ODE problem, the time integrator, and some options such as the time step (`dt`), whether to save every step (`save_everystep`), and the callbacks.
+
+```@example PERK-example-1
+# Solve the ODE problem using PERK2
+sol = Trixi.solve(ode, ode_algorithm,
+                  dt = 1.0, # overwritten by `stepsize_callback`
+                  save_everystep = false, callback = callbacks)
+```
+
+##### Advanced constructors:
+There are two additional constructors for the `PairedExplicitRK2` method besides the one taking in a semidiscretization `semi`:
+  - `PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString)` constructs a `num_stages`-stage method from the given optimal monomial_coeffs $\boldsymbol \alpha$.
+  These are expected to be present in the provided directory in the form of a `gamma_<S>.txt` file, where `<S>` is the number of stages `num_stages`.
+  This constructor is useful when the optimal coefficients cannot be obtained using the optimization routine by Ketcheson and Ahmadia, possibly due to a large number of stages $S$.
+  - `PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64})` constructs a `num_stages`-stage using the optimization approach by Ketcheson and Ahmadia for the (reduced) spectrum `eig_vals`.
+  The use-case for this constructor would be a large system, for which the computation of all eigenvalues is infeasible.
+
+#### Automatic computation of a stable CFL Number
+
+In the previous tutorial the CFL number was set manually to $2.5$.
+To avoid the manual trial-and-error process behind this, instantiations of `AbstractPairedExplicitRK` methods can automatically compute the stable CFL number for a given simulation setup using the [`Trixi.calculate_cfl`](@ref) function.
+When constructing the time integrator from a semidiscretization `semi`, 
+```@example PERK-example-1
+# Construct third-order paired-explicit Runge-Kutta method with 8 stages for given simulation setup.                               
+ode_algorithm = Trixi.PairedExplicitRK3(8, tspan, semi)
+```
+the maximum timestep `dt` is stored by the `ode_algorithm`.
+This can then be used to compute the stable CFL number for the given simulation setup:
+```@example PERK-example-1
+cfl_number = Trixi.calculate_cfl(ode_algorithm, ode)
+```
+For nonlinear problems, the spectrum will in general change over the course of the simulation.
+Thus, it is often necessary to reduce the optimal `cfl_number` by a safety factor:
+```@example PERK-example-1
+# For non-linear problems, the CFL number should be reduced by a safety factor
+# as the spectrum changes (in general) over the course of a simulation
+stepsize_callback = StepsizeCallback(cfl = 0.85 * cfl_number)
+```
+If the optimal monomial coefficients are precomputed, the user needs to provide the obtained maximum timestep `dt_opt` from the optimization at construction stage.
+The corresponding constructor has signature
+```julia
+PairedExplicitRK3(num_stages, base_path_a_coeffs::AbstractString,
+                  dt_opt = nothing; cS2 = 1.0f0)
+```
+Then, the stable CFL number can be computed as described above.
+
+#### Currently implemented PERK methods
+
+##### Single/Standalone methods
+
+- [`Trixi.PairedExplicitRK2`](@ref): Second-order PERK method with at least two stages.
+- [`Trixi.PairedExplicitRK3`](@ref): Third-order PERK method with at least three stages.
\ No newline at end of file
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock.jl
new file mode 100644
index 00000000000..feec3d534d0
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock.jl
@@ -0,0 +1,183 @@
+using OrdinaryDiffEq
+using Trixi
+
+# This is the classic 1D viscous shock wave problem with analytical solution 
+# for a special value of the Prandtl number.
+# The original references are:
+#
+# - R. Becker (1922)
+#   Stoßwelle und Detonation.
+#   [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605)
+#
+#   English translations:
+#   Impact waves and detonation. Part I.
+#   https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf
+#   Impact waves and detonation. Part II.
+#   https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf
+#
+# - M. Morduchow, P. A. Libby (1949)
+#   On a Complete Solution of the One-Dimensional Flow Equations 
+#   of a Viscous, Head-Conducting, Compressible Gas
+#   [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882)
+#
+#
+# The particular problem considered here is described in
+# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) 
+#   Entropy in self-similar shock profiles
+#   [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+
+### Fixed parameters ###
+
+# Special value for which nonlinear solver can be omitted
+# Corresponds essentially to fixing the Mach number
+alpha = 0.5
+# We want kappa = cp * mu = mu_bar to ensure constant enthalpy
+prandtl_number() = 3 / 4
+
+### Free choices: ###
+gamma() = 5 / 3
+
+# In Margolin et al., the Navier-Stokes equations are given for an 
+# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu 
+mu_isotropic() = 0.15
+mu_bar() = mu_isotropic() / (gamma() - 1) # Re-scaled viscosity
+
+rho_0() = 1
+v() = 1 # Shock speed
+
+domain_length = 4.0
+
+### Derived quantities ###
+
+Ma() = 2 / sqrt(3 - gamma()) # Mach number for alpha = 0.5
+c_0() = v() / Ma() # Speed of sound ahead of the shock
+
+# From constant enthalpy condition
+p_0() = c_0()^2 * rho_0() / gamma()
+
+l() = mu_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale
+
+"""
+    initial_condition_viscous_shock(x, t, equations)
+
+Classic 1D viscous shock wave problem with analytical solution 
+for a special value of the Prandtl number.
+The version implemented here is described in
+- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
+  Entropy in self-similar shock profiles
+  [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+"""
+function initial_condition_viscous_shock(x, t, equations)
+    y = x[1] - v() * t # Translated coordinate
+
+    # Coordinate transformation. See eq. (33) in Margolin et al. (2017)
+    chi = 2 * exp(y / (2 * l()))
+
+    w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2))
+
+    rho = rho_0() / w
+    u = v() * (1 - w)
+    p = p_0() * 1 / w * (1 + (gamma() - 1) / 2 * Ma()^2 * (1 - w^2))
+
+    return prim2cons(SVector(rho, u, 0, p), equations)
+end
+initial_condition = initial_condition_viscous_shock
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+equations = CompressibleEulerEquations2D(gamma())
+
+# Trixi implements the stress tensor in deviatoric form, thus we need to 
+# convert the "isotropic viscosity" to the "deviatoric viscosity"
+mu_deviatoric() = mu_bar() * 3 / 4
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu_deviatoric(),
+                                                          Prandtl = prandtl_number(),
+                                                          gradient_variables = GradientVariablesPrimitive())
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hlle)
+
+coordinates_min = (-domain_length / 2, -domain_length / 2)
+coordinates_max = (domain_length / 2, domain_length / 2)
+
+trees_per_dimension = (8, 2)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 0,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (false, true))
+
+### Inviscid boundary conditions ###
+
+# Prescribe pure influx based on initial conditions
+function boundary_condition_inflow(u_inner, normal_direction::AbstractVector, x, t,
+                                   surface_flux_function,
+                                   equations::CompressibleEulerEquations2D)
+    u_cons = initial_condition_viscous_shock(x, t, equations)
+    flux = Trixi.flux(u_cons, normal_direction, equations)
+
+    return flux
+end
+
+# Completely free outflow
+function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
+                                    surface_flux_function,
+                                    equations::CompressibleEulerEquations2D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, normal_direction, equations)
+
+    return flux
+end
+
+boundary_conditions = Dict(:x_neg => boundary_condition_inflow,
+                           :x_pos => boundary_condition_outflow)
+
+### Viscous boundary conditions ###
+# For the viscous BCs, we use the known analytical solution
+velocity_bc = NoSlip() do x, t, equations_parabolic
+    Trixi.velocity(initial_condition_viscous_shock(x,
+                                                   t,
+                                                   equations_parabolic),
+                   equations_parabolic)
+end
+
+heat_bc = Isothermal() do x, t, equations_parabolic
+    Trixi.temperature(initial_condition_viscous_shock(x,
+                                                      t,
+                                                      equations_parabolic),
+                      equations_parabolic)
+end
+
+boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc)
+
+boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_parabolic,
+                                     :x_pos => boundary_condition_parabolic)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions_parabolic))
+
+###############################################################################
+# 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-3, ode_default_options()..., callback = callbacks)
+
+summary_callback() # print the timer summary
diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_viscous_shock.jl
new file mode 100644
index 00000000000..43cdafa3753
--- /dev/null
+++ b/examples/p4est_3d_dgsem/elixir_navierstokes_viscous_shock.jl
@@ -0,0 +1,183 @@
+using OrdinaryDiffEq
+using Trixi
+
+# This is the classic 1D viscous shock wave problem with analytical solution 
+# for a special value of the Prandtl number.
+# The original references are:
+#
+# - R. Becker (1922)
+#   Stoßwelle und Detonation.
+#   [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605)
+#
+#   English translations:
+#   Impact waves and detonation. Part I.
+#   https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf
+#   Impact waves and detonation. Part II.
+#   https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf
+#
+# - M. Morduchow, P. A. Libby (1949)
+#   On a Complete Solution of the One-Dimensional Flow Equations 
+#   of a Viscous, Head-Conducting, Compressible Gas
+#   [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882)
+#
+#
+# The particular problem considered here is described in
+# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) 
+#   Entropy in self-similar shock profiles
+#   [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+
+### Fixed parameters ###
+
+# Special value for which nonlinear solver can be omitted
+# Corresponds essentially to fixing the Mach number
+alpha = 0.5
+# We want kappa = cp * mu = mu_bar to ensure constant enthalpy
+prandtl_number() = 3 / 4
+
+### Free choices: ###
+gamma() = 5 / 3
+
+# In Margolin et al., the Navier-Stokes equations are given for an 
+# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu 
+mu_isotropic() = 0.15
+mu_bar() = mu_isotropic() / (gamma() - 1) # Re-scaled viscosity
+
+rho_0() = 1
+v() = 1 # Shock speed
+
+domain_length = 4.0
+
+### Derived quantities ###
+
+Ma() = 2 / sqrt(3 - gamma()) # Mach number for alpha = 0.5
+c_0() = v() / Ma() # Speed of sound ahead of the shock
+
+# From constant enthalpy condition
+p_0() = c_0()^2 * rho_0() / gamma()
+
+l() = mu_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale
+
+"""
+    initial_condition_viscous_shock(x, t, equations)
+
+Classic 1D viscous shock wave problem with analytical solution 
+for a special value of the Prandtl number.
+The version implemented here is described in
+- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
+  Entropy in self-similar shock profiles
+  [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+"""
+function initial_condition_viscous_shock(x, t, equations)
+    y = x[1] - v() * t # Translated coordinate
+
+    # Coordinate transformation. See eq. (33) in Margolin et al. (2017)
+    chi = 2 * exp(y / (2 * l()))
+
+    w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2))
+
+    rho = rho_0() / w
+    u = v() * (1 - w)
+    p = p_0() * 1 / w * (1 + (gamma() - 1) / 2 * Ma()^2 * (1 - w^2))
+
+    return prim2cons(SVector(rho, u, 0, 0, p), equations)
+end
+initial_condition = initial_condition_viscous_shock
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+equations = CompressibleEulerEquations3D(gamma())
+
+# Trixi implements the stress tensor in deviatoric form, thus we need to 
+# convert the "isotropic viscosity" to the "deviatoric viscosity"
+mu_deviatoric() = mu_bar() * 3 / 4
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu_deviatoric(),
+                                                          Prandtl = prandtl_number(),
+                                                          gradient_variables = GradientVariablesPrimitive())
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hlle)
+
+coordinates_min = (-domain_length / 2, -domain_length / 2, -domain_length / 2)
+coordinates_max = (domain_length / 2, domain_length / 2, domain_length / 2)
+
+trees_per_dimension = (8, 2, 2)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 0,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (false, true, true))
+
+### Inviscid boundary conditions ###
+
+# Prescribe pure influx based on initial conditions
+function boundary_condition_inflow(u_inner, normal_direction::AbstractVector, x, t,
+                                   surface_flux_function,
+                                   equations::CompressibleEulerEquations3D)
+    u_cons = initial_condition_viscous_shock(x, t, equations)
+    flux = Trixi.flux(u_cons, normal_direction, equations)
+
+    return flux
+end
+
+# Completely free outflow
+function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
+                                    surface_flux_function,
+                                    equations::CompressibleEulerEquations3D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, normal_direction, equations)
+
+    return flux
+end
+
+boundary_conditions = Dict(:x_neg => boundary_condition_inflow,
+                           :x_pos => boundary_condition_outflow)
+
+### Viscous boundary conditions ###
+# For the viscous BCs, we use the known analytical solution
+velocity_bc = NoSlip() do x, t, equations_parabolic
+    Trixi.velocity(initial_condition_viscous_shock(x,
+                                                   t,
+                                                   equations_parabolic),
+                   equations_parabolic)
+end
+
+heat_bc = Isothermal() do x, t, equations_parabolic
+    Trixi.temperature(initial_condition_viscous_shock(x,
+                                                      t,
+                                                      equations_parabolic),
+                      equations_parabolic)
+end
+
+boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc)
+
+boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_parabolic,
+                                     :x_pos => boundary_condition_parabolic)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions_parabolic))
+
+###############################################################################
+# 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-3, ode_default_options()..., callback = callbacks)
+
+summary_callback() # print the timer summary
diff --git a/examples/structured_1d_dgsem/elixir_advection_float128.jl b/examples/structured_1d_dgsem/elixir_advection_float128.jl
new file mode 100644
index 00000000000..7f72eba58e3
--- /dev/null
+++ b/examples/structured_1d_dgsem/elixir_advection_float128.jl
@@ -0,0 +1,60 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+using Quadmath
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+# See https://github.com/JuliaMath/Quadmath.jl
+RealT = Float128
+
+advection_velocity = 4 / 3 # Does not need to be in higher precision
+equations = LinearScalarAdvectionEquation1D(advection_velocity)
+
+solver = DGSEM(RealT = RealT, polydeg = 13, surface_flux = flux_lax_friedrichs)
+
+# CARE: Important to use higher precision datatype for coordinates
+# as these are used for type promotion of the mesh (points etc.)
+coordinates_min = (-one(RealT),)
+coordinates_max = (one(RealT),)
+cells_per_dimension = (1,)
+
+# `StructuredMesh` infers datatype from coordinates
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# CARE: Important to use higher precision datatype in specification of final time
+tspan = (zero(RealT), one(RealT))
+
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+analysis_callback = AnalysisCallback(semi, interval = 100,
+                                     extra_analysis_errors = (:conservation_error,))
+
+# cfl does not need to be in higher precision
+stepsize_callback = StepsizeCallback(cfl = 0.25)
+
+callbacks = CallbackSet(summary_callback,
+                        stepsize_callback,
+                        analysis_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, Feagin14(),
+            # Turn off adaptivity to avoid setting very small tolerances
+            adaptive = false,
+            dt = 42, # `dt` does not need to be in higher precision
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/structured_1d_dgsem/elixir_burgers_perk3.jl b/examples/structured_1d_dgsem/elixir_burgers_perk3.jl
index bf91fde74ea..5ff3aff3678 100644
--- a/examples/structured_1d_dgsem/elixir_burgers_perk3.jl
+++ b/examples/structured_1d_dgsem/elixir_burgers_perk3.jl
@@ -1,9 +1,9 @@
 # Convex and ECOS are imported because they are used for finding the optimal time step and optimal 
-# monomial coefficients in the stability polynomial of P-ERK time integrators.
+# monomial coefficients in the stability polynomial of PERK time integrators.
 using Convex, ECOS
 
 # NLsolve is imported to solve the system of nonlinear equations to find the coefficients
-# in the Butcher tableau in the third order P-ERK time integrator.
+# in the Butcher tableau in the third order PERK time integrator.
 using NLsolve
 
 using OrdinaryDiffEq
diff --git a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
index 3cb4735e132..6b7dd007128 100644
--- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
+++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
@@ -19,7 +19,8 @@ mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
 # Green light that at x = 0 which switches at t = 0 from red to green.
 # To the left there are cars bumper to bumper, to the right there are no cars.
 function initial_condition_greenlight(x, t, equation::TrafficFlowLWREquations1D)
-    scalar = x[1] < 0.0 ? 1.0 : 0.0
+    RealT = eltype(x)
+    scalar = x[1] < 0 ? one(RealT) : zero(RealT)
 
     return SVector(scalar)
 end
@@ -29,7 +30,8 @@ end
 
 # Assume that there are always cars waiting at the left 
 function inflow(x, t, equations::TrafficFlowLWREquations1D)
-    return initial_condition_greenlight(coordinates_min, t, equations)
+    # -1.0 = coordinates_min
+    return initial_condition_greenlight(-1.0, t, equations)
 end
 boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
 
diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl
index 6ef227d882b..92abc0b765c 100644
--- a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl
+++ b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl
@@ -152,7 +152,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl
index d19b9230c5c..97250074e31 100644
--- a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl
+++ b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl
@@ -93,7 +93,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
             save_everystep = false, callback = callbacks);
 
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_advection_basic.jl b/examples/t8code_2d_dgsem/elixir_advection_basic.jl
index 648661898ec..36fececcc7a 100644
--- a/examples/t8code_2d_dgsem/elixir_advection_basic.jl
+++ b/examples/t8code_2d_dgsem/elixir_advection_basic.jl
@@ -55,7 +55,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary
 summary_callback()
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_advection_extended.jl b/examples/t8code_2d_dgsem/elixir_advection_extended.jl
index 27f271e9870..95793a309ef 100644
--- a/examples/t8code_2d_dgsem/elixir_advection_extended.jl
+++ b/examples/t8code_2d_dgsem/elixir_advection_extended.jl
@@ -1,4 +1,3 @@
-
 using OrdinaryDiffEq
 using Trixi
 
@@ -83,7 +82,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary
 summary_callback()
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl
index 66e38555d75..76b82a3a712 100644
--- a/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl
+++ b/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl
@@ -90,7 +90,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary
 summary_callback()
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_advection_restart.jl b/examples/t8code_2d_dgsem/elixir_advection_restart.jl
index 4f5d29de318..bd6f4c663de 100644
--- a/examples/t8code_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/t8code_2d_dgsem/elixir_advection_restart.jl
@@ -1,4 +1,3 @@
-
 using OrdinaryDiffEq
 using Trixi
 
@@ -41,7 +40,3 @@ load_timestep!(integrator, restart_filename)
 
 sol = solve!(integrator)
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_advection_restart_amr.jl b/examples/t8code_2d_dgsem/elixir_advection_restart_amr.jl
index 56a2e5b2abc..3e1a35581d5 100644
--- a/examples/t8code_2d_dgsem/elixir_advection_restart_amr.jl
+++ b/examples/t8code_2d_dgsem/elixir_advection_restart_amr.jl
@@ -1,4 +1,3 @@
-
 using OrdinaryDiffEq
 using Trixi
 
@@ -60,7 +59,3 @@ load_timestep!(integrator, restart_filename)
 
 sol = solve!(integrator)
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl
index f13f0fae05a..aed6c755d5c 100644
--- a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl
+++ b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl
@@ -76,7 +76,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary.
 summary_callback()
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl
index 349678f1637..b2d49e3ccfe 100644
--- a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl
+++ b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl
@@ -91,7 +91,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_euler_sedov.jl b/examples/t8code_2d_dgsem/elixir_euler_sedov.jl
index e8ec3a7a9e6..fae7d818ad8 100644
--- a/examples/t8code_2d_dgsem/elixir_euler_sedov.jl
+++ b/examples/t8code_2d_dgsem/elixir_euler_sedov.jl
@@ -98,7 +98,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl b/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl
index a7d10006593..2a0c4a6ee20 100644
--- a/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl
+++ b/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl
@@ -70,7 +70,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl
index 9f9ad16dc9b..19502f1ce0e 100644
--- a/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl
+++ b/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl
@@ -94,7 +94,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
index a3366caa317..b8dbf94ac0c 100644
--- a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
+++ b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -1,4 +1,3 @@
-
 using OrdinaryDiffEq
 using Trixi
 
@@ -108,7 +107,3 @@ 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
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl
index 2fdb624593c..0fe308701b5 100644
--- a/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl
+++ b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl
@@ -81,7 +81,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
             save_everystep = false, callback = callbacks);
 summary_callback() # print the timer summary
 println("Number of gravity subcycles: ", semi.gravity_counter.ncalls_since_readout)
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl
index 04c36dd8642..e184cb3fd05 100644
--- a/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl
+++ b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl
@@ -57,7 +57,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl
index fa8d27a01b4..5dcea4d3150 100644
--- a/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl
+++ b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl
@@ -137,7 +137,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl
index 0a937664493..688ddb2dbb5 100644
--- a/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl
+++ b/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl
@@ -59,7 +59,3 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav
 sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8,
             ode_default_options()..., callback = callbacks);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_advection_amr.jl b/examples/t8code_3d_dgsem/elixir_advection_amr.jl
index 7db6bbc3366..44c796163c1 100644
--- a/examples/t8code_3d_dgsem/elixir_advection_amr.jl
+++ b/examples/t8code_3d_dgsem/elixir_advection_amr.jl
@@ -77,7 +77,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl
index 45daaf0fbb2..fdea156e222 100644
--- a/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl
+++ b/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl
@@ -112,7 +112,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_advection_basic.jl b/examples/t8code_3d_dgsem/elixir_advection_basic.jl
index 415a9e78ba2..ae97a73d182 100644
--- a/examples/t8code_3d_dgsem/elixir_advection_basic.jl
+++ b/examples/t8code_3d_dgsem/elixir_advection_basic.jl
@@ -65,7 +65,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary
 summary_callback()
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_advection_nonconforming.jl b/examples/t8code_3d_dgsem/elixir_advection_nonconforming.jl
index f12c760bb1b..0755a76ef45 100644
--- a/examples/t8code_3d_dgsem/elixir_advection_nonconforming.jl
+++ b/examples/t8code_3d_dgsem/elixir_advection_nonconforming.jl
@@ -87,7 +87,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary
 summary_callback()
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_advection_restart.jl b/examples/t8code_3d_dgsem/elixir_advection_restart.jl
index 83290c6014a..4c72a0b3ee3 100644
--- a/examples/t8code_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/t8code_3d_dgsem/elixir_advection_restart.jl
@@ -38,7 +38,3 @@ load_timestep!(integrator, restart_filename)
 
 sol = solve!(integrator)
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl
index 9b8e8825672..400a3cca9c0 100644
--- a/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl
+++ b/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl
@@ -96,7 +96,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
 
 # Print the timer summary
 summary_callback()
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_euler_ec.jl b/examples/t8code_3d_dgsem/elixir_euler_ec.jl
index 2e13f7edd57..e75b0f69636 100644
--- a/examples/t8code_3d_dgsem/elixir_euler_ec.jl
+++ b/examples/t8code_3d_dgsem/elixir_euler_ec.jl
@@ -87,7 +87,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl
index daf431a86ef..d45de658cc0 100644
--- a/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl
+++ b/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl
@@ -114,7 +114,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl
index 19cdea7ee3c..d24512a4cdd 100644
--- a/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl
+++ b/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl
@@ -102,7 +102,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), #maxiters=1
             dt = 1.0, # 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
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_euler_sedov.jl b/examples/t8code_3d_dgsem/elixir_euler_sedov.jl
index 9f42ad24ac7..f897249ed2e 100644
--- a/examples/t8code_3d_dgsem/elixir_euler_sedov.jl
+++ b/examples/t8code_3d_dgsem/elixir_euler_sedov.jl
@@ -100,7 +100,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl
index 1346f603608..4b87b646df9 100644
--- a/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl
+++ b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl
@@ -116,7 +116,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl
index 612b1640239..ce1662c8e50 100644
--- a/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl
+++ b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl
@@ -66,7 +66,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl
index 106b4821144..83641f085c1 100644
--- a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl
+++ b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -1,4 +1,3 @@
-
 using OrdinaryDiffEq
 using Trixi
 
@@ -110,7 +109,3 @@ 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);
 summary_callback() # print the timer summary
-
-# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh)
diff --git a/examples/tree_1d_dgsem/elixir_advection_doublefloat.jl b/examples/tree_1d_dgsem/elixir_advection_doublefloat.jl
new file mode 100644
index 00000000000..8405d0125d0
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_advection_doublefloat.jl
@@ -0,0 +1,63 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+using DoubleFloats
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+# See https://github.com/JuliaMath/DoubleFloats.jl
+RealT = Double64
+
+advection_velocity = 4 / 3 # Does not need to be in higher precision
+equations = LinearScalarAdvectionEquation1D(advection_velocity)
+
+solver = DGSEM(RealT = RealT, polydeg = 7, surface_flux = flux_lax_friedrichs)
+
+# CARE: Important to use higher precision datatype for coordinates
+# as these are used for type promotion of the mesh (points etc.)
+coordinates_min = -one(RealT) # minimum coordinate
+coordinates_max = one(RealT) # maximum coordinate
+
+# For `TreeMesh` the datatype has to be specified explicitly,
+# i.e., is not inferred from the coordinates.
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 30_000,
+                RealT = RealT)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# CARE: Important to use higher precision datatype in specification of final time
+tspan = (zero(RealT), one(RealT))
+
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+analysis_callback = AnalysisCallback(semi, interval = 100,
+                                     extra_analysis_errors = (:conservation_error,))
+
+# cfl does not need to be in higher precision
+stepsize_callback = StepsizeCallback(cfl = 1.4)
+
+callbacks = CallbackSet(summary_callback,
+                        stepsize_callback,
+                        analysis_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, DP8(),
+            # Turn off adaptivity to avoid setting very small tolerances
+            adaptive = false,
+            dt = 42, # `dt` does not need to be in higher precision
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2.jl b/examples/tree_1d_dgsem/elixir_advection_perk2.jl
index 69b587b42ec..86b173d6564 100644
--- a/examples/tree_1d_dgsem/elixir_advection_perk2.jl
+++ b/examples/tree_1d_dgsem/elixir_advection_perk2.jl
@@ -1,6 +1,6 @@
 
 # Convex and ECOS are imported because they are used for finding the optimal time step and optimal 
-# monomial coefficients in the stability polynomial of P-ERK time integrators.
+# monomial coefficients in the stability polynomial of PERK time integrators.
 using Convex, ECOS
 
 using OrdinaryDiffEq
diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl
new file mode 100644
index 00000000000..df6a453ac83
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl
@@ -0,0 +1,176 @@
+using OrdinaryDiffEq
+using Trixi
+
+# This is the classic 1D viscous shock wave problem with analytical solution 
+# for a special value of the Prandtl number.
+# The original references are:
+#
+# - R. Becker (1922)
+#   Stoßwelle und Detonation.
+#   [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605)
+#
+#   English translations:
+#   Impact waves and detonation. Part I.
+#   https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf
+#   Impact waves and detonation. Part II.
+#   https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf
+#
+# - M. Morduchow, P. A. Libby (1949)
+#   On a Complete Solution of the One-Dimensional Flow Equations 
+#   of a Viscous, Head-Conducting, Compressible Gas
+#   [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882)
+#
+#
+# The particular problem considered here is described in
+# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017) 
+#   Entropy in self-similar shock profiles
+#   [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+
+### Fixed parameters ###
+
+# Special value for which nonlinear solver can be omitted
+# Corresponds essentially to fixing the Mach number
+alpha = 0.5
+# We want kappa = cp * mu = mu_bar to ensure constant enthalpy
+prandtl_number() = 1
+
+### Free choices: ###
+gamma() = 5 / 3
+
+mu() = 0.15
+mu_bar() = mu() / (gamma() - 1) # Re-scaled viscosity
+
+rho_0() = 1
+v() = 1 # Shock speed
+
+domain_length = 4.0
+
+### Derived quantities ###
+
+Ma() = 2 / sqrt(3 - gamma()) # Mach number for alpha = 0.5
+c_0() = v() / Ma() # Speed of sound ahead of the shock
+
+# From constant enthalpy condition
+p_0() = c_0()^2 * rho_0() / gamma()
+
+l() = mu_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale
+
+"""
+    initial_condition_viscous_shock(x, t, equations)
+
+Classic 1D viscous shock wave problem with analytical solution 
+for a special value of the Prandtl number.
+The version implemented here is described in
+- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
+  Entropy in self-similar shock profiles
+  [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+"""
+function initial_condition_viscous_shock(x, t, equations)
+    y = x[1] - v() * t # Translated coordinate
+
+    # Coordinate transformation. See eq. (33) in Margolin et al. (2017)
+    chi = 2 * exp(y / (2 * l()))
+
+    w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2))
+
+    rho = rho_0() / w
+    u = v() * (1 - w)
+    p = p_0() * 1 / w * (1 + (gamma() - 1) / 2 * Ma()^2 * (1 - w^2))
+
+    return prim2cons(SVector(rho, u, p), equations)
+end
+initial_condition = initial_condition_viscous_shock
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+equations = CompressibleEulerEquations1D(gamma())
+equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu = mu_bar(),
+                                                          Prandtl = prandtl_number(),
+                                                          gradient_variables = GradientVariablesPrimitive())
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hlle)
+
+coordinates_min = -domain_length / 2
+coordinates_max = domain_length / 2
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                periodicity = false,
+                n_cells_max = 30_000)
+
+### Inviscid boundary conditions ###
+
+# Prescribe pure influx based on initial conditions
+function boundary_condition_inflow(u_inner, orientation::Integer, normal_direction, x, t,
+                                   surface_flux_function,
+                                   equations::CompressibleEulerEquations1D)
+    u_cons = initial_condition_viscous_shock(x, t, equations)
+    flux = Trixi.flux(u_cons, orientation, equations)
+
+    return flux
+end
+
+# Completely free outflow
+function boundary_condition_outflow(u_inner, orientation::Integer, normal_direction, x, t,
+                                    surface_flux_function,
+                                    equations::CompressibleEulerEquations1D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, orientation, equations)
+
+    return flux
+end
+
+boundary_conditions = (; x_neg = boundary_condition_inflow,
+                       x_pos = boundary_condition_outflow)
+
+### Viscous boundary conditions ###
+# For the viscous BCs, we use the known analytical solution
+velocity_bc = NoSlip() do x, t, equations_parabolic
+    Trixi.velocity(initial_condition_viscous_shock(x,
+                                                   t,
+                                                   equations_parabolic),
+                   equations_parabolic)
+end
+
+heat_bc = Isothermal() do x, t, equations_parabolic
+    Trixi.temperature(initial_condition_viscous_shock(x,
+                                                      t,
+                                                      equations_parabolic),
+                      equations_parabolic)
+end
+
+boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc)
+
+boundary_conditions_parabolic = (; x_neg = boundary_condition_parabolic,
+                                 x_pos = boundary_condition_parabolic)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions_parabolic))
+
+###############################################################################
+# 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 = 100)
+
+analysis_interval = 1000
+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-3, ode_default_options()..., callback = callbacks)
+
+summary_callback() # print the timer summary
diff --git a/examples/tree_3d_dgsem/elixir_euler_convergence_pure_fv.jl b/examples/tree_3d_dgsem/elixir_euler_convergence_pure_fv.jl
index 063053c1255..4789b46dacc 100644
--- a/examples/tree_3d_dgsem/elixir_euler_convergence_pure_fv.jl
+++ b/examples/tree_3d_dgsem/elixir_euler_convergence_pure_fv.jl
@@ -53,8 +53,3 @@ 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);
 summary_callback() # print the timer summary
-
-# This elixir is also tested with `T8codeMesh` which needs to be finalized
-# explicitly. Call finalize to make sure MPI related objects in t8code are
-# released before `MPI` finalizes.
-!isinteractive() && finalize(mesh);
diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl
index 9b897436dbb..a83ac0a524f 100644
--- a/ext/TrixiConvexECOSExt.jl
+++ b/ext/TrixiConvexECOSExt.jl
@@ -1,7 +1,7 @@
 # Package extension for adding Convex-based features to Trixi.jl
 module TrixiConvexECOSExt
 
-# Required for coefficient optimization in P-ERK scheme integrators
+# Required for coefficient optimization in PERK scheme integrators
 if isdefined(Base, :get_extension)
     using Convex: MOI, solve!, Variable, minimize, evaluate
     using ECOS: Optimizer
diff --git a/ext/TrixiNLsolveExt.jl b/ext/TrixiNLsolveExt.jl
index fa188d04c71..2e6bca0ba7b 100644
--- a/ext/TrixiNLsolveExt.jl
+++ b/ext/TrixiNLsolveExt.jl
@@ -2,7 +2,7 @@
 module TrixiNLsolveExt
 
 # Required for finding coefficients in Butcher tableau in the third order of 
-# P-ERK scheme integrators
+# PERK scheme integrators
 if isdefined(Base, :get_extension)
     using NLsolve: nlsolve
 else
diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl
index dfda8ece687..5a1de036808 100644
--- a/src/auxiliary/precompile.jl
+++ b/src/auxiliary/precompile.jl
@@ -329,8 +329,10 @@ function _precompile_manual_()
     @assert Base.precompile(Tuple{typeof(Trixi.gauss_nodes_weights), Int})
     @assert Base.precompile(Tuple{typeof(Trixi.calc_forward_upper), Int})
     @assert Base.precompile(Tuple{typeof(Trixi.calc_forward_lower), Int})
-    @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_upper), Int, Val{:gauss}})
-    @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_lower), Int, Val{:gauss}})
+    @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_upper), Int,
+                                  Val{:gauss}})
+    @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_lower), Int,
+                                  Val{:gauss}})
     @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_upper), Int,
                                   Val{:gauss_lobatto}})
     @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_lower), Int,
diff --git a/src/auxiliary/special_elixirs.jl b/src/auxiliary/special_elixirs.jl
index d71a27aa96a..3e960945e7a 100644
--- a/src/auxiliary/special_elixirs.jl
+++ b/src/auxiliary/special_elixirs.jl
@@ -6,10 +6,11 @@
 #! format: noindent
 
 """
-    convergence_test([mod::Module=Main,] elixir::AbstractString, iterations; kwargs...)
+    convergence_test([mod::Module=Main,] elixir::AbstractString, iterations, RealT = Float64; kwargs...)
 
 Run `iterations` Trixi.jl simulations using the setup given in `elixir` and compute
 the experimental order of convergence (EOC) in the ``L^2`` and ``L^\\infty`` norm.
+Use `RealT` as the data type to represent the errors.
 In each iteration, the resolution of the respective mesh will be doubled.
 Additional keyword arguments `kwargs...` and the optional module `mod` are passed directly
 to [`trixi_include`](@ref).
@@ -18,12 +19,14 @@ This function assumes that the spatial resolution is set via the keywords
 `initial_refinement_level` (an integer) or `cells_per_dimension` (a tuple of
 integers, one per spatial dimension).
 """
-function convergence_test(mod::Module, elixir::AbstractString, iterations; kwargs...)
+function convergence_test(mod::Module, elixir::AbstractString, iterations,
+                          RealT = Float64;
+                          kwargs...)
     @assert(iterations>1,
             "Number of iterations must be bigger than 1 for a convergence analysis")
 
     # Types of errors to be calculated
-    errors = Dict(:l2 => Float64[], :linf => Float64[])
+    errors = Dict(:l2 => RealT[], :linf => RealT[])
 
     initial_resolution = extract_initial_resolution(elixir, kwargs)
 
@@ -105,7 +108,7 @@ function analyze_convergence(errors, iterations,
         println("")
 
         # Print mean EOCs
-        mean_values = zeros(nvariables)
+        mean_values = zeros(eltype(errors[:l2]), nvariables)
         for v in 1:nvariables
             mean_values[v] = sum(eocs[kind][:, v]) ./ length(eocs[kind][:, v])
             @printf("%-10s", "mean")
@@ -119,8 +122,9 @@ function analyze_convergence(errors, iterations,
     return eoc_mean_values
 end
 
-function convergence_test(elixir::AbstractString, iterations; kwargs...)
-    convergence_test(Main, elixir::AbstractString, iterations; kwargs...)
+function convergence_test(elixir::AbstractString, iterations, RealT = Float64;
+                          kwargs...)
+    convergence_test(Main, elixir::AbstractString, iterations, RealT; kwargs...)
 end
 
 # Helper methods used in the functions defined above
diff --git a/src/auxiliary/t8code.jl b/src/auxiliary/t8code.jl
index ddc0a2f1399..ffb81fd4689 100644
--- a/src/auxiliary/t8code.jl
+++ b/src/auxiliary/t8code.jl
@@ -32,10 +32,12 @@ function init_t8code()
             T8code.Libt8.p4est_init(C_NULL, LOG_LEVEL)
         end
 
+        # Clean up t8code before MPI shuts down.
         MPI.add_finalize_hook!() do
+            T8code.clean_up()
             status = T8code.Libt8.sc_finalize_noabort()
             if status != 0
-                @warn("Inconsistent state detected after finalizing t8code. Have you finalized all `T8codeMesh` objects and/or properly freed/un-referenced all t8code related objects?")
+                @warn("Inconsistent state detected after finalizing t8code.")
             end
         end
 
@@ -96,8 +98,8 @@ end
 # \return greater zero if the first entry in `elements` should be refined,
 #         smaller zero if the family `elements` shall be coarsened,
 #         zero else.
-function adapt_callback(forest,
-                        forest_from,
+function adapt_callback(forest::Ptr{t8_forest},
+                        forest_from::Ptr{t8_forest},
                         which_tree,
                         lelement_id,
                         ts,
@@ -189,7 +191,7 @@ function trixi_t8_adapt!(mesh, indicators)
 
     differences = trixi_t8_get_difference(old_levels, new_levels, 2^ndims(mesh))
 
-    mesh.forest = forest_cached
+    update_forest!(mesh, forest_cached)
 
     return differences
 end
diff --git a/src/callbacks_step/analysis_dg1d.jl b/src/callbacks_step/analysis_dg1d.jl
index d2613c325be..b4d49a39d5a 100644
--- a/src/callbacks_step/analysis_dg1d.jl
+++ b/src/callbacks_step/analysis_dg1d.jl
@@ -61,16 +61,17 @@ function calc_error_norms(func, u, t, analyzer,
                                        inv.(view(inverse_jacobian, :, element)))
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
-
         for i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i), t,
                                         equations)
             diff = func(u_exact, equations) -
                    func(get_node_vars(u_local, equations, dg, i), equations)
-            l2_error += diff .^ 2 * (weights[i] * jacobian_local[i])
+            # We take absolute value as we need the Jacobian here for the volume calculation
+            abs_jacobian_local_i = abs(jacobian_local[i])
+
+            l2_error += diff .^ 2 * (weights[i] * abs_jacobian_local_i)
             linf_error = @. max(linf_error, abs(diff))
-            total_volume += weights[i] * jacobian_local[i]
+            total_volume += weights[i] * abs_jacobian_local_i
         end
     end
 
diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl
index e089133fa17..e33533738ec 100644
--- a/src/callbacks_step/analysis_dg2d.jl
+++ b/src/callbacks_step/analysis_dg2d.jl
@@ -161,16 +161,17 @@ function calc_error_norms(func, u, t, analyzer,
                                        jacobian_tmp1)
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
-
         for j in eachnode(analyzer), i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j),
                                         t, equations)
             diff = func(u_exact, equations) -
                    func(get_node_vars(u_local, equations, dg, i, j), equations)
-            l2_error += diff .^ 2 * (weights[i] * weights[j] * jacobian_local[i, j])
+            # We take absolute value as we need the Jacobian here for the volume calculation
+            abs_jacobian_local_ij = abs(jacobian_local[i, j])
+
+            l2_error += diff .^ 2 * (weights[i] * weights[j] * abs_jacobian_local_ij)
             linf_error = @. max(linf_error, abs(diff))
-            total_volume += weights[i] * weights[j] * jacobian_local[i, j]
+            total_volume += weights[i] * weights[j] * abs_jacobian_local_ij
         end
     end
 
diff --git a/src/callbacks_step/analysis_dg2d_parallel.jl b/src/callbacks_step/analysis_dg2d_parallel.jl
index 5b3ae858ab7..10b6e30f95c 100644
--- a/src/callbacks_step/analysis_dg2d_parallel.jl
+++ b/src/callbacks_step/analysis_dg2d_parallel.jl
@@ -114,16 +114,17 @@ function calc_error_norms(func, u, t, analyzer,
                                        jacobian_tmp1)
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
-
         for j in eachnode(analyzer), i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j),
                                         t, equations)
             diff = func(u_exact, equations) -
                    func(get_node_vars(u_local, equations, dg, i, j), equations)
-            l2_error += diff .^ 2 * (weights[i] * weights[j] * jacobian_local[i, j])
+            # We take absolute value as we need the Jacobian here for the volume calculation
+            abs_jacobian_local_ij = abs(jacobian_local[i, j])
+
+            l2_error += diff .^ 2 * (weights[i] * weights[j] * abs_jacobian_local_ij)
             linf_error = @. max(linf_error, abs(diff))
-            volume += weights[i] * weights[j] * jacobian_local[i, j]
+            volume += weights[i] * weights[j] * abs_jacobian_local_ij
         end
     end
 
diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl
index 3c9c3a69f5e..efc1dc3bb3c 100644
--- a/src/callbacks_step/analysis_dg3d.jl
+++ b/src/callbacks_step/analysis_dg3d.jl
@@ -187,18 +187,19 @@ function calc_error_norms(func, u, t, analyzer,
                                        jacobian_tmp1, jacobian_tmp2)
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
-
         for k in eachnode(analyzer), j in eachnode(analyzer), i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j,
                                                         k), t, equations)
             diff = func(u_exact, equations) -
                    func(get_node_vars(u_local, equations, dg, i, j, k), equations)
+            # We take absolute value as we need the Jacobian here for the volume calculation
+            abs_jacobian_local_ijk = abs(jacobian_local[i, j, k])
+
             l2_error += diff .^ 2 *
-                        (weights[i] * weights[j] * weights[k] * jacobian_local[i, j, k])
+                        (weights[i] * weights[j] * weights[k] * abs_jacobian_local_ijk)
             linf_error = @. max(linf_error, abs(diff))
-            total_volume += weights[i] * weights[j] * weights[k] *
-                            jacobian_local[i, j, k]
+            total_volume += (weights[i] * weights[j] * weights[k] *
+                             abs_jacobian_local_ijk)
         end
     end
 
diff --git a/src/callbacks_step/analysis_dg3d_parallel.jl b/src/callbacks_step/analysis_dg3d_parallel.jl
index 285f0f62de6..ab9ba6a0055 100644
--- a/src/callbacks_step/analysis_dg3d_parallel.jl
+++ b/src/callbacks_step/analysis_dg3d_parallel.jl
@@ -31,17 +31,18 @@ function calc_error_norms(func, u, t, analyzer,
                                        jacobian_tmp1, jacobian_tmp2)
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
-
         for k in eachnode(analyzer), j in eachnode(analyzer), i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j,
                                                         k), t, equations)
             diff = func(u_exact, equations) -
                    func(get_node_vars(u_local, equations, dg, i, j, k), equations)
+            # We take absolute value as we need the Jacobian here for the volume calculation
+            abs_jacobian_local_ijk = abs(jacobian_local[i, j, k])
+
             l2_error += diff .^ 2 *
-                        (weights[i] * weights[j] * weights[k] * jacobian_local[i, j, k])
+                        (weights[i] * weights[j] * weights[k] * abs_jacobian_local_ijk)
             linf_error = @. max(linf_error, abs(diff))
-            volume += weights[i] * weights[j] * weights[k] * jacobian_local[i, j, k]
+            volume += weights[i] * weights[j] * weights[k] * abs_jacobian_local_ijk
         end
     end
 
diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl
index 024ce2d7e87..fe79a7e0590 100644
--- a/src/equations/compressible_navier_stokes_1d.jl
+++ b/src/equations/compressible_navier_stokes_1d.jl
@@ -148,7 +148,7 @@ end
 function flux(u, gradients, orientation::Integer,
               equations::CompressibleNavierStokesDiffusion1D)
     # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`.
-    rho, v1, _ = convert_transformed_to_primitive(u, equations)
+    _, v1, _ = convert_transformed_to_primitive(u, equations)
     # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T)
     # either computed directly or reverse engineered from the gradient of the entropy variables
     # by way of the `convert_gradient_variables` function.
@@ -257,6 +257,12 @@ end
     return T
 end
 
+@inline function velocity(u, equations::CompressibleNavierStokesDiffusion1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    return v1
+end
+
 @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                         <:Adiabatic})(flux_inner,
                                                                                       u_inner,
diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl
index 5708abb4e00..977dae18a43 100644
--- a/src/equations/compressible_navier_stokes_2d.jl
+++ b/src/equations/compressible_navier_stokes_2d.jl
@@ -148,7 +148,7 @@ end
 function flux(u, gradients, orientation::Integer,
               equations::CompressibleNavierStokesDiffusion2D)
     # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`.
-    rho, v1, v2, _ = convert_transformed_to_primitive(u, equations)
+    _, v1, v2, _ = convert_transformed_to_primitive(u, equations)
     # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T)
     # either computed directly or reverse engineered from the gradient of the entropy variables
     # by way of the `convert_gradient_variables` function.
@@ -281,6 +281,13 @@ end
     return T
 end
 
+@inline function velocity(u, equations::CompressibleNavierStokesDiffusion2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
 @inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion2D)
     # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v
 
diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl
index 9622882fc1a..551462e7595 100644
--- a/src/equations/compressible_navier_stokes_3d.jl
+++ b/src/equations/compressible_navier_stokes_3d.jl
@@ -148,7 +148,7 @@ end
 function flux(u, gradients, orientation::Integer,
               equations::CompressibleNavierStokesDiffusion3D)
     # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`.
-    rho, v1, v2, v3, _ = convert_transformed_to_primitive(u, equations)
+    _, v1, v2, v3, _ = convert_transformed_to_primitive(u, equations)
     # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, v3, T)
     # either computed directly or reverse engineered from the gradient of the entropy variables
     # by way of the `convert_gradient_variables` function.
@@ -309,6 +309,14 @@ end
     return T
 end
 
+@inline function velocity(u, equations::CompressibleNavierStokesDiffusion3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
 @inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion3D)
     # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v
 
diff --git a/src/equations/lattice_boltzmann_2d.jl b/src/equations/lattice_boltzmann_2d.jl
index bad417d9c84..920f862f436 100644
--- a/src/equations/lattice_boltzmann_2d.jl
+++ b/src/equations/lattice_boltzmann_2d.jl
@@ -58,7 +58,7 @@ The main sources for the base implementation were
    Raum**, Master Thesis, University of Cologne, 2018.
 3. Dieter Hänel, **Molekulare Gasdynamik**, Springer-Verlag Berlin Heidelberg, 2004
    [doi:10.1007/3-540-35047-0](https://doi.org/10.1007/3-540-35047-0)
-4. Dieter Krüger et al., **The Lattice Boltzmann Method**, Springer International Publishing, 2017
+4. Timm Krüger et al., **The Lattice Boltzmann Method**, Springer International Publishing, 2017
    [doi:10.1007/978-3-319-44649-3](https://doi.org/10.1007/978-3-319-44649-3)
 """
 struct LatticeBoltzmannEquations2D{RealT <: Real, CollisionOp} <:
diff --git a/src/equations/lattice_boltzmann_3d.jl b/src/equations/lattice_boltzmann_3d.jl
index bf4c365e0fe..1451eb7cbc1 100644
--- a/src/equations/lattice_boltzmann_3d.jl
+++ b/src/equations/lattice_boltzmann_3d.jl
@@ -97,7 +97,7 @@ The main sources for the base implementation were
    Raum**, Master Thesis, University of Cologne, 2018.
 3. Dieter Hänel, **Molekulare Gasdynamik**, Springer-Verlag Berlin Heidelberg, 2004
    [doi:10.1007/3-540-35047-0](https://doi.org/10.1007/3-540-35047-0)
-4. Dieter Krüger et al., **The Lattice Boltzmann Method**, Springer International Publishing, 2017
+4. Timm Krüger et al., **The Lattice Boltzmann Method**, Springer International Publishing, 2017
    [doi:10.1007/978-3-319-44649-3](https://doi.org/10.1007/978-3-319-44649-3)
 """
 struct LatticeBoltzmannEquations3D{RealT <: Real, CollisionOp} <:
diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl
index 0d16101257f..fa5b1d8f81a 100644
--- a/src/meshes/t8code_mesh.jl
+++ b/src/meshes/t8code_mesh.jl
@@ -7,7 +7,7 @@ to manage trees and mesh refinement.
 """
 mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <:
                AbstractMesh{NDIMS}
-    forest      :: Ptr{t8_forest} # cpointer to forest
+    forest      :: T8code.ForestWrapper
     is_parallel :: IsParallel
 
     # This specifies the geometry interpolation for each tree.
@@ -34,7 +34,7 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <:
                                RealT = Float64) where {NDIMS}
         is_parallel = mpi_isparallel() ? True() : False()
 
-        mesh = new{NDIMS, RealT, typeof(is_parallel), NDIMS + 2, length(nodes)}(forest,
+        mesh = new{NDIMS, RealT, typeof(is_parallel), NDIMS + 2, length(nodes)}(T8code.ForestWrapper(forest),
                                                                                 is_parallel)
 
         mesh.nodes = nodes
@@ -44,20 +44,37 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <:
         mesh.unsaved_changes = true
 
         finalizer(mesh) do mesh
-            # When finalizing, `forest`, `scheme`, `cmesh`, and `geometry` are
-            # also cleaned up from within `t8code`. The cleanup code for
-            # `cmesh` does some MPI calls for deallocating shared memory
-            # arrays. Due to garbage collection in Julia the order of shutdown
-            # is not deterministic. Hence, "manual" finalization might be
-            # necessary in order to avoid MPI-related error output when closing
-            # the Julia program/session.
-            t8_forest_unref(Ref(mesh.forest))
+            # In serial mode we can finalize the forest right away. In parallel
+            # mode we keep the forest till Trixi shuts down or the user
+            # finalizes the forest wrapper explicitly.
+            if !mpi_isparallel()
+                finalize(mesh.forest)
+            end
         end
 
         return mesh
     end
 end
 
+"""
+    Trixi.update_forest!(mesh::T8codeMesh, new_forest::Ptr{t8_forest})
+
+Update the `mesh` object with the `new_forest`. Ownership of the old forest
+goes to caller.
+
+# Arguments
+- `mesh::T8codeMesh`: Initialized mesh.
+- `new_forest::Ptr{t8_forest}`: New forest.
+
+Returns `nothing`.
+"""
+function update_forest!(mesh::T8codeMesh, new_forest::Ptr{t8_forest})
+    # `mesh.forest` must not be overwritten. Its lifetime is attached to the
+    # `mesh` lifetime. Thus, this setter function.
+    mesh.forest.pointer = new_forest
+    return nothing
+end
+
 const SerialT8codeMesh{NDIMS} = T8codeMesh{NDIMS, <:Real, <:False}
 const ParallelT8codeMesh{NDIMS} = T8codeMesh{NDIMS, <:Real, <:True}
 @inline mpi_parallel(mesh::SerialT8codeMesh) = False()
@@ -232,13 +249,13 @@ function T8codeMesh(ndims, ntrees, nelements, tree_node_coordinates, nodes,
     # The adapt callback refines the forest according to the `levels` array.
     # For each tree the callback recursively increases the refinement level
     # till it matches with the associated section in `levels.
-    forest = adapt(forest, adapt_callback; recursive = true, balance = false,
-                   partition = false, ghost = false, user_data = C_NULL)
+    forest = adapt_forest(forest, adapt_callback; recursive = true, balance = false,
+                          partition = false, ghost = false, user_data = C_NULL)
 
     @assert t8_forest_get_global_num_elements(forest) == nelements
 
     if mpi_isparallel()
-        forest = partition(forest)
+        forest = partition_forest(forest)
     end
 
     return T8codeMesh{ndims}(forest, tree_node_coordinates, nodes, boundary_names, "")
@@ -807,42 +824,27 @@ function adapt_callback_wrapper(forest,
 end
 
 """
-    Trixi.adapt!(mesh::T8codeMesh, adapt_callback; kwargs...)
+    Trixi.adapt_forest(forest::Ptr{t8_forest}, adapt_callback; kwargs...)
 
-Adapt a `T8codeMesh` according to a user-defined `adapt_callback`.
+Adapt a `T8codeMesh` according to a user-defined `adapt_callback`. This
+function is primarily for internal use. See `Trixi.adapt!(mesh::T8codeMesh,
+adapt_callback; kwargs...)` for a more detailed documentation.
 
 # Arguments
-- `mesh::T8codeMesh`: Initialized mesh object.
+- `forest::Ptr{t8_forest}`: New forest.
 - `adapt_callback`: A user-defined callback which tells the adaption routines
                     if an element should be refined, coarsened or stay unchanged.
+- `kwargs`: Refer to `Trixi.adapt!(mesh::T8codeMesh, adapt_callback; kwargs...)`.
 
-    The expected callback signature is as follows:
-
-      `adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, user_data)`
-        # Arguments
-        - `forest`: Pointer to the analyzed forest.
-        - `ltreeid`: Local index of the current tree where the analyzed elements are part of.
-        - `eclass_scheme`: Element class of `elements`.
-        - `lelemntid`: Local index of the first element in `elements`.
-        - `elements`: Array of elements. If consecutive elements form a family
-                      they are passed together, otherwise `elements` consists of just one element.
-        - `is_family`: Boolean signifying if `elements` represents a family or not.
-        - `user_data`: Void pointer to some arbitrary user data. Default value is `C_NULL`.
-        # Returns
-          -1 : Coarsen family of elements.
-           0 : Stay unchanged.
-           1 : Refine element.
+Note that the old forest usually gets deallocated within t8code. Call
+`t8_forest_ref(Ref(mesh.forest))` beforehand to prevent that.
 
-- `kwargs`:
-    - `recursive = true`: Adapt the forest recursively. If true the caller must ensure that the callback
-                          returns 0 for every analyzed element at some point to stop the recursion.
-    - `balance = true`: Make sure the adapted forest is 2^(NDIMS-1):1 balanced.
-    - `partition = true`: Partition the forest to redistribute elements evenly among MPI ranks.
-    - `ghost = true`: Create a ghost layer for MPI data exchange.
-    - `user_data = C_NULL`: Pointer to some arbitrary user-defined data.
+Returns a `Ptr{t8_forest}` to a new forest.
 """
-function adapt(forest::Ptr{t8_forest}, adapt_callback; recursive = true, balance = true,
-               partition = true, ghost = true, user_data = C_NULL)
+function adapt_forest(forest::Union{T8code.ForestWrapper, Ptr{t8_forest}}, adapt_callback;
+                      recursive = true,
+                      balance = true,
+                      partition = true, ghost = true, user_data = C_NULL)
     # Check that forest is a committed, that is valid and usable, forest.
     @assert t8_forest_is_committed(forest) != 0
 
@@ -885,33 +887,79 @@ function adapt(forest::Ptr{t8_forest}, adapt_callback; recursive = true, balance
     return new_forest
 end
 
+"""
+    Trixi.adapt!(mesh::T8codeMesh, adapt_callback; kwargs...)
+
+Adapt a `T8codeMesh` according to a user-defined `adapt_callback`.
+
+# Arguments
+- `mesh::T8codeMesh`: Initialized mesh object.
+- `adapt_callback`: A user-defined callback which tells the adaption routines
+                    if an element should be refined, coarsened or stay unchanged.
+
+    The expected callback signature is as follows:
+
+      `adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, user_data)`
+        # Arguments
+        - `forest`: Pointer to the analyzed forest.
+        - `ltreeid`: Local index of the current tree where the analyzed elements are part of.
+        - `eclass_scheme`: Element class of `elements`.
+        - `lelemntid`: Local index of the first element in `elements`.
+        - `elements`: Array of elements. If consecutive elements form a family
+                      they are passed together, otherwise `elements` consists of just one element.
+        - `is_family`: Boolean signifying if `elements` represents a family or not.
+        - `user_data`: Void pointer to some arbitrary user data. Default value is `C_NULL`.
+        # Returns
+          -1 : Coarsen family of elements.
+           0 : Stay unchanged.
+           1 : Refine element.
+
+- `kwargs`:
+    - `recursive = true`: Adapt the forest recursively. If true the caller must ensure that the callback
+                          returns 0 for every analyzed element at some point to stop the recursion.
+    - `balance = true`: Make sure the adapted forest is 2^(NDIMS-1):1 balanced.
+    - `partition = true`: Partition the forest to redistribute elements evenly among MPI ranks.
+    - `ghost = true`: Create a ghost layer for MPI data exchange.
+    - `user_data = C_NULL`: Pointer to some arbitrary user-defined data.
+
+Returns `nothing`.
+"""
 function adapt!(mesh::T8codeMesh, adapt_callback; kwargs...)
     # Call `t8_forest_ref(Ref(mesh.forest))` to keep it.
-    mesh.forest = adapt(mesh.forest, adapt_callback; kwargs...)
+    update_forest!(mesh, adapt_forest(mesh.forest, adapt_callback; kwargs...))
+    return nothing
 end
 
-"""
-    Trixi.balance!(mesh::T8codeMesh)
-
-Balance a `T8codeMesh` to ensure 2^(NDIMS-1):1 face neighbors.
-"""
-function balance!(mesh::T8codeMesh)
+function balance_forest(forest::Union{T8code.ForestWrapper, Ptr{t8_forest}})
     new_forest_ref = Ref{t8_forest_t}()
     t8_forest_init(new_forest_ref)
     new_forest = new_forest_ref[]
 
-    let set_from = mesh.forest, no_repartition = 1, do_ghost = 1
+    let set_from = forest, no_repartition = 1, do_ghost = 1
         t8_forest_set_balance(new_forest, set_from, no_repartition)
         t8_forest_set_ghost(new_forest, do_ghost, T8_GHOST_FACES)
         t8_forest_commit(new_forest)
     end
 
-    mesh.forest = new_forest
+    return new_forest
+end
 
+"""
+    Trixi.balance!(mesh::T8codeMesh)
+
+Balance a `T8codeMesh` to ensure 2^(NDIMS-1):1 face neighbors.
+
+# Arguments
+- `mesh::T8codeMesh`: Initialized mesh object.
+
+Returns `nothing`.
+"""
+function balance!(mesh::T8codeMesh)
+    update_forest!(mesh, balance_forest(mesh.forest))
     return nothing
 end
 
-function partition(forest::Ptr{t8_forest})
+function partition_forest(forest::Union{T8code.ForestWrapper, Ptr{t8_forest}})
     new_forest_ref = Ref{t8_forest_t}()
     t8_forest_init(new_forest_ref)
     new_forest = new_forest_ref[]
@@ -932,9 +980,11 @@ Partition a `T8codeMesh` in order to redistribute elements evenly among MPI rank
 
 # Arguments
 - `mesh::T8codeMesh`: Initialized mesh object.
+
+Returns `nothing`.
 """
 function partition!(mesh::T8codeMesh)
-    mesh.forest = partition(mesh.forest)
+    update_forest!(mesh, partition_forest(mesh.forest))
     return nothing
 end
 
@@ -951,7 +1001,7 @@ function count_interfaces(mesh::T8codeMesh)
     return count_interfaces(mesh.forest, ndims(mesh))
 end
 
-function count_interfaces(forest::Ptr{t8_forest}, ndims)
+function count_interfaces(forest, ndims)
     @assert t8_forest_is_committed(forest) != 0
 
     num_local_elements = t8_forest_get_local_num_elements(forest)
diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl
index 3ca429f63b6..294cc69f471 100644
--- a/src/semidiscretization/semidiscretization_euler_gravity.jl
+++ b/src/semidiscretization/semidiscretization_euler_gravity.jl
@@ -19,8 +19,8 @@ struct ParametersEulerGravity{RealT <: Real, TimestepGravity}
     background_density     :: RealT # aka rho0
     gravitational_constant :: RealT # aka G
     cfl                    :: RealT
-    resid_tol              :: RealT
-    n_iterations_max       :: Int
+    resid_tol              :: RealT # Hyp.-Diff. Eq. steady state tolerance
+    n_iterations_max       :: Int   # Max. number of iterations of the pseudo-time gravity solver
     timestep_gravity       :: TimestepGravity
 end
 
@@ -124,6 +124,7 @@ function SemidiscretizationEulerGravity(semi_euler::SemiEuler,
           SemidiscretizationHyperbolic{Mesh, <:AbstractHyperbolicDiffusionEquations}}
     u_ode = compute_coefficients(zero(real(semi_gravity)), semi_gravity)
     du_ode = similar(u_ode)
+    # Registers for gravity solver, tailored to the 2N and 3S* methods implemented below
     u_tmp1_ode = similar(u_ode)
     u_tmp2_ode = similar(u_ode)
     cache = (; u_ode, du_ode, u_tmp1_ode, u_tmp2_ode)
@@ -217,6 +218,9 @@ end
     calc_error_norms(func, u, t, analyzer, semi.semi_euler, cache_analysis)
 end
 
+# Coupled Euler and gravity solver at each Runge-Kutta stage, 
+# corresponding to Algorithm 2 in Schlottke-Lakemper et al. (2020),
+# https://dx.doi.org/10.1016/j.jcp.2021.110467
 function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t)
     @unpack semi_euler, semi_gravity, cache = semi
 
@@ -275,33 +279,33 @@ function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode)
     finalstep = false
     @unpack n_iterations_max, cfl, resid_tol, timestep_gravity = parameters
     iter = 0
-    t = zero(real(semi_gravity.solver))
+    tau = zero(real(semi_gravity.solver)) # Pseudo-time
 
     # iterate gravity solver until convergence or maximum number of iterations are reached
     @unpack equations = semi_gravity
     while !finalstep
-        dt = @trixi_timeit timer() "calculate dt" begin
-            cfl * max_dt(u_gravity, t, semi_gravity.mesh,
+        dtau = @trixi_timeit timer() "calculate dtau" begin
+            cfl * max_dt(u_gravity, tau, semi_gravity.mesh,
                    have_constant_speed(equations), equations,
                    semi_gravity.solver, semi_gravity.cache)
         end
 
         # evolve solution by one pseudo-time step
         time_start = time_ns()
-        timestep_gravity(cache, u_euler, t, dt, parameters, semi_gravity)
+        timestep_gravity(cache, u_euler, tau, dtau, parameters, semi_gravity)
         runtime = time_ns() - time_start
         put!(gravity_counter, runtime)
 
         # update iteration counter
         iter += 1
-        t += dt
+        tau += dtau
 
         # check if we reached the maximum number of iterations
         if n_iterations_max > 0 && iter >= n_iterations_max
             @warn "Max iterations reached: Gravity solver failed to converge!" residual=maximum(abs,
                                                                                                 @views du_gravity[1,
                                                                                                                   ..,
-                                                                                                                  :]) t=t dt=dt
+                                                                                                                  :]) tau=tau dtau=dtau
             finalstep = true
         end
 
@@ -315,34 +319,37 @@ function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode)
 end
 
 # Integrate gravity solver for 2N-type low-storage schemes
-function timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
+function timestep_gravity_2N!(cache, u_euler, tau, dtau, gravity_parameters,
+                              semi_gravity,
                               a, b, c)
     G = gravity_parameters.gravitational_constant
     rho0 = gravity_parameters.background_density
     grav_scale = -4.0 * pi * G
 
+    # Note that `u_ode` is `u_gravity` in `rhs!` above
     @unpack u_ode, du_ode, u_tmp1_ode = cache
     u_tmp1_ode .= zero(eltype(u_tmp1_ode))
     du_gravity = wrap_array(du_ode, semi_gravity)
     for stage in eachindex(c)
-        t_stage = t + dt * c[stage]
+        tau_stage = tau + dtau * c[stage]
 
         # rhs! has the source term for the harmonic problem
         # We don't need a `@trixi_timeit timer() "rhs!"` here since that's already
         # included in the `rhs!` call.
-        rhs!(du_ode, u_ode, semi_gravity, t_stage)
+        rhs!(du_ode, u_ode, semi_gravity, tau_stage)
 
         # Source term: Jeans instability OR coupling convergence test OR blast wave
         # put in gravity source term proportional to Euler density
         # OBS! subtract off the background density ρ_0 (spatial mean value)
+        # Note: Adding to `du_gravity` is essentially adding to `du_ode`!
         @views @. du_gravity[1, .., :] += grav_scale * (u_euler[1, .., :] - rho0)
 
         a_stage = a[stage]
-        b_stage_dt = b[stage] * dt
+        b_stage_dtau = b[stage] * dtau
         @trixi_timeit timer() "Runge-Kutta step" begin
             @threaded for idx in eachindex(u_ode)
                 u_tmp1_ode[idx] = du_ode[idx] - u_tmp1_ode[idx] * a_stage
-                u_ode[idx] += u_tmp1_ode[idx] * b_stage_dt
+                u_ode[idx] += u_tmp1_ode[idx] * b_stage_dtau
             end
         end
     end
@@ -350,7 +357,7 @@ function timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gr
     return nothing
 end
 
-function timestep_gravity_carpenter_kennedy_erk54_2N!(cache, u_euler, t, dt,
+function timestep_gravity_carpenter_kennedy_erk54_2N!(cache, u_euler, tau, dtau,
                                                       gravity_parameters, semi_gravity)
     # Coefficients for Carpenter's 5-stage 4th-order low-storage Runge-Kutta method
     a = SVector(0.0, 567301805773.0 / 1357537059087.0,
@@ -363,47 +370,50 @@ function timestep_gravity_carpenter_kennedy_erk54_2N!(cache, u_euler, t, dt,
                 2526269341429.0 / 6820363962896.0,
                 2006345519317.0 / 3224310063776.0, 2802321613138.0 / 2924317926251.0)
 
-    timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gravity, a, b,
-                         c)
+    timestep_gravity_2N!(cache, u_euler, tau, dtau, gravity_parameters, semi_gravity,
+                         a, b, c)
 end
 
 # Integrate gravity solver for 3S*-type low-storage schemes
-function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
+function timestep_gravity_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
                                   semi_gravity,
                                   gamma1, gamma2, gamma3, beta, delta, c)
     G = gravity_parameters.gravitational_constant
     rho0 = gravity_parameters.background_density
     grav_scale = -4 * G * pi
 
+    # Note that `u_ode` is `u_gravity` in `rhs!` above
     @unpack u_ode, du_ode, u_tmp1_ode, u_tmp2_ode = cache
     u_tmp1_ode .= zero(eltype(u_tmp1_ode))
     u_tmp2_ode .= u_ode
     du_gravity = wrap_array(du_ode, semi_gravity)
     for stage in eachindex(c)
-        t_stage = t + dt * c[stage]
+        tau_stage = tau + dtau * c[stage]
 
         # rhs! has the source term for the harmonic problem
         # We don't need a `@trixi_timeit timer() "rhs!"` here since that's already
         # included in the `rhs!` call.
-        rhs!(du_ode, u_ode, semi_gravity, t_stage)
+        rhs!(du_ode, u_ode, semi_gravity, tau_stage)
 
         # Source term: Jeans instability OR coupling convergence test OR blast wave
         # put in gravity source term proportional to Euler density
         # OBS! subtract off the background density ρ_0 around which the Jeans instability is perturbed
+        # Note: Adding to `du_gravity` is essentially adding to `du_ode`!
         @views @. du_gravity[1, .., :] += grav_scale * (u_euler[1, .., :] - rho0)
 
         delta_stage = delta[stage]
         gamma1_stage = gamma1[stage]
         gamma2_stage = gamma2[stage]
         gamma3_stage = gamma3[stage]
-        beta_stage_dt = beta[stage] * dt
+        beta_stage_dtau = beta[stage] * dtau
         @trixi_timeit timer() "Runge-Kutta step" begin
             @threaded for idx in eachindex(u_ode)
+                # See Algorithm 1 (3S* method) in Schlottke-Lakemper et al. (2020)
                 u_tmp1_ode[idx] += delta_stage * u_ode[idx]
                 u_ode[idx] = (gamma1_stage * u_ode[idx] +
                               gamma2_stage * u_tmp1_ode[idx] +
                               gamma3_stage * u_tmp2_ode[idx] +
-                              beta_stage_dt * du_ode[idx])
+                              beta_stage_dtau * du_ode[idx])
             end
         end
     end
@@ -411,7 +421,8 @@ function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
     return nothing
 end
 
-function timestep_gravity_erk51_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
+# First-order, 5-stage, 3S*-storage optimized method
+function timestep_gravity_erk51_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
                                         semi_gravity)
     # New 3Sstar coefficients optimized for polynomials of degree polydeg=3
     # and examples/parameters_hypdiff_lax_friedrichs.toml
@@ -434,11 +445,13 @@ function timestep_gravity_erk51_3Sstar!(cache, u_euler, t, dt, gravity_parameter
     c = SVector(0.0000000000000000E+00, 1.9189497208340553E-01, 1.9580448818599061E-01,
                 2.4241635859769023E-01, 5.0728347557552977E-01)
 
-    timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
+    timestep_gravity_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
+                             semi_gravity,
                              gamma1, gamma2, gamma3, beta, delta, c)
 end
 
-function timestep_gravity_erk52_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
+# Second-order, 5-stage, 3S*-storage optimized method
+function timestep_gravity_erk52_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
                                         semi_gravity)
     # New 3Sstar coefficients optimized for polynomials of degree polydeg=3
     # and examples/parameters_hypdiff_lax_friedrichs.toml
@@ -461,11 +474,13 @@ function timestep_gravity_erk52_3Sstar!(cache, u_euler, t, dt, gravity_parameter
     c = SVector(0.0000000000000000E+00, 4.5158640252832094E-01, 1.0221535725056414E+00,
                 1.4280257701954349E+00, 7.1581334196229851E-01)
 
-    timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
+    timestep_gravity_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
+                             semi_gravity,
                              gamma1, gamma2, gamma3, beta, delta, c)
 end
 
-function timestep_gravity_erk53_3Sstar!(cache, u_euler, t, dt, gravity_parameters,
+# Third-order, 5-stage, 3S*-storage optimized method
+function timestep_gravity_erk53_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
                                         semi_gravity)
     # New 3Sstar coefficients optimized for polynomials of degree polydeg=3
     # and examples/parameters_hypdiff_lax_friedrichs.toml
@@ -488,7 +503,8 @@ function timestep_gravity_erk53_3Sstar!(cache, u_euler, t, dt, gravity_parameter
     c = SVector(0.0000000000000000E+00, 8.4476964977404881E-02, 2.8110631488732202E-01,
                 5.7093842145029405E-01, 7.2999896418559662E-01)
 
-    timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, semi_gravity,
+    timestep_gravity_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters,
+                             semi_gravity,
                              gamma1, gamma2, gamma3, beta, delta, c)
 end
 
diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl
index cac1dba9c74..3ea668c0264 100644
--- a/src/solvers/dgsem/basis_lobatto_legendre.jl
+++ b/src/solvers/dgsem/basis_lobatto_legendre.jl
@@ -6,7 +6,7 @@
 #! format: noindent
 
 """
-    LobattoLegendreBasis([RealT=Float64,] polydeg::Integer)
+    LobattoLegendreBasis([RealT = Float64,] polydeg::Integer)
 
 Create a nodal Lobatto-Legendre basis for polynomials of degree `polydeg`.
 
@@ -37,15 +37,14 @@ end
 function LobattoLegendreBasis(RealT, polydeg::Integer)
     nnodes_ = polydeg + 1
 
-    # compute everything using `Float64` by default
-    nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_)
+    nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_, RealT)
     inverse_weights_ = inv.(weights_)
 
-    _, inverse_vandermonde_legendre_ = vandermonde_legendre(nodes_)
+    _, inverse_vandermonde_legendre_ = vandermonde_legendre(nodes_, RealT)
 
-    boundary_interpolation_ = zeros(nnodes_, 2)
-    boundary_interpolation_[:, 1] = calc_lhat(-1.0, nodes_, weights_)
-    boundary_interpolation_[:, 2] = calc_lhat(1.0, nodes_, weights_)
+    boundary_interpolation_ = zeros(RealT, nnodes_, 2)
+    boundary_interpolation_[:, 1] = calc_lhat(-one(RealT), nodes_, weights_)
+    boundary_interpolation_[:, 2] = calc_lhat(one(RealT), nodes_, weights_)
 
     derivative_matrix_ = polynomial_derivative_matrix(nodes_)
     derivative_split_ = calc_dsplit(nodes_, weights_)
@@ -79,7 +78,6 @@ function LobattoLegendreBasis(RealT, polydeg::Integer)
                                                            derivative_split_transpose,
                                                            derivative_dhat)
 end
-
 LobattoLegendreBasis(polydeg::Integer) = LobattoLegendreBasis(Float64, polydeg)
 
 function Base.show(io::IO, basis::LobattoLegendreBasis)
@@ -165,11 +163,10 @@ function MortarL2(basis::LobattoLegendreBasis)
     RealT = real(basis)
     nnodes_ = nnodes(basis)
 
-    # compute everything using `Float64` by default
-    forward_upper_ = calc_forward_upper(nnodes_)
-    forward_lower_ = calc_forward_lower(nnodes_)
-    reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss))
-    reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss))
+    forward_upper_ = calc_forward_upper(nnodes_, RealT)
+    forward_lower_ = calc_forward_lower(nnodes_, RealT)
+    reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss), RealT)
+    reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss), RealT)
 
     # type conversions to get the requested real type and enable possible
     # optimizations of runtime performance and latency
@@ -228,10 +225,10 @@ end
 # end
 
 # function MortarEC(basis::LobattoLegendreBasis{RealT}, surface_flux)
-#   forward_upper   = calc_forward_upper(n_nodes)
-#   forward_lower   = calc_forward_lower(n_nodes)
-#   l2reverse_upper = calc_reverse_upper(n_nodes, Val(:gauss_lobatto))
-#   l2reverse_lower = calc_reverse_lower(n_nodes, Val(:gauss_lobatto))
+#   forward_upper   = calc_forward_upper(n_nodes, RealT)
+#   forward_lower   = calc_forward_lower(n_nodes, RealT)
+#   l2reverse_upper = calc_reverse_upper(n_nodes, Val(:gauss_lobatto), RealT)
+#   l2reverse_lower = calc_reverse_lower(n_nodes, Val(:gauss_lobatto), RealT)
 
 #   # type conversions to make use of StaticArrays etc.
 #   nnodes_ = nnodes(basis)
@@ -263,8 +260,7 @@ function SolutionAnalyzer(basis::LobattoLegendreBasis;
     RealT = real(basis)
     nnodes_ = analysis_polydeg + 1
 
-    # compute everything using `Float64` by default
-    nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_)
+    nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_, RealT)
     vandermonde_ = polynomial_interpolation_matrix(get_nodes(basis), nodes_)
 
     # type conversions to get the requested real type and enable possible
@@ -322,11 +318,10 @@ end
 function AdaptorL2(basis::LobattoLegendreBasis{RealT}) where {RealT}
     nnodes_ = nnodes(basis)
 
-    # compute everything using `Float64` by default
-    forward_upper_ = calc_forward_upper(nnodes_)
-    forward_lower_ = calc_forward_lower(nnodes_)
-    reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss))
-    reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss))
+    forward_upper_ = calc_forward_upper(nnodes_, RealT)
+    forward_lower_ = calc_forward_lower(nnodes_, RealT)
+    reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss), RealT)
+    reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss), RealT)
 
     # type conversions to get the requested real type and enable possible
     # optimizations of runtime performance and latency
@@ -408,7 +403,7 @@ end
 # This implements algorithm 37 "PolynomialDerivativeMatrix" from Kopriva's book.
 function polynomial_derivative_matrix(nodes)
     n_nodes = length(nodes)
-    d = zeros(n_nodes, n_nodes)
+    d = zeros(eltype(nodes), n_nodes, n_nodes)
     wbary = barycentric_weights(nodes)
 
     for i in 1:n_nodes, j in 1:n_nodes
@@ -481,7 +476,7 @@ For details, see (especially Section 3)
 """
 function barycentric_weights(nodes)
     n_nodes = length(nodes)
-    weights = ones(n_nodes)
+    weights = ones(eltype(nodes), n_nodes)
 
     for j in 2:n_nodes, k in 1:(j - 1)
         weights[k] *= nodes[k] - nodes[j]
@@ -528,7 +523,7 @@ For details, see e.g. Section 2 of
 """
 function lagrange_interpolating_polynomials(x, nodes, wbary)
     n_nodes = length(nodes)
-    polynomials = zeros(n_nodes)
+    polynomials = zeros(eltype(nodes), n_nodes)
 
     for i in 1:n_nodes
         # Avoid division by zero when `x` is close to node by using 
@@ -553,7 +548,7 @@ function lagrange_interpolating_polynomials(x, nodes, wbary)
 end
 
 """
-    gauss_lobatto_nodes_weights(n_nodes::Integer)
+    gauss_lobatto_nodes_weights(n_nodes::Integer, RealT = Float64)
 
 Computes nodes ``x_j`` and weights ``w_j`` for the (Legendre-)Gauss-Lobatto quadrature.
 This implements algorithm 25 "GaussLobattoNodesAndWeights" from the book
@@ -563,15 +558,14 @@ This implements algorithm 25 "GaussLobattoNodesAndWeights" from the book
   Algorithms for scientists and engineers. 
   [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
 """
-# From FLUXO (but really from blue book by Kopriva)
-function gauss_lobatto_nodes_weights(n_nodes::Integer)
-    # From Kopriva's book
-    n_iterations = 10
-    tolerance = 1e-15
+function gauss_lobatto_nodes_weights(n_nodes::Integer, RealT = Float64)
+    # From FLUXO (but really from blue book by Kopriva)
+    n_iterations = 20
+    tolerance = 2 * eps(RealT) # Relative tolerance for Newton iteration
 
     # Initialize output
-    nodes = zeros(n_nodes)
-    weights = zeros(n_nodes)
+    nodes = zeros(RealT, n_nodes)
+    weights = zeros(RealT, n_nodes)
 
     # Special case for polynomial degree zero (first order finite volume)
     if n_nodes == 1
@@ -584,15 +578,15 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer)
     N = n_nodes - 1
 
     # Calculate values at boundary
-    nodes[1] = -1.0
-    nodes[end] = 1.0
-    weights[1] = 2 / (N * (N + 1))
+    nodes[1] = -1
+    nodes[end] = 1
+    weights[1] = RealT(2) / (N * (N + 1))
     weights[end] = weights[1]
 
     # Calculate interior values
     if N > 1
-        cont1 = pi / N
-        cont2 = 3 / (8 * N * pi)
+        cont1 = convert(RealT, pi) / N
+        cont2 = 3 / (8 * N * convert(RealT, pi))
 
         # Use symmetry -> only left side is computed
         for i in 1:(div(N + 1, 2) - 1)
@@ -608,6 +602,10 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer)
                 if abs(dx) < tolerance * abs(nodes[i + 1])
                     break
                 end
+
+                if k == n_iterations
+                    @warn "`gauss_lobatto_nodes_weights` Newton iteration did not converge"
+                end
             end
 
             # Calculate weight
@@ -622,8 +620,8 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer)
 
     # If odd number of nodes, set center node to origin (= 0.0) and calculate weight
     if n_nodes % 2 == 1
-        _, _, L = calc_q_and_l(N, 0)
-        nodes[div(N, 2) + 1] = 0.0
+        _, _, L = calc_q_and_l(N, zero(RealT))
+        nodes[div(N, 2) + 1] = 0
         weights[div(N, 2) + 1] = weights[1] / L^2
     end
 
@@ -631,11 +629,13 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer)
 end
 
 # From FLUXO (but really from blue book by Kopriva, algorithm 24)
-function calc_q_and_l(N::Integer, x::Float64)
-    L_Nm2 = 1.0
+function calc_q_and_l(N::Integer, x::Real)
+    RealT = typeof(x)
+
+    L_Nm2 = one(RealT)
     L_Nm1 = x
-    Lder_Nm2 = 0.0
-    Lder_Nm1 = 1.0
+    Lder_Nm2 = zero(RealT)
+    Lder_Nm1 = one(RealT)
 
     local L
     for i in 2:N
@@ -652,10 +652,9 @@ function calc_q_and_l(N::Integer, x::Float64)
 
     return q, qder, L
 end
-calc_q_and_l(N::Integer, x::Real) = calc_q_and_l(N, convert(Float64, x))
 
 """
-    gauss_nodes_weights(n_nodes::Integer)
+    gauss_nodes_weights(n_nodes::Integer, RealT = Float64)
 
 Computes nodes ``x_j`` and weights ``w_j`` for the Gauss-Legendre quadrature.
 This implements algorithm 23 "LegendreGaussNodesAndWeights" from the book
@@ -665,31 +664,30 @@ This implements algorithm 23 "LegendreGaussNodesAndWeights" from the book
   Algorithms for scientists and engineers. 
   [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
 """
-function gauss_nodes_weights(n_nodes::Integer)
-    # From Kopriva's book
-    n_iterations = 10
-    tolerance = 1e-15
+function gauss_nodes_weights(n_nodes::Integer, RealT = Float64)
+    n_iterations = 20
+    tolerance = 2 * eps(RealT) # Relative tolerance for Newton iteration
 
     # Initialize output
-    nodes = ones(n_nodes) * 1000
-    weights = zeros(n_nodes)
+    nodes = ones(RealT, n_nodes)
+    weights = zeros(RealT, n_nodes)
 
     # Get polynomial degree for convenience
     N = n_nodes - 1
     if N == 0
-        nodes .= 0.0
-        weights .= 2.0
+        nodes .= 0
+        weights .= 2
         return nodes, weights
     elseif N == 1
-        nodes[1] = -sqrt(1 / 3)
+        nodes[1] = -sqrt(one(RealT) / 3)
         nodes[end] = -nodes[1]
-        weights .= 1.0
+        weights .= 1
         return nodes, weights
     else # N > 1
         # Use symmetry property of the roots of the Legendre polynomials
         for i in 0:(div(N + 1, 2) - 1)
             # Starting guess for Newton method
-            nodes[i + 1] = -cos(pi / (2 * N + 2) * (2 * i + 1))
+            nodes[i + 1] = -cos(convert(RealT, pi) / (2 * N + 2) * (2 * i + 1))
 
             # Newton iteration to find root of Legendre polynomial (= integration node)
             for k in 0:n_iterations
@@ -699,6 +697,10 @@ function gauss_nodes_weights(n_nodes::Integer)
                 if abs(dx) < tolerance * abs(nodes[i + 1])
                     break
                 end
+
+                if k == n_iterations
+                    @warn "`gauss_nodes_weights` Newton iteration did not converge"
+                end
             end
 
             # Calculate weight
@@ -712,8 +714,8 @@ function gauss_nodes_weights(n_nodes::Integer)
 
         # If odd number of nodes, set center node to origin (= 0.0) and calculate weight
         if n_nodes % 2 == 1
-            poly, deriv = legendre_polynomial_and_derivative(N + 1, 0.0)
-            nodes[div(N, 2) + 1] = 0.0
+            poly, deriv = legendre_polynomial_and_derivative(N + 1, zero(RealT))
+            nodes[div(N, 2) + 1] = 0
             weights[div(N, 2) + 1] = (2 * N + 3) / deriv^2
         end
 
@@ -733,20 +735,21 @@ This implements algorithm 22 "LegendrePolynomialAndDerivative" from the book
   [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
 """
 function legendre_polynomial_and_derivative(N::Int, x::Real)
+    RealT = typeof(x)
     if N == 0
-        poly = 1.0
-        deriv = 0.0
+        poly = one(RealT)
+        deriv = zero(RealT)
     elseif N == 1
-        poly = convert(Float64, x)
-        deriv = 1.0
+        poly = x
+        deriv = one(RealT)
     else
-        poly_Nm2 = 1.0
-        poly_Nm1 = convert(Float64, x)
-        deriv_Nm2 = 0.0
-        deriv_Nm1 = 1.0
+        poly_Nm2 = one(RealT)
+        poly_Nm1 = copy(x)
+        deriv_Nm2 = zero(RealT)
+        deriv_Nm1 = one(RealT)
 
-        poly = 0.0
-        deriv = 0.0
+        poly = zero(RealT)
+        deriv = zero(RealT)
         for i in 2:N
             poly = ((2 * i - 1) * x * poly_Nm1 - (i - 1) * poly_Nm2) / i
             deriv = deriv_Nm2 + (2 * i - 1) * poly_Nm1
@@ -765,10 +768,10 @@ function legendre_polynomial_and_derivative(N::Int, x::Real)
 end
 
 # Calculate Legendre vandermonde matrix and its inverse
-function vandermonde_legendre(nodes, N)
+function vandermonde_legendre(nodes, N::Integer, RealT = Float64)
     n_nodes = length(nodes)
     n_modes = N + 1
-    vandermonde = zeros(n_nodes, n_modes)
+    vandermonde = zeros(RealT, n_nodes, n_modes)
 
     for i in 1:n_nodes
         for m in 1:n_modes
@@ -779,5 +782,7 @@ function vandermonde_legendre(nodes, N)
     inverse_vandermonde = inv(vandermonde)
     return vandermonde, inverse_vandermonde
 end
-vandermonde_legendre(nodes) = vandermonde_legendre(nodes, length(nodes) - 1)
+vandermonde_legendre(nodes, RealT = Float64) = vandermonde_legendre(nodes,
+                                                                    length(nodes) - 1,
+                                                                    RealT)
 end # @muladd
diff --git a/src/solvers/dgsem/l2projection.jl b/src/solvers/dgsem/l2projection.jl
index 0bb46f5ca15..436c1b92032 100644
--- a/src/solvers/dgsem/l2projection.jl
+++ b/src/solvers/dgsem/l2projection.jl
@@ -23,9 +23,9 @@
 # Calculate forward projection matrix for discrete L2 projection from large to upper
 #
 # Note: This is actually an interpolation.
-function calc_forward_upper(n_nodes)
+function calc_forward_upper(n_nodes, RealT = Float64)
     # Calculate nodes, weights, and barycentric weights
-    nodes, weights = gauss_lobatto_nodes_weights(n_nodes)
+    nodes, _ = gauss_lobatto_nodes_weights(n_nodes, RealT)
     wbary = barycentric_weights(nodes)
 
     # Calculate projection matrix (actually: interpolation)
@@ -43,9 +43,9 @@ end
 # Calculate forward projection matrix for discrete L2 projection from large to lower
 #
 # Note: This is actually an interpolation.
-function calc_forward_lower(n_nodes)
+function calc_forward_lower(n_nodes, RealT = Float64)
     # Calculate nodes, weights, and barycentric weights
-    nodes, weights = gauss_lobatto_nodes_weights(n_nodes)
+    nodes, _ = gauss_lobatto_nodes_weights(n_nodes, RealT)
     wbary = barycentric_weights(nodes)
 
     # Calculate projection matrix (actually: interpolation)
@@ -64,9 +64,9 @@ end
 #
 # Note: To not make the L2 projection exact, first convert to Gauss nodes,
 # perform projection, and convert back to Gauss-Lobatto.
-function calc_reverse_upper(n_nodes, ::Val{:gauss})
+function calc_reverse_upper(n_nodes, ::Val{:gauss}, RealT = Float64)
     # Calculate nodes, weights, and barycentric weights for Legendre-Gauss
-    gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes)
+    gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes, RealT)
     gauss_wbary = barycentric_weights(gauss_nodes)
 
     # Calculate projection matrix (actually: discrete L2 projection with errors)
@@ -80,7 +80,7 @@ function calc_reverse_upper(n_nodes, ::Val{:gauss})
     end
 
     # Calculate Vandermondes
-    lobatto_nodes, lobatto_weights = gauss_lobatto_nodes_weights(n_nodes)
+    lobatto_nodes, _ = gauss_lobatto_nodes_weights(n_nodes, RealT)
     gauss2lobatto = polynomial_interpolation_matrix(gauss_nodes, lobatto_nodes)
     lobatto2gauss = polynomial_interpolation_matrix(lobatto_nodes, gauss_nodes)
 
@@ -91,9 +91,9 @@ end
 #
 # Note: To not make the L2 projection exact, first convert to Gauss nodes,
 # perform projection, and convert back to Gauss-Lobatto.
-function calc_reverse_lower(n_nodes, ::Val{:gauss})
+function calc_reverse_lower(n_nodes, ::Val{:gauss}, RealT = Float64)
     # Calculate nodes, weights, and barycentric weights for Legendre-Gauss
-    gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes)
+    gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes, RealT)
     gauss_wbary = barycentric_weights(gauss_nodes)
 
     # Calculate projection matrix (actually: discrete L2 projection with errors)
@@ -107,7 +107,7 @@ function calc_reverse_lower(n_nodes, ::Val{:gauss})
     end
 
     # Calculate Vandermondes
-    lobatto_nodes, lobatto_weights = gauss_lobatto_nodes_weights(n_nodes)
+    lobatto_nodes, _ = gauss_lobatto_nodes_weights(n_nodes, RealT)
     gauss2lobatto = polynomial_interpolation_matrix(gauss_nodes, lobatto_nodes)
     lobatto2gauss = polynomial_interpolation_matrix(lobatto_nodes, gauss_nodes)
 
@@ -116,9 +116,9 @@ end
 
 # Calculate reverse projection matrix for discrete L2 projection from upper to large (Gauss-Lobatto
 # version)
-function calc_reverse_upper(n_nodes, ::Val{:gauss_lobatto})
+function calc_reverse_upper(n_nodes, ::Val{:gauss_lobatto}, RealT = Float64)
     # Calculate nodes, weights, and barycentric weights
-    nodes, weights = gauss_lobatto_nodes_weights(n_nodes)
+    nodes, weights = gauss_lobatto_nodes_weights(n_nodes, RealT)
     wbary = barycentric_weights(nodes)
 
     # Calculate projection matrix (actually: discrete L2 projection with errors)
@@ -135,9 +135,9 @@ end
 
 # Calculate reverse projection matrix for discrete L2 projection from lower to large (Gauss-Lobatto
 # version)
-function calc_reverse_lower(n_nodes, ::Val{:gauss_lobatto})
+function calc_reverse_lower(n_nodes, ::Val{:gauss_lobatto}, RealT = Float64)
     # Calculate nodes, weights, and barycentric weights
-    nodes, weights = gauss_lobatto_nodes_weights(n_nodes)
+    nodes, weights = gauss_lobatto_nodes_weights(n_nodes, RealT)
     wbary = barycentric_weights(nodes)
 
     # Calculate projection matrix (actually: discrete L2 projection with errors)
diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl
index 1b38eaa2d40..6a7bbbb71d7 100644
--- a/src/solvers/dgsem_p4est/dg_2d.jl
+++ b/src/solvers/dgsem_p4est/dg_2d.jl
@@ -384,7 +384,7 @@ end
 function prolong2mortars!(cache, u,
                           mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations,
                           mortar_l2::LobattoLegendreMortarL2,
-                          surface_integral, dg::DGSEM)
+                          dg::DGSEM)
     @unpack neighbor_ids, node_indices = cache.mortars
     index_range = eachnode(dg)
 
diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index baa9f7bf516..9f8306031c1 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -33,7 +33,7 @@ for `P4estMesh`es we call
 
     ```
     prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic,
-                                dg.mortar, dg.surface_integral, dg)
+                                dg.mortar, dg)
 
     calc_mortar_flux_divergence!(cache_parabolic.elements.surface_flux_values,
                                  mesh, equations_parabolic, dg.mortar,
@@ -42,11 +42,11 @@ for `P4estMesh`es we call
 instead of
     ```
     prolong2mortars!(cache, flux_viscous, mesh, equations_parabolic,
-                     dg.mortar, dg.surface_integral, dg)
+                     dg.mortar, dg)
 
     calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, mesh,
                       equations_parabolic,
-                      dg.mortar, dg.surface_integral, dg, cache)
+                      dg.mortar, dg, cache)
     ```
 =#
 function rhs_parabolic!(du, u, t, mesh::Union{P4estMesh{2}, P4estMesh{3}},
@@ -126,7 +126,7 @@ function rhs_parabolic!(du, u, t, mesh::Union{P4estMesh{2}, P4estMesh{3}},
     # !!! NOTE: we reuse the hyperbolic cache here since it contains "mortars" and "u_threaded". See https://github.com/trixi-framework/Trixi.jl/issues/1674 for a discussion
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic,
-                                    dg.mortar, dg.surface_integral, dg)
+                                    dg.mortar, dg)
     end
 
     # Calculate mortar fluxes (specialized for AbstractEquationsParabolic)
@@ -253,7 +253,7 @@ function calc_gradient!(gradients, u_transformed, t,
     # Prolong solution to mortars. This resues the hyperbolic version of `prolong2mortars`
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars!(cache, u_transformed, mesh, equations_parabolic,
-                         dg.mortar, dg.surface_integral, dg)
+                         dg.mortar, dg)
     end
 
     # Calculate mortar fluxes. This reuses the hyperbolic version of `calc_mortar_flux`,
@@ -676,7 +676,6 @@ function prolong2mortars_divergence!(cache, flux_viscous::Vector{Array{uEltype,
                                      mesh::Union{P4estMesh{2}, T8codeMesh{2}},
                                      equations,
                                      mortar_l2::LobattoLegendreMortarL2,
-                                     surface_integral,
                                      dg::DGSEM) where {uEltype <: Real}
     @unpack neighbor_ids, node_indices = cache.mortars
     @unpack contravariant_vectors = cache.elements
diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl
index 633c6f1d067..5d5225c6f9e 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl
@@ -137,7 +137,7 @@ function prolong2mpimortars!(cache, u,
                              mesh::Union{ParallelP4estMesh{2}, ParallelT8codeMesh{2}},
                              equations,
                              mortar_l2::LobattoLegendreMortarL2,
-                             surface_integral, dg::DGSEM)
+                             dg::DGSEM)
     @unpack node_indices = cache.mpi_mortars
     index_range = eachnode(dg)
 
diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl
index 0632e3d3d8e..ac87f339d7c 100644
--- a/src/solvers/dgsem_p4est/dg_3d.jl
+++ b/src/solvers/dgsem_p4est/dg_3d.jl
@@ -418,7 +418,7 @@ end
 function prolong2mortars!(cache, u,
                           mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations,
                           mortar_l2::LobattoLegendreMortarL2,
-                          surface_integral, dg::DGSEM)
+                          dg::DGSEM)
     @unpack fstar_tmp_threaded = cache
     @unpack neighbor_ids, node_indices = cache.mortars
     index_range = eachnode(dg)
diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl
index b1204e5a15c..f79d599364d 100644
--- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl
@@ -148,7 +148,7 @@ function calc_gradient!(gradients, u_transformed, t,
     # !!! should we have a separate mortars/u_threaded in cache_parabolic?
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars!(cache, u_transformed, mesh, equations_parabolic,
-                         dg.mortar, dg.surface_integral, dg)
+                         dg.mortar, dg)
     end
 
     # Calculate mortar fluxes. These should reuse the hyperbolic version of `calc_mortar_flux`,
@@ -647,7 +647,7 @@ function prolong2mortars_divergence!(cache, flux_viscous,
                                      mesh::Union{P4estMesh{3}, T8codeMesh{3}},
                                      equations,
                                      mortar_l2::LobattoLegendreMortarL2,
-                                     surface_integral, dg::DGSEM)
+                                     dg::DGSEM)
     @unpack neighbor_ids, node_indices = cache.mortars
     @unpack fstar_tmp_threaded = cache
     @unpack contravariant_vectors = cache.elements
diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel.jl b/src/solvers/dgsem_p4est/dg_3d_parallel.jl
index e51836e8a51..cd64e77e739 100644
--- a/src/solvers/dgsem_p4est/dg_3d_parallel.jl
+++ b/src/solvers/dgsem_p4est/dg_3d_parallel.jl
@@ -20,7 +20,7 @@ function rhs!(du, u, t,
     # Prolong solution to MPI mortars
     @trixi_timeit timer() "prolong2mpimortars" begin
         prolong2mpimortars!(cache, u, mesh, equations,
-                            dg.mortar, dg.surface_integral, dg)
+                            dg.mortar, dg)
     end
 
     # Start to send MPI data
@@ -64,7 +64,7 @@ function rhs!(du, u, t,
     # Prolong solution to mortars
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars!(cache, u, mesh, equations,
-                         dg.mortar, dg.surface_integral, dg)
+                         dg.mortar, dg)
     end
 
     # Calculate mortar fluxes
@@ -305,7 +305,7 @@ function prolong2mpimortars!(cache, u,
                              mesh::Union{ParallelP4estMesh{3}, ParallelT8codeMesh{3}},
                              equations,
                              mortar_l2::LobattoLegendreMortarL2,
-                             surface_integral, dg::DGSEM)
+                             dg::DGSEM)
     @unpack node_indices = cache.mpi_mortars
     index_range = eachnode(dg)
 
diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl
index 069d4f09bc4..98bd0889206 100644
--- a/src/solvers/dgsem_tree/dg_2d.jl
+++ b/src/solvers/dgsem_tree/dg_2d.jl
@@ -152,7 +152,7 @@ function rhs!(du, u, t,
     # Prolong solution to mortars
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars!(cache, u, mesh, equations,
-                         dg.mortar, dg.surface_integral, dg)
+                         dg.mortar, dg)
     end
 
     # Calculate mortar fluxes
@@ -798,7 +798,7 @@ end
 
 function prolong2mortars!(cache, u,
                           mesh::TreeMesh{2}, equations,
-                          mortar_l2::LobattoLegendreMortarL2, surface_integral,
+                          mortar_l2::LobattoLegendreMortarL2,
                           dg::DGSEM)
     @threaded for mortar in eachmortar(dg, cache)
         large_element = cache.mortars.neighbor_ids[3, mortar]
diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
index 22c29eb4725..2123eae5459 100644
--- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
@@ -89,7 +89,7 @@ function rhs_parabolic!(du, u, t, mesh::Union{TreeMesh{2}, TreeMesh{3}},
     # Prolong solution to mortars
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars!(cache, flux_viscous, mesh, equations_parabolic,
-                         dg.mortar, dg.surface_integral, dg)
+                         dg.mortar, dg)
     end
 
     # Calculate mortar fluxes
@@ -519,7 +519,7 @@ end
 function prolong2mortars!(cache, flux_viscous::Vector{Array{uEltype, 4}},
                           mesh::TreeMesh{2},
                           equations_parabolic::AbstractEquationsParabolic,
-                          mortar_l2::LobattoLegendreMortarL2, surface_integral,
+                          mortar_l2::LobattoLegendreMortarL2,
                           dg::DGSEM) where {uEltype <: Real}
     flux_viscous_x, flux_viscous_y = flux_viscous
     @threaded for mortar in eachmortar(dg, cache)
@@ -830,7 +830,7 @@ function calc_gradient!(gradients, u_transformed, t,
     # NOTE: This re-uses the implementation for hyperbolic terms in "dg_2d.jl"
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars!(cache, u_transformed, mesh, equations_parabolic,
-                         dg.mortar, dg.surface_integral, dg)
+                         dg.mortar, dg)
     end
 
     # Calculate mortar fluxes
diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl
index b21933f0a5f..b77308c32fe 100644
--- a/src/solvers/dgsem_tree/dg_2d_parallel.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl
@@ -463,7 +463,7 @@ function rhs!(du, u, t,
     # Prolong solution to MPI mortars
     @trixi_timeit timer() "prolong2mpimortars" begin
         prolong2mpimortars!(cache, u, mesh, equations,
-                            dg.mortar, dg.surface_integral, dg)
+                            dg.mortar, dg)
     end
 
     # Start to send MPI data
@@ -510,7 +510,7 @@ function rhs!(du, u, t,
     # Prolong solution to mortars
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars!(cache, u, mesh, equations,
-                         dg.mortar, dg.surface_integral, dg)
+                         dg.mortar, dg)
     end
 
     # Calculate mortar fluxes
@@ -597,7 +597,7 @@ end
 
 function prolong2mpimortars!(cache, u,
                              mesh::ParallelTreeMesh{2}, equations,
-                             mortar_l2::LobattoLegendreMortarL2, surface_integral,
+                             mortar_l2::LobattoLegendreMortarL2,
                              dg::DGSEM)
     @unpack mpi_mortars = cache
 
diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl
index 381534d2079..43f88bc9198 100644
--- a/src/solvers/dgsem_tree/dg_3d.jl
+++ b/src/solvers/dgsem_tree/dg_3d.jl
@@ -200,7 +200,7 @@ function rhs!(du, u, t,
     # Prolong solution to mortars
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars!(cache, u, mesh, equations,
-                         dg.mortar, dg.surface_integral, dg)
+                         dg.mortar, dg)
     end
 
     # Calculate mortar fluxes
@@ -864,7 +864,7 @@ end
 function prolong2mortars!(cache, u,
                           mesh::TreeMesh{3}, equations,
                           mortar_l2::LobattoLegendreMortarL2,
-                          surface_integral, dg::DGSEM)
+                          dg::DGSEM)
     # temporary buffer for projections
     @unpack fstar_tmp1_threaded = cache
 
diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
index 657f39bbffb..445a40e3a48 100644
--- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
@@ -496,7 +496,7 @@ function prolong2mortars!(cache,
                           mesh::TreeMesh{3},
                           equations_parabolic::AbstractEquationsParabolic,
                           mortar_l2::LobattoLegendreMortarL2,
-                          surface_integral, dg::DGSEM) where {uEltype <: Real}
+                          dg::DGSEM) where {uEltype <: Real}
     # temporary buffer for projections
     @unpack fstar_tmp1_threaded = cache
 
@@ -966,7 +966,7 @@ function calc_gradient!(gradients, u_transformed, t,
     # NOTE: This re-uses the implementation for hyperbolic terms in "dg_3d.jl"
     @trixi_timeit timer() "prolong2mortars" begin
         prolong2mortars!(cache, u_transformed, mesh, equations_parabolic,
-                         dg.mortar, dg.surface_integral, dg)
+                         dg.mortar, dg)
     end
 
     # Calculate mortar fluxes
diff --git a/src/solvers/fdsbp_tree/fdsbp.jl b/src/solvers/fdsbp_tree/fdsbp.jl
index 11b09c6df9c..5824c6595fa 100644
--- a/src/solvers/fdsbp_tree/fdsbp.jl
+++ b/src/solvers/fdsbp_tree/fdsbp.jl
@@ -51,7 +51,7 @@ create_cache(mesh, equations, mortar::Nothing, uEltype) = NamedTuple()
 nmortars(mortar::Nothing) = 0
 
 function prolong2mortars!(cache, u, mesh, equations, mortar::Nothing,
-                          surface_integral, dg::DG)
+                          dg::DG)
     @assert isempty(eachmortar(dg, cache))
 end
 
diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl
index 1143e4ecb93..fb2847e1c91 100644
--- a/src/time_integration/methods_2N.jl
+++ b/src/time_integration/methods_2N.jl
@@ -189,13 +189,15 @@ function step!(integrator::SimpleIntegrator2N)
     integrator.iter += 1
     integrator.t += integrator.dt
 
-    # handle callbacks
-    if callbacks isa CallbackSet
-        foreach(callbacks.discrete_callbacks) do cb
-            if cb.condition(integrator.u, integrator.t, integrator)
-                cb.affect!(integrator)
+    @trixi_timeit timer() "Step-Callbacks" begin
+        # handle callbacks
+        if callbacks isa CallbackSet
+            foreach(callbacks.discrete_callbacks) do cb
+                if cb.condition(integrator.u, integrator.t, integrator)
+                    cb.affect!(integrator)
+                end
+                return nothing
             end
-            return nothing
         end
     end
 
diff --git a/src/time_integration/methods_3Sstar.jl b/src/time_integration/methods_3Sstar.jl
index a3ab023b64b..f197ef60902 100644
--- a/src/time_integration/methods_3Sstar.jl
+++ b/src/time_integration/methods_3Sstar.jl
@@ -265,13 +265,15 @@ function step!(integrator::SimpleIntegrator3Sstar)
     integrator.iter += 1
     integrator.t += integrator.dt
 
-    # handle callbacks
-    if callbacks isa CallbackSet
-        foreach(callbacks.discrete_callbacks) do cb
-            if cb.condition(integrator.u, integrator.t, integrator)
-                cb.affect!(integrator)
+    @trixi_timeit timer() "Step-Callbacks" begin
+        # handle callbacks
+        if callbacks isa CallbackSet
+            foreach(callbacks.discrete_callbacks) do cb
+                if cb.condition(integrator.u, integrator.t, integrator)
+                    cb.affect!(integrator)
+                end
+                return nothing
             end
-            return nothing
         end
     end
 
diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl
index 5184ada1374..c6533aa2005 100644
--- a/src/time_integration/methods_SSP.jl
+++ b/src/time_integration/methods_SSP.jl
@@ -39,7 +39,7 @@ struct SimpleSSPRK33{StageCallbacks} <: SimpleAlgorithmSSP
         c = SVector(0.0, 1.0, 1 / 2)
 
         # Butcher tableau
-        #   c |       a
+        #   c |       A
         #   0 |
         #   1 |   1
         # 1/2 | 1/4  1/4
@@ -210,11 +210,14 @@ function solve!(integrator::SimpleIntegratorSSP)
         integrator.iter += 1
         integrator.t += integrator.dt
 
-        # handle callbacks
-        if callbacks isa CallbackSet
-            foreach(callbacks.discrete_callbacks) do cb
-                if cb.condition(integrator.u, integrator.t, integrator)
-                    cb.affect!(integrator)
+        @trixi_timeit timer() "Step-Callbacks" begin
+            # handle callbacks
+            if callbacks isa CallbackSet
+                foreach(callbacks.discrete_callbacks) do cb
+                    if cb.condition(integrator.u, integrator.t, integrator)
+                        cb.affect!(integrator)
+                    end
+                    return nothing
                 end
             end
         end
diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl
index 2451680a505..37d994dc9bf 100644
--- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl
+++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl
@@ -34,8 +34,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan,
     # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency)
     coeffs_max = num_stages - 2
 
-    a_matrix = zeros(coeffs_max, 2)
-    a_matrix[:, 1] = c[3:end]
+    a_matrix = zeros(2, coeffs_max)
+    a_matrix[1, :] = c[3:end]
 
     consistency_order = 2
 
@@ -56,8 +56,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan,
         num_monomial_coeffs = length(monomial_coeffs)
         @assert num_monomial_coeffs == coeffs_max
         A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs)
-        a_matrix[:, 1] -= A
-        a_matrix[:, 2] = A
+        a_matrix[1, :] -= A
+        a_matrix[2, :] = A
     end
 
     return a_matrix, c, dt_opt
@@ -78,8 +78,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages,
     # - 2 Since first entry of A is always zero (explicit method) and second is given by c_2 (consistency)
     coeffs_max = num_stages - 2
 
-    a_matrix = zeros(coeffs_max, 2)
-    a_matrix[:, 1] = c[3:end]
+    a_matrix = zeros(2, coeffs_max)
+    a_matrix[1, :] = c[3:end]
 
     path_monomial_coeffs = joinpath(base_path_monomial_coeffs,
                                     "gamma_" * string(num_stages) * ".txt")
@@ -91,8 +91,8 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages,
     @assert num_monomial_coeffs == coeffs_max
     A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs)
 
-    a_matrix[:, 1] -= A
-    a_matrix[:, 2] = A
+    a_matrix[1, :] -= A
+    a_matrix[2, :] = A
 
     return a_matrix, c
 end
@@ -131,16 +131,17 @@ optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solve
 
 Note: To use this integrator, the user must import the `Convex` and `ECOS` packages.
 """
-mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle
-    const num_stages::Int
+struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle
+    num_stages::Int
 
     a_matrix::Matrix{Float64}
     c::Vector{Float64}
     b1::Float64
     bS::Float64
     cS::Float64
+
     dt_opt::Union{Float64, Nothing}
-end # struct PairedExplicitRK2
+end
 
 # Constructor that reads the coefficients from a file
 function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString,
@@ -200,9 +201,8 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F,
     finalstep::Bool # added for convenience
     dtchangeable::Bool
     force_stepfail::Bool
-    # PairedExplicitRK2 stages:
+    # Additional PERK register
     k1::uType
-    k_higher::uType
 end
 
 function init(ode::ODEProblem, alg::PairedExplicitRK2;
@@ -211,9 +211,7 @@ function init(ode::ODEProblem, alg::PairedExplicitRK2;
     du = zero(u0)
     u_tmp = zero(u0)
 
-    # PairedExplicitRK2 stages
-    k1 = zero(u0)
-    k_higher = zero(u0)
+    k1 = zero(u0) # Additional PERK register
 
     t0 = first(ode.tspan)
     tdir = sign(ode.tspan[end] - ode.tspan[1])
@@ -226,7 +224,7 @@ function init(ode::ODEProblem, alg::PairedExplicitRK2;
                                                                      ode.tspan;
                                                                      kwargs...),
                                              false, true, false,
-                                             k1, k_higher)
+                                             k1)
 
     # initialize callbacks
     if callback isa CallbackSet
@@ -262,48 +260,21 @@ function step!(integrator::PairedExplicitRK2Integrator)
     end
 
     @trixi_timeit timer() "Paired Explicit Runge-Kutta ODE integration step" begin
-        # k1
-        integrator.f(integrator.du, integrator.u, prob.p, integrator.t)
-        @threaded for i in eachindex(integrator.du)
-            integrator.k1[i] = integrator.du[i] * integrator.dt
-        end
-
-        # Construct current state
-        @threaded for i in eachindex(integrator.u)
-            integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i]
-        end
-        # k2
-        integrator.f(integrator.du, integrator.u_tmp, prob.p,
-                     integrator.t + alg.c[2] * integrator.dt)
-
-        @threaded for i in eachindex(integrator.du)
-            integrator.k_higher[i] = integrator.du[i] * integrator.dt
-        end
+        # First and second stage are identical across all single/standalone PERK methods
+        PERK_k1!(integrator, prob.p)
+        PERK_k2!(integrator, prob.p, alg)
 
         # Higher stages
         for stage in 3:(alg.num_stages)
-            # Construct current state
-            @threaded for i in eachindex(integrator.u)
-                integrator.u_tmp[i] = integrator.u[i] +
-                                      alg.a_matrix[stage - 2, 1] *
-                                      integrator.k1[i] +
-                                      alg.a_matrix[stage - 2, 2] *
-                                      integrator.k_higher[i]
-            end
-
-            integrator.f(integrator.du, integrator.u_tmp, prob.p,
-                         integrator.t + alg.c[stage] * integrator.dt)
-
-            @threaded for i in eachindex(integrator.du)
-                integrator.k_higher[i] = integrator.du[i] * integrator.dt
-            end
+            PERK_ki!(integrator, prob.p, alg, stage)
         end
 
         @threaded for i in eachindex(integrator.u)
-            integrator.u[i] += alg.b1 * integrator.k1[i] +
-                               alg.bS * integrator.k_higher[i]
+            integrator.u[i] += integrator.dt *
+                               (alg.b1 * integrator.k1[i] +
+                                alg.bS * integrator.du[i])
         end
-    end # PairedExplicitRK2 step
+    end
 
     integrator.iter += 1
     integrator.t += integrator.dt
diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl
index 02bc3eeba36..5be1d498d59 100644
--- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl
+++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl
@@ -59,11 +59,11 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan,
                                                     monomial_coeffs, c;
                                                     verbose)
     end
-    # Fill A-matrix in P-ERK style
-    a_matrix = zeros(num_stages - 2, 2)
-    a_matrix[:, 1] = c[3:end]
-    a_matrix[:, 1] -= a_unknown
-    a_matrix[:, 2] = a_unknown
+    # Fill A-matrix in PERK style
+    a_matrix = zeros(2, num_stages - 2)
+    a_matrix[1, :] = c[3:end]
+    a_matrix[1, :] -= a_unknown
+    a_matrix[2, :] = a_unknown
 
     return a_matrix, c, dt_opt
 end
@@ -80,8 +80,8 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages,
     # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency)
     a_coeffs_max = num_stages - 2
 
-    a_matrix = zeros(a_coeffs_max, 2)
-    a_matrix[:, 1] = c[3:end]
+    a_matrix = zeros(2, a_coeffs_max)
+    a_matrix[1, :] = c[3:end]
 
     path_a_coeffs = joinpath(base_path_a_coeffs,
                              "a_" * string(num_stages) * ".txt")
@@ -91,9 +91,9 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages,
     num_a_coeffs = size(a_coeffs, 1)
 
     @assert num_a_coeffs == a_coeffs_max
-    # Fill A-matrix in P-ERK style
-    a_matrix[:, 1] -= a_coeffs
-    a_matrix[:, 2] = a_coeffs
+    # Fill A-matrix in PERK style
+    a_matrix[1, :] -= a_coeffs
+    a_matrix[2, :] = a_coeffs
 
     return a_matrix, c
 end
@@ -107,7 +107,7 @@ end
                       verbose = false, cS2 = 1.0f0)
 
     Parameters:
-    - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (P-ERK) method.
+    - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (PERK) method.
     - `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in 
       the Butcher tableau of the Runge Kutta method.
       The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages).txt")` and separated by line breaks.
@@ -122,7 +122,7 @@ end
       s is the number of stages, default is 1.0f0.
 
 The following structures and methods provide an implementation of
-the third-order paired explicit Runge-Kutta (P-ERK) method
+the third-order paired explicit Runge-Kutta (PERK) method
 optimized for a certain simulation setup (PDE, IC & BC, Riemann Solver, DG Solver).
 The original paper is
 - Nasab, Vermeire (2022)
@@ -135,13 +135,14 @@ Multirate Time-Integration based on Dynamic ODE Partitioning through Adaptively
 
 Note: To use this integrator, the user must import the `Convex`, `ECOS`, and `NLsolve` packages.
 """
-mutable struct PairedExplicitRK3 <: AbstractPairedExplicitRKSingle
-    const num_stages::Int # S
+struct PairedExplicitRK3 <: AbstractPairedExplicitRKSingle
+    num_stages::Int # S
 
     a_matrix::Matrix{Float64}
     c::Vector{Float64}
+
     dt_opt::Union{Float64, Nothing}
-end # struct PairedExplicitRK3
+end
 
 # Constructor for previously computed A Coeffs
 function PairedExplicitRK3(num_stages, base_path_a_coeffs::AbstractString,
@@ -195,9 +196,9 @@ mutable struct PairedExplicitRK3Integrator{RealT <: Real, uType, Params, Sol, F,
     finalstep::Bool # added for convenience
     dtchangeable::Bool
     force_stepfail::Bool
-    # PairedExplicitRK stages:
+    # Additional PERK3 registers
     k1::uType
-    k_higher::uType
+    kS1::uType
 end
 
 function init(ode::ODEProblem, alg::PairedExplicitRK3;
@@ -206,9 +207,9 @@ function init(ode::ODEProblem, alg::PairedExplicitRK3;
     du = zero(u0)
     u_tmp = zero(u0)
 
-    # PairedExplicitRK stages
+    # Additional PERK3 registers
     k1 = zero(u0)
-    k_higher = zero(u0)
+    kS1 = zero(u0)
 
     t0 = first(ode.tspan)
     tdir = sign(ode.tspan[end] - ode.tspan[1])
@@ -221,7 +222,7 @@ function init(ode::ODEProblem, alg::PairedExplicitRK3;
                                                                      ode.tspan;
                                                                      kwargs...),
                                              false, true, false,
-                                             k1, k_higher)
+                                             k1, kS1)
 
     # initialize callbacks
     if callback isa CallbackSet
@@ -258,72 +259,42 @@ function step!(integrator::PairedExplicitRK3Integrator)
     end
 
     @trixi_timeit timer() "Paired Explicit Runge-Kutta ODE integration step" begin
-        # k1
-        integrator.f(integrator.du, integrator.u, prob.p, integrator.t)
-        @threaded for i in eachindex(integrator.du)
-            integrator.k1[i] = integrator.du[i] * integrator.dt
-        end
-
-        # Construct current state
-        @threaded for i in eachindex(integrator.du)
-            integrator.u_tmp[i] = integrator.u[i] + alg.c[2] * integrator.k1[i]
-        end
-        # k2
-        integrator.f(integrator.du, integrator.u_tmp, prob.p,
-                     integrator.t + alg.c[2] * integrator.dt)
-
-        @threaded for i in eachindex(integrator.du)
-            integrator.k_higher[i] = integrator.du[i] * integrator.dt
-        end
+        # First and second stage are identical across all single/standalone PERK methods
+        PERK_k1!(integrator, prob.p)
+        PERK_k2!(integrator, prob.p, alg)
 
-        # Higher stages
         for stage in 3:(alg.num_stages - 1)
-            # Construct current state
-            @threaded for i in eachindex(integrator.du)
-                integrator.u_tmp[i] = integrator.u[i] +
-                                      alg.a_matrix[stage - 2, 1] *
-                                      integrator.k1[i] +
-                                      alg.a_matrix[stage - 2, 2] *
-                                      integrator.k_higher[i]
-            end
-
-            integrator.f(integrator.du, integrator.u_tmp, prob.p,
-                         integrator.t + alg.c[stage] * integrator.dt)
-
-            @threaded for i in eachindex(integrator.du)
-                integrator.k_higher[i] = integrator.du[i] * integrator.dt
-            end
+            PERK_ki!(integrator, prob.p, alg, stage)
         end
 
-        # Last stage
-        @threaded for i in eachindex(integrator.du)
-            integrator.u_tmp[i] = integrator.u[i] +
-                                  alg.a_matrix[alg.num_stages - 2, 1] *
-                                  integrator.k1[i] +
-                                  alg.a_matrix[alg.num_stages - 2, 2] *
-                                  integrator.k_higher[i]
+        # We need to store `du` of the S-1 stage in `kS1` for the final update:
+        @threaded for i in eachindex(integrator.u)
+            integrator.kS1[i] = integrator.du[i]
         end
 
-        integrator.f(integrator.du, integrator.u_tmp, prob.p,
-                     integrator.t + alg.c[alg.num_stages] * integrator.dt)
+        PERK_ki!(integrator, prob.p, alg, alg.num_stages)
 
         @threaded for i in eachindex(integrator.u)
             # "Own" PairedExplicitRK based on SSPRK33.
-            # Note that 'k_higher' carries the values of K_{S-1}
+            # Note that 'kS1' carries the values of K_{S-1}
             # and that we construct 'K_S' "in-place" from 'integrator.du'
-            integrator.u[i] += (integrator.k1[i] + integrator.k_higher[i] +
-                                4.0 * integrator.du[i] * integrator.dt) / 6.0
+            integrator.u[i] += integrator.dt *
+                               (integrator.k1[i] + integrator.kS1[i] +
+                                4.0 * integrator.du[i]) / 6.0
         end
-    end # PairedExplicitRK step timer
+    end
 
     integrator.iter += 1
     integrator.t += integrator.dt
 
-    # handle callbacks
-    if callbacks isa CallbackSet
-        for cb in callbacks.discrete_callbacks
-            if cb.condition(integrator.u, integrator.t, integrator)
-                cb.affect!(integrator)
+    @trixi_timeit timer() "Step-Callbacks" begin
+        # handle callbacks
+        if callbacks isa CallbackSet
+            foreach(callbacks.discrete_callbacks) do cb
+                if cb.condition(integrator.u, integrator.t, integrator)
+                    cb.affect!(integrator)
+                end
+                return nothing
             end
         end
     end
@@ -334,4 +305,13 @@ function step!(integrator::PairedExplicitRK3Integrator)
         terminate!(integrator)
     end
 end
+
+function Base.resize!(integrator::PairedExplicitRK3Integrator, new_size)
+    resize!(integrator.u, new_size)
+    resize!(integrator.du, new_size)
+    resize!(integrator.u_tmp, new_size)
+
+    resize!(integrator.k1, new_size)
+    resize!(integrator.kS1, new_size)
+end
 end # @muladd
diff --git a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl
index c606326738f..87f067fc39d 100644
--- a/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl
+++ b/src/time_integration/paired_explicit_runge_kutta/paired_explicit_runge_kutta.jl
@@ -108,13 +108,42 @@ function solve!(integrator::AbstractPairedExplicitRKIntegrator)
 
     @trixi_timeit timer() "main loop" while !integrator.finalstep
         step!(integrator)
-    end # "main loop" timer
+    end
 
     return TimeIntegratorSolution((first(prob.tspan), integrator.t),
                                   (prob.u0, integrator.u),
                                   integrator.sol.prob)
 end
 
+# Function that computes the first stage of a general PERK method
+@inline function PERK_k1!(integrator::AbstractPairedExplicitRKIntegrator, p)
+    integrator.f(integrator.k1, integrator.u, p, integrator.t)
+end
+
+@inline function PERK_k2!(integrator::AbstractPairedExplicitRKSingleIntegrator, p, alg)
+    @threaded for i in eachindex(integrator.du)
+        integrator.u_tmp[i] = integrator.u[i] +
+                              alg.c[2] * integrator.dt * integrator.k1[i]
+    end
+
+    integrator.f(integrator.du, integrator.u_tmp, p,
+                 integrator.t + alg.c[2] * integrator.dt)
+end
+
+@inline function PERK_ki!(integrator::AbstractPairedExplicitRKSingleIntegrator, p, alg,
+                          stage)
+    # Construct current state
+    @threaded for i in eachindex(integrator.u)
+        integrator.u_tmp[i] = integrator.u[i] +
+                              integrator.dt *
+                              (alg.a_matrix[1, stage - 2] * integrator.k1[i] +
+                               alg.a_matrix[2, stage - 2] * integrator.du[i])
+    end
+
+    integrator.f(integrator.du, integrator.u_tmp, p,
+                 integrator.t + alg.c[stage] * integrator.dt)
+end
+
 # used for AMR (Adaptive Mesh Refinement)
 function Base.resize!(integrator::AbstractPairedExplicitRKIntegrator, new_size)
     resize!(integrator.u, new_size)
@@ -122,7 +151,6 @@ function Base.resize!(integrator::AbstractPairedExplicitRKIntegrator, new_size)
     resize!(integrator.u_tmp, new_size)
 
     resize!(integrator.k1, new_size)
-    resize!(integrator.k_higher, new_size)
 end
 
 # get a cache where the RHS can be stored
diff --git a/test/Project.toml b/test/Project.toml
index 64c20ea4c37..f9a6241421b 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,6 +3,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
@@ -14,6 +15,7 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -23,6 +25,7 @@ Aqua = "0.8"
 CairoMakie = "0.10"
 Convex = "0.16"
 DelimitedFiles = "1"
+DoubleFloats = "1.4.0"
 Downloads = "1"
 ECOS = "1.1.2"
 ExplicitImports = "1.0.1"
@@ -34,6 +37,7 @@ NLsolve = "4.5.1"
 OrdinaryDiffEq = "6.49.1"
 Plots = "1.19"
 Printf = "1"
+Quadmath = "0.5.10"
 Random = "1"
 StableRNGs = "1.0.2"
 Test = "1"
diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl
index 7a1a3f3ac09..0e70282e77c 100644
--- a/test/test_p4est_2d.jl
+++ b/test/test_p4est_2d.jl
@@ -231,10 +231,10 @@ end
                             0.3507514f0
                         ],
                         linf=[
-                            0.2942562f0,
+                            0.2930063f0,
                             0.4079147f0,
                             0.3972956f0,
-                            1.0810697f0
+                            1.0764117f0
                         ],
                         tspan=(0.0f0, 1.0f0),
                         rtol=10 * sqrt(eps(Float32)), # to make CI pass
diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl
index 7de57a45178..c640be3c8b8 100644
--- a/test/test_parabolic_1d.jl
+++ b/test/test_parabolic_1d.jl
@@ -228,6 +228,29 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "TreeMesh1D: elixir_navierstokes_viscous_shock.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
+                                 "elixir_navierstokes_viscous_shock.jl"),
+                        l2=[
+                            0.00025762354103445303,
+                            0.0001433692781569829,
+                            0.00017369861968287976
+                        ],
+                        linf=[
+                            0.0016731940030498826,
+                            0.0010638575921477766,
+                            0.0011495207677434394
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 # Clean up afterwards: delete Trixi output directory
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index dc3e776154f..62142e1fc3b 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -760,6 +760,31 @@ end
     @test isapprox(drag_f, 1.5427441885921553, atol = 1e-13)
     @test isapprox(lift_f, 0.005621910087395724, atol = 1e-13)
 end
+
+@trixi_testset "P4estMesh2D: elixir_navierstokes_viscous_shock.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
+                                 "elixir_navierstokes_viscous_shock.jl"),
+                        l2=[
+                            0.0002576236264053728,
+                            0.00014336949098706463,
+                            7.189100338239794e-17,
+                            0.00017369905124642074
+                        ],
+                        linf=[
+                            0.0016731940983241156,
+                            0.0010638640749656147,
+                            5.59044079947959e-16,
+                            0.001149532023891009
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory
diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index 99ad05c7f70..9b7cb8a75ce 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -503,6 +503,33 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "P4estMesh3D: elixir_navierstokes_viscous_shock.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_3d_dgsem",
+                                 "elixir_navierstokes_viscous_shock.jl"),
+                        l2=[
+                            0.0002576235461250765,
+                            0.0001433693418567713,
+                            1.5583069105517042e-16,
+                            1.257551423107977e-16,
+                            0.00017369872990116004
+                        ],
+                        linf=[
+                            0.0016731930282756213,
+                            0.0010638586882356638,
+                            2.738015991633e-15,
+                            3.281831854493919e-15,
+                            0.0011495231318404686
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory
diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl
index 3061257d508..a5b561ed8fe 100644
--- a/test/test_structured_1d.jl
+++ b/test/test_structured_1d.jl
@@ -58,6 +58,20 @@ end
     end
 end
 
+@trixi_testset "elixir_advection_float128.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_float128.jl"),
+                        l2=Float128[6.49879312655540217059228636803492411e-09],
+                        linf=Float128[5.35548407857266390181158920649552284e-08])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
 # Testing the third-order paired explicit Runge-Kutta (PERK) method with its optimal CFL number
 @trixi_testset "elixir_burgers_perk3.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_perk3.jl"),
@@ -72,6 +86,15 @@ end
         du_ode = similar(u_ode)
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000
     end
+
+    # Test `resize!`
+    integrator = init(ode, ode_algorithm, dt = 42.0, callback = callbacks)
+    resize!(integrator, 42)
+    @test length(integrator.u) == 42
+    @test length(integrator.du) == 42
+    @test length(integrator.u_tmp) == 42
+    @test length(integrator.k1) == 42
+    @test length(integrator.kS1) == 42
 end
 
 # Testing the third-order paired explicit Runge-Kutta (PERK) method without stepsize callback
@@ -82,8 +105,8 @@ end
                         save_solution=SaveSolutionCallback(dt = 0.1 + 1.0e-8), # Adding a small epsilon to avoid floating-point precision issues
                         callbacks=CallbackSet(summary_callback, save_solution,
                                               analysis_callback, alive_callback),
-                        l2=[5.726144786001842e-7],
-                        linf=[3.430730019182704e-6])
+                        l2=[5.726144824784944e-7],
+                        linf=[3.43073006914274e-6])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl
index f061e2e1c30..0665a1be105 100644
--- a/test/test_tree_1d_advection.jl
+++ b/test/test_tree_1d_advection.jl
@@ -142,6 +142,20 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000
     end
 end
+
+@trixi_testset "elixir_advection_doublefloat.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_doublefloat.jl"),
+                        l2=Double64[6.80895929885700039832943251427357703e-11],
+                        linf=Double64[5.82834770064525291688100323411704252e-10])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 end # module
diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl
index b7152cb8b75..5d52f4a5d42 100644
--- a/test/test_tree_2d_mhd.jl
+++ b/test/test_tree_2d_mhd.jl
@@ -150,24 +150,24 @@ end
 
 @trixi_testset "elixir_mhd_ec_float32.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_float32.jl"),
-                        l2=Float32[0.036355644,
+                        l2=Float32[0.03635566,
+                                   0.042947732,
                                    0.042947736,
-                                   0.04294775,
-                                   0.025748005,
-                                   0.16211236,
-                                   0.017452478,
-                                   0.017452475,
-                                   0.026877576,
-                                   2.0951437f-7],
-                        linf=Float32[0.22100884,
-                                     0.28798944,
-                                     0.2879896,
-                                     0.20858234,
-                                     0.8812654,
-                                     0.09208202,
-                                     0.09208143,
-                                     0.14795381,
-                                     2.0912607f-6])
+                                   0.025748001,
+                                   0.16211228,
+                                   0.01745248,
+                                   0.017452491,
+                                   0.026877586,
+                                   2.417227f-7],
+                        linf=Float32[0.2210092,
+                                     0.28798974,
+                                     0.28799006,
+                                     0.20858109,
+                                     0.8812673,
+                                     0.09208107,
+                                     0.09208131,
+                                     0.14795369,
+                                     2.2078211f-6])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
diff --git a/test/test_trixi.jl b/test/test_trixi.jl
index be3521ed16a..024f3654ea0 100644
--- a/test/test_trixi.jl
+++ b/test/test_trixi.jl
@@ -15,6 +15,7 @@ are compared approximately against these reference values, using `atol, rtol`
 as absolute/relative tolerance.
 """
 macro test_trixi_include(elixir, args...)
+    # Note: The variables below are just Symbols, not actual errors/types
     local l2 = get_kwarg(args, :l2, nothing)
     local linf = get_kwarg(args, :linf, nothing)
     local RealT = get_kwarg(args, :RealT, :Float64)
@@ -24,6 +25,12 @@ macro test_trixi_include(elixir, args...)
     elseif RealT === :Float32
         atol_default = 500 * eps(Float32)
         rtol_default = sqrt(eps(Float32))
+    elseif RealT === :Float128
+        atol_default = 500 * eps(Float128)
+        rtol_default = sqrt(eps(Float128))
+    elseif RealT === :Double64
+        atol_default = 500 * eps(Double64)
+        rtol_default = sqrt(eps(Double64))
     end
     local atol = get_kwarg(args, :atol, atol_default)
     local rtol = get_kwarg(args, :rtol, rtol_default)
@@ -167,7 +174,10 @@ macro test_nowarn_mod(expr, additional_ignore_content = String[])
                                      r"WARNING: Method definition .* in module .* at .* overwritten .*.\n",
                                      # Warnings from third party packages
                                      r"┌ Warning: Problem status ALMOST_INFEASIBLE; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n",
-                                     r"┌ Warning: Problem status ALMOST_OPTIMAL; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n"]
+                                     r"┌ Warning: Problem status ALMOST_OPTIMAL; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n",
+                                     # Warnings for higher-precision floating data types
+                                     r"┌ Warning: #= /home/runner/work/Trixi.jl/Trixi.jl/src/solvers/dgsem/interpolation.jl:118 =#:\n│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.\n│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.\n└ @ Trixi ~/.julia/packages/LoopVectorization/tIJUA/src/condense_loopset.jl:1166\n",
+                                     r"┌ Warning: #= /home/runner/work/Trixi.jl/Trixi.jl/src/solvers/dgsem/interpolation.jl:136 =#:\n│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.\n│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.\n└ @ Trixi ~/.julia/packages/LoopVectorization/tIJUA/src/condense_loopset.jl:1166\n"]
                 append!(ignore_content, $additional_ignore_content)
                 for pattern in ignore_content
                     stderr_content = replace(stderr_content, pattern => "")
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 641c7244efe..97e32726ca7 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -268,6 +268,12 @@ end
 @timed_testset "interpolation" begin
     @testset "nodes and weights" begin
         @test Trixi.gauss_nodes_weights(1) == ([0.0], [2.0])
+
+        @test Trixi.gauss_nodes_weights(2)[1] ≈ [-1 / sqrt(3), 1 / sqrt(3)]
+        @test Trixi.gauss_nodes_weights(2)[2] == [1.0, 1.0]
+
+        @test Trixi.gauss_nodes_weights(3)[1] ≈ [-sqrt(3 / 5), 0.0, sqrt(3 / 5)]
+        @test Trixi.gauss_nodes_weights(3)[2] ≈ [5 / 9, 8 / 9, 5 / 9]
     end
 
     @testset "multiply_dimensionwise" begin
@@ -1700,7 +1706,7 @@ end
 
     ode_algorithm = Trixi.PairedExplicitRK2(6, path_coeff_file)
 
-    @test isapprox(ode_algorithm.a_matrix,
+    @test isapprox(transpose(ode_algorithm.a_matrix),
                    [0.12405417889682908 0.07594582110317093
                     0.16178873711001726 0.13821126288998273
                     0.16692313960864164 0.2330768603913584
@@ -1713,7 +1719,7 @@ end
     tspan = (0.0, 1.0)
     ode_algorithm = Trixi.PairedExplicitRK2(12, tspan, vec(eig_vals))
 
-    @test isapprox(ode_algorithm.a_matrix,
+    @test isapprox(transpose(ode_algorithm.a_matrix),
                    [0.06453812656711647 0.02637096434197444
                     0.09470601372274887 0.041657622640887494
                     0.12332877820069793 0.058489403617483886
@@ -1733,7 +1739,7 @@ end
 
     ode_algorithm = Trixi.PairedExplicitRK3(8, path_coeff_file)
 
-    @test isapprox(ode_algorithm.a_matrix,
+    @test isapprox(transpose(ode_algorithm.a_matrix),
                    [0.33551678438002486 0.06448322158043965
                     0.49653494442225443 0.10346507941960345
                     0.6496890912144586 0.15031092070647037
@@ -1748,7 +1754,7 @@ end
     tspan = (0.0, 1.0)
     ode_algorithm = Trixi.PairedExplicitRK3(13, tspan, vec(eig_vals))
 
-    @test isapprox(ode_algorithm.a_matrix,
+    @test isapprox(transpose(ode_algorithm.a_matrix),
                    [0.19121164778938382 0.008788355190848427
                     0.28723462747227385 0.012765384448655121
                     0.38017717196008227 0.019822834000382223
diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl
index c433338a238..fed71e49f0b 100644
--- a/test/test_unstructured_2d.jl
+++ b/test/test_unstructured_2d.jl
@@ -335,8 +335,8 @@ end
                         ],
                         linf=[
                             Float32(2.776782342826098),
-                            3.2162943f0, # this needs to be adapted
-                            3.6683278f0, # this needed to be adapted
+                            3.2116454f0, # this needs to be adapted
+                            3.6616623f0, # this needed to be adapted
                             Float32(2.052861364219655)
                         ],
                         tspan=(0.0f0, 0.25f0),