diff --git a/src/main/c/includes/electrosphere_FluidSim.h b/src/main/c/includes/electrosphere_FluidSim.h index 4d2531b..7c5adb8 100644 --- a/src/main/c/includes/electrosphere_FluidSim.h +++ b/src/main/c/includes/electrosphere_FluidSim.h @@ -14,7 +14,7 @@ extern "C" { #undef electrosphere_FluidSim_VISCOSITY_CONSTANT #define electrosphere_FluidSim_VISCOSITY_CONSTANT 0.0f #undef electrosphere_FluidSim_LINEARSOLVERTIMES -#define electrosphere_FluidSim_LINEARSOLVERTIMES 20L +#define electrosphere_FluidSim_LINEARSOLVERTIMES 10L #undef electrosphere_FluidSim_GRAVITY #define electrosphere_FluidSim_GRAVITY -100.0f /* diff --git a/src/main/c/velocitystep.c b/src/main/c/velocitystep.c index fcafbf3..e731d2b 100644 --- a/src/main/c/velocitystep.c +++ b/src/main/c/velocitystep.c @@ -196,59 +196,59 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection for(k=1; kN-1){ - // for(i=i-8; i < N-1; i++){ - // p[IX(i,j,k)] = (div[IX(i,j,k)] + a*(p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c; - // } - // } - for(i=1; i < N-1; i++){ - p[IX(i,j,k)] = (div[IX(i,j,k)] + a*(p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c; + for(i = 1; i < N-1; i=i+8){ + __m256 vector = _mm256_loadu_ps(&p[IX(i-1,j,k)]); + vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i+1,j,k)])); + vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j-1,k)])); + vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j+1,k)])); + vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j,k-1)])); + vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j,k+1)])); + // vector = _mm256_mul_ps(vector,aScalar); + vector = _mm256_add_ps(vector,_mm256_loadu_ps(&div[IX(i,j,k)])); + vector = _mm256_div_ps(vector,cScalar); + _mm256_storeu_ps(&p[IX(i,j,k)],vector); } + //If there is any leftover, perform manual solving + if(i>N-1){ + for(i=i-8; i < N-1; i++){ + p[IX(i,j,k)] = (div[IX(i,j,k)] + a*(p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c; + } + } + // for(i=1; i < N-1; i++){ + // p[IX(i,j,k)] = (div[IX(i,j,k)] + a*(p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c; + // } } } } @@ -347,58 +347,58 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection //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,xScalar); - // 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,xScalar); - // 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,yScalar); - // 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,yScalar); - // 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,zScalar); - // 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,zScalar); - // vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(9,j,k)]),vector); - // _mm256_storeu_ps(&w[IX(9,j,k)],vector); - for(i = 1; i < N-1; i++){ - u[IX(i,j,k)] = u[IX(i,j,k)] - 0.5 * (p[IX(i+1,j,k)] - p[IX(i-1,j,k)]) / h; - v[IX(i,j,k)] = v[IX(i,j,k)] - 0.5 * (p[IX(i,j+1,k)] - p[IX(i,j-1,k)]) / h; - w[IX(i,j,k)] = w[IX(i,j,k)] - 0.5 * (p[IX(i,j,k+1)] - p[IX(i,j,k-1)]) / h; - } + 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,xScalar); + 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,xScalar); + 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,yScalar); + 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,yScalar); + 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,zScalar); + 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,zScalar); + vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(9,j,k)]),vector); + _mm256_storeu_ps(&w[IX(9,j,k)],vector); + // for(i = 1; i < N-1; i++){ + // u[IX(i,j,k)] = u[IX(i,j,k)] - 0.5 * (p[IX(i+1,j,k)] - p[IX(i-1,j,k)]) / h; + // v[IX(i,j,k)] = v[IX(i,j,k)] - 0.5 * (p[IX(i,j+1,k)] - p[IX(i,j-1,k)]) / h; + // w[IX(i,j,k)] = w[IX(i,j,k)] - 0.5 * (p[IX(i,j,k+1)] - p[IX(i,j,k-1)]) / h; + // } } } } @@ -663,14 +663,6 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, j t0*u1*d0[IX(i1,j0,k1)]+ t1*u1*d0[IX(i1,j1,k1)] ); - // if(i == 1 && j == 1 && k == 1 && m == 2){ - // printf("%d %d %d\n",m,n,o); - // printf("%d %d %d %d %d %d\n",i0,i1,j0,j1,k0,k1); - // printf("%.2f vs\n",d0[IX(i,j,k)]); - // printf("%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n\n", - // d0[IX(i0,j0,k0)],d0[IX(i1,j0,k0)],d0[IX(i0,j1,k0)],d0[IX(i1,j1,k0)], - // d0[IX(i0,j0,k1)],d0[IX(i1,j0,k1)],d0[IX(i0,j1,k1)],d0[IX(i1,j1,k1)]); - // } } } } @@ -686,35 +678,16 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setBoundsToNeighbors int DIM = N; float * target = GET_ARR(env,neighborArray,CENTER_LOC); float * source; - // if(ARR_EXISTS(chunk_mask,0,1,1)){ - // source = GET_ARR(env,neighborArray,CK(0,1,1)); - // for(int x=1; x < DIM-1; x++){ - // for(int y = 1; y < DIM-1; y++){ - // target[IX(0,x,y)] = source[IX(DIM-2,x,y)]; - // } - // } - // } else { - for(int x=1; x < DIM-1; x++){ - for(int y = 1; y < DIM-1; y++){ - target[IX(0,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(1,x,y)] : target[IX(1,x,y)]; - } + for(int x=1; x < DIM-1; x++){ + for(int y = 1; y < DIM-1; y++){ + target[IX(0,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(1,x,y)] : target[IX(1,x,y)]; } - // } - - // if(ARR_EXISTS(chunk_mask,2,1,1)){ - // source = GET_ARR(env,neighborArray,CK(2,1,1)); - // for(int x=1; x < DIM-1; x++){ - // for(int y = 1; y < DIM-1; y++){ - // target[IX(DIM-1,x,y)] = source[IX(1,x,y)]; - // } - // } - // } else { - for(int x=1; x < DIM-1; x++){ - for(int y = 1; y < DIM-1; y++){ - target[IX(DIM-1,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)]; - } + } + for(int x=1; x < DIM-1; x++){ + for(int y = 1; y < DIM-1; y++){ + target[IX(DIM-1,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)]; } - // } + } for(int x=1; x < DIM-1; x++){ for(int y = 1; y < DIM-1; y++){ //((x)+(DIM)*(y) + (DIM)*(DIM)*(z)) diff --git a/src/main/java/electrosphere/FluidSim.java b/src/main/java/electrosphere/FluidSim.java index 6fb915c..862dfff 100644 --- a/src/main/java/electrosphere/FluidSim.java +++ b/src/main/java/electrosphere/FluidSim.java @@ -175,34 +175,13 @@ public class FluidSim { // //Vector stage solveChunkMask(simArray); - // System.out.println("Prior to add"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); addVectorSources(simArray, timestep); - // System.out.println("after add"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); swapAllVectorFields(simArray, timestep); - // System.out.println("after swap 11"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); solveVectorDiffusion(simArray, timestep); - // System.out.println("after diffuse"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); solveProjection(simArray, step, timestep); - // System.out.println("after proj 1"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); swapAllVectorFields(simArray, timestep); - // System.out.println("after swap 2"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); advectVectorsAcrossBoundaries(simArray, timestep); - // System.out.println("after advect"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); solveProjection(simArray, step, timestep); - // System.out.println("after proj 2"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); - // System.out.println("\n\n\n\n\n"); - // if(step == 7){ - // System.out.println(step); - // // System.exit(1); - // } // //Density stage @@ -213,18 +192,6 @@ public class FluidSim { advectDensity(simArray, timestep); // mirrorNeighborDensities(simArray, timestep); - // for(int x = 0; x < simArray.length; x++){ - // for(int y = 0; y < simArray[0].length; y++){ - // for(int z = 0; z < simArray[0][0].length; z++){ - // // - // //Reads out the results of the fluid sim - // // - // simArray[x][y][z].readDataIntoArrays(); - // } - // } - // } - // System.out.println("End density stage"); - // System.out.println(sumAllDensity(simArray)); @@ -455,41 +422,6 @@ public class FluidSim { } } } - // System.out.println("after swap in proj"); - // System.out.println(sumAllV(simArray) + " " + sumAllV0(simArray)); - // simArray[0][0][0].copyNeighborsWrapper(2, 0, simArray[0][0][0].vAdditionVector); - // simArray[1][0][0].copyNeighborsWrapper(2, 1, simArray[1][0][0].vAdditionVector); - // simArray[0][0][0].readDataIntoArrays(); - // simArray[1][0][0].readDataIntoArrays(); - // System.out.println("\n\n\n"); - // System.out.println("0 0 0"); - // for(int i = 0; i < 18; i++){ - // for(int j = 0; j < 18; j++){ - // System.out.print(simArray[0][0][0].vArrayView[IX(16,i,j)] + " "); - // } - // System.out.println(); - // } - // System.out.println("\n\n\n"); - // System.out.println("1 0 0"); - // for(int i = 0; i < 18; i++){ - // for(int j = 0; j < 18; j++){ - // System.out.print(simArray[1][0][0].vArrayView[IX(0,i,j)] + " "); - // } - // System.out.println(); - // } - // System.out.println("\n\n\n"); - // for(int i = 0; i < 18; i++){ - // for(int j = 0; j < 18; j++){ - // System.out.print(simArray[1][0][0].vArrayView[IX(0,i,j)] - simArray[0][0][0].vArrayView[IX(16,i,j)] + " "); - // } - // System.out.println(); - // } - // float value = simArray[1][0][0].vArrayView[IX(0,1,1)] - simArray[0][0][0].vArrayView[IX(16,1,1)]; - // System.out.println(simArray[0][0][0].vArrayView[IX(17,1,1)] + " " + simArray[1][0][0].vArrayView[IX(1,1,1)]); - // System.out.println("\n\n\n"); - // if(step == 1){ - // System.exit(1); - // } for(int x = 0; x < simArray.length; x++){ for(int y = 0; y < simArray[0].length; y++){ for(int z = 0; z < simArray[0][0].length; z++){ @@ -502,8 +434,6 @@ public class FluidSim { } } } - // System.out.println("after setup proj"); - // System.out.println(sumAllV(simArray) + " " + sumAllV0(simArray)); for(int x = 0; x < simArray.length; x++){ for(int y = 0; y < simArray[0].length; y++){ for(int z = 0; z < simArray[0][0].length; z++){ @@ -514,8 +444,6 @@ public class FluidSim { } } } - // System.out.println("after bound set 1 in proj"); - // System.out.println(sumAllV(simArray) + " " + sumAllV0(simArray)); //samples u0, v0 //sets u0 //these should have just been mirrored in the above @@ -540,11 +468,7 @@ public class FluidSim { } } } - // System.out.println("after proj iteration"); - // System.out.println(sumAllV(simArray) + " " + sumAllV0(simArray)); } - // System.out.println("after proj solver"); - // System.out.println(sumAllV(simArray) + " " + sumAllV0(simArray)); //samples u,v,w,u0 //sets u,v,w //Finalize projection @@ -557,8 +481,6 @@ public class FluidSim { } } } - // System.out.println("after finalize proj"); - // System.out.println(sumAllV(simArray) + " " + sumAllV0(simArray)); //set boundaries a final time for u,v,w //... for(int x = 0; x < simArray.length; x++){ @@ -656,8 +578,6 @@ public class FluidSim { } } } - // System.out.println("after first bound swap in advect"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); //samples u,v,w,u0,v0,w0 for(int x = 0; x < simArray.length; x++){ for(int y = 0; y < simArray[0].length; y++){ @@ -670,8 +590,6 @@ public class FluidSim { } } } - // System.out.println("before bound swap in advect"); - // System.out.println(sumAllU(simArray) + " " + sumAllU0(simArray)); //mirror neighbor data for(int x = 0; x < simArray.length; x++){ for(int y = 0; y < simArray[0].length; y++){ @@ -1018,16 +936,18 @@ public class FluidSim { FloatBuffer u0FloatView = uAdditionVector[13].asFloatBuffer(); FloatBuffer v0FloatView = vAdditionVector[13].asFloatBuffer(); FloatBuffer w0FloatView = wAdditionVector[13].asFloatBuffer(); + int index = 0; for(int i = 0; i < DIM; i++){ for(int j = 0; j < DIM; j++){ for(int k = 0; k < DIM; k++){ - densityArrayView[IX(i,j,k)] = xFloatView.get(); - uArrayView[IX(i,j,k)] = uFloatView.get(); - vArrayView[IX(i,j,k)] = vFloatView.get(); - wArrayView[IX(i,j,k)] = wFloatView.get(); - u0ArrayView[IX(i,j,k)] = u0FloatView.get(); - v0ArrayView[IX(i,j,k)] = v0FloatView.get(); - w0ArrayView[IX(i,j,k)] = w0FloatView.get(); + index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k)); + densityArrayView[index] = xFloatView.get(); + uArrayView[index] = uFloatView.get(); + vArrayView[index] = vFloatView.get(); + wArrayView[index] = wFloatView.get(); + u0ArrayView[index] = u0FloatView.get(); + v0ArrayView[index] = v0FloatView.get(); + w0ArrayView[index] = w0FloatView.get(); } } } @@ -1053,13 +973,15 @@ public class FluidSim { FloatBuffer u0FloatView = uAdditionVector[13].asFloatBuffer(); FloatBuffer v0FloatView = vAdditionVector[13].asFloatBuffer(); FloatBuffer w0FloatView = wAdditionVector[13].asFloatBuffer(); + int index = 0; for(int i = 0; i < DIM; i++){ for(int j = 0; j < DIM; j++){ for(int k = 0; k < DIM; k++){ - x0FloatView.put(density0ArrayView[IX(i,j,k)]); - u0FloatView.put(u0ArrayView[IX(i,j,k)]); - v0FloatView.put(v0ArrayView[IX(i,j,k)]); - w0FloatView.put(w0ArrayView[IX(i,j,k)]); + index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k)); + x0FloatView.put(density0ArrayView[index]); + u0FloatView.put(u0ArrayView[index]); + v0FloatView.put(v0ArrayView[index]); + w0FloatView.put(w0ArrayView[index]); } } } @@ -1069,12 +991,14 @@ public class FluidSim { * Adds gravity to the simulation */ private void addGravity(){ + int index = 0; for(int i = 0; i < DIM; i++){ for(int j = 0; j < DIM; j++){ for(int k = 0; k < DIM; k++){ - u0ArrayView[IX(i,j,k)] = 0; - v0ArrayView[IX(i,j,k)] = densityArrayView[IX(i,j,k)] * GRAVITY; - w0ArrayView[IX(i,j,k)] = 0; + index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k)); + u0ArrayView[index] = 0; + v0ArrayView[index] = densityArrayView[index] * GRAVITY; + w0ArrayView[index] = 0; } } } diff --git a/src/main/java/electrosphere/Main.java b/src/main/java/electrosphere/Main.java index ec418fb..0d7440b 100644 --- a/src/main/java/electrosphere/Main.java +++ b/src/main/java/electrosphere/Main.java @@ -70,8 +70,8 @@ public class Main { //redraw GLFWContext.redraw(meshArray); i++; - if(i == 1000){ - System.out.println(time / 1000.0); + if(i == 100){ + System.out.println(time / 100.0); } if(i > 3){ // scan.next();