diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
index 8518d202..41a71d06 100644
--- a/.JuliaFormatter.toml
+++ b/.JuliaFormatter.toml
@@ -6,3 +6,5 @@ yas_style_nesting = true
 
 # Align struct fields for better readability of large struct definitions
 align_struct_field = true
+
+always_use_return = true
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index f43fe2be..6eef7699 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -23,6 +23,7 @@ for elixir in elixirs
     println("Running $benchname...")
     redirect_stdout(devnull) do
         trixi_include(@__MODULE__, elixir, tspan = (0.0, 1e-11))
+        return nothing
     end
     SUITE[benchname] = @benchmarkable DispersiveShallowWater.rhs!($(similar(sol.u[end])),
                                                                   $(copy(sol.u[end])),
diff --git a/src/boundary_conditions.jl b/src/boundary_conditions.jl
index a28dc225..a76d0d4d 100644
--- a/src/boundary_conditions.jl
+++ b/src/boundary_conditions.jl
@@ -18,4 +18,5 @@ const boundary_condition_reflecting = BoundaryConditionReflecting()
 
 function Base.show(io::IO, ::BoundaryConditionReflecting)
     print(io, "boundary_condition_reflecting")
+    return nothing
 end
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 7f044c03..5dec2226 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -37,6 +37,7 @@ function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:AnalysisCallback})
     analysis_callback = cb.affect!
     @unpack interval = analysis_callback
     print(io, "AnalysisCallback(interval=", interval, ")")
