From 6f07dcf93f81cf89f060f065bcc187307774c38a Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 16:08:22 +0200
Subject: [PATCH 01/24] use VectorOfArray from RecursiveArrayTools.jl as
 datastructure for time integration

---
 Project.toml                                  |  2 +
 src/DispersiveShallowWater.jl                 |  1 +
 src/callbacks_step/analysis.jl                | 32 +++---
 src/callbacks_step/relaxation.jl              | 40 ++++----
 src/equations/bbm_bbm_1d.jl                   | 26 ++---
 .../bbm_bbm_variable_bathymetry_1d.jl         | 30 +++---
 src/equations/svaerd_kalisch_1d.jl            | 29 +++---
 src/semidiscretization.jl                     | 97 +++++++------------
 src/solver.jl                                 | 29 +++---
 src/visualization.jl                          | 14 +--
 10 files changed, 121 insertions(+), 179 deletions(-)

diff --git a/Project.toml b/Project.toml
index 3c918980..d5427cf8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,6 +11,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PolynomialBases = "c74db56a-226d-5e98-8bb0-a6049094aeea"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -29,6 +30,7 @@ LinearAlgebra = "1"
 PolynomialBases = "0.4.15"
 Printf = "1"
 RecipesBase = "1.1"
+RecursiveArrayTools = "3"
 Reexport = "1.0"
 Roots = "2.0.17"
 SciMLBase = "1.93, 2"
diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl
index a3835cb0..568a51ad 100644
--- a/src/DispersiveShallowWater.jl
+++ b/src/DispersiveShallowWater.jl
@@ -7,6 +7,7 @@ using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, lu, cholesky, ch
 using PolynomialBases: PolynomialBases
 using Printf: @printf, @sprintf
 using RecipesBase: RecipesBase, @recipe, @series
+using RecursiveArrayTools: VectorOfArray
 using Reexport: @reexport
 using Roots: AlefeldPotraShi, find_zero
 
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index b7f72282..ddf55a6b 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -149,10 +149,10 @@ function integrals(cb::DiscreteCallback{Condition,
     return (; zip(names, integrals)...)
 end
 
-function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, t,
+function initialize!(cb::DiscreteCallback{Condition, Affect!}, q, t,
                      integrator) where {Condition, Affect! <: AnalysisCallback}
     semi = integrator.p
-    initial_state_integrals = integrate(u_ode, semi)
+    initial_state_integrals = integrate(q, semi)
 
     analysis_callback = cb.affect!
     analysis_callback.initial_state_integrals = initial_state_integrals
@@ -184,8 +184,8 @@ end
 
 # This method is just called internally from `(analysis_callback::AnalysisCallback)(integrator)`
 # and serves as a function barrier. Additionally, it makes the code easier to profile and optimize.
-function (analysis_callback::AnalysisCallback)(io, u_ode, integrator, semi)
-    _, equations, solver = mesh_equations_solver(semi)
+function (analysis_callback::AnalysisCallback)(io, q, integrator, semi)
+    equations = semi.equations
     @unpack analysis_errors, analysis_integrals, tstops, errors, integrals = analysis_callback
     @unpack t, dt = integrator
     push!(tstops, t)
@@ -247,7 +247,7 @@ function (analysis_callback::AnalysisCallback)(io, u_ode, integrator, semi)
         println(io)
 
         # Calculate L2/Linf errors, which are also returned
-        l2_error, linf_error = calc_error_norms(u_ode, t, semi)
+        l2_error, linf_error = calc_error_norms(q, t, semi)
         current_errors = zeros(real(semi), (length(analysis_errors), nvariables(equations)))
         current_errors[1, :] = l2_error
         current_errors[2, :] = linf_error
@@ -266,7 +266,7 @@ function (analysis_callback::AnalysisCallback)(io, u_ode, integrator, semi)
         # Conservation error
         if :conservation_error in analysis_errors
             @unpack initial_state_integrals = analysis_callback
-            state_integrals = integrate(u_ode, semi)
+            state_integrals = integrate(q, semi)
             current_errors[3, :] = abs.(state_integrals - initial_state_integrals)
             print(io, " |∫q - ∫q₀|:  ")
             for v in eachvariable(equations)
@@ -282,7 +282,7 @@ function (analysis_callback::AnalysisCallback)(io, u_ode, integrator, semi)
             println(io, " Integrals:    ")
         end
         current_integrals = zeros(real(semi), length(analysis_integrals))
-        analyze_integrals!(io, current_integrals, 1, analysis_integrals, u_ode, t, semi)
+        analyze_integrals!(io, current_integrals, 1, analysis_integrals, q, t, semi)
         push!(integrals, current_integrals)
 
         println(io, "─"^100)
@@ -292,26 +292,25 @@ end
 
 # Iterate over tuples of analysis integrals in a type-stable way using "lispy tuple programming".
 function analyze_integrals!(io, current_integrals, i, analysis_integrals::NTuple{N, Any},
-                            u_ode,
-                            t, semi) where {N}
+                            q, t, semi) where {N}
 
     # Extract the first analysis integral and process it; keep the remaining to be processed later
     quantity = first(analysis_integrals)
     remaining_quantities = Base.tail(analysis_integrals)
 
-    res = analyze(quantity, u_ode, t, semi)
+    res = analyze(quantity, q, t, semi)
     current_integrals[i] = res
     @printf(io, " %-12s:", pretty_form_utf(quantity))
     @printf(io, "  % 10.8e", res)
     println(io)
 
     # Recursively call this method with the unprocessed integrals
-    analyze_integrals!(io, current_integrals, i + 1, remaining_quantities, u_ode, t, semi)
+    analyze_integrals!(io, current_integrals, i + 1, remaining_quantities, q, t, semi)
     return nothing
 end
 
 # terminate the type-stable iteration over tuples
-function analyze_integrals!(io, current_integrals, i, analysis_integrals::Tuple{}, u_ode, t,
+function analyze_integrals!(io, current_integrals, i, analysis_integrals::Tuple{}, q, t,
                             semi)
     nothing
 end
@@ -320,7 +319,6 @@ end
 function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition,
                                                                 Affect! <:
                                                                 AnalysisCallback}
-    analysis_callback = cb.affect!
     semi = sol.prob.p
 
     l2_error, linf_error = calc_error_norms(sol.u[end], sol.t[end], semi)
@@ -328,14 +326,14 @@ function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition,
     return (; l2 = l2_error, linf = linf_error)
 end
 
-function analyze(quantity, u_ode, t, semi::Semidiscretization)
-    integrate_quantity(u_ode -> quantity(u_ode, semi.equations), u_ode, semi)
+function analyze(quantity, q, t, semi::Semidiscretization)
+    integrate_quantity(q -> quantity(q, semi.equations), q, semi)
 end
 
 # modified entropy from Svärd-Kalisch equations need to take the whole vector `u` for every point in space
 function analyze(quantity::Union{typeof(energy_total_modified), typeof(entropy_modified)},
-                 u_ode, t, semi::Semidiscretization)
-    integrate_quantity(quantity, u_ode, semi)
+                 q, t, semi::Semidiscretization)
+    integrate_quantity(quantity, q, semi)
 end
 
 pretty_form_utf(::typeof(waterheight_total)) = "∫η"
diff --git a/src/callbacks_step/relaxation.jl b/src/callbacks_step/relaxation.jl
index 097b77a7..2562e62c 100644
--- a/src/callbacks_step/relaxation.jl
+++ b/src/callbacks_step/relaxation.jl
@@ -59,45 +59,40 @@ end
 @inline function (relaxation_callback::RelaxationCallback)(integrator)
     semi = integrator.p
     told = integrator.tprev
-    uold_ode = integrator.uprev
-    uold = wrap_array(uold_ode, semi)
+    qold = integrator.uprev
     tnew = integrator.t
-    unew_ode = integrator.u
-    unew = wrap_array(unew_ode, semi)
+    qnew = integrator.u
 
-    gamma = one(tnew)
     terminate_integration = false
-    gamma_lo = one(gamma) / 2
-    gamma_hi = 3 * one(gamma) / 2
+    gamma_lo = one(tnew) / 2
+    gamma_hi = 3 * one(tnew) / 2
 
-    function relaxation_functional(u, semi)
+    function relaxation_functional(q, semi)
         @unpack tmp1 = semi.cache
-        # modified entropy from Svärd-Kalisch equations need to take the whole vector `u` for every point in space
+        # modified entropy from Svärd-Kalisch equations need to take the whole vector `q` for every point in space
         if relaxation_callback.invariant isa
            Union{typeof(energy_total_modified), typeof(entropy_modified)}
-            return integrate_quantity!(tmp1, relaxation_callback.invariant, u, semi;
-                                       wrap = false)
+            return integrate_quantity!(tmp1, relaxation_callback.invariant, q, semi)
         else
             return integrate_quantity!(tmp1,
-                                       u_ode -> relaxation_callback.invariant(u_ode,
-                                                                              semi.equations),
-                                       u, semi; wrap = false)
+                                       q -> relaxation_callback.invariant(q, semi.equations),
+                                       q, semi)
         end
     end
 
-    function convex_combination(gamma, uold, unew)
-        @. uold + gamma * (unew - uold)
+    function convex_combination(gamma, old, new)
+        @. old + gamma * (new - old)
     end
-    energy_old = relaxation_functional(uold, semi)
+    energy_old = relaxation_functional(qold, semi)
 
     @trixi_timeit timer() "relaxation" begin
