From 72f9fcdcbaefeb78a43341c5d12ee13ed42d65ec Mon Sep 17 00:00:00 2001 From: austin Date: Tue, 10 Dec 2024 19:12:47 -0500 Subject: [PATCH] grid2 work --- .vscode/settings.json | 3 +- .../includes/fluid/sim/grid2/solver_consts.h | 11 +- .../c/includes/math/ode/conjugate_gradient.h | 14 +- src/main/c/src/fluid/sim/grid2/density.c | 15 ++ src/main/c/src/fluid/sim/grid2/fluidsim.c | 6 +- src/main/c/src/fluid/sim/grid2/velocity.c | 28 ++- src/main/c/src/math/ode/conjugate_gradient.c | 196 +++++++++++++++--- src/main/c/src/math/ode/multigrid.c | 8 + src/test/c/CMakeLists.txt | 2 + .../fluid/sim/grid2/advect_projection_tests.c | 12 +- .../sim/grid2/finalize_projection_tests.c | 6 +- .../fluid/sim/grid2/setup_projection_tests.c | 18 +- .../fluid/sim/grid2/solve_projection_tests.c | 6 +- 13 files changed, 259 insertions(+), 66 deletions(-) diff --git a/.vscode/settings.json b/.vscode/settings.json index 9140e023..e612e93c 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -45,6 +45,7 @@ "math.h": "c", "multigrid.h": "c", "gauss_seidel.h": "c", - "ode_utils.h": "c" + "ode_utils.h": "c", + "util.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 77788bcd..08cb95d2 100644 --- a/src/main/c/includes/fluid/sim/grid2/solver_consts.h +++ b/src/main/c/includes/fluid/sim/grid2/solver_consts.h @@ -5,23 +5,28 @@ /** * The number of times to relax most solvers */ -#define FLUID_GRID2_LINEARSOLVERTIMES 3 +#define FLUID_GRID2_LINEARSOLVERTIMES 2 /** * The number of times to relax the vector diffusion solver */ -#define FLUID_GRID2_VECTOR_DIFFUSE_TIMES 3 +#define FLUID_GRID2_VECTOR_DIFFUSE_TIMES 2 /** * Width of a single grid cell */ -#define FLUID_GRID2_H (1.0/(DIM-2)) +#define FLUID_GRID2_H (1.0/DIM) /** * Timestep to simulate by */ #define FLUID_GRID2_SIM_STEP 0.01f +/** + * Const to multiply the advection stages by to offset numeric instability + */ +#define FLUID_GRID2_CORRECTION_CONST 1.0005f + /** * Really small value used for something */ diff --git a/src/main/c/includes/math/ode/conjugate_gradient.h b/src/main/c/includes/math/ode/conjugate_gradient.h index 3147fe46..a9844639 100644 --- a/src/main/c/includes/math/ode/conjugate_gradient.h +++ b/src/main/c/includes/math/ode/conjugate_gradient.h @@ -8,8 +8,9 @@ * @param phi0 The phi array from the last frame * @param a The a const * @param c The c const + * @return 1 if the system has already been solved, 0 otherwise */ -void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c); +int solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c); /** * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method @@ -18,6 +19,15 @@ void 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_iterate(float * phi, float * phi0, float a, float c); +void solver_conjugate_gradient_iterate_parallel(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 + * @param phi0 The phi array from the last frame + * @param a The a const + * @param c The c const + */ +void solver_conjugate_gradient_solve_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/density.c b/src/main/c/src/fluid/sim/grid2/density.c index 07a1e932..dc3c241e 100644 --- a/src/main/c/src/fluid/sim/grid2/density.c +++ b/src/main/c/src/fluid/sim/grid2/density.c @@ -8,6 +8,7 @@ #include "fluid/queue/chunk.h" #include "fluid/sim/grid2/solver_consts.h" #include "fluid/sim/grid2/utilities.h" +#include "math/ode/multigrid.h" /** @@ -54,6 +55,10 @@ LIBRARY_API void fluid_grid2_solveDiffuseDensity( __m256 aScalar = _mm256_set1_ps(a); __m256 cScalar = _mm256_set1_ps(c); + // for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){ + // solver_multigrid_iterate(x,x0,a,c); + // } + //transform u direction for(k=1; k DIM-1 || j0 > DIM-1 || k0 > DIM-1 || + i1 < 0 || j1 < 0 || k1 < 0 || + i1 > DIM-1 || j1 > DIM-1 || k1 > DIM-1 + ){ + printf("advect dens: %d %d %d %d %d %d --- %f %f %f\n", i0, j0, k0, i1, j1, k1, x, y, z); + fflush(stdout); + } center_d[IX(i,j,k)] = s0*( diff --git a/src/main/c/src/fluid/sim/grid2/fluidsim.c b/src/main/c/src/fluid/sim/grid2/fluidsim.c index 362c29be..7bc7848a 100644 --- a/src/main/c/src/fluid/sim/grid2/fluidsim.c +++ b/src/main/c/src/fluid/sim/grid2/fluidsim.c @@ -98,10 +98,8 @@ void fluid_grid2_simulate( //these should have just been mirrored in the above // //Perform main projection solver - for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){ - fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,timestep); - fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0); - } + fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,timestep); + fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0); //samples u,v,w,u0 //sets u,v,w //Finalize projection diff --git a/src/main/c/src/fluid/sim/grid2/velocity.c b/src/main/c/src/fluid/sim/grid2/velocity.c index b8316d7e..60eb9a22 100644 --- a/src/main/c/src/fluid/sim/grid2/velocity.c +++ b/src/main/c/src/fluid/sim/grid2/velocity.c @@ -239,12 +239,18 @@ LIBRARY_API void fluid_grid2_solveProjection( float * p = GET_ARR_RAW(jru0,CENTER_LOC); float * div = GET_ARR_RAW(jrv0,CENTER_LOC); - //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); + //perform iteration of v cycle multigrid method + for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){ + solver_multigrid_iterate(p,div,a,c); + } + + // solver_conjugate_gradient_solve_serial(p,div,a,c); + + //finest grain iteration with conjugate gradient method + // int shouldSolve = solver_conjugate_gradient_init(p,div,a,c); + // if(shouldSolve < 1){ + // solver_conjugate_gradient_iterate(p,div,a,c); + // } } /** @@ -419,6 +425,16 @@ void fluid_grid2_advect_velocity(int b, float ** jrd, float ** jrd0, float * u, u1 = z-k0; u0 = 1-u1; + + if( + i0 < 0 || j0 < 0 || k0 < 0 || + i0 > DIM-1 || j0 > DIM-1 || k0 > DIM-1 || + i1 < 0 || j1 < 0 || k1 < 0 || + i1 > DIM-1 || j1 > DIM-1 || k1 > DIM-1 + ){ + printf("advect vel: %d %d %d %d %d %d --- %f %f %f\n", i0, j0, k0, i1, j1, k1, x, y, z); + fflush(stdout); + } d[IX(i,j,k)] = diff --git a/src/main/c/src/math/ode/conjugate_gradient.c b/src/main/c/src/math/ode/conjugate_gradient.c index 11d235d2..09002bce 100644 --- a/src/main/c/src/math/ode/conjugate_gradient.c +++ b/src/main/c/src/math/ode/conjugate_gradient.c @@ -3,18 +3,33 @@ #include #include "fluid/queue/chunk.h" +#include "fluid/sim/grid2/solver_consts.h" #include "math/ode/ode_utils.h" +/** + * Used for preventing divide by 0's + */ +static float CONJUGATE_GRADIENT_EPSILON = 1e-6; + +static float CONJUGATE_GRADIENT_CONVERGENCE_THRESHOLD = 0.0001f; + /** * The search direction array */ -float * p = NULL; +static float * p = NULL; + +/** + * Maxmum frames to iterate for + */ +static int max_frames = 100; /** * The residual array */ -float * r = NULL; +static float * r = NULL; + +static float * A = NULL; /** @@ -23,14 +38,18 @@ float * r = NULL; * @param phi0 The phi array from the last frame * @param a The a const * @param c The c const + * @return 1 if the system has already been solved, 0 otherwise */ -void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c){ +int 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)); } + if(A == NULL){ + A = (float *)calloc(1,DIM*DIM*DIM*sizeof(float)); + } int i, j, k; @@ -41,7 +60,17 @@ void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c) //iniitalize the r (residual) and p (search direction) arrays for(k=1; k CONJUGATE_GRADIENT_CONVERGENCE_THRESHOLD && framecount < max_frames){ + //solve Ap + for(k=1; k -CONJUGATE_GRADIENT_EPSILON){ + printf("Divide by 0! %f \n", denominator); + fflush(stdout); + break; + } + float alpha = convergence / denominator; + for(k=1; k= 100 || convergence > CONJUGATE_GRADIENT_CONVERGENCE_THRESHOLD){ + printf("Failed to converge! %f \n", convergence); + fflush(stdout); + } + +} + diff --git a/src/main/c/src/math/ode/multigrid.c b/src/main/c/src/math/ode/multigrid.c index 3f381893..97023394 100644 --- a/src/main/c/src/math/ode/multigrid.c +++ b/src/main/c/src/math/ode/multigrid.c @@ -36,6 +36,10 @@ static float * quarterGridPhi0; */ void solver_multigrid_iterate(float * phi, float * phi0, float a, float c){ + // + //gauss-seidel at highest res + solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,DIM); + // //half res for(int x = 1; x < halfDim - 2; x++){ @@ -89,6 +93,10 @@ void solver_multigrid_iterate(float * phi, float * phi0, float a, float c){ } } } + + // + //gauss-seidel at highest res + solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,DIM); } diff --git a/src/test/c/CMakeLists.txt b/src/test/c/CMakeLists.txt index 0cf62235..6a770d81 100644 --- a/src/test/c/CMakeLists.txt +++ b/src/test/c/CMakeLists.txt @@ -58,6 +58,8 @@ foreach (TEST_FILE ${TEST_SOURCES}) endif() endforeach () +target_compile_options(test_runner PRIVATE -m64 -mavx -mavx2) + # make test runner depend on library add_dependencies(test_runner StormEngine) 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 b8eff790..48f1ea2b 100644 --- a/src/test/c/fluid/sim/grid2/advect_projection_tests.c +++ b/src/test/c/fluid/sim/grid2/advect_projection_tests.c @@ -62,10 +62,8 @@ int fluid_sim_grid2_advect_projection_test1(){ 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_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 @@ -129,10 +127,8 @@ int fluid_sim_grid2_advect_projection_compute_error_over_time(){ 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_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 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 aca05a76..ebe2c429 100644 --- a/src/test/c/fluid/sim/grid2/finalize_projection_tests.c +++ b/src/test/c/fluid/sim/grid2/finalize_projection_tests.c @@ -36,10 +36,8 @@ int fluid_sim_grid2_finalize_projection_test1(){ Chunk * currentChunk = queue[0]; currentChunk->u[CENTER_LOC][IX(3,3,3)] = 1.0f; fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP); - 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_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP); + fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0); //finalize fluid_grid2_finalizeProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP); 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 ff223095..b0509d33 100644 --- a/src/test/c/fluid/sim/grid2/setup_projection_tests.c +++ b/src/test/c/fluid/sim/grid2/setup_projection_tests.c @@ -41,10 +41,10 @@ int fluid_sim_grid2_setup_projection_test1(){ fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP); //test the result - rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],-0.5f / DIM,"Divergence of the vector field at 3,3,3 should be -0.5/DIM! %f %f \n"); - rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(3,3,3)],0,"Divergence of the vector field at 3,3,3 should be 0! %f %f \n"); - rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],0.5f / DIM,"Divergence of the vector field at 4,3,3 should be 0.5/DIM! %f %f \n"); - rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(5,3,3)],0,"Divergence of the vector field at 5,3,3 should be 0! %f %f \n"); + rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],-0.5f * FLUID_GRID2_H,"Divergence of the vector field at 3,3,3 should be -0.5 * h! actual: %f expected: %f \n"); + rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(3,3,3)],0,"Divergence of the vector field at 3,3,3 should be 0! actual: %f expected: %f \n"); + rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],0.5f * FLUID_GRID2_H,"Divergence of the vector field at 4,3,3 should be 0.5 * h! actual: %f expected: %f \n"); + rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(5,3,3)],0,"Divergence of the vector field at 5,3,3 should be 0! actual: %f expected: %f \n"); return rVal; } @@ -53,7 +53,7 @@ int fluid_sim_grid2_setup_projection_test1(){ * Testing velocity advection */ int fluid_sim_grid2_setup_projection_test2(){ - printf("fluid_sim_grid2_setup_projection_test1\n"); + printf("fluid_sim_grid2_setup_projection_test2\n"); int rVal = 0; Environment * env = fluid_environment_create(); Chunk ** queue = NULL; @@ -69,9 +69,9 @@ int fluid_sim_grid2_setup_projection_test2(){ fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP); //test the result - rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],0.5f / DIM,"Divergence of the vector field at 3,3,3 should be 0.5/DIM! %f %f \n"); + rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],0.5f * FLUID_GRID2_H,"Divergence of the vector field at 3,3,3 should be 0.5 * h! actual: %f expected: %f \n"); rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(3,3,3)],0,"Divergence of the vector field at 3,3,3 should be 0! %f %f \n"); - rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],-0.5f / DIM,"Divergence of the vector field at 4,3,3 should be -0.5/DIM! %f %f \n"); + rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],-0.5f * FLUID_GRID2_H,"Divergence of the vector field at 4,3,3 should be -0.5 * h! actual: %f expected: %f \n"); return rVal; } @@ -84,8 +84,8 @@ int fluid_sim_grid2_setup_projection_tests(int argc, char **argv){ solver_multigrid_allocate(); - // rVal += fluid_sim_grid2_setup_projection_test1(); - // rVal += fluid_sim_grid2_setup_projection_test2(); + 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 c1294484..a97f1cca 100644 --- a/src/test/c/fluid/sim/grid2/solve_projection_tests.c +++ b/src/test/c/fluid/sim/grid2/solve_projection_tests.c @@ -38,10 +38,8 @@ int fluid_sim_grid2_solve_projection_test1(){ fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP); //actually solve - 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_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP); + fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0); //test the result float expected, actual;