updates

Reviewed-on: https://git.austinwhoover.com/studiorailgun/fluid-sim/pulls/1

.gitignore

@@ -5,6 +5,6 @@
 /.settings
 /.classpath
 /.project
-/.vscode
 /shared-folder
-/shared-folder/**
+/shared-folder/**
+/chunks

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "src/main/c/lib/stb"]
+	path = src/main/c/lib/stb
+	url = https://github.com/nothings/stb.git

.vscode/c_cpp_properties.json

@@ -0,0 +1,22 @@
+{
+    "configurations": [
+        {
+            "name": "Win32",
+            "includePath": [
+                "${workspaceFolder}/**",
+                "C:/Program Files/Eclipse Adoptium/jdk-21.0.5.11-hotspot/include/**"
+            ],
+            "defines": [
+                "_DEBUG",
+                "UNICODE",
+                "_UNICODE"
+            ],
+            "windowsSdkVersion": "10.0.19041.0",
+            "compilerPath": "C:/ProgramData/mingw64/mingw64/bin/gcc.exe",
+            "cStandard": "c17",
+            "cppStandard": "c++17",
+            "intelliSenseMode": "windows-gcc-x64"
+        }
+    ],
+    "version": 4
+}

.vscode/extensions.json

@@ -0,0 +1,6 @@
+{
+    "recommendations": [
+        "ms-vscode.cpptools-extension-pack",
+        "vscjava.vscode-java-pack"
+    ]
+}

.vscode/launch.json

@@ -0,0 +1,20 @@
+{
+    "configurations": [
+        {
+            "type": "java",
+            "name": "Launch Java Program",
+            "request": "launch",
+            "mainClass": "electrosphere.Main",
+            "preLaunchTask": "compileNoSteps",
+            "vmArgs": "-Djava.library.path=./shared-folder"
+        },
+        {
+            "type": "java",
+            "name": "Launch Java Program (SAVE CHUNKS TO DISK)",
+            "request": "launch",
+            "mainClass": "electrosphere.Main",
+            "preLaunchTask": "compileWithSteps",
+            "vmArgs": "-Djava.library.path=./shared-folder"
+        }
+    ]
+}

.vscode/settings.json

@@ -0,0 +1,21 @@
+{
+    "java.configuration.updateBuildConfiguration": "automatic",
+    "files.associations": {
+        "stdio.h": "c",
+        "jni.h": "c",
+        "immintrin.h": "c",
+        "stdint.h": "c",
+        "utilities.h": "c",
+        "chunkmask.h": "c",
+        "pthread.h": "c",
+        "semaphore.h": "c",
+        "libfluidsim.h": "c",
+        "stb_ds.h": "c",
+        "threadpool.h": "c",
+        "electrosphere_fluidsim.h": "c",
+        "type_traits": "c",
+        "mainfunctions.h": "c",
+        "chunk.h": "c",
+        "simulation.h": "c"
+    }
+}

.vscode/tasks.json

@@ -0,0 +1,15 @@
+{
+    "version": "2.0.0",
+    "tasks": [
+        {
+            "label": "compileWithSteps",
+            "type": "shell",
+            "command": "SAVE_STEPS=1 ./src/main/c/compile.sh",
+        },
+        {
+            "label": "compileNoSteps",
+            "type": "shell",
+            "command": "SAVE_STEPS=0 ./src/main/c/compile.sh",
+        }
+    ]
+}

Jenkinsfile

@@ -0,0 +1,17 @@
+pipeline {
+    agent any
+    tools {
+        maven '3.9.6'
+    }
+    stages {
+        stage('Build') {
+            steps {
+                sh 'mvn --version'
+                sh 'java -version'
+                sh 'git submodule init'
+                sh 'git submodule update'
+                sh 'mvn clean generate-resources package'
+            }
+        }
+    }
+}

README.md

@@ -0,0 +1,11 @@
+# Fluid Sim
+
+
+
+## How to run
+
+Clone with `git clone --recurse-submodules git@git.austinwhoover.com:studiorailgun/fluid-sim.git`
+
+Then run with VSCode debug plugin
+
+

compare.sh

@@ -0,0 +1,9 @@
+search_dir=./chunks
+for entry in "$search_dir"/*
+do
+    if ! diff "/c/Users/satellite/Documents/fluid-sim/$entry" "/c/Users/satellite/Documents/fluid-sim-the-golden-code/$entry" > /dev/null;
+    then
+        echo "“$entry differs"
+    fi
+done
+echo "done!"

dependency-reduced-pom.xml

@@ -73,6 +73,9 @@
             </goals>
             <configuration>
               <executable>bash</executable>
+              <environmentVariables>
+                <SAVE_STEPS>0</SAVE_STEPS>
+              </environmentVariables>
               <arguments>
                 <argument>./src/main/c/compile.sh</argument>
               </arguments>
@@ -80,6 +83,10 @@
           </execution>
         </executions>
       </plugin>
+      <plugin>
+        <artifactId>maven-surefire-plugin</artifactId>
+        <version>3.2.5</version>
+      </plugin>
     </plugins>
   </build>
   <profiles>
@@ -103,6 +110,56 @@
       </properties>
     </profile>
   </profiles>
+  <dependencies>
+    <dependency>
+      <groupId>org.junit.jupiter</groupId>
+      <artifactId>junit-jupiter-engine</artifactId>
+      <version>5.10.2</version>
+      <scope>test</scope>
+      <exclusions>
+        <exclusion>
+          <artifactId>junit-platform-engine</artifactId>
+          <groupId>org.junit.platform</groupId>
+        </exclusion>
+        <exclusion>
+          <artifactId>junit-jupiter-api</artifactId>
+          <groupId>org.junit.jupiter</groupId>
+        </exclusion>
+        <exclusion>
+          <artifactId>apiguardian-api</artifactId>
+          <groupId>org.apiguardian</groupId>
+        </exclusion>
+      </exclusions>
+    </dependency>
+    <dependency>
+      <groupId>org.junit.platform</groupId>
+      <artifactId>junit-platform-runner</artifactId>
+      <version>1.10.2</version>
+      <scope>test</scope>
+      <exclusions>
+        <exclusion>
+          <artifactId>junit</artifactId>
+          <groupId>junit</groupId>
+        </exclusion>
+        <exclusion>
+          <artifactId>junit-platform-launcher</artifactId>
+          <groupId>org.junit.platform</groupId>
+        </exclusion>
+        <exclusion>
+          <artifactId>junit-platform-suite-api</artifactId>
+          <groupId>org.junit.platform</groupId>
+        </exclusion>
+        <exclusion>
+          <artifactId>junit-platform-suite-commons</artifactId>
+          <groupId>org.junit.platform</groupId>
+        </exclusion>
+        <exclusion>
+          <artifactId>apiguardian-api</artifactId>
+          <groupId>org.apiguardian</groupId>
+        </exclusion>
+      </exclusions>
+    </dependency>
+  </dependencies>
   <dependencyManagement>
     <dependencies>
       <dependency>

optimizationideas.md

@@ -0,0 +1,26 @@
+'sleep' chunks if there are no massive velocity swings for ~100 frames
+ - unsleep as soon as new density or velocity is added
+
+ VVV These two may not work because the projection phase currently iterates over all cells with avx instructions
+      it also may generally not work as fast moving water basically works by having empty cells query area that already has water (so maybe no sleeping air)
+'sleep' individual cells that are full of water and fully surrounded by water
+'sleep' terrain cells 
+'sleep' "inactive" cells
+ - When reading data back into java, while iterating over each cell check if they've changed value
+ - If not, add 1 to an accumulating array alongside the other ones (d, u, etc)
+ - In C, if the accumulating array is >100 or some threshold, skip the cell
+ - If a force gets added java side that is greater than some threshold, set all accumulating array to 0
+
+multigrid method for projection code
+
+intelligent neighbor awakening --
+If a neighbor is asleep (say a neighbor full of air above this cell),
+do not simulate that neighbor at all
+Instead poll the currently awake cell for large velocity spikes in the direction of the neighbor
+If there is a large spike, awake the neighbor cell
+if a neighbor is awoken, the simulator should go back through and simulate for that chunk as well (idk if this is possible)
+
+
+use avx bitmasking for handling terrain in projection code ???
+
+

pom.xml

@@ -64,6 +64,7 @@
 
 
     <dependencies>
+        <!--lwjgl-->
         <dependency>
             <groupId>org.lwjgl</groupId>
             <artifactId>lwjgl</artifactId>
@@ -109,6 +110,19 @@
             <artifactId>lwjgl-opengl</artifactId>
             <classifier>${lwjgl.natives}</classifier>
         </dependency>
+        <!--junit-->
+        <dependency>
+            <groupId>org.junit.jupiter</groupId>
+            <artifactId>junit-jupiter-engine</artifactId>
+            <version>5.10.2</version>
+            <scope>test</scope>
+        </dependency>
+        <dependency>
+            <groupId>org.junit.platform</groupId>
+            <artifactId>junit-platform-runner</artifactId>
+            <version>1.10.2</version>
+            <scope>test</scope>
+        </dependency>
     </dependencies>
     
     <build>
@@ -181,6 +195,9 @@
                         </goals>
                         <configuration>
                             <executable>bash</executable>
+                            <environmentVariables>
+                                <SAVE_STEPS>0</SAVE_STEPS>
+                            </environmentVariables>
                             <arguments>
                                 <argument>./src/main/c/compile.sh</argument>
                             </arguments>
@@ -188,6 +205,11 @@
                     </execution>
                 </executions>
             </plugin>
+            <plugin>
+                <groupId>org.apache.maven.plugins</groupId>
+                <artifactId>maven-surefire-plugin</artifactId>
+                <version>3.2.5</version>
+            </plugin>
         </plugins>
     </build>
 

src/main/c/compile.sh

@@ -1,5 +1,7 @@
 cd ./src/main/c
 
+echo "SAVE STEPS $SAVE_STEPS"
+
 LIB_ENDING=".so"
 BASE_INCLUDE_DIR=""
 OS_INCLUDE_DIR=""
@@ -41,8 +43,28 @@ rm -f ./*.dll
 
 
 #compile object files
-COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
-INPUT_FILES="./fluidsim.c"
+# COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -march=native -Ofast -msse -msse2 -msse3 -mmmx -m3dnow"
+# INPUT_FILES="./src/densitystep.c"
+# OUTPUT_FILE="./densitystep.o"
+# gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
+
+# COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -march=native -Ofast -msse -msse2 -msse3 -mmmx -m3dnow"
+# INPUT_FILES="./src/velocitystep.c"
+# OUTPUT_FILE="./velocitystep.o"
+# gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
+
+# COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -march=native -Ofast -msse -msse2 -msse3 -mmmx -m3dnow"
+# INPUT_FILES="./src/chunkmask.c"
+# OUTPUT_FILE="./chunkmask.o"
+# gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
+
+COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -march=native -Ofast -DSAVE_STEPS=$SAVE_STEPS"
+INPUT_FILES="./src/javainterface.c"
+OUTPUT_FILE="./javainterface.o"
+gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
+
+COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -march=native -Ofast -DSAVE_STEPS=$SAVE_STEPS"
+INPUT_FILES="./src/fluidsim.c"
 OUTPUT_FILE="./fluidsim.o"
 gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
 
@@ -51,7 +73,7 @@ gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OU
 #compile shared object file
 OUTPUT_FILE="libfluidsim$LIB_ENDING"
 COMPILE_FLAGS="-shared"
-INPUT_FILES="fluidsim.o"
+INPUT_FILES="fluidsim.o javainterface.o"
 gcc $COMPILE_FLAGS $INPUT_FILES -o $OUTPUT_FILE
 
 #move to resources

src/main/c/fluidsim.c

@@ -1,382 +0,0 @@
-#include <jni.h>
-#include <stdio.h>
-#include <immintrin.h>
-
-
-#define SWAP(x0,x) {float *tmp=x0;x0=x;x=tmp;}
-#define IX(i,j,k) ((i)+(N)*(j)+(N*N)*(k))
-
-#define LINEARSOLVERTIMES 10
-
-
-
-void diffuse(int N, int b, float * x, float * x0, float diff, float dt);
-void advect(int N, int b, float * d, float * d0, float * u, float * v, float * w, float dt);
-void project(int N, float * u, float * v, float * w, float * p, float * div);
-void set_bnd(int N, int b, float * x);
-void dens_step(int N, float * x, float * x0, float * u, float * v, float * w, float diff, float dt);
-void vel_step(int N, float * u, float * v, float * w, float * u0, float * v0, float * w0, float visc, float dt);
-void lin_solve(int N, int b, float* x, float* x0, float a, float c);
-
-/**
- * The core simulation function
-*/
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
-    JNIEnv * env,
-    jobject this,
-    jint DIM_X,
-    jint DIM_Y,
-    jint DIM_Z,
-    jobject jx,
-    jobject jx0,
-    jobject ju,
-    jobject jv,
-    jobject jw,
-    jobject ju0,
-    jobject jv0,
-    jobject jw0,
-    jfloat DIFFUSION_RATE,
-    jfloat VISCOSITY_RATE,
-    jfloat timestep){
-    jboolean isCopy;
-    float * x = (*env)->GetDirectBufferAddress(env,jx);
-    float * x0 = (*env)->GetDirectBufferAddress(env,jx0);
-    float * u = (*env)->GetDirectBufferAddress(env,ju);
-    float * v = (*env)->GetDirectBufferAddress(env,jv);
-    float * w = (*env)->GetDirectBufferAddress(env,jw);
-    float * u0 = (*env)->GetDirectBufferAddress(env,ju0);
-    float * v0 = (*env)->GetDirectBufferAddress(env,jv0);
-    float * w0 = (*env)->GetDirectBufferAddress(env,jw0);
-    int N = DIM_X;
-    int i,j,k;
-    vel_step(DIM_X, u, v, w, u0, v0, w0, VISCOSITY_RATE, timestep);
-    dens_step(DIM_X, x, x0, u, v, w, DIFFUSION_RATE, timestep);
-}
-
-/**
- * Adds values from a source array to a current frame array (eg more density to the main density array)
-*/
-void add_source(int N, float * x, float * s, float dt){
-	int i;
-    int size=N*N*N;
-	for(i=0; i<size; i++){
-        x[i] += dt*s[i];
-    }
-}
-
-/**
- * Diffuses a given array by a diffusion constant
-*/
-void diffuse(int N, int b, float * x, float * x0, float diff, float dt){
-    float a=dt*diff*N*N*N;
-    lin_solve(N, b, x, x0, a, 1+6*a);
-}
-
-/**
- * Advects a given array based on the force vectors in the simulation
-*/
-void advect(int N, 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*N;
-
-    for(k=1; k<N-1; k++){
-        for(j=1; j<N-1; j++){
-            for(i=1; i<N-1; i++){
-                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>N+0.5f) x=N+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>N+0.5f) z=N+0.5f; k0=(int)z; k1=k0+1;
-
-                s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; u1 = z-k0; u0 = 1-u1;
-                if(i0 >= N){
-                    i0 = N - 1;
-                }
-                // if(i0 < 0){
-                //     i0 = 0;
-                // }
-                if(j0 >= N){
-                    j0 = N - 1;
-                }
-                // if(j0 < 0){
-                //     j0 = 0;
-                // }
-                if(k0 >= N){
-                    k0 = N - 1;
-                }
-                // if(k0 < 0){
-                //     k0 = 0;
-                // }
-                if(i1 >= N){
-                    i1 = N - 1;
-                }
-                // if(i1 < 0){
-                //     i1 = 0;
-                // }
-                if(j1 >= N){
-                    j1 = N - 1;
-                }
-                // if(j1 < 0){
-                //     j1 = 0;
-                // }
-                if(k1 >= N){
-                    k1 = N - 1;
-                }
-                // if(k1 < 0){
-                //     k1 = 0;
-                // }
-                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)]);
-            }
-        }
-    }
-    set_bnd(N, b, d);
-}
-
-/**
- * The main density step function
-*/
-void dens_step(int N, float * x, float * x0, float * u, float * v, float * w, float diff, float dt){
-    add_source(N, x, x0, dt);
-    SWAP(x0, x);
-    diffuse(N, 0, x, x0, diff, dt);
-    SWAP(x0, x);
-    advect(N, 0, x, x0, u, v, w, dt);
-}
-
-/**
- * The main velocity step function
-*/
-void vel_step(int N, float * u, float * v, float * w, float * u0, float * v0, float * w0, float visc, float dt){
-    add_source(N, u, u0, dt);
-    add_source(N, v, v0, dt);
-    add_source(N, w, w0, dt);
-    SWAP(u0, u);
-    diffuse(N, 1, u, u0, visc, dt);
-    SWAP(v0, v);
-    diffuse(N, 2, v, v0, visc, dt);
-    SWAP(w0, w);
-    diffuse(N, 3, w, w0, visc, dt);
-    project(N, u, v, w, u0, v0);
-    SWAP(u0, u);
-    SWAP(v0, v);
-    SWAP(w0, w);
-    advect(N, 1, u, u0, u0, v0, w0, dt);
-    advect(N, 2, v, v0, u0, v0, w0, dt);
-    advect(N, 3, w, w0, u0, v0, w0, dt);
-    project(N, u, v, w, u0, v0);
-}
-
-//used for temporary vector storage when appropriate
-float container[16];
-
-/**
- * Projects a given array based on force vectors
-*/
-void project(int N, float * u, float * v, float * w, float * p, float * div){
-    int i, j, k;
-
-    __m256 nVector = _mm256_set1_ps(N);
-    __m256 constScalar = _mm256_set1_ps(-1.0/3.0);
-    __m256 zeroVec = _mm256_set1_ps(0);
-    __m256 vector, vector2, vector3;
-
-    for(k=1; k<N-1; k++){
-        for(j=1; j<N-1; j++){
-            i = 1;
-            //
-            //lower
-            //
-            //first part
-            vector = _mm256_loadu_ps(&u[IX(i+1,j,k)]);
-            vector = _mm256_sub_ps(vector,_mm256_loadu_ps(&u[IX(i-1,j,k)]));
-            vector = _mm256_div_ps(vector,nVector);
-            //second part
-            vector2 = _mm256_loadu_ps(&v[IX(i,j+1,k)]);
-            vector2 = _mm256_sub_ps(vector2,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
-            vector2 = _mm256_div_ps(vector2,nVector);
-            //third part
-            vector3 = _mm256_loadu_ps(&w[IX(i,j,k+1)]);
-            vector3 = _mm256_sub_ps(vector3,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
-            vector3 = _mm256_div_ps(vector3,nVector);
-            //multiply and finalize
-            vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
-            vector = _mm256_mul_ps(vector,constScalar);
-            //store
-            _mm256_storeu_ps(&div[IX(i,j,k)],vector);
-            _mm256_storeu_ps(&p[IX(i,j,k)],zeroVec);
-            i = 9;
-            //
-            //upper
-            //
-            //first part
-            vector = _mm256_loadu_ps(&u[IX(i+1,j,k)]);
-            vector = _mm256_sub_ps(vector,_mm256_loadu_ps(&u[IX(i-1,j,k)]));
-            vector = _mm256_div_ps(vector,nVector);
-            //second part
-            vector2 = _mm256_loadu_ps(&v[IX(i,j+1,k)]);
-            vector2 = _mm256_sub_ps(vector2,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
-            vector2 = _mm256_div_ps(vector2,nVector);
-            //third part
-            vector3 = _mm256_loadu_ps(&w[IX(i,j,k+1)]);
-            vector3 = _mm256_sub_ps(vector3,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
-            vector3 = _mm256_div_ps(vector3,nVector);
-            //multiply and finalize
-            vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
-            vector = _mm256_mul_ps(vector,constScalar);
-            //store
-            _mm256_storeu_ps(&div[IX(i,j,k)],vector);
-            _mm256_storeu_ps(&p[IX(i,j,k)],zeroVec);
-        }
-    }
-    
-    set_bnd(N, 0, div);
-    set_bnd(N, 0, p);
-
-    lin_solve(N, 0, p, div, 1, 6);
-
-
-    constScalar = _mm256_set1_ps(0.5f*N);
-    for ( k=1 ; k<N-1 ; k++ ) {
-        for ( j=1 ; j<N-1 ; j++ ) {
-            //
-            //v
-            //
-            //lower
-            vector =  _mm256_loadu_ps(&p[IX(1+1,j,k)]);
-            vector2 = _mm256_loadu_ps(&p[IX(1-1,j,k)]);
-            vector =  _mm256_sub_ps(vector,vector2);
-            vector = _mm256_mul_ps(vector,constScalar);
-            vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(1,j,k)]),vector);
-            _mm256_storeu_ps(&u[IX(1,j,k)],vector);
-            //upper
-            vector =  _mm256_loadu_ps(&p[IX(9+1,j,k)]);
-            vector2 = _mm256_loadu_ps(&p[IX(9-1,j,k)]);
-            vector =  _mm256_sub_ps(vector,vector2);
-            vector = _mm256_mul_ps(vector,constScalar);
-            vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(9,j,k)]),vector);
-            _mm256_storeu_ps(&u[IX(9,j,k)],vector);
-            //
-            //v
-            //
-            //lower
-            vector =  _mm256_loadu_ps(&p[IX(1,j+1,k)]);
-            vector2 = _mm256_loadu_ps(&p[IX(1,j-1,k)]);
-            vector =  _mm256_sub_ps(vector,vector2);
-            vector = _mm256_mul_ps(vector,constScalar);
-            vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(1,j,k)]),vector);
-            _mm256_storeu_ps(&v[IX(1,j,k)],vector);
-            //upper
-            vector =  _mm256_loadu_ps(&p[IX(9,j+1,k)]);
-            vector2 = _mm256_loadu_ps(&p[IX(9,j-1,k)]);
-            vector =  _mm256_sub_ps(vector,vector2);
-            vector = _mm256_mul_ps(vector,constScalar);
-            vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(9,j,k)]),vector);
-            _mm256_storeu_ps(&v[IX(9,j,k)],vector);
-            //
-            //w
-            //
-            //lower
-            vector =  _mm256_loadu_ps(&p[IX(1,j,k+1)]);
-            vector2 = _mm256_loadu_ps(&p[IX(1,j,k-1)]);
-            vector =  _mm256_sub_ps(vector,vector2);
-            vector = _mm256_mul_ps(vector,constScalar);
-            vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(1,j,k)]),vector);
-            _mm256_storeu_ps(&w[IX(1,j,k)],vector);
-            //upper
-            vector =  _mm256_loadu_ps(&p[IX(9,j,k+1)]);
-            vector2 = _mm256_loadu_ps(&p[IX(9,j,k-1)]);
-            vector =  _mm256_sub_ps(vector,vector2);
-            vector = _mm256_mul_ps(vector,constScalar);
-            vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(9,j,k)]),vector);
-            _mm256_storeu_ps(&w[IX(9,j,k)],vector);
-        }
-    }
-    
-    set_bnd(N, 1, u);
-    set_bnd(N, 2, v);
-    set_bnd(N, 3, w);
-}
-
-/**
- * Solves a linear system of equations in a vectorized manner
-*/
-void lin_solve(int N, int b, float* x, float* x0, float a, float c){
-    int i, j, k, l, m;
-    __m256 aScalar = _mm256_set1_ps(a);
-    __m256 cScalar = _mm256_set1_ps(c);
-    // iterate the solver
-    for ( l=0 ; l<LINEARSOLVERTIMES ; l++ ) {
-        // update for each cell
-        for(k=1; k<N-1; k++){
-            for(j=1; j<N-1; j++){
-                int n = 0;
-                //solve as much as possible vectorized
-                for(i = 1; i < N-1; i=i+8){
-                    __m256 vector = _mm256_loadu_ps(&x[IX(i-1,j,k)]);
-                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i+1,j,k)]));
-                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j-1,k)]));
-                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j+1,k)]));
-                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k-1)]));
-                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k+1)]));
-                    vector = _mm256_mul_ps(vector,aScalar);
-                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x0[IX(i,j,k)]));
-                    vector = _mm256_div_ps(vector,cScalar);
-                    _mm256_storeu_ps(&x[IX(i,j,k)],vector);
-                }
-                //If there is any leftover, perform manual solving
-                if(i>N-1){
-                    for(i=i-8; i < N-1; i++){
-                        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;
-                    }
-                }
-            }
-        }
-        set_bnd(N, b, x);
-    }
-}
-
-/**
- * Sets the bounds of the simulation
-*/
-void set_bnd(int N, int b, float * target){
-    int DIM = N;
-    for(int x=1; x < DIM-1; x++){
-        for(int y = 1; y < DIM-1; y++){
-            //((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[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[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[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[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);
-}

src/main/c/includes/chunk.h

@@ -0,0 +1,19 @@
+#ifndef CHUNK_H
+#define CHUNK_H
+
+/**
+ * A chunk
+*/
+typedef struct {
+    float * d[27];
+    float * d0[27];
+    float * u[27];
+    float * v[27];
+    float * w[27];
+    float * u0[27];
+    float * v0[27];
+    float * w0[27];
+    int chunkMask;
+} Chunk;
+
+#endif

src/main/c/includes/chunkmask.h

@@ -0,0 +1,48 @@
+#include <stdint.h>
+
+#ifndef CHUNKMASK_H
+#define CHUNKMASK_H
+
+#define CENTER_LOC 13
+
+#define CHUNK_222 1
+#define CHUNK_122 2
+#define CHUNK_022 4
+#define CHUNK_212 8
+#define CHUNK_112 16
+#define CHUNK_012 32
+#define CHUNK_202 64
+#define CHUNK_102 128
+#define CHUNK_002 256
+
+#define CHUNK_221 512
+#define CHUNK_121 1024
+#define CHUNK_021 2048
+#define CHUNK_211 4096
+#define CHUNK_111 8192
+#define CHUNK_011 16384
+#define CHUNK_201 32768
+#define CHUNK_101 65536
+#define CHUNK_001 131072
+
+#define CHUNK_220 262144
+#define CHUNK_120 524288
+#define CHUNK_020 1048576
+#define CHUNK_210 2097152
+#define CHUNK_110 4194304
+#define CHUNK_010 8388608
+#define CHUNK_200 16777216
+#define CHUNK_100 33554432
+#define CHUNK_000 67108864
+
+extern const uint32_t CHUNK_INDEX_ARR[];
+
+
+//control offsetting the advect sampler location if a valid neighbor chunk is hit
+extern const char CHUNK_NORMALIZE_U[];
+
+extern const char CHUNK_NORMALIZE_V[];
+
+extern const char CHUNK_NORMALIZE_W[];
+
+#endif

src/main/c/includes/electrosphere_FluidSim.h

@@ -20,10 +20,10 @@ extern "C" {
 /*
  * Class:     electrosphere_FluidSim
  * Method:    simulate
- * Signature: (IIILjava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
+ * Signature: (Ljava/util/List;F)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate
-  (JNIEnv *, jobject, jint, jint, jint, jobject, jobject, jobject, jobject, jobject, jobject, jobject, jobject, jfloat, jfloat, jfloat);
+  (JNIEnv *, jclass, jobject, jfloat);
 
 #ifdef __cplusplus
 }

src/main/c/includes/mainFunctions.h

@@ -0,0 +1,101 @@
+#ifndef MAINFUNC
+#define MAINFUNC
+
+
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    addSourceToVectors
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
+ */
+static inline void addSourceToVectors
+  (int, int, float **, float **, float **, float **, float **, float **, float, float, float);
+
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    solveVectorDiffuse
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
+ */
+static inline void solveVectorDiffuse
+  (int, int, float **, float **, float **, float **, float **, float **, float, float, float);
+
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    setupProjection
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
+ */
+static inline void setupProjection
+  (
+  int N,
+  int chunk_mask,
+  float ** ur,
+  float ** vr,
+  float ** wr,
+  float ** pr,
+  float ** divr,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt);
+
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    solveProjection
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
+ */
+static inline void solveProjection
+  (int, int, float **, float **, float **, float **, float **, float **, float, float, float);
+
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    finalizeProjection
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
+ */
+static inline void finalizeProjection
+  (int, int,  float **, float **, float **, float **, float **, float **, float, float, float);
+
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    advectVectors
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
+ */
+static inline void advectVectors
+  (int, int, float **, float **, float **, float **, float **, float **, float, float, float);
+
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    addDensity
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;F)V
+ */
+static inline void addDensity
+  (int, int, float **, float **, float);
+
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    solveDiffuseDensity
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
+ */
+static inline void solveDiffuseDensity
+  (int, int, float **, float **, float **, float **, float **, float, float, float);
+
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    advectDensity
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
+ */
+static inline void advectDensity(uint32_t chunk_mask, int N, float ** d, float ** d0, float ** ur, float ** vr, float ** wr, float dt);
+
+static inline void setBoundsToNeighborsRaw
+  (
+  int N,
+  int chunk_mask,
+  int vector_dir,
+  float ** neighborArray);
+
+static inline void copyNeighborsRaw
+  (
+  int N,
+  int chunk_mask,
+  int cx,
+  int vector_dir,
+  float ** neighborArray);
+
+#endif

src/main/c/includes/simulation.h

@@ -0,0 +1,8 @@
+#ifndef SIMULATION_H
+#define SIMULATION_H
+
+#include "./chunk.h"
+
+void simulate(int numChunks, Chunk ** passedInChunks, float timestep);
+
+#endif

src/main/c/includes/utilities.h

@@ -0,0 +1,12 @@
+#include <stdint.h>
+
+#ifndef UTILITIES_H
+#define UTILITIES_H
+
+#define SWAP(x0,x) {float *tmp=x0;x0=x;x=tmp;}
+#define IX(i,j,k) ((i)+(N)*(j)+(N*N)*(k))
+#define CK(m,n,o) ((m)+(n)*(3)+(o)*(3)*(3))
+#define GET_ARR_RAW(src,i) src[i]
+#define ARR_EXISTS(chunk_mask,m,n,o) (chunk_mask & CHUNK_INDEX_ARR[CK(m,n,o)]) > 0
+
+#endif

src/main/c/lib/stb

@@ -0,0 +1 @@
+Subproject commit ae721c50eaf761660b4f90cc590453cdb0c2acd0

src/main/c/src/chunkmask.c

@@ -0,0 +1,61 @@
+#include <stdint.h>
+#include "../includes/utilities.h"
+#include "../includes/chunkmask.h"
+
+const uint32_t CHUNK_INDEX_ARR[] = {
+    CHUNK_000, CHUNK_100, CHUNK_200,
+    CHUNK_010, CHUNK_110, CHUNK_210,
+    CHUNK_020, CHUNK_120, CHUNK_220,
+
+    CHUNK_001, CHUNK_101, CHUNK_201,
+    CHUNK_011, CHUNK_111, CHUNK_211,
+    CHUNK_021, CHUNK_121, CHUNK_221,
+
+    CHUNK_002, CHUNK_102, CHUNK_202,
+    CHUNK_012, CHUNK_112, CHUNK_212,
+    CHUNK_022, CHUNK_122, CHUNK_222,
+};
+
+
+//control offsetting the advect sampler location if a valid neighbor chunk is hit
+const char CHUNK_NORMALIZE_U[] = {
+    1, 0, -1,
+    1, 0, -1,
+    1, 0, -1,
+
+    1, 0, -1,
+    1, 0, -1,
+    1, 0, -1,
+
+    1, 0, -1,
+    1, 0, -1,
+    1, 0, -1,
+};
+
+const char CHUNK_NORMALIZE_V[] = {
+     1,  1,  1,
+     0,  0,  0,
+    -1, -1, -1,
+
+    1,  1,  1,
+     0,  0,  0,
+    -1, -1, -1,
+
+    1,  1,  1,
+     0,  0,  0,
+    -1, -1, -1,
+};
+
+const char CHUNK_NORMALIZE_W[] = {
+    1, 1, 1,
+    1, 1, 1,
+    1, 1, 1,
+
+    0, 0, 0,
+    0, 0, 0,
+    0, 0, 0,
+
+    -1, -1, -1,
+    -1, -1, -1,
+    -1, -1, -1,
+};

src/main/c/src/densitystep.c

@@ -0,0 +1,230 @@
+#include <stdio.h>
+#include <immintrin.h>
+#include <stdint.h>
+
+#include "../includes/utilities.h"
+#include "../includes/chunkmask.h"
+
+
+/**
+ * Adds density to the density array
+*/
+static inline void addDensity(
+  int N,
+  int chunk_mask,
+  float ** d,
+  float ** d0,
+  float dt
+){
+    int i;
+    int size=N*N*N;
+    float * x = GET_ARR_RAW(d,CENTER_LOC);
+    float * s = GET_ARR_RAW(d0,CENTER_LOC);
+    for(i=0; i<size; i++){
+        x[i] += dt*s[i];
+    }
+}
+
+/*
+ * A single iteration of the jacobi to solve density diffusion
+ */
+static inline void solveDiffuseDensity(
+  int N,
+  int chunk_mask,
+  float ** d,
+  float ** d0,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+){
+    float a=dt*DIFFUSION_CONST*N*N*N;
+    float c=1+6*a;
+    int i, j, k, l, m;
+    float * x = GET_ARR_RAW(d,CENTER_LOC);
+    float * x0 = GET_ARR_RAW(d0,CENTER_LOC);
+    
+    __m256 aScalar = _mm256_set1_ps(a);
+    __m256 cScalar = _mm256_set1_ps(c);
+
+    //transform u direction
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&x[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x0[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&x[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    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;
+                }
+            }
+        }
+    }
+}
+
+/**
+ * Advects the density based on the vectors
+*/
+static inline void advectDensity(uint32_t chunk_mask, int N, float ** d, float ** d0, float ** ur, float ** vr, float ** wr, float dt){
+    int i, j, k, i0, j0, k0, i1, j1, k1;
+    int m,n,o;
+    float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz;
+    
+    dtx=dty=dtz=dt*N;
+
+    float * center_d = GET_ARR_RAW(d,CENTER_LOC);
+    float * center_d0 = GET_ARR_RAW(d0,CENTER_LOC);
+
+    float * u = GET_ARR_RAW(ur,CENTER_LOC);
+    float * v = GET_ARR_RAW(vr,CENTER_LOC);
+    float * w = GET_ARR_RAW(wr,CENTER_LOC);
+
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            for(i=1; i<N-1; i++){
+                center_d0 = GET_ARR_RAW(d0,CENTER_LOC);
+                //calculate location to pull from
+                x = i-dtx*u[IX(i,j,k)];
+                y = j-dty*v[IX(i,j,k)];
+                z = k-dtz*w[IX(i,j,k)];
+
+                m = n = o = 1;
+
+                if(x < 1){ m -= 1; }
+                if(x >= N-1){ m += 1; }
+                if(y < 1){ n -= 1; }
+                if(y >= N-1){ n += 1; }
+                if(z < 1){ o -= 1; }
+                if(z >= N-1){ o += 1; }
+
+                //If the out of bounds coordinate is in bounds for a neighbor chunk, use that chunk as source instead
+                // if(CK(m,n,o) != CENTER_LOC){
+                //     printf("Looking in border chunk\n");
+                // }
+                // if(x > 16){
+                //     printf("%f %d %d %d\n",m,n,o);
+                // }
+                // if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
+                //     // printf("Hit other chunk\n");
+                //     d0 = GET_ARR(env,jrd0,CK(m,n,o));
+                //     x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * (N-1);
+                //     y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * (N-1);
+                //     z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * (N-1);
+                // }
+
+                if(x < 0.001f){
+                    //cases to consider:
+                    //m = 0, x = -10
+                    //m = 2, x = 0.01
+                    x=0.001f;
+                    i0=(int)0;
+                    i1=1;
+                    s0 = 0.999f;
+                    s1 = 0.001f;
+                } else if(x >= N - 1){
+                    //cases to consider:
+                    //m = 0, x = 17.01
+                    //m = 2, x = 20
+                    x = N-1;
+                    i0=(int)N-2;
+                    i1=N-1;
+                    s0 = 0.001f;
+                    s1 = 0.999f;
+                } else {
+                    i0=(int)x;
+                    i1=i0+1;
+                    s1 = x-i0;
+                    s0 = 1-s1;
+                }
+
+                //clamp location within chunk
+                // if (x<0.5f) x=0.5f;
+                // if (x>N+0.5f) x=N+0.5f;
+                if (y<0.5f) y=0.5f;
+                if (y>N+0.5f) y=N+0.5f;
+                if (z<0.5f) z=0.5f;
+                if (z>N+0.5f) z=N+0.5f;
+
+                //get actual indices
+                // i0=(int)x;
+                // i1=i0+1;
+                j0=(int)y;
+                j1=j0+1;
+                k0=(int)z;
+                k1=k0+1;
+
+                //calculate percentage of each index
+                // s1 = x-i0;
+                // s0 = 1-s1;
+                t1 = y-j0;
+                t0 = 1-t1;
+                u1 = z-k0;
+                u0 = 1-u1;
+
+                if(i0 >= N){
+                    i0 = N - 1;
+                }
+                // if(i0 < 0){
+                //     i0 = 0;
+                // }
+                if(j0 >= N){
+                    j0 = N - 1;
+                }
+                // if(j0 < 0){
+                //     j0 = 0;
+                // }
+                if(k0 >= N){
+                    k0 = N - 1;
+                }
+                // if(k0 < 0){
+                //     k0 = 0;
+                // }
+                if(i1 >= N){
+                    i1 = N - 1;
+                }
+                // if(i1 < 0){
+                //     i1 = 0;
+                // }
+                if(j1 >= N){
+                    j1 = N - 1;
+                }
+                // if(j1 < 0){
+                //     j1 = 0;
+                // }
+                if(k1 >= N){
+                    k1 = N - 1;
+                }
+                // if(k1 < 0){
+                //     k1 = 0;
+                // }
+                center_d[IX(i,j,k)] = 
+                s0*(
+                    t0*u0*center_d0[IX(i0,j0,k0)]+
+                    t1*u0*center_d0[IX(i0,j1,k0)]+
+                    t0*u1*center_d0[IX(i0,j0,k1)]+
+                    t1*u1*center_d0[IX(i0,j1,k1)]
+                )+
+                s1*(
+                    t0*u0*center_d0[IX(i1,j0,k0)]+
+                    t1*u0*center_d0[IX(i1,j1,k0)]+
+                    t0*u1*center_d0[IX(i1,j0,k1)]+
+                    t1*u1*center_d0[IX(i1,j1,k1)]
+                );
+            }
+        }
+    }
+}

