From 58fc2b05d701cfae3054e21238490cc10157a1e9 Mon Sep 17 00:00:00 2001
From: Niklas Neher <73897120+LasNikas@users.noreply.github.com>
Date: Tue, 6 Feb 2024 14:56:27 +0100
Subject: [PATCH 1/4] `WendlandC2Kernel` for `TotalLagrangianSPHSystem` (#371)

* change kernel for tlsph

* implement suggestions

* implement suggestions
---
 examples/fsi/dam_break_2d.jl          |  4 ++--
 examples/fsi/dam_break_gate_2d.jl     |  4 ++--
 examples/solid/oscillating_beam_2d.jl | 30 +++++++++++++++++++++++----
 3 files changed, 30 insertions(+), 8 deletions(-)

diff --git a/examples/fsi/dam_break_2d.jl b/examples/fsi/dam_break_2d.jl
index 1e886aff6..9e6a58506 100644
--- a/examples/fsi/dam_break_2d.jl
+++ b/examples/fsi/dam_break_2d.jl
@@ -88,8 +88,8 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Solid
-solid_smoothing_length = sqrt(2) * solid_particle_spacing
-solid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+solid_smoothing_length = 2 * sqrt(2) * solid_particle_spacing
+solid_smoothing_kernel = WendlandC2Kernel{2}()
 
 # For the FSI we need the hydrodynamic masses and densities in the solid boundary model
 hydrodynamic_densites = fluid_density * ones(size(solid.density))
diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl
index 0300cd529..d3a3e9d3a 100644
--- a/examples/fsi/dam_break_gate_2d.jl
+++ b/examples/fsi/dam_break_gate_2d.jl
@@ -117,8 +117,8 @@ boundary_system_gate = BoundarySPHSystem(gate, boundary_model_gate, movement=gat
 
 # ==========================================================================================
 # ==== Solid
-solid_smoothing_length = sqrt(2) * solid_particle_spacing
-solid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+solid_smoothing_length = 2 * sqrt(2) * solid_particle_spacing
+solid_smoothing_kernel = WendlandC2Kernel{2}()
 
 # For the FSI we need the hydrodynamic masses and densities in the solid boundary model
 hydrodynamic_densites = fluid_density * ones(size(solid.density))
diff --git a/examples/solid/oscillating_beam_2d.jl b/examples/solid/oscillating_beam_2d.jl
index 3791197ff..fc1ffc1c9 100644
--- a/examples/solid/oscillating_beam_2d.jl
+++ b/examples/solid/oscillating_beam_2d.jl
@@ -45,9 +45,10 @@ solid = union(beam, fixed_particles)
 
 # ==========================================================================================
 # ==== Solid
-
-smoothing_length = sqrt(2) * particle_spacing
-smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+# The kernel in the reference uses a differently scaled smoothing length,
+# so this is equivalent to the smoothing length of `sqrt(2) * particle_spacing` used in the paper.
+smoothing_length = 2 * sqrt(2) * particle_spacing
+smoothing_kernel = WendlandC2Kernel{2}()
 
 solid_system = TotalLagrangianSPHSystem(solid,
                                         smoothing_kernel, smoothing_length,
@@ -62,7 +63,28 @@ semi = Semidiscretization(solid_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)
-saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+
+# Track the position of the particle in the middle of the tip of the beam.
+particle_id = Int(n_particles_per_dimension[1] * (n_particles_per_dimension[2] + 1) / 2)
+
+shift_x = beam.coordinates[1, particle_id]
+shift_y = beam.coordinates[2, particle_id]
+
+function x_deflection(v, u, t, system)
+    particle_position = TrixiParticles.current_coords(u, system, particle_id)
+
+    return particle_position[1] - shift_x
+end
+
+function y_deflection(v, u, t, system)
+    particle_position = TrixiParticles.current_coords(u, system, particle_id)
+
+    return particle_position[2] - shift_y
+end
+
+saving_callback = SolutionSavingCallback(dt=0.02, prefix="",
+                                         x_deflection=x_deflection,
+                                         y_deflection=y_deflection)
 
 callbacks = CallbackSet(info_callback, saving_callback)
 

From f220877c7b333b1383ee6cc8f6bfdd9ff644ef7f Mon Sep 17 00:00:00 2001
From: Sven Berger <berger.sven@gmail.com>
Date: Tue, 6 Feb 2024 16:56:57 +0100
Subject: [PATCH 2/4] Save additional output for TLSPH (#348)

* write missing values

* replace incorrect nomenclature

* calc von mises stress

* rm comment

* format

* change docs to F

* add test

* format

* also save the cauchy stress components

* format

* implement suggestions

* format

* update

* update

* format

* add citation

---------

Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>
---
 .gitignore                                    |  1 +
 examples/solid/oscillating_beam_2d.jl         |  8 +--
 .../solid/total_lagrangian_sph/system.jl      | 66 +++++++++++++++----
 src/visualization/write2vtk.jl                | 14 ++++
 test/systems/solid_system.jl                  | 34 ++++++++++
 5 files changed, 107 insertions(+), 16 deletions(-)

diff --git a/.gitignore b/.gitignore
index be4f7441f..f7151b63f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,3 +18,4 @@ run
 *.vtu
 *.pvd
 *.mp4
+*.csv
diff --git a/examples/solid/oscillating_beam_2d.jl b/examples/solid/oscillating_beam_2d.jl
index fc1ffc1c9..8d541debb 100644
--- a/examples/solid/oscillating_beam_2d.jl
+++ b/examples/solid/oscillating_beam_2d.jl
@@ -50,12 +50,10 @@ solid = union(beam, fixed_particles)
 smoothing_length = 2 * sqrt(2) * particle_spacing
 smoothing_kernel = WendlandC2Kernel{2}()
 
-solid_system = TotalLagrangianSPHSystem(solid,
-                                        smoothing_kernel, smoothing_length,
-                                        E, nu,
+solid_system = TotalLagrangianSPHSystem(solid, smoothing_kernel, smoothing_length,
+                                        E, nu, nothing, #No boundary model
                                         n_fixed_particles=nparticles(fixed_particles),
-                                        acceleration=(0.0, -gravity),
-                                        nothing) # No boundary model
+                                        acceleration=(0.0, -gravity))
 
 # ==========================================================================================
 # ==== Simulation
diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl
index 85d1cc45c..e75ec6a45 100644
--- a/src/schemes/solid/total_lagrangian_sph/system.jl
+++ b/src/schemes/solid/total_lagrangian_sph/system.jl
@@ -33,15 +33,15 @@ The zero subscript on quantities denotes that the quantity is to be measured in
 The difference in the initial coordinates is denoted by $\bm{X}_{ab} = \bm{X}_a - \bm{X}_b$,
 the difference in the current coordinates is denoted by $\bm{x}_{ab} = \bm{x}_a - \bm{x}_b$.
 
-For the computation of the PK1 stress tensor, the deformation gradient $\bm{J}$ is computed per particle as
+For the computation of the PK1 stress tensor, the deformation gradient $\bm{F}$ is computed per particle as
 ```math
-\bm{J}_a = \sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ba} (\bm{L}_{0a}\nabla_{0a} W(\bm{X}_{ab}))^T \\
+\bm{F}_a = \sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ba} (\bm{L}_{0a}\nabla_{0a} W(\bm{X}_{ab}))^T \\
     \qquad  = -\left(\sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ab} (\nabla_{0a} W(\bm{X}_{ab}))^T \right) \bm{L}_{0a}^T
 ```
 with $1 \leq i,j \leq d$.
 From the deformation gradient, the Green-Lagrange strain
 ```math
-\bm{E} = \frac{1}{2}(\bm{J}^T\bm{J} - \bm{I})
+\bm{E} = \frac{1}{2}(\bm{F}^T\bm{F} - \bm{I})
 ```
 and the second Piola-Kirchhoff stress tensor
 ```math
@@ -49,7 +49,7 @@ and the second Piola-Kirchhoff stress tensor
 ```
 are computed to obtain the PK1 stress tensor as
 ```math
-\bm{P} = \bm{J}\bm{S}.
+\bm{P} = \bm{F}\bm{S}.
 ```
 
 Here,
@@ -285,8 +285,8 @@ end
     calc_deformation_grad!(deformation_grad, neighborhood_search, system)
 
     @threaded for particle in eachparticle(system)
-        J_particle = deformation_gradient(system, particle)
-        pk1_particle = pk1_stress_tensor(J_particle, system)
+        F_particle = deformation_gradient(system, particle)
+        pk1_particle = pk1_stress_tensor(F_particle, system)
         pk1_particle_corrected = pk1_particle * correction_matrix(system, particle)
 
         @inbounds for j in 1:ndims(system), i in 1:ndims(system)
@@ -332,18 +332,18 @@ end
 end
 
 # First Piola-Kirchhoff stress tensor
-@inline function pk1_stress_tensor(J, system)
-    S = pk2_stress_tensor(J, system)
+@inline function pk1_stress_tensor(F, system)
+    S = pk2_stress_tensor(F, system)
 
-    return J * S
+    return F * S
 end
 
 # Second Piola-Kirchhoff stress tensor
-@inline function pk2_stress_tensor(J, system)
+@inline function pk2_stress_tensor(F, system)
     (; lame_lambda, lame_mu) = system
 
     # Compute the Green-Lagrange strain
-    E = 0.5 * (transpose(J) * J - I)
+    E = 0.5 * (transpose(F) * F - I)
 
     return lame_lambda * tr(E) * I + 2 * lame_mu * E
 end
@@ -411,3 +411,47 @@ end
 function viscosity_model(system::TotalLagrangianSPHSystem)
     return system.boundary_model.viscosity
 end
+
+# An explanation of these equation can be found in
+# J. Lubliner, 2008. Plasticity theory.
+# See here below Equation 5.3.21 for the equation for the equivalent stress.
+# The von-Mises stress is one form of equivalent stress, where sigma is the deviatoric stress.
+# See pages 32 and 123.
+function von_mises_stress(system::TotalLagrangianSPHSystem)
+    von_mises_stress = zeros(eltype(system.pk1_corrected), nparticles(system))
+
+    @threaded for particle in each_moving_particle(system)
+        F = deformation_gradient(system, particle)
+        J = det(F)
+        P = pk1_corrected(system, particle)
+        sigma = (1.0 / J) * P * F'
+
+        # Calculate deviatoric stress tensor
+        s = sigma - (1.0 / 3.0) * tr(sigma) * I
+
+        von_mises_stress[particle] = sqrt(3.0 / 2.0 * sum(s .^ 2))
+    end
+
+    return von_mises_stress
+end
+
+# An explanation of these equation can be found in
+# J. Lubliner, 2008. Plasticity theory.
+# See here page 473 for the relation between the `pk1`, the first Piola-Kirchhoff tensor,
+# and the Cauchy stress.
+function cauchy_stress(system::TotalLagrangianSPHSystem)
+    NDIMS = ndims(system)
+
+    cauchy_stress_tensors = zeros(eltype(system.pk1_corrected), NDIMS, NDIMS,
+                                  nparticles(system))
+
+    @threaded for particle in each_moving_particle(system)
+        F = deformation_gradient(system, particle)
+        J = det(F)
+        P = pk1_corrected(system, particle)
+        sigma = (1.0 / J) * P * F'
+        cauchy_stress_tensors[:, :, particle] = sigma
+    end
+
+    return cauchy_stress_tensors
+end
diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl
index 11b88524a..bad071243 100644
--- a/src/visualization/write2vtk.jl
+++ b/src/visualization/write2vtk.jl
@@ -216,11 +216,25 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d
 
     vtk["velocity"] = hcat(view(v, 1:ndims(system), :),
                            zeros(ndims(system), n_fixed_particles))
+    vtk["jacobian"] = [det(deformation_gradient(system, particle))
+                       for particle in eachparticle(system)]
+
+    vtk["von_mises_stress"] = von_mises_stress(system)
+
+    sigma = cauchy_stress(system)
+    vtk["sigma_11"] = sigma[1, 1, :]
+    vtk["sigma_22"] = sigma[2, 2, :]
+    if ndims(system) == 3
+        vtk["sigma_33"] = sigma[3, 3, :]
+    end
+
     vtk["material_density"] = system.material_density
     vtk["young_modulus"] = system.young_modulus
     vtk["poisson_ratio"] = system.poisson_ratio
     vtk["lame_lambda"] = system.lame_lambda
     vtk["lame_mu"] = system.lame_mu
+    vtk["smoothing_kernel"] = type2string(system.smoothing_kernel)
+    vtk["smoothing_length"] = system.smoothing_length
 
     write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data)
 end
diff --git a/test/systems/solid_system.jl b/test/systems/solid_system.jl
index a486aac66..7e69f0041 100644
--- a/test/systems/solid_system.jl
+++ b/test/systems/solid_system.jl
@@ -335,4 +335,38 @@
 
         @test v0 == velocity
     end
+
+    @testset verbose=true "compute_von_mises_stress" begin
+        # System setup
+        coordinates = [1.0 2.0; 1.0 2.0]
+        velocity = zero(coordinates)
+        mass = [1.25, 1.5]
+        material_densities = [990.0, 1000.0]
+        smoothing_kernel = Val(:smoothing_kernel)
+        smoothing_length = 0.362
+        nu = 0.25  # Poisson's ratio
+        E = 2.5    # Young's modulus
+        boundary_model = Val(:boundary_model)
+
+        initial_condition = InitialCondition(; coordinates, velocity, mass,
+                                             density=material_densities)
+        system = TotalLagrangianSPHSystem(initial_condition, smoothing_kernel,
+                                          smoothing_length, E, nu, boundary_model)
+
+        # Initialize deformation_grad and pk1_corrected with arbitrary values
+        for particle in TrixiParticles.eachparticle(system)
+            system.deformation_grad[:, :, particle] = [1.0 0.2; 0.2 1.0]
+            system.pk1_corrected[:, :, particle] = [1.0 0.5; 0.5 1.0]
+        end
+
+        von_mises_stress = TrixiParticles.von_mises_stress(system)
+        cauchy_stress = TrixiParticles.cauchy_stress(system)
+
+        reference_stress_tensor = [1.145833 0.729167; 0.729167 1.145833;;;
+                                   1.145833 0.729167; 0.729167 1.145833]
+
+        # Verify against calculation by hand
+        @test isapprox(von_mises_stress[1], 1.4257267477533202, atol=1e-14)
+        @test isapprox(reference_stress_tensor, cauchy_stress, atol=1e-6)
+    end
 end

From 878c2b49694fd7a10bc5fac731f5dd44bd6ea9ee Mon Sep 17 00:00:00 2001
From: Sven Berger <berger.sven@gmail.com>
Date: Wed, 7 Feb 2024 09:48:57 +0100
Subject: [PATCH 3/4] Switch integration method for fluid examples (#358)

* switch to RDPK3SpFSAL35

* rm

---------

Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
---
 examples/fluid/accelerated_tank_2d.jl      | 2 +-
 examples/fluid/dam_break_2d.jl             | 2 +-
 examples/fluid/dam_break_3d.jl             | 2 +-
 examples/fluid/falling_water_column_2d.jl  | 2 +-
 examples/fluid/moving_wall_2d.jl           | 2 +-
 examples/fluid/periodic_channel_2d.jl      | 2 +-
 examples/fluid/rectangular_tank_2d.jl      | 2 +-
 examples/fluid/rectangular_tank_edac_2d.jl | 2 +-
 8 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/examples/fluid/accelerated_tank_2d.jl b/examples/fluid/accelerated_tank_2d.jl
index 2cc797fbe..67c0da7ef 100644
--- a/examples/fluid/accelerated_tank_2d.jl
+++ b/examples/fluid/accelerated_tank_2d.jl
@@ -77,7 +77,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
 # Sometimes, the method fails to do so because forces become extremely large when
 # fluid particles are very close to boundary particles, and the time integration method
 # interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL49(),
+sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
             reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
             dtmax=1e-2, # Limit stepsize to prevent crashing
diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl
index 7f696b324..867911f5e 100644
--- a/examples/fluid/dam_break_2d.jl
+++ b/examples/fluid/dam_break_2d.jl
@@ -84,7 +84,7 @@ callbacks = CallbackSet(info_callback, saving_callback, density_reinit_cb)
 # Sometimes, the method fails to do so because forces become extremely large when
 # fluid particles are very close to boundary particles, and the time integration method
 # interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL49(),
+sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
             reltol=1e-5, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
             dtmax=1e-3, # Limit stepsize to prevent crashing
diff --git a/examples/fluid/dam_break_3d.jl b/examples/fluid/dam_break_3d.jl
index 8ee08f4b5..ad01dde42 100644
--- a/examples/fluid/dam_break_3d.jl
+++ b/examples/fluid/dam_break_3d.jl
@@ -70,7 +70,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
 # Sometimes, the method fails to do so because forces become extremely large when
 # fluid particles are very close to boundary particles, and the time integration method
 # interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL49(),
+sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
             reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
             dtmax=1e-2, # Limit stepsize to prevent crashing
diff --git a/examples/fluid/falling_water_column_2d.jl b/examples/fluid/falling_water_column_2d.jl
index 5346a70d5..8174ce117 100644
--- a/examples/fluid/falling_water_column_2d.jl
+++ b/examples/fluid/falling_water_column_2d.jl
@@ -71,7 +71,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
 # Sometimes, the method fails to do so because forces become extremely large when
 # fluid particles are very close to boundary particles, and the time integration method
 # interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL49(),
+sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
             reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
             dtmax=1e-2, # Limit stepsize to prevent crashing
diff --git a/examples/fluid/moving_wall_2d.jl b/examples/fluid/moving_wall_2d.jl
index 553c3e929..0720b76c0 100644
--- a/examples/fluid/moving_wall_2d.jl
+++ b/examples/fluid/moving_wall_2d.jl
@@ -92,7 +92,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
 # Sometimes, the method fails to do so because forces become extremely large when
 # fluid particles are very close to boundary particles, and the time integration method
 # interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL49(),
+sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
             reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
             dtmax=1e-2, # Limit stepsize to prevent crashing
diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl
index 65152f7f7..cd033717a 100644
--- a/examples/fluid/periodic_channel_2d.jl
+++ b/examples/fluid/periodic_channel_2d.jl
@@ -68,7 +68,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
 # Sometimes, the method fails to do so because forces become extremely large when
 # fluid particles are very close to boundary particles, and the time integration method
 # interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL49(),
+sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-8, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
             reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
             dtmax=1e-2, # Limit stepsize to prevent crashing
diff --git a/examples/fluid/rectangular_tank_2d.jl b/examples/fluid/rectangular_tank_2d.jl
index 37d81f1bb..35bc7acac 100644
--- a/examples/fluid/rectangular_tank_2d.jl
+++ b/examples/fluid/rectangular_tank_2d.jl
@@ -78,7 +78,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
 # Sometimes, the method fails to do so because forces become extremely large when
 # fluid particles are very close to boundary particles, and the time integration method
 # interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL49(),
+sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
             reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
             dtmax=1e-2, # Limit stepsize to prevent crashing
diff --git a/examples/fluid/rectangular_tank_edac_2d.jl b/examples/fluid/rectangular_tank_edac_2d.jl
index 9b1a63746..f2af73179 100644
--- a/examples/fluid/rectangular_tank_edac_2d.jl
+++ b/examples/fluid/rectangular_tank_edac_2d.jl
@@ -64,7 +64,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
 # Sometimes, the method fails to do so because forces become extremely large when
 # fluid particles are very close to boundary particles, and the time integration method
 # interprets this as an instability.
-sol = solve(ode, RDPK3SpFSAL49(),
+sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
             reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
             dtmax=1e-2, # Limit stepsize to prevent crashing

From 58fcf856ab089b9e0bea5dc08752337d5a110a47 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 7 Feb 2024 12:26:56 +0100
Subject: [PATCH 4/4] Use `trixi_include` from TrixiBase.jl (#375)

* Use `trixi_include` from TrixiBase.jl

* Add compat entry for TrixiBase.jl

* Add TrixiBase.jl API reference

* Add TrixiBase.jl API reference

* Add TrixiBase to docs dependencies

* Reformat

* Add `@test_nowarn_mod`

* Reformat

---------

Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>
---
 Project.toml                    |   2 +
 docs/Project.toml               |   5 +
 docs/make.jl                    |   7 +-
 docs/src/reference-trixibase.md |   9 ++
 src/TrixiParticles.jl           |   3 +-
 src/util.jl                     | 131 -------------------------
 test/examples/examples.jl       | 168 +++++++++++++++++---------------
 test/general/general.jl         |   2 -
 test/general/util.jl            |  89 -----------------
 test/test_util.jl               |  51 +++++++++-
 10 files changed, 160 insertions(+), 307 deletions(-)
 create mode 100644 docs/src/reference-trixibase.md
 delete mode 100644 test/general/util.jl

diff --git a/Project.toml b/Project.toml
index d80ab2193..7a9a127a5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,6 +19,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
 ThreadingUtilities = "8290d209-cae3-49c0-8002-c8c24d57dab5"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
@@ -34,4 +35,5 @@ StaticArrays = "1"
 StrideArrays = "0.1"
 ThreadingUtilities = "0.5"
 TimerOutputs = "0.5"
+TrixiBase = "0.1"
 WriteVTK = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 2ca447130..5644118db 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,3 +1,8 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 TrixiParticles = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
+
+[compat]
+Documenter = "1"
+TrixiBase = "0.1"
diff --git a/docs/make.jl b/docs/make.jl
index 1f43f51c6..838425328 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,4 +1,6 @@
-using Documenter, TrixiParticles
+using Documenter
+using TrixiParticles
+using TrixiBase
 
 # Get TrixiParticles.jl root directory
 trixiparticles_root_dir = dirname(@__DIR__)
@@ -67,8 +69,9 @@ makedocs(sitename="TrixiParticles.jl",
                                                                         "total_lagrangian_sph.md"),
                      "Boundary" => joinpath("systems", "boundary.md"),
                  ],
-                 "Time integration" => "time_integration.md",
+                 "Time Integration" => "time_integration.md",
                  "Callbacks" => "callbacks.md",
+                 "TrixiBase.jl API Reference" => "reference-trixibase.md",
              ],
              "Authors" => "authors.md",
              "Contributing" => "contributing.md",
diff --git a/docs/src/reference-trixibase.md b/docs/src/reference-trixibase.md
new file mode 100644
index 000000000..c7a970f88
--- /dev/null
+++ b/docs/src/reference-trixibase.md
@@ -0,0 +1,9 @@
+# TrixiBase.jl API
+
+```@meta
+CurrentModule = TrixiBase
+```
+
+```@autodocs
+Modules = [TrixiBase]
+```
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index f232a664b..ab0ade4f2 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -5,6 +5,7 @@ using Reexport: @reexport
 using Dates
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 using FastPow: @fastpow
+using ForwardDiff: ForwardDiff
 using LinearAlgebra: norm, dot, I, tr, inv, pinv, det
 using Morton: cartesian2morton
 using MuladdMacro: @muladd
@@ -17,8 +18,8 @@ using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
 using ThreadingUtilities
 using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!
+using TrixiBase: trixi_include
 using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save
-using ForwardDiff
 
 # util needs to be first because of macro @trixi_timeit
 include("util.jl")
diff --git a/src/util.jl b/src/util.jl
index c1e4f5bad..1a24ec23b 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -136,137 +136,6 @@ readdir(examples_dir())
 """
 examples_dir() = pkgdir(TrixiParticles, "examples")
 
-# Note: We can't call the method below `TrixiParticles.include` since that is created automatically
-# inside `module TrixiParticles` to `include` source files and evaluate them within the global scope
-# of `TrixiParticles`. However, users will want to evaluate in the global scope of `Main` or something
-# similar to manage dependencies on their own.
-"""
-    trixi_include([mod::Module=Main,] example::AbstractString; kwargs...)
-
-`include` the file `example` and evaluate its content in the global scope of module `mod`.
-You can override specific assignments in `example` by supplying keyword arguments.
-It's basic purpose is to make it easier to modify some parameters while running TrixiParticles from the
-REPL. Additionally, this is used in tests to reduce the computational burden for CI while still
-providing examples with sensible default values for users.
-
-Before replacing assignments in `example`, the keyword argument `maxiters` is inserted
-into calls to `solve` and `TrixiParticles.solve` with it's default value used in the SciML ecosystem
-for ODEs, see https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Miscellaneous.
-
-Copied from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl).
-
-# Examples
-
-```jldoctest
-julia> redirect_stdout(devnull) do
-           trixi_include(@__MODULE__, joinpath(examples_dir(), "dam_break_2d.jl"),
-                         tspan=(0.0, 0.1))
-           sol.t[end]
-       end
-0.1
-```
-"""
-function trixi_include(mod::Module, elixir::AbstractString; kwargs...)
-    # Check that all kwargs exist as assignments
-    code = read(elixir, String)
-    expr = Meta.parse("begin \n$code \nend")
-    expr = insert_maxiters(expr)
-
-    for (key, val) in kwargs
-        # This will throw an error when `key` is not found
-        find_assignment(expr, key)
-    end
-
-    # Include `elixir` and replace assignments
-    Base.include(ex -> replace_assignments(insert_maxiters(ex); kwargs...), mod, elixir)
-end
-
-trixi_include(elixir::AbstractString; kwargs...) = trixi_include(Main, elixir; kwargs...)
-
-# Helper methods used in the functions defined above, also copied from Trixi.jl
-
-# Apply the function `f` to `expr` and all sub-expressions recursively.
-walkexpr(f, expr::Expr) = f(Expr(expr.head, (walkexpr(f, arg) for arg in expr.args)...))
-walkexpr(f, x) = f(x)
-
-# Insert the keyword argument `maxiters` into calls to `solve` and `TrixiParticles.solve`
-# with default value `10^5` if it is not already present.
-function insert_maxiters(expr)
-    maxiters_default = 10^5
-
-    expr = walkexpr(expr) do x
-        if x isa Expr
-            is_plain_solve = x.head === Symbol("call") && x.args[1] === Symbol("solve")
-            is_trixi_solve = (x.head === Symbol("call") && x.args[1] isa Expr &&
-                              x.args[1].head === Symbol(".") &&
-                              x.args[1].args[1] === Symbol("TrixiParticles") &&
-                              x.args[1].args[2] isa QuoteNode &&
-                              x.args[1].args[2].value === Symbol("solve"))
-
-            if is_plain_solve || is_trixi_solve
-                # Do nothing if `maxiters` is already set as keyword argument...
-                for arg in x.args
-                    if arg isa Expr && arg.head === Symbol("kw") &&
-                       arg.args[1] === Symbol("maxiters")
-                        return x
-                    end
-                end
-
-                # ...and insert it otherwise.
-                push!(x.args, Expr(Symbol("kw"), Symbol("maxiters"), maxiters_default))
-            end
-        end
-        return x
-    end
-
-    return expr
-end
-
-# Replace assignments to `key` in `expr` by `key = val` for all `(key,val)` in `kwargs`.
-function replace_assignments(expr; kwargs...)
-    # replace explicit and keyword assignments
-    expr = walkexpr(expr) do x
-        if x isa Expr
-            for (key, val) in kwargs
-                if (x.head === Symbol("=") || x.head === :kw) && x.args[1] === Symbol(key)
-                    x.args[2] = :($val)
-                    # dump(x)
-                end
-            end
-        end
-        return x
-    end
-
-    return expr
-end
-
-# find a (keyword or common) assignment to `destination` in `expr`
-# and return the assigned value
-function find_assignment(expr, destination)
-    # declare result to be able to assign to it in the closure
-    local result
-    found = false
-
-    # find explicit and keyword assignments
-    walkexpr(expr) do x
-        if x isa Expr
-            if (x.head === Symbol("=") || x.head === :kw) &&
-               x.args[1] === Symbol(destination)
-                result = x.args[2]
-                found = true
-                # dump(x)
-            end
-        end
-        return x
-    end
-
-    if !found
-        throw(ArgumentError("assignment $destination not found in expression"))
-    end
-
-    return result
-end
-
 """
     @autoinfiltrate
     @autoinfiltrate condition::Bool
diff --git a/test/examples/examples.jl b/test/examples/examples.jl
index ac4f6218f..cf2301674 100644
--- a/test/examples/examples.jl
+++ b/test/examples/examples.jl
@@ -3,177 +3,185 @@
 @testset verbose=true "Examples" begin
     @testset verbose=true "Fluid" begin
         @trixi_testset "fluid/oscillating_drop_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid",
-                                                "oscillating_drop_2d.jl"))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "oscillating_drop_2d.jl"))
             @test sol.retcode == ReturnCode.Success
             # This is the error on an Apple M2 Pro. We need this tolerance to make CI pass.
             @test isapprox(error_A, 0.0001717690010767381, atol=1e-8)
         end
 
         @trixi_testset "fluid/rectangular_tank_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid",
