multigrid working with diffusion solvers

.vscode/settings.json

@@ -48,6 +48,7 @@
         "ode_utils.h": "c",
         "util.h": "c",
         "conjugate_gradient.h": "c",
-        "flux.h": "c"
+        "flux.h": "c",
+        "diffusion_ode.h": "c"
     }
 }

src/main/c/includes/fluid/env/environment.h

@@ -4,6 +4,7 @@
 #include <jni.h>
 #include "public.h"
 #include "fluid/queue/chunk.h"
+#include "math/ode/diffusion_ode.h"
 
 /**
  * The List lookup table
@@ -73,6 +74,7 @@ typedef struct {
 typedef struct {
     //density data
     double densityTotal;
+    double densityMaintenance;
     double densityDiffuse;
     double densityAdvect;
 
@@ -115,6 +117,11 @@ typedef struct {
     float * fluid_grid2_neighborArr_bounds;
     float * fluid_grid2_neighborArr_divergenceCache;
     float * fluid_grid2_neighborArr_scalarCache;
+    
+    /**
+     * Data for computing diffusion ODEs
+     */
+    OdeDiffuseData diffuseData;
 } FluidGrid2State;
 
 /**

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

@@ -7,6 +7,11 @@
  */
 #define FLUID_GRID2_LINEARSOLVERTIMES 10
 
+/**
+ * Convergence threshold for density diffusion
+ */
+#define FLUID_GRID2_DENSITY_DIFFUSE_THRESHOLD 0.001f
+
 /**
  * The number of times to relax most solvers
  */

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

@@ -1,7 +1,15 @@
 #ifndef MATH_CONJUGATE_GRADIENT_H
 #define MATH_CONJUGATE_GRADIENT_H
 
+#include "math/ode/ode.h"
 
+
+
+//
+//
+//   NAVIER STOKES SPECIFIC
+//
+//
 /**
  * Iniitalizes the conjugate gradient solver with the phi values
  * @param phi The phi array
@@ -19,7 +27,7 @@ int 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_init_serial(float * phi, float * phi0, float a, float c);
+void solver_conjugate_gradient_init_navier_stokes_serial(float * phi, float * phi0, float a, float c);
 
 /**
  * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method serially
@@ -29,7 +37,7 @@ void solver_conjugate_gradient_init_serial(float * phi, float * phi0, float a, f
  * @param c The c const
  * @return The residual
  */
-float solver_conjugate_gradient_iterate_serial(float * phi, float * phi0, float a, float c);
+float solver_conjugate_gradient_iterate_navier_stokes_serial(float * phi, float * phi0, float a, float c);
 
 /**
  * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method
@@ -41,4 +49,36 @@ float solver_conjugate_gradient_iterate_serial(float * phi, float * phi0, float
  */
 float solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float a, float c);
 
+
+
+//
+//
+//    GENERIC ODES
+//
+//
+
+
+/**
+ * Computes the stencil for the conjugate gradient solver
+ */
+typedef float (* ode_cg_search_direction_stencil)(float * phi, int x, int y, int z, OdeData * data);
+
+/**
+ * Initializes the conjugate gradient solver
+ * @param phi The phi array
+ * @param phi0 The phi array from the last frame
+ */
+void solver_conjugate_gradient_init_serial(float * phi, float * phi0);
+
+
+/**
+ * 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 stencil_func The stencil to compute the search direction
+ * @param odeData The ode data
+ * @return The residual
+ */
+float solver_conjugate_gradient_iterate_serial(float * phi, float * phi0, ode_cg_search_direction_stencil stencil_func, OdeData * odeData);
+
 #endif

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

@@ -0,0 +1,27 @@
+#ifndef MATH_ODE_DIFFUSION_ODE_H
+#define MATH_ODE_DIFFUSION_ODE_H
+
+#include "math/ode/ode.h"
+#include "fluid/env/utilities.h"
+#include "fluid/queue/chunk.h"
+#include "fluid/sim/grid2/solver_consts.h"
+
+/**
+ * Data for computing the diffusion ode
+ */
+typedef struct {
+    /**
+     * Simulation timestep
+     */
+    float dt;
+} OdeDiffuseData;
+
+
+/**
+ * Computes the residual of a given position in a diffusion ode
+ */
+float ode_diffusion_cg_stencil(float * phi, int x, int y, int z, OdeData * data);
+
+
+
+#endif

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