src/main/c/src/fluidsim.c

@@ -0,0 +1,542 @@
+#include <stdint.h>
+#include "../includes/utilities.h"
+#include "../includes/chunkmask.h"
+#include "../includes/electrosphere_FluidSim.h"
+#include "../includes/mainFunctions.h"
+#include "../includes/chunk.h"
+#include "../includes/simulation.h"
+
+#include "./chunkmask.c"
+#include "./velocitystep.c"
+#include "./densitystep.c"
+
+#ifndef SAVE_STEPS
+#define SAVE_STEPS 0
+#endif
+
+
+#define DIM 18
+#define LINEARSOLVERTIMES 20
+#define REALLY_SMALL_VALUE 0.00001
+
+#define DIFFUSION_CONSTANT 0.00001
+#define VISCOSITY_CONSTANT 0.00001
+
+char fileNameBuff[50];
+
+//all chunks
+Chunk ** chunks = NULL;
+
+//jni help:
+//https://stackoverflow.com/questions/39823375/clarification-about-getfieldid
+
+static inline void saveStep(float * values, const char * name);
+
+void simulate(
+    int numChunks,
+    Chunk ** passedInChunks,
+    jfloat timestep
+){
+    chunks = passedInChunks;
+
+    // printf("%p\n",chunks[0].d);
+    saveStep(chunks[0]->u[CENTER_LOC], "./chunks/beginU");
+    saveStep(chunks[0]->v[CENTER_LOC], "./chunks/beginV");
+    saveStep(chunks[0]->w[CENTER_LOC], "./chunks/beginW");
+    saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/beginU0");
+    saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/beginV0");
+    saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/beginW0");
+
+    //solve chunk mask
+    for(int i = 0; i < numChunks; i++){
+        Chunk * currentChunk = chunks[i];
+        addSourceToVectors(
+            DIM,
+            currentChunk->chunkMask,
+            currentChunk->u,
+            currentChunk->v,
+            currentChunk->w,
+            currentChunk->u0,
+            currentChunk->v0,
+            currentChunk->w0,
+            DIFFUSION_CONSTANT,
+            VISCOSITY_CONSTANT,
+            timestep
+        );
+        saveStep(currentChunk->u[CENTER_LOC], "./chunks/addSrcU");
+        saveStep(currentChunk->v[CENTER_LOC], "./chunks/addSrcV");
+        saveStep(currentChunk->w[CENTER_LOC], "./chunks/addSrcW");
+        saveStep(currentChunk->u0[CENTER_LOC], "./chunks/addSrcU0");
+        saveStep(currentChunk->v0[CENTER_LOC], "./chunks/addSrcV0");
+        saveStep(currentChunk->w0[CENTER_LOC], "./chunks/addSrcW0");
+    }
+    //swap all vector fields
+    {
+        //swap vector fields
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+
+            float * tmpArr;
+            for(int j = 0; j < 27; j++){
+                tmpArr = currentChunk->u[j];
+                currentChunk->u[j] = currentChunk->u0[j];
+                currentChunk->u0[j] = tmpArr;
+            }
+            for(int j = 0; j < 27; j++){
+                tmpArr = currentChunk->v[j];
+                currentChunk->v[j] = currentChunk->v0[j];
+                currentChunk->v0[j] = tmpArr;
+            }
+            for(int j = 0; j < 27; j++){
+                tmpArr = currentChunk->w[j];
+                currentChunk->w[j] = currentChunk->w0[j];
+                currentChunk->w0[j] = tmpArr;
+            }
+        }
+        //copy neighbors
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w0);
+        }
+    }
+    saveStep(chunks[0]->u[CENTER_LOC], "./chunks/swapU");
+    saveStep(chunks[0]->v[CENTER_LOC], "./chunks/swapV");
+    saveStep(chunks[0]->w[CENTER_LOC], "./chunks/swapW");
+    saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/swapU0");
+    saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/swapV0");
+    saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/swapW0");
+    // printf("after swap vecs u\n");
+    // printLayer(chunks[0]->u[CENTER_LOC],targetLayer);
+    // printf("after swap vecs u0\n");
+    // printLayer(chunks[0]->u0[CENTER_LOC],targetLayer);
+    //solve vector diffusion
+    {
+        for(int l = 0; l < LINEARSOLVERTIMES; l++){
+            //solve vector diffusion
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                solveVectorDiffuse(DIM,currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
+            }
+            if(SAVE_STEPS){
+                sprintf(fileNameBuff, "./chunks/diffuseUStep%dx", l);
+                saveStep(chunks[0]->u[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseUStep%dx0", l);
+                saveStep(chunks[0]->u0[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseVStep%dx", l);
+                saveStep(chunks[0]->v[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseVStep%dx0", l);
+                saveStep(chunks[0]->v0[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseWStep%dx", l);
+                saveStep(chunks[0]->w[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseWStep%dx0", l);
+                saveStep(chunks[0]->w0[CENTER_LOC], fileNameBuff);
+            }
+            //update array for vectors
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u);
+                setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v);
+                setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,3,currentChunk->w);
+                // setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u0);
+                // setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v0);
+                // setBoundsToNeighborsRaw(DIM,chunkMask,3,currentChunk->w0);
+                copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->u);
+                copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->v);
+                copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->w);
+                copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->u0);
+                copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->v0);
+                copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->w0);
+            }
+            if(SAVE_STEPS){
+                sprintf(fileNameBuff, "./chunks/diffuseUStep%dxBnd", l);
+                saveStep(chunks[0]->u[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseUStep%dx0Bnd", l);
+                saveStep(chunks[0]->u0[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseVStep%dxBnd", l);
+                saveStep(chunks[0]->v[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseVStep%dx0Bnd", l);
+                saveStep(chunks[0]->v0[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseWStep%dxBnd", l);
+                saveStep(chunks[0]->w[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/diffuseWStep%dx0Bnd", l);
+                saveStep(chunks[0]->w0[CENTER_LOC], fileNameBuff);
+            }
+        }
+    }
+    saveStep(chunks[0]->u[CENTER_LOC], "./chunks/diffuseU");
+    saveStep(chunks[0]->v[CENTER_LOC], "./chunks/diffuseV");
+    saveStep(chunks[0]->w[CENTER_LOC], "./chunks/diffuseW");
+    saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/diffuseU0");
+    saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/diffuseV0");
+    saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/diffuseW0");
+    //solve projection
+    {
+        //update array for vectors
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            // setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u);
+            // setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v);
+            // setBoundsToNeighborsRaw(DIM,chunkMask,3,currentChunk->w);
+            // setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u0);
+            // setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v0);
+        }
+        //setup projection
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setupProjection(DIM,currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
+        }
+
+        saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/setupProj1Div");
+        saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/setupProj1P");
+
+        //update array for vectors
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,0,currentChunk->u0);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,0,currentChunk->v0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->u0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->v0);
+        }
+
+        saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/setupProj1DivBnd");
+        saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/setupProj1PBnd");
+
+        //samples u0, v0
+        //sets u0
+        //these should have just been mirrored in the above
+        //
+        //Perform main projection solver
+        for(int l = 0; l < LINEARSOLVERTIMES; l++){
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                solveProjection(DIM,currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
+            }
+            if(SAVE_STEPS){
+                sprintf(fileNameBuff, "./chunks/proj1Step%dx", l);
+                saveStep(chunks[0]->u0[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/proj1Step%dx0", l);
+                saveStep(chunks[0]->v0[CENTER_LOC], fileNameBuff);
+            }
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,0,currentChunk->u0);
+                copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->u0);
+            }
+            if(SAVE_STEPS){
+                sprintf(fileNameBuff, "./chunks/proj1Step%dxBnd", l);
+                saveStep(chunks[0]->u0[CENTER_LOC], fileNameBuff);
+                sprintf(fileNameBuff, "./chunks/proj1Step%dx0Bnd", l);
+                saveStep(chunks[0]->v0[CENTER_LOC], fileNameBuff);
+            }
+        }
+        //samples u,v,w,u0
+        //sets u,v,w
+        //Finalize projection
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            finalizeProjection(DIM,currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
+        }
+
+        saveStep(chunks[0]->u[CENTER_LOC], "./chunks/finalizeProj1U");
+        saveStep(chunks[0]->v[CENTER_LOC], "./chunks/finalizeProj1V");
+        saveStep(chunks[0]->w[CENTER_LOC], "./chunks/finalizeProj1W");
+        // exit(0);
+        //set boundaries a final time for u,v,w
+        //...
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,3,currentChunk->w);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u0);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v0);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,3,currentChunk->w0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w0);
+        }
+    }
+    saveStep(chunks[0]->u[CENTER_LOC], "./chunks/projU");
+    saveStep(chunks[0]->v[CENTER_LOC], "./chunks/projV");
+    saveStep(chunks[0]->w[CENTER_LOC], "./chunks/projW");
+    saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/projU0");
+    saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/projV0");
+    saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/projW0");
+    // exit(0);
+    //swap all vector fields
+    {
+        //swap vector fields
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+
+            float * tmpArr;
+            for(int j = 0; j < 27; j++){
+                tmpArr = currentChunk->u[j];
+                currentChunk->u[j] = currentChunk->u0[j];
+                currentChunk->u0[j] = tmpArr;
+            }
+            for(int j = 0; j < 27; j++){
+                tmpArr = currentChunk->v[j];
+                currentChunk->v[j] = currentChunk->v0[j];
+                currentChunk->v0[j] = tmpArr;
+            }
+            for(int j = 0; j < 27; j++){
+                tmpArr = currentChunk->w[j];
+                currentChunk->w[j] = currentChunk->w0[j];
+                currentChunk->w0[j] = tmpArr;
+            }
+        }
+        //copy neighbors
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w0);
+        }
+    }
+    saveStep(chunks[0]->u[CENTER_LOC], "./chunks/swap2U");
+    saveStep(chunks[0]->v[CENTER_LOC], "./chunks/swap2V");
+    saveStep(chunks[0]->w[CENTER_LOC], "./chunks/swap2W");
+    saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/swap2U0");
+    saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/swap2V0");
+    saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/swap2W0");
+    //advect vectors across boundaries
+    {
+        //update border arrs
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,3,currentChunk->w);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u0);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v0);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,3,currentChunk->w0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w0);
+        }
+        //advect
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            advectVectors(DIM,currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
+        }
+        //update neighbor arr
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,3,currentChunk->w);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w);
+        }
+    }
+    saveStep(chunks[0]->u[CENTER_LOC], "./chunks/advectU");
+    saveStep(chunks[0]->v[CENTER_LOC], "./chunks/advectV");
+    saveStep(chunks[0]->w[CENTER_LOC], "./chunks/advectW");
+    saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/advectU0");
+    saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/advectV0");
+    saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/advectW0");
+    //solve projection
+    {
+        //update array for vectors
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,3,currentChunk->w);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u0);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v0);
+        }
+        //setup projection
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setupProjection(DIM,currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
+        }
+        //update array for vectors
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u0);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v0);
+        }
+        //samples u0, v0
+        //sets u0
+        //these should have just been mirrored in the above
+        //
+        //Perform main projection solver
+        for(int l = 0; l < LINEARSOLVERTIMES; l++){
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                solveProjection(DIM,currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
+            }
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,0,currentChunk->u0);
+                copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->u0);
+            }
+        }
+        //samples u,v,w,u0
+        //sets u,v,w
+        //Finalize projection
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            finalizeProjection(DIM,currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
+        }
+        //set boundaries a final time for u,v,w
+        //...
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,3,currentChunk->w);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,1,currentChunk->u0);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,2,currentChunk->v0);
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,3,currentChunk->w0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,1,currentChunk->u0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,2,currentChunk->v0);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,3,currentChunk->w0);
+        }
+    }
+
+
+
+
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
+
+
+
+    //add density
+    {
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            addDensity(DIM,currentChunk->chunkMask,currentChunk->d,currentChunk->d0,timestep);
+        }
+    }
+    //swap all density arrays
+    {
+        //swap vector fields
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+
+            float * tmpArr;
+            for(int j = 0; j < 27; j++){
+                tmpArr = currentChunk->d[j];
+                currentChunk->d[j] = currentChunk->d0[j];
+                currentChunk->d0[j] = tmpArr;
+            }
+        }
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->d);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->d0);
+        }
+    }
+    //diffuse density
+    {
+        for(int l = 0; l < LINEARSOLVERTIMES; l++){
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                solveDiffuseDensity(DIM,currentChunk->chunkMask,currentChunk->d,currentChunk->d0,currentChunk->u,currentChunk->v,currentChunk->w,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
+            }
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,0,currentChunk->d);
+            }
+        }
+    }
+    //swap all density arrays
+    {
+        //swap vector fields
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            float * tmpArr;
+            for(int j = 0; j < 27; j++){
+                tmpArr = currentChunk->d[j];
+                currentChunk->d[j] = currentChunk->d0[j];
+                currentChunk->d0[j] = tmpArr;
+            }
+        }
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->d);
+            copyNeighborsRaw(DIM,currentChunk->chunkMask,0,0,currentChunk->d0);
+        }
+    }
+    //advect density
+    {
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            advectDensity(currentChunk->chunkMask,DIM,currentChunk->d,currentChunk->d0,currentChunk->u,currentChunk->v,currentChunk->w,timestep);
+        }
+    }
+    //mirror densities
+    {
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            setBoundsToNeighborsRaw(DIM,currentChunk->chunkMask,0,currentChunk->d);
+        }
+    }
+}
+
+
+static inline void saveStep(float * values, const char * name){
+    if(SAVE_STEPS){
+        FILE *fp;
+        int N = DIM;
+
+        // ... fill the array somehow ...
+
+        fp = fopen(name, "w");
+        // check for error here
+
+        for(int x = 0; x < DIM; x++){
+            for(int y = 0; y < DIM; y++){
+                for(int z = 0; z < DIM; z++){
+                    float val = values[IX(x,y,z)];
+                    if(val < REALLY_SMALL_VALUE && val > -REALLY_SMALL_VALUE){
+                        val = 0;
+                    }
+                    fprintf(fp, "%f\t", val);
+                }
+                fprintf(fp, "\n");
+            }
+            fprintf(fp, "\n");
+        }
+
+        fclose(fp);
+    }
+}