-                                                "rectangular_tank_2d.jl"), tspan=(0.0, 0.1))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "rectangular_tank_2d.jl"),
+                                           tspan=(0.0, 0.1))
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fluid/rectangular_tank_2d.jl with SummationDensity" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid",
-                                                "rectangular_tank_2d.jl"), tspan=(0.0, 0.1),
-                                       fluid_density_calculator=SummationDensity(),
-                                       clip_negative_pressure=true)
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "rectangular_tank_2d.jl"),
+                                           tspan=(0.0, 0.1),
+                                           fluid_density_calculator=SummationDensity(),
+                                           clip_negative_pressure=true)
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fluid/rectangular_tank_3d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid",
-                                                "rectangular_tank_3d.jl"), tspan=(0.0, 0.1))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "rectangular_tank_3d.jl"),
+                                           tspan=(0.0, 0.1))
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fluid/rectangular_tank_3d.jl with SummationDensity" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid",
-                                                "rectangular_tank_3d.jl"), tspan=(0.0, 0.1),
-                                       fluid_density_calculator=SummationDensity(),
-                                       clip_negative_pressure=true)
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "rectangular_tank_3d.jl"),
+                                           tspan=(0.0, 0.1),
+                                           fluid_density_calculator=SummationDensity(),
+                                           clip_negative_pressure=true)
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fluid/rectangular_tank_edac_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid",
-                                                "rectangular_tank_edac_2d.jl"),
-                                       tspan=(0.0, 0.1))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "rectangular_tank_edac_2d.jl"),
+                                           tspan=(0.0, 0.1))
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fluid/dam_break_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
-                                       tspan=(0.0, 0.1))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "dam_break_2d.jl"),
+                                           tspan=(0.0, 0.1))
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fluid/dam_break_2d_corrections.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid",
-                                                "dam_break_2d_corrections.jl"),
-                                       tspan=(0.0, 0.1))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "dam_break_2d_corrections.jl"),
+                                           tspan=(0.0, 0.1))
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fluid/dam_break_3d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid", "dam_break_3d.jl"),
-                                       tspan=(0.0, 0.1), fluid_particle_spacing=0.1)
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "dam_break_3d.jl"),
+                                           tspan=(0.0, 0.1), fluid_particle_spacing=0.1)
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fluid/falling_water_column_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid",
-                                                "falling_water_column_2d.jl"),
-                                       tspan=(0.0, 0.4))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "falling_water_column_2d.jl"),
+                                           tspan=(0.0, 0.4))
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fluid/periodic_channel_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fluid",
-                                                "periodic_channel_2d.jl"),
-                                       tspan=(0.0, 0.4))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fluid",
+                                                    "periodic_channel_2d.jl"),
+                                           tspan=(0.0, 0.4))
             @test sol.retcode == ReturnCode.Success
         end
     end
 
     @testset verbose=true "Solid" begin
         @trixi_testset "solid/oscillating_beam_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "solid",
