我目前正在学习我需要在我的论文后面使用的平滑粒子流体动力学方法。我已经研究了该方法背后的数学，我想编写一个简单的示例来说明经典 SPH 如何不精确地再现常量字段。所以，如果我们考虑例如
我想知道如何获得以下近似值
我对 Fortran 很熟悉，所以我更愿意在这个例子中使用它。
对此的任何帮助将不胜感激，并在此先感谢您。
一个简单的 SPH 示例的实现
计算科学
正则
sph
2021-12-15 19:29:31
2个回答
您问题中函数的 SPH 实现如下：
一个非常有用的模拟工具可能是开源的Aboria，它的构建是为了促进像 SPH 这样的粒子方法的应用。尽管它需要 C++ 编程技能，但这对于您的问题来说是一个不错的选择。您可以在 GitHub 上找到 Aboria，在此处找到相关的 SPH 示例。
只要不需要编程，研究 SPH 的这个和其他属性的最简单方法可能是使用Nauicle，它也是一个开源工具，其目标类似于 Aboria。Nauicle 代码可在此处获得示例和用户指南。但是，如果您更喜欢手动编程 SPH，则不能选择此工具。
使用以下简单配置文件运行 Nauicle
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
给出你要求的结果。下面第四行的方程实现了 SPH 微分算子
Paraview 中结果的可视化：
这是一个示例 SPH 代码，它是 Python 中的 dam break 的一个版本，在 3d 中。要运行，请 pip 安装PyOpenGL==3.1.0包，并在 Ubuntu 上运行
sudo apt-get install mesa-common-dev libgl1-mesa-dev libglu1-mesa-dev freeglut3-dev
基于这个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()
```
其它你可能感兴趣的问题