diff --git a/examples/t8code_2d_fv/elixir_advection_basic.jl b/examples/t8code_2d_fv/elixir_advection_basic.jl
new file mode 100644
index 00000000000..12f9363892f
--- /dev/null
+++ b/examples/t8code_2d_fv/elixir_advection_basic.jl
@@ -0,0 +1,52 @@
+#!/usr/bin/env julia
+# Has to be at least Julia v1.9.0.
+
+# Running in parallel:
+#   ${JULIA_DEPOT_PATH}/.julia/bin/mpiexecjl --project=. -n 3 julia hybrid-t8code-mesh.jl
+#
+# More information: https://juliaparallel.org/MPI.jl/stable/usage/
+using OrdinaryDiffEq
+using Trixi
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+solver = FV(order = 1, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (8, 8)
+
+mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                  initial_refinement_level = 1)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+ode = semidiscretize(semi, (0.0, 1.0));
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback, save_solution, analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback()
diff --git a/src/meshes/t8code_fv_mesh.jl b/src/meshes/t8code_fv_mesh.jl
index f941212e749..14597104d10 100644
--- a/src/meshes/t8code_fv_mesh.jl
+++ b/src/meshes/t8code_fv_mesh.jl
@@ -111,19 +111,6 @@ function Base.show(io::IO, ::MIME"text/plain", mesh::T8codeFVMesh)
     end
 end
 
-function create_cache(mesh::T8codeFVMesh, equations,
-                      solver, RealT, uEltype)
-    elements = init_elements(mesh, RealT, uEltype)
-
-    interfaces = init_interfaces(mesh, equations, elements)
-
-    u_ = init_solution!(mesh, equations)
-
-    cache = (; elements, interfaces, u_)
-
-    return cache
-end
-
 # Write the forest as vtu and also write the element's volumes in the file.
 #
 # t8code supports writing element based data to vtu as long as its stored
diff --git a/src/solvers/fv_t8code/containers.jl b/src/solvers/fv_t8code/containers.jl
index 95b9b768c50..48cbb64c013 100644
--- a/src/solvers/fv_t8code/containers.jl
+++ b/src/solvers/fv_t8code/containers.jl
@@ -52,17 +52,7 @@ function Base.show(container::T8codeElementContainer)
     println("neighbor_faces     = ", container.neighbor_faces[1:num_faces])
 end
 
-function init_elements(mesh::T8codeFVMesh, RealT, uEltype)
-    # Initialize container
-    elements = init_elements(mesh)
-
-    # Exchange the neighboring data at MPI process boundaries.
-    exchange_ghost_data(mesh, elements)
-
-    return elements
-end
-
-function init_elements(mesh::T8codeFVMesh)
+function init_elements(mesh::T8codeFVMesh, uEltype)
     @unpack forest = mesh
     # Check that the forest is a committed.
     @assert(t8_forest_is_committed(forest)==1)
@@ -81,6 +71,19 @@ function init_elements(mesh::T8codeFVMesh)
                                                                         num_local_elements +
                                                                         num_ghost_elements)
 
+    init_elements!(elements, mesh)
+
+    # Exchange the neighboring data at MPI process boundaries.
+    exchange_ghost_data(mesh, elements)
+
+    return elements
+end
+
+function init_elements!(elements, mesh::T8codeFVMesh)
+    @unpack forest = mesh
+    n_dims = ndims(mesh)
+    @unpack max_number_faces = mesh
+
     midpoint = Vector{Cdouble}(undef, n_dims)
 
     face_midpoints = Matrix{Cdouble}(undef, 3, max_number_faces) # Need NDIMS=3 for t8code API. Also, consider that Julia is column major.
@@ -175,7 +178,7 @@ function init_elements(mesh::T8codeFVMesh)
         end
     end
 
-    return elements
+    return nothing
 end
 
 mutable struct T8codeInterfaceContainer{NDIMS, uEltype <: Real} <: AbstractContainer
diff --git a/src/solvers/fv_t8code/fv_2d.jl b/src/solvers/fv_t8code/fv_2d.jl
index 76e3d9918ad..8e894e7a54d 100644
--- a/src/solvers/fv_t8code/fv_2d.jl
+++ b/src/solvers/fv_t8code/fv_2d.jl
@@ -45,6 +45,19 @@ Base.summary(io::IO, solver::FV) = print(io, "FV(order=$(solver.order))")
     ndofs(mesh, solver, cache)
 end
 
+function create_cache(mesh::T8codeFVMesh, equations,
+                      solver, RealT, uEltype)
+    elements = init_elements(mesh, uEltype)
+
+    interfaces = init_interfaces(mesh, equations, elements)
+
+    u_ = init_solution!(mesh, equations)
+
+    cache = (; elements, interfaces, u_)
+
+    return cache
+end
+
 function compute_coefficients!(u, func, t, mesh::T8codeFVMesh, equations,
                                solver::FV, cache)
     for element in eachelement(mesh, solver, cache)