Add tutorial on visualization with Plots.jl and interpolation

NEWS.md

@@ -6,6 +6,11 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 ## Version 0.5.4
 
+### Features
+
+- Added a plotting recipe to plot particles of individual systems (#1298).
+- Added a tutorial on plotting particles and interpolated fields (#1298).
+
 ### Important Bugfixes
 
 - Fixed mathematical inconsistencies in the SPH documentation, including incorrect

docs/literate/src/tut_visualization.jl

@@ -0,0 +1,93 @@
+# # [Visualizing particle data with Plots.jl](@id tut_visualization)
+
+# In this tutorial, we run the two-dimensional vortex street from
+# [`examples/fluid/vortex_street_2d.jl`](https://github.com/trixi-framework/TrixiParticles.jl/blob/main/examples/fluid/vortex_street_2d.jl)
+# and visualize the particle data with [`Plots.jl`](https://github.com/juliaplots/plots.jl).
+
+using TrixiParticles
+using Plots
+#src # Reset GR's process-wide color table, which can be exhausted by earlier tutorials.
+Plots.closeall() # hide
+
+# The example defines the particle spacing as `particle_spacing_factor * cylinder_diameter`.
+# We deliberately use a very coarse particle resolution. This makes the distinction between
+# the discrete particles and the interpolated field in the next section clear.
+# To remove visual clutter, we disable the info callback.
+# Since we visualize with Plots.jl, we also disable the saving callback.
+trixi_include(@__MODULE__,
+              joinpath(examples_dir(), "fluid", "vortex_street_2d.jl");
+              particle_spacing_factor=0.2,
+              info_callback=nothing, saving_callback=nothing);
+nothing # hide
+
+# ## Visualizing discrete particles
+
+# SPH stores the solution on moving particles. The standard plotting recipe provides the
+# quickest way to inspect their distribution at the final time. We color the fluid particles
+# by the magnitude of the velocity stored on each particle.
+v_ode, _ = sol.u[end].x
+v_fluid = TrixiParticles.wrap_v(v_ode, fluid_system, semi)
+
+active_particles = TrixiParticles.eachparticle(fluid_system)
+particle_velocity = TrixiParticles.current_velocity(v_fluid,
+                                                    fluid_system)[:, active_particles]
+particle_velocity_magnitude = vec(sqrt.(sum(abs2, particle_velocity; dims=1)))
+
+particle_plot = plot(fluid_system, sol; zcolor=particle_velocity_magnitude, color=:viridis,
+                     xlims=(0.25, 1.8), ylims=(0.1, 0.9), legend=false,
+                     xlabel="x", ylabel="y", colorbar=true, colorbar_title="|v|",
+                     size=(900, 450))
+plot!(particle_plot; dpi=200) # hide
+savefig(particle_plot, "tut_visualization_particles.png") # hide
+nothing # hide
+
+# ![Particle visualization of the vortex street](tut_visualization_particles.png)
+
+# ## Interpolating particle data onto a regular grid
+
+# Smoothed particle hydrodynamics (SPH) represents a continuous (smoothed) field
+# by a discrete set of particles. While visualizing individual particles is straightforward
+# and often sufficient, in order to visualize the actual field approximation,
+# the particle data must be interpolated.
+#
+# Importantly, interpolation does not add physical resolution: features that are not resolved
+# by the particles cannot be recovered by choosing a finer interpolation grid.
+# It simply visualizes the SPH approximation instead of only the interpolation points.
+
+# [`interpolate_plane_2d`](@ref) constructs regularly spaced sample points between two corners
+# and uses the SPH kernel to reconstruct the requested fields there. The interpolation spacing
+# is one quarter of the particle spacing, so the plot contains many more pixels than
+# the simulation contains particles.
+interpolation_min = [0.0, 0.0]
+interpolation_max = domain_size
+interpolation_spacing = particle_spacing / 4
+
+interpolated = interpolate_plane_2d(interpolation_min, interpolation_max,
+                                    interpolation_spacing, semi, fluid_system, sol)
+interpolated_velocity_magnitude = vec(sqrt.(sum(abs2, interpolated.velocity; dims=1)))
+nothing # hide
+
+# The returned named tuple also contains `pressure`, `density`, `neighbor_count`, and
+# `computed_density`. Here we visualize the magnitude of the interpolated velocity.
+interpolated_plot = scatter(interpolated.point_coords[1, :],
+                            interpolated.point_coords[2, :];
+                            marker_z=interpolated_velocity_magnitude,
+                            color=:viridis,
+                            marker=:square, markerstrokewidth=0, markersize=2.5,
+                            aspect_ratio=:equal, size=(900, 450),
+                            xlims=(0.25, 1.8), ylims=(0.1, 0.9), xlabel="x", ylabel="y",
+                            label=nothing, colorbar_title="|v|")
+plot!(interpolated_plot; dpi=200) # hide
+savefig(interpolated_plot, "tut_visualization_interpolated_velocity.png") # hide
+nothing # hide
+
+# ![Interpolated velocity magnitude](tut_visualization_interpolated_velocity.png)
+
+# Compared with the visibly discrete particle distribution, the interpolated field shows
+# much more detail, representing the continuous SPH approximation of the solution.
+
+# To write the same reconstruction as a VTI image for ParaView, replace the interpolation call
+# above with [`interpolate_plane_2d_vtk`](@ref):
+interpolate_plane_2d_vtk(interpolation_min, interpolation_max, interpolation_spacing,
+                         semi, fluid_system, sol; filename="vortex_street_velocity")
+nothing # hide

docs/make.jl

@@ -34,14 +34,20 @@ function copy_file(filename, replaces...;
     write(new_file, content)
 end
 
+tutorials_output_dir = joinpath("docs", "src", "tutorials")
+rm(tutorials_output_dir; recursive=true, force=true)
+mkpath(tutorials_output_dir)
+
 Literate.markdown(joinpath("docs", "literate", "src", "tut_setup.jl"),
-                  joinpath("docs", "src", "tutorials"))
+                  tutorials_output_dir)
 Literate.markdown(joinpath("docs", "literate", "src", "tut_custom_kernel.jl"),
-                  joinpath("docs", "src", "tutorials"))
+                  tutorials_output_dir)
 Literate.markdown(joinpath("docs", "literate", "src", "tut_rigid_body_fsi.jl"),
-                  joinpath("docs", "src", "tutorials"))
+                  tutorials_output_dir)
 Literate.markdown(joinpath("docs", "literate", "src", "tut_packing.jl"),
-                  joinpath("docs", "src", "tutorials"))
+                  tutorials_output_dir)
+Literate.markdown(joinpath("docs", "literate", "src", "tut_visualization.jl"),
+                  tutorials_output_dir)
 
 copy_file("AUTHORS.md",
           "in the [LICENSE.md](LICENSE.md) file" => "under [License](@ref)")
