Drafty draft

NEWS.md

@@ -4,25 +4,34 @@ TrixiParticles.jl follows the interpretation of
 [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
-## Version 0.5.4
-
-### API Changes
-
-- Replaced the experimental `MechanicalWorkCalculatorCallback` with
-  `MechanicalWorkCalculator`, which can be passed as a custom quantity to
-  `PostprocessCallback` (#1228).
-
-### Important Bugfixes
-
-- Fixed mathematical inconsistencies in the SPH documentation, including incorrect
-  formulas, inconsistent force-vs-acceleration notation, and wrong LaTeX text-mode
-  commands (#1086).
-- Fixed restarting with EDAC from solution objects (#1213) and from VTK files (#1297).
-- Fixed the custom quantities `kinetic_energy`, `total_mass`, `max_pressure`, `min_pressure`,
-  `avg_pressure`, `max_density`, `min_density` and `avg_density` to only take active particles
-  into account (#1184).
-
-## Version 0.5.3
+## Version 0.5.4
+
+### API Changes
+
+- Replaced the experimental `MechanicalWorkCalculatorCallback` with
+  `MechanicalWorkCalculator`, which can be passed as a custom quantity to
+  `PostprocessCallback` (#1228).
+
+### Important Bugfixes
+
+- Fixed mathematical inconsistencies in the SPH documentation, including incorrect
+  formulas, inconsistent force-vs-acceleration notation, and wrong LaTeX text-mode
+  commands (#1086).
+- Fixed restarting with EDAC from solution objects (#1213) and from VTK files (#1297).
+- Fixed the custom quantities `kinetic_energy`, `total_mass`, `max_pressure`, `min_pressure`,
+  `avg_pressure`, `max_density`, `min_density` and `avg_density` to only take active particles
+  into account (#1184).
+
+- Hardened surface tension model configuration by validating coefficients and surface-normal
+  thresholds, avoiding unnecessary normal allocation for `CohesionForceAkinci`, and stabilizing
+  Akinci cohesion and adhesion kernels across floating-point scales.
+
+### Features
+
+- Added an optional Makie recipe for rendering two- and three-dimensional particle systems
+  with `plot`, `plot!`, `trixi2makie`, and `trixi2makie!`.
+
+## Version 0.5.3
 
 ### Features
 

Project.toml

@@ -35,6 +35,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d"
 OrdinaryDiffEqSymplecticRK = "fa646aed-7ef9-47eb-84c4-9443fc8cbfa8"
@@ -43,6 +44,7 @@ Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 [extensions]
 TrixiParticlesOrdinaryDiffEqSymplecticRKExt = ["OrdinaryDiffEqSymplecticRK", "OrdinaryDiffEqCore"]
 TrixiParticlesCUDAExt = "CUDA"
+TrixiParticlesMakieExt = "Makie"
 
 [compat]
 Accessors = "0.1.43"
@@ -60,6 +62,7 @@ GPUArraysCore = "0.2"
 JSON = "1"
 KernelAbstractions = "0.9"
 LinearAlgebra = "1"
+Makie = "0.24"
 OrdinaryDiffEqLowStorageRK = "3"
 OrdinaryDiffEqCore = "4"
 OrdinaryDiffEqSymplecticRK = "2"

README.md

@@ -62,24 +62,24 @@ It offers intuitive configuration, robust pre- and post-processing, and vendor-a
 We provide several example simulation setups in the `examples` folder (which can be accessed from Julia via `examples_dir()`).
 
 <table align="center" border="0">
-  <tr>
-  </tr>
   <tr>
     <td align="center">
-      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/683e9363-5705-49cc-9a5c-3b47d73ea4b8" style="width: 80% !important;"/><br><figcaption>2D Dam Break</figcaption>
+      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/683e9363-5705-49cc-9a5c-3b47d73ea4b8" style="width: 80% !important;"/><br />
+      <figcaption>2D Dam Break</figcaption>
     </td>
     <td align="center">
-      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/c10faddf-0400-47c9-b225-f5d286a8ecb8" style="width: 80% !important;"/><br><figcaption>Moving Wall</figcaption>
+      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/c10faddf-0400-47c9-b225-f5d286a8ecb8" style="width: 80% !important;"/><br />
+      <figcaption>Moving Wall</figcaption>
     </td>
   </tr>
-  <tr>
-  </tr>
   <tr>
     <td align="center">
-      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/e05ace63-e330-441a-a391-eda3d2764074" style="width: 80% !important;"/><br><figcaption>Oscillating Beam</figcaption>
+      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/e05ace63-e330-441a-a391-eda3d2764074" style="width: 80% !important;"/><br />
+      <figcaption>Oscillating Beam</figcaption>
     </td>
     <td align="center">
-      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/ada0d554-e0ba-44ed-923d-2b77ef252258" style="width: 80% !important;"/><br><figcaption>Dam Break with Elastic Plate</figcaption>
+      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/ada0d554-e0ba-44ed-923d-2b77ef252258" style="width: 80% !important;"/><br />
+      <figcaption>Dam Break with Elastic Plate</figcaption>
     </td>
   </tr>
 </table>

docs/literate/src/tut_2d_geometry.jl

@@ -0,0 +1,150 @@
+# # [Setting up a 2D simulation from geometry files](@id tut_2d_geometry)
+
+# In this tutorial, we build two 2D setups from geometry files:
+# 1. a curved pipe, where one geometry file defines the outer wall and another
+#    defines the channel cut out of it,
+# 2. a dam-break basin with a coastline profile, where one geometry file defines
+#    the coastline wall and the seawall on the right.
+#
+# We use 2D geometry formats such as `.asc` or `.dxf` for 2D setups.
+# STL files describe surfaces and are therefore better suited to 3D setups.
+
+# First, we import TrixiParticles.jl together with
+# `OrdinaryDiffEqLowStorageRK` of
+# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
+# and [Plots.jl](https://docs.juliaplots.org/stable/).
+using TrixiParticles
+using OrdinaryDiffEqLowStorageRK
+using Plots
+
+# ## Resolution
+
+# We use the same particle spacing for the fluid and for the wall geometries.
+particle_spacing = 0.03
+fluid_density = 1000.0
+gravity = 9.81
+sound_speed = 10.0
+state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                   exponent=7)
+nothing # hide
+
+# ## Loading 2D geometry files
+
+# The following helper loads a closed 2D geometry file and samples particles in its interior:
+# 1. load the polygon with [`load_geometry`](@ref),
+# 2. fill the polygon with [`ComplexShape`](@ref).
+#
+# This creates a filled 2D solid region rather than particles only along the polygon edges.
+function solid_from_geometry_file(file; particle_spacing, density)
+    geometry = load_geometry(file)
+    solid = ComplexShape(geometry; particle_spacing, density,
+                         grid_offset=0.5particle_spacing)
+
+    return (; geometry, solid)
+end
+
+# ## A curved pipe from two filled geometries
+
+# The pipe wall is an L-shaped solid region with a channel cut out of it:
+# 1. one geometry file describes the outer pipe envelope,
+# 2. one geometry file describes the empty channel,
+# 3. `setdiff` subtracts the channel from the solid envelope.
+pipe_outer_file = pkgdir(TrixiParticles, "examples", "preprocessing", "data",
+                         "curved_pipe_outer_2d.asc")
+pipe_channel_file = pkgdir(TrixiParticles, "examples", "preprocessing", "data",
+                           "curved_pipe_channel_2d.asc")
+
+pipe_outer = solid_from_geometry_file(pipe_outer_file; particle_spacing,
+                                      density=fluid_density)
+pipe_channel = load_geometry(pipe_channel_file)
+
+pipe_setup = (; wall=setdiff(pipe_outer.solid, pipe_channel),
+              outer_geometry=pipe_outer.geometry,
+              channel_geometry=pipe_channel)
+
+# ## A dam-break basin with a coastline profile
+
+# In the second setup, one 2D geometry file defines the coastline wall:
+# the beach profile, a finite wall thickness below it, and the seawall on the right.
+coast_file = pkgdir(TrixiParticles, "examples", "preprocessing", "data",
+                    "coastline_profile_2d.asc")
+coast = solid_from_geometry_file(coast_file; particle_spacing, density=fluid_density)
+
+# The geometry file defines the coastline bed and right wall as a solid region.
+# We add the left wall as a rectangular particle block and place a rectangular
+# dam-break water column beside it.
+left_wall = RectangularShape(particle_spacing, (5, 50), (0.0, -0.12),
+                             density=fluid_density)
+reservoir = RectangularShape(particle_spacing, (28, 42), (0.15, 0.03),
+                             acceleration=(0.0, -gravity),
+                             state_equation=state_equation)
+coast_setup = (; geometry=coast.geometry,
+               wall=union(coast.solid, left_wall),
+               fluid=setdiff(reservoir, coast.geometry))
+
+p_pipe = plot(pipe_setup.wall, label="wall", title="Curved pipe",
+              markerstrokewidth=0, markersize=4)
+plot!(p_pipe, showaxis=false, aspect_ratio=:equal,
+      xlims=(-0.03, 1.23), ylims=(-0.03, 1.23))
+
+p_coast = plot(coast_setup.fluid, coast_setup.wall,
+               labels=["fluid" "wall"], title="Coastline dam break",
+               markerstrokewidth=0, markersize=3)
+plot!(p_coast, showaxis=false, aspect_ratio=:equal,
+      xlims=(0.0, 2.75), ylims=(-0.15, 1.35))
+
+plot(p_pipe, p_coast, layout=(1, 2), size=(900, 360))
+savefig("tut_2d_geometry_plot.png"); # hide
+# ![2D geometry based initial conditions](tut_2d_geometry_plot.png)
+
+# ## Building the simulation systems
+
+# We continue with the coastline setup. The remaining steps are the same as for
+# other 2D simulations.
+setup = coast_setup
+tspan = (0.0, 0.03)
+nothing # hide
+
+# We define the state equation, smoothing kernel, and viscosity for the
+# weakly compressible SPH simulation.
+smoothing_length = 1.2 * particle_spacing
+smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
+
+fluid_density_calculator = ContinuityDensity()
+density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
+
+fluid_system = WeaklyCompressibleSPHSystem(setup.fluid;
+                                           density_calculator=fluid_density_calculator,
+                                           state_equation, smoothing_kernel,
+                                           smoothing_length, viscosity=viscosity,
+                                           density_diffusion=density_diffusion,
+                                           acceleration=(0.0, -gravity))
+nothing # hide
+
+# For the wall, we reuse the combined solid wall particles created above. The high-level
+# constructor obtains the smoothing kernel, smoothing length, and state equation from the fluid.
+boundary_model = BoundaryModelDummyParticles(setup.wall; fluid_system)
+boundary_system = WallBoundarySystem(setup.wall, boundary_model)
+nothing # hide
+
+# ## Semidiscretization
+
+# We construct the [`Semidiscretization`](@ref TrixiParticles.Semidiscretization)
+# from the fluid and boundary systems.
+semi = Semidiscretization(fluid_system, boundary_system)
+ode = semidiscretize(semi, tspan)
+nothing # hide
+
+# ## Time integration
+
+# We can now solve the problem. An [`InfoCallback`](@ref) prints progress during
+# the simulation.
+callbacks = CallbackSet(InfoCallback(interval=10))
+nothing # hide
+
+sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks) #!md
+
+# For more accurate body-fitted particles around sharper features, you can also
+# apply the [particle packing workflow](@ref tut_packing) to the 2D geometry files
+# before starting the simulation.

docs/literate/src/tut_packing.jl

@@ -75,8 +75,8 @@ plot!(right_margin=5Plots.mm) #hide
 # ## Creating an initial configuration of boundary particles
 
 # To create the initial configuration of the boundary particles,
-# we use the sampled points of the SDF whose signed distance lies between 0
-# and `boundary_thickness`.
+# we use the sampled points of the SDF whose signed distance lies between the
+# geometry offset implied by `place_on_shell` and `boundary_thickness`.
 # Here, we need to specify the `density` of the boundary particles.
 # As an example, we choose `1.0` for all particles.
 # This gives us an [`InitialCondition`](@ref InitialCondition) for the boundary particles.
@@ -125,7 +125,8 @@ plot!(geometry, linestyle=:dash, label=nothing, showaxis=false, color=:black,
 # ## Particle packing
 
 # In the following, we will essentially follow the same steps described in the fluid tutorials.
-# That means we will generate systems that are then passed to the [`Semidiscretization`](@ref).
+# That means we will generate systems that are then passed to the
+# [`Semidiscretization`](@ref TrixiParticles.Semidiscretization).
 # The difference from a typical physical simulation is that we use [`ParticlePackingSystem`](@ref),
 # which does not represent any physical law. Instead, we only use the simulation framework to time-integrate
 # the packing process.
@@ -211,7 +212,7 @@ plot!(geometry, seriestype=:path, color=:black, label=nothing, linewidth=2)
 boundary_system = ParticlePackingSystem(boundary_sampled; is_boundary=true,
                                         smoothing_kernel, smoothing_length,
                                         boundary_compress_factor=0.7, signed_distance_field,
-                                        background_pressure)
+                                        boundary_thickness, background_pressure)
 
 # We can now couple the boundary system with the interior system:
 semi = Semidiscretization(packing_system, boundary_system)
@@ -251,10 +252,12 @@ fixed_system = ParticlePackingSystem(packed_ic; smoothing_kernel, smoothing_leng
 
 # Now we define a rectangular domain that we want to pack.
 # In practice, you could create any `InitialCondition` that encloses your complex geometry.
-tank_domain = RectangularTank(particle_spacing, (4, 4), (0, 0), min_coordinates=(-1, -2),
-                              density)
+domain_size = (4, 4)
+n_particles_per_dimension = round.(Int, domain_size ./ particle_spacing)
+tank_domain = RectangularShape(particle_spacing, n_particles_per_dimension, (-1, -2);
+                               density)
 
-sampled_outer_domain = setdiff(tank_domain.fluid, packed_ic)
+sampled_outer_domain = setdiff(tank_domain, packed_ic)
 
 # If we plot these two `InitialCondition`s, we can see
 # that the geometry interface is not properly represented yet.