From a26ed3a0015e7fee2ebe9592a0c2aaa2f598e287 Mon Sep 17 00:00:00 2001 From: austin Date: Sun, 15 Dec 2024 00:30:19 -0500 Subject: [PATCH] work on pressurecell solver --- src/main/c/includes/fluid/queue/chunk.h | 26 +++ .../includes/fluid/sim/pressurecell/bounds.h | 25 +++ .../fluid/sim/pressurecell/normalization.h | 20 ++ .../fluid/sim/pressurecell/solver_consts.h | 29 ++- src/main/c/src/fluid/dispatch/dispatcher.c | 4 + src/main/c/src/fluid/queue/chunk.c | 17 ++ .../c/src/fluid/sim/pressurecell/bounds.c | 140 ++++++++++++- .../c/src/fluid/sim/pressurecell/density.c | 31 +-- .../fluid/sim/pressurecell/normalization.c | 29 +++ .../c/src/fluid/sim/pressurecell/pressure.c | 84 +++++++- .../src/fluid/sim/pressurecell/pressurecell.c | 43 +++- .../c/src/fluid/sim/pressurecell/velocity.c | 189 +++++++++++------- .../client/script/ScriptClientVoxelUtils.java | 2 +- .../fluid/manager/ServerFluidChunk.java | 17 ++ .../fluid/manager/ServerFluidManager.java | 1 + .../simulator/FluidAcceleratedSimulator.java | 2 +- .../c/fluid/sim/pressurecell/diffuse_tests.c | 30 +-- .../fluid/sim/pressurecell/divergence_tests.c | 78 ++++++++ .../c/fluid/sim/pressurecell/pressure_tests.c | 23 ++- 19 files changed, 665 insertions(+), 125 deletions(-) create mode 100644 src/main/c/includes/fluid/sim/pressurecell/normalization.h create mode 100644 src/main/c/src/fluid/sim/pressurecell/normalization.c create mode 100644 src/test/c/fluid/sim/pressurecell/divergence_tests.c diff --git a/src/main/c/includes/fluid/queue/chunk.h b/src/main/c/includes/fluid/queue/chunk.h index fead4f85..99a7bde4 100644 --- a/src/main/c/includes/fluid/queue/chunk.h +++ b/src/main/c/includes/fluid/queue/chunk.h @@ -36,6 +36,27 @@ */ #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 */ @@ -120,6 +141,11 @@ typedef struct { */ int simLOD; + /** + * Octree that stores whether a location is of interest or not + */ + char ** interestTree; + /** * The convergence of this chunk */ diff --git a/src/main/c/includes/fluid/sim/pressurecell/bounds.h b/src/main/c/includes/fluid/sim/pressurecell/bounds.h index 5aa83572..5e16e3c6 100644 --- a/src/main/c/includes/fluid/sim/pressurecell/bounds.h +++ b/src/main/c/includes/fluid/sim/pressurecell/bounds.h @@ -6,9 +6,34 @@ #include "fluid/queue/chunk.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 */ 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 \ No newline at end of file diff --git a/src/main/c/includes/fluid/sim/pressurecell/normalization.h b/src/main/c/includes/fluid/sim/pressurecell/normalization.h new file mode 100644 index 00000000..8c7fc223 --- /dev/null +++ b/src/main/c/includes/fluid/sim/pressurecell/normalization.h @@ -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 \ No newline at end of file diff --git a/src/main/c/includes/fluid/sim/pressurecell/solver_consts.h b/src/main/c/includes/fluid/sim/pressurecell/solver_consts.h index 39426586..37aa1bb6 100644 --- a/src/main/c/includes/fluid/sim/pressurecell/solver_consts.h +++ b/src/main/c/includes/fluid/sim/pressurecell/solver_consts.h @@ -7,6 +7,21 @@ */ #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 */ @@ -15,17 +30,27 @@ /** * Diffusion constant */ -#define FLUID_PRESSURECELL_DIFFUSION_CONSTANT 0.0001f +#define FLUID_PRESSURECELL_DIFFUSION_CONSTANT 0.01f /** * Viscosity constant */ -#define FLUID_PRESSURECELL_VISCOSITY_CONSTANT 0.0001f +#define FLUID_PRESSURECELL_VISCOSITY_CONSTANT 0.01f /** * Amount that density contributes to the pressure */ #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 \ No newline at end of file diff --git a/src/main/c/src/fluid/dispatch/dispatcher.c b/src/main/c/src/fluid/dispatch/dispatcher.c index d746d8bf..9eecdd9b 100644 --- a/src/main/c/src/fluid/dispatch/dispatcher.c +++ b/src/main/c/src/fluid/dispatch/dispatcher.c @@ -28,6 +28,10 @@ LIBRARY_API void fluid_dispatch(int numReadIn, Chunk ** chunkViewC, Environment if(countCurrent > 0){ 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){ diff --git a/src/main/c/src/fluid/queue/chunk.c b/src/main/c/src/fluid/queue/chunk.c index bfe6d07d..e4d2ed7c 100644 --- a/src/main/c/src/fluid/queue/chunk.c +++ b/src/main/c/src/fluid/queue/chunk.c @@ -13,5 +13,22 @@ LIBRARY_API Chunk * chunk_create(){ rVal->vTempCache = (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)); + + /** + * 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; } diff --git a/src/main/c/src/fluid/sim/pressurecell/bounds.c b/src/main/c/src/fluid/sim/pressurecell/bounds.c index 2a9b579e..df96acd5 100644 --- a/src/main/c/src/fluid/sim/pressurecell/bounds.c +++ b/src/main/c/src/fluid/sim/pressurecell/bounds.c @@ -1,9 +1,147 @@ #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 */ 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]); -} \ No newline at end of file + 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; + } + } + } + } +} diff --git a/src/main/c/src/fluid/sim/pressurecell/density.c b/src/main/c/src/fluid/sim/pressurecell/density.c index b3f80219..9b00f1b4 100644 --- a/src/main/c/src/fluid/sim/pressurecell/density.c +++ b/src/main/c/src/fluid/sim/pressurecell/density.c @@ -1,3 +1,4 @@ +#include #include "fluid/sim/pressurecell/density.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; float * densityArr = chunk->d[CENTER_LOC]; float * densityTemp = chunk->dTempCache; + float D = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt; for(z = 1; z < DIM-1; z++){ for(y = 1; y < DIM-1; y++){ for(x = 1; x < DIM-1; x++){ densityTemp[IX(x,y,z)] = densityArr[IX(x,y,z)] + - densityArr[IX(x,y,z)] * -6 * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + - densityArr[IX(x-1,y,z)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + - densityArr[IX(x+1,y,z)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + - densityArr[IX(x,y-1,z)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + - densityArr[IX(x,y+1,z)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + - densityArr[IX(x,y,z-1)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + - densityArr[IX(x,y,z+1)] * FLUID_PRESSURECELL_DIFFUSION_CONSTANT + densityArr[IX(x,y,z)] * -6 * D + + densityArr[IX(x-1,y,z)] * D + + densityArr[IX(x+1,y,z)] * D + + densityArr[IX(x,y-1,z)] * D + + densityArr[IX(x,y+1,z)] * D + + densityArr[IX(x,y,z-1)] * D + + densityArr[IX(x,y,z+1)] * D ; } } @@ -61,9 +63,9 @@ LIBRARY_API void pressurecell_advect_density(Environment * environment, Chunk * 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; - yp = y - v[IX(x,y,z)] * environment->consts.dt; - zp = z - w[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 * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; + zp = z - w[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; //clamp to border x0 = xp; @@ -91,7 +93,12 @@ LIBRARY_API void pressurecell_advect_density(Environment * environment, Chunk * x1 < 0 || y1 < 0 || z1 < 0 || 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); } @@ -109,7 +116,7 @@ LIBRARY_API void pressurecell_advect_density(Environment * environment, Chunk * 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)); } } } diff --git a/src/main/c/src/fluid/sim/pressurecell/normalization.c b/src/main/c/src/fluid/sim/pressurecell/normalization.c new file mode 100644 index 00000000..25ee5ef6 --- /dev/null +++ b/src/main/c/src/fluid/sim/pressurecell/normalization.c @@ -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){ + +} + + diff --git a/src/main/c/src/fluid/sim/pressurecell/pressure.c b/src/main/c/src/fluid/sim/pressurecell/pressure.c index 4b0a2ade..d7208334 100644 --- a/src/main/c/src/fluid/sim/pressurecell/pressure.c +++ b/src/main/c/src/fluid/sim/pressurecell/pressure.c @@ -13,21 +13,63 @@ LIBRARY_API void pressurecell_approximate_pressure(Environment * environment, Ch //values stored across frames float * presureCache = chunk->pressureCache[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 + float * pressureCache = chunk->pressureCache[CENTER_LOC]; float * pressureTemp = chunk->pressureTempCache; - float gridSpacing = 1; + //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+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 = - (divArr[IX(x+1,y,z)] - divArr[IX(x-1,y,z)]) / (2.0f * gridSpacing) + - (divArr[IX(x+1,y,z)] - divArr[IX(x-1,y,z)]) / (2.0f * gridSpacing) + - (divArr[IX(x+1,y,z)] - divArr[IX(x-1,y,z)]) / (2.0f * gridSpacing) - ; - newPressure = newPressure / (gridSpacing * 2); + 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; } } @@ -47,18 +89,38 @@ LIBRARY_API void pressurecell_approximate_divergence(Environment * environment, float * presureCache = chunk->pressureCache[CENTER_LOC]; //temporary caches float * pressureTemp = chunk->pressureTempCache; - float gridSpacing = 1; + float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING; float du, dv, dw; float newDivergence; + 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 caused by static outflow due to diffusion + float outflowDiv = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 6 * dt; + //compute divergence 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); - dw = (wArr[IX(x+1,y,z)] - wArr[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,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; - 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 presureCache[IX(x,y,z)] = pressureTemp[IX(x,y,z)]; diff --git a/src/main/c/src/fluid/sim/pressurecell/pressurecell.c b/src/main/c/src/fluid/sim/pressurecell/pressurecell.c index c70d9e11..e1d1589c 100644 --- a/src/main/c/src/fluid/sim/pressurecell/pressurecell.c +++ b/src/main/c/src/fluid/sim/pressurecell/pressurecell.c @@ -9,11 +9,18 @@ #include "fluid/env/utilities.h" #include "fluid/queue/chunkmask.h" #include "fluid/queue/chunk.h" +#include "fluid/sim/pressurecell/bounds.h" #include "fluid/sim/pressurecell/density.h" #include "fluid/sim/pressurecell/pressure.h" #include "fluid/sim/pressurecell/solver_consts.h" #include "fluid/sim/pressurecell/velocity.h" + + + +static inline void fluid_pressurecell_clearArr(float * d); + + /** * Used for storing timings */ @@ -31,6 +38,16 @@ LIBRARY_API void fluid_pressurecell_simulate( //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 // @@ -96,6 +113,23 @@ LIBRARY_API void fluid_pressurecell_simulate( Chunk * currentChunk = chunks[i]; 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; + } +} diff --git a/src/main/c/src/fluid/sim/pressurecell/velocity.c b/src/main/c/src/fluid/sim/pressurecell/velocity.c index 03146801..6f1a0b31 100644 --- a/src/main/c/src/fluid/sim/pressurecell/velocity.c +++ b/src/main/c/src/fluid/sim/pressurecell/velocity.c @@ -1,4 +1,6 @@ +#include +#include "fluid/queue/chunk.h" #include "fluid/sim/pressurecell/velocity.h" #include "fluid/sim/pressurecell/solver_consts.h" @@ -9,7 +11,6 @@ LIBRARY_API void pressurecell_add_gravity(Environment * environment, Chunk * chunk){ int x, y, z; float * dArr = chunk->d[CENTER_LOC]; - float * vArr = chunk->v[CENTER_LOC]; float * vSourceArr = chunk->v0[CENTER_LOC]; for(z = 1; z < DIM-1; z++){ 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){ int x, y, z; + float * dArr = chunk->d[CENTER_LOC]; 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 D = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt; for(z = 1; z < DIM-1; z++){ for(y = 1; y < DIM-1; y++){ for(x = 1; x < DIM-1; x++){ - //diffuse u uTemp[IX(x,y,z)] = uArr[IX(x,y,z)] + - uArr[IX(x,y,z)] * -6 * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - uArr[IX(x-1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - uArr[IX(x+1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - uArr[IX(x,y-1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - uArr[IX(x,y+1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - uArr[IX(x,y,z-1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - uArr[IX(x,y,z+1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + uArr[IX(x,y,z)] * -6 * D + + uArr[IX(x-1,y,z)] * D + + uArr[IX(x+1,y,z)] * D + + uArr[IX(x,y-1,z)] * D + + uArr[IX(x,y+1,z)] * D + + uArr[IX(x,y,z-1)] * D + + 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)] + - vArr[IX(x,y,z)] * -6 * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - vArr[IX(x-1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - vArr[IX(x+1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - vArr[IX(x,y-1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - vArr[IX(x,y+1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - vArr[IX(x,y,z-1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - vArr[IX(x,y,z+1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + vArr[IX(x,y,z)] * -6 * D + + vArr[IX(x-1,y,z)] * D + + vArr[IX(x+1,y,z)] * D + + vArr[IX(x,y-1,z)] * D + + vArr[IX(x,y+1,z)] * D + + vArr[IX(x,y,z-1)] * D + + 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)] + - wArr[IX(x,y,z)] * -6 * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - wArr[IX(x-1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - wArr[IX(x+1,y,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - wArr[IX(x,y-1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - wArr[IX(x,y+1,z)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - wArr[IX(x,y,z-1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + - wArr[IX(x,y,z+1)] * FLUID_PRESSURECELL_VISCOSITY_CONSTANT + wArr[IX(x,y,z)] * -6 * D + + wArr[IX(x-1,y,z)] * D + + wArr[IX(x+1,y,z)] * D + + wArr[IX(x,y-1,z)] * D + + wArr[IX(x,y+1,z)] * D + + wArr[IX(x,y,z-1)] * D + + 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; float xp, yp, zp; float s0, s1, t0, t1, u0, u1; - float interpolated; + float interpolatedU, interpolatedV, interpolatedW; + float magnitude; 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 - uTemp[IX(x,y,z)] * environment->consts.dt; - yp = y - vTemp[IX(x,y,z)] * environment->consts.dt; - zp = z - wTemp[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 * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; + zp = z - wTemp[IX(x,y,z)] * environment->consts.dt * FLUID_PRESSURECELL_TRANSFERENCE_MULTIPLIER; //clamp to border x0 = xp; @@ -138,21 +150,7 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk * u0 = 1-u1; - 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", x0, y0, z0); - printf("%d %d %d\n ", x1, y1, z1); - printf("%f %f %f\n", x, y, z); - printf("\n"); - fflush(stdout); - } - - interpolated = + interpolatedU = s0*( t0*u0*uTemp[IX(x0,y0,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)]+ t1*u1*uTemp[IX(x1,y1,z1)] ); - uArr[IX(x,y,z)] = interpolated; - - interpolated = + interpolatedV = s0*( t0*u0*vTemp[IX(x0,y0,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)]+ t1*u1*vTemp[IX(x1,y1,z1)] ); - vArr[IX(x,y,z)] = interpolated; - - interpolated = + interpolatedW = s0*( t0*u0*wTemp[IX(x0,y0,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)]+ 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){ int x, y, z; - float * presureCache = chunk->pressureCache[CENTER_LOC]; + float * presureCache = chunk->pressureTempCache; float * uArr = chunk->u[CENTER_LOC]; float * vArr = chunk->v[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 * wTemp = chunk->wTempCache; //temporary caches - float * pressureTemp = chunk->pressureTempCache; - float gridSpacing = 1; + float gridSpacing = FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING; float du, dv, dw; float pressureDivergence; - for(z = 1; z < DIM-1; z++){ - for(y = 1; y < DIM-1; y++){ + float magnitude; + for(y = 1; y < DIM-1; y++){ + for(z = 1; z < DIM-1; z++){ 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 - uTemp[IX(x,y,z)] = uTemp[IX(x,y,z)] - pressureDivergence; - vTemp[IX(x,y,z)] = vTemp[IX(x,y,z)] - pressureDivergence; - wTemp[IX(x,y,z)] = wTemp[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)] - (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); + + // 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 - 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)]; + // if(magnitude > 1.0f){ + // uArr[IX(x,y,z)] = uTemp[IX(x,y,z)] / magnitude; + // 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"); + // } } } } diff --git a/src/main/java/electrosphere/client/script/ScriptClientVoxelUtils.java b/src/main/java/electrosphere/client/script/ScriptClientVoxelUtils.java index cab695fe..ba13b577 100644 --- a/src/main/java/electrosphere/client/script/ScriptClientVoxelUtils.java +++ b/src/main/java/electrosphere/client/script/ScriptClientVoxelUtils.java @@ -66,7 +66,7 @@ public class ScriptClientVoxelUtils { Vector3d centerPos = new Vector3d(CameraEntityUtils.getCameraCenter(camera)); Vector3d cursorPos = collisionEngine.rayCastPosition(new Vector3d(centerPos), new Vector3d(eyePos).mul(-1.0), CollisionEngine.DEFAULT_INTERACT_DISTANCE); 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); Vector3i worldPos = new Vector3i( diff --git a/src/main/java/electrosphere/server/fluid/manager/ServerFluidChunk.java b/src/main/java/electrosphere/server/fluid/manager/ServerFluidChunk.java index 32c39d34..0b649413 100644 --- a/src/main/java/electrosphere/server/fluid/manager/ServerFluidChunk.java +++ b/src/main/java/electrosphere/server/fluid/manager/ServerFluidChunk.java @@ -106,6 +106,11 @@ public class ServerFluidChunk { */ FloatBuffer bounds; + /** + * The pressure cache + */ + FloatBuffer pressureCache; + /** * 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.velocityZ = this.bVelocityZ[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); } + /** + * 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 * @param x The x position diff --git a/src/main/java/electrosphere/server/fluid/manager/ServerFluidManager.java b/src/main/java/electrosphere/server/fluid/manager/ServerFluidManager.java index 7cc41a73..2c290ea3 100644 --- a/src/main/java/electrosphere/server/fluid/manager/ServerFluidManager.java +++ b/src/main/java/electrosphere/server/fluid/manager/ServerFluidManager.java @@ -275,6 +275,7 @@ public class ServerFluidManager { ServerFluidChunk fluidChunk = this.getChunk(worldPos.x, worldPos.y, worldPos.z); fluidChunk.setAsleep(false); fluidChunk.setWeight(voxelPos.x, voxelPos.y, voxelPos.z, weight); + fluidChunk.setPressure(voxelPos.x, voxelPos.y, voxelPos.z, weight); lock.unlock(); } diff --git a/src/main/java/electrosphere/server/fluid/simulator/FluidAcceleratedSimulator.java b/src/main/java/electrosphere/server/fluid/simulator/FluidAcceleratedSimulator.java index 238415bf..1d9c01d5 100644 --- a/src/main/java/electrosphere/server/fluid/simulator/FluidAcceleratedSimulator.java +++ b/src/main/java/electrosphere/server/fluid/simulator/FluidAcceleratedSimulator.java @@ -29,7 +29,7 @@ public class FluidAcceleratedSimulator implements ServerFluidSimulator { /** * The gravity constant */ - public static final float GRAVITY_CONST = -1000f; + public static final float GRAVITY_CONST = -100f; /** * Load fluid sim library diff --git a/src/test/c/fluid/sim/pressurecell/diffuse_tests.c b/src/test/c/fluid/sim/pressurecell/diffuse_tests.c index e0a42c0f..f7ddd449 100644 --- a/src/test/c/fluid/sim/pressurecell/diffuse_tests.c +++ b/src/test/c/fluid/sim/pressurecell/diffuse_tests.c @@ -50,7 +50,7 @@ int fluid_sim_pressurecell_diffuse_test1(){ // // 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)]; 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"); @@ -60,37 +60,37 @@ int fluid_sim_pressurecell_diffuse_test1(){ // // 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)]; 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"); } - 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)]; - 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"); } - 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)]; 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"); } - 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)]; 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"); } - 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)]; 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"); } - 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)]; 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"); @@ -126,7 +126,7 @@ int fluid_sim_pressurecell_diffuse_test2(){ // // 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)]; 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"); @@ -136,37 +136,37 @@ int fluid_sim_pressurecell_diffuse_test2(){ // // 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)]; 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"); } - 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)]; 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"); } - 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)]; 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"); } - 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)]; 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"); } - 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)]; 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"); } - 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)]; 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"); diff --git a/src/test/c/fluid/sim/pressurecell/divergence_tests.c b/src/test/c/fluid/sim/pressurecell/divergence_tests.c new file mode 100644 index 00000000..46607263 --- /dev/null +++ b/src/test/c/fluid/sim/pressurecell/divergence_tests.c @@ -0,0 +1,78 @@ +#include + +#include "stb/stb_ds.h" + +#include "fluid/queue/boundsolver.h" +#include "fluid/queue/chunkmask.h" +#include "fluid/queue/chunk.h" +#include "fluid/env/environment.h" +#include "fluid/env/utilities.h" +#include "fluid/sim/pressurecell/pressure.h" +#include "fluid/sim/pressurecell/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; +} \ No newline at end of file diff --git a/src/test/c/fluid/sim/pressurecell/pressure_tests.c b/src/test/c/fluid/sim/pressurecell/pressure_tests.c index b2cee159..747bb7f3 100644 --- a/src/test/c/fluid/sim/pressurecell/pressure_tests.c +++ b/src/test/c/fluid/sim/pressurecell/pressure_tests.c @@ -39,11 +39,14 @@ int fluid_sim_pressurecell_pressure_test1(){ //setup chunk values float deltaDensity = 0.01f; Chunk * currentChunk = queue[0]; - currentChunk->d[CENTER_LOC][IX(4,4,4)] = MAX_FLUID_VALUE - deltaDensity; - currentChunk->d0[CENTER_LOC][IX(4,4,4)] = deltaDensity; - currentChunk->u0[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY; - currentChunk->v0[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY; - currentChunk->w0[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY; + 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; //actually simulate pressurecell_approximate_pressure(env,currentChunk); @@ -54,11 +57,11 @@ int fluid_sim_pressurecell_pressure_test1(){ // // cell that originall had values // - // expected = MAX_FLUID_VALUE - FLUID_PRESSURECELL_MAX_VELOCITY * env->consts.dt * MAX_FLUID_VALUE; - // actual = currentChunk->pressureCache[CENTER_LOC][IX(4,4,4)]; - // 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"); - // } + expected = -FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER * currentChunk->divergenceCache[CENTER_LOC][IX(4,4,4)]; + actual = currentChunk->pressureTempCache[IX(4,4,4)]; + 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"); + } return rVal; }