grid2 work

.vscode/settings.json

@@ -45,6 +45,7 @@
         "math.h": "c",
         "multigrid.h": "c",
         "gauss_seidel.h": "c",
-        "ode_utils.h": "c"
+        "ode_utils.h": "c",
+        "util.h": "c"
     }
 }

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

@@ -5,23 +5,28 @@
 /**
  * The number of times to relax most solvers
  */
-#define FLUID_GRID2_LINEARSOLVERTIMES 3
+#define FLUID_GRID2_LINEARSOLVERTIMES 2
 
 /**
  * The number of times to relax the vector diffusion solver
  */
-#define FLUID_GRID2_VECTOR_DIFFUSE_TIMES 3
+#define FLUID_GRID2_VECTOR_DIFFUSE_TIMES 2
 
 /**
  * Width of a single grid cell
  */
-#define FLUID_GRID2_H (1.0/(DIM-2))
+#define FLUID_GRID2_H (1.0/DIM)
 
 /**
  * Timestep to simulate by
  */
 #define FLUID_GRID2_SIM_STEP 0.01f
 
+/**
+ * Const to multiply the advection stages by to offset numeric instability
+ */
+#define FLUID_GRID2_CORRECTION_CONST 1.0005f
+
 /**
  * Really small value used for something
  */

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

@@ -8,8 +8,9 @@
  * @param phi0 The phi array from the last frame
  * @param a The a const
  * @param c The c const
+ * @return 1 if the system has already been solved, 0 otherwise
  */
-void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c);
+int solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c);
 
 /**
  * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method
@@ -18,6 +19,15 @@ void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c)
  * @param a The a const
  * @param c The c const
  */
-void solver_conjugate_gradient_iterate(float * phi, float * phi0, float a, float c);
+void solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float a, float c);
+
+/**
+ * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method serially
+ * @param phi The phi array
+ * @param phi0 The phi array from the last frame
+ * @param a The a const
+ * @param c The c const
+ */
+void solver_conjugate_gradient_solve_serial(float * phi, float * phi0, float a, float c);
 
 #endif

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

@@ -8,6 +8,7 @@
 #include "fluid/queue/chunk.h"
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/utilities.h"
+#include "math/ode/multigrid.h"
 
 
 /**
@@ -54,6 +55,10 @@ LIBRARY_API void fluid_grid2_solveDiffuseDensity(
     __m256 aScalar = _mm256_set1_ps(a);
     __m256 cScalar = _mm256_set1_ps(c);
 
+    // for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+    //     solver_multigrid_iterate(x,x0,a,c);
+    // }
+
     //transform u direction
     for(k=1; k<DIM-1; k++){
         for(j=1; j<DIM-1; j++){
@@ -158,6 +163,16 @@ LIBRARY_API void fluid_grid2_advectDensity(float ** d, float ** d0, float ** ur,
 
                 u1 = z-k0;
                 u0 = 1-u1;
+
+                if(
+                    i0 < 0 || j0 < 0 || k0 < 0 ||
+                    i0 > DIM-1 || j0 > DIM-1 || k0 > DIM-1 ||
+                    i1 < 0 || j1 < 0 || k1 < 0 ||
+                    i1 > DIM-1 || j1 > DIM-1 || k1 > DIM-1
+                ){
+                    printf("advect dens: %d %d %d   %d %d %d ---  %f %f %f\n", i0, j0, k0, i1, j1, k1, x, y, z);
+                    fflush(stdout);
+                }
                 
                 center_d[IX(i,j,k)] = 
                 s0*(

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

@@ -98,10 +98,8 @@ void fluid_grid2_simulate(
         //these should have just been mirrored in the above
         //
         //Perform main projection solver
-        for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+        fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,timestep);
-            fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,timestep);
+        fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-            fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-        }
         //samples u,v,w,u0
         //sets u,v,w
         //Finalize projection

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

@@ -239,12 +239,18 @@ LIBRARY_API void fluid_grid2_solveProjection(
     float * p = GET_ARR_RAW(jru0,CENTER_LOC);
     float * div = GET_ARR_RAW(jrv0,CENTER_LOC);
 
-    //transform u direction
+    //perform iteration of v cycle multigrid method
-    solver_gauss_seidel_iterate_parallel(p,div,a,c,DIM);
+    for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
-    solver_multigrid_iterate(p,div,a,c);
+        solver_multigrid_iterate(p,div,a,c);
-    solver_gauss_seidel_iterate_parallel(p,div,a,c,DIM);
+    }
-    solver_conjugate_gradient_init(p,div,a,c);
+
-    solver_conjugate_gradient_iterate(p,div,a,c);
+    // solver_conjugate_gradient_solve_serial(p,div,a,c);
+
+    //finest grain iteration with conjugate gradient method
+    // int shouldSolve = solver_conjugate_gradient_init(p,div,a,c);
+    // if(shouldSolve < 1){
+    //     solver_conjugate_gradient_iterate(p,div,a,c);
+    // }
 }
 
 /**
@@ -419,6 +425,16 @@ void fluid_grid2_advect_velocity(int b, float ** jrd, float ** jrd0, float * u,
 
                 u1 = z-k0;
                 u0 = 1-u1;
+
+                if(
+                    i0 < 0 || j0 < 0 || k0 < 0 ||
+                    i0 > DIM-1 || j0 > DIM-1 || k0 > DIM-1 ||
+                    i1 < 0 || j1 < 0 || k1 < 0 ||
+                    i1 > DIM-1 || j1 > DIM-1 || k1 > DIM-1
+                ){
+                    printf("advect vel: %d %d %d   %d %d %d ---  %f %f %f\n", i0, j0, k0, i1, j1, k1, x, y, z);
+                    fflush(stdout);
+                }
                 
 
                 d[IX(i,j,k)] = 

src/main/c/src/math/ode/conjugate_gradient.c

@@ -3,18 +3,33 @@
 #include <immintrin.h>
 
 #include "fluid/queue/chunk.h"
+#include "fluid/sim/grid2/solver_consts.h"
 #include "math/ode/ode_utils.h"
 
 
+/**
+ * Used for preventing divide by 0's
+ */
+static float CONJUGATE_GRADIENT_EPSILON = 1e-6;
+
+static float CONJUGATE_GRADIENT_CONVERGENCE_THRESHOLD = 0.0001f;
+
 /**
  * The search direction array
  */
-float * p = NULL;
+static float * p = NULL;
+
+/**
+ * Maxmum frames to iterate for
+ */
+static int max_frames = 100;
 
 /**
  * The residual array
  */
-float * r = NULL;
+static float * r = NULL;
+
+static float * A = NULL;
 
 
 /**
@@ -23,14 +38,18 @@ float * r = NULL;
  * @param phi0 The phi array from the last frame
  * @param a The a const
  * @param c The c const
+ * @return 1 if the system has already been solved, 0 otherwise
  */
-void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c){
+int solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c){
     if(p == NULL){
         p = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
     }
     if(r == NULL){
         r = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
     }
+    if(A == NULL){
+        A = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
+    }
 
     int i, j, k;
     
@@ -41,7 +60,17 @@ void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c)
     //iniitalize the r (residual) and p (search direction) arrays
     for(k=1; k<DIM-1; k++){
         for(j=1; j<DIM-1; j++){
-            int n = 0;
+            //set borders
+            i = 0;
+            _mm256_storeu_ps(&r[ode_index(i,j,k,DIM)],_mm256_setzero_ps());
+            _mm256_storeu_ps(&p[ode_index(i,j,k,DIM)],_mm256_setzero_ps());
+            i = 8;
+            _mm256_storeu_ps(&r[ode_index(i,j,k,DIM)],_mm256_setzero_ps());
+            _mm256_storeu_ps(&p[ode_index(i,j,k,DIM)],_mm256_setzero_ps());
+            i = DIM-9;
+            _mm256_storeu_ps(&r[ode_index(i,j,k,DIM)],_mm256_setzero_ps());
+            _mm256_storeu_ps(&p[ode_index(i,j,k,DIM)],_mm256_setzero_ps());
+
             //solve as much as possible vectorized
             for(i = 1; i < DIM-1; i=i+8){
                 //calculate the stencil applied to phi
@@ -71,27 +100,35 @@ void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c)
             }
         }
     }
+    float sampleResidual = r[ode_index(3,3,3,DIM)];
+    if(sampleResidual <= CONJUGATE_GRADIENT_EPSILON){
+        return 1;
+    }
+    printf("sampleResidual: %f\n",sampleResidual);
+    return 0;
 }
 
 /**
- * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method
+ * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method parallelly
  * @param phi The phi array
  * @param phi0 The phi array from the last frame
  * @param a The a const
  * @param c The c const
  */