@@ -87,7 +93,9 @@ makedocs(sitename="TrixiParticles.jl",
                      "Setting up your simulation from scratch" => joinpath("tutorials",
                                                                            "tut_setup.md"),
                      "Modifying or extending components of TrixiParticles.jl within a simulation file" => joinpath("tutorials",
-                                                                                                                   "tut_custom_kernel.md")
+                                                                                                                   "tut_custom_kernel.md"),
+                     "Visualizing particle data with Plots.jl" => joinpath("tutorials",
+                                                                           "tut_visualization.md")
                  ],
                  "Fluid-Structure Interaction" => [
                      "Fluid-structure interaction with rigid bodies" => joinpath("tutorials",

docs/src/tutorial.md

@@ -6,7 +6,8 @@
 
 1. [Setting up your simulation from scratch](tutorials/tut_setup.md): learn the structure of a simulation file and run a complete WCSPH example.
 2. [Modifying or extending components of TrixiParticles.jl within a simulation file](tutorials/tut_custom_kernel.md): replace selected parts of an existing setup without cloning the package.
-3. [Particle packing tutorial](tutorials/tut_packing.md): build a body-fitted particle configuration for complex geometries.
+3. [Visualizing particle data with Plots.jl](tutorials/tut_visualization.md): compare discrete particle plots with interpolated fields on a Cartesian grid.
+4. [Particle packing tutorial](tutorials/tut_packing.md): build a body-fitted particle configuration for complex geometries.
 
 ## Tutorials
 
@@ -38,6 +39,20 @@ directly in the file you run.
 - Focus: `trixi_include`, custom kernels, rapid iteration
 - Choose this if: you want to prototype changes without cloning and modifying the package
 
+### [Visualizing particle data with Plots.jl](tutorials/tut_visualization.md)
+
+```@raw html
+<img src="../tutorials/tut_visualization_interpolated_velocity.png"
+     alt="Velocity from a coarse vortex-street simulation interpolated onto a regular grid"
+     style="max-width: 360px; width: 100%; border-radius: 12px;" />
+```
+
+Compare the velocity carried by moving SPH particles with its kernel reconstruction on a
+Cartesian grid in a deliberately coarse vortex-street simulation.
+
+- Focus: particle plots, plane interpolation, regular-grid visualization
+- Choose this if: you want to postprocess particle data or export a smooth field for plotting
+
 ### [Particle packing tutorial](tutorials/tut_packing.md)
 
 ```@raw html

docs/src/visualization.md

@@ -1,5 +1,8 @@
 # Visualization
 
+For instructions on how to visualize simulation data in Julia, see the tutorial on
+[visualizing particle data with Plots.jl](@ref tut_visualization).
+
 ## Export VTK files
 You can export particle data as VTK files by using the [`SolutionSavingCallback`](@ref).
 All [predefined examples](examples.md) already use this callback to export VTK files to the `out`

src/visualization/recipes_plots.jl

@@ -13,55 +13,67 @@ RecipesBase.@recipe function f(sol::TrixiParticlesODESolution)
     return sol.u[end].x..., sol.prob.p.semi
 end
 
-# GPU version
-RecipesBase.@recipe function f(v_ode::AbstractGPUArray, u_ode::AbstractGPUArray,
-                               semi::Semidiscretization)
-    # Move GPU data to the CPU
-    v_ode_, u_ode_, semi_ = transfer2cpu(v_ode, u_ode, semi)
+RecipesBase.@recipe function f(system::AbstractSystem, sol::TrixiParticlesODESolution)
+    # Redirect everything to the single-system recipe
+    return system, sol.u[end].x..., sol.prob.p.semi
+end
 
-    # Redirect everything to the next recipe
-    return v_ode_, u_ode_, semi_
+function get_system_plot_data(u_ode, system, semi, particle_spacing)
+    u = wrap_u(u_ode, system, semi)
+    periodic_box = get_neighborhood_search(system, semi).periodic_box
+    coordinates = PointNeighbors.periodic_coords(active_coordinates(u, system),
+                                                 periodic_box)
+
+    x = collect(coordinates[1, :])
+    y = collect(coordinates[2, :])
+
+    if particle_spacing < 0
+        particle_spacing = 0.0
+    end
+
+    x_min, x_max = extrema(x)
+    y_min, y_max = extrema(y)
+
+    # Add one particle radius around the center of the particles
+    # to obtain the domain size.
+    x_min -= 0.5particle_spacing
+    x_max += 0.5particle_spacing
+    y_min -= 0.5particle_spacing
+    y_max += 0.5particle_spacing
+
+    return (; x, y, x_min, x_max, y_min, y_max, particle_spacing,
+            label=timer_name(system))
+end
+
+RecipesBase.@recipe function f(system::AbstractSystem, v_ode::AbstractArray,
+                               u_ode::AbstractArray, semi::Semidiscretization)
+    # Move data to the CPU if on the GPU.
+    # `transfer2cpu` also validates that `system` is in `semi.systems`, even on the CPU.
+    v_ode_, u_ode_, system_, semi_ = transfer2cpu(v_ode, u_ode, system, semi)
+
+    particle_spacing = system_.initial_condition.particle_spacing
+    system_data = get_system_plot_data(u_ode_, system_, semi_, particle_spacing)
+
+    return (system_, system_data)
 end
 
 RecipesBase.@recipe function f(v_ode::AbstractArray, u_ode::AbstractArray,
                                semi::Semidiscretization;
                                particle_spacings=TrixiParticles.particle_spacings(semi),
                                size=(600, 400), # Default size
                                xlims=(-Inf, Inf), ylims=(-Inf, Inf))
-    # We need to split this in two recipes in order to find the minimum and maximum
-    # coordinates across all systems.
-    # In this first recipe, we collect the data for each system,
-    # and then pass it to the next recipe.
-    systems_data = map(enumerate(semi.systems)) do (i, system)
-        u = wrap_u(u_ode, system, semi)
-        periodic_box = get_neighborhood_search(system, semi).periodic_box
-        coordinates = PointNeighbors.periodic_coords(active_coordinates(u, system),
-                                                     periodic_box)
-
-        x = collect(coordinates[1, :])
-        y = collect(coordinates[2, :])
-
-        particle_spacing = particle_spacings[i]
-        if particle_spacing < 0
-            particle_spacing = 0.0
-        end
-
-        x_min, x_max = extrema(x)
-        y_min, y_max = extrema(y)
-
-        # Add one particle radius around the center of the particles
-        # to obtain the domain size.
-        x_min -= 0.5particle_spacing
-        x_max += 0.5particle_spacing
-        y_min -= 0.5particle_spacing
-        y_max += 0.5particle_spacing
+    # Move data to the CPU if on the GPU.
+    # `transfer2cpu` also validates that all systems are in `semi.systems`, even on the CPU.
+    v_ode_, u_ode_, semi_ = transfer2cpu(v_ode, u_ode, semi)
 
-        return (; x, y, x_min, x_max, y_min, y_max, particle_spacing,
-                label=timer_name(system))
+    # Find the minimum and maximum coordinates across all systems,
+    # and then pass it to the next recipe.
+    systems_data = map(enumerate(semi_.systems)) do (i, system)
+        get_system_plot_data(u_ode_, system, semi_, particle_spacings[i])
     end
 
     # Pass the semidiscretization and the collected data to the next recipe
-    return (semi, systems_data...)
+    return (semi_, systems_data...)
 end
 
 function particle_spacings(semi::Semidiscretization)
@@ -93,7 +105,8 @@ RecipesBase.@recipe function f(initial_conditions::InitialCondition...)
     return (first(initial_conditions), ics...)
 end
 
-RecipesBase.@recipe function f(::Union{InitialCondition, Semidiscretization},
+RecipesBase.@recipe function f(::Union{AbstractSystem, InitialCondition,
+                                       Semidiscretization},
                                data...; size=(600, 400), xlims=(Inf, Inf), ylims=(Inf, Inf))
     # `data` is a tuple of named tuples, passed from the recipe above.
     # Each named tuple contains coordinates and metadata for a system or initial condition.