@@ -191,7 +191,8 @@ static inline void solver_multigrid_parallel_store_residual(float * phi, float *
         return;
     }
     __m256 laplacian;
-    __m256 constVec = _mm256_set1_ps(6);
+    __m256 constVecA = _mm256_set1_ps(a);
+    __m256 constVecC = _mm256_set1_ps(c);
     //calculate residual
     int i, j, k;
     for(k=1; k<GRIDDIM-1; k++){
@@ -201,23 +202,26 @@ static inline void solver_multigrid_parallel_store_residual(float * phi, float *
                 _mm256_sub_ps(
                     _mm256_mul_ps(
                         _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)]),
-                        constVec
+                        constVecC
                     ),
-                    _mm256_add_ps(
+                    _mm256_mul_ps(
                         _mm256_add_ps(
                             _mm256_add_ps(
-                                _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,GRIDDIM)]),
+                                _mm256_add_ps(
-                                _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,GRIDDIM)])
+                                    _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,GRIDDIM)]),
+                                    _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,GRIDDIM)])
+                                ),
+                                _mm256_add_ps(
+                                    _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,GRIDDIM)]),
+                                    _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,GRIDDIM)])
+                                )
                             ),
                             _mm256_add_ps(
-                                _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,GRIDDIM)]),
+                                _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,GRIDDIM)]),
-                                _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,GRIDDIM)])
+                                _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,GRIDDIM)])
                             )
                         ),
