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

This commit is contained in:
austin 2024-12-09 22:20:38 -05:00
parent d0593f4e88
commit 19bd6d365d
7 changed files with 178 additions and 319 deletions

View File

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

View File

@ -45,8 +45,7 @@ LIBRARY_API void fluid_grid2_solveDiffuseDensity(
float DIFFUSION_CONST,
float dt
){
float h = FLUID_GRID2_H;
float a=dt*DIFFUSION_CONST/(h*h);
float a=dt*DIFFUSION_CONST/(FLUID_GRID2_H*FLUID_GRID2_H);
float c=1+6*a;
int i, j, k, l, m;
float * x = GET_ARR_RAW(d,CENTER_LOC);
@ -89,10 +88,9 @@ 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;
float h = FLUID_GRID2_H;
dtx = dt*h;
dty = dt*h;
dtz = dt*h;
dtx = dt*FLUID_GRID2_H;
dty = dt*FLUID_GRID2_H;
dtz = dt*FLUID_GRID2_H;
float * center_d = GET_ARR_RAW(d,CENTER_LOC);
float * center_d0 = GET_ARR_RAW(d0,CENTER_LOC);
@ -110,116 +108,57 @@ LIBRARY_API void fluid_grid2_advectDensity(float ** d, float ** d0, float ** ur,
y = j-dty*v[IX(i,j,k)];
z = k-dtz*w[IX(i,j,k)];
m = n = o = 1;
if(x < 1){ m -= 1; }
if(x >= DIM-1){ m += 1; }
if(y < 1){ n -= 1; }
if(y >= DIM-1){ n += 1; }
if(z < 1){ o -= 1; }
if(z >= DIM-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){
// 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");
// 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);
// }
if(x < 0.001f){
//cases to consider:
//m = 0, x = -10
//m = 2, x = 0.01
x=0.001f;
i0=(int)0;
//clamp location within chunk
//get indices, and calculate percentage to pull from each index
if(x < 0.5f){
x=0.5f;
i0=0;
i1=1;
s0 = 0.999f;
s1 = 0.001f;
} else if(x >= DIM - 1){
//cases to consider:
//m = 0, x = 17.01
//m = 2, x = 20
x = DIM-1;
i0=(int)DIM-2;
} else if(x > (DIM - 2) + 0.5f){
x = (DIM - 2) + 0.5f;
i0=DIM-2;
i1=DIM-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 (y<0.5f) y=0.5f;
if (y>DIM+0.5f) y=DIM+0.5f;
if (z<0.5f) z=0.5f;
if (z>DIM+0.5f) z=DIM+0.5f;
if(y < 0.5f){
y=0.5f;
j0=(int)0;
j1=1;
} else if(y > (DIM - 2) + 0.5f){
y = (DIM - 2) + 0.5f;
j0=DIM-2;
j1=DIM-1;
} else {
j0=(int)y;
j1=j0+1;
}
//get actual indices
// 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;
if(z < 0.5f){
z=0.5f;
k0=(int)0;
k1=1;
} else if(z > (DIM - 2) + 0.5f){
z = (DIM - 2) + 0.5f;
k0=DIM-2;
k1=DIM-1;
} else {
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 >= DIM){
i0 = DIM - 1;
}
// if(i0 < 0){
// i0 = 0;
// }
if(j0 >= DIM){
j0 = DIM - 1;
}
// if(j0 < 0){
// j0 = 0;
// }
if(k0 >= DIM){
k0 = DIM - 1;
}
// if(k0 < 0){
// k0 = 0;
// }
if(i1 >= DIM){
i1 = DIM - 1;
}
// if(i1 < 0){
// i1 = 0;
// }
if(j1 >= DIM){
j1 = DIM - 1;
}
// if(j1 < 0){
// j1 = 0;
// }
if(k1 >= DIM){
k1 = DIM - 1;
}
// if(k1 < 0){
// k1 = 0;
// }
center_d[IX(i,j,k)] =
s0*(
t0*u0*center_d0[IX(i0,j0,k0)]+

View File

@ -45,8 +45,7 @@ LIBRARY_API void fluid_grid2_solveVectorDiffuse (
float ** jrw0,
float dt
){
float h = FLUID_GRID2_H;
float a=dt*FLUID_GRID2_VISCOSITY_CONSTANT/(h*h);
float a=dt*FLUID_GRID2_VISCOSITY_CONSTANT/(FLUID_GRID2_H*FLUID_GRID2_H);
float c=1+6*a;
int i, j, k, l, m;
float * u = GET_ARR_RAW(jru,CENTER_LOC);
@ -161,9 +160,7 @@ LIBRARY_API void fluid_grid2_setupProjection(
){
int i, j, k;
float h = FLUID_GRID2_H;
__m256 nVector = _mm256_set1_ps(1);
__m256 constScalar = _mm256_set1_ps(-0.5 * h);
__m256 constScalar = _mm256_set1_ps(-0.5 * FLUID_GRID2_H);
__m256 zeroVec = _mm256_set1_ps(0);
__m256 vector, vector2, vector3;
@ -283,8 +280,7 @@ LIBRARY_API void fluid_grid2_finalizeProjection(
float dt
){
int i, j, k;
float h = FLUID_GRID2_H;
__m256 constScalar = _mm256_set1_ps(2.0f*h);
__m256 constScalar = _mm256_set1_ps(2.0f*FLUID_GRID2_H);
__m256 vector, vector2, vector3;
float * u = GET_ARR_RAW(jru,CENTER_LOC);
@ -391,210 +387,59 @@ void fluid_grid2_advect_velocity(int b, float ** jrd, float ** jrd0, float * u,
y = j-dty*v[IX(i,j,k)];
z = k-dtz*w[IX(i,j,k)];
m = n = o = 1;
if(x < 0){ m += 1; }
else if(x >= DIM){ m -= 1; }
if(y < 0){ n += 1; }
else if(y >= DIM){ n -= 1; }
if(z < 0){ o += 1; }
else if(z >= DIM){ 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 && GET_ARR_RAW(jrd,CK(m,n,o)) != NULL){
// // 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_RAW(jrd0,CK(m,n,o));
// x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * (DIM-2);
// // printf("%d => %f\n",m,x);
// y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * (DIM-2);
// z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * (DIM-2);
// }
//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;
if(x < 0.5f){
x=0.5f;
i0=0;
i1=1;
s0 = 0.999f;
s1 = 0.001f;
} else if(x > DIM - 1){
//cases to consider:
//m = 0, x = 17.01
//m = 2, x = 20
x = DIM-1;
i0=(int)DIM-2;
} else if(x > (DIM - 2) + 0.5f){
x = (DIM - 2) + 0.5f;
i0=DIM-2;
i1=DIM-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;
if(y < 0.5f){
y=0.5f;
j0=(int)0;
j1=1;
t0 = 0.999f;
t1 = 0.001f;
} else if(y > DIM - 1){
//cases to consider:
//m = 0, x = 17.01
//m = 2, x = 20
y = DIM-1;
j0=(int)DIM-2;
} else if(y > (DIM - 2) + 0.5f){
y = (DIM - 2) + 0.5f;
j0=DIM-2;
j1=DIM-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;
if(z < 0.5f){
z=0.5f;
k0=(int)0;
k1=1;
u0 = 0.999f;
u1 = 0.001f;
} else if(z > DIM - 1){
//cases to consider:
//m = 0, x = 17.01
//m = 2, x = 20
z = DIM-1;
k0=(int)DIM-2;
} else if(z > (DIM - 2) + 0.5f){
z = (DIM - 2) + 0.5f;
k0=DIM-2;
k1=DIM-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;
//get actual indices
// i0=(int)x;
// i1=i0+1;
// j0=(int)y;
// j1=j0+1;
// k0=(int)z;
// k1=k0+1;
s1 = x-i0;
s0 = 1-s1;
//calculate percentage of each index
// s1 = x-i0;
// s0 = 1-s1;
// t1 = y-j0;
// t0 = 1-t1;
// u1 = z-k0;
// u0 = 1-u1;
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 >= DIM){
j0 = DIM - 1;
}
// if(j0 < 0){
// j0 = 0;
// }
if(k0 >= DIM){
k0 = DIM - 1;
}
// if(k0 < 0){
// k0 = 0;
// }
// if(i1 >= N){
// i1 = N - 1;
// }
// if(i1 < 0){
// i1 = 0;
// }
if(j1 >= DIM){
j1 = DIM - 1;
}
// if(j1 < 0){
// j1 = 0;
// }
if(k1 >= DIM){
k1 = DIM - 1;
}
// if(k1 < 0){
// k1 = 0;
// }
d[IX(i,j,k)] =
s0*(
t0*u0*d0[IX(i0,j0,k0)]+

View File

@ -19,6 +19,10 @@
*/
#define FLUID_GRID2_PROJECTION_CELL_CENTER 24
/**
* Error margin for tests
*/
#define FLUID_GRID2_PROJECTION_ERROR_MARGIN 0.00001f
/**
* Testing velocity advection
@ -34,11 +38,9 @@ int fluid_sim_grid2_advect_projection_test1(){
//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);
currentChunk->d[CENTER_LOC][IX(4,4,4)] = MAX_FLUID_VALUE;
currentChunk->u[CENTER_LOC][IX(4,4,4)] = MAX_FLUID_VALUE;
float beforeSum = chunk_queue_sum_density(queue);
//actually simulate
int frameCount = 1;
@ -46,7 +48,7 @@ int fluid_sim_grid2_advect_projection_test1(){
int chunkCount = arrlen(queue);
for(int chunkIndex = 0; chunkIndex < 1; chunkIndex++){
currentChunk = queue[chunkIndex];
//advect
//advect velocity
fluid_grid2_flip_arrays(currentChunk->u,currentChunk->u0);
fluid_grid2_flip_arrays(currentChunk->v,currentChunk->v0);
fluid_grid2_flip_arrays(currentChunk->w,currentChunk->w0);
@ -64,21 +66,17 @@ int fluid_sim_grid2_advect_projection_test1(){
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);
//advect density
fluid_grid2_flip_arrays(currentChunk->d,currentChunk->d0);
fluid_grid2_advectDensity(currentChunk->d,currentChunk->d0,currentChunk->u,currentChunk->v,currentChunk->w,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");
float afterSum = chunk_queue_sum_density(queue);
if(fabs(beforeSum - afterSum) > FLUID_GRID2_PROJECTION_ERROR_MARGIN){
rVal += assertEqualsFloat(beforeSum,afterSum,"Advection changed density! %f %f \n");
}
return rVal;

View File

@ -0,0 +1,85 @@
#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"
/**
* Error margin for tests
*/
#define FLUID_GRID2_PROJECTION_ERROR_MARGIN 0.00001f
/**
* Testing velocity advection
*/
int fluid_sim_grid2_finalize_projection_test1(){
printf("fluid_sim_grid2_finalize_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);
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);
}
//finalize
fluid_grid2_finalizeProjection(currentChunk->u,currentChunk->v,currentChunk->w,currentChunk->u0,currentChunk->v0,FLUID_GRID2_SIM_STEP);
//test the result
float expected, actual;
{
float xVel_at_2_3_3 = 0;
//2,3,3
expected = currentChunk->u0[CENTER_LOC][IX(3,3,3)] - currentChunk->u0[CENTER_LOC][IX(1,3,3)];
expected = expected / (2.0f * FLUID_GRID2_H);
expected = xVel_at_2_3_3 - expected;
actual = currentChunk->u[CENTER_LOC][IX(2,3,3)];
if(fabs(expected - actual) > FLUID_GRID2_PROJECTION_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual," - Conservative velocity at 2,3,3 is above error margin! expected: %f actual: %f \n");
}
}
{
float xVel_at_3_3_3 = 1.0f;
//3,3,3
expected = currentChunk->u0[CENTER_LOC][IX(4,3,3)] - currentChunk->u0[CENTER_LOC][IX(2,3,3)];
expected = expected / (2.0f * FLUID_GRID2_H);
expected = xVel_at_3_3_3 - expected;
actual = currentChunk->u[CENTER_LOC][IX(3,3,3)];
if(fabs(expected - actual) > FLUID_GRID2_PROJECTION_ERROR_MARGIN){
rVal += assertEqualsFloat(expected,actual," - Conservative velocity at 3,3,3 is above error margin! expected: %f actual: %f \n");
}
}
return rVal;
}
/**
* Testing velocity advection
*/
int fluid_sim_grid2_finalize_projection_tests(int argc, char **argv){
int rVal = 0;
rVal += fluid_sim_grid2_finalize_projection_test1();
return rVal;
}

View File

@ -14,18 +14,13 @@
#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
* Testing gradient approximation
*/
int fluid_sim_grid2_solve_projection_test1(){
printf("fluid_sim_grid2_solve_projection_test1\n");
@ -41,16 +36,13 @@ int fluid_sim_grid2_solve_projection_test1(){
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
//actually solve
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;
{
@ -110,7 +102,7 @@ int fluid_sim_grid2_solve_projection_test1(){
/**
* Testing velocity advection
* Testing gradient approximation
*/
int fluid_sim_grid2_solve_projection_tests(int argc, char **argv){
int rVal = 0;

View File

@ -148,12 +148,12 @@ int fluid_sim_grid2_velocity_advection_test3(){
int frameCount = 50;
for(int frame = 0; frame < frameCount; frame++){
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);
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);
}
//test the result