Allow disabling reinterpolation in 1D plots (#2258)

JoshuaLampert · DanielDoehring · ranocha · web-flow · commit 810984592e06 · 2025-02-06T09:15:24.000+01:00
* allow disabling reinterpolation in 1D plots

* Apply suggestions from code review

Co-authored-by: Daniel Doehring &lt;daniel.doehring@rwth-aachen.de&gt;

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha &lt;ranocha@users.noreply.github.com&gt;

* use reinterpolate = false for FDSBP

* format

* simplify for reinterpolate = false (don't do anything)

* clarify docstring

* reintroduce blank line

* Update src/visualization/types.jl

Co-authored-by: Hendrik Ranocha &lt;ranocha@users.noreply.github.com&gt;

* Update src/visualization/utilities.jl

Co-authored-by: Daniel Doehring &lt;daniel.doehring@rwth-aachen.de&gt;

---------

Co-authored-by: Daniel Doehring &lt;daniel.doehring@rwth-aachen.de&gt;
Co-authored-by: Hendrik Ranocha &lt;ranocha@users.noreply.github.com&gt;
diff --git a/src/visualization/recipes_plots.jl b/src/visualization/recipes_plots.jl
@@ -190,11 +190,13 @@ end
 RecipesBase.@recipe function f(u, semi::SemidiscretizationHyperbolic{<:TreeMesh};
                                solution_variables = nothing,
                                grid_lines = true, max_supported_level = 11,
-                               nvisnodes = nothing, slice = :xy,
-                               point = (0.0, 0.0, 0.0), curve = nothing)
+                               nvisnodes = nothing,
+                               reinterpolate = default_reinterpolate(semi.solver),
+                               slice = :xy, point = (0.0, 0.0, 0.0), curve = nothing)
     # Create a PlotData1D or PlotData2D object depending on the dimension.
     if ndims(semi) == 1
-        return PlotData1D(u, semi; solution_variables, nvisnodes, slice, point, curve)
+        return PlotData1D(u, semi; solution_variables, nvisnodes, reinterpolate,
+                          slice, point, curve)
     else
         return PlotData2D(u, semi;
                           solution_variables, grid_lines, max_supported_level,
diff --git a/src/visualization/types.jl b/src/visualization/types.jl
@@ -516,7 +516,9 @@ end
 
 """
     PlotData1D(u, semi [or mesh, equations, solver, cache];
-               solution_variables=nothing, nvisnodes=nothing)
+               solution_variables=nothing, nvisnodes=nothing,
+               reinterpolate=Trixi.default_reinterpolate(solver),
+               slice=:x, point=(0.0, 0.0, 0.0), curve=nothing)
 
 Create a new `PlotData1D` object that can be used for visualizing 1D DGSEM solution data array
 `u` with `Plots.jl`. All relevant geometrical information is extracted from the
@@ -531,7 +533,10 @@ e.g., to visualize an analytical solution.
 
 `nvisnodes` specifies the number of visualization nodes to be used. If it is `nothing`,
 twice the number of solution DG nodes are used for visualization, and if set to `0`,
-exactly the number of nodes in the DG elements are used.
+exactly the number of nodes as in the DG elements are used. The solution is interpolated to these
+nodes for plotting. If `reinterpolate` is `false`, the solution is not interpolated to the `visnodes`,
+i.e., the solution is plotted at the solution nodes. In this case `nvisnodes` is ignored. By default,
+`reinterpolate` is set to `false` for `FDSBP` approximations and to `true` otherwise.
 
 When visualizing data from a two-dimensional simulation, a 1D slice is extracted for plotting.
 `slice` specifies the axis along which the slice is extracted and may be `:x`, or `:y`.
@@ -560,6 +565,7 @@ end
 
 function PlotData1D(u, mesh::TreeMesh, equations, solver, cache;
                     solution_variables = nothing, nvisnodes = nothing,
+                    reinterpolate = default_reinterpolate(solver),
                     slice = :x, point = (0.0, 0.0, 0.0), curve = nothing,
                     variable_names = nothing)
     solution_variables_ = digest_solution_variables(equations, solution_variables)
@@ -574,7 +580,7 @@ function PlotData1D(u, mesh::TreeMesh, equations, solver, cache;
 
     if ndims(mesh) == 1
         x, data, mesh_vertices_x = get_data_1d(original_nodes, unstructured_data,
-                                               nvisnodes)
+                                               nvisnodes, reinterpolate)
         orientation_x = 1
 
         # Special care is required for first-order FV approximations since the nodes are the
@@ -608,7 +614,8 @@ function PlotData1D(u, mesh::TreeMesh, equations, solver, cache;
         else
             x, data, mesh_vertices_x = unstructured_2d_to_1d(original_nodes,
                                                              unstructured_data,
-                                                             nvisnodes, slice, point)
+                                                             nvisnodes, reinterpolate,
+                                                             slice, point)
         end
     else # ndims(mesh) == 3
         if curve !== nothing
@@ -619,7 +626,8 @@ function PlotData1D(u, mesh::TreeMesh, equations, solver, cache;
         else
             x, data, mesh_vertices_x = unstructured_3d_to_1d(original_nodes,
                                                              unstructured_data,
-                                                             nvisnodes, slice, point)
+                                                             nvisnodes, reinterpolate,
+                                                             slice, point)
         end
     end
 
@@ -629,6 +637,7 @@ end
 
 function PlotData1D(u, mesh, equations, solver, cache;
                     solution_variables = nothing, nvisnodes = nothing,
+                    reinterpolate = default_reinterpolate(solver),
                     slice = :x, point = (0.0, 0.0, 0.0), curve = nothing,
                     variable_names = nothing)
     solution_variables_ = digest_solution_variables(equations, solution_variables)
@@ -643,7 +652,7 @@ function PlotData1D(u, mesh, equations, solver, cache;
 
     if ndims(mesh) == 1
         x, data, mesh_vertices_x = get_data_1d(original_nodes, unstructured_data,
-                                               nvisnodes)
+                                               nvisnodes, reinterpolate)
         orientation_x = 1
     elseif ndims(mesh) == 2
         # Create a 'PlotData2DTriangulated' object so a triangulation can be used when extracting relevant data.
diff --git a/src/visualization/utilities.jl b/src/visualization/utilities.jl
@@ -391,54 +391,66 @@ function get_data_2d(center_level_0, length_level_0, leaf_cells, coordinates, le
     return xs, ys, node_centered_data, mesh_vertices_x, mesh_vertices_y
 end
 
+# For finite difference methods (FDSBP), we do not want to reinterpolate the data, but use the same
+# nodes as input nodes. For DG methods, we usually want to reinterpolate the data on an equidistant grid.
+default_reinterpolate(solver) = solver isa FDSBP ? false : true
+
 # Extract data from a 1D DG solution and prepare it for visualization as a line plot.
 # This returns a tuple with
 # - x coordinates
 # - nodal 1D data
 #
 # Note: This is a low-level function that is not considered as part of Trixi's interface and may
-#       thus be changed in future releases.
-function get_data_1d(original_nodes, unstructured_data, nvisnodes)
+# thus be changed in future releases.
+function get_data_1d(original_nodes, unstructured_data, nvisnodes, reinterpolate)
     # Get the dimensions of u; where n_vars is the number of variables, n_nodes the number of nodal values per element and n_elements the total number of elements.
     n_nodes, n_elements, n_vars = size(unstructured_data)
 
-    # Set the amount of nodes visualized according to nvisnodes.
-    if nvisnodes === nothing
-        max_nvisnodes = 2 * n_nodes
-    elseif nvisnodes == 0
-        max_nvisnodes = n_nodes
+    # If `reinterpolate` is `false`, we do nothing.
+    # If `reinterpolate` is `true`, the output nodes are equidistantly spaced.
+    if reinterpolate == false
+        interpolated_nodes = original_nodes
+        interpolated_data = unstructured_data
     else
-        @assert nvisnodes>=2 "nvisnodes must be zero or >= 2"
-        max_nvisnodes = nvisnodes
-    end
+        # Set the amount of nodes visualized according to nvisnodes.
+        if nvisnodes === nothing
+            max_nvisnodes = 2 * n_nodes
+        elseif nvisnodes == 0
+            max_nvisnodes = n_nodes
+        else
+            @assert nvisnodes>=2 "nvisnodes must be zero or >= 2"
+            max_nvisnodes = nvisnodes
+        end
 
-    interpolated_nodes = Array{eltype(original_nodes), 2}(undef, max_nvisnodes,
-                                                          n_elements)
-    interpolated_data = Array{eltype(unstructured_data), 3}(undef, max_nvisnodes,
-                                                            n_elements, n_vars)
+        interpolated_nodes = Array{eltype(original_nodes), 2}(undef, max_nvisnodes,
+                                                              n_elements)
+        interpolated_data = Array{eltype(unstructured_data), 3}(undef, max_nvisnodes,
+                                                                n_elements, n_vars)
 
-    for j in 1:n_elements
-        # Interpolate on an equidistant grid.
-        interpolated_nodes[:, j] .= range(original_nodes[1, 1, j],
-                                          original_nodes[1, end, j],
-                                          length = max_nvisnodes)
-    end
+        for j in 1:n_elements
+            # Interpolate on an equidistant grid.
+            interpolated_nodes[:, j] .= range(original_nodes[1, 1, j],
+                                              original_nodes[1, end, j],
+                                              length = max_nvisnodes)
+        end
 
-    nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes)
-    nodes_out = collect(range(-1, 1, length = max_nvisnodes))
-
-    # Calculate vandermonde matrix for interpolation.
-    vandermonde = polynomial_interpolation_matrix(nodes_in, nodes_out)
-
-    # Iterate over all variables.
-    for v in 1:n_vars
-        # Interpolate data for each element.
-        for element in 1:n_elements
-            multiply_scalar_dimensionwise!(@view(interpolated_data[:, element, v]),
-                                           vandermonde,
-                                           @view(unstructured_data[:, element, v]))
+        nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes)
+        nodes_out = collect(range(-1, 1, length = max_nvisnodes))
+
+        # Calculate vandermonde matrix for interpolation.
+        vandermonde = polynomial_interpolation_matrix(nodes_in, nodes_out)
+
+        # Iterate over all variables.
+        for v in 1:n_vars
+            # Interpolate data for each element.
+            for element in 1:n_elements
+                multiply_scalar_dimensionwise!(@view(interpolated_data[:, element, v]),
+                                               vandermonde,
+                                               @view(unstructured_data[:, element, v]))
+            end
         end
     end
+
     # Return results after data is reshaped
     return vec(interpolated_nodes), reshape(interpolated_data, :, n_vars),
            vcat(original_nodes[1, 1, :], original_nodes[1, end, end])
@@ -675,8 +687,8 @@ function unstructured_3d_to_2d(unstructured_data, coordinates, levels,
 end
 
 # Convert 2d unstructured data to 1d slice and interpolate them.
-function unstructured_2d_to_1d(original_nodes, unstructured_data, nvisnodes, slice,
-                               point)
+function unstructured_2d_to_1d(original_nodes, unstructured_data, nvisnodes,
+                               reinterpolate, slice, point)
     if slice === :x
         slice_dimension = 2
         other_dimension = 1
@@ -749,7 +761,7 @@ function unstructured_2d_to_1d(original_nodes, unstructured_data, nvisnodes, sli
     end
 
     return get_data_1d(reshape(new_nodes[:, 1:new_id], 1, n_nodes_in, new_id),
-                       new_unstructured_data[:, 1:new_id, :], nvisnodes)
+                       new_unstructured_data[:, 1:new_id, :], nvisnodes, reinterpolate)
 end
 
 # Calculate the arc length of a curve given by ndims x npoints point coordinates (piece-wise linear approximation)
@@ -1085,8 +1097,8 @@ function get_value_at_point(point, nodes, data)
 end
 
 # Convert 3d unstructured data to 1d slice and interpolate them.
-function unstructured_3d_to_1d(original_nodes, unstructured_data, nvisnodes, slice,
-                               point)
+function unstructured_3d_to_1d(original_nodes, unstructured_data, nvisnodes,
+                               reinterpolate, slice, point)
     if slice === :x
         slice_dimension = 1
         other_dimensions = [2, 3]
@@ -1178,7 +1190,7 @@ function unstructured_3d_to_1d(original_nodes, unstructured_data, nvisnodes, sli
     end
 
     return get_data_1d(reshape(new_nodes[:, 1:new_id], 1, n_nodes_in, new_id),
-                       new_unstructured_data[:, 1:new_id, :], nvisnodes)
+                       new_unstructured_data[:, 1:new_id, :], nvisnodes, reinterpolate)
 end
 
 # Interpolate unstructured DG data to structured data (cell-centered)
diff --git a/test/test_visualization.jl b/test/test_visualization.jl
@@ -185,6 +185,7 @@ end
         pd = PlotData1D(sol)
 
         @test_nowarn_mod Plots.plot(sol)
+        @test_nowarn_mod Plots.plot(sol, reinterpolate = false)
         @test_nowarn_mod Plots.plot(pd)
         @test_nowarn_mod Plots.plot(pd["p"])
         @test_nowarn_mod Plots.plot(getmesh(pd))
