diff --git a/src/main/c/includes/fluid/sim/grid2/solver_consts.h b/src/main/c/includes/fluid/sim/grid2/solver_consts.h index 4eefdff4..d54e7476 100644 --- a/src/main/c/includes/fluid/sim/grid2/solver_consts.h +++ b/src/main/c/includes/fluid/sim/grid2/solver_consts.h @@ -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 diff --git a/src/main/c/includes/fluid/sim/grid2/utilities.h b/src/main/c/includes/fluid/sim/grid2/utilities.h index 344d0de9..7bfefe42 100644 --- a/src/main/c/includes/fluid/sim/grid2/utilities.h +++ b/src/main/c/includes/fluid/sim/grid2/utilities.h @@ -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 ); diff --git a/src/main/c/includes/fluid/sim/grid2/velocity.h b/src/main/c/includes/fluid/sim/grid2/velocity.h index 55f30c14..b33e694d 100644 --- a/src/main/c/includes/fluid/sim/grid2/velocity.h +++ b/src/main/c/includes/fluid/sim/grid2/velocity.h @@ -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, diff --git a/src/main/c/src/fluid/sim/grid2/density.c b/src/main/c/src/fluid/sim/grid2/density.c index a8d3c24a..5defac5a 100644 --- a/src/main/c/src/fluid/sim/grid2/density.c +++ b/src/main/c/src/fluid/sim/grid2/density.c @@ -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); diff --git a/src/main/c/src/fluid/sim/grid2/fluidsim.c b/src/main/c/src/fluid/sim/grid2/fluidsim.c index 22ddab4d..362c29be 100644 --- a/src/main/c/src/fluid/sim/grid2/fluidsim.c +++ b/src/main/c/src/fluid/sim/grid2/fluidsim.c @@ -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 diff --git a/src/main/c/src/fluid/sim/grid2/utilities.c b/src/main/c/src/fluid/sim/grid2/utilities.c index 804232f1..64f97003 100644 --- a/src/main/c/src/fluid/sim/grid2/utilities.c +++ b/src/main/c/src/fluid/sim/grid2/utilities.c @@ -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++){ diff --git a/src/main/c/src/fluid/sim/grid2/velocity.c b/src/main/c/src/fluid/sim/grid2/velocity.c index 2b927373..ee6748d7 100644 --- a/src/main/c/src/fluid/sim/grid2/velocity.c +++ b/src/main/c/src/fluid/sim/grid2/velocity.c @@ -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, diff --git a/src/test/c/fluid/sim/grid2/advect_projection_tests.c b/src/test/c/fluid/sim/grid2/advect_projection_tests.c new file mode 100644 index 00000000..fc148dbf --- /dev/null +++ b/src/test/c/fluid/sim/grid2/advect_projection_tests.c @@ -0,0 +1,96 @@ +#include + +#include "stb/stb_ds.h" + +#include "fluid/queue/boundsolver.h" +#include "fluid/queue/chunkmask.h" +#include "fluid/queue/chunk.h" +#include "fluid/env/environment.h" +#include "fluid/env/utilities.h" +#include "fluid/sim/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; +} \ No newline at end of file diff --git a/src/test/c/fluid/sim/grid2/density_diffuse_tests.c b/src/test/c/fluid/sim/grid2/density_diffuse_tests.c index fa881f2c..506baead 100644 --- a/src/test/c/fluid/sim/grid2/density_diffuse_tests.c +++ b/src/test/c/fluid/sim/grid2/density_diffuse_tests.c @@ -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 diff --git a/src/test/c/fluid/sim/grid2/setup_projection_tests.c b/src/test/c/fluid/sim/grid2/setup_projection_tests.c new file mode 100644 index 00000000..bc14acce --- /dev/null +++ b/src/test/c/fluid/sim/grid2/setup_projection_tests.c @@ -0,0 +1,88 @@ +#include + +#include "stb/stb_ds.h" + +#include "fluid/queue/boundsolver.h" +#include "fluid/queue/chunkmask.h" +#include "fluid/queue/chunk.h" +#include "fluid/env/environment.h" +#include "fluid/env/utilities.h" +#include "fluid/sim/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; +} \ No newline at end of file diff --git a/src/test/c/fluid/sim/grid2/solve_projection_tests.c b/src/test/c/fluid/sim/grid2/solve_projection_tests.c new file mode 100644 index 00000000..ff47b26e --- /dev/null +++ b/src/test/c/fluid/sim/grid2/solve_projection_tests.c @@ -0,0 +1,121 @@ +#include + +#include "stb/stb_ds.h" + +#include "fluid/queue/boundsolver.h" +#include "fluid/queue/chunkmask.h" +#include "fluid/queue/chunk.h" +#include "fluid/env/environment.h" +#include "fluid/env/utilities.h" +#include "fluid/sim/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; +} \ No newline at end of file diff --git a/src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c b/src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c index 7a26603f..31f3d148 100644 --- a/src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c +++ b/src/test/c/fluid/sim/grid2/velocity_diffuse_tests.c @@ -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