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

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

Browse Source
This commit is contained in:
austin 2024-12-15 00:30:19 -05:00
parent eb845a3891
commit a26ed3a001
19 changed files with 665 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/main/c/includes/fluid/queue/chunk.h
View File

@ -36,6 +36,27 @@
*/ */
#define CHUNK_SPACING 16 #define CHUNK_SPACING 16
/**
* The maximum level of the interest tree
*/
#define CHUNK_MAX_INTEREST_LEVEL 4
/**
* The dimensions of the interest tree at each level
*/
static char INTEREST_MODIFIER_DIMS[CHUNK_MAX_INTEREST_LEVEL+1] = {
16,
8,
4,
2,
1,
};
/**
* Gets the interest tree at a position and level
*/
#define INTEREST(tree,level,x,y,z) tree[level][x/INTEREST_MODIFIER_DIMS[level]*INTEREST_MODIFIER_DIMS[level]*INTEREST_MODIFIER_DIMS[level]+(y/INTEREST_MODIFIER_DIMS[level])*INTEREST_MODIFIER_DIMS[level]+(z/INTEREST_MODIFIER_DIMS[level])]
/** /**
* A chunk * A chunk
*/ */
@ -120,6 +141,11 @@ typedef struct {
*/ */
int simLOD; int simLOD;
/**
* Octree that stores whether a location is of interest or not
*/
char ** interestTree;
/** /**
* The convergence of this chunk * The convergence of this chunk
*/ */

25
src/main/c/includes/fluid/sim/pressurecell/bounds.h
View File

@ -6,9 +6,34 @@
#include "fluid/queue/chunk.h" #include "fluid/queue/chunk.h"
#include "fluid/queue/chunkmask.h" #include "fluid/queue/chunkmask.h"
/**
* Used for signaling the bounds setting method to not use adjacent cells when evaluating borders
*/
#define FLUID_PRESSURECELL_BOUND_NO_DIR 0
/**
* Used for signaling the bounds setting method to use adjacent cells when evaluating x axis borders
*/
#define FLUID_PRESSURECELL_DIRECTION_U 1
/**
* Used for signaling the bounds setting method to use adjacent cells when evaluating y axis borders
*/
#define FLUID_PRESSURECELL_DIRECTION_V 2
/**
* Used for signaling the bounds setting method to use adjacent cells when evaluating z axis borders
*/
#define FLUID_PRESSURECELL_DIRECTION_W 3
/** /**
* Updates the bounds of the chunk based on its neighbors * Updates the bounds of the chunk based on its neighbors
*/ */
LIBRARY_API void pressurecell_update_bounds(Environment * environment, Chunk * chunk); LIBRARY_API void pressurecell_update_bounds(Environment * environment, Chunk * chunk);
/**
* Updates the interest tree for this chunk
*/
LIBRARY_API void pressurecell_update_interest(Environment * environment, Chunk * chunk);
#endif #endif

20
src/main/c/includes/fluid/sim/pressurecell/normalization.h Normal file
View File

@ -0,0 +1,20 @@
#ifndef FLUID_PRESSURECELL_NORMALIZATION_H
#define FLUID_PRESSURECELL_NORMALIZATION_H
#include "public.h"
#include "fluid/env/environment.h"
#include "fluid/queue/chunk.h"
#include "fluid/queue/chunkmask.h"
/**
* Calculates the ratio to normalize the chunk by
*/
LIBRARY_API void fluid_pressurecell_calculate_normalization_ratio(Environment * env, Chunk * chunk);
/**
* Normalizes the chunk
*/
LIBRARY_API void fluid_pressurecell_normalize_chunk(Environment * env, Chunk * chunk);
#endif

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

@ -7,6 +7,21 @@
*/ */
#define FLUID_PRESSURECELL_SIM_STEP 0.01f #define FLUID_PRESSURECELL_SIM_STEP 0.01f
/**
* Multiplier applied to advection to allow more mass/velocity to transfer per frame
*/
#define FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER 1.0f
/**
* Spacing of cells
*/
#define FLUID_PRESSURECELL_SPACING 1.0f
/**
* Multiplier applied to pressure calculations to encourage advection
*/
#define FLUID_PRESSURECELL_PRESSURE_MULTIPLIER 1.0f
/** /**
* Maximum allowed velocity of the pressurecell simulator * Maximum allowed velocity of the pressurecell simulator
*/ */
@ -15,17 +30,27 @@
/** /**
* Diffusion constant * Diffusion constant
*/ */
#define FLUID_PRESSURECELL_DIFFUSION_CONSTANT 0.0001f #define FLUID_PRESSURECELL_DIFFUSION_CONSTANT 0.01f
/** /**
* Viscosity constant * Viscosity constant
*/ */
#define FLUID_PRESSURECELL_VISCOSITY_CONSTANT 0.0001f #define FLUID_PRESSURECELL_VISCOSITY_CONSTANT 0.01f
/** /**
* Amount that density contributes to the pressure * Amount that density contributes to the pressure
*/ */
#define FLUID_PRESSURECELL_DENSITY_CONST 0.001f #define FLUID_PRESSURECELL_DENSITY_CONST 0.001f
/**
* Amount of the residual to add to the pressure field each frame
*/
#define FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER 0.1f
/**
* Constant applied to divergence to get new pressure
*/
#define FLUID_PRESSURECELL_DIV_PRESSURE_CONST 5.0f
#endif #endif

4
src/main/c/src/fluid/dispatch/dispatcher.c
View File

@ -28,6 +28,10 @@ LIBRARY_API void fluid_dispatch(int numReadIn, Chunk ** chunkViewC, Environment
if(countCurrent > 0){ if(countCurrent > 0){
arrdeln(environment->queue.grid2Queue,0,countCurrent); arrdeln(environment->queue.grid2Queue,0,countCurrent);
} }
countCurrent = arrlen(environment->queue.pressurecellQueue);
if(countCurrent > 0){
arrdeln(environment->queue.pressurecellQueue,0,countCurrent);
}
if(override == FLUID_DISPATCHER_OVERRIDE_CELLULAR){ if(override == FLUID_DISPATCHER_OVERRIDE_CELLULAR){

17
src/main/c/src/fluid/queue/chunk.c
View File

@ -13,5 +13,22 @@ LIBRARY_API Chunk * chunk_create(){
rVal->vTempCache = (float *)calloc(1,DIM*DIM*DIM*sizeof(float)); rVal->vTempCache = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
rVal->wTempCache = (float *)calloc(1,DIM*DIM*DIM*sizeof(float)); rVal->wTempCache = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
rVal->pressureTempCache = (float *)calloc(1,DIM*DIM*DIM*sizeof(float)); rVal->pressureTempCache = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
/**
* Should store 5 tables:
* 1x1x1
* 2x2x2
* 4x4x4
* 8x8x8
* 16x16x16
*/
rVal->interestTree = (char **)calloc(1,sizeof(char *) * 6);
//allocating extra data for ghost cells even though we will not evaluate the ghost cells
rVal->interestTree[0] = (char *)calloc(3*3*3,sizeof(char));
rVal->interestTree[1] = (char *)calloc(4*4*4,sizeof(char));
rVal->interestTree[2] = (char *)calloc(6*6*6,sizeof(char));
rVal->interestTree[3] = (char *)calloc(10*10*10,sizeof(char));
rVal->interestTree[4] = (char *)calloc(18*18*18,sizeof(char));
return rVal; return rVal;
} }

140
src/main/c/src/fluid/sim/pressurecell/bounds.c
View File

@ -1,9 +1,147 @@
#include "fluid/sim/pressurecell/bounds.h" #include "fluid/sim/pressurecell/bounds.h"
/**
* Sets the bounds reflecting off hard borders and otherwise assuming continuity
*/
void fluid_pressurecell_set_bounds_legacy(
Environment * environment,
int vector_dir,
float * target
){
//set the boundary planes
for(int x = 1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
//x-direction boundary planes
if(vector_dir == FLUID_PRESSURECELL_DIRECTION_U){
target[IX(0,x,y)] = -target[IX(1,x,y)];
target[IX(DIM-1,x,y)] = -target[IX(DIM-2,x,y)];
} else {
target[IX(0,x,y)] = target[IX(1,x,y)];
target[IX(DIM-1,x,y)] = target[IX(DIM-2,x,y)];
}
//y-direction boundary planes
if(vector_dir == FLUID_PRESSURECELL_DIRECTION_V){
target[IX(x,0,y)] = -target[IX(x,1,y)];
target[IX(x,DIM-1,y)] = -target[IX(x,DIM-2,y)];
} else {
target[IX(x,0,y)] = target[IX(x,1,y)];
target[IX(x,DIM-1,y)] = target[IX(x,DIM-2,y)];
}
//z-direction boundary planes
if(vector_dir == FLUID_PRESSURECELL_DIRECTION_W){
target[IX(x,y,0)] = -target[IX(x,y,1)];
target[IX(x,y,DIM-1)] = -target[IX(x,y,DIM-2)];
} else {
target[IX(x,y,0)] = target[IX(x,y,1)];
target[IX(x,y,DIM-1)] = target[IX(x,y,DIM-2)];
}
}
}
//sets the edges of the chunk
//this should logically follow from how we're treating the boundary planes
for(int x = 1; x < DIM-1; x++){
target[IX(x,0,0)] = (float)(0.5f * (target[IX(x,1,0)] + target[IX(x,0,1)]));
target[IX(x,DIM-1,0)] = (float)(0.5f * (target[IX(x,DIM-2,0)] + target[IX(x,DIM-1,1)]));
target[IX(x,0,DIM-1)] = (float)(0.5f * (target[IX(x,1,DIM-1)] + target[IX(x,0,DIM-2)]));
target[IX(x,DIM-1,DIM-1)] = (float)(0.5f * (target[IX(x,DIM-2,DIM-1)] + target[IX(x,DIM-1,DIM-2)]));
target[IX(0,x,0)] = (float)(0.5f * (target[IX(1,x,0)] + target[IX(0,x,1)]));
target[IX(DIM-1,x,0)] = (float)(0.5f * (target[IX(DIM-2,x,0)] + target[IX(DIM-1,x,1)]));
target[IX(0,x,DIM-1)] = (float)(0.5f * (target[IX(1,x,DIM-1)] + target[IX(0,x,DIM-2)]));
target[IX(DIM-1,x,DIM-1)] = (float)(0.5f * (target[IX(DIM-2,x,DIM-1)] + target[IX(DIM-1,x,DIM-2)]));
target[IX(0,0,x)] = (float)(0.5f * (target[IX(1,0,x)] + target[IX(0,1,x)]));
target[IX(DIM-1,0,x)] = (float)(0.5f * (target[IX(DIM-2,0,x)] + target[IX(DIM-1,1,x)]));
target[IX(0,DIM-1,x)] = (float)(0.5f * (target[IX(1,DIM-1,x)] + target[IX(0,DIM-2,x)]));
target[IX(DIM-1,DIM-1,x)] = (float)(0.5f * (target[IX(DIM-2,DIM-1,x)] + target[IX(DIM-1,DIM-2,x)]));
}
//sets the corners of the chunk
//this should logically follow from how we're treating the boundary planes
target[IX(0,0,0)] = (float)((target[IX(1,0,0)]+target[IX(0,1,0)]+target[IX(0,0,1)])/3.0);
target[IX(DIM-1,0,0)] = (float)((target[IX(DIM-2,0,0)]+target[IX(DIM-1,1,0)]+target[IX(DIM-1,0,1)])/3.0);
target[IX(0,DIM-1,0)] = (float)((target[IX(1,DIM-1,0)]+target[IX(0,DIM-2,0)]+target[IX(0,DIM-1,1)])/3.0);
target[IX(0,0,DIM-1)] = (float)((target[IX(0,0,DIM-2)]+target[IX(1,0,DIM-1)]+target[IX(0,1,DIM-1)])/3.0);
target[IX(DIM-1,DIM-1,0)] = (float)((target[IX(DIM-2,DIM-1,0)]+target[IX(DIM-1,DIM-2,0)]+target[IX(DIM-1,DIM-1,1)])/3.0);
target[IX(0,DIM-1,DIM-1)] = (float)((target[IX(1,DIM-1,DIM-1)]+target[IX(0,DIM-2,DIM-1)]+target[IX(0,DIM-1,DIM-2)])/3.0);
target[IX(DIM-1,0,DIM-1)] = (float)((target[IX(DIM-1,0,DIM-2)]+target[IX(DIM-2,0,DIM-1)]+target[IX(DIM-1,1,DIM-1)])/3.0);
target[IX(DIM-1,DIM-1,DIM-1)] = (float)((target[IX(DIM-1,DIM-1,DIM-2)]+target[IX(DIM-1,DIM-2,DIM-1)]+target[IX(DIM-1,DIM-1,DIM-2)])/3.0);
}
/** /**
* Updates the bounds of the chunk based on its neighbors * Updates the bounds of the chunk based on its neighbors
*/ */
LIBRARY_API void pressurecell_update_bounds(Environment * environment, Chunk * chunk){ LIBRARY_API void pressurecell_update_bounds(Environment * environment, Chunk * chunk){
// fluid_pressurecell_set_bounds_legacy(environment,FLUID_PRESSURECELL_BOUND_NO_DIR,chunk->d[CENTER_LOC]);
// fluid_pressurecell_set_bounds_legacy(environment,FLUID_PRESSURECELL_BOUND_NO_DIR,chunk->d0[CENTER_LOC]);
} fluid_pressurecell_set_bounds_legacy(environment,FLUID_PRESSURECELL_DIRECTION_U,chunk->u[CENTER_LOC]);
fluid_pressurecell_set_bounds_legacy(environment,FLUID_PRESSURECELL_DIRECTION_V,chunk->v[CENTER_LOC]);
fluid_pressurecell_set_bounds_legacy(environment,FLUID_PRESSURECELL_DIRECTION_W,chunk->w[CENTER_LOC]);
// fluid_pressurecell_set_bounds_legacy(environment,FLUID_PRESSURECELL_DIRECTION_U,chunk->u0[CENTER_LOC]);
// fluid_pressurecell_set_bounds_legacy(environment,FLUID_PRESSURECELL_DIRECTION_V,chunk->v0[CENTER_LOC]);
// fluid_pressurecell_set_bounds_legacy(environment,FLUID_PRESSURECELL_DIRECTION_W,chunk->w0[CENTER_LOC]);
}
/**
* Updates the interest tree for this chunk
*/
LIBRARY_API void pressurecell_update_interest(Environment * environment, Chunk * chunk){
int level, x, y, z;
int dim;
char ** interestTree = chunk->interestTree;
float * densityArr = chunk->d[CENTER_LOC];
float * densitSrcArr = chunk->d0[CENTER_LOC];
// for(level = 0; level < CHUNK_MAX_INTEREST_LEVEL; level++){
// dim = pow(2,level);
// for(x = 0; x < dim; x++){
// for(y = 0; y < dim; y++){
// for(z = 0; z < dim; z++){
// if(
// densityArr[IX(x,y,z)] > 0 ||
// densityArr[IX(x+1,y,z)] > 0 ||
// densityArr[IX(x-1,y,z)] > 0 ||
// densityArr[IX(x,y+1,z)] > 0 ||
// densityArr[IX(x,y-1,z)] > 0 ||
// densityArr[IX(x,y,z+1)] > 0 ||
// densityArr[IX(x,y,z-1)] > 0
// ){
// }
// }
// }
// }
// }
for(x = 1; x < DIM-1; x++){
for(y = 1; y < DIM-1; y++){
for(z = 1; z < DIM-1; z++){
if(
densityArr[IX(x,y,z)] > MIN_FLUID_VALUE ||
densityArr[IX(x+1,y,z)] > MIN_FLUID_VALUE ||
densityArr[IX(x-1,y,z)] > MIN_FLUID_VALUE ||
densityArr[IX(x,y+1,z)] > MIN_FLUID_VALUE ||
densityArr[IX(x,y-1,z)] > MIN_FLUID_VALUE ||
densityArr[IX(x,y,z+1)] > MIN_FLUID_VALUE ||
densityArr[IX(x,y,z-1)] > MIN_FLUID_VALUE ||
densitSrcArr[IX(x,y,z)] > MIN_FLUID_VALUE
){
INTEREST(interestTree,0,x,y,z) = 1;
INTEREST(interestTree,1,x,y,z) = 1;
INTEREST(interestTree,2,x,y,z) = 1;
INTEREST(interestTree,3,x,y,z) = 1;
INTEREST(interestTree,4,x,y,z) = 1;
}
}
}
}
}

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

@ -1,3 +1,4 @@
#include <math.h>
#include "fluid/sim/pressurecell/density.h" #include "fluid/sim/pressurecell/density.h"
#include "fluid/sim/pressurecell/solver_consts.h" #include "fluid/sim/pressurecell/solver_consts.h"
@ -25,17 +26,18 @@ LIBRARY_API void pressurecell_diffuse_density(Environment * environment, Chunk *
int x, y, z; int x, y, z;
float * densityArr = chunk->d[CENTER_LOC]; float * densityArr = chunk->d[CENTER_LOC];
float * densityTemp = chunk->dTempCache; float * densityTemp = chunk->dTempCache;
float D = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt;
for(z = 1; z < DIM-1; z++){ for(z = 1; z < DIM-1; z++){
for(y = 1; y < DIM-1; y++){ for(y = 1; y < DIM-1; y++){
for(x = 1; x < DIM-1; x++){ for(x = 1; x < DIM-1; x++){
densityTemp[IX(x,y,z)] = densityArr[IX(x,y,z)] + densityTemp[IX(x,y,z)] = densityArr[IX(x,y,z)] +
densityArr[IX(x,y,z)] * -6 * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + densityArr[IX(x,y,z)] * -6 * D +
densityArr[IX(x-1,y,z)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + densityArr[IX(x-1,y,z)] * D +
densityArr[IX(x+1,y,z)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + densityArr[IX(x+1,y,z)] * D +
densityArr[IX(x,y-1,z)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + densityArr[IX(x,y-1,z)] * D +
densityArr[IX(x,y+1,z)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + densityArr[IX(x,y+1,z)] * D +
densityArr[IX(x,y,z-1)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + densityArr[IX(x,y,z-1)] * D +
densityArr[IX(x,y,z+1)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT densityArr[IX(x,y,z+1)] * D
; ;
} }
} }
@ -61,9 +63,9 @@ LIBRARY_API void pressurecell_advect_density(Environment * environment, Chunk *
for(z = 1; z < DIM-1; z++){ for(z = 1; z < DIM-1; z++){
for(x = 1; x < DIM-1; x++){ for(x = 1; x < DIM-1; x++){
//calculate the real (float) position we are at //calculate the real (float) position we are at
xp = x - u[IX(x,y,z)] * environment->consts.dt; xp = x - u[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER;
yp = y - v[IX(x,y,z)] * environment->consts.dt; yp = y - v[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER;
zp = z - w[IX(x,y,z)] * environment->consts.dt; zp = z - w[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER;
//clamp to border //clamp to border
x0 = xp; x0 = xp;
@ -91,7 +93,12 @@ LIBRARY_API void pressurecell_advect_density(Environment * environment, Chunk *
x1 < 0 || y1 < 0 || z1 < 0 || x1 < 0 || y1 < 0 || z1 < 0 ||
x1 > DIM-1 || y1 > DIM-1 || z1 > DIM-1 x1 > DIM-1 || y1 > DIM-1 || z1 > DIM-1
){ ){
printf("advect dens: %d %d %d %d %d %d --- %f %f %f\n", x0, y0, z0, x1, y1, z1, x, y, z); printf("advect dens:\n");
printf("%d %d %d\n", x0, y0, z0);
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\n", environment->consts.dt);
fflush(stdout); fflush(stdout);
} }
@ -109,7 +116,7 @@ LIBRARY_API void pressurecell_advect_density(Environment * environment, Chunk *
t1*u1*densityTemp[IX(x1,y1,z1)] t1*u1*densityTemp[IX(x1,y1,z1)]
); );
densityArr[IX(x,y,z)] = interpolated; densityArr[IX(x,y,z)] = fmax(MIN_FLUID_VALUE,fmin(MAX_FLUID_VALUE,interpolated));
} }
} }
} }

29
src/main/c/src/fluid/sim/pressurecell/normalization.c Normal file
View File

@ -0,0 +1,29 @@
#include "fluid/sim/pressurecell/normalization.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
*/
LIBRARY_API void fluid_pressurecell_calculate_normalization_ratio(Environment * env, Chunk * chunk){
}
/**
* Normalizes the chunk
*/
LIBRARY_API void fluid_pressurecell_normalize_chunk(Environment * env, Chunk * chunk){
}

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

@ -13,21 +13,63 @@ LIBRARY_API void pressurecell_approximate_pressure(Environment * environment, Ch
//values stored across frames //values stored across frames
float * presureCache = chunk->pressureCache[CENTER_LOC]; float * presureCache = chunk->pressureCache[CENTER_LOC];
float * divArr = chunk->divergenceCache[CENTER_LOC]; float * divArr = chunk->divergenceCache[CENTER_LOC];
// float * uArr = chunk->u[CENTER_LOC];
// float * vArr = chunk->v[CENTER_LOC];
// float * wArr = chunk->w[CENTER_LOC];
float * uArr = chunk->uTempCache;
float * vArr = chunk->vTempCache;
float * wArr = chunk->wTempCache;
//temporary caches //temporary caches
float * pressureCache = chunk->pressureCache[CENTER_LOC];
float * pressureTemp = chunk->pressureTempCache; float * pressureTemp = chunk->pressureTempCache;
float gridSpacing = 1; //consts/vars
float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING;
float du, dv, dw; float du, dv, dw;
float newPressure; float newPressure;
float dt = environment->consts.dt;
for(z = 1; z < DIM-1; z++){ for(z = 1; z < DIM-1; z++){
for(y = 1; y < DIM-1; y++){ for(y = 1; y < DIM-1; y++){
for(x = 1; x < DIM-1; x++){ 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 //compute the new pressure value
newPressure = newPressure = pressureCache[IX(x,y,z)] + residual;
(divArr[IX(x+1,y,z)] - divArr[IX(x-1,y,z)]) / (2.0f * gridSpacing) + // newPressure = -FLUID_PRESSURECELL_DIV_PRESSURE_CONST * (newDiv + outflowDiv);
(divArr[IX(x+1,y,z)] - divArr[IX(x-1,y,z)]) / (2.0f * gridSpacing) + // if(newPressure > 0){
(divArr[IX(x+1,y,z)] - divArr[IX(x-1,y,z)]) / (2.0f * gridSpacing) // printf("%f \n",newPressure);
; // }
newPressure = newPressure / (gridSpacing * 2);
pressureTemp[IX(x,y,z)] = newPressure; pressureTemp[IX(x,y,z)] = newPressure;
} }
} }
@ -47,18 +89,38 @@ LIBRARY_API void pressurecell_approximate_divergence(Environment * environment,
float * presureCache = chunk->pressureCache[CENTER_LOC]; float * presureCache = chunk->pressureCache[CENTER_LOC];
//temporary caches //temporary caches
float * pressureTemp = chunk->pressureTempCache; float * pressureTemp = chunk->pressureTempCache;
float gridSpacing = 1; float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING;
float du, dv, dw; float du, dv, dw;
float newDivergence; float newDivergence;
float dt = environment->consts.dt;
for(z = 1; z < DIM-1; z++){ for(z = 1; z < DIM-1; z++){
for(y = 1; y < DIM-1; y++){ for(y = 1; y < DIM-1; y++){
for(x = 1; x < DIM-1; x++){ for(x = 1; x < DIM-1; x++){
//divergence caused by static outflow due to diffusion
float outflowDiv = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 6 * dt;
//compute divergence //compute divergence
du = (uArr[IX(x+1,y,z)] - uArr[IX(x-1,y,z)]) / (2.0f * gridSpacing); du = (uArr[IX(x+1,y,z)] - uArr[IX(x-1,y,z)]) / (2.0f * gridSpacing);
dv = (vArr[IX(x+1,y,z)] - vArr[IX(x-1,y,z)]) / (2.0f * gridSpacing); dv = (vArr[IX(x,y+1,z)] - vArr[IX(x,y-1,z)]) / (2.0f * gridSpacing);
dw = (wArr[IX(x+1,y,z)] - wArr[IX(x-1,y,z)]) / (2.0f * gridSpacing); dw = (wArr[IX(x,y,z+1)] - wArr[IX(x,y,z-1)]) / (2.0f * gridSpacing);
// if(x == 1){
// du = 0;
// } else if(x == DIM-2){
// du = 0;
// }
// if(y == 1){
// dv = 0;
// } else if(y == DIM-2){
// dv = 0;
// }
// if(z == 1){
// dw = 0;
// } else if(z == DIM-2){
// dw = 0;
// }
newDivergence = du+dv+dw; newDivergence = du+dv+dw;
divArr[IX(x,y,z)] = newDivergence; divArr[IX(x,y,z)] = divArr[IX(x,y,z)] + newDivergence - FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER * divArr[IX(x,y,z)] + outflowDiv;
//store pressure value from this frame //store pressure value from this frame
presureCache[IX(x,y,z)] = pressureTemp[IX(x,y,z)]; presureCache[IX(x,y,z)] = pressureTemp[IX(x,y,z)];

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

@ -9,11 +9,18 @@
#include "fluid/env/utilities.h" #include "fluid/env/utilities.h"
#include "fluid/queue/chunkmask.h" #include "fluid/queue/chunkmask.h"
#include "fluid/queue/chunk.h" #include "fluid/queue/chunk.h"
#include "fluid/sim/pressurecell/bounds.h"
#include "fluid/sim/pressurecell/density.h" #include "fluid/sim/pressurecell/density.h"
#include "fluid/sim/pressurecell/pressure.h" #include "fluid/sim/pressurecell/pressure.h"
#include "fluid/sim/pressurecell/solver_consts.h" #include "fluid/sim/pressurecell/solver_consts.h"
#include "fluid/sim/pressurecell/velocity.h" #include "fluid/sim/pressurecell/velocity.h"
static inline void fluid_pressurecell_clearArr(float * d);
/** /**
* Used for storing timings * Used for storing timings
*/ */
@ -31,6 +38,16 @@ LIBRARY_API void fluid_pressurecell_simulate(
//This is the section of non-parallel code //This is the section of non-parallel code
// //
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
pressurecell_update_bounds(environment,currentChunk);
// pressurecell_update_interest(environment,currentChunk);
}
// //
//This is the section of parallel code //This is the section of parallel code
// //
@ -96,6 +113,23 @@ LIBRARY_API void fluid_pressurecell_simulate(
Chunk * currentChunk = chunks[i]; Chunk * currentChunk = chunks[i];
pressurecell_advect_density(environment,currentChunk); pressurecell_advect_density(environment,currentChunk);
} }
//
// Setup for next iteration
//
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
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->uTempCache);
fluid_pressurecell_clearArr(currentChunk->vTempCache);
fluid_pressurecell_clearArr(currentChunk->wTempCache);
}
} }
@ -103,7 +137,14 @@ LIBRARY_API void fluid_pressurecell_simulate(
/**
* Clears an array
*/
static inline void fluid_pressurecell_clearArr(float * d){
for(int j = 0; j < DIM * DIM * DIM; j++){
d[j] = 0;
}
}

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

@ -1,4 +1,6 @@
#include <math.h>
#include "fluid/queue/chunk.h"
#include "fluid/sim/pressurecell/velocity.h" #include "fluid/sim/pressurecell/velocity.h"
#include "fluid/sim/pressurecell/solver_consts.h" #include "fluid/sim/pressurecell/solver_consts.h"
@ -9,7 +11,6 @@
LIBRARY_API void pressurecell_add_gravity(Environment * environment, Chunk * chunk){ LIBRARY_API void pressurecell_add_gravity(Environment * environment, Chunk * chunk){
int x, y, z; int x, y, z;
float * dArr = chunk->d[CENTER_LOC]; float * dArr = chunk->d[CENTER_LOC];
float * vArr = chunk->v[CENTER_LOC];
float * vSourceArr = chunk->v0[CENTER_LOC]; float * vSourceArr = chunk->v0[CENTER_LOC];
for(z = 1; z < DIM-1; z++){ for(z = 1; z < DIM-1; z++){
for(y = 1; y < DIM-1; y++){ for(y = 1; y < DIM-1; y++){
@ -48,46 +49,57 @@ LIBRARY_API void pressurecell_add_velocity(Environment * environment, Chunk * ch
*/ */
LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk * chunk){ LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk * chunk){
int x, y, z; int x, y, z;
float * dArr = chunk->d[CENTER_LOC];
float * uArr = chunk->u[CENTER_LOC]; float * uArr = chunk->u[CENTER_LOC];
float * vArr = chunk->v[CENTER_LOC]; float * vArr = chunk->v[CENTER_LOC];
float * wArr = chunk->w[CENTER_LOC]; float * wArr = chunk->w[CENTER_LOC];
float * uTemp = chunk->uTempCache; float * uTemp = chunk->uTempCache;
float * vTemp = chunk->vTempCache; float * vTemp = chunk->vTempCache;
float * wTemp = chunk->wTempCache; float * wTemp = chunk->wTempCache;
float D = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt;
for(z = 1; z < DIM-1; z++){ for(z = 1; z < DIM-1; z++){
for(y = 1; y < DIM-1; y++){ for(y = 1; y < DIM-1; y++){
for(x = 1; x < DIM-1; x++){ for(x = 1; x < DIM-1; x++){
//diffuse u
uTemp[IX(x,y,z)] = uArr[IX(x,y,z)] + uTemp[IX(x,y,z)] = uArr[IX(x,y,z)] +
uArr[IX(x,y,z)] * -6 * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + uArr[IX(x,y,z)] * -6 * D +
uArr[IX(x-1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + uArr[IX(x-1,y,z)] * D +
uArr[IX(x+1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + uArr[IX(x+1,y,z)] * D +
uArr[IX(x,y-1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + uArr[IX(x,y-1,z)] * D +
uArr[IX(x,y+1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + uArr[IX(x,y+1,z)] * D +
uArr[IX(x,y,z-1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + uArr[IX(x,y,z-1)] * D +
uArr[IX(x,y,z+1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT uArr[IX(x,y,z+1)] * D
; ;
}
}
}
//diffuse v for(z = 1; z < DIM-1; z++){
for(y = 1; y < DIM-1; y++){
for(x = 1; x < DIM-1; x++){
vTemp[IX(x,y,z)] = vArr[IX(x,y,z)] + vTemp[IX(x,y,z)] = vArr[IX(x,y,z)] +
vArr[IX(x,y,z)] * -6 * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + vArr[IX(x,y,z)] * -6 * D +
vArr[IX(x-1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + vArr[IX(x-1,y,z)] * D +
vArr[IX(x+1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + vArr[IX(x+1,y,z)] * D +
vArr[IX(x,y-1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + vArr[IX(x,y-1,z)] * D +
vArr[IX(x,y+1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + vArr[IX(x,y+1,z)] * D +
vArr[IX(x,y,z-1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + vArr[IX(x,y,z-1)] * D +
vArr[IX(x,y,z+1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT vArr[IX(x,y,z+1)] * D
; ;
}
}
}
//diffuse w for(z = 1; z < DIM-1; z++){
for(y = 1; y < DIM-1; y++){
for(x = 1; x < DIM-1; x++){
wTemp[IX(x,y,z)] = wArr[IX(x,y,z)] + wTemp[IX(x,y,z)] = wArr[IX(x,y,z)] +
wArr[IX(x,y,z)] * -6 * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + wArr[IX(x,y,z)] * -6 * D +
wArr[IX(x-1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + wArr[IX(x-1,y,z)] * D +
wArr[IX(x+1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + wArr[IX(x+1,y,z)] * D +
wArr[IX(x,y-1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + wArr[IX(x,y-1,z)] * D +
wArr[IX(x,y+1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + wArr[IX(x,y+1,z)] * D +
wArr[IX(x,y,z-1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + wArr[IX(x,y,z-1)] * D +
wArr[IX(x,y,z+1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT wArr[IX(x,y,z+1)] * D
; ;
} }
} }
@ -108,15 +120,15 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk *
int x0, x1, y0, y1, z0, z1; int x0, x1, y0, y1, z0, z1;
float xp, yp, zp; float xp, yp, zp;
float s0, s1, t0, t1, u0, u1; float s0, s1, t0, t1, u0, u1;
float interpolated; float interpolatedU, interpolatedV, interpolatedW;
float magnitude;
for(y = 1; y < DIM-1; y++){ 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(z = 1; z < DIM-1; z++){
for(x = 1; x < DIM-1; x++){ for(x = 1; x < DIM-1; x++){
//calculate the real (float) position we are at //calculate the real (float) position we are at
xp = x - uTemp[IX(x,y,z)] * environment->consts.dt; xp = x - uTemp[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER;
yp = y - vTemp[IX(x,y,z)] * environment->consts.dt; yp = y - vTemp[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER;
zp = z - wTemp[IX(x,y,z)] * environment->consts.dt; zp = z - wTemp[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER;
//clamp to border //clamp to border
x0 = xp; x0 = xp;
@ -138,21 +150,7 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk *
u0 = 1-u1; u0 = 1-u1;
if( interpolatedU =
x0 < 0 || y0 < 0 || z0 < 0 ||
x0 > DIM-1 || y0 > DIM-1 || z0 > DIM-1 ||
x1 < 0 || y1 < 0 || z1 < 0 ||
x1 > DIM-1 || y1 > DIM-1 || z1 > DIM-1
){
printf("advect vel: out of bounds \n");
printf("%d %d %d\n", x0, y0, z0);
printf("%d %d %d\n ", x1, y1, z1);
printf("%f %f %f\n", x, y, z);
printf("\n");
fflush(stdout);
}
interpolated =
s0*( s0*(
t0*u0*uTemp[IX(x0,y0,z0)]+ t0*u0*uTemp[IX(x0,y0,z0)]+
t1*u0*uTemp[IX(x0,y1,z0)]+ t1*u0*uTemp[IX(x0,y1,z0)]+
@ -165,9 +163,7 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk *
t0*u1*uTemp[IX(x1,y0,z1)]+ t0*u1*uTemp[IX(x1,y0,z1)]+
t1*u1*uTemp[IX(x1,y1,z1)] t1*u1*uTemp[IX(x1,y1,z1)]
); );
uArr[IX(x,y,z)] = interpolated; interpolatedV =
interpolated =
s0*( s0*(
t0*u0*vTemp[IX(x0,y0,z0)]+ t0*u0*vTemp[IX(x0,y0,z0)]+
t1*u0*vTemp[IX(x0,y1,z0)]+ t1*u0*vTemp[IX(x0,y1,z0)]+
@ -180,9 +176,7 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk *
t0*u1*vTemp[IX(x1,y0,z1)]+ t0*u1*vTemp[IX(x1,y0,z1)]+
t1*u1*vTemp[IX(x1,y1,z1)] t1*u1*vTemp[IX(x1,y1,z1)]
); );
vArr[IX(x,y,z)] = interpolated; interpolatedW =
interpolated =
s0*( s0*(
t0*u0*wTemp[IX(x0,y0,z0)]+ t0*u0*wTemp[IX(x0,y0,z0)]+
t1*u0*wTemp[IX(x0,y1,z0)]+ t1*u0*wTemp[IX(x0,y1,z0)]+
@ -195,7 +189,41 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk *
t0*u1*wTemp[IX(x1,y0,z1)]+ t0*u1*wTemp[IX(x1,y0,z1)]+
t1*u1*wTemp[IX(x1,y1,z1)] t1*u1*wTemp[IX(x1,y1,z1)]
); );
wArr[IX(x,y,z)] = interpolated;
if(
x0 < 0 || y0 < 0 || z0 < 0 ||
x0 > DIM-1 || y0 > DIM-1 || z0 > DIM-1 ||
x1 < 0 || y1 < 0 || z1 < 0 ||
x1 > DIM-1 || y1 > DIM-1 || z1 > DIM-1
){
printf("advect vel: out of bounds \n");
printf("%d %d %d\n", x, y, z);
printf("%d %d %d\n", x0, y0, z0);
printf("%d %d %d\n", x1, y1, z1);
printf("percentages: \n");
printf("%f %f %f\n", s0, t0, u0);
printf("%f %f %f\n", s1, t1, u1);
printf("%f %f %f\n", xp, yp, zp);
printf("interpolated:\n");
printf("%f %f %f\n", interpolatedU, interpolatedV, interpolatedW);
printf("%f %f %f\n", uTemp[IX(x,y,z)], vTemp[IX(x,y,z)], wTemp[IX(x,y,z)]);
printf("%f\n", environment->consts.dt);
printf("\n");
fflush(stdout);
}
// magnitude = sqrt(interpolatedU * interpolatedU + interpolatedV * interpolatedV + interpolatedW * interpolatedW);
// if(magnitude > 1){
// interpolatedU = interpolatedU / magnitude;
// interpolatedV = interpolatedV / magnitude;
// interpolatedW = interpolatedW / magnitude;
// }
uArr[IX(x,y,z)] = interpolatedU;
vArr[IX(x,y,z)] = interpolatedV;
wArr[IX(x,y,z)] = interpolatedW;
} }
} }
} }
@ -206,7 +234,7 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk *
*/ */
LIBRARY_API void pressurecell_interpolate_velocity(Environment * environment, Chunk * chunk){ LIBRARY_API void pressurecell_interpolate_velocity(Environment * environment, Chunk * chunk){
int x, y, z; int x, y, z;
float * presureCache = chunk->pressureCache[CENTER_LOC]; float * presureCache = chunk->pressureTempCache;
float * uArr = chunk->u[CENTER_LOC]; float * uArr = chunk->u[CENTER_LOC];
float * vArr = chunk->v[CENTER_LOC]; float * vArr = chunk->v[CENTER_LOC];
float * wArr = chunk->w[CENTER_LOC]; float * wArr = chunk->w[CENTER_LOC];
@ -214,30 +242,49 @@ LIBRARY_API void pressurecell_interpolate_velocity(Environment * environment, Ch
float * vTemp = chunk->vTempCache; float * vTemp = chunk->vTempCache;
float * wTemp = chunk->wTempCache; float * wTemp = chunk->wTempCache;
//temporary caches //temporary caches
float * pressureTemp = chunk->pressureTempCache; float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING;
float gridSpacing = 1;
float du, dv, dw; float du, dv, dw;
float pressureDivergence; float pressureDivergence;
for(z = 1; z < DIM-1; z++){ float magnitude;
for(y = 1; y < DIM-1; y++){ for(y = 1; y < DIM-1; y++){
for(z = 1; z < DIM-1; z++){
for(x = 1; x < DIM-1; x++){ for(x = 1; x < DIM-1; x++){
//compute the new pressure value
pressureDivergence =
(presureCache[IX(x+1,y,z)] - presureCache[IX(x-1,y,z)]) / (2.0f * gridSpacing) +
(presureCache[IX(x+1,y,z)] - presureCache[IX(x-1,y,z)]) / (2.0f * gridSpacing) +
(presureCache[IX(x+1,y,z)] - presureCache[IX(x-1,y,z)]) / (2.0f * gridSpacing)
;
pressureDivergence = pressureDivergence / (gridSpacing * 2);
//project the pressure gradient onto the velocity field //project the pressure gradient onto the velocity field
uTemp[IX(x,y,z)] = uTemp[IX(x,y,z)] - pressureDivergence; 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)] - pressureDivergence; 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)] - pressureDivergence; 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);
// 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)]);
//map the new velocity field onto the old one //map the new velocity field onto the old one
uArr[IX(x,y,z)] = uTemp[IX(x,y,z)]; // if(magnitude > 1.0f){
vArr[IX(x,y,z)] = vTemp[IX(x,y,z)]; // uArr[IX(x,y,z)] = uTemp[IX(x,y,z)] / magnitude;
wArr[IX(x,y,z)] = wTemp[IX(x,y,z)]; // vArr[IX(x,y,z)] = vTemp[IX(x,y,z)] / magnitude;
// wArr[IX(x,y,z)] = wTemp[IX(x,y,z)] / magnitude;
// } else if(magnitude > 0.0f){
// 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)];
// } else {
// uArr[IX(x,y,z)] = 0;
// vArr[IX(x,y,z)] = 0;
// wArr[IX(x,y,z)] = 0;
// }
// 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
// ){
// 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("\n");
// }
} }
} }
} }

2
src/main/java/electrosphere/client/script/ScriptClientVoxelUtils.java
View File

@ -66,7 +66,7 @@ public class ScriptClientVoxelUtils {
Vector3d centerPos = new Vector3d(CameraEntityUtils.getCameraCenter(camera)); Vector3d centerPos = new Vector3d(CameraEntityUtils.getCameraCenter(camera));
Vector3d cursorPos = collisionEngine.rayCastPosition(new Vector3d(centerPos), new Vector3d(eyePos).mul(-1.0), CollisionEngine.DEFAULT_INTERACT_DISTANCE); Vector3d cursorPos = collisionEngine.rayCastPosition(new Vector3d(centerPos), new Vector3d(eyePos).mul(-1.0), CollisionEngine.DEFAULT_INTERACT_DISTANCE);
if(cursorPos == null){ if(cursorPos == null){
cursorPos = new Vector3d(centerPos).add(new Vector3d(eyePos).mul(-CollisionEngine.DEFAULT_INTERACT_DISTANCE)); cursorPos = new Vector3d(centerPos).add(new Vector3d(eyePos).normalize().mul(-CollisionEngine.DEFAULT_INTERACT_DISTANCE));
} }
cursorPos = cursorPos.add(cursorVerticalOffset); cursorPos = cursorPos.add(cursorVerticalOffset);
Vector3i worldPos = new Vector3i( Vector3i worldPos = new Vector3i(

17
src/main/java/electrosphere/server/fluid/manager/ServerFluidChunk.java
View File

@ -106,6 +106,11 @@ public class ServerFluidChunk {
*/ */
FloatBuffer bounds; FloatBuffer bounds;
/**
* The pressure cache
*/
FloatBuffer pressureCache;
/** /**
* The array of all adjacent weight buffers for the fluid sim * The array of all adjacent weight buffers for the fluid sim
*/ */
@ -244,6 +249,7 @@ public class ServerFluidChunk {
this.velocityY = this.bVelocityY[CENTER_BUFF].asFloatBuffer(); this.velocityY = this.bVelocityY[CENTER_BUFF].asFloatBuffer();
this.velocityZ = this.bVelocityZ[CENTER_BUFF].asFloatBuffer(); this.velocityZ = this.bVelocityZ[CENTER_BUFF].asFloatBuffer();
this.bounds = this.bBounds[CENTER_BUFF].asFloatBuffer(); this.bounds = this.bBounds[CENTER_BUFF].asFloatBuffer();
this.pressureCache = this.bPressureCache[CENTER_BUFF].asFloatBuffer();
} }
/** /**
@ -488,6 +494,17 @@ public class ServerFluidChunk {
this.bounds.put(this.IX(x,y,z),bound); this.bounds.put(this.IX(x,y,z),bound);
} }
/**
* Sets a pressure
* @param x The x coordinate
* @param y The y coordinate
* @param z The z coordinate
* @param pressure The pressure
*/
public void setPressure(int x, int y, int z, float pressure){
pressureCache.put(this.IX(x,y,z),pressure);
}
/** /**
* Gets the inddex into the buffer * Gets the inddex into the buffer
* @param x The x position * @param x The x position

1
src/main/java/electrosphere/server/fluid/manager/ServerFluidManager.java
View File

@ -275,6 +275,7 @@ public class ServerFluidManager {
ServerFluidChunk fluidChunk = this.getChunk(worldPos.x, worldPos.y, worldPos.z); ServerFluidChunk fluidChunk = this.getChunk(worldPos.x, worldPos.y, worldPos.z);
fluidChunk.setAsleep(false); fluidChunk.setAsleep(false);
fluidChunk.setWeight(voxelPos.x, voxelPos.y, voxelPos.z, weight); fluidChunk.setWeight(voxelPos.x, voxelPos.y, voxelPos.z, weight);
fluidChunk.setPressure(voxelPos.x, voxelPos.y, voxelPos.z, weight);
lock.unlock(); lock.unlock();
} }

2
src/main/java/electrosphere/server/fluid/simulator/FluidAcceleratedSimulator.java
View File

@ -29,7 +29,7 @@ public class FluidAcceleratedSimulator implements ServerFluidSimulator {
/** /**
* The gravity constant * The gravity constant
*/ */
public static final float GRAVITY_CONST = -1000f; public static final float GRAVITY_CONST = -100f;
/** /**
* Load fluid sim library * Load fluid sim library

30
src/test/c/fluid/sim/pressurecell/diffuse_tests.c
View File

@ -50,7 +50,7 @@ int fluid_sim_pressurecell_diffuse_test1(){
// //
// cell that originall had values // cell that originall had values
// //
expected = MAX_FLUID_VALUE - FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 6 * MAX_FLUID_VALUE; expected = MAX_FLUID_VALUE - FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 6 * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->dTempCache[IX(4,4,4)]; actual = currentChunk->dTempCache[IX(4,4,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,4)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,4)! expected: %f actual: %f \n");
@ -60,37 +60,37 @@ int fluid_sim_pressurecell_diffuse_test1(){
// //
// neighbors // neighbors
// //
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->dTempCache[IX(3,4,4)]; actual = currentChunk->dTempCache[IX(3,4,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (3,4,4)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (3,4,4)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->dTempCache[IX(5,4,4)]; actual = currentChunk->dTempCache[IX(5,4,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN * env->consts.dt){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (5,4,4)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (5,4,4)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->dTempCache[IX(4,3,4)]; actual = currentChunk->dTempCache[IX(4,3,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,3,4)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,3,4)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->dTempCache[IX(4,5,4)]; actual = currentChunk->dTempCache[IX(4,5,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,5,4)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,5,4)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->dTempCache[IX(4,4,3)]; actual = currentChunk->dTempCache[IX(4,4,3)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,3)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,3)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->dTempCache[IX(4,4,5)]; actual = currentChunk->dTempCache[IX(4,4,5)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,5)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,5)! expected: %f actual: %f \n");
@ -126,7 +126,7 @@ int fluid_sim_pressurecell_diffuse_test2(){
// //
// cell that originall had values // cell that originall had values
// //
expected = MAX_FLUID_VALUE - FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 6 * MAX_FLUID_VALUE; expected = MAX_FLUID_VALUE - FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 6 * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->uTempCache[IX(4,4,4)]; actual = currentChunk->uTempCache[IX(4,4,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ 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"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,4)! expected: %f actual: %f \n");
@ -136,37 +136,37 @@ int fluid_sim_pressurecell_diffuse_test2(){
// //
// neighbors // neighbors
// //
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->uTempCache[IX(3,4,4)]; actual = currentChunk->uTempCache[IX(3,4,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (3,4,4)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (3,4,4)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->uTempCache[IX(5,4,4)]; actual = currentChunk->uTempCache[IX(5,4,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (5,4,4)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (5,4,4)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->uTempCache[IX(4,3,4)]; actual = currentChunk->uTempCache[IX(4,3,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,3,4)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,3,4)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->uTempCache[IX(4,5,4)]; actual = currentChunk->uTempCache[IX(4,5,4)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,5,4)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,5,4)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->uTempCache[IX(4,4,3)]; actual = currentChunk->uTempCache[IX(4,4,3)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,3)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,3)! expected: %f actual: %f \n");
} }
expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE; expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
actual = currentChunk->uTempCache[IX(4,4,5)]; actual = currentChunk->uTempCache[IX(4,4,5)];
if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,5)! expected: %f actual: %f \n"); rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,5)! expected: %f actual: %f \n");

78
src/test/c/fluid/sim/pressurecell/divergence_tests.c Normal file
View File

@ -0,0 +1,78 @@
#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/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.00001f
/**
* Number of chunks
*/
#define CHUNK_DIM 4
/**
* Testing pressure values
*/
int fluid_sim_pressurecell_divergence_test1(){
printf("fluid_sim_pressurecell_divergence_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)] = 0;
currentChunk->u[CENTER_LOC][IX(3,4,4)] = 1;
currentChunk->u[CENTER_LOC][IX(5,4,4)] = -1;
currentChunk->v[CENTER_LOC][IX(4,3,4)] = 1;
currentChunk->v[CENTER_LOC][IX(4,5,4)] = -1;
currentChunk->w[CENTER_LOC][IX(4,4,3)] = 1;
currentChunk->w[CENTER_LOC][IX(4,4,5)] = -1;
//actually simulate
pressurecell_approximate_divergence(env,currentChunk);
//test the result
float expected, actual;
//
// cell that originall had values
//
expected = -3 * env->consts.dt;
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");
}
return rVal;
}
/**
* Testing pressure values
*/
int fluid_sim_pressurecell_divergence_tests(int argc, char **argv){
int rVal = 0;
rVal += fluid_sim_pressurecell_divergence_test1();
return rVal;
}

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

@ -39,11 +39,14 @@ int fluid_sim_pressurecell_pressure_test1(){
//setup chunk values //setup chunk values
float deltaDensity = 0.01f; float deltaDensity = 0.01f;
Chunk * currentChunk = queue[0]; Chunk * currentChunk = queue[0];
currentChunk->d[CENTER_LOC][IX(4,4,4)] = MAX_FLUID_VALUE - deltaDensity; currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)] = 3;
currentChunk->d0[CENTER_LOC][IX(4,4,4)] = deltaDensity;
currentChunk->u0[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY; currentChunk->pressureCache[CENTER_LOC][IX(3,4,4)] = 1;
currentChunk->v0[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY; currentChunk->pressureCache[CENTER_LOC][IX(5,4,4)] = -1;
currentChunk->w0[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY; 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;
//actually simulate //actually simulate
pressurecell_approximate_pressure(env,currentChunk); pressurecell_approximate_pressure(env,currentChunk);
@ -54,11 +57,11 @@ int fluid_sim_pressurecell_pressure_test1(){
// //
// cell that originall had values // cell that originall had values
// //
// expected = MAX_FLUID_VALUE - FLUID_PRESSURECELL_MAX_VELOCITY * env->consts.dt * MAX_FLUID_VALUE; expected = -FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER * currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)];
// actual = currentChunk->pressureCache[CENTER_LOC][IX(4,4,4)]; actual = currentChunk->pressureTempCache[IX(4,4,4)];
// if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){ 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"); rVal += assertEqualsFloat(expected,actual,"Failed to advect density correctly (4,4,4)! expected: %f actual: %f \n");
// } }
return rVal; return rVal;
} }