From b43ed2d139f4ccdadc4f38698009b618b46bacba Mon Sep 17 00:00:00 2001
From: patrickersing <patrick.ersing@liu.se>
Date: Wed, 13 Mar 2024 08:15:42 +0100
Subject: [PATCH 1/9] add TimeSeriesCallback support for TreeMesh1D

---
 .../elixir_euler_source_terms.jl              |  5 ++-
 src/callbacks_step/time_series_dg2d.jl        | 32 ++++++++++++++++++-
 test/test_tree_1d_euler.jl                    |  5 +++
 3 files changed, 40 insertions(+), 2 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_source_terms.jl
index 555910f69f0..bf6bf7488a9 100644
--- a/examples/tree_1d_dgsem/elixir_euler_source_terms.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_source_terms.jl
@@ -44,10 +44,13 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 stepsize_callback = StepsizeCallback(cfl = 0.8)
 
+time_series = TimeSeriesCallback(semi, [0.0, 0.33, 1.0], interval = 10)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
                         save_solution,
-                        stepsize_callback)
+                        stepsize_callback,
+                        time_series)
 
 ###############################################################################
 # run the simulation
diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl
index c15945d6e16..4da8daf403c 100644
--- a/src/callbacks_step/time_series_dg2d.jl
+++ b/src/callbacks_step/time_series_dg2d.jl
@@ -6,7 +6,7 @@
 #! format: noindent
 
 # Creates cache for time series callback
-function create_cache_time_series(point_coordinates, mesh::TreeMesh{2}, dg, cache)
+function create_cache_time_series(point_coordinates, mesh::TreeMesh, dg, cache)
     # Determine element ids for point coordinates
     element_ids = get_elements_by_coordinates(point_coordinates, mesh, dg, cache)
 
@@ -119,6 +119,36 @@ function calc_interpolating_polynomials(coordinates, element_ids, mesh::TreeMesh
 end
 
 # Record the solution variables at each given point
+function record_state_at_points!(point_data, u, solution_variables,
+                                 n_solution_variables,
+                                 mesh::TreeMesh{1}, equations, dg::DG,
+                                 time_series_cache)
+    @unpack element_ids, interpolating_polynomials = time_series_cache
+    old_length = length(first(point_data))
+    new_length = old_length + n_solution_variables
+
+    # Loop over all points/elements that should be recorded
+    for index in 1:length(element_ids)
+        # Extract data array and element id
+        data = point_data[index]
+        element_id = element_ids[index]
+
+        # Make room for new data to be recorded
+        resize!(data, new_length)
+        data[(old_length + 1):new_length] .= zero(eltype(data))
+
+        # Loop over all nodes to compute their contribution to the interpolated values
+        for i in eachnode(dg)
+            u_node = solution_variables(get_node_vars(u, equations, dg, i,
+                                                      element_id), equations)
+
+            for v in 1:length(u_node)
+                data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index])
+            end
+        end
+    end
+end
+
 function record_state_at_points!(point_data, u, solution_variables,
                                  n_solution_variables,
                                  mesh::TreeMesh{2}, equations, dg::DG,
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index f26500b411c..2d36ef3bc17 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -30,6 +30,11 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
         du_ode = similar(u_ode)
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
+    # Extra test to make sure the "TimeSeriesCallback" made correct data.
+    # Extracts data at all points from the first step of the time series and compares it to the exact solution.
+    point_data = [getindex(time_series.affect!.point_data[i], 1:3) for i=1:3]
+    ref_data = [initial_condition_convergence_test(time_series.affect!.point_coordinates[i], time_series.affect!.time[1], equations) for i=1:3]
+    @test point_data ≈ ref_data atol=1e-6
 end
 
 @trixi_testset "elixir_euler_convergence_pure_fv.jl" begin

From 1c347a622a4d14d45a826d3e97fcdf4a8251dda4 Mon Sep 17 00:00:00 2001
From: patrickersing <patrick.ersing@liu.se>
Date: Wed, 13 Mar 2024 09:45:03 +0100
Subject: [PATCH 2/9] add TimeSeriesCallback support for TreeMesh3D

---
 .../elixir_euler_source_terms.jl              |  7 +++-
 src/callbacks_step/time_series_dg2d.jl        | 36 ++++++++++++++++++-
 2 files changed, 41 insertions(+), 2 deletions(-)

diff --git a/examples/tree_3d_dgsem/elixir_euler_source_terms.jl b/examples/tree_3d_dgsem/elixir_euler_source_terms.jl
index f0246c30490..2d57ac279bd 100644
--- a/examples/tree_3d_dgsem/elixir_euler_source_terms.jl
+++ b/examples/tree_3d_dgsem/elixir_euler_source_terms.jl
@@ -41,10 +41,15 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 stepsize_callback = StepsizeCallback(cfl = 0.6)
 
+time_series = TimeSeriesCallback(semi,
+                                 [(0.0, 0.0, 0.0), (0.33, 0.33, 0.33), (1.0, 1.0, 1.0)],
+                                 interval = 100)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
                         save_solution,
-                        stepsize_callback)
+                        stepsize_callback,
+                        time_series)
 
 ###############################################################################
 # run the simulation
diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl
index 4da8daf403c..d85b32b8489 100644
--- a/src/callbacks_step/time_series_dg2d.jl
+++ b/src/callbacks_step/time_series_dg2d.jl
@@ -143,7 +143,8 @@ function record_state_at_points!(point_data, u, solution_variables,
                                                       element_id), equations)
 
             for v in 1:length(u_node)
-                data[old_length + v] += (u_node[v] * interpolating_polynomials[i, 1, index])
+                data[old_length + v] += (u_node[v] *
+                                         interpolating_polynomials[i, 1, index])
             end
         end
     end
@@ -180,4 +181,37 @@ function record_state_at_points!(point_data, u, solution_variables,
         end
     end
 end
+
+function record_state_at_points!(point_data, u, solution_variables,
+                                 n_solution_variables,
+                                 mesh::TreeMesh{3}, equations, dg::DG,
+                                 time_series_cache)
+    @unpack element_ids, interpolating_polynomials = time_series_cache
+    old_length = length(first(point_data))
+    new_length = old_length + n_solution_variables
+
+    # Loop over all points/elements that should be recorded
+    for index in 1:length(element_ids)
+        # Extract data array and element id
+        data = point_data[index]
+        element_id = element_ids[index]
+
+        # Make room for new data to be recorded
+        resize!(data, new_length)
+        data[(old_length + 1):new_length] .= zero(eltype(data))
+
+        # Loop over all nodes to compute their contribution to the interpolated values
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_node = solution_variables(get_node_vars(u, equations, dg, i, j, k,
+                                                      element_id), equations)
+
+            for v in 1:length(u_node)
+                data[old_length + v] += (u_node[v]
+                                         * interpolating_polynomials[i, 1, index]
+                                         * interpolating_polynomials[j, 2, index]
+                                         * interpolating_polynomials[k, 3, index])
+            end
+        end
+    end
+end
 end # @muladd

From 91e1287e0d0e0ff76ba294748fb90e0009265ed7 Mon Sep 17 00:00:00 2001
From: patrickersing <patrick.ersing@liu.se>
Date: Wed, 13 Mar 2024 09:48:46 +0100
Subject: [PATCH 3/9] add testing

---
 test/test_tree_1d_euler.jl | 15 +++++++++++----
 test/test_tree_3d_euler.jl | 32 ++++++++++++++++++++++++++++++++
 2 files changed, 43 insertions(+), 4 deletions(-)

diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index 2d36ef3bc17..5b34f077d98 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -31,10 +31,17 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
     # Extra test to make sure the "TimeSeriesCallback" made correct data.
-    # Extracts data at all points from the first step of the time series and compares it to the exact solution.
-    point_data = [getindex(time_series.affect!.point_data[i], 1:3) for i=1:3]
-    ref_data = [initial_condition_convergence_test(time_series.affect!.point_coordinates[i], time_series.affect!.time[1], equations) for i=1:3]
-    @test point_data ≈ ref_data atol=1e-6
+    # Extracts data at all points from the first step of the time series and compares it to the 
+    # exact solution and an interpolated reference solution
+    point_data = [getindex(time_series.affect!.point_data[i], 1:3) for i in 1:3]
+    exact_data = [initial_condition_convergence_test(time_series.affect!.point_coordinates[i],
+                                                     time_series.affect!.time[1],
+                                                     equations) for i in 1:3]
+    ref_data = [[1.968279088772251, 1.9682791565395945, 3.874122958278797],
+        [2.0654816955822017, 2.0654817326611883, 4.26621471136323],
+        [2.0317209235018936, 2.0317209516429506, 4.127889808862571]]
+    @test point_data≈exact_data atol=1e-6
+    @test point_data ≈ ref_data
 end
 
 @trixi_testset "elixir_euler_convergence_pure_fv.jl" begin
diff --git a/test/test_tree_3d_euler.jl b/test/test_tree_3d_euler.jl
index e9e2b82fec5..3b59bbe2c9f 100644
--- a/test/test_tree_3d_euler.jl
+++ b/test/test_tree_3d_euler.jl
@@ -34,6 +34,38 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_dgsem")
         du_ode = similar(u_ode)
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
+    # Extra test to make sure the "TimeSeriesCallback" made correct data.
+    # Extracts data at all points from the first step of the time series and compares it to the 
+    # exact solution and an interpolated reference solution
+    point_data = [getindex(time_series.affect!.point_data[i], 1:5) for i in 1:3]
+    exact_data = [initial_condition_convergence_test(time_series.affect!.point_coordinates[:,
+                                                                                           i],
+                                                     time_series.affect!.time[1],
+                                                     equations) for i in 1:3]
+    ref_data = [
+        [
+            1.9053441545123802,
+            1.9053448755267233,
+            1.905344875526724,
+            1.9053448755267266,
+            3.6303259410202022,
+        ],
+        [
+            2.0956374376888744,
+            2.095637084878075,
+            2.0956370848780743,
+            2.095637084878072,
+            4.391698734498684,
+        ],
+        [
+            2.0946764789674774,
+            2.0946781952809355,
+            2.094678195280935,
+            2.0946781952809324,
+            4.3876634952412505,
+        ]]
+    @test point_data≈exact_data atol=1e-1
+    @test point_data ≈ ref_data
 end
 
 @trixi_testset "elixir_euler_convergence_pure_fv.jl" begin

From 1305daa2798b42158786bf72853b79b56cb505b1 Mon Sep 17 00:00:00 2001
From: patrickersing <patrick.ersing@liu.se>
Date: Fri, 15 Mar 2024 10:05:24 +0100
Subject: [PATCH 4/9] update tests

---
 .../elixir_euler_source_terms.jl              |  2 +-
 test/test_tree_1d_euler.jl                    |  5 ++-
 test/test_tree_3d_euler.jl                    | 31 +++++--------------
 3 files changed, 13 insertions(+), 25 deletions(-)

diff --git a/examples/tree_3d_dgsem/elixir_euler_source_terms.jl b/examples/tree_3d_dgsem/elixir_euler_source_terms.jl
index 2d57ac279bd..0cf615fda10 100644
--- a/examples/tree_3d_dgsem/elixir_euler_source_terms.jl
+++ b/examples/tree_3d_dgsem/elixir_euler_source_terms.jl
@@ -43,7 +43,7 @@ stepsize_callback = StepsizeCallback(cfl = 0.6)
 
 time_series = TimeSeriesCallback(semi,
                                  [(0.0, 0.0, 0.0), (0.33, 0.33, 0.33), (1.0, 1.0, 1.0)],
-                                 interval = 100)
+                                 interval = 10)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index 5b34f077d98..784d123128e 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -21,7 +21,10 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
                             1.6205433861493646e-7,
                             1.465427772462391e-7,
                             5.372255111879554e-7,
