Approaches to conservative Smoothed Particle Hydrodynamics with entropy

Hamiltonian continuum mechanics on the Lagrangian (reference) manifold can be derived by means of the principle of stationary action. The Lagrangian is a function of the Lagrange$\rightarrow$Euler mapping $\mathbf{x}(t,\mathbf{X})$ and its time-derivative, where $t$ is the time, $\mathbf{X}$ the Lagrangian coordinates of a point, and $\mathbf{x}$ its Eulerian coordinate (current configuration) [17]. The principle of stationary action then gives the Lagrange-Euler equation, which can be rewritten (by the Legendre transformation) to the Hamiltonian setting with the fields of position and momentum density, $x^{i}(\mathbf{X})$ and $M_{i}(\mathbf{X})$, as the state variables. These state variables are equipped with the canonical Poisson bracket

$$\{F,G\}^{(\mathrm{Lagrangian})}=\int d\mathbf{X}\left(F_{x^{i}}G_{M_{i}}-G_{x^{i}}F_{M_{i}}\right),$$

(1)

and Hamiltonian evolution of any functional of the state variables, $F(\mathbf{x},\mathbf{M})$, is given by the Poisson bracket of the functional and the energy $E(\mathbf{x},\mathbf{M})$,

$$\dot{F}=\{F,E\}.$$

(2)

In particular, if $F$ is chosen as the state variables, we obtain the Hamiltonian evolution equations of the positions and momenta,

	$\displaystyle\dot{\mathbf{x}}$	$\displaystyle=E_{\mathbf{M}}$		(3a)
	$\displaystyle\dot{\mathbf{M}}$	$\displaystyle=-E_{\mathbf{x}},$		(3b)
where the subscripts denote derivatives (more precisely, the functional derivatives [43]).

Equations (3) represent the canonical Hamiltonian equations for the Lagrangian continuum mechanics, which can be used as dynamics of elastic bodies [31, 45].

However, to close the equations, we have to prescribe the energy $E(\mathbf{x},\mathbf{M})$. If we know the internal energy density per the Eulerian volume $\epsilon(\rho)$, the Lagrangian energy is the sum of the kinetic energy and the internal energy transformed to the Lagrangian manifold,

$$E^{Lagrangian}=\int d\mathbf{X}\left(\frac{\mathbf{M}^{2}}{2\rho_{0}}+\det\frac{\partial\mathbf{x}}{\partial\mathbf{X}}\epsilon(\rho(\mathbf{x}(\mathbf{X})))\right).$$

(4)

Here, $\rho$ is the Eulerian mass density, which is related to the Lagrangian density $\rho_{0}(\mathbf{X})$ (mass density per the Lagrangian volume) through

$$\rho(\mathbf{x})=\rho_{0}(\mathbf{X}(\mathbf{x}))\det\frac{\partial\mathbf{X}}{\partial\mathbf{x}}.$$

(5)

Without going to more details on the Hamiltonian continuum mechanics on the Lagrangian manifold, let us go to the Hamiltonian formulation of the standard SPH method.

The standard SPH [15, 55] expresses evolution of positions and velocities of the smoothed particles, and the reversible part of the evolution (disregarding viscosity) is also a Hamiltonian system. To write down the classical SPH formulation, we need two approximations:

1. 1.

   smoothed state variables (positions and momenta of the particles)

2. 2.

   differential operators (appearing in the evolution equations like $\mathrm{div}$ or $\mathrm{grad}$)

Different operator discretization may yield different results as clearly observed in the discussion of density evolution [55, Chap 5.3.1]. Here, we propose quite a different approach, employing the least unnecessary knowledge for deriving discrete operators, thus reducing ambiguity in the formulation of SPH.

The mapping

	$\displaystyle\mathbf{x}_{\alpha}$	$\displaystyle=\int d\mathbf{X}\bar{\chi}_{\alpha}(\mathbf{X})\mathbf{x}(\mathbf{X})$		(7a)
	$\displaystyle\mathbf{M}_{\alpha}$	$\displaystyle=\int d\mathbf{X}V_{0\alpha}\bar{\chi}_{\alpha}(\mathbf{X})\mathbf{M}(\mathbf{X})$		(7b)