-                                                "oscillating_beam_2d.jl"), tspan=(0.0, 0.1))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "solid",
+                                                    "oscillating_beam_2d.jl"),
+                                           tspan=(0.0, 0.1))
             @test sol.retcode == ReturnCode.Success
         end
     end
 
     @testset verbose=true "FSI" begin
         @trixi_testset "fsi/falling_water_column_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fsi",
-                                                "falling_water_column_2d.jl"),
-                                       tspan=(0.0, 0.4))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fsi",
+                                                    "falling_water_column_2d.jl"),
+                                           tspan=(0.0, 0.4))
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fsi/dam_break_2d.jl" begin
             # Use rounded dimensions to avoid warnings
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fsi", "dam_break_2d.jl"),
-                                       initial_fluid_size=(0.15, 0.29),
-                                       tspan=(0.0, 0.4),
-                                       dtmax=1e-3)
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fsi",
+                                                    "dam_break_2d.jl"),
+                                           initial_fluid_size=(0.15, 0.29),
+                                           tspan=(0.0, 0.4),
+                                           dtmax=1e-3)
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fsi/dam_break_gate_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fsi",
-                                                "dam_break_gate_2d.jl"),
-                                       tspan=(0.0, 0.4),
-                                       dtmax=1e-3)
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fsi",
+                                                    "dam_break_gate_2d.jl"),
+                                           tspan=(0.0, 0.4),
+                                           dtmax=1e-3)
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "fsi/falling_spheres_2d.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "fsi",
-                                                "falling_spheres_2d.jl"),
-                                       tspan=(0.0, 1.0))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "fsi",
+                                                    "falling_spheres_2d.jl"),
+                                           tspan=(0.0, 1.0))
             @test sol.retcode == ReturnCode.Success
         end
     end
 
     @testset verbose=true "N-Body" begin
         @trixi_testset "n_body/n_body_solar_system.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "n_body",