+    return nothing
 end
 
 function Base.show(io::IO, ::MIME"text/plain",
@@ -62,7 +63,7 @@ end
 
 function AnalysisCallback(semi::Semidiscretization; kwargs...)
     mesh, equations, solver = mesh_equations_solver(semi)
-    AnalysisCallback(mesh, equations, solver; kwargs...)
+    return AnalysisCallback(mesh, equations, solver; kwargs...)
 end
 
 function AnalysisCallback(mesh, equations::AbstractEquations, solver;
@@ -96,9 +97,9 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver;
                                          Vector{Vector{real(solver)}}(),
                                          io)
 
-    DiscreteCallback(condition, analysis_callback,
-                     save_positions = (false, false),
-                     initialize = initialize!)
+    return DiscreteCallback(condition, analysis_callback,
+                            save_positions = (false, false),
+                            initialize = initialize!)
 end
 
 """
@@ -318,7 +319,7 @@ end
 # terminate the type-stable iteration over tuples
 function analyze_integrals!(io, current_integrals, i, analysis_integrals::Tuple{}, q, t,
                             semi)
-    nothing
+    return nothing
 end
 
 # used for error checks and EOC analysis
@@ -333,7 +334,7 @@ function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition,
 end
 
 function analyze!(semi::Semidiscretization, quantity, q, t)
-    integrate_quantity!(semi.cache.tmp1, quantity, q, semi)
+    return integrate_quantity!(semi.cache.tmp1, quantity, q, semi)
 end
 
 pretty_form_utf(::typeof(waterheight)) = "∫h"
diff --git a/src/callbacks_step/relaxation.jl b/src/callbacks_step/relaxation.jl
index 6a36b33b..e73ce42e 100644
--- a/src/callbacks_step/relaxation.jl
+++ b/src/callbacks_step/relaxation.jl
@@ -21,6 +21,7 @@ function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:RelaxationCallback})
     relaxation_callback = cb.affect!
     @unpack invariant = relaxation_callback
     print(io, "RelaxationCallback(invariant=", string(nameof(invariant)), ")")
+    return nothing
 end
 
 function Base.show(io::IO, ::MIME"text/plain",
@@ -40,9 +41,9 @@ end
 function RelaxationCallback(; invariant)
     relaxation_callback = RelaxationCallback(invariant)
 
-    DiscreteCallback(relaxation_callback, relaxation_callback, # the first one is the condition, the second the affect!
-                     save_positions = (false, false),
-                     initialize = initialize!)
+    return DiscreteCallback(relaxation_callback, relaxation_callback, # the first one is the condition, the second the affect!
+                            save_positions = (false, false),
+                            initialize = initialize!)
 end
 
 function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t,
diff --git a/src/callbacks_step/summary.jl b/src/callbacks_step/summary.jl
index 8bb4b60a..e528424f 100644
--- a/src/callbacks_step/summary.jl
+++ b/src/callbacks_step/summary.jl
@@ -9,19 +9,19 @@ struct SummaryCallback
 
     function SummaryCallback(io::IO = stdout)
         function initialize(cb, u, t, integrator)
-            initialize_summary_callback(cb, u, t, integrator)
+            return initialize_summary_callback(cb, u, t, integrator)
         end
         # At the end of the simulation, the timer is printed
         function finalize(cb, u, t, integrator)
-            finalize_summary_callback(cb, u, t, integrator)
+            return finalize_summary_callback(cb, u, t, integrator)
         end
         summary_callback = new(io)
         # SummaryCallback is never called during the simulation
         condition = (u, t, integrator) -> false
-        DiscreteCallback(condition, summary_callback,
-                         save_positions = (false, false),
-                         initialize = initialize,
-                         finalize = finalize)
+        return DiscreteCallback(condition, summary_callback,
+                                save_positions = (false, false),
+                                initialize = initialize,
+                                finalize = finalize)
     end
 end
 
@@ -29,6 +29,7 @@ function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SummaryCallback})
     @nospecialize cb # reduce precompilation time
 
     print(io, "SummaryCallback")
+    return nothing
 end
 
 function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator)
diff --git a/src/dispersion_relation.jl b/src/dispersion_relation.jl
index 8a509685..49d025af 100644
--- a/src/dispersion_relation.jl
+++ b/src/dispersion_relation.jl
@@ -15,6 +15,7 @@ end
 
 function Base.show(io::IO, disp_rel::LinearDispersionRelation)
     print(io, "LinearDispersionRelation(h0 = ", disp_rel.ref_height, ")")
+    return nothing
 end
 
 Base.broadcastable(disp_rel::LinearDispersionRelation) = (disp_rel,)
diff --git a/src/equations/bbm_1d.jl b/src/equations/bbm_1d.jl
index 6fd4b672..c7010165 100644
--- a/src/equations/bbm_1d.jl
+++ b/src/equations/bbm_1d.jl
@@ -50,7 +50,7 @@ end
 
 function BBMEquation1D(; gravity, D = 1.0, eta0 = 0.0, split_form = true)
     eta0 == 0.0 || @warn "The still-water surface needs to be 0 for the BBM equations"
-    BBMEquation1D(gravity, D, eta0, split_form)
+    return BBMEquation1D(gravity, D, eta0, split_form)
 end
 
 function check_solver(::BBMEquation1D, solver, boundary_conditions)
diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl
index d5172620..817ef25f 100644
--- a/src/equations/bbm_bbm_1d.jl
+++ b/src/equations/bbm_bbm_1d.jl
@@ -61,7 +61,7 @@ end
 function BBMBBMEquations1D(; bathymetry_type = bathymetry_variable,
                            gravity, eta0 = 0.0)
     eta0 == 0.0 || @warn "The still-water surface needs to be 0 for the BBM-BBM equations"
-    BBMBBMEquations1D(bathymetry_type, gravity, eta0)
+    return BBMBBMEquations1D(bathymetry_type, gravity, eta0)
 end
 
 function check_solver(::BBMBBMEquations1D, solver, boundary_conditions)
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index cebe699f..4224b45b 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -569,15 +569,16 @@ function solve_system_matrix!(dv, system_matrix, rhs,
                                                SerreGreenNaghdiEquations1D},
                               D1, cache, ::BoundaryConditionPeriodic)
     scale_by_mass_matrix!(rhs, D1)
-    solve_system_matrix!(dv, system_matrix, rhs, equations, D1, cache)
+    return solve_system_matrix!(dv, system_matrix, rhs, equations, D1, cache)
 end
 
 function solve_system_matrix!(dv, system_matrix, rhs,
                               equations::SvaerdKalischEquations1D,
                               D1, cache, ::BoundaryConditionReflecting)
     scale_by_mass_matrix!(rhs, D1)
-    solve_system_matrix!(dv, system_matrix, (@view rhs[(begin + 1):(end - 1)]), equations,
-                         D1, cache)
+    return solve_system_matrix!(dv, system_matrix, (@view rhs[(begin + 1):(end - 1)]),
+                                equations,
+                                D1, cache)
 end
 
 function solve_system_matrix!(dv, system_matrix, rhs,
@@ -585,7 +586,7 @@ function solve_system_matrix!(dv, system_matrix, rhs,
                               D1, cache, ::BoundaryConditionReflecting)
     scale_by_mass_matrix!(rhs, D1)
     rhs[begin] = rhs[end] = 0
-    solve_system_matrix!(dv, system_matrix, rhs, equations, D1, cache)
+    return solve_system_matrix!(dv, system_matrix, rhs, equations, D1, cache)
 end
 
 function solve_system_matrix!(dv, system_matrix, rhs,
@@ -611,7 +612,7 @@ function solve_system_matrix!(dv, system_matrix, rhs,
 end
 
 function solve_system_matrix!(dv, system_matrix, ::Union{BBMEquation1D, BBMBBMEquations1D})
-    ldiv!(system_matrix, dv)
+    return ldiv!(system_matrix, dv)
 end
 
 # To match the experimental data from Dingemans, the initial condition needs to be shifted
diff --git a/src/equations/hyperbolic_serre_green_naghdi_1d.jl b/src/equations/hyperbolic_serre_green_naghdi_1d.jl
index 57ddc4ac..34f2840f 100644
--- a/src/equations/hyperbolic_serre_green_naghdi_1d.jl
+++ b/src/equations/hyperbolic_serre_green_naghdi_1d.jl
@@ -98,7 +98,7 @@ function HyperbolicSerreGreenNaghdiEquations1D(; bathymetry_type = bathymetry_mi
                                                gravity,
                                                eta0 = 0.0,
                                                lambda)
-    HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type, gravity, eta0, lambda)
+    return HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type, gravity, eta0, lambda)
 end
 
 function varnames(::typeof(prim2prim), ::HyperbolicSerreGreenNaghdiEquations1D)
diff --git a/src/equations/kdv_1d.jl b/src/equations/kdv_1d.jl
index bb9e6ca3..3de80075 100644
--- a/src/equations/kdv_1d.jl
+++ b/src/equations/kdv_1d.jl
@@ -50,7 +50,7 @@ end
 
 function KdVEquation1D(; gravity, D = 1.0, eta0 = 0.0)
     eta0 == 0.0 || @warn "The still-water surface needs to be 0 for the KdV equations"
-    KdVEquation1D(gravity, D, eta0)
+    return KdVEquation1D(gravity, D, eta0)
 end
 
 function check_solver(::KdVEquation1D, solver, boundary_conditions)
diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl
index 22bcdd46..971b42db 100644
--- a/src/equations/serre_green_naghdi_1d.jl
+++ b/src/equations/serre_green_naghdi_1d.jl
@@ -78,7 +78,7 @@ end
 
 function SerreGreenNaghdiEquations1D(; bathymetry_type = bathymetry_variable,
                                      gravity, eta0 = 0.0)
-    SerreGreenNaghdiEquations1D(bathymetry_type, gravity, eta0)
+    return SerreGreenNaghdiEquations1D(bathymetry_type, gravity, eta0)
 end
 
 """
diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
index 3f501f7a..9128b286 100644
--- a/src/equations/svaerd_kalisch_1d.jl
+++ b/src/equations/svaerd_kalisch_1d.jl
@@ -72,7 +72,7 @@ const SvärdKalischEquations1D = SvaerdKalischEquations1D
 
 function SvaerdKalischEquations1D(; bathymetry_type = bathymetry_variable, gravity,
                                   eta0 = 0.0, alpha = 0.0, beta = 1 / 3, gamma = 0.0)
-    SvaerdKalischEquations1D(bathymetry_type, gravity, eta0, alpha, beta, gamma)
+    return SvaerdKalischEquations1D(bathymetry_type, gravity, eta0, alpha, beta, gamma)
 end
 
 function check_solver(::SvaerdKalischEquations1D, solver, ::BoundaryConditionPeriodic)
diff --git a/src/mesh.jl b/src/mesh.jl
index 093d2af5..830a4cf2 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -11,7 +11,7 @@ struct Mesh1D{RealT}
     function Mesh1D{RealT}(xmin::RealT, xmax::RealT, N::Int) where {RealT}
         @assert xmin < xmax
         @assert N > 0
-        new(xmin, xmax, N)
+        return new(xmin, xmax, N)
     end
 end
 
@@ -22,11 +22,12 @@ Create a simple homogeneous one-dimensional mesh from `xmin` to `xmax` with `N`
 """
 function Mesh1D(xmin, xmax, N)
     xmin, xmax = promote(xmin, xmax)
-    Mesh1D{typeof(xmin)}(xmin, xmax, N)
+    return Mesh1D{typeof(xmin)}(xmin, xmax, N)
 end
 
 function Base.show(io::IO, mesh::Mesh1D{RealT}) where {RealT}
     print(io, "Mesh1D{", RealT, "} with length ", mesh.N)
+    return nothing
 end
 
 function Base.show(io::IO, ::MIME"text/plain", mesh::Mesh1D{RealT}) where {RealT}
diff --git a/src/semidiscretization.jl b/src/semidiscretization.jl
index 6cd442a2..c8f56f34 100644
--- a/src/semidiscretization.jl
+++ b/src/semidiscretization.jl
@@ -35,8 +35,8 @@ struct Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions,
         @assert xmax(mesh) == xmax(solver.D1)
         @assert nnodes(mesh) == length(grid(solver))
 
-        new(mesh, equations, initial_condition, boundary_conditions, source_terms, solver,
-            cache)
+        return new(mesh, equations, initial_condition, boundary_conditions, source_terms,
+                   solver, cache)
     end
 end
 
@@ -82,20 +82,22 @@ function Semidiscretization(mesh, equations, initial_condition, solver;
                           RealT, uEltype)...,
              initial_cache...)
 
-    Semidiscretization{typeof(mesh), typeof(equations), typeof(initial_condition),
-                       typeof(boundary_conditions), typeof(source_terms),
-                       typeof(solver), typeof(cache)}(mesh, equations, initial_condition,
-                                                      boundary_conditions, source_terms,
-                                                      solver, cache)
+    return Semidiscretization{typeof(mesh), typeof(equations), typeof(initial_condition),
+                              typeof(boundary_conditions), typeof(source_terms),
+                              typeof(solver), typeof(cache)}(mesh, equations,
+                                                             initial_condition,
+                                                             boundary_conditions,
+                                                             source_terms,
+                                                             solver, cache)
 end
 
 """
     check_solver(equations, solver, boundary_conditions)
 
-Check that the `solver` is compatible with the given `equations` and 
+Check that the `solver` is compatible with the given `equations` and
 `boundary_conditions`. The default implementation performs no checks.
-Specific equation types can override this method to validate that 
-required derivative operators are present (e.g., some equations 
+Specific equation types can override this method to validate that
+required derivative operators are present (e.g., some equations
 require `D2` or `D3` to be non-`nothing`).
 
 Throws an `ArgumentError` if the solver is incompatible.
@@ -118,6 +120,7 @@ function Base.show(io::IO, semi::Semidiscretization)
         print(io, key)
     end
     print(io, "))")
+    return nothing
 end
 
 function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization)
@@ -158,22 +161,22 @@ function PolynomialBases.integrate(func, q::ArrayPartition, semi::Semidiscretiza
     return integrals
 end
 function PolynomialBases.integrate(func, quantity, semi::Semidiscretization)
-    integrate(func, quantity, semi.solver.D1)
+    return integrate(func, quantity, semi.solver.D1)
 end
 function PolynomialBases.integrate(q, semi::Semidiscretization)
-    integrate(identity, q, semi)
+    return integrate(identity, q, semi)
 end
 
 function integrate_quantity(func, q, semi)
     quantity = zeros(eltype(q), nnodes(semi))
-    integrate_quantity!(quantity, func, q, semi)
+    return integrate_quantity!(quantity, func, q, semi)
 end
 
 function integrate_quantity!(quantity, func, q, semi)
     for i in eachnode(semi)
         quantity[i] = func(get_node_vars(q, semi.equations, i), semi.equations)
     end
-    integrate(quantity, semi)
+    return integrate(quantity, semi)
 end
 
 # Obtain the function, which has an additional `!` appended to the name
@@ -188,7 +191,7 @@ function integrate_quantity!(quantity,
                                          typeof(hamiltonian)}, q,
                              semi)
     inplace_version(func)(quantity, q, semi.equations, semi.cache)
-    integrate(quantity, semi)
+    return integrate(quantity, semi)
 end
 
 @inline function mesh_equations_solver(semi)
@@ -202,7 +205,7 @@ end
 end
 
 function calc_error_norms(q, t, semi)
-    calc_error_norms(q, t, semi.initial_condition, mesh_equations_solver(semi)...)
+    return calc_error_norms(q, t, semi.initial_condition, mesh_equations_solver(semi)...)
 end
 
 function rhs!(dq, q, semi, t)
@@ -224,9 +227,9 @@ end
 function compute_coefficients!(q, func, t, semi)
     # Call `compute_coefficients` defined by the solver
     mesh, equations, solver = mesh_equations_solver(semi)
-    compute_coefficients!(q, func, t, mesh,
-                          is_hyperbolic_appproximation(equations), equations,
-                          solver)
+    return compute_coefficients!(q, func, t, mesh,
+                                 is_hyperbolic_appproximation(equations), equations,
+                                 solver)
 end
 
 check_bathymetry(equations, q0) = nothing
@@ -271,6 +274,7 @@ function jacobian(semi::Semidiscretization;
                   q0 = compute_coefficients(semi.initial_condition, t, semi))
     J = ForwardDiff.jacobian(similar(q0), q0) do dq, q
         DispersiveShallowWater.rhs!(dq, q, semi, t)
+        return nothing
     end
     return J
 end
diff --git a/src/solver.jl b/src/solver.jl
index 73fb052c..6d8c291f 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -46,7 +46,7 @@ struct Solver{RealT <: Real,
             @assert derivative_order(D3) == 3
         end
 
-        new(D1, D2, D3)
+        return new(D1, D2, D3)
     end
 end
 
@@ -69,7 +69,7 @@ function Solver(mesh, accuracy_order)
     D2 = periodic_derivative_operator(2, accuracy_order, mesh.xmin, mesh.xmax, mesh.N)
     D3 = periodic_derivative_operator(3, accuracy_order + 2, mesh.xmin, mesh.xmax, mesh.N)
     @assert real(D1) == real(D2)
-    Solver{real(D1), typeof(D1), typeof(D2), typeof(D3)}(D1, D2, D3)
+    return Solver{real(D1), typeof(D1), typeof(D2), typeof(D3)}(D1, D2, D3)
 end
 
 # Also allow to pass custom SBP operators (for convenience without explicitly specifying the type)
@@ -89,11 +89,12 @@ function Solver(D1::AbstractDerivativeOperator{RealT},
                           Nothing} = nothing,
                 D3::Union{AbstractDerivativeOperator{RealT}, AbstractMatrix{RealT},
                           Nothing} = nothing) where {RealT}
-    Solver{RealT, typeof(D1), typeof(D2), typeof(D3)}(D1, D2, D3)
+    return Solver{RealT, typeof(D1), typeof(D2), typeof(D3)}(D1, D2, D3)
 end
 
 function Base.show(io::IO, solver::Solver{RealT}) where {RealT}
     print(io, "Solver{", RealT, "}")
+    return nothing
 end
 
 function Base.show(io::IO, ::MIME"text/plain", solver::Solver{RealT}) where {RealT}
@@ -125,7 +126,7 @@ grid(solver::Solver) = grid(solver.D1)
     # 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))))
+    return SVector(ntuple(@inline(v->q.x[v][indices...]), Val(nvariables(equations))))
 end
 
 @inline function set_node_vars!(q, q_node, equations, indices...)
diff --git a/src/util.jl b/src/util.jl
index 06854caa..db7df191 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -44,7 +44,8 @@ See also [`examples_dir`](@ref) and [`get_examples`](@ref).
 Copied from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl).
 """
 function default_example()
-    joinpath(examples_dir(), "serre_green_naghdi_1d", "serre_green_naghdi_soliton.jl")
+    return joinpath(examples_dir(), "serre_green_naghdi_1d",
+                    "serre_green_naghdi_soliton.jl")
 end
 
 """
