fully recursive multigrid solver

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

@@ -14,57 +14,6 @@ static inline int solver_gauss_seidel_get_index(int x, int y, int z, int N){
     return (x + (N * y) +  (N * N * z));
 }
 
-/**
- * Relaxes an ODE matrix by 1 iteration of gauss seidel parallelized
- * @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 dimensions of the grid
- */
-static inline void solver_gauss_seidel_iterate_parallel(float * phi, float * phi0, float a, float c, int gridDim){
-    int i, j, k, l, m;
-    
-    __m256 aScalar = _mm256_set1_ps(a);
-    __m256 cScalar = _mm256_set1_ps(c);
-
-    //transform u direction
-    for(k=1; k<gridDim-1; k++){
-        for(j=1; j<gridDim-1; j++){
-            int n = 0;
-            //solve as much as possible vectorized
-            for(i = 1; i < gridDim-1; i=i+8){
-                __m256 vector = _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,gridDim)]);
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,gridDim)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,gridDim)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,gridDim)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,gridDim)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,gridDim)]));
-                vector = _mm256_mul_ps(vector,aScalar);
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,gridDim)]));
-                vector = _mm256_div_ps(vector,cScalar);
-                _mm256_storeu_ps(&phi[solver_gauss_seidel_get_index(i,j,k,gridDim)],vector);
-            }
-            //If there is any leftover, perform manual solving
-            if(i>gridDim-1){
-                for(i=i-8; i < gridDim-1; i++){
-                    phi[solver_gauss_seidel_get_index(i,j,k,gridDim)] =
-                    (
-                        phi0[solver_gauss_seidel_get_index(i,j,k,gridDim)] + 
-                        a * (
-                            phi[solver_gauss_seidel_get_index(i-1,j,k,gridDim)]+
-                            phi[solver_gauss_seidel_get_index(i+1,j,k,gridDim)]+
-                            phi[solver_gauss_seidel_get_index(i,j-1,k,gridDim)]+
-                            phi[solver_gauss_seidel_get_index(i,j+1,k,gridDim)]+
-                            phi[solver_gauss_seidel_get_index(i,j,k-1,gridDim)]+
-                            phi[solver_gauss_seidel_get_index(i,j,k+1,gridDim)]
-                        )
-                    ) / c;
-                }
-            }
-        }
-    }
-}
 
 /**
  * Relaxes an ODE matrix by 1 iteration of gauss seidel serially
@@ -77,9 +26,6 @@ static inline void solver_gauss_seidel_iterate_parallel(float * phi, float * phi
 static inline void solver_gauss_seidel_iterate_serial(float * phi, float * phi0, float a, float c, int gridDim){
     int i, j, k, l, m;
     
-    __m256 aScalar = _mm256_set1_ps(a);
-    __m256 cScalar = _mm256_set1_ps(c);
-
     //transform u direction
     for(k=1; k<gridDim-1; k++){
         for(j=1; j<gridDim-1; j++){
@@ -103,6 +49,99 @@ static inline void solver_gauss_seidel_iterate_serial(float * phi, float * phi0,
     }
 }
 
+/**
+ * Relaxes an ODE matrix by 1 iteration of gauss seidel parallelized
+ * @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 dimensions of the grid
+ */
+static inline void solver_gauss_seidel_iterate_parallel(float * phi, float * phi0, float a, float c, int gridDim){
+    int i, j, k, l, m;
+    if(gridDim >= 10){
+        __m256 aScalar = _mm256_set1_ps(a);
+        __m256 cScalar = _mm256_set1_ps(c);
+        //transform u direction
+        for(k=1; k<gridDim-1; k++){
+            for(j=1; j<gridDim-1; j++){
+                int n = 0;
+                //solve as much as possible vectorized
+                for(i = 1; i < gridDim-1; i=i+8){
+                    __m256 vector = _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,gridDim)]);
+                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,gridDim)]));
+                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,gridDim)]));
+                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,gridDim)]));
+                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,gridDim)]));
+                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,gridDim)]));
+                    vector = _mm256_mul_ps(vector,aScalar);
+                    vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,gridDim)]));
+                    vector = _mm256_div_ps(vector,cScalar);
+                    _mm256_storeu_ps(&phi[solver_gauss_seidel_get_index(i,j,k,gridDim)],vector);
+                }
+                //If there is any leftover, perform manual solving
+                if(i>gridDim-1){
+                    for(i=i-8; i < gridDim-1; i++){
+                        phi[solver_gauss_seidel_get_index(i,j,k,gridDim)] =
+                        (
+                            phi0[solver_gauss_seidel_get_index(i,j,k,gridDim)] + 
+                            a * (
+                                phi[solver_gauss_seidel_get_index(i-1,j,k,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i+1,j,k,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i,j-1,k,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i,j+1,k,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i,j,k-1,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i,j,k+1,gridDim)]
+                            )
+                        ) / c;
+                    }
+                }
+            }
+        }
+    } else if(gridDim >= 6){
+        __m128 aScalar = _mm_set1_ps(a);
+        __m128 cScalar = _mm_set1_ps(c);
+        //transform u direction
+        for(k=1; k<gridDim-1; k++){
+            for(j=1; j<gridDim-1; j++){
+                int n = 0;
+                //solve as much as possible vectorized
+                for(i = 1; i < gridDim-1; i=i+8){
+                    __m128 vector = _mm_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,gridDim)]);
+                    vector = _mm_add_ps(vector,_mm_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,gridDim)]));
+                    vector = _mm_add_ps(vector,_mm_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,gridDim)]));
+                    vector = _mm_add_ps(vector,_mm_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,gridDim)]));
+                    vector = _mm_add_ps(vector,_mm_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,gridDim)]));
+                    vector = _mm_add_ps(vector,_mm_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,gridDim)]));
+                    vector = _mm_mul_ps(vector,aScalar);
+                    vector = _mm_add_ps(vector,_mm_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,gridDim)]));
+                    vector = _mm_div_ps(vector,cScalar);
+                    _mm_storeu_ps(&phi[solver_gauss_seidel_get_index(i,j,k,gridDim)],vector);
+                }
+                //If there is any leftover, perform manual solving
+                if(i>gridDim-1){
+                    for(i=i-8; i < gridDim-1; i++){
+                        phi[solver_gauss_seidel_get_index(i,j,k,gridDim)] =
+                        (
+                            phi0[solver_gauss_seidel_get_index(i,j,k,gridDim)] + 
+                            a * (
+                                phi[solver_gauss_seidel_get_index(i-1,j,k,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i+1,j,k,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i,j-1,k,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i,j+1,k,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i,j,k-1,gridDim)]+
+                                phi[solver_gauss_seidel_get_index(i,j,k+1,gridDim)]
+                            )
+                        ) / c;
+                    }
+                }
+            }
+        }
+    } else {
+        solver_gauss_seidel_iterate_serial(phi,phi0,a,c,gridDim);
+    }
+}
+
 
 
 #endif

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

