grid2 work

docs/src/progress/renderertodo.md

@@ -1273,6 +1273,13 @@ More cellular determinism verification
 reintroduce sleeping code
 Cellular bounds desync fixes
 
+(12/08/2024)
+cellular stability work
+
+(12/09/2024)
+cellular stability work
+Add grid2 sim files
+
 
 # TODO
 

src/main/c/includes/fluid/dispatch/dispatcher.h

@@ -6,14 +6,19 @@
 #include "fluid/env/environment.h"
 
 
+#define FLUID_DISPATCHER_OVERRIDE_NONE 0
+
+#define FLUID_DISPATCHER_OVERRIDE_CELLULAR 1
+
 
 /**
  * Dispatches chunks to different simulation queues based on the chunk's properties
  * @param numReadIn The number of chunks
  * @param chunkViewC The array of chunks
  * @param environment The environment storing the simulation queues
+ * @param override Overrides the queueing system to force all chunks into a given queue (ie for testing)
  */
-LIBRARY_API void fluid_dispatch(int numReadIn, Chunk ** chunkViewC, Environment * environment);
+LIBRARY_API void fluid_dispatch(int numReadIn, Chunk ** chunkViewC, Environment * environment, int override);
 
 
 

src/main/c/includes/fluid/env/environment.h

