up max multigrid iteration

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

@@ -5,12 +5,12 @@
 /**
  * The number of times to relax most solvers
  */
-#define FLUID_GRID2_LINEARSOLVERTIMES 2
+#define FLUID_GRID2_LINEARSOLVERTIMES 10
 
 /**
  * Width of a single grid cell
  */
-#define FLUID_GRID2_H (1.0/DIM)
+#define FLUID_GRID2_H (1.0/(DIM-2))
 
 /**
  * Timestep to simulate by

src/main/c/includes/math/ode/multigrid.h

@@ -13,14 +13,9 @@
  * @param phi0 The phi array from the last frame
  * @param a The a const
  * @param c The c const
+ * @return The residual
  */
-void solver_multigrid_iterate(float * phi, float * phi0, float a, float c);
-
-
-/**
- * Allocates the data for the multigrid solver
- */
-LIBRARY_API void solver_multigrid_allocate();
+float solver_multigrid_iterate(float * phi, float * phi0, float a, float c);
 
 
 

src/main/c/includes/util/vector.h

@@ -0,0 +1,14 @@
+#ifndef UTIL_VECTOR_H
+#define UTIL_VECTOR_H
+
+#include <immintrin.h>
+
+#include "public.h"
+
+
+/**
+ * Sums a float vector's elements
+ */
+LIBRARY_API float vec_sum(__m256 x);
+
+#endif

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

@@ -144,7 +144,6 @@ JNIEXPORT void JNICALL Java_electrosphere_server_fluid_simulator_FluidAccelerate
 
     //init grid2 sim
     fluid_grid2_allocate_arrays();
-    solver_multigrid_allocate();
 }
 
 /**

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

@@ -130,10 +130,8 @@ LIBRARY_API void fluid_grid2_simulate(
     for(int i = 0; i < numChunks; i++){
         Chunk * currentChunk = chunks[i];
 
-        environment->state.existingDensity = environment->state.existingDensity + fluid_grid2_calculateSum(currentChunk->d);
-        environment->state.newDensity = environment->state.newDensity + fluid_grid2_calculateSum(currentChunk->d0);
-
         fluid_grid2_addDensity(environment,currentChunk->d,currentChunk->d0,timestep);
+        environment->state.existingDensity = environment->state.existingDensity + fluid_grid2_calculateSum(currentChunk->d);
 
         //swap all density arrays
         fluid_grid2_flip_arrays(currentChunk->d,currentChunk->d0);

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

@@ -178,9 +178,12 @@ LIBRARY_API void fluid_grid2_solveProjection(
     float * div = GET_ARR_RAW(jrv0,CENTER_LOC);
 
     //perform iteration of v cycle multigrid method
-    for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
-        solver_multigrid_iterate(p,div,a,c);
+    float residual = 1;
+    int iteration = 0;
+    while(iteration < FLUID_GRID2_LINEARSOLVERTIMES && residual > FLUID_GRID2_REALLY_SMALL_VALUE){
+        residual = solver_multigrid_iterate(p,div,a,c);
         fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
+        iteration++;
     }
 
     // solver_conjugate_gradient_solve_serial(p,div,a,c);

src/main/c/src/math/ode/multigrid.c

@@ -5,6 +5,7 @@
 #include "fluid/sim/grid2/utilities.h"
 #include "math/ode/gauss_seidel.h"
 #include "math/ode/multigrid.h"
+#include "util/vector.h"
 
 /**
  * Dimension of the half resolution grid
@@ -19,14 +20,14 @@ static int quarterDim = ((DIM - 2) / 4) + 2;
 /**
  * The half resolution grids
  */
-static float * halfGridPhi;
-static float * halfGridPhi0;
+static float * halfGridPhi = NULL;
+static float * halfGridPhi0 = NULL;
 
 /**
  * The quarter resolution grids
  */
-static float * quarterGridPhi;
-static float * quarterGridPhi0;
+static float * quarterGridPhi = NULL;
+static float * quarterGridPhi0 = NULL;
 
 /**
  * Relaxes an ODE matrix by 1 iteration of multigrid method
@@ -34,8 +35,21 @@ static float * quarterGridPhi0;
  * @param phi0 The phi array from the last frame
  * @param a The a const
  * @param c The c const
+ * @return The residual
  */
