diff --git a/.vscode/settings.json b/.vscode/settings.json
index fdd3769d..9140e023 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -42,6 +42,9 @@
         "velocity.h": "c",
         "density.h": "c",
         "grid2.h": "c",
-        "math.h": "c"
+        "math.h": "c",
+        "multigrid.h": "c",
+        "gauss_seidel.h": "c",
+        "ode_utils.h": "c"
     }
 }
\ No newline at end of file
diff --git a/docs/src/progress/renderertodo.md b/docs/src/progress/renderertodo.md
index 24a81a47..7d9dc2c2 100644
--- a/docs/src/progress/renderertodo.md
+++ b/docs/src/progress/renderertodo.md
@@ -1279,6 +1279,12 @@ cellular stability work
 (12/09/2024)
 cellular stability work
 Add grid2 sim files
+grid2 functioning + lots of tests
+multigrid implementation
+
+(12/10/2024)
+Implement conjugate gradient solver
+
 
 
 # TODO
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 bf168bcc..77788bcd 100644
--- a/src/main/c/includes/fluid/sim/grid2/solver_consts.h
+++ b/src/main/c/includes/fluid/sim/grid2/solver_consts.h
@@ -5,17 +5,17 @@
 /**
  * The number of times to relax most solvers
  */
-#define FLUID_GRID2_LINEARSOLVERTIMES 20
+#define FLUID_GRID2_LINEARSOLVERTIMES 3
 
 /**
  * The number of times to relax the vector diffusion solver
  */
-#define FLUID_GRID2_VECTOR_DIFFUSE_TIMES 20
+#define FLUID_GRID2_VECTOR_DIFFUSE_TIMES 3
 
 /**
  * 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
diff --git a/src/main/c/includes/math/ode/conjugate_gradient.h b/src/main/c/includes/math/ode/conjugate_gradient.h
new file mode 100644
index 00000000..3147fe46
--- /dev/null
+++ b/src/main/c/includes/math/ode/conjugate_gradient.h
@@ -0,0 +1,23 @@
+#ifndef MATH_CONJUGATE_GRADIENT_H
+#define MATH_CONJUGATE_GRADIENT_H
+
+
+/**
+ * Iniitalizes the conjugate gradient solver with the phi values
+ * @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(float * phi, float * phi0, float a, float c);
+
+/**
+ * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient 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
+ */
+void solver_conjugate_gradient_iterate(float * phi, float * phi0, float a, float c);
+
+#endif
\ No newline at end of file
diff --git a/src/main/c/includes/math/ode/gauss_seidel.h b/src/main/c/includes/math/ode/gauss_seidel.h
new file mode 100644
index 00000000..eaf6eea0
--- /dev/null
+++ b/src/main/c/includes/math/ode/gauss_seidel.h
@@ -0,0 +1,108 @@
+
+#ifndef MATH_GAUSS_SEIDEL_H
+#define MATH_GAUSS_SEIDEL_H
+
+#include "immintrin.h"
+
+#include "fluid/queue/chunk.h"
+#include "fluid/queue/chunkmask.h"
+
+/**
+ * Gets the index into the array
+ */
+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
+ * @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_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++){
+            int n = 0;
+            //If there is any leftover, perform manual solving
+            for(i=1; 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;
+            }
+        }
+    }
+}
+
+
+
+#endif
\ No newline at end of file
diff --git a/src/main/c/includes/math/ode/multigrid.h b/src/main/c/includes/math/ode/multigrid.h
new file mode 100644
index 00000000..c53fee1d
--- /dev/null
+++ b/src/main/c/includes/math/ode/multigrid.h
@@ -0,0 +1,27 @@
+
+#ifndef MATH_SOLVER_MULTIGRID_H
+#define MATH_SOLVER_MULTIGRID_H
+
+#include "immintrin.h"
+
+#include "fluid/queue/chunk.h"
+#include "fluid/queue/chunkmask.h"
+
+/**
+ * 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
+ */
+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();
+
+
+
+#endif
\ No newline at end of file
diff --git a/src/main/c/includes/math/ode/ode_utils.h b/src/main/c/includes/math/ode/ode_utils.h
new file mode 100644
index 00000000..9f79523c
--- /dev/null
+++ b/src/main/c/includes/math/ode/ode_utils.h
@@ -0,0 +1,18 @@
+
+#ifndef MATH_ODE_UTILS_H
+#define MATH_ODE_UTILS_H
+
+
+/**
+ * Gets the index into the array
+ */
+static inline int ode_index(int x, int y, int z, int N){
+    return (x + (N * y) +  (N * N * z));
+}
+
+/**
+ * Computes the stencil of a given source array at a given position
+ */
+__m256 ode_poisson_stencil_parallel(float * source, int x, int y, int z);
+
+#endif
\ No newline at end of file
diff --git a/src/main/c/src/fluid/queue/javainterface.c b/src/main/c/src/fluid/queue/javainterface.c
index 5e40e154..f9e917f0 100644
--- a/src/main/c/src/fluid/queue/javainterface.c
+++ b/src/main/c/src/fluid/queue/javainterface.c
@@ -16,6 +16,7 @@
 #include "fluid/sim/grid2/grid2.h"
 #include "fluid/sim/simulator.h"
 #include "fluid/dispatch/dispatcher.h"
