Add functionality for TimeSeries callback on UnstructuredMesh2D

examples/unstructured_2d_dgsem/elixir_euler_basic.jl

@@ -58,11 +58,15 @@ save_solution = SaveSolutionCallback(interval = 10,
 
 stepsize_callback = StepsizeCallback(cfl = 0.9)
 
+time_series = TimeSeriesCallback(semi, [(2.0, -0.5), (0.28, 1 / 3), (1.87, 2 / 3)];
+                                 interval = 10)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
                         alive_callback,
                         save_restart,
                         save_solution,
+                        time_series,
                         stepsize_callback)
 
 ###############################################################################

src/callbacks_step/time_series_dg.jl

@@ -5,8 +5,10 @@
 @muladd begin
 #! format: noindent
 
-# Store time series file for a TreeMesh with a DG solver
-function save_time_series_file(time_series_callback, mesh::TreeMesh, equations, dg::DG)
+# Store time series file for a DG solver
+function save_time_series_file(time_series_callback,
+                               mesh::Union{TreeMesh, UnstructuredMesh2D},
+                               equations, dg::DG)
     @unpack (interval, solution_variables, variable_names,
     output_directory, filename, point_coordinates,
     point_data, time, step, time_series_cache) = time_series_callback

src/callbacks_step/time_series_dg2d.jl

@@ -6,7 +6,9 @@
 #! 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::Union{TreeMesh{2}, UnstructuredMesh2D},
+                                  dg, cache)
     # Determine element ids for point coordinates
     element_ids = get_elements_by_coordinates(point_coordinates, mesh, dg, cache)
 
