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

rearching CG solver

Browse Source
This commit is contained in:
austin 2024-12-11 21:48:39 -05:00
parent 7dd7dcfc4d
commit 2247a5d9dd
5 changed files with 128 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.vscode/settings.json vendored
View File

@ -46,6 +46,7 @@
"multigrid.h": "c", "multigrid.h": "c",
"gauss_seidel.h": "c", "gauss_seidel.h": "c",
"ode_utils.h": "c", "ode_utils.h": "c",
"util.h": "c" "util.h": "c",
"conjugate_gradient.h": "c"
} }
} }

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

@ -5,7 +5,27 @@
/** /**
* The number of times to relax most solvers * The number of times to relax most solvers
*/ */
#define FLUID_GRID2_LINEARSOLVERTIMES 20 #define FLUID_GRID2_LINEARSOLVERTIMES 10
/**
* The number of times to relax most solvers
*/
#define FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS 20
/**
* Tolerance to target for multigrid precomputing in projection
*/
#define FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE 0.005f
/**
* The number of times to relax most solvers
*/
#define FLUID_GRID2_SOLVER_CG_MAX_ITERATIONS 10
/**
* Tolerance to target for conjugate gradient precomputing in projection
*/
#define FLUID_GRID2_SOLVER_CG_TOLERANCE 0.001f
/** /**
* Width of a single grid cell * Width of a single grid cell

11
src/main/c/includes/math/ode/conjugate_gradient.h
View File

@ -21,6 +21,15 @@ int solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c);
*/ */
void solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float a, float c); void solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float a, float c);
/**
* Initializes the conjugate gradient solver
* @param phi The phi array
* @param phi0 The phi array from the last frame
* @param a The a const
* @param c The c const
*/
void solver_conjugate_gradient_init_serial(float * phi, float * phi0, float a, float c);
/** /**
* Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method serially * Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method serially
* @param phi The phi array * @param phi The phi array
@ -28,6 +37,6 @@ void solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float
* @param a The a const * @param a The a const
* @param c The c const * @param c The c const
*/ */
void solver_conjugate_gradient_solve_serial(float * phi, float * phi0, float a, float c); float solver_conjugate_gradient_iterate_serial(float * phi, float * phi0, float a, float c);
#endif #endif

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