@@ -15,7 +15,17 @@
  * @param c The c const
  * @return The residual
  */
-float solver_multigrid_iterate(float * phi, float * phi0, float a, float c);
+float solver_multigrid_iterate_serial(float * phi, float * phi0, float a, float c);
+
+/**
+ * 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
+ * @return The residual
+ */
+float solver_multigrid_iterate_parallel(float * phi, float * phi0, float a, float c);
 
 
 

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

@@ -181,7 +181,7 @@ LIBRARY_API void fluid_grid2_solveProjection(
     float residual = 1;
     int iteration = 0;
     while(iteration < FLUID_GRID2_LINEARSOLVERTIMES && (residual > FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE || residual < -FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE)){
-        residual = solver_multigrid_iterate(p,div,a,c);
+        residual = solver_multigrid_iterate_serial(p,div,a,c);
         fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
         iteration++;
     }

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

@@ -18,6 +18,7 @@ static int halfDim = ((DIM - 2) / 2) + 2;
  * Dimension of the quarter resolution grid
  */
 static int quarterDim = ((DIM - 2) / 4) + 2;
+static int LOWEST_DIM = ((DIM - 2) / 4) + 2;
 
 /**
  * The full resolution grids
@@ -38,10 +39,21 @@ static float * quarterGridPhi = NULL;
 static float * quarterGridPhi0 = NULL;
 static float * quarterGridResidual = NULL;
 
+/**
+ * The eighth resolution grids
+ */
+static float * eighthGridPhi = NULL;
+static float * eighthGridPhi0 = NULL;
+static float * eighthGridResidual = NULL;
+
 
 float solver_multigrid_calculate_residual_norm_serial(float * phi, float * phi0, float a, float c, int GRIDDIM);
 float solver_multigrid_calculate_residual_parallel(float * phi, float * phi0, float a, float c, int GRIDDIM);
 float solver_multigrid_store_residual_serial(float * phi, float * phi0, float * residualGrid, float a, float c, int GRIDDIM);
+float * get_current_phi(int dim);
+float * get_current_phi0(int dim);
+float * get_current_residual(int dim);
+void initialization_check();
 
 /**
  * Relaxes an ODE matrix by 1 iteration of multigrid method
@@ -51,40 +63,171 @@ float solver_multigrid_store_residual_serial(float * phi, float * phi0, float *
  * @param c The c const
  * @return The residual
  */
-float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
+float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float a, float c, int GRIDDIM){
+    initialization_check();
+    float residual;
+    printf("Current dim: %d \n",GRIDDIM);
+    fflush(stdout);
 
+    int LOWERDIM = ((GRIDDIM - 2) / 2) + 2;
+    float * currResidual = get_current_residual(GRIDDIM);
+    float * lowerPhi = get_current_phi(LOWERDIM);
+    float * lowerPhi0 = get_current_phi0(LOWERDIM);
+    float * lowerResidual = get_current_residual(LOWERDIM);
+
+    printf("Solved arrays: %p  %p  %p  %p  \n",currResidual,lowerPhi, lowerPhi0, lowerResidual);
+    fflush(stdout);
+
+    printf("smooth\n");
+    fflush(stdout);
+    //
+    //gauss-seidel at highest res
+    solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM);
+    if(GRIDDIM == DIM){
+        fluid_grid2_set_bounds(0,phi);
+    }
+    //compute residuals
+    residual = solver_multigrid_store_residual_serial(phi,phi0,currResidual,a,c,GRIDDIM);
+
+    printf("restrict\n");
+    fflush(stdout);
+    //restrict
+    //(current operator is injection -- inject r^2 from this grid at phi0 of the smaller grid)
+    for(int x = 0; x < LOWERDIM; x++){
+        for(int y = 0; y < LOWERDIM; y++){
+            for(int z = 0; z < LOWERDIM; z++){
+                lowerPhi[solver_gauss_seidel_get_index(x,y,z,LOWERDIM)] = 0;
+                lowerPhi0[solver_gauss_seidel_get_index(x,y,z,LOWERDIM)] = 0;
+            }
+        }
+    }
+    //populate grid
+    for(int x = 1; x < LOWERDIM - 1; x++){
+        for(int y = 1; y < LOWERDIM - 1; y++){
+            for(int z = 1; z < LOWERDIM - 1; z++){
+                //direct transfer operator (faster, lower accuracy)
+                lowerPhi0[solver_gauss_seidel_get_index(x,y,z,LOWERDIM)] = currResidual[solver_gauss_seidel_get_index(x*2,y*2,z*2,GRIDDIM)];
+            }
+        }
+    }
+
+    printf("recurse\n");
+    fflush(stdout);
+    //solve next-coarsest grid
+    if(GRIDDIM <= LOWEST_DIM){
+        //smooth
+        solver_gauss_seidel_iterate_parallel(quarterGridPhi,quarterGridPhi0,a,c,quarterDim);
+        //compute residual
+        residual = solver_multigrid_store_residual_serial(quarterGridPhi,quarterGridPhi0,quarterGridResidual,a,c,quarterDim);
+    } else {
+        printf("solve\n");
+        fflush(stdout);
+        float solution = 
+            (
+                phi0[solver_gauss_seidel_get_index(1,1,1,GRIDDIM)] + 
+                phi0[solver_gauss_seidel_get_index(2,1,1,GRIDDIM)] + 
+                phi0[solver_gauss_seidel_get_index(1,2,1,GRIDDIM)] + 
+                phi0[solver_gauss_seidel_get_index(1,1,2,GRIDDIM)] +
+                phi0[solver_gauss_seidel_get_index(1,2,2,GRIDDIM)] + 
+                phi0[solver_gauss_seidel_get_index(2,1,2,GRIDDIM)] + 
+                phi0[solver_gauss_seidel_get_index(2,2,1,GRIDDIM)] + 
+                phi0[solver_gauss_seidel_get_index(2,2,2,GRIDDIM)]
+            ) / 8.0f
+        ;
+        
+        printf("interpolate solution\n");
+        fflush(stdout);
+        //interpolate from the lower grid
+        for(int x = 1; x < GRIDDIM - 1; x++){
+            for(int y = 1; y < GRIDDIM - 1; y++){
+                for(int z = 1; z < GRIDDIM - 1; z++){
+                    //direct transfer operator (faster, lower accuracy)
+                    phi[solver_gauss_seidel_get_index(x,y,z,GRIDDIM)] = phi[solver_gauss_seidel_get_index(x,y,z,GRIDDIM)] + solution;
+                }
+            }
+        }
+    }
+
+    printf("interpolate\n");
+    fflush(stdout);
+    //interpolate from the lower grid
+    for(int x = 1; x < GRIDDIM - 1; x++){
+        for(int y = 1; y < GRIDDIM - 1; y++){
+            for(int z = 1; z < GRIDDIM - 1; z++){
+                //direct transfer operator (faster, lower accuracy)
+                phi[solver_gauss_seidel_get_index(x,y,z,GRIDDIM)] = 
+                phi[solver_gauss_seidel_get_index(x,y,z,GRIDDIM)] + 
+                lowerPhi[solver_gauss_seidel_get_index(  x/2+0,  y/2+0,  z/2+0  ,LOWERDIM)]
+                ;
+
+                //interpolation operator (slower, better accuracy)
+                // phi[solver_gauss_seidel_get_index(x,y,z,DIM)] = 
+                // phi[solver_gauss_seidel_get_index(x,y,z,DIM)] + 
+                // (
+                //     halfGridPhi[solver_gauss_seidel_get_index(  x/2+0,  y/2+0,  z/2+0  ,halfDim)] +
+                //     halfGridPhi[solver_gauss_seidel_get_index(  x/2+0,  y/2+0,  z/2+1  ,halfDim)] +
+                //     halfGridPhi[solver_gauss_seidel_get_index(  x/2+0,  y/2+1,  z/2+0  ,halfDim)] +
+                //     halfGridPhi[solver_gauss_seidel_get_index(  x/2+0,  y/2+1,  z/2+1  ,halfDim)] +
+                //     halfGridPhi[solver_gauss_seidel_get_index(  x/2+1,  y/2+0,  z/2+0  ,halfDim)] +
+                //     halfGridPhi[solver_gauss_seidel_get_index(  x/2+1,  y/2+0,  z/2+1  ,halfDim)] +
+                //     halfGridPhi[solver_gauss_seidel_get_index(  x/2+1,  y/2+1,  z/2+0  ,halfDim)] +
+                //     halfGridPhi[solver_gauss_seidel_get_index(  x/2+1,  y/2+1,  z/2+1  ,halfDim)]
+                // )
+                // ;
+            }
+        }
+    }
+
+
+
+
+
+
+    printf("smooth\n");
+    fflush(stdout);
+    //
     // full res