+#include "math/ode/multigrid.h"
 
 
 //defines
@@ -143,6 +144,7 @@ JNIEXPORT void JNICALL Java_electrosphere_server_fluid_simulator_FluidAccelerate
 
     //init grid2 sim
     fluid_grid2_allocate_arrays();
+    solver_multigrid_allocate();
 }
 
 /**
diff --git a/src/main/c/src/fluid/sim/grid2/velocity.c b/src/main/c/src/fluid/sim/grid2/velocity.c
index 3152f906..b8316d7e 100644
--- a/src/main/c/src/fluid/sim/grid2/velocity.c
+++ b/src/main/c/src/fluid/sim/grid2/velocity.c
@@ -8,6 +8,9 @@
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/velocity.h"
 #include "fluid/sim/grid2/utilities.h"
+#include "math/ode/gauss_seidel.h"
+#include "math/ode/multigrid.h"
+#include "math/ode/conjugate_gradient.h"
 
 #define SET_BOUND_IGNORE 0
 #define SET_BOUND_USE_NEIGHBOR 1
@@ -227,6 +230,7 @@ LIBRARY_API void fluid_grid2_solveProjection(
   float ** jrv0,
   float dt
 ){
+    float a = 1;
     float c = 6;
     int i, j, k, l, m;
     __m256 cScalar = _mm256_set1_ps(c);
@@ -234,36 +238,13 @@ LIBRARY_API void fluid_grid2_solveProjection(
 
     float * p = GET_ARR_RAW(jru0,CENTER_LOC);
     float * div = GET_ARR_RAW(jrv0,CENTER_LOC);
-    // update for each cell
-    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){
 
-                //compute Q
-                vector =                      _mm256_loadu_ps(&p[IX(i-1,j,k)]);
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i+1,j,k)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j-1,k)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j+1,k)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j,k-1)]));
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j,k+1)]));
-
-                //add A
-                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&div[IX(i,j,k)]));
-
-                //divide by 6
-                vector = _mm256_div_ps(vector,cScalar);
-                _mm256_storeu_ps(&p[IX(i,j,k)],vector);
-            }
-            //If there is any leftover, perform manual solving
-            if(i>DIM-1){
-                for(i=i-8; i < DIM-1; i++){
-                    p[IX(i,j,k)] = (div[IX(i,j,k)] + (p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c;
-                }
-            }
-        }
-    }
+    //transform u direction
+    solver_gauss_seidel_iterate_parallel(p,div,a,c,DIM);
+    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);
 }
 
 /**
diff --git a/src/main/c/src/math/ode/conjugate_gradient.c b/src/main/c/src/math/ode/conjugate_gradient.c
new file mode 100644
index 00000000..11d235d2
--- /dev/null
+++ b/src/main/c/src/math/ode/conjugate_gradient.c
@@ -0,0 +1,143 @@
+#include <stdlib.h>
+#include <math.h>
+#include <immintrin.h>
+
+#include "fluid/queue/chunk.h"
+#include "math/ode/ode_utils.h"
+
+
+/**
+ * The search direction array
+ */
+float * p = NULL;
+
+/**
+ * The residual array
+ */
+float * r = NULL;
+
+
+/**
+ * Iniitalizes the conjugate gradient solver with the phi values
+ * @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(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));
+    }
+
+    int i, j, k;
+    
+    __m256 aScalar = _mm256_set1_ps(a);
+    __m256 cScalar = _mm256_set1_ps(c);
+    __m256 f, residual, stencil;
+
+    //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 phi
+                //get values from neighbors
+                stencil =                       _mm256_loadu_ps(&phi[ode_index( i-1, j,   k,   DIM )]);
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i+1, j,   k,   DIM )]));
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i,   j-1, k,   DIM )]));
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i,   j+1, k,   DIM )]));
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i,   j,   k-1, DIM )]));
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i,   j,   k+1, DIM )]));
+                //multiply by a
+                stencil = _mm256_mul_ps(stencil,aScalar);
+                //add previous value
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index(i,j,k,DIM)]));
+                //divide by 6
+                stencil = _mm256_div_ps(stencil,cScalar);
+
+                //grab the f value (the value of phi0)
+                f = _mm256_loadu_ps(&phi0[ode_index(i,j,k,DIM)]);
+
+                //subtract stencil from f
+                residual = _mm256_sub_ps(f,stencil);
+
+                //store in residual and searchDir
+                _mm256_storeu_ps(&r[ode_index(i,j,k,DIM)],residual);
+                _mm256_storeu_ps(&p[ode_index(i,j,k,DIM)],residual);
+            }
+        }
+    }
+}
+
+/**
+ * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient 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
+ */
+void solver_conjugate_gradient_iterate(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 denominator, r0_squared, r1_squared, a_k, r_delta, beta, p0, p1, p_new, phi_old, phi1, r0, r1;
+
+    //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
+                //get values from neighbors
+                stencil =                       _mm256_loadu_ps(&p[ode_index( i-1, j,   k,   DIM )]);
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i+1, j,   k,   DIM )]));
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i,   j-1, k,   DIM )]));
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i,   j+1, k,   DIM )]));
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i,   j,   k-1, DIM )]));
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i,   j,   k+1, DIM )]));
+                //multiply by a
+                stencil = _mm256_mul_ps(stencil,aScalar);
+                //add previous value
+                stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index(i,j,k,DIM)]));
+                //divide by 6
+                stencil = _mm256_div_ps(stencil,cScalar);
+
+                //load p0
+                p0 = _mm256_loadu_ps(&p[ode_index(i,j,k,DIM)]);
+
+                //calculate ak
+                denominator = _mm256_mul_ps(stencil,p0);
+                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)]));
+                a_k = _mm256_div_ps(r0_squared,denominator);
+
+                //update phi
+                phi_old = _mm256_loadu_ps(&phi[ode_index(i,j,k,DIM)]);
+                phi1 = _mm256_mul_ps(p1,a_k);
+                phi1 = _mm256_add_ps(phi1,phi_old);
+                _mm256_storeu_ps(&phi[ode_index(i,j,k,DIM)],phi1);
+
+                //update residual
+                r_delta = _mm256_add_ps(stencil,a_k);
+                r0 = _mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]);
+                r1 = _mm256_sub_ps(r0,r_delta);
+                _mm256_storeu_ps(&r[ode_index(i,j,k,DIM)],r1);
+
+                //calculate beta
+                r1_squared = _mm256_mul_ps(r1,r1);
+                beta = _mm256_div_ps(r1_squared,r0_squared);
+
+                //update p (search dir)
+                p_new = _mm256_mul_ps(beta,p0);
+                p1 = _mm256_add_ps(r1,p_new);
+                _mm256_storeu_ps(&p[ode_index(i,j,k,DIM)],p1);
+            }
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/main/c/src/math/ode/multigrid.c b/src/main/c/src/math/ode/multigrid.c
new file mode 100644
index 00000000..3f381893
--- /dev/null
+++ b/src/main/c/src/math/ode/multigrid.c
@@ -0,0 +1,105 @@
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "fluid/queue/chunk.h"
+#include "math/ode/gauss_seidel.h"
+#include "math/ode/multigrid.h"
+
+/**
+ * Dimension of the half resolution grid
+ */
+static int halfDim = ((DIM - 2) / 2) + 2;
+
+/**
+ * Dimension of the quarter resolution grid
+ */
+static int quarterDim = ((DIM - 2) / 4) + 2;
+
+/**
+ * The half resolution grids
+ */
+static float * halfGridPhi;
+static float * halfGridPhi0;
+
+/**
+ * The quarter resolution grids
+ */
+static float * quarterGridPhi;
+static float * quarterGridPhi0;
+
+/**
+ * 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
+ */
+void solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
+
+    //
+    //half res
+    for(int x = 1; x < halfDim - 2; x++){
+        for(int y = 1; y < halfDim - 2; y++){
+            for(int z = 1; z < halfDim - 2; z++){
+                halfGridPhi[solver_gauss_seidel_get_index(x,y,z,halfDim)] = phi[solver_gauss_seidel_get_index(x*2,y*2,z*2,DIM)];
+                halfGridPhi0[solver_gauss_seidel_get_index(x,y,z,halfDim)] = phi0[solver_gauss_seidel_get_index(x*2,y*2,z*2,DIM)];
+            }
+        }
+    }
+    solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim);
+
+
+
+    //
+    //quarter res
+    //
+    //half res
+    for(int x = 1; x < quarterDim - 2; x++){
+        for(int y = 1; y < quarterDim - 2; y++){
+            for(int z = 1; z < quarterDim - 2; z++){
+                quarterGridPhi[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = halfGridPhi[solver_gauss_seidel_get_index(x*2,y*2,z*2,halfDim)];
+                quarterGridPhi0[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = halfGridPhi0[solver_gauss_seidel_get_index(x*2,y*2,z*2,halfDim)];
+            }
+        }
+    }
+    solver_gauss_seidel_iterate_parallel(quarterGridPhi,quarterGridPhi0,a,c,quarterDim);
+    for(int x = 1; x < halfDim - 2; x++){
+        for(int y = 1; y < halfDim - 2; y++){
+            for(int z = 1; z < halfDim - 2; z++){
+                halfGridPhi[solver_gauss_seidel_get_index(x,y,z,halfDim)] = quarterGridPhi[solver_gauss_seidel_get_index(x/2,y/2,z/2,quarterDim)];
+                halfGridPhi0[solver_gauss_seidel_get_index(x,y,z,halfDim)] = quarterGridPhi0[solver_gauss_seidel_get_index(x/2,y/2,z/2,quarterDim)];
+            }
+        }
+    }
+
+
+
+
+
+
+
+    //
+    //half res
+    solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim);
+    for(int x = 1; x < DIM - 2; x++){
+        for(int y = 1; y < DIM - 2; y++){
+            for(int z = 1; z < DIM - 2; z++){
+                phi[solver_gauss_seidel_get_index(x,y,z,DIM)] = halfGridPhi[solver_gauss_seidel_get_index(x/2,y/2,z/2,halfDim)];
+                phi0[solver_gauss_seidel_get_index(x,y,z,DIM)] = halfGridPhi0[solver_gauss_seidel_get_index(x/2,y/2,z/2,halfDim)];
+            }
+        }
+    }
+}
+
+
+/**
+ * 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));
+}
+
+
diff --git a/src/main/c/src/math/ode/ode_utils.c b/src/main/c/src/math/ode/ode_utils.c
new file mode 100644
index 00000000..718040a5
--- /dev/null
+++ b/src/main/c/src/math/ode/ode_utils.c
@@ -0,0 +1,23 @@
+#include "immintrin.h"
+
+#include "fluid/queue/chunk.h"
+#include "math/ode/ode_utils.h"
+
+// /**
+//  * Computes the stencil of a given source array at a given position
+//  */
+// __m256 ode_poisson_stencil_parallel(float * source, float * sourcePrev, int x, int y, int z){
+//     __m256 POISSON_STENCIL_C_SCALAR = _mm256_set1_ps(6);
+//     //get values from neighbors
+//     __m256 stencil =                _mm256_loadu_ps(&source[ode_index( x-1, y,   z,   DIM )]);
+//     stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x+1, y,   z,   DIM )]));
+//     stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x,   y-1, z,   DIM )]));
+//     stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x,   y+1, z,   DIM )]));
+//     stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x,   y,   z-1, DIM )]));
+//     stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x,   y,   z+1, DIM )]));
+//     //add previous value
+//     stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&sourcePrev[ode_index(x,y,z,DIM)]));
+//     //divide by 6
+//     stencil = _mm256_div_ps(stencil,POISSON_STENCIL_C_SCALAR);
+//     _mm256_storeu_ps(&source[ode_index(x,y,z,DIM)],stencil);
+// }
\ No newline at end of file
diff --git a/src/test/c/fluid/sim/grid2/advect_projection_tests.c b/src/test/c/fluid/sim/grid2/advect_projection_tests.c
index 2f459db4..b8eff790 100644
--- a/src/test/c/fluid/sim/grid2/advect_projection_tests.c
+++ b/src/test/c/fluid/sim/grid2/advect_projection_tests.c
@@ -11,6 +11,7 @@
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/utilities.h"
 #include "fluid/sim/grid2/velocity.h"