-        if (relaxation_functional(convex_combination(gamma_lo, uold, unew), semi) -
+        if (relaxation_functional(convex_combination(gamma_lo, qold, qnew), semi) -
             energy_old) *
-           (relaxation_functional(convex_combination(gamma_hi, uold, unew), semi) -
+           (relaxation_functional(convex_combination(gamma_hi, qold, qnew), semi) -
             energy_old) > 0
             terminate_integration = true
         else
-            gamma = find_zero(g -> relaxation_functional(convex_combination(g, uold, unew),
+            gamma = find_zero(g -> relaxation_functional(convex_combination(g, qold, qnew),
                                                          semi) -
                                    energy_old, (gamma_lo, gamma_hi), AlefeldPotraShi())
         end
@@ -106,9 +101,8 @@ end
             terminate_integration = true
         end
 
-        unew .= convex_combination(gamma, uold, unew)
-        unew_ode = reshape(unew, prod(size(unew)))
-        DiffEqBase.set_u!(integrator, unew_ode)
+        qnew .= convex_combination(gamma, qold, qnew)
+        DiffEqBase.set_u!(integrator, qnew)
         if !isapprox(tnew, first(integrator.opts.tstops))
             tgamma = convex_combination(gamma, told, tnew)
             DiffEqBase.set_t!(integrator, tgamma)
diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl
index 0f0fcda2..639c500e 100644
--- a/src/equations/bbm_bbm_1d.jl
+++ b/src/equations/bbm_bbm_1d.jl
@@ -173,17 +173,14 @@ end
 # - Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
 #   A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
 #   [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
-function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_condition,
+function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
               ::BoundaryConditionPeriodic, source_terms, solver, cache)
     @unpack invImD2, tmp1, tmp2 = cache
 
-    q = wrap_array(u_ode, mesh, equations, solver)
-    dq = wrap_array(du_ode, mesh, equations, solver)
-
-    eta = view(q, 1, :)
-    v = view(q, 2, :)
-    deta = view(dq, 1, :)
-    dv = view(dq, 2, :)
+    eta = q.u[1]
+    v = q.u[2]
+    deta = dq.u[1]
+    dv = dq.u[2]
 
     D = equations.D
     # energy and mass conservative semidiscretization
@@ -226,17 +223,14 @@ end
 # - Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
 #   A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
 #   [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
-function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_condition,
+function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
               ::BoundaryConditionReflecting, source_terms, solver, cache)
     @unpack invImD2d, invImD2n, tmp1, tmp2, tmp3 = cache
 
-    q = wrap_array(u_ode, mesh, equations, solver)
-    dq = wrap_array(du_ode, mesh, equations, solver)
-
-    eta = view(q, 1, :)
-    v = view(q, 2, :)
-    deta = view(dq, 1, :)
-    dv = view(dq, 2, :)
+    eta = q.u[1]
+    v = q.u[2]
+    deta = dq.u[1]
+    dv = dq.u[2]
 
     D = equations.D
     # energy and mass conservative semidiscretization
diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
index b677d2fc..0b1dea1a 100644
--- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl
+++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
@@ -257,19 +257,16 @@ end
 # - Joshua Lampert and Hendrik Ranocha (2024)
 #   Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations
 #   [DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669)
-function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
+function rhs!(dq, q, t, mesh, equations::BBMBBMVariableEquations1D,
               initial_condition, ::BoundaryConditionPeriodic, source_terms,
               solver, cache)
     @unpack invImDKD, invImD2K, D, tmp1, tmp2 = cache
 
-    q = wrap_array(u_ode, mesh, equations, solver)
-    dq = wrap_array(du_ode, mesh, equations, solver)
-
-    eta = view(q, 1, :)
-    v = view(q, 2, :)
-    deta = view(dq, 1, :)
-    dv = view(dq, 2, :)
-    dD = view(dq, 3, :)
+    eta = q.u[1]
+    v = q.u[2]
+    deta = dq.u[1]
+    dv = dq.u[2]
+    dD = dq.u[3]
     fill!(dD, zero(eltype(dD)))
 
     if solver.D1 isa PeriodicDerivativeOperator ||
@@ -305,19 +302,16 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
     return nothing
 end
 
-function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
+function rhs!(dq, q, t, mesh, equations::BBMBBMVariableEquations1D,
               initial_condition, ::BoundaryConditionReflecting, source_terms,
               solver, cache)
     @unpack invImDKDn, invImD2Kd, D, tmp1, tmp2, tmp3 = cache
 
-    q = wrap_array(u_ode, mesh, equations, solver)
-    dq = wrap_array(du_ode, mesh, equations, solver)
-
-    eta = view(q, 1, :)
-    v = view(q, 2, :)
-    deta = view(dq, 1, :)
-    dv = view(dq, 2, :)
-    dD = view(dq, 3, :)
+    eta = q.u[1]
+    v = q.u[2]
+    deta = dq.u[1]
+    dv = dq.u[2]
+    dD = dq.u[3]
     fill!(dD, zero(eltype(dD)))
 
     # energy and mass conservative semidiscretization
diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
index 4d9699c2..84a78ad0 100644
--- a/src/equations/svaerd_kalisch_1d.jl
+++ b/src/equations/svaerd_kalisch_1d.jl
@@ -199,18 +199,16 @@ end
 # - Joshua Lampert and Hendrik Ranocha (2024)
 #   Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations
 #   [DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669)
-function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D,
+function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D,
               initial_condition, ::BoundaryConditionPeriodic, source_terms,
               solver, cache)
     @unpack factorization, D1betaD1, D, h, hv, alpha_hat, gamma_hat, tmp1, tmp2, D1_central, M = cache
-    q = wrap_array(u_ode, mesh, equations, solver)
-    dq = wrap_array(du_ode, mesh, equations, solver)
-
-    eta = view(q, 1, :)
-    v = view(q, 2, :)
-    deta = view(dq, 1, :)
-    dv = view(dq, 2, :)
-    dD = view(dq, 3, :)
+
+    eta = q.u[1]
+    v = q.u[2]
+    deta = dq.u[1]
+    dv = dq.u[2]
+    dD = dq.u[3]
     fill!(dD, zero(eltype(dD)))
 
     @trixi_timeit timer() "deta hyperbolic" begin
@@ -332,11 +330,10 @@ number of nodes as length of the second dimension.
 `cache` needs to hold the first-derivative SBP operator `D1`.
 """
 @inline function energy_total_modified(q, equations::SvaerdKalischEquations1D, cache)
-    e_modified = zeros(eltype(q), size(q, 2))
+    e_modified = zeros(eltype(q), size(q, 1))
     # Need to compute new beta_hat, do not use the old one from the `cache`
-    eta = view(q, 1, :)
-    v = view(q, 2, :)
-    D = view(q, 3, :)
+    v = q.u[2]
+    D = q.u[3]
     beta_hat = equations.beta * D .^ 3
     if cache.D1 isa PeriodicDerivativeOperator ||
        cache.D1 isa UniformPeriodicCoupledOperator
@@ -344,10 +341,10 @@ number of nodes as length of the second dimension.
     elseif cache.D1 isa PeriodicUpwindOperators
         tmp = 0.5 * beta_hat .* ((cache.D1.minus * v) .^ 2)
     else
-        @error "unknown type of first-derivative operator"
+        @error "unknown type of first-derivative operator: $(typeof(cache.D1))"
     end
-    for i in 1:size(q, 2)
-        e_modified[i] = energy_total(view(q, :, i), equations) + tmp[i]
+    for i in 1:size(q, 1)
+        e_modified[i] = energy_total(view(q, i, :), equations) + tmp[i]
     end
     return e_modified
 end
diff --git a/src/semidiscretization.jl b/src/semidiscretization.jl
index 83aae58f..f532b0e6 100644
--- a/src/semidiscretization.jl
+++ b/src/semidiscretization.jl
@@ -117,71 +117,45 @@ Get the grid of a semidiscretization.
 """
 grid(semi::Semidiscretization) = grid(semi.solver)
 
-function PolynomialBases.integrate(func, u_ode, semi::Semidiscretization; wrap = true)
-    if wrap == true
-        q = wrap_array(u_ode, semi)
-        integrals = zeros(real(semi), nvariables(semi))
-        for v in eachvariable(semi.equations)
-            integrals[v] = integrate(func, q[v, :], semi.solver.D1)
-        end
-        return integrals
-    else
-        integrate(func, u_ode, semi.solver.D1)
+function PolynomialBases.integrate(func, q::VectorOfArray, semi::Semidiscretization)
+    integrals = zeros(real(semi), nvariables(semi))
+    for v in eachvariable(semi.equations)
+        integrals[v] = integrate(func, q.u[v], semi.solver.D1)
     end
+    return integrals
+end
+function PolynomialBases.integrate(func, quantity, semi::Semidiscretization)
+    integrate(func, quantity, semi.solver.D1)
 end
-function PolynomialBases.integrate(u_ode, semi::Semidiscretization; wrap = true)
-    integrate(identity, u_ode, semi; wrap = wrap)
+function PolynomialBases.integrate(q, semi::Semidiscretization)
+    integrate(identity, q, semi)
 end
 
-function integrate_quantity(func, u_ode, semi::Semidiscretization; wrap = true)
-    if wrap == true
-        q = wrap_array(u_ode, semi)
-    else
-        q = u_ode
-    end
-    quantity = zeros(eltype(q), size(q, 2))
-    for i in 1:size(q, 2)
-        quantity[i] = func(view(q, :, i))
-    end
-    integrate(quantity, semi; wrap = false)
+function integrate_quantity(func, q, semi::Semidiscretization)
+    quantity = zeros(eltype(q), nnodes(semi))
+    integrate_quantity!(quantity, func, q, semi)
 end
 
-function integrate_quantity!(quantity, func, u_ode, semi::Semidiscretization; wrap = true)
-    if wrap == true
-        q = wrap_array(u_ode, semi)
-    else
-        q = u_ode
-    end
-    for i in 1:size(q, 2)
-        quantity[i] = func(view(q, :, i))
+function integrate_quantity!(quantity, func, q, semi::Semidiscretization)
+    for i in eachnode(semi)
+        quantity[i] = func([q.u[v][i] for v in eachvariable(semi.equations)])
     end
-    integrate(quantity, semi; wrap = false)
+    integrate(quantity, semi)
 end
 
 # modified entropy from Svärd-Kalisch equations need to take the whole vector `q` for every point in space
 function integrate_quantity(func::Union{typeof(energy_total_modified),
-                                        typeof(entropy_modified)}, u_ode,
-                            semi::Semidiscretization; wrap = true)
-    if wrap == true
-        q = wrap_array(u_ode, semi)
-    else
-        q = u_ode
-    end
+                                        typeof(entropy_modified)}, q,
+                            semi::Semidiscretization)
     quantity = func(q, semi.equations, semi.cache)
