density not crashing, but broken

src/main/c/compile.sh

@@ -61,12 +61,17 @@ INPUT_FILES="./src/threadpool.c"
 OUTPUT_FILE="./threadpool.o"
 gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
 
+COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
+INPUT_FILES="./src/libfluidsim.c"
+OUTPUT_FILE="./libfluidsim.o"
+gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
+
 
 
 #compile shared object file
 OUTPUT_FILE="libfluidsim$LIB_ENDING"
 COMPILE_FLAGS="-shared"
-INPUT_FILES="densitystep.o velocitystep.o chunkmask.o threadpool.o"
+INPUT_FILES="densitystep.o velocitystep.o chunkmask.o threadpool.o libfluidsim.o"
 gcc $COMPILE_FLAGS $INPUT_FILES -o $OUTPUT_FILE
 
 #move to resources

src/main/c/includes/electrosphere_FluidSim.h

@@ -16,7 +16,7 @@ extern "C" {
 #undef electrosphere_FluidSim_LINEARSOLVERTIMES
 #define electrosphere_FluidSim_LINEARSOLVERTIMES 20L
 #undef electrosphere_FluidSim_GRAVITY
-#define electrosphere_FluidSim_GRAVITY -100.0f
+#define electrosphere_FluidSim_GRAVITY 0.0f
 /*
  * Class:     electrosphere_FluidSim
  * Method:    createThreadpool

src/main/c/includes/libfluidsim.h

@@ -2,5 +2,15 @@
 #ifndef LIB_FLUID_SIM
 #define LIB_FLUID_SIM
 
+#include "./threadpool.h"
+
+void simulate(LibraryContext * context, int threadIndex);
+
+void addDensity(int N, int chunk_mask, float * jrx, float * x0, float dt);
+
+void diffuseDensity(int N, int chunk_mask, float * jrx, float * jrx0, float DIFFUSION_CONST, float VISCOSITY_CONST, float dt);
+
+void advectDensity(uint32_t chunk_mask, int N, float * jrd, float * jrd0, float * u, float * v, float * w, float dt);
+
 //close include guard
 #endif

src/main/c/includes/threadpool.h

@@ -1,6 +1,8 @@
 #ifndef THREADPOOL
 #define THREADPOOL
 
+#include <pthread.h>
+
 
 /**
  * A threadpool used to distribute work
@@ -18,6 +20,7 @@ typedef struct {
 */
 typedef struct {
     int chunkMask;
+    float timestep;
     float * d;
     float * u;
     float * v;
@@ -26,6 +29,8 @@ typedef struct {
     float * u0;
     float * v0;
     float * w0;
+    float diffuseConst;
+    float viscosityConst;
 } Chunk;
 
 

src/main/c/src/densitystep.c

@@ -1,4 +1,3 @@
-#include <jni.h>
 #include <stdio.h>
 #include <immintrin.h>
 #include <stdint.h>
@@ -7,53 +6,35 @@
 #include "../includes/utilities.h"
 #include "../includes/chunkmask.h"
 
-void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, jobjectArray d0, float * u, float * v, float * w, float dt);
 
-/*
+void addDensity(
- * Class:     electrosphere_FluidSim
+  int N,
- * Method:    addDensity
+  int chunk_mask,
- * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;F)V
+  float * x,
- */
+  float * s,
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_addDensity
+  float dt){
-  (JNIEnv * env,
-  jobject this,
-  jint N,
-  jint chunk_mask,
-  jobjectArray jrx,
-  jobjectArray x0,
-  jfloat dt){
     int i;
     int size=N*N*N;
-    float * x = GET_ARR(env,jrx,CENTER_LOC);
-    float * s = GET_ARR(env,x0,CENTER_LOC);
     for(i=0; i<size; i++){
         x[i] += dt*s[i];
     }
 }
 
-/*
+
- * Class:     electrosphere_FluidSim
+void diffuseDensity(
- * Method:    solveDiffuseDensity
+  int N,
- * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
+  int chunk_mask,
- */
+  float * jrx,
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
+  float * jrx0,
-  (JNIEnv * env,
+  float DIFFUSION_CONST,
-  jobject this,
+  float VISCOSITY_CONST,
-  jint N,
+  float dt){
-  jint chunk_mask,
-  jobjectArray jrx,
-  jobjectArray jrx0,
-  jobjectArray jru,
-  jobjectArray jrv,
-  jobjectArray jrw,
-  jfloat DIFFUSION_CONST,
-  jfloat VISCOSITY_CONST,
-  jfloat dt){
     float a=dt*DIFFUSION_CONST*N*N*N;
     float c=1+6*a;
     int i, j, k, l, m;
-    float * x = GET_ARR(env,jrx,CENTER_LOC);
+    float * x = jrx;
-    float * x0 = GET_ARR(env,jrx0,CENTER_LOC);
+    float * x0 = jrx0;
+
     
     __m256 aScalar = _mm256_set1_ps(a);
     __m256 cScalar = _mm256_set1_ps(c);
@@ -85,42 +66,21 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
     }
 }
 
