Renderer/src/main/c/src/fluid/sim/grid2/fluidsim.c
austin 1b36112e24
Some checks failed
studiorailgun/Renderer/pipeline/head There was a failure building this commit
refactoring grid2 headers
2024-12-09 18:01:27 -05:00

455 lines
19 KiB
C

#include <stdint.h>
#include <stdlib.h>
#include <immintrin.h>
//native interfaces
#include "native/electrosphere_server_fluid_simulator_FluidAcceleratedSimulator.h"
//fluid lib
#include "fluid/env/utilities.h"
#include "fluid/queue/chunkmask.h"
#include "fluid/queue/chunk.h"
#include "fluid/sim/grid2/grid2.h"
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/velocity.h"
#include "fluid/sim/grid2/density.h"
#include "fluid/sim/grid2/utilities.h"
#ifndef SAVE_STEPS
#define SAVE_STEPS 0
#endif
#define FLUID_GRID2_REALLY_SMALL_VALUE 0.00001
#define FLUID_GRID2_DIFFUSION_CONSTANT 0.00001
#define FLUID_GRID2_VISCOSITY_CONSTANT 0.00001
static inline void fluid_grid2_saveStep(float * values, const char * name);
static inline void fluid_grid2_applyGravity(Chunk * currentChunk, Environment * environment);
static inline void fluid_grid2_clearArr(float ** d);
static inline void fluid_grid2_rewrite_bounds(Chunk * chunk);
/**
* A grid that stores a mask of the border locations
*/
float * fluid_grid2_border_mask;
float * fluid_grid2_border_mask_inverted;
/**
* An array that stores values for the neighbor's arrays
*/
float * fluid_grid2_neighborArr_d;
float * fluid_grid2_neighborArr_d0;
float * fluid_grid2_neighborArr_u;
float * fluid_grid2_neighborArr_v;
float * fluid_grid2_neighborArr_w;
float * fluid_grid2_neighborArr_u0;
float * fluid_grid2_neighborArr_v0;
float * fluid_grid2_neighborArr_w0;
float * fluid_grid2_neighborArr_bounds;
void fluid_grid2_simulate(
int numChunks,
Chunk ** passedInChunks,
Environment * environment,
jfloat timestep
){
Chunk ** chunks = passedInChunks;
//solve chunk mask
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
//update the bounds arrays
fluid_grid2_rewrite_bounds(currentChunk);
//add velocity
fluid_grid2_applyGravity(currentChunk,environment);
fluid_grid2_addSourceToVectors(
currentChunk->chunkMask,
currentChunk->u,
currentChunk->v,
currentChunk->w,
currentChunk->u0,
currentChunk->v0,
currentChunk->w0,
FLUID_GRID2_DIFFUSION_CONSTANT,
FLUID_GRID2_VISCOSITY_CONSTANT,
timestep
);
//swap all vector fields
//swap vector fields
float * tmpArr;
for(int j = 0; j < 27; j++){
tmpArr = currentChunk->u[j];
currentChunk->u[j] = currentChunk->u0[j];
currentChunk->u0[j] = tmpArr;
}
for(int j = 0; j < 27; j++){
tmpArr = currentChunk->v[j];
currentChunk->v[j] = currentChunk->v0[j];
currentChunk->v0[j] = tmpArr;
}
for(int j = 0; j < 27; j++){
tmpArr = currentChunk->w[j];
currentChunk->w[j] = currentChunk->w0[j];
currentChunk->w0[j] = tmpArr;
}
//solve vector diffusion
for(int l = 0; l < FLUID_GRID2_VECTOR_DIFFUSE_TIMES; l++){
//solve vector diffusion
fluid_grid2_solveVectorDiffuse(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
//update array for vectors
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
}
//setup projection
fluid_grid2_setupProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
//update array for vectors
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_NO_DIR,currentChunk->v0);
//samples u0, v0
//sets u0
//these should have just been mirrored in the above
//
//Perform main projection solver
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
}
//samples u,v,w,u0
//sets u,v,w
//Finalize projection
fluid_grid2_finalizeProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
//set boundaries a final time for u,v,w
//...
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
//swap all vector fields
//swap vector fields
for(int j = 0; j < 27; j++){
tmpArr = currentChunk->u[j];
currentChunk->u[j] = currentChunk->u0[j];
currentChunk->u0[j] = tmpArr;
}
for(int j = 0; j < 27; j++){
tmpArr = currentChunk->v[j];
currentChunk->v[j] = currentChunk->v0[j];
currentChunk->v0[j] = tmpArr;
}
for(int j = 0; j < 27; j++){
tmpArr = currentChunk->w[j];
currentChunk->w[j] = currentChunk->w0[j];
currentChunk->w0[j] = tmpArr;
}
//advect vectors across boundaries
//advect
fluid_grid2_advectVectors(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
//update neighbor arr
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
//setup projection
fluid_grid2_setupProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
//update array for vectors
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_U,currentChunk->u0);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_V,currentChunk->v0);
//samples u0, v0
//sets u0
//these should have just been mirrored in the above
//
//Perform main projection solver
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
}
//samples u,v,w,u0
//sets u,v,w
//Finalize projection
fluid_grid2_finalizeProjection(currentChunk->chunkMask,currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
//set boundaries a final time for u,v,w
//...
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
}
///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
///------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//add density
double deltaDensity = 0;
environment->state.existingDensity = 0;
environment->state.newDensity = 0;
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
fluid_grid2_addDensity(environment,currentChunk->chunkMask,currentChunk->d,currentChunk->d0,timestep);
//swap all density arrays
//swap vector fields
float * tmpArr;
for(int j = 0; j < 27; j++){
tmpArr = currentChunk->d[j];
currentChunk->d[j] = currentChunk->d0[j];
currentChunk->d0[j] = tmpArr;
}
//diffuse density
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveDiffuseDensity(currentChunk->chunkMask,currentChunk->d,currentChunk->d0,currentChunk->u,currentChunk->v,currentChunk->w,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_VISCOSITY_CONSTANT,timestep);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->d);
}
//swap all density arrays
//swap vector fields
for(int j = 0; j < 27; j++){
tmpArr = currentChunk->d[j];
currentChunk->d[j] = currentChunk->d0[j];
currentChunk->d0[j] = tmpArr;
}
//advect density
fluid_grid2_advectDensity(currentChunk->chunkMask,currentChunk->d,currentChunk->d0,currentChunk->u,currentChunk->v,currentChunk->w,timestep);
}
//mirror densities
{
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->d);
}
}
//normalize densities
{
double transformedDensity = 0;
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
transformedDensity = transformedDensity + fluid_grid2_calculateSum(currentChunk->chunkMask,currentChunk->d);
}
float normalizationRatio = 0;
if(transformedDensity != 0){
normalizationRatio = (environment->state.existingDensity + environment->state.newDensity) / transformedDensity;
environment->state.normalizationRatio = normalizationRatio;
}
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
fluid_grid2_normalizeDensity(currentChunk->d,normalizationRatio);
}
}
//clear delta arrays
{
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
fluid_grid2_clearArr(currentChunk->d0);
fluid_grid2_clearArr(currentChunk->u0);
fluid_grid2_clearArr(currentChunk->v0);
fluid_grid2_clearArr(currentChunk->w0);
}
}
}
static inline void fluid_grid2_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 < FLUID_GRID2_REALLY_SMALL_VALUE && val > -FLUID_GRID2_REALLY_SMALL_VALUE){
val = 0;
}
fprintf(fp, "%f\t", val);
}
fprintf(fp, "\n");
}
fprintf(fp, "\n");
}
fclose(fp);
}
}
/**
* Applies gravity to the chunk
* @param currentChunk The chunk to apply on
* @param environment The environment data of the world
*/
static inline void fluid_grid2_applyGravity(Chunk * currentChunk, Environment * environment){
int N = DIM;
for(int x = 0; x < DIM; x++){
for(int y = 0; y < DIM; y++){
for(int z = 0; z < DIM; z++){
GET_ARR_RAW(currentChunk->v0,CENTER_LOC)[IX(x,y,z)] = GET_ARR_RAW(currentChunk->v0,CENTER_LOC)[IX(x,y,z)] + GET_ARR_RAW(currentChunk->d,CENTER_LOC)[IX(x,y,z)] * environment->consts.gravity;
}
}
}
}
/**
* Clears an array
*/
static inline void fluid_grid2_clearArr(float ** d){
float * x = GET_ARR_RAW(d,CENTER_LOC);
for(int j = 0; j < DIM * DIM * DIM; j++){
x[j] = 0;
}
}
/**
* Allocates the arrays necessary for grid2 simulation
*/
void fluid_grid2_allocate_arrays(){
fluid_grid2_border_mask = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
fluid_grid2_border_mask_inverted = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
for(int x = 0; x < DIM; x++){
for(int y = 0; y < DIM; y++){
for(int z = 0; z < DIM; z++){
if(x == 0 || x == DIM-1 || y == 0 || y == DIM-1 || z == 0 || z == DIM-1){
fluid_grid2_border_mask[IX(x,y,z)] = 1;
fluid_grid2_border_mask_inverted[IX(x,y,z)] = 0;
} else {
fluid_grid2_border_mask[IX(x,y,z)] = 0;
fluid_grid2_border_mask_inverted[IX(x,y,z)] = 1;
}
}
}
}
fluid_grid2_neighborArr_d = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
fluid_grid2_neighborArr_d0 = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
fluid_grid2_neighborArr_u = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
fluid_grid2_neighborArr_v = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
fluid_grid2_neighborArr_w = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
fluid_grid2_neighborArr_u0 = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
fluid_grid2_neighborArr_v0 = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
fluid_grid2_neighborArr_w0 = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
fluid_grid2_neighborArr_bounds = (float *)calloc(1,DIM * DIM * DIM * sizeof(float));
}
/**
* Applies a bounds array to a source array
*/
static inline void fluid_grid2_apply_bounds_mask(float * realArr, float * boundsArr){
__m256 maskedBounds, realVal, invertedMask, realPart, finalVec;
int x;
for(int z = 0; z < 18; z++){
for(int y = 0; y < 18; y++){
//lower part
x = 0;
//border part
maskedBounds = _mm256_loadu_ps(&fluid_grid2_border_mask[IX(x,y,z)]);
//real part
realVal = _mm256_loadu_ps(&realArr[IX(x,y,z)]);
invertedMask = _mm256_loadu_ps(&fluid_grid2_border_mask_inverted[IX(x,y,z)]);
realPart = _mm256_mul_ps(realVal,invertedMask);
finalVec = _mm256_add_ps(realPart,maskedBounds);
_mm256_storeu_ps(&realArr[IX(x,y,z)],finalVec);
//upper part
x = 2;
//border part
maskedBounds = _mm256_loadu_ps(&fluid_grid2_border_mask[IX(x,y,z)]);
//real part
realVal = _mm256_loadu_ps(&realArr[IX(x,y,z)]);
invertedMask = _mm256_loadu_ps(&fluid_grid2_border_mask_inverted[IX(x,y,z)]);
realPart = _mm256_mul_ps(realVal,invertedMask);
finalVec = _mm256_add_ps(realPart,maskedBounds);
_mm256_storeu_ps(&realArr[IX(x,y,z)],finalVec);
}
}
}
/**
* Applies the bounds mask to the current chunk
*/
static inline void fluid_grid2_set_bounds(Chunk * chunk){
fluid_grid2_apply_bounds_mask(chunk->d[CENTER_LOC], fluid_grid2_neighborArr_d);
fluid_grid2_apply_bounds_mask(chunk->d0[CENTER_LOC], fluid_grid2_neighborArr_d0);
fluid_grid2_apply_bounds_mask(chunk->u[CENTER_LOC], fluid_grid2_neighborArr_u);
fluid_grid2_apply_bounds_mask(chunk->v[CENTER_LOC], fluid_grid2_neighborArr_v);
fluid_grid2_apply_bounds_mask(chunk->w[CENTER_LOC], fluid_grid2_neighborArr_w);
fluid_grid2_apply_bounds_mask(chunk->u0[CENTER_LOC], fluid_grid2_neighborArr_u0);
fluid_grid2_apply_bounds_mask(chunk->v0[CENTER_LOC], fluid_grid2_neighborArr_v0);
fluid_grid2_apply_bounds_mask(chunk->w0[CENTER_LOC], fluid_grid2_neighborArr_w0);
fluid_grid2_apply_bounds_mask(chunk->bounds[CENTER_LOC], fluid_grid2_neighborArr_bounds);
}
/**
* Quickly masks a chunk's arrays
*/
static inline void fluid_grid2_populate_masked_arr(float * sourceArr, float * workingArr){
__m256 arrVal, maskVal, masked;
for(int z = 0; z < 18; z++){
for(int y = 0; y < 18; y++){
//lower part
arrVal = _mm256_loadu_ps(&sourceArr[IX(0,y,z)]);
maskVal = _mm256_loadu_ps(&fluid_grid2_border_mask[IX(0,y,z)]);
masked = _mm256_mul_ps(arrVal,maskVal);
_mm256_storeu_ps(&workingArr[IX(0,y,z)],masked);
//upper part
arrVal = _mm256_loadu_ps(&sourceArr[IX(2,y,z)]);
maskVal = _mm256_loadu_ps(&fluid_grid2_border_mask[IX(0,y,z)]);
masked = _mm256_mul_ps(arrVal,maskVal);
_mm256_storeu_ps(&workingArr[IX(0,y,z)],masked);
}
}
}
/**
* Rewrites the bounds arrays to contain the bounds of the current chunk
*/
static inline void fluid_grid2_rewrite_bounds(Chunk * chunk){
fluid_grid2_populate_masked_arr(chunk->d[CENTER_LOC], fluid_grid2_neighborArr_d);
fluid_grid2_populate_masked_arr(chunk->d0[CENTER_LOC], fluid_grid2_neighborArr_d0);
fluid_grid2_populate_masked_arr(chunk->u[CENTER_LOC], fluid_grid2_neighborArr_u);
fluid_grid2_populate_masked_arr(chunk->v[CENTER_LOC], fluid_grid2_neighborArr_v);
fluid_grid2_populate_masked_arr(chunk->w[CENTER_LOC], fluid_grid2_neighborArr_w);
fluid_grid2_populate_masked_arr(chunk->u0[CENTER_LOC], fluid_grid2_neighborArr_u0);
fluid_grid2_populate_masked_arr(chunk->v0[CENTER_LOC], fluid_grid2_neighborArr_v0);
fluid_grid2_populate_masked_arr(chunk->w0[CENTER_LOC], fluid_grid2_neighborArr_w0);
fluid_grid2_populate_masked_arr(chunk->bounds[CENTER_LOC], fluid_grid2_neighborArr_bounds);
}