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

Compare commits

base: studiorailgun:eae31e2ec3b6df438489c7a2f9eb3adc2db6d08b
Branches Tags
studiorailgun:master
studiorailgun:jenkinsAutoTesting
studiorailgun:migrate-multichunk-to-c
studiorailgun:migration-successful
studiorailgun:multichunking-success
studiorailgun:c-single-thread
..
compare: studiorailgun:37a9afa06934527e93e10457b3ead1f2fd63d8d1
Branches Tags
studiorailgun:master
studiorailgun:jenkinsAutoTesting
studiorailgun:migrate-multichunk-to-c
studiorailgun:migration-successful
studiorailgun:multichunking-success
studiorailgun:c-single-thread

No commits in common. "eae31e2ec3b6df438489c7a2f9eb3adc2db6d08b" and "37a9afa06934527e93e10457b3ead1f2fd63d8d1" have entirely different histories.
eae31e2ec3 ... 37a9afa069

14 changed files with 680 additions and 945 deletions
Show Stats Expand all files Collapse all files

1
.gitignore vendored
View File

@ -8,4 +8,3 @@
/.vscode
/shared-folder
/shared-folder/**
/src/main/c/lib/**

3
.gitmodules vendored
View File

@ -1,3 +0,0 @@
[submodule "src/main/c/lib/stb"]
path = src/main/c/lib/stb
url = https://github.com/nothings/stb.git

17
Jenkinsfile vendored
View File

@ -1,17 +0,0 @@
pipeline {
agent any
tools {
maven '3.9.6'
}
stages {
stage('Build') {
steps {
sh 'mvn --version'
sh 'java -version'
sh 'apt-get update'
sh 'apt-get install gcc'
sh 'mvn clean generate-resources package'
}
}
}
}

47
src/main/c/src/chunkmask.c → src/main/c/chunkmask.c
View File

@ -1,8 +1,7 @@
#include <jni.h>
#include <stdint.h>
#include "../includes/libfluidsim.h"
#include "../includes/utilities.h"
#include "../includes/chunkmask.h"
#include "includes/utilities.h"
#include "includes/chunkmask.h"
uint32_t matrix_transform(JNIEnv * env, jobjectArray jrx);
@ -61,43 +60,43 @@ const uint32_t CHUNK_INDEX_ARR[] = {
//control offsetting the advect sampler location if a valid neighbor chunk is hit
const char CHUNK_NORMALIZE_U[] = {
1, 0, -1,
1, 0, -1,
1, 0, -1,
-1, 0, 1,
-1, 0, 1,
-1, 0, 1,
1, 0, -1,
1, 0, -1,
1, 0, -1,
-1, 0, 1,
-1, 0, 1,
-1, 0, 1,
1, 0, -1,
1, 0, -1,
1, 0, -1,
-1, 0, 1,
-1, 0, 1,
-1, 0, 1,
};
const char CHUNK_NORMALIZE_V[] = {
-1, -1, -1,
0, 0, 0,
1, 1, 1,
0, 0, 0,
-1, -1, -1,
1, 1, 1,
0, 0, 0,
-1, -1, -1,
0, 0, 0,
1, 1, 1,
1, 1, 1,
0, 0, 0,
-1, -1, -1,
0, 0, 0,
1, 1, 1,
};
const char CHUNK_NORMALIZE_W[] = {
1, 1, 1,
1, 1, 1,
1, 1, 1,
-1, -1, -1,
-1, -1, -1,
-1, -1, -1,
0, 0, 0,
0, 0, 0,
0, 0, 0,
-1, -1, -1,
-1, -1, -1,
-1, -1, -1,
1, 1, 1,
1, 1, 1,
1, 1, 1,
};

11
src/main/c/compile.sh
View File

@ -42,17 +42,22 @@ rm -f ./*.dll
#compile object files
COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
INPUT_FILES="./src/densitystep.c"
INPUT_FILES="./fluidsim.c"
OUTPUT_FILE="./fluidsim.o"
gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
INPUT_FILES="./densitystep.c"
OUTPUT_FILE="./densitystep.o"
gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
INPUT_FILES="./src/velocitystep.c"
INPUT_FILES="./velocitystep.c"
OUTPUT_FILE="./velocitystep.o"
gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE
COMPILE_FLAGS="-c -fPIC -m64 -mavx -mavx2 -O1"
INPUT_FILES="./src/chunkmask.c"
INPUT_FILES="./chunkmask.c"
OUTPUT_FILE="./chunkmask.o"
gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OUTPUT_FILE

100
src/main/c/src/densitystep.c → src/main/c/densitystep.c
View File

@ -2,13 +2,11 @@
#include <stdio.h>
#include <immintrin.h>
#include <stdint.h>
#include <pthread.h>
#include "../includes/libfluidsim.h"
#include "../includes/utilities.h"
#include "../includes/chunkmask.h"
#include "includes/utilities.h"
#include "includes/chunkmask.h"
void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, jobjectArray d0, float * u, float * v, float * w, float dt);
void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt);
/*
* Class: electrosphere_FluidSim
@ -21,12 +19,12 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_addDensity
jint N,
jint chunk_mask,
jobjectArray jrx,
jobjectArray x0,
jobject x0,
jfloat dt){
int i;
int size=N*N*N;
float * x = GET_ARR(env,jrx,CENTER_LOC);
float * s = GET_ARR(env,x0,CENTER_LOC);
float * s = (*env)->GetDirectBufferAddress(env,x0);
for(i=0; i<size; i++){
x[i] += dt*s[i];
}
@ -43,7 +41,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
jint N,
jint chunk_mask,
jobjectArray jrx,
jobjectArray jrx0,
jobject jrx0,
jobjectArray jru,
jobjectArray jrv,
jobjectArray jrw,
@ -54,7 +52,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
float c=1+6*a;
int i, j, k, l, m;
float * x = GET_ARR(env,jrx,CENTER_LOC);
float * x0 = GET_ARR(env,jrx0,CENTER_LOC);
float * x0 = (*env)->GetDirectBufferAddress(env,jrx0);
__m256 aScalar = _mm256_set1_ps(a);
__m256 cScalar = _mm256_set1_ps(c);
@ -97,17 +95,18 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectDensity
jint N,
jint chunk_mask,
jobjectArray jrx,
jobjectArray jrx0,
jobject jrx0,
jobjectArray jru,
jobjectArray jrv,
jobjectArray jrw,
jfloat DIFFUSION_CONST,
jfloat VISCOSITY_CONST,
jfloat dt){
advectDensity(env,chunk_mask,N,3,jrx,jrx0,GET_ARR(env,jru,CENTER_LOC),GET_ARR(env,jrv,CENTER_LOC),GET_ARR(env,jrw,CENTER_LOC),dt);
float * x0 = (*env)->GetDirectBufferAddress(env,jrx0);
advectDensity(env,chunk_mask,N,3,jrx,x0,GET_ARR(env,jru,CENTER_LOC),GET_ARR(env,jrv,CENTER_LOC),GET_ARR(env,jrw,CENTER_LOC),dt);
}
void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, jobjectArray jrd0, float * u, float * v, float * w, float dt){
void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt){
int i, j, k, i0, j0, k0, i1, j1, k1;
int m,n,o;
float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz;
@ -116,12 +115,12 @@ void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray
float * d = GET_ARR(env,jrd,CENTER_LOC);
float * d0 = GET_ARR(env,jrd0,CENTER_LOC);
float * sampleArr = d0;
for(k=1; k<N-1; k++){
for(j=1; j<N-1; j++){
for(i=1; i<N-1; i++){
d0 = GET_ARR(env,jrd0,CENTER_LOC);
sampleArr = d0;
//calculate location to pull from
x = i-dtx*u[IX(i,j,k)];
y = j-dty*v[IX(i,j,k)];
@ -129,12 +128,12 @@ void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray
m = n = o = 1;
// if(x < 1){ m -= 1; }
// if(x >= N-1){ m += 1; }
// if(y < 1){ n -= 1; }
// if(y >= N-1){ n += 1; }
// if(z < 1){ o -= 1; }
// if(z >= N-1){ o += 1; }
if(x < 1){ m -= 1; }
if(x >= N-1){ m += 1; }
if(y < 1){ n -= 1; }
if(y >= N-1){ n += 1; }
if(z < 1){ o -= 1; }
if(z >= N-1){ o += 1; }
//If the out of bounds coordinate is in bounds for a neighbor chunk, use that chunk as source instead
// if(CK(m,n,o) != CENTER_LOC){
@ -145,56 +144,31 @@ void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray
// }
// if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
// // printf("Hit other chunk\n");
// d0 = GET_ARR(env,jrd0,CK(m,n,o));
// x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * (N-1);
// y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * (N-1);
// z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * (N-1);
// sampleArr = GET_ARR(env,jrd,CK(m,n,o));
// x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * N;
// y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * N;
// z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * N;
// }
if(x < 0.001f){
//cases to consider:
//m = 0, x = -10
//m = 2, x = 0.01
x=0.001f;
i0=(int)0;
i1=1;
s0 = 0.999f;
s1 = 0.001f;
} else if(x >= N - 1){
//cases to consider:
//m = 0, x = 17.01
//m = 2, x = 20
x = N-1;
i0=(int)N-2;
i1=N-1;
s0 = 0.001f;
s1 = 0.999f;
} else {
i0=(int)x;
i1=i0+1;
s1 = x-i0;
s0 = 1-s1;
}
//clamp location within chunk
// if (x<0.5f) x=0.5f;
// if (x>N+0.5f) x=N+0.5f;
if (x<0.5f) x=0.5f;
if (x>N+0.5f) x=N+0.5f;
if (y<0.5f) y=0.5f;
if (y>N+0.5f) y=N+0.5f;
if (z<0.5f) z=0.5f;
if (z>N+0.5f) z=N+0.5f;
//get actual indices
// i0=(int)x;
// i1=i0+1;
i0=(int)x;
i1=i0+1;
j0=(int)y;
j1=j0+1;
k0=(int)z;
k1=k0+1;
//calculate percentage of each index
// s1 = x-i0;
// s0 = 1-s1;
s1 = x-i0;
s0 = 1-s1;
t1 = y-j0;
t0 = 1-t1;
u1 = z-k0;
@ -238,16 +212,16 @@ void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray
// }
d[IX(i,j,k)] =
s0*(
t0*u0*d0[IX(i0,j0,k0)]+
t1*u0*d0[IX(i0,j1,k0)]+
t0*u1*d0[IX(i0,j0,k1)]+
t1*u1*d0[IX(i0,j1,k1)]
t0*u0*sampleArr[IX(i0,j0,k0)]+
t1*u0*sampleArr[IX(i0,j1,k0)]+
t0*u1*sampleArr[IX(i0,j0,k1)]+
t1*u1*sampleArr[IX(i0,j1,k1)]
)+
s1*(
t0*u0*d0[IX(i1,j0,k0)]+
t1*u0*d0[IX(i1,j1,k0)]+
t0*u1*d0[IX(i1,j0,k1)]+
t1*u1*d0[IX(i1,j1,k1)]
t0*u0*sampleArr[IX(i1,j0,k0)]+
t1*u0*sampleArr[IX(i1,j1,k0)]+
t0*u1*sampleArr[IX(i1,j0,k1)]+
t1*u1*sampleArr[IX(i1,j1,k1)]
);
}
}

396
src/main/c/fluidsim.c Normal file
View File

@ -0,0 +1,396 @@
#include <jni.h>
#include <stdio.h>
#include <immintrin.h>
#include <stdint.h>
#include "includes/utilities.h"
#include "includes/chunkmask.h"
#define LINEARSOLVERTIMES 10
void diffuse(JNIEnv * env, uint32_t chunk_mask, int N, int b, float * x, float * x0, float diff, float dt);
void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt);
void project(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float * p, float * div);
void set_bnd(JNIEnv * env, uint32_t chunk_mask, int N, int b, float * x);
void dens_step(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jrx, float * x0, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float diff, float dt);
void vel_step(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float * u0, float * v0, float * w0, float visc, float dt);
void lin_solve(JNIEnv * env, uint32_t chunk_mask, int N, int b, float* x, float* x0, float a, float c);
/**
* The core simulation function
*/
JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate(
JNIEnv * env,
jobject this,
jint DIM_X,
jint chunk_mask,
jobject jx,
jobject jx0,
jobject ju,
jobject jv,
jobject jw,
jobject ju0,
jobject jv0,
jobject jw0,
jfloat DIFFUSION_RATE,
jfloat VISCOSITY_RATE,
jfloat timestep){
jboolean isCopy;
// float * x = (*env)->GetDirectBufferAddress(env,jx);
float * x0 = (*env)->GetDirectBufferAddress(env,jx0);
// float * u = (*env)->GetDirectBufferAddress(env,ju);
// float * v = (*env)->GetDirectBufferAddress(env,jv);
// float * w = (*env)->GetDirectBufferAddress(env,jw);
float * u0 = (*env)->GetDirectBufferAddress(env,ju0);
float * v0 = (*env)->GetDirectBufferAddress(env,jv0);
float * w0 = (*env)->GetDirectBufferAddress(env,jw0);
int N = DIM_X;
int i,j,k;
vel_step(env, chunk_mask, DIM_X, ju, jv, jw, u0, v0, w0, VISCOSITY_RATE, timestep);
dens_step(env, chunk_mask, DIM_X, jx, x0, ju, jv, jw, DIFFUSION_RATE, timestep);
}
/**
* Adds values from a source array to a current frame array (eg more density to the main density array)
*/
void add_source(int N, float * x, float * s, float dt){
int i;
int size=N*N*N;
for(i=0; i<size; i++){
x[i] += dt*s[i];
}
}
/**
* Diffuses a given array by a diffusion constant
*/
void diffuse(JNIEnv * env, uint32_t chunk_mask, int N, int b, float * x, float * x0, float diff, float dt){
float a=dt*diff*N*N*N;
lin_solve(env, chunk_mask, N, b, x, x0, a, 1+6*a);
}
/**
* Advects a given array based on the force vectors in the simulation
*/
void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt){
int i, j, k, i0, j0, k0, i1, j1, k1;
float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz;
dtx=dty=dtz=dt*N;
float * d = GET_ARR(env,jrd,CENTER_LOC);
for(k=1; k<N-1; k++){
for(j=1; j<N-1; j++){
for(i=1; i<N-1; i++){
x = i-dtx*u[IX(i,j,k)]; y = j-dty*v[IX(i,j,k)]; z = k-dtz*w[IX(i,j,k)];
if (x<0.5f) x=0.5f; if (x>N+0.5f) x=N+0.5f; i0=(int)x; i1=i0+1;
if (y<0.5f) y=0.5f; if (y>N+0.5f) y=N+0.5f; j0=(int)y; j1=j0+1;
if (z<0.5f) z=0.5f; if (z>N+0.5f) z=N+0.5f; k0=(int)z; k1=k0+1;
s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; u1 = z-k0; u0 = 1-u1;
if(i0 >= N){
i0 = N - 1;
}
// if(i0 < 0){
// i0 = 0;
// }
if(j0 >= N){
j0 = N - 1;
}
// if(j0 < 0){
// j0 = 0;
// }
if(k0 >= N){
k0 = N - 1;
}
// if(k0 < 0){
// k0 = 0;
// }
if(i1 >= N){
i1 = N - 1;
}
// if(i1 < 0){
// i1 = 0;
// }
if(j1 >= N){
j1 = N - 1;
}
// if(j1 < 0){
// j1 = 0;
// }
if(k1 >= N){
k1 = N - 1;
}
// if(k1 < 0){
// k1 = 0;
// }
d[IX(i,j,k)] = s0*(t0*u0*d0[IX(i0,j0,k0)]+t1*u0*d0[IX(i0,j1,k0)]+t0*u1*d0[IX(i0,j0,k1)]+t1*u1*d0[IX(i0,j1,k1)])+
s1*(t0*u0*d0[IX(i1,j0,k0)]+t1*u0*d0[IX(i1,j1,k0)]+t0*u1*d0[IX(i1,j0,k1)]+t1*u1*d0[IX(i1,j1,k1)]);
}
}
}
set_bnd(env, chunk_mask, N, b, d);
}
/**
* The main density step function
*/
void dens_step(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jrx, float * x0, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float diff, float dt){
float * x = GET_ARR(env,jrx,CENTER_LOC);
float * u = GET_ARR(env,jru,CENTER_LOC);
float * v = GET_ARR(env,jrv,CENTER_LOC);
float * w = GET_ARR(env,jrw,CENTER_LOC);
add_source(N, x, x0, dt);
SWAP(x0, x);
diffuse(env, chunk_mask, N, 0, x, x0, diff, dt);
SWAP(x0, x);
advect(env, chunk_mask, N, 0, jrx, x0, u, v, w, dt);
}
/**
* The main velocity step function
*/
void vel_step(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float * u0, float * v0, float * w0, float visc, float dt){
float * u = GET_ARR(env,jru,CENTER_LOC);
float * v = GET_ARR(env,jrv,CENTER_LOC);
float * w = GET_ARR(env,jrw,CENTER_LOC);
add_source(N, u, u0, dt);
add_source(N, v, v0, dt);
add_source(N, w, w0, dt);
SWAP(u0, u);
diffuse(env, chunk_mask, N, 1, u, u0, visc, dt);
SWAP(v0, v);
diffuse(env, chunk_mask, N, 2, v, v0, visc, dt);
SWAP(w0, w);
diffuse(env, chunk_mask, N, 3, w, w0, visc, dt);
project(env, chunk_mask, N, jru, jrv, jrw, u0, v0);
SWAP(u0, u);
SWAP(v0, v);
SWAP(w0, w);
advect(env, chunk_mask, N, 1, jru, u0, u0, v0, w0, dt);
advect(env, chunk_mask, N, 2, jrv, v0, u0, v0, w0, dt);
advect(env, chunk_mask, N, 3, jrw, w0, u0, v0, w0, dt);
project(env, chunk_mask, N, jru, jrv, jrw, u0, v0);
}
//used for temporary vector storage when appropriate
float container[16];
/**
* Projects a given array based on force vectors
*/
void project(JNIEnv * env, uint32_t chunk_mask, int N, jobjectArray jru, jobjectArray jrv, jobjectArray jrw, float * p, float * div){
int i, j, k;
__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(env,jru,CENTER_LOC);
float * v = GET_ARR(env,jrv,CENTER_LOC);
float * w = GET_ARR(env,jrw,CENTER_LOC);
for(k=1; k<N-1; k++){
for(j=1; j<N-1; j++){
i = 1;
//
//lower
//
//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,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,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,nVector);
//multiply and finalize
vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
vector = _mm256_mul_ps(vector,constScalar);
//store
_mm256_storeu_ps(&div[IX(i,j,k)],vector);
_mm256_storeu_ps(&p[IX(i,j,k)],zeroVec);
i = 9;
//
//upper
//
//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,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,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,nVector);
//multiply and finalize
vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
vector = _mm256_mul_ps(vector,constScalar);
//store
_mm256_storeu_ps(&div[IX(i,j,k)],vector);
_mm256_storeu_ps(&p[IX(i,j,k)],zeroVec);
}
}
set_bnd(env, chunk_mask, N, 0, div);
set_bnd(env, chunk_mask, N, 0, p);
lin_solve(env, chunk_mask, N, 0, p, div, 1, 6);
constScalar = _mm256_set1_ps(0.5f*N);
for ( k=1 ; k<N-1 ; k++ ) {
for ( j=1 ; j<N-1 ; j++ ) {
//
//v
//
//lower
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,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,constScalar);
vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(9,j,k)]),vector);
_mm256_storeu_ps(&u[IX(9,j,k)],vector);
//
//v
//
//lower
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,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,constScalar);
vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(9,j,k)]),vector);
_mm256_storeu_ps(&v[IX(9,j,k)],vector);
//
//w
//
//lower
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,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,constScalar);
vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(9,j,k)]),vector);
_mm256_storeu_ps(&w[IX(9,j,k)],vector);
}
}
set_bnd(env, chunk_mask, N, 1, u);
set_bnd(env, chunk_mask, N, 2, v);
set_bnd(env, chunk_mask, N, 3, w);
}
/**
* Solves a linear system of equations in a vectorized manner
*/
void lin_solve(JNIEnv * env, uint32_t chunk_mask, int N, int b, float* x, float* x0, float a, float c){
int i, j, k, l, m;
__m256 aScalar = _mm256_set1_ps(a);
__m256 cScalar = _mm256_set1_ps(c);
// iterate the solver
for ( l=0 ; l<LINEARSOLVERTIMES ; l++ ) {
// update for each cell
for(k=1; k<N-1; k++){
for(j=1; j<N-1; j++){
int n = 0;
//solve as much as possible vectorized
for(i = 1; i < N-1; i=i+8){
__m256 vector = _mm256_loadu_ps(&x[IX(i-1,j,k)]);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i+1,j,k)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j-1,k)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j+1,k)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k-1)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k+1)]));
vector = _mm256_mul_ps(vector,aScalar);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x0[IX(i,j,k)]));
vector = _mm256_div_ps(vector,cScalar);
_mm256_storeu_ps(&x[IX(i,j,k)],vector);
}
//If there is any leftover, perform manual solving
if(i>N-1){
for(i=i-8; i < N-1; i++){
x[IX(i,j,k)] = (x0[IX(i,j,k)] + a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]+x[IX(i,j,k-1)]+x[IX(i,j,k+1)]))/c;
}
}
}
}
set_bnd(env, chunk_mask, N, b, x);
}
}
/**
* Sets the bounds of the simulation
*/
void set_bnd(JNIEnv * env, uint32_t chunk_mask, int N, int b, float * target){
int DIM = N;
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
//((x)+(DIM)*(y) + (DIM)*(DIM)*(z))
target[0 + DIM * x + DIM * DIM * y] = b==1 ? -target[1 + DIM * x + DIM * DIM * y] : target[1 + DIM * x + DIM * DIM * y];
target[IX(DIM-1,x,y)] = b==1 ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)];
target[IX(x,0,y)] = b==2 ? -target[IX(x,1,y)] : target[IX(x,1,y)];
target[IX(x,DIM-1,y)] = b==2 ? -target[IX(x,DIM-2,y)] : target[IX(x,DIM-2,y)];
target[IX(x,y,0)] = b==3 ? -target[IX(x,y,1)] : target[IX(x,y,1)];
target[IX(x,y,DIM-1)] = b==3 ? -target[IX(x,y,DIM-2)] : target[IX(x,y,DIM-2)];
}
}
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)]));
}
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);
}

