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

projection solver tests
Some checks failed
studiorailgun/Renderer/pipeline/head There was a failure building this commit
Details

Browse Source
This commit is contained in:
austin 2024-12-09 21:40:39 -05:00
parent 305ec1d600
commit d0593f4e88
12 changed files with 355 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/main/c/includes/fluid/sim/grid2/solver_consts.h
View File

@ -15,7 +15,7 @@
/**
* Width of a single grid cell
*/
#define FLUID_GRID2_H (1)
#define FLUID_GRID2_H 1.0/DIM
/**
* Timestep to simulate by

1
src/main/c/includes/fluid/sim/grid2/utilities.h
View File

@ -17,7 +17,6 @@ void fluid_grid2_add_source(float * x, float * s, float dt);
* Sets the bounds of this cube to those of its neighbor
*/
LIBRARY_API void fluid_grid2_setBoundsToNeighborsRaw(
int chunk_mask,
int vector_dir,
float ** neighborArray
);

6
src/main/c/includes/fluid/sim/grid2/velocity.h
View File

@ -55,7 +55,7 @@ void fluid_grid2_addSourceToVectors(
* @param divr The grid that will be zeroed out in preparation of the solver
* @param dt The timestep for the simulation
*/
void fluid_grid2_setupProjection(
LIBRARY_API void fluid_grid2_setupProjection(
float ** ur,
float ** vr,
float ** wr,
@ -72,7 +72,7 @@ void fluid_grid2_setupProjection(
* @param jru0 The gradient field
* @param jrv0 The first derivative field
*/
void fluid_grid2_solveProjection(
LIBRARY_API void fluid_grid2_solveProjection(
float ** jru0,
float ** jrv0,
float dt
@ -84,7 +84,7 @@ void fluid_grid2_solveProjection(
* This subtracts the difference delta along the approximated gradient field.
* Thus we are left with an approximately mass-conserved field.
*/
void fluid_grid2_finalizeProjection(
LIBRARY_API void fluid_grid2_finalizeProjection(
float ** jru,
float ** jrv,
float ** jrw,

6
src/main/c/src/fluid/sim/grid2/density.c
View File

@ -89,8 +89,10 @@ LIBRARY_API void fluid_grid2_advectDensity(float ** d, float ** d0, float ** ur,
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;
dtx=dty=dtz=dt*DIM;
float h = FLUID_GRID2_H;
dtx = dt*h;
dty = dt*h;
dtz = dt*h;
float * center_d = GET_ARR_RAW(d,CENTER_LOC);
float * center_d0 = GET_ARR_RAW(d0,CENTER_LOC);

40
src/main/c/src/fluid/sim/grid2/fluidsim.c
View File

@ -82,16 +82,16 @@ void fluid_grid2_simulate(
//solve vector diffusion
fluid_grid2_solveVectorDiffuse(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,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);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
}
//setup projection
fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,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);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->v0);
//samples u0, v0
//sets u0
@ -100,7 +100,7 @@ void fluid_grid2_simulate(
//Perform main projection solver
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,timestep);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
}
//samples u,v,w,u0
//sets u,v,w
@ -109,9 +109,9 @@ void fluid_grid2_simulate(
//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);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
//swap all vector fields
//swap vector fields
@ -123,14 +123,14 @@ void fluid_grid2_simulate(
//advect
fluid_grid2_advectVectors(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,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);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
//setup projection
fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,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);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u0);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v0);
//samples u0, v0
//sets u0
//these should have just been mirrored in the above
@ -138,7 +138,7 @@ void fluid_grid2_simulate(
//Perform main projection solver
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,timestep);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
}
//samples u,v,w,u0
//sets u,v,w
@ -146,9 +146,9 @@ void fluid_grid2_simulate(
fluid_grid2_finalizeProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,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);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
}
@ -177,7 +177,7 @@ void fluid_grid2_simulate(
//diffuse density
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveDiffuseDensity(currentChunk->d,currentChunk->d0,FLUID_GRID2_DIFFUSION_CONSTANT,timestep);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->d);
fluid_grid2_setBoundsToNeighborsRaw(0,currentChunk->d);
}
//swap all density arrays
//swap vector fields
@ -190,7 +190,7 @@ void fluid_grid2_simulate(
{
for(int i = 0; i < numChunks; i++){
Chunk * currentChunk = chunks[i];
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->d);
fluid_grid2_setBoundsToNeighborsRaw(0,currentChunk->d);
}
}
//normalize densities

2
src/main/c/src/fluid/sim/grid2/utilities.c
View File

@ -28,12 +28,10 @@ void fluid_grid2_add_source(float * x, float * s, float dt){
* Sets the bounds of this cube to those of its neighbor
*/
LIBRARY_API void fluid_grid2_setBoundsToNeighborsRaw(
int chunk_mask,
int vector_dir,
float ** neighborArray
){
float * target = GET_ARR_RAW(neighborArray,CENTER_LOC);
float * source;
//set the faces bounds
for(int x=1; x < DIM-1; x++){
for(int y = 1; y < DIM-1; y++){

30
src/main/c/src/fluid/sim/grid2/velocity.c
View File

@ -151,7 +151,7 @@ LIBRARY_API void fluid_grid2_solveVectorDiffuse (
* @param VISCOSITY_CONST The viscosity constant
* @param dt The timestep for the simulation
*/
void fluid_grid2_setupProjection(
LIBRARY_API void fluid_grid2_setupProjection(
float ** ur,
float ** vr,
float ** wr,
@ -163,7 +163,7 @@ void fluid_grid2_setupProjection(
float h = FLUID_GRID2_H;
__m256 nVector = _mm256_set1_ps(1);
__m256 constScalar = _mm256_set1_ps(-1.0/(2.0*h));
__m256 constScalar = _mm256_set1_ps(-0.5 * h);
__m256 zeroVec = _mm256_set1_ps(0);
__m256 vector, vector2, vector3;
@ -183,15 +183,12 @@ void fluid_grid2_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,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);
@ -205,15 +202,12 @@ void fluid_grid2_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,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);
@ -231,16 +225,15 @@ void fluid_grid2_setupProjection(
* @param jru0 The gradient field
* @param jrv0 The first derivative field
*/
void fluid_grid2_solveProjection(
LIBRARY_API void fluid_grid2_solveProjection(
float ** jru0,
float ** jrv0,
float dt
){
int a = 1;
int c = 6;
float c = 6;
int i, j, k, l, m;
__m256 aScalar = _mm256_set1_ps(a);
__m256 cScalar = _mm256_set1_ps(c);
__m256 vector;
float * p = GET_ARR_RAW(jru0,CENTER_LOC);
float * div = GET_ARR_RAW(jrv0,CENTER_LOC);
@ -250,21 +243,26 @@ void fluid_grid2_solveProjection(
int n = 0;
//solve as much as possible vectorized
for(i = 1; i < DIM-1; i=i+8){
__m256 vector = _mm256_loadu_ps(&p[IX(i-1,j,k)]);
//compute Q
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);
//add A
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&div[IX(i,j,k)]));
//divide by 6
vector = _mm256_div_ps(vector,cScalar);
_mm256_storeu_ps(&p[IX(i,j,k)],vector);
}
//If there is any leftover, perform manual solving
if(i>DIM-1){
for(i=i-8; i < DIM-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;
p[IX(i,j,k)] = (div[IX(i,j,k)] + (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;
}
}
}
@ -276,7 +274,7 @@ void fluid_grid2_solveProjection(
* This subtracts the difference delta along the approximated gradient field.
* Thus we are left with an approximately mass-conserved field.
*/
void fluid_grid2_finalizeProjection(
LIBRARY_API void fluid_grid2_finalizeProjection(
float ** jru,
float ** jrv,
float ** jrw,

96
src/test/c/fluid/sim/grid2/advect_projection_tests.c Normal file
View File

@ -0,0 +1,96 @@
#include <math.h>
#include "stb/stb_ds.h"
#include "fluid/queue/boundsolver.h"
#include "fluid/queue/chunkmask.h"
#include "fluid/queue/chunk.h"
#include "fluid/env/environment.h"
#include "fluid/env/utilities.h"
#include "fluid/sim/grid2/density.h"
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
/**
* Center of the advection cell
*/
#define FLUID_GRID2_PROJECTION_CELL_CENTER 24
/**
* Testing velocity advection
*/
int fluid_sim_grid2_advect_projection_test1(){
printf("fluid_sim_grid2_advect_projection_test1\n");
int rVal = 0;
Environment * env = fluid_environment_create();
Chunk ** queue = NULL;
queue = createChunkGrid(env,3,3,3);
//setup chunk values
Chunk * currentChunk = queue[0];
currentChunk->d[CENTER_LOC][IX(2,2,2)] = MAX_FLUID_VALUE;
advection_setup_convection_cell(queue, FLUID_GRID2_PROJECTION_CELL_CENTER);
float beforeSumX = chunk_queue_sum_velocity(queue,FLUID_GRID2_BOUND_DIR_U);
float beforeSumY = chunk_queue_sum_velocity(queue,FLUID_GRID2_BOUND_DIR_V);
float beforeSumZ = chunk_queue_sum_velocity(queue,FLUID_GRID2_BOUND_DIR_W);
//actually simulate
int frameCount = 1;
for(int frame = 0; frame < frameCount; frame++){
int chunkCount = arrlen(queue);
for(int chunkIndex = 0; chunkIndex < 1; chunkIndex++){
currentChunk = queue[chunkIndex];
//advect
fluid_grid2_flip_arrays(currentChunk->u,currentChunk->u0);
fluid_grid2_flip_arrays(currentChunk->v,currentChunk->v0);
fluid_grid2_flip_arrays(currentChunk->w,currentChunk->w0);
fluid_grid2_advectVectors(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_SIM_STEP);
fluid_grid2_flip_arrays(currentChunk->u,currentChunk->u0);
fluid_grid2_flip_arrays(currentChunk->v,currentChunk->v0);
fluid_grid2_flip_arrays(currentChunk->w,currentChunk->w0);
////project
fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u0);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v0);
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
}
fluid_grid2_finalizeProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
}
}
//test the result
float afterSumX = chunk_queue_sum_velocity(queue,FLUID_GRID2_BOUND_DIR_U);
float afterSumY = chunk_queue_sum_velocity(queue,FLUID_GRID2_BOUND_DIR_V);
float afterSumZ = chunk_queue_sum_velocity(queue,FLUID_GRID2_BOUND_DIR_W);
if(fabs(beforeSumX - afterSumX) > FLUID_GRID2_REALLY_SMALL_VALUE){
rVal += assertEqualsFloat(beforeSumX,afterSumX,"Velocity advection step changed x-velocity sum! %f %f \n");
}
if(fabs(beforeSumY - afterSumY) > FLUID_GRID2_REALLY_SMALL_VALUE){
rVal += assertEqualsFloat(beforeSumX,afterSumX,"Velocity advection step changed y-density sum! %f %f \n");
}
if(fabs(beforeSumZ - afterSumZ) > FLUID_GRID2_REALLY_SMALL_VALUE){
rVal += assertEqualsFloat(beforeSumX,afterSumX,"Velocity advection step changed z-density sum! %f %f \n");
}
return rVal;
}
/**
* Testing velocity advection
*/
int fluid_sim_grid2_advect_projection_tests(int argc, char **argv){
int rVal = 0;
// rVal += fluid_sim_grid2_advect_projection_test1();
return rVal;
}

4
src/test/c/fluid/sim/grid2/density_diffuse_tests.c
View File

@ -76,7 +76,7 @@ int fluid_sim_grid2_density_diffuse_test1(){
//diffuse density
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveDiffuseDensity(currentChunk->d,currentChunk->d0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_SIM_STEP);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->d);
fluid_grid2_setBoundsToNeighborsRaw(0,currentChunk->d);
}
//swap all density arrays
//swap vector fields
@ -134,7 +134,7 @@ int fluid_sim_grid2_density_diffuse_test2(){
//diffuse density
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveDiffuseDensity(currentChunk->d,currentChunk->d0,FLUID_GRID2_DIFFUSION_CONSTANT,FLUID_GRID2_SIM_STEP);
fluid_grid2_setBoundsToNeighborsRaw(currentChunk->chunkMask,0,currentChunk->d);
fluid_grid2_setBoundsToNeighborsRaw(0,currentChunk->d);
}
//swap all density arrays
//swap vector fields

88
src/test/c/fluid/sim/grid2/setup_projection_tests.c Normal file
View File

@ -0,0 +1,88 @@
#include <math.h>
#include "stb/stb_ds.h"
#include "fluid/queue/boundsolver.h"
#include "fluid/queue/chunkmask.h"
#include "fluid/queue/chunk.h"
#include "fluid/env/environment.h"
#include "fluid/env/utilities.h"
#include "fluid/sim/grid2/density.h"
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
/**
* Center of the advection cell
*/
#define FLUID_GRID2_PROJECTION_CELL_CENTER 24
/**
* Testing velocity advection
*/
int fluid_sim_grid2_setup_projection_test1(){
printf("fluid_sim_grid2_setup_projection_test1\n");
int rVal = 0;
Environment * env = fluid_environment_create();
Chunk ** queue = NULL;
queue = createChunkGrid(env,3,3,3);
//setup chunk values
Chunk * currentChunk = queue[0];
currentChunk->u[CENTER_LOC][IX(3,3,3)] = 1.0f;
//actually simulate
fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
//test the result
rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],-0.5f / DIM,"Divergence of the vector field at 3,3,3 should be -0.5/DIM! %f %f \n");
rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(3,3,3)],0,"Divergence of the vector field at 3,3,3 should be 0! %f %f \n");
rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],0.5f / DIM,"Divergence of the vector field at 4,3,3 should be 0.5/DIM! %f %f \n");
rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(5,3,3)],0,"Divergence of the vector field at 5,3,3 should be 0! %f %f \n");
return rVal;
}
/**
* Testing velocity advection
*/
int fluid_sim_grid2_setup_projection_test2(){
printf("fluid_sim_grid2_setup_projection_test1\n");
int rVal = 0;
Environment * env = fluid_environment_create();
Chunk ** queue = NULL;
queue = createChunkGrid(env,3,3,3);
//setup chunk values
Chunk * currentChunk = queue[0];
currentChunk->u[CENTER_LOC][IX(3,3,3)] = -1.0f;
//actually simulate
fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
//test the result
rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(2,3,3)],0.5f / DIM,"Divergence of the vector field at 3,3,3 should be 0.5/DIM! %f %f \n");
rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(3,3,3)],0,"Divergence of the vector field at 3,3,3 should be 0! %f %f \n");
rVal += assertEqualsFloat(currentChunk->v0[CENTER_LOC][IX(4,3,3)],-0.5f / DIM,"Divergence of the vector field at 4,3,3 should be -0.5/DIM! %f %f \n");
return rVal;
}
/**
* Testing velocity advection
*/
int fluid_sim_grid2_setup_projection_tests(int argc, char **argv){
int rVal = 0;
rVal += fluid_sim_grid2_setup_projection_test1();
rVal += fluid_sim_grid2_setup_projection_test2();
return rVal;
}

