Globally time-reversible fluid simulations with smoothed particle hydrodynamics

by   Ondrej Kincl, et al.

This paper describes an energy-preserving and globally time-reversible code for weakly compressible smoothed particle hydrodynamics (SPH). We do not add any additional dynamics to the Monaghan's original SPH scheme at the level of ordinary differential equation, but we show how to discretize the equations by using a corrected expression for density and by invoking a symplectic integrator. Moreover, to achieve the global-in-time reversibility, we have to correct the initial state, implement a conservative fluid-wall interaction, and use the fixed-point arithmetic. Although the numerical scheme is reversible globally in time (solvable backwards in time while recovering the initial conditions), we observe thermalization of the particle velocities and growth of the Boltzmann entropy. In other words, when we do not see all the possible details, as in the Boltzmann entropy, which depends only on the one-particle distribution function, we observe the emergence of the second law of thermodynamics (irreversible behavior) from purely reversible dynamics.


page 4

page 8

page 12


Non-reversible sampling schemes on submanifolds

Calculating averages with respect to probability measures on submanifold...

An Exact Bitwise Reversible Integrator

At a fundamental level most physical equations are time reversible. In t...

An entropic method for discrete systems with Gibbs entropy

We consider general systems of ordinary differential equations with mono...

A Deep Dive into the Distribution Function: Understanding Phase Space Dynamics with Continuum Vlasov-Maxwell Simulations

In collisionless and weakly collisional plasmas, the particle distributi...

The Fictitious Domain Method Based on Navier Slip Boundary Condition for Simulation of Flow-Particle Interaction

In this article, we develop a least–squares/fictitious domain method for...

Numerical verification of the microscopic time reversibility of Newton's equations of motion: Fighting exponential divergence

Numerical solutions to Newton's equations of motion for chaotic self gra...

Explicit, time-reversible and symplectic integrator for Hamiltonians in isotropic uniformly curved geometries

The kinetic term of the N-body Hamiltonian system defined on the surface...

1 Introduction

The original SPH scheme, now called weakly compressible SPH (WCSPH), was developed in 1960s by R.Gingold and J. Monaghan [1]. The scheme used the leap-frog time integration and relied on an artificial term stabilizing the fluid in the presence of shocks. Many other stabilization strategies were conceived, including -SPH [2] and methods based on Riemann solvers [3]. These approaches often provide results in reasonable agreement with experimental data or other numerical methods. Unfortunately, they often come at the cost of increased complexity and necessity to fine-tune additional numerical parameters. Failure in such fine-tuning can result in unrealistic dissipation, or large energy oscillations unless a convenient limiter is used [4].

There are essentially two different ways, how pressure can be calculated in WCSPH [5]

. A commonly used approach is to update the density iteratively, using the velocities obtained from the discretized balance of momentum. Alternatively, one can compute density directly as an SPH interpolation of particle masses. The latter approach is less common, but it leads to a symplectic structure of the system of SPH differential equations, which makes the WCSPH discretization completely stable even in the absence of viscosity or stabilization. Moreover, the WCSPH numerical scheme becomes globally reversible in time when also the fixed-point arithmetic is used. For instance, the breaking dam simulation can be reversed to recover the initial conditions. The purpose of this paper is to show key advantages of the symplectic WCSPH compared to the standard approach.

However, as pointed out by D. Violeau, symplectic WCSPH shows unphysical behavior at free surfaces, where constant functions are not correctly interpolated. We suggest two ways how to correct this issue. The first method modifies the density by an appropriate integration constant, which is easy to implement. The second and more complex methods renormalizes the initial state by solving a certain nonlinear problem using Newton’s method. We also demonstrate that a careful implementation of the symplectic WCSPH scheme, using the corrected treatment of free surfaces and the fixed-point (as opposed to the more usual floating-point) arithmetic leads to bit-perfect global time reversibility. This can be considered an additional symmetry-preserving feature alongside the conservation of energy, mass, momentum and angular momentum.

