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?

  • $\begingroup$ 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. $\endgroup$
    – Bill Barth
    Commented Nov 29, 2014 at 16:50
  • $\begingroup$ @BillBarth I was indeed missing a nabla, could you explain me the derivation? $\endgroup$
    – BRabbit27
    Commented Nov 29, 2014 at 16:55

3 Answers 3


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]$.

  • $\begingroup$ When you say interpolate between the points $r_b$ in the neighborhood you mean interpolate between $r$ and all $r_b$ in the neighborhood? $\endgroup$
    – BRabbit27
    Commented Dec 3, 2014 at 14:09
  • $\begingroup$ 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$. $\endgroup$ Commented Dec 5, 2014 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} $$

  • $\begingroup$ Ok I see, then the question that's still bugging me is how does the divergence estimator is derived for SPH? $\endgroup$
    – BRabbit27
    Commented Nov 29, 2014 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


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}$

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.