28
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
#define electrosphere_FluidSim_LINEARSOLVERTIMES 10L
#undef electrosphere_FluidSim_GRAVITY
#define electrosphere_FluidSim_GRAVITY -100.0f
#define electrosphere_FluidSim_GRAVITY -1000.0f
/*
* Class: electrosphere_FluidSim
* Method: simulate
@ -84,26 +84,26 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectVectors
/*
* Class: electrosphere_FluidSim
* Method: addDensity
* Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;F)V
* Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;F)V
*/
JNIEXPORT void JNICALL Java_electrosphere_FluidSim_addDensity
(JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jfloat);
(JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jfloat);
/*
* Class: electrosphere_FluidSim
* Method: solveDiffuseDensity
* Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
* Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
*/
JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
(JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
(JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
/*
* Class: electrosphere_FluidSim
* Method: advectDensity
* Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
* Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
*/
JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectDensity
(JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
(JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
/*
* Class: electrosphere_FluidSim
@ -113,14 +113,6 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectDensity
JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setBoundsToNeighbors
(JNIEnv *, jobject, jint, jint, jint, jobjectArray);
/*
* Class: electrosphere_FluidSim
* Method: copyNeighbors
* Signature: (IIII[Ljava/nio/ByteBuffer;)V
*/
JNIEXPORT void JNICALL Java_electrosphere_FluidSim_copyNeighbors
(JNIEnv *, jobject, jint, jint, jint, jint, jobjectArray);
#ifdef __cplusplus
}
#endif

10
src/main/c/includes/libfluidsim.h
View File

@ -1,10 +0,0 @@
//include guard
#ifndef LIB_FLUID_SIM
#define LIB_FLUID_SIM
//include stb ds
#define STB_IMAGE_IMPLEMENTATION
#include "../lib/stb/stb_ds.h"
//close include guard
#endif

1
src/main/c/lib/stb

@ -1 +0,0 @@
Subproject commit ae721c50eaf761660b4f90cc590453cdb0c2acd0

584
src/main/c/src/velocitystep.c → src/main/c/velocitystep.c
View File

@ -3,9 +3,8 @@
#include <immintrin.h>
#include <stdint.h>
#include "../includes/libfluidsim.h"
#include "../includes/utilities.h"
#include "../includes/chunkmask.h"
#include "includes/utilities.h"
#include "includes/chunkmask.h"
#define BOUND_NO_DIR 0
@ -17,7 +16,7 @@
#define SET_BOUND_USE_NEIGHBOR 1
void add_source(int N, float * x, float * s, float dt);
void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, jobjectArray jrd0, float * u, float * v, float * w, float dt);
void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt);
/*
@ -67,7 +66,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveVectorDiffuse
jfloat DIFFUSION_CONST,
jfloat VISCOSITY_CONST,
jfloat dt){
float a=dt*VISCOSITY_CONST*N*N*N;
float a=dt*DIFFUSION_CONST*N*N*N;
float c=1+6*a;
int i, j, k, l, m;
float * u = GET_ARR(env,jru,CENTER_LOC);
@ -178,9 +177,7 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection
jfloat 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;
@ -192,9 +189,6 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection
float * p = GET_ARR(env,jru0,CENTER_LOC);
float * div = GET_ARR(env,jrv0,CENTER_LOC);
float scalar = 1.0/3.0;
float h = 1.0/N;
for(k=1; k<N-1; k++){
for(j=1; j<N-1; j++){
i = 1;
@ -204,15 +198,15 @@ JNIEXPORT void JNICALL 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);
@ -226,30 +220,21 @@ JNIEXPORT void JNICALL 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);
//store
_mm256_storeu_ps(&div[IX(i,j,k)],vector);
_mm256_storeu_ps(&p[IX(i,j,k)],zeroVec);
// for(i = 1; i < N - 1; i++){
// div[IX(i,j,k)] =
// -scalar*h*(u[IX(i+1,j,k)]-u[IX(i-1,j,k)]+
// v[IX(i,j+1,k)]-v[IX(i,j-1,k)]+
// w[IX(i,j,k+1)]-w[IX(i,j,k-1)]);
// p[IX(i,j,k)] = 0;
// }
}
}
}
@ -277,34 +262,31 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveProjection
__m256 aScalar = _mm256_set1_ps(a);
__m256 cScalar = _mm256_set1_ps(c);
float * p = GET_ARR(env,jru0,CENTER_LOC);
float * div = GET_ARR(env,jrv0,CENTER_LOC);
float * x = GET_ARR(env,jru0,CENTER_LOC);
float * x0 = GET_ARR(env,jrv0,CENTER_LOC);
// update for each cell
for(k=1; k<N-1; k++){
for(j=1; j<N-1; j++){
int n = 0;
//solve as much as possible vectorized
for(i = 1; i < N-1; i=i+8){
__m256 vector = _mm256_loadu_ps(&p[IX(i-1,j,k)]);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&p[IX(i+1,j,k)]));
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+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_add_ps(vector,_mm256_loadu_ps(&div[IX(i,j,k)]));
__m256 vector = _mm256_loadu_ps(&x[IX(i-1,j,k)]);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i+1,j,k)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j-1,k)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j+1,k)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k-1)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k+1)]));
vector = _mm256_mul_ps(vector,aScalar);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x0[IX(i,j,k)]));
vector = _mm256_div_ps(vector,cScalar);
_mm256_storeu_ps(&p[IX(i,j,k)],vector);
_mm256_storeu_ps(&x[IX(i,j,k)],vector);
}
//If there is any leftover, perform manual solving
if(i>N-1){
for(i=i-8; i < N-1; i++){
p[IX(i,j,k)] = (div[IX(i,j,k)] + a*(p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c;
x[IX(i,j,k)] = (x0[IX(i,j,k)] + a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]+x[IX(i,j,k-1)]+x[IX(i,j,k+1)]))/c;
}
}
// for(i=1; i < N-1; i++){
// p[IX(i,j,k)] = (div[IX(i,j,k)] + a*(p[IX(i-1,j,k)]+p[IX(i+1,j,k)]+p[IX(i,j-1,k)]+p[IX(i,j+1,k)]+p[IX(i,j,k-1)]+p[IX(i,j,k+1)]))/c;
// }
}
}
}
@ -327,10 +309,9 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
jfloat VISCOSITY_CONST,
jfloat 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 nVector = _mm256_set1_ps(N);
__m256 constScalar = _mm256_set1_ps(0.5f*N);
__m256 zeroVec = _mm256_set1_ps(0);
__m256 vector, vector2, vector3;
float * u = GET_ARR(env,jru,CENTER_LOC);
@ -340,8 +321,6 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
float * p = GET_ARR(env,jru0,CENTER_LOC);
float * div = GET_ARR(env,jrv0,CENTER_LOC);
float h = 1.0 / N;
for ( k=1 ; k<N-1 ; k++ ) {
for ( j=1 ; j<N-1 ; j++ ) {
//
@ -351,14 +330,14 @@ JNIEXPORT void JNICALL 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);
//
@ -368,14 +347,14 @@ JNIEXPORT void JNICALL 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);
//
@ -385,21 +364,16 @@ JNIEXPORT void JNICALL 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;
// }
}
}
}
@ -421,13 +395,13 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectVectors
jfloat DIFFUSION_CONST,
jfloat VISCOSITY_CONST,
jfloat dt){
advect(env,chunk_mask,N,1,jru,jru0,GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
advect(env,chunk_mask,N,2,jrv,jrv0,GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
advect(env,chunk_mask,N,3,jrw,jrw0,GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
advect(env,chunk_mask,N,1,jru,GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
advect(env,chunk_mask,N,2,jrv,GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
advect(env,chunk_mask,N,3,jrw,GET_ARR(env,jrw0,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
}
void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, jobjectArray jrd0, float * u, float * v, float * w, float dt){
void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt){
int i, j, k, i0, j0, k0, i1, j1, k1;
int m,n,o;
float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz;
@ -436,12 +410,12 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, j
float * d = GET_ARR(env,jrd,CENTER_LOC);
float * d0 = GET_ARR(env,jrd0,CENTER_LOC);
float * sampleArr = d0;
for(k=1; k<N-1; k++){
for(j=1; j<N-1; j++){
for(i=1; i<N-1; i++){
d0 = GET_ARR(env,jrd0,CENTER_LOC);
sampleArr = d0;
//calculate location to pull from
x = i-dtx*u[IX(i,j,k)];
y = j-dty*v[IX(i,j,k)];
@ -449,175 +423,55 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, j
m = n = o = 1;
if(x < 0){ m += 1; }
else if(x >= N){ m -= 1; }
if(y < 0){ n += 1; }
else if(y >= N){ n -= 1; }
if(z < 0){ o += 1; }
else if(z >= N){ o -= 1; }
if(x < 1){ m -= 1; }
if(x >= N-1){ m += 1; }
if(y < 1){ n -= 1; }
if(y >= N-1){ n += 1; }
if(z < 1){ o -= 1; }
if(z >= N-1){ o += 1; }
//If the out of bounds coordinate is in bounds for a neighbor chunk, use that chunk as source instead
if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
// if(i == 1 && j == 1 && k == 1){
// printf("\narr indices: %d %d %d\n\n",m,n,o);
// }
//cases:
//if x = 17.01, m = 2
// 17 in current array is 1 in neighbor
// 18 in current array is 2 in neighbor
// 19 in current array is 3 in neighbor
//want to sample neighbor array at 1 & 2
//x becomes 1.01, sampling new array (keep in mind that 0 in the new array should contain the current array values)
//modification: subtract 16
//cases:
//if x = 16.99, m = 2
// 16 in current array is 0 in neighbor
// 17 in current array is 1 in neighbor
// 18 in current array is 2 in neighbor
// 19 in current array is 3 in neighbor
//want to sample current array still
//x becomes 1.01, sampling new array (keep in mind that 0 in the new array should contain the current array values)
//modification: no modification
//if x = 0.01, m = 0
// 0 in current array is 16 in neighbor
//-1 in current array is 15 in neighbor
//-2 in current array is 14 in neighbor
//want to sample current array still
//x becomes 15.01, sampling new array (keep in mind that 17 in the new array should contain the current array)
//modification: no modification
//if x = -0.01, m = 0
// 0 in current array is 16 in neighbor
//-1 in current array is 15 in neighbor
//-2 in current array is 14 in neighbor
//want to sample -1 & 0, so i0 becomes 15
//x becomes 15.99, sampling new array (keep in mind that 17 in the new array should contain the current array)
//modification: add 16
//if x = -2, m = 0
// 0 in current array is 16 in neighbor
//-1 in current array is 15 in neighbor
//-2 in current array is 14 in neighbor
//x becomes 14, sampling new array (keep in mind that 17 in the new array should contain the current array)
//modification: add 16
// printf("Hit other chunk\n");
d0 = GET_ARR(env,jrd0,CK(m,n,o));
x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * (N-2);
// printf("%d => %f\n",m,x);
y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * (N-2);
z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * (N-2);
}
// if(CK(m,n,o) != CENTER_LOC){
// printf("Looking in border chunk\n");
// }
// if(x > 16){
// printf("%f %d %d %d\n",m,n,o);
// }
// if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
// // printf("Hit other chunk\n");
// sampleArr = GET_ARR(env,jrd,CK(m,n,o));
// x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * N;
// y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * N;
// z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * N;
// }
//clamp location within chunk
//get indices, and calculate percentage to pull from each index
if(x < 0.001f){
//cases to consider:
//m = 0, x = -10
//m = 2, x = 0.01
x=0.001f;
i0=(int)0;
i1=1;
s0 = 0.999f;
s1 = 0.001f;
} else if(x > N - 1){
//cases to consider:
//m = 0, x = 17.01
//m = 2, x = 20
x = N-1;
i0=(int)N-2;
i1=N-1;
s0 = 0.001f;
s1 = 0.999f;
} else {
i0=(int)x;
i1=i0+1;
s1 = x-i0;
s0 = 1-s1;
}
if(y < 0.001f){
//cases to consider:
//m = 0, x = -10
//m = 2, x = 0.01
y=0.001f;
j0=(int)0;
j1=1;
t0 = 0.999f;
t1 = 0.001f;
} else if(y > N - 1){
//cases to consider:
//m = 0, x = 17.01
//m = 2, x = 20
y = N-1;
j0=(int)N-2;
j1=N-1;
t0 = 0.001f;
t1 = 0.999f;
} else {
j0=(int)y;
j1=j0+1;
t1 = y-j0;
t0 = 1-t1;
}
if(z < 0.001f){
//cases to consider:
//m = 0, x = -10
//m = 2, x = 0.01
z=0.001f;
k0=(int)0;
k1=1;
u0 = 0.999f;
u1 = 0.001f;
} else if(z > N - 1){
//cases to consider:
//m = 0, x = 17.01
//m = 2, x = 20
z = N-1;
k0=(int)N-2;
k1=N-1;
u0 = 0.001f;
u1 = 0.999f;
} else {
k0=(int)z;
k1=k0+1;
u1 = z-k0;
u0 = 1-u1;
}
// if (x<0.001f) x=0.001f;
// if (x>N+0.5f) x=N+0.5f;
// if (y<0.001f) y=0.001f;
// if (y>N+0.5f) y=N+0.5f;
// if (z<0.001f) z=0.001f;
// if (z>N+0.5f) z=N+0.5f;
if (x<0.5f) x=0.5f;
if (x>N+0.5f) x=N+0.5f;
if (y<0.5f) y=0.5f;
if (y>N+0.5f) y=N+0.5f;
if (z<0.5f) z=0.5f;
if (z>N+0.5f) z=N+0.5f;
//get actual indices
// i0=(int)x;
// i1=i0+1;
// j0=(int)y;
// j1=j0+1;
// k0=(int)z;
// k1=k0+1;
i0=(int)x;
i1=i0+1;
j0=(int)y;
j1=j0+1;
k0=(int)z;
k1=k0+1;
//calculate percentage of each index
// s1 = x-i0;
// s0 = 1-s1;
// t1 = y-j0;
// t0 = 1-t1;
// u1 = z-k0;
// u0 = 1-u1;
s1 = x-i0;
s0 = 1-s1;
t1 = y-j0;
t0 = 1-t1;
u1 = z-k0;
u0 = 1-u1;
// if(i0 >= N){
// i0 = N - 1;
// }
if(i0 >= N){
i0 = N - 1;
}
// if(i0 < 0){
// i0 = 0;
// }
@ -633,9 +487,9 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, j
// if(k0 < 0){
// k0 = 0;
// }
// if(i1 >= N){
// i1 = N - 1;
// }
if(i1 >= N){
i1 = N - 1;
}
// if(i1 < 0){
// i1 = 0;
// }
@ -653,16 +507,16 @@ void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, j
// }
d[IX(i,j,k)] =
s0*(
t0*u0*d0[IX(i0,j0,k0)]+
t1*u0*d0[IX(i0,j1,k0)]+
t0*u1*d0[IX(i0,j0,k1)]+
t1*u1*d0[IX(i0,j1,k1)]
t0*u0*sampleArr[IX(i0,j0,k0)]+
t1*u0*sampleArr[IX(i0,j1,k0)]+
t0*u1*sampleArr[IX(i0,j0,k1)]+
t1*u1*sampleArr[IX(i0,j1,k1)]
)+
s1*(
t0*u0*d0[IX(i1,j0,k0)]+
t1*u0*d0[IX(i1,j1,k0)]+
t0*u1*d0[IX(i1,j0,k1)]+
t1*u1*d0[IX(i1,j1,k1)]
t0*u0*sampleArr[IX(i1,j0,k0)]+
t1*u0*sampleArr[IX(i1,j1,k0)]+
t0*u1*sampleArr[IX(i1,j0,k1)]+
t1*u1*sampleArr[IX(i1,j1,k1)]
);
}
}
@ -679,16 +533,35 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setBoundsToNeighbors
int DIM = N;
float * target = GET_ARR(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)];
// if(ARR_EXISTS(chunk_mask,0,1,1)){
// source = GET_ARR(env,neighborArray,CK(0,1,1));
// for(int x=1; x < DIM-1; x++){
// for(int y = 1; y < DIM-1; y++){
// target[IX(0,x,y)] = source[IX(DIM-2,x,y)];
// }
// }
// } else {
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)];
// }
// if(ARR_EXISTS(chunk_mask,2,1,1)){
// source = GET_ARR(env,neighborArray,CK(2,1,1));
// for(int x=1; x < DIM-1; x++){
// for(int y = 1; y < DIM-1; y++){
// target[IX(DIM-1,x,y)] = source[IX(1,x,y)];
// }
// }
// } else {
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++){
//((x)+(DIM)*(y) + (DIM)*(DIM)*(z))
@ -725,228 +598,3 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setBoundsToNeighbors
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);
}
/**
* This exclusively copies neighbors to make sure zeroing out stuff doesn't break sim
*/
JNIEXPORT void JNICALL Java_electrosphere_FluidSim_copyNeighbors
(JNIEnv * env,
jobject this,
jint N,
jint chunk_mask,
jint cx,
jint vector_dir,
jobjectArray neighborArray){
int DIM = N;
float * target = GET_ARR(env,neighborArray,CENTER_LOC);
float * source;
//
//
// PLANES
//
//
if(ARR_EXISTS(chunk_mask,0,1,1)){
source = GET_ARR(env,neighborArray,CK(0,1,1));
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
target[IX(0,x,y)] = source[IX(DIM-2,x,y)];
}
}
}
if(ARR_EXISTS(chunk_mask,2,1,1)){
source = GET_ARR(env,neighborArray,CK(2,1,1));
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
target[IX(DIM-1,x,y)] = source[IX(1,x,y)];
}
}
}
if(ARR_EXISTS(chunk_mask,1,0,1)){
source = GET_ARR(env,neighborArray,CK(1,0,1));
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
target[IX(x,0,y)] = source[IX(x,DIM-2,y)];
}
}
}
if(ARR_EXISTS(chunk_mask,1,2,1)){
source = GET_ARR(env,neighborArray,CK(1,2,1));
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
target[IX(x,DIM-1,y)] = source[IX(x,1,y)];
}
}
}
if(ARR_EXISTS(chunk_mask,1,1,0)){
source = GET_ARR(env,neighborArray,CK(1,1,0));
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
target[IX(x,y,0)] = source[IX(x,y,DIM-2)];
}
}
}
if(ARR_EXISTS(chunk_mask,1,1,2)){
source = GET_ARR(env,neighborArray,CK(1,1,2));
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){
target[IX(x,y,DIM-1)] = source[IX(x,y,1)];
}
}
}
//
//
// EDGES
//
//
if(ARR_EXISTS(chunk_mask,0,0,1)){
source = GET_ARR(env,neighborArray,CK(0,0,1));
for(int x=1; x < DIM-1; x++){
target[IX(0,0,x)] = source[IX(DIM-2,DIM-2,x)];
}
}
if(ARR_EXISTS(chunk_mask,2,0,1)){
source = GET_ARR(env,neighborArray,CK(2,0,1));
for(int x=1; x < DIM-1; x++){
target[IX(DIM-1,0,x)] = source[IX(1,DIM-2,x)];
}
}
if(ARR_EXISTS(chunk_mask,0,2,1)){
source = GET_ARR(env,neighborArray,CK(0,2,1));
for(int x=1; x < DIM-1; x++){
target[IX(0,DIM-1,x)] = source[IX(DIM-2,1,x)];
}
}
if(ARR_EXISTS(chunk_mask,2,2,1)){
source = GET_ARR(env,neighborArray,CK(2,2,1));
for(int x=1; x < DIM-1; x++){
target[IX(DIM-1,DIM-1,x)] = source[IX(1,1,x)];
}
}
//
//
if(ARR_EXISTS(chunk_mask,0,1,0)){
source = GET_ARR(env,neighborArray,CK(0,1,0));
for(int x=1; x < DIM-1; x++){
target[IX(0,x,0)] = source[IX(DIM-2,x,DIM-2)];
}
}
if(ARR_EXISTS(chunk_mask,2,1,0)){
source = GET_ARR(env,neighborArray,CK(2,1,0));
for(int x=1; x < DIM-1; x++){
target[IX(DIM-1,x,0)] = source[IX(1,x,DIM-2)];
}
}
if(ARR_EXISTS(chunk_mask,0,1,2)){
source = GET_ARR(env,neighborArray,CK(0,1,2));
for(int x=1; x < DIM-1; x++){
target[IX(0,x,DIM-1)] = source[IX(DIM-2,x,1)];
}
}
if(ARR_EXISTS(chunk_mask,2,1,2)){
source = GET_ARR(env,neighborArray,CK(2,1,2));
for(int x=1; x < DIM-1; x++){
target[IX(DIM-1,x,DIM-1)] = source[IX(1,x,1)];
}
}
//
//
if(ARR_EXISTS(chunk_mask,1,0,0)){
source = GET_ARR(env,neighborArray,CK(1,0,0));
for(int x=1; x < DIM-1; x++){
target[IX(x,0,0)] = source[IX(x,DIM-2,DIM-2)];
}
}
if(ARR_EXISTS(chunk_mask,1,2,0)){
source = GET_ARR(env,neighborArray,CK(1,2,0));
for(int x=1; x < DIM-1; x++){
target[IX(x,DIM-1,0)] = source[IX(x,1,DIM-2)];
}
}
if(ARR_EXISTS(chunk_mask,1,0,2)){
source = GET_ARR(env,neighborArray,CK(1,0,2));
for(int x=1; x < DIM-1; x++){
target[IX(x,0,DIM-1)] = source[IX(x,DIM-2,1)];
}
}
if(ARR_EXISTS(chunk_mask,1,2,2)){
source = GET_ARR(env,neighborArray,CK(1,2,2));
for(int x=1; x < DIM-1; x++){
target[IX(x,DIM-1,DIM-1)] = source[IX(x,1,1)];
}
}
//
//
// CORNERS
//
//
if(ARR_EXISTS(chunk_mask,0,0,0)){
source = GET_ARR(env,neighborArray,CK(0,0,0));
target[IX(0,0,0)] = source[IX(DIM-2,DIM-2,DIM-2)];
}
if(ARR_EXISTS(chunk_mask,2,0,0)){
source = GET_ARR(env,neighborArray,CK(2,0,0));
target[IX(DIM-1,0,0)] = source[IX(1,DIM-2,DIM-2)];
}
if(ARR_EXISTS(chunk_mask,0,2,0)){
source = GET_ARR(env,neighborArray,CK(0,2,0));
target[IX(0,DIM-1,0)] = source[IX(DIM-2,1,DIM-2)];
}
if(ARR_EXISTS(chunk_mask,2,2,0)){
source = GET_ARR(env,neighborArray,CK(2,2,0));
target[IX(DIM-1,DIM-1,0)] = source[IX(1,1,DIM-2)];
}
//
//
if(ARR_EXISTS(chunk_mask,0,0,2)){
source = GET_ARR(env,neighborArray,CK(0,0,2));
target[IX(0,0,DIM-1)] = source[IX(DIM-2,DIM-2,1)];
}
if(ARR_EXISTS(chunk_mask,2,0,2)){
source = GET_ARR(env,neighborArray,CK(2,0,2));
target[IX(DIM-1,0,DIM-1)] = source[IX(1,DIM-2,1)];
}
if(ARR_EXISTS(chunk_mask,0,2,2)){
source = GET_ARR(env,neighborArray,CK(0,2,2));
target[IX(0,DIM-1,DIM-1)] = source[IX(DIM-2,1,1)];
}
if(ARR_EXISTS(chunk_mask,2,2,2)){
source = GET_ARR(env,neighborArray,CK(2,2,2));
target[IX(DIM-1,DIM-1,DIM-1)] = source[IX(1,1,1)];
}
}

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