Although the resulting code is energy-preserving and reversible, it can still be considered dissipative in the sense that the Boltzmann entropy of the SPH particles grows. In other words, we observe the emergence of the second law of thermodynamics from purely reversible dynamics.

We demonstrate the various versions of WCSPH on dam-break and Gresho vortex benchmarks.

All SPH codes used to make this paper are freely available within the new SmoothedParticles.jl package [6] written in the Julia programming language [7].

2 WCSPH for inviscid fluids

The standard weakly compressible spatial semi-discretization of Euler fluid equations in a uniform gravitational field is [8]:


where the pressure is a function of density via the baratropic formula,


Constant stands for the referential density, is the gravitational acceleration, and is the numerical speed of sound (typically chosen as ten times the characteristic flow speed ). It is well known that the system of ordinary differential equations (1) conserves the total energy in the form



Moreover, for , the momentum

and angular momentum

are conserved as well. Naturally, gravitational force can deliver some momentum to the fluid, while the equal opposite reaction of the fluid on Earth is considered grounded.

The WCSPH equations (1) represent a Hamiltonian system, which can be checked for instance using program [9]. But it is not a symplectic system because the mass density is treated as a state variable. If, on the other hand, the density would always be calculated from the current positions of the particles, then only the positions and momenta would be state variables and the system of equations would be symplectic. Symplecticity would make it easier to choose a geometric structure-preserving numerical integrator, for instance using the Verlet scheme.

A symplectic form of WCSPH is obtained by solving the ordinary differential equation for density (1a),


instead of updating iteratively by 1a. However, this closed expression for density (4) is rarely used in practice because it leads to questionable behavior near free surfaces [5]. This happens because free boundary particles are under-occupied and, according to (4), they have . Therefore, they are equipped with large internal energy and, affected by its negative gradient, they start to vibrate. A remedy is provided in the following section, which leads to the possibility to use the symplectic form of WCSPH.

3 Treatment of free surfaces

3.1 Solving the equation for with an appropriate integration constant

The issue of unphysical behavior at free surfaces can be corrected by choosing an appropriate integration constant when solving the differential equation for ,



Indeed, since WCSPH based on (1a) typically assumes zero internal forces at , it is actually (5) and not (4) that is equivalent to the standard WCSPH formulation, since (4) leads to non-zero density gradients near the boundary. Although schemes based on (1a) and (5) are completely equivalent at the ODE level, they become very much different once time discretization is considered. It appears that (1a) is currently used by most researchers, despite the advantages of Equation (5):

  1. Inferring from positions using Equation (5) does not accumulate the density error.

  2. The system of SPH evolution equations becomes symplectic (and not only Hamiltonian).

  3. Therefore, symplectic structure-preserving integrators can be employed, which leads to simulations globally reversible in time.

However, even if we use the correct closed formula for density (5) and use a symplectic integrator, are facing another obstacle. For problems involving free surfaces, integration constant is typically much bigger for boundary particles and, therefore, whenever such particle submerges into interior of the fluid body, it generates a small void pocket around itself. This unwanted numerical artifact is removed in the following section.

3.2 Initial state correction (ISC)

The numerical artifact that particles that are originally at the boundary are treated differently than the rest of the particles, even if the submerge into the in the interior of the fluid, is rooted in the difference in the constants at the beginning. Therefore, w´e shall remove that artifact by correcting the initial state.

A solution is to correct the initial positions by a vector field

such that

These non-linear equations can be solved by a strategy based on the Newton’s method, details of which follow. Linearization with respect to leads to

On the left hand side, we identify the discrete divergence operator

see [5] for details. This linear problem for unknowns is underdetermined, so shall be restricted to in the form of a discrete gradient
which can be also rewritten as

