Skip to content

Instantly share code, notes, and snippets.

@kishida

kishida/ParticleParallel.java

Created April 30, 2023 02:57
Show Gist options
  • Save kishida/e0330adeebc750987e3722d3e252b048 to your computer and use it in GitHub Desktop.
Save kishida/e0330adeebc750987e3722d3e252b048 to your computer and use it in GitHub Desktop.
Download ZIP
並列化した粒子法での流体シミュレーション
Raw
ParticleParallel.java
import java.awt.*;
import java.awt.image.BufferedImage;
import java.util.stream.IntStream;
import javax.swing.*;
public class ParticleParallel {
enum Type{
WATER(new Color(0, 0, 255, 96)), OUTWALL(Color.BLACK), INWALL(Color.GREEN);
Color color;
Type(Color c){
color = c;
}
}
static class Elem{
double x;
double y;
double vx;
double vy;
Type type;
boolean surface = false;
public Elem(double x, double y, Type type) {
this.x = x;
this.y = y;
this.type = type;
}
Color getColor(){
return surface ? new Color(255, 255, 255, 96) : type.color;
}
}
static Elem[] elements = new Elem[1500];
static int elementCount;
static JLabel label;
static Image img = new BufferedImage(500, 400, BufferedImage.TYPE_INT_RGB);
static final double parcSize = 8e-3;//1.5e-2;//粒子サイズ
static final double minlen = 0.5;//粒子の最低距離(粒子サイズに対する)
static final double re = 2.1;//近傍粒子の範囲
static final double relap = 4.0;//ラプラシアン用の近傍粒子の範囲
static double deltat = 0.0001;//時間の単位(粒子の最高速度によって調整される)
static final double grav = 9.8;//重力加速度(m/s^2)
static final double dim = 2;//次元数
static final double nyu = 1.0038E-6;//動粘性係数(水 1.0038E-6 m^2/s)
static final double rho = 1.0e3;//粒子の比重(kg/m^3)
static final double deltatmax = 1.0e-4;//時間の単位の最大値(本来は1e-3)
static double n0;//基準粒子密度
static double lambda;//ラプラシアン用の係数
static final double beta = 0.97;//自由表面条件
static double time = 0;
static final int PARALLEL = 8;
public static void main(String[] args){
//左右の外壁
double maxy = 0;
double oldparc = 1.5e-2;
double ysize = 0.5 * parcSize / oldparc;
double xsize = 0.8 * parcSize / oldparc;
for(double y = 0; y < ysize; y += parcSize){
elements[elementCount++] = new Elem(0, y, Type.INWALL);
elements[elementCount++] = new Elem(-parcSize, y, Type.OUTWALL);
elements[elementCount++] = new Elem(-parcSize * 2, y, Type.OUTWALL);
elements[elementCount++] = new Elem(xsize, y, Type.INWALL);
elements[elementCount++] = new Elem(xsize+parcSize, y, Type.OUTWALL);
elements[elementCount++] = new Elem(xsize+parcSize * 2, y, Type.OUTWALL);
maxy = y;
}
//下の壁
for(double x = -parcSize * 2; x < xsize + parcSize * 2; x += parcSize){
elements[elementCount++] = new Elem(x, maxy + parcSize,
(x < 0 || x > xsize) ? Type.OUTWALL : Type.INWALL);
elements[elementCount++] = new Elem(x, maxy + parcSize * 2, Type.OUTWALL);
elements[elementCount++] = new Elem(x, maxy + parcSize * 3, Type.OUTWALL);
}
//水柱
int sample = 0;
double waterx = .3 * parcSize / oldparc;
double watercxmin = .2 * parcSize / oldparc;
double watercxmax = .25 * parcSize / oldparc;
double watercymin = .1 * parcSize / oldparc;
double watercymax = .15 * parcSize / oldparc;
for(double x = parcSize; x < waterx; x+= parcSize){
for(double y = parcSize * 3; y < ysize; y += parcSize){//.5まで
if(y > watercxmin && y < watercxmax && x > watercymin && x < watercymax){
//内部になる点を一点記憶しておく
sample = elementCount;
}
elements[elementCount++] = new Elem(x, y, Type.WATER);
}
}
//基本の粒子密度n0とラプラシアン用の係数rambda
n0 = 0;
double a = 0;
for(int i = 0; i < elementCount; ++i){
if(i == sample) continue;
double d = Math.sqrt(
(elements[sample].x - elements[i].x) * (elements[sample].x - elements[i].x) +
(elements[sample].y - elements[i].y) * (elements[sample].y - elements[i].y));
double w = (d < re * parcSize) ? parcSize * re / d - 1 : 0;
n0 += w;
a += (elements[sample].x - elements[i].x) * (elements[sample].x - elements[i].x) +
(elements[sample].y - elements[i].y) * (elements[sample].y - elements[i].y) * w;
}
lambda = a / n0;
//ウィンドウの表示
JFrame f = new JFrame("粒子法");
f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
f.setLocation(400, 300);
f.setSize(500, 450);
label = new JLabel(new ImageIcon(img));
draw(1000000);
f.add(label);
f.setVisible(true);
//処理スレッド
new Thread(() -> {
try{
for(;;){
long start = System.nanoTime();
calc();
draw(System.nanoTime() - start);
}
}catch(InterruptedException e){}
}).start();
}
static void calc() throws InterruptedException{
final double[][] w = new double[elementCount][elementCount];
//粒子密度を計算しておく
IntStream.range(0, PARALLEL).parallel().forEach(p -> {
int range = Math.min(elementCount * (p + 1) / PARALLEL, elementCount);
for(int i = p * elementCount / PARALLEL; i < range; ++i){
for(int j = 0; j < i; ++j){
double d = Math.sqrt(
(elements[i].x - elements[j].x) * (elements[i].x - elements[j].x) +
(elements[i].y - elements[j].y) * (elements[i].y - elements[j].y));
double wd = (d < relap * parcSize) ? parcSize * relap / d - 1 : 0;
w[i][j] = wd;
w[j][i] = wd;
}
}
});
//時間間隔
double velolim = 0.2 * parcSize / (deltatmax * 0.02);
double svelolim = velolim * velolim;
double maxvelo = IntStream.range(0, PARALLEL).parallel().mapToDouble(p -> {
double maxvtemp = 0;
int range = Math.min(elementCount * (p + 1) / PARALLEL, elementCount);
for(int i = p * elementCount / PARALLEL; i < range; ++i){
if(elements[i].type != Type.WATER) continue;
double v = elements[i].vx * elements[i].vx + elements[i].vy * elements[i].vy;
if(v > velolim){
//速すぎる場合、遅くしちゃう
double sv = Math.sqrt(v);
elements[i].vx *= svelolim / sv;
elements[i].vy *= svelolim / sv;
v = velolim;
}
if(v > maxvtemp) maxvtemp = v;
}
return maxvtemp;
}).max().getAsDouble();
maxvelo = Math.sqrt(maxvelo);
double dtlim = deltat * 1.1;
if(maxvelo == 0){
deltat = dtlim;
}else{
deltat = 0.2 * parcSize / maxvelo;
}
if(deltat > dtlim) deltat = dtlim;
if(deltat > deltatmax) deltat = deltatmax;
//仮の速度・位置を計算
final double oldx[] = new double[elementCount];
final double oldy[] = new double[elementCount];
for(int i = 0; i < elementCount; ++i){
oldx[i] = elements[i].vx;
oldy[i] = elements[i].vy;
}
IntStream.range(0, PARALLEL).parallel().forEach(p -> {
int range = Math.min(elementCount * (p + 1) / PARALLEL, elementCount);
for(int i = p * elementCount / PARALLEL; i < range; ++i){
if(elements[i].type != Type.WATER){
continue;
}
//速度のラプラシアンを計算
double laplasux = 0;
double laplasuy = 0;
for(int j = 0; j < elementCount; ++j){
if(j == i) continue;
if(w[i][j] == 0) continue;
laplasux += (oldx[j] - elements[i].vx) * w[i][j];
laplasuy += (oldy[j] - elements[i].vy) * w[i][j];
}
laplasux *= 2 * dim / n0 / lambda;
laplasuy *= 2 * dim / n0 / lambda;
//仮の位置・速度を計算
//
elements[i].vx += + deltat * (nyu * laplasux);
elements[i].vy += deltat * (nyu * laplasuy + grav);
elements[i].x += elements[i].vx * deltat;
elements[i].y += elements[i].vy * deltat;
}
});
//粒子が近づきすぎてないか
double r2lim = parcSize * parcSize * minlen * minlen;
for(int i = 0; i < elementCount; ++i){
for(int j = 0; j < i; ++j){
if(w[i][j] == 0) continue;
double x = elements[j].x - elements[i].x;
double y = elements[j].y - elements[i].y;
double rr = x * x + y * y;
if(rr < r2lim){
double r = Math.sqrt(rr);
//double m2 = rho + rho;
double vgx = (elements[i].vx + elements[j].vx) / 2;//本来は速度を粘度で按分
double vgy = (elements[i].vy + elements[j].vy) / 2;
double vrx = rho * (elements[i].vx - vgx);
double vry = rho * (elements[i].vy - vgy);
double vabs = (vrx * x + vry * y) / r;
if(vabs < 0) continue;
final double colrat = 0.15;
double vrat = 1 + colrat;
double vmx = vrat * vabs * x / r;
double vmy = vrat * vabs * y / r;
if(elements[i].type == Type.WATER){
elements[i].vx -= vmx / rho;
elements[i].vy -= vmy / rho;
elements[i].x -= deltat * vmx / rho;
elements[i].y -= deltat * vmy / rho;
}
if(elements[j].type == Type.WATER){
elements[j].vx += vmx / rho;
elements[j].vy += vmy / rho;
elements[j].x += deltat * vmx / rho;
elements[j].y += deltat * vmy / rho;
}
}
}
}
//仮の粒子密度を計算
final double[][] wasta = new double[elementCount][elementCount];
final double[][] rsqu = new double[elementCount][elementCount];
final int[] neighborcount = new int[elementCount];
final int[][] neighbor = new int[elementCount][elementCount];
for(int i = 0; i < elementCount; ++i){
for(int j = 0; j < i; ++j){
rsqu[i][j] = rsqu[j][i] =
(elements[i].x - elements[j].x) * (elements[i].x - elements[j].x) +
(elements[i].y - elements[j].y) * (elements[i].y - elements[j].y);
double d = Math.sqrt(rsqu[i][j]);
double rp = re * parcSize;
if(d < rp){
double wd = rp / d - 1;
wasta[i][j] = wd;
wasta[j][i] = wd;
neighbor[i][neighborcount[i]++] = j;
neighbor[j][neighborcount[j]++] = i;
}
}
}
//粒子密度から自由表面になるところを検出
final double[] nasta = new double[elementCount];
for(int i = 0; i < elementCount; ++i){
for(int j = 0; j < neighborcount[i]; ++j){
nasta[i] += wasta[i][neighbor[i][j]];
}
if(elements[i].type == Type.WATER){
elements[i].surface = nasta[i] < n0 * beta;
}
}
//よくわからんけど。
final int[] check = new int[elementCount];
for(int i = 0; i < elementCount; ++i){
if(elements[i].type != Type.WATER){
check[i] = 2;
}else if(elements[i].surface){
check[i] = 1;
}else{
check[i] = 0;
}
}
for(boolean cont = true; !cont;){
cont = false;
for(int i = 0; i < elementCount; ++i){
if(check[i] == 1){
for(int j = 0; j < neighborcount[i]; ++j){
if(check[neighbor[i][j]] == 0) check[neighbor[i][j]] = 1;
}
check[i] = 2;
cont = true;
}
}
}
//連立方程式の係数行列を作成
final double[] rk = new double[elementCount];
final double[][] coef = new double[elementCount][elementCount];
final double[] val = new double[elementCount];
//exec = Executors.newFixedThreadPool(threadLimit);
for(int i = 0; i < elementCount; ++i){
if(elements[i].surface){
val[i] = 0;
}else{
val[i] = -rho / deltat / deltat * (nasta[i] - n0) * lambda / 2 / dim;
}
for(int j = 0; j < neighborcount[i]; ++j){
coef[i][j + 1] = wasta[i][neighbor[i][j]];
coef[i][0] -= wasta[i][neighbor[i][j]];
}
coef[i][0] -= 1.0e-7 / deltat / deltat * rho * n0;//よくわからないけどソースから
if(check[i] == 0) coef[i][0] *= 2;
}
//連立方程式を解く(とりあえずガウスザイデル法。CG法がいいらしい)
for(int n = 0; n < 1000; ++n){
double div = 0;
for(int i = 0; i < elementCount; ++i){
double sum = 0;
for(int j = 0; j < neighborcount[i]; ++j){
sum += coef[i][j + 1] * rk[neighbor[i][j]];
}
double old = rk[i];
rk[i] = (val[i] - sum) / coef[i][0];
div += Math.abs(rk[i] - old);
}
//収束判定
if(div < 0.0001) {
break;
}
}
for(int i = 0; i < elementCount; ++i){
//圧力が負の場合は0に
if(rk[i] < 0) rk[i] = 0;
}
//速度・場所を再計算
for(int i = 0; i < elementCount; ++i){
oldx[i] = elements[i].x;
oldy[i] = elements[i].y;
}
IntStream.range(0, PARALLEL).parallel().forEach(p -> {
int range = Math.min(elementCount * (p + 1) / PARALLEL, elementCount);
for(int i = p * elementCount / PARALLEL; i < range; ++i){
// final int i = idx;
if(elements[i].type != Type.WATER) continue;//場所がかわるのは水だけ
double udashx = 0;
double udashy = 0;
//圧力の最小値を得る
double hatp = rk[i];
for(int j = 0; j < neighborcount[i]; ++j){
if(hatp > rk[neighbor[i][j]]) hatp = rk[neighbor[i][j]];
}
//速度の修正量を計算
for(int j = 0; j < neighborcount[i]; ++j){
int n = neighbor[i][j];
udashx += (rk[n] - hatp) / rsqu[i][n] * (oldx[n] - elements[i].x) * wasta[i][n];
udashy += (rk[n] - hatp) / rsqu[i][n] * (oldy[n] - elements[i].y) * wasta[i][n];
}
udashx = -udashx * deltat / rho * dim / n0;
udashy = -udashy * deltat / rho * dim / n0;
//速度の再計算
elements[i].vx += udashx;
elements[i].vy += udashy;
//場所を再計算
elements[i].x += deltat * udashx;
elements[i].y += deltat * udashy;
}
});
time += deltat;
}
static void draw(long delta){
Graphics2D g = (Graphics2D) img.getGraphics();
g.setColor(new Color(128, 212, 255));
g.fillRect(0, 0, 500, 400);
int size = 8;
double s = size / parcSize;
for(int i = 0; i < elementCount; ++i){
Elem el = elements[i];
g.setColor(el.getColor());
g.fillOval((int)(el.x * s) + 30, (int)(el.y * s) + 20, size * 3/2, size * 3/2);
}
g.setColor(Color.BLACK);
g.drawString(String.format("%2.4f (%1.4fs/m)", time, deltat / delta * 1000 * 1000 * 1000 * 60), 10, 15);
g.dispose();
label.repaint();
}
}
@kishida
Copy link
Author

kishida commented Apr 30, 2023

10x
https://user-images.githubusercontent.com/1178746/235333375-5c4a1bf9-a623-453a-b8d6-e8892cfdbdbb.mp4

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment