parallel multigrid performance work

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

@@ -5,7 +5,7 @@
 /**
  * The number of times to relax most solvers
  */
-#define FLUID_GRID2_LINEARSOLVERTIMES 10
+#define FLUID_GRID2_LINEARSOLVERTIMES 20
 
 /**
  * Width of a single grid cell

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

@@ -4,7 +4,7 @@
 #include "math/ode/gauss_seidel.h"
 #include "math/ode/multigrid.h"
 
-#define MULTIGRID_PARALLEL_LOWEST_PARALLEL_DIM ((DIM - 2) / 1) + 2
+#define MULTIGRID_PARALLEL_LOWEST_PARALLEL_DIM ((DIM - 2) / 2) + 2
 
 
 
@@ -246,7 +246,7 @@ static inline void solver_multigrid_parallel_store_residual(float * phi, float *
  * @param GRIDDIM The dimension of the phi grid
  * @return The residual
  */
-static inline void solver_multigrid_parallel_iterate_recurse(float * phi, float * phi0, float a, float c, int GRIDDIM){
+static inline void solver_multigrid_parallel_iterate_recurse_v_cycle(float * phi, float * phi0, float a, float c, int GRIDDIM){
     int LOWERDIM = ((GRIDDIM - 2) / 2) + 2;
     float * currResidual = solver_multigrid_get_current_residual(GRIDDIM);
     float * lowerPhi = solver_multigrid_get_current_phi(LOWERDIM);
@@ -266,7 +266,47 @@ static inline void solver_multigrid_parallel_iterate_recurse(float * phi, float
     if(GRIDDIM <= MULTIGRID_PARALLEL_LOWEST_PARALLEL_DIM){
         solver_multigrid_iterate_serial_recursive(lowerPhi,lowerPhi0,a,c,LOWERDIM);
     } else {
-        solver_multigrid_parallel_iterate_recurse(lowerPhi,lowerPhi0,a,c,LOWERDIM);
+        solver_multigrid_parallel_iterate_recurse_v_cycle(lowerPhi,lowerPhi0,a,c,LOWERDIM);
+    }
+
+    //interpolate from the lower grid
+    prolongate_parallel(phi,GRIDDIM,lowerPhi,LOWERDIM);
+
+    //smooth
+    solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM);
+}
+
+
+/**
+ * Relaxes an ODE matrix by 1 iteration of multigrid method
+ * @param phi The phi array
+ * @param phi0 The phi array from the last frame
+ * @param a The a const
+ * @param c The c const
+ * @param GRIDDIM The dimension of the phi grid
+ * @return The residual
+ */
+static inline void solver_multigrid_parallel_iterate_recurse_f_cycle(float * phi, float * phi0, float a, float c, int GRIDDIM){
+    int LOWERDIM = ((GRIDDIM - 2) / 2) + 2;
+    float * currResidual = solver_multigrid_get_current_residual(GRIDDIM);
+    float * lowerPhi = solver_multigrid_get_current_phi(LOWERDIM);
+    float * lowerPhi0 = solver_multigrid_get_current_phi0(LOWERDIM);
+    float * lowerResidual = solver_multigrid_get_current_residual(LOWERDIM);
+
+    //smooth
+    solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM);
+
+    //compute residuals
+    solver_multigrid_parallel_store_residual(phi,phi0,currResidual,a,c,GRIDDIM);
+
+    //restrict
+    restrict_parallel(currResidual,GRIDDIM,lowerPhi,lowerPhi0,LOWERDIM);
+
+    //solve next-coarsest grid
+    if(GRIDDIM <= MULTIGRID_PARALLEL_LOWEST_PARALLEL_DIM){
+        solver_multigrid_iterate_serial_recursive(lowerPhi,lowerPhi0,a,c,LOWERDIM);
+    } else {
+        solver_multigrid_parallel_iterate_recurse_v_cycle(lowerPhi,lowerPhi0,a,c,LOWERDIM);
     }
 
     //interpolate from the lower grid
@@ -286,7 +326,7 @@ static inline void solver_multigrid_parallel_iterate_recurse(float * phi, float
  */
 static inline float solver_multigrid_parallel_iterate(float * phi, float * phi0, float a, float c){
     solver_multigrid_initialization_check();
-    solver_multigrid_parallel_iterate_recurse(phi,phi0,a,c,DIM);
+    solver_multigrid_parallel_iterate_recurse_v_cycle(phi,phi0,a,c,DIM);
     return solver_multigrid_calculate_residual_norm_serial(phi,phi0,a,c,DIM);
 }
 

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

@@ -186,7 +186,7 @@ LIBRARY_API void fluid_grid2_solveProjection(
         iteration++;
     }
     if(residual > FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE || residual < -FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE){
-        printf("Projection residual didn't converge!   %f  \n",residual);
+        // printf("Projection residual didn't converge!   %f  \n",residual);
     }
 
     // for(i = 0; i < 100; i++){

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

@@ -484,7 +484,7 @@ float solver_multigrid_calculate_residual_norm_serial(float * phi, float * phi0,
         }
     }
     residual_norm = (float)sqrt(residual_norm);
-    printf("residual_norm: %lf\n",residual_norm);
+    // printf("residual_norm: %lf\n",residual_norm);
     return residual_norm;
 }