-                                                "n_body_solar_system.jl"))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "n_body",
+                                                    "n_body_solar_system.jl"))
             @test sol.retcode == ReturnCode.Success
         end
 
         @trixi_testset "n_body/n_body_benchmark_trixi.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "n_body",
-                                                "n_body_benchmark_trixi.jl"))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "n_body",
+                                                    "n_body_benchmark_trixi.jl"))
         end
 
         @trixi_testset "n_body/n_body_benchmark_reference.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "n_body",
-                                                "n_body_benchmark_reference.jl"))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "n_body",
+                                                    "n_body_benchmark_reference.jl"))
         end
 
         @trixi_testset "n_body/n_body_benchmark_reference_faster.jl" begin
-            @test_nowarn trixi_include(joinpath(examples_dir(), "n_body",
-                                                "n_body_benchmark_reference_faster.jl"))
+            @test_nowarn_mod trixi_include(joinpath(examples_dir(), "n_body",
+                                                    "n_body_benchmark_reference_faster.jl"))
         end
     end
 
     @testset verbose=true "Postprocessing" begin
         @trixi_testset "postprocessing/interpolation_plane.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "postprocessing",
-                                                "interpolation_plane.jl"),
-                                       tspan=(0.0, 0.01))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "postprocessing",
+                                                    "interpolation_plane.jl"),
+                                           tspan=(0.0, 0.01))
             @test sol.retcode == ReturnCode.Success
         end
         @trixi_testset "postprocessing/interpolation_point_line.jl" begin
