From 1d75ad9af9d46eeb8c92f1cbb8569d153e63cff8 Mon Sep 17 00:00:00 2001 From: austin Date: Wed, 11 Dec 2024 18:16:41 -0500 Subject: [PATCH] fully recursive multigrid solver --- src/main/c/includes/math/ode/gauss_seidel.h | 147 ++++++---- src/main/c/includes/math/ode/multigrid.h | 12 +- src/main/c/src/fluid/sim/grid2/velocity.c | 2 +- src/main/c/src/math/ode/multigrid.c | 296 ++++++++++++++++++-- 4 files changed, 373 insertions(+), 84 deletions(-) diff --git a/src/main/c/includes/math/ode/gauss_seidel.h b/src/main/c/includes/math/ode/gauss_seidel.h index eaf6eea0..bd48b608 100644 --- a/src/main/c/includes/math/ode/gauss_seidel.h +++ b/src/main/c/includes/math/ode/gauss_seidel.h @@ -14,57 +14,6 @@ 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; kgridDim-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 @@ -77,9 +26,6 @@ static inline void solver_gauss_seidel_iterate_parallel(float * phi, float * phi 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= 10){ + __m256 aScalar = _mm256_set1_ps(a); + __m256 cScalar = _mm256_set1_ps(c); + //transform u direction + for(k=1; kgridDim-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; + } + } + } + } + } else if(gridDim >= 6){ + __m128 aScalar = _mm_set1_ps(a); + __m128 cScalar = _mm_set1_ps(c); + //transform u direction + for(k=1; kgridDim-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; + } + } + } + } + } else { + solver_gauss_seidel_iterate_serial(phi,phi0,a,c,gridDim); + } +} + #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 index a36ad121..9e808778 100644 --- a/src/main/c/includes/math/ode/multigrid.h +++ b/src/main/c/includes/math/ode/multigrid.h @@ -15,7 +15,17 @@ * @param c The c const * @return The residual */ -float solver_multigrid_iterate(float * phi, float * phi0, float a, float c); +float solver_multigrid_iterate_serial(float * phi, float * phi0, float a, float c); + +/** + * 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 + * @return The residual + */ +float solver_multigrid_iterate_parallel(float * phi, float * phi0, float a, float c); diff --git a/src/main/c/src/fluid/sim/grid2/velocity.c b/src/main/c/src/fluid/sim/grid2/velocity.c index ea2bf452..f6274540 100644 --- a/src/main/c/src/fluid/sim/grid2/velocity.c +++ b/src/main/c/src/fluid/sim/grid2/velocity.c @@ -181,7 +181,7 @@ LIBRARY_API void fluid_grid2_solveProjection( float residual = 1; int iteration = 0; while(iteration < FLUID_GRID2_LINEARSOLVERTIMES && (residual > FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE || residual < -FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE)){ - residual = solver_multigrid_iterate(p,div,a,c); + residual = solver_multigrid_iterate_serial(p,div,a,c); fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p); iteration++; } diff --git a/src/main/c/src/math/ode/multigrid.c b/src/main/c/src/math/ode/multigrid.c index 47377d40..7440a06d 100644 --- a/src/main/c/src/math/ode/multigrid.c +++ b/src/main/c/src/math/ode/multigrid.c @@ -18,6 +18,7 @@ static int halfDim = ((DIM - 2) / 2) + 2; * Dimension of the quarter resolution grid */ static int quarterDim = ((DIM - 2) / 4) + 2; +static int LOWEST_DIM = ((DIM - 2) / 4) + 2; /** * The full resolution grids @@ -38,10 +39,21 @@ static float * quarterGridPhi = NULL; static float * quarterGridPhi0 = NULL; static float * quarterGridResidual = NULL; +/** + * The eighth resolution grids + */ +static float * eighthGridPhi = NULL; +static float * eighthGridPhi0 = NULL; +static float * eighthGridResidual = NULL; + float solver_multigrid_calculate_residual_norm_serial(float * phi, float * phi0, float a, float c, int GRIDDIM); float solver_multigrid_calculate_residual_parallel(float * phi, float * phi0, float a, float c, int GRIDDIM); float solver_multigrid_store_residual_serial(float * phi, float * phi0, float * residualGrid, float a, float c, int GRIDDIM); +float * get_current_phi(int dim); +float * get_current_phi0(int dim); +float * get_current_residual(int dim); +void initialization_check(); /** * Relaxes an ODE matrix by 1 iteration of multigrid method @@ -51,40 +63,171 @@ float solver_multigrid_store_residual_serial(float * phi, float * phi0, float * * @param c The c const * @return The residual */ -float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){ +float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float a, float c, int GRIDDIM){ + initialization_check(); + float residual; + printf("Current dim: %d \n",GRIDDIM); + fflush(stdout); + int LOWERDIM = ((GRIDDIM - 2) / 2) + 2; + float * currResidual = get_current_residual(GRIDDIM); + float * lowerPhi = get_current_phi(LOWERDIM); + float * lowerPhi0 = get_current_phi0(LOWERDIM); + float * lowerResidual = get_current_residual(LOWERDIM); + + printf("Solved arrays: %p %p %p %p \n",currResidual,lowerPhi, lowerPhi0, lowerResidual); + fflush(stdout); + + printf("smooth\n"); + fflush(stdout); + // + //gauss-seidel at highest res + solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM); + if(GRIDDIM == DIM){ + fluid_grid2_set_bounds(0,phi); + } + //compute residuals + residual = solver_multigrid_store_residual_serial(phi,phi0,currResidual,a,c,GRIDDIM); + + printf("restrict\n"); + fflush(stdout); + //restrict + //(current operator is injection -- inject r^2 from this grid at phi0 of the smaller grid) + for(int x = 0; x < LOWERDIM; x++){ + for(int y = 0; y < LOWERDIM; y++){ + for(int z = 0; z < LOWERDIM; z++){ + lowerPhi[solver_gauss_seidel_get_index(x,y,z,LOWERDIM)] = 0; + lowerPhi0[solver_gauss_seidel_get_index(x,y,z,LOWERDIM)] = 0; + } + } + } + //populate grid + for(int x = 1; x < LOWERDIM - 1; x++){ + for(int y = 1; y < LOWERDIM - 1; y++){ + for(int z = 1; z < LOWERDIM - 1; z++){ + //direct transfer operator (faster, lower accuracy) + lowerPhi0[solver_gauss_seidel_get_index(x,y,z,LOWERDIM)] = currResidual[solver_gauss_seidel_get_index(x*2,y*2,z*2,GRIDDIM)]; + } + } + } + + printf("recurse\n"); + fflush(stdout); + //solve next-coarsest grid + if(GRIDDIM <= LOWEST_DIM){ + //smooth + solver_gauss_seidel_iterate_parallel(quarterGridPhi,quarterGridPhi0,a,c,quarterDim); + //compute residual + residual = solver_multigrid_store_residual_serial(quarterGridPhi,quarterGridPhi0,quarterGridResidual,a,c,quarterDim); + } else { + printf("solve\n"); + fflush(stdout); + float solution = + ( + phi0[solver_gauss_seidel_get_index(1,1,1,GRIDDIM)] + + phi0[solver_gauss_seidel_get_index(2,1,1,GRIDDIM)] + + phi0[solver_gauss_seidel_get_index(1,2,1,GRIDDIM)] + + phi0[solver_gauss_seidel_get_index(1,1,2,GRIDDIM)] + + phi0[solver_gauss_seidel_get_index(1,2,2,GRIDDIM)] + + phi0[solver_gauss_seidel_get_index(2,1,2,GRIDDIM)] + + phi0[solver_gauss_seidel_get_index(2,2,1,GRIDDIM)] + + phi0[solver_gauss_seidel_get_index(2,2,2,GRIDDIM)] + ) / 8.0f + ; + + printf("interpolate solution\n"); + fflush(stdout); + //interpolate from the lower grid + for(int x = 1; x < GRIDDIM - 1; x++){ + for(int y = 1; y < GRIDDIM - 1; y++){ + for(int z = 1; z < GRIDDIM - 1; z++){ + //direct transfer operator (faster, lower accuracy) + phi[solver_gauss_seidel_get_index(x,y,z,GRIDDIM)] = phi[solver_gauss_seidel_get_index(x,y,z,GRIDDIM)] + solution; + } + } + } + } + + printf("interpolate\n"); + fflush(stdout); + //interpolate from the lower grid + for(int x = 1; x < GRIDDIM - 1; x++){ + for(int y = 1; y < GRIDDIM - 1; y++){ + for(int z = 1; z < GRIDDIM - 1; z++){ + //direct transfer operator (faster, lower accuracy) + phi[solver_gauss_seidel_get_index(x,y,z,GRIDDIM)] = + phi[solver_gauss_seidel_get_index(x,y,z,GRIDDIM)] + + lowerPhi[solver_gauss_seidel_get_index( x/2+0, y/2+0, z/2+0 ,LOWERDIM)] + ; + + //interpolation operator (slower, better accuracy) + // phi[solver_gauss_seidel_get_index(x,y,z,DIM)] = + // phi[solver_gauss_seidel_get_index(x,y,z,DIM)] + + // ( + // halfGridPhi[solver_gauss_seidel_get_index( x/2+0, y/2+0, z/2+0 ,halfDim)] + + // halfGridPhi[solver_gauss_seidel_get_index( x/2+0, y/2+0, z/2+1 ,halfDim)] + + // halfGridPhi[solver_gauss_seidel_get_index( x/2+0, y/2+1, z/2+0 ,halfDim)] + + // halfGridPhi[solver_gauss_seidel_get_index( x/2+0, y/2+1, z/2+1 ,halfDim)] + + // halfGridPhi[solver_gauss_seidel_get_index( x/2+1, y/2+0, z/2+0 ,halfDim)] + + // halfGridPhi[solver_gauss_seidel_get_index( x/2+1, y/2+0, z/2+1 ,halfDim)] + + // halfGridPhi[solver_gauss_seidel_get_index( x/2+1, y/2+1, z/2+0 ,halfDim)] + + // halfGridPhi[solver_gauss_seidel_get_index( x/2+1, y/2+1, z/2+1 ,halfDim)] + // ) + // ; + } + } + } + + + + + + + printf("smooth\n"); + fflush(stdout); + // // full res - if(fullGridResidual == NULL){ - fullGridResidual = (float *)calloc(1,DIM * DIM * DIM * sizeof(float)); + // + // + //smooth + solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM); + if(GRIDDIM == DIM){ + fluid_grid2_set_bounds(0,phi); } - // half res - if(halfGridPhi == NULL){ - halfGridPhi = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float)); - } - if(halfGridPhi0 == NULL){ - halfGridPhi0 = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float)); - } - if(halfGridResidual == NULL){ - halfGridResidual = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float)); - } + + return solver_multigrid_calculate_residual_norm_serial(phi,phi0,a,c,GRIDDIM); +} - // quarter res - if(quarterGridPhi == NULL){ - quarterGridPhi = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float)); - } - if(quarterGridPhi0 == NULL){ - quarterGridPhi0 = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float)); - } - if(quarterGridResidual == NULL){ - quarterGridResidual = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float)); - } +/** + * 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 + * @return The residual + */ +float solver_multigrid_iterate_serial(float * phi, float * phi0, float a, float c){ + return solver_multigrid_iterate_serial_recursive(phi,phi0,a,c,DIM); +} + + +/** + * 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 + * @return The residual + */ +float solver_multigrid_iterate_parallel(float * phi, float * phi0, float a, float c){ + initialization_check(); float residual; // //gauss-seidel at highest res - solver_gauss_seidel_iterate_serial(phi,phi0,a,c,DIM); + solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,DIM); fluid_grid2_set_bounds(0,phi); //compute residuals residual = solver_multigrid_store_residual_serial(phi,phi0,fullGridResidual,a,c,DIM); @@ -117,7 +260,7 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){ //half res // //smooth - solver_gauss_seidel_iterate_serial(halfGridPhi,halfGridPhi0,a,c,halfDim); + solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim); //compute residual residual = solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim); @@ -151,7 +294,7 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){ } //smooth - solver_gauss_seidel_iterate_serial(quarterGridPhi,quarterGridPhi0,a,c,quarterDim); + solver_gauss_seidel_iterate_parallel(quarterGridPhi,quarterGridPhi0,a,c,quarterDim); //compute residual residual = solver_multigrid_store_residual_serial(quarterGridPhi,quarterGridPhi0,quarterGridResidual,a,c,quarterDim); @@ -193,7 +336,7 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){ // //smooth - solver_gauss_seidel_iterate_serial(halfGridPhi,halfGridPhi0,a,c,halfDim); + solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim); //compute residual residual = solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim); @@ -236,7 +379,7 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){ // // //smooth - solver_gauss_seidel_iterate_serial(phi,phi0,a,c,DIM); + solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,DIM); fluid_grid2_set_bounds(0,phi); @@ -245,6 +388,103 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){ +/** + * Gets the phi to use for the current level of the multigrid solver + */ +float * get_current_phi(int dim){ + if(dim == ((DIM-2)/2) + 2){ + return halfGridPhi; + } else if(dim == ((DIM-2)/4) + 2){ + return quarterGridPhi; + } else if(dim == ((DIM-2)/8) + 2){ + return eighthGridPhi; + } else { + printf("[get_current_phi] Invalid dim detected! %d\n",dim); + fflush(stdout); + return NULL; + } +} + +/** + * Gets the phi0 to use for the current level of the multigrid solver + */ +float * get_current_phi0(int dim){ + if(dim == ((DIM-2)/2) + 2){ + return halfGridPhi0; + } else if(dim == ((DIM-2)/4) + 2){ + return quarterGridPhi0; + } else if(dim == ((DIM-2)/8) + 2){ + return eighthGridPhi0; + } else { + printf("[get_current_phi0] Invalid dim detected! %d\n",dim); + fflush(stdout); + return NULL; + } +} + +/** + * Gets the residual to use for the current level of the multigrid solver + */ +float * get_current_residual(int dim){ + if(dim == 18){ + return fullGridResidual; + } else if(dim == ((DIM-2)/2) + 2){ + return halfGridResidual; + } else if(dim == ((DIM-2)/4) + 2){ + return quarterGridResidual; + } else if(dim == ((DIM-2)/8) + 2){ + return eighthGridResidual; + } else { + printf("[get_current_residual] Invalid dim detected! %d\n",dim); + fflush(stdout); + return NULL; + } +} + +/** + * Verifies that all grids are allocated + */ +void initialization_check(){ + // full res + if(fullGridResidual == NULL){ + fullGridResidual = (float *)calloc(1,DIM * DIM * DIM * sizeof(float)); + } + + // half res + if(halfGridPhi == NULL){ + halfGridPhi = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float)); + } + if(halfGridPhi0 == NULL){ + halfGridPhi0 = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float)); + } + if(halfGridResidual == NULL){ + halfGridResidual = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float)); + } + + // quarter res + if(quarterGridPhi == NULL){ + quarterGridPhi = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float)); + } + if(quarterGridPhi0 == NULL){ + quarterGridPhi0 = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float)); + } + if(quarterGridResidual == NULL){ + quarterGridResidual = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float)); + } + + // quarter res + int eighthResDim = ((DIM-2)/8) + 2; + if(eighthGridPhi == NULL){ + eighthGridPhi = (float *)calloc(1,eighthResDim * eighthResDim * eighthResDim * sizeof(float)); + } + if(eighthGridPhi0 == NULL){ + eighthGridPhi0 = (float *)calloc(1,eighthResDim * eighthResDim * eighthResDim * sizeof(float)); + } + if(eighthGridResidual == NULL){ + eighthGridResidual = (float *)calloc(1,eighthResDim * eighthResDim * eighthResDim * sizeof(float)); + } +} + /** * Calculates the residual of the grid