then defines the position and momentum of the SPH particle $\alpha$, which play the role of the SPH state variables. Later, in Section 3, we will also add entropy to the state variables, but first we show how to derive the standard SPH by reduction from the Lagrangian continuum mechanics.

In order to reduce the continuum Poisson bracket (1) to a Poisson bracket for the SPH state variables ($\mathbf{x}_{\alpha}$ and $\mathbf{M}_{\alpha}$), we plug functionals dependent on the SPH variables into the Poisson bracket. Because the Poisson bracket contains derivatives of the functionals, we need to calculate derivatives of the SPH state variables with respect to the continuum state variables,

$$\frac{\delta x^{i}_{\alpha}}{\delta x^{j}(\mathbf{X})}=\delta^{i}_{j}\bar{\chi}_{\alpha}(\mathbf{X})\quad\text{and}\quad\frac{\delta M_{\alpha i}}{\delta M_{j}(\mathbf{X})}=\delta^{j}_{i}V_{0\alpha}\bar{\chi}_{\alpha}(\mathbf{X})$$

(8)

while the remaining derivatives are zero. Derivatives of an arbitrary functional of the SPH state variables, $F(\mathbf{x}_{\alpha},\mathbf{M}_{\alpha})$, are then

$$\frac{\delta F}{\delta x^{i}(\mathbf{X})}=\sum_{\alpha}\frac{\partial F}{\partial x^{i}_{\alpha}}\bar{\chi}_{\alpha}(\mathbf{X})\quad\text{and}\quad\frac{\delta F}{\delta M_{i}(\mathbf{X})}=\sum_{\beta}\frac{\partial F}{\partial M_{\beta i}}V_{0\beta}\bar{\chi}_{\beta}(\mathbf{X}).$$

(9)

Plugging two arbitrary functionals dependent on the discrete state variables to Poisson bracket (1) leads to

	$\displaystyle\{F,G\}^{(\mathrm{SPH})}$	$\displaystyle=\sum_{\alpha}\sum_{\beta}\left(F_{x^{i}_{\alpha}}G_{M_{\beta i}}-G_{x^{i}_{\alpha}}F_{M_{\beta i}}\right)\int d\mathbf{X}V_{0\beta}\bar{\chi}_{\alpha}(\mathbf{X})\bar{\chi}_{\beta}(\mathbf{X})$
		$\displaystyle=\sum_{\alpha}\left(F_{x^{i}_{\alpha}}G_{M_{\alpha i}}-G_{x^{i}_{\alpha}}F_{M_{\alpha i}}\right),$		(10)

which is the SPH Poisson bracket expressing kinematics of the SPH state variables. This Poisson bracket can be seen as the Hamiltonian counterpart of the Lagrangian form of SPH [55].

The reversible Hamiltonian evolution equations implied by this bracket are

	$\displaystyle\dot{\mathbf{x}}_{\alpha}$	$\displaystyle=E_{\mathbf{M}_{\alpha}}$		(11a)
	$\displaystyle\dot{\mathbf{M}}_{\alpha}$	$\displaystyle=-E_{\mathbf{x}_{\alpha}},$		(11b)
which are the Hamilton canonical equations for the SPH state variables.

However, to close the equations, we need to express the energy in terms of the SPH state variables.

To find the dependence of the energy on the SPH state variables, we need to approximate the exact energy of the continuous system (4), which depends also on the full Lagrange$\rightarrow$Euler mapping $\mathbf{x}(\mathbf{X})$. Knowledge of discrete particle positions is not enough for the precise reconstruction of the exact continuum energy, and hence we turn to the standard SPH smoothing using the weighing functions $W$:

$$\overline{A}_{\alpha}=\sum_{\beta}A_{\beta}W(|\mathbf{x}_{\alpha}-\mathbf{x}_{\beta}|),$$

(12)