@@ -68,6 +70,58 @@ function get_elements_by_coordinates!(element_ids, coordinates, mesh::TreeMesh,
     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.
+# We use the barycenter of each element to determine if a given coordinate (x,y)
+# lies within said element. The shortest distance between the point and all the
+# barycenters "wins".
+function get_elements_by_coordinates!(element_ids, coordinates,
+                                      mesh::UnstructuredMesh2D,
+                                      dg, cache)
+    if length(element_ids) != size(coordinates, 2)
+        throw(DimensionMismatch("storage length for element ids does not match the number of coordinates"))
+    end
+
+    # Reset element ids - 0 indicates "not (yet) found"
+    element_ids .= 0
+
+    # Compute and save the barycentric coordinate on each element
+    bary_centers = zeros(eltype(mesh.corners), 2, mesh.n_elements)
+    calc_bary_centers!(bary_centers, dg, cache)
+
+    # Iterate over coordinates
+    for index in 1:length(element_ids)
+        point = SVector(coordinates[1, index], 
+                        coordinates[2, index])
+
+        # Search for the element with minimal distance between the point
+        # we are interested in and the barycenter of the element.
+        element_ids[index] = argmin(eachelement(dg, cache)) do element
+            bary_center = SVector(bary_centers[1, element],
+                                  bary_centers[2, element])
+            # Compute the squared Euclidean distance between the point we are
+            # interested in and the barycenter of the element.
+            # We compute the squared norm to avoid computing computationally
+            # more expensive square roots.
+            sum(abs2, point - bary_center)
+        end
+    end
+
+    return element_ids
+end
+
+# Use the available `node_coordinates` on each element to compute and save the barycenter.
+@inline function calc_bary_centers!(bary_centers, dg, cache)
+    n = nnodes(dg)
+    @views for element in eachelement(dg, cache)
+        bary_centers[1, element] = sum(cache.elements.node_coordinates[1, :, :,
+                                                                       element]) / n^2
+        bary_centers[2, element] = sum(cache.elements.node_coordinates[2, :, :,
+                                                                       element]) / n^2
+    end
+    return nothing
+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)
@@ -106,8 +160,139 @@ function calc_interpolating_polynomials!(interpolating_polynomials, coordinates,
     return interpolating_polynomials
 end
 
-function calc_interpolating_polynomials(coordinates, element_ids, mesh::TreeMesh, dg,
-                                        cache)
+function calc_interpolating_polynomials!(interpolating_polynomials, coordinates,
+                                         element_ids,
+                                         mesh::UnstructuredMesh2D, dg::DGSEM, cache)
+    @unpack nodes = dg.basis
+
+    wbary = barycentric_weights(nodes)
+
+    # Helper array for a straight-sided quadrilateral element
+    corners = zeros(eltype(mesh.corners), 4, 2)
+
+    for index in 1:length(element_ids)
+        # Construct point
+        x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh)))
+
+        # Convert to unit coordinates; procedure differs for straight-sided
+        # versus curvilinear elements
+        element = element_ids[index]
+        if !mesh.element_is_curved[element]
+            for j in 1:2, i in 1:4
+                # Pull the (x,y) values of the element corners from the global corners array
+                corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]]
+            end
+            # Compute coordinates in reference system
+            unit_coordinates = invert_bilinear_interpolation(mesh, x, corners)
+
+            # Sanity check that the computed `unit_coordinates` indeed recover the desired point `x`
+            x_check = straight_side_quad_map(unit_coordinates[1], unit_coordinates[2],
+                                             corners)
+            if !isapprox(x[1], x_check[1], atol = 1e-13) ||
+               !isapprox(x[2], x_check[2], atol = 1e-13)
+                error("failed to compute computational coordinates for the time series point $(x), closest candidate was $(x_check)")
+            end
+        else # mesh.element_is_curved[element]
+            unit_coordinates = invert_transfinite_interpolation(mesh, x,
+                                                                view(mesh.surface_curves,
+                                                                     :, element))
+
+            # Sanity check that the computed `unit_coordinates` indeed recover the desired point `x`
+            x_check = transfinite_quad_map(unit_coordinates[1], unit_coordinates[2],
+                                           view(mesh.surface_curves, :, element))
+            if !isapprox(x[1], x_check[1], atol = 1e-13) ||
+               !isapprox(x[2], x_check[2], atol = 1e-13)
+                error("failed to compute computational coordinates for the time series point $(x), closest candidate was $(x_check)")
+            end
+        end
+
+        # 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
+
+# Use a Newton iteration to determine the computational coordinates
+# (xi, eta) of given (x,y) `point` that is given in physical coordinates
+# by inverting the transformation. For straight-sided elements this
+# amounts to inverting a bi-linear interpolation. For curved
+# elements we invert the transfinite interpolation with linear blending.
+# The residual function for the Newton iteration is
+#    r(xi,eta) = X(xi,eta) - point
+# and the Jacobian entries are computed accordingly from either
+# `straight_side_quad_map_metrics` or `transfinite_quad_map_metrics`.
+# We exploit the 2x2 nature of the problem and directly compute the matrix
+# inverse to make things faster. The implementations below are inspired by
+# an answer on Stack Overflow (https://stackoverflow.com/a/18332009) where
+# the author explicitly states that their code is released to the public domain.
+@inline function invert_bilinear_interpolation(mesh::UnstructuredMesh2D, point,
+                                               element_corners)
+    # Initial guess for the point (center of the reference element)
+    xi = zero(eltype(point))
+    eta = zero(eltype(point))
+    for k in 1:5 # Newton's method should converge quickly
+        # Compute current x and y coordinate and the Jacobian matrix
+        # J = (X_xi, X_eta; Y_xi, Y_eta)
+        x, y = straight_side_quad_map(xi, eta, element_corners)
+        J11, J12, J21, J22 = straight_side_quad_map_metrics(xi, eta, element_corners)
+
+        # Compute residuals for the Newton teration for the current (x, y) coordinate
+        r1 = x - point[1]
+        r2 = y - point[2]
+
+        # Newton update that directly applies the inverse of the 2x2 Jacobian matrix
+        inv_detJ = inv(J11 * J22 - J12 * J21)
+
+        # Update with explicitly inverted Jacobian
+        xi = xi - inv_detJ * (J22 * r1 - J12 * r2)
+        eta = eta - inv_detJ * (-J21 * r1 + J11 * r2)
+
+        # Ensure updated point is in the reference element
+        xi = min(max(xi, -1), 1)
+        eta = min(max(eta, -1), 1)
+    end
+
+    return SVector(xi, eta)
+end
+
+@inline function invert_transfinite_interpolation(mesh::UnstructuredMesh2D, point,
+                                                  surface_curves::AbstractVector{<:CurvedSurface})
+    # Initial guess for the point (center of the reference element)
+    xi = zero(eltype(point))
+    eta = zero(eltype(point))
+    for k in 1:5 # Newton's method should converge quickly
+        # Compute current x and y coordinate and the Jacobian matrix
+        # J = (X_xi, X_eta; Y_xi, Y_eta)
+        x, y = transfinite_quad_map(xi, eta, surface_curves)
+        J11, J12, J21, J22 = transfinite_quad_map_metrics(xi, eta, surface_curves)
+
+        # Compute residuals for the Newton teration for the current (x,y) coordinate
+        r1 = x - point[1]
+        r2 = y - point[2]
+
+        # Newton update that directly applies the inverse of the 2x2 Jacobian matrix
+        inv_detJ = inv(J11 * J22 - J12 * J21)
+
+        # Update with explicitly inverted Jacobian
+        xi = xi - inv_detJ * (J22 * r1 - J12 * r2)
+        eta = eta - inv_detJ * (-J21 * r1 + J11 * r2)
+
+        # Ensure updated point is in the reference element
+        xi = min(max(xi, -1), 1)
+        eta = min(max(eta, -1), 1)
+    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))
@@ -121,8 +306,8 @@ end
 # Record the solution variables at each given point
 function record_state_at_points!(point_data, u, solution_variables,
                                  n_solution_variables,
-                                 mesh::TreeMesh{2}, equations, dg::DG,
-                                 time_series_cache)
+                                 mesh::Union{TreeMesh{2}, UnstructuredMesh2D},
+                                 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