-                        ])
+                        ],
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time series to check against.
+                        coverage_override=(maxiters = 20,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
diff --git a/test/test_tree_3d_euler.jl b/test/test_tree_3d_euler.jl
index 3b59bbe2c9f..98549f2a5ab 100644
--- a/test/test_tree_3d_euler.jl
+++ b/test/test_tree_3d_euler.jl
@@ -25,7 +25,10 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_dgsem")
                             0.032179231640894645,
                             0.032179231640895534,
                             0.0655408023333299,
-                        ])
+                        ],
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time series to check against.
+                        coverage_override=(maxiters = 20,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -43,28 +46,10 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_dgsem")
                                                      time_series.affect!.time[1],
                                                      equations) for i in 1:3]
     ref_data = [
-        [
-            1.9053441545123802,
-            1.9053448755267233,
-            1.905344875526724,
-            1.9053448755267266,
-            3.6303259410202022,
-        ],
-        [
-            2.0956374376888744,
-            2.095637084878075,
-            2.0956370848780743,
-            2.095637084878072,
-            4.391698734498684,
-        ],
-        [
-            2.0946764789674774,
-            2.0946781952809355,
-            2.094678195280935,
-            2.0946781952809324,
-            4.3876634952412505,
-        ]]
-    @test point_data≈exact_data atol=1e-1
+        [2.0908436034061966, 2.0907469131413534, 2.090746913141344, 2.0907469131413516, 4.376022861751483],
+        [1.9121470846408657, 1.911738389880639, 1.9117383898806395, 1.9117383898806386, 3.639369160978031],
+        [1.8960587498057229, 1.896520160078096, 1.8965201600780932, 1.896520160078095, 3.614473843601274]]
+    @test point_data ≈ exact_data atol=1e-1
     @test point_data ≈ ref_data
 end
 

From ec8b3effc2f7e858f5ab0eae4b7da71953726bbf Mon Sep 17 00:00:00 2001
From: patrickersing <patrick.ersing@liu.se>
Date: Fri, 15 Mar 2024 10:11:50 +0100
Subject: [PATCH 5/9] reorganize files structure for TimeSeriesCallback

---
 src/callbacks_step/time_series.jl         |   1 +
 src/callbacks_step/time_series_dg.jl      |  37 +++++
 src/callbacks_step/time_series_dg2d.jl    | 183 +---------------------
 src/callbacks_step/time_series_dg_tree.jl | 183 ++++++++++++++++++++++
 test/test_tree_3d_euler.jl                |  26 ++-
 5 files changed, 244 insertions(+), 186 deletions(-)
 create mode 100644 src/callbacks_step/time_series_dg_tree.jl

diff --git a/src/callbacks_step/time_series.jl b/src/callbacks_step/time_series.jl
index f6d76f0fb15..cf90a47900c 100644
--- a/src/callbacks_step/time_series.jl
+++ b/src/callbacks_step/time_series.jl
@@ -218,5 +218,6 @@ function (time_series_callback::TimeSeriesCallback)(integrator)
 end
 
 include("time_series_dg.jl")
+include("time_series_dg_tree.jl")
 include("time_series_dg2d.jl")
 end # @muladd