@@ -61,7 +62,7 @@ readdir(data_dir())
 data_dir() = pkgdir(DispersiveShallowWater, "data")::String
 
 function convergence_test(example::AbstractString, iterations_or_Ns; kwargs...)
-    convergence_test(Main, example::AbstractString, iterations_or_Ns; kwargs...)
+    return convergence_test(Main, example::AbstractString, iterations_or_Ns; kwargs...)
 end
 
 """
@@ -84,7 +85,7 @@ function convergence_test(mod::Module, example::AbstractString, iterations; io::
 
     initial_N = extract_initial_N(example, kwargs)
     Ns = initial_N * 2 .^ (0:(iterations - 1))
-    convergence_test(mod, example, Ns; io = io, kwargs...)
+    return convergence_test(mod, example, Ns; io = io, kwargs...)
 end
 
 function convergence_test(mod::Module, example::AbstractString, Ns::AbstractVector;
@@ -112,7 +113,7 @@ function convergence_test(mod::Module, example::AbstractString, Ns::AbstractVect
     end
 
     # Use raw error values to compute EOC
-    analyze_convergence(io, errors, iterations, (@invokelatest mod.semi), Ns)
+    return analyze_convergence(io, errors, iterations, (@invokelatest mod.semi), Ns)
 end
 
 # Analyze convergence for any semidiscretization
@@ -120,7 +121,7 @@ end
 function analyze_convergence(io, errors, iterations, semi::Semidiscretization, Ns)
     _, equations, _, _ = mesh_equations_solver_cache(semi)
     variablenames = varnames(prim2prim, equations)
-    analyze_convergence(io, errors, iterations, variablenames, Ns)
+    return analyze_convergence(io, errors, iterations, variablenames, Ns)
 end
 
 # This method is called with the collected error values to actually compute and print the EOC
diff --git a/src/visualization.jl b/src/visualization.jl
index f5889671..52958f62 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -158,23 +158,25 @@ end
 @recipe function f(semisol::Pair{<:Semidiscretization, <:ODESolution}; plot_initial = false,
                    plot_analytical = false, plot_bathymetry = true, conversion = prim2phys,
                    step = -1)
-    PlotData(semisol, plot_initial, plot_analytical, plot_bathymetry, conversion, step)
+    return PlotData(semisol, plot_initial, plot_analytical, plot_bathymetry, conversion,
+                    step)
 end
 
 @recipe function f(semi::Semidiscretization, sol::ODESolution; plot_initial = false,
                    plot_analytical = false, plot_bathymetry = true, conversion = prim2phys,
                    step = -1)
-    PlotData(semi => sol, plot_initial, plot_analytical, plot_bathymetry, conversion, step)
+    return PlotData(semi => sol, plot_initial, plot_analytical, plot_bathymetry, conversion,
+                    step)
 end
 
 @recipe function f(semisol::Pair{<:Semidiscretization, <:ODESolution}, x_value;
                    conversion = prim2phys)
-    PlotDataOverTime(semisol, x_value, conversion)
+    return PlotDataOverTime(semisol, x_value, conversion)
 end
 
 @recipe function f(semi::Semidiscretization, sol::ODESolution, x_value;
                    conversion = prim2phys)
-    PlotDataOverTime(semi => sol, x_value, conversion)
+    return PlotDataOverTime(semi => sol, x_value, conversion)
 end
 
 function pretty_form_utf(name)