src/main/c/src/javainterface.c

@@ -0,0 +1,170 @@
+#include <jni.h>
+
+//library includes
+//include stb ds
+#define STB_DS_IMPLEMENTATION
+#include "../lib/stb/stb_ds.h"
+
+//local includes
+#include "../includes/chunk.h"
+#include "../includes/chunkmask.h"
+#include "../includes/utilities.h"
+#include "../includes/simulation.h"
+
+
+//defines
+
+
+//function defines
+#define getChunk(i) (*env)->CallObjectMethod(env,chunkList,jListGet,i)
+#define getBuffArr(buffId) (*env)->GetObjectField(env,chunkJRaw,buffId)
+#define setBuffArr(buffId,value) (*env)->SetObjectField(env,chunkJRaw,buffId,value)
+#define GET_ARR(env,src,i) (*env)->GetDirectBufferAddress(env,(*env)->GetObjectArrayElement(env,src,i))
+
+
+
+//declarations
+void readInChunks(JNIEnv * env, jobject chunkList);
+int calculateChunkMask(JNIEnv * env, jobjectArray jrx);
+uint32_t matrix_transform(JNIEnv * env, jobjectArray jrx);
+
+
+
+//the list of chunks
+Chunk ** javaChunkView = NULL;
+//the number of chunks
+int numChunks = 0;
+
+
+JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
+    JNIEnv * env,
+    jclass fluidSimClass,
+    jobject chunkList,
+    jfloat dt
+    ){
+    readInChunks(env,chunkList);
+    simulate(numChunks,javaChunkView,dt);
+}
+
+
+/**
+ * Reads chunks into the dynamic array
+*/
+void readInChunks(JNIEnv * env, jobject chunkList){
+    jclass listClass = (*env)->FindClass(env,"java/util/List");
+    jclass fluidSimClass = (*env)->FindClass(env,"electrosphere/FluidSim");
+    //JNIEnv *env, jclass clazz, const char *name, const char *sig
+    jmethodID jListSize = (*env)->GetMethodID(env, listClass, "size", "()I");
+    jmethodID jListGet  = (*env)->GetMethodID(env, listClass, "get", "(I)Ljava/lang/Object;");
+    jmethodID jListAdd  = (*env)->GetMethodID(env, listClass, "add", "(Ljava/lang/Object;)Z");
+
+    //ByteBuffer[]
+    jfieldID dJId = (*env)->GetFieldID(env,fluidSimClass,"density","[Ljava/nio/ByteBuffer;");
+    jfieldID d0JId = (*env)->GetFieldID(env,fluidSimClass,"densityAddition","[Ljava/nio/ByteBuffer;");
+    jfieldID uJId = (*env)->GetFieldID(env,fluidSimClass,"uVector","[Ljava/nio/ByteBuffer;");
+    jfieldID vJId = (*env)->GetFieldID(env,fluidSimClass,"vVector","[Ljava/nio/ByteBuffer;");
+    jfieldID wJId = (*env)->GetFieldID(env,fluidSimClass,"wVector","[Ljava/nio/ByteBuffer;");
+    jfieldID u0JId = (*env)->GetFieldID(env,fluidSimClass,"uAdditionVector","[Ljava/nio/ByteBuffer;");
+    jfieldID v0JId = (*env)->GetFieldID(env,fluidSimClass,"vAdditionVector","[Ljava/nio/ByteBuffer;");
+    jfieldID w0JId = (*env)->GetFieldID(env,fluidSimClass,"wAdditionVector","[Ljava/nio/ByteBuffer;");
+    jfieldID chunkmaskJId = (*env)->GetFieldID(env,fluidSimClass,"chunkMask","I");
+
+    //the number of chunks
+    numChunks = (*env)->CallIntMethod(env,chunkList,jListSize);
+
+    //current chunk (this)
+    jobject chunkJRaw;
+    //current chunk fields
+    jobjectArray jd;
+    jobjectArray jd0;
+    jobjectArray u;
+    jobjectArray v;
+    jobjectArray w;
+    jobjectArray u0;
+    jobjectArray v0;
+    jobjectArray w0;
+    int chunkMask;
+
+    //solve chunk mask
+    for(int i = 0; i < numChunks; i++){
+        chunkJRaw = getChunk(i);
+        chunkMask = calculateChunkMask(env,getBuffArr(dJId));
+        (*env)->SetIntField(env,chunkJRaw,chunkmaskJId,chunkMask);
+
+        Chunk * newChunk;
+        if(i >= stbds_arrlen(javaChunkView)){
+            // printf("allocate chunk %d\n",i);
+            // fflush(stdout);
+            newChunk = (Chunk *)malloc(sizeof(Chunk));
+            // printf("new chunk %p\n",newChunk);
+            // fflush(stdout);
+            stbds_arrput(javaChunkView,newChunk);
+            // printf("new chunk %p\n",chunks[i]);
+            // fflush(stdout);
+        } else {
+            newChunk = javaChunkView[i];
+            // printf("get chunk %d: %p\n",i,newChunk);
+            // fflush(stdout);
+        }
+        jd = (*env)->GetObjectField(env,chunkJRaw,dJId);
+        jd0 = (*env)->GetObjectField(env,chunkJRaw,d0JId);
+        u = (*env)->GetObjectField(env,chunkJRaw,uJId);
+        v = (*env)->GetObjectField(env,chunkJRaw,vJId);
+        w = (*env)->GetObjectField(env,chunkJRaw,wJId);
+        u0 = (*env)->GetObjectField(env,chunkJRaw,u0JId);
+        v0 = (*env)->GetObjectField(env,chunkJRaw,v0JId);
+        w0 = (*env)->GetObjectField(env,chunkJRaw,w0JId);
+        newChunk->chunkMask = chunkMask;
+        for(int j = 0; j < 27; j++){
+            if((chunkMask & CHUNK_INDEX_ARR[j]) > 0){
+                newChunk->d[j] = GET_ARR(env,jd,j);
+                newChunk->d0[j] = GET_ARR(env,jd0,j);
+                newChunk->u[j] = GET_ARR(env,u,j);
+                newChunk->v[j] = GET_ARR(env,v,j);
+                newChunk->w[j] = GET_ARR(env,w,j);
+                newChunk->u0[j] = GET_ARR(env,u0,j);
+                newChunk->v0[j] = GET_ARR(env,v0,j);
+                newChunk->w0[j] = GET_ARR(env,w0,j);
+            }
+        }
+    }
+}
+
+
+/**
+ * Calculates the bitmask for available chunks for the provided chunk's neighbor array
+*/
+int calculateChunkMask(JNIEnv * env, jobjectArray jrx){
+    return matrix_transform(env,jrx);
+}
+
+/**
+ * Calculates a mask that represents all nearby chunks that are actually accessible and exist
+*/
+uint32_t matrix_transform(JNIEnv * env, jobjectArray jrx){
+
+    //The returned value, an availability mask that contains the availability of each neighbor chunk
+    uint32_t rVal = 0;
+
+    //Add to maks for initial chunks
+    for(int i = 0; i < CENTER_LOC; i++){
+        if((*env)->GetObjectArrayElement(env,jrx,i)!=NULL){
+            rVal = rVal + 1;
+        }
+        rVal = rVal << 1;
+    }
+    //add 1 for center chunk because we already have that
+    rVal = rVal + 1;
+    rVal = rVal << 1;
+    //continue on for remaining chunks
+    for(int i = CENTER_LOC+1; i < 27; i++){
+        if((*env)->GetObjectArrayElement(env,jrx,i)!=NULL){
+            rVal = rVal + 1;
+        }
+        if(i < 26){
+            rVal = rVal << 1;
+        }
+    }
+
+    return rVal;
+}