diff --git a/src/callbacks_step/time_series_dg.jl b/src/callbacks_step/time_series_dg.jl
index ae394afbbfd..3781a10662d 100644
--- a/src/callbacks_step/time_series_dg.jl
+++ b/src/callbacks_step/time_series_dg.jl
@@ -34,4 +34,41 @@ function save_time_series_file(time_series_callback,
         end
     end
 end
+
+# Creates cache for time series callback
+function create_cache_time_series(point_coordinates,
+                                  mesh::Union{TreeMesh, UnstructuredMesh2D},
+                                  dg, cache)
+    # Determine element ids for point coordinates
+    element_ids = get_elements_by_coordinates(point_coordinates, mesh, dg, cache)
+
+    # Calculate & store Lagrange interpolation polynomials
+    interpolating_polynomials = calc_interpolating_polynomials(point_coordinates,
+                                                               element_ids, mesh,
+                                                               dg, cache)
+
+    time_series_cache = (; element_ids, interpolating_polynomials)
+
+    return time_series_cache
+end
+
+function get_elements_by_coordinates(coordinates, mesh, dg, cache)
+    element_ids = Vector{Int}(undef, size(coordinates, 2))
+    get_elements_by_coordinates!(element_ids, coordinates, mesh, dg, cache)
+
+    return element_ids
+end
+
+function calc_interpolating_polynomials(coordinates, element_ids,
+                                        mesh::Union{TreeMesh, UnstructuredMesh2D},
+                                        dg, cache)
+    interpolating_polynomials = Array{real(dg), 3}(undef,
+                                                   nnodes(dg), ndims(mesh),
+                                                   length(element_ids))
+    calc_interpolating_polynomials!(interpolating_polynomials, coordinates, element_ids,
+                                    mesh, dg,
+                                    cache)
+
+    return interpolating_polynomials
+end
 end # @muladd
diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg2d.jl
index b8c1f6b0418..f6d1bb48f24 100644
--- a/src/callbacks_step/time_series_dg2d.jl
+++ b/src/callbacks_step/time_series_dg2d.jl
@@ -5,71 +5,6 @@
 @muladd begin
 #! format: noindent
 
-# Creates cache for time series callback
-function create_cache_time_series(point_coordinates,
-                                  mesh::Union{TreeMesh, UnstructuredMesh2D},
-                                  dg, cache)
-    # Determine element ids for point coordinates
-    element_ids = get_elements_by_coordinates(point_coordinates, mesh, dg, cache)
-
-    # Calculate & store Lagrange interpolation polynomials
-    interpolating_polynomials = calc_interpolating_polynomials(point_coordinates,
-                                                               element_ids, mesh,
-                                                               dg, cache)
-
-    time_series_cache = (; element_ids, interpolating_polynomials)
-
-    return time_series_cache
-end
-
-# Find element ids containing coordinates given as a matrix [ndims, npoints]
-function get_elements_by_coordinates!(element_ids, coordinates, mesh::TreeMesh, dg,
-                                      cache)
-    if length(element_ids) != size(coordinates, 2)
-        throw(DimensionMismatch("storage length for element ids does not match the number of coordinates"))
-    end
-
-    @unpack cell_ids = cache.elements
-    @unpack tree = mesh
-
-    # Reset element ids - 0 indicates "not (yet) found"
-    element_ids .= 0
-    found_elements = 0
-
-    # Iterate over all elements
-    for element in eachelement(dg, cache)
-        # Get cell id
-        cell_id = cell_ids[element]
-
-        # Iterate over coordinates
-        for index in 1:length(element_ids)
-            # Skip coordinates for which an element has already been found
-            if element_ids[index] > 0
-                continue
-            end
-
-            # Construct point
-            x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh)))
-
-            # Skip if point is not in cell
-            if !is_point_in_cell(tree, x, cell_id)
-                continue
-            end
-
-            # Otherwise point is in cell and thus in element
-            element_ids[index] = element
-            found_elements += 1
-        end
-
-        # Exit loop if all elements have already been found
-        if found_elements == length(element_ids)
-            break
-        end
-    end
-
-    return element_ids
-end
-
 # Elements on an `UnstructuredMesh2D` are possibly curved. Assume that each
 # element is convex, i.e., all interior angles are less than 180 degrees.
 # This routine computes the shortest distance from a given point to each element