+#include "math/ode/multigrid.h"
 #include "../../../util/chunk_test_utils.h"
 #include "../../../util/test.h"
 
@@ -82,13 +83,87 @@ int fluid_sim_grid2_advect_projection_test1(){
     return rVal;
 }
 
+/**
+ * Testing velocity advection
+ */
+int fluid_sim_grid2_advect_projection_compute_error_over_time(){
+    printf("fluid_sim_grid2_advect_projection_compute_error_over_time\n");
+    int rVal = 0;
+    Environment * env = fluid_environment_create();
+    Chunk ** queue = NULL;
+    queue = createChunkGrid(env,3,3,3);
+
+
+
+    for(int complexity = 2; complexity < DIM-1; complexity++){
+        //actually simulate
+        for(int frameCount = 1; frameCount < 50; frameCount++){
+            //setup chunk values
+            Chunk * currentChunk = queue[0];
+            chunk_fill(currentChunk,0);
+            for(int x = 1; x < complexity; x++){
+                for(int y = 1; y < complexity; y++){
+                    for(int z = 1; z < complexity; z++){
+                        currentChunk->d[CENTER_LOC][IX(x,y,z)] = MAX_FLUID_VALUE;
+                        currentChunk->u[CENTER_LOC][IX(x,y,z)] = MAX_FLUID_VALUE;
+                    }
+                }
+            }
+            float beforeSum = chunk_queue_sum_density(queue);
+
+            //solve
+            for(int frame = 0; frame < frameCount; frame++){
+                int chunkCount = arrlen(queue);
+                for(int chunkIndex = 0; chunkIndex < 1; chunkIndex++){
+                    currentChunk = queue[chunkIndex];
+                    //advect velocity
+                    fluid_grid2_flip_arrays(currentChunk->u,currentChunk->u0);
+                    fluid_grid2_flip_arrays(currentChunk->v,currentChunk->v0);
+                    fluid_grid2_flip_arrays(currentChunk->w,currentChunk->w0);
+                    fluid_grid2_advectVectors(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_SIM_STEP);
+                    fluid_grid2_flip_arrays(currentChunk->u,currentChunk->u0);
+                    fluid_grid2_flip_arrays(currentChunk->v,currentChunk->v0);
+                    fluid_grid2_flip_arrays(currentChunk->w,currentChunk->w0);
+
+                    ////project
+                    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_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
+                    fluid_grid2_flip_arrays(currentChunk->d,currentChunk->d0);
+                    fluid_grid2_advectDensity(currentChunk->d,currentChunk->d0,currentChunk->u,currentChunk->v,currentChunk->w,FLUID_GRID2_SIM_STEP);
+                }
+            }
+            //test the result
+            float afterSum = chunk_queue_sum_density(queue);
+            if(fabs(beforeSum - afterSum) > FLUID_GRID2_PROJECTION_ERROR_MARGIN){
+                printf("%f,",fabs(beforeSum - afterSum));
+                // rVal += assertEqualsFloat(beforeSum,afterSum,"Advection changed density!  %f %f  \n");
+            }
+        }
+        printf("\n");
+    }
+
+    rVal++;
+    return rVal;
+}
+
 /**
  * Testing velocity advection
  */
 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();
 
     return rVal;
 }