which is an approximate form of mollification. It can be shown [55] that if the weighing function satisfies certain normalization property, it provides an approximation precise up to the second-order. The normalization requirement can be translated into a condition that $W$ is a homogeneous function of the distance of a degree $-n$, that is $W(\alpha|\mathbf{x}|)=\alpha^{-n}W(|\mathbf{x}|)$, where $n$ denotes the spatial dimension. From this homogeneous property, one can also see that any scaling of the local volume (as in bulk expansion) is exactly captured by the weighing function $W$. Finally, as the dimension of $W$ is $1/m^{n}$, we may define any density variable discretization in Eulerian frame $\overline{A}_{\alpha}$ from its Lagrangian counterpart $A_{\alpha}$ via the above expression (12). We shall use such smoothing to define the SPH mass density $\rho_{\alpha}$. In particular, $\rho_{\alpha}$ is the Eulerian mass density of particle $\alpha$, $m_{\alpha}=\int d\mathbf{X}\rho_{0}(\mathbf{X})\bar{\chi}_{\alpha}(\mathbf{X})$ is the mass of the Lagrangian particle and hence

$$\rho_{\alpha}=\sum_{\beta}m_{\beta}W_{\alpha\beta},$$

(13)

where $W_{\alpha\beta}\stackrel{{\scriptstyle def}}{{=}}W(|\mathbf{x}_{\alpha}-\mathbf{x}_{\beta}|)$.

Now we can close the Hamilton canonical equations (11) by supplying an energy functional that depends only on the SPH state variables. This is actually the step that depends on material properties, since the Poisson bracket (2.2) is already fixed by the SPH state variables. The continuum energy (4) can be approximated by

	$\displaystyle E^{SPH}$	$\displaystyle=\sum_{\alpha}V_{0\alpha}\frac{\mathbf{M}^{2}_{\alpha}}{2\rho_{0\alpha}}+\sum_{\alpha}V_{0\alpha}\frac{\rho_{0\alpha}}{\rho_{\alpha}}\epsilon(\rho_{\alpha})$
		$\displaystyle=\sum_{\alpha}\frac{\mathbf{M}^{2}_{\alpha}}{2m_{\alpha}}+\sum_{\alpha}V_{\alpha}\epsilon(\rho_{\alpha}),$		(14)

where $\rho_{0\alpha}=m_{\alpha}/V_{0\alpha}$ is its density in the Lagrangian space. This approximate energy functional can also be obtained by means of the principle of maximum entropy (MaxEnt), which is shown in Appendix A. We keep the Eulerian volume $V_{\alpha}$ unspecified (unrelated to Eulerian density, for example), as it is a key step subjected to further in Section 3, where five possible definitions are shown.

The purpose of this Section is to collect the above results on the SPH Poisson bracket and SPH energy and to show that they indeed lead to the standard formulation of SPH.

When energy (2.3) is plugged into the Hamilton canonical equations (11), we need to take its derivatives with respect to the SPH state variables. In particular, derivative of $\rho_{\alpha}$ (defined in Equation (13)) with respect to the positions reads

	$\displaystyle\frac{\partial\rho_{\beta}}{\partial x^{i}_{\alpha}}$	$\displaystyle=\sum_{\gamma}m_{\gamma}W^{\prime}(|\mathbf{x}_{\beta}-\mathbf{x}_{\gamma}|)\frac{\partial\sqrt{(x^{j}_{\beta}-x^{j}_{\gamma})(x^{j}_{\beta}-x^{j}_{\gamma})}}{\partial x^{i}_{\alpha}}$
		$\displaystyle=\sum_{\gamma}m_{\gamma}W^{\prime}(|\mathbf{x}_{\beta}-\mathbf{x}_{\gamma}|)\frac{x^{j}_{\beta}-x^{j}_{\gamma}}{|\mathbf{x}_{\beta}-\mathbf{x}_{\gamma}|}\delta^{j}_{i}\left(\delta_{\beta\alpha}-\delta_{\gamma\alpha}\right)$
		$\displaystyle=\sum_{\gamma}m_{\gamma}W^{\prime}_{\beta\gamma}e_{i\beta\gamma}\left(\delta_{\beta\alpha}-\delta_{\gamma\alpha}\right),$		(15)