-    if(fullGridResidual == NULL){
+    //
-        fullGridResidual = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
+    //
+    //smooth
+    solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM);
+    if(GRIDDIM == DIM){
+        fluid_grid2_set_bounds(0,phi);
     }
 
-    // half res
+    
-    if(halfGridPhi == NULL){
+    return solver_multigrid_calculate_residual_norm_serial(phi,phi0,a,c,GRIDDIM);
-        halfGridPhi = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
+}
-    }
-    if(halfGridPhi0 == NULL){
-        halfGridPhi0 = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
-    }
-    if(halfGridResidual == NULL){
-        halfGridResidual = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
-    }
 
-    // quarter res
-    if(quarterGridPhi == NULL){
-        quarterGridPhi = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
-    }
-    if(quarterGridPhi0 == NULL){
-        quarterGridPhi0 = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
-    }
-    if(quarterGridResidual == NULL){
-        quarterGridResidual = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
-    }
 
+/**
+ * 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
+ * @return The residual
+ */
+float solver_multigrid_iterate_serial(float * phi, float * phi0, float a, float c){
+    return solver_multigrid_iterate_serial_recursive(phi,phi0,a,c,DIM);
+}
+
+
+/**
+ * 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
+ * @return The residual
+ */
+float solver_multigrid_iterate_parallel(float * phi, float * phi0, float a, float c){
+    initialization_check();
     float residual;
 
     //
     //gauss-seidel at highest res
-    solver_gauss_seidel_iterate_serial(phi,phi0,a,c,DIM);
+    solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,DIM);
     fluid_grid2_set_bounds(0,phi);
     //compute residuals
     residual = solver_multigrid_store_residual_serial(phi,phi0,fullGridResidual,a,c,DIM);