-/*
+void advectDensity(uint32_t chunk_mask, int N, float * jrd, float * jrd0, float * u, float * v, float * w, float dt){
- * Class:     electrosphere_FluidSim
- * Method:    advectDensity
- * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
- */
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectDensity
-  (JNIEnv * env,
-  jobject this,
-  jint N,
-  jint chunk_mask,
-  jobjectArray jrx,
-  jobjectArray jrx0,
-  jobjectArray jru,
-  jobjectArray jrv,
-  jobjectArray jrw,
-  jfloat DIFFUSION_CONST,
-  jfloat VISCOSITY_CONST,
-  jfloat dt){
-    advectDensity(env,chunk_mask,N,3,jrx,jrx0,GET_ARR(env,jru,CENTER_LOC),GET_ARR(env,jrv,CENTER_LOC),GET_ARR(env,jrw,CENTER_LOC),dt);
-}
-
-void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, jobjectArray 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(env,jrd,CENTER_LOC);
+    float * d = jrd;
 
-    float * d0 = GET_ARR(env,jrd0,CENTER_LOC);
+    float * d0 = jrd0;
 
     for(k=1; k<N-1; k++){
         for(j=1; j<N-1; j++){
             for(i=1; i<N-1; i++){
-                d0 = GET_ARR(env,jrd0,CENTER_LOC);
+                d0 = jrd0;
                 //calculate location to pull from
                 x = i-dtx*u[IX(i,j,k)];
                 y = j-dty*v[IX(i,j,k)];

src/main/c/src/libfluidsim.c

@@ -0,0 +1,77 @@
+#include <jni.h>
+#include <stdio.h>
+#include <immintrin.h>
+#include <stdint.h>
+#include <unistd.h>
+#include <pthread.h>
+#include <semaphore.h>
+
+//git libs
+#include "../lib/stb/stb_ds.h"
+
+//all other headers
+#include "../includes/libfluidsim.h"
+#include "../includes/threadpool.h"
+#include "../includes/chunkmask.h"
+#include "../includes/utilities.h"
+
+#define DIM 18
+#define LINEARSOLVERTIMES 20
+
+
+//call this if the child threads need to wait on one another
+void sync(LibraryContext * context){
+    pthread_barrier_wait(context->threadpool->barrierMain);
+}
+
+/**
+ * Main simulation function
+*/
+void simulate(LibraryContext * context, int threadIndex){
+    int numThreads = context->threadpool->numThreads;
+    int numChunks = context->numChunks;
+
+    //add densities to chunks
+    for(int i = threadIndex; i < numChunks; i = i + numThreads){
+        Chunk chunk = context->chunks[i];
+        addDensity(DIM,chunk.chunkMask,chunk.d,chunk.d0,chunk.timestep);
+    }
+    sync(context);
+
+    //diffuse density
+    for(int l = 0; l < LINEARSOLVERTIMES; l++){
+        for(int i = threadIndex; i < numChunks; i = i + numThreads){
+            Chunk chunk = context->chunks[i];
+            SWAP(chunk.d,chunk.d0);
+            diffuseDensity(
+                DIM,
+                chunk.chunkMask,
+                chunk.d,
+                chunk.d0,
+                chunk.diffuseConst,
+                chunk.viscosityConst,
+                chunk.timestep
+            );
+            SWAP(chunk.d,chunk.d0);
+        }
+        sync(context);
+    }
+
+    //advect density
+    for(int i = threadIndex; i < numChunks; i = i + numThreads){
+        Chunk chunk = context->chunks[i];
+        advectDensity(
+            chunk.chunkMask,
+            DIM,
+            chunk.d,
+            chunk.d0,
+            chunk.u,
+            chunk.v,
+            chunk.w,
+            chunk.timestep
+        );
+    }
+
+    //call this if the child threads need to wait on one another
+    sync(context);
+}

src/main/c/src/threadpool.c

@@ -17,7 +17,7 @@
 #include "../includes/utilities.h"
 
 typedef struct {
-  ThreadPool * threadPool;
+  LibraryContext * libraryContext;
   int threadIndex; //the index for this thread in all threadpool objects
 } ThreadData;
 
@@ -39,6 +39,7 @@ JNIEXPORT jlong JNICALL Java_electrosphere_FluidSim_createThreadpool
 
   LibraryContext * libraryContext = (LibraryContext *)malloc(sizeof(LibraryContext));
   libraryContext->numChunks = 0;
+  libraryContext->chunks = NULL; //must be tied to NULL at the start
 
   //init threadpool
   libraryContext->threadpool = (ThreadPool *)malloc(sizeof(ThreadPool));
@@ -71,7 +72,7 @@ JNIEXPORT jlong JNICALL Java_electrosphere_FluidSim_createThreadpool
 
     //create data to be sent to thread
     ThreadData * threadData = (ThreadData *)malloc(sizeof(ThreadData));
-    threadData->threadPool = pool;
+    threadData->libraryContext = libraryContext;
     threadData->threadIndex = i;
 
     //create thread
@@ -124,7 +125,6 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_queueChunk(
   //if the array is full, malloc a new chunk, otherwise overwrite an existing one
   int arrayCapacity = stbds_arrlen(libraryContext->chunks);
   if(libraryContext->numChunks >= arrayCapacity){
-    printf("%d > %d\n",libraryContext->numChunks,arrayCapacity);
     Chunk * newChunk = (Chunk *)malloc(sizeof(Chunk));
     newChunk->chunkMask = chunkMask;
     newChunk->u = GET_ARR(env,ur,CENTER_LOC);
@@ -134,8 +134,11 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_queueChunk(
     newChunk->v0 = GET_ARR(env,v0r,CENTER_LOC);
     newChunk->w0 = GET_ARR(env,w0r,CENTER_LOC);
     newChunk->d = GET_ARR(env,dr,CENTER_LOC);
-    stbds_arrput(libraryContext->chunks,newChunk[0]);
     newChunk->d0 = GET_ARR(env,d0r,CENTER_LOC);
+    newChunk->timestep = timestep;
+    newChunk->diffuseConst = DIFFUSION_CONSTANT;
+    newChunk->viscosityConst = VISCOSITY_CONSTANT;
+    stbds_arrput(libraryContext->chunks,newChunk[0]);
   } else {
     Chunk * currentChunk = &libraryContext->chunks[libraryContext->numChunks];
     currentChunk->chunkMask = chunkMask;
@@ -147,10 +150,12 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_queueChunk(
     currentChunk->w0 = GET_ARR(env,w0r,CENTER_LOC);
     currentChunk->d = GET_ARR(env,dr,CENTER_LOC);
     currentChunk->d0 = GET_ARR(env,d0r,CENTER_LOC);
+    currentChunk->timestep = timestep;
+    currentChunk->diffuseConst = DIFFUSION_CONSTANT;
+    currentChunk->viscosityConst = VISCOSITY_CONSTANT;
   }
 
   libraryContext->numChunks++;
-  // arrput(libraryContext->chunks,newChunk);
 }
 
 /*
@@ -191,30 +196,26 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_fetchWork(
 */
 void * mainThreadLoop(void * dataRaw){
   ThreadData * threadData = (ThreadData *)dataRaw;
-  ThreadPool * threadPool = threadData->threadPool;
+  LibraryContext * libraryContext = threadData->libraryContext;
   int threadIndex = threadData->threadIndex;
 
-  printf("thread work\n");
 
   //main thread loop
   int running = 1;
-  while(running){
+  while(running==1){
 
     //begin work
     //wait to begin work until parent signals its ready
-    pthread_barrier_wait(threadPool->barrierWithParentThread);
+    pthread_barrier_wait(libraryContext->threadpool->barrierWithParentThread);
 
     //do main work
-
+    simulate(libraryContext,threadData->threadIndex);
-
-    //call this if the child threads need to wait on one another
-    pthread_barrier_wait(threadPool->barrierMain);
 
 
 
     //finalize work
     //the parent thread needs to call the barrier as well with whatever method it uses to grab data from the sim back to java
-    pthread_barrier_wait(threadPool->barrierWithParentThread);
+    pthread_barrier_wait(libraryContext->threadpool->barrierWithParentThread);
 
 
   }

src/main/java/electrosphere/FluidSim.java

@@ -79,7 +79,8 @@ public class FluidSim {
 
     static final int LINEARSOLVERTIMES = 20;
 
-    static final float GRAVITY = -100f;
+    static final float GRAVITY = 0.0f;
+    // static final float GRAVITY = -100f; //<-- original value
     
     public void setup(Vector3i offset){
         //allocate buffers for this chunk
@@ -713,6 +714,162 @@ public class FluidSim {
 
 
 
+
+    /**
+     *  OLD NATIVE DECLARATIONS
+     */
+
+
+
+    // /**
+    //  * 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 chunkMask,
+    //     ByteBuffer[] x,
+    //     ByteBuffer x0,
+    //     ByteBuffer[] u,
+    //     ByteBuffer[] v,
+    //     ByteBuffer[] w,
+    //     ByteBuffer[] u0,
+    //     ByteBuffer[] v0,
+    //     ByteBuffer[] w0,
+    //     float DIFFUSION_CONSTANT,
+    //     float VISCOSITY_CONSTANT,
+    //     float timestep
+    // );
+
+    // private void calculateChunkMaskWrapper(){
+    //     this.chunkMask = this.calculateChunkMask(density);
+    // }
+    // /**
+    //  * Calculates the mask of chunk neighbors
+    //  * @param densityBuffers The neighbor array
+    //  * @return The mask
+    //  */
+    // private native int calculateChunkMask(ByteBuffer[] densityBuffers);
+
+    // /**
+    //  * Add vector values to u, v, and w all at once
+    //  */
+    // private void addSourceToVectorsWrapper(float timestep){
+    //     addSourceToVectors(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
+    // }
+    // private native void addSourceToVectors(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+
+    // /**
+    //  * Solves u, v, and w diffusion systems of equations
+    //  */
+    // private void solveVectorDiffuseWrapper(float timestep){
+    //     solveVectorDiffuse(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
+    // }
+    // private native void solveVectorDiffuse(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+
+    // /**
+    //  * Setup projection system
+    //  */
+    // private void setupProjectionWrapper(float timestep){
+    //     setupProjection(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
+    // }
+    // private native void setupProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+
+    // /**
+    //  * Solve projection system
+    //  */
+    // private void solveProjectionWrapper(float timestep){
+    //     solveProjection(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
+    // }
+    // private native void solveProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+
+    // /**
+    //  * Does work like subtracting curl from vector field, setting boundaries, etc
+    //  */
+    // private void finalizeProjectionWrapper(float timestep){
+    //     finalizeProjection(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
+    // }
+    // private native void finalizeProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+
+    // /**
+    //  * Advects vectors
+    //  */
+    // private void advectVectorsWrapper(float timestep){
+    //     advectVectors(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
+    // }
+    // private native void advectVectors(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+
+    // /**
+    //  * Adds density to the simulation
+    //  */
+    // private void addDensityWrapper(float timestep){
+    //     addDensity(DIM, chunkMask, density, densityAddition, timestep);
+    // }
+    // private native void addDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer[] x0, float timestep);
+
+    // /**
+    //  * Solve density diffusion
+    //  */
+    // private void solveDiffuseDensityWrapper(float timestep){
+    //     solveDiffuseDensity(DIM, chunkMask, density, densityAddition, uVector, vVector, wVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
+    // }
+    // private native void solveDiffuseDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer[] x0, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+
+    // /**
+    //  * Solve density diffusion
+    //  */
+    // private void advectDensityWrapper(float timestep){
+    //     advectDensity(DIM, chunkMask, density, densityAddition, uVector, vVector, wVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
+    // }
+    // private native void advectDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer[] x0, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+
+
+    // /**
+    //  * Sets the bounds of the neighbormap to neighbor values if available
+    //  */
+    // private void setBoundsToNeighborsWrapper(int vectorDir, ByteBuffer[] neighborMap){
+    //     setBoundsToNeighbors(DIM, chunkMask, vectorDir, neighborMap);
+    // }
+    // private native void setBoundsToNeighbors(int DIM_X, int chunkMask, int vectorDir, ByteBuffer[] neighborMap);
+
+    // /**
+    //  * Sets the bounds of the neighbormap to neighbor values if available, otherwise doesn't mess with them.
+    //  * This is to make sure zeroing out doesn't mess up the sim
+    //  */
+    // private void copyNeighborsWrapper(int vectorDir, int x, ByteBuffer[] neighborMap){
+    //     copyNeighbors(DIM, chunkMask, x, vectorDir, neighborMap);
+    // }
+    // private native void copyNeighbors(int DIM_X, int chunkMask, int x, int vectorDir, ByteBuffer[] neighborMap);
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 
 
@@ -957,7 +1114,7 @@ public class FluidSim {
                 for(int k = 0; k < DIM; k++){
                     index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
                     u0ArrayView[index] = 0;
-                    v0ArrayView[index] = densityArrayView[index] * GRAVITY;
+                    v0ArrayView[index] = 0;
                     w0ArrayView[index] = 0;
                 }
             }

src/main/java/electrosphere/Main.java

@@ -56,7 +56,7 @@ public class Main {
                 //Simulate
                 //
                 FluidSim.unlockThreadsWrapper();
-                FluidSim.simChunks(simArray,i,0.01f);
+                FluidSim.simChunks(simArray,i,0.0001f);
                 FluidSim.unlockThreadsWrapper();
                 time = time + (System.currentTimeMillis() - lastTime);
                 //