-void solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
+float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
+    if(halfGridPhi == NULL){
+        halfGridPhi = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
+    }
+    if(halfGridPhi0 == NULL){
+        halfGridPhi0 = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
+    }
+    if(quarterGridPhi == NULL){
+        quarterGridPhi = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
+    }
+    if(quarterGridPhi0 == NULL){
+        quarterGridPhi0 = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
+    }
 
     //
     //gauss-seidel at highest res
@@ -73,8 +87,8 @@ void solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
     for(int x = 0; x < quarterDim - 1; x++){
         for(int y = 0; y < quarterDim - 1; y++){
             for(int z = 0; z < quarterDim - 1; z++){
-                halfGridPhi[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = 0;
-                halfGridPhi0[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = 0;
+                quarterGridPhi[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = 0;
+                quarterGridPhi0[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = 0;
             }
         }
     }
@@ -135,17 +149,44 @@ void solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
     //
     //gauss-seidel at highest res
     solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,DIM);
+
+    //calculate residual
+    __m256 aScalar = _mm256_set1_ps(a);
+    __m256 cScalar = _mm256_set1_ps(c);
+    __m256 collector = _mm256_setzero_ps();
+    __m256 vector;
+    float residual = 1;
+    //transform u direction
+    // int i, j, k;
+    // for(k=1; k<DIM-2; k++){
+    //     for(j=1; j<DIM-2; j++){
+    //         //lower
+    //         i = 1;
+    //         vector = _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,DIM)]);
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,DIM)]));
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,DIM)]));
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,DIM)]));
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,DIM)]));
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,DIM)]));
+    //         vector = _mm256_mul_ps(vector,aScalar);
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,DIM)]));
+    //         vector = _mm256_div_ps(vector,cScalar);
+    //         collector = _mm256_add_ps(collector,vector);
+
+    //         //upper
+    //         i = 9;
+    //         vector = _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,DIM)]);
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,DIM)]));
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,DIM)]));
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,DIM)]));
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,DIM)]));
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,DIM)]));
+    //         vector = _mm256_mul_ps(vector,aScalar);
+    //         vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,DIM)]));
+    //         vector = _mm256_div_ps(vector,cScalar);
+    //         collector = _mm256_add_ps(collector,vector);
+    //     }
+    // }
+    // residual = vec_sum(collector);
+    return residual;
 }
-
-
-/**
- * Allocates the data for the multigrid solver
- */
-LIBRARY_API void solver_multigrid_allocate(){
-    halfGridPhi = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
-    halfGridPhi0 = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
-    quarterGridPhi = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
-    quarterGridPhi0 = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
-}
-
-

src/main/c/src/util/vector.c

@@ -0,0 +1,29 @@
+#include "util/vector.h"
+
+
+
+
+/**
+ * Sums a float vector's elements
+ */
+LIBRARY_API float vec_sum(__m256 x) {
+    // hiQuad = ( x7, x6, x5, x4 )
+    const __m128 hiQuad = _mm256_extractf128_ps(x, 1);
+    // loQuad = ( x3, x2, x1, x0 )
+    const __m128 loQuad = _mm256_castps256_ps128(x);
+    // sumQuad = ( x3 + x7, x2 + x6, x1 + x5, x0 + x4 )
+    const __m128 sumQuad = _mm_add_ps(loQuad, hiQuad);
+    // loDual = ( -, -, x1 + x5, x0 + x4 )
+    const __m128 loDual = sumQuad;
+    // hiDual = ( -, -, x3 + x7, x2 + x6 )
+    const __m128 hiDual = _mm_movehl_ps(sumQuad, sumQuad);
+    // sumDual = ( -, -, x1 + x3 + x5 + x7, x0 + x2 + x4 + x6 )
+    const __m128 sumDual = _mm_add_ps(loDual, hiDual);
+    // lo = ( -, -, -, x0 + x2 + x4 + x6 )
+    const __m128 lo = sumDual;
+    // hi = ( -, -, -, x1 + x3 + x5 + x7 )
+    const __m128 hi = _mm_shuffle_ps(sumDual, sumDual, 0x1);
+    // sum = ( -, -, -, x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 )
+    const __m128 sum = _mm_add_ss(lo, hi);
+    return _mm_cvtss_f32(sum);
+}

src/test/c/fluid/sim/grid2/add_dens_tests.c