@@ -117,7 +260,7 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
     //half res
     //
     //smooth
-    solver_gauss_seidel_iterate_serial(halfGridPhi,halfGridPhi0,a,c,halfDim);
+    solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim);
     //compute residual
     residual = solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim);
 
@@ -151,7 +294,7 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
     }
 
     //smooth
-    solver_gauss_seidel_iterate_serial(quarterGridPhi,quarterGridPhi0,a,c,quarterDim);
+    solver_gauss_seidel_iterate_parallel(quarterGridPhi,quarterGridPhi0,a,c,quarterDim);
     //compute residual
     residual = solver_multigrid_store_residual_serial(quarterGridPhi,quarterGridPhi0,quarterGridResidual,a,c,quarterDim);
 
@@ -193,7 +336,7 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
     //
 
     //smooth
-    solver_gauss_seidel_iterate_serial(halfGridPhi,halfGridPhi0,a,c,halfDim);
+    solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim);
     //compute residual
     residual = solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim);
 
@@ -236,7 +379,7 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
     //
     //
     //smooth
-    solver_gauss_seidel_iterate_serial(phi,phi0,a,c,DIM);
+    solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,DIM);
     fluid_grid2_set_bounds(0,phi);
 
     
@@ -245,6 +388,103 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
 
 
 
