# How is the divergence and curl approximated in Smoothed Particle Hydrodynamics?

I was reading through Monaghan, SPH, 1992 paper and he states that any quantity $A$ at $r_b$ can be approximated by $A(r) = \sum_b m_b \frac{A_b}{\rho_b}W(r-r_b,h)$ and so the gradient of that quantity can be computed as $\nabla A(r) = \sum_b m_b \frac{A_b}{\rho_b} \nabla W(r-r_b,h)$.

So if we want to approximate the density we have $\rho (r) = \sum_b m_b W(r-r_b,h)$, so far so good.

Then there is the equation to approximate the divergence of the velocity as $\nabla \cdot v = \sum_b m_b v_b \cdot \nabla W(r-r_b,h)$ which I don't understand exactly where it comes from. Is the velocity $v$ playing the role of the approximated quantity, $A_b = v$ ?

Now, then the

second golden rule of SPH is to rewrite formulae with density placed inside operators

which gives $\nabla \cdot v = [\nabla \cdot (\rho v) - v \cdot \nabla \rho] / \rho$. How is this obtained? why is the density appearing in the formula and how is related to the above formula for the divergence of the velocity?

I guess something similar happens when deriving the estimator for the curl of the velocity.

Can someone explain a little more in details how are the formulas derived?

• Are you maybe missing a gradient of $W$ in your expansion of the divergence of $v$? My searching seems to indicate that it should be there, and then the derivation should be straightforward. – Bill Barth Nov 29 '14 at 16:50
• @BillBarth I was indeed missing a nabla, could you explain me the derivation? – BRabbit27 Nov 29 '14 at 16:55

When you approximate a quantity, for example $A(r) = \sum_b m_b \frac{A_b}{\rho_b}W(r-r_b,h)$, then what you are saying is that you are assigning a scalar value $A_b$ to every nodal point $b$. At any point $r$, this formula then simply provides a way to interpolate between the points $r_b$ in the neighborhood. Because the masses and nodal values $A_b$ are just fixed scalars, the gradient of your approximation is simply $\nabla A(r) = \sum_b m_b \frac{A_b}{\rho_b}\nabla W(r-r_b,h)$.

You can do the same construction with vector fields. If you assign a vector $v_b$ to every nodal point $b$, then your approximation becomes $\mathbf v(r) = \sum_b m_b \frac{\mathbf v_b}{\rho_b}W(r-r_b,h)$ and its divergence -- because the $\mathbf v_b$ are just fixed vectors, independent of $r$ -- is $\nabla \cdot \mathbf v(r) = \sum_b m_b \frac{\mathbf v_b}{\rho_b} \cdot \nabla W(r-r_b,h)$. The same can of course be done to the curl: $\nabla \times \mathbf v(r) = \sum_b m_b \nabla \times \left[\frac{\mathbf v_b}{\rho_b} W(r-r_b,h)\right] = \sum_b m_b \left[\frac{\mathbf v_b}{\rho_b} \times \nabla W(r-r_b,h)\right]$.

• When you say interpolate between the points $r_b$ in the neighborhood you mean interpolate between $r$ and all $r_b$ in the neighborhood? – BRabbit27 Dec 3 '14 at 14:09
• Yes. The phrase "interpolate between the points $r_b$" is intended to mean that we find some function that can be evaluated everywhere (e.g., at the point $r$), not just at the points $r_b$. – Wolfgang Bangerth Dec 5 '14 at 2:36

Your second question is answered by expanding the divergence of $\rho v$ and rearranging terms: $$\nabla\cdot({\rho v}) = \nabla\rho\cdot v + \rho \nabla \cdot v$$ so $$\nabla \cdot v = \frac{\nabla\cdot({\rho v}) - \nabla\rho\cdot v}{\rho}$$

• Ok I see, then the question that's still bugging me is how does the divergence estimator is derived for SPH? – BRabbit27 Nov 29 '14 at 16:20

Given the fundamental idea of the SPH interpolant with an arbitrary function $$A$$ (which could be the velocity $$\textbf{v}$$):

$$A(r)=\int_{-\infty}^{\infty} A(r')\delta(r-r')dr'$$,

the convolution of an arbitrary $$A(r)$$ with the Dirac delta $$\delta$$. Obviously, is not applicable in numerical simulations due to two reasons:

