In this tutorial, we will guide you through the general structure of simulation files. We will set up a simulation similar to the example simulation examples/fluid/hydrostatic_water_column_2d.jl, which is one of our simplest example simulations. In the second part of this tutorial, we will show how to replace components of TrixiParticles.jl by custom implementations from within a simulation file, without ever cloning the repository.
For different setups and physics, have a look at our other example files.
First, we import TrixiParticles.jl and OrdinaryDiffEq.jl, which we will use at the very end for the time integration.
using TrixiParticles
+using OrdinaryDiffEq
Now, we define the particle spacing, which is our numerical resolution. We usually call the variable fluid_particle_spacing, so that we can easily change the resolution of an example file by overwriting this variable with trixi_include. In 2D, the number of particles will grow quadratically, in 3D cubically with the spacing.
We also set the number of boundary layers, which need to be sufficiently large, depending on the smoothing kernel and smoothing length, so that the compact support of the smoothing kernel is fully sampled with particles for a fluid particle close to a boundary. In particular, we require boundary_layers >= compact_support. The value for the compact support for each kernel can be found in the smoothing kernel overview.
fluid_particle_spacing = 0.05
+boundary_layers = 3
We want to simulate a water column resting under hydrostatic pressure inside a rectangular tank. 
First, we define the physical parameters gravitational acceleration, simulation time, initial fluid size, tank size and fluid density.
gravity = 9.81
tspan = (0.0, 1.0)
initial_fluid_size = (1.0, 0.9)
tank_size = (1.0, 1.0)
-fluid_density = 1000.0
In order to have the initial particle mass and density correspond to the hydrostatic pressure gradient, we need to define a state equation, which relates the fluid density to pressure. Note that we could also skip this part here and define the state equation later when we define the fluid system, but then the fluid would be initialized with constant density, which would cause it to oscillate under gravity.
sound_speed = 10.0
+fluid_density = 1000.0
In order to have the initial particle mass and density correspond to the hydrostatic pressure gradient, we need to define a state equation, which relates the fluid density to pressure. Note that we could also skip this part here and define the state equation later when we define the fluid system, but then the fluid would be initialized with constant density, which would cause it to oscillate under gravity.
sound_speed = 10.0
state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
exponent=7)
The speed of sound here is numerical and not physical. We artificially lower the speed of sound, since the physical speed of sound in water would lead to prohibitively small time steps. The speed of sound in Weakly Compressible SPH should be chosen as small as possible for numerical efficiency, but large enough to limit density fluctuations to about 1%.
TrixiParticles.jl requires the initial particle positions and quantities in form of an InitialCondition. Instead of manually defining particle positions, you can work with our pre-defined setups. Among others, we provide setups for rectangular shapes, circles, and spheres. Initial conditions can also be combined with common set operations. See this page for a list of pre-defined setups and details on set operations on initial conditions.
Here, we use the RectangularTank setup, which generates a rectangular fluid inside a rectangular tank, and supports a hydrostatic pressure gradient by passing a gravitational acceleration and a state equation (see above).
tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size,
fluid_density, n_layers=boundary_layers,
- acceleration=(0.0, -gravity), state_equation=state_equation)
To model the water column, we use the Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH) method. This method requires a smoothing kernel and a corresponding smoothing length, which should be chosen in relation to the particle spacing.
smoothing_length = 1.2 * fluid_particle_spacing
-smoothing_kernel = SchoenbergCubicSplineKernel{2}()
You can find an overview over smoothing kernels and corresponding smoothing lengths here.
For stability, we need numerical dissipation in form of an artificial viscosity term.
viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
We choose the parameters as small as possible to avoid visible viscosity, but as large as possible to stabilize the simulation.
The WCSPH method can either compute the particle density by a kernel summation over all neighboring particles (see SummationDensity) or by making the particle density a variable in the ODE system and integrating it over time. We choose the latter approach here by using the density calculator ContinuityDensity.
fluid_density_calculator = ContinuityDensity()
+ acceleration=(0.0, -gravity), state_equation=state_equation)
To model the water column, we use the Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH) method. This method requires a smoothing kernel and a corresponding smoothing length, which should be chosen in relation to the particle spacing.
smoothing_length = 1.2 * fluid_particle_spacing
+smoothing_kernel = SchoenbergCubicSplineKernel{2}()
You can find an overview over smoothing kernels and corresponding smoothing lengths here.
For stability, we need numerical dissipation in form of an artificial viscosity term. Other viscosity models offer a physical approach based on the kinematic viscosity of the fluid.
viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
We choose the parameters as small as possible to avoid non-physical behavior, but as large as possible to stabilize the simulation.
The WCSPH method can either compute the particle density directly with a kernel summation over all neighboring particles (see SummationDensity) or by making the particle density a variable in the ODE system and integrating its change over time. We choose the latter approach here by using the density calculator ContinuityDensity, which is more efficient and handles free surfaces without the need for additional correction terms.
fluid_density_calculator = ContinuityDensity()
fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
state_equation, smoothing_kernel,
smoothing_length, viscosity=viscosity,
- acceleration=(0.0, -gravity))
In order to define a boundary system, we first have to choose a boundary model, which defines how the fluid interacts with boundary particles. We will use the BoundaryModelDummyParticles with AdamiPressureExtrapolation, which generally produces the best results of the implemented methods. See here for a comprehensive overview over boundary models.
boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
+ acceleration=(0.0, -gravity))
To model the boundary, we use particle-based boundary conditions, in which particles are sampled in the boundary that interact with the fluid particles to avoid penetration. In order to define a boundary system, we first have to choose a boundary model, which defines how the fluid interacts with boundary particles. We will use the BoundaryModelDummyParticles with AdamiPressureExtrapolation. See here for a comprehensive overview over boundary models.
boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
state_equation=state_equation,
AdamiPressureExtrapolation(),
smoothing_kernel, smoothing_length)
boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
The key component of every simulation is the Semidiscretization, which couples all systems of the simulation. All methods in TrixiParticles.jl are semidiscretizations, which discretize the equations in time to provide an ordinary differential equation that still has to be solved in time. By providing a simulation time span, we can call semidiscretize, which returns an ODEProblem that can be solved with a time integration method.
semi = Semidiscretization(fluid_system, boundary_system)
-ode = semidiscretize(semi, tspan)
We use the methods provided by OrdinaryDiffEq.jl, but note that other packages or custom implementations can also be used.
OrdinaryDiffEq.jl supports callbacks, which are executed during the simulation. For this simulation, we use the InfoCallback, which prints information about the simulation setup at the beginning of the simulation, information about the current simulation time and runtime during the simulation, and a performance summary at the end of the simulation. We also want to save the current solution in regular intervals in terms of simulation time as VTK, so that we can look at the solution in ParaView. The SolutionSavingCallback provides this functionality. To pass the callbacks to OrdinaryDiffEq.jl, we have to bundle them into a CallbackSet.
info_callback = InfoCallback(interval=50)
+ode = semidiscretize(semi, tspan)
We use the methods provided by OrdinaryDiffEq.jl, but note that other packages or custom implementations can also be used.
OrdinaryDiffEq.jl supports callbacks, which are executed during the simulation. For this simulation, we use the InfoCallback, which prints information about the simulation setup at the beginning of the simulation, information about the current simulation time and runtime during the simulation, and a performance summary at the end of the simulation. We also want to save the current solution in regular intervals in terms of simulation time as VTK, so that we can look at the solution in ParaView. The SolutionSavingCallback provides this functionality. To pass the callbacks to OrdinaryDiffEq.jl, we have to bundle them into a CallbackSet.
info_callback = InfoCallback(interval=50)
saving_callback = SolutionSavingCallback(dt=0.02)
-callbacks = CallbackSet(info_callback, saving_callback)
Finally, we can start the simulation by solving the ODEProblem. We use the method RDPK3SpFSAL35 of OrdinaryDiffEq.jl, which is a Runge-Kutta method with automatic (error based) time step size control. This method is usually a good choice for prototyping, since we do not have to worry about choosing a stable step size and can just run the simulation. For better performance, it might be beneficial to tweak the tolerances of this method or choose a different method that is more efficient for the respective simulation. You can find both approaches in our example files. Here, we just use the method with the default parameters, and only disable save_everystep to avoid expensive saving of the solution in every time step.
sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks);
See Visualization for how to visualize the solution.
If we would like to use an implementation of a component that is not available in TrixiParticles.jl, we can implement it ourselves within the simulation file, without ever cloning the TrixiParticles.jl repository. A good starting point is to check out the available implementations in TrixiParticles.jl, then copy the relevant functions to the simulation file and modify them as needed.
To implement a custom smoothing kernel, we define a struct extending TrixiParticles.SmoothingKernel. This abstract struct has a type parameter for the number of dimensions, which we set to 2 in this case.
struct MyGaussianKernel <: TrixiParticles.SmoothingKernel{2} end
This kernel is going to be an implementation of the Gaussian kernel with a cutoff for compact support. Note that the same kernel in a more optimized version is already implemented in TrixiParticles.jl as GaussianKernel.
In order to use our new kernel, we have to define three functions. TrixiParticles.kernel, which is the kernel function itself, TrixiParticles.kernel_deriv, which is the derivative of the kernel function, and TrixiParticles.compact_support, which defines the compact support of the kernel in relation to the smoothing length. The latter is relevant for determining the search radius of the neighborhood search.
function TrixiParticles.kernel(kernel::MyGaussianKernel, r, h)
+callbacks = CallbackSet(info_callback, saving_callback)
Finally, we can start the simulation by solving the ODEProblem. We use the method RDPK3SpFSAL35 of OrdinaryDiffEq.jl, which is a Runge-Kutta method with automatic (error based) time step size control. This method is usually a good choice for prototyping, since we do not have to worry about choosing a stable step size and can just run the simulation. For better performance, it might be beneficial to tweak the tolerances of this method or choose a different method that is more efficient for the respective simulation. You can find both approaches in our example files. Here, we just use the method with the default parameters, and only disable save_everystep to avoid expensive saving of the solution in every time step.
sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks);
+████████╗██████╗ ██╗██╗ ██╗██╗██████╗ █████╗ ██████╗ ████████╗██╗ ██████╗██╗ ███████╗███████╗
+╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║██╔══██╗██╔══██╗██╔══██╗╚══██╔══╝██║██╔════╝██║ ██╔════╝██╔════╝
+ ██║ ██████╔╝██║ ╚███╔╝ ██║██████╔╝███████║██████╔╝ ██║ ██║██║ ██║ █████╗ ███████╗
+ ██║ ██╔══██╗██║ ██╔██╗ ██║██╔═══╝ ██╔══██║██╔══██╗ ██║ ██║██║ ██║ ██╔══╝ ╚════██║
+ ██║ ██║ ██║██║██╔╝ ██╗██║██║ ██║ ██║██║ ██║ ██║ ██║╚██████╗███████╗███████╗███████║
+ ╚═╝ ╚═╝ ╚═╝╚═╝╚═╝ ╚═╝╚═╝╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝ ╚═╝ ╚═╝ ╚═════╝╚══════╝╚══════╝╚══════╝
+
+
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ Semidiscretization │
+│ ══════════════════ │
+│ #spatial dimensions: ………………………… 2 │
+│ #systems: ……………………………………………………… 2 │
+│ neighborhood search: ………………………… GridNeighborhoodSearch │
+│ total #particles: ………………………………… 636 │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ WeaklyCompressibleSPHSystem{2} │
+│ ══════════════════════════════ │
+│ #particles: ………………………………………………… 360 │
+│ density calculator: …………………………… ContinuityDensity │
+│ correction method: ……………………………… Nothing │
+│ state equation: ……………………………………… StateEquationCole │
+│ smoothing kernel: ………………………………… SchoenbergCubicSplineKernel │
+│ viscosity: …………………………………………………… ArtificialViscosityMonaghan{Float64}(0.02, 0.0, 0.01) │
+│ density diffusion: ……………………………… nothing │
+│ surface tension: …………………………………… nothing │
+│ acceleration: …………………………………………… [0.0, -9.81] │
+│ source terms: …………………………………………… Nothing │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ BoundarySPHSystem{2} │
+│ ════════════════════ │
+│ #particles: ………………………………………………… 276 │
+│ boundary model: ……………………………………… BoundaryModelDummyParticles(AdamiPressureExtrapolation, Nothing) │
+│ movement function: ……………………………… nothing │
+│ adhesion coefficient: ……………………… 0.0 │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ SolutionSavingCallback │
+│ ══════════════════════ │
+│ dt: ……………………………………………………………………… 0.02 │
+│ custom quantities: ……………………………… nothing │
+│ save initial solution: …………………… yes │
+│ save final solution: ………………………… yes │
+│ output directory: ………………………………… /home/runner/work/TrixiParticles…les.jl/docs/build/tutorials/out │
+│ prefix: …………………………………………………………… │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ Time integration │
+│ ════════════════ │
+│ Start time: ………………………………………………… 0.0 │
+│ Final time: ………………………………………………… 1.0 │
+│ time integrator: …………………………………… RDPK3SpFSAL35 │
+│ adaptive: ……………………………………………………… true │
+│ abstol: …………………………………………………………… 1.0e-6 │
+│ reltol: …………………………………………………………… 0.001 │
+│ controller: ………………………………………………… PIDController(beta=[0.7, -0.23, …iter=default_dt_factor_limiter) │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ Environment information │
+│ ═══════════════════════ │
+│ #threads: ……………………………………………………… 1 │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+
+#timesteps: 50 │ Δt: 1.4965e-03 │ sim. time: 1.8236e-01 (18.236%) │ run time: 1.2462e-01 s
+#timesteps: 100 │ Δt: 5.3979e-03 │ sim. time: 4.2540e-01 (42.540%) │ run time: 2.6655e-01 s
+#timesteps: 150 │ Δt: 5.7273e-03 │ sim. time: 6.7662e-01 (67.662%) │ run time: 3.9337e-01 s
+#timesteps: 200 │ Δt: 4.8634e-03 │ sim. time: 9.1219e-01 (91.219%) │ run time: 5.2495e-01 s
+────────────────────────────────────────────────────────────────────────────────────────────────────
+Trixi simulation finished. Final time: 1.0 Time steps: 222 (accepted), 222 (total)
+────────────────────────────────────────────────────────────────────────────────────────────────────
+
+────────────────────────────────────────────────────────────────────────────────
+ TrixiParticles.jl Time Allocations
+ ─────────────────────── ────────────────────────
+ Tot / % measured: 576ms / 98.0% 95.4MiB / 98.4%
+
+Section ncalls time %tot avg alloc %tot avg
+────────────────────────────────────────────────────────────────────────────────
+kick! 1.11k 358ms 63.5% 322μs 3.77MiB 4.0% 3.47KiB
+ system interaction 1.11k 286ms 50.7% 257μs 1.05MiB 1.1% 985B
+ fluid1-fluid1 1.11k 226ms 40.0% 203μs 0.00B 0.0% 0.00B
+ fluid1-boundary2 1.11k 57.0ms 10.1% 51.2μs 0.00B 0.0% 0.00B
+ ~system intera... 1.11k 3.07ms 0.5% 2.76μs 1.05MiB 1.1% 985B
+ boundary2-fluid1 1.11k 32.7μs 0.0% 29.4ns 0.00B 0.0% 0.00B
+ boundary2-boun... 1.11k 32.6μs 0.0% 29.2ns 0.00B 0.0% 0.00B
+ update systems a... 1.11k 71.4ms 12.7% 64.1μs 2.72MiB 2.9% 2.50KiB
+ compute bounda... 1.11k 48.8ms 8.6% 43.8μs 0.00B 0.0% 0.00B
+ update nhs 1.11k 11.4ms 2.0% 10.2μs 2.72MiB 2.9% 2.50KiB
+ inverse state ... 1.11k 6.14ms 1.1% 5.52μs 0.00B 0.0% 0.00B
+ ~update system... 1.11k 5.06ms 0.9% 4.55μs 1.55KiB 0.0% 1.42B
+ update density... 1.11k 32.9μs 0.0% 29.6ns 0.00B 0.0% 0.00B
+ source terms 1.11k 425μs 0.1% 381ns 0.00B 0.0% 0.00B
+ ~kick!~ 1.11k 407μs 0.1% 366ns 1.55KiB 0.0% 1.42B
+ reset ∂v/∂t 1.11k 151μs 0.0% 135ns 0.00B 0.0% 0.00B
+save solution 50 205ms 36.3% 4.09ms 90.0MiB 96.0% 1.80MiB
+ write to vtk 100 176ms 31.3% 1.76ms 86.7MiB 92.5% 888KiB
+ ~save solution~ 50 25.2ms 4.5% 504μs 3.15MiB 3.4% 64.6KiB
+ update systems 50 3.09ms 0.5% 61.8μs 127KiB 0.1% 2.53KiB
+ compute bounda... 50 2.10ms 0.4% 41.9μs 0.00B 0.0% 0.00B
+ update nhs 50 494μs 0.1% 9.89μs 125KiB 0.1% 2.50KiB
+ inverse state ... 50 269μs 0.0% 5.38μs 0.00B 0.0% 0.00B
+ ~update systems~ 50 229μs 0.0% 4.58μs 1.55KiB 0.0% 31.7B
+ update density... 50 1.50μs 0.0% 30.0ns 0.00B 0.0% 0.00B
+drift! 1.11k 873μs 0.2% 784ns 976B 0.0% 0.88B
+ velocity 1.11k 526μs 0.1% 473ns 0.00B 0.0% 0.00B
+ ~drift!~ 1.11k 221μs 0.0% 198ns 976B 0.0% 0.88B
+ reset ∂u/∂t 1.11k 126μs 0.0% 113ns 0.00B 0.0% 0.00B
+compute boundary p... 1 59.3μs 0.0% 59.3μs 0.00B 0.0% 0.00B
+update nhs 1 54.9μs 0.0% 54.9μs 2.50KiB 0.0% 2.50KiB
+inverse state equa... 1 7.96μs 0.0% 7.96μs 0.00B 0.0% 0.00B
+update density dif... 1 29.0ns 0.0% 29.0ns 0.00B 0.0% 0.00B
+────────────────────────────────────────────────────────────────────────────────
See Visualization for how to visualize the solution. For the simplest visualization, we can use Plots.jl:
using Plots
+plot(sol)
GKS: cannot open display - headless operation mode active
If we would like to use an implementation of a component that is not available in TrixiParticles.jl, we can implement it ourselves within the simulation file, without ever cloning the TrixiParticles.jl repository. A good starting point is to check out the available implementations in TrixiParticles.jl, then copy the relevant functions to the simulation file and modify them as needed.
To implement a custom smoothing kernel, we define a struct extending TrixiParticles.SmoothingKernel. This abstract struct has a type parameter for the number of dimensions, which we set to 2 in this case.
struct MyGaussianKernel <: TrixiParticles.SmoothingKernel{2} end
This kernel is going to be an implementation of the Gaussian kernel with a cutoff for compact support. Note that the same kernel in a more optimized version is already implemented in TrixiParticles.jl as GaussianKernel.
In order to use our new kernel, we have to define three functions. TrixiParticles.kernel, which is the kernel function itself, TrixiParticles.kernel_deriv, which is the derivative of the kernel function, and TrixiParticles.compact_support, which defines the compact support of the kernel in relation to the smoothing length. The latter is relevant for determining the search radius of the neighborhood search.
function TrixiParticles.kernel(kernel::MyGaussianKernel, r, h)
q = r / h
if q < 2
@@ -42,5 +161,140 @@
return 0.0
end
-TrixiParticles.compact_support(::MyGaussianKernel, h) = 3 * h
This is all we need to use our custom kernel implementation in a simulation. We only need to replace the definition above by
smoothing_kernel = MyGaussianKernel()
and run the simulation file again.
In order to use our kernel in a pre-defined example file, we can use the function trixi_include to replace the definition of the variable smoothing_kernel. The following will run the example simulation examples/fluid/hydrostatic_water_column_2d.jl with our custom kernel.
julia> trixi_include(joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"),
- smoothing_kernel=MyGaussianKernel());