studiorailgun/fluid-sim
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

stable-ish multichunk solver
All checks were successful
studiorailgun/fluid-sim/pipeline/head This commit looks good
Details

Browse Source
This commit is contained in:
unknown 2024-03-15 16:56:09 -04:00
parent d4d5d2c94e
commit 7a1da2d29c
9 changed files with 332 additions and 123 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@ -7,4 +7,5 @@
/.project
/shared-folder
/shared-folder/**
/src/main/c/lib/**
/src/main/c/lib/**
/chunks

9
compare.sh Normal file
View File

@ -0,0 +1,9 @@
search_dir=./chunks
for entry in "$search_dir"/*
do
if ! diff "/c/Users/satellite/Documents/fluid-sim/$entry" "/c/Users/satellite/Documents/fluid-sim-the-golden-code/$entry" > /dev/null;
then
echo "“$entry differs"
fi
done
echo "done!"

26
optimizationideas.md Normal file
View File

@ -0,0 +1,26 @@
'sleep' chunks if there are no massive velocity swings for ~100 frames
- unsleep as soon as new density or velocity is added
VVV These two may not work because the projection phase currently iterates over all cells with avx instructions
it also may generally not work as fast moving water basically works by having empty cells query area that already has water (so maybe no sleeping air)
'sleep' individual cells that are full of water and fully surrounded by water
'sleep' terrain cells
'sleep' "inactive" cells
- When reading data back into java, while iterating over each cell check if they've changed value
- If not, add 1 to an accumulating array alongside the other ones (d, u, etc)
- In C, if the accumulating array is >100 or some threshold, skip the cell
- If a force gets added java side that is greater than some threshold, set all accumulating array to 0
multigrid method for projection code
intelligent neighbor awakening --
If a neighbor is asleep (say a neighbor full of air above this cell),
do not simulate that neighbor at all
Instead poll the currently awake cell for large velocity spikes in the direction of the neighbor
If there is a large spike, awake the neighbor cell
if a neighbor is awoken, the simulator should go back through and simulate for that chunk as well (idk if this is possible)
use avx bitmasking for handling terrain in projection code ???

182
src/main/c/fluidsim.c
View File

@ -15,9 +15,13 @@
#define DIM 18
#define LINEARSOLVERTIMES 20
#define SAVE_STEPS 0
#define REALLY_SMALL_VALUE 0.00001
#define DIFFUSION_CONSTANT 0.000001
#define VISCOSITY_CONSTANT 0.000001
#define DIFFUSION_CONSTANT 0.00001
#define VISCOSITY_CONSTANT 0.00001
char fileNameBuff[50];
/**
* A chunk
@ -44,6 +48,8 @@ Chunk ** chunks = NULL;
//jni help:
//https://stackoverflow.com/questions/39823375/clarification-about-getfieldid
void saveStep(float * values, const char * name);
JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
JNIEnv * env,
jclass class,
@ -139,6 +145,12 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
}
// printf("%p\n",chunks[0].d);
saveStep(chunks[0]->u[CENTER_LOC], "./chunks/beginU");
saveStep(chunks[0]->v[CENTER_LOC], "./chunks/beginV");
saveStep(chunks[0]->w[CENTER_LOC], "./chunks/beginW");
saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/beginU0");
saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/beginV0");
saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/beginW0");
//solve chunk mask
for(int i = 0; i < numChunks; i++){
@ -157,6 +169,12 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
VISCOSITY_CONSTANT,
timestep
);
saveStep(currentChunk->u[CENTER_LOC], "./chunks/addSrcU");
saveStep(currentChunk->v[CENTER_LOC], "./chunks/addSrcV");
saveStep(currentChunk->w[CENTER_LOC], "./chunks/addSrcW");
saveStep(currentChunk->u0[CENTER_LOC], "./chunks/addSrcU0");
saveStep(currentChunk->v0[CENTER_LOC], "./chunks/addSrcV0");
saveStep(currentChunk->w0[CENTER_LOC], "./chunks/addSrcW0");
}
//swap all vector fields
{
@ -194,6 +212,16 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
copyNeighborsRaw(DIM,chunkMask,0,3,currentChunk->w0);
}
}
saveStep(chunks[0]->u[CENTER_LOC], "./chunks/swapU");
saveStep(chunks[0]->v[CENTER_LOC], "./chunks/swapV");
saveStep(chunks[0]->w[CENTER_LOC], "./chunks/swapW");
saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/swapU0");
saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/swapV0");
saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/swapW0");
// printf("after swap vecs u\n");
// printLayer(chunks[0]->u[CENTER_LOC],targetLayer);
// printf("after swap vecs u0\n");
// printLayer(chunks[0]->u0[CENTER_LOC],targetLayer);
//solve vector diffusion
{
for(int l = 0; l < LINEARSOLVERTIMES; l++){
@ -203,6 +231,20 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
chunkMask = currentChunk->chunkMask;
Java_electrosphere_FluidSim_solveVectorDiffuse(DIM,chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
}
if(SAVE_STEPS){
sprintf(fileNameBuff, "./chunks/diffuseUStep%dx", l);
saveStep(chunks[0]->u[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseUStep%dx0", l);
saveStep(chunks[0]->u0[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseVStep%dx", l);
saveStep(chunks[0]->v[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseVStep%dx0", l);
saveStep(chunks[0]->v0[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseWStep%dx", l);
saveStep(chunks[0]->w[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseWStep%dx0", l);
saveStep(chunks[0]->w0[CENTER_LOC], fileNameBuff);
}
//update array for vectors
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
@ -210,23 +252,49 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u);
setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v);
setBoundsToNeighborsRaw(DIM,chunkMask,3,currentChunk->w);
copyNeighborsRaw(DIM,chunkMask,0,1,currentChunk->u);
copyNeighborsRaw(DIM,chunkMask,0,2,currentChunk->v);
copyNeighborsRaw(DIM,chunkMask,0,3,currentChunk->w);
// setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u0);
// setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v0);
// setBoundsToNeighborsRaw(DIM,chunkMask,3,currentChunk->w0);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->u);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->v);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->w);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->u0);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->v0);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->w0);
}
if(SAVE_STEPS){
sprintf(fileNameBuff, "./chunks/diffuseUStep%dxBnd", l);
saveStep(chunks[0]->u[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseUStep%dx0Bnd", l);
saveStep(chunks[0]->u0[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseVStep%dxBnd", l);
saveStep(chunks[0]->v[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseVStep%dx0Bnd", l);
saveStep(chunks[0]->v0[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseWStep%dxBnd", l);
saveStep(chunks[0]->w[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/diffuseWStep%dx0Bnd", l);
saveStep(chunks[0]->w0[CENTER_LOC], fileNameBuff);
}
}
}
saveStep(chunks[0]->u[CENTER_LOC], "./chunks/diffuseU");
saveStep(chunks[0]->v[CENTER_LOC], "./chunks/diffuseV");
saveStep(chunks[0]->w[CENTER_LOC], "./chunks/diffuseW");
saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/diffuseU0");
saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/diffuseV0");
saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/diffuseW0");
//solve projection
{
//update array for vectors
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
chunkMask = currentChunk->chunkMask;
setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u);
setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v);
setBoundsToNeighborsRaw(DIM,chunkMask,3,currentChunk->w);
setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u0);
setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v0);
// setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u);
// setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v);
// setBoundsToNeighborsRaw(DIM,chunkMask,3,currentChunk->w);
// setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u0);
// setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v0);
copyNeighborsRaw(DIM,chunkMask,0,1,currentChunk->u);
copyNeighborsRaw(DIM,chunkMask,0,2,currentChunk->v);
copyNeighborsRaw(DIM,chunkMask,0,3,currentChunk->w);
@ -237,17 +305,25 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
chunkMask = currentChunk->chunkMask;
Java_electrosphere_FluidSim_setupProjection(DIM,chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
Java_electrosphere_FluidSim_setupProjection(DIM,chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
}
saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/setupProj1Div");
saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/setupProj1P");
//update array for vectors
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
chunkMask = currentChunk->chunkMask;
setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u0);
setBoundsToNeighborsRaw(DIM,chunkMask,2,currentChunk->v0);
copyNeighborsRaw(DIM,chunkMask,0,1,currentChunk->u0);
copyNeighborsRaw(DIM,chunkMask,0,2,currentChunk->v0);
setBoundsToNeighborsRaw(DIM,chunkMask,0,currentChunk->u0);
setBoundsToNeighborsRaw(DIM,chunkMask,0,currentChunk->v0);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->u0);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->v0);
}
saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/setupProj1DivBnd");
saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/setupProj1PBnd");
//samples u0, v0
//sets u0
//these should have just been mirrored in the above
@ -259,11 +335,23 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
chunkMask = currentChunk->chunkMask;
Java_electrosphere_FluidSim_solveProjection(DIM,chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
}
if(SAVE_STEPS){
sprintf(fileNameBuff, "./chunks/proj1Step%dx", l);
saveStep(chunks[0]->u0[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/proj1Step%dx0", l);
saveStep(chunks[0]->v0[CENTER_LOC], fileNameBuff);
}
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
chunkMask = currentChunk->chunkMask;
setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u0);
copyNeighborsRaw(DIM,chunkMask,0,1,currentChunk->u0);
setBoundsToNeighborsRaw(DIM,chunkMask,0,currentChunk->u0);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->u0);
}
if(SAVE_STEPS){
sprintf(fileNameBuff, "./chunks/proj1Step%dxBnd", l);
saveStep(chunks[0]->u0[CENTER_LOC], fileNameBuff);
sprintf(fileNameBuff, "./chunks/proj1Step%dx0Bnd", l);
saveStep(chunks[0]->v0[CENTER_LOC], fileNameBuff);
}
}
//samples u,v,w,u0
@ -274,6 +362,11 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
chunkMask = currentChunk->chunkMask;
Java_electrosphere_FluidSim_finalizeProjection(DIM,chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
}
saveStep(chunks[0]->u[CENTER_LOC], "./chunks/finalizeProj1U");
saveStep(chunks[0]->v[CENTER_LOC], "./chunks/finalizeProj1V");
saveStep(chunks[0]->w[CENTER_LOC], "./chunks/finalizeProj1W");
// exit(0);
//set boundaries a final time for u,v,w
//...
for(int i = 0; i < numChunks; i++){
@ -293,6 +386,13 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
copyNeighborsRaw(DIM,chunkMask,0,3,currentChunk->w0);
}
}
saveStep(chunks[0]->u[CENTER_LOC], "./chunks/projU");
saveStep(chunks[0]->v[CENTER_LOC], "./chunks/projV");
saveStep(chunks[0]->w[CENTER_LOC], "./chunks/projW");
saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/projU0");
saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/projV0");
saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/projW0");
// exit(0);
//swap all vector fields
{
//swap vector fields
@ -329,6 +429,12 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
copyNeighborsRaw(DIM,chunkMask,0,3,currentChunk->w0);
}
}
saveStep(chunks[0]->u[CENTER_LOC], "./chunks/swap2U");
saveStep(chunks[0]->v[CENTER_LOC], "./chunks/swap2V");
saveStep(chunks[0]->w[CENTER_LOC], "./chunks/swap2W");
saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/swap2U0");
saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/swap2V0");
saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/swap2W0");
//advect vectors across boundaries
{
//update border arrs
@ -366,6 +472,12 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
copyNeighborsRaw(DIM,chunkMask,0,3,currentChunk->w);
}
}
saveStep(chunks[0]->u[CENTER_LOC], "./chunks/advectU");
saveStep(chunks[0]->v[CENTER_LOC], "./chunks/advectV");
saveStep(chunks[0]->w[CENTER_LOC], "./chunks/advectW");
saveStep(chunks[0]->u0[CENTER_LOC], "./chunks/advectU0");
saveStep(chunks[0]->v0[CENTER_LOC], "./chunks/advectV0");
saveStep(chunks[0]->w0[CENTER_LOC], "./chunks/advectW0");
//solve projection
{
//update array for vectors
@ -387,7 +499,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
chunkMask = currentChunk->chunkMask;
Java_electrosphere_FluidSim_setupProjection(DIM,chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
Java_electrosphere_FluidSim_setupProjection(DIM,chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,DIFFUSION_CONSTANT,VISCOSITY_CONSTANT,timestep);
}
//update array for vectors
for(int i = 0; i < numChunks; i++){
@ -412,8 +524,8 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
chunkMask = currentChunk->chunkMask;
setBoundsToNeighborsRaw(DIM,chunkMask,1,currentChunk->u0);
copyNeighborsRaw(DIM,chunkMask,0,1,currentChunk->u0);
setBoundsToNeighborsRaw(DIM,chunkMask,0,currentChunk->u0);
copyNeighborsRaw(DIM,chunkMask,0,0,currentChunk->u0);
}
}
//samples u,v,w,u0
@ -534,4 +646,34 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
setBoundsToNeighborsRaw(DIM,chunkMask,0,currentChunk->d);
}
}
}
void saveStep(float * values, const char * name){
if(SAVE_STEPS){
FILE *fp;
int N = DIM;
// ... fill the array somehow ...
fp = fopen(name, "w");
// check for error here
for(int x = 0; x < DIM; x++){
for(int y = 0; y < DIM; y++){
for(int z = 0; z < DIM; z++){
float val = values[IX(x,y,z)];
if(val < REALLY_SMALL_VALUE && val > -REALLY_SMALL_VALUE){
val = 0;
}
fprintf(fp, "%f\t", val);
}
fprintf(fp, "\n");
}
fprintf(fp, "\n");
}
fclose(fp);
}
}

6
src/main/c/includes/electrosphere_FluidSim.h
View File

@ -10,13 +10,13 @@ extern "C" {
#undef electrosphere_FluidSim_DIM
#define electrosphere_FluidSim_DIM 18L
#undef electrosphere_FluidSim_DIFFUSION_CONSTANT
#define electrosphere_FluidSim_DIFFUSION_CONSTANT 0.0f
#define electrosphere_FluidSim_DIFFUSION_CONSTANT 1.0E-5f
#undef electrosphere_FluidSim_VISCOSITY_CONSTANT
#define electrosphere_FluidSim_VISCOSITY_CONSTANT 0.0f
#define electrosphere_FluidSim_VISCOSITY_CONSTANT 1.0E-5f
#undef electrosphere_FluidSim_LINEARSOLVERTIMES
#define electrosphere_FluidSim_LINEARSOLVERTIMES 20L
#undef electrosphere_FluidSim_GRAVITY
#define electrosphere_FluidSim_GRAVITY -100.0f
#define electrosphere_FluidSim_GRAVITY -1000.0f
/*
* Class: electrosphere_FluidSim
* Method: simulate

12
src/main/c/includes/mainFunctions.h
View File

@ -33,7 +33,17 @@ void Java_electrosphere_FluidSim_solveVectorDiffuse
* Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
*/
void Java_electrosphere_FluidSim_setupProjection
(int, int, float **, float **, float **, float **, float **, float **, float, float, float);
(
int N,
int chunk_mask,
float ** ur,
float ** vr,
float ** wr,
float ** pr,
float ** divr,
float DIFFUSION_CONST,
float VISCOSITY_CONST,
float dt);
/*
* Class: electrosphere_FluidSim

88
src/main/c/velocitystep.c
View File

@ -14,6 +14,8 @@
#define SET_BOUND_IGNORE 0
#define SET_BOUND_USE_NEIGHBOR 1
#define LINEARSOLVERTIMES 20
void add_source(int N, float * x, float * s, float dt);
void advect(uint32_t chunk_mask, int N, int b, float ** jrd, float ** jrd0, float * u, float * v, float * w, float dt);
@ -162,30 +164,27 @@ void Java_electrosphere_FluidSim_setupProjection
(
int N,
int chunk_mask,
float ** jru,
float ** jrv,
float ** jrw,
float ** jru0,
float ** jrv0,
float ** jrw0,
float ** ur,
float ** vr,
float ** wr,
float ** pr,
float ** divr,
float DIFFUSION_CONST,
float VISCOSITY_CONST,
float dt){
int i, j, k;
__m256 xVector = _mm256_set1_ps(N);
__m256 yVector = _mm256_set1_ps(N);
__m256 zVector = _mm256_set1_ps(N);
__m256 nVector = _mm256_set1_ps(N);
__m256 constScalar = _mm256_set1_ps(-1.0/3.0);
__m256 zeroVec = _mm256_set1_ps(0);
__m256 vector, vector2, vector3;
float * u = GET_ARR_RAW(env,jru,CENTER_LOC);
float * v = GET_ARR_RAW(env,jrv,CENTER_LOC);
float * w = GET_ARR_RAW(env,jrw,CENTER_LOC);
float * u = GET_ARR_RAW(env,ur,CENTER_LOC);
float * v = GET_ARR_RAW(env,vr,CENTER_LOC);
float * w = GET_ARR_RAW(env,wr,CENTER_LOC);
float * p = GET_ARR_RAW(env,jru0,CENTER_LOC);
float * div = GET_ARR_RAW(env,jrv0,CENTER_LOC);
float * p = GET_ARR_RAW(env,pr,CENTER_LOC);
float * div = GET_ARR_RAW(env,divr,CENTER_LOC);
float scalar = 1.0/3.0;
float h = 1.0/N;
@ -199,15 +198,15 @@ void Java_electrosphere_FluidSim_setupProjection
//first part
vector = _mm256_loadu_ps(&u[IX(i+1,j,k)]);
vector = _mm256_sub_ps(vector,_mm256_loadu_ps(&u[IX(i-1,j,k)]));
vector = _mm256_div_ps(vector,xVector);
vector = _mm256_div_ps(vector,nVector);
//second part
vector2 = _mm256_loadu_ps(&v[IX(i,j+1,k)]);
vector2 = _mm256_sub_ps(vector2,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
vector2 = _mm256_div_ps(vector2,yVector);
vector2 = _mm256_div_ps(vector2,nVector);
//third part
vector3 = _mm256_loadu_ps(&w[IX(i,j,k+1)]);
vector3 = _mm256_sub_ps(vector3,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
vector3 = _mm256_div_ps(vector3,zVector);
vector3 = _mm256_div_ps(vector3,nVector);
//multiply and finalize
vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
vector = _mm256_mul_ps(vector,constScalar);
@ -221,15 +220,15 @@ void Java_electrosphere_FluidSim_setupProjection
//first part
vector = _mm256_loadu_ps(&u[IX(i+1,j,k)]);
vector = _mm256_sub_ps(vector,_mm256_loadu_ps(&u[IX(i-1,j,k)]));
vector = _mm256_div_ps(vector,xVector);
vector = _mm256_div_ps(vector,nVector);
//second part
vector2 = _mm256_loadu_ps(&v[IX(i,j+1,k)]);
vector2 = _mm256_sub_ps(vector2,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
vector2 = _mm256_div_ps(vector2,yVector);
vector2 = _mm256_div_ps(vector2,nVector);
//third part
vector3 = _mm256_loadu_ps(&w[IX(i,j,k+1)]);
vector3 = _mm256_sub_ps(vector3,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
vector3 = _mm256_div_ps(vector3,zVector);
vector3 = _mm256_div_ps(vector3,nVector);
//multiply and finalize
vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
vector = _mm256_mul_ps(vector,constScalar);
@ -285,7 +284,7 @@ void Java_electrosphere_FluidSim_solveProjection
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j+1,k)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j,k-1)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i,j,k+1)]));
// vector = _mm256_mul_ps(vector,aScalar);
vector = _mm256_mul_ps(vector,aScalar);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&div[IX(i,j,k)]));
vector = _mm256_div_ps(vector,cScalar);
_mm256_storeu_ps(&p[IX(i,j,k)],vector);
@ -320,10 +319,7 @@ void Java_electrosphere_FluidSim_finalizeProjection
float VISCOSITY_CONST,
float dt){
int i, j, k;
// __m256 constScalar = _mm256_set1_ps(0.5f*N);
__m256 xScalar = _mm256_set1_ps(0.5*N);
__m256 yScalar = _mm256_set1_ps(0.5*N);
__m256 zScalar = _mm256_set1_ps(0.5*N);
__m256 constScalar = _mm256_set1_ps(0.5f*N);
__m256 vector, vector2, vector3;
float * u = GET_ARR_RAW(env,jru,CENTER_LOC);
@ -333,8 +329,6 @@ void Java_electrosphere_FluidSim_finalizeProjection
float * p = GET_ARR_RAW(env,jru0,CENTER_LOC);
float * div = GET_ARR_RAW(env,jrv0,CENTER_LOC);
float h = 1.0 / N;
for ( k=1 ; k<N-1 ; k++ ) {
for ( j=1 ; j<N-1 ; j++ ) {
//
@ -344,14 +338,14 @@ void Java_electrosphere_FluidSim_finalizeProjection
vector = _mm256_loadu_ps(&p[IX(1+1,j,k)]);
vector2 = _mm256_loadu_ps(&p[IX(1-1,j,k)]);
vector = _mm256_sub_ps(vector,vector2);
vector = _mm256_mul_ps(vector,xScalar);
vector = _mm256_mul_ps(vector,constScalar);
vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(1,j,k)]),vector);
_mm256_storeu_ps(&u[IX(1,j,k)],vector);
//upper
vector = _mm256_loadu_ps(&p[IX(9+1,j,k)]);
vector2 = _mm256_loadu_ps(&p[IX(9-1,j,k)]);
vector = _mm256_sub_ps(vector,vector2);
vector = _mm256_mul_ps(vector,xScalar);
vector = _mm256_mul_ps(vector,constScalar);
vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(9,j,k)]),vector);
_mm256_storeu_ps(&u[IX(9,j,k)],vector);
//
@ -361,14 +355,14 @@ void Java_electrosphere_FluidSim_finalizeProjection
vector = _mm256_loadu_ps(&p[IX(1,j+1,k)]);
vector2 = _mm256_loadu_ps(&p[IX(1,j-1,k)]);
vector = _mm256_sub_ps(vector,vector2);
vector = _mm256_mul_ps(vector,yScalar);
vector = _mm256_mul_ps(vector,constScalar);
vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(1,j,k)]),vector);
_mm256_storeu_ps(&v[IX(1,j,k)],vector);
//upper
vector = _mm256_loadu_ps(&p[IX(9,j+1,k)]);
vector2 = _mm256_loadu_ps(&p[IX(9,j-1,k)]);
vector = _mm256_sub_ps(vector,vector2);
vector = _mm256_mul_ps(vector,yScalar);
vector = _mm256_mul_ps(vector,constScalar);
vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(9,j,k)]),vector);
_mm256_storeu_ps(&v[IX(9,j,k)],vector);
//
@ -378,21 +372,16 @@ void Java_electrosphere_FluidSim_finalizeProjection
vector = _mm256_loadu_ps(&p[IX(1,j,k+1)]);
vector2 = _mm256_loadu_ps(&p[IX(1,j,k-1)]);
vector = _mm256_sub_ps(vector,vector2);
vector = _mm256_mul_ps(vector,zScalar);
vector = _mm256_mul_ps(vector,constScalar);
vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(1,j,k)]),vector);
_mm256_storeu_ps(&w[IX(1,j,k)],vector);
//upper
vector = _mm256_loadu_ps(&p[IX(9,j,k+1)]);
vector2 = _mm256_loadu_ps(&p[IX(9,j,k-1)]);
vector = _mm256_sub_ps(vector,vector2);
vector = _mm256_mul_ps(vector,zScalar);
vector = _mm256_mul_ps(vector,constScalar);
vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(9,j,k)]),vector);
_mm256_storeu_ps(&w[IX(9,j,k)],vector);
// for(i = 1; i < N-1; i++){
// u[IX(i,j,k)] = u[IX(i,j,k)] - 0.5 * (p[IX(i+1,j,k)] - p[IX(i-1,j,k)]) / h;
// v[IX(i,j,k)] = v[IX(i,j,k)] - 0.5 * (p[IX(i,j+1,k)] - p[IX(i,j-1,k)]) / h;
// w[IX(i,j,k)] = w[IX(i,j,k)] - 0.5 * (p[IX(i,j,k+1)] - p[IX(i,j,k-1)]) / h;
// }
}
}
}
@ -670,19 +659,22 @@ void setBoundsToNeighborsRaw
int DIM = N;
float * target = GET_ARR_RAW(env,neighborArray,CENTER_LOC);
float * source;
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
target[IX(0,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(1,x,y)] : target[IX(1,x,y)];
}
}
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
target[IX(DIM-1,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)];
}
}
// for(int x=1; x < DIM-1; x++){
// for(int y = 1; y < DIM-1; y++){
// target[IX(0,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(1,x,y)] : target[IX(1,x,y)];
// }
// }
// for(int x=1; x < DIM-1; x++){
// for(int y = 1; y < DIM-1; y++){
// target[IX(DIM-1,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)];
// }
// }
//set the faces bounds
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
//((x)+(DIM)*(y) + (DIM)*(DIM)*(z))
target[IX(0,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(1,x,y)] : target[IX(1,x,y)];
target[IX(DIM-1,x,y)] = vector_dir==BOUND_DIR_U ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)];
target[IX(x,0,y)] = vector_dir==BOUND_DIR_V ? -target[IX(x,1,y)] : target[IX(x,1,y)];
target[IX(x,DIM-1,y)] = vector_dir==BOUND_DIR_V ? -target[IX(x,DIM-2,y)] : target[IX(x,DIM-2,y)];
target[IX(x,y,0)] = vector_dir==BOUND_DIR_W ? -target[IX(x,y,1)] : target[IX(x,y,1)];

99
src/main/java/electrosphere/FluidSim.java
View File

@ -73,12 +73,12 @@ public class FluidSim {
public int chunkMask = 0;
static final float DIFFUSION_CONSTANT = 0.0f;
static final float VISCOSITY_CONSTANT = 0.0f;
static final float DIFFUSION_CONSTANT = 0.00001f;
static final float VISCOSITY_CONSTANT = 0.00001f;
static final int LINEARSOLVERTIMES = 20;
static final float GRAVITY = -100f;
static final float GRAVITY = -1000f;
public void setup(Vector3i offset){
//allocate buffers for this chunk
@ -111,45 +111,47 @@ public class FluidSim {
for(int i = 0; i < DIM; i++){
for(int j = 0; j < DIM; j++){
for(int k = 0; k < DIM; k++){
if(offset.x == 0 && offset.y == 0 && offset.z == 0){
// if(offset.x == 0 && offset.y == 0 && offset.z == 0){
if(
Math.abs(16 - i) < 5 &&
Math.abs(j) < 5 &&
Math.abs(16 - k) < 5 &&
i < 17 && i > 0 &&
j < 17 && j > 0 &&
k < 17 && k > 0
Math.abs(16 - i) < 4 &&
Math.abs(8 - j) < 4 &&
Math.abs(10 - k) < 4
// &&
// i < 17 && i > 0 &&
// j < 17 && j > 0 &&
// k < 17 && k > 0
){
xf.put(1);
uf.put(50);
vf.put(0);
wf.put(0);
} else {
uf.put(0.1f);
vf.put(-1);
wf.put(0.1f);
}
else {
xf.put(0);
uf.put(0);
vf.put(0);
wf.put(0);
}
} else {
if(
Math.abs(0 - i) < 5 &&
Math.abs(j) < 5 &&
Math.abs(0 - k) < 5 &&
i < 17 && i > 0 &&
j < 17 && j > 0 &&
k < 17 && k > 0
){
// xf.put(1);
// uf.put(50);
// vf.put(0);
// wf.put(rand.nextFloat() * 0.1f);
} else {
xf.put(0);
uf.put(0);
vf.put(0);
wf.put(0);
}
}
// } else {
// if(
// Math.abs(0 - i) < 5 &&
// Math.abs(j) < 5 &&
// Math.abs(0 - k) < 5 &&
// i < 17 && i > 0 &&
// j < 17 && j > 0 &&
// k < 17 && k > 0
// ){
// // xf.put(1);
// // uf.put(50);
// // vf.put(0);
// // wf.put(rand.nextFloat() * 0.1f);
// } else {
// xf.put(0);
// uf.put(0);
// vf.put(0);
// wf.put(0);
// }
// }
}
}
}
@ -160,9 +162,9 @@ public class FluidSim {
static double lastTime = 0;
public static void simChunks(FluidSim[][][] simArray, int step, float timestep){
simArray[0][1][2].density0ArrayView[IX(10,10,3)] = 3.0f;
simArray[2][1][2].density0ArrayView[IX(10,10,3)] = 3.0f;
simArray[2][1][0].density0ArrayView[IX(10,10,3)] = 3.0f;
// simArray[0][1][2].density0ArrayView[IX(10,10,3)] = 3.0f;
// simArray[2][1][2].density0ArrayView[IX(10,10,3)] = 3.0f;
// simArray[2][1][0].density0ArrayView[IX(10,10,3)] = 3.0f;
List<FluidSim> chunksToSim = new LinkedList<FluidSim>();
//
@ -189,7 +191,7 @@ public class FluidSim {
lastTime = GLFW.glfwGetTime();
//
//simulate
simulateWrapper(chunksToSim,timestep);
simulateWrapper(chunksToSim,0.01f);
//clock
time = time + (GLFW.glfwGetTime() - lastTime);
i++;
@ -197,8 +199,6 @@ public class FluidSim {
System.out.println(time / 100.0 * 1000.0);
}
//
//Read out of buffers back into accessible arrays
@ -212,6 +212,25 @@ public class FluidSim {
}
}
}
if(Main.endStep == 0){
System.exit(0);
}
}
static void printLayer(FluidSim[][][] simArray, int step, int layer){
if(step == 0){
for(int x = 0; x < DIM; x++){
for(int y = 0; y < DIM; y++){
System.out.printf("%.2f\t",simArray[0][0][0].uArrayView[IX(x,layer,y)]);
}
System.out.println();
}
}
}
private static double sumAllDensity(FluidSim[][][] simArray){

30
src/main/java/electrosphere/Main.java
View File

@ -17,7 +17,8 @@ import electrosphere.render.Mesh;
*/
public class Main {
static boolean render = true;
public static int endStep = -1;
public static void main(String args[]){
@ -31,15 +32,20 @@ public class Main {
try {
GLFWContext.init();
if(render){
GLFWContext.init();
//init shader program
Mesh.initShaderProgram();
//init shader program
Mesh.initShaderProgram();
}
FluidSim[][][] simArray = initFluidSim(dim,vdim,dim);
Mesh[][][] meshArray = initMeshes(dim,vdim,dim,simArray);
Mesh[][][] meshArray = null;
if(render){
meshArray = initMeshes(dim,vdim,dim,simArray);
}
@ -60,16 +66,20 @@ public class Main {
//
//Remesh
//
for(int x = 0; x < simArray.length; x++){
for(int y = 0; y < simArray[0].length; y++){
for(int z = 0; z < simArray[0][0].length; z++){
meshArray[x][y][z].remesh();
if(render){
for(int x = 0; x < simArray.length; x++){
for(int y = 0; y < simArray[0].length; y++){
for(int z = 0; z < simArray[0][0].length; z++){
meshArray[x][y][z].remesh();
}
}
}
}
time = time + (System.currentTimeMillis() - lastTime);
//redraw
GLFWContext.redraw(meshArray);
if(render){
GLFWContext.redraw(meshArray);
}
i++;
if(i == 100){
System.out.println("overall time: " + time / 100.0);