Sparse linear system (6) for and has to be solved iteratively, updating and reevaluating according to (4), until with satisfactory precision. This initial state correction (ISC) procedure has been tested in our numerical examples, and the procedure converges quadratically except for problems with too much symmetry (for instance rectangle with particles arranged in square grid nearly always diverges)111For these cases, we recommend to initially disrupt particle positions by a random noise.. The whole procedure may be then summarized as follows:

before ISC
after ISC
Figure 1: Applying ISC to particles in square arranged regularly. Density (4) is shown in color (red to blue). Initially, it varies due to boundary effect but becomes everywhere constant after 7 iterations of ISC (up to relative error less than ).

Initial state correction (ISC) set randomly for all repeat until sufficient precision is acquired: for all recompute neighbor list find for all solve linear problem

This computation is expensive, but is performed only once per simulation. In place of the linear system with block matrix, one could also solve directly for :

However, the product leads to a matrix which is much more dense. A matrix-free iterative method could be used without evaluating the product directly, is a future possibility.

Let us note that after ISC, all constant fields will be correctly interpolated, at least initially. Weak compressibility helps to preserve this trait approximately in a way that does not accumulate error provided (4) is used. In this sense, ISC can be considered an alternative to (zeroth order) operator renormalization with the benefit that it does not break the symmetry of smoothing kernel. Unfortunately, this technique has two drawbacks. Firstly, it may slightly deform geometry in an unwanted way. Secondly, adapting it to a concrete setting may be difficult — for instance, when an anti-clump term is included.

Both symplecticity of the system of equations and ISC improve the numerical solutions, but in order to obtain a globally reversible behavior we have to tackle a further problem caused at the walls of the simulation box, as in the following section.

3.3 Wall modelling

In WCSPH, walls can be modeled by a layer of dummy particles which behave as if they belonged to the fluid, except that their positions and velocities are fixed. In this approach, it is hoped that the forces preventing compression will deflect any incoming particle. This method, of course, violates balance of momentum, but this is plausible, since we can imagine that any momentum yielded by the fluid is grounded. More importantly though, dummy particles violate conservation of energy. To see this, consider a wave stopped by a wall. The fluid acts on the solid with equal and opposite force. Therefore, a wall particle is subject to acceleration, which should change its velocity after time step to

To preserve wall integrity, this velocity must be annulled, reducing the total energy by a small kinetic contribution. From this, we see that dummy particles are essentially dissipative.

One could argue that this is still physically reasonable because every fluid-solid collision creates a sound wave, leading to loss of energy. However, it is dubious whether such effect is correctly modeled by dummy particles. Moreover, it complicates efforts to make a strictly energy-conservative scheme.

An alternative approach uses a fluid-solid force inspired by the Lennard-Jones potential

by which a wall particle acts onto a fluid particle . Here,

and are numerical parameters [10]. Except for this interaction, wall particles are neither considered in the balance of mass (1b) nor in the calculation of density (5). Unlike the dummy particle approach, this solid wall modeling is conservative, with a modified energy


This makes it potentially advantageous for problems involving inviscid flows.

In our case, when we seek a globally reversible-in-time simulation, the conservativeness of the fluid-wall interaction becomes important. However, even if we use the correct close formula for density (1a), the initial state correction, and the conservative fluid-wall interaction, the symplectic Verlet scheme is not reversible globally in time because of the errors caused by the floating-point arithmetic. The following section contains the final ingredient necessary for that reversibility, the symplectic Verlet integrator with fixed-point arithmetic.

4 Symplectic integrator and fixed-point arithmetic

The symplectic character of evolution equations (1b) and (1c) is preserved in symplectic numerical schemes, for instance the classical Verlet scheme [11]. Consequence are for instance the long-time stability of the trajectories and conservation of of integrals of motion (energy, momentum, and angular momentum, among others). However, even application of the Verlet scheme in SPH does not lead to globally in-time reversible simulations due to the floating-point errors, but a remedy is in the fixed-point arithmetic, as we show below.

The usual discretization of the SPH equations using the Verlet scheme is

