# implementation of a simple SPH example

I am currently learning the Smoothed Particle Hydrodynamics method that I need to use later in my thesis. I have studied the mathematics behind the method and I want to code a simple example to show how the classical SPH does not reproduce a constant field exactly. So, if we consider for example

$$f(x,y)=\frac{1}{2}\left(\frac{\partial x}{\partial x}+\frac{\partial y}{\partial y}\right)=1$$

I would like to know how to get the following approximation of $f$

I am familiar with Fortran, so I would prefer to use it for this example.

Any help with this would be appreciated and thank you in advance.

• I think that your question is not clear enough. Are you asking how to modify the classical SPH to obtain a constant field? – nicoguaro Feb 1 '18 at 21:18
• No, i am just asking how to create a code that produce the linked image (to show that sph does not reproduce constant field exactly) i don't know how to do it – hamza boulahia Feb 1 '18 at 21:20
• So, you are asking how to program the method? – nicoguaro Feb 1 '18 at 21:21
• Yes if you want, but even some instruction to start wloud be very helpfull – hamza boulahia Feb 1 '18 at 21:23
• Have you read an introductory text on the topic? – nicoguaro Feb 1 '18 at 21:31

$f(x,y)=\frac{1}{2}\big(\frac{\partial x}{\partial x}+\frac{\partial y}{\partial y}\big)=\frac{1}{2}\nabla\cdot\textbf{r}\approx\frac{1}{2}\sum_j^n(r_j-r_i)\frac{m_j}{\rho_j}\nabla W_{ij}=\langle f_i\rangle.\qquad(*)$

A very useful simulation tool could be the open-source Aboria, constructed to facilitate the application of particle methods like SPH. Although it requires C++ programming skills, this is a good choice for your problem. You can find Aboria on GitHub here and a related SPH example here.

As long as programming is not a requirement, probably the easiest way to investigate this and other properties of SPH is using Nauticle, which is also an open-source tool with objectives similar to that of Aboria. The Nauticle code is available here with examples and user-guide. However, this tool is not an option if you prefer programming SPH manually.

Running Nauticle with the following simple configuration file

simulation:
case:
workspace:
constants:
- L: 1
- rho0: 1000
- csize: L/2
- dx: csize/2.5
- mass: dx^2*rho0
variables:
- dt: 1
particle_system:
domain:
cell_size: csize|csize
minimum: 0|0
maximum: L/csize|L/csize
boundary: 2|2
grid:
gpos: 0|0
gsize: L|L
goffset: 0|0
gip_dist: dx|dx
fields:
- f: 0
equations:
- eq: f=sph_D00(r,mass,rho0,Wp52220,csize)/2
parameter_space:
simulated_time: dt
print_interval: dt


gives the result you requested. The equation in the fourth line from below implements the SPH differential operator $(*)$

Visualisation of the result in Paraview:

• Any suggested documentations on how to visualize SPH results on Paraview ? – hamza boulahia Mar 2 '18 at 16:48
• Most of the SPH simulation tools generate .vtk result files, which can be simply read in Paraview. You can read about the usage through the link below. This material is not related to particle data but the procedure is similar. cac.cornell.edu/education/Training/data10/… – BalazsToth Mar 2 '18 at 21:21

Here is a sample SPH code, a version of dam break in Python, in 3d. To run, pip install PyOpenGL==3.1.0 package, and on Ubuntu

sudo apt-get install mesa-common-dev libgl1-mesa-dev libglu1-mesa-dev freeglut3-dev


The simulation code based on this repo

https://github.com/cerrno/mueller-sph

from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *
from random import random
from PIL import Image
from PIL import ImageOps
from collections import defaultdict
import numpy as np, datetime
import sys, numpy.linalg as lin

p1,p2,p3 = 73856093, 19349663, 83492791
G = np.array([0.0, 0.0, -9.8*2])

B = 10 # top
l = 0.2 # bolec kutu buyuklugu
n = B*20 # bolec sozluk buyuklugu

REST_DENS = 10.0
GAS_CONST = 0.5
MASS = 100.0
VISC = 20.0
DT = 0.1
H = 0.2 # kernel radius
HSQ = H*H # radius^2 for optimization
POLY6 = 315.0/(65.0*np.pi*np.power(H, 9.));
VISC_LAP = 45.0/(np.pi*np.power(H, 6.));
EPS = 0.05
BOUND_DAMPING = -0.5
img = True

def spatial_hash(x):
"""
x = [x0,x1,x2] uc boyutlu kordinatlari icin bir bolec (hash) degeri uret
"""
ix,iy,iz = np.floor((x[0]+2.0)/l), np.floor((x[1]+2.0)/l), np.floor((x[2]+2.0)/l)
return (int(ix*p1) ^ int(iy*p2) ^ int(iz*p3)) % n

class Simulation:
def __init__(self):
self.geo_hash_list = None
self.i = 0
self.r   = 0.05
self.cor = 0.5
self.balls = []
self.tm  = 0.0
self.th  = 200.0
self.mmax =  1.0-self.r
self.mmin = -1.0+self.r
self.right = False
self.left = False