@ -181,21 +181,31 @@ LIBRARY_API void fluid_grid2_solveProjection(
//perform iteration of v cycle multigrid method //perform iteration of v cycle multigrid method
chunk->projectionResidual = 1; chunk->projectionResidual = 1;
chunk->projectionIterations = 0; chunk->projectionIterations = 0;
while(chunk->projectionIterations < FLUID_GRID2_LINEARSOLVERTIMES && (chunk->projectionResidual > FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE)){ while(chunk->projectionIterations < FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS && (chunk->projectionResidual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){
chunk->projectionResidual = solver_multigrid_parallel_iterate(p,div,a,c); chunk->projectionResidual = solver_multigrid_parallel_iterate(p,div,a,c);
fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p); fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
chunk->projectionIterations++; chunk->projectionIterations++;
} }
if(chunk->projectionResidual > FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_PROJECTION_CONVERGENCE_TOLERANCE){ // if(chunk->projectionResidual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE){
// printf("Projection residual didn't converge! %f \n",residual); // // printf("Projection residual didn't converge! %f \n",residual);
} // }
// for(i = 0; i < 100; i++){ // for(i = 0; i < 100; i++){
// solver_gauss_seidel_iterate_parallel(p,div,a,c,DIM); // solver_gauss_seidel_iterate_parallel(p,div,a,c,DIM);
// fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p); // fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
// } // }
// solver_conjugate_gradient_solve_serial(p,div,a,c); //init CG solver
// solver_conjugate_gradient_init_serial(p,div,a,c);
// //solve with CG
// while(chunk->projectionIterations < FLUID_GRID2_SOLVER_CG_MAX_ITERATIONS && (chunk->projectionResidual > FLUID_GRID2_SOLVER_CG_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_SOLVER_CG_TOLERANCE)){
// chunk->projectionResidual = solver_conjugate_gradient_iterate_serial(p,div,a,c);;
// fluid_grid2_set_bounds(FLUID_GRID2_BOUND_NO_DIR,p);
// chunk->projectionIterations++;
// }
// if(chunk->projectionResidual > FLUID_GRID2_SOLVER_CG_TOLERANCE || chunk->projectionResidual < -FLUID_GRID2_SOLVER_CG_TOLERANCE){
// // printf("Projection residual didn't converge! %f \n",residual);
// }
//finest grain iteration with conjugate gradient method //finest grain iteration with conjugate gradient method
// int shouldSolve = solver_conjugate_gradient_init(p,div,a,c); // int shouldSolve = solver_conjugate_gradient_init(p,div,a,c);

154
src/main/c/src/math/ode/conjugate_gradient.c
View File

@ -190,8 +190,87 @@ void solver_conjugate_gradient_iterate_parallel(float * phi, float * phi0, float
* @param phi0 The phi array from the last frame * @param phi0 The phi array from the last frame
* @param a The a const * @param a The a const
* @param c The c const * @param c The c const
* @return The residual
*/ */
void solver_conjugate_gradient_solve_serial(float * phi, float * phi0, float a, float c){ float solver_conjugate_gradient_iterate_serial(float * phi, float * phi0, float a, float c){
int i, j, k;
float convergence, denominator;
float laplacian, alpha, r_new_dot, beta;
//solve Ap
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
laplacian =
(
6 * p[ode_index (i , j , k ,DIM)] +
- (
p[ode_index (i+1 , j , k ,DIM)] +
p[ode_index (i-1 , j , k ,DIM)] +
p[ode_index (i , j+1 , k ,DIM)] +
p[ode_index (i , j-1 , k ,DIM)] +
p[ode_index (i , j , k+1 ,DIM)] +
p[ode_index (i , j , k-1 ,DIM)]
)
);
A[ode_index(i,j,k,DIM)] = laplacian;
}
}
}
convergence = 0;
denominator = CONJUGATE_GRADIENT_EPSILON;
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
convergence = convergence + r[ode_index(i,j,k,DIM)] * r[ode_index(i,j,k,DIM)];
denominator = denominator + p[ode_index(i,j,k,DIM)] * A[ode_index(i,j,k,DIM)];
}
}
}
if(denominator < CONJUGATE_GRADIENT_EPSILON && denominator > -CONJUGATE_GRADIENT_EPSILON){
printf("Divide by 0! %f \n", denominator);
fflush(stdout);
}
alpha = convergence / denominator;
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
phi[ode_index(i,j,k,DIM)] = phi[ode_index(i,j,k,DIM)] + alpha * p[ode_index(i,j,k,DIM)];
r[ode_index(i,j,k,DIM)] = r[ode_index(i,j,k,DIM)] - alpha * A[ode_index(i,j,k,DIM)];
}
}
}
r_new_dot = 0;
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
r_new_dot = r_new_dot + r[ode_index(i,j,k,DIM)] * r[ode_index(i,j,k,DIM)];
}
}
}
beta = r_new_dot / convergence;
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
p[ode_index(i,j,k,DIM)] = r[ode_index(i,j,k,DIM)] + beta * p[ode_index(i,j,k,DIM)];
}
}
}
return (float)sqrt(convergence);
}
/**
* Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method serially
* @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
*/
void solver_conjugate_gradient_init_serial(float * phi, float * phi0, float a, float c){
int i, j, k; int i, j, k;
if(p == NULL){ if(p == NULL){
p = (float *)calloc(1,DIM*DIM*DIM*sizeof(float)); p = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
@ -212,78 +291,5 @@ void solver_conjugate_gradient_solve_serial(float * phi, float * phi0, float a,
} }
} }
} }
float convergence = 1;
int framecount = 0;
while(convergence > CONJUGATE_GRADIENT_CONVERGENCE_THRESHOLD && framecount < max_frames){
//solve Ap
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
A[ode_index(i,j,k,DIM)] =
(
p[ode_index (i , j , k ,DIM)] +
a * (
p[ode_index (i+1 , j , k ,DIM)] +
p[ode_index (i-1 , j , k ,DIM)] +
p[ode_index (i , j+1 , k ,DIM)] +
p[ode_index (i , j-1 , k ,DIM)] +
p[ode_index (i , j , k+1 ,DIM)] +
p[ode_index (i , j , k-1 ,DIM)]
)
) / c;
}
}
}
convergence = 0;
float denominator = CONJUGATE_GRADIENT_EPSILON;
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
convergence = convergence + r[ode_index(i,j,k,DIM)] * r[ode_index(i,j,k,DIM)];
denominator = denominator + p[ode_index(i,j,k,DIM)] * A[ode_index(i,j,k,DIM)];
}
}
}
if(denominator < CONJUGATE_GRADIENT_EPSILON && denominator > -CONJUGATE_GRADIENT_EPSILON){
printf("Divide by 0! %f \n", denominator);
fflush(stdout);
break;
}
float alpha = convergence / denominator;
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
phi[ode_index(i,j,k,DIM)] = phi[ode_index(i,j,k,DIM)] + alpha * p[ode_index(i,j,k,DIM)];
r[ode_index(i,j,k,DIM)] = r[ode_index(i,j,k,DIM)] - alpha * A[ode_index(i,j,k,DIM)];
}
}
}
float r_new_dot = 0;
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
r_new_dot = r_new_dot + r[ode_index(i,j,k,DIM)] * r[ode_index(i,j,k,DIM)];
}
}
}
float beta = r_new_dot / convergence;
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
for(i = 1; i < DIM-1; i++){
p[ode_index(i,j,k,DIM)] = r[ode_index(i,j,k,DIM)] + beta * p[ode_index(i,j,k,DIM)];
}
}
}
framecount++;
}
if(framecount >= 100 || convergence > CONJUGATE_GRADIENT_CONVERGENCE_THRESHOLD){
printf("Failed to converge! %f \n", convergence);
fflush(stdout);
}
} }