\ No newline at end of file
diff --git a/src/test/c/fluid/sim/grid2/density_advection_tests.c b/src/test/c/fluid/sim/grid2/density_advection_tests.c
index 9d70d323..26110022 100644
--- a/src/test/c/fluid/sim/grid2/density_advection_tests.c
+++ b/src/test/c/fluid/sim/grid2/density_advection_tests.c
@@ -10,6 +10,7 @@
 #include "fluid/sim/grid2/density.h"
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/utilities.h"
+#include "math/ode/multigrid.h"
 #include "../../../util/chunk_test_utils.h"
 #include "../../../util/test.h"
 
@@ -227,6 +228,8 @@ 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();
diff --git a/src/test/c/fluid/sim/grid2/density_diffuse_tests.c b/src/test/c/fluid/sim/grid2/density_diffuse_tests.c
index 506baead..340abd2a 100644
--- a/src/test/c/fluid/sim/grid2/density_diffuse_tests.c
+++ b/src/test/c/fluid/sim/grid2/density_diffuse_tests.c
@@ -9,6 +9,7 @@
 #include "fluid/sim/grid2/density.h"
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/utilities.h"
+#include "math/ode/multigrid.h"
 #include "../../../util/chunk_test_utils.h"
 #include "../../../util/test.h"
 
