diff --git a/.vscode/settings.json b/.vscode/settings.json
index e612e93c..e2324420 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -46,6 +46,7 @@
         "multigrid.h": "c",
         "gauss_seidel.h": "c",
         "ode_utils.h": "c",
-        "util.h": "c"
+        "util.h": "c",
+        "conjugate_gradient.h": "c"
     }
 }
\ No newline at end of file
diff --git a/src/main/c/includes/fluid/sim/grid2/solver_consts.h b/src/main/c/includes/fluid/sim/grid2/solver_consts.h
index dd1d7cde..20ac3225 100644
--- a/src/main/c/includes/fluid/sim/grid2/solver_consts.h
+++ b/src/main/c/includes/fluid/sim/grid2/solver_consts.h
@@ -5,7 +5,27 @@
 /**
  * The number of times to relax most solvers
  */
-#define FLUID_GRID2_LINEARSOLVERTIMES 20
+#define FLUID_GRID2_LINEARSOLVERTIMES 10
+
+/**
+ * The number of times to relax most solvers
+ */
+#define FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS 20
+
+/**
+ * Tolerance to target for multigrid precomputing in projection
+ */
+#define FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE 0.005f
+
+/**
+ * The number of times to relax most solvers
+ */
+#define FLUID_GRID2_SOLVER_CG_MAX_ITERATIONS 10
+
+/**
+ * Tolerance to target for conjugate gradient precomputing in projection
+ */
+#define FLUID_GRID2_SOLVER_CG_TOLERANCE 0.001f
 
 /**
  * Width of a single grid cell
diff --git a/src/main/c/includes/math/ode/conjugate_gradient.h b/src/main/c/includes/math/ode/conjugate_gradient.h
index a9844639..1ea491bf 100644
--- a/src/main/c/includes/math/ode/conjugate_gradient.h
+++ b/src/main/c/includes/math/ode/conjugate_gradient.h
@@ -21,6 +21,15 @@ int solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c);
  */
 void solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float a, float c);
 
+/**
+ * Initializes the conjugate gradient solver
+ * @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_init_serial(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
@@ -28,6 +37,6 @@ void solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float
  * @param a The a const
  * @param c The c const
  */
-void solver_conjugate_gradient_solve_serial(float * phi, float * phi0, float a, float c);
+float solver_conjugate_gradient_iterate_serial(float * phi, float * phi0, float a, float c);
 
 #endif
\ No newline at end of file
diff --git a/src/main/c/src/fluid/sim/grid2/velocity.c b/src/main/c/src/fluid/sim/grid2/velocity.c
index 44c57fbd..2647ee3e 100644
--- a/src/main/c/src/fluid/sim/grid2/velocity.c
+++ b/src/main/c/src/fluid/sim/grid2/velocity.c
@@ -181,21 +181,31 @@ LIBRARY_API void fluid_grid2_solveProjection(
     //perform iteration of v cycle multigrid method
     chunk->projectionResidual = 1;
     chunk->projectionIterations = 0;
-    while(chunk->projectionIterations < FLUID_GRID2_LINEARSOLVERTIMES && (chunk->projectionResidual > FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE)){
+    while(chunk->projectionIterations < FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS && (chunk->projectionResidual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){
         chunk->projectionResidual = solver_multigrid_parallel_iterate(p,div,a,c);
         fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
         chunk->projectionIterations++;
     }
-    if(chunk->projectionResidual > FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE){
-        // printf("Projection residual didn't converge!   %f  \n",residual);
-    }
+    // if(chunk->projectionResidual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE){
+    //     // printf("Projection residual didn't converge!   %f  \n",residual);
+    // }
 
     // for(i = 0; i < 100; i++){
     //     solver_gauss_seidel_iterate_parallel(p,div,a,c,DIM);
     //     fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
     // }
 
-    // solver_conjugate_gradient_solve_serial(p,div,a,c);
+    //init CG solver
+    // solver_conjugate_gradient_init_serial(p,div,a,c);
+    // //solve with CG
+    // while(chunk->projectionIterations < FLUID_GRID2_SOLVER_CG_MAX_ITERATIONS && (chunk->projectionResidual > FLUID_GRID2_SOLVER_CG_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_SOLVER_CG_TOLERANCE)){
+    //     chunk->projectionResidual = solver_conjugate_gradient_iterate_serial(p,div,a,c);;
+    //     fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
+    //     chunk->projectionIterations++;
+    // }
+    // if(chunk->projectionResidual > FLUID_GRID2_SOLVER_CG_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_SOLVER_CG_TOLERANCE){
+    //     // printf("Projection residual didn't converge!   %f  \n",residual);
+    // }
 
     //finest grain iteration with conjugate gradient method
     // int shouldSolve = solver_conjugate_gradient_init(p,div,a,c);
diff --git a/src/main/c/src/math/ode/conjugate_gradient.c b/src/main/c/src/math/ode/conjugate_gradient.c
index 09002bce..8d46705e 100644
--- a/src/main/c/src/math/ode/conjugate_gradient.c
+++ b/src/main/c/src/math/ode/conjugate_gradient.c
@@ -190,8 +190,87 @@ void solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float
  * @param phi0 The phi array from the last frame
  * @param a The a const
  * @param c The c const
+ * @return The residual
  */
-void solver_conjugate_gradient_solve_serial(float * phi, float * phi0, float a, float c){
+float solver_conjugate_gradient_iterate_serial(float * phi, float * phi0, float a, float c){
+    int i, j, k;
+    float convergence, denominator;
+    float laplacian, alpha, r_new_dot, beta;
+    //solve Ap
+    for(k=1; k<DIM-1; k++){
+        for(j=1; j<DIM-1; j++){
+            for(i = 1; i < DIM-1; i++){
+                laplacian = 
+                (
+                    6 * p[ode_index (i   , j   , k   ,DIM)] +
+                    - (
+                        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)]
+                    )
+                );
+                A[ode_index(i,j,k,DIM)] = laplacian;
+            }
+        }
+    }
+    convergence = 0;
+    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);
+    }
+     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)];
+            }
+        }
+    }
+
+    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)];
+            }
+        }
+    }
+
+    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)];
+            }
+        }
+    }
+    return (float)sqrt(convergence);
+}
+
+
+/**
+ * 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
+ * @return The residual
+ */
+void solver_conjugate_gradient_init_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));
@@ -212,78 +291,5 @@ void solver_conjugate_gradient_solve_serial(float * phi, float * phi0, float a,
             }
         }
     }
-
-    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);
-    }
-
 }