diff --git a/Project.toml b/Project.toml
index 285d29ce639..e2d1873af79 100644
--- a/Project.toml
+++ b/Project.toml
@@ -22,7 +22,6 @@ LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
-Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
@@ -82,7 +81,6 @@ LoopVectorization = "0.12.151"
 MPI = "0.20"
 Makie = "0.19, 0.20"
 MuladdMacro = "0.2.2"
-Octavian = "0.3.21"
 OffsetArrays = "1.12"
 P4est = "0.4.9"
 Polyester = "0.7.10"
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 0cedec782df..f5d3fca48a0 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -53,7 +53,6 @@ using LinearMaps: LinearMap
 using LoopVectorization: LoopVectorization, @turbo, indices
 using StaticArrayInterface: static_length # used by LoopVectorization
 using MuladdMacro: @muladd
-using Octavian: Octavian, matmul!
 using Polyester: Polyester, @batch # You know, the cheapest threads you can find...
 using OffsetArrays: OffsetArray, OffsetVector
 using P4est
diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index 2d588c5c79d..6063be9e12a 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -6,7 +6,7 @@
 #! format: noindent
 
 # out <- A*x
-mul_by!(A) = @inline (out, x) -> matmul!(out, A, x)
+mul_by!(A) = @inline (out, x) -> mul!(out, A, x)
 mul_by!(A::T) where {T <: SimpleKronecker} = @inline (out, x) -> mul!(out, A, x)
 mul_by!(A::AbstractSparseMatrix) = @inline (out, x) -> mul!(out, A, x)
 function mul_by!(A::LinearAlgebra.AdjOrTrans{T, S}) where {T, S <: AbstractSparseMatrix}
@@ -14,7 +14,7 @@ function mul_by!(A::LinearAlgebra.AdjOrTrans{T, S}) where {T, S <: AbstractSpars
 end
 
 #  out <- out + α * A * x
-mul_by_accum!(A, α) = @inline (out, x) -> matmul!(out, A, x, α, One())
+mul_by_accum!(A, α) = @inline (out, x) -> mul!(out, A, x, α, One())
 function mul_by_accum!(A::AbstractSparseMatrix, α)
     @inline (out, x) -> mul!(out, A, x, α, One())
 end
@@ -33,6 +33,11 @@ mul_by_accum!(A::UniformScaling) = MulByAccumUniformScaling()
     StructArrays.foreachfield(f, args...)
 end
 
+# Matrix{<:SVector} fallback
+@inline function apply_to_each_field(f::F, args::Vararg{Any, N}) where {F, N}
+    f(args...)
+end
+
 # specialize for UniformScaling types: works for either StructArray{SVector} or Matrix{SVector}
 # solution storage formats.
 @inline apply_to_each_field(f::MulByUniformScaling, out, x, args...) = copy!(out, x)
@@ -136,12 +141,8 @@ function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys, ValueT
     return boundary_conditions
 end
 
-# Allocate nested array type for DGMulti solution storage.
 function allocate_nested_array(uEltype, nvars, array_dimensions, dg)
-    # store components as separate arrays, combine via StructArrays
-    return StructArray{SVector{nvars, uEltype}}(ntuple(_ -> zeros(uEltype,
-                                                                  array_dimensions...),
-                                                       nvars))
+    return zeros(SVector{nvars, uEltype}, array_dimensions...)
 end
 
 function reset_du!(du, dg::DGMulti, other_args...)
diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl
index 88f06607019..aca950ecf2d 100644
--- a/src/solvers/dgmulti/flux_differencing.jl
+++ b/src/solvers/dgmulti/flux_differencing.jl
@@ -288,12 +288,6 @@ function compute_flux_differencing_SBP_matrices(dg::DGMultiFluxDiffSBP,
     return Qrst_skew
 end
 
-# For flux differencing SBP-type approximations, store solutions in Matrix{SVector{nvars}}.
-# This results in a slight speedup for `calc_volume_integral!`.
-function allocate_nested_array(uEltype, nvars, array_dimensions, dg::DGMultiFluxDiffSBP)
-    return zeros(SVector{nvars, uEltype}, array_dimensions...)
-end
-
 function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiffSBP, RealT,
                       uEltype)
     rd = dg.basis