src/main/c/src/velocitystep.c

@@ -0,0 +1,936 @@
+#include <stdio.h>
+#include <immintrin.h>
+#include <stdint.h>
+
+#include "../includes/utilities.h"
+#include "../includes/chunkmask.h"
+
+
+#define BOUND_NO_DIR 0
+#define BOUND_DIR_U 1
+#define BOUND_DIR_V 2
+#define BOUND_DIR_W 3
+
+#define SET_BOUND_IGNORE 0
+#define SET_BOUND_USE_NEIGHBOR 1
+
+#define LINEARSOLVERTIMES 20
+
+static inline void add_source(int N, float * x, float * s, float dt);
+static inline void advect(uint32_t chunk_mask, int N, int b, float ** jrd, float ** jrd0, float * u, float * v, float * w, float dt);
+
+
+/*
+ * Adds force to all vectors
+ */
+static inline void addSourceToVectors
+  (
+  int N,
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt){
+    add_source(N,GET_ARR_RAW(jru,CENTER_LOC),GET_ARR_RAW(jru0,CENTER_LOC),dt);
+    add_source(N,GET_ARR_RAW(jrv,CENTER_LOC),GET_ARR_RAW(jrv0,CENTER_LOC),dt);
+    add_source(N,GET_ARR_RAW(jrw,CENTER_LOC),GET_ARR_RAW(jrw0,CENTER_LOC),dt);
+}
+
+/**
+ * Adds from a source array to a destination array
+*/
+static inline void add_source(int N, float * x, float * s, float dt){
+	int i;
+    int size=N*N*N;
+	for(i=0; i<size; i++){
+        x[i] += dt*s[i];
+    }
+}
+
+/*
+ * Solves vector diffusion along all axis
+ */
+static inline void solveVectorDiffuse (
+  int N,
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+){
+    float a=dt*VISCOSITY_CONST*N*N*N;
+    float c=1+6*a;
+    int i, j, k, l, m;
+    float * u = GET_ARR_RAW(jru,CENTER_LOC);
+    float * v = GET_ARR_RAW(jrv,CENTER_LOC);
+    float * w = GET_ARR_RAW(jrw,CENTER_LOC);
+    float * u0 = GET_ARR_RAW(jru0,CENTER_LOC);
+    float * v0 = GET_ARR_RAW(jrv0,CENTER_LOC);
+    float * w0 = GET_ARR_RAW(jrw0,CENTER_LOC);
+    
+    __m256 aScalar = _mm256_set1_ps(a);
+    __m256 cScalar = _mm256_set1_ps(c);
+
+    //transform u direction
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&u[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u0[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&u[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    u[IX(i,j,k)] = (u0[IX(i,j,k)] + a*(u[IX(i-1,j,k)]+u[IX(i+1,j,k)]+u[IX(i,j-1,k)]+u[IX(i,j+1,k)]+u[IX(i,j,k-1)]+u[IX(i,j,k+1)]))/c;
+                }
+            }
+        }
+    }
+
+    //transform v direction
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&v[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v0[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&v[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    v[IX(i,j,k)] = (v0[IX(i,j,k)] + a*(v[IX(i-1,j,k)]+v[IX(i+1,j,k)]+v[IX(i,j-1,k)]+v[IX(i,j+1,k)]+v[IX(i,j,k-1)]+v[IX(i,j,k+1)]))/c;
+                }
+            }
+        }
+    }
+
+    //transform w direction
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&w[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w0[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&w[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    w[IX(i,j,k)] = (w0[IX(i,j,k)] + a*(w[IX(i-1,j,k)]+w[IX(i+1,j,k)]+w[IX(i,j-1,k)]+w[IX(i,j+1,k)]+w[IX(i,j,k-1)]+w[IX(i,j,k+1)]))/c;
+                }
+            }
+        }
+    }
+}
+
+/*
+ * Sets up a projection system of equations
+ */
+static inline void setupProjection(
+  int N,
+  int chunk_mask,
+  float ** ur,
+  float ** vr,
+  float ** wr,
+  float ** pr,
+  float ** divr,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+){
+    int i, j, k;
+
+    __m256 nVector = _mm256_set1_ps(N);
+    __m256 constScalar = _mm256_set1_ps(-1.0/3.0);
+    __m256 zeroVec = _mm256_set1_ps(0);
+    __m256 vector, vector2, vector3;
+
+    float * u = GET_ARR_RAW(ur,CENTER_LOC);
+    float * v = GET_ARR_RAW(vr,CENTER_LOC);
+    float * w = GET_ARR_RAW(wr,CENTER_LOC);
+
+    float * p = GET_ARR_RAW(pr,CENTER_LOC);
+    float * div = GET_ARR_RAW(divr,CENTER_LOC);
+
+    float scalar = 1.0/3.0;
+    float h = 1.0/N;
+
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            i = 1;
+            //
+            //lower
+            //
+            //first part
+            vector = _mm256_loadu_ps(&u[IX(i+1,j,k)]);
+            vector = _mm256_sub_ps(vector,_mm256_loadu_ps(&u[IX(i-1,j,k)]));
+            vector = _mm256_div_ps(vector,nVector);
+            //second part
+            vector2 = _mm256_loadu_ps(&v[IX(i,j+1,k)]);
+            vector2 = _mm256_sub_ps(vector2,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
+            vector2 = _mm256_div_ps(vector2,nVector);
+            //third part
+            vector3 = _mm256_loadu_ps(&w[IX(i,j,k+1)]);
+            vector3 = _mm256_sub_ps(vector3,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
+            vector3 = _mm256_div_ps(vector3,nVector);
+            //multiply and finalize
+            vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
+            vector = _mm256_mul_ps(vector,constScalar);
+            //store
+            _mm256_storeu_ps(&div[IX(i,j,k)],vector);
+            _mm256_storeu_ps(&p[IX(i,j,k)],zeroVec);
+            i = 9;
+            //
+            //upper
+            //
+            //first part
+            vector = _mm256_loadu_ps(&u[IX(i+1,j,k)]);
+            vector = _mm256_sub_ps(vector,_mm256_loadu_ps(&u[IX(i-1,j,k)]));
+            vector = _mm256_div_ps(vector,nVector);
+            //second part
+            vector2 = _mm256_loadu_ps(&v[IX(i,j+1,k)]);
+            vector2 = _mm256_sub_ps(vector2,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
+            vector2 = _mm256_div_ps(vector2,nVector);
+            //third part
+            vector3 = _mm256_loadu_ps(&w[IX(i,j,k+1)]);
+            vector3 = _mm256_sub_ps(vector3,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
+            vector3 = _mm256_div_ps(vector3,nVector);
+            //multiply and finalize
+            vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
+            vector = _mm256_mul_ps(vector,constScalar);
+            //store
+            _mm256_storeu_ps(&div[IX(i,j,k)],vector);
+            _mm256_storeu_ps(&p[IX(i,j,k)],zeroVec);
+            
+        }
+    }
+}
+
+/*
+ * Solves a projection system of equations
+ */
+static inline void solveProjection(
+  int N,
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+){
+    int a = 1;
+    int c = 6;
+    int i, j, k, l, m;
+    __m256 aScalar = _mm256_set1_ps(a);
+    __m256 cScalar = _mm256_set1_ps(c);
+
+    float * p = GET_ARR_RAW(jru0,CENTER_LOC);
+    float * div = GET_ARR_RAW(jrv0,CENTER_LOC);
+    // update for each cell
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&p[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&div[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&p[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    p[IX(i,j,k)] = (div[IX(i,j,k)] + a*(p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c;
+                }
+            }
+        }
+    }
+}
+
+/*
+ * Finalizes a projection (subtract curl, set bounds, etc)
+ */
+static inline void finalizeProjection(
+  int N,
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+){
+    int i, j, k;
+    __m256 constScalar = _mm256_set1_ps(0.5f*N);
+    __m256 vector, vector2, vector3;
+
+    float * u = GET_ARR_RAW(jru,CENTER_LOC);
+    float * v = GET_ARR_RAW(jrv,CENTER_LOC);
+    float * w = GET_ARR_RAW(jrw,CENTER_LOC);
+
+    float * p = GET_ARR_RAW(jru0,CENTER_LOC);
+    float * div = GET_ARR_RAW(jrv0,CENTER_LOC);
+
+    for ( k=1 ; k<N-1 ; k++ ) {
+        for ( j=1 ; j<N-1 ; j++ ) {
+            //
+            //v
+            //
+            //lower
+            vector =  _mm256_loadu_ps(&p[IX(1+1,j,k)]);
+            vector2 = _mm256_loadu_ps(&p[IX(1-1,j,k)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(1,j,k)]),vector);
+            _mm256_storeu_ps(&u[IX(1,j,k)],vector);
+            //upper
+            vector =  _mm256_loadu_ps(&p[IX(9+1,j,k)]);
+            vector2 = _mm256_loadu_ps(&p[IX(9-1,j,k)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(9,j,k)]),vector);
+            _mm256_storeu_ps(&u[IX(9,j,k)],vector);
+            //
+            //v
+            //
+            //lower
+            vector =  _mm256_loadu_ps(&p[IX(1,j+1,k)]);
+            vector2 = _mm256_loadu_ps(&p[IX(1,j-1,k)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(1,j,k)]),vector);
+            _mm256_storeu_ps(&v[IX(1,j,k)],vector);
+            //upper
+            vector =  _mm256_loadu_ps(&p[IX(9,j+1,k)]);
+            vector2 = _mm256_loadu_ps(&p[IX(9,j-1,k)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(9,j,k)]),vector);
+            _mm256_storeu_ps(&v[IX(9,j,k)],vector);
+            //
+            //w
+            //
+            //lower
+            vector =  _mm256_loadu_ps(&p[IX(1,j,k+1)]);
+            vector2 = _mm256_loadu_ps(&p[IX(1,j,k-1)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(1,j,k)]),vector);
+            _mm256_storeu_ps(&w[IX(1,j,k)],vector);
+            //upper
+            vector =  _mm256_loadu_ps(&p[IX(9,j,k+1)]);
+            vector2 = _mm256_loadu_ps(&p[IX(9,j,k-1)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(9,j,k)]),vector);
+            _mm256_storeu_ps(&w[IX(9,j,k)],vector);
+        }
+    }
+}
+
+/*
+ * Advects u, v, and w
+ */
+static inline void advectVectors(
+  int N,
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+){
+    advect(chunk_mask,N,1,jru,jru0,GET_ARR_RAW(jru0,CENTER_LOC),GET_ARR_RAW(jrv0,CENTER_LOC),GET_ARR_RAW(jrw0,CENTER_LOC),dt);
+    advect(chunk_mask,N,2,jrv,jrv0,GET_ARR_RAW(jru0,CENTER_LOC),GET_ARR_RAW(jrv0,CENTER_LOC),GET_ARR_RAW(jrw0,CENTER_LOC),dt);
+    advect(chunk_mask,N,3,jrw,jrw0,GET_ARR_RAW(jru0,CENTER_LOC),GET_ARR_RAW(jrv0,CENTER_LOC),GET_ARR_RAW(jrw0,CENTER_LOC),dt);
+}
+
+/**
+ * Actually performs the advection
+*/
+static inline void advect(uint32_t chunk_mask, int N, int b, float ** jrd, float ** jrd0, float * u, float * v, float * w, float dt){
+    int i, j, k, i0, j0, k0, i1, j1, k1;
+    int m,n,o;
+    float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz;
+    
+    dtx=dty=dtz=dt*N;
+
+    float * d = GET_ARR_RAW(jrd,CENTER_LOC);
+
+    float * d0 = GET_ARR_RAW(jrd0,CENTER_LOC);
+
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            for(i=1; i<N-1; i++){
+                d0 = GET_ARR_RAW(jrd0,CENTER_LOC);
+                //calculate location to pull from
+                x = i-dtx*u[IX(i,j,k)];
+                y = j-dty*v[IX(i,j,k)];
+                z = k-dtz*w[IX(i,j,k)];
+
+                m = n = o = 1;
+
+                if(x < 0){ m += 1; }
+                else if(x >= N){ m -= 1; }
+                if(y < 0){ n += 1; }
+                else if(y >= N){ n -= 1; }
+                if(z < 0){ o += 1; }
+                else if(z >= N){ o -= 1; }
+
+                //If the out of bounds coordinate is in bounds for a neighbor chunk, use that chunk as source instead
+                if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
+
+                    // if(i == 1 && j == 1 && k == 1){
+                    //     printf("\narr indices: %d %d %d\n\n",m,n,o);
+                    // }
+
+                    //cases:
+                    //if x = 17.01, m = 2
+                    // 17 in current array is 1 in neighbor
+                    // 18 in current array is 2 in neighbor
+                    // 19 in current array is 3 in neighbor
+                    //want to sample neighbor array at 1 & 2
+                    //x becomes 1.01, sampling new array (keep in mind that 0 in the new array should contain the current array values)
+                    //modification: subtract 16
+
+                    //cases:
+                    //if x = 16.99, m = 2
+                    // 16 in current array is 0 in neighbor
+                    // 17 in current array is 1 in neighbor
+                    // 18 in current array is 2 in neighbor
+                    // 19 in current array is 3 in neighbor
+                    //want to sample current array still
+                    //x becomes 1.01, sampling new array (keep in mind that 0 in the new array should contain the current array values)
+                    //modification: no modification
+
+                    //if x = 0.01, m = 0
+                    // 0 in current array is 16 in neighbor
+                    //-1 in current array is 15 in neighbor
+                    //-2 in current array is 14 in neighbor
+                    //want to sample current array still
+                    //x becomes 15.01, sampling new array (keep in mind that 17 in the new array should contain the current array)
+                    //modification: no modification
+
+                    //if x = -0.01, m = 0
+                    // 0 in current array is 16 in neighbor
+                    //-1 in current array is 15 in neighbor
+                    //-2 in current array is 14 in neighbor
+                    //want to sample -1 & 0, so i0 becomes 15
+                    //x becomes 15.99, sampling new array (keep in mind that 17 in the new array should contain the current array)
+                    //modification: add 16
+
+                    //if x = -2, m = 0
+                    // 0 in current array is 16 in neighbor
+                    //-1 in current array is 15 in neighbor
+                    //-2 in current array is 14 in neighbor
+                    //x becomes 14, sampling new array (keep in mind that 17 in the new array should contain the current array)
+                    //modification: add 16
+
+
+                    // printf("Hit other chunk\n");
+                    d0 = GET_ARR_RAW(jrd0,CK(m,n,o));
+                    x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * (N-2);
+                    // printf("%d => %f\n",m,x);
+                    y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * (N-2);
+                    z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * (N-2);
+                }
+
+                //clamp location within chunk
+                //get indices, and calculate percentage to pull from each index
+                if(x < 0.001f){
+                    //cases to consider:
+                    //m = 0, x = -10
+                    //m = 2, x = 0.01
+                    x=0.001f;
+                    i0=(int)0;
+                    i1=1;
+                    s0 = 0.999f;
+                    s1 = 0.001f;
+                } else if(x > N - 1){
+                    //cases to consider:
+                    //m = 0, x = 17.01
+                    //m = 2, x = 20
+                    x = N-1;
+                    i0=(int)N-2;
+                    i1=N-1;
+                    s0 = 0.001f;
+                    s1 = 0.999f;
+                } else {
+                    i0=(int)x;
+                    i1=i0+1;
+                    s1 = x-i0;
+                    s0 = 1-s1;
+                }
+
+                if(y < 0.001f){
+                    //cases to consider:
+                    //m = 0, x = -10
+                    //m = 2, x = 0.01
+                    y=0.001f;
+                    j0=(int)0;
+                    j1=1;
+                    t0 = 0.999f;
+                    t1 = 0.001f;
+                } else if(y > N - 1){
+                    //cases to consider:
+                    //m = 0, x = 17.01
+                    //m = 2, x = 20
+                    y = N-1;
+                    j0=(int)N-2;
+                    j1=N-1;
+                    t0 = 0.001f;
+                    t1 = 0.999f;
+                } else {
+                    j0=(int)y;
+                    j1=j0+1;
+                    t1 = y-j0;
+                    t0 = 1-t1;
+                }
+
+
+                if(z < 0.001f){
+                    //cases to consider:
+                    //m = 0, x = -10
+                    //m = 2, x = 0.01
+                    z=0.001f;
+                    k0=(int)0;
+                    k1=1;
+                    u0 = 0.999f;
+                    u1 = 0.001f;
+                } else if(z > N - 1){
+                    //cases to consider:
+                    //m = 0, x = 17.01
+                    //m = 2, x = 20
+                    z = N-1;
+                    k0=(int)N-2;
+                    k1=N-1;
+                    u0 = 0.001f;
+                    u1 = 0.999f;
+                } else {
+                    k0=(int)z;
+                    k1=k0+1;
+                    u1 = z-k0;
+                    u0 = 1-u1;
+                }
+
+                // if (x<0.001f) x=0.001f;
+                // if (x>N+0.5f) x=N+0.5f;
+                // if (y<0.001f) y=0.001f;
+                // if (y>N+0.5f) y=N+0.5f;
+                // if (z<0.001f) z=0.001f;
+                // if (z>N+0.5f) z=N+0.5f;
+
+                //get actual indices
+                // i0=(int)x;
+                // i1=i0+1;
+                // j0=(int)y;
+                // j1=j0+1;
+                // k0=(int)z;
+                // k1=k0+1;
+
+                //calculate percentage of each index
+                // s1 = x-i0;
+                // s0 = 1-s1;
+                // t1 = y-j0;
+                // t0 = 1-t1;
+                // u1 = z-k0;
+                // u0 = 1-u1;
+
+                // if(i0 >= N){
+                //     i0 = N - 1;
+                // }
+                // if(i0 < 0){
+                //     i0 = 0;
+                // }
+                if(j0 >= N){
+                    j0 = N - 1;
+                }
+                // if(j0 < 0){
+                //     j0 = 0;
+                // }
+                if(k0 >= N){
+                    k0 = N - 1;
+                }
+                // if(k0 < 0){
+                //     k0 = 0;
+                // }
+                // if(i1 >= N){
+                //     i1 = N - 1;
+                // }
+                // if(i1 < 0){
+                //     i1 = 0;
+                // }
+                if(j1 >= N){
+                    j1 = N - 1;
+                }
+                // if(j1 < 0){
+                //     j1 = 0;
+                // }
+                if(k1 >= N){
+                    k1 = N - 1;
+                }
+                // if(k1 < 0){
+                //     k1 = 0;
+                // }
+                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)]
+                );
+            }
+        }
+    }
+}
+
+/**
+ * Sets the bounds of this cube to those of its neighbor
+*/
+static inline void setBoundsToNeighborsRaw(
+  int N,
+  int chunk_mask,
+  int vector_dir,
+  float ** neighborArray
+){
+    int DIM = N;
+    float * target = GET_ARR_RAW(neighborArray,CENTER_LOC);
+    float * source;
+    //set the faces bounds
+    for(int x=1; x < DIM-1; x++){
+        for(int y = 1; y < DIM-1; y++){
+            //((x)+(DIM)*(y) + (DIM)*(DIM)*(z))
+            target[IX(0,x,y)] =         vector_dir==BOUND_DIR_U ? -target[IX(1,x,y)] : target[IX(1,x,y)];
+            target[IX(DIM-1,x,y)] =     vector_dir==BOUND_DIR_U ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)];
+            target[IX(x,0,y)] =         vector_dir==BOUND_DIR_V ? -target[IX(x,1,y)] : target[IX(x,1,y)];
+            target[IX(x,DIM-1,y)] =     vector_dir==BOUND_DIR_V ? -target[IX(x,DIM-2,y)] : target[IX(x,DIM-2,y)];
+            target[IX(x,y,0)] =         vector_dir==BOUND_DIR_W ? -target[IX(x,y,1)] : target[IX(x,y,1)];
+            target[IX(x,y,DIM-1)] =     vector_dir==BOUND_DIR_W ? -target[IX(x,y,DIM-2)] : target[IX(x,y,DIM-2)];
+        }
+    }
+    //sets the edges of the chunk
+    for(int x = 1; x < DIM-1; x++){
+        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[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[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)]));
+
+    }
+    //sets the corners of the chunk
+    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);
+}
+
+
+/**
+ * This exclusively copies neighbors to make sure zeroing out stuff doesn't break sim
+*/
+static inline void copyNeighborsRaw(
+  int N,
+  int chunk_mask,
+  int cx,
+  int vector_dir,
+  float ** neighborArray
+){
+    int DIM = N;
+    float * target = GET_ARR_RAW(neighborArray,CENTER_LOC);
+    float * source;
+
+
+    //
+    //
+    //   PLANES
+    //
+    //
+    // __m512 transferVector;// = _mm512_set1_ps(0.5*N);
+
+    //__m256 vector = _mm256_loadu_ps(&p[IX(i-1,j,k)]);
+    //vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i+1,j,k)]));
+    //vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j-1,k)]));
+    //_mm256_storeu_ps(&p[IX(i,j,k)],vector);
+    //__m256
+    //_mm256_loadu_ps
+    //_mm256_storeu_ps
+    if(ARR_EXISTS(chunk_mask,0,1,1)){
+        source = GET_ARR_RAW(neighborArray,CK(0,1,1));
+        for(int x=1; x < DIM-1; x++){
+            // transferVector = _mm512_loadu_ps(&source[IX(DIM-2,x,1)]);
+            // _mm512_storeu_ps(&target[IX(0,x,1)],_mm512_loadu_ps(&source[IX(DIM-2,x,1)]));
+            for(int y = 1; y < DIM-1; y++){
+                target[IX(0,x,y)] = source[IX(DIM-2,x,y)];
+            }
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,2,1,1)){
+        source = GET_ARR_RAW(neighborArray,CK(2,1,1));
+        for(int x=1; x < DIM-1; x++){
+            // _mm512_storeu_ps(&target[IX(DIM-1,x,1)],_mm512_loadu_ps(&source[IX(1,x,1)]));
+            for(int y = 1; y < DIM-1; y++){
+                target[IX(DIM-1,x,y)] = source[IX(1,x,y)];
+            }
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,1,0,1)){
+        source = GET_ARR_RAW(neighborArray,CK(1,0,1));
+        for(int x=1; x < DIM-1; x++){
+            for(int y = 1; y < DIM-1; y++){
+                target[IX(x,0,y)] = source[IX(x,DIM-2,y)];
+            }
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,1,2,1)){
+        source = GET_ARR_RAW(neighborArray,CK(1,2,1));
+        for(int x=1; x < DIM-1; x++){
+            for(int y = 1; y < DIM-1; y++){
+                target[IX(x,DIM-1,y)] = source[IX(x,1,y)];
+            }
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,1,1,0)){
+        source = GET_ARR_RAW(neighborArray,CK(1,1,0));
+        for(int x=1; x < DIM-1; x++){
+            for(int y = 1; y < DIM-1; y++){
+                target[IX(x,y,0)] = source[IX(x,y,DIM-2)];
+            }
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,1,1,2)){
+        source = GET_ARR_RAW(neighborArray,CK(1,1,2));
+        for(int x=1; x < DIM-1; x++){
+            for(int y = 1; y < DIM-1; y++){
+                target[IX(x,y,DIM-1)] = source[IX(x,y,1)];
+            }
+        }
+    }
+
+
+    //
+    //
+    //    EDGES
+    //
+    //
+    if(ARR_EXISTS(chunk_mask,0,0,1)){
+        source = GET_ARR_RAW(neighborArray,CK(0,0,1));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(0,0,x)] = source[IX(DIM-2,DIM-2,x)];
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,2,0,1)){
+        source = GET_ARR_RAW(neighborArray,CK(2,0,1));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(DIM-1,0,x)] = source[IX(1,DIM-2,x)];
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,0,2,1)){
+        source = GET_ARR_RAW(neighborArray,CK(0,2,1));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(0,DIM-1,x)] = source[IX(DIM-2,1,x)];
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,2,2,1)){
+        source = GET_ARR_RAW(neighborArray,CK(2,2,1));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(DIM-1,DIM-1,x)] = source[IX(1,1,x)];
+        }
+    }
+
+    //
+    //
+
+    if(ARR_EXISTS(chunk_mask,0,1,0)){
+        source = GET_ARR_RAW(neighborArray,CK(0,1,0));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(0,x,0)] = source[IX(DIM-2,x,DIM-2)];
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,2,1,0)){
+        source = GET_ARR_RAW(neighborArray,CK(2,1,0));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(DIM-1,x,0)] = source[IX(1,x,DIM-2)];
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,0,1,2)){
+        source = GET_ARR_RAW(neighborArray,CK(0,1,2));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(0,x,DIM-1)] = source[IX(DIM-2,x,1)];
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,2,1,2)){
+        source = GET_ARR_RAW(neighborArray,CK(2,1,2));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(DIM-1,x,DIM-1)] = source[IX(1,x,1)];
+        }
+    }
+
+    //
+    //
+
+    if(ARR_EXISTS(chunk_mask,1,0,0)){
+        source = GET_ARR_RAW(neighborArray,CK(1,0,0));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(x,0,0)] = source[IX(x,DIM-2,DIM-2)];
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,1,2,0)){
+        source = GET_ARR_RAW(neighborArray,CK(1,2,0));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(x,DIM-1,0)] = source[IX(x,1,DIM-2)];
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,1,0,2)){
+        source = GET_ARR_RAW(neighborArray,CK(1,0,2));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(x,0,DIM-1)] = source[IX(x,DIM-2,1)];
+        }
+    }
+
+    if(ARR_EXISTS(chunk_mask,1,2,2)){
+        source = GET_ARR_RAW(neighborArray,CK(1,2,2));
+        for(int x=1; x < DIM-1; x++){
+                target[IX(x,DIM-1,DIM-1)] = source[IX(x,1,1)];
+        }
+    }
+
+
+    //
+    //
+    //     CORNERS
+    //
+    //
+
+    if(ARR_EXISTS(chunk_mask,0,0,0)){
+        source = GET_ARR_RAW(neighborArray,CK(0,0,0));
+        target[IX(0,0,0)] = source[IX(DIM-2,DIM-2,DIM-2)];
+    }
+
+    if(ARR_EXISTS(chunk_mask,2,0,0)){
+        source = GET_ARR_RAW(neighborArray,CK(2,0,0));
+        target[IX(DIM-1,0,0)] = source[IX(1,DIM-2,DIM-2)];
+    }
+
+    if(ARR_EXISTS(chunk_mask,0,2,0)){
+        source = GET_ARR_RAW(neighborArray,CK(0,2,0));
+        target[IX(0,DIM-1,0)] = source[IX(DIM-2,1,DIM-2)];
+    }
+
+    if(ARR_EXISTS(chunk_mask,2,2,0)){
+        source = GET_ARR_RAW(neighborArray,CK(2,2,0));
+        target[IX(DIM-1,DIM-1,0)] = source[IX(1,1,DIM-2)];
+    }
+
+    //
+    //
+
+    if(ARR_EXISTS(chunk_mask,0,0,2)){
+        source = GET_ARR_RAW(neighborArray,CK(0,0,2));
+        target[IX(0,0,DIM-1)] = source[IX(DIM-2,DIM-2,1)];
+    }
+
+    if(ARR_EXISTS(chunk_mask,2,0,2)){
+        source = GET_ARR_RAW(neighborArray,CK(2,0,2));
+        target[IX(DIM-1,0,DIM-1)] = source[IX(1,DIM-2,1)];
+    }
+
+    if(ARR_EXISTS(chunk_mask,0,2,2)){
+        source = GET_ARR_RAW(neighborArray,CK(0,2,2));
+        target[IX(0,DIM-1,DIM-1)] = source[IX(DIM-2,1,1)];
+    }
+
+    if(ARR_EXISTS(chunk_mask,2,2,2)){
+        source = GET_ARR_RAW(neighborArray,CK(2,2,2));
+        target[IX(DIM-1,DIM-1,DIM-1)] = source[IX(1,1,1)];
+    }
+
+
+
+}

