From 81e083f6c2231e1e9cf67661f1824574ea0115f5 Mon Sep 17 00:00:00 2001 From: austin Date: Wed, 11 Dec 2024 18:29:58 -0500 Subject: [PATCH] multigrid architecting --- src/main/c/src/math/ode/multigrid.c | 134 ++++++++++++---------------- 1 file changed, 57 insertions(+), 77 deletions(-) diff --git a/src/main/c/src/math/ode/multigrid.c b/src/main/c/src/math/ode/multigrid.c index 6212eab1..23325276 100644 --- a/src/main/c/src/math/ode/multigrid.c +++ b/src/main/c/src/math/ode/multigrid.c @@ -49,11 +49,13 @@ 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); +void 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(); +void restrict_serial(float * currResidual, int GRIDDIM, float * lowerPhi, float * lowerPhi0, int LOWERDIM); +void prolongate_serial(float * phi, int GRIDDIM, float * lowerPhi, int LOWERDIM); /** * Relaxes an ODE matrix by 1 iteration of multigrid method @@ -65,7 +67,6 @@ void initialization_check(); */ float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float a, float c, int GRIDDIM){ initialization_check(); - float residual; int LOWERDIM = ((GRIDDIM - 2) / 2) + 2; float * currResidual = get_current_residual(GRIDDIM); @@ -80,34 +81,17 @@ float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float fluid_grid2_set_bounds(0,phi); } //compute residuals - residual = solver_multigrid_store_residual_serial(phi,phi0,currResidual,a,c,GRIDDIM); + solver_multigrid_store_residual_serial(phi,phi0,currResidual,a,c,GRIDDIM); //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)]; - } - } - } + restrict_serial(currResidual,GRIDDIM,lowerPhi,lowerPhi0,LOWERDIM); //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); + solver_multigrid_store_residual_serial(quarterGridPhi,quarterGridPhi0,quarterGridResidual,a,c,quarterDim); } else { float solution = ( @@ -122,7 +106,7 @@ float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float ) / 8.0f ; - //interpolate from the lower grid + //interpolate solution to this 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++){ @@ -134,42 +118,8 @@ float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float } //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)] - ; + prolongate_serial(phi,GRIDDIM,lowerPhi,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)] - // ) - // ; - } - } - } - - - - - - - // - // full res - // - // //smooth solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM); if(GRIDDIM == DIM){ @@ -211,7 +161,7 @@ float solver_multigrid_iterate_parallel(float * phi, float * phi0, float a, floa 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); + solver_multigrid_store_residual_serial(phi,phi0,fullGridResidual,a,c,DIM); //restrict @@ -243,7 +193,7 @@ float solver_multigrid_iterate_parallel(float * phi, float * phi0, float a, floa //smooth solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim); //compute residual - residual = solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim); + solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim); @@ -277,7 +227,7 @@ float solver_multigrid_iterate_parallel(float * phi, float * phi0, float a, floa //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); + solver_multigrid_store_residual_serial(quarterGridPhi,quarterGridPhi0,quarterGridResidual,a,c,quarterDim); //interpolate into phi of the higher grid for(int x = 1; x < halfDim - 1; x++){ @@ -319,7 +269,7 @@ float solver_multigrid_iterate_parallel(float * phi, float * phi0, float a, floa //smooth solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim); //compute residual - residual = solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim); + solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim); //interpolate phi of the higher grid for(int x = 1; x < DIM - 1; x++){ @@ -422,6 +372,50 @@ float * get_current_residual(int dim){ } } +/** + * Serially restricts the current residual into the lower phi grid + */ +void restrict_serial(float * currResidual, int GRIDDIM, float * lowerPhi, float * lowerPhi0, int LOWERDIM){ + int x, y, z; + //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)]; + } + } + } +} + +/** + * Prolongates a lower grid into a higher grid + */ +void prolongate_serial(float * phi, int GRIDDIM, float * lowerPhi, int LOWERDIM){ + int x, y, z; + 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)] + ; + } + } + } +} + /** * Verifies that all grids are allocated */ @@ -557,12 +551,9 @@ float solver_multigrid_calculate_residual_parallel(float * phi, float * phi0, fl /** * Calculates the residual of the grid */ -float solver_multigrid_store_residual_serial(float * phi, float * phi0, float * residualGrid, float a, float c, int GRIDDIM){ +void solver_multigrid_store_residual_serial(float * phi, float * phi0, float * residualGrid, float a, float c, int GRIDDIM){ //calculate residual int i, j, k; - int increment = 0; - float h = FLUID_GRID2_H; - float residual_norm = 0; for(k=1; k GRIDDIM*GRIDDIM*GRIDDIM){ - printf("INCREMENT FAILURE %d %d %d %d \n",i,j,k,increment); - return -1; - } } } } - residual_norm = (float)sqrt(residual_norm); - return residual_norm; } \ No newline at end of file