@ -42,7 +42,7 @@ public class FluidSim {
//Buffers that contain density for current frame
ByteBuffer[] density = new ByteBuffer[27];
//Buffers that contain new density to add to the simulation
ByteBuffer[] densityAddition = new ByteBuffer[27];
ByteBuffer densityAddition;
//Buffers that contain u vector directions
ByteBuffer[] uVector = new ByteBuffer[27];
//Buffers that contain v vector directions
@ -66,23 +66,23 @@ public class FluidSim {
float[] wArrayView = new float[DIM * DIM * DIM];
//these should be set to the
float[] u0ArrayView = new float[DIM * DIM * DIM];
public float[] v0ArrayView = new float[DIM * DIM * DIM];
float[] v0ArrayView = new float[DIM * DIM * DIM];
float[] w0ArrayView = new float[DIM * DIM * DIM];
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 int LINEARSOLVERTIMES = 10;
static final float GRAVITY = -100f;
static final float GRAVITY = -1000f;
public void setup(Vector3i offset){
//allocate buffers for this chunk
density[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
densityAddition[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
densityAddition = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
uVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
vVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
wVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
@ -91,7 +91,7 @@ public class FluidSim {
wAdditionVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
//order endian-ness
density[13].order(ByteOrder.LITTLE_ENDIAN);
densityAddition[13].order(ByteOrder.LITTLE_ENDIAN);
densityAddition.order(ByteOrder.LITTLE_ENDIAN);
uVector[13].order(ByteOrder.LITTLE_ENDIAN);
vVector[13].order(ByteOrder.LITTLE_ENDIAN);
wVector[13].order(ByteOrder.LITTLE_ENDIAN);
@ -110,19 +110,19 @@ 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 == 1){
if(
Math.abs(16 - i) < 5 &&
Math.abs(j) < 5 &&
Math.abs(16 - k) < 5 &&
Math.abs(5 - i) < 4 &&
Math.abs(j) < 4 &&
Math.abs(5 - 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);
uf.put(1);
vf.put(-1f);
wf.put(rand.nextFloat() * 0.1f);
} else {
xf.put(0);
uf.put(0);
@ -131,17 +131,17 @@ public class FluidSim {
}
} else {
if(
Math.abs(0 - i) < 5 &&
Math.abs(j) < 5 &&
Math.abs(0 - k) < 5 &&
Math.abs(8 - i) < 4 &&
Math.abs(j) < 4 &&
Math.abs(16 - 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(rand.nextFloat() * 0.1f);
xf.put(1);
uf.put(1);
vf.put(-1f);
wf.put(rand.nextFloat() * 0.1f);
} else {
xf.put(0);
uf.put(0);
@ -178,23 +178,20 @@ public class FluidSim {
addVectorSources(simArray, timestep);
swapAllVectorFields(simArray, timestep);
solveVectorDiffusion(simArray, timestep);
solveProjection(simArray, step, timestep);
solveProjection(simArray, timestep);
swapAllVectorFields(simArray, timestep);
advectVectorsAcrossBoundaries(simArray, timestep);
solveProjection(simArray, step, timestep);
solveProjection(simArray, timestep);
//
//Density stage
addDensity(simArray, timestep);
swapAllDensityArrays(simArray, timestep);
diffuseDensity(simArray, timestep);
swapAllDensityArrays(simArray, timestep);
advectDensity(simArray, timestep);
// mirrorNeighborDensities(simArray, timestep);
//
//Read out of buffers back into accessible arrays
for(int x = 0; x < simArray.length; x++){
@ -209,128 +206,25 @@ public class FluidSim {
}
}
private static double sumAllDensity(FluidSim[][][] simArray){
double rVal = 0;
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++){
rVal = rVal + simArray[x][y][z].sumDensity();
}
}
}
return rVal;
}
private double sumDensity(){
double rVal = 0;
for(int x = 1; x < DIM - 1; x++){
for(int y = 1; y < DIM - 1; y++){
for(int z = 1; z < DIM - 1; z++){
rVal = rVal + densityArrayView[IX(x,y,z)];
}
}
}
return rVal;
}
private double sumU(){
double rVal = 0;
for(int x = 1; x < DIM - 1; x++){
for(int y = 1; y < DIM - 1; y++){
for(int z = 1; z < DIM - 1; z++){
rVal = rVal + Math.abs(uArrayView[IX(x,y,z)]);
}
}
}
return rVal;
}
private static double sumAllU(FluidSim[][][] simArray){
double rVal = 0;
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++){
simArray[x][y][z].readDataIntoArrays();
rVal = rVal + simArray[x][y][z].sumU();
}
}
}
return rVal;
}
private double sumU0(){
double rVal = 0;
for(int x = 1; x < DIM - 1; x++){
for(int y = 1; y < DIM - 1; y++){
for(int z = 1; z < DIM - 1; z++){
rVal = rVal + Math.abs(u0ArrayView[IX(x,y,z)]);
}
}
}
return rVal;
}
private static double sumAllU0(FluidSim[][][] simArray){
double rVal = 0;
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++){
simArray[x][y][z].readDataIntoArrays();
rVal = rVal + simArray[x][y][z].sumU0();
}
}
}
return rVal;
}
private double sumV(){
double rVal = 0;
for(int x = 1; x < DIM - 1; x++){
for(int y = 1; y < DIM - 1; y++){
for(int z = 1; z < DIM - 1; z++){
rVal = rVal + Math.abs(vArrayView[IX(x,y,z)]);
}
}
}
return rVal;
}
private static double sumAllV(FluidSim[][][] simArray){
double rVal = 0;
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++){
simArray[x][y][z].readDataIntoArrays();
rVal = rVal + simArray[x][y][z].sumV();
}
}
}
return rVal;
}
private double sumV0(){
double rVal = 0;
for(int x = 1; x < DIM - 1; x++){
for(int y = 1; y < DIM - 1; y++){
for(int z = 1; z < DIM - 1; z++){
rVal = rVal + Math.abs(v0ArrayView[IX(x,y,z)]);
}
}
}
return rVal;
}
private static double sumAllV0(FluidSim[][][] simArray){
double rVal = 0;
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++){
simArray[x][y][z].readDataIntoArrays();
rVal = rVal + simArray[x][y][z].sumV0();
}
}
}
return rVal;
/**
* Runs a frame of the fluid simulation
*/
private void simulate(int step, float timestep){
simulate(
DIM,
0,
density,
densityAddition,
uVector,
vVector,
wVector,
uAdditionVector,
vAdditionVector,
wAdditionVector,
DIFFUSION_CONSTANT,
VISCOSITY_CONSTANT,
timestep
);
}
private static void solveChunkMask(FluidSim[][][] simArray){
@ -362,20 +256,6 @@ public class FluidSim {
}
private static void solveVectorDiffusion(FluidSim[][][] simArray, float timestep){
//samples u,v,w,u0,v0,w0
//sets u,v,w
// 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++){
// simArray[x][y][z].copyNeighborsWrapper(1, simArray[x][y][z].uVector);
// simArray[x][y][z].copyNeighborsWrapper(2, simArray[x][y][z].vVector);
// simArray[x][y][z].copyNeighborsWrapper(3, simArray[x][y][z].wVector);
// simArray[x][y][z].copyNeighborsWrapper(1, simArray[x][y][z].uAdditionVector);
// simArray[x][y][z].copyNeighborsWrapper(2, simArray[x][y][z].vAdditionVector);
// simArray[x][y][z].copyNeighborsWrapper(3, simArray[x][y][z].wAdditionVector);
// }
// }
// }
for(int l = 0; l < LINEARSOLVERTIMES; l++){
for(int x = 0; x < simArray.length; x++){
for(int y = 0; y < simArray[0].length; y++){
@ -394,38 +274,16 @@ public class FluidSim {
simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wVector);
simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
}
}
}
}
}
private static void solveProjection(FluidSim[][][] simArray, int step, float timestep){
//samples u,v,w
//sets u0,v0
private static void solveProjection(FluidSim[][][] simArray, float timestep){
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++){
simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vAdditionVector);
}
}
}
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++){
// System.out.println("Setup " + x + " " + y + " " + z);
//setup projection across boundaries
//...
//set boundaries appropriately
@ -439,14 +297,9 @@ public class FluidSim {
for(int z = 0; z < simArray[0][0].length; z++){
simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].vAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].vAdditionVector);
}
}
}
//samples u0, v0
//sets u0
//these should have just been mirrored in the above
//
//Perform main projection solver
for(int l = 0; l < LINEARSOLVERTIMES; l++){
@ -464,13 +317,11 @@ public class FluidSim {
for(int y = 0; y < simArray[0].length; y++){
for(int z = 0; z < simArray[0][0].length; z++){
simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].uAdditionVector);
// simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].vAdditionVector);
}
}
}
}
//samples u,v,w,u0
//sets u,v,w
//Finalize projection
for(int x = 0; x < simArray.length; x++){
for(int y = 0; y < simArray[0].length; y++){
@ -489,15 +340,6 @@ public class FluidSim {
simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wVector);
simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vAdditionVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].vAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].wAdditionVector);
}
}
}
@ -521,19 +363,6 @@ public class FluidSim {
}
}
}
//then need to mirror each array as relevant
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++){
simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].wAdditionVector);
}
}
}
}
private void swapVectorFields(){
@ -543,42 +372,21 @@ public class FluidSim {
// densityAddition = density[13];
// density[13] = tmp;
//swap u0 <-> u
for(int i = 0; i < 27; i++){
tmp = uAdditionVector[i];
uAdditionVector[i] = uVector[i];
uVector[i] = tmp;
//swap v0 <-> v
tmp = vAdditionVector[i];
vAdditionVector[i] = vVector[i];
vVector[i] = tmp;
//swap w0 <-> w
tmp = wAdditionVector[i];
wAdditionVector[i] = wVector[i];
wVector[i] = tmp;
}
tmp = uAdditionVector[13];
uAdditionVector[13] = uVector[13];
uVector[13] = tmp;
//swap v0 <-> v
tmp = vAdditionVector[13];
vAdditionVector[13] = vVector[13];
vVector[13] = tmp;
//swap w0 <-> w
tmp = wAdditionVector[13];
wAdditionVector[13] = wVector[13];
wVector[13] = tmp;
//...
}
private static void advectVectorsAcrossBoundaries(FluidSim[][][] simArray, float timestep){
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++){
simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vAdditionVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vAdditionVector);
simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wAdditionVector);
}
}
}
//samples u,v,w,u0,v0,w0
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++){
@ -597,9 +405,6 @@ public class FluidSim {
simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vVector);
simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wVector);
simArray[x][y][z].copyNeighborsWrapper(1, x, simArray[x][y][z].uVector);
simArray[x][y][z].copyNeighborsWrapper(2, x, simArray[x][y][z].vVector);
simArray[x][y][z].copyNeighborsWrapper(3, x, simArray[x][y][z].wVector);
}
}
}
@ -613,36 +418,16 @@ public class FluidSim {
simArray[x][y][z].addDensityWrapper(timestep);
//swap x <=> x0
//swap arrays in java side...
// simArray[x][y][z].swapDensityArrays();
}
}
}
}
private static void swapAllDensityArrays(FluidSim[][][] simArray, float timestep){
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++){
simArray[x][y][z].swapDensityArrays();
}
}
}
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++){
simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].density);
simArray[x][y][z].copyNeighborsWrapper(0, x, simArray[x][y][z].densityAddition);
}
}
}
}
private void swapDensityArrays(){
for(int i = 0; i < 27; i++){
ByteBuffer tmp = density[i];
density[i] = densityAddition[i];
densityAddition[i] = tmp;
}
ByteBuffer tmp = density[13];
density[13] = densityAddition;
densityAddition = tmp;
}
private static void diffuseDensity(FluidSim[][][] simArray, float timestep){
@ -672,7 +457,7 @@ public class FluidSim {
for(int y = 0; y < simArray[0].length; y++){
for(int z = 0; z < simArray[0][0].length; z++){
//swap x <=> x0 again
// simArray[x][y][z].swapDensityArrays();
simArray[x][y][z].swapDensityArrays();
//advect density
simArray[x][y][z].advectDensityWrapper(timestep);
}
@ -790,6 +575,14 @@ public class FluidSim {
}
private native void solveProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
/**
* Solve projection system
*/
// private void setProjectionBordersWrapper(float timestep){
// setProjectionBorders(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
// }
// private native void setProjectionBorders(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
/**
* Does work like subtracting curl from vector field, setting boundaries, etc
*/
@ -812,7 +605,7 @@ public class FluidSim {
private void addDensityWrapper(float timestep){
addDensity(DIM, chunkMask, density, densityAddition, timestep);
}
private native void addDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer[] x0, float timestep);
private native void addDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer x0, float timestep);
/**
* Solve density diffusion
@ -820,7 +613,7 @@ public class FluidSim {
private void solveDiffuseDensityWrapper(float timestep){
solveDiffuseDensity(DIM, chunkMask, density, densityAddition, uVector, vVector, wVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
}
private native void solveDiffuseDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer[] x0, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
private native void solveDiffuseDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer x0, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
/**
* Solve density diffusion
@ -828,7 +621,7 @@ public class FluidSim {
private void advectDensityWrapper(float timestep){
advectDensity(DIM, chunkMask, density, densityAddition, uVector, vVector, wVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
}
private native void advectDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer[] x0, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
private native void advectDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer x0, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
/**
@ -839,16 +632,6 @@ public class FluidSim {
}
private native void setBoundsToNeighbors(int DIM_X, int chunkMask, int vectorDir, ByteBuffer[] neighborMap);
/**
* Sets the bounds of the neighbormap to neighbor values if available, otherwise doesn't mess with them.
* This is to make sure zeroing out doesn't mess up the sim
*/
private void copyNeighborsWrapper(int vectorDir, int x, ByteBuffer[] neighborMap){
copyNeighbors(DIM, chunkMask, x, vectorDir, neighborMap);
}
private native void copyNeighbors(int DIM_X, int chunkMask, int x, int vectorDir, ByteBuffer[] neighborMap);
@ -920,34 +703,17 @@ public class FluidSim {
if(wVector[13].position() > 0){
wVector[13].position(0);
}
if(uAdditionVector[13].position() > 0){
uAdditionVector[13].position(0);
}
if(vAdditionVector[13].position() > 0){
vAdditionVector[13].position(0);
}
if(wAdditionVector[13].position() > 0){
wAdditionVector[13].position(0);
}
FloatBuffer xFloatView = density[13].asFloatBuffer();
FloatBuffer uFloatView = uVector[13].asFloatBuffer();
FloatBuffer vFloatView = vVector[13].asFloatBuffer();
FloatBuffer wFloatView = wVector[13].asFloatBuffer();
FloatBuffer u0FloatView = uAdditionVector[13].asFloatBuffer();
FloatBuffer v0FloatView = vAdditionVector[13].asFloatBuffer();
FloatBuffer w0FloatView = wAdditionVector[13].asFloatBuffer();
int index = 0;
for(int i = 0; i < DIM; i++){
for(int j = 0; j < DIM; j++){
for(int k = 0; k < DIM; k++){
index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
densityArrayView[index] = xFloatView.get();
uArrayView[index] = uFloatView.get();
vArrayView[index] = vFloatView.get();
wArrayView[index] = wFloatView.get();
u0ArrayView[index] = u0FloatView.get();
v0ArrayView[index] = v0FloatView.get();
w0ArrayView[index] = w0FloatView.get();
densityArrayView[IX(i,j,k)] = xFloatView.get();
uArrayView[IX(i,j,k)] = uFloatView.get();
vArrayView[IX(i,j,k)] = vFloatView.get();
wArrayView[IX(i,j,k)] = wFloatView.get();
}
}
}
@ -957,8 +723,8 @@ public class FluidSim {
* Writes data from the java-side arrays into buffers that get passed into c-side
*/
private void writeNewStateIntoBuffers(){
if(densityAddition[13].position() > 0){
densityAddition[13].position(0);
if(densityAddition.position() > 0){
densityAddition.position(0);
}
if(uAdditionVector[13].position() > 0){
uAdditionVector[13].position(0);
@ -969,19 +735,17 @@ public class FluidSim {
if(wAdditionVector[13].position() > 0){
wAdditionVector[13].position(0);
}
FloatBuffer x0FloatView = densityAddition[13].asFloatBuffer();
FloatBuffer x0FloatView = densityAddition.asFloatBuffer();
FloatBuffer u0FloatView = uAdditionVector[13].asFloatBuffer();
FloatBuffer v0FloatView = vAdditionVector[13].asFloatBuffer();
FloatBuffer w0FloatView = wAdditionVector[13].asFloatBuffer();
int index = 0;
for(int i = 0; i < DIM; i++){
for(int j = 0; j < DIM; j++){
for(int k = 0; k < DIM; k++){
index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
x0FloatView.put(density0ArrayView[index]);
u0FloatView.put(u0ArrayView[index]);
v0FloatView.put(v0ArrayView[index]);
w0FloatView.put(w0ArrayView[index]);
x0FloatView.put(density0ArrayView[IX(i,j,k)]);
u0FloatView.put(u0ArrayView[IX(i,j,k)]);
v0FloatView.put(v0ArrayView[IX(i,j,k)]);
w0FloatView.put(w0ArrayView[IX(i,j,k)]);
}
}
}
@ -991,14 +755,10 @@ public class FluidSim {
* Adds gravity to the simulation
*/
private void addGravity(){
int index = 0;
for(int i = 0; i < DIM; i++){
for(int j = 0; j < DIM; j++){
for(int k = 0; k < DIM; k++){
index = ((i)+(DIM)*(j) + (DIM)*(DIM)*(k));
u0ArrayView[index] = 0;
v0ArrayView[index] = densityArrayView[index] * GRAVITY;
w0ArrayView[index] = 0;
v0ArrayView[IX(i,j,k)] = densityArrayView[IX(i,j,k)] * GRAVITY;
}
}
}
@ -1063,10 +823,6 @@ public class FluidSim {
return density[getNeighborIndex(1,1,1)];
}
public ByteBuffer getDensityAdditionBuffer(){
return densityAddition[getNeighborIndex(1,1,1)];
}
public ByteBuffer getUBuffer(){
return uVector[getNeighborIndex(1,1,1)];
}
@ -1093,7 +849,6 @@ public class FluidSim {
public void setNeighbor(int x, int y, int z, FluidSim neighbor){
density[getNeighborIndex(x,y,z)] = neighbor.getDensityBuffer();
densityAddition[getNeighborIndex(x,y,z)] = neighbor.getDensityAdditionBuffer();
uVector[getNeighborIndex(x,y,z)] = neighbor.getUBuffer();
vVector[getNeighborIndex(x,y,z)] = neighbor.getVBuffer();
wVector[getNeighborIndex(x,y,z)] = neighbor.getWBuffer();

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

@ -22,7 +22,7 @@ public class Main {
public static void main(String args[]){
int dim = 10;
int dim = 2;
int i = 0;
long time = 0;
long lastTime = 0;
@ -55,7 +55,7 @@ public class Main {
//
//Simulate
//
FluidSim.simChunks(simArray,i,0.01f);
FluidSim.simChunks(simArray,i,0.001f);
time = time + (System.currentTimeMillis() - lastTime);
//
//Remesh
@ -70,8 +70,8 @@ public class Main {
//redraw
GLFWContext.redraw(meshArray);
i++;
if(i == 100){
System.out.println(time / 100.0);
if(i == 1000){
System.out.println(time / 1000.0);
}
if(i > 3){
// scan.next();

6
src/main/java/electrosphere/render/Mesh.java
View File

@ -2,13 +2,11 @@ package electrosphere.render;
import java.io.BufferedInputStream;
import java.io.BufferedReader;
import java.io.File;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.nio.FloatBuffer;
import java.nio.IntBuffer;
import java.nio.file.Files;
import java.util.HashMap;
import java.util.LinkedList;
import java.util.List;
@ -243,7 +241,7 @@ public class Mesh {
public static void initShaderProgram(){
String vsSrc = "";
ClassLoader classloader = Thread.currentThread().getContextClassLoader();
try (BufferedReader is = new BufferedReader(Files.newBufferedReader(new File("C:\\Users\\satellite\\Documents\\fluid-sim\\src\\main\\resources\\shader.vs").toPath()))){
try (BufferedReader is = new BufferedReader(new InputStreamReader(classloader.getResourceAsStream("shader.vs")))){
String temp;
while((temp = is.readLine())!=null){
vsSrc = vsSrc + temp + "\n";
@ -254,7 +252,7 @@ public class Mesh {
}
String fsSrc = "";
try (BufferedReader is = new BufferedReader(Files.newBufferedReader(new File("C:\\Users\\satellite\\Documents\\fluid-sim\\src\\main\\resources\\shader.fs").toPath()))){
try (BufferedReader is = new BufferedReader(new InputStreamReader(classloader.getResourceAsStream("shader.fs")))){
String temp;
while((temp = is.readLine())!=null){
fsSrc = fsSrc + temp + "\n";