From 072d4f7784d7b1faeb6d618e66625239cd1bbfc8 Mon Sep 17 00:00:00 2001
From: unknown <>
Date: Thu, 6 Jul 2023 20:45:37 -0400
Subject: [PATCH] SIMD support

---
 pom.xml                                      |   2 +-
 src/main/java/electrosphere/FluidSim.java    | 289 +++++++++++++------
 src/main/java/electrosphere/Main.java        |  17 +-
 src/main/java/electrosphere/render/Mesh.java |  14 +-
 src/main/resources/Test1.cl                  |  76 -----
 5 files changed, 217 insertions(+), 181 deletions(-)
 delete mode 100644 src/main/resources/Test1.cl

diff --git a/pom.xml b/pom.xml
index 2b5933d..b1e600d 100644
--- a/pom.xml
+++ b/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>
diff --git a/src/main/java/electrosphere/FluidSim.java b/src/main/java/electrosphere/FluidSim.java
index 385af5a..3eddf66 100644
--- a/src/main/java/electrosphere/FluidSim.java
+++ b/src/main/java/electrosphere/FluidSim.java
@@ -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;
+    }
+
 }
diff --git a/src/main/java/electrosphere/Main.java b/src/main/java/electrosphere/Main.java
index c3bd4a6..2fdaf01 100644
--- a/src/main/java/electrosphere/Main.java
+++ b/src/main/java/electrosphere/Main.java
@@ -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);
+            }
         }
     }
 
diff --git a/src/main/java/electrosphere/render/Mesh.java b/src/main/java/electrosphere/render/Mesh.java
index 933346e..1e6fe62 100644
--- a/src/main/java/electrosphere/render/Mesh.java
+++ b/src/main/java/electrosphere/render/Mesh.java
@@ -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);
diff --git a/src/main/resources/Test1.cl b/src/main/resources/Test1.cl
deleted file mode 100644
index af10ae7..0000000
--- a/src/main/resources/Test1.cl
+++ /dev/null
@@ -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;
-    // }
-}
\ No newline at end of file