SIMD support
This commit is contained in:
unknown
parent
7ab4168a5c
commit
072d4f7784
2
pom.xml
2
pom.xml
@ -2,7 +2,7 @@
|
||||
<project xmlns="http://maven.apache.org/POM/4.0.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/xsd/maven-4.0.0.xsd">
|
||||
<modelVersion>4.0.0</modelVersion>
|
||||
<groupId>electrosphere</groupId>
|
||||
<artifactId>sample-cli-java</artifactId>
|
||||
<artifactId>fluid-sim</artifactId>
|
||||
<version>1.0-SNAPSHOT</version>
|
||||
<packaging>jar</packaging>
|
||||
<properties>
|
||||
|
||||
@ -8,7 +8,12 @@ import java.util.LinkedList;
|
||||
import java.util.List;
|
||||
import java.util.Random;
|
||||
import java.util.Set;
|
||||
// import jdk.incubator.vector.FloatVector;
|
||||
|
||||
import jdk.incubator.foreign.MemorySegment;
|
||||
import jdk.incubator.vector.FloatVector;
|
||||
import jdk.incubator.vector.VectorMask;
|
||||
import jdk.incubator.vector.VectorOperators;
|
||||
import jdk.incubator.vector.VectorOperators.Associative;
|
||||
|
||||
import org.joml.Vector2i;
|
||||
import org.lwjgl.BufferUtils;
|
||||
@ -36,16 +41,15 @@ public class FluidSim {
|
||||
|
||||
public static final int DIM = 18;
|
||||
|
||||
// float[][][][] data = new float[4][DIM][DIM][DIM];
|
||||
float[][][] x = new float[DIM][DIM][DIM];
|
||||
float[][][] x0 = new float[DIM][DIM][DIM];
|
||||
float[][][] u = new float[DIM][DIM][DIM];
|
||||
float[][][] v = new float[DIM][DIM][DIM];
|
||||
float[][][] w = new float[DIM][DIM][DIM];
|
||||
float[][][] u0 = new float[DIM][DIM][DIM];
|
||||
float[][][] v0 = new float[DIM][DIM][DIM];
|
||||
float[][][] w0 = new float[DIM][DIM][DIM];
|
||||
// float[][][][] oldData = new float[4][DIM][DIM][DIM];
|
||||
float[] x = new float[DIM * DIM * DIM];
|
||||
float[] x0 = new float[DIM * DIM * DIM];
|
||||
float[] u = new float[DIM * DIM * DIM];
|
||||
float[] v = new float[DIM * DIM * DIM];
|
||||
float[] w = new float[DIM * DIM * DIM];
|
||||
float[] u0 = new float[DIM * DIM * DIM];
|
||||
float[] v0 = new float[DIM * DIM * DIM];
|
||||
float[] w0 = new float[DIM * DIM * DIM];
|
||||
float[] jacobiAltArray = new float[DIM * DIM * DIM];
|
||||
|
||||
static final float DIFFUSION_CONSTANT = 0.0001f;
|
||||
static final float VISCOSITY_CONSTANT = 0.0001f;
|
||||
@ -65,10 +69,10 @@ public class FluidSim {
|
||||
Math.abs(8 - j) < 4 &&
|
||||
Math.abs(10 - k) < 4
|
||||
){
|
||||
x[i][j][k] = 1;
|
||||
u[i][j][k] = rand.nextFloat() * 0.1f;
|
||||
v[i][j][k] = -1f;
|
||||
w[i][j][k] = rand.nextFloat() * 0.1f;
|
||||
x[IX(i,j,k)] = 1;
|
||||
u[IX(i,j,k)] = rand.nextFloat() * 0.1f;
|
||||
v[IX(i,j,k)] = -1f;
|
||||
w[IX(i,j,k)] = rand.nextFloat() * 0.1f;
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -108,13 +112,12 @@ public class FluidSim {
|
||||
if(step % 1500 == 0){
|
||||
for(int i = 0; i < 4; i++){
|
||||
for(int j = 0; j < 4; j++){
|
||||
x[15 - i][15][15 - j] = 1;
|
||||
v[15 - i][15][15 - j] = -1;
|
||||
x[IX(15 - i,15,15 - j)] = 1;
|
||||
v[IX(15 - i,15,15 - j)] = -1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// x[16][16][16] = 1;
|
||||
|
||||
vel_step(DIM, DIM, DIM, u, v, w, u0, v0, w0, VISCOSITY_CONSTANT, timestep);
|
||||
|
||||
@ -124,73 +127,74 @@ public class FluidSim {
|
||||
}
|
||||
|
||||
|
||||
void lin_solve ( int M, int N, int O, int b, float[][][] x, float[][][] x0, float a, float c ){
|
||||
void lin_solve(int M, int N, int O, int b, float[] x, float[] x0, float a, float c){
|
||||
int i, j, k, l;
|
||||
|
||||
// iterate the solver
|
||||
for ( l=0 ; l<LINEARSOLVERTIMES ; l++ ) {
|
||||
// update for each cell
|
||||
for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
|
||||
x[i][j][k] = (x0[i][j][k] + a*(x[i-1][j][k]+x[i+1][j][k]+x[i][j-1][k]+x[i][j+1][k]+x[i][j][k-1]+x[i][j][k+1]))/c;
|
||||
x[IX(i,j,k)] = (x0[IX(i,j,k)] + a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]+x[IX(i,j,k-1)]+x[IX(i,j,k+1)]))/c;
|
||||
}}}
|
||||
setBounds(x, b);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void diffuse(int M, int N, int O, int b, float[][][] x, float[][][] x0, float diff, float dt){
|
||||
void diffuse(int M, int N, int O, int b, float[] x, float[] x0, float diff, float dt){
|
||||
int max = Math.max(Math.max(M, N), Math.max(N, O));
|
||||
float a=dt*diff*max*max*max;
|
||||
lin_solve ( M, N, O, b, x, x0, a, 1+6*a );
|
||||
linSolveJacobian(M, N, O, b, x, x0, a, 1+6*a);
|
||||
}
|
||||
|
||||
void setBounds(float[][][] target, int b){
|
||||
void setBounds(float[] target, int b){
|
||||
for(int x=1; x < DIM-1; x++){
|
||||
for(int y = 1; y < DIM-1; y++){
|
||||
target[0][x][y] = b==1 ? -target[1][x][y] : target[1][x][y];
|
||||
target[DIM-1][x][y] = b==1 ? -target[DIM-2][x][y] : target[DIM-2][x][y];
|
||||
target[x][0][y] = b==2 ? -target[x][1][y] : target[x][1][y];
|
||||
target[x][DIM-1][y] = b==2 ? -target[x][DIM-2][y] : target[x][DIM-2][y];
|
||||
target[x][y][0] = b==3 ? -target[x][y][1] : target[x][y][1];
|
||||
target[x][y][DIM-1] = b==3 ? -target[x][y][DIM-2] : target[x][y][DIM-2];
|
||||
//((x)+(DIM)*(y) + (DIM)*(DIM)*(z))
|
||||
target[0 + DIM * x + DIM * DIM * y] = b==1 ? -target[1 + DIM * x + DIM * DIM * y] : target[1 + DIM * x + DIM * DIM * y];
|
||||
target[IX(DIM-1,x,y)] = b==1 ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)];
|
||||
target[IX(x,0,y)] = b==2 ? -target[IX(x,1,y)] : target[IX(x,1,y)];
|
||||
target[IX(x,DIM-1,y)] = b==2 ? -target[IX(x,DIM-2,y)] : target[IX(x,DIM-2,y)];
|
||||
target[IX(x,y,0)] = b==3 ? -target[IX(x,y,1)] : target[IX(x,y,1)];
|
||||
target[IX(x,y,DIM-1)] = b==3 ? -target[IX(x,y,DIM-2)] : target[IX(x,y,DIM-2)];
|
||||
}
|
||||
}
|
||||
for(int x = 1; x < DIM-1; x++){
|
||||
target[x][0][0] = (float)(0.5f * (target[x][1][0] + target[x][0][1]));
|
||||
target[x][DIM-1][0] = (float)(0.5f * (target[x][DIM-2][0] + target[x][DIM-1][1]));
|
||||
target[x][0][DIM-1] = (float)(0.5f * (target[x][1][DIM-1] + target[x][0][DIM-2]));
|
||||
target[x][DIM-1][DIM-1] = (float)(0.5f * (target[x][DIM-2][DIM-1] + target[x][DIM-1][DIM-2]));
|
||||
target[IX(x,0,0)] = (float)(0.5f * (target[IX(x,1,0)] + target[IX(x,0,1)]));
|
||||
target[IX(x,DIM-1,0)] = (float)(0.5f * (target[IX(x,DIM-2,0)] + target[IX(x,DIM-1,1)]));
|
||||
target[IX(x,0,DIM-1)] = (float)(0.5f * (target[IX(x,1,DIM-1)] + target[IX(x,0,DIM-2)]));
|
||||
target[IX(x,DIM-1,DIM-1)] = (float)(0.5f * (target[IX(x,DIM-2,DIM-1)] + target[IX(x,DIM-1,DIM-2)]));
|
||||
|
||||
target[0][x][0] = (float)(0.5f * (target[1][x][0] + target[0][x][1]));
|
||||
target[DIM-1][x][0] = (float)(0.5f * (target[DIM-2][x][0] + target[DIM-1][x][1]));
|
||||
target[0][x][DIM-1] = (float)(0.5f * (target[1][x][DIM-1] + target[0][x][DIM-2]));
|
||||
target[DIM-1][x][DIM-1] = (float)(0.5f * (target[DIM-2][x][DIM-1] + target[DIM-1][x][DIM-2]));
|
||||
target[IX(0,x,0)] = (float)(0.5f * (target[IX(1,x,0)] + target[IX(0,x,1)]));
|
||||
target[IX(DIM-1,x,0)] = (float)(0.5f * (target[IX(DIM-2,x,0)] + target[IX(DIM-1,x,1)]));
|
||||
target[IX(0,x,DIM-1)] = (float)(0.5f * (target[IX(1,x,DIM-1)] + target[IX(0,x,DIM-2)]));
|
||||
target[IX(DIM-1,x,DIM-1)] = (float)(0.5f * (target[IX(DIM-2,x,DIM-1)] + target[IX(DIM-1,x,DIM-2)]));
|
||||
|
||||
|
||||
target[0][0][x] = (float)(0.5f * (target[1][0][x] + target[0][1][x]));
|
||||
target[DIM-1][0][x] = (float)(0.5f * (target[DIM-2][0][x] + target[DIM-1][1][x]));
|
||||
target[0][DIM-1][x] = (float)(0.5f * (target[1][DIM-1][x] + target[0][DIM-2][x]));
|
||||
target[DIM-1][DIM-1][x] = (float)(0.5f * (target[DIM-2][DIM-1][x] + target[DIM-1][DIM-2][x]));
|
||||
target[IX(0,0,x)] = (float)(0.5f * (target[IX(1,0,x)] + target[IX(0,1,x)]));
|
||||
target[IX(DIM-1,0,x)] = (float)(0.5f * (target[IX(DIM-2,0,x)] + target[IX(DIM-1,1,x)]));
|
||||
target[IX(0,DIM-1,x)] = (float)(0.5f * (target[IX(1,DIM-1,x)] + target[IX(0,DIM-2,x)]));
|
||||
target[IX(DIM-1,DIM-1,x)] = (float)(0.5f * (target[IX(DIM-2,DIM-1,x)] + target[IX(DIM-1,DIM-2,x)]));
|
||||
|
||||
}
|
||||
target[0][0][0] = (float)((target[1][0][0]+target[0][1][0]+target[0][0][1])/3.0);
|
||||
target[DIM-1][0][0] = (float)((target[DIM-2][0][0]+target[DIM-1][1][0]+target[DIM-1][0][1])/3.0);
|
||||
target[0][DIM-1][0] = (float)((target[1][DIM-1][0]+target[0][DIM-2][0]+target[0][DIM-1][1])/3.0);
|
||||
target[0][0][DIM-1] = (float)((target[0][0][DIM-2]+target[1][0][DIM-1]+target[0][1][DIM-1])/3.0);
|
||||
target[DIM-1][DIM-1][0] = (float)((target[DIM-2][DIM-1][0]+target[DIM-1][DIM-2][0]+target[DIM-1][DIM-1][1])/3.0);
|
||||
target[0][DIM-1][DIM-1] = (float)((target[1][DIM-1][DIM-1]+target[0][DIM-2][DIM-1]+target[0][DIM-1][DIM-2])/3.0);
|
||||
target[DIM-1][0][DIM-1] = (float)((target[DIM-1][0][DIM-2]+target[DIM-2][0][DIM-1]+target[DIM-1][1][DIM-1])/3.0);
|
||||
target[DIM-1][DIM-1][DIM-1] = (float)((target[DIM-1][DIM-1][DIM-2]+target[DIM-1][DIM-2][DIM-1]+target[DIM-1][DIM-1][DIM-2])/3.0);
|
||||
target[IX(0,0,0)] = (float)((target[IX(1,0,0)]+target[IX(0,1,0)]+target[IX(0,0,1)])/3.0);
|
||||
target[IX(DIM-1,0,0)] = (float)((target[IX(DIM-2,0,0)]+target[IX(DIM-1,1,0)]+target[IX(DIM-1,0,1)])/3.0);
|
||||
target[IX(0,DIM-1,0)] = (float)((target[IX(1,DIM-1,0)]+target[IX(0,DIM-2,0)]+target[IX(0,DIM-1,1)])/3.0);
|
||||
target[IX(0,0,DIM-1)] = (float)((target[IX(0,0,DIM-2)]+target[IX(1,0,DIM-1)]+target[IX(0,1,DIM-1)])/3.0);
|
||||
target[IX(DIM-1,DIM-1,0)] = (float)((target[IX(DIM-2,DIM-1,0)]+target[IX(DIM-1,DIM-2,0)]+target[IX(DIM-1,DIM-1,1)])/3.0);
|
||||
target[IX(0,DIM-1,DIM-1)] = (float)((target[IX(1,DIM-1,DIM-1)]+target[IX(0,DIM-2,DIM-1)]+target[IX(0,DIM-1,DIM-2)])/3.0);
|
||||
target[IX(DIM-1,0,DIM-1)] = (float)((target[IX(DIM-1,0,DIM-2)]+target[IX(DIM-2,0,DIM-1)]+target[IX(DIM-1,1,DIM-1)])/3.0);
|
||||
target[IX(DIM-1,DIM-1,DIM-1)] = (float)((target[IX(DIM-1,DIM-1,DIM-2)]+target[IX(DIM-1,DIM-2,DIM-1)]+target[IX(DIM-1,DIM-1,DIM-2)])/3.0);
|
||||
}
|
||||
|
||||
void advect(int M, int N, int O, int b, float[][][] d, float[][][] d0, float[][][] u, float[][][] v, float[][][] w, float dt){
|
||||
void advect(int M, int N, int O, int b, float[] d, float[] d0, float[] u, float[] v, float[] w, float dt){
|
||||
int i, j, k, i0, j0, k0, i1, j1, k1;
|
||||
float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz;
|
||||
|
||||
dtx=dty=dtz=dt*Math.max(Math.max(M, N), Math.max(N, O));
|
||||
|
||||
for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
|
||||
x = i-dtx*u[i][j][k]; y = j-dty*v[i][j][k]; z = k-dtz*w[i][j][k];
|
||||
x = i-dtx*u[IX(i,j,k)]; y = j-dty*v[IX(i,j,k)]; z = k-dtz*w[IX(i,j,k)];
|
||||
if (x<0.5f) x=0.5f; if (x>M+0.5f) x=M+0.5f; i0=(int)x; i1=i0+1;
|
||||
if (y<0.5f) y=0.5f; if (y>N+0.5f) y=N+0.5f; j0=(int)y; j1=j0+1;
|
||||
if (z<0.5f) z=0.5f; if (z>O+0.5f) z=O+0.5f; k0=(int)z; k1=k0+1;
|
||||
@ -214,43 +218,39 @@ public class FluidSim {
|
||||
if(k1 >= DIM){
|
||||
k1 = DIM - 1;
|
||||
}
|
||||
d[i][j][k] = s0*(t0*u0*d0[i0][j0][k0]+t1*u0*d0[i0][j1][k0]+t0*u1*d0[i0][j0][k1]+t1*u1*d0[i0][j1][k1])+
|
||||
s1*(t0*u0*d0[i1][j0][k0]+t1*u0*d0[i1][j1][k0]+t0*u1*d0[i1][j0][k1]+t1*u1*d0[i1][j1][k1]);
|
||||
d[IX(i,j,k)] = s0*(t0*u0*d0[IX(i0,j0,k0)]+t1*u0*d0[IX(i0,j1,k0)]+t0*u1*d0[IX(i0,j0,k1)]+t1*u1*d0[IX(i0,j1,k1)])+
|
||||
s1*(t0*u0*d0[IX(i1,j0,k0)]+t1*u0*d0[IX(i1,j1,k0)]+t0*u1*d0[IX(i1,j0,k1)]+t1*u1*d0[IX(i1,j1,k1)]);
|
||||
}}}
|
||||
|
||||
setBounds(d, b);
|
||||
}
|
||||
|
||||
void moveData(float[][][] set1, float[][][] set2){
|
||||
void moveData(float[] set1, float[] set2){
|
||||
float tmp = 0;
|
||||
for(int x = 0; x < DIM; x++){
|
||||
for(int y = 0; y < DIM; y++){
|
||||
for(int z = 0; z < DIM; z++){
|
||||
tmp = set1[x][y][z];
|
||||
set1[x][y][z] = set2[x][y][z];
|
||||
set2[x][y][z] = tmp;
|
||||
}
|
||||
}
|
||||
for(int i = 0; i < set1.length; i++){
|
||||
tmp = set1[i];
|
||||
set1[i] = set2[i];
|
||||
set2[i] = tmp;
|
||||
}
|
||||
}
|
||||
|
||||
void project(int M, int N, int O, float[][][] u, float[][][] v, float[][][] w, float[][][] p, float[][][] div){
|
||||
void project(int M, int N, int O, float[] u, float[] v, float[] w, float[] p, float[] div){
|
||||
int i, j, k;
|
||||
|
||||
for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
|
||||
div[i][j][k] = (float)(-1.0/3.0*((u[i+1][j][k]-u[i-1][j][k])/M+(v[i][j+1][k]-v[i][j-1][k])/M+(w[i][j][k+1]-w[i][j][k-1])/M));
|
||||
p[i][j][k] = 0;
|
||||
div[IX(i,j,k)] = (float)(-1.0/3.0*((u[IX(i+1,j,k)]-u[IX(i-1,j,k)])/M+(v[IX(i,j+1,k)]-v[IX(i,j-1,k)])/M+(w[IX(i,j,k+1)]-w[IX(i,j,k-1)])/M));
|
||||
p[IX(i,j,k)] = 0;
|
||||
}}}
|
||||
|
||||
setBounds(div, 0);
|
||||
setBounds(p, 0);
|
||||
|
||||
lin_solve ( M, N, O, 0, p, div, 1, 6 );
|
||||
linSolveJacobian(M, N, O, 0, p, div, 1, 6);
|
||||
|
||||
for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
|
||||
u[i][j][k] -= 0.5f*M*(p[i+1][j][k]-p[i-1][j][k]);
|
||||
v[i][j][k] -= 0.5f*M*(p[i][j+1][k]-p[i][j-1][k]);
|
||||
w[i][j][k] -= 0.5f*M*(p[i][j][k+1]-p[i][j][k-1]);
|
||||
u[IX(i,j,k)] -= 0.5f*M*(p[IX(i+1,j,k)]-p[IX(i-1,j,k)]);
|
||||
v[IX(i,j,k)] -= 0.5f*M*(p[IX(i,j+1,k)]-p[IX(i,j-1,k)]);
|
||||
w[IX(i,j,k)] -= 0.5f*M*(p[IX(i,j,k+1)]-p[IX(i,j,k-1)]);
|
||||
}}}
|
||||
|
||||
setBounds(u, 1);
|
||||
@ -258,7 +258,7 @@ public class FluidSim {
|
||||
setBounds(w, 3);
|
||||
}
|
||||
|
||||
void dens_step(int M, int N, int O, float[][][] x, float[][][] x0, float[][][] u, float[][][] v, float[][][] w, float diff, float dt){
|
||||
void dens_step(int M, int N, int O, float[] x, float[] x0, float[] u, float[] v, float[] w, float diff, float dt){
|
||||
// add_source ( M, N, O, x, x0, dt );
|
||||
moveData ( x0, x ); diffuse ( M, N, O, 0, x, x0, diff, dt );
|
||||
moveData ( x0, x ); advect ( M, N, O, 0, x, x0, u, v, w, dt );
|
||||
@ -266,25 +266,32 @@ public class FluidSim {
|
||||
|
||||
|
||||
|
||||
void add_source(int M, int N, int O, float[][][] x, float[][][] density, float dt){
|
||||
void add_source(int M, int N, int O, float[] x, float[] density, float dt){
|
||||
int i, j, k;
|
||||
for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
|
||||
if(density[i][j][k] > 0.0f){
|
||||
x[i][j][k] += dt * GRAVITY * density[i][j][k];
|
||||
if(density[IX(i,j,k)] > 0.0f){
|
||||
x[IX(i,j,k)] += dt * GRAVITY * density[IX(i,j,k)];
|
||||
}
|
||||
}}}
|
||||
}
|
||||
|
||||
void vel_step(int M, int N, int O, float[][][] u, float[][][] v, float[][][] w, float[][][] u0, float[][][] v0, float[][][] w0, float visc, float dt ){
|
||||
void vel_step(int M, int N, int O, float[] u, float[] v, float[] w, float[] u0, float[] v0, float[] w0, float visc, float dt ){
|
||||
// add_source ( M, N, O, u, u0, dt ); add_source ( M, N, O, v, v0, dt );add_source ( M, N, O, w, w0, dt );
|
||||
add_source ( M, N, O, v, x, dt );
|
||||
moveData ( u0, u ); diffuse ( M, N, O, 1, u, u0, visc, dt );
|
||||
moveData ( v0, v ); diffuse ( M, N, O, 2, v, v0, visc, dt );
|
||||
moveData ( w0, w ); diffuse ( M, N, O, 3, w, w0, visc, dt );
|
||||
project ( M, N, O, u, v, w, u0, v0 );
|
||||
moveData ( u0, u ); moveData ( v0, v );moveData ( w0, w );
|
||||
advect ( M, N, O, 1, u, u0, u0, v0, w0, dt ); advect ( M, N, O, 2, v, v0, u0, v0, w0, dt );advect ( M, N, O, 3, w, w0, u0, v0, w0, dt );
|
||||
project ( M, N, O, u, v, w, u0, v0 );
|
||||
add_source(M, N, O, v, x, dt);
|
||||
moveData(u0, u);
|
||||
diffuse(M, N, O, 1, u, u0, visc, dt);
|
||||
moveData(v0, v);
|
||||
diffuse(M, N, O, 2, v, v0, visc, dt);
|
||||
moveData(w0, w);
|
||||
diffuse(M, N, O, 3, w, w0, visc, dt);
|
||||
project(M, N, O, u, v, w, u0, v0);
|
||||
moveData(u0, u);
|
||||
moveData(v0, v);
|
||||
moveData(w0, w);
|
||||
advect(M, N, O, 1, u, u0, u0, v0, w0, dt);
|
||||
advect(M, N, O, 2, v, v0, u0, v0, w0, dt);
|
||||
advect(M, N, O, 3, w, w0, u0, v0, w0, dt);
|
||||
project(M, N, O, u, v, w, u0, v0);
|
||||
}
|
||||
|
||||
void sum(){
|
||||
@ -292,18 +299,126 @@ public class FluidSim {
|
||||
for(int i = 0; i < DIM; i++){
|
||||
for(int j = 0; j < DIM; j++){
|
||||
for(int k = 0; k < DIM; k++){
|
||||
sum[0] = sum[0] + x[i][j][k];
|
||||
sum[1] = sum[1] + u[i][j][k];
|
||||
sum[2] = sum[2] + v[i][j][k];
|
||||
sum[3] = sum[3] + w[i][j][k];
|
||||
sum[0] = sum[0] + x[IX(i,j,k)];
|
||||
sum[1] = sum[1] + u[IX(i,j,k)];
|
||||
sum[2] = sum[2] + v[IX(i,j,k)];
|
||||
sum[3] = sum[3] + w[IX(i,j,k)];
|
||||
}
|
||||
}
|
||||
}
|
||||
// System.out.println(sum[0] + " " + sum[1] + " " + sum[2] + " " + sum[3]);
|
||||
}
|
||||
|
||||
public float[][][] getData(){
|
||||
public float[] getData(){
|
||||
return x;
|
||||
}
|
||||
|
||||
public static final int IX(int x, int y, int z){
|
||||
return ((x)+(DIM)*(y) + (DIM)*(DIM)*(z));
|
||||
}
|
||||
|
||||
static final int[] vectorOffsetMap = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)];
|
||||
static final int[] vectorOffsetMapXN = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)];
|
||||
static final int[] vectorOffsetMapXP = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)];
|
||||
static final int[] vectorOffsetMapYN = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)];
|
||||
static final int[] vectorOffsetMapYP = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)];
|
||||
static final int[] vectorOffsetMapZN = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)];
|
||||
static final int[] vectorOffsetMapZP = new int[(DIM - 2) * (DIM - 2) * (DIM - 2)];
|
||||
static {
|
||||
int i = 0;
|
||||
for(int x = 1; x < DIM - 1; x++){
|
||||
for(int y = 1; y < DIM - 1; y++){
|
||||
for(int z = 1; z < DIM - 1; z++){
|
||||
// if(x != y || y != z){
|
||||
vectorOffsetMap[i] = IX(x,y,z);
|
||||
vectorOffsetMapXN[i] = IX(x-1, y, z);
|
||||
vectorOffsetMapXP[i] = IX(x+1, y, z);
|
||||
vectorOffsetMapYN[i] = IX( x,y-1, z);
|
||||
vectorOffsetMapYP[i] = IX( x,y+1, z);
|
||||
vectorOffsetMapZN[i] = IX( x, y,z-1);
|
||||
vectorOffsetMapZP[i] = IX( x, y,z+1);
|
||||
i++;
|
||||
// }
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static final int LANE_WIDTH = 8;
|
||||
void linSolveJacobian(int M, int N, int O, int b, float[] x, float x0[], float a, float c){
|
||||
//x contains the solution to the problem
|
||||
//x0 contains A
|
||||
|
||||
// float[] a = new float[]{
|
||||
// 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16
|
||||
// };
|
||||
// FloatVector vec = FloatVector.fromArray(FloatVector.SPECIES_128, a, 0);
|
||||
// FloatVector finalVec = vec.add(3);
|
||||
// float[] b = finalVec.toArray();
|
||||
// for(int i = 0; i < b.length; i++){
|
||||
// System.out.println(b[i]);
|
||||
// }
|
||||
float[] vectorArrView;
|
||||
float[] srcArr, destArr;
|
||||
for(int k = 0; k < LINEARSOLVERTIMES; k++){
|
||||
if(k / 2 == 0){
|
||||
srcArr = x;
|
||||
destArr = jacobiAltArray;
|
||||
} else {
|
||||
srcArr = jacobiAltArray;
|
||||
destArr = x;
|
||||
}
|
||||
for(int i = 0; i < vectorOffsetMap.length; i += LANE_WIDTH){
|
||||
if(i + LANE_WIDTH < vectorOffsetMap.length){
|
||||
|
||||
// for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
|
||||
// x[IX(i,j,k)] = (x0[IX(i,j,k)] + a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]+x[IX(i,j,k-1)]+x[IX(i,j,k+1)]))/c;
|
||||
// }}}
|
||||
|
||||
// adderVector = FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapXN,i);
|
||||
// x0Vector = FloatVector.fromArray(FloatVector.SPECIES_256,x0,0,vectorOffsetMap,i);
|
||||
vectorArrView = FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapXN,i)
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapXP,i))
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapYN,i))
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapYP,i))
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapZN,i))
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapZP,i))
|
||||
.mul(a)
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x0,0,vectorOffsetMap,i))
|
||||
.div(c).toArray();
|
||||
|
||||
for(int l = 0; l < LANE_WIDTH; l++){
|
||||
x[vectorOffsetMap[i+l]] = vectorArrView[l];
|
||||
}
|
||||
} else {
|
||||
int numLeft = vectorOffsetMap.length - i;
|
||||
|
||||
VectorMask<Float> mask = FloatVector.SPECIES_256.indexInRange(O, numLeft);
|
||||
|
||||
// for ( i=1 ; i<M-1 ; i++ ) { for ( j=1 ; j<N-1 ; j++ ) { for ( k=1 ; k<O-1 ; k++ ) {
|
||||
// x[IX(i,j,k)] = (x0[IX(i,j,k)] + a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]+x[IX(i,j,k-1)]+x[IX(i,j,k+1)]))/c;
|
||||
// }}}
|
||||
|
||||
// adderVector = FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapXN,i,mask);
|
||||
// x0Vector = FloatVector.fromArray(FloatVector.SPECIES_256,x0,0,vectorOffsetMap,i,mask);
|
||||
vectorArrView = FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapXN,i,mask)
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapXP,i))
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapYN,i))
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapYP,i))
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapZN,i))
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x,0,vectorOffsetMapZP,i))
|
||||
.mul(a)
|
||||
.add(FloatVector.fromArray(FloatVector.SPECIES_256,x0,0,vectorOffsetMap,i,mask))
|
||||
.div(c).toArray();
|
||||
|
||||
for(int l = 0; l < numLeft; l++){
|
||||
x[vectorOffsetMap[i+l]] = vectorArrView[l];
|
||||
}
|
||||
}
|
||||
}
|
||||
setBounds(destArr, b);
|
||||
}
|
||||
//x[IX(i,j,k)] = (x0[IX(i,j,k)] + a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]+x[IX(i,j,k-1)]+x[IX(i,j,k+1)]))/c;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -4,7 +4,6 @@ import java.io.File;
|
||||
import java.io.IOException;
|
||||
import java.nio.file.Files;
|
||||
import java.nio.file.Paths;
|
||||
// import jdk.incubator.vector.FloatVector;
|
||||
import java.util.concurrent.TimeUnit;
|
||||
|
||||
import electrosphere.render.GLFWContext;
|
||||
@ -16,30 +15,28 @@ import electrosphere.render.Mesh;
|
||||
public class Main {
|
||||
|
||||
public static void main(String args[]){
|
||||
// float[] a = new float[]{
|
||||
// 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16
|
||||
// };
|
||||
// FloatVector vec = FloatVector.fromArray(FloatVector.SPECIES_128, a, 0);
|
||||
// FloatVector finalVec = vec.add(3);
|
||||
// float[] b = finalVec.toArray();
|
||||
// for(int i = 0; i < b.length; i++){
|
||||
// System.out.println(b[i]);
|
||||
// }
|
||||
GLFWContext.init();
|
||||
FluidSim sim = new FluidSim();
|
||||
sim.setup();
|
||||
Mesh.meshInitially(sim);
|
||||
int i = 0;
|
||||
long time = 0;
|
||||
long lastTime = 0;
|
||||
while(true){
|
||||
try {
|
||||
TimeUnit.MILLISECONDS.sleep(1);
|
||||
} catch (InterruptedException e) {
|
||||
e.printStackTrace();
|
||||
}
|
||||
lastTime = System.currentTimeMillis();
|
||||
sim.simulate(i,0.001f);
|
||||
time = time + (System.currentTimeMillis() - lastTime);
|
||||
Mesh.remesh(sim);
|
||||
GLFWContext.redraw();
|
||||
i++;
|
||||
if(i == 1000){
|
||||
System.out.println(time / 1000.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -38,7 +38,7 @@ public class Mesh {
|
||||
static Matrix4f projectionMatrix = new Matrix4f();
|
||||
static Matrix4f model = new Matrix4f().identity();
|
||||
|
||||
static MemoryStack stack = MemoryStack.create(4 * 1000 * 1000);
|
||||
static MemoryStack stack = MemoryStack.create(16 * 1000 * 1000);
|
||||
|
||||
public static void meshInitially(FluidSim sim){
|
||||
//create and bind vao
|
||||
@ -820,7 +820,7 @@ public class Mesh {
|
||||
return new Vector3f(x,y,z);
|
||||
}
|
||||
|
||||
public static TerrainChunkData generateTerrainChunkData(float[][][] data){
|
||||
public static TerrainChunkData generateTerrainChunkData(float[] data){
|
||||
|
||||
// 5 6
|
||||
// +-------------+ +-----5-------+ ^ Y
|
||||
@ -849,15 +849,15 @@ public class Mesh {
|
||||
|
||||
|
||||
|
||||
for(int x = 0; x < data.length - 1; x++){
|
||||
for(int y = 0; y < data[0].length - 1; y++){
|
||||
for(int z = 0; z < data[0][0].length - 1; z++){
|
||||
for(int x = 0; x < FluidSim.DIM - 1; x++){
|
||||
for(int y = 0; y < FluidSim.DIM - 1; y++){
|
||||
for(int z = 0; z < FluidSim.DIM - 1; z++){
|
||||
//push the current cell's values into the gridcell
|
||||
currentCell.setValues(
|
||||
new Vector3f(x+0,y+0,z+0), new Vector3f(x+0,y+0,z+1), new Vector3f(x+1,y+0,z+1), new Vector3f(x+1,y+0,z+0),
|
||||
new Vector3f(x+0,y+1,z+0), new Vector3f(x+0,y+1,z+1), new Vector3f(x+1,y+1,z+1), new Vector3f(x+1,y+1,z+0),
|
||||
data[x+0][y+0][z+0], data[x+0][y+0][z+1], data[x+1][y+0][z+1], data[x+1][y+0][z+0],
|
||||
data[x+0][y+1][z+0], data[x+0][y+1][z+1], data[x+1][y+1][z+1], data[x+1][y+1][z+0]
|
||||
data[FluidSim.IX(x+0,y+0,z+0)], data[FluidSim.IX(x+0,y+0,z+1)], data[FluidSim.IX(x+1,y+0,z+1)], data[FluidSim.IX(x+1,y+0,z+0)],
|
||||
data[FluidSim.IX(x+0,y+1,z+0)], data[FluidSim.IX(x+0,y+1,z+1)], data[FluidSim.IX(x+1,y+1,z+1)], data[FluidSim.IX(x+1,y+1,z+0)]
|
||||
);
|
||||
//polygonize the current gridcell
|
||||
polygonize(currentCell, 0.01, triangles, vertMap, verts, normals, trianglesSharingVert);
|
||||
|
||||
@ -1,76 +0,0 @@
|
||||
// #ifdef DOUBLE_FP
|
||||
// #ifdef AMD_FP
|
||||
// #pragma OPENCL EXTENSION cl_amd_fp64 : enable
|
||||
// #else
|
||||
// #ifndef CL_VERSION_1_2
|
||||
// #pragma OPENCL EXTENSION cl_khr_fp64 : enable
|
||||
// #endif
|
||||
// #endif
|
||||
// #define varfloat double
|
||||
// #define _255 255.0
|
||||
// #else
|
||||
// #define varfloat float
|
||||
// #define _255 255.0f
|
||||
// #endif
|
||||
|
||||
// #ifdef USE_TEXTURE
|
||||
// #define OUTPUT_TYPE __write_only image2d_t
|
||||
// #else
|
||||
// #define OUTPUT_TYPE global uint *
|
||||
// #endif
|
||||
|
||||
|
||||
__kernel void test(
|
||||
__read_only image2d_t input,
|
||||
__write_only image2d_t output
|
||||
) {
|
||||
unsigned int ix = get_global_id(0);
|
||||
unsigned int iy = get_global_id(1);
|
||||
|
||||
|
||||
write_imagef(output, (int2)(ix, iy), (float4)127 - read_imagef(input, (int2)(ix, iy)));
|
||||
|
||||
// varfloat r = x0 + ix * rangeX / width;
|
||||
// varfloat i = y0 + iy * rangeY / height;
|
||||
|
||||
// varfloat x = 0;
|
||||
// varfloat y = 0;
|
||||
|
||||
// varfloat magnitudeSquared = 0;
|
||||
// int iteration = 0;
|
||||
|
||||
// while ( magnitudeSquared < 4 && iteration < maxIterations ) {
|
||||
// varfloat x2 = x*x;
|
||||
// varfloat y2 = y*y;
|
||||
// y = 2 * x * y + i;
|
||||
// x = x2 - y2 + r;
|
||||
// magnitudeSquared = x2+y2;
|
||||
// iteration++;
|
||||
// }
|
||||
|
||||
// if ( iteration == maxIterations ) {
|
||||
// #ifdef USE_TEXTURE
|
||||
// write_imagef(output, (int2)(ix, iy), (float4)0);
|
||||
// #else
|
||||
// output[iy * width + ix] = 0;
|
||||
// #endif
|
||||
// } else {
|
||||
// float alpha = (float)iteration / maxIterations;
|
||||
// int colorIndex = (int)(alpha * colorMapSize);
|
||||
// #ifdef USE_TEXTURE
|
||||
// // We could have changed colorMap to a texture + sampler, but the
|
||||
// // unpacking below has minimal overhead and it's kinda interesting.
|
||||
// uint c = colorMap[colorIndex];
|
||||
// write_imageui(output, (int2)(ix, iy), (uint4)(
|
||||
// (c >> 0) & 0xFF,
|
||||
// (c >> 8) & 0xFF,
|
||||
// (c >> 16) & 0xFF,
|
||||
// 0xFF
|
||||
// ));
|
||||
// #else
|
||||
// output[iy * width + ix] = colorMap[colorIndex];
|
||||
// #endif
|
||||
// // monochrom
|
||||
// //output[iy * width + ix] = 255*iteration/maxIterations;
|
||||
// }
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user