+/**
+ * Gets the phi to use for the current level of the multigrid solver
+ */
+float * get_current_phi(int dim){
+    if(dim == ((DIM-2)/2) + 2){
+        return halfGridPhi;
+    } else if(dim == ((DIM-2)/4) + 2){
+        return quarterGridPhi;
+    } else if(dim == ((DIM-2)/8) + 2){
+        return eighthGridPhi;
+    } else {
+        printf("[get_current_phi] Invalid dim detected! %d\n",dim);
+        fflush(stdout);
+        return NULL;
+    }
+}
+
+/**
+ * Gets the phi0 to use for the current level of the multigrid solver
+ */
+float * get_current_phi0(int dim){
+    if(dim == ((DIM-2)/2) + 2){
+        return halfGridPhi0;
+    } else if(dim == ((DIM-2)/4) + 2){
+        return quarterGridPhi0;
+    } else if(dim == ((DIM-2)/8) + 2){
+        return eighthGridPhi0;
+    } else {
+        printf("[get_current_phi0] Invalid dim detected! %d\n",dim);
+        fflush(stdout);
+        return NULL;
+    }
+}
+
+/**
+ * Gets the residual to use for the current level of the multigrid solver
+ */
+float * get_current_residual(int dim){
+    if(dim == 18){
+        return fullGridResidual;
+    } else if(dim == ((DIM-2)/2) + 2){
+        return halfGridResidual;
+    } else if(dim == ((DIM-2)/4) + 2){
+        return quarterGridResidual;
+    } else if(dim == ((DIM-2)/8) + 2){
+        return eighthGridResidual;
+    } else {
+        printf("[get_current_residual] Invalid dim detected! %d\n",dim);
+        fflush(stdout);
+        return NULL;
+    }
+}
+
+/**
+ * Verifies that all grids are allocated
+ */
+void initialization_check(){
+    // full res
+    if(fullGridResidual == NULL){
+        fullGridResidual = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
+    }
+
+    // half res
+    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(halfGridResidual == NULL){
+        halfGridResidual = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
+    }
+
+    // quarter res
+    if(quarterGridPhi == NULL){
+        quarterGridPhi = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
+    }
+    if(quarterGridPhi0 == NULL){
+        quarterGridPhi0 = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
+    }
+    if(quarterGridResidual == NULL){
+        quarterGridResidual = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
+    }
+
+    // quarter res
+    int eighthResDim = ((DIM-2)/8) + 2;
+    if(eighthGridPhi == NULL){
+        eighthGridPhi = (float *)calloc(1,eighthResDim * eighthResDim * eighthResDim * sizeof(float));
+    }
+    if(eighthGridPhi0 == NULL){
+        eighthGridPhi0 = (float *)calloc(1,eighthResDim * eighthResDim * eighthResDim * sizeof(float));
+    }
+    if(eighthGridResidual == NULL){
+        eighthGridResidual = (float *)calloc(1,eighthResDim * eighthResDim * eighthResDim * sizeof(float));
+    }
+}
+
 
 /**
  * Calculates the residual of the grid