-void solver_conjugate_gradient_iterate(float * phi, float * phi0, float a, float c){
+void solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float a, float c){
     int i, j, k;
     
     __m256 aScalar = _mm256_set1_ps(a);
     __m256 cScalar = _mm256_set1_ps(c);
-    __m256 f, residual, stencil;
+    __m256 epsilonScalar = _mm256_set1_ps(CONJUGATE_GRADIENT_EPSILON);
-    __m256 denominator, r0_squared, r1_squared, a_k, r_delta, beta, p0, p1, p_new, phi_old, phi1, r0, r1;
+    __m256 stencil;
+    __m256 denominator, r_old_squared, r_new_squared, a_k, r_delta, beta, p_old, p_new, phi_old, phi_new, r_old, r_new;
+    __m256 nanMask, nanComp;
+    __m256i iMask0, iMask1;
 
     //iniitalize the r (residual) and p (search direction) arrays
     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){
                 //calculate the stencil applied to p
@@ -110,34 +147,143 @@ void solver_conjugate_gradient_iterate(float * phi, float * phi0, float a, float
                 stencil = _mm256_div_ps(stencil,cScalar);
 
                 //load p0
-                p0 = _mm256_loadu_ps(&p[ode_index(i,j,k,DIM)]);
+                p_old = _mm256_loadu_ps(&p[ode_index(i,j,k,DIM)]);
 
                 //calculate ak
-                denominator = _mm256_mul_ps(stencil,p0);
+                denominator = _mm256_mul_ps(stencil,p_old);
-                r0_squared = _mm256_mul_ps(_mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]),_mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]));
+                denominator = _mm256_add_ps(denominator,epsilonScalar);
-                a_k = _mm256_div_ps(r0_squared,denominator);
+                r_old_squared = _mm256_mul_ps(_mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]),_mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]));
+                a_k = _mm256_div_ps(r_old_squared,denominator);
+
 
                 //update phi
                 phi_old = _mm256_loadu_ps(&phi[ode_index(i,j,k,DIM)]);
-                phi1 = _mm256_mul_ps(p1,a_k);
+                phi_new = _mm256_mul_ps(p_old,a_k);
-                phi1 = _mm256_add_ps(phi1,phi_old);
+                phi_new = _mm256_add_ps(phi_new,phi_old);
-                _mm256_storeu_ps(&phi[ode_index(i,j,k,DIM)],phi1);
+
+                //store new phi
+                _mm256_storeu_ps(&phi[ode_index(i,j,k,DIM)],phi_new);
 
                 //update residual
                 r_delta = _mm256_add_ps(stencil,a_k);
-                r0 = _mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]);
+                r_old = _mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]);
-                r1 = _mm256_sub_ps(r0,r_delta);
+                r_new = _mm256_sub_ps(r_old,r_delta);
-                _mm256_storeu_ps(&r[ode_index(i,j,k,DIM)],r1);
+                _mm256_storeu_ps(&r[ode_index(i,j,k,DIM)],r_new);
 
                 //calculate beta
-                r1_squared = _mm256_mul_ps(r1,r1);
+                r_new_squared = _mm256_mul_ps(r_new,r_new);
-                beta = _mm256_div_ps(r1_squared,r0_squared);
+                denominator = _mm256_add_ps(r_old_squared,epsilonScalar);
+                beta = _mm256_div_ps(r_new_squared,denominator);
 
                 //update p (search dir)
-                p_new = _mm256_mul_ps(beta,p0);
+                p_new = _mm256_mul_ps(beta,p_old);
-                p1 = _mm256_add_ps(r1,p_new);
+                p_new = _mm256_add_ps(r_new,p_new);
-                _mm256_storeu_ps(&p[ode_index(i,j,k,DIM)],p1);
+                _mm256_storeu_ps(&p[ode_index(i,j,k,DIM)],p_new);
             }
         }
     }