-    integrate(quantity, semi; wrap = false)
+    integrate(quantity, semi)
 end
 
 function integrate_quantity!(quantity,
                              func::Union{typeof(energy_total_modified),
-                                         typeof(entropy_modified)}, u_ode,
-                             semi::Semidiscretization; wrap = true)
-    if wrap == true
-        q = wrap_array(u_ode, semi)
-    else
-        q = u_ode
-    end
-    quantity = func(q, semi.equations, semi.cache)
-    integrate(quantity, semi; wrap = false)
+                                         typeof(entropy_modified)}, q,
+                             semi::Semidiscretization)
+    integrate_quantity(func, q, semi)
 end
 
 @inline function mesh_equations_solver(semi::Semidiscretization)
@@ -194,18 +168,14 @@ end
     return mesh, equations, solver, cache
 end
 
-function calc_error_norms(u, t, semi::Semidiscretization)
-    calc_error_norms(u, t, semi.initial_condition, mesh_equations_solver(semi)...)
-end
-
-function wrap_array(u_ode, semi::Semidiscretization)
-    wrap_array(u_ode, mesh_equations_solver(semi)...)
+function calc_error_norms(q, t, semi::Semidiscretization)
+    calc_error_norms(q, t, semi.initial_condition, mesh_equations_solver(semi)...)
 end
 
-function rhs!(du_ode, u_ode, semi::Semidiscretization, t)
+function rhs!(dq, q, semi::Semidiscretization, t)
     @unpack mesh, equations, initial_condition, boundary_conditions, solver, source_terms, cache = semi
 