where is the acceleration composed of internal, gravitational, and wall forces. The scheme is of second order and symplectic, and therefore it conserves energy with error

that does not depend explicitly on time, provided that is sufficiently small [12]222Here, we rely on the fact that using the closed expression (5).. If we took as a separate state variable with its own evolution equation (1a) (usual WCSPH), symplecticity of the system of equations would be violated. The usual WCSPH equations are still Hamiltonian (generated by a Poisson bracket and a Hamiltonian), as checked by program [9], but they are not symplectic and they thus the (symplectic) Verlet scheme does not preserve their geometric structure, which causes energy to grow exponentially, unless some stabilization is used - see Figure 2.333In [8], it is suggested to correct this error by a mid-step extrapolation of density. However, it should be noted this approach adds additional inaccuracy to the scheme.

Figure 2: Energy growth in standard WCSPH (STD) and symplectic WCSPH (SYM) in dambreak simulation with zero viscosity and parameters from Table 1. Without a stabilization of some sort, STD diverges, whereas SYM is stable.

Besides being symplectic, Verlet scheme is also time reversible, which means that after changing the sign of velocities, the scheme returns back in time. The scheme thus preserves the time-reversal symmetry of the original system (1a), (1b), (1c). We can verify this by substituting

which produces the same set of equations

but in the reversed order. Unfortunately, the time reversibility fails in the floating-point arithmetic (FloPA), where addition is not associative. In particular, addition and subtraction of a value to a float does not recover . For instance


Although this error is usually very small, it tends to accumulate in the presence of non-linearities and destroys the reversibility in longer simulations.

There is, however, an easy solution suggested by [13] in the context of N-body simulations of Solar System. It converts vectors , , and to the fixed-point arithmetic (FixPA) just before they are used in the time-step evaluation. Since addition is associative in FixPA, this method allows for bit-precise time reversibility and long-time conservation of energy – see Figure 3. Note that the intermediate computations of the sums in (5) and (1b) can be still performed in FloPA and only then converted to FixPA. However, due to the non-associativity of FloPA, summation of more than two elements in FloPA is order-dependent, and to ensure reversibility, the order of the summands must be chosen in a way that does not depend on the current state of the simulation (for instance, it is possible to perform summations always in the order of the particle indices).

(a) FixPA
(b) FloPA
Figure 3: Results for dam break test with parameters from Table 1 using symplectic WCSPH. At , we reset the time counter and revert all velocities. In FixPA, the simulation is completely reversible and returns to the initial state. In FloPA, this fails due to round-off errors illustrated in (8).
Figure 4: Comparison of symplectic WCSPH (SYM) in dam-break simulation with experimental data [14] (EXP) and numerical result from Violeau’s book [5] (VIO). Here, is a dimensionless time and , where is the -coordinate of the wave’s leading edge.
density 1000 kg/
spatial step 0.005 m
num. sound speed 120 m/s
gravitational constant 9.8 m/
water column width 1m
water column height 2m
box width 4m
box height 3m
kernel support radius
Lennard-Jones radius
Lennard-Jones energy
time step
Table 1: Parameters used for dam-break test.

In summary, symplectic WCSPH with the correct closed formula for density, the initial state correction, the conservative wall-fluid interaction, and the symplectic numerical Verlet scheme in the fixed-point arithmetic finally lead to globally-in-time reversible simulations. We have demonstrated the reversibility for instance on the dam break benchmark, and Appendix A contains the Gresho vortex benchmark. Although the simulation is symmetric, we can still observe the growth of entropy if we choose not see all the details of the simulation, which is the matter of the following section.

5 Emergence of the second law of thermodynamics

The second law of thermodynamics tells that the entropy of each isolated system grows until the system reaches the thermodynamic equilibrium. Where does this irreversible behavior come from when the underlying evolution equations for particle dynamics (here SPH) are reversible? In this section we illustrate that the emergence of the second law from completely reversible dynamics is caused by ignoring details of the dynamics, similarly as in [15].

