From 757a4aa409233de7d450b7fba79a9294c78f0f06 Mon Sep 17 00:00:00 2001 From: austin Date: Sun, 15 Dec 2024 22:16:04 -0500 Subject: [PATCH] more pressurecell work --- .../c/src/fluid/sim/pressurecell/density.c | 58 ++++++--- .../c/src/fluid/sim/pressurecell/velocity.c | 110 ++++++++++++++---- 2 files changed, 130 insertions(+), 38 deletions(-) diff --git a/src/main/c/src/fluid/sim/pressurecell/density.c b/src/main/c/src/fluid/sim/pressurecell/density.c index 8d3f2f59..4b49c931 100644 --- a/src/main/c/src/fluid/sim/pressurecell/density.c +++ b/src/main/c/src/fluid/sim/pressurecell/density.c @@ -2,6 +2,7 @@ #include "fluid/sim/pressurecell/density.h" #include "fluid/sim/pressurecell/solver_consts.h" +#include "math/ode/multigrid_parallel.h" /** * Adds density from the delta buffer to this chunk @@ -26,29 +27,56 @@ LIBRARY_API void pressurecell_diffuse_density(Environment * environment, Chunk * int x, y, z; float * densityArr = chunk->d[CENTER_LOC]; float * densityTemp = chunk->dTempCache; - float D = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt; + float a = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt / (FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING); + float c = 1+6*a; + // float residual = 1; + // int iterations = 0; + // while(iterations < FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){ + // residual = solver_multigrid_parallel_iterate(densityTemp,densityArr,a,c); + // // fluid_grid2_set_bounds(environment,BOUND_SET_DENSITY_PHI,x); + // for(x = 1; x < DIM-1; x++){ + // for(y = 1; y < DIM-1; y++){ + // //diffuse back into the grid + // densityTemp[IX(0,x,y)] = 0; + // densityTemp[IX(DIM-1,x,y)] = 0; + // densityTemp[IX(x,0,y)] = 0; + // densityTemp[IX(x,DIM-1,y)] = 0; + // densityTemp[IX(x,y,0)] = 0; + // densityTemp[IX(x,y,DIM-1)] = 0; + // } + // } + // iterations++; + // } for(z = 1; z < DIM-1; z++){ for(y = 1; y < DIM-1; y++){ for(x = 1; x < DIM-1; x++){ densityTemp[IX(x,y,z)] = densityArr[IX(x,y,z)] + - densityArr[IX(x,y,z)] * -6 * D + - densityArr[IX(x-1,y,z)] * D + - densityArr[IX(x+1,y,z)] * D + - densityArr[IX(x,y-1,z)] * D + - densityArr[IX(x,y+1,z)] * D + - densityArr[IX(x,y,z-1)] * D + - densityArr[IX(x,y,z+1)] * D + ( + densityArr[IX(x,y,z)] * -6 + + ( + densityArr[IX(x-1,y,z)] + + densityArr[IX(x+1,y,z)] + + densityArr[IX(x,y-1,z)] + + densityArr[IX(x,y+1,z)] + + densityArr[IX(x,y,z-1)] + + densityArr[IX(x,y,z+1)] + ) + ) * a ; } } } - // for(z = 1; z < DIM-1; z++){ - // for(y = 1; y < DIM-1; y++){ - // for(x = 1; x < DIM-1; x++){ - // densityArr[IX(x,y,z)] = densityTemp[IX(x,y,z)]; - // } - // } - // } + for(x = 1; x < DIM-1; x++){ + for(y = 1; y < DIM-1; y++){ + //diffuse back into the grid + densityTemp[IX(1,x,y)] = densityTemp[IX(1,x,y)] + densityTemp[IX(0,x,y)]; + densityTemp[IX(DIM-2,x,y)] = densityTemp[IX(DIM-2,x,y)] + densityTemp[IX(DIM-1,x,y)]; + densityTemp[IX(x,1,y)] = densityTemp[IX(x,1,y)] + densityTemp[IX(x,0,y)]; + densityTemp[IX(x,DIM-2,y)] = densityTemp[IX(x,DIM-2,y)] + densityTemp[IX(x,DIM-1,y)]; + densityTemp[IX(x,y,1)] = densityTemp[IX(x,y,1)] + densityTemp[IX(x,y,0)]; + densityTemp[IX(x,y,DIM-2)] = densityTemp[IX(x,y,DIM-2)] + densityTemp[IX(x,y,DIM-1)]; + } + } } /** diff --git a/src/main/c/src/fluid/sim/pressurecell/velocity.c b/src/main/c/src/fluid/sim/pressurecell/velocity.c index 94542dc1..bdeb6cbf 100644 --- a/src/main/c/src/fluid/sim/pressurecell/velocity.c +++ b/src/main/c/src/fluid/sim/pressurecell/velocity.c @@ -3,6 +3,7 @@ #include "fluid/queue/chunk.h" #include "fluid/sim/pressurecell/velocity.h" #include "fluid/sim/pressurecell/solver_consts.h" +#include "math/ode/multigrid_parallel.h" /** @@ -59,18 +60,73 @@ LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk float * uTemp = chunk->uTempCache; float * vTemp = chunk->vTempCache; float * wTemp = chunk->wTempCache; - float D = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt / (FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING); + float a = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt / (FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING); + float c = 1+6*a; + // float residual = 1; + // int iterations = 0; + // while(iterations < FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){ + // residual = solver_multigrid_parallel_iterate(uArr,uTemp,a,c); + // // fluid_grid2_set_bounds(environment,BOUND_SET_DENSITY_PHI,x); + // for(x = 1; x < DIM-1; x++){ + // for(y = 1; y < DIM-1; y++){ + // //diffuse back into the grid + // uArr[IX(0,x,y)] = 0; + // uArr[IX(DIM-1,x,y)] = 0; + // uArr[IX(x,0,y)] = 0; + // uArr[IX(x,DIM-1,y)] = 0; + // uArr[IX(x,y,0)] = 0; + // uArr[IX(x,y,DIM-1)] = 0; + // } + // } + // iterations++; + // } + // while(iterations < FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){ + // residual = solver_multigrid_parallel_iterate(vArr,vTemp,a,c); + // for(x = 1; x < DIM-1; x++){ + // for(y = 1; y < DIM-1; y++){ + // //diffuse back into the grid + // vArr[IX(0,x,y)] = 0; + // vArr[IX(DIM-1,x,y)] = 0; + // vArr[IX(x,0,y)] = 0; + // vArr[IX(x,DIM-1,y)] = 0; + // vArr[IX(x,y,0)] = 0; + // vArr[IX(x,y,DIM-1)] = 0; + // } + // } + // // fluid_grid2_set_bounds(environment,BOUND_SET_DENSITY_PHI,x); + // iterations++; + // } + // while(iterations < FLUID_GRID2_SOLVER_MULTIGRID_MAX_ITERATIONS && (residual > FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE || residual < -FLUID_GRID2_SOLVER_MULTIGRID_TOLERANCE)){ + // residual = solver_multigrid_parallel_iterate(wArr,wTemp,a,c); + // // fluid_grid2_set_bounds(environment,BOUND_SET_DENSITY_PHI,x); + // for(x = 1; x < DIM-1; x++){ + // for(y = 1; y < DIM-1; y++){ + // //diffuse back into the grid + // wArr[IX(0,x,y)] = 0; + // wArr[IX(DIM-1,x,y)] = 0; + // wArr[IX(x,0,y)] = 0; + // wArr[IX(x,DIM-1,y)] = 0; + // wArr[IX(x,y,0)] = 0; + // wArr[IX(x,y,DIM-1)] = 0; + // } + // } + // iterations++; + // } for(z = 1; z < DIM-1; z++){ for(y = 1; y < DIM-1; y++){ for(x = 1; x < DIM-1; x++){ uArr[IX(x,y,z)] = uTemp[IX(x,y,z)] + - uTemp[IX(x,y,z)] * -6 * D + - uTemp[IX(x-1,y,z)] * D + - uTemp[IX(x+1,y,z)] * D + - uTemp[IX(x,y-1,z)] * D + - uTemp[IX(x,y+1,z)] * D + - uTemp[IX(x,y,z-1)] * D + - uTemp[IX(x,y,z+1)] * D + ( + uTemp[IX(x,y,z)] * -6 + + ( + uTemp[IX(x-1,y,z)] + + uTemp[IX(x+1,y,z)] + + uTemp[IX(x,y-1,z)] + + uTemp[IX(x,y+1,z)] + + uTemp[IX(x,y,z-1)] + + uTemp[IX(x,y,z+1)] + ) + ) * c ; } } @@ -80,13 +136,17 @@ LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk for(y = 1; y < DIM-1; y++){ for(x = 1; x < DIM-1; x++){ vArr[IX(x,y,z)] = vTemp[IX(x,y,z)] + - vTemp[IX(x,y,z)] * -6 * D + - vTemp[IX(x-1,y,z)] * D + - vTemp[IX(x+1,y,z)] * D + - vTemp[IX(x,y-1,z)] * D + - vTemp[IX(x,y+1,z)] * D + - vTemp[IX(x,y,z-1)] * D + - vTemp[IX(x,y,z+1)] * D + ( + vTemp[IX(x,y,z)] * -6 + + ( + vTemp[IX(x-1,y,z)] + + vTemp[IX(x+1,y,z)] + + vTemp[IX(x,y-1,z)] + + vTemp[IX(x,y+1,z)] + + vTemp[IX(x,y,z-1)] + + vTemp[IX(x,y,z+1)] + ) + ) * c ; } } @@ -96,13 +156,17 @@ LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk for(y = 1; y < DIM-1; y++){ for(x = 1; x < DIM-1; x++){ wArr[IX(x,y,z)] = wTemp[IX(x,y,z)] + - wTemp[IX(x,y,z)] * -6 * D + - wTemp[IX(x-1,y,z)] * D + - wTemp[IX(x+1,y,z)] * D + - wTemp[IX(x,y-1,z)] * D + - wTemp[IX(x,y+1,z)] * D + - wTemp[IX(x,y,z-1)] * D + - wTemp[IX(x,y,z+1)] * D + ( + wTemp[IX(x,y,z)] * -6 + + ( + wTemp[IX(x-1,y,z)] + + wTemp[IX(x+1,y,z)] + + wTemp[IX(x,y-1,z)] + + wTemp[IX(x,y+1,z)] + + wTemp[IX(x,y,z-1)] + + wTemp[IX(x,y,z+1)] + ) + ) * c ; } } @@ -125,7 +189,7 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk * float s0, s1, t0, t1, u0, u1; float interpolatedU, interpolatedV, interpolatedW; float magnitude; - float interpConst = environment->consts.dt / FLUID_PRESSURECELL_SPACING; + float interpConst = environment->consts.dt / (FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING); for(y = 1; y < DIM-1; y++){ for(z = 1; z < DIM-1; z++){ for(x = 1; x < DIM-1; x++){