227 lines
7.0 KiB
Python
227 lines
7.0 KiB
Python
import sys
|
|
import numpy as np
|
|
import cupy as cp
|
|
import pygame
|
|
|
|
# Minimaler, GPU-beschleunigter 2D-Fluid-Solver mit persistenter Flüssigkeit und Gravitation.
|
|
# Installation:
|
|
# pip install cupy-cudaXX pygame numpy
|
|
# XX = passende CUDA-Version, z.B. cupy-cuda112 für CUDA 11.2
|
|
|
|
# ==============================
|
|
# EINSTELLUNGEN
|
|
# ==============================
|
|
N = 64 # Gittergröße (N x N)
|
|
ITER = 2 # Iterationen; verringern für Geschwindigkeit
|
|
DT = 0.1 # Zeitschritt
|
|
DIFF = 0.0001 # Diffusionskoeffizient (klein, aber nicht null für Stabilität)
|
|
VISC = 0.0001 # Viskosität (klein, aber nicht null für Stabilität)
|
|
ZOOM = 8 # Pixel pro Gitterzelle in der Anzeige
|
|
GRAVITY = 0.1 # Gravitationskraft (anpassbar)
|
|
|
|
WIDTH, HEIGHT = N * ZOOM, N * ZOOM
|
|
|
|
# ==============================
|
|
# HILFSFUNKTIONEN
|
|
# ==============================
|
|
def set_bnd(b, x):
|
|
"""
|
|
Randbehandlung für die Simulation.
|
|
b = 0: Dichte
|
|
b = 1: u-Geschwindigkeit
|
|
b = 2: v-Geschwindigkeit
|
|
"""
|
|
if b == 1:
|
|
# u-Geschwindigkeit: links und rechts auf 0 setzen (Wände)
|
|
x[0, 1:-1] = 0
|
|
x[-1, 1:-1] = 0
|
|
elif b == 2:
|
|
# v-Geschwindigkeit: oben und unten auf 0 setzen (Wände)
|
|
x[1:-1, 0] = 0
|
|
x[1:-1, -1] = 0
|
|
else:
|
|
# Dichte: Spiegeln an den Rändern
|
|
x[1:-1, 0] = x[1:-1, 1]
|
|
x[1:-1, -1] = x[1:-1, -2]
|
|
x[0, 1:-1] = x[1, 1:-1]
|
|
x[-1, 1:-1] = x[-2, 1:-1]
|
|
|
|
# Ecken setzen
|
|
x[0, 0] = 0.5 * (x[1, 0] + x[0, 1])
|
|
x[0, -1] = 0.5 * (x[1, -1] + x[0, -2])
|
|
x[-1, 0] = 0.5 * (x[-2, 0] + x[-1, 1])
|
|
x[-1, -1] = 0.5 * (x[-2, -1] + x[-1, -2])
|
|
|
|
def add_source(x, s):
|
|
"""
|
|
Fügt eine Quelle zur Simulation hinzu.
|
|
"""
|
|
x += DT * s
|
|
|
|
def diffuse(b, x, x0, diff):
|
|
"""
|
|
Diffusionsschritt der Simulation.
|
|
"""
|
|
a = DT * diff * (N - 2) * (N - 2)
|
|
for _ in range(ITER):
|
|
x[1:-1, 1:-1] = (x0[1:-1, 1:-1] + a * (
|
|
x[0:-2, 1:-1] + x[2:, 1:-1] +
|
|
x[1:-1, 0:-2] + x[1:-1, 2:]
|
|
)) / (1 + 4*a)
|
|
set_bnd(b, x)
|
|
|
|
def advect(b, d, d0, u, v):
|
|
"""
|
|
Advektionsschritt der Simulation.
|
|
"""
|
|
dt0 = DT * (N - 2)
|
|
|
|
# Rückverfolgte Position berechnen
|
|
X, Y = cp.meshgrid(cp.arange(N), cp.arange(N), indexing='ij')
|
|
X = X[1:-1, 1:-1] - dt0 * u[1:-1, 1:-1]
|
|
Y = Y[1:-1, 1:-1] - dt0 * v[1:-1, 1:-1]
|
|
|
|
X = cp.clip(X, 0.5, N - 1.5)
|
|
Y = cp.clip(Y, 0.5, N - 1.5)
|
|
|
|
i0 = X.astype(int)
|
|
i1 = i0 + 1
|
|
j0 = Y.astype(int)
|
|
j1 = j0 + 1
|
|
|
|
s1 = X - i0
|
|
s0 = 1 - s1
|
|
t1 = Y - j0
|
|
t0 = 1 - t1
|
|
|
|
d_new = (s0 * (t0 * d0[i0, j0] + t1 * d0[i0, j1]) +
|
|
s1 * (t0 * d0[i1, j0] + t1 * d0[i1, j1]))
|
|
|
|
d[1:-1, 1:-1] = d_new
|
|
set_bnd(b, d)
|
|
|
|
def project(u, v, p, div):
|
|
"""
|
|
Projektionsschritt zur Inkompressibilität der Flüssigkeit.
|
|
"""
|
|
# Divergenz berechnen
|
|
div[1:-1,1:-1] = -0.5 * (
|
|
(u[2:,1:-1] - u[0:-2,1:-1]) +
|
|
(v[1:-1,2:] - v[1:-1,0:-2])
|
|
) / N
|
|
p.fill(0)
|
|
set_bnd(0, div)
|
|
set_bnd(0, p)
|
|
|
|
# Poisson-Gleichung lösen
|
|
for _ in range(ITER):
|
|
p[1:-1,1:-1] = (div[1:-1,1:-1] + p[0:-2,1:-1] + p[2:,1:-1] +
|
|
p[1:-1,0:-2] + p[1:-1,2:]) / 4
|
|
set_bnd(0, p)
|
|
|
|
# Geschwindigkeiten aktualisieren
|
|
u[1:-1,1:-1] -= 0.5 * (p[2:,1:-1] - p[0:-2,1:-1]) * N
|
|
v[1:-1,1:-1] -= 0.5 * (p[1:-1,2:] - p[1:-1,0:-2]) * N
|
|
set_bnd(1, u)
|
|
set_bnd(2, v)
|
|
|
|
# ==============================
|
|
# FLUID-KLASSE
|
|
# ==============================
|
|
class Fluid:
|
|
def __init__(self):
|
|
self.dens = cp.zeros((N, N), dtype=cp.float32)
|
|
self.dens_prev = cp.zeros((N, N), dtype=cp.float32)
|
|
self.u = cp.zeros((N, N), dtype=cp.float32)
|
|
self.u_prev = cp.zeros((N, N), dtype=cp.float32)
|
|
self.v = cp.zeros((N, N), dtype=cp.float32)
|
|
self.v_prev = cp.zeros((N, N), dtype=cp.float32)
|
|
|
|
def step(self, gravity):
|
|
"""
|
|
Führt einen Simulationsschritt durch, einschließlich der Anwendung von Gravitation.
|
|
"""
|
|
# Gravitation hinzufügen (nur zur v-Komponente)
|
|
self.v_prev += gravity * DT
|
|
|
|
# Dichte
|
|
add_source(self.dens, self.dens_prev)
|
|
diffuse(0, self.dens, self.dens_prev, DIFF)
|
|
advect(0, self.dens, self.dens_prev, self.u, self.v)
|
|
|
|
# Geschwindigkeit
|
|
add_source(self.u, self.u_prev)
|
|
add_source(self.v, self.v_prev)
|
|
diffuse(1, self.u, self.u_prev, VISC)
|
|
diffuse(2, self.v, self.v_prev, VISC)
|
|
project(self.u, self.v, self.u_prev, self.v_prev)
|
|
|
|
advect(1, self.u, self.u_prev, self.u, self.v)
|
|
advect(2, self.v, self.v_prev, self.u, self.v)
|
|
project(self.u, self.v, self.u_prev, self.v_prev)
|
|
|
|
# Quellen zurücksetzen
|
|
self.dens_prev.fill(0)
|
|
self.u_prev.fill(0)
|
|
self.v_prev.fill(0)
|
|
|
|
# ==============================
|
|
# HAUPTPROGRAMM
|
|
# ==============================
|
|
def main():
|
|
pygame.init()
|
|
screen = pygame.display.set_mode((WIDTH, HEIGHT))
|
|
pygame.display.set_caption("GPU Fluid Simulation mit Gravitation")
|
|
clock = pygame.time.Clock()
|
|
|
|
fluid = Fluid()
|
|
|
|
# Zentrum für kontinuierliche Dichtequelle
|
|
center = (N//2, N//2)
|
|
|
|
running = True
|
|
while running:
|
|
clock.tick(60) # 60 FPS anvisiert
|
|
|
|
for event in pygame.event.get():
|
|
if event.type == pygame.QUIT:
|
|
running = False
|
|
|
|
# Mausinteraktion: Linke Maustaste fügt Dichte und zufällige Velocity hinzu
|
|
if pygame.mouse.get_pressed()[0]:
|
|
mx, my = pygame.mouse.get_pos()
|
|
i = int(mx // ZOOM)
|
|
j = int(my // ZOOM)
|
|
if 1 <= i < N-1 and 1 <= j < N-1:
|
|
fluid.dens_prev[j, i] += 100
|
|
fluid.u_prev[j, i] += cp.random.uniform(-2, 2)
|
|
fluid.v_prev[j, i] += cp.random.uniform(-2, 2)
|
|
|
|
# Immer ein klein wenig Dichte in der Mitte
|
|
fluid.dens_prev[center[1], center[0]] += 10 # Hinweis: [j, i] entspricht [y, x]
|
|
|
|
# Ein Schritt Fluid-Simulation mit Gravitation
|
|
fluid.step(GRAVITY)
|
|
|
|
# Dichte von GPU holen für Pygame
|
|
dens_host = cp.asnumpy(fluid.dens)
|
|
# In Graustufen umwandeln und skalieren
|
|
dens_normalized = np.clip(dens_host * 255, 0, 255).astype(np.uint8)
|
|
color_img = np.stack((dens_normalized, dens_normalized, dens_normalized), axis=-1)
|
|
|
|
# Skalieren des Bildes auf die Bildschirmgröße
|
|
color_img_scaled = np.repeat(np.repeat(color_img, ZOOM, axis=0), ZOOM, axis=1)
|
|
|
|
# Erstellen eines Pygame-Surfaces aus dem Array
|
|
surface = pygame.surfarray.make_surface(color_img_scaled.swapaxes(0, 1))
|
|
|
|
# Blitten der Surface auf den Bildschirm
|
|
screen.blit(surface, (0, 0))
|
|
pygame.display.flip()
|
|
|
|
pygame.quit()
|
|
sys.exit()
|
|
|
|
if __name__ == "__main__":
|
|
main()
|