@@ -60,8 +60,9 @@ int fluid_sim_grid2_add_dens_test1(){
 
     //test the result
     float afterSum = chunk_queue_sum_density(queue);
-    if(fabs(beforeSum - afterSum) > FLUID_GRID2_PROJECTION_ERROR_MARGIN){
-        rVal += assertEqualsFloat(beforeSum,afterSum,"Advection changed density!  %f %f  \n");
+    float expectedSum = additionFrameCutoff * MAX_FLUID_VALUE * FLUID_GRID2_SIM_STEP;
+    if(fabs(expectedSum - afterSum) > FLUID_GRID2_PROJECTION_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expectedSum,afterSum,"Simulation did not properly add density! expected: %f     actual: %f  \n");
     }
 
     return rVal;
@@ -73,10 +74,9 @@ int fluid_sim_grid2_add_dens_test1(){
 int fluid_sim_grid2_add_dens_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
     fluid_grid2_allocate_arrays();
 
-    // rVal += fluid_sim_grid2_add_dens_test1();
+    rVal += fluid_sim_grid2_add_dens_test1();
 
     return rVal;
 }

src/test/c/fluid/sim/grid2/advect_projection_tests.c

@@ -156,8 +156,6 @@ int fluid_sim_grid2_advect_projection_compute_error_over_time(){
 int fluid_sim_grid2_advect_projection_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
-
     // rVal += fluid_sim_grid2_advect_projection_test1();
     // rVal += fluid_sim_grid2_advect_projection_compute_error_over_time();
 

src/test/c/fluid/sim/grid2/density_advection_tests.c

@@ -228,8 +228,6 @@ int fluid_sim_grid2_density_advection_test5(){
 int fluid_sim_grid2_density_advection_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
-
     rVal += fluid_sim_grid2_density_advection_test1();
     rVal += fluid_sim_grid2_density_advection_test2();
     rVal += fluid_sim_grid2_density_advection_test3();

src/test/c/fluid/sim/grid2/density_diffuse_tests.c

@@ -162,8 +162,6 @@ int fluid_sim_grid2_density_diffuse_test2(){
 int fluid_sim_grid2_density_diffuse_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
-
     rVal += fluid_sim_grid2_density_diffuse_test1();
     rVal += fluid_sim_grid2_density_diffuse_test2();
 

src/test/c/fluid/sim/grid2/finalize_projection_tests.c

@@ -78,8 +78,6 @@ int fluid_sim_grid2_finalize_projection_test1(){
 int fluid_sim_grid2_finalize_projection_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
-
     rVal += fluid_sim_grid2_finalize_projection_test1();
 
     return rVal;

src/test/c/fluid/sim/grid2/full_sim_tests.c

@@ -137,7 +137,6 @@ int fluid_sim_grid2_full_sim_test3(){
 int fluid_sim_grid2_full_sim_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
     fluid_grid2_allocate_arrays();
 
     rVal += fluid_sim_grid2_full_sim_test1();

src/test/c/fluid/sim/grid2/setup_projection_tests.c

@@ -82,8 +82,6 @@ int fluid_sim_grid2_setup_projection_test2(){
 int fluid_sim_grid2_setup_projection_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
-
     rVal += fluid_sim_grid2_setup_projection_test1();
     rVal += fluid_sim_grid2_setup_projection_test2();
 

src/test/c/fluid/sim/grid2/solve_projection_tests.c

@@ -39,7 +39,6 @@ int fluid_sim_grid2_solve_projection_test1(){
 
     //actually solve
     fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
-    fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0[CENTER_LOC]);
 
     //test the result
     float expected, actual;
@@ -106,8 +105,6 @@ int fluid_sim_grid2_solve_projection_test1(){
 int fluid_sim_grid2_solve_projection_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
-
     rVal += fluid_sim_grid2_solve_projection_test1();
 
     return rVal;

src/test/c/fluid/sim/grid2/velocity_advection_tests.c

@@ -183,8 +183,6 @@ int fluid_sim_grid2_velocity_advection_test3(){
 int fluid_sim_grid2_velocity_advection_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
-
     rVal += fluid_sim_grid2_velocity_advection_test1();
     rVal += fluid_sim_grid2_velocity_advection_test2();
     rVal += fluid_sim_grid2_velocity_advection_test3();

src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c

@@ -150,8 +150,6 @@ int fluid_sim_grid2_velocity_diffuse_test2(){
 int fluid_sim_grid2_velocity_diffuse_tests(int argc, char **argv){
     int rVal = 0;
 
-    solver_multigrid_allocate();
-
     rVal += fluid_sim_grid2_velocity_diffuse_test1();
 
     return rVal;