-    @trixi_timeit timer() "rhs!" rhs!(du_ode, u_ode, t, mesh, equations, initial_condition,
+    @trixi_timeit timer() "rhs!" rhs!(dq, q, t, mesh, equations, initial_condition,
                                       boundary_conditions, source_terms, solver, cache)
 
     return nothing
@@ -213,13 +183,12 @@ end
 
 function compute_coefficients(func, t, semi::Semidiscretization)
     @unpack mesh, equations, solver = semi
-    u_ode = allocate_coefficients(mesh_equations_solver(semi)...)
-    compute_coefficients!(u_ode, func, t, semi)
-    return u_ode
+    q = allocate_coefficients(mesh_equations_solver(semi)...)
+    compute_coefficients!(q, func, t, semi)
+    return q
 end
 
-function compute_coefficients!(u_ode, func, t, semi::Semidiscretization)
-    q = wrap_array(u_ode, semi)
+function compute_coefficients!(q, func, t, semi::Semidiscretization)
     # Call `compute_coefficients` defined by the solver
     compute_coefficients!(q, func, t, mesh_equations_solver(semi)...)
 end
@@ -231,7 +200,7 @@ Wrap the semidiscretization `semi` as an ODE problem in the time interval `tspan
 that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/).
 """
 function semidiscretize(semi::Semidiscretization, tspan)
-    u0_ode = compute_coefficients(semi.initial_condition, first(tspan), semi)
+    q0 = compute_coefficients(semi.initial_condition, first(tspan), semi)
     iip = true # is-inplace, i.e., we modify a vector when calling rhs!
-    return ODEProblem{iip}(rhs!, u0_ode, tspan, semi)
+    return ODEProblem{iip}(rhs!, q0, tspan, semi)
 end
diff --git a/src/solver.jl b/src/solver.jl
index 6070ac19..6206fede 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -76,37 +76,30 @@ end
 grid(solver::Solver) = grid(solver.D1)
 
 @inline eachnode(solver::Solver) = Base.OneTo(length(grid(solver)))
-@inline Base.real(solver::Solver{RealT}) where {RealT} = RealT
+@inline Base.real(::Solver{RealT}) where {RealT} = RealT
 
-@inline function set_node_vars!(u, u_node, equations, indices...)
+@inline function set_node_vars!(q, q_node, equations, indices...)
     for v in eachvariable(equations)
-        u[v, indices...] = u_node[v]
+        q.u[v][indices...] = q_node[v]
     end
     return nothing
 end
 
 function allocate_coefficients(mesh::Mesh1D, equations, solver::AbstractSolver)
-    # cf. wrap_array
-    zeros(real(solver), nvariables(equations) * nnodes(mesh)^ndims(mesh))
+    return VectorOfArray([zeros(real(solver), nnodes(mesh)) for _ in eachvariable(equations)])
 end
 
-function wrap_array(u_ode, mesh, equations, solver)
-    unsafe_wrap(Array{eltype(u_ode), ndims(mesh) + 1}, pointer(u_ode),
-                (nvariables(equations), ntuple(_ -> nnodes(mesh), ndims(mesh))...))
-end
-
-function compute_coefficients!(u, func, t, mesh::Mesh1D, equations, solver::AbstractSolver)
+function compute_coefficients!(q, func, t, mesh::Mesh1D, equations, solver::AbstractSolver)
     x = grid(solver)
     for i in eachnode(solver)
-        u_node = func(x[i], t, equations, mesh)
-        set_node_vars!(u, u_node, equations, i)
+        q_node = func(x[i], t, equations, mesh)
+        set_node_vars!(q, q_node, equations, i)
     end
 end
 
-function calc_error_norms(u_ode, t, initial_condition, mesh::Mesh1D, equations,
+function calc_error_norms(q, t, initial_condition, mesh::Mesh1D, equations,
                           solver::AbstractSolver)
     x = grid(solver)
-    q = wrap_array(u_ode, mesh, equations, solver)
     q_exact = zeros(real(solver), (nvariables(equations), nnodes(mesh)))
     for i in eachnode(solver)
         q_exact[:, i] = initial_condition(x[i], t, equations, mesh)
@@ -114,7 +107,7 @@ function calc_error_norms(u_ode, t, initial_condition, mesh::Mesh1D, equations,
     l2_error = zeros(real(solver), nvariables(equations))
     linf_error = similar(l2_error)
     for v in eachvariable(equations)
-        @views diff = q[v, :] - q_exact[v, :]
+        @views diff = q.u[v] - q_exact[v, :]
         l2_error[v] = integrate(q -> q^2, diff, solver.D1) |> sqrt
         linf_error[v] = maximum(abs.(diff))
     end
@@ -130,9 +123,9 @@ function calc_sources!(dq, q, t, source_terms,
                        equations::AbstractEquations{1}, solver::Solver)
     x = grid(solver)
     for i in eachnode(solver)
-        local_source = source_terms(view(q, :, i), x[i], t, equations)
+        local_source = source_terms(view(q, i, :), x[i], t, equations)
         for v in eachvariable(equations)
-            dq[v, i] += local_source[v]
+            dq.u[v][i] += local_source[v]
         end
     end
     return nothing
diff --git a/src/visualization.jl b/src/visualization.jl
index 1f29df24..cffd2003 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -33,22 +33,22 @@ end
     t = sol.t[step]
 
     if plot_initial == true
-        q_exact = wrap_array(compute_coefficients(initial_condition, t, semi), semi)
+        q_exact = compute_coefficients(initial_condition, t, semi)
     end
 
-    q = wrap_array(sol.u[step], semi)
+    q = sol.u[step]
     data = zeros(nvars, nnodes(semi))
     if plot_bathymetry == true
         bathy = zeros(nnodes(semi))
     end
     for j in eachnode(semi)
         if plot_bathymetry == true
-            bathy[j] = bathymetry(view(q, :, j), equations)
+            bathy[j] = bathymetry(view(q, j, :), equations)
         end
         if plot_initial == true
-            q_exact[:, j] .= conversion(view(q_exact, :, j), equations)
+            q_exact[j, :] .= conversion(q_exact[j, :], equations)
         end
-        data[:, j] .= conversion(view(q, :, j), equations)
+        data[:, j] .= conversion(view(q, j, :), equations)
     end
 
     plot_title --> "$(get_name(semi.equations)) at t = $(round(t, digits=5))"
@@ -63,7 +63,7 @@ end
                 subplot --> i
                 linestyle := :solid
                 label := "initial $(names[i])"
-                grid(semi), q_exact[i, :]
+                grid(semi), q_exact[:, i]
             end
         end
 
@@ -112,7 +112,7 @@ end
             # Allow that the spatial value `x` is not on the grid. Thus, interpolate the given values to the provided `x`
             # with a linear spline.
             solution[i, k] = linear_interpolation(grid(semi),
-                                                  view(wrap_array(sol.u[k], semi), i, :))(x)
+                                                  view(sol.u[k], :, i))(x)
         end
     end
 

From f8fa9d558fdac463fd2794bcff6b32cda954dc9a Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 16:47:54 +0200
Subject: [PATCH 02/24] remove temporary vectors

---
 src/equations/bbm_bbm_1d.jl                   | 29 +++++------------
 .../bbm_bbm_variable_bathymetry_1d.jl         | 31 +++++--------------
 2 files changed, 15 insertions(+), 45 deletions(-)

diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl
index 639c500e..284424b1 100644
--- a/src/equations/bbm_bbm_1d.jl
+++ b/src/equations/bbm_bbm_1d.jl
@@ -133,8 +133,7 @@ function create_cache(mesh, equations::BBMBBMEquations1D,
                       RealT, uEltype)
     D = equations.D
     invImD2 = lu(I - 1 / 6 * D^2 * sparse(solver.D2))
-    tmp2 = Array{RealT}(undef, nnodes(mesh))
-    return (invImD2 = invImD2, tmp2 = tmp2)
+    return (invImD2 = invImD2)
 end
 
 function create_cache(mesh, equations::BBMBBMEquations1D,
@@ -164,9 +163,7 @@ function create_cache(mesh, equations::BBMBBMEquations1D,
     else
         @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
     end
-    tmp2 = Array{RealT}(undef, nnodes(mesh))
-    tmp3 = Array{RealT}(undef, nnodes(mesh) - 2)
-    return (invImD2d = invImD2d, invImD2n = invImD2n, tmp2 = tmp2, tmp3 = tmp3)
+    return (invImD2d = invImD2d, invImD2n = invImD2n)
 end
 
 # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see
@@ -204,17 +201,11 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
     @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,
                                                        solver)
 
-    # To use the in-place version `ldiv!` instead of `\`, we need temporary arrays
-    # since `deta` and `dv` are not stored contiguously
     @trixi_timeit timer() "deta elliptic" begin
-        tmp1[:] = deta
-        ldiv!(tmp2, invImD2, tmp1)
-        deta[:] = tmp2
+        ldiv!(invImD2, deta)
     end
     @trixi_timeit timer() "dv elliptic" begin
-        tmp2[:] = dv
-        ldiv!(tmp1, invImD2, tmp2)
-        dv[:] = tmp1
+        ldiv!(invImD2, dv)
     end
     return nothing
 end
@@ -225,7 +216,7 @@ end
 #   [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
 function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
               ::BoundaryConditionReflecting, source_terms, solver, cache)
-    @unpack invImD2d, invImD2n, tmp1, tmp2, tmp3 = cache
+    @unpack invImD2d, invImD2n = cache
 
     eta = q.u[1]
     v = q.u[2]
@@ -251,18 +242,12 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
     @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,
                                                        solver)
 
-    # To use the in-place version `ldiv!` instead of `\`, we need temporary arrays
-    # since `deta` and `dv` are not stored contiguously
     @trixi_timeit timer() "deta elliptic" begin
-        tmp1[:] = deta
-        ldiv!(tmp2, invImD2n, tmp1)
-        deta[:] = tmp2
+        ldiv!(invImD2n, deta)
     end
     @trixi_timeit timer() "dv elliptic" begin
-        tmp2[:] = dv
-        ldiv!(tmp3, invImD2d, tmp2[2:(end - 1)])
+        ldiv!(invImD2d, dv[2:(end - 1)])
         dv[1] = dv[end] = zero(eltype(dv))
-        dv[2:(end - 1)] = tmp3
     end
 
     return nothing
diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
index 0b1dea1a..bb1e17c0 100644
--- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl
+++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
@@ -210,8 +210,7 @@ function create_cache(mesh, equations::BBMBBMVariableEquations1D,
         @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
     end
     invImD2K = lu(I - 1 / 6 * sparse(solver.D2) * K)
-    tmp2 = Array{RealT}(undef, nnodes(mesh))
-    return (invImDKD = invImDKD, invImD2K = invImD2K, D = D, tmp2 = tmp2)
+    return (invImDKD = invImDKD, invImD2K = invImD2K, D = D)
 end
 
 function create_cache(mesh, equations::BBMBBMVariableEquations1D,
@@ -248,9 +247,7 @@ function create_cache(mesh, equations::BBMBBMVariableEquations1D,
     else
         @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
     end
-    tmp2 = Array{RealT}(undef, nnodes(mesh))
-    tmp3 = Array{RealT}(undef, nnodes(mesh) - 2)
-    return (invImD2Kd = invImD2Kd, invImDKDn = invImDKDn, D = D, tmp2 = tmp2, tmp3 = tmp3)
+    return (invImD2Kd = invImD2Kd, invImDKDn = invImDKDn, D = D)
 end
 
 # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see
@@ -260,7 +257,7 @@ end
 function rhs!(dq, q, t, mesh, equations::BBMBBMVariableEquations1D,
               initial_condition, ::BoundaryConditionPeriodic, source_terms,
               solver, cache)
-    @unpack invImDKD, invImD2K, D, tmp1, tmp2 = cache
+    @unpack invImDKD, invImD2K, D = cache
 
     eta = q.u[1]
     v = q.u[2]
@@ -286,17 +283,11 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMVariableEquations1D,
     @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,
                                                        solver)
 
-    # To use the in-place version `ldiv!` instead of `\`, we need temporary arrays
-    # since `deta` and `dv` are not stored contiguously
     @trixi_timeit timer() "deta elliptic" begin
-        tmp1[:] = deta
-        ldiv!(tmp2, invImDKD, tmp1)
-        deta[:] = tmp2
+        ldiv!(invImDKD, deta)
     end
     @trixi_timeit timer() "dv elliptic" begin
-        tmp2[:] = dv
-        ldiv!(tmp1, invImD2K, tmp2)
-        dv[:] = tmp1
+        ldiv!(invImD2K, dv)
     end
 
     return nothing
@@ -305,7 +296,7 @@ end
 function rhs!(dq, q, t, mesh, equations::BBMBBMVariableEquations1D,
               initial_condition, ::BoundaryConditionReflecting, source_terms,
               solver, cache)
-    @unpack invImDKDn, invImD2Kd, D, tmp1, tmp2, tmp3 = cache
+    @unpack invImDKDn, invImD2Kd, D = cache
 
     eta = q.u[1]
     v = q.u[2]
@@ -332,18 +323,12 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMVariableEquations1D,
     @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,
                                                        solver)
 
-    # To use the in-place version `ldiv!` instead of `\`, we need temporary arrays
-    # since `deta` and `dv` are not stored contiguously
     @trixi_timeit timer() "deta elliptic" begin
-        tmp1[:] = deta
-        ldiv!(tmp2, invImDKDn, tmp1)
-        deta[:] = tmp2
+        ldiv!(invImDKDn, deta)
     end
     @trixi_timeit timer() "dv elliptic" begin
-        tmp2[:] = dv
-        ldiv!(tmp3, invImD2Kd, tmp2[2:(end - 1)])
+        ldiv!(invImD2Kd, dv[2:(end - 1)])
         dv[1] = dv[end] = zero(eltype(dv))
-        dv[2:(end - 1)] = tmp3
     end
 
     return nothing

From e4ab265784d7197e8f48783ac07c92bbad44164e Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 17:08:29 +0200
Subject: [PATCH 03/24] some fixes

---
 src/equations/bbm_bbm_1d.jl                     | 6 +++---
 src/equations/bbm_bbm_variable_bathymetry_1d.jl | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl
index 284424b1..6c3efaec 100644
--- a/src/equations/bbm_bbm_1d.jl
+++ b/src/equations/bbm_bbm_1d.jl
@@ -133,7 +133,7 @@ function create_cache(mesh, equations::BBMBBMEquations1D,
                       RealT, uEltype)
     D = equations.D
     invImD2 = lu(I - 1 / 6 * D^2 * sparse(solver.D2))
-    return (invImD2 = invImD2)
+    return (invImD2 = invImD2,)
 end
 
 function create_cache(mesh, equations::BBMBBMEquations1D,
@@ -172,7 +172,7 @@ end
 #   [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
 function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
               ::BoundaryConditionPeriodic, source_terms, solver, cache)
-    @unpack invImD2, tmp1, tmp2 = cache
+    @unpack invImD2 = cache
 
     eta = q.u[1]
     v = q.u[2]
@@ -246,7 +246,7 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
         ldiv!(invImD2n, deta)
     end
     @trixi_timeit timer() "dv elliptic" begin
-        ldiv!(invImD2d, dv[2:(end - 1)])
+        ldiv!(invImD2d, (@view dv[2:(end - 1)]))
         dv[1] = dv[end] = zero(eltype(dv))
     end
 
diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
index bb1e17c0..e3fb894e 100644
--- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl
+++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
@@ -327,7 +327,7 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMVariableEquations1D,
         ldiv!(invImDKDn, deta)
     end
     @trixi_timeit timer() "dv elliptic" begin
-        ldiv!(invImD2Kd, dv[2:(end - 1)])
+        ldiv!(invImD2Kd, (@view dv[2:(end - 1)]))
         dv[1] = dv[end] = zero(eltype(dv))
     end
 

From 06057590f9a4b553fd39efc7a5f68f27b437246a Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 17:19:07 +0200
Subject: [PATCH 04/24] use view instead of creating vector

---
 src/semidiscretization.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/semidiscretization.jl b/src/semidiscretization.jl
index f532b0e6..a05247e9 100644
--- a/src/semidiscretization.jl
+++ b/src/semidiscretization.jl
@@ -138,7 +138,7 @@ end
 
 function integrate_quantity!(quantity, func, q, semi::Semidiscretization)
     for i in eachnode(semi)
-        quantity[i] = func([q.u[v][i] for v in eachvariable(semi.equations)])
+        quantity[i] = func(view(q, i, :))
     end
     integrate(quantity, semi)
 end

From fd8d2dcd42c2e7758229a4a4075592cb069a9e4a Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 17:36:33 +0200
Subject: [PATCH 05/24] format

---
 src/solver.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/solver.jl b/src/solver.jl
index 6206fede..42c5d165 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -86,7 +86,8 @@ grid(solver::Solver) = grid(solver.D1)
 end
 
 function allocate_coefficients(mesh::Mesh1D, equations, solver::AbstractSolver)
-    return VectorOfArray([zeros(real(solver), nnodes(mesh)) for _ in eachvariable(equations)])
+    return VectorOfArray([zeros(real(solver), nnodes(mesh))
+                          for _ in eachvariable(equations)])
 end
 
 function compute_coefficients!(q, func, t, mesh::Mesh1D, equations, solver::AbstractSolver)

From 4e34bc40315688e6aebff57c788d1e02f431ed7c Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 17:43:04 +0200
Subject: [PATCH 06/24] bump lower compat of RecursiveArrayTools

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index d5427cf8..fc4d37a7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -30,7 +30,7 @@ LinearAlgebra = "1"
 PolynomialBases = "0.4.15"
 Printf = "1"
 RecipesBase = "1.1"
-RecursiveArrayTools = "3"
+RecursiveArrayTools = "3.22"
 Reexport = "1.0"
 Roots = "2.0.17"
 SciMLBase = "1.93, 2"

From 633738d1e71722994d72d9e6fb9a0ae5b98b9aa4 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 17:52:53 +0200
Subject: [PATCH 07/24] set lower compat of RecursiveArrayTools to 3.3

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index fc4d37a7..048873ce 100644
--- a/Project.toml
+++ b/Project.toml
@@ -30,7 +30,7 @@ LinearAlgebra = "1"
 PolynomialBases = "0.4.15"
 Printf = "1"
 RecipesBase = "1.1"
-RecursiveArrayTools = "3.22"
+RecursiveArrayTools = "3.3"
 Reexport = "1.0"
 Roots = "2.0.17"
 SciMLBase = "1.93, 2"

From c9dbd508edca4c85f3b14f8a84040aed0c3a9c64 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 17:56:39 +0200
Subject: [PATCH 08/24] remove compat for SciMLBase v1.93

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 048873ce..f80f7478 100644
--- a/Project.toml
+++ b/Project.toml
@@ -33,7 +33,7 @@ RecipesBase = "1.1"
 RecursiveArrayTools = "3.3"
 Reexport = "1.0"
 Roots = "2.0.17"
-SciMLBase = "1.93, 2"
+SciMLBase = "2"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 StaticArrays = "1"

From e0b3988544ee5d801c612ef7ab7ea2d2f901541c Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 18:05:03 +0200
Subject: [PATCH 09/24] bump lower compat of SciMLBase to 2.11

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index f80f7478..03867e17 100644
--- a/Project.toml
+++ b/Project.toml
@@ -33,7 +33,7 @@ RecipesBase = "1.1"
 RecursiveArrayTools = "3.3"
 Reexport = "1.0"
 Roots = "2.0.17"
-SciMLBase = "2"
+SciMLBase = "2.11"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 StaticArrays = "1"

From be754646a68fc05169e9a0c2538a6dc19e8cae4f Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 18:07:03 +0200
Subject: [PATCH 10/24] bump lower compat of DiffEqBase to 6.130

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 03867e17..1745bf0c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,7 +24,7 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 
 [compat]
 BandedMatrices = "0.17, 1"
-DiffEqBase = "6.128"
+DiffEqBase = "6.130"
 Interpolations = "0.14.2, 0.15"
 LinearAlgebra = "1"
 PolynomialBases = "0.4.15"

From 81e828dd150a308048d441b783b541d53dfb808e Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 18:12:01 +0200
Subject: [PATCH 11/24] bump lower compat of SummationByPartsOperators to
 0.5.50

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 1745bf0c..9c8f5a67 100644
--- a/Project.toml
+++ b/Project.toml
@@ -37,7 +37,7 @@ SciMLBase = "2.11"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 StaticArrays = "1"
-SummationByPartsOperators = "0.5.41"
+SummationByPartsOperators = "0.5.50"
 TimerOutputs = "0.5.7"
 TrixiBase = "0.1.3"
 julia = "1.9"

From d030e2356434e759c638a847f9300dc5053cddc3 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 18:14:44 +0200
Subject: [PATCH 12/24] bump lower compat of DiffEqBase to 6.143

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 9c8f5a67..767cb18b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,7 +24,7 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 
 [compat]
 BandedMatrices = "0.17, 1"
-DiffEqBase = "6.130"
+DiffEqBase = "6.143"
 Interpolations = "0.14.2, 0.15"
 LinearAlgebra = "1"
 PolynomialBases = "0.4.15"

From 22d492b2706eb1d2e19679a45121f84011db0237 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 18:18:34 +0200
Subject: [PATCH 13/24] bump lower compat of SummationByPartsOperators to
 0.5.52

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 767cb18b..73fd80d6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -37,7 +37,7 @@ SciMLBase = "2.11"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 StaticArrays = "1"
-SummationByPartsOperators = "0.5.50"
+SummationByPartsOperators = "0.5.52"
 TimerOutputs = "0.5.7"
 TrixiBase = "0.1.3"
 julia = "1.9"

From c3c0348bfe78a1c7dcf06310c01d04824bbce22c Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 18:21:12 +0200
Subject: [PATCH 14/24] remove compat 0.17 of BandedMatrices

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 73fd80d6..f9f9a342 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,7 +23,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 
 [compat]
-BandedMatrices = "0.17, 1"
+BandedMatrices = "1"
 DiffEqBase = "6.143"
 Interpolations = "0.14.2, 0.15"
 LinearAlgebra = "1"

From f6617db1e840faa23c07ee68529429e45c04d3b5 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 18:23:22 +0200
Subject: [PATCH 15/24] bump lower compat in tests of SummationByPartsOperators
 to 0.5.52

---
 test/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/Project.toml b/test/Project.toml
index 29b2aa9b..6ad4fe04 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -13,5 +13,5 @@ ExplicitImports = "1.0.1"
 OrdinaryDiffEq = "6.49.1"
 Plots = "1.23"
 SparseArrays = "1"
-SummationByPartsOperators = "0.5.41"
+SummationByPartsOperators = "0.5.52"
 Test = "1"

From fc04fee190fa721f9b3e5faa954f4b93d31dd7bb Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Mon, 17 Jun 2024 18:26:28 +0200
Subject: [PATCH 16/24] bump more compats

---
 test/Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/Project.toml b/test/Project.toml
index 6ad4fe04..c1f377de 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,8 +10,8 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 Aqua = "0.7, 0.8"
 ExplicitImports = "1.0.1"
-OrdinaryDiffEq = "6.49.1"
-Plots = "1.23"
+OrdinaryDiffEq = "6.56"
+Plots = "1.25"
 SparseArrays = "1"
 SummationByPartsOperators = "0.5.52"
 Test = "1"

From 6044989982a7487711ccbcb1edb5b90b12fbda2b Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Tue, 18 Jun 2024 11:10:44 +0200
Subject: [PATCH 17/24] use ArrayPartition instead of VectorOfArrays

---
 src/DispersiveShallowWater.jl                 |  2 +-
 src/equations/bbm_bbm_1d.jl                   | 16 +++++------
 .../bbm_bbm_variable_bathymetry_1d.jl         | 20 ++++++-------
 src/equations/svaerd_kalisch_1d.jl            | 21 +++++++-------
 src/semidiscretization.jl                     |  6 ++--
 src/solver.jl                                 | 28 +++++++++++++++----
 src/visualization.jl                          | 16 +++++------
 7 files changed, 63 insertions(+), 46 deletions(-)

diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl
index 568a51ad..19b94484 100644
--- a/src/DispersiveShallowWater.jl
+++ b/src/DispersiveShallowWater.jl
@@ -7,7 +7,7 @@ using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, lu, cholesky, ch
 using PolynomialBases: PolynomialBases
 using Printf: @printf, @sprintf
 using RecipesBase: RecipesBase, @recipe, @series
-using RecursiveArrayTools: VectorOfArray
+using RecursiveArrayTools: ArrayPartition
 using Reexport: @reexport
 using Roots: AlefeldPotraShi, find_zero
 
diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl
index 6c3efaec..298ddd7c 100644
--- a/src/equations/bbm_bbm_1d.jl
+++ b/src/equations/bbm_bbm_1d.jl
@@ -174,10 +174,10 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
               ::BoundaryConditionPeriodic, source_terms, solver, cache)
     @unpack invImD2 = cache
 
-    eta = q.u[1]
-    v = q.u[2]
-    deta = dq.u[1]
-    dv = dq.u[2]
+    eta = q.x[1]
+    v = q.x[2]
+    deta = dq.x[1]
+    dv = dq.x[2]
 
     D = equations.D
     # energy and mass conservative semidiscretization
@@ -218,10 +218,10 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
               ::BoundaryConditionReflecting, source_terms, solver, cache)
     @unpack invImD2d, invImD2n = cache
 
-    eta = q.u[1]
-    v = q.u[2]
-    deta = dq.u[1]
-    dv = dq.u[2]
+    eta = q.x[1]
+    v = q.x[2]
+    deta = dq.x[1]
+    dv = dq.x[2]
 
     D = equations.D
     # energy and mass conservative semidiscretization
diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
index e3fb894e..4e9bbd97 100644
--- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl
+++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
@@ -259,11 +259,11 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMVariableEquations1D,
               solver, cache)
     @unpack invImDKD, invImD2K, D = cache
 
-    eta = q.u[1]
-    v = q.u[2]
-    deta = dq.u[1]
-    dv = dq.u[2]
-    dD = dq.u[3]
+    eta = q.x[1]
+    v = q.x[2]
+    deta = dq.x[1]
+    dv = dq.x[2]
+    dD = dq.x[3]
     fill!(dD, zero(eltype(dD)))
 
     if solver.D1 isa PeriodicDerivativeOperator ||
@@ -298,11 +298,11 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMVariableEquations1D,
               solver, cache)
     @unpack invImDKDn, invImD2Kd, D = cache
 
-    eta = q.u[1]
-    v = q.u[2]
-    deta = dq.u[1]
-    dv = dq.u[2]
-    dD = dq.u[3]
+    eta = q.x[1]
+    v = q.x[2]
+    deta = dq.x[1]
+    dv = dq.x[2]
+    dD = dq.x[3]
     fill!(dD, zero(eltype(dD)))
 
     # energy and mass conservative semidiscretization
diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
index 84a78ad0..1b8e6cb4 100644
--- a/src/equations/svaerd_kalisch_1d.jl
+++ b/src/equations/svaerd_kalisch_1d.jl
@@ -204,11 +204,11 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D,
               solver, cache)
     @unpack factorization, D1betaD1, D, h, hv, alpha_hat, gamma_hat, tmp1, tmp2, D1_central, M = cache
 
-    eta = q.u[1]
-    v = q.u[2]
-    deta = dq.u[1]
-    dv = dq.u[2]
-    dD = dq.u[3]
+    eta = q.x[1]
+    v = q.x[2]
+    deta = dq.x[1]
+    dv = dq.x[2]
+    dD = dq.x[3]
     fill!(dD, zero(eltype(dD)))
 
     @trixi_timeit timer() "deta hyperbolic" begin
@@ -330,10 +330,11 @@ number of nodes as length of the second dimension.
 `cache` needs to hold the first-derivative SBP operator `D1`.
 """
 @inline function energy_total_modified(q, equations::SvaerdKalischEquations1D, cache)
-    e_modified = zeros(eltype(q), size(q, 1))
     # Need to compute new beta_hat, do not use the old one from the `cache`
-    v = q.u[2]
-    D = q.u[3]
+    v = q.x[2]
+    D = q.x[3]
+    N = length(v)
+    e_modified = zeros(eltype(q), N)
     beta_hat = equations.beta * D .^ 3
     if cache.D1 isa PeriodicDerivativeOperator ||
        cache.D1 isa UniformPeriodicCoupledOperator
@@ -343,8 +344,8 @@ number of nodes as length of the second dimension.
     else
         @error "unknown type of first-derivative operator: $(typeof(cache.D1))"
     end
-    for i in 1:size(q, 1)
-        e_modified[i] = energy_total(view(q, i, :), equations) + tmp[i]
+    for i in 1:N
+        e_modified[i] = energy_total(get_node_vars(q, equations, i), equations) + tmp[i]
     end
     return e_modified
 end
diff --git a/src/semidiscretization.jl b/src/semidiscretization.jl
index a05247e9..dbf4ef52 100644
--- a/src/semidiscretization.jl
+++ b/src/semidiscretization.jl
@@ -117,10 +117,10 @@ Get the grid of a semidiscretization.
 """
 grid(semi::Semidiscretization) = grid(semi.solver)
 
-function PolynomialBases.integrate(func, q::VectorOfArray, semi::Semidiscretization)
+function PolynomialBases.integrate(func, q::ArrayPartition, semi::Semidiscretization)
     integrals = zeros(real(semi), nvariables(semi))
     for v in eachvariable(semi.equations)
-        integrals[v] = integrate(func, q.u[v], semi.solver.D1)
+        integrals[v] = integrate(func, q.x[v], semi.solver.D1)
     end
     return integrals
 end
@@ -138,7 +138,7 @@ end
 
 function integrate_quantity!(quantity, func, q, semi::Semidiscretization)
     for i in eachnode(semi)
-        quantity[i] = func(view(q, i, :))
+        quantity[i] = func(get_node_vars(q, semi.equations, i))
     end
     integrate(quantity, semi)
 end
diff --git a/src/solver.jl b/src/solver.jl
index 42c5d165..0b9bd9e4 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -78,16 +78,32 @@ grid(solver::Solver) = grid(solver.D1)
 @inline eachnode(solver::Solver) = Base.OneTo(length(grid(solver)))
 @inline Base.real(::Solver{RealT}) where {RealT} = RealT
 
+# Adapted from Trixi.jl
+# https://github.com/trixi-framework/Trixi.jl/blob/75d8c67629562efd24b2a04e46d22b0a1f4f572c/src/solvers/dg.jl#L539
+@inline function get_node_vars(q, equations, indices...)
+    # There is a cut-off at `n == 10` inside of the method
+    # `ntuple(f::F, n::Integer) where F` in Base at ntuple.jl:17
+    # in Julia `v1.5`, leading to type instabilities if
+    # more than ten variables are used. That's why we use
+    # `Val(...)` below.
+    # We use `@inline` to make sure that the `getindex` calls are
+    # really inlined, which might be the default choice of the Julia
+    # compiler for standard `Array`s but not necessarily for more
+    # advanced array types such as `PtrArray`s, cf.
+    # https://github.com/JuliaSIMD/VectorizationBase.jl/issues/55
+    SVector(ntuple(@inline(v->q.x[v][indices...]), Val(nvariables(equations))))
+end
+
 @inline function set_node_vars!(q, q_node, equations, indices...)
     for v in eachvariable(equations)
-        q.u[v][indices...] = q_node[v]
+        q.x[v][indices...] = q_node[v]
     end
     return nothing
 end
 
 function allocate_coefficients(mesh::Mesh1D, equations, solver::AbstractSolver)
-    return VectorOfArray([zeros(real(solver), nnodes(mesh))
-                          for _ in eachvariable(equations)])
+    return ArrayPartition([zeros(real(solver), nnodes(mesh))
+                           for _ in eachvariable(equations)]...)
 end
 
 function compute_coefficients!(q, func, t, mesh::Mesh1D, equations, solver::AbstractSolver)
@@ -108,7 +124,7 @@ function calc_error_norms(q, t, initial_condition, mesh::Mesh1D, equations,
     l2_error = zeros(real(solver), nvariables(equations))
     linf_error = similar(l2_error)
     for v in eachvariable(equations)
-        @views diff = q.u[v] - q_exact[v, :]
+        @views diff = q.x[v] - q_exact[v, :]
         l2_error[v] = integrate(q -> q^2, diff, solver.D1) |> sqrt
         linf_error[v] = maximum(abs.(diff))
     end
@@ -124,9 +140,9 @@ function calc_sources!(dq, q, t, source_terms,
                        equations::AbstractEquations{1}, solver::Solver)
     x = grid(solver)
     for i in eachnode(solver)
-        local_source = source_terms(view(q, i, :), x[i], t, equations)
+        local_source = source_terms(get_node_vars(q, equations, i), x[i], t, equations)
         for v in eachvariable(equations)
-            dq.u[v][i] += local_source[v]
+            dq.x[v][i] += local_source[v]
         end
     end
     return nothing
diff --git a/src/visualization.jl b/src/visualization.jl
index cffd2003..de627a70 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -43,15 +43,15 @@ end
     end
     for j in eachnode(semi)
         if plot_bathymetry == true
-            bathy[j] = bathymetry(view(q, j, :), equations)
+            bathy[j] = bathymetry(get_node_vars(q, equations, j), equations)
         end
         if plot_initial == true
-            q_exact[j, :] .= conversion(q_exact[j, :], equations)
+            set_node_vars!(q_exact, conversion(get_node_vars(q_exact, equations, j), equations), equations, j)
         end
-        data[:, j] .= conversion(view(q, j, :), equations)
+        data[:, j] .= conversion(get_node_vars(q, equations, j), equations)
     end
 
-    plot_title --> "$(get_name(semi.equations)) at t = $(round(t, digits=5))"
+    plot_title --> "$(get_name(equations)) at t = $(round(t, digits=5))"
     layout --> nsubplots
 
     for i in 1:nsubplots
@@ -63,7 +63,7 @@ end
                 subplot --> i
                 linestyle := :solid
                 label := "initial $(names[i])"
-                grid(semi), q_exact[:, i]
+                grid(semi), q_exact.x[i]
             end
         end
 
@@ -107,12 +107,12 @@ end
 
     solution = zeros(nvariables(semi), length(sol.t))
     data = zeros(nvars, length(sol.t))
-    for i in 1:nvariables(semi)
+    for v in 1:nvariables(semi)
         for k in 1:length(sol.t)
             # Allow that the spatial value `x` is not on the grid. Thus, interpolate the given values to the provided `x`
             # with a linear spline.
-            solution[i, k] = linear_interpolation(grid(semi),
-                                                  view(sol.u[k], :, i))(x)
+            solution[v, k] = linear_interpolation(grid(semi),
+                                                  sol.u[k].x[v])(x)
         end
     end
 

From cc4e3e54b265755a28ee3e9308723b0ca1de2387 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Tue, 18 Jun 2024 11:12:38 +0200
Subject: [PATCH 18/24] introduce eachvariable(semi)

---
 src/semidiscretization.jl | 3 ++-
 src/visualization.jl      | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/semidiscretization.jl b/src/semidiscretization.jl
index dbf4ef52..b13599e7 100644
--- a/src/semidiscretization.jl
+++ b/src/semidiscretization.jl
@@ -106,6 +106,7 @@ end
 
 @inline Base.ndims(semi::Semidiscretization) = ndims(semi.mesh)
 @inline nvariables(semi::Semidiscretization) = nvariables(semi.equations)
+@inline eachvariable(semi::Semidiscretization) = eachvariable(semi.equations)
 @inline nnodes(semi::Semidiscretization) = nnodes(semi.mesh)
 @inline eachnode(semi::Semidiscretization) = eachnode(semi.mesh)
 @inline Base.real(semi::Semidiscretization) = real(semi.solver)
@@ -119,7 +120,7 @@ grid(semi::Semidiscretization) = grid(semi.solver)
 
 function PolynomialBases.integrate(func, q::ArrayPartition, semi::Semidiscretization)
     integrals = zeros(real(semi), nvariables(semi))
-    for v in eachvariable(semi.equations)
+    for v in eachvariable(semi)
         integrals[v] = integrate(func, q.x[v], semi.solver.D1)
     end
     return integrals
diff --git a/src/visualization.jl b/src/visualization.jl
index de627a70..a9e8ff9f 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -107,7 +107,7 @@ end
 
     solution = zeros(nvariables(semi), length(sol.t))
     data = zeros(nvars, length(sol.t))
-    for v in 1:nvariables(semi)
+    for v in eachvariable(semi)
         for k in 1:length(sol.t)
             # Allow that the spatial value `x` is not on the grid. Thus, interpolate the given values to the provided `x`
             # with a linear spline.

From 859f45c2f310ef768958a6042bd2ca02f9375b68 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Tue, 18 Jun 2024 11:14:54 +0200
Subject: [PATCH 19/24] format

---
 src/visualization.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/visualization.jl b/src/visualization.jl
index a9e8ff9f..48ed2179 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -46,7 +46,9 @@ end
             bathy[j] = bathymetry(get_node_vars(q, equations, j), equations)
         end
         if plot_initial == true
-            set_node_vars!(q_exact, conversion(get_node_vars(q_exact, equations, j), equations), equations, j)
+            set_node_vars!(q_exact,
+                           conversion(get_node_vars(q_exact, equations, j), equations),
+                           equations, j)
         end
         data[:, j] .= conversion(get_node_vars(q, equations, j), equations)
     end

From ba147c66595f60088aa45ae73cfb2554b4eb6def Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Tue, 18 Jun 2024 11:21:41 +0200
Subject: [PATCH 20/24] fix plotting with conversion

---
 src/visualization.jl | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/visualization.jl b/src/visualization.jl
index 48ed2179..e4a51d44 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -34,6 +34,7 @@ end
 
     if plot_initial == true
         q_exact = compute_coefficients(initial_condition, t, semi)
+        data_exact = zeros(nvars, nnodes(semi))
     end
 
     q = sol.u[step]
@@ -46,9 +47,7 @@ end
             bathy[j] = bathymetry(get_node_vars(q, equations, j), equations)
         end
         if plot_initial == true
-            set_node_vars!(q_exact,
-                           conversion(get_node_vars(q_exact, equations, j), equations),
-                           equations, j)
+            data_exact[:, j] .= conversion(get_node_vars(q_exact, equations, j), equations)
         end
         data[:, j] .= conversion(get_node_vars(q, equations, j), equations)
     end
@@ -65,7 +64,7 @@ end
                 subplot --> i
                 linestyle := :solid
                 label := "initial $(names[i])"
-                grid(semi), q_exact.x[i]
+                grid(semi), data_exact[i, :]
             end
         end
 

From cb0c62f29b964af886a1a36c9e493b0d482a0ef5 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Tue, 18 Jun 2024 11:27:41 +0200
Subject: [PATCH 21/24] bump lower compat of RecipesBase to 1.2

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index f9f9a342..ae068e56 100644
--- a/Project.toml
+++ b/Project.toml
@@ -29,7 +29,7 @@ Interpolations = "0.14.2, 0.15"
 LinearAlgebra = "1"
 PolynomialBases = "0.4.15"
 Printf = "1"
-RecipesBase = "1.1"
+RecipesBase = "1.2"
 RecursiveArrayTools = "3.3"
 Reexport = "1.0"
 Roots = "2.0.17"

From 5f672fe2aeb2eccd3afb6a1769670f3f6d301f8a Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Tue, 18 Jun 2024 11:31:21 +0200
Subject: [PATCH 22/24] bump lower compat in test of OrdinaryDiffEq to 6.62

---
 test/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/Project.toml b/test/Project.toml
index c1f377de..9c9acb41 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,7 +10,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 Aqua = "0.7, 0.8"
 ExplicitImports = "1.0.1"
-OrdinaryDiffEq = "6.56"
+OrdinaryDiffEq = "6.62"
 Plots = "1.25"
 SparseArrays = "1"
 SummationByPartsOperators = "0.5.52"

From e94e833718bf214ccc6d8087d5a0ad290b651494 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Wed, 19 Jun 2024 13:10:43 +0200
Subject: [PATCH 23/24] add allocation tests

---
 test/test_bbm_bbm_1d.jl                     | 14 ++++++++++++++
 test/test_bbm_bbm_variable_bathymetry_1d.jl | 20 +++++++++++++++++++-
 test/test_svaerd_kalisch_1d.jl              | 12 ++++++++++++
 test/test_util.jl                           | 15 +++++++++++++++
 4 files changed, 60 insertions(+), 1 deletion(-)

diff --git a/test/test_bbm_bbm_1d.jl b/test/test_bbm_bbm_1d.jl
index e862b807..ea0229cc 100644
--- a/test/test_bbm_bbm_1d.jl
+++ b/test/test_bbm_bbm_1d.jl
@@ -19,6 +19,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d")
                             change_entropy=0.00019552920957721653,
                             atol_ints=1e-10) # in order to make CI pass
 
+        @test_allocations(semi, sol)
+
         # test upwind operators
         using SummationByPartsOperators: upwind_operators, periodic_derivative_operator
         using SparseArrays: sparse
@@ -47,6 +49,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d")
         for (linf_expected, linf_actual) in zip(linf, linf_measured)
             @test isapprox(linf_expected, linf_actual, atol = atol, rtol = rtol)
         end
+
+        @test_allocations(semi, sol)
     end
 
     @trixi_testset "bbm_bbm_1d_dg" begin
@@ -58,6 +62,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d")
                             change_waterheight=1.0543424823256011e-14,
                             change_velocity=-3.552713678800501e-15,
                             change_entropy=-0.021843627302246205)
+
+        @test_allocations(semi, sol, allocs = 60000)
     end
 
     @trixi_testset "bbm_bbm_1d_relaxation" begin
@@ -69,6 +75,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d")
                             change_waterheight=2.1746572188776938e-13,
                             change_velocity=-4.547473508864641e-13,
                             change_entropy=0.0)
+
+        @test_allocations(semi, sol)
     end
 
     @trixi_testset "bbm_bbm_1d_manufactured" begin
@@ -82,6 +90,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d")
                             change_entropy=17.387441847193436,
                             atol=1e-10,
                             atol_ints=1e-10) # in order to make CI pass
+
+        @test_allocations(semi, sol)
     end
 
     @trixi_testset "bbm_bbm_1d_basic_reflecting" begin
@@ -94,6 +104,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d")
                             change_velocity=0.5469460931577768,
                             change_entropy=130.69415963528348)
 
+        @test_allocations(semi, sol)
+
         # test upwind operators
         using SummationByPartsOperators: upwind_operators, Mattsson2017
         using SparseArrays: sparse
@@ -123,6 +135,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d")
         for (linf_expected, linf_actual) in zip(linf, linf_measured)
             @test isapprox(linf_expected, linf_actual, atol = atol, rtol = rtol)
         end
+
+        @test_allocations(semi, sol)
     end
 end
 
diff --git a/test/test_bbm_bbm_variable_bathymetry_1d.jl b/test/test_bbm_bbm_variable_bathymetry_1d.jl
index 9a9f4f35..836630c0 100644
--- a/test/test_bbm_bbm_variable_bathymetry_1d.jl
+++ b/test/test_bbm_bbm_variable_bathymetry_1d.jl
@@ -19,6 +19,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                             change_velocity=-3.410605131648481e-13,
                             change_entropy=0.0006197172569955001,
                             atol_ints=1e-10) # in order to make CI pass
+
+        @test_allocations(semi, sol)
     end
 
     @trixi_testset "bbm_bbm_variable_bathymetry_1d_relaxation" begin
@@ -31,6 +33,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                             change_waterheight=-2.1271873151818e-13,
                             change_velocity=2.7874259106214573e-13,
                             change_entropy=0.0)
+
+        @test_allocations(semi, sol)
     end
 
     @trixi_testset "bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation" begin
@@ -43,6 +47,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                             change_waterheight=3.1796787425264483e-13,
                             change_velocity=-1.0584519372081047e-12,
                             change_entropy=1.4210854715202004e-14)
+
+        @test_allocations(semi, sol, allocs = 60000)
     end
 
     @trixi_testset "bbm_bbm_variable_bathymetry_1d_upwind_relaxation" begin
@@ -55,7 +61,9 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                             change_waterheight=-1.5765166949677223e-13,
                             change_velocity=1.8791999206735355e-13,
                             change_entropy=2.1316282072803006e-14,
-                            atol=1e-11) # in order to make CI pass)
+                            atol=1e-11) # in order to make CI pass
+
+        @test_allocations(semi, sol)
     end
 
     @trixi_testset "bbm_bbm_variable_bathymetry_1d_manufactured" begin
@@ -70,6 +78,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                             change_entropy=17.817012269328853,
                             atol=1e-10,
                             atol_ints=1e-10) # in order to make CI pass
+
+        @test_allocations(semi, sol)
     end
 
     @trixi_testset "bbm_bbm_variable_bathymetry_1d_well_balanced" begin
@@ -83,6 +93,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                             change_velocity=-5.728935765975126e-15,
                             change_entropy=0.0,
                             lake_at_rest=6.725397799854067e-15)
+
+        @test_allocations(semi, sol)
     end
 
     @trixi_testset "bbm_bbm_variable_bathymetry_1d_dingemans" begin
@@ -96,6 +108,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                             change_waterheight=-1.4210854715202004e-13,
                             change_velocity=-3.1478183774857893e-15,
                             change_entropy=3.417533493699221e-7)
+
+        @test_allocations(semi, sol)
     end
 
     @trixi_testset "bbm_bbm_variable_bathymetry_1d_basic_reflecting" begin
@@ -111,6 +125,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                             atol=1e-10,
                             atol_ints=1e-10) # in order to make CI pass
 
+        @test_allocations(semi, sol)
+
         # test upwind operators
         using SummationByPartsOperators: upwind_operators, Mattsson2017
         using SparseArrays: sparse
@@ -140,6 +156,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
         for (linf_expected, linf_actual) in zip(linf, linf_measured)
             @test isapprox(linf_expected, linf_actual, atol = atol, rtol = rtol)
         end
+
+        @test_allocations(semi, sol)
     end
 end
 
diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl
index f120f3f2..164cef94 100644
--- a/test/test_svaerd_kalisch_1d.jl
+++ b/test/test_svaerd_kalisch_1d.jl
@@ -18,6 +18,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_waterheight=-2.42861286636753e-16,
                             change_entropy=0.1868146724821993,
                             atol=1e-9) # in order to make CI pass
+
+    @test_allocations(semi, sol, allocs = 130_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans" begin
@@ -31,6 +33,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_waterheight=-3.979039320256561e-13,
                             change_entropy=-0.00024362648639453255,
                             change_entropy_modified=-6.311893230304122e-9)
+
+    @test_allocations(semi, sol, allocs = 500_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans_cg" begin
@@ -43,6 +47,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_waterheight=-1.4210854715202004e-13,
                             change_entropy=-0.0002425303440531934,
                             change_entropy_modified=-2.6815314413397573e-9)
+
+    @test_allocations(semi, sol, allocs = 1_000_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans_upwind" begin
@@ -55,6 +61,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_waterheight=-1.1368683772161603e-13,
                             change_entropy=-0.00023645232727176335,
                             change_entropy_modified=-6.654090611846186e-9)
+
+    @test_allocations(semi, sol, allocs = 600_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans_relaxation" begin
@@ -68,6 +76,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_waterheight=-3.979039320256561e-13,
                             change_entropy=-0.00024362054875837202,
                             change_entropy_modified=0.0)
+
+    @test_allocations(semi, sol, allocs = 500_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_well_balanced" begin
@@ -81,6 +91,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_momentum=1.5679986322667355e-14,
                             change_entropy=0.0,
                             lake_at_rest=0.0)
+
+    @test_allocations(semi, sol, allocs = 200_000)
     end
 end
 
diff --git a/test/test_util.jl b/test/test_util.jl
index fc40501a..eba0f81d 100644
--- a/test/test_util.jl
+++ b/test/test_util.jl
@@ -190,6 +190,21 @@ macro trixi_testset(name, expr)
     end
 end
 
+"""
+    @test_allocations(semi, sol, allocs = 50000)
+
+Test that the memory allocations of `DispersiveShallowWater.rhs!` are below `allocs`
+(e.g., from type instabilities).
+"""
+macro test_allocations(semi, sol, allocs = 50000)
+    quote
+        t = $sol.t[end]
+        q = $sol.u[end]
+        dq = similar(q)
+        @test (@allocated DispersiveShallowWater.rhs!(dq, q, $semi, t)) < $allocs
+    end
+end
+
 # Copied from TrixiBase.jl. See https://github.com/trixi-framework/TrixiBase.jl/issues/9.
 """
     @test_nowarn_mod expr

From 0f5182c4724f55156cae77d531ed679e5d4fa016 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Wed, 19 Jun 2024 13:17:14 +0200
Subject: [PATCH 24/24] format

---
 test/test_bbm_bbm_1d.jl                     |  2 +-
 test/test_bbm_bbm_variable_bathymetry_1d.jl |  2 +-
 test/test_svaerd_kalisch_1d.jl              | 12 ++++++------
 3 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/test/test_bbm_bbm_1d.jl b/test/test_bbm_bbm_1d.jl
index ea0229cc..236acafe 100644
--- a/test/test_bbm_bbm_1d.jl
+++ b/test/test_bbm_bbm_1d.jl
@@ -63,7 +63,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d")
                             change_velocity=-3.552713678800501e-15,
                             change_entropy=-0.021843627302246205)
 
-        @test_allocations(semi, sol, allocs = 60000)
+        @test_allocations(semi, sol, allocs=60000)
     end
 
     @trixi_testset "bbm_bbm_1d_relaxation" begin
diff --git a/test/test_bbm_bbm_variable_bathymetry_1d.jl b/test/test_bbm_bbm_variable_bathymetry_1d.jl
index 836630c0..3a271c81 100644
--- a/test/test_bbm_bbm_variable_bathymetry_1d.jl
+++ b/test/test_bbm_bbm_variable_bathymetry_1d.jl
@@ -48,7 +48,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                             change_velocity=-1.0584519372081047e-12,
                             change_entropy=1.4210854715202004e-14)
 
-        @test_allocations(semi, sol, allocs = 60000)
+        @test_allocations(semi, sol, allocs=60000)
     end
 
     @trixi_testset "bbm_bbm_variable_bathymetry_1d_upwind_relaxation" begin
diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl
index 164cef94..bd98adb0 100644
--- a/test/test_svaerd_kalisch_1d.jl
+++ b/test/test_svaerd_kalisch_1d.jl
@@ -19,7 +19,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_entropy=0.1868146724821993,
                             atol=1e-9) # in order to make CI pass
 
-    @test_allocations(semi, sol, allocs = 130_000)
+        @test_allocations(semi, sol, allocs=130_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans" begin
@@ -34,7 +34,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_entropy=-0.00024362648639453255,
                             change_entropy_modified=-6.311893230304122e-9)
 
-    @test_allocations(semi, sol, allocs = 500_000)
+        @test_allocations(semi, sol, allocs=500_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans_cg" begin
@@ -48,7 +48,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_entropy=-0.0002425303440531934,
                             change_entropy_modified=-2.6815314413397573e-9)
 
-    @test_allocations(semi, sol, allocs = 1_000_000)
+        @test_allocations(semi, sol, allocs=1_000_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans_upwind" begin
@@ -62,7 +62,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_entropy=-0.00023645232727176335,
                             change_entropy_modified=-6.654090611846186e-9)
 
-    @test_allocations(semi, sol, allocs = 600_000)
+        @test_allocations(semi, sol, allocs=600_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans_relaxation" begin
@@ -77,7 +77,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_entropy=-0.00024362054875837202,
                             change_entropy_modified=0.0)
 
-    @test_allocations(semi, sol, allocs = 500_000)
+        @test_allocations(semi, sol, allocs=500_000)
     end
 
     @trixi_testset "svaerd_kalisch_1d_well_balanced" begin
@@ -92,7 +92,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             change_entropy=0.0,
                             lake_at_rest=0.0)
 
-    @test_allocations(semi, sol, allocs = 200_000)
+        @test_allocations(semi, sol, allocs=200_000)
     end
 end