The emergence of the second law is actually expected due to the result of Lanford [16] and following works [17], where a system of classical particles with a short-range potential is shown to obey the Boltzmann equation. The collision term in the Boltzmann equation then causes the growth of the Boltzmann entropy. Figure 5 shows the growth of Boltzmann entropy in a reversible SPH breaking-dam simulation.

Figure 5: Left: Histogram of the velocities of the particles at time

(blue) in comparison with the equilibrium Maxwell-Boltzmann distribution (orange) (

19) with fitted temperature . Right: Evolution of Boltzmann entropy (see Appendix B) in the reversible breaking-dam simulation (blue). Despite the reversibility of the dynamics, the entropy increases and approaches equilibrium entropy corresponding to the temperature of the measured Maxwell-Boltzmann distribution (blue line) from (21).

But how to interpret the Lanford’s mathematical result and the observed growth of Boltzmann entropy of the SPH particles from the physical point of view? The SPH particles obey reversible Hamiltonian dynamics. Another possibility to describe their motion is the Liouville equation, which expresses reversible evolution of the N-particle distribution function, . Indeed, Liouville equation follows from the particle mechanics, and, vice versa, if we set the distribution function to be the product of Dirac -distributions in the positions and momenta of all particles, we recover Hamilton canonical equations. Liouville equation also conserves the Liouville entropy


which thus remains constant. This is consistent with the reversibility of the underlying Hamiltonian particle mechanics.

Now we decide not to see all the positions and momenta of the individual particles, observing for instance only the probability distribution of momenta at given space and time. Such one-particle distribution function

indeed does not contain the knowledge of positions and momenta of all particles and thus the entropy of the system expressed in terms of must be higher than the Liouville entropy (functional of ). The former is the Boltzmann entropy and it is obtained by maximization of the Liouville entropy subject to the constraint that the one-particle distribution function is known [18]. Plugging the resulting back into the Liouville entropy, we obtain the Boltzmann entropy


see [18]. Because the Boltzmann entropy is the Liouville entropy evaluated at the point where it is maximal (keeping the knowledge of ), the Boltzmann entropy is higher than the Liouville entropy. Although the Liouville entropy remains constant in the dynamics, the Boltzmann entropy can grow in the dynamics following the Liouville equation.

The irreversibility of the dynamics of the one-particle distribution function can be illustrated on the reversible breaking-dam simulation with fixed-point arithmetic. Section 4 contains details of the simulation, but let us also show the histogram of the particle velocity distribution in the middle of the simulation, just before the velocities are inverted. Figure 5 shows the measured histogram of the velocities of the SPH particles in comparison with the equilibrium Maxwell-Boltzmann distribution function (19) obtained by fitting the effective temperature444Note that the effective temperature is not the physical temperature of the system because it is calculated from the overall motion of the macroscopic SPH particles. To see this, consider the extreme case of only one SPH particle with weight one kilogram, where the effective temperature obtained from the motion of that macroscopic particle does not coincide with the physical temperature of that particle.. Although all velocities were initially zero, their distribution approaches the equilibrium distribution. This approach towards equilibrium is reflected in the growth of Boltzmann entropy during the long-time reversible simulation.

In other words, when the particles are described by the Boltzmann equation, our knowledge is incomplete, since the positions and momenta of all the particles remain unknown, in contrast to the description by means of Hamilton canonical equations or Liouville equation. And this lack of knowledge makes the evolution of the one-particle distribution function irreversible. The irreversibility is then explicitly expressed by the collision integral in the Boltzmann equation. The second law of thermodynamics thus emerges from completely reversible dynamics when our description is incomplete (not seeing all positions and momenta of the particles).

In summary, if we see all the positions and momenta in the SPH simulation, we can not see the second law of thermodynamics. Indeed, the simulation is reversible and the Liouville entropy remains constant. However, when we only focus on the one-particle distribution function, we can see the growth of Boltzmann entropy and thus irreversible behavior. The second law of thermodynamics emerges when we our knowledge about the precise state of the system is incomplete [15].

