From e9c79046d510155753f183c7d3f0ea9e81291e2a Mon Sep 17 00:00:00 2001 From: austin Date: Sun, 15 Dec 2024 17:05:11 -0500 Subject: [PATCH] fluid work --- .vscode/settings.json | 3 +- .../fluid/sim/pressurecell/normalization.h | 5 + .../fluid/sim/pressurecell/solver_consts.h | 20 + .../fluid/sim/pressurecell/velocity.h | 7 +- .../c/src/fluid/sim/pressurecell/density.c | 35 +- .../fluid/sim/pressurecell/normalization.c | 8 +- .../c/src/fluid/sim/pressurecell/pressure.c | 194 +++++-- .../src/fluid/sim/pressurecell/pressurecell.c | 42 +- .../c/src/fluid/sim/pressurecell/velocity.c | 205 ++++++- .../fluid/sim/pressurecell/add_source_tests.c | 6 +- .../fluid/sim/pressurecell/advection_tests.c | 2 +- .../fluid/sim/pressurecell/divergence_tests.c | 124 +++- .../c/fluid/sim/pressurecell/pressure_tests.c | 541 +++++++++++++++++- .../fluid/sim/pressurecell/projection_tests.c | 443 ++++++++++++++ 14 files changed, 1535 insertions(+), 100 deletions(-) create mode 100644 src/test/c/fluid/sim/pressurecell/projection_tests.c diff --git a/.vscode/settings.json b/.vscode/settings.json index 9661ac3b..7a90885a 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -52,6 +52,7 @@ "diffusion_ode.h": "c", "pressurecell.h": "c", "pressure.h": "c", - "tracking.h": "c" + "tracking.h": "c", + "multigrid_navier_stokes.h": "c" } } \ No newline at end of file diff --git a/src/main/c/includes/fluid/sim/pressurecell/normalization.h b/src/main/c/includes/fluid/sim/pressurecell/normalization.h index 8c7fc223..3fc1f39f 100644 --- a/src/main/c/includes/fluid/sim/pressurecell/normalization.h +++ b/src/main/c/includes/fluid/sim/pressurecell/normalization.h @@ -6,6 +6,11 @@ #include "fluid/queue/chunk.h" #include "fluid/queue/chunkmask.h" +/** + * Calculates the expected density and pressure + */ +LIBRARY_API void fluid_pressurecell_calculate_expected_intake(Environment * env, Chunk * chunk); + /** * Calculates the ratio to normalize the chunk by */ diff --git a/src/main/c/includes/fluid/sim/pressurecell/solver_consts.h b/src/main/c/includes/fluid/sim/pressurecell/solver_consts.h index 37aa1bb6..e1f70587 100644 --- a/src/main/c/includes/fluid/sim/pressurecell/solver_consts.h +++ b/src/main/c/includes/fluid/sim/pressurecell/solver_consts.h @@ -52,5 +52,25 @@ */ #define FLUID_PRESSURECELL_DIV_PRESSURE_CONST 5.0f +/** + * The number of times to relax most solvers + */ +#define FLUID_PRESSURECELL_SOLVER_MULTIGRID_MAX_ITERATIONS 5 + +/** + * The tolerance to shoot for when approximating pressure + */ +#define FLUID_PRESSURECELL_PROJECTION_CONVERGENCE_TOLERANCE 0.01f + +/** + * The minimum pressure allowed + */ +// #define FLUID_PRESSURECELL_MIN_PRESSURE -100000.0f + +/** + * The maximum pressure allowed + */ +// #define FLUID_PRESSURECELL_MAX_PRESSURE 100000.0f + #endif \ No newline at end of file diff --git a/src/main/c/includes/fluid/sim/pressurecell/velocity.h b/src/main/c/includes/fluid/sim/pressurecell/velocity.h index 0a90691c..5f7b6c5d 100644 --- a/src/main/c/includes/fluid/sim/pressurecell/velocity.h +++ b/src/main/c/includes/fluid/sim/pressurecell/velocity.h @@ -31,7 +31,12 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk * /** * Interpolates between the advected velocity and the previous frame's velocity by the pressure divergence amount */ -LIBRARY_API void pressurecell_interpolate_velocity(Environment * environment, Chunk * chunk); +LIBRARY_API double pressurecell_project_velocity(Environment * environment, Chunk * chunk); + +/** + * Copy temp velocities to next frame +*/ +LIBRARY_API void pressurecell_copy_for_next_frame(Environment * environment, Chunk * chunk); #endif \ No newline at end of file diff --git a/src/main/c/src/fluid/sim/pressurecell/density.c b/src/main/c/src/fluid/sim/pressurecell/density.c index 9b00f1b4..4a8ce307 100644 --- a/src/main/c/src/fluid/sim/pressurecell/density.c +++ b/src/main/c/src/fluid/sim/pressurecell/density.c @@ -42,6 +42,13 @@ LIBRARY_API void pressurecell_diffuse_density(Environment * environment, Chunk * } } } + // 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)]; + // } + // } + // } } /** @@ -58,14 +65,35 @@ LIBRARY_API void pressurecell_advect_density(Environment * environment, Chunk * float xp, yp, zp; float s0, s1, t0, t1, u0, u1; float interpolated; + float vecU, vecV, vecW; for(y = 1; y < DIM-1; y++){ //TODO: eventually skip y levels if there is no density to advect for(z = 1; z < DIM-1; z++){ for(x = 1; x < DIM-1; x++){ //calculate the real (float) position we are at - xp = x - u[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; - yp = y - v[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; - zp = z - w[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; + vecU = u[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; + vecV = v[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; + vecW = w[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; + if(vecU > 0.999f){ + vecU = 0.999f; + } else if(vecU < -0.999f){ + vecU = -0.999f; + } + if(vecV > 0.999f){ + vecV = 0.999f; + } else if(vecV < -0.999f){ + vecV = -0.999f; + } + if(vecW > 0.999f){ + vecW = 0.999f; + } else if(vecW < -0.999f){ + vecW = -0.999f; + } + + //calculate percentage to pull from existing vs other + xp = x - vecU; + yp = y - vecV; + zp = z - vecW; //clamp to border x0 = xp; @@ -98,6 +126,7 @@ LIBRARY_API void pressurecell_advect_density(Environment * environment, Chunk * printf("%d %d %d\n", x1, y1, z1); printf("%f %f %f\n", xp, yp, zp); printf("%f %f %f\n", u[IX(x,y,z)], v[IX(x,y,z)], w[IX(x,y,z)]); + printf("%f %f %f\n", vecU, vecV, vecW); printf("%f\n", environment->consts.dt); fflush(stdout); } diff --git a/src/main/c/src/fluid/sim/pressurecell/normalization.c b/src/main/c/src/fluid/sim/pressurecell/normalization.c index 8c9c8b57..3c154139 100644 --- a/src/main/c/src/fluid/sim/pressurecell/normalization.c +++ b/src/main/c/src/fluid/sim/pressurecell/normalization.c @@ -38,8 +38,12 @@ LIBRARY_API void fluid_pressurecell_calculate_normalization_ratio(Environment * } double expected = chunk->pressureCellData.densitySum; env->state.existingDensity = env->state.existingDensity + expected; - double normalizationRatio = expected / sum; - chunk->pressureCellData.densitySum = normalizationRatio; + if(sum > 0){ + double normalizationRatio = expected / sum; + chunk->pressureCellData.densitySum = normalizationRatio; + } else { + chunk->pressureCellData.densitySum = expected; + } } /** diff --git a/src/main/c/src/fluid/sim/pressurecell/pressure.c b/src/main/c/src/fluid/sim/pressurecell/pressure.c index d7208334..2b31cdbb 100644 --- a/src/main/c/src/fluid/sim/pressurecell/pressure.c +++ b/src/main/c/src/fluid/sim/pressurecell/pressure.c @@ -1,7 +1,8 @@ - +#include #include "fluid/sim/pressurecell/pressure.h" #include "fluid/sim/pressurecell/solver_consts.h" +#include "math/ode/multigrid_parallel.h" @@ -22,58 +23,155 @@ LIBRARY_API void pressurecell_approximate_pressure(Environment * environment, Ch //temporary caches float * pressureCache = chunk->pressureCache[CENTER_LOC]; float * pressureTemp = chunk->pressureTempCache; + float * phi0 = chunk->dTempCache; //consts/vars float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING; float du, dv, dw; float newPressure; float dt = environment->consts.dt; + // for(z = 1; z < DIM-1; z++){ + // for(y = 1; y < DIM-1; y++){ + // for(x = 1; x < DIM-1; x++){ + // //divergence part + // //positive value means inflow + // //negative value means outflow + // float newDiv = + // ( + // uArr[IX(x+1,y,z)] - uArr[IX(x-1,y,z)] + + // vArr[IX(x,y+1,z)] - vArr[IX(x,y-1,z)] + + // wArr[IX(x,y,z+1)] - wArr[IX(x,y,z-1)] + // ) / (2 * gridSpacing); + // float finalDiv = divArr[IX(x,y,z)]; + + + // //poisson stencil + // float stencil = ( + // -6 * pressureCache[IX(x,y,z)] + + // ( + // pressureCache[IX(x+1,y,z)] + + // pressureCache[IX(x-1,y,z)] + + // pressureCache[IX(x,y+1,z)] + + // pressureCache[IX(x,y-1,z)] + + // pressureCache[IX(x,y,z+1)] + + // pressureCache[IX(x,y,z-1)] + // ) + // ) / gridSpacing; + // float residual = stencil - FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER * finalDiv; + // // if(residual > 0){ + // // printf("%f \n", finalDiv); + // // printf("%f \n", stencil); + // // printf("%f \n", residual); + // // } + + // //divergence caused by static outflow due to diffusion + // float outflowDiv = -FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 6 * dt; + + // //compute the new pressure value + // // newPressure = pressureCache[IX(x,y,z)] + residual; + // newPressure = FLUID_PRESSURECELL_DIV_PRESSURE_CONST * (newDiv + outflowDiv); + // // if(newPressure > 0){ + // // printf("%f \n",newPressure); + // // } + // pressureTemp[IX(x,y,z)] = pressureCache[IX(x,y,z)] + newPressure; + // pressureTemp[IX(x,y,z)] = fmax(FLUID_PRESSURECELL_MIN_PRESSURE,fmin(FLUID_PRESSURECELL_MAX_PRESSURE,pressureTemp[IX(x,y,z)])); + // } + // } + // } + + //setup multigrid for(z = 1; z < DIM-1; z++){ for(y = 1; y < DIM-1; y++){ for(x = 1; x < DIM-1; x++){ - //divergence part - //positive value means inflow - //negative value means outflow - float newDiv = - ( - uArr[IX(x+1,y,z)] - uArr[IX(x-1,y,z)] + - vArr[IX(x+1,y,z)] - vArr[IX(x-1,y,z)] + - wArr[IX(x+1,y,z)] - wArr[IX(x-1,y,z)] - ) / (2 * gridSpacing); - float finalDiv = divArr[IX(x,y,z)]; - - - //poisson stencil - float stencil = ( - -6 * pressureCache[IX(x,y,z)] + - ( - pressureCache[IX(x+1,y,z)] + - pressureCache[IX(x-1,y,z)] + - pressureCache[IX(x,y+1,z)] + - pressureCache[IX(x,y-1,z)] + - pressureCache[IX(x,y,z+1)] + - pressureCache[IX(x,y,z-1)] - ) - ) / gridSpacing; - float residual = stencil - FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER * finalDiv; - // if(residual > 0){ - // printf("%f \n", finalDiv); - // printf("%f \n", stencil); - // printf("%f \n", residual); - // } - - //divergence caused by static outflow due to diffusion - float outflowDiv = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 6 * dt; - - //compute the new pressure value - newPressure = pressureCache[IX(x,y,z)] + residual; - // newPressure = -FLUID_PRESSURECELL_DIV_PRESSURE_CONST * (newDiv + outflowDiv); - // if(newPressure > 0){ - // printf("%f \n",newPressure); - // } - pressureTemp[IX(x,y,z)] = newPressure; + phi0[IX(x,y,z)] = divArr[IX(x,y,z)]; + pressureTemp[IX(x,y,z)] = 0; } } } + for(x = 0; x < DIM; x++){ + for(y = 0; y < DIM; y++){ + //pressure borders + //essentially, if the pressure builds up next to an edge, we don't have to have a 0 pressure area right next to it on the edge itself + //there are two values that should potentially be set to here + //either, same pressure as voxel in normal direction if this edge is actually an edge + //otherwise, set to the pressure of the neighboring chunk + pressureTemp[IX(0,x,y)] = pressureCache[IX(0,x,y)]; + pressureTemp[IX(DIM-1,x,y)] = pressureCache[IX(DIM-1,x,y)]; + pressureTemp[IX(x,0,y)] = pressureCache[IX(x,0,y)]; + pressureTemp[IX(x,DIM-1,y)] = pressureCache[IX(x,DIM-1,y)]; + pressureTemp[IX(x,y,0)] = pressureCache[IX(x,y,0)]; + pressureTemp[IX(x,y,DIM-1)] = pressureCache[IX(x,y,DIM-1)]; + + //divergence borders + phi0[IX(0,x,y)] = divArr[IX(0,x,y)]; + phi0[IX(DIM-1,x,y)] = divArr[IX(DIM-1,x,y)]; + phi0[IX(x,0,y)] = divArr[IX(x,0,y)]; + phi0[IX(x,DIM-1,y)] = divArr[IX(x,DIM-1,y)]; + phi0[IX(x,y,0)] = divArr[IX(x,y,0)]; + phi0[IX(x,y,DIM-1)] = divArr[IX(x,y,DIM-1)]; + } + } + float a = 1; + float c = 6; + chunk->projectionIterations = 0; + chunk->projectionResidual = 1; + while(chunk->projectionIterations < FLUID_PRESSURECELL_SOLVER_MULTIGRID_MAX_ITERATIONS && (chunk->projectionResidual > FLUID_PRESSURECELL_PROJECTION_CONVERGENCE_TOLERANCE || chunk->projectionResidual < -FLUID_PRESSURECELL_PROJECTION_CONVERGENCE_TOLERANCE)){ + chunk->projectionResidual = solver_multigrid_parallel_iterate(pressureTemp,phi0,a,c); + // for(z = 1; z < DIM-1; z++){ + // for(y = 1; y < DIM-1; y++){ + // for(x = 1; x < DIM-1; x++){ + // pressureTemp[IX(x,y,z)] = fmax(FLUID_PRESSURECELL_MIN_PRESSURE,fmin(FLUID_PRESSURECELL_MAX_PRESSURE,pressureTemp[IX(x,y,z)])); + // } + // } + // } + //essentially, if the pressure builds up next to an edge, we don't have to have a 0 pressure area right next to it on the edge itself + //there are two values that should potentially be set to here + //either, same pressure as voxel in normal direction if this edge is actually an edge + //otherwise, set to the pressure of the neighboring chunk + for(x = 0; x < DIM; x++){ + for(y = 0; y < DIM; y++){ + //pressure borders + pressureTemp[IX(0,x,y)] = 0; + pressureTemp[IX(DIM-1,x,y)] = 0; + pressureTemp[IX(x,0,y)] = 0; + pressureTemp[IX(x,DIM-1,y)] = 0; + pressureTemp[IX(x,y,0)] = 0; + pressureTemp[IX(x,y,DIM-1)] = 0; + // pressureTemp[IX(0,x,y)] = pressureTemp[IX(1,x,y)]; + // pressureTemp[IX(DIM-1,x,y)] = pressureTemp[IX(DIM-2,x,y)]; + // pressureTemp[IX(x,0,y)] = pressureTemp[IX(x,1,y)]; + // pressureTemp[IX(x,DIM-1,y)] = pressureTemp[IX(x,DIM-2,y)]; + // pressureTemp[IX(x,y,0)] = pressureTemp[IX(x,y,1)]; + // pressureTemp[IX(x,y,DIM-1)] = pressureTemp[IX(x,y,DIM-2)]; + } + } + chunk->projectionIterations++; + } + // for(z = 1; z < DIM-1; z++){ + // for(y = 1; y < DIM-1; y++){ + // for(x = 1; x < DIM-1; x++){ + // //check for NaNs + // if(pressureTemp[IX(x,y,z)] != pressureTemp[IX(x,y,z)]){ + // pressureTemp[IX(x,y,z)] = 0; + // } + // // pressureTemp[IX(x,y,z)] = pressureTemp[IX(x,y,z)] / 10.0f; + // pressureTemp[IX(x,y,z)] = fmax(FLUID_PRESSURECELL_MIN_PRESSURE,fmin(FLUID_PRESSURECELL_MAX_PRESSURE,pressureTemp[IX(x,y,z)])); + // } + // } + // } + + //do NOT zero out pressure on edges + //this will cause the fluid to advect into the walls + // for(x = 0; x < DIM; x++){ + // for(y = 0; y < DIM; y++){ + // //pressure borders + // pressureTemp[IX(0,x,y)] = 0; + // pressureTemp[IX(DIM-1,x,y)] = 0; + // pressureTemp[IX(x,0,y)] = 0; + // pressureTemp[IX(x,DIM-1,y)] = 0; + // pressureTemp[IX(x,y,0)] = 0; + // pressureTemp[IX(x,y,DIM-1)] = 0; + // } + // } } /** @@ -120,7 +218,17 @@ LIBRARY_API void pressurecell_approximate_divergence(Environment * environment, // dw = 0; // } newDivergence = du+dv+dw; - divArr[IX(x,y,z)] = divArr[IX(x,y,z)] + newDivergence - FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER * divArr[IX(x,y,z)] + outflowDiv; + // divArr[IX(x,y,z)] = divArr[IX(x,y,z)] + newDivergence - FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER * divArr[IX(x,y,z)] + outflowDiv; + divArr[IX(x,y,z)] = newDivergence; + if(newDivergence > 3 || newDivergence < -3){ + printf("Invalid divergence! \n"); + printf("%f \n",newDivergence); + printf("%f %f \n", uArr[IX(x+1,y,z)], uArr[IX(x-1,y,z)]); + printf("%f %f \n", vArr[IX(x,y+1,z)], vArr[IX(x,y-1,z)]); + printf("%f %f \n", wArr[IX(x,y,z+1)], wArr[IX(x,y,z-1)]); + printf("%f %f %f \n", du, dv, dw); + printf("\n"); + } //store pressure value from this frame presureCache[IX(x,y,z)] = pressureTemp[IX(x,y,z)]; diff --git a/src/main/c/src/fluid/sim/pressurecell/pressurecell.c b/src/main/c/src/fluid/sim/pressurecell/pressurecell.c index 05660c92..72620cff 100644 --- a/src/main/c/src/fluid/sim/pressurecell/pressurecell.c +++ b/src/main/c/src/fluid/sim/pressurecell/pressurecell.c @@ -13,6 +13,7 @@ #include "fluid/sim/pressurecell/density.h" #include "fluid/sim/pressurecell/pressure.h" #include "fluid/sim/pressurecell/solver_consts.h" +#include "fluid/sim/pressurecell/normalization.h" #include "fluid/sim/pressurecell/velocity.h" @@ -56,37 +57,48 @@ LIBRARY_API void fluid_pressurecell_simulate( // // Velocity phase // + + //approximate pressure first using the values from last frame + //this allows us to guarantee that we're using the divergence from neighbors (eventually) + for(int i = 0; i < numChunks; i++){ + Chunk * currentChunk = chunks[i]; + pressurecell_approximate_pressure(environment,currentChunk); + } + for(int i = 0; i < numChunks; i++){ + Chunk * currentChunk = chunks[i]; + pressurecell_project_velocity(environment,currentChunk); + } + // printf("grav\n"); + // fflush(stdout); for(int i = 0; i < numChunks; i++){ Chunk * currentChunk = chunks[i]; pressurecell_add_gravity(environment,currentChunk); } + // printf("+vel\n"); + // fflush(stdout); for(int i = 0; i < numChunks; i++){ Chunk * currentChunk = chunks[i]; pressurecell_add_velocity(environment,currentChunk); } + // printf("diff vel\n"); + // fflush(stdout); for(int i = 0; i < numChunks; i++){ Chunk * currentChunk = chunks[i]; pressurecell_diffuse_velocity(environment,currentChunk); } + // printf("adv vel\n"); + // fflush(stdout); for(int i = 0; i < numChunks; i++){ Chunk * currentChunk = chunks[i]; pressurecell_advect_velocity(environment,currentChunk); } - for(int i = 0; i < numChunks; i++){ - Chunk * currentChunk = chunks[i]; - pressurecell_approximate_pressure(environment,currentChunk); - } - - for(int i = 0; i < numChunks; i++){ - Chunk * currentChunk = chunks[i]; - pressurecell_interpolate_velocity(environment,currentChunk); - } - + // printf("approx div\n"); + // fflush(stdout); for(int i = 0; i < numChunks; i++){ Chunk * currentChunk = chunks[i]; pressurecell_approximate_divergence(environment,currentChunk); @@ -98,18 +110,22 @@ LIBRARY_API void fluid_pressurecell_simulate( // // Density phase // - - + // printf("+dens\n"); + // fflush(stdout); for(int i = 0; i < numChunks; i++){ Chunk * currentChunk = chunks[i]; pressurecell_add_density(environment,currentChunk); } + // printf("diff dens\n"); + // fflush(stdout); for(int i = 0; i < numChunks; i++){ Chunk * currentChunk = chunks[i]; pressurecell_diffuse_density(environment,currentChunk); } + // printf("adv dens\n"); + // fflush(stdout); for(int i = 0; i < numChunks; i++){ Chunk * currentChunk = chunks[i]; pressurecell_advect_density(environment,currentChunk); @@ -126,11 +142,13 @@ LIBRARY_API void fluid_pressurecell_simulate( for(int i = 0; i < numChunks; i++){ Chunk * currentChunk = chunks[i]; fluid_pressurecell_normalize_chunk(environment,currentChunk); + pressurecell_copy_for_next_frame(environment,currentChunk); fluid_pressurecell_clearArr(currentChunk->d0[CENTER_LOC]); fluid_pressurecell_clearArr(currentChunk->u0[CENTER_LOC]); fluid_pressurecell_clearArr(currentChunk->v0[CENTER_LOC]); fluid_pressurecell_clearArr(currentChunk->w0[CENTER_LOC]); fluid_pressurecell_clearArr(currentChunk->pressureTempCache); + fluid_pressurecell_clearArr(currentChunk->dTempCache); fluid_pressurecell_clearArr(currentChunk->uTempCache); fluid_pressurecell_clearArr(currentChunk->vTempCache); fluid_pressurecell_clearArr(currentChunk->wTempCache); diff --git a/src/main/c/src/fluid/sim/pressurecell/velocity.c b/src/main/c/src/fluid/sim/pressurecell/velocity.c index 6f1a0b31..adb6d9c6 100644 --- a/src/main/c/src/fluid/sim/pressurecell/velocity.c +++ b/src/main/c/src/fluid/sim/pressurecell/velocity.c @@ -30,15 +30,18 @@ LIBRARY_API void pressurecell_add_velocity(Environment * environment, Chunk * ch float * uArr = chunk->u[CENTER_LOC]; float * vArr = chunk->v[CENTER_LOC]; float * wArr = chunk->w[CENTER_LOC]; + float * uTemp = chunk->uTempCache; + float * vTemp = chunk->vTempCache; + float * wTemp = chunk->wTempCache; float * uSourceArr = chunk->u0[CENTER_LOC]; float * vSourceArr = chunk->v0[CENTER_LOC]; float * wSourceArr = chunk->w0[CENTER_LOC]; 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)] = uArr[IX(x,y,z)] + uSourceArr[IX(x,y,z)]; - vArr[IX(x,y,z)] = vArr[IX(x,y,z)] + vSourceArr[IX(x,y,z)]; - wArr[IX(x,y,z)] = wArr[IX(x,y,z)] + wSourceArr[IX(x,y,z)]; + uTemp[IX(x,y,z)] = uArr[IX(x,y,z)] + uSourceArr[IX(x,y,z)]; + vTemp[IX(x,y,z)] = vArr[IX(x,y,z)] + vSourceArr[IX(x,y,z)]; + wTemp[IX(x,y,z)] = wArr[IX(x,y,z)] + wSourceArr[IX(x,y,z)]; } } } @@ -221,9 +224,9 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk * // interpolatedW = interpolatedW / magnitude; // } - uArr[IX(x,y,z)] = interpolatedU; - vArr[IX(x,y,z)] = interpolatedV; - wArr[IX(x,y,z)] = interpolatedW; + uTemp[IX(x,y,z)] = interpolatedU; + vTemp[IX(x,y,z)] = interpolatedV; + wTemp[IX(x,y,z)] = interpolatedW; } } } @@ -232,27 +235,115 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk * /** * Interpolates between the advected velocity and the previous frame's velocity by the pressure divergence amount */ -LIBRARY_API void pressurecell_interpolate_velocity(Environment * environment, Chunk * chunk){ +LIBRARY_API double pressurecell_project_velocity(Environment * environment, Chunk * chunk){ int x, y, z; - float * presureCache = chunk->pressureTempCache; + float * pressureTemp = chunk->pressureTempCache; float * uArr = chunk->u[CENTER_LOC]; float * vArr = chunk->v[CENTER_LOC]; float * wArr = chunk->w[CENTER_LOC]; - float * uTemp = chunk->uTempCache; - float * vTemp = chunk->vTempCache; - float * wTemp = chunk->wTempCache; + // float * uTemp = chunk->uTempCache; + // float * vTemp = chunk->vTempCache; + // float * wTemp = chunk->wTempCache; //temporary caches float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING; - float du, dv, dw; float pressureDivergence; float magnitude; + float pressureDifferenceX, pressureDifferenceY, pressureDifferenceZ; + double maxMagnitude = 0; + //project for(y = 1; y < DIM-1; y++){ for(z = 1; z < DIM-1; z++){ for(x = 1; x < DIM-1; x++){ + /* + + lets say pressure is like: + 1 0 0 + x-1 x x+1 + + The pressure gradient becomes + -0.5 + + this pressure gradient is subtracted from the existing x velocity + so if the existing velocity was + 1 + + it is now + 1 - (-0.5) = 1.5 + + the higher pressure value has pushed the vector away from itself + + lets say pressure is like: + -300 500 1000 + x-1 x x+1 + + the pressure gradient becomes + -650 + + so when we modify the velocity of this field, we will now get + -649 + + Obviously this is really high, so we account for this by normalizing the velocity + if the vector was originally pushing up along y, now it is almost exclusively pushing along x + + */ + pressureDifferenceX = (pressureTemp[IX(x+1,y,z)] - pressureTemp[IX(x-1,y,z)]) / (gridSpacing * 2.0f); + pressureDifferenceY = (pressureTemp[IX(x,y+1,z)] - pressureTemp[IX(x,y-1,z)]) / (gridSpacing * 2.0f); + pressureDifferenceZ = (pressureTemp[IX(x,y,z+1)] - pressureTemp[IX(x,y,z-1)]) / (gridSpacing * 2.0f); + + //check for NaNs + if(pressureDifferenceX != pressureDifferenceX){ + pressureDifferenceX = 0; + } + if(pressureDifferenceY != pressureDifferenceY){ + pressureDifferenceY = 0; + } + if(pressureDifferenceZ != pressureDifferenceZ){ + pressureDifferenceZ = 0; + } + + //make sure the pressure gradient does not push the velocity into walls + if(x == 1 && pressureDifferenceX > 0){ + pressureDifferenceX = 0; + } + if(x == DIM-2 && pressureDifferenceX < 0){ + pressureDifferenceX = 0; + } + if(y == 1 && pressureDifferenceY > 0){ + pressureDifferenceY = 0; + } + if(y == DIM-2 && pressureDifferenceY < 0){ + pressureDifferenceY = 0; + } + if(z == 1 && pressureDifferenceZ > 0){ + pressureDifferenceZ = 0; + } + if(z == DIM-2 && pressureDifferenceZ < 0){ + pressureDifferenceZ = 0; + } + //project the pressure gradient onto the velocity field - uTemp[IX(x,y,z)] = uTemp[IX(x,y,z)] - (presureCache[IX(x+1,y,z)] - presureCache[IX(x-1,y,z)]) / (gridSpacing * 2.0f); - vTemp[IX(x,y,z)] = vTemp[IX(x,y,z)] - (presureCache[IX(x,y+1,z)] - presureCache[IX(x,y-1,z)]) / (gridSpacing * 2.0f); - wTemp[IX(x,y,z)] = wTemp[IX(x,y,z)] - (presureCache[IX(x,y,z+1)] - presureCache[IX(x,y,z-1)]) / (gridSpacing * 2.0f); + uArr[IX(x,y,z)] = uArr[IX(x,y,z)] - pressureDifferenceX; + vArr[IX(x,y,z)] = vArr[IX(x,y,z)] - pressureDifferenceY; + wArr[IX(x,y,z)] = wArr[IX(x,y,z)] - pressureDifferenceZ; + + float magnitude = sqrt(uArr[IX(x,y,z)] * uArr[IX(x,y,z)] + vArr[IX(x,y,z)] * vArr[IX(x,y,z)] + wArr[IX(x,y,z)] * wArr[IX(x,y,z)]); + // if(maxMagnitude < magnitude){ + // maxMagnitude = magnitude; + // } + + if(magnitude != magnitude || magnitude > 1000000){ + uArr[IX(x,y,z)] = 0; + vArr[IX(x,y,z)] = 0; + wArr[IX(x,y,z)] = 0; + } + + //normalize if the projection has pushed us wayyy out of bounds + //ie, large pressure differentials can create huge imbalances + if(magnitude > 1.0f){ + uArr[IX(x,y,z)] = uArr[IX(x,y,z)] / magnitude; + vArr[IX(x,y,z)] = vArr[IX(x,y,z)] / magnitude; + wArr[IX(x,y,z)] = wArr[IX(x,y,z)] / magnitude; + } // magnitude = sqrt(uTemp[IX(x,y,z)] * uTemp[IX(x,y,z)] + vTemp[IX(x,y,z)] * vTemp[IX(x,y,z)] + wTemp[IX(x,y,z)] * wTemp[IX(x,y,z)]); @@ -272,20 +363,90 @@ LIBRARY_API void pressurecell_interpolate_velocity(Environment * environment, Ch // } // if( - // uTemp[x,y,z] < -100.0f || uTemp[x,y,z] > 100.0f || - // vTemp[x,y,z] < -100.0f || vTemp[x,y,z] > 100.0f || - // wTemp[x,y,z] < -100.0f || wTemp[x,y,z] > 100.0f || - // magnitude < -1000 || magnitude > 1000 || - // pressureDivergence < -1000 || pressureDivergence > 1000 + // uArr[x,y,z] < -100.0f || uArr[x,y,z] > 100.0f || + // vArr[x,y,z] < -100.0f || vArr[x,y,z] > 100.0f || + // wArr[x,y,z] < -100.0f || wArr[x,y,z] > 100.0f + // || magnitude < -1000 || magnitude > 1000 + // // pressureDivergence < -1000 || pressureDivergence > 1000 // ){ // printf("pressure divergence thing is off!!\n"); // printf("%f \n", magnitude); // printf("%f \n", pressureDivergence); - // printf("%f %f %f \n", uTemp[IX(x,y,z)], vTemp[IX(x,y,z)], wTemp[IX(x,y,z)]); // printf("%f %f %f \n", uArr[IX(x,y,z)], vArr[IX(x,y,z)], wArr[IX(x,y,z)]); + // printf("%f %f \n", presureCache[IX(x+1,y,z)], presureCache[IX(x-1,y,z)]); + // printf("%f %f \n", presureCache[IX(x,y+1,z)], presureCache[IX(x,y-1,z)]); + // printf("%f %f \n", presureCache[IX(x,y,z+1)], presureCache[IX(x,y,z-1)]); // printf("\n"); // } } } } -} \ No newline at end of file + //normalize vector field + if(maxMagnitude > 1){ + for(y = 1; y < DIM-1; y++){ + for(z = 1; z < DIM-1; z++){ + for(x = 1; x < DIM-1; x++){ + + //project the pressure gradient onto the velocity field + uArr[IX(x,y,z)] = uArr[IX(x,y,z)] / maxMagnitude; + vArr[IX(x,y,z)] = vArr[IX(x,y,z)] / maxMagnitude; + wArr[IX(x,y,z)] = wArr[IX(x,y,z)] / maxMagnitude; + + //check for NaNs + if(uArr[IX(x,y,z)] != uArr[IX(x,y,z)]){ + uArr[IX(x,y,z)] = 0; + } + if(vArr[IX(x,y,z)] != vArr[IX(x,y,z)]){ + vArr[IX(x,y,z)] = 0; + } + if(wArr[IX(x,y,z)] != wArr[IX(x,y,z)]){ + wArr[IX(x,y,z)] = 0; + } + + if( + uArr[x,y,z] < -100.0f || uArr[x,y,z] > 100.0f || + vArr[x,y,z] < -100.0f || vArr[x,y,z] > 100.0f || + wArr[x,y,z] < -100.0f || wArr[x,y,z] > 100.0f + // || magnitude < -1000 || magnitude > 1000 + // pressureDivergence < -1000 || pressureDivergence > 1000 + ){ + printf("pressure divergence thing is off!!\n"); + printf("%f \n", magnitude); + printf("%f \n", pressureDivergence); + printf("%f %f %f \n", uArr[IX(x,y,z)], vArr[IX(x,y,z)], wArr[IX(x,y,z)]); + printf("%f %f \n", pressureTemp[IX(x+1,y,z)], pressureTemp[IX(x-1,y,z)]); + printf("%f %f \n", pressureTemp[IX(x,y+1,z)], pressureTemp[IX(x,y-1,z)]); + printf("%f %f \n", pressureTemp[IX(x,y,z+1)], pressureTemp[IX(x,y,z-1)]); + printf("\n"); + } + } + } + } + } + return maxMagnitude; +} + +/** + * Copy temp velocities to next frame +*/ +LIBRARY_API void pressurecell_copy_for_next_frame(Environment * environment, Chunk * chunk){ + int x, y, z; + float * presureCache = chunk->pressureTempCache; + float * uArr = chunk->u[CENTER_LOC]; + float * vArr = chunk->v[CENTER_LOC]; + float * wArr = chunk->w[CENTER_LOC]; + float * uTemp = chunk->uTempCache; + float * vTemp = chunk->vTempCache; + float * wTemp = chunk->wTempCache; + for(y = 1; y < DIM-1; y++){ + for(z = 1; z < DIM-1; z++){ + for(x = 1; x < DIM-1; x++){ + //project the pressure gradient onto the velocity field + uArr[IX(x,y,z)] = uTemp[IX(x,y,z)]; + vArr[IX(x,y,z)] = vTemp[IX(x,y,z)]; + wArr[IX(x,y,z)] = wTemp[IX(x,y,z)]; + } + } + } +} + diff --git a/src/test/c/fluid/sim/pressurecell/add_source_tests.c b/src/test/c/fluid/sim/pressurecell/add_source_tests.c index db785a2d..16ba8f78 100644 --- a/src/test/c/fluid/sim/pressurecell/add_source_tests.c +++ b/src/test/c/fluid/sim/pressurecell/add_source_tests.c @@ -81,17 +81,17 @@ int fluid_sim_pressurecell_add_source_test2(){ //test the result float expected, actual; expected = FLUID_PRESSURECELL_MAX_VELOCITY; - actual = currentChunk->u[CENTER_LOC][IX(4,4,4)]; + actual = currentChunk->uTempCache[IX(4,4,4)]; if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ rVal += assertEqualsFloat(expected,actual,"Failed to add u0! expected: %f actual: %f \n"); } expected = FLUID_PRESSURECELL_MAX_VELOCITY; - actual = currentChunk->v[CENTER_LOC][IX(4,4,4)]; + actual = currentChunk->vTempCache[IX(4,4,4)]; if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ rVal += assertEqualsFloat(expected,actual,"Failed to add v0! expected: %f actual: %f \n"); } expected = FLUID_PRESSURECELL_MAX_VELOCITY; - actual = currentChunk->w[CENTER_LOC][IX(4,4,4)]; + actual = currentChunk->wTempCache[IX(4,4,4)]; if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ rVal += assertEqualsFloat(expected,actual,"Failed to add w0! expected: %f actual: %f \n"); } diff --git a/src/test/c/fluid/sim/pressurecell/advection_tests.c b/src/test/c/fluid/sim/pressurecell/advection_tests.c index 34cc1e42..9230fca2 100644 --- a/src/test/c/fluid/sim/pressurecell/advection_tests.c +++ b/src/test/c/fluid/sim/pressurecell/advection_tests.c @@ -98,7 +98,7 @@ int fluid_sim_pressurecell_advection_test2(){ // cell that originall had values // expected = FLUID_PRESSURECELL_MAX_VELOCITY - FLUID_PRESSURECELL_MAX_VELOCITY * env->consts.dt; - actual = currentChunk->u0[CENTER_LOC][IX(4,4,4)]; + actual = currentChunk->uTempCache[IX(4,4,4)]; if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,4)! expected: %f actual: %f \n"); } diff --git a/src/test/c/fluid/sim/pressurecell/divergence_tests.c b/src/test/c/fluid/sim/pressurecell/divergence_tests.c index 46607263..22026184 100644 --- a/src/test/c/fluid/sim/pressurecell/divergence_tests.c +++ b/src/test/c/fluid/sim/pressurecell/divergence_tests.c @@ -57,7 +57,7 @@ int fluid_sim_pressurecell_divergence_test1(){ // // cell that originall had values // - expected = -3 * env->consts.dt; + expected = -3; actual = currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)]; if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ rVal += assertEqualsFloat(expected,actual,"Failed to calculate divergence correctly (4,4,4)! expected: %f actual: %f \n"); @@ -66,6 +66,126 @@ int fluid_sim_pressurecell_divergence_test1(){ return rVal; } +/** + * Testing pressure values + */ +int fluid_sim_pressurecell_divergence_test2(){ + printf("fluid_sim_pressurecell_divergence_test2\n"); + int rVal = 0; + Environment * env = fluid_environment_create(); + Chunk ** queue = NULL; + queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM); + int chunkCount = arrlen(queue); + + + + //setup chunk values + float deltaDensity = 0.01f; + Chunk * currentChunk = queue[0]; + + currentChunk->u[CENTER_LOC][IX(0,1,1)] = 0; + currentChunk->u[CENTER_LOC][IX(2,1,1)] = 1; + currentChunk->v[CENTER_LOC][IX(1,0,1)] = 0; + currentChunk->v[CENTER_LOC][IX(1,2,1)] = -1; + currentChunk->w[CENTER_LOC][IX(1,1,0)] = 0; + currentChunk->w[CENTER_LOC][IX(1,1,2)] = -1; + + //actually simulate + pressurecell_approximate_divergence(env,currentChunk); + + //test the result + float expected, actual; + + // + // cell that originall had values + // + expected = -0.5; + actual = currentChunk->divergenceCache[CENTER_LOC][IX(1,1,1)]; + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to calculate divergence correctly (1,1,1)! expected: %f actual: %f \n"); + } + + return rVal; +} + +/** + * Testing pressure values + */ +int fluid_sim_pressurecell_divergence_test3(){ + printf("fluid_sim_pressurecell_divergence_test3\n"); + int rVal = 0; + Environment * env = fluid_environment_create(); + Chunk ** queue = NULL; + queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM); + int chunkCount = arrlen(queue); + + + + //setup chunk values + float deltaDensity = 0.01f; + Chunk * currentChunk = queue[0]; + + currentChunk->u[CENTER_LOC][IX(1,1,1)] = 1; + currentChunk->v[CENTER_LOC][IX(1,1,1)] = 1; + currentChunk->w[CENTER_LOC][IX(1,1,1)] = 1; + + currentChunk->u[CENTER_LOC][IX(2,1,1)] = 1; + currentChunk->v[CENTER_LOC][IX(1,2,1)] = 1; + currentChunk->w[CENTER_LOC][IX(1,1,2)] = 1; + + //actually simulate + pressurecell_approximate_divergence(env,currentChunk); + + //test the result + float expected, actual; + int cx, cy, cz; + float u, v, w; + float sum; + + // + // cell that originall had values + // + cx = 1; + cy = 1; + cz = 1; + expected = currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy,cz)]; + if(expected < 1.5){ //we expect 1.5 velocity to leave this cell in one iteration + rVal++; + printf("Divergence calc failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("%f \n",expected); + printf("div at[%d,%d,%d] %f \n", cx,cy,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy,cz)]); + printf("div at[%d,%d,%d] %f \n", cx-1,cy,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx-1,cy,cz)]); + printf("div at[%d,%d,%d] %f \n", cx+1,cy,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx+1,cy,cz)]); + printf("div at[%d,%d,%d] %f \n", cx,cy-1,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy-1,cz)]); + printf("div at[%d,%d,%d] %f \n", cx,cy+1,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy+1,cz)]); + printf("div at[%d,%d,%d] %f \n", cx,cy,cz-1, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy,cz-1)]); + printf("div at[%d,%d,%d] %f \n", cx,cy,cz+1, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy,cz+1)]); + printf("\n"); + } + + cx = 2; + cy = 1; + cz = 1; + expected = currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy,cz)]; + if(expected != -0.5){//we expect 0.5 velocity to enter this cell in one iteration + rVal++; + printf("Divergence calc failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("%f \n",expected); + printf("div at[%d,%d,%d] %f \n", cx,cy,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy,cz)]); + printf("div at[%d,%d,%d] %f \n", cx-1,cy,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx-1,cy,cz)]); + printf("div at[%d,%d,%d] %f \n", cx+1,cy,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx+1,cy,cz)]); + printf("div at[%d,%d,%d] %f \n", cx,cy-1,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy-1,cz)]); + printf("div at[%d,%d,%d] %f \n", cx,cy+1,cz, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy+1,cz)]); + printf("div at[%d,%d,%d] %f \n", cx,cy,cz-1, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy,cz-1)]); + printf("div at[%d,%d,%d] %f \n", cx,cy,cz+1, currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy,cz+1)]); + printf("\n"); + } + + return rVal; +} + /** * Testing pressure values */ @@ -73,6 +193,8 @@ int fluid_sim_pressurecell_divergence_tests(int argc, char **argv){ int rVal = 0; rVal += fluid_sim_pressurecell_divergence_test1(); + rVal += fluid_sim_pressurecell_divergence_test2(); + rVal += fluid_sim_pressurecell_divergence_test3(); return rVal; } \ No newline at end of file diff --git a/src/test/c/fluid/sim/pressurecell/pressure_tests.c b/src/test/c/fluid/sim/pressurecell/pressure_tests.c index 747bb7f3..acd33769 100644 --- a/src/test/c/fluid/sim/pressurecell/pressure_tests.c +++ b/src/test/c/fluid/sim/pressurecell/pressure_tests.c @@ -16,7 +16,7 @@ /** * Error margin for tests */ -#define FLUID_PRESSURE_CELL_ERROR_MARGIN 0.00001f +#define FLUID_PRESSURE_CELL_ERROR_MARGIN 0.01f /** * Number of chunks @@ -39,14 +39,8 @@ int fluid_sim_pressurecell_pressure_test1(){ //setup chunk values float deltaDensity = 0.01f; Chunk * currentChunk = queue[0]; - currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)] = 3; - - currentChunk->pressureCache[CENTER_LOC][IX(3,4,4)] = 1; - currentChunk->pressureCache[CENTER_LOC][IX(5,4,4)] = -1; - currentChunk->pressureCache[CENTER_LOC][IX(4,3,4)] = 1; - currentChunk->pressureCache[CENTER_LOC][IX(4,5,4)] = -1; - currentChunk->pressureCache[CENTER_LOC][IX(4,4,3)] = 1; - currentChunk->pressureCache[CENTER_LOC][IX(4,4,5)] = -1; + currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)] = 1; + currentChunk->u[CENTER_LOC][IX(4,4,4)] = 1; //actually simulate pressurecell_approximate_pressure(env,currentChunk); @@ -57,10 +51,531 @@ int fluid_sim_pressurecell_pressure_test1(){ // // cell that originall had values // - expected = -FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER * currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)]; - actual = currentChunk->pressureTempCache[IX(4,4,4)]; + expected = currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)]; + actual = currentChunk->pressureTempCache[IX(4,4,4)] * 6 - ( + currentChunk->pressureTempCache[IX(3,4,4)] + + currentChunk->pressureTempCache[IX(5,4,4)] + + currentChunk->pressureTempCache[IX(4,3,4)] + + currentChunk->pressureTempCache[IX(4,5,4)] + + currentChunk->pressureTempCache[IX(4,4,3)] + + currentChunk->pressureTempCache[IX(4,4,5)] + ); if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ rVal += assertEqualsFloat(expected,actual,"Failed to advect density correctly (4,4,4)! expected: %f actual: %f \n"); + printf("%f \n", currentChunk->pressureTempCache[IX(4,4,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(3,4,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(5,4,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,3,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,5,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,4,3)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,4,5)]); + } + + + expected = currentChunk->divergenceCache[CENTER_LOC][IX(4,3,4)]; + actual = currentChunk->pressureTempCache[IX(4,3,4)] * 6 - ( + currentChunk->pressureTempCache[IX(3,3,4)] + + currentChunk->pressureTempCache[IX(5,3,4)] + + currentChunk->pressureTempCache[IX(4,2,4)] + + currentChunk->pressureTempCache[IX(4,4,4)] + + currentChunk->pressureTempCache[IX(4,3,3)] + + currentChunk->pressureTempCache[IX(4,3,5)] + ); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to advect density correctly (4,3,4)! expected: %f actual: %f \n"); + printf("%f \n", currentChunk->pressureTempCache[IX(4,3,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(3,3,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(5,3,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,2,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,4,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,3,3)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,3,5)]); + } + + return rVal; +} + +/** + * Testing pressure values + */ +int fluid_sim_pressurecell_pressure_test2(){ + printf("fluid_sim_pressurecell_pressure_test2\n"); + int rVal = 0; + Environment * env = fluid_environment_create(); + Chunk ** queue = NULL; + queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM); + int chunkCount = arrlen(queue); + + + + //setup chunk values + float deltaDensity = 0.01f; + Chunk * currentChunk = queue[0]; + currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)] = 1; + currentChunk->u[CENTER_LOC][IX(4,4,4)] = 1; + + //actually simulate + pressurecell_approximate_pressure(env,currentChunk); + + //test the result + float expected, actual; + + // + // cell that originall had values + // + int cx, cy, cz; + cx = 1; + cy = 1; + cz = 1; + expected = currentChunk->divergenceCache[CENTER_LOC][IX(cx,cy,cz)]; + actual = currentChunk->pressureTempCache[IX(cx,cy,cz)] * 6 - ( + currentChunk->pressureTempCache[IX(cx-1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx+1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx,cy-1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy+1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy,cz-1)] + + currentChunk->pressureTempCache[IX(cx,cy,cz+1)] + ); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to approximate pressure correctly! expected: %f actual: %f \n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("%f \n", currentChunk->pressureTempCache[IX(cx,cy,cz)]); + printf("%f \n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("%f \n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf("%f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("%f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf("%f \n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("%f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + } + + return rVal; +} + +/** + * Testing pressure values + */ +int fluid_sim_pressurecell_pressure_test3(){ + printf("fluid_sim_pressurecell_pressure_test3\n"); + int rVal = 0; + Environment * env = fluid_environment_create(); + Chunk ** queue = NULL; + queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM); + int chunkCount = arrlen(queue); + + + + //setup chunk values + float deltaDensity = 0.01f; + Chunk * currentChunk = queue[0]; + currentChunk->u[CENTER_LOC][IX(1,1,1)] = 1; + currentChunk->v[CENTER_LOC][IX(1,1,1)] = 1; + currentChunk->w[CENTER_LOC][IX(1,1,1)] = 1; + + currentChunk->u[CENTER_LOC][IX(2,1,1)] = 1; + currentChunk->v[CENTER_LOC][IX(1,2,1)] = 1; + currentChunk->w[CENTER_LOC][IX(1,1,2)] = 1; + + //divergence at 1,1,1 should be 1.5 + pressurecell_approximate_divergence(env,currentChunk); + //divergence at 2,1,1 should be -0.5 + + + + //actually simulate + pressurecell_approximate_pressure(env,currentChunk); + + //test the result + float expected, actual; + + // + // cell that originall had values + // + int cx, cy, cz; + cx = 1; + cy = 1; + cz = 1; + //essentially this is stating that we expect neighbors of this cell to give out 0.21f more velocity than they currently will this frame + //ergo, we should subtract that from this cell in order to prevent velocity compressibility + expected = 0.21; + actual = currentChunk->pressureTempCache[IX(cx,cy,cz)]; + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to approximate pressure correctly! expected: %f actual: %f \n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("Residual: %f\n", currentChunk->projectionResidual); + printf("\n\n\n"); + } + + + + cx = 2; + cy = 1; + cz = 1; + //essentially this is stating that we expect neighbors of this cell to give out -0.5 less velocity than they currently will this frame + //ergo, we should subtract that from this cell in order to prevent velocity compressibility + expected = -0.5f; + actual = currentChunk->pressureTempCache[IX(cx,cy,cz)] * 6 - ( + currentChunk->pressureTempCache[IX(cx-1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx+1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx,cy-1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy+1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy,cz-1)] + + currentChunk->pressureTempCache[IX(cx,cy,cz+1)] + ); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to approximate pressure correctly! expected: %f actual: %f \n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("Residual: %f\n", currentChunk->projectionResidual); + printf("\n\n\n"); + } + + cx = 3; + cy = 1; + cz = 1; + //essentially this is stating that we expect neighbors of this cell to give out -0.462 more velocity than they currently will this frame + //ergo, we should subtract that from this cell in order to prevent velocity compressibility + expected = -0.5; + actual = currentChunk->pressureTempCache[IX(cx,cy,cz)] * 6 - ( + currentChunk->pressureTempCache[IX(cx-1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx+1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx,cy-1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy+1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy,cz-1)] + + currentChunk->pressureTempCache[IX(cx,cy,cz+1)] + ); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to approximate pressure correctly! expected: %f actual: %f \n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("Residual: %f\n", currentChunk->projectionResidual); + printf("\n\n\n"); + } + + cx = 4; + cy = 1; + cz = 1; + //essentially this is stating that we expect neighbors of this cell to give out 0.03 more velocity than they currently will this frame + //ergo, we should subtract that from this cell in order to prevent velocity compressibility + expected = 0.0; + actual = currentChunk->pressureTempCache[IX(cx,cy,cz)] * 6 - ( + currentChunk->pressureTempCache[IX(cx-1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx+1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx,cy-1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy+1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy,cz-1)] + + currentChunk->pressureTempCache[IX(cx,cy,cz+1)] + ); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to approximate pressure correctly! expected: %f actual: %f \n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("Residual: %f\n", currentChunk->projectionResidual); + printf("\n\n\n"); + } + + return rVal; +} + +/** + * Testing pressure values + */ +int fluid_sim_pressurecell_pressure_test4(){ + printf("fluid_sim_pressurecell_pressure_test4\n"); + int rVal = 0; + Environment * env = fluid_environment_create(); + Chunk ** queue = NULL; + queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM); + int chunkCount = arrlen(queue); + + + + //setup chunk values + float deltaDensity = 0.01f; + Chunk * currentChunk = queue[0]; + currentChunk->u[CENTER_LOC][IX(1,1,1)] = 1; + currentChunk->v[CENTER_LOC][IX(1,1,1)] = 1; + currentChunk->u[CENTER_LOC][IX(2,1,1)] = 1; + + + + //actually simulate + int frameCount = 5; + for(int frame = 0; frame < frameCount; frame++){ + //simulate velocity transfering + currentChunk->u[CENTER_LOC][IX(1,1,1)] = currentChunk->u[CENTER_LOC][IX(1,1,1)] - 0.1f; + currentChunk->u[CENTER_LOC][IX(2,1,1)] = currentChunk->u[CENTER_LOC][IX(2,1,1)] - 0.1f; + currentChunk->u[CENTER_LOC][IX(3,1,1)] = currentChunk->u[CENTER_LOC][IX(3,1,1)] + 0.1f; + pressurecell_approximate_divergence(env,currentChunk); + pressurecell_approximate_pressure(env,currentChunk); + } + + //test the result + float expected, actual; + + // State of the world (velocity) + // + // ^ + // | + // | + // | + // | + // 0.5f -- 0.5f --- 0.5f --- > + + // + // cell that originall had values + // + int cx, cy, cz; + cx = 1; + cy = 1; + cz = 1; + //essentially this is stating that we expect neighbors of this cell to give out 0.02 more velocity than they should this frame + //ergo, we should subtract that from this cell in order to prevent velocity compressibility + expected = 0.02; + actual = currentChunk->pressureTempCache[IX(cx,cy,cz)]; + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to approximate pressure correctly! expected: %f actual: %f \n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("Residual: %f\n", currentChunk->projectionResidual); + printf("\n\n\n"); + } + + + + cx = 2; + cy = 1; + cz = 1; + //essentially this is stating that we expect neighbors of this cell to give out 0.02 more velocity than they should this frame + //ergo, we should subtract that from this cell in order to prevent velocity compressibility + expected = 0.02; + actual = currentChunk->pressureTempCache[IX(cx,cy,cz)] * 6 - ( + currentChunk->pressureTempCache[IX(cx-1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx+1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx,cy-1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy+1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy,cz-1)] + + currentChunk->pressureTempCache[IX(cx,cy,cz+1)] + ); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to approximate pressure correctly! expected: %f actual: %f \n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("Residual: %f\n", currentChunk->projectionResidual); + printf("\n\n\n"); + } + + cx = 3; + cy = 1; + cz = 1; + //essentially this is stating that we expect neighbors of this cell to give out -0.221 less velocity than they should this frame + //ergo, we should subtract that from this cell in order to prevent velocity compressibility + expected = -0.221; + actual = currentChunk->pressureTempCache[IX(cx,cy,cz)] * 6 - ( + currentChunk->pressureTempCache[IX(cx-1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx+1,cy,cz)] + + currentChunk->pressureTempCache[IX(cx,cy-1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy+1,cz)] + + currentChunk->pressureTempCache[IX(cx,cy,cz-1)] + + currentChunk->pressureTempCache[IX(cx,cy,cz+1)] + ); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to approximate pressure correctly! expected: %f actual: %f \n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("Residual: %f\n", currentChunk->projectionResidual); + printf("\n\n\n"); + } + + return rVal; +} + + +/** + * Testing pressure values + */ +int fluid_sim_pressurecell_pressure_test5(){ + printf("fluid_sim_pressurecell_pressure_test5\n"); + int rVal = 0; + Environment * env = fluid_environment_create(); + Chunk ** queue = NULL; + queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM); + int chunkCount = arrlen(queue); + + + + //setup chunk values + float deltaDensity = 0.01f; + Chunk * currentChunk = queue[0]; + int x, y, z; + //if everything is pulling away from 1,1,1 with maximum speed, calculate the pressure at 1,1,1 + for(x = 0; x < DIM; x++){ + for(y = 0; y < DIM; y++){ + for(z = 0; z < DIM; z++){ + currentChunk->u[CENTER_LOC][IX(x,y,z)] = 1; + currentChunk->v[CENTER_LOC][IX(x,y,z)] = 1; + currentChunk->u[CENTER_LOC][IX(x,y,z)] = 1; + if( + x == 0 || x == DIM-1 || + y == 0 || y == DIM-1 || + z == 0 || z == DIM-1 + ){ + currentChunk->u[CENTER_LOC][IX(x,y,z)] = 0; + currentChunk->v[CENTER_LOC][IX(x,y,z)] = 0; + currentChunk->u[CENTER_LOC][IX(x,y,z)] = 0; + } + } + } + } + + + + + //actually simulate + pressurecell_approximate_divergence(env,currentChunk); + pressurecell_approximate_pressure(env,currentChunk); + + //test the result + float expected, actual; + + // State of the world (velocity) + // + // ^ + // | + // | + // | + // | + // 0.5f -- 0.5f --- 0.5f --- > + + // + // cell that originall had values + // + int cx, cy, cz; + cx = 1; + cy = 1; + cz = 1; + //essentially this is stating that we expect neighbors of this cell to give out 0.02 more velocity than they should this frame + //ergo, we should subtract that from this cell in order to prevent velocity compressibility + expected = 0.02; + actual = currentChunk->pressureTempCache[IX(cx,cy,cz)]; + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal += assertEqualsFloat(expected,actual,"Failed to approximate pressure correctly! expected: %f actual: %f \n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("Residual: %f\n", currentChunk->projectionResidual); + printf("\n\n\n"); } return rVal; @@ -73,6 +588,10 @@ int fluid_sim_pressurecell_pressure_tests(int argc, char **argv){ int rVal = 0; rVal += fluid_sim_pressurecell_pressure_test1(); + rVal += fluid_sim_pressurecell_pressure_test2(); + rVal += fluid_sim_pressurecell_pressure_test3(); + // rVal += fluid_sim_pressurecell_pressure_test4(); + // rVal += fluid_sim_pressurecell_pressure_test5(); return rVal; } \ No newline at end of file diff --git a/src/test/c/fluid/sim/pressurecell/projection_tests.c b/src/test/c/fluid/sim/pressurecell/projection_tests.c new file mode 100644 index 00000000..712d88dc --- /dev/null +++ b/src/test/c/fluid/sim/pressurecell/projection_tests.c @@ -0,0 +1,443 @@ +#include + +#include "stb/stb_ds.h" + +#include "fluid/queue/boundsolver.h" +#include "fluid/queue/chunkmask.h" +#include "fluid/queue/chunk.h" +#include "fluid/env/environment.h" +#include "fluid/env/utilities.h" +#include "fluid/sim/pressurecell/pressure.h" +#include "fluid/sim/pressurecell/velocity.h" +#include "fluid/sim/pressurecell/solver_consts.h" +#include "math/ode/multigrid.h" +#include "../../../util/chunk_test_utils.h" +#include "../../../util/test.h" + +/** + * Error margin for tests + */ +#define FLUID_PRESSURE_CELL_ERROR_MARGIN 0.01f + +/** + * Number of chunks + */ +#define CHUNK_DIM 4 + +/** + * Testing projection values + */ +int fluid_sim_pressurecell_projection_test1(){ + printf("fluid_sim_pressurecell_projection_test1\n"); + int rVal = 0; + Environment * env = fluid_environment_create(); + Chunk ** queue = NULL; + queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM); + int chunkCount = arrlen(queue); + + + + //setup chunk values + float deltaDensity = 0.01f; + Chunk * currentChunk = queue[0]; + currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY; + currentChunk->u[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY; + float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING; + pressurecell_approximate_pressure(env,currentChunk); + + //actually simulate + pressurecell_project_velocity(env,currentChunk); + + //test the result + float expected, actual; + + // + // cell that originall had values + // + expected = FLUID_PRESSURECELL_MAX_VELOCITY - (currentChunk->pressureTempCache[IX(5,4,4)] - currentChunk->pressureTempCache[IX(3,4,4)]) / (gridSpacing * 2.0f); + actual = currentChunk->u[CENTER_LOC][IX(4,4,4)]; + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal++; + printf("%f \n", currentChunk->pressureTempCache[IX(4,4,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(3,4,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(5,4,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,3,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,5,4)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,4,3)]); + printf("%f \n", currentChunk->pressureTempCache[IX(4,4,5)]); + } + + return rVal; +} + +/** + * Testing projection mass conservation + */ +int fluid_sim_pressurecell_projection_test2(){ + printf("fluid_sim_pressurecell_projection_test2\n"); + int rVal = 0; + Environment * env = fluid_environment_create(); + Chunk ** queue = NULL; + queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM); + int chunkCount = arrlen(queue); + + + + //setup chunk values + float deltaDensity = 0.01f; + Chunk * currentChunk = queue[0]; + currentChunk->u[CENTER_LOC][IX(2,1,1)] = 1; + currentChunk->v[CENTER_LOC][IX(1,2,1)] = 1; + currentChunk->w[CENTER_LOC][IX(1,1,2)] = 1; + float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING; + pressurecell_approximate_divergence(env,currentChunk); + pressurecell_approximate_pressure(env,currentChunk); + + //actually simulate + pressurecell_project_velocity(env,currentChunk); + + //test the result + float expected, actual; + + // + // cell that originall had values + // + int cx, cy, cz; + float u, v, w; + float sum; + cx = 2; + cy = 1; + cz = 1; + u = currentChunk->u[CENTER_LOC][IX(cx,cy,cz)]; + v = currentChunk->v[CENTER_LOC][IX(cx,cy,cz)]; + w = currentChunk->w[CENTER_LOC][IX(cx,cy,cz)]; + sum = u + v + w; + if(u < 0 || v < 0 || w < 0 || sum <= 0.1f){ + rVal++; + printf("Projection failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("%f %f %f \n", u, v, w); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz, currentChunk->pressureTempCache[IX(cx,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx-1,cy,cz, currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx+1,cy,cz, currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy-1,cz, currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy+1,cz, currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz-1, currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz+1, currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf("pdiv x %f \n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)] - currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pdiv y %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)] - currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pdiv z %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)] - currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("\n"); + } + + // + // + // + cx = 1; + cy = 2; + cz = 1; + u = currentChunk->u[CENTER_LOC][IX(cx,cy,cz)]; + v = currentChunk->v[CENTER_LOC][IX(cx,cy,cz)]; + w = currentChunk->w[CENTER_LOC][IX(cx,cy,cz)]; + sum = u + v + w; + if(u < 0 || v < 0 || w < 0 || sum <= 0.1f){ + rVal++; + printf("Projection failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("%f %f %f \n", u, v, w); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz, currentChunk->pressureTempCache[IX(cx,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx-1,cy,cz, currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx+1,cy,cz, currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy-1,cz, currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy+1,cz, currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz-1, currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz+1, currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf("pdiv x %f \n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)] - currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pdiv y %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)] - currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pdiv z %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)] - currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("\n"); + } + + // + // + // + cx = 1; + cy = 1; + cz = 2; + u = currentChunk->u[CENTER_LOC][IX(cx,cy,cz)]; + v = currentChunk->v[CENTER_LOC][IX(cx,cy,cz)]; + w = currentChunk->w[CENTER_LOC][IX(cx,cy,cz)]; + sum = u + v + w; + if(u < 0 || v < 0 || w < 0 || sum <= 0.1f){ + rVal++; + printf("Projection failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("%f %f %f \n", u, v, w); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz, currentChunk->pressureTempCache[IX(cx,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx-1,cy,cz, currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx+1,cy,cz, currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy-1,cz, currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy+1,cz, currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz-1, currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz+1, currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf("pdiv x %f \n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)] - currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pdiv y %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)] - currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pdiv z %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)] - currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("\n"); + } + + // + // + // + cx = 1; + cy = 1; + cz = 1; + u = currentChunk->u[CENTER_LOC][IX(cx,cy,cz)]; + v = currentChunk->v[CENTER_LOC][IX(cx,cy,cz)]; + w = currentChunk->w[CENTER_LOC][IX(cx,cy,cz)]; + sum = u + v + w; + if(u < 0 || v < 0 || w < 0 || sum <= 0.1f){ + rVal++; + printf("Projection failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("%f %f %f \n", u, v, w); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz, currentChunk->pressureTempCache[IX(cx,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx-1,cy,cz, currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx+1,cy,cz, currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy-1,cz, currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy+1,cz, currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz-1, currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("pressure at[%d,%d,%d] %f \n", cx,cy,cz+1, currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf("pdiv x %f \n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)] - currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pdiv y %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)] - currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pdiv z %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)] - currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("\n"); + } + + return rVal; +} + + +/** + * Testing projection mass conservation + */ +int fluid_sim_pressurecell_projection_test3(){ + printf("fluid_sim_pressurecell_projection_test3\n"); + int rVal = 0; + Environment * env = fluid_environment_create(); + Chunk ** queue = NULL; + queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM); + int chunkCount = arrlen(queue); + + + + //setup chunk values + float deltaDensity = 0.01f; + Chunk * currentChunk = queue[0]; + currentChunk->u[CENTER_LOC][IX(1,1,1)] = 1; + currentChunk->v[CENTER_LOC][IX(1,1,1)] = 1; + currentChunk->w[CENTER_LOC][IX(1,1,1)] = 1; + currentChunk->u[CENTER_LOC][IX(2,1,1)] = 1; + currentChunk->v[CENTER_LOC][IX(1,2,1)] = 1; + currentChunk->w[CENTER_LOC][IX(1,1,2)] = 1; + float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING; + pressurecell_approximate_divergence(env,currentChunk); + pressurecell_approximate_pressure(env,currentChunk); + + //actually simulate + double maxMagnitude = pressurecell_project_velocity(env,currentChunk); + + //test the result + float expected, actual; + + // + // cell that originall had values + // + int cx, cy, cz; + float u, v, w; + float sum; + + + + cx = 1; + cy = 1; + cz = 1; + u = currentChunk->u[CENTER_LOC][IX(cx,cy,cz)]; + v = currentChunk->v[CENTER_LOC][IX(cx,cy,cz)]; + w = currentChunk->w[CENTER_LOC][IX(cx,cy,cz)]; + expected = 1.0f; + actual = sqrt(u*u + v*v + w*w); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal++; + printf("Projection failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("new force: <%f,%f,%f> \n", u, v, w); + printf("expected: %f\n", expected); + printf("actual: %f\n", actual); + printf("max magnitude: %lf\n",maxMagnitude); + printf("\n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("pdiv x %f \n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)] - currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pdiv y %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)] - currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pdiv z %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)] - currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("\n\n\n"); + } + + + + + cx = 2; + cy = 1; + cz = 1; + u = currentChunk->u[CENTER_LOC][IX(cx,cy,cz)]; + v = currentChunk->v[CENTER_LOC][IX(cx,cy,cz)]; + w = currentChunk->w[CENTER_LOC][IX(cx,cy,cz)]; + //essentially, the pressure difference between the previous point and the next point around this one + //are pulling on this velocity such that it slows down + expected = 0.596777f; + actual = sqrt(u*u + v*v + w*w); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal++; + printf("Projection failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("new force: <%f,%f,%f> \n", u, v, w); + printf("expected: %f\n", expected); + printf("actual: %f\n", actual); + printf("max magnitude: %lf\n",maxMagnitude); + printf("\n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("pdiv x %f \n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)] - currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pdiv y %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)] - currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pdiv z %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)] - currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("\n\n\n"); + } + + // + // + // + cx = 1; + cy = 2; + cz = 1; + u = currentChunk->u[CENTER_LOC][IX(cx,cy,cz)]; + v = currentChunk->v[CENTER_LOC][IX(cx,cy,cz)]; + w = currentChunk->w[CENTER_LOC][IX(cx,cy,cz)]; + //essentially, the pressure difference between the previous point and the next point around this one + //are pulling on this velocity such that it slows down + expected = 0.596777f; + actual = sqrt(u*u + v*v + w*w); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal++; + printf("Projection failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("new force: <%f,%f,%f> \n", u, v, w); + printf("expected: %f\n", expected); + printf("actual: %f\n", actual); + printf("max magnitude: %lf\n",maxMagnitude); + printf("\n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("pdiv x %f \n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)] - currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pdiv y %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)] - currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pdiv z %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)] - currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("\n\n\n"); + } + + // + // + // + cx = 1; + cy = 1; + cz = 2; + u = currentChunk->u[CENTER_LOC][IX(cx,cy,cz)]; + v = currentChunk->v[CENTER_LOC][IX(cx,cy,cz)]; + w = currentChunk->w[CENTER_LOC][IX(cx,cy,cz)]; + //essentially, the pressure difference between the previous point and the next point around this one + //are pulling on this velocity such that it slows down + expected = 0.596777f; + actual = sqrt(u*u + v*v + w*w); + if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ + rVal++; + printf("Projection failed\n"); + printf("at point (%d,%d,%d) \n", cx, cy, cz); + printf("new force: <%f,%f,%f> \n", u, v, w); + printf("expected: %f\n", expected); + printf("actual: %f\n", actual); + printf("max magnitude: %lf\n",maxMagnitude); + printf("\n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)]); + printf(" +Y | \n"); + printf(" | %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)]); + printf(" | / +Z \n"); + printf(" -X | / \n"); + printf(" B---------A--------------%f\n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)]); + printf(" / | +X \n"); + printf(" / | \n"); + printf(" -Z / | \n"); + printf(" C | -Y \n"); + printf(" %f \n", currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf(" \n"); + printf("A: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz)] ); + printf("B: %f\n", currentChunk->pressureTempCache[IX(cx-1,cy,cz)] ); + printf("C: %f\n", currentChunk->pressureTempCache[IX(cx,cy,cz-1)] ); + printf("pdiv x %f \n", currentChunk->pressureTempCache[IX(cx+1,cy,cz)] - currentChunk->pressureTempCache[IX(cx-1,cy,cz)]); + printf("pdiv y %f \n", currentChunk->pressureTempCache[IX(cx,cy+1,cz)] - currentChunk->pressureTempCache[IX(cx,cy-1,cz)]); + printf("pdiv z %f \n", currentChunk->pressureTempCache[IX(cx,cy,cz+1)] - currentChunk->pressureTempCache[IX(cx,cy,cz-1)]); + printf("\n\n\n"); + } + + return rVal; +} + +/** + * Testing projection values + */ +int fluid_sim_pressurecell_projection_tests(int argc, char **argv){ + int rVal = 0; + + rVal += fluid_sim_pressurecell_projection_test1(); + rVal += fluid_sim_pressurecell_projection_test2(); + rVal += fluid_sim_pressurecell_projection_test3(); + + return rVal; +} \ No newline at end of file