diff --git a/.vscode/settings.json b/.vscode/settings.json index b6676953..d805572a 100644 --- a/.vscode/settings.json +++ b/.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" } } \ No newline at end of file diff --git a/src/main/c/includes/fluid/env/environment.h b/src/main/c/includes/fluid/env/environment.h index d40048ed..67be4901 100644 --- a/src/main/c/includes/fluid/env/environment.h +++ b/src/main/c/includes/fluid/env/environment.h @@ -4,6 +4,7 @@ #include #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; /** 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 b5ce07da..c12c8887 100644 --- a/src/main/c/includes/fluid/sim/grid2/solver_consts.h +++ b/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 */ diff --git a/src/main/c/includes/math/ode/conjugate_gradient.h b/src/main/c/includes/math/ode/conjugate_gradient.h index e80f1acc..fbc579be 100644 --- a/src/main/c/includes/math/ode/conjugate_gradient.h +++ b/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 \ No newline at end of file diff --git a/src/main/c/includes/math/ode/diffusion_ode.h b/src/main/c/includes/math/ode/diffusion_ode.h new file mode 100644 index 00000000..b964adbe --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/main/c/includes/math/ode/multigrid_parallel.h b/src/main/c/includes/math/ode/multigrid_parallel.h index 9f037bea..436afd18 100644 --- a/src/main/c/includes/math/ode/multigrid_parallel.h +++ b/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 FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){ - // residual = solver_multigrid_parallel_iterate(x,x0,a,c); - // fluid_grid2_set_bounds(BOUND_SET_DENSITY_PHI,x); - // iterations++; - // } - - for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){ - //iterate - solver_gauss_seidel_iterate_parallel(x,x0,a,c,DIM); - - //set bounds + //about ~40% faster than gauss seidel + float residual = 1; + int iterations = 0; + while(iterations < FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){ + residual = solver_multigrid_parallel_iterate(x,x0,a,c); 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); + // } } /** diff --git a/src/main/c/src/fluid/sim/grid2/grid2.c b/src/main/c/src/fluid/sim/grid2/grid2.c index beccd4a3..4addfe8a 100644 --- a/src/main/c/src/fluid/sim/grid2/grid2.c +++ b/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; - // - //Velocity step - // + //maintenance 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 - // gettimeofday(&tv,NULL); - // end = 1000000.0 * tv.tv_sec + tv.tv_usec; - // perMilli = (end - start) / 1000.0f; - // environment->state.timeTracking.velocityDiffuse = 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.velocityDiffuse = perMilli; + start = end; - // for(int i = 0; i < numChunks; i++){ - // Chunk * currentChunk = chunks[i]; - // //update the bounds arrays - // fluid_grid2_rewrite_bounds(currentChunk); + + //project + for(int i = 0; i < numChunks; i++){ + 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 - // gettimeofday(&tv,NULL); - // end = 1000000.0 * tv.tv_sec + tv.tv_usec; - // perMilli = (end - start) / 1000.0f; - // environment->state.timeTracking.velocityProject = 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.velocityProject = perMilli; + start = end; - // for(int i = 0; i < numChunks; i++){ - // Chunk * currentChunk = chunks[i]; - // Chunk * currentChunk = chunks[i]; - // //update the bounds arrays - // fluid_grid2_rewrite_bounds(currentChunk); + //advect + for(int i = 0; i < numChunks; i++){ + Chunk * currentChunk = chunks[i]; + //update the bounds arrays + 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 - // gettimeofday(&tv,NULL); - // end = 1000000.0 * tv.tv_sec + tv.tv_usec; - // perMilli = (end - start) / 1000.0f; - // environment->state.timeTracking.velocityAdvect = 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.velocityAdvect = perMilli; + start = end; - // for(int i = 0; i < numChunks; i++){ - // Chunk * currentChunk = chunks[i]; - // //update the bounds arrays - // fluid_grid2_rewrite_bounds(currentChunk); + //project again + for(int i = 0; i < numChunks; i++){ + Chunk * currentChunk = chunks[i]; + //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++){ - // Chunk * currentChunk = chunks[i]; - // //update the bounds arrays - // fluid_grid2_rewrite_bounds(environment, currentChunk); + //advect + for(int i = 0; i < numChunks; i++){ + Chunk * currentChunk = chunks[i]; + //update the bounds arrays + fluid_grid2_rewrite_bounds(environment, currentChunk); //advect density diff --git a/src/main/c/src/fluid/sim/grid2/velocity.c b/src/main/c/src/fluid/sim/grid2/velocity.c index 9435b1ec..c02b4549 100644 --- a/src/main/c/src/fluid/sim/grid2/velocity.c +++ b/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++){ - //transform u direction - solver_gauss_seidel_iterate_parallel(u,u0,a,c,DIM); + // //about ~30% faster + // for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){ + // //transform u direction + // solver_gauss_seidel_iterate_parallel(u,u0,a,c,DIM); - //transform v direction - solver_gauss_seidel_iterate_parallel(v,v0,a,c,DIM); + // //transform v direction + // solver_gauss_seidel_iterate_parallel(v,v0,a,c,DIM); - //transform w direction - solver_gauss_seidel_iterate_parallel(w,w0,a,c,DIM); + // //transform w direction + // 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); - fluid_grid2_set_bounds(environment,BOUND_SET_VECTOR_DIFFUSE_PHI_W,w); + 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(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 diff --git a/src/main/c/src/math/ode/conjugate_gradient.c b/src/main/c/src/math/ode/conjugate_gradient.c index 020f2b70..39876a6c 100644 --- a/src/main/c/src/math/ode/conjugate_gradient.c +++ b/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 1000000){ + // printf("%f\n",laplacian); + // } + } + } + } + convergence = 0; + denominator = CONJUGATE_GRADIENT_EPSILON; + for(k=1; k 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; kdt*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)] + ) + ; +} + diff --git a/src/test/c/fluid/sim/grid2/speed_tests.c b/src/test/c/fluid/sim/grid2/speed_tests.c index 84850e9a..0b615457 100644 --- a/src/test/c/fluid/sim/grid2/speed_tests.c +++ b/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_test2(); -// rVal += fluid_sim_grid2_speed_test3(); + rVal += fluid_sim_grid2_speed_test1(); + rVal += fluid_sim_grid2_speed_test2(); + rVal += fluid_sim_grid2_speed_test3(); return rVal; } \ No newline at end of file