1
$\begingroup$

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$

enter image description here

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.

$\endgroup$
11
  • $\begingroup$ I think that your question is not clear enough. Are you asking how to modify the classical SPH to obtain a constant field? $\endgroup$
    – nicoguaro
    Feb 1, 2018 at 21:18
  • $\begingroup$ 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 $\endgroup$
    – Varazda
    Feb 1, 2018 at 21:20
  • 1
    $\begingroup$ So, you are asking how to program the method? $\endgroup$
    – nicoguaro
    Feb 1, 2018 at 21:21
  • $\begingroup$ Yes if you want, but even some instruction to start wloud be very helpfull $\endgroup$
    – Varazda
    Feb 1, 2018 at 21:23
  • 1
    $\begingroup$ Have you read an introductory text on the topic? $\endgroup$
    – nicoguaro
    Feb 1, 2018 at 21:31

2 Answers 2

5
$\begingroup$

The SPH implementation of the function in your question reads as:

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

enter image description here

$\endgroup$
2
  • $\begingroup$ Any suggested documentations on how to visualize SPH results on Paraview ? $\endgroup$
    – Varazda
    Mar 2, 2018 at 16:48
  • $\begingroup$ 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/… $\endgroup$
    – BalazsToth
    Mar 2, 2018 at 21:21
0
$\begingroup$

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.));
SPIKY_GRAD = -45.0/(np.pi*np.power(H, 6.));
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)
        glLoadIdentity()
        gluPerspective(60.0,1.0,1.0,50.0)
        glTranslatef(0.0,0.0,-3.5)
        glMatrixMode(GL_MODELVIEW)
        glLoadIdentity()

    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'])
                        tmp2 = SPIKY_GRAD*np.power(H-r,2.0)
                        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()


```
$\endgroup$

Your Answer

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

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