-                        _mm256_add_ps(
+                        constVecA
-                            _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,GRIDDIM)]),
-                            _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,GRIDDIM)])
-                        )
                     )
                 );
                 _mm256_storeu_ps(

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

@@ -0,0 +1,21 @@
+#ifndef MATH_ODE_H
+#define MATH_ODE_H
+
+
+/**
+ * Data for computing the ode (ie could hold timestep for instance)
+ */
+typedef void * OdeData;
+
+/**
+ * Computes the residual of a given position in an ode
+ */
+typedef float (* ode_approximate_stencil)(float * phi, float * phi0, int x, int y, int z, OdeData * data);
+
+/**
+ * Computes the residual of a given position in an ode
+ */
+typedef float (* ode_residual_stencil)(float * phi, float * phi0, int x, int y, int z, OdeData * data);
+
+
+#endif

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

@@ -68,21 +68,23 @@ LIBRARY_API void fluid_grid2_solveDiffuseDensity(
     float * x = GET_ARR_RAW(d,CENTER_LOC);
     float * x0 = GET_ARR_RAW(d0,CENTER_LOC);
 
-    // float residual = 1;
+    //about ~40% faster than gauss seidel
-    // int iterations = 0;
+    float residual = 1;
-    // while(iterations < FLUID_GRID2_LINEARSOLVERTIMES && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){
+    int iterations = 0;
-    //     residual = solver_multigrid_parallel_iterate(x,x0,a,c);
+    while(iterations < FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){
-    //     fluid_grid2_set_bounds(BOUND_SET_DENSITY_PHI,x);
+        residual = solver_multigrid_parallel_iterate(x,x0,a,c);
-    //     iterations++;
-    // }
-    
-    for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
-        //iterate
-        solver_gauss_seidel_iterate_parallel(x,x0,a,c,DIM);
-
-        //set bounds
         fluid_grid2_set_bounds(environment,BOUND_SET_DENSITY_PHI,x);
+        iterations++;
     }
+    
+    // //about ~40% slower than multigrid
+    // for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+    //     //iterate
+    //     solver_gauss_seidel_iterate_parallel(x,x0,a,c,DIM);
+
+    //     //set bounds
+    //     fluid_grid2_set_bounds(environment,BOUND_SET_DENSITY_PHI,x);
+    // }
 }
 
 /**

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

@@ -39,12 +39,13 @@ LIBRARY_API void fluid_grid2_simulate(
     Chunk ** chunks = passedInChunks;
     double start, end, perMilli;
 
+    //update ODE solver data
+    environment->state.grid2.diffuseData.dt = timestep;
+
 
     gettimeofday(&tv,NULL);
     start = 1000000.0 * tv.tv_sec + tv.tv_usec;
-    //
+    //maintenance
-    //Velocity step
-    //
     for(int i = 0; i < numChunks; i++){
         Chunk * currentChunk = chunks[i];
         //update the bounds arrays
@@ -69,6 +70,16 @@ LIBRARY_API void fluid_grid2_simulate(
         fluid_grid2_flip_arrays(currentChunk->u,currentChunk->u0);
         fluid_grid2_flip_arrays(currentChunk->v,currentChunk->v0);
         fluid_grid2_flip_arrays(currentChunk->w,currentChunk->w0);
+    }
+
+    gettimeofday(&tv,NULL);
+    start = 1000000.0 * tv.tv_sec + tv.tv_usec;
+
+    //diffuse velocity
+    for(int i = 0; i < numChunks; i++){
+        Chunk * currentChunk = chunks[i];
+        //update the bounds arrays
+        fluid_grid2_rewrite_bounds(environment,currentChunk);
 
         //solve vector diffusion
         fluid_grid2_solveVectorDiffuse(
@@ -82,19 +93,21 @@ LIBRARY_API void fluid_grid2_simulate(
             timestep
         );
         
-    // }
+    }
 
-    // //time tracking
+    //time tracking
-    // gettimeofday(&tv,NULL);
+    gettimeofday(&tv,NULL);
-    // end = 1000000.0 * tv.tv_sec + tv.tv_usec;
+    end = 1000000.0 * tv.tv_sec + tv.tv_usec;
-    // perMilli = (end - start) / 1000.0f;
+    perMilli = (end - start) / 1000.0f;
-    // environment->state.timeTracking.velocityDiffuse = perMilli;
+    environment->state.timeTracking.velocityDiffuse = perMilli;
-    // start = end;
+    start = end;
 
-    // for(int i = 0; i < numChunks; i++){
+
-    //     Chunk * currentChunk = chunks[i];
+    //project
-    //     //update the bounds arrays
+    for(int i = 0; i < numChunks; i++){
-    //     fluid_grid2_rewrite_bounds(currentChunk);
+        Chunk * currentChunk = chunks[i];
+        //update the bounds arrays
+        fluid_grid2_rewrite_bounds(environment,currentChunk);
 
         // setup projection
         fluid_grid2_setupProjection(
@@ -120,37 +133,38 @@ LIBRARY_API void fluid_grid2_simulate(
         fluid_grid2_flip_arrays(currentChunk->v,currentChunk->v0);
         fluid_grid2_flip_arrays(currentChunk->w,currentChunk->w0);
 
-    // }
+    }
 
-    // //time tracking
+    //time tracking
-    // gettimeofday(&tv,NULL);
+    gettimeofday(&tv,NULL);
-    // end = 1000000.0 * tv.tv_sec + tv.tv_usec;
+    end = 1000000.0 * tv.tv_sec + tv.tv_usec;
-    // perMilli = (end - start) / 1000.0f;
+    perMilli = (end - start) / 1000.0f;
-    // environment->state.timeTracking.velocityProject = perMilli;
+    environment->state.timeTracking.velocityProject = perMilli;
-    // start = end;
+    start = end;
 
-    // for(int i = 0; i < numChunks; i++){
+    //advect
-    //     Chunk * currentChunk = chunks[i];
+    for(int i = 0; i < numChunks; i++){
-    //     Chunk * currentChunk = chunks[i];
+        Chunk * currentChunk = chunks[i];
-    //     //update the bounds arrays
+        //update the bounds arrays
-    //     fluid_grid2_rewrite_bounds(currentChunk);
+        fluid_grid2_rewrite_bounds(environment,currentChunk);
     
         // advect
         fluid_grid2_advectVectors(environment,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,timestep);
 
-    // }
+    }
 
-    // //time tracking
+    //time tracking
-    // gettimeofday(&tv,NULL);
+    gettimeofday(&tv,NULL);
-    // end = 1000000.0 * tv.tv_sec + tv.tv_usec;
+    end = 1000000.0 * tv.tv_sec + tv.tv_usec;
-    // perMilli = (end - start) / 1000.0f;
+    perMilli = (end - start) / 1000.0f;
-    // environment->state.timeTracking.velocityAdvect = perMilli;
+    environment->state.timeTracking.velocityAdvect = perMilli;
-    // start = end;
+    start = end;
 
-    // for(int i = 0; i < numChunks; i++){
+    //project again
-    //     Chunk * currentChunk = chunks[i];
+    for(int i = 0; i < numChunks; i++){
-    //     //update the bounds arrays
+        Chunk * currentChunk = chunks[i];
-    //     fluid_grid2_rewrite_bounds(currentChunk);
+        //update the bounds arrays
+        fluid_grid2_rewrite_bounds(environment,currentChunk);
 
         //setup projection
         fluid_grid2_setupProjection(environment,currentChunk,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,timestep);
@@ -185,32 +199,58 @@ LIBRARY_API void fluid_grid2_simulate(
     environment->state.newDensity = 0;
     for(int i = 0; i < numChunks; i++){
         Chunk * currentChunk = chunks[i];
-        //update the bounds arrays
-        fluid_grid2_rewrite_bounds(environment, currentChunk);
         //add density
         fluid_grid2_addDensity(environment,currentChunk->d,currentChunk->d0,timestep);
         environment->state.existingDensity = environment->state.existingDensity + fluid_grid2_calculateSum(currentChunk->d);
         //swap all density arrays
         fluid_grid2_flip_arrays(currentChunk->d,currentChunk->d0);
+    }
+
+    //time tracking
+    gettimeofday(&tv,NULL);
+    end = 1000000.0 * tv.tv_sec + tv.tv_usec;
+    perMilli = (end - start) / 1000.0f;
+    environment->state.timeTracking.densityMaintenance = perMilli;
+    start = end;
+
+
+    //solve density diffusion
+    for(int i = 0; i < numChunks; i++){
+        Chunk * currentChunk = chunks[i];
+        //update the bounds arrays
+        // fluid_grid2_rewrite_bounds(environment, currentChunk); //33% more time than just diffusion step
         //diffuse density
         fluid_grid2_solveDiffuseDensity(environment,currentChunk->d,currentChunk->d0,timestep);
+    }
+
+    //time tracking
+    gettimeofday(&tv,NULL);
+    end = 1000000.0 * tv.tv_sec + tv.tv_usec;
+    perMilli = (end - start) / 1000.0f;
+    environment->state.timeTracking.densityDiffuse = perMilli;
+    start = end;
+
+
+    //flip arrays
+    for(int i = 0; i < numChunks; i++){
+        Chunk * currentChunk = chunks[i];
         //swap all density arrays
         fluid_grid2_flip_arrays(currentChunk->d,currentChunk->d0);
+    }
 
-        
+    //time tracking
-    // }
+    gettimeofday(&tv,NULL);
+    end = 1000000.0 * tv.tv_sec + tv.tv_usec;
+    perMilli = (end - start) / 1000.0f;
+    environment->state.timeTracking.densityMaintenance = environment->state.timeTracking.densityMaintenance + perMilli;
+    start = end;
 
-    // //time tracking
-    // gettimeofday(&tv,NULL);
-    // end = 1000000.0 * tv.tv_sec + tv.tv_usec;
-    // perMilli = (end - start) / 1000.0f;
-    // environment->state.timeTracking.densityDiffuse = perMilli;
-    // start = end;
     
-    // for(int i = 0; i < numChunks; i++){
+    //advect
-    //     Chunk * currentChunk = chunks[i];
+    for(int i = 0; i < numChunks; i++){
-    //     //update the bounds arrays
+        Chunk * currentChunk = chunks[i];
-    //     fluid_grid2_rewrite_bounds(environment, currentChunk);
+        //update the bounds arrays
+        fluid_grid2_rewrite_bounds(environment, currentChunk);
 
 
         //advect density

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

@@ -11,6 +11,7 @@
 #include "math/ode/gauss_seidel.h"
 #include "math/ode/multigrid_parallel.h"
 #include "math/ode/conjugate_gradient.h"
+#include "math/ode/diffusion_ode.h"
 #include "util/matrix.h"
 
 #define SET_BOUND_IGNORE 0
@@ -58,21 +59,87 @@ LIBRARY_API void fluid_grid2_solveVectorDiffuse(
     float * v0 = GET_ARR_RAW(jrv0,CENTER_LOC);
     float * w0 = GET_ARR_RAW(jrw0,CENTER_LOC);
     
-    for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
+    // //about ~30% faster
-        //transform u direction
+    // for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
-        solver_gauss_seidel_iterate_parallel(u,u0,a,c,DIM);
+    //     //transform u direction
+    //     solver_gauss_seidel_iterate_parallel(u,u0,a,c,DIM);
 
-        //transform v direction
+    //     //transform v direction
-        solver_gauss_seidel_iterate_parallel(v,v0,a,c,DIM);
+    //     solver_gauss_seidel_iterate_parallel(v,v0,a,c,DIM);
 
-        //transform w direction
+    //     //transform w direction
-        solver_gauss_seidel_iterate_parallel(w,w0,a,c,DIM);
+    //     solver_gauss_seidel_iterate_parallel(w,w0,a,c,DIM);
 
-        //set bounds
+    //     //set bounds
+    //     fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_U,u);
+    //     fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_V,v);
+    //     fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_W,w);
+    // }
+
+    float residual;
+    int iterations;
+
+    residual = 1;
+    iterations = 0;
+    while(iterations < FLUID_GRID2_LINEARSOLVERTIMES && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){
+        residual = solver_multigrid_parallel_iterate(u,u0,a,c);
         fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_U,u);
-        fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_V,v);
+        iterations++;
-        fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_W,w);
     }
+
+    residual = 1;
+    iterations = 0;
+    while(iterations < FLUID_GRID2_LINEARSOLVERTIMES && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){
+        residual = solver_multigrid_parallel_iterate(v,v0,a,c);
+        fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_V,v);
+        iterations++;
+    }
+
+    residual = 1;
+    iterations = 0;
+    while(iterations < FLUID_GRID2_LINEARSOLVERTIMES && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){
+        residual = solver_multigrid_parallel_iterate(w,w0,a,c);
+        fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_W,w);
+        iterations++;
+    }
+
+
+
+
+
+    // //init CG solver
+    // solver_conjugate_gradient_init_serial(u,u0);
+    // residual = 1;
+    // iterations = 0;
+    // //solve with CG
+    // while(iterations < FLUID_GRID2_SOLVER_CG_MAX_ITERATIONS && (residual > FLUID_GRID2_SOLVER_CG_TOLERANCE || residual < -FLUID_GRID2_SOLVER_CG_TOLERANCE)){
+    //     residual = solver_conjugate_gradient_iterate_serial(u,u0,ode_diffusion_cg_stencil, (OdeData *)&(environment->state.grid2.diffuseData));
+    //     fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_U,u);
+    //     iterations++;
+    // }
+
+    // //init CG solver
+    // solver_conjugate_gradient_init_serial(v,v0);
+    // residual = 1;
+    // iterations = 0;
+    // //solve with CG
+    // while(iterations < FLUID_GRID2_SOLVER_CG_MAX_ITERATIONS && (residual > FLUID_GRID2_SOLVER_CG_TOLERANCE || residual < -FLUID_GRID2_SOLVER_CG_TOLERANCE)){
+    //     residual = solver_conjugate_gradient_iterate_parallel(v,v0,a,c);
+    //     residual = solver_conjugate_gradient_iterate_serial(v,v0,ode_diffusion_cg_stencil, (OdeData *)&(environment->state.grid2.diffuseData));
+    //     fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_V,v);
+    //     iterations++;
+    // }
+
+    // //init CG solver
+    // solver_conjugate_gradient_init_serial(w,w0);
+    // residual = 1;
+    // iterations = 0;
+    // //solve with CG
+    // while(iterations < FLUID_GRID2_SOLVER_CG_MAX_ITERATIONS && (residual > FLUID_GRID2_SOLVER_CG_TOLERANCE || residual < -FLUID_GRID2_SOLVER_CG_TOLERANCE)){
+    //     residual = solver_conjugate_gradient_iterate_serial(w,w0,ode_diffusion_cg_stencil, (OdeData *)&(environment->state.grid2.diffuseData));
+    //     fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_W,w);
+    //     iterations++;
+    // }
 }
 
 /**
@@ -201,7 +268,7 @@ LIBRARY_API void fluid_grid2_solveProjection(
     // }
 
     //init CG solver
-    solver_conjugate_gradient_init_serial(p,div,a,c);
+    solver_conjugate_gradient_init_serial(p,div);
     chunk->projectionIterations = 0;
     //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)){
@@ -210,7 +277,7 @@ LIBRARY_API void fluid_grid2_solveProjection(
         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",chunk->projectionResidual);
+        // printf("Projection residual didn't converge!   %f  \n",chunk->projectionResidual);
     }
 
     //store scalar potential in cache

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

@@ -4,8 +4,10 @@
 
 #include "fluid/queue/chunk.h"
 #include "fluid/sim/grid2/solver_consts.h"
+#include "math/ode/conjugate_gradient.h"
 #include "math/ode/ode_utils.h"
 #include "math/ode/gauss_seidel.h"
+#include "math/ode/ode.h"
 
 
 /**
@@ -330,7 +332,7 @@ float solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, floa
  * @param c The c const
  * @return The residual
  */
-float solver_conjugate_gradient_iterate_serial(float * phi, float * phi0, float a, float c){
+float solver_conjugate_gradient_iterate_navier_stokes_serial(float * phi, float * phi0, float a, float c){
     int i, j, k;
     float convergence, denominator;
     float laplacian, alpha, r_new_dot, beta;
@@ -398,9 +400,8 @@ float solver_conjugate_gradient_iterate_serial(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_init_serial(float * phi, float * phi0, float a, float c){
+void solver_conjugate_gradient_init_navier_stokes_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));
@@ -423,3 +424,118 @@ void solver_conjugate_gradient_init_serial(float * phi, float * phi0, float a, f
     }
 }
 
+/**
+ * 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
+ */
+void solver_conjugate_gradient_init_serial(float * phi, float * phi0){
+    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)];
+            }
+        }
+    }
+}
+
+
+
+/**
+ * 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 approximation_func The function to approximate the ode
+ * @param residual_func The function to compute the residual
+ * @param odeData The ode data
+ * @return The residual
+ */
+float solver_conjugate_gradient_iterate_serial(float * phi, float * phi0, ode_cg_search_direction_stencil stencil_func, OdeData * odeData){
+    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 = stencil_func(phi,i,j,k,odeData);
+                A[ode_index(i,j,k,DIM)] = laplacian;
+                // if(fabs(laplacian) > 1000000){
+                //     printf("%f\n",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(fabs(p[ode_index(i,j,k,DIM)] * A[ode_index(i,j,k,DIM)]) > 1000){
+                //     printf("convergence: %f       denominator: %f      \n",
+                //         (r[ode_index(i,j,k,DIM)] * r[ode_index(i,j,k,DIM)]),
+                //         (p[ode_index(i,j,k,DIM)] * A[ode_index(i,j,k,DIM)])
+                //     );
+                //     printf("A: %f      \n",
+                //        A[ode_index(i,j,k,DIM)]
+                //     );
+                //     printf("r: %f      \n",
+                //        r[ode_index(i,j,k,DIM)]
+                //     );
+                //     printf("p: %f      \n",
+                //        p[ode_index(i,j,k,DIM)]
+                //     );
+                //     printf("\n");
+                //     fflush(stdout);
+                // }
+            }
+        }
+    }
+
+    //have hit the desired level of convergence
+    if(convergence < CONJUGATE_GRADIENT_EPSILON && convergence > -CONJUGATE_GRADIENT_EPSILON){
+        return 0.0f;
+    }
+    if(denominator < CONJUGATE_GRADIENT_EPSILON && denominator > -CONJUGATE_GRADIENT_EPSILON){
+        printf("Divide by 0! %f \n", denominator);
+        printf("Convergence:  %f  \n",convergence);
+        fflush(stdout);
+    }
+    alpha = convergence / denominator;
+    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++){
+                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 = 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);
+}
+