where $W^{\prime}_{\beta\gamma}=W^{\prime}(|\mathbf{x}_{\beta}-\mathbf{x}_{\gamma}|)$ is the derivative of the SPH kernel $W$ evaluated on the distance between particles $\beta$ and $\gamma$ and where $\mathbf{e}_{\beta\gamma}=\frac{\mathbf{x}_{\beta}-\mathbf{x}_{\gamma}}{|\mathbf{x}_{\beta}-\mathbf{x}_{\gamma}|}$ is the unit vector¹¹1Actually, it is the covector expressing gradient of the distance between the particles while the vector is dual to the covector. But in Euclidean space we do not need to distinguish them [14]. pointing from $\mathbf{x}_{\gamma}$ to $\mathbf{x}_{\beta}$.

Now, we are finally in position to write down the resulting ordinary differential equations expressing the Hamiltonian evolution of the SPH state variables,

	$\displaystyle\dot{x}^{i}_{\alpha}$	$\displaystyle=\frac{M^{i}_{\alpha}}{m_{\alpha}}$		(16a)
	$\displaystyle\dot{M}_{\alpha i}$	$\displaystyle=-\sum_{\beta}\left(-\frac{m_{\beta}}{\rho^{2}_{\beta}}\epsilon+\frac{m_{\beta}}{\rho_{\beta}}\frac{\partial\epsilon}{\partial\rho_{\beta}}\right)\sum_{\gamma}m_{\gamma}W^{\prime}_{\beta\gamma}e_{i\beta\gamma}(\delta_{\beta\alpha}-\delta_{\gamma\alpha})$
		$\displaystyle=-\sum_{\beta}\left(\frac{m_{\alpha}m_{\beta}}{\rho^{2}_{\alpha}}p_{\alpha}+\frac{m_{\alpha}m_{\beta}}{\rho^{2}_{\beta}}p_{\beta}\right)W^{\prime}_{\alpha\beta}e_{i\alpha\beta},$		(16b)

where the pressure of particle $\alpha$ is defined by the usual barotropic relation

$$p_{\alpha}=-\epsilon+\rho_{\alpha}\frac{\partial\epsilon}{\partial\rho_{\alpha}},$$

(17)

see [7, 43]. When the dependence of the internal energy on the mass density is prescribed, which can be integrated from an equation of state, the system of equations is completely specified. See [55] for some choices of internal energy typical in SPH.

Equations (16) represent a symplectic SPH discretization of fluid mechanics where the energy depends on the mass density [55]. This system of equations forms a Hamiltonian system because it is constructed from a Poisson bracket and energy. Combined with a symplectic time integrator, such as the Verlet scheme, the numerical solution has energy error bounded uniformly with respect to the simulation time $t$, and the numerical schemes can even be made time reversible by using fixed-point arithmetic [26].

However, it is customary in the SPH literature [15] to use the discretized continuity equation instead of the closed expression (13). Then, Equations (16) can be generalized, using the definition of $\rho_{a}$, to another Hamiltonian system,

	$\displaystyle\dot{x}^{i}_{\alpha}$	$\displaystyle=\frac{M^{i}_{\alpha}}{m_{\alpha}}$		(18a)
	$\displaystyle\dot{\rho}_{\alpha}$	$\displaystyle=\sum_{\beta}m_{\beta}W^{\prime}_{\alpha\beta}\mathbf{e}_{\alpha\beta}\left(\frac{\mathbf{M}_{\alpha}}{m_{\alpha}}-\frac{\mathbf{M}_{\beta}}{m_{\beta}}\right)$		(18b)
	$\displaystyle\dot{M}_{\alpha i}$	$\displaystyle=-\sum_{\beta}\left(\frac{m_{\alpha}m_{\beta}}{\rho^{2}_{\alpha}}p_{\alpha}+\frac{m_{\alpha}m_{\beta}}{\rho^{2}_{\beta}}p_{\beta}\right)W^{\prime}_{\alpha\beta}e_{i\alpha\beta},$		(18c)
where the pressure is again given by Equation (17). Equations (18) are the standard SPH equations [55].