@@ -54,6 +54,7 @@ typedef struct {
 typedef struct {
     Chunk ** cellularQueue;
     Chunk ** gridQueue;
+    Chunk ** grid2Queue;
 } FluidSimQueue;
 
 /**

src/main/c/includes/fluid/sim/grid2/mainFunctions.h

@@ -0,0 +1,250 @@
+#include <stdint.h>
+#include "fluid/env/environment.h"
+
+#ifndef FLUID_GRID2_MAINFUNC
+#define FLUID_GRID2_MAINFUNC
+
+void fluid_grid2_simulate(
+    int numChunks,
+    Chunk ** passedInChunks,
+    Environment * environment,
+    jfloat timestep
+);
+
+void fluid_grid2_addSourceToVectors
+  (
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt);
+
+/**
+ * Adds from a source array to a destination array
+*/
+void fluid_grid2_add_source(float * x, float * s, float dt);
+
+
+
+/*
+ * Solves vector diffusion along all axis
+ */
+void fluid_grid2_solveVectorDiffuse (
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+);
+
+
+
+/**
+ * Sets up a projection system of equations
+ * It stores the first derivative of the field in pr, and zeroes out divr.
+ * This allows you to calculate the second derivative into divr using the derivative stored in pr.
+ * @param ur The x velocity grid
+ * @param vr The y velocity grid
+ * @param wr The z velocity grid
+ * @param pr The grid that will contain the first derivative
+ * @param divr The grid that will be zeroed out in preparation of the solver
+ * @param DIFFUSION_CONST The diffusion constant
+ * @param VISCOSITY_CONST The viscosity constant
+ * @param dt The timestep for the simulation
+ */
+void fluid_grid2_setupProjection(
+  int chunk_mask,
+  float ** ur,
+  float ** vr,
+  float ** wr,
+  float ** pr,
+  float ** divr,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+);
+
+
+
+/**
+ * Solves a projection system of equations.
+ * This performs a single iteration across a the p grid to approximate the gradient field.
+ * @param jru0 The gradient field
+ * @param jrv0 The first derivative field
+ */
+void fluid_grid2_solveProjection(
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+);
+
+
+/**
+ * Finalizes a projection.
+ * This subtracts the difference delta along the approximated gradient field.
+ * Thus we are left with an approximately mass-conserved field. 
+ */
+void fluid_grid2_finalizeProjection(
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+);
+
+
+/*
+ * Advects u, v, and w
+ */
+void fluid_grid2_advectVectors(
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+);
+
+
+/**
+ * Actually performs the advection
+*/
+void fluid_grid2_advect(uint32_t chunk_mask, int b, float ** jrd, float ** jrd0, float * u, float * v, float * w, float dt);
+
+
+/**
+ * Sets the bounds of this cube to those of its neighbor
+*/
+void fluid_grid2_setBoundsToNeighborsRaw(
+  int chunk_mask,
+  int vector_dir,
+  float ** neighborArray
+);
+
+
+/**
+ * This exclusively copies neighbors to make sure zeroing out stuff doesn't break sim
+*/
+void fluid_grid2_copyNeighborsRaw(
+  int chunk_mask,
+  int cx,
+  int vector_dir,
+  float ** neighborArray
+);
+
+
+
+
+
+
+
+/**
+ * Adds density to the density array
+ * @return The change in density within this chunk for this frame
+*/
+void fluid_grid2_addDensity(
+    Environment * environment,
+    int chunk_mask,
+    float ** d,
+    float ** d0,
+    float dt
+);
+
+
+
+
+/*
+ * A single iteration of the jacobi to solve density diffusion
+ */
+void fluid_grid2_solveDiffuseDensity(
+  int chunk_mask,
+  float ** d,
+  float ** d0,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+);
+
+
+
+
+
+/**
+ * Advects the density based on the vectors
+*/
+void fluid_grid2_advectDensity(uint32_t chunk_mask, float ** d, float ** d0, float ** ur, float ** vr, float ** wr, float dt);
+
+
+
+
+
+
+/**
+ * Sums the density of the chunk
+ */
+double fluid_grid2_calculateSum(uint32_t chunk_mask, float ** d);
+
+
+
+
+
+
+/**
+ * Normalizes the density array with a given ratio
+ */
+void fluid_grid2_normalizeDensity(float ** d, float ratio);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#endif

src/main/c/includes/fluid/sim/grid2/simulation.h

@@ -0,0 +1,16 @@
+#ifndef FLUID_GRID2_SIMULATION_H
+#define FLUID_GRID2_SIMULATION_H
+
+#include "fluid/queue/chunk.h"
+#include "fluid/env/environment.h"
+
+/**
+ * Performs the main simulation
+ * @param numChunks The number of chunks
+ * @param passedInChunks The chunks to simulate
+ * @param environment The environment data
+ * @param timestep The timestep to simulate by
+ */
+void fluid_grid2_simulate(int numChunks, Chunk ** passedInChunks, Environment * environment, float timestep);
+
+#endif

src/main/c/includes/fluid/sim/grid2/solver_consts.h

@@ -0,0 +1,7 @@
+#ifndef FLUID_GRID2_SOLVER_CONSTS_H
+#define FLUID_GRID2_SOLVER_CONSTS_H
+
+#define FLUID_GRID2_LINEARSOLVERTIMES 3
+#define FLUID_GRID2_VECTOR_DIFFUSE_TIMES 1
+
+#endif

src/main/c/src/fluid/dispatch/dispatcher.c

@@ -11,8 +11,9 @@
  * @param numReadIn The number of chunks
  * @param chunkViewC The array of chunks
  * @param environment The environment storing the simulation queues
+ * @param override Overrides the queueing system to force all chunks into a given queue (ie for testing)
  */
-LIBRARY_API void fluid_dispatch(int numReadIn, Chunk ** chunkViewC, Environment * environment){
+LIBRARY_API void fluid_dispatch(int numReadIn, Chunk ** chunkViewC, Environment * environment, int override){
     //clear queues
     int countCurrent;
     countCurrent = arrlen(environment->queue.cellularQueue);
@@ -23,12 +24,27 @@ LIBRARY_API void fluid_dispatch(int numReadIn, Chunk ** chunkViewC, Environment
     if(countCurrent > 0){
         arrdeln(environment->queue.gridQueue,0,countCurrent);
     }
+    countCurrent = arrlen(environment->queue.grid2Queue);
+    if(countCurrent > 0){
+        arrdeln(environment->queue.grid2Queue,0,countCurrent);
+    }
+
+
+    if(override == FLUID_DISPATCHER_OVERRIDE_CELLULAR){
+        //queue new chunks
+        for(int i = 0; i < numReadIn; i++){
+            Chunk * currentChunk = chunkViewC[i];
+            //TODO: conditionally add to queues based on some values (ie lod, spatial loc, etc)
+            arrput(environment->queue.cellularQueue,currentChunk);
+        }
+        return;
+    }
 
     //queue new chunks
     for(int i = 0; i < numReadIn; i++){
         Chunk * currentChunk = chunkViewC[i];
         //TODO: conditionally add to queues based on some values (ie lod, spatial loc, etc)
-        arrput(environment->queue.cellularQueue,currentChunk);
+        arrput(environment->queue.grid2Queue,currentChunk);
     }
 }
 

src/main/c/src/fluid/queue/javainterface.c

@@ -59,7 +59,7 @@ JNIEXPORT void JNICALL Java_electrosphere_server_fluid_simulator_FluidAccelerate
     environment->consts.dt = dt;
     int numReadIn = readInChunks(env,chunkList,environment);
     fluid_solve_bounds(numReadIn,chunkViewC,environment);
-    fluid_dispatch(numReadIn,chunkViewC,environment);
+    fluid_dispatch(numReadIn,chunkViewC,environment,FLUID_DISPATCHER_OVERRIDE_NONE);
     fluid_simulate(environment);
     updateMetadata(env,numReadIn,chunkViewC,environment);
     

src/main/c/src/fluid/sim/grid2/densitystep.c

@@ -0,0 +1,269 @@
+#include <stdio.h>
+#include <immintrin.h>
+#include <stdint.h>
+#include <jni.h>
+
+#include "fluid/env/utilities.h"
+#include "fluid/queue/chunkmask.h"
+#include "fluid/env/environment.h"
+#include "fluid/queue/chunk.h"
+
+
+/**
+ * Adds density to the density array
+ * @return The change in density within this chunk for this frame
+*/
+void fluid_grid2_addDensity(
+    Environment * environment,
+    int chunk_mask,
+    float ** d,
+    float ** d0,
+    float dt
+){
+    int i;
+    int size=DIM*DIM*DIM;
+    float * x = GET_ARR_RAW(d,CENTER_LOC);
+    float * s = GET_ARR_RAW(d0,CENTER_LOC);
+    for(i=0; i<size; i++){
+        environment->state.newDensity = environment->state.newDensity + dt * s[i];
+        environment->state.existingDensity = environment->state.existingDensity + x[i];
+        x[i] += dt*s[i];
+        if(x[i] < MIN_FLUID_VALUE){
+            x[i] = MIN_FLUID_VALUE;
+        } else if(x[i] > MAX_FLUID_VALUE){
+            x[i] = MAX_FLUID_VALUE;
+        }
+    }
+}
+
+/*
+ * A single iteration of the jacobi to solve density diffusion
+ */
+void fluid_grid2_solveDiffuseDensity(
+  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*DIM*DIM*DIM;
+    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<DIM-1; k++){
+        for(j=1; j<DIM-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < DIM-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>DIM-1){
+                for(i=i-8; i < DIM-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
+*/
+void fluid_grid2_advectDensity(uint32_t chunk_mask, 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*DIM;
+
+    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<DIM-1; k++){
+        for(j=1; j<DIM-1; j++){
+            for(i=1; i<DIM-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 >= DIM-1){ m += 1; }
+                if(y < 1){ n -= 1; }
+                if(y >= DIM-1){ n += 1; }
+                if(z < 1){ o -= 1; }
+                if(z >= DIM-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 >= DIM - 1){
+                    //cases to consider:
+                    //m = 0, x = 17.01
+                    //m = 2, x = 20
+                    x = DIM-1;
+                    i0=(int)DIM-2;
+                    i1=DIM-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>DIM+0.5f) y=DIM+0.5f;
+                if (z<0.5f) z=0.5f;
+                if (z>DIM+0.5f) z=DIM+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 >= DIM){
+                    i0 = DIM - 1;
+                }
+                // if(i0 < 0){
+                //     i0 = 0;
+                // }
+                if(j0 >= DIM){
+                    j0 = DIM - 1;
+                }
+                // if(j0 < 0){
+                //     j0 = 0;
+                // }
+                if(k0 >= DIM){
+                    k0 = DIM - 1;
+                }
+                // if(k0 < 0){
+                //     k0 = 0;
+                // }
+                if(i1 >= DIM){
+                    i1 = DIM - 1;
+                }
+                // if(i1 < 0){
+                //     i1 = 0;
+                // }
+                if(j1 >= DIM){
+                    j1 = DIM - 1;
+                }
+                // if(j1 < 0){
+                //     j1 = 0;
+                // }
+                if(k1 >= DIM){
+                    k1 = DIM - 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)]
+                );
+            }
+        }
+    }
+}
+
+
+/**
+ * Sums the density of the chunk
+ */
+double fluid_grid2_calculateSum(uint32_t chunk_mask, float ** d){
+    int j;
+    int size=DIM*DIM*DIM;
+    float * x = GET_ARR_RAW(d,CENTER_LOC);
+    double rVal = 0;
+    for(j=0; j<size; j++){
+        rVal = rVal + x[j];
+    }
+    return rVal;
+}
+
+/**
+ * Normalizes the density array with a given ratio
+ */
+void fluid_grid2_normalizeDensity(float ** d, float ratio){
+    int j;
+    int size=DIM*DIM*DIM;
+    float * x = GET_ARR_RAW(d,CENTER_LOC);
+    for(j=0; j<size; j++){
+        float value = x[j];
+        value = value * ratio;
+        x[j] = value;
+    }
+}

src/main/c/src/fluid/sim/grid2/fluidsim.c

@@ -0,0 +1,592 @@
+#include <stdint.h>
+
+//native interfaces
+#include "native/electrosphere_server_fluid_simulator_FluidAcceleratedSimulator.h"
+
+//fluid lib
+#include "fluid/env/utilities.h"
+#include "fluid/queue/chunkmask.h"
+#include "fluid/sim/grid2/mainFunctions.h"
+#include "fluid/queue/chunk.h"
+#include "fluid/sim/grid2/simulation.h"
+#include "fluid/sim/grid2/solver_consts.h"
+
+#ifndef SAVE_STEPS
+#define SAVE_STEPS 0
+#endif
+
+
+#define FLUID_GRID2_REALLY_SMALL_VALUE 0.00001
+
+#define FLUID_GRID2_DIFFUSION_CONSTANT 0.00001
+#define FLUID_GRID2_VISCOSITY_CONSTANT 0.00001
+
+static inline void fluid_grid2_saveStep(float * values, const char * name);
+static inline void fluid_grid2_applyGravity(Chunk * currentChunk, Environment * environment);
+static inline void fluid_grid2_clearArr(float ** d);
+
+void fluid_grid2_simulate(
+    int numChunks,
+    Chunk ** passedInChunks,
+    Environment * environment,
+    jfloat timestep
+){
+    Chunk ** 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];
+        fluid_grid2_applyGravity(currentChunk,environment);
+        fluid_grid2_addSourceToVectors(
+            currentChunk->chunkMask,
+            currentChunk->u,
+            currentChunk->v,
+            currentChunk->w,
+            currentChunk->u0,
+            currentChunk->v0,
+            currentChunk->w0,
+            FLUID_GRID2_DIFFUSION_CONSTANT,
+            FLUID_GRID2_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];
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,3,currentChunk->w);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v0);
+            fluid_grid2_copyNeighborsRaw(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 < FLUID_GRID2_VECTOR_DIFFUSE_TIMES; l++){
+            //solve vector diffusion
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                fluid_grid2_solveVectorDiffuse(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_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];
+                fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u);
+                fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v);
+                fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,3,currentChunk->w);
+                // setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u0);
+                // setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v0);
+                // setBoundsToNeighborsRaw(DIM,chunkMask,3,currentChunk->w0);
+                fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->u);
+                fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->v);
+                fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->w);
+                fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->u0);
+                fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->v0);
+                fluid_grid2_copyNeighborsRaw(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);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,3,currentChunk->w);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v0);
+        }
+        //setup projection
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_setupProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_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];
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->u0);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->v0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->u0);
+            fluid_grid2_copyNeighborsRaw(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 < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                fluid_grid2_solveProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_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];
+                fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->u0);
+                fluid_grid2_copyNeighborsRaw(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];
+            fluid_grid2_finalizeProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_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];
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,3,currentChunk->w);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u0);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v0);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,3,currentChunk->w0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,3,currentChunk->w);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v0);
+            fluid_grid2_copyNeighborsRaw(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];
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,3,currentChunk->w);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v0);
+            fluid_grid2_copyNeighborsRaw(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];
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,3,currentChunk->w);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u0);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v0);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,3,currentChunk->w0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,3,currentChunk->w);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,3,currentChunk->w0);
+        }
+        //advect
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_advectVectors(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
+        }
+        //update neighbor arr
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,3,currentChunk->w);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v);
+            fluid_grid2_copyNeighborsRaw(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];
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,3,currentChunk->w);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u0);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,3,currentChunk->w);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v0);
+        }
+        //setup projection
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_setupProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
+        }
+        //update array for vectors
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u0);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u0);
+            fluid_grid2_copyNeighborsRaw(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 < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                fluid_grid2_solveProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
+            }
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->u0);
+                fluid_grid2_copyNeighborsRaw(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];
+            fluid_grid2_finalizeProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
+        }
+        //set boundaries a final time for u,v,w
+        //...
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,3,currentChunk->w);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,1,currentChunk->u0);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,2,currentChunk->v0);
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,3,currentChunk->w0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,3,currentChunk->w);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,1,currentChunk->u0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,2,currentChunk->v0);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,3,currentChunk->w0);
+        }
+    }
+
+
+
+
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
+
+
+
+    //add density
+    {
+        double deltaDensity = 0;
+        environment->state.existingDensity = 0;
+        environment->state.newDensity = 0;
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_addDensity(environment,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];
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->d);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->d0);
+        }
+    }
+    //diffuse density
+    {
+        for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                fluid_grid2_solveDiffuseDensity(currentChunk->chunkMask,currentChunk->d,currentChunk->d0,currentChunk->u,currentChunk->v,currentChunk->w,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
+            }
+            for(int i = 0; i < numChunks; i++){
+                Chunk * currentChunk = chunks[i];
+                fluid_grid2_setBoundsToNeighborsRaw(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];
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->d);
+            fluid_grid2_copyNeighborsRaw(currentChunk->chunkMask,0,0,currentChunk->d0);
+        }
+    }
+    //advect density
+    {
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_advectDensity(currentChunk->chunkMask,currentChunk->d,currentChunk->d0,currentChunk->u,currentChunk->v,currentChunk->w,timestep);
+        }
+    }
+    //mirror densities
+    {
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->d);
+        }
+    }
+    //normalize densities
+    {
+        double transformedDensity = 0;
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            transformedDensity = transformedDensity + fluid_grid2_calculateSum(currentChunk->chunkMask,currentChunk->d);
+        }
+        float normalizationRatio = 0;
+        if(transformedDensity != 0){
+            normalizationRatio = (environment->state.existingDensity + environment->state.newDensity) / transformedDensity;
+            environment->state.normalizationRatio = normalizationRatio;
+        }
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_normalizeDensity(currentChunk->d,normalizationRatio);
+        }
+    }
+    //clear delta arrays
+    {
+        for(int i = 0; i < numChunks; i++){
+            Chunk * currentChunk = chunks[i];
+            fluid_grid2_clearArr(currentChunk->d0);
+            fluid_grid2_clearArr(currentChunk->u0);
+            fluid_grid2_clearArr(currentChunk->v0);
+            fluid_grid2_clearArr(currentChunk->w0);
+        }
+    }
+}
+
+
+static inline void fluid_grid2_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 < FLUID_GRID2_REALLY_SMALL_VALUE && val > -FLUID_GRID2_REALLY_SMALL_VALUE){
+                        val = 0;
+                    }
+                    fprintf(fp, "%f\t", val);
+                }
+                fprintf(fp, "\n");
+            }
+            fprintf(fp, "\n");
+        }
+
+        fclose(fp);
+    }
+}
+
+/**
+ * Applies gravity to the chunk
+ * @param currentChunk The chunk to apply on
+ * @param environment The environment data of the world
+ */
+static inline void fluid_grid2_applyGravity(Chunk * currentChunk, Environment * environment){
+    int N = DIM;
+    for(int x = 0; x < DIM; x++){
+        for(int y = 0; y < DIM; y++){
+            for(int z = 0; z < DIM; z++){
+                GET_ARR_RAW(currentChunk->v0,CENTER_LOC)[IX(x,y,z)] = GET_ARR_RAW(currentChunk->v0,CENTER_LOC)[IX(x,y,z)] + GET_ARR_RAW(currentChunk->d,CENTER_LOC)[IX(x,y,z)] * environment->consts.gravity;
+            }
+        }
+    }
+}
+
+/**
+ * Clears an array
+ */
+static inline void fluid_grid2_clearArr(float ** d){
+    float * x = GET_ARR_RAW(d,CENTER_LOC);
+    for(int j = 0; j < DIM * DIM * DIM; j++){
+        x[j] = 0;
+    }
+}

