studiorailgun/Renderer
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

parallel residual storage
Some checks failed
studiorailgun/Renderer/pipeline/head There was a failure building this commit
Details

Browse Source
This commit is contained in:
austin 2024-12-11 19:16:18 -05:00
parent 37a6180be8
commit 99e89c9d88
1 changed files with 61 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
src/main/c/src/math/ode/multigrid.c
View File

@ -61,6 +61,7 @@ void prolongate_serial(float * phi, int GRIDDIM, float * lowerPhi, int LOWERDIM)
//parallelized operations
void restrict_parallel(float * currResidual, int GRIDDIM, float * lowerPhi, float * lowerPhi0, int LOWERDIM);
void prolongate_parallel(float * phi, int GRIDDIM, float * lowerPhi, int LOWERDIM);
void solver_multigrid_store_residual_parallel(float * phi, float * phi0, float * residualGrid, float a, float c, int GRIDDIM);
/**
* Relaxes an ODE matrix by 1 iteration of multigrid method
@ -71,7 +72,7 @@ void prolongate_parallel(float * phi, int GRIDDIM, float * lowerPhi, int LOWERDI
* @param GRIDDIM The dimension of the phi grid
* @return The residual
*/
float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float a, float c, int GRIDDIM){
void solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float a, float c, int GRIDDIM){
int LOWERDIM = ((GRIDDIM - 2) / 2) + 2;
float * currResidual = get_current_residual(GRIDDIM);
float * lowerPhi = get_current_phi(LOWERDIM);
@ -82,10 +83,10 @@ float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float
solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM);
//compute residuals
solver_multigrid_store_residual_serial(phi,phi0,currResidual,a,c,GRIDDIM);
solver_multigrid_store_residual_parallel(phi,phi0,currResidual,a,c,GRIDDIM);
//restrict
restrict_serial(currResidual,GRIDDIM,lowerPhi,lowerPhi0,LOWERDIM);
restrict_parallel(currResidual,GRIDDIM,lowerPhi,lowerPhi0,LOWERDIM);
//solve next-coarsest grid
if(GRIDDIM <= LOWEST_DIM){
@ -116,12 +117,10 @@ float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float
}
//interpolate from the lower grid
prolongate_serial(phi,GRIDDIM,lowerPhi,LOWERDIM);
prolongate_parallel(phi,GRIDDIM,lowerPhi,LOWERDIM);
//smooth
solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM);
return solver_multigrid_calculate_residual_norm_serial(phi,phi0,a,c,GRIDDIM);
}
@ -135,7 +134,8 @@ float solver_multigrid_iterate_serial_recursive(float * phi, float * phi0, float
*/
float solver_multigrid_iterate_serial(float * phi, float * phi0, float a, float c){
initialization_check();
return solver_multigrid_iterate_serial_recursive(phi,phi0,a,c,DIM);
solver_multigrid_iterate_serial_recursive(phi,phi0,a,c,DIM);
return solver_multigrid_calculate_residual_norm_serial(phi,phi0,a,c,DIM);
}
@ -148,7 +148,7 @@ float solver_multigrid_iterate_serial(float * phi, float * phi0, float a, float
* @param GRIDDIM The dimension of the phi grid
* @return The residual
*/
float solver_multigrid_iterate_parallel_recursive(float * phi, float * phi0, float a, float c, int GRIDDIM){
void solver_multigrid_iterate_parallel_recursive(float * phi, float * phi0, float a, float c, int GRIDDIM){
int LOWERDIM = ((GRIDDIM - 2) / 2) + 2;
float * currResidual = get_current_residual(GRIDDIM);
float * lowerPhi = get_current_phi(LOWERDIM);
@ -159,7 +159,7 @@ float solver_multigrid_iterate_parallel_recursive(float * phi, float * phi0, flo
solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM);
//compute residuals
solver_multigrid_store_residual_serial(phi,phi0,currResidual,a,c,GRIDDIM);
solver_multigrid_store_residual_parallel(phi,phi0,currResidual,a,c,GRIDDIM);
//restrict
restrict_parallel(currResidual,GRIDDIM,lowerPhi,lowerPhi0,LOWERDIM);
@ -176,8 +176,6 @@ float solver_multigrid_iterate_parallel_recursive(float * phi, float * phi0, flo
//smooth
solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,GRIDDIM);
return solver_multigrid_calculate_residual_norm_serial(phi,phi0,a,c,GRIDDIM);
}
/**
@ -190,7 +188,8 @@ float solver_multigrid_iterate_parallel_recursive(float * phi, float * phi0, flo
*/
float solver_multigrid_iterate_parallel(float * phi, float * phi0, float a, float c){
initialization_check();
return solver_multigrid_iterate_parallel_recursive(phi,phi0,a,c,DIM);
solver_multigrid_iterate_parallel_recursive(phi,phi0,a,c,DIM);
return solver_multigrid_calculate_residual_norm_serial(phi,phi0,a,c,DIM);
}
@ -566,4 +565,54 @@ void solver_multigrid_store_residual_serial(float * phi, float * phi0, float * r
}
}
}
}
/**
* Calculates the residual of the grid
*/
void solver_multigrid_store_residual_parallel(float * phi, float * phi0, float * residualGrid, float a, float c, int GRIDDIM){
if(GRIDDIM < 10){
solver_multigrid_store_residual_serial(phi,phi0,residualGrid,a,c,GRIDDIM);
return;
}
__m256 laplacian;
__m256 constVec = _mm256_set1_ps(6);
//calculate residual
int i, j, k;
for(k=1; k<GRIDDIM-1; k++){
for(j=1; j<GRIDDIM-1; j++){
for(i=1; i<GRIDDIM-1; i=i+8){
laplacian =
_mm256_sub_ps(
_mm256_mul_ps(
_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)]),
constVec
),
_mm256_add_ps(
_mm256_add_ps(
_mm256_add_ps(
_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,GRIDDIM)]),
_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,GRIDDIM)])
),
_mm256_add_ps(
_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,GRIDDIM)]),
_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,GRIDDIM)])
)
),
_mm256_add_ps(
_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,GRIDDIM)]),
_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,GRIDDIM)])
)
)
);
_mm256_storeu_ps(
&residualGrid[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)],
_mm256_sub_ps(
_mm256_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)]),
laplacian
)
);
}
}
}
}