@@ -161,6 +162,8 @@ 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();
 
diff --git a/src/test/c/fluid/sim/grid2/finalize_projection_tests.c b/src/test/c/fluid/sim/grid2/finalize_projection_tests.c
index 0d5e572f..aca05a76 100644
--- a/src/test/c/fluid/sim/grid2/finalize_projection_tests.c
+++ b/src/test/c/fluid/sim/grid2/finalize_projection_tests.c
@@ -11,6 +11,7 @@
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/utilities.h"
 #include "fluid/sim/grid2/velocity.h"
+#include "math/ode/multigrid.h"
 #include "../../../util/chunk_test_utils.h"
 #include "../../../util/test.h"
 
@@ -79,6 +80,8 @@ 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;
diff --git a/src/test/c/fluid/sim/grid2/setup_projection_tests.c b/src/test/c/fluid/sim/grid2/setup_projection_tests.c
index bc14acce..ff223095 100644
--- a/src/test/c/fluid/sim/grid2/setup_projection_tests.c
+++ b/src/test/c/fluid/sim/grid2/setup_projection_tests.c
@@ -11,6 +11,7 @@
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/utilities.h"
 #include "fluid/sim/grid2/velocity.h"
+#include "math/ode/multigrid.h"
 #include "../../../util/chunk_test_utils.h"
 #include "../../../util/test.h"
 