src/main/java/electrosphere/FluidSim.java

@@ -2,40 +2,63 @@ package electrosphere;
 
 import java.io.File;
 import java.io.IOException;
+import java.io.OutputStream;
 import java.nio.ByteBuffer;
 import java.nio.ByteOrder;
 import java.nio.FloatBuffer;
 import java.nio.file.Files;
 import java.util.LinkedList;
 import java.util.List;
-import java.util.Random;
-import java.util.Set;
 
-import org.joml.Vector2i;
-import org.lwjgl.BufferUtils;
-import org.lwjgl.PointerBuffer;
-import org.lwjgl.system.MemoryUtil;
+import org.joml.Vector3i;
+import org.lwjgl.glfw.GLFW;
 
 /**
  * Simulates a fluid via opencl
  */
 public class FluidSim {
 
+    /**
+     * Load fluid sim library
+     */
     static {
-        System.out.println(System.getProperty("user.dir"));
-        System.load(System.getProperty("user.dir") + "/shared-folder/libfluidsim.dll");
+        String osName = System.getProperty("os.name").toLowerCase();
+        System.out.println(osName);
+        if(osName.contains("win")){
+            System.load(new File("./shared-folder/libfluidsim.dll").toPath().toAbsolutePath().toString());
+        } else {
+            System.load(new File("./shared-folder/libfluidsim.so").toPath().toAbsolutePath().toString());
+        }
     }
 
     public static final int DIM = 18;
 
-    ByteBuffer x = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-    ByteBuffer x0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-    ByteBuffer u = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-    ByteBuffer v = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-    ByteBuffer w = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-    ByteBuffer u0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-    ByteBuffer v0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-    ByteBuffer w0 = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+    //                           
+    //              +-------------+             ^ Y                
+    //            / |           / |             |    _             
+    //          /   |         /   |             |    /\  Z         
+    //(0,2,0) +-----+-------+     |             |   /              
+    //        |     +-------+-----+             |  /               
+    //        |   /         |   /               | /                
+    //        | /           | /                 |/                 
+    //        +-------------+ (2,0,0)           +---------------> X
+
+    //Buffers that contain density for current frame
+    public ByteBuffer[] density = new ByteBuffer[27];
+    //Buffers that contain new density to add to the simulation
+    public ByteBuffer[] densityAddition = new ByteBuffer[27];
+    //Buffers that contain u vector directions
+    public ByteBuffer[] uVector = new ByteBuffer[27];
+    //Buffers that contain v vector directions
+    public ByteBuffer[] vVector = new ByteBuffer[27];
+    //Buffers that contain w vector directions
+    public ByteBuffer[] wVector = new ByteBuffer[27];
+    //Buffers that contain u vector directions to add to the simulation
+    public ByteBuffer[] uAdditionVector = new ByteBuffer[27];
+    //Buffers that contain v vector directions to add to the simulation
+    public ByteBuffer[] vAdditionVector = new ByteBuffer[27];
+    //Buffers that contain w vector directions to add to the simulation
+    public ByteBuffer[] wAdditionVector = new ByteBuffer[27];
 
     //The densities for every voxel for the current frame
     float[] densityArrayView = new float[DIM * DIM * DIM];
@@ -47,8 +70,10 @@ public class FluidSim {
     float[] wArrayView = new float[DIM * DIM * DIM];
     //these should be set to the 
     float[] u0ArrayView = new float[DIM * DIM * DIM];
-    float[] v0ArrayView = new float[DIM * DIM * DIM];
+    public float[] v0ArrayView = new float[DIM * DIM * DIM];
     float[] w0ArrayView = new float[DIM * DIM * DIM];
+
+    public int chunkMask = 0;
     
 
     static final float DIFFUSION_CONSTANT = 0.00001f;
@@ -58,92 +83,434 @@ public class FluidSim {
 
     static final float GRAVITY = -1000f;
     
-    public void setup(){
-        x.order(ByteOrder.LITTLE_ENDIAN);
-        x0.order(ByteOrder.LITTLE_ENDIAN);
-        u.order(ByteOrder.LITTLE_ENDIAN);
-        v.order(ByteOrder.LITTLE_ENDIAN);
-        w.order(ByteOrder.LITTLE_ENDIAN);
-        u0.order(ByteOrder.LITTLE_ENDIAN);
-        v0.order(ByteOrder.LITTLE_ENDIAN);
-        w0.order(ByteOrder.LITTLE_ENDIAN);
-        FloatBuffer xf = x.asFloatBuffer();
-        FloatBuffer uf = u.asFloatBuffer();
-        FloatBuffer vf = v.asFloatBuffer();
-        FloatBuffer wf = w.asFloatBuffer();
+    public void setup(Vector3i offset){
+        //allocate buffers for this chunk
+        density[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        densityAddition[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        uVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        vVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        wVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        uAdditionVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        vAdditionVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        wAdditionVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        //order endian-ness
+        density[13].order(ByteOrder.LITTLE_ENDIAN);
+        densityAddition[13].order(ByteOrder.LITTLE_ENDIAN);
+        uVector[13].order(ByteOrder.LITTLE_ENDIAN);
+        vVector[13].order(ByteOrder.LITTLE_ENDIAN);
+        wVector[13].order(ByteOrder.LITTLE_ENDIAN);
+        uAdditionVector[13].order(ByteOrder.LITTLE_ENDIAN);
+        vAdditionVector[13].order(ByteOrder.LITTLE_ENDIAN);
+        wAdditionVector[13].order(ByteOrder.LITTLE_ENDIAN);
+
+        //set initial values for chunk
+        FloatBuffer xf = density[13].asFloatBuffer();
+        FloatBuffer uf = uVector[13].asFloatBuffer();
+        FloatBuffer vf = vVector[13].asFloatBuffer();
+        FloatBuffer wf = wVector[13].asFloatBuffer();
 
-        Random rand = new Random(1);
         //make a cube of water in the center
         for(int i = 0; i < DIM; i++){
             for(int j = 0; j < DIM; j++){
                 for(int k = 0; k < DIM; k++){
-                    if(
-                        Math.abs(16 - i) < 4 &&
-                        Math.abs(8 - j) < 4 &&
-                        Math.abs(10 - k) < 4
-                    ){
-                        xf.put(1);
-                        uf.put(rand.nextFloat() * 0.1f);
-                        vf.put(-1f);
-                        wf.put(rand.nextFloat() * 0.1f);
-                    } else {
-                        xf.put(0);
-                        uf.put(0);
-                        vf.put(0);
-                        wf.put(0);
-                    }
+                    // if(offset.x == 0 && offset.y == 0 && offset.z == 0){
+                        if(
+                            Math.abs(16 - i) < 4 &&
+                            Math.abs(8 - j) < 4 &&
+                            Math.abs(10 - k) < 4
+                            // &&
+                            // i < 17 && i > 0 &&
+                            // j < 17 && j > 0 &&
+                            // k < 17 && k > 0
+                        ){
+                            xf.put(1);
+                            uf.put(0.1f);
+                            vf.put(-1);
+                            wf.put(0.1f);
+                        }
+                         else {
+                            xf.put(0);
+                            uf.put(0);
+                            vf.put(0);
+                            wf.put(0);
+                        }
+                    // } else {
+                    //     if(
+                    //         Math.abs(0 - i) < 5 &&
+                    //         Math.abs(j) < 5 &&
+                    //         Math.abs(0 - k) < 5 &&
+                    //         i < 17 && i > 0 &&
+                    //         j < 17 && j > 0 &&
+                    //         k < 17 && k > 0
+                    //     ){
+                    //         // xf.put(1);
+                    //         // uf.put(50);
+                    //         // vf.put(0);
+                    //         // wf.put(rand.nextFloat() * 0.1f);
+                    //     } else {
+                    //         xf.put(0);
+                    //         uf.put(0);
+                    //         vf.put(0);
+                    //         wf.put(0);
+                    //     }
+                    // }
                 }
             }
         }
     }
 
-    /**
-     * Runs a frame of the fluid simulation
-     */
-    public void simulate(int step, float timestep){
+    static int i = 0;
+    static double time = 0;
+    static double lastTime = 0;
+    public static void simChunks(FluidSim[][][] simArray, int step, float timestep){
+
+        // simArray[0][1][2].density0ArrayView[IX(10,10,3)] = 3.0f;
+        // simArray[2][1][2].density0ArrayView[IX(10,10,3)] = 3.0f;
+        // simArray[2][1][0].density0ArrayView[IX(10,10,3)] = 3.0f;
+
+        List<FluidSim> chunksToSim = new LinkedList<FluidSim>();
+        //
+        //init data for upcoming frame
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    //
+                    // Add forces and density here
+                    //
+                    simArray[x][y][z].addGravity();
+                    //
+                    //Performs main fluid simulation logic
+                    //
+                    simArray[x][y][z].writeNewStateIntoBuffers();
+                    //
+                    // add to queue
+                    //
+                    chunksToSim.add(simArray[x][y][z]);
+                }
+            }
+        }
+
+        lastTime = GLFW.glfwGetTime();
+        //
+        //simulate
+        simulateWrapper(chunksToSim,0.01f);
+        //clock
+        time = time + (GLFW.glfwGetTime() - lastTime);
+        i++;
+        if(i == 100){
+            System.out.println(time / 100.0 * 1000.0);
+        }
+
 
         //
-        // Add forces and density here
-        //
-        addGravity();
+        //Read out of buffers back into accessible arrays
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    //
+                    //Reads out the results of the fluid sim
+                    //
+                    simArray[x][y][z].readDataIntoArrays();
+                }
+            }
+        }
 
-        //
-        //Performs main fluid simulation logic
-        //
-        writeNewStateIntoBuffers();
 
-        if(x.position() > 0){
-            x.position(0);
+        if(Main.endStep == 0){
+            System.exit(0);
         }
-        if(x0.position() > 0){
-            x0.position(0);
-        }
-        if(u.position() > 0){
-            u.position(0);
-        }
-        if(v.position() > 0){
-            v.position(0);
-        }
-        if(w.position() > 0){
-            w.position(0);
-        }
-        if(u0.position() > 0){
-            u0.position(0);
-        }
-        if(v0.position() > 0){
-            v0.position(0);
-        }
-        if(w0.position() > 0){
-            w0.position(0);
-        }
-        simulate(DIM, DIM, DIM, x, x0, u, v, w, u0, v0, w0, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
 
-        //
-        //Reads out the results of the fluid sim
-        //
-        readDataIntoArrays();
+        
 
     }
+
+
+
+    //
+    //
+    //Management functions
+    //
+    //
+    /**
+     * Creates a fluid simulation array based on a given dimension set
+     * @param dimx the x dimension
+     * @param dimy the y dimension
+     * @param dimz the z dimension
+     * @return The array of fluid sim chunks
+     */
+    public static FluidSim[][][] initFluidSim(int dimx, int dimy, int dimz){
+        FluidSim[][][] simArray = new FluidSim[dimx][dimy][dimz];
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    simArray[x][y][z] = new FluidSim();
+                    simArray[x][y][z].setup(new Vector3i(x,y,z));
+                }
+            }
+        }
+        //set sim adjacencies
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    FluidSim current = simArray[x][y][z];
+
+                    for(int i = -1; i < 2; i++){
+                        for(int j = -1; j < 2; j++){
+                            for(int k = -1; k < 2; k++){
+                                if(i == j && j == k && k == 0){
+                                    continue;
+                                }
+                                if(
+                                    0 <= x + i && x + i < simArray.length &&
+                                    0 <= y + j && y + j < simArray[0].length &&
+                                    0 <= z + k && z + k < simArray[0][0].length
+                                ){
+                                    current.setNeighbor(i+1,j+1,k+1,simArray[x+i][y+j][z+k]);
+                                }
+                            }
+                        }
+                    }
+
+                }
+            }
+        }
+        return simArray;
+    }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+    static void printLayer(FluidSim[][][] simArray, int step, int layer){
+        if(step == 0){
+            for(int x = 0; x < DIM; x++){
+                for(int y = 0; y < DIM; y++){
+                    System.out.printf("%.2f\t",simArray[0][0][0].uArrayView[IX(x,layer,y)]);
+                }
+                System.out.println();
+            }
+        }
+    }
+
+    private static double sumAllDensity(FluidSim[][][] simArray){
+        double rVal = 0;
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    rVal = rVal + simArray[x][y][z].sumDensity();
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private double sumDensity(){
+        double rVal = 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++){
+                    rVal = rVal + densityArrayView[IX(x,y,z)];
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private double sumU(){
+        double rVal = 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++){
+                    rVal = rVal + Math.abs(uArrayView[IX(x,y,z)]);
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private static double sumAllU(FluidSim[][][] simArray){
+        double rVal = 0;
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    simArray[x][y][z].readDataIntoArrays();
+                    rVal = rVal + simArray[x][y][z].sumU();
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private double sumU0(){
+        double rVal = 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++){
+                    rVal = rVal + Math.abs(u0ArrayView[IX(x,y,z)]);
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private static double sumAllU0(FluidSim[][][] simArray){
+        double rVal = 0;
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    simArray[x][y][z].readDataIntoArrays();
+                    rVal = rVal + simArray[x][y][z].sumU0();
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private double sumV(){
+        double rVal = 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++){
+                    rVal = rVal + Math.abs(vArrayView[IX(x,y,z)]);
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private static double sumAllV(FluidSim[][][] simArray){
+        double rVal = 0;
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    simArray[x][y][z].readDataIntoArrays();
+                    rVal = rVal + simArray[x][y][z].sumV();
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private double sumV0(){
+        double rVal = 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++){
+                    rVal = rVal + Math.abs(v0ArrayView[IX(x,y,z)]);
+                }
+            }
+        }
+        return rVal;
+    }
+
+    private static double sumAllV0(FluidSim[][][] simArray){
+        double rVal = 0;
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    simArray[x][y][z].readDataIntoArrays();
+                    rVal = rVal + simArray[x][y][z].sumV0();
+                }
+            }
+        }
+        return rVal;
+    }
+
+
+
+    /**
+     * Main simulation function
+     * @param timestep
+     */
+    private static void simulateWrapper(List<FluidSim> chunks, float timestep){
+        simulate(chunks,timestep);
+    }
+    private static native void simulate(List<FluidSim> chunks, float timestep);
+
+
+
+    /**
+     * Dumps this chunk's density values to disk
+     * @param path the path to dump to
+     */
+    public void dumpToDisk(String path){
+        OutputStream fileOS = null;
+        try {
+            fileOS = Files.newOutputStream(new File(path).toPath());
+            ByteBuffer density = this.density[13];
+            while(density.hasRemaining()){
+                fileOS.write(density.get());
+            }
+            fileOS.close();
+        } catch (IOException e) {
+            e.printStackTrace();
+        }
+    }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 
     /**
@@ -153,25 +520,46 @@ public class FluidSim {
         //
         //Read data back into arrays
         //
-        if(x.position() > 0){
-            x.position(0);
+        if(density[13].position() > 0){
+            density[13].position(0);
         }
-        if(u.position() > 0){
-            u.position(0);
+        if(uVector[13].position() > 0){
+            uVector[13].position(0);
         }
-        if(v.position() > 0){
-            v.position(0);
+        if(vVector[13].position() > 0){
+            vVector[13].position(0);
         }
-        if(w.position() > 0){
-            w.position(0);
+        if(wVector[13].position() > 0){
+            wVector[13].position(0);
         }
+        if(uAdditionVector[13].position() > 0){
+            uAdditionVector[13].position(0);
+        }
+        if(vAdditionVector[13].position() > 0){
+            vAdditionVector[13].position(0);
+        }
+        if(wAdditionVector[13].position() > 0){
+            wAdditionVector[13].position(0);
+        }
+        FloatBuffer xFloatView = density[13].asFloatBuffer();
+        FloatBuffer uFloatView = uVector[13].asFloatBuffer();
+        FloatBuffer vFloatView = vVector[13].asFloatBuffer();
+        FloatBuffer wFloatView = wVector[13].asFloatBuffer();
+        FloatBuffer u0FloatView = uAdditionVector[13].asFloatBuffer();
+        FloatBuffer v0FloatView = vAdditionVector[13].asFloatBuffer();
+        FloatBuffer w0FloatView = wAdditionVector[13].asFloatBuffer();
+        int index = 0;
         for(int i = 0; i < DIM; i++){
             for(int j = 0; j < DIM; j++){
                 for(int k = 0; k < DIM; k++){
-                    densityArrayView[IX(i,j,k)] = x.getFloat();
-                    uArrayView[IX(i,j,k)] = u.getFloat();
-                    vArrayView[IX(i,j,k)] = v.getFloat();
-                    wArrayView[IX(i,j,k)] = w.getFloat();
+                    index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
+                    densityArrayView[index] = xFloatView.get();
+                    uArrayView[index] = uFloatView.get();
+                    vArrayView[index] = vFloatView.get();
+                    wArrayView[index] = wFloatView.get();
+                    u0ArrayView[index] = u0FloatView.get();
+                    v0ArrayView[index] = v0FloatView.get();
+                    w0ArrayView[index] = w0FloatView.get();
                 }
             }
         }
@@ -181,29 +569,31 @@ public class FluidSim {
      * Writes data from the java-side arrays into buffers that get passed into c-side
      */
     private void writeNewStateIntoBuffers(){
-        if(x0.position() > 0){
-            x0.position(0);
+        if(densityAddition[13].position() > 0){
+            densityAddition[13].position(0);
         }
-        if(u0.position() > 0){
-            u0.position(0);
+        if(uAdditionVector[13].position() > 0){
+            uAdditionVector[13].position(0);
         }
-        if(v0.position() > 0){
-            v0.position(0);
+        if(vAdditionVector[13].position() > 0){
+            vAdditionVector[13].position(0);
         }
-        if(w0.position() > 0){
-            w0.position(0);
+        if(wAdditionVector[13].position() > 0){
+            wAdditionVector[13].position(0);
         }
-        FloatBuffer x0FloatView = x0.asFloatBuffer();
-        FloatBuffer u0FloatView = u0.asFloatBuffer();
-        FloatBuffer v0FloatView = v0.asFloatBuffer();
-        FloatBuffer w0FloatView = w0.asFloatBuffer();
+        FloatBuffer x0FloatView = densityAddition[13].asFloatBuffer();
+        FloatBuffer u0FloatView = uAdditionVector[13].asFloatBuffer();
+        FloatBuffer v0FloatView = vAdditionVector[13].asFloatBuffer();
+        FloatBuffer w0FloatView = wAdditionVector[13].asFloatBuffer();
+        int index = 0;
         for(int i = 0; i < DIM; i++){
             for(int j = 0; j < DIM; j++){
                 for(int k = 0; k < DIM; k++){
-                    x0FloatView.put(density0ArrayView[IX(i,j,k)]);
-                    u0FloatView.put(u0ArrayView[IX(i,j,k)]);
-                    v0FloatView.put(v0ArrayView[IX(i,j,k)]);
-                    w0FloatView.put(w0ArrayView[IX(i,j,k)]);
+                    index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
+                    x0FloatView.put(density0ArrayView[index]);
+                    u0FloatView.put(u0ArrayView[index]);
+                    v0FloatView.put(v0ArrayView[index]);
+                    w0FloatView.put(w0ArrayView[index]);
                 }
             }
         }
@@ -213,48 +603,61 @@ public class FluidSim {
      * Adds gravity to the simulation
      */
     private void addGravity(){
+        int index = 0;
         for(int i = 0; i < DIM; i++){
             for(int j = 0; j < DIM; j++){
                 for(int k = 0; k < DIM; k++){
-                    v0ArrayView[IX(i,j,k)] = densityArrayView[IX(i,j,k)] * GRAVITY;
+                    index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
+                    u0ArrayView[index] = 0;
+                    v0ArrayView[index] = densityArrayView[index] * GRAVITY;
+                    w0ArrayView[index] = 0;
                 }
             }
         }
     }
 
-    /**
-     * The native function call to simulate a frame of fluid
-     * @param DIM_X
-     * @param DIM_Y
-     * @param DIM_Z
-     * @param x
-     * @param x0
-     * @param u
-     * @param v
-     * @param w
-     * @param u0
-     * @param v0
-     * @param w0
-     * @param DIFFUSION_CONSTANT
-     * @param VISCOSITY_CONSTANT
-     * @param timestep
-     */
-    private native void simulate(
-        int DIM_X,
-        int DIM_Y,
-        int DIM_Z,
-        ByteBuffer x,
-        ByteBuffer x0,
-        ByteBuffer u,
-        ByteBuffer v,
-        ByteBuffer w,
-        ByteBuffer u0,
-        ByteBuffer v0,
-        ByteBuffer w0,
-        float DIFFUSION_CONSTANT,
-        float VISCOSITY_CONSTANT,
-        float timestep
-    );
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
     public float[] getData(){
         return densityArrayView;
@@ -264,4 +667,51 @@ public class FluidSim {
         return ((x)+(DIM)*(y) + (DIM)*(DIM)*(z));
     }
 
+    public static final int getNeighborIndex(int x, int y, int z){
+        return x + y * 3 + z * 3 * 3;
+    }
+
+    public ByteBuffer getDensityBuffer(){
+        return density[getNeighborIndex(1,1,1)];
+    }
+
+    public ByteBuffer getDensityAdditionBuffer(){
+        return densityAddition[getNeighborIndex(1,1,1)];
+    }
+
+    public ByteBuffer getUBuffer(){
+        return uVector[getNeighborIndex(1,1,1)];
+    }
+
+    public ByteBuffer getVBuffer(){
+        return vVector[getNeighborIndex(1,1,1)];
+    }
+
+    public ByteBuffer getWBuffer(){
+        return wVector[getNeighborIndex(1,1,1)];
+    }
+
+    public ByteBuffer getUAddBuffer(){
+        return uAdditionVector[getNeighborIndex(1,1,1)];
+    }
+
+    public ByteBuffer getVAddBuffer(){
+        return vAdditionVector[getNeighborIndex(1,1,1)];
+    }
+
+    public ByteBuffer getWAddBuffer(){
+        return wAdditionVector[getNeighborIndex(1,1,1)];
+    }
+
+    public void setNeighbor(int x, int y, int z, FluidSim neighbor){
+        density[getNeighborIndex(x,y,z)] = neighbor.getDensityBuffer();
+        densityAddition[getNeighborIndex(x,y,z)] = neighbor.getDensityAdditionBuffer();
+        uVector[getNeighborIndex(x,y,z)] = neighbor.getUBuffer();
+        vVector[getNeighborIndex(x,y,z)] = neighbor.getVBuffer();
+        wVector[getNeighborIndex(x,y,z)] = neighbor.getWBuffer();
+        uAdditionVector[getNeighborIndex(x,y,z)] = neighbor.getUAddBuffer();
+        vAdditionVector[getNeighborIndex(x,y,z)] = neighbor.getVAddBuffer();
+        wAdditionVector[getNeighborIndex(x,y,z)] = neighbor.getWAddBuffer();
+    }
+
 }

src/main/java/electrosphere/Main.java

@@ -2,10 +2,13 @@ package electrosphere;
 
 import java.io.File;
 import java.io.IOException;
+import java.nio.ByteBuffer;
 import java.nio.file.Files;
-import java.nio.file.Paths;
+import java.nio.file.Path;
 import java.util.concurrent.TimeUnit;
 
+import org.joml.Vector3i;
+
 import electrosphere.render.GLFWContext;
 import electrosphere.render.Mesh;
 
@@ -14,29 +17,78 @@ import electrosphere.render.Mesh;
  */
 public class Main {
 
+    static boolean render = false;
+    public static int endStep = -1;
+
+    public static final float TIMESTEP = 0.001f;
+    
+
     public static void main(String args[]){
+
+        int dim = 3;
+        int vdim = 3;
+        int i = 0;
+        long time = 0;
+        long lastTime = 0;
+
         try {
-            GLFWContext.init();
-            FluidSim sim = new FluidSim();
-            sim.setup();
-            Mesh.meshInitially(sim);
-            int i = 0;
-            long time = 0;
-            long lastTime = 0;
+
+            if(render){
+                GLFWContext.init();
+
+                //init shader program
+                Mesh.initShaderProgram();
+            }
+
+
+            FluidSim[][][] simArray = FluidSim.initFluidSim(dim,vdim,dim);
+            
+            Mesh[][][] meshArray = null;
+            if(render){
+                meshArray = initMeshes(dim,vdim,dim,simArray);
+            }
+
+
+            //uncomment this to generate test data
+            // generateTestData();
+
             while(true){
                 try {
                     TimeUnit.MILLISECONDS.sleep(2);
                 } catch (InterruptedException e) {
                     e.printStackTrace();
                 }
+                //
+                //Transfer data
+                //
                 lastTime = System.currentTimeMillis();
-                sim.simulate(i,0.001f);
+                //
+                //Simulate
+                //
+                FluidSim.simChunks(simArray,i,TIMESTEP);
+                //
+                //Remesh
+                //
+                if(render){
+                    for(int x = 0; x < simArray.length; x++){
+                        for(int y = 0; y < simArray[0].length; y++){
+                            for(int z = 0; z < simArray[0][0].length; z++){
+                                meshArray[x][y][z].remesh();
+                            }
+                        }
+                    }
+                }
                 time = time + (System.currentTimeMillis() - lastTime);
-                Mesh.remesh(sim);
-                GLFWContext.redraw();
+                //redraw
+                if(render){
+                    GLFWContext.redraw(meshArray);
+                }
                 i++;
-                if(i == 1000){
-                    System.out.println(time / 1000.0);
+                if(i == 100){
+                    System.out.println("overall time: " + time / 100.0);
+                }
+                if(i > 3){
+                    // scan.next();
                 }
             }
         } catch(Throwable ex){
@@ -44,4 +96,139 @@ public class Main {
         }
     }
 
+    private static Mesh[][][] initMeshes(int dimx, int dimy, int dimz, FluidSim[][][] simArray){
+        Mesh[][][] meshArray = new Mesh[dimx][dimy][dimz];
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    meshArray[x][y][z] = new Mesh(simArray[x][y][z], new Vector3i(z * 16, y * 16, x * 16));
+                    meshArray[x][y][z].meshInitially();
+                }
+            }
+        }
+        return meshArray;
+    }
+
+    /**
+     * Generates test data on disk in the test resources folder
+     */
+    private static void generateTestData(){
+
+        //variables that will be used in generating data
+        int i = 0;
+        int dim = 1;
+        int length = 500;
+        FluidSim[][][] simArray;
+        int[] dims = new int[]{1,3,5};
+        int[] lengths = new int[]{1,50,100,500};
+        render = false;
+
+        //clean up existing data
+        //this is scary recursive file logic, but the innermost callback has a guard to guarantee the file is underneath the testdata path
+        while(Files.exists(new File("./src/test/resources/testdata").toPath())){
+            try {
+                Files.walk(new File("./src/test/resources/testdata").toPath()).forEach(path -> {
+                    Path parent = new File("./src/test/resources/testdata").toPath().toAbsolutePath();
+                    Path child = path.toAbsolutePath();
+                    try {
+                        TimeUnit.MILLISECONDS.sleep(1);
+                    } catch (InterruptedException e) {
+                        // TODO Auto-generated catch block
+                        e.printStackTrace();
+                    }
+                    if(child.startsWith(parent)){
+                        try {
+                            System.out.println("Deleting " + path);
+                            Files.delete(path);
+                        } catch (IOException e) {
+                            System.out.println("Failed to delete " + path);
+                            e.printStackTrace();
+                        }
+                    }
+                });
+            } catch (IOException e) {
+                // TODO Auto-generated catch block
+                e.printStackTrace();
+            }
+        }
+
+        //create folders to contain test data
+        try {
+            Files.createDirectory(new File("./src/test/resources/testdata").toPath());
+        } catch (IOException e) {
+            // TODO Auto-generated catch block
+            e.printStackTrace();
+        }
+        for(int dimIterator = 0; dimIterator < dims.length; dimIterator++){
+            dim = dims[dimIterator];
+            try {
+                Files.createDirectory(new File("./src/test/resources/testdata/" + dim + "by" + dim).toPath());
+            } catch (IOException e) {
+                // TODO Auto-generated catch block
+                e.printStackTrace();
+            }
+            for(int lengthIterator = 0; lengthIterator < lengths.length; lengthIterator++){
+                length = lengths[lengthIterator];
+                try {
+                    Files.createDirectory(new File("./src/test/resources/testdata/" + dim + "by" + dim + "/" + length + "steps").toPath());
+                } catch (IOException e) {
+                    // TODO Auto-generated catch block
+                    e.printStackTrace();
+                }
+            }
+        }
+
+        //main data generation loop
+        for(int dimIterator = 0; dimIterator < dims.length; dimIterator++){
+            for(int lengthIterator = 0; lengthIterator < lengths.length; lengthIterator++){
+
+
+                //actual main routine to generate data starts here
+                dim = dims[dimIterator];
+                length = lengths[lengthIterator];
+                System.out.println("Generating " + dim + "x" + dim + "x" + dim + " for " + length + " steps");
+
+                i = 0;
+                simArray = FluidSim.initFluidSim(dim,dim,dim);
+                for(i = 0; i < length; i++){
+                    FluidSim.simChunks(simArray,i,TIMESTEP);
+                }
+
+                for(int x = 0; x < dim; x++){
+                    for(int y = 0; y < dim; y++){
+                        for(int z = 0; z < dim; z++){
+                            simArray[x][y][z].dumpToDisk("./src/test/resources/testdata/" + dim + "by" + dim + "/" + length + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + length + "Step.data");
+                        }
+                    }
+                }
+
+                for(int x = 0; x < dim; x++){
+                    for(int y = 0; y < dim; y++){
+                        for(int z = 0; z < dim; z++){
+                            byte[] bytes;
+                            try {
+                                bytes = Files.readAllBytes(new File("./src/test/resources/testdata/" + dim + "by" + dim + "/" + length + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + length + "Step.data").toPath());
+                                ByteBuffer densityBytes = simArray[0][0][0].getDensityBuffer();
+                                i = 0;
+                                while(densityBytes.hasRemaining()){
+                                    boolean pass = bytes[i] == densityBytes.get();
+                                    if(!pass){
+                                        System.err.println("failed to pass!");
+                                    }
+                                    i++;
+                                }
+                            } catch (IOException e) {
+                                // TODO Auto-generated catch block
+                                e.printStackTrace();
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+
+        System.exit(0);
+    }
+
 }

src/main/java/electrosphere/render/GLFWContext.java

@@ -2,7 +2,6 @@ package electrosphere.render;
 
 import java.nio.IntBuffer;
 
-import org.lwjgl.BufferUtils;
 import org.lwjgl.glfw.GLFW;
 import org.lwjgl.opengl.GL;
 import org.lwjgl.opengl.GL44;
@@ -38,9 +37,6 @@ public class GLFWContext {
         GLFW.glfwMakeContextCurrent(window);
         //Maximize it
         GLFW.glfwMaximizeWindow(window);
-        //grab actual framebuffer 
-        int bufferWidth = 0;
-        int bufferHeight = 0;
 
         try(MemoryStack stack = MemoryStack.stackPush()){
             IntBuffer xBuffer = MemoryUtil.memAllocInt(1);
@@ -48,9 +44,6 @@ public class GLFWContext {
             GLFW.glfwGetFramebufferSize(window, xBuffer, yBuffer);
             MemoryUtil.memFree(xBuffer);
             MemoryUtil.memFree(yBuffer);
-            
-            bufferWidth = xBuffer.get();
-            bufferHeight = yBuffer.get();
         }
 
         //
@@ -82,13 +75,19 @@ public class GLFWContext {
         GLFW.glfwSwapInterval(0);
     }
 
-    public static void redraw(){
+    public static void redraw(Mesh[][][] meshes){
         //clear screen
         GL44.glClearColor(0.5f, 0.5f, 0.5f, 1.0f); 
         GL44.glClear(GL44.GL_COLOR_BUFFER_BIT | GL44.
         GL_DEPTH_BUFFER_BIT);
         //draw here
-        Mesh.draw();
+        for(int x = 0; x < meshes.length; x++){
+            for(int y = 0; y < meshes[0].length; y++){
+                for(int z = 0; z < meshes[0][0].length; z++){
+                    meshes[x][y][z].draw();
+                }
+            }
+        }
         //ending stuff
         GLFW.glfwSwapBuffers(window);
         GLFW.glfwPollEvents();

src/main/java/electrosphere/render/Mesh.java

@@ -1,12 +1,11 @@
 package electrosphere.render;
 
-import java.io.BufferedInputStream;
 import java.io.BufferedReader;
+import java.io.File;
 import java.io.IOException;
-import java.io.InputStream;
-import java.io.InputStreamReader;
 import java.nio.FloatBuffer;
 import java.nio.IntBuffer;
+import java.nio.file.Files;
 import java.util.HashMap;
 import java.util.LinkedList;
 import java.util.List;
@@ -14,33 +13,37 @@ import java.util.Map;
 
 import org.joml.Matrix4f;
 import org.joml.Vector3f;
-import org.lwjgl.BufferUtils;
-import org.lwjgl.opengl.GL33;
+import org.joml.Vector3i;
 import org.lwjgl.opengl.GL44;
 import org.lwjgl.opengl.GL45;
-import org.lwjgl.system.MemoryStack;
 import org.lwjgl.system.MemoryUtil;
 
 import electrosphere.FluidSim;
 
 public class Mesh {
     
-    public static int vaoPtr;
-    static int vertBufferPtr;
-    static int normBufferPtr;
-    static int faceBufferPtr;
+    public int vaoPtr;
+    int vertBufferPtr;
+    int normBufferPtr;
+    int faceBufferPtr;
 
     static int shaderProgram;
 
-    public static int elementCount = 0;
+    public int elementCount = 0;
 
-    static Matrix4f viewMatrix = new Matrix4f();
-    static Matrix4f projectionMatrix = new Matrix4f();
-    static Matrix4f model = new Matrix4f().identity();
+    Matrix4f viewMatrix = new Matrix4f();
+    Matrix4f projectionMatrix = new Matrix4f();
+    Matrix4f model = new Matrix4f().identity();
 
-    static MemoryStack stack = MemoryStack.create(32 * 1000 * 1000);
+    FluidSim sim;
+    Vector3i offset;
 
-    public static void meshInitially(FluidSim sim){
+    public Mesh(FluidSim sim, Vector3i offset){
+        this.sim = sim;
+        this.offset = offset;
+    }
+
+    public void meshInitially(){
         //create and bind vao
         vaoPtr = GL44.glGenVertexArrays();
         GL44.glBindVertexArray(vaoPtr);
@@ -101,7 +104,6 @@ public class Mesh {
             FloatBuffer NormalArrayBufferData;
             if(normalCount > 0){
                 NormalArrayBufferData = MemoryUtil.memAllocFloat(normalCount * 3);
-                float[] temp = new float[3];
                 for(float normalValue : data.normals){
                     NormalArrayBufferData.put(normalValue);
                 }
@@ -115,93 +117,34 @@ public class Mesh {
             ex.printStackTrace();
         }
 
-        String vsSrc = "";
-        ClassLoader classloader = Thread.currentThread().getContextClassLoader();
-        try (BufferedReader is = new BufferedReader(new InputStreamReader(classloader.getResourceAsStream("shader.vs")))){
-            String temp;
-            while((temp = is.readLine())!=null){
-                vsSrc = vsSrc + temp + "\n";
-            }
-        } catch (IOException e) {
-            // TODO Auto-generated catch block
-            e.printStackTrace();
-        }
-
-        String fsSrc = "";
-        try (BufferedReader is = new BufferedReader(new InputStreamReader(classloader.getResourceAsStream("shader.fs")))){
-            String temp;
-            while((temp = is.readLine())!=null){
-                fsSrc = fsSrc + temp + "\n";
-            }
-        } catch (IOException e) {
-            // TODO Auto-generated catch block
-            e.printStackTrace();
-        }
-        //
-        //Create shader
-        //
-        int vertexShader = GL45.glCreateShader(GL45.GL_VERTEX_SHADER);
-        GL45.glShaderSource(vertexShader, vsSrc);
-        //Compiles the source for the vertex shader object
-        GL45.glCompileShader(vertexShader);
-        //The following tests if the vertex shader compiles successfully
-        int success;
-        success = GL45.glGetShaderi(vertexShader, GL45.GL_COMPILE_STATUS);
-        if (success != GL45.GL_TRUE) {
-            System.err.println("Vertex Shader failed to compile!");
-            System.err.println("Source is: ");
-            System.err.println(GL45.glGetShaderSource(vertexShader));
-            System.err.println(GL45.glGetShaderInfoLog(vertexShader));
-        }
-        //Creates and opengl object for a fragment shader and assigns its 'pointer' to the integer fragmentShader
-        int fragmentShader = GL45.glCreateShader(GL45.GL_FRAGMENT_SHADER);
-        //This points the opengl shadder object to its proper source
-        GL45.glShaderSource(fragmentShader, fsSrc);
-        //This compiles the shader object
-        GL45.glCompileShader(fragmentShader);
-        //This tests for the success of the compile attempt
-        success = GL45.glGetShaderi(fragmentShader, GL45.GL_COMPILE_STATUS);
-        if (success != GL45.GL_TRUE) {
-            System.err.println("Fragment Shader failed to compile!");
-            System.err.println("Source is: ");
-            System.err.println(GL45.glGetShaderSource(fragmentShader));
-            System.err.println(GL45.glGetShaderInfoLog(fragmentShader));
-        }
-        //This creates a shader program opengl object and assigns its 'pointer' to the integer shaderProgram
-        shaderProgram = GL45.glCreateProgram();
-        //This attaches the vertex and fragment shaders to the program
-        GL45.glAttachShader(shaderProgram, vertexShader);
-        GL45.glAttachShader(shaderProgram, fragmentShader);
-        //This links the program to the GPU (I think its to the GPU anyway)
-        GL45.glLinkProgram(shaderProgram);
-        //Tests for the success of the shader program creation
-        success = GL45.glGetProgrami(shaderProgram, GL45.GL_LINK_STATUS);
-        if (success != GL45.GL_TRUE) {
-            throw new RuntimeException(GL45.glGetProgramInfoLog(shaderProgram));
-        }        
         
-        //Deletes the individual shader objects to free up memory
-        GL45.glDeleteShader(vertexShader);
-        GL45.glDeleteShader(fragmentShader);
-
-        //bind shader
-        GL45.glUseProgram(shaderProgram);
 
         viewMatrix = new Matrix4f();
-        viewMatrix.setLookAt(new Vector3f(-8,20,0), new Vector3f(8,8,8), new Vector3f(0,1,0));
+        viewMatrix.setLookAt(new Vector3f(0,20,16), new Vector3f(8,8,16), new Vector3f(0,1,0));
         projectionMatrix = new Matrix4f();
         projectionMatrix.setPerspective(90,1920.0f/1080.0f,0.0001f,1000f);
-        model = new Matrix4f().identity();
+        model = new Matrix4f().identity().translate(offset.x, offset.y, offset.z);
 
         GL45.glUniformMatrix4fv(GL45.glGetUniformLocation(shaderProgram, "view"), false, viewMatrix.get(new float[16]));
         GL45.glUniformMatrix4fv(GL45.glGetUniformLocation(shaderProgram, "projection"), false, projectionMatrix.get(new float[16]));
         GL45.glUniformMatrix4fv(GL45.glGetUniformLocation(shaderProgram, "model"), false, model.get(new float[16]));
     }
 
-    public static void remesh(FluidSim sim){
+
+    public void remesh(){
         //generate verts
         TerrainChunkData data = generateTerrainChunkData(sim.getData());
 
+        GL44.glBindVertexArray(vaoPtr);
+
+        // for(int i = 0; i < data.normals.size(); i+=3){
+        //     // if(i > 100 && data.normals.get(i+1) < 0.7){
+        //         data.normals.set(i+1, 1f);
+        //         data.normals.set(i+2, 0f);
+        //         System.out.printf("%.2f %.2f %.2f,\n",data.normals.get(i+0),data.normals.get(i+1),data.normals.get(i+2));
+        //     // }
+        // }
+
 
         //
         //Buffer data to GPU
@@ -267,9 +210,9 @@ public class Mesh {
 
     static float angle = 0;
 
-    public static void draw(){
+    public void draw(){
         GL45.glUseProgram(shaderProgram);
-        GL45.glBindFramebuffer(GL45.GL_FRAMEBUFFER, 0);
+        // GL45.glBindFramebuffer(GL45.GL_FRAMEBUFFER, 0);
         GL44.glBindVertexArray(vaoPtr);
         // angle = angle + 0.01f;
         // viewMatrix.identity().setLookAt(
@@ -282,7 +225,85 @@ public class Mesh {
         GL45.glUniformMatrix4fv(GL45.glGetUniformLocation(shaderProgram, "view"), false, viewMatrix.get(new float[16]));
         GL45.glUniformMatrix4fv(GL45.glGetUniformLocation(shaderProgram, "projection"), false, projectionMatrix.get(new float[16]));
         GL45.glUniformMatrix4fv(GL45.glGetUniformLocation(shaderProgram, "model"), false, model.get(new float[16]));
-        GL44.glDrawElements(GL44.GL_TRIANGLES, Mesh.elementCount, GL44.GL_UNSIGNED_INT, 0);
+        GL44.glDrawElements(GL44.GL_TRIANGLES, elementCount, GL44.GL_UNSIGNED_INT, 0);
+    }
+
+
+
+
+
+
+    public static void initShaderProgram(){
+        String vsSrc = "";
+        try (BufferedReader is = new BufferedReader(Files.newBufferedReader(new File("./src/main/resources/shader.vs").toPath()))){
+            String temp;
+            while((temp = is.readLine())!=null){
+                vsSrc = vsSrc + temp + "\n";
+            }
+        } catch (IOException e) {
+            // TODO Auto-generated catch block
+            e.printStackTrace();
+        }
+
+        String fsSrc = "";
+        try (BufferedReader is = new BufferedReader(Files.newBufferedReader(new File("./src/main/resources/shader.fs").toPath()))){
+            String temp;
+            while((temp = is.readLine())!=null){
+                fsSrc = fsSrc + temp + "\n";
+            }
+        } catch (IOException e) {
+            // TODO Auto-generated catch block
+            e.printStackTrace();
+        }
+        //
+        //Create shader
+        //
+        int vertexShader = GL45.glCreateShader(GL45.GL_VERTEX_SHADER);
+        GL45.glShaderSource(vertexShader, vsSrc);
+        //Compiles the source for the vertex shader object
+        GL45.glCompileShader(vertexShader);
+        //The following tests if the vertex shader compiles successfully
+        int success;
+        success = GL45.glGetShaderi(vertexShader, GL45.GL_COMPILE_STATUS);
+        if (success != GL45.GL_TRUE) {
+            System.err.println("Vertex Shader failed to compile!");
+            System.err.println("Source is: ");
+            System.err.println(GL45.glGetShaderSource(vertexShader));
+            System.err.println(GL45.glGetShaderInfoLog(vertexShader));
+        }
+        //Creates and opengl object for a fragment shader and assigns its 'pointer' to the integer fragmentShader
+        int fragmentShader = GL45.glCreateShader(GL45.GL_FRAGMENT_SHADER);
+        //This points the opengl shadder object to its proper source
+        GL45.glShaderSource(fragmentShader, fsSrc);
+        //This compiles the shader object
+        GL45.glCompileShader(fragmentShader);
+        //This tests for the success of the compile attempt
+        success = GL45.glGetShaderi(fragmentShader, GL45.GL_COMPILE_STATUS);
+        if (success != GL45.GL_TRUE) {
+            System.err.println("Fragment Shader failed to compile!");
+            System.err.println("Source is: ");
+            System.err.println(GL45.glGetShaderSource(fragmentShader));
+            System.err.println(GL45.glGetShaderInfoLog(fragmentShader));
+        }
+        //This creates a shader program opengl object and assigns its 'pointer' to the integer shaderProgram
+        shaderProgram = GL45.glCreateProgram();
+        //This attaches the vertex and fragment shaders to the program
+        GL45.glAttachShader(shaderProgram, vertexShader);
+        GL45.glAttachShader(shaderProgram, fragmentShader);
+        //This links the program to the GPU (I think its to the GPU anyway)
+        GL45.glLinkProgram(shaderProgram);
+        //Tests for the success of the shader program creation
+        success = GL45.glGetProgrami(shaderProgram, GL45.GL_LINK_STATUS);
+        if (success != GL45.GL_TRUE) {
+            throw new RuntimeException(GL45.glGetProgramInfoLog(shaderProgram));
+        }        
+        
+        //Deletes the individual shader objects to free up memory
+        GL45.glDeleteShader(vertexShader);
+        GL45.glDeleteShader(fragmentShader);
+
+        //bind shader
+        GL45.glUseProgram(shaderProgram);
     }
 
 
@@ -612,6 +633,37 @@ public class Mesh {
         (byte)0x13
     };
 
+    //allows reverse lookup where vertList index returns the points that vert is concerned with
+    static final int[][] vertIndexInverseTable = new int[][]{
+        {0,1},
+        {1,2},
+        {2,3},
+        {3,0},
+        {4,5},
+        {5,6},
+        {6,7},
+        {7,4},
+        {0,4},
+        {1,5},
+        {2,6},
+        {3,7},
+    };
+
+    //Lookup for the direction vectors for every vertex pair
+    static final Vector3f[] indexVectorLookupTable = new Vector3f[]{
+        new Vector3f(0,0,1),
+        new Vector3f(1,0,0),
+        new Vector3f(0,0,-1),
+        new Vector3f(-1,0,0),
+        new Vector3f(0,0,1),
+        new Vector3f(1,0,0),
+        new Vector3f(0,0,-1),
+        new Vector3f(-1,0,0),
+        new Vector3f(0,1,0),
+        new Vector3f(0,1,0),
+        new Vector3f(0,1,0),
+        new Vector3f(0,1,0),
+    };
 
 
     protected static int polygonize(
@@ -721,34 +773,76 @@ public class Mesh {
             Vector3f v = verts.get(index2).sub(verts.get(index1), new Vector3f());
             Vector3f n = new Vector3f(u.y * v.z - u.z * v.y, u.z * v.x - u.x * v.z, u.x * v.y - u.y * v.x).normalize();
 
+            // Vector3f alignmentVector = new Vector3f();
 
+            // for(int j = 0; j < 3; j++){
+            //     int[] points = vertIndexInverseTable[triTable[cubeIndex][i+j]];
+            //     float factor = 1.0f / Math.abs(points[1] - points[0]);
+            //     //if point two is greater than point one, flip vector
+            //     if(points[0] < points[1]){
+            //         alignmentVector.add(
+            //             -indexVectorLookupTable[triTable[cubeIndex][i+j]].x * factor,
+            //             -indexVectorLookupTable[triTable[cubeIndex][i+j]].y * factor,
+            //             -indexVectorLookupTable[triTable[cubeIndex][i+j]].z * factor
+            //         );
+            //     } else {
+            //         alignmentVector.add(
+            //             indexVectorLookupTable[triTable[cubeIndex][i+j]].z * factor,
+            //             indexVectorLookupTable[triTable[cubeIndex][i+j]].y * factor,
+            //             indexVectorLookupTable[triTable[cubeIndex][i+j]].z * factor
+            //         );
+            //     }
+            // }
+            // n.x = 0;
+            // n.y = 1;
+            // // n.z = 0;
+            // n.normalize();
+            // // n.set(alignmentVector.normalize().mul(-1));
+            // // if(n.dot(alignmentVector.normalize()) >= 0){
+            // //     n.add(alignmentVector.normalize()).div(2.0f);
+            // // }
+
+
+            float oldProportion;
+            float newProportion;
+            Vector3f currentNormal;
 
             //for each vertex, average the new normal with the normals that are already there
             int trianglesSharingIndex0 = trianglesSharingVert.get(index0);
-            //calculate proportion of each normal
-            float oldProportion = trianglesSharingIndex0 / (float)(trianglesSharingIndex0 + 1);
-            float newProportion = 1.0f / (float)(trianglesSharingIndex0 + 1);
-            //increment number of triangles sharing vert
-            trianglesSharingVert.set(index0, trianglesSharingIndex0 + 1);
+            if(trianglesSharingIndex0 > 0){
+                //calculate proportion of each normal
+                oldProportion = trianglesSharingIndex0 / (float)(trianglesSharingIndex0 + 1);
+                newProportion = 1.0f / (float)(trianglesSharingIndex0 + 1);
+                //increment number of triangles sharing vert
+                trianglesSharingVert.set(index0, trianglesSharingIndex0 + 1);
 
-            Vector3f currentNormal = normals.get(index0);
-            currentNormal = averageNormals(currentNormal,oldProportion,n,newProportion);
-            normals.get(index0).set(currentNormal);
+                currentNormal = normals.get(index0);
+                currentNormal = averageNormals(currentNormal,oldProportion,n,newProportion);
+                normals.get(index0).set(currentNormal);
+            } else {
+                trianglesSharingVert.set(index0, trianglesSharingIndex0 + 1);
+                normals.get(index0).set(n);
+            }
 
 
 
 
 
             int trianglesSharingIndex1 = trianglesSharingVert.get(index1);
-            //calculate proportion of each normal
-            oldProportion = trianglesSharingIndex1 / (float)(trianglesSharingIndex1 + 1);
-            newProportion = 1.0f / (float)(trianglesSharingIndex1 + 1);
-            //increment number of triangles sharing vert
-            trianglesSharingVert.set(index1, trianglesSharingIndex1 + 1);
+            if(trianglesSharingIndex1 > 0){
+                //calculate proportion of each normal
+                oldProportion = trianglesSharingIndex1 / (float)(trianglesSharingIndex1 + 1);
+                newProportion = 1.0f / (float)(trianglesSharingIndex1 + 1);
+                //increment number of triangles sharing vert
+                trianglesSharingVert.set(index1, trianglesSharingIndex1 + 1);
 
-            currentNormal = normals.get(index1);
-            currentNormal = averageNormals(currentNormal,oldProportion,n,newProportion);
-            normals.get(index1).set(currentNormal);
+                currentNormal = normals.get(index1);
+                currentNormal = averageNormals(currentNormal,oldProportion,n,newProportion);
+                normals.get(index1).set(currentNormal);
+            } else {
+                trianglesSharingVert.set(index1, trianglesSharingIndex1 + 1);
+                normals.get(index1).set(n);
+            }
 
 
 
@@ -758,15 +852,20 @@ public class Mesh {
 
 
             int trianglesSharingIndex2 = trianglesSharingVert.get(index2);
-            //calculate proportion of each normal
-            oldProportion = trianglesSharingIndex2 / (float)(trianglesSharingIndex2 + 1);
-            newProportion = 1.0f / (float)(trianglesSharingIndex2 + 1);
-            //increment number of triangles sharing vert
-            trianglesSharingVert.set(index2, trianglesSharingIndex2 + 1);
+            if(trianglesSharingIndex2 > 0){
+                //calculate proportion of each normal
+                oldProportion = trianglesSharingIndex2 / (float)(trianglesSharingIndex2 + 1);
+                newProportion = 1.0f / (float)(trianglesSharingIndex2 + 1);
+                //increment number of triangles sharing vert
+                trianglesSharingVert.set(index2, trianglesSharingIndex2 + 1);
 
-            currentNormal = normals.get(index2);
-            currentNormal = averageNormals(currentNormal,oldProportion,n,newProportion);
-            normals.get(index2).set(currentNormal);
+                currentNormal = normals.get(index2);
+                currentNormal = averageNormals(currentNormal,oldProportion,n,newProportion);
+                normals.get(index2).set(currentNormal);
+            } else {
+                trianglesSharingVert.set(index2, trianglesSharingIndex2 + 1);
+                normals.get(index2).set(n);
+            }
 
         }
 
@@ -821,9 +920,9 @@ public class Mesh {
  
  
  
-        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++){
+        for(int x = 1; x < FluidSim.DIM - 1; x++){
+            for(int y = 1; y < FluidSim.DIM - 1; y++){
+                for(int z = 1; 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),
@@ -889,7 +988,16 @@ public class Mesh {
 
     private static Vector3f averageNormals(Vector3f normal0, float proportion0, Vector3f normal1, float proportion1){
         Vector3f rVal = new Vector3f(normal0);
-        rVal = rVal.mul(proportion0).add(new Vector3f(normal1).mul(proportion1));
+        rVal = rVal.mul(proportion0).add(new Vector3f(normal1).mul(proportion1)).normalize();
+        if(Float.isNaN(rVal.x)){
+            rVal.x = 0;
+        }
+        if(Float.isNaN(rVal.y)){
+            rVal.y = 0;
+        }
+        if(Float.isNaN(rVal.z)){
+            rVal.z = 0;
+        }
         return rVal;
     }
 

src/main/resources/shader.fs

@@ -10,6 +10,7 @@ vec3 dLDiffuse = vec3(0.5,0.5,0.5);
 
 in vec3 FragPos;
 in vec3 Normal;
+in vec3 rawPos;
 
 
 uniform vec3 viewPos;
@@ -26,7 +27,15 @@ void main(){
     
     vec3 lightAmount = CalcDirLight(norm, viewDir);
 
-    vec3 color = vec3(0.3,0.7,0.9) * lightAmount;
+    vec3 color = vec3(0.3,0.7,0.9);
+    if(
+        rawPos.x < 2 || rawPos.x > 16 ||
+        rawPos.y < 2 || rawPos.y > 16 ||
+        rawPos.z < 2 || rawPos.z > 16
+    ){
+        color = vec3(0.9,0.7,0.3);
+    }
+    color = color * lightAmount;
 
     //this final calculation is for transparency
     FragColor = vec4(color,1.0);

src/main/resources/shader.vs

@@ -16,6 +16,7 @@ uniform mat4 projection;
 //output buffers
 out vec3 Normal;
 out vec3 FragPos;
+out vec3 rawPos;
 
 
 
@@ -29,6 +30,7 @@ void main() {
     //push frag, normal, and texture positions to fragment shader
     FragPos = vec3(model * FinalVertex);
     Normal = mat3(transpose(inverse(model))) * aNormal;
+    rawPos = aPos;
 
     
     //set final position with opengl space

src/test/java/LongRunTests.java

@@ -0,0 +1,185 @@
+import java.io.IOException;
+import java.io.InputStream;
+import java.nio.ByteBuffer;
+import java.nio.ByteOrder;
+
+import org.junit.Test;
+
+import electrosphere.FluidSim;
+import electrosphere.Main;
+
+public class LongRunTests {
+
+    // @Test
+    // public void test5by5Chunk1Step(){
+
+    //     int dim = 5;
+    //     int maxTimestep = 1;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer fromDiskBuffer = ByteBuffer.allocate(FluidSim.DIM * FluidSim.DIM * FluidSim.DIM * 4);
+    //                     fromDiskBuffer.order(ByteOrder.LITTLE_ENDIAN);
+    //                     fromDiskBuffer.put(bytes);
+    //                     fromDiskBuffer.flip();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = fromDiskBuffer.get() == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+    // @Test
+    // public void test5by5Chunk50Step(){
+
+    //     int dim = 5;
+    //     int maxTimestep = 50;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = bytes[i] == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+    // @Test
+    // public void test5by5Chunk100Step(){
+
+    //     int dim = 5;
+    //     int maxTimestep = 100;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = bytes[i] == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+    // @Test
+    // public void test5by5Chunk500Step(){
+
+    //     int dim = 5;
+    //     int maxTimestep = 500;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = bytes[i] == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+    
+}

src/test/java/MediumRunTests.java

@@ -0,0 +1,197 @@
+import java.io.IOException;
+import java.io.InputStream;
+import java.nio.ByteBuffer;
+import java.nio.ByteOrder;
+
+import org.junit.Test;
+
+import electrosphere.FluidSim;
+import electrosphere.Main;
+
+public class MediumRunTests {
+
+    // @Test
+    // public void test3by3Chunk1Step(){
+
+    //     int dim = 3;
+    //     int maxTimestep = 1;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer fromDiskBuffer = ByteBuffer.allocate(FluidSim.DIM * FluidSim.DIM * FluidSim.DIM * 4);
+    //                     fromDiskBuffer.order(ByteOrder.LITTLE_ENDIAN);
+    //                     fromDiskBuffer.put(bytes);
+    //                     fromDiskBuffer.flip();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = fromDiskBuffer.get() == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+    // @Test
+    // public void test3by3Chunk50Step(){
+
+    //     int dim = 3;
+    //     int maxTimestep = 50;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer fromDiskBuffer = ByteBuffer.allocate(FluidSim.DIM * FluidSim.DIM * FluidSim.DIM * 4);
+    //                     fromDiskBuffer.order(ByteOrder.LITTLE_ENDIAN);
+    //                     fromDiskBuffer.put(bytes);
+    //                     fromDiskBuffer.flip();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = fromDiskBuffer.get() == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+    // @Test
+    // public void test3by3Chunk100Step(){
+
+    //     int dim = 3;
+    //     int maxTimestep = 100;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer fromDiskBuffer = ByteBuffer.allocate(FluidSim.DIM * FluidSim.DIM * FluidSim.DIM * 4);
+    //                     fromDiskBuffer.order(ByteOrder.LITTLE_ENDIAN);
+    //                     fromDiskBuffer.put(bytes);
+    //                     fromDiskBuffer.flip();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = fromDiskBuffer.get() == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+    // @Test
+    // public void test3by3Chunk500Step(){
+
+    //     int dim = 3;
+    //     int maxTimestep = 500;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer fromDiskBuffer = ByteBuffer.allocate(FluidSim.DIM * FluidSim.DIM * FluidSim.DIM * 4);
+    //                     fromDiskBuffer.order(ByteOrder.LITTLE_ENDIAN);
+    //                     fromDiskBuffer.put(bytes);
+    //                     fromDiskBuffer.flip();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = fromDiskBuffer.get() == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+    
+}

src/test/java/ShortRunTest.java

@@ -0,0 +1,202 @@
+
+
+import java.io.IOException;
+import java.io.InputStream;
+import java.nio.ByteBuffer;
+import java.nio.ByteOrder;
+
+import org.junit.Test;
+
+import electrosphere.FluidSim;
+import electrosphere.Main;
+
+/**
+ * Tests stepping a single chunk by a single frame
+ */
+public class ShortRunTest {
+
+    // @Test
+    // public void test1by1Chunk1Step(){
+
+    //     int dim = 1;
+    //     int maxTimestep = 1;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer fromDiskBuffer = ByteBuffer.allocate(FluidSim.DIM * FluidSim.DIM * FluidSim.DIM * 4);
+    //                     fromDiskBuffer.order(ByteOrder.LITTLE_ENDIAN);
+    //                     fromDiskBuffer.put(bytes);
+    //                     fromDiskBuffer.flip();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = fromDiskBuffer.get() == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+    // @Test
+    // public void test1by1Chunk50Step(){
+
+    //     int dim = 1;
+    //     int maxTimestep = 50;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer fromDiskBuffer = ByteBuffer.allocate(FluidSim.DIM * FluidSim.DIM * FluidSim.DIM * 4);
+    //                     fromDiskBuffer.order(ByteOrder.LITTLE_ENDIAN);
+    //                     fromDiskBuffer.put(bytes);
+    //                     fromDiskBuffer.flip();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = fromDiskBuffer.get() == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+    // @Test
+    // public void test1by1Chunk100Step(){
+
+    //     int dim = 1;
+    //     int maxTimestep = 100;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer fromDiskBuffer = ByteBuffer.allocate(FluidSim.DIM * FluidSim.DIM * FluidSim.DIM * 4);
+    //                     fromDiskBuffer.order(ByteOrder.LITTLE_ENDIAN);
+    //                     fromDiskBuffer.put(bytes);
+    //                     fromDiskBuffer.flip();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = fromDiskBuffer.get() == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+    // @Test
+    // public void test1by1Chunk500Step(){
+
+    //     int dim = 1;
+    //     int maxTimestep = 500;
+
+    //     System.out.println("TEST: " + dim + "x" + dim + "x" + dim + " for " + maxTimestep + " steps");
+
+    //     //init chunk array
+    //     FluidSim[][][] simArray = FluidSim.initFluidSim(dim,dim,dim);
+
+    //     //simulate the chunk
+    //     for(int i = 0; i < maxTimestep; i++){
+    //         FluidSim.simChunks(simArray, i, Main.TIMESTEP);
+    //     }
+
+    //     for(int x = 0; x < dim; x++){
+    //         for(int y = 0; y < dim; y++){
+    //             for(int z = 0; z < dim; z++){
+    //                 InputStream testFileIS = this.getClass().getResourceAsStream("./testdata/" + dim + "by" + dim + "/" + maxTimestep + "steps/chunk_" + x + "_" + y + "_" + z + "_" + dim + "by" + dim + "Chunk" + maxTimestep + "Step.data");
+    //                 byte[] bytes;
+    //                 try {
+    //                     bytes = testFileIS.readAllBytes();
+    //                     ByteBuffer fromDiskBuffer = ByteBuffer.allocate(FluidSim.DIM * FluidSim.DIM * FluidSim.DIM * 4);
+    //                     fromDiskBuffer.order(ByteOrder.LITTLE_ENDIAN);
+    //                     fromDiskBuffer.put(bytes);
+    //                     fromDiskBuffer.flip();
+    //                     ByteBuffer densityBytes = simArray[x][y][z].getDensityBuffer();
+    //                     int i = 0;
+    //                     while(densityBytes.hasRemaining()){
+    //                         boolean pass = fromDiskBuffer.get() == densityBytes.get();
+    //                         assert(pass);
+    //                         i++;
+    //                     }
+    //                 } catch (IOException e) {
+    //                     e.printStackTrace();
+    //                     assert(false);
+    //                 }
+    //             }
+    //         }
+    //     }
+
+    //     System.out.println("PASSED");
+
+    // }
+
+}

src/test/resources/shader.fs

@@ -0,0 +1,57 @@
+#version 330 core
+
+out vec4 FragColor;
+
+
+vec3 color = vec3(0.3,0.7,0.9);
+
+vec3 dLDirection = vec3(0.1,-0.9,0.1);
+vec3 dLDiffuse = vec3(0.5,0.5,0.5);
+
+in vec3 FragPos;
+in vec3 Normal;
+in vec3 rawPos;
+
+
+uniform vec3 viewPos;
+// uniform DirLight dirLight;
+// uniform PointLight pointLights[NR_POINT_LIGHTS];
+// uniform SpotLight spotLight;
+
+// function prototypes
+vec3 CalcDirLight(vec3 normal, vec3 viewDir);
+
+void main(){
+    vec3 norm = normalize(Normal);
+    vec3 viewDir = normalize(viewPos - FragPos);
+    
+    vec3 lightAmount = CalcDirLight(norm, viewDir);
+
+    vec3 color = vec3(0.3,0.7,0.9);
+    if(
+        rawPos.x < 2 || rawPos.x > 16 ||
+        rawPos.y < 2 || rawPos.y > 16 ||
+        rawPos.z < 2 || rawPos.z > 16
+    ){
+        color = vec3(0.9,0.7,0.3);
+    }
+    color = color * lightAmount;
+
+    //this final calculation is for transparency
+    FragColor = vec4(color,1.0);
+}
+
+vec3 CalcDirLight(vec3 normal, vec3 viewDir){
+    vec3 lightDir = normalize(-dLDirection);
+    // diffuse shading
+    float diff = max(dot(normal, lightDir), 0.0);
+    // specular shading
+    vec3 reflectDir = reflect(-lightDir, normal);
+    float spec = pow(max(dot(viewDir, reflectDir), 0.0), 0.6);
+    // combine results
+    vec3 diffuse = dLDiffuse * diff;
+
+    vec3 specular = spec * color;
+    
+    return diffuse + specular;
+}

src/test/resources/shader.vs

@@ -0,0 +1,38 @@
+//Vertex Shader
+#version 330 core
+
+//input buffers
+layout (location = 0) in vec3 aPos;
+layout (location = 1) in vec3 aNormal;
+
+
+//coordinate space transformation matrices
+uniform mat4 model;
+uniform mat4 view;
+uniform mat4 projection;
+
+
+
+//output buffers
+out vec3 Normal;
+out vec3 FragPos;
+out vec3 rawPos;
+
+
+
+
+void main() {
+    //normalize posiiton and normal
+    vec4 FinalVertex = vec4(aPos, 1.0);
+    vec4 FinalNormal = vec4(aNormal, 1.0);
+
+
+    //push frag, normal, and texture positions to fragment shader
+    FragPos = vec3(model * FinalVertex);
+    Normal = mat3(transpose(inverse(model))) * aNormal;
+    rawPos = aPos;
+
+    
+    //set final position with opengl space
+    gl_Position = projection * view * model * FinalVertex;
+}