multigrid architecting
Some checks failed
studiorailgun/Renderer/pipeline/head There was a failure building this commit
Some checks failed
studiorailgun/Renderer/pipeline/head There was a failure building this commit
This commit is contained in:
austin
parent
96e03bea33
commit
81e083f6c2
@ -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-1; k++){
|
||||
for(j=1; j<GRIDDIM-1; j++){
|
||||
for(i=1; i<GRIDDIM-1; i++){
|
||||
@ -579,18 +570,7 @@ float solver_multigrid_store_residual_serial(float * phi, float * phi0, float *
|
||||
)
|
||||
);
|
||||
residualGrid[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)] = phi0[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)] - laplacian;
|
||||
// printf("%f - %f \n",phi0[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)],laplacian);
|
||||
residual_norm = residual_norm + (
|
||||
residualGrid[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)] * residualGrid[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)]
|
||||
);
|
||||
increment++;
|
||||
if(increment > 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;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user