This system of equations is Hamiltonian because it is generated by a Poisson bracket

	$\displaystyle\{F,G\}^{(\mathrm{SPH},\rho_{\alpha})}$	$\displaystyle=\{F,G\}^{(\mathrm{SPH})}$
		$\displaystyle+\sum_{\alpha}\sum_{\beta}m_{\beta}W^{\prime}_{\alpha\beta}e_{i\alpha\beta}\left(F_{\rho_{\alpha}}(G_{M_{\alpha i}}-G_{M_{\beta i}})-G_{\rho_{\alpha}}(F_{M_{\alpha i}}-F_{M_{\beta i}})\right),$		(19)

which can be obtained by projection of the continuum Poisson bracket (1) to the state variables $\mathbf{x}_{\alpha}$, $\mathbf{M}_{\alpha}$, and $\rho_{\alpha}$, see Appendix B for details (here we dropped $s_{\alpha}$ from state variables considered in Appendix B).

Poisson bracket (2.4) is different from the bracket (2.2) and leads to Hamiltonian but non-symplectic mechanics [36]. From a numerical perspective, it is advantageous to keep the symplecticity (non-degenerate Poisson bracket) because, for instance, the energy error is then capped by a constant decreasing with the time step when a symplectic integrator is used [20, 33]. Therefore, we make the dynamics symplectic by directly evaluating the density from its definition (13) at each time step while keeping the evolution equations for $\mathbf{x}_{\alpha}$ and $\mathbf{M}_{\alpha}$. This results, for example, in a globally reversible scheme for SPH [26]. On the other hand, this direct evaluation of the density within the boundaries may cause other numerical artifacts [2] if not properly treated [26].

In summary, the standard SPH equations (18) can be derived from Lagrangian continuum mechanics by reducing the continuum Poisson bracket (1) to the SPH Poisson bracket (2.4). The reduction relies on the projection from the continuous fields of positions and momenta to the SPH positions and momenta (7), where the Lagrangian manifold is partitioned into discrete particles. Then, the SPH energy (2.3) is found as an approximation of the continuum energy (4), and this approximation can be seen as an application of the principle of the MaxEnt principle (Sect. A). When a concrete dependence of the internal energy on the mass density is chosen (a material is chosen), the SPH equations (18) can be solved as long as one chooses an appropriate interpretation of the particle volume.

However, what if the energy also depends on the entropy of the material, which means that the fluid is not barotropic? Then, the formulation of SPH has to be enriched to include also the entropy of the SPH particles, which is the purpose of the following Section.

The fundamental idea behind SPH discretization lies in approximative formulas for the Dirac delta function and convolution. In particular, if we take a classical integrable and normalized function, that is, $\int_{\mathbb{R}}W=1$, and define a rescaled function as $W_{h}(x)=\frac{1}{h}W(x/h)$, then $\int_{\mathbb{R}}W_{h}=1$ and $\lim_{h\to 0}W_{h}=\delta(x)$. As $f\star\delta=f$, we get $f\approx f\star W_{h}$. In particular, we have

$$f(x_{\alpha})\approx\int_{\mathbb{R}^{n}}f(y)W_{h}(x-y)\approx\sum_{\beta}V_{\beta}f(x_{\beta})W_{h}(x_{\alpha}-x_{\beta}).$$

With particular choices of the approximated function $f$, we obtain different approximation formulas for volume and different corresponding errors.

For the Wendland kernel, one can show that the relative error in 1D is proportional to the gradient of a logarithm at that point [55, Eq. (5.115)]

$$\frac{f(x)-f(x_{\alpha})}{f(x)}=C(h,r)\frac{d\ln f}{dx},$$

(22)

where $C(h,r)$ is a constant independent of the function $f$ but dependent of the kernel width $h$ and particle size $r$.

Using the above approximations of convolution with the Dirac delta function, we may propose several distinct definitions of particle volumes. Namely, the choice $f=\rho$ gives

$$\rho(x_{\alpha})=\rho_{\alpha}\approx\sum_{\beta}V_{\beta}\rho(x_{\beta})W_{\alpha\beta}=\sum_{\beta}m_{\beta}W_{\alpha\beta},$$