Interestingly, since our discrete system is deterministic, reversible and can exist in only finitely many states, it is recurrent [19]. In other words, if the system evolves long enough, it comes back to the initial state exactly (provided that the particles are not allowed to escape to infinity). However, since the phase space of the system contains thousands of SPH particles, it has enormous amount of states, which makes it highly unlikely to observe the recurrence theorem in practice.

6 Conclusion

In this paper we have turned the usual weakly compressible smoothed particle hydrodynamics (WCSPH) to a symplectic form by finding the correct closed formula for the mass density. Because the closed formula for density leads to different treatment of particles that are initially near the boundary, the method of initial state correction (ICS) was introduced, due to which all particles are then treated in the same way. This leads to stable SPH simulations in the presence of free surfaces without any further stabilization.

In order to get simulations that preserve the energy, we use a conservative fluid-wall interaction. A symplectic (Verlet) scheme, which is suitable for the symplectic WCSPH, implemented in the fixed-point arithmetic then leads to globally-in-time reversible SPH simulations. This is demonstrated on the dam break benchmark, where inversion of velocities at a later stage of the simulation eventually leads the system back to the initial state. The simulations are available in a new Julia package SmoothedParticles.jl [6].

Despite the reversibility of the simulations, we observe thermodynamic behavior when we do not use all the details of the simulation. The Boltzmann entropy, which depends only on the one-particle distribution function and not on positions of individual particles, grows in the dam-break simulation and approaches the equilibrium value. In other words, we observe the emergence of the second law of thermodynamics from purely reversible dynamics, caused by reduction of our knowledge about the system.

In future, we would like to extend the WCSPH framework to non-isothermal fluids and solids while keeping the Hamiltonianity of the equations and numerical schemes.


  • [1] R. Gingold and J. Monaghan, “Smoothed particle hydrodynamics: theory and application to non-spherical stars,” Mon. Not. R. Astron. Soc., vol. 181, no. 3, pp. 375–389, 1977.
  • [2] M. Antuono, A. Colagrossi, S. Marrone, and D. Molteni, “Free-surface flows solved by means of sph schemes with numerical diffusive terms,” Computer Physics Communications, vol. 181, no. 3, pp. 532–549, 2010.
  • [3] C. Zhang, X. Y. Hu, and N. A. Adams, “A weakly compressible sph method based on a low-dissipation riemann solver,” J. Comput. Phys., vol. 335, pp. 605–620, 2017.
  • [4] Z.-F. Meng, A.-M. Zhang, P.-P. Wang, and F.-R. Ming, “A shock-capturing scheme with a novel limiter for compressible flows solved by smoothed particle hydrodynamics,” Computer Methods in Applied Mechanics and Engineering, vol. 386, p. 114082, 2021.
  • [5] D. Violeau, Fluid Mechanics and the SPH Method: Theory and Applications. Oxford University Press, Oxford, UK, 2012.
  • [6] O. Kincl and M. Pavelka, “SmoothedParticles.jl,” 2021,
  • [7] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM Review, vol. 59, no. 1, pp. 65–98, 2017.
  • [8] J. Monaghan, “Smoothed particle hydrodynamics,” Reports on Progress in Physics, vol. 68, p. 1703, 2005.
  • [9] M. Kroeger and M. Huetter, “Automated symbolic calculations in nonequilibrium thermodynamics,” Comput. Phys. Commun., vol. 181, p. 2149–2157, 2010.
  • [10] J. J. Monaghan, “Simulating free surface flows with sph,” Journal of computational physics, vol. 110, no. 2, pp. 399–406, 1994.
  • [11] B. Leimkuhler, S. Reich, and C. U. Press, Simulating Hamiltonian Dynamics. Cambridge Monographs on Applie, Cambridge University Press, 2004.
  • [12] J. Candy and W. Rozmus, “A symplectic integration algorithm for separable hamiltonian functions,” Journal of Computational Physics, vol. 92, no. 1, pp. 230 – 256, 1991.
  • [13] H. Rein and D. Tamayo, “JANUS: a bit-wise reversible integrator for N-body dynamics,” Monthly Notices of the Royal Astronomical Society, vol. 473, pp. 3351–3357, 2017.
  • [14] S. Koshizuka and Y. Oka, “Moving-particle semi-implicit method for fragmentation of incompressible fluid,” Nuclear science and engineering, vol. 123, no. 3, pp. 421–434, 1996.
  • [15] M. Pavelka, V. Klika, and M. Grmela, “Generalization of the dynamical lack-of-fit reduction,” Journal of Statistical Physics, vol. 181, no. 1, pp. 19–52, 2020.
  • [16] O. E. Lanford, “Time evolution of large classical systems,” in Dynamical systems, theory and applications (J. Moser, ed.), vol. 38 of Lecture Notes in Physics, pp. 1–111, Springer–Verlag, Berlin, 1975.
  • [17] M. Pulvirenti, C. Saffirio, and S. Simonella, “On the validity of the boltzmann equation for short range potentials,” Reviews in Mathematical Physics, vol. 26, no. 02, p. 1450001, 2014.
  • [18] M. Pavelka, V. Klika, and M. Grmela, Multiscale Thermo-Dynamics. Berlin: de Gruyter, 2018.
  • [19] D. Wallace, “Recurrence theorems: a unified account,” Journal of Mathematical Physics, vol. 56, no. 2, p. 022105, 2015.
  • [20] H. Grimm-Strele, F. Kupka, and H. Muthsam, “Curvilinear grids for weno methods in astrophysical simulations,” Computer Physics Communications, vol. 185, no. 3, pp. 764–776, 2014.