-}
+}
+
+/**
+ * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method serially
+ * @param phi The phi array
+ * @param phi0 The phi array from the last frame
+ * @param a The a const
+ * @param c The c const
+ */
+void solver_conjugate_gradient_solve_serial(float * phi, float * phi0, float a, float c){
+    int i, j, k;
+    if(p == NULL){
+        p = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
+    }
+    if(r == NULL){
+        r = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
+    }
+    if(A == NULL){
+        A = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
+    }
+    
+    //iniitalize the r (residual) and p (search direction) arrays
+    for(k=1; k<DIM-1; k++){
+        for(j=1; j<DIM-1; j++){
+            for(i = 1; i < DIM-1; i++){
+                r[ode_index(i,j,k,DIM)] = phi0[ode_index(i,j,k,DIM)];
+                p[ode_index(i,j,k,DIM)] = r[ode_index(i,j,k,DIM)];
+            }
+        }
+    }
+
+    float convergence = 1;
+    int framecount = 0;
+    while(convergence > CONJUGATE_GRADIENT_CONVERGENCE_THRESHOLD && framecount < max_frames){
+        //solve Ap
+        for(k=1; k<DIM-1; k++){
+            for(j=1; j<DIM-1; j++){
+                for(i = 1; i < DIM-1; i++){
+                    A[ode_index(i,j,k,DIM)] = 
+                    (
+                        p[ode_index (i   , j   , k   ,DIM)] + 
+                        a * (
+                            p[ode_index (i+1 , j   , k   ,DIM)] +
+                            p[ode_index (i-1 , j   , k   ,DIM)] +
+                            p[ode_index (i   , j+1 , k   ,DIM)] +
+                            p[ode_index (i   , j-1 , k   ,DIM)] +
+                            p[ode_index (i   , j   , k+1 ,DIM)] +
+                            p[ode_index (i   , j   , k-1 ,DIM)]
+                        )
+                    ) / c;
+                }
+            }
+        }
+        convergence = 0;
+        float denominator = CONJUGATE_GRADIENT_EPSILON;
+        for(k=1; k<DIM-1; k++){
+            for(j=1; j<DIM-1; j++){
+                for(i = 1; i < DIM-1; i++){
+                    convergence = convergence + r[ode_index(i,j,k,DIM)] * r[ode_index(i,j,k,DIM)];
+                    denominator = denominator + p[ode_index(i,j,k,DIM)] * A[ode_index(i,j,k,DIM)];
+                }
+            }
+        }
+        if(denominator < CONJUGATE_GRADIENT_EPSILON && denominator > -CONJUGATE_GRADIENT_EPSILON){
+            printf("Divide by 0! %f \n", denominator);
+            fflush(stdout);
+            break;
+        }
+        float alpha = convergence / denominator;
+        for(k=1; k<DIM-1; k++){
+            for(j=1; j<DIM-1; j++){
+                for(i = 1; i < DIM-1; i++){
+                    phi[ode_index(i,j,k,DIM)] = phi[ode_index(i,j,k,DIM)] + alpha * p[ode_index(i,j,k,DIM)];
+                    r[ode_index(i,j,k,DIM)] = r[ode_index(i,j,k,DIM)] - alpha * A[ode_index(i,j,k,DIM)];
+                }
+            }
+        }
+
+        float r_new_dot = 0;
+        for(k=1; k<DIM-1; k++){
+            for(j=1; j<DIM-1; j++){
+                for(i = 1; i < DIM-1; i++){
+                    r_new_dot = r_new_dot + r[ode_index(i,j,k,DIM)] * r[ode_index(i,j,k,DIM)];
+                }
+            }
+        }
+
+        float beta = r_new_dot / convergence;
+
+        for(k=1; k<DIM-1; k++){
+            for(j=1; j<DIM-1; j++){
+                for(i = 1; i < DIM-1; i++){
+                    p[ode_index(i,j,k,DIM)] = r[ode_index(i,j,k,DIM)] + beta * p[ode_index(i,j,k,DIM)];
+                }
+            }
+        }
+        framecount++;
+    }
+    if(framecount >= 100 || convergence > CONJUGATE_GRADIENT_CONVERGENCE_THRESHOLD){
+        printf("Failed to converge! %f  \n", convergence);
+        fflush(stdout);
+    }
+
+}
+

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

@@ -36,6 +36,10 @@ static float * quarterGridPhi0;
  */
 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);
+
     //
     //half res
     for(int x = 1; x < halfDim - 2; x++){
@@ -89,6 +93,10 @@ 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);
 }
 
 

src/test/c/CMakeLists.txt

@@ -58,6 +58,8 @@ foreach (TEST_FILE ${TEST_SOURCES})
   endif()
 endforeach ()
 