and hence entails the above mass density update expression Eq. (13) or, generally, Eq. (12). From this expression, one may define a particle volume (based on mass density) as

$$(V^{m}_{\alpha})^{-1}=\frac{\rho_{\alpha}}{m_{\alpha}}\approx\sum_{\beta}\frac{m_{\beta}}{m_{\alpha}}W_{\alpha\beta}.$$

(23)

However, the introduction of another volume density variable, as is the case with entropy density, gives another possibility

$$(V^{s}_{\alpha})^{-1}=\frac{s_{\alpha}}{S_{\alpha}}\approx\sum_{\beta}\frac{S_{\beta}}{S_{\alpha}}W_{\alpha\beta},$$

(24)

where we chose $f=s$ and where $S_{\sigma}$ is the entropy of particle $\alpha$.

At the same time, the choice $f:x\mapsto 1/V_{\alpha}$, $x\in V_{\alpha}$ gives directly an estimate of the particle volume as

$$(V^{d}_{\alpha})^{-1}\approx\sum_{\beta}W_{\alpha\beta}.$$

(25)

Finally, the choice $f=\chi_{\Omega}$ being the characteristic function of a given domain $\Omega$ yields an implicit relation for particle volumes

$$1\approx\sum_{\beta}\tilde{V}^{I}_{\beta}W_{\alpha\beta},$$

where $\tilde{V}^{I}_{\alpha}$ is the volume of the part of the particle that is within the domain $\Omega$. This expression can be rearranged into

$$\tilde{V}_{\alpha}^{I}=\sum_{\beta}W_{\alpha\beta}^{-1},$$

(26)

where $W_{\alpha\beta}^{-1}$ denotes the inverse to the matrix $W_{\alpha\beta}$.

Let us make few observations regarding the errors and convergence. All four definitions of volumes are based on Eulerian quantities whose evolution we aim to calculate. Hence, the error upper bound is only a rough guidance to the suitable volume definition selection as it is dependent on the systems evolution. Nevertheless, we expect from Eq. (22) that the mass-based volume should give a good particle volume representations in systems, where the spatial gradient of density is small and entropy-based volume should similarly yield a well-represented system with shallow entropy gradients. The direct definition of particle volume, $V^{d}_{\alpha}$, has errors proportional to the spatial variation of particle volumes and hence is dependent on the heterogeneity of discretization itself. Finally, the implicit definition of volume, $V^{I}_{\alpha}$, can be expected to have the smallest discretization error in the bulk, however, is expected to introduce errors in the boundary particles where there is a jump in the discretized function – the characteristic function of the domain $\Omega$. In addition, the determination of implicit volume $V^{I}_{\alpha}$ comes at a greater computational cost as it involves inverting the large matrix $W_{\alpha\beta}$.

The precision of the implicit particle volumes $V^{I}_{\alpha}$ for the bulk integration can be perceived from the opposite perspective: we are looking for particle volumes such that the integral is exactly represented by the sum. Namely, the particle volumes $V^{I}_{\alpha}$ can be equivalently defined as weights of a Gaussian quadrature rule

$$\int f\mathrm{d}\mathbf{x}\approx Q[f]=\sum_{\alpha}V^{I}_{\alpha}f(\mathbf{x}_{\alpha}),$$

which is exact for all functions in the form

$$f(\mathbf{x})=\sum_{\beta}f_{\beta}W(\mathbf{x}-\mathbf{x}_{\beta})$$

(27)

(where $f_{\beta}\in\mathbb{R}$ are arbitrary coefficients). Indeed, from the definition (26):

$$Q[f]=\sum_{\alpha}V^{I}_{\alpha}f(\mathbf{x}_{\alpha})=\sum_{\alpha}\sum_{\beta}V^{I}_{\alpha}f_{\beta}W_{\alpha\beta}=\sum_{\beta}f_{\beta}=\sum_{\beta}f_{\beta}\int W(\mathbf{x}-\mathbf{x}_{\beta})\mathrm{d}\mathbf{x}=\int f(\mathbf{x})\mathrm{d}\mathbf{x}.$$

