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

fix multigrid solver
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 17:31:52 -05:00
parent 02e7cd2374
commit e24332df83
5 changed files with 250 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
src/main/c/includes/fluid/sim/grid2/solver_consts.h
View File

@ -27,6 +27,11 @@
*/
#define FLUID_GRID2_REALLY_SMALL_VALUE 0.00001f
/**
* The tolerance for convergence of the projection operator
*/
#define FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE 0.01f
/**
* Diffusion constant
*/

1
src/main/c/src/fluid/sim/grid2/density.c
View File

@ -8,6 +8,7 @@
#include "fluid/queue/chunk.h"
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "math/ode/multigrid.h"
#include "math/ode/gauss_seidel.h"

9
src/main/c/src/fluid/sim/grid2/velocity.c
View File

@ -180,15 +180,20 @@ LIBRARY_API void fluid_grid2_solveProjection(
//perform iteration of v cycle multigrid method
float residual = 1;
int iteration = 0;
while(iteration < FLUID_GRID2_LINEARSOLVERTIMES && (residual > FLUID_GRID2_REALLY_SMALL_VALUE || residual < -FLUID_GRID2_REALLY_SMALL_VALUE)){
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);
fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
iteration++;
}
if(residual > FLUID_GRID2_REALLY_SMALL_VALUE || residual < -FLUID_GRID2_REALLY_SMALL_VALUE){
if(residual > FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE || residual < -FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE){
printf("Projection residual didn't converge! %f \n",residual);
}
// for(i = 0; i < 100; i++){
// solver_gauss_seidel_iterate_parallel(p,div,a,c,DIM);
// fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
// }
// solver_conjugate_gradient_solve_serial(p,div,a,c);
//finest grain iteration with conjugate gradient method

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

@ -1,8 +1,10 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "fluid/queue/chunk.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/solver_consts.h"
#include "math/ode/gauss_seidel.h"
#include "math/ode/multigrid.h"
#include "util/vector.h"
@ -17,20 +19,29 @@ static int halfDim = ((DIM - 2) / 2) + 2;
*/
static int quarterDim = ((DIM - 2) / 4) + 2;
/**
* The full resolution grids
*/
static float * fullGridResidual = NULL;
/**
* The half resolution grids
*/
static float * halfGridPhi = NULL;
static float * halfGridPhi0 = NULL;
static float * halfGridResidual = NULL;
/**
* The quarter resolution grids
*/
static float * quarterGridPhi = NULL;
static float * quarterGridPhi0 = NULL;
static float * quarterGridResidual = NULL;
float solver_multigrid_calculate_residual(float * phi, float * phi0, float a, float c);
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);
/**
* Relaxes an ODE matrix by 1 iteration of multigrid method
@ -41,27 +52,46 @@ float solver_multigrid_calculate_residual(float * phi, float * phi0, float a, fl
* @return The residual
*/
float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
// 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));
}
float residual;
//
//gauss-seidel at highest res
solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,DIM);
solver_gauss_seidel_iterate_serial(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);
//
//half res
//clear grid
//restrict
//(current operator is injection -- inject r^2 from this grid at phi0 of the smaller grid)
for(int x = 0; x < halfDim - 1; x++){
for(int y = 0; y < halfDim - 1; y++){
for(int z = 0; z < halfDim - 1; z++){
@ -74,19 +104,34 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
for(int x = 1; x < halfDim - 1; x++){
for(int y = 1; y < halfDim - 1; y++){
for(int z = 1; z < halfDim - 1; z++){
halfGridPhi[solver_gauss_seidel_get_index(x,y,z,halfDim)] = phi[solver_gauss_seidel_get_index(x*2,y*2,z*2,DIM)];
halfGridPhi0[solver_gauss_seidel_get_index(x,y,z,halfDim)] = phi0[solver_gauss_seidel_get_index(x*2,y*2,z*2,DIM)];
//direct transfer operator (faster, lower accuracy)
halfGridPhi0[solver_gauss_seidel_get_index(x,y,z,halfDim)] = fullGridResidual[solver_gauss_seidel_get_index(x*2,y*2,z*2,DIM)];
}
}
}
solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim);
//
//half res
//
//smooth
solver_gauss_seidel_iterate_serial(halfGridPhi,halfGridPhi0,a,c,halfDim);
//compute residual
residual = solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim);
//
//quarter res
//
//clear grid
//restrict
//(current operator is injection -- inject r^2 from this grid at phi0 of the smaller grid)
for(int x = 0; x < quarterDim - 1; x++){
for(int y = 0; y < quarterDim - 1; y++){
for(int z = 0; z < quarterDim - 1; z++){
@ -99,12 +144,46 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
for(int x = 1; x < quarterDim - 1; x++){
for(int y = 1; y < quarterDim - 1; y++){
for(int z = 1; z < quarterDim - 1; z++){
quarterGridPhi[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = halfGridPhi[solver_gauss_seidel_get_index(x*2,y*2,z*2,halfDim)];
quarterGridPhi0[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = halfGridPhi0[solver_gauss_seidel_get_index(x*2,y*2,z*2,halfDim)];
//direct transfer operator (faster, lower accuracy)
quarterGridPhi0[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = halfGridResidual[solver_gauss_seidel_get_index(x*2,y*2,z*2,halfDim)];
}
}
}
solver_gauss_seidel_iterate_parallel(quarterGridPhi,quarterGridPhi0,a,c,quarterDim);
//smooth
solver_gauss_seidel_iterate_serial(quarterGridPhi,quarterGridPhi0,a,c,quarterDim);
//compute residual
residual = 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++){
for(int y = 1; y < halfDim - 1; y++){
for(int z = 1; z < halfDim - 1; z++){
//direct transfer operator (faster, lower accuracy)
halfGridPhi[solver_gauss_seidel_get_index(x,y,z,halfDim)] =
halfGridPhi[solver_gauss_seidel_get_index(x,y,z,halfDim)] +
quarterGridPhi[solver_gauss_seidel_get_index( x/2+0, y/2+0, z/2+0 ,quarterDim)]
;
//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)]
// )
// ;
}
}
}
@ -112,92 +191,136 @@ float solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
//
//half res
//
//clear grid
for(int x = 0; x < halfDim - 1; x++){
for(int y = 0; y < halfDim - 1; y++){
for(int z = 0; z < halfDim - 1; z++){
halfGridPhi[solver_gauss_seidel_get_index(x,y,z,halfDim)] = 0;
halfGridPhi0[solver_gauss_seidel_get_index(x,y,z,halfDim)] = 0;
}
}
}
//populate grid
for(int x = 1; x < halfDim - 1; x++){
for(int y = 1; y < halfDim - 1; y++){
for(int z = 1; z < halfDim - 1; z++){
halfGridPhi[solver_gauss_seidel_get_index(x,y,z,halfDim)] = quarterGridPhi[solver_gauss_seidel_get_index(x/2,y/2,z/2,quarterDim)];
halfGridPhi0[solver_gauss_seidel_get_index(x,y,z,halfDim)] = quarterGridPhi0[solver_gauss_seidel_get_index(x/2,y/2,z/2,quarterDim)];
}
}
}
solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim);
//smooth
solver_gauss_seidel_iterate_serial(halfGridPhi,halfGridPhi0,a,c,halfDim);
//compute residual
residual = solver_multigrid_store_residual_serial(halfGridPhi,halfGridPhi0,halfGridResidual,a,c,halfDim);
//
//full res
//interpolate phi of the higher grid
for(int x = 1; x < DIM - 1; x++){
for(int y = 1; y < DIM - 1; y++){
for(int z = 1; z < DIM - 1; z++){
phi[solver_gauss_seidel_get_index(x,y,z,DIM)] = halfGridPhi[solver_gauss_seidel_get_index(x/2,y/2,z/2,halfDim)];
phi0[solver_gauss_seidel_get_index(x,y,z,DIM)] = halfGridPhi0[solver_gauss_seidel_get_index(x/2,y/2,z/2,halfDim)];
//direct transfer operator (faster, lower 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)]
;
//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)]
// )
// ;
}
}
}
//
//gauss-seidel at highest res
solver_gauss_seidel_iterate_parallel(phi,phi0,a,c,DIM);
// full res
//
//
//smooth
solver_gauss_seidel_iterate_serial(phi,phi0,a,c,DIM);
fluid_grid2_set_bounds(0,phi);
return solver_multigrid_calculate_residual(phi,phi0,a,c);
return solver_multigrid_calculate_residual_norm_serial(phi,phi0,a,c,DIM);
}
/**
* Calculates the residual of the grid
* @return Returns the residual norm of the grid
*/
float solver_multigrid_calculate_residual_norm_serial(float * phi, float * phi0, float a, float c, int GRIDDIM){
//calculate residual
float residual_norm = 0;
int i, j, k;
int increment = 0;
float h = FLUID_GRID2_H;
float residual;
for(k=1; k<GRIDDIM-1; k++){
for(j=1; j<GRIDDIM-1; j++){
for(i=1; i<GRIDDIM-1; i++){
float laplacian =
(
6 * phi[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)] +
- (
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)]
)
);
residual = phi0[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)] - laplacian;
residual_norm = residual_norm + (residual * residual);
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);
printf("residual_norm: %lf\n",residual_norm);
return residual_norm;
}
/**
* Calculates the residual of the grid
*/
float solver_multigrid_calculate_residual(float * phi, float * phi0, float a, float c){
float solver_multigrid_calculate_residual_parallel(float * phi, float * phi0, float a, float c, int GRIDDIM){
//calculate residual
__m256 aScalar = _mm256_set1_ps(a);
__m256 cScalar = _mm256_set1_ps(c);
__m256 collector = _mm256_setzero_ps();
__m256 vector;
float vec_sum_storage[8];
float residual = 1;
float residual = 0;
//transform u direction
int i, j, k;
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
//lower
i = 1;
vector = _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,DIM)]);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,DIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,DIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,DIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,DIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,DIM)]));
vector = _mm256_mul_ps(vector,aScalar);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,DIM)]));
vector = _mm256_div_ps(vector,cScalar);
collector = _mm256_add_ps(collector,vector);
//upper
i = 9;
vector = _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,DIM)]);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,DIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,DIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,DIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,DIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,DIM)]));
vector = _mm256_mul_ps(vector,aScalar);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,DIM)]));
vector = _mm256_div_ps(vector,cScalar);
collector = _mm256_add_ps(collector,vector);
int increment = 0;
for(k=1; k<GRIDDIM-1; k++){
for(j=1; j<GRIDDIM-1; j++){
for(i=1; i<GRIDDIM-1; i=i+8){
vector = _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,GRIDDIM)]);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,GRIDDIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,GRIDDIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,GRIDDIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,GRIDDIM)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,GRIDDIM)]));
vector = _mm256_mul_ps(vector,aScalar);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)]));
vector = _mm256_div_ps(vector,cScalar);
collector = _mm256_add_ps(collector,vector);
increment++;
if(increment > GRIDDIM*GRIDDIM*GRIDDIM){
printf("INCREMENT FAILURE %d %d %d %d \n",i,j,k,increment);
return -1;
}
}
}
}
_mm256_storeu_ps(vec_sum_storage,collector);
@ -208,4 +331,45 @@ float solver_multigrid_calculate_residual(float * phi, float * phi0, float a, fl
vec_sum_storage[6] + vec_sum_storage[7]
;
return residual;
}
/**
* Calculates the residual of the grid
*/
float 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++){
float laplacian =
(
6 * phi[solver_gauss_seidel_get_index(i,j,k,GRIDDIM)] +
- (
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)]
)
);
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;
}

2
src/test/c/fluid/sim/grid2/solve_projection_tests.c
View File

@ -18,7 +18,7 @@
/**
* Error margin for tests
*/
#define FLUID_GRId2_PROJECTION_ERROR_MARGIN 0.00001f
#define FLUID_GRId2_PROJECTION_ERROR_MARGIN 0.001f
/**
* Testing gradient approximation