@@ -208,44 +143,6 @@ function calc_minimum_surface_distance(point, node_coordinates,
     return min_distance2, indices
 end
 
-function get_elements_by_coordinates(coordinates, mesh, dg, cache)
-    element_ids = Vector{Int}(undef, size(coordinates, 2))
-    get_elements_by_coordinates!(element_ids, coordinates, mesh, dg, cache)
-
-    return element_ids
-end
-
-# Calculate the interpolating polynomials to extract data at the given coordinates
-# The coordinates are known to be located in the respective element in `element_ids`
-function calc_interpolating_polynomials!(interpolating_polynomials, coordinates,
-                                         element_ids,
-                                         mesh::TreeMesh, dg::DGSEM, cache)
-    @unpack tree = mesh
-    @unpack nodes = dg.basis
-
-    wbary = barycentric_weights(nodes)
-
-    for index in 1:length(element_ids)
-        # Construct point
-        x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh)))
-
-        # Convert to unit coordinates
-        cell_id = cache.elements.cell_ids[element_ids[index]]
-        cell_coordinates_ = cell_coordinates(tree, cell_id)
-        cell_length = length_at_cell(tree, cell_id)
-        unit_coordinates = (x .- cell_coordinates_) * 2 / cell_length
-
-        # Calculate interpolating polynomial for each dimension, making use of tensor product structure
-        for d in 1:ndims(mesh)
-            interpolating_polynomials[:, d, index] .= lagrange_interpolating_polynomials(unit_coordinates[d],
-                                                                                         nodes,
-                                                                                         wbary)
-        end
-    end
-
-    return interpolating_polynomials
-end
-
 function calc_interpolating_polynomials!(interpolating_polynomials, coordinates,
                                          element_ids,
                                          mesh::UnstructuredMesh2D, dg::DGSEM, cache)
@@ -374,54 +271,9 @@ end
     return SVector(xi, eta)
 end
 
-function calc_interpolating_polynomials(coordinates, element_ids,
-                                        mesh::Union{TreeMesh, UnstructuredMesh2D},
-                                        dg, cache)
-    interpolating_polynomials = Array{real(dg), 3}(undef,
-                                                   nnodes(dg), ndims(mesh),
-                                                   length(element_ids))
-    calc_interpolating_polynomials!(interpolating_polynomials, coordinates, element_ids,
-                                    mesh, dg,
-                                    cache)
-
-    return interpolating_polynomials
-end
-
-# Record the solution variables at each given point
 function record_state_at_points!(point_data, u, solution_variables,
                                  n_solution_variables,
-                                 mesh::TreeMesh{1}, equations, dg::DG,
-                                 time_series_cache)
-    @unpack element_ids, interpolating_polynomials = time_series_cache
-    old_length = length(first(point_data))
-    new_length = old_length + n_solution_variables
-
-    # Loop over all points/elements that should be recorded
-    for index in 1:length(element_ids)
-        # Extract data array and element id
-        data = point_data[index]
-        element_id = element_ids[index]
-
-        # Make room for new data to be recorded
-        resize!(data, new_length)
-        data[(old_length + 1):new_length] .= zero(eltype(data))
-
-        # Loop over all nodes to compute their contribution to the interpolated values
-        for i in eachnode(dg)
-            u_node = solution_variables(get_node_vars(u, equations, dg, i,
-                                                      element_id), equations)
-
-            for v in 1:length(u_node)
-                data[old_length + v] += (u_node[v] *
-                                         interpolating_polynomials[i, 1, index])
-            end
-        end
-    end
-end
-
-function record_state_at_points!(point_data, u, solution_variables,
-                                 n_solution_variables,
-                                 mesh::Union{TreeMesh{2}, UnstructuredMesh2D},
+                                 mesh::UnstructuredMesh2D,
                                  equations, dg::DG, time_series_cache)
     @unpack element_ids, interpolating_polynomials = time_series_cache
     old_length = length(first(point_data))