• $$\delta$$ has infinitesimally small radius, and infinite height
• The integral can be evaluated only analytically.

To solve both issues, first we change $$\delta$$ to a mollifier function (kernel) to achieve a finite (or sometimes infinitely large) influence radius:

$$A(r)\approx\int_{\Omega} A(r')W(r-r',h)dr'$$,

where $$h$$ is analogous to the standard deviation in case of Gaussian function and $$\Omega$$ is the domain of the kernel. $$W(r-r',h)$$ is usually a compact bell shaped function of unit volume. Secondly, change the integral to the summation form:

$$\langle A(r_i)\rangle = \sum_{j}{A(r_j)\frac{m_j}{\rho_j}W(r_i-r_j,h)}$$.

Now $$r_i$$ is corresponding to the $$i^{th}$$ discrete point often referred to as particle. This is can be considered as a discrete convolution.

Now apply the continuum convolution to the derivative of A:

$$\nabla_{r} A(r)=\int_{\Omega} \nabla_{r'} A(r')W(r-r',h)dr'$$,

where $$\nabla_{r}$$ denotes the derivative at $$r$$.

Using the rules of differentiation, it also can be written as:

$$\nabla_{r} A(r)\approx\int_{\Omega} \nabla_{r'} \big[A(r')W(r-r',h)\big]dr'-\int_{\Omega}A(r')\nabla_{r'} W(r-r',h)dr'$$.

Transform the first term in the RHS from volume integral to surface integral using the Gauss-Ostrogradsky theorem:

$$\nabla A(r)\approx\int_{\partial\Omega}A(r')W(r-r',h)dS-\int_{\Omega}A(r')\nabla_{r'} W(r-r',h)dr'$$.

Since $$W(r-r',h)$$ is compact and has a zero value at its boundary, the first term on RHS is now zero, so:

$$\nabla A(r)\approx-\int_{\Omega}A(r')\nabla_{r'} W(r-r',h)dr'$$,

where we can use the identity

$$\nabla_{r'}W(r-r',h)=-\nabla_{r}W(r-r',h)$$,

to obtain the discrete form

$$\langle \nabla A_i\rangle = \sum_{j}{A_j\frac{m_j}{\rho_j}\nabla_i W(r_i-r_j,h)}$$.

Unfortunately, this gives a poor estimation on the derivative, because it is extremely sensitive to the particle distribution. To overcome this problem, you can apply the discrete convolution to the formula you've shown:

$$\langle\nabla A_i\rangle=\frac{\langle \nabla(\rho_i A_i)\rangle-A_i\langle\nabla \rho_i\rangle}{\rho_i} = \frac{1}{\rho_i}\sum_{j}{(\rho_jA_j)\frac{m_j}{\rho_j}\nabla_i W(r_i-r_j,h)}-\frac{1}{\rho_i}\sum_{j}{(\rho_jA_i)\frac{m_j}{\rho_j}\nabla_i W(r_i-r_j,h)}$$,

which is finally:

$$\langle\nabla A_i\rangle=\frac{1}{\rho_i}\sum_{j}{(A_j-A_i)m_j \nabla_i W(r_i-r_j,h)}$$.

It is important, that this is only one of the several possible forms of differential operators in SPH. They have different consistency and conservation properties forming a trade-off. Which to use? It depends on the equation you want to solve.

Nevertheless, you are free to change the product to use the obtained operator for curl:

$$\langle\nabla\times \textbf{v}_i\rangle=\frac{1}{\rho_i}\sum_{j}{(\textbf{v}_j-\textbf{v}_i)\times\nabla_i W(r_i-r_j,h)m_j}$$

• A question: where did the minus sign go in the discretization after the application of the Gauss-Ostrogradsky theorem? – Michele Feb 24 '20 at 14:44
• The negative sign comes from the rule of differentiation of function products below the line starting with "Using the rules of differentiation,...". – BalazsToth Feb 24 '20 at 19:42
• That I know, but why did it disappear after the line "which can be discretized:"? – Michele Feb 25 '20 at 11:33
• Then this is a good question, thank you for noticing. I edited the answer to clarify this issue. – BalazsToth Feb 27 '20 at 9:23