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 5950477e..013681b3 100644 --- a/src/main/c/includes/fluid/sim/grid2/solver_consts.h +++ b/src/main/c/includes/fluid/sim/grid2/solver_consts.h @@ -5,7 +5,7 @@ /** * The number of times to relax most solvers */ -#define FLUID_GRID2_LINEARSOLVERTIMES 10 +#define FLUID_GRID2_LINEARSOLVERTIMES 20 /** * Width of a single grid cell diff --git a/src/main/c/includes/math/ode/multigrid_parallel.h b/src/main/c/includes/math/ode/multigrid_parallel.h index 8607ffaf..9f037bea 100644 --- a/src/main/c/includes/math/ode/multigrid_parallel.h +++ b/src/main/c/includes/math/ode/multigrid_parallel.h @@ -4,7 +4,7 @@ #include "math/ode/gauss_seidel.h" #include "math/ode/multigrid.h" -#define MULTIGRID_PARALLEL_LOWEST_PARALLEL_DIM ((DIM - 2) / 1) + 2 +#define MULTIGRID_PARALLEL_LOWEST_PARALLEL_DIM ((DIM - 2) / 2) + 2 @@ -246,7 +246,7 @@ static inline void solver_multigrid_parallel_store_residual(float * phi, float * * @param GRIDDIM The dimension of the phi grid * @return The residual */ -static inline void solver_multigrid_parallel_iterate_recurse(float * phi, float * phi0, float a, float c, int GRIDDIM){ +static inline void solver_multigrid_parallel_iterate_recurse_v_cycle(float * phi, float * phi0, float a, float c, int GRIDDIM){ int LOWERDIM = ((GRIDDIM - 2) / 2) + 2; float * currResidual = solver_multigrid_get_current_residual(GRIDDIM); float * lowerPhi = solver_multigrid_get_current_phi(LOWERDIM); @@ -266,7 +266,47 @@ static inline void solver_multigrid_parallel_iterate_recurse(float * phi, float if(GRIDDIM <= MULTIGRID_PARALLEL_LOWEST_PARALLEL_DIM){ solver_multigrid_iterate_serial_recursive(lowerPhi,lowerPhi0,a,c,LOWERDIM); } else { - solver_multigrid_parallel_iterate_recurse(lowerPhi,lowerPhi0,a,c,LOWERDIM); + solver_multigrid_parallel_iterate_recurse_v_cycle(lowerPhi,lowerPhi0,a,c,LOWERDIM); + } + + //interpolate from the lower grid + prolongate_parallel(phi,GRIDDIM,lowerPhi,LOWERDIM); + + //smooth + solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM); +} + + +/** + * 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 + * @param GRIDDIM The dimension of the phi grid + * @return The residual + */ +static inline void solver_multigrid_parallel_iterate_recurse_f_cycle(float * phi, float * phi0, float a, float c, int GRIDDIM){ + int LOWERDIM = ((GRIDDIM - 2) / 2) + 2; + float * currResidual = solver_multigrid_get_current_residual(GRIDDIM); + float * lowerPhi = solver_multigrid_get_current_phi(LOWERDIM); + float * lowerPhi0 = solver_multigrid_get_current_phi0(LOWERDIM); + float * lowerResidual = solver_multigrid_get_current_residual(LOWERDIM); + + //smooth + solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM); + + //compute residuals + solver_multigrid_parallel_store_residual(phi,phi0,currResidual,a,c,GRIDDIM); + + //restrict + restrict_parallel(currResidual,GRIDDIM,lowerPhi,lowerPhi0,LOWERDIM); + + //solve next-coarsest grid + if(GRIDDIM <= MULTIGRID_PARALLEL_LOWEST_PARALLEL_DIM){ + solver_multigrid_iterate_serial_recursive(lowerPhi,lowerPhi0,a,c,LOWERDIM); + } else { + solver_multigrid_parallel_iterate_recurse_v_cycle(lowerPhi,lowerPhi0,a,c,LOWERDIM); } //interpolate from the lower grid @@ -286,7 +326,7 @@ static inline void solver_multigrid_parallel_iterate_recurse(float * phi, float */ static inline float solver_multigrid_parallel_iterate(float * phi, float * phi0, float a, float c){ solver_multigrid_initialization_check(); - solver_multigrid_parallel_iterate_recurse(phi,phi0,a,c,DIM); + solver_multigrid_parallel_iterate_recurse_v_cycle(phi,phi0,a,c,DIM); return solver_multigrid_calculate_residual_norm_serial(phi,phi0,a,c,DIM); } diff --git a/src/main/c/src/fluid/sim/grid2/velocity.c b/src/main/c/src/fluid/sim/grid2/velocity.c index a4cf0a1f..86663480 100644 --- a/src/main/c/src/fluid/sim/grid2/velocity.c +++ b/src/main/c/src/fluid/sim/grid2/velocity.c @@ -186,7 +186,7 @@ LIBRARY_API void fluid_grid2_solveProjection( iteration++; } if(residual > FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE || residual < -FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE){ - printf("Projection residual didn't converge! %f \n",residual); + // printf("Projection residual didn't converge! %f \n",residual); } // for(i = 0; i < 100; i++){ diff --git a/src/main/c/src/math/ode/multigrid.c b/src/main/c/src/math/ode/multigrid.c index e851f87e..fb6af8df 100644 --- a/src/main/c/src/math/ode/multigrid.c +++ b/src/main/c/src/math/ode/multigrid.c @@ -484,7 +484,7 @@ float solver_multigrid_calculate_residual_norm_serial(float * phi, float * phi0, } } residual_norm = (float)sqrt(residual_norm); - printf("residual_norm: %lf\n",residual_norm); + // printf("residual_norm: %lf\n",residual_norm); return residual_norm; }