+target_compile_options(test_runner PRIVATE -m64 -mavx -mavx2)
+
 # make test runner depend on library
 add_dependencies(test_runner StormEngine)
 

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

@@ -62,10 +62,8 @@ int fluid_sim_grid2_advect_projection_test1(){
             fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
             fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u0);
             fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v0);
-            for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+            fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
-                fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
+            fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-                fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-            }
             fluid_grid2_finalizeProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
 
             //advect density
@@ -129,10 +127,8 @@ int fluid_sim_grid2_advect_projection_compute_error_over_time(){
                     fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
                     fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u0);
                     fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v0);
-                    for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+                    fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
-                        fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
+                    fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-                        fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-                    }
                     fluid_grid2_finalizeProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
 
                     //advect density

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

@@ -36,10 +36,8 @@ int fluid_sim_grid2_finalize_projection_test1(){
     Chunk * currentChunk = queue[0];
     currentChunk->u[CENTER_LOC][IX(3,3,3)] = 1.0f;
     fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
-    for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+    fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
-        fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
+    fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-        fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-    }
 
     //finalize
     fluid_grid2_finalizeProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);

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

@@ -41,10 +41,10 @@ int fluid_sim_grid2_setup_projection_test1(){
     fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
 
     //test the result
-    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],-0.5f / DIM,"Divergence of the vector field at 3,3,3 should be -0.5/DIM!   %f %f  \n");
+    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],-0.5f * FLUID_GRID2_H,"Divergence of the vector field at 3,3,3 should be -0.5 * h!   actual: %f    expected: %f  \n");
-    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(3,3,3)],0,"Divergence of the vector field at 3,3,3 should be 0!   %f %f  \n");
+    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(3,3,3)],0,"Divergence of the vector field at 3,3,3 should be 0!   actual: %f    expected: %f  \n");
-    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],0.5f / DIM,"Divergence of the vector field at 4,3,3 should be 0.5/DIM!   %f %f  \n");
+    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],0.5f * FLUID_GRID2_H,"Divergence of the vector field at 4,3,3 should be 0.5 * h!   actual: %f    expected: %f  \n");
-    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(5,3,3)],0,"Divergence of the vector field at 5,3,3 should be 0!   %f %f  \n");
+    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(5,3,3)],0,"Divergence of the vector field at 5,3,3 should be 0!   actual: %f    expected: %f  \n");
 
     return rVal;
 }
@@ -53,7 +53,7 @@ int fluid_sim_grid2_setup_projection_test1(){
  * Testing velocity advection
  */
 int fluid_sim_grid2_setup_projection_test2(){
-    printf("fluid_sim_grid2_setup_projection_test1\n");
+    printf("fluid_sim_grid2_setup_projection_test2\n");
     int rVal = 0;
     Environment * env = fluid_environment_create();
     Chunk ** queue = NULL;
@@ -69,9 +69,9 @@ int fluid_sim_grid2_setup_projection_test2(){
     fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
 
     //test the result
-    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],0.5f / DIM,"Divergence of the vector field at 3,3,3 should be 0.5/DIM!   %f %f  \n");
+    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],0.5f * FLUID_GRID2_H,"Divergence of the vector field at 3,3,3 should be 0.5 * h!   actual: %f    expected: %f  \n");
     rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(3,3,3)],0,"Divergence of the vector field at 3,3,3 should be 0!   %f %f  \n");
-    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],-0.5f / DIM,"Divergence of the vector field at 4,3,3 should be -0.5/DIM!   %f %f  \n");
+    rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],-0.5f * FLUID_GRID2_H,"Divergence of the vector field at 4,3,3 should be -0.5 * h!   actual: %f    expected: %f  \n");
 
     return rVal;
 }
@@ -84,8 +84,8 @@ int fluid_sim_grid2_setup_projection_tests(int argc, char **argv){
 
     solver_multigrid_allocate();
 
-    // rVal += fluid_sim_grid2_setup_projection_test1();
+    rVal += fluid_sim_grid2_setup_projection_test1();
-    // rVal += fluid_sim_grid2_setup_projection_test2();
+    rVal += fluid_sim_grid2_setup_projection_test2();
 
     return rVal;
 }

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

@@ -38,10 +38,8 @@ int fluid_sim_grid2_solve_projection_test1(){
     fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
 
     //actually solve
-    for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+    fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
-        fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
+    fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-        fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
-    }
 
     //test the result
     float expected, actual;