Appendix A Gresho vortex benchmark

The Gresho vortex benchmark [20] prescribes a tangential velocity component of Eulerian fluid in polar coordinates as

The vortex is confined in a box with no-slip wall (though different variants can be found in literature). Theoretically, in the absence of viscosity, the vortex should be stationary, but obtaining this in a numerical scheme is challenging. In our numerical experiment, we simulate the flow in time-interval and measure the error using a discretized Sobolev norm


The meaning of the factor is that the zero velocity field corresponds to error . In order to prevent formation of void space in the center of vortex, we had to include an anti-clump force term in the form

Figure 6: Simulation result with Vogel spiral + ISC arrangement at (left) and (right).
no filter passive filter active filter
square 28.16% 12.20% 53.07%
hexagonal 27.84% 11.77% 55.79%
Vogel 29.85% 13.31% 78.06%
Vogel + ISC 28.19% 11.26% 55.25%
Table 2: Error of Gresho vortex benchmark for different grid types (square, hexagonal, Vogel spiral, Vogel spiral + ISC) and with different types of noise filtering. The error is computed according to formula (11). It appears from data that passive (a posteriori) filter is better than active filtering (which dissipates energy). Compared to WENO method [20], numerical dissipation in symplectic WCSPH is quite strong.
density 1
spatial step 1e-2
num. sound speed 20
kernel support radius
box size 1
anti-clump pressure 10
time step
filter frequency 30
Table 3: Paramers used for Gresho vortex simulation. This benchmark is dimensionless and with no gravity. Walls were implemented by two layers of dummy particles.

where is a smoothing kernel with support radius . Additionally, we consider noise reduction by the Shepard filter


This filter can be applied either in post-processing (passive filter), or by setting every time steps (active filter). Results are summarized in Table 2 and Figure 6. Simulation parameters are listed in Table 3.

Appendix B The Maxwell-Boltzmann entropy

In this section we derive the formula for Boltzmann entropy in terms of the Maxwell-Boltzmann distribution and, subsequently, to find its equilibrium value.

b.1 The reduced Boltzmann entropy in terms of the Maxwell-Boltzmann distribution