@@ -81,8 +82,10 @@ int fluid_sim_grid2_setup_projection_test2(){
 int fluid_sim_grid2_setup_projection_tests(int argc, char **argv){
     int rVal = 0;
 
-    rVal += fluid_sim_grid2_setup_projection_test1();
-    rVal += fluid_sim_grid2_setup_projection_test2();
+    solver_multigrid_allocate();
+
+    // rVal += fluid_sim_grid2_setup_projection_test1();
+    // rVal += fluid_sim_grid2_setup_projection_test2();
 
     return rVal;
 }
\ No newline at end of file
diff --git a/src/test/c/fluid/sim/grid2/solve_projection_tests.c b/src/test/c/fluid/sim/grid2/solve_projection_tests.c
index b3fb3dcf..c1294484 100644
--- a/src/test/c/fluid/sim/grid2/solve_projection_tests.c
+++ b/src/test/c/fluid/sim/grid2/solve_projection_tests.c
@@ -11,6 +11,7 @@
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/utilities.h"
 #include "fluid/sim/grid2/velocity.h"
+#include "math/ode/multigrid.h"
 #include "../../../util/chunk_test_utils.h"
 #include "../../../util/test.h"
 
@@ -107,6 +108,8 @@ 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;
diff --git a/src/test/c/fluid/sim/grid2/speed_tests.c b/src/test/c/fluid/sim/grid2/speed_tests.c
new file mode 100644
index 00000000..3c072a3f
--- /dev/null
+++ b/src/test/c/fluid/sim/grid2/speed_tests.c
@@ -0,0 +1,45 @@
+
+
+#define TARGET_SIM_RADIUS 3.0
+#define TARGET_SIM_DIAMETER (TARGET_SIM_RADIUS * 2) + 1
+
+/**
+ * Ie probably shouldn't be simulating the sky
+ */
+#define PERCENT_ACTIVE_IN_RADIUS 0.5
+#define CHUNKS_PER_PLAYER ((TARGET_SIM_DIAMETER * TARGET_SIM_DIAMETER * TARGET_SIM_DIAMETER) * PERCENT_ACTIVE_IN_RADIUS)
+
+/**
+ * Say 8 players are logged in
+ */
+#define PLAYERS_PER_SERVER 8.0
+
+//worst case scenario is they're all on different boats trying to reach one another
+
+#define TARGET_CHUNKS (PLAYERS_PER_SERVER * CHUNKS_PER_PLAYER)
+
+/**
+ * Number of threads on server to dedicate to fluid sim
+ */
+#define SIM_THREADS 4.0
+
+
+/**
+ * Allow 8 milli seconds per thread to simulate
+ */
+#define TOTAL_ALLOWED_TIME_IN_MILLI_SECONDS 8.0
+
+/**
+ * Allowed time per chunk
+ */
+#define TIME_PER_CHUNK (TOTAL_ALLOWED_TIME_IN_MILLI_SECONDS / (TARGET_CHUNKS / SIM_THREADS))
+
+
+
+int fluid_sim_grid2_speed_tests(){
+    int rVal = 0;
+
+    int itmePErChunk = TIME_PER_CHUNK;
+
+    return rVal;
+}
\ No newline at end of file
diff --git a/src/test/c/fluid/sim/grid2/velocity_advection_tests.c b/src/test/c/fluid/sim/grid2/velocity_advection_tests.c
index a75678a5..b852eae9 100644
--- a/src/test/c/fluid/sim/grid2/velocity_advection_tests.c
+++ b/src/test/c/fluid/sim/grid2/velocity_advection_tests.c
@@ -11,6 +11,7 @@
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/utilities.h"
 #include "fluid/sim/grid2/velocity.h"
+#include "math/ode/multigrid.h"
 #include "../../../util/chunk_test_utils.h"
 #include "../../../util/test.h"
 
@@ -182,6 +183,8 @@ 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();
diff --git a/src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c b/src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c
index 31f3d148..9993d40c 100644
--- a/src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c
+++ b/src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c
@@ -9,6 +9,7 @@
 #include "fluid/sim/grid2/solver_consts.h"
 #include "fluid/sim/grid2/utilities.h"
 #include "fluid/sim/grid2/velocity.h"
+#include "math/ode/multigrid.h"
 #include "../../../util/chunk_test_utils.h"
 #include "../../../util/test.h"
 
@@ -149,6 +150,8 @@ 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;