121
src/test/c/fluid/sim/grid2/solve_projection_tests.c Normal file
View File

@ -0,0 +1,121 @@
#include <math.h>
#include "stb/stb_ds.h"
#include "fluid/queue/boundsolver.h"
#include "fluid/queue/chunkmask.h"
#include "fluid/queue/chunk.h"
#include "fluid/env/environment.h"
#include "fluid/env/utilities.h"
#include "fluid/sim/grid2/density.h"
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
/**
* Center of the advection cell
*/
#define FLUID_GRID2_PROJECTION_CELL_CENTER 24
/**
* Error margin for tests
*/
#define FLUID_GRId2_PROJECTION_ERROR_MARGIN 0.00001f
/**
* Testing velocity advection
*/
int fluid_sim_grid2_solve_projection_test1(){
printf("fluid_sim_grid2_solve_projection_test1\n");
int rVal = 0;
Environment * env = fluid_environment_create();
Chunk ** queue = NULL;
queue = createChunkGrid(env,3,3,3);
//setup chunk values
Chunk * currentChunk = queue[0];
currentChunk->u[CENTER_LOC][IX(3,3,3)] = 1.0f;
fluid_grid2_setupProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
//actually simulate
for(int l = 0; l < FLUID_GRID2_LINEARSOLVERTIMES; l++){
fluid_grid2_solveProjection(currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_NO_DIR,currentChunk->u0);
}
//test the result
//divergence of the gradient should be ___ above and below
// rVal += assertEqualsFloat(currentChunk->u0[CENTER_LOC][IX(3,2,3)],0,"First derivative of the scalar at 3,2,3 should be 0! %f %f \n");
float expected, actual;
{
//2,3,3
expected =
currentChunk->u0[CENTER_LOC][IX(1,3,3)] +
currentChunk->u0[CENTER_LOC][IX(3,3,3)] +
currentChunk->u0[CENTER_LOC][IX(2,2,3)] +
currentChunk->u0[CENTER_LOC][IX(2,4,3)] +
currentChunk->u0[CENTER_LOC][IX(2,3,2)] +
currentChunk->u0[CENTER_LOC][IX(2,3,4)]
;
expected = expected + currentChunk->v0[CENTER_LOC][IX(2,3,3)];
expected = expected / 6.0f;
actual = currentChunk->u0[CENTER_LOC][IX(2,3,3)];
if(fabs(expected - actual) > FLUID_GRId2_PROJECTION_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual," - Scalar potential at 2,3,3 is above error margin! expected: %f actual: %f \n");
}
}
{
//3,3,3
expected =
currentChunk->u0[CENTER_LOC][IX(2,3,3)] +
currentChunk->u0[CENTER_LOC][IX(4,3,3)] +
currentChunk->u0[CENTER_LOC][IX(3,2,3)] +
currentChunk->u0[CENTER_LOC][IX(3,4,3)] +
currentChunk->u0[CENTER_LOC][IX(3,3,2)] +
currentChunk->u0[CENTER_LOC][IX(3,3,4)]
;
expected = expected + currentChunk->v0[CENTER_LOC][IX(3,3,3)]; //this should be -0.5 * h * divergence of vector field
expected = expected / 6.0f;
actual = currentChunk->u0[CENTER_LOC][IX(3,3,3)];
if(fabs(expected - actual) > FLUID_GRId2_PROJECTION_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual," - Scalar potential at 3,3,3 is above error margin! expected: %f actual: %f \n");
}
}
{
//4,3,3
expected =
currentChunk->u0[CENTER_LOC][IX(3,3,3)] +
currentChunk->u0[CENTER_LOC][IX(5,3,3)] +
currentChunk->u0[CENTER_LOC][IX(4,2,3)] +
currentChunk->u0[CENTER_LOC][IX(4,4,3)] +
currentChunk->u0[CENTER_LOC][IX(4,3,2)] +
currentChunk->u0[CENTER_LOC][IX(4,3,4)]
;
expected = expected + currentChunk->v0[CENTER_LOC][IX(4,3,3)]; //this should be -0.5 * h * divergence of vector field
expected = expected / 6.0f;
actual = currentChunk->u0[CENTER_LOC][IX(4,3,3)];
if(fabs(expected - actual) > FLUID_GRId2_PROJECTION_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual," - Scalar potential at 4,3,3 is above error margin! expected: %f actual: %f \n");
}
}
return rVal;
}
/**
* Testing velocity advection
*/
int fluid_sim_grid2_solve_projection_tests(int argc, char **argv){
int rVal = 0;
rVal += fluid_sim_grid2_solve_projection_test1();
return rVal;
}

12
src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c
View File

@ -56,9 +56,9 @@ int fluid_sim_grid2_velocity_diffuse_test1(){
//solve vector diffusion
fluid_grid2_solveVectorDiffuse(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_SIM_STEP);
//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);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
}
//swap all density arrays
//swap vector fields
@ -116,9 +116,9 @@ int fluid_sim_grid2_velocity_diffuse_test2(){
//solve vector diffusion
fluid_grid2_solveVectorDiffuse(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,currentChunk->w0,FLUID_GRID2_SIM_STEP);
//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);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_U,currentChunk->u);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_V,currentChunk->v);
fluid_grid2_setBoundsToNeighborsRaw(FLUID_GRID2_BOUND_DIR_W,currentChunk->w);
}
//swap all density arrays
//swap vector fields