@@ -450,37 +302,4 @@ function record_state_at_points!(point_data, u, solution_variables,
         end
     end
 end
-
-function record_state_at_points!(point_data, u, solution_variables,
-                                 n_solution_variables,
-                                 mesh::TreeMesh{3}, equations, dg::DG,
-                                 time_series_cache)
-    @unpack element_ids, interpolating_polynomials = time_series_cache
-    old_length = length(first(point_data))
-    new_length = old_length + n_solution_variables
-
-    # Loop over all points/elements that should be recorded
-    for index in 1:length(element_ids)
-        # Extract data array and element id
-        data = point_data[index]
-        element_id = element_ids[index]
-
-        # Make room for new data to be recorded
-        resize!(data, new_length)
-        data[(old_length + 1):new_length] .= zero(eltype(data))
-
-        # Loop over all nodes to compute their contribution to the interpolated values
-        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
-            u_node = solution_variables(get_node_vars(u, equations, dg, i, j, k,
-                                                      element_id), equations)
-
-            for v in 1:length(u_node)
-                data[old_length + v] += (u_node[v]
-                                         * interpolating_polynomials[i, 1, index]
-                                         * interpolating_polynomials[j, 2, index]
-                                         * interpolating_polynomials[k, 3, index])
-            end
-        end
-    end
-end
 end # @muladd
diff --git a/src/callbacks_step/time_series_dg_tree.jl b/src/callbacks_step/time_series_dg_tree.jl
new file mode 100644
index 00000000000..f6e844021e2
--- /dev/null
+++ b/src/callbacks_step/time_series_dg_tree.jl
@@ -0,0 +1,183 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+# Find element ids containing coordinates given as a matrix [ndims, npoints]
+function get_elements_by_coordinates!(element_ids, coordinates, mesh::TreeMesh, dg,
+                                      cache)
+    if length(element_ids) != size(coordinates, 2)
+        throw(DimensionMismatch("storage length for element ids does not match the number of coordinates"))
+    end
+
+    @unpack cell_ids = cache.elements
+    @unpack tree = mesh
+
+    # Reset element ids - 0 indicates "not (yet) found"
+    element_ids .= 0
+    found_elements = 0
+
+    # Iterate over all elements
+    for element in eachelement(dg, cache)
+        # Get cell id
+        cell_id = cell_ids[element]
+
+        # Iterate over coordinates
+        for index in 1:length(element_ids)
+            # Skip coordinates for which an element has already been found
+            if element_ids[index] > 0
+                continue
+            end
+
+            # Construct point
+            x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh)))
+
+            # Skip if point is not in cell
+            if !is_point_in_cell(tree, x, cell_id)
+                continue
+            end
+
+            # Otherwise point is in cell and thus in element
+            element_ids[index] = element
+            found_elements += 1
+        end
+
+        # Exit loop if all elements have already been found
+        if found_elements == length(element_ids)
+            break
+        end
+    end
+
+    return element_ids
+end
+
+# Calculate the interpolating polynomials to extract data at the given coordinates
+# The coordinates are known to be located in the respective element in `element_ids`
+function calc_interpolating_polynomials!(interpolating_polynomials, coordinates,
+                                         element_ids,
+                                         mesh::TreeMesh, dg::DGSEM, cache)
+    @unpack tree = mesh
+    @unpack nodes = dg.basis
+
+    wbary = barycentric_weights(nodes)
+
+    for index in 1:length(element_ids)
+        # Construct point
+        x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh)))
+
+        # Convert to unit coordinates
+        cell_id = cache.elements.cell_ids[element_ids[index]]
+        cell_coordinates_ = cell_coordinates(tree, cell_id)
+        cell_length = length_at_cell(tree, cell_id)
+        unit_coordinates = (x .- cell_coordinates_) * 2 / cell_length
+
+        # Calculate interpolating polynomial for each dimension, making use of tensor product structure
+        for d in 1:ndims(mesh)
+            interpolating_polynomials[:, d, index] .= lagrange_interpolating_polynomials(unit_coordinates[d],
+                                                                                         nodes,
+                                                                                         wbary)
+        end
+    end
+
+    return interpolating_polynomials
+end
+
+# Record the solution variables at each given point
+function record_state_at_points!(point_data, u, solution_variables,
+                                 n_solution_variables,
+                                 mesh::TreeMesh{1}, equations, dg::DG,
+                                 time_series_cache)
+    @unpack element_ids, interpolating_polynomials = time_series_cache
+    old_length = length(first(point_data))
+    new_length = old_length + n_solution_variables
+
+    # Loop over all points/elements that should be recorded
+    for index in 1:length(element_ids)
+        # Extract data array and element id
+        data = point_data[index]
+        element_id = element_ids[index]
+
+        # Make room for new data to be recorded
+        resize!(data, new_length)
+        data[(old_length + 1):new_length] .= zero(eltype(data))
+
+        # Loop over all nodes to compute their contribution to the interpolated values
+        for i in eachnode(dg)
+            u_node = solution_variables(get_node_vars(u, equations, dg, i,
+                                                      element_id), equations)
+
+            for v in 1:length(u_node)
+                data[old_length + v] += (u_node[v] *
+                                         interpolating_polynomials[i, 1, index])
+            end
+        end
+    end
+end
+
+function record_state_at_points!(point_data, u, solution_variables,
+                                 n_solution_variables,
+                                 mesh::TreeMesh{2},
+                                 equations, dg::DG, time_series_cache)
+    @unpack element_ids, interpolating_polynomials = time_series_cache
+    old_length = length(first(point_data))
+    new_length = old_length + n_solution_variables
+
+    # Loop over all points/elements that should be recorded
+    for index in 1:length(element_ids)
+        # Extract data array and element id
+        data = point_data[index]
+        element_id = element_ids[index]
+
+        # Make room for new data to be recorded
+        resize!(data, new_length)
+        data[(old_length + 1):new_length] .= zero(eltype(data))
+
+        # Loop over all nodes to compute their contribution to the interpolated values
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = solution_variables(get_node_vars(u, equations, dg, i, j,
+                                                      element_id), equations)
+
+            for v in 1:length(u_node)
+                data[old_length + v] += (u_node[v]
+                                         * interpolating_polynomials[i, 1, index]
+                                         * interpolating_polynomials[j, 2, index])
+            end
+        end
+    end
+end
+
+function record_state_at_points!(point_data, u, solution_variables,
+                                 n_solution_variables,
+                                 mesh::TreeMesh{3}, equations, dg::DG,
+                                 time_series_cache)
+    @unpack element_ids, interpolating_polynomials = time_series_cache
+    old_length = length(first(point_data))
+    new_length = old_length + n_solution_variables
+
+    # Loop over all points/elements that should be recorded
+    for index in 1:length(element_ids)
+        # Extract data array and element id
+        data = point_data[index]
+        element_id = element_ids[index]
+
+        # Make room for new data to be recorded
+        resize!(data, new_length)
+        data[(old_length + 1):new_length] .= zero(eltype(data))
+
+        # Loop over all nodes to compute their contribution to the interpolated values
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_node = solution_variables(get_node_vars(u, equations, dg, i, j, k,
+                                                      element_id), equations)
+
+            for v in 1:length(u_node)
+                data[old_length + v] += (u_node[v]
+                                         * interpolating_polynomials[i, 1, index]
+                                         * interpolating_polynomials[j, 2, index]
+                                         * interpolating_polynomials[k, 3, index])
+            end
+        end
+    end
+end
+end # @muladd
diff --git a/test/test_tree_3d_euler.jl b/test/test_tree_3d_euler.jl
index 98549f2a5ab..1a8c9b1ef71 100644
--- a/test/test_tree_3d_euler.jl
+++ b/test/test_tree_3d_euler.jl
@@ -46,10 +46,28 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_dgsem")
                                                      time_series.affect!.time[1],
                                                      equations) for i in 1:3]
     ref_data = [
-        [2.0908436034061966, 2.0907469131413534, 2.090746913141344, 2.0907469131413516, 4.376022861751483],
-        [1.9121470846408657, 1.911738389880639, 1.9117383898806395, 1.9117383898806386, 3.639369160978031],
-        [1.8960587498057229, 1.896520160078096, 1.8965201600780932, 1.896520160078095, 3.614473843601274]]
-    @test point_data ≈ exact_data atol=1e-1
+        [
+            2.0908436034061966,
+            2.0907469131413534,
+            2.090746913141344,
+            2.0907469131413516,
+            4.376022861751483,
+        ],
+        [
+            1.9121470846408657,
+            1.911738389880639,
+            1.9117383898806395,
+            1.9117383898806386,
+            3.639369160978031,
+        ],
+        [
+            1.8960587498057229,
+            1.896520160078096,
+            1.8965201600780932,
+            1.896520160078095,
+            3.614473843601274,
+        ]]
+    @test point_data≈exact_data atol=1e-1
     @test point_data ≈ ref_data
 end
 

