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

fluid work
Some checks failed
studiorailgun/Renderer/pipeline/head There was a failure building this commit
Details

Browse Source
This commit is contained in:
austin 2024-12-15 17:05:11 -05:00
parent 5efb1b7ea2
commit e9c79046d5
14 changed files with 1535 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.vscode/settings.json vendored
View File

@ -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"
}
}

5
src/main/c/includes/fluid/sim/pressurecell/normalization.h
View File

@ -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
*/

20
src/main/c/includes/fluid/sim/pressurecell/solver_consts.h
View File

@ -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

7
src/main/c/includes/fluid/sim/pressurecell/velocity.h
View File

@ -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

35
src/main/c/src/fluid/sim/pressurecell/density.c
View File

@ -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);
}

8
src/main/c/src/fluid/sim/pressurecell/normalization.c
View File

@ -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;
}
}
/**

194
src/main/c/src/fluid/sim/pressurecell/pressure.c
View File

@ -1,7 +1,8 @@
#include <math.h>
#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)];

42
src/main/c/src/fluid/sim/pressurecell/pressurecell.c
View File

@ -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);

205
src/main/c/src/fluid/sim/pressurecell/velocity.c
View File

@ -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");
// }
}
}
}
}
//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)];
}
}
}
}

6
src/test/c/fluid/sim/pressurecell/add_source_tests.c
View File

@ -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");
}

2
src/test/c/fluid/sim/pressurecell/advection_tests.c
View File

@ -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");
}

124
src/test/c/fluid/sim/pressurecell/divergence_tests.c
View File

@ -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;
}

541
src/test/c/fluid/sim/pressurecell/pressure_tests.c
View File

@ -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;
}

443
src/test/c/fluid/sim/pressurecell/projection_tests.c Normal file
View File

@ -0,0 +1,443 @@
#include <math.h>
#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;
}