-            @test_nowarn trixi_include(@__MODULE__,
-                                       joinpath(examples_dir(), "postprocessing",
-                                                "interpolation_point_line.jl"))
+            @test_nowarn_mod trixi_include(@__MODULE__,
+                                           joinpath(examples_dir(), "postprocessing",
+                                                    "interpolation_point_line.jl"))
             @test sol.retcode == ReturnCode.Success
         end
     end
diff --git a/test/general/general.jl b/test/general/general.jl
index ac84de771..b59cff789 100644
--- a/test/general/general.jl
+++ b/test/general/general.jl
@@ -1,5 +1,3 @@
-# `util.jl` is added here to avoid confusion with `test_util.jl`
-include("util.jl")
 include("initial_condition.jl")
 include("smoothing_kernels.jl")
 include("density_calculator.jl")
diff --git a/test/general/util.jl b/test/general/util.jl
deleted file mode 100644
index 5da8700c2..000000000
--- a/test/general/util.jl
+++ /dev/null
@@ -1,89 +0,0 @@
-@testset verbose=true "trixi_include" begin
-    @trixi_testset "Basic" begin
-        example = """
-            x = 4
-            """
-
-        filename = tempname()
-        try
-            open(filename, "w") do file
-                write(file, example)
-            end
-
-            # Use `@trixi_testset`, which wraps code in a temporary module, and call
-            # `trixi_include` with `@__MODULE__` in order to isolate this test.
-            @test_nowarn trixi_include(@__MODULE__, filename)
-            @test @isdefined x
-            @test x == 4
-
-            @test_nowarn trixi_include(@__MODULE__, filename, x=7)
-            @test x == 7
-
-            @test_throws "assignment y not found in expression" trixi_include(@__MODULE__,
-                                                                              filename, y=3)
-        finally
-            rm(filename, force=true)
-        end
-    end
-
-    @trixi_testset "With `solve` Without `maxiters`" begin
-        example = """
-            solve() = 0
-            x = solve()
-            """
-
-        filename = tempname()
-        try
-            open(filename, "w") do file
-                write(file, example)
-            end
-
-            # Use `@trixi_testset`, which wraps code in a temporary module, and call
-            # `trixi_include` with `@__MODULE__` in order to isolate this test.
-            @test_throws "no method matching solve(; maxiters::Int64)" trixi_include(@__MODULE__,
-                                                                                     filename)
-
-            @test_throws "no method matching solve(; maxiters::Int64)" trixi_include(@__MODULE__,
-                                                                                     filename,
-                                                                                     maxiters=3)
-        finally
-            rm(filename, force=true)
-        end
-    end
-
-    @trixi_testset "With `solve` with `maxiters`" begin
-        # We need another example file that we include with `Base.include` first, in order to
-        # define the `solve` method without `trixi_include` trying to insert `maxiters` kwargs.
-        example1 = """
-            solve(; maxiters=0) = maxiters
-            """
-
-        example2 = """
-            x = solve()
-            """
-
-        filename1 = tempname()
-        filename2 = tempname()
-        try
-            open(filename1, "w") do file
-                write(file, example1)
-            end
-            open(filename2, "w") do file
-                write(file, example2)
-            end
-
-            # Use `@trixi_testset`, which wraps code in a temporary module, and call
-            # `Base.include` and `trixi_include` with `@__MODULE__` in order to isolate this test.
-            Base.include(@__MODULE__, filename1)
-            @test_nowarn trixi_include(@__MODULE__, filename2)
-            @test @isdefined x
-            @test x == 10^5
-
-            @test_nowarn trixi_include(@__MODULE__, filename2, maxiters=7)
-            @test x == 7
-        finally
-            rm(filename1, force=true)
-            rm(filename2, force=true)
-        end
-    end
-end
diff --git a/test/test_util.jl b/test/test_util.jl
index 68d6b48d9..55a2edda8 100644
--- a/test/test_util.jl
+++ b/test/test_util.jl
@@ -21,12 +21,16 @@ macro trixi_testset(name, expr)
         println("═"^100)
         println($name)
 