src/main/c/src/fluid/sim/grid2/velocitystep.c

@@ -0,0 +1,941 @@
+#include <stdio.h>
+#include <immintrin.h>
+#include <stdint.h>
+
+#include "fluid/env/utilities.h"
+#include "fluid/queue/chunkmask.h"
+#include "fluid/queue/chunk.h"
+#include "fluid/sim/grid2/solver_consts.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
+
+void fluid_grid2_add_source(float * x, float * s, float dt);
+void fluid_grid2_advect(uint32_t chunk_mask, int b, float ** jrd, float ** jrd0, float * u, float * v, float * w, float dt);
+
+
+/*
+ * Adds force to all vectors
+ */
+void fluid_grid2_addSourceToVectors
+  (
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt){
+    fluid_grid2_add_source(GET_ARR_RAW(jru,CENTER_LOC),GET_ARR_RAW(jru0,CENTER_LOC),dt);
+    fluid_grid2_add_source(GET_ARR_RAW(jrv,CENTER_LOC),GET_ARR_RAW(jrv0,CENTER_LOC),dt);
+    fluid_grid2_add_source(GET_ARR_RAW(jrw,CENTER_LOC),GET_ARR_RAW(jrw0,CENTER_LOC),dt);
+}
+
+/**
+ * Adds from a source array to a destination array
+*/
+void fluid_grid2_add_source(float * x, float * s, float dt){
+	int i;
+    int size=DIM*DIM*DIM;
+	for(i=0; i<size; i++){
+        x[i] += dt*s[i];
+    }
+}
+
+/*
+ * Solves vector diffusion along all axis
+ */
+void fluid_grid2_solveVectorDiffuse (
+  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*DIM*DIM*DIM;
+    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<DIM-1; k++){
+        for(j=1; j<DIM-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < DIM-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>DIM-1){
+                for(i=i-8; i < DIM-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<DIM-1; k++){
+        for(j=1; j<DIM-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < DIM-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>DIM-1){
+                for(i=i-8; i < DIM-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<DIM-1; k++){
+        for(j=1; j<DIM-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < DIM-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>DIM-1){
+                for(i=i-8; i < DIM-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
+ * It stores the first derivative of the field in pr, and zeroes out divr.
+ * This allows you to calculate the second derivative into divr using the derivative stored in pr.
+ * @param ur The x velocity grid
+ * @param vr The y velocity grid
+ * @param wr The z velocity grid
+ * @param pr The grid that will contain the first derivative
+ * @param divr The grid that will be zeroed out in preparation of the solver
+ * @param DIFFUSION_CONST The diffusion constant
+ * @param VISCOSITY_CONST The viscosity constant
+ * @param dt The timestep for the simulation
+ */
+void fluid_grid2_setupProjection(
+  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(DIM);
+    __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/DIM;
+
+    for(k=1; k<DIM-1; k++){
+        for(j=1; j<DIM-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.
+ * This performs a single iteration across a the p grid to approximate the gradient field.
+ * @param jru0 The gradient field
+ * @param jrv0 The first derivative field
+ */
+void fluid_grid2_solveProjection(
+  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<DIM-1; k++){
+        for(j=1; j<DIM-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < DIM-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>DIM-1){
+                for(i=i-8; i < DIM-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.
+ * This subtracts the difference delta along the approximated gradient field.
+ * Thus we are left with an approximately mass-conserved field. 
+ */
+void fluid_grid2_finalizeProjection(
+  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*DIM);
+    __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<DIM-1 ; k++ ) {
+        for ( j=1 ; j<DIM-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
+ */
+void fluid_grid2_advectVectors(
+  int chunk_mask,
+  float ** jru,
+  float ** jrv,
+  float ** jrw,
+  float ** jru0,
+  float ** jrv0,
+  float ** jrw0,
+  float DIFFUSION_CONST,
+  float VISCOSITY_CONST,
+  float dt
+){
+    fluid_grid2_advect(chunk_mask,1,jru,jru0,GET_ARR_RAW(jru0,CENTER_LOC),GET_ARR_RAW(jrv0,CENTER_LOC),GET_ARR_RAW(jrw0,CENTER_LOC),dt);
+    fluid_grid2_advect(chunk_mask,2,jrv,jrv0,GET_ARR_RAW(jru0,CENTER_LOC),GET_ARR_RAW(jrv0,CENTER_LOC),GET_ARR_RAW(jrw0,CENTER_LOC),dt);
+    fluid_grid2_advect(chunk_mask,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
+*/
+void fluid_grid2_advect(uint32_t chunk_mask, 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*DIM;
+
+    float * d = GET_ARR_RAW(jrd,CENTER_LOC);
+
+    float * d0 = GET_ARR_RAW(jrd0,CENTER_LOC);
+
+    for(k=1; k<DIM-1; k++){
+        for(j=1; j<DIM-1; j++){
+            for(i=1; i<DIM-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 >= DIM){ m -= 1; }
+                if(y < 0){ n += 1; }
+                else if(y >= DIM){ n -= 1; }
+                if(z < 0){ o += 1; }
+                else if(z >= DIM){ 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)] * (DIM-2);
+                    // printf("%d => %f\n",m,x);
+                    y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * (DIM-2);
+                    z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * (DIM-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 > DIM - 1){
+                    //cases to consider:
+                    //m = 0, x = 17.01
+                    //m = 2, x = 20
+                    x = DIM-1;
+                    i0=(int)DIM-2;
+                    i1=DIM-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 > DIM - 1){
+                    //cases to consider:
+                    //m = 0, x = 17.01
+                    //m = 2, x = 20
+                    y = DIM-1;
+                    j0=(int)DIM-2;
+                    j1=DIM-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 > DIM - 1){
+                    //cases to consider:
+                    //m = 0, x = 17.01
+                    //m = 2, x = 20
+                    z = DIM-1;
+                    k0=(int)DIM-2;
+                    k1=DIM-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 >= DIM){
+                    j0 = DIM - 1;
+                }
+                // if(j0 < 0){
+                //     j0 = 0;
+                // }
+                if(k0 >= DIM){
+                    k0 = DIM - 1;
+                }
+                // if(k0 < 0){
+                //     k0 = 0;
+                // }
+                // if(i1 >= N){
+                //     i1 = N - 1;
+                // }
+                // if(i1 < 0){
+                //     i1 = 0;
+                // }
+                if(j1 >= DIM){
+                    j1 = DIM - 1;
+                }
+                // if(j1 < 0){
+                //     j1 = 0;
+                // }
+                if(k1 >= DIM){
+                    k1 = DIM - 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
+*/
+void fluid_grid2_setBoundsToNeighborsRaw(
+  int chunk_mask,
+  int vector_dir,
+  float ** neighborArray
+){
+    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
+*/
+void fluid_grid2_copyNeighborsRaw(
+  int chunk_mask,
+  int cx,
+  int vector_dir,
+  float ** neighborArray
+){
+    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/c/src/fluid/sim/simulator.c

@@ -5,6 +5,7 @@
 #include "fluid/dispatch/dispatcher.h"
 #include "fluid/sim/simulator.h"
 #include "fluid/sim/grid/simulation.h"
+#include "fluid/sim/grid2/simulation.h"
 #include "fluid/sim/cellular/cellular.h"
 #include "fluid/queue/chunk.h"
 #include "fluid/env/environment.h"
@@ -19,12 +20,27 @@ LIBRARY_API void fluid_simulate(Environment * environment){
     int currentCount, i;
 
     //cellular sim
-    fluid_cellular_simulate(environment);
+    {
+        currentCount = stbds_arrlen(queue.cellularQueue);
+        if(currentCount > 0){
+            fluid_cellular_simulate(environment);
+        }
+    }
 
     //grid sim
     {
         currentCount = stbds_arrlen(queue.gridQueue);
-        fluid_grid_simulate(currentCount,queue.gridQueue,environment,environment->consts.dt);
+        if(currentCount > 0){
+            fluid_grid_simulate(currentCount,queue.gridQueue,environment,environment->consts.dt);
+        }
+    }
+
+    //grid2 sim
+    {
+        currentCount = stbds_arrlen(queue.grid2Queue);
+        if(currentCount > 0){
+            fluid_grid2_simulate(currentCount,queue.grid2Queue,environment,environment->consts.dt);
+        }
     }
 }
 

src/test/c/fluid/dispatch/dispatcher_tests.c

@@ -14,7 +14,7 @@ int fluid_dispatch_dispatcher_tests(){
 
     int queueSize = 10;
     Chunk ** queue = chunk_create_queue(queueSize);
-    fluid_dispatch(queueSize,queue,env);
+    fluid_dispatch(queueSize,queue,env,FLUID_DISPATCHER_OVERRIDE_NONE);
 
     int gridChunksFound = stbds_arrlen(env->queue.gridQueue) +
         stbds_arrlen(env->queue.cellularQueue)

src/test/c/fluid/sim/cellular/cellular_tests.c

@@ -283,7 +283,7 @@ int fluid_sim_cellular_bounds_test1(){
     }
 
     //dispatch and simulate
-    fluid_dispatch(chunkCount,queue,env);
+    fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
     fluid_simulate(env);
 
     //solve bounds afterwards to properly push data back into real arrays
@@ -344,7 +344,7 @@ int fluid_sim_cellular_bounds_test2(){
     }
 
     //dispatch and simulate
-    fluid_dispatch(chunkCount,queue,env);
+    fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
     fluid_simulate(env);
 
     //solve bounds afterwards to properly push data back into real arrays
@@ -408,7 +408,7 @@ int fluid_sim_cellular_stability_test1(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -477,7 +477,7 @@ int fluid_sim_cellular_stability_test2(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -549,7 +549,7 @@ int fluid_sim_cellular_stability_test3(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -620,7 +620,7 @@ int fluid_sim_cellular_stability_test4(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -698,7 +698,7 @@ int fluid_sim_cellular_stability_test5(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -776,7 +776,7 @@ int fluid_sim_cellular_stability_test6(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -847,7 +847,7 @@ int fluid_sim_cellular_stability_test7(){
         float currentSum = chunk_queue_sum(queue);
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         float postSum = chunk_queue_sum(queue);
@@ -934,7 +934,7 @@ int fluid_sim_cellular_stability_test8(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -1013,7 +1013,7 @@ int fluid_sim_cellular_stability_test9(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -1091,7 +1091,7 @@ int fluid_sim_cellular_stability_test10(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -1161,7 +1161,7 @@ int fluid_sim_cellular_stability_test11(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -1235,7 +1235,7 @@ int fluid_sim_cellular_stability_test12(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -1314,7 +1314,7 @@ int fluid_sim_cellular_stability_test13(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         // printf("\n");
@@ -1398,7 +1398,7 @@ int fluid_sim_cellular_stability_test14(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         env->state.frame++;
@@ -1486,7 +1486,7 @@ int fluid_sim_cellular_stability_test15(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         env->state.frame++;
@@ -1569,7 +1569,7 @@ int fluid_sim_cellular_stability_test16(){
         }
         // printf("frame: %d   --- %f   \n", frameCounter,currentSum);
         fluid_solve_bounds(chunkCount,queue,env);
-        fluid_dispatch(chunkCount,queue,env);
+        fluid_dispatch(chunkCount,queue,env,FLUID_DISPATCHER_OVERRIDE_CELLULAR);
         fluid_simulate(env);
         fluid_solve_bounds(chunkCount,queue,env);
         env->state.frame++;

src/test/c/fluid/sim/simulator_tests.c

@@ -20,7 +20,7 @@ int fluid_sim_simulator_tests(int argc, char **argv){
 
     int queueSize = 10;
     Chunk ** queue = chunk_create_queue(queueSize);
-    fluid_dispatch(queueSize,queue,env);
+    fluid_dispatch(queueSize,queue,env,FLUID_DISPATCHER_OVERRIDE_NONE);
 
     fluid_simulate(env);