From c3e17db1375ed88e6ff6baafe109a958db24fa09 Mon Sep 17 00:00:00 2001
From: patrickersing <patrick.ersing@liu.se>
Date: Fri, 15 Mar 2024 11:14:54 +0100
Subject: [PATCH 6/9] update test values for julia 1.9

---
 test/test_tree_3d_euler.jl | 30 +++++++++++++++---------------
 1 file changed, 15 insertions(+), 15 deletions(-)

diff --git a/test/test_tree_3d_euler.jl b/test/test_tree_3d_euler.jl
index 1a8c9b1ef71..47669dce2fb 100644
--- a/test/test_tree_3d_euler.jl
+++ b/test/test_tree_3d_euler.jl
@@ -47,25 +47,25 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_dgsem")
                                                      equations) for i in 1:3]
     ref_data = [
         [
-            2.0908436034061966,
-            2.0907469131413534,
-            2.090746913141344,
-            2.0907469131413516,
-            4.376022861751483,
+            1.951156832316166,
+            1.952073047561595,
+            1.9520730475615966,
+            1.9520730475615953,
+            3.814390510967551,
         ],
         [
-            1.9121470846408657,
-            1.911738389880639,
-            1.9117383898806395,
-            1.9117383898806386,
-            3.639369160978031,
+            2.0506452262144363,
+            2.050727319703708,
+            2.0507273197037073,
+            2.0507273197037077,
+            4.203653999433724,
         ],
         [
-            1.8960587498057229,
-            1.896520160078096,
-            1.8965201600780932,
-            1.896520160078095,
-            3.614473843601274,
+            2.046982357537558,
+            2.0463728824399654,
+            2.0463728824399654,
+            2.0463728824399645,
+            4.190033459318115,
         ]]
     @test point_data≈exact_data atol=1e-1
     @test point_data ≈ ref_data

From 4ccc47379456f48be788c9dd3830b54c1a33e673 Mon Sep 17 00:00:00 2001
From: patrickersing <patrick.ersing@liu.se>
Date: Fri, 15 Mar 2024 14:38:55 +0100
Subject: [PATCH 7/9] update news.md

---
 NEWS.md | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index 5b08d51ab89..022252e61a9 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,7 +7,8 @@ for human readability.
 ## Changes in the v0.7 lifecycle
 
 #### Added
-- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`.
+- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension
+  to 1D and 3D on `TreeMesh`.
 
 
 ## Changes when updating to v0.7 from v0.6.x

From d71bee51eb150f1c9d8cdfd50f4dedd358d80ddc Mon Sep 17 00:00:00 2001
From: patrickersing <patrick.ersing@liu.se>
Date: Fri, 15 Mar 2024 14:40:17 +0100
Subject: [PATCH 8/9] rename time_series_dg2d.jl, add comments from code review

---
 examples/tree_1d_dgsem/elixir_euler_source_terms.jl         | 6 +++---
 examples/tree_3d_dgsem/elixir_euler_source_terms.jl         | 4 ++--
 src/callbacks_step/time_series.jl                           | 5 ++---
 src/callbacks_step/time_series_dg_tree.jl                   | 4 +++-
 .../{time_series_dg2d.jl => time_series_dg_unstructured.jl} | 0
 5 files changed, 10 insertions(+), 9 deletions(-)
 rename src/callbacks_step/{time_series_dg2d.jl => time_series_dg_unstructured.jl} (100%)

diff --git a/examples/tree_1d_dgsem/elixir_euler_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_source_terms.jl
index bf6bf7488a9..017f921b635 100644
--- a/examples/tree_1d_dgsem/elixir_euler_source_terms.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_source_terms.jl
@@ -48,9 +48,9 @@ time_series = TimeSeriesCallback(semi, [0.0, 0.33, 1.0], interval = 10)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
-                        stepsize_callback,
-                        time_series)
+                        save_solution, 
+                        time_series,
+                        stepsize_callback)
 
 ###############################################################################
 # run the simulation
diff --git a/examples/tree_3d_dgsem/elixir_euler_source_terms.jl b/examples/tree_3d_dgsem/elixir_euler_source_terms.jl
index 0cf615fda10..021fd09f316 100644
--- a/examples/tree_3d_dgsem/elixir_euler_source_terms.jl
+++ b/examples/tree_3d_dgsem/elixir_euler_source_terms.jl
@@ -48,8 +48,8 @@ time_series = TimeSeriesCallback(semi,
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
                         save_solution,
-                        stepsize_callback,
-                        time_series)
+                        time_series,
+                        stepsize_callback)
 
 ###############################################################################
 # run the simulation
diff --git a/src/callbacks_step/time_series.jl b/src/callbacks_step/time_series.jl
index cf90a47900c..ae18c85700d 100644
--- a/src/callbacks_step/time_series.jl
+++ b/src/callbacks_step/time_series.jl
@@ -23,8 +23,7 @@ After the last time step, the results are stored in an HDF5 file `filename` in d
 The real data type `RealT` and data type for solution variables `uEltype` default to the respective
 types used in the solver and the cache.
 
-Currently this callback is only implemented for [`TreeMesh`](@ref) in 2D
-and [`UnstructuredMesh2D`](@ref).
+Currently this callback is only implemented for [`TreeMesh`](@ref) and [`UnstructuredMesh2D`](@ref).
 """
 mutable struct TimeSeriesCallback{RealT <: Real, uEltype <: Real, SolutionVariables,
                                   VariableNames, Cache}
@@ -219,5 +218,5 @@ end
 
 include("time_series_dg.jl")
 include("time_series_dg_tree.jl")
-include("time_series_dg2d.jl")
+include("time_series_dg_unstructured.jl")
 end # @muladd
diff --git a/src/callbacks_step/time_series_dg_tree.jl b/src/callbacks_step/time_series_dg_tree.jl
index f6e844021e2..37d4e6ea705 100644
--- a/src/callbacks_step/time_series_dg_tree.jl
+++ b/src/callbacks_step/time_series_dg_tree.jl
@@ -84,7 +84,7 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates,
     return interpolating_polynomials
 end
 
-# Record the solution variables at each given point
+# Record the solution variables at each given point for the 1D case
 function record_state_at_points!(point_data, u, solution_variables,
                                  n_solution_variables,
                                  mesh::TreeMesh{1}, equations, dg::DG,
@@ -116,6 +116,7 @@ function record_state_at_points!(point_data, u, solution_variables,
     end
 end
 
+# Record the solution variables at each given point for the 2D case
 function record_state_at_points!(point_data, u, solution_variables,
                                  n_solution_variables,
                                  mesh::TreeMesh{2},
@@ -148,6 +149,7 @@ function record_state_at_points!(point_data, u, solution_variables,
     end
 end
 
+# Record the solution variables at each given point for the 3D case
 function record_state_at_points!(point_data, u, solution_variables,
                                  n_solution_variables,
                                  mesh::TreeMesh{3}, equations, dg::DG,
diff --git a/src/callbacks_step/time_series_dg2d.jl b/src/callbacks_step/time_series_dg_unstructured.jl
similarity index 100%
rename from src/callbacks_step/time_series_dg2d.jl
rename to src/callbacks_step/time_series_dg_unstructured.jl

From 29fbefde48331b51140d661f63146e9468f14ce7 Mon Sep 17 00:00:00 2001
From: patrickersing <patrick.ersing@liu.se>
Date: Fri, 15 Mar 2024 14:44:50 +0100
Subject: [PATCH 9/9] apply formatter

---
 examples/tree_1d_dgsem/elixir_euler_source_terms.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_source_terms.jl
index 017f921b635..cb8a09057d9 100644
--- a/examples/tree_1d_dgsem/elixir_euler_source_terms.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_source_terms.jl
@@ -48,7 +48,7 @@ time_series = TimeSeriesCallback(semi, [0.0, 0.33, 1.0], interval = 10)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution, 
+                        save_solution,
                         time_series,
                         stepsize_callback)