diff --git a/NEWS.md b/NEWS.md index 5e9256aecd..fc2a9e1082 100644 --- a/NEWS.md +++ b/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 diff --git a/docs/literate/src/tut_visualization.jl b/docs/literate/src/tut_visualization.jl new file mode 100644 index 0000000000..296482e09e --- /dev/null +++ b/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 diff --git a/docs/make.jl b/docs/make.jl index ee7d9f9c8e..656513af69 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 4225faf486..5d0614aca3 100644 --- a/docs/src/tutorial.md +++ b/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 +Velocity from a coarse vortex-street simulation interpolated onto a regular grid +``` + +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 diff --git a/docs/src/visualization.md b/docs/src/visualization.md index 9dda94bfba..effbc90f57 100644 --- a/docs/src/visualization.md +++ b/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` diff --git a/src/visualization/recipes_plots.jl b/src/visualization/recipes_plots.jl index c1378497a1..8ce2a13491 100644 --- a/src/visualization/recipes_plots.jl +++ b/src/visualization/recipes_plots.jl @@ -13,14 +13,48 @@ 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, @@ -28,40 +62,18 @@ RecipesBase.@recipe function f(v_ode::AbstractArray, u_ode::AbstractArray, 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.