src/main/c/src/math/ode/diffusion_ode.c

@@ -0,0 +1,20 @@
+
+#include "math/ode/diffusion_ode.h"
+
+/**
+ * Computes the residual of a given position in a diffusion ode
+ */
+float ode_diffusion_cg_stencil(float * phi, int x, int y, int z, OdeData * data){
+    OdeDiffuseData * diffuseData = (OdeDiffuseData *)data;
+    float a = diffuseData->dt*FLUID_GRID2_VISCOSITY_CONSTANT/(FLUID_GRID2_H*FLUID_GRID2_H);
+    float c = 1+6*a;
+    return
+        6 * phi[IX(x,y,z)] -
+        1 * (
+            phi[IX(x+1,y,z)] + phi[IX(x-1,y,z)] + 
+            phi[IX(x,y+1,z)] + phi[IX(x,y-1,z)] + 
+            phi[IX(x,y,z+1)] + phi[IX(x,y,z-1)]
+        )
+    ;
+}
+

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

@@ -58,6 +58,11 @@
  */
 #define FLUID_GRID2_PROJECTION_ERROR_MARGIN 0.00001f
 
+/**
+ * Target number of fluid frames/second
+ */
+#define TARGET_FPS 60
+
 /**
  * Used for storing timings
  */
@@ -73,7 +78,7 @@ int fluid_sim_grid2_speed_test1(){
     int rVal = 0;
     Environment * env = fluid_environment_create();
     Chunk ** queue = NULL;
-    queue = createChunkGrid(env,3,3,3);
+    queue = createChunkGrid(env,TARGET_SIM_DIAMETER,TARGET_SIM_DIAMETER,TARGET_SIM_DIAMETER);
     int chunkCount = arrlen(queue);
 
 
@@ -117,6 +122,7 @@ int fluid_sim_grid2_speed_test1(){
         printf("Density time (milli): %f  \n",env->state.timeTracking.densityTotal);
         printf(" - Advect (milli): %f  \n",env->state.timeTracking.densityAdvect);
         printf(" - Diffuse (milli): %f  \n",env->state.timeTracking.densityDiffuse);
+        printf(" - Maintenance (milli): %f  \n",env->state.timeTracking.densityMaintenance);
         printf("\n");
         rVal++;
     }
@@ -132,7 +138,7 @@ int fluid_sim_grid2_speed_test2(){
     int rVal = 0;
     Environment * env = fluid_environment_create();
     Chunk ** queue = NULL;
-    queue = createChunkGrid(env,3,3,3);
+    queue = createChunkGrid(env,TARGET_SIM_DIAMETER,TARGET_SIM_DIAMETER,TARGET_SIM_DIAMETER);
     int chunkCount = arrlen(queue);
 
 
@@ -144,7 +150,7 @@ int fluid_sim_grid2_speed_test2(){
     float beforeSum = chunk_queue_sum_density(queue);
 
     //actually simulate
-    int frameCount = 50;
+    int frameCount = TARGET_FPS;
     double frameTimeAccumulator = 0;
     for(int frame = 0; frame < frameCount; frame++){
         //get time at start
@@ -176,6 +182,7 @@ int fluid_sim_grid2_speed_test2(){
         printf("Density time (milli): %f  \n",env->state.timeTracking.densityTotal);
         printf(" - Advect (milli): %f  \n",env->state.timeTracking.densityAdvect);
         printf(" - Diffuse (milli): %f  \n",env->state.timeTracking.densityDiffuse);
+        printf(" - Maintenance (milli): %f  \n",env->state.timeTracking.densityMaintenance);
         printf("\n");
         rVal++;
     }
@@ -191,7 +198,7 @@ int fluid_sim_grid2_speed_test3(){
     int rVal = 0;
     Environment * env = fluid_environment_create();
     Chunk ** queue = NULL;
-    queue = createChunkGrid(env,5,5,5);
+    queue = createChunkGrid(env,TARGET_SIM_DIAMETER,TARGET_SIM_DIAMETER,TARGET_SIM_DIAMETER);
     int chunkCount = arrlen(queue);
 
 
@@ -203,7 +210,7 @@ int fluid_sim_grid2_speed_test3(){
     float beforeSum = chunk_queue_sum_density(queue);
 
     //actually simulate
-    int frameCount = 50;
+    int frameCount = TARGET_FPS;
     double frameTimeAccumulator = 0;
     for(int frame = 0; frame < frameCount; frame++){
         //get time at start
@@ -235,6 +242,7 @@ int fluid_sim_grid2_speed_test3(){
         printf("Density time (milli): %f  \n",env->state.timeTracking.densityTotal);
         printf(" - Advect (milli): %f  \n",env->state.timeTracking.densityAdvect);
         printf(" - Diffuse (milli): %f  \n",env->state.timeTracking.densityDiffuse);
+        printf(" - Maintenance (milli): %f  \n",env->state.timeTracking.densityMaintenance);
         printf("\n");
         rVal++;
     }
@@ -248,9 +256,9 @@ int fluid_sim_grid2_speed_test3(){
 int fluid_sim_grid2_speed_tests(){
     int rVal = 0;
 
-// rVal += fluid_sim_grid2_speed_test1();
+    rVal += fluid_sim_grid2_speed_test1();
-// rVal += fluid_sim_grid2_speed_test2();
+    rVal += fluid_sim_grid2_speed_test2();
-// rVal += fluid_sim_grid2_speed_test3();
+    rVal += fluid_sim_grid2_speed_test3();
 
     return rVal;
 }