-        local time_start = time_ns()
-
         @eval module $mod
         using Test
         using TrixiParticles
 
+        # We also include this file again to provide the definition of
+        # the other testing macros. This allows to use `@trixi_testset`
+        # in a nested fashion and also call `@test_nowarn_mod` from
+        # there.
+        include(@__FILE__)
+
         @testset verbose=true $name $expr
         end
 
@@ -34,6 +38,49 @@ macro trixi_testset(name, expr)
     end
 end
 
+# Copied from TrixiBase.jl. See https://github.com/trixi-framework/TrixiBase.jl/issues/9.
+"""
+    @test_nowarn_mod expr
+
+Modified version of `@test_nowarn expr` that prints the content of `stderr` when
+it is not empty and ignores some common info statements printed in Trixi.jl
+uses.
+"""
+macro test_nowarn_mod(expr, additional_ignore_content=String[])
+    quote
+        let fname = tempname()
+            try
+                ret = open(fname, "w") do f
+                    redirect_stderr(f) do
+                        $(esc(expr))
+                    end
+                end
+                stderr_content = read(fname, String)
+                if !isempty(stderr_content)
+                    println("Content of `stderr`:\n", stderr_content)
+                end
+
+                # Patterns matching the following ones will be ignored. Additional patterns
+                # passed as arguments can also be regular expressions, so we just use the
+                # type `Any` for `ignore_content`.
+                ignore_content = Any["[ Info: You just called `trixi_include`. Julia may now compile the code, please be patient.\n"]
+                append!(ignore_content, $additional_ignore_content)
+                for pattern in ignore_content
+                    stderr_content = replace(stderr_content, pattern => "")
+                end
+
+                # We also ignore simple module redefinitions for convenience. Thus, we
+                # check whether every line of `stderr_content` is of the form of a
+                # module replacement warning.
+                @test occursin(r"^(WARNING: replacing module .+\.\n)*$", stderr_content)
+                ret
+            finally
+                rm(fname, force=true)
+            end
+        end
+    end
+end
+
 function perturb!(data, amplitude)
     for i in eachindex(data)
         # Perturbation in the interval (-amplitude, amplitude)