To obtain the reverse implication, substitute $f(\mathbf{x})=W(\mathbf{x}-\mathbf{x}_{\alpha})$ into $\int f\mathrm{d}\mathbf{x}=Q[f]$. Note that again the issue with boundary elements appear, here in the requirement of the function form (27), which effectively requires $\Omega$ to be exactly discretized into SPH particles.

Let us, moreover, comment on the convergence of all the four above definitions of volume. The discretization of the convolution integral and Dirac delta function are improving, and hence the error being reduced, if the discretization step is decreased. At the same time, it is known that there is a tradeoff between the particle size and the SPH kernel width, reaching an optimum and preventing the truncation error to be arbitrarily small, [48].

Yet another possibility is to define both mass and entropy densities as the SPH averages,

	$\displaystyle\rho_{\alpha}=$	$\displaystyle\sum_{\beta}m_{\beta}W_{\alpha\beta}$		(28a)
	$\displaystyle s_{\alpha}=$	$\displaystyle\sum_{\beta}S_{\beta}W_{\alpha\beta}.$		(28b)

This choice, called the mixed-volume approach, however, does not specify how a single particle volume is defined and instead needs two particle volumes at once. Despite this loss of elegance, it shows good results in the presence of discontinuities, see Section 3.5.

We shall now turn to the discussion of the suitability of the particle volume definitions from a different perspective – the validity of discretized conservation laws.

Using the mass-based volume $V^{m}_{\alpha}$, Eq. (23), has the advantage of conserving the total mass, $M=\sum_{\alpha}V^{m}_{\alpha}\rho_{\alpha}=\sum_{\alpha}m_{\alpha}$, as particle masses $m_{\alpha}$ are constant. Similarly, the total entropy, $S=\sum_{\alpha}V^{m}_{\alpha}s_{\alpha}=\sum_{\alpha}S_{\alpha}$, is conserved by the reversible part of the evolution (the Hamiltonian part, which does not affect $S_{\alpha}$).

The expression for total energy, Eq. (3), however, differs according to the choice of the particle volume definition. In the mass-based volume case, it reads

$$E^{SPH}=\sum_{\alpha}\frac{\mathbf{M}^{2}_{\alpha}}{2m_{\alpha}}+\sum_{\alpha}\frac{m_{\alpha}}{\rho_{\alpha}}\epsilon(\rho_{\alpha},s_{\alpha}).$$

If we consider the entropy-based particle volumes, the energy discretization takes on a slightly different form:

$$E^{SPH}=\sum_{\alpha}\frac{\mathbf{M}^{2}_{\alpha}}{2m_{\alpha}}+\sum_{\alpha}\frac{S_{\alpha}}{s_{\alpha}}\epsilon(\rho_{\alpha},s_{\alpha}).$$

In the rest of this paper, we shall explore these four choices of particle volume definitions both in the reversible and irreversible processes, and discuss the conservation of energy.

As a purely mathematical comparison, let us investigate the accuracy of numerical integration

$$\int f_{1}(\mathbf{x})\,\mathrm{d}\mathbf{x}\approx\sum_{a}V_{a}f(\mathbf{x}_{a})$$

(29)

for the different definitions of particle volume $V_{a}$. To this end, we consider functions

$$\begin{split}f_{1}&=e^{-x-y}\\ f_{2}&=1+\cos 8x+\sin 8x\\ f_{3}&=1-xy+4x^{2}y^{2}\end{split}$$

(30)

in a square $[0,1]\times[0,1]$. The exact integral can be easily evaluated. We approximate this value by $N$ particles arranged in a hexagonal pattern. Let us define the spatial resolution of the problem as

$$\mathrm{d}r=\frac{1}{\sqrt{N}}.$$

(31)

The particle masses $m_{\alpha}$ and entropies $S_{\alpha}$ were all set equal to 1 and Wendland’s quintic kernel with support radius $h=3\mathrm{d}r$ was used. The result is shown in Figure 2. We found that the implicit volume is more precise in most cases (at the expense of a longer computation time). Therefore, we suggest that the implicit-volume quadrature should be used in a situation where an integral has to be computed from data in a post-process.