def init(self):
i = 0
for xs in np.linspace(-0.3, 0.3, 10):
for ys in np.linspace(-0.3, 0.3, 10):
for zs in np.linspace(0.0, 0.2, 4):
v = np.array([0.0, 0.0, 0.0])
f = np.array([0,0,0])
x = np.array([xs, ys, zs])
d = {'x': x, 'f':f, 'v': v, 'i': i, 'rho': 0.0, 'p': 0.0}
self.balls.append(d)
i += 1

tm = 0.0

glEnable(GL_LIGHTING)
glEnable(GL_LIGHT0)
glEnable(GL_DEPTH_TEST)
glClearColor(1.0,1.0,1.0,1.0)

glMatrixMode(GL_PROJECTION)
gluPerspective(60.0,1.0,1.0,50.0)
glTranslatef(0.0,0.0,-3.5)
glMatrixMode(GL_MODELVIEW)

def hash_balls(self):
self.geo_hash_list = defaultdict(list)
for i,b in enumerate(self.balls):
self.geo_hash_list[spatial_hash(self.balls[i]['x'])].append(self.balls[i])

def computeDensityPressure(self):
for i,pi in enumerate(self.balls):
pi['rho'] = 0.0
h = spatial_hash(self.balls[i]['x']) # su anki topun boleci
if (len(self.geo_hash_list[h])>1): # yakinda top var mi
otherList = self.geo_hash_list[h] # varsa isle
for j,pj in enumerate(otherList):
r2 = lin.norm(pj['x']-pi['x'])**2
if  r2 < HSQ:
pi['rho'] += MASS*POLY6*np.power(HSQ-r2, 3.0)
pi['p'] = GAS_CONST*(pi['rho'] - REST_DENS)

def computeForces(self):
for i,pi in enumerate(self.balls):
fpress = np.array([0.0, 0.0, 0.0])
fvisc = np.array([0.0, 0.0, 0.0])
h = spatial_hash(self.balls[i]['x']) # su anki topun boleci
if (len(self.geo_hash_list[h])>1): # yakinda top var mi
otherList = self.geo_hash_list[h] # varsa isle
for j,pj in enumerate(otherList):
if pj['i'] == pi['i']: continue
rij = pi['x']-pj['x']
r = lin.norm(rij)
if r < H:
if np.sum(rij)>0.0: rij = rij / r
tmp1 = -rij*MASS*(pi['p'] + pj['p']) / (2.0 * pj['rho'])
fpress += (tmp1 * tmp2)
tmp1 = VISC*MASS*(pj['v'] - pi['v'])
tmp2 = pj['rho'] * VISC_LAP*(H-r)
fvisc += (tmp1 / tmp2)
fgrav = G * pi['rho']
pi['f'] = fpress + fvisc + fgrav

def integrate(self):
for j,p in enumerate(self.balls):
if p['rho'] > 0.0:
p['v'] += DT*p['f']/p['rho']
p['x'] += DT*p['v']

if p['x'][0]-EPS < -1.0:
p['v'][0] *= BOUND_DAMPING
p['x'][0] = -1.0
if p['x'][0]+EPS > 1.0:
p['v'][0] *= BOUND_DAMPING
p['x'][0] = 1.0-EPS

if p['x'][1]-EPS < -1.0:
p['v'][1] *= BOUND_DAMPING
p['x'][1] = -1.0
if p['x'][1]+EPS > 1.0:
p['v'][1] *= BOUND_DAMPING
p['x'][1] = 1.0-EPS

if p['x'][2]-EPS < -1.0:
p['v'][2] *= BOUND_DAMPING
p['x'][2] = -1.0
if p['x'][2]+EPS > 1.0:
p['v'][2] *= BOUND_DAMPING
p['x'][2] = 1.0-EPS

self.hash_balls()

def update(self):
self.hash_balls()
self.computeDensityPressure()
self.computeForces()
self.integrate()
if self.i > 40: exit()
glutPostRedisplay()

def display(self):
glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
glPushMatrix()
glRotatef(self.th,0.0,1.0,0.0)
glRotatef(90.0,-1.0,0.0,0.0)
glutWireCube(2.0)
for j,b in enumerate(self.balls):
glPushMatrix()
glTranslatef(b['x'][0],b['x'][1],b['x'][2])
glMaterialfv(GL_FRONT, GL_DIFFUSE, [0.0, 0.0, 1.0, 1.0])
glutSolidSphere(self.r,50,50)
glPopMatrix()
glPopMatrix()
glutSwapBuffers()

if img and self.i % 2 == 0:
width,height = 480,480
data = glReadPixels(0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE)
image = Image.frombytes("RGBA", (width, height), data)
image = ImageOps.flip(image)
image.save('/tmp/glut/glutout-%03d.png' % self.i, 'PNG')

self.i += 1

if __name__ == '__main__':
if (os.path.exists("/tmp/glut") == False): os.mkdir("/tmp/glut")
s = Simulation()
glutInit(())
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH)
glutInitWindowSize(500,500)
glutCreateWindow("GLUT Bouncing Ball in Python")
glutDisplayFunc(s.display)
glutIdleFunc(s.update)
s.init()
glutMainLoop()

$$$$
`