Let us start with the Boltzmann entropy in two dimensions expressed in terms of the one-particle distribution function ,


where the position is constrained to a box with volume and where

is the Planck constant. This formula can be obtained for instance by the principle of maximum entropy from the Liouville entropy, from where it also follows that the one-particle distribution function is normalized to the number of particles, see


Assuming that the distribution function depends only on the magnitude of the momentum (isotropic dependence), the Boltzmann entropy can be rewritten as


where is the norm of the momentum. This dependence on the one-particle distribution function has to be converted to a dependence on the Maxwell-Boltzmann distribution function

, which tells the probability that the norm of velocity of a particle

is in the interval . This is done using the normalization,


which leads to the expression of the two-dimensional Maxwell-Boltzmann distribution function in terms of the one-particle distribution,


which is normalized to unity. Note that we assume that the one-particle distribution is homogeneous in space (the total volume being ) and isotropic in momentum. These assumptions are valid in the thermodynamic equilibrium, but they approximately hold also before the equilibrium is reached, and for the purpose of showing that the Boltzmann entropy grow in our simulations, the approximation is satisfactory. Finally, the Boltzmann entropy in terms of the Maxwell-Boltzmann distribution function becomes


Since the number of particles is conserved in our simulation, we only evaluate the part of Boltzmann entropy that varies when changes,


called reduced Boltzmann entropy. Note that the integral converges near the origin because the Maxwell-Boltzmann distribution typically has linear growth there.

In our simulations, this reduced Boltzmann entropy is numerically integrated using the approximate Maxwell-Boltzmann distribution function obtained from a histogram of the particle velocities at given time instant.

b.2 Equilibrium Boltzmann entropy

What is the final value of the reduced Boltzmann entropy when the system of particles reaches the thermodynamic equilibrium? This question can be answered in two steps. First, the equilibrium distribution function is calculated by the principle of maximum entropy (MaxEnt). Second, the equilibrium distribution function is plugged into the formula for the reduced Boltzmann entropy. Note, however, that the MaxEnt step depends on the energy (Hamiltonian) of the system and it depends on the complexity of the Hamiltonian whether the calculation can proceed purely analytically (without numerical solutions). Therefore, we use the simplest Hamiltonian approximating the true Hamiltonian of our SPH particles, consisting only of the kinetic energy of the particles. Maximization of the Boltzmann entropy subject to the constraints given by the total energy and total number of particles leads to the equilibrium one-particle distribution function


where and are the Lagrange multipliers corresponding to the two constraints (number of particles and total energy). From the normalization it follows that , while the other Lagrange multiplier can be interpreted as the inverse temperature, .

The equilibrium Maxwell-Boltzmann distribution function then becomes


As time proceeds in our simulations, the histogram of particle velocities approaches the equilibrium Maxwell-Boltzmann distribution.

When this equilibrium distribution function is plugged back into the Boltzmann entropy, the reduced Boltzmann entropy becomes


This equilibrium entropy, which depends on , can be calculated once the temperature is obtained by fitting the histogram of particle velocities to the equilibrium Maxwell-Boltzmann distribution function (19). This makes sense, however, only at later stages in the simulations, when the histogram approaches the equilibrium distribution.

Instead of using the equilibrium temperature, which makes sense only in later stages of the simulations, we can express the equilibrium entropy in terms of energy, which can be measured anytime. The total kinetic energy of the two-dimensional system of particles is equal to , which makes it possible to write the equilibrium entropy in terms of the energy,


This value of equilibrium entropy is close to the Boltzmann entropy obtained directly from the approximated Maxwell-Boltzmann distribution, while the entropy based on the equilibrium temperature is only approached at later stages of the simulations.


OK was supported by project No. START/SCI/053 of Charles University Research program. MP was supported by project No. UNCE/SCI/023 of Charles University Research program. OK and MP were also supported by the Czech Science Foundation (project no. 20-22092S)