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

implement conjugate gradient solver
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-10 15:56:39 -05:00
parent 19bd6d365d
commit c8c2ba7ad7
21 changed files with 615 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
.vscode/settings.json vendored
View File

@ -42,6 +42,9 @@
"velocity.h": "c",
"density.h": "c",
"grid2.h": "c",
"math.h": "c"
"math.h": "c",
"multigrid.h": "c",
"gauss_seidel.h": "c",
"ode_utils.h": "c"
}
}

6
docs/src/progress/renderertodo.md
View File

@ -1279,6 +1279,12 @@ cellular stability work
(12/09/2024)
cellular stability work
Add grid2 sim files
grid2 functioning + lots of tests
multigrid implementation
(12/10/2024)
Implement conjugate gradient solver
# TODO

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

@ -5,17 +5,17 @@
/**
* The number of times to relax most solvers
*/
#define FLUID_GRID2_LINEARSOLVERTIMES 20
#define FLUID_GRID2_LINEARSOLVERTIMES 3
/**
* The number of times to relax the vector diffusion solver
*/
#define FLUID_GRID2_VECTOR_DIFFUSE_TIMES 20
#define FLUID_GRID2_VECTOR_DIFFUSE_TIMES 3
/**
* Width of a single grid cell
*/
#define FLUID_GRID2_H (1.0/DIM)
#define FLUID_GRID2_H (1.0/(DIM-2))
/**
* Timestep to simulate by

23
src/main/c/includes/math/ode/conjugate_gradient.h Normal file
View File

@ -0,0 +1,23 @@
#ifndef MATH_CONJUGATE_GRADIENT_H
#define MATH_CONJUGATE_GRADIENT_H
/**
* Iniitalizes the conjugate gradient solver with the phi values
* @param phi The phi array
* @param phi0 The phi array from the last frame
* @param a The a const
* @param c The c const
*/
void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c);
/**
* Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method
* @param phi The phi array
* @param phi0 The phi array from the last frame
* @param a The a const
* @param c The c const
*/
void solver_conjugate_gradient_iterate(float * phi, float * phi0, float a, float c);
#endif

108
src/main/c/includes/math/ode/gauss_seidel.h Normal file
View File

@ -0,0 +1,108 @@
#ifndef MATH_GAUSS_SEIDEL_H
#define MATH_GAUSS_SEIDEL_H
#include "immintrin.h"
#include "fluid/queue/chunk.h"
#include "fluid/queue/chunkmask.h"
/**
* Gets the index into the array
*/
static inline int solver_gauss_seidel_get_index(int x, int y, int z, int N){
return (x + (N * y) + (N * N * z));
}
/**
* Relaxes an ODE matrix by 1 iteration of gauss seidel parallelized
* @param phi The phi array
* @param phi0 The phi array from the last frame
* @param a The a const
* @param c The c const
* @param gridDim The dimensions of the grid
*/
static inline void solver_gauss_seidel_iterate_parallel(float * phi, float * phi0, float a, float c, int gridDim){
int i, j, k, l, m;
__m256 aScalar = _mm256_set1_ps(a);
__m256 cScalar = _mm256_set1_ps(c);
//transform u direction
for(k=1; k<gridDim-1; k++){
for(j=1; j<gridDim-1; j++){
int n = 0;
//solve as much as possible vectorized
for(i = 1; i < gridDim-1; i=i+8){
__m256 vector = _mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i-1,j,k,gridDim)]);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i+1,j,k,gridDim)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j-1,k,gridDim)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j+1,k,gridDim)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k-1,gridDim)]));
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi[solver_gauss_seidel_get_index(i,j,k+1,gridDim)]));
vector = _mm256_mul_ps(vector,aScalar);
vector = _mm256_add_ps(vector,_mm256_loadu_ps(&phi0[solver_gauss_seidel_get_index(i,j,k,gridDim)]));
vector = _mm256_div_ps(vector,cScalar);
_mm256_storeu_ps(&phi[solver_gauss_seidel_get_index(i,j,k,gridDim)],vector);
}
//If there is any leftover, perform manual solving
if(i>gridDim-1){
for(i=i-8; i < gridDim-1; i++){
phi[solver_gauss_seidel_get_index(i,j,k,gridDim)] =
(
phi0[solver_gauss_seidel_get_index(i,j,k,gridDim)] +
a * (
phi[solver_gauss_seidel_get_index(i-1,j,k,gridDim)]+
phi[solver_gauss_seidel_get_index(i+1,j,k,gridDim)]+
phi[solver_gauss_seidel_get_index(i,j-1,k,gridDim)]+
phi[solver_gauss_seidel_get_index(i,j+1,k,gridDim)]+
phi[solver_gauss_seidel_get_index(i,j,k-1,gridDim)]+
phi[solver_gauss_seidel_get_index(i,j,k+1,gridDim)]
)
) / c;
}
}
}
}
}
/**
* Relaxes an ODE matrix by 1 iteration of gauss seidel serially
* @param phi The phi array
* @param phi0 The phi array from the last frame
* @param a The a const
* @param c The c const
* @param gridDim The dimensions of the grid
*/
static inline void solver_gauss_seidel_iterate_serial(float * phi, float * phi0, float a, float c, int gridDim){
int i, j, k, l, m;
__m256 aScalar = _mm256_set1_ps(a);
__m256 cScalar = _mm256_set1_ps(c);
//transform u direction
for(k=1; k<gridDim-1; k++){
for(j=1; j<gridDim-1; j++){
int n = 0;
//If there is any leftover, perform manual solving
for(i=1; i < gridDim-1; i++){
phi[solver_gauss_seidel_get_index(i,j,k,gridDim)] =
(
phi0[solver_gauss_seidel_get_index(i,j,k,gridDim)] +
a * (
phi[solver_gauss_seidel_get_index(i-1,j,k,gridDim)]+
phi[solver_gauss_seidel_get_index(i+1,j,k,gridDim)]+
phi[solver_gauss_seidel_get_index(i,j-1,k,gridDim)]+
phi[solver_gauss_seidel_get_index(i,j+1,k,gridDim)]+
phi[solver_gauss_seidel_get_index(i,j,k-1,gridDim)]+
phi[solver_gauss_seidel_get_index(i,j,k+1,gridDim)]
)
) / c;
}
}
}
}
#endif

27
src/main/c/includes/math/ode/multigrid.h Normal file
View File

@ -0,0 +1,27 @@
#ifndef MATH_SOLVER_MULTIGRID_H
#define MATH_SOLVER_MULTIGRID_H
#include "immintrin.h"
#include "fluid/queue/chunk.h"
#include "fluid/queue/chunkmask.h"
/**
* Relaxes an ODE matrix by 1 iteration of multigrid method
* @param phi The phi array
* @param phi0 The phi array from the last frame
* @param a The a const
* @param c The c const
*/
void solver_multigrid_iterate(float * phi, float * phi0, float a, float c);
/**
* Allocates the data for the multigrid solver
*/
LIBRARY_API void solver_multigrid_allocate();
#endif

18
src/main/c/includes/math/ode/ode_utils.h Normal file
View File

@ -0,0 +1,18 @@
#ifndef MATH_ODE_UTILS_H
#define MATH_ODE_UTILS_H
/**
* Gets the index into the array
*/
static inline int ode_index(int x, int y, int z, int N){
return (x + (N * y) + (N * N * z));
}
/**
* Computes the stencil of a given source array at a given position
*/
__m256 ode_poisson_stencil_parallel(float * source, int x, int y, int z);
#endif

2
src/main/c/src/fluid/queue/javainterface.c
View File

@ -16,6 +16,7 @@
#include "fluid/sim/grid2/grid2.h"
#include "fluid/sim/simulator.h"
#include "fluid/dispatch/dispatcher.h"
#include "math/ode/multigrid.h"
//defines
@ -143,6 +144,7 @@ JNIEXPORT void JNICALL Java_electrosphere_server_fluid_simulator_FluidAccelerate
//init grid2 sim
fluid_grid2_allocate_arrays();
solver_multigrid_allocate();
}
/**

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

@ -8,6 +8,9 @@
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/velocity.h"
#include "fluid/sim/grid2/utilities.h"
#include "math/ode/gauss_seidel.h"
#include "math/ode/multigrid.h"
#include "math/ode/conjugate_gradient.h"
#define SET_BOUND_IGNORE 0
#define SET_BOUND_USE_NEIGHBOR 1
@ -227,6 +230,7 @@ LIBRARY_API void fluid_grid2_solveProjection(
float ** jrv0,
float dt
){
float a = 1;
float c = 6;
int i, j, k, l, m;
__m256 cScalar = _mm256_set1_ps(c);
@ -234,36 +238,13 @@ LIBRARY_API void fluid_grid2_solveProjection(
float * p = GET_ARR_RAW(jru0,CENTER_LOC);
float * div = GET_ARR_RAW(jrv0,CENTER_LOC);
// update for each cell
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
int n = 0;
//solve as much as possible vectorized
for(i = 1; i < DIM-1; i=i+8){
//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)]));
//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)] + (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;
}
}
}
}
//transform u direction
solver_gauss_seidel_iterate_parallel(p,div,a,c,DIM);
solver_multigrid_iterate(p,div,a,c);
solver_gauss_seidel_iterate_parallel(p,div,a,c,DIM);
solver_conjugate_gradient_init(p,div,a,c);
solver_conjugate_gradient_iterate(p,div,a,c);
}
/**

143
src/main/c/src/math/ode/conjugate_gradient.c Normal file
View File

@ -0,0 +1,143 @@
#include <stdlib.h>
#include <math.h>
#include <immintrin.h>
#include "fluid/queue/chunk.h"
#include "math/ode/ode_utils.h"
/**
* The search direction array
*/
float * p = NULL;
/**
* The residual array
*/
float * r = NULL;
/**
* Iniitalizes the conjugate gradient solver with the phi values
* @param phi The phi array
* @param phi0 The phi array from the last frame
* @param a The a const
* @param c The c const
*/
void solver_conjugate_gradient_init(float * phi, float * phi0, float a, float c){
if(p == NULL){
p = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
}
if(r == NULL){
r = (float *)calloc(1,DIM*DIM*DIM*sizeof(float));
}
int i, j, k;
__m256 aScalar = _mm256_set1_ps(a);
__m256 cScalar = _mm256_set1_ps(c);
__m256 f, residual, stencil;
//iniitalize the r (residual) and p (search direction) arrays
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
int n = 0;
//solve as much as possible vectorized
for(i = 1; i < DIM-1; i=i+8){
//calculate the stencil applied to phi
//get values from neighbors
stencil = _mm256_loadu_ps(&phi[ode_index( i-1, j, k, DIM )]);
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i+1, j, k, DIM )]));
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i, j-1, k, DIM )]));
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i, j+1, k, DIM )]));
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i, j, k-1, DIM )]));
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index( i, j, k+1, DIM )]));
//multiply by a
stencil = _mm256_mul_ps(stencil,aScalar);
//add previous value
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&phi[ode_index(i,j,k,DIM)]));
//divide by 6
stencil = _mm256_div_ps(stencil,cScalar);
//grab the f value (the value of phi0)
f = _mm256_loadu_ps(&phi0[ode_index(i,j,k,DIM)]);
//subtract stencil from f
residual = _mm256_sub_ps(f,stencil);
//store in residual and searchDir
_mm256_storeu_ps(&r[ode_index(i,j,k,DIM)],residual);
_mm256_storeu_ps(&p[ode_index(i,j,k,DIM)],residual);
}
}
}
}
/**
* Iteratively solves an ODE matrix by 1 iteration of conjugate gradient method
* @param phi The phi array
* @param phi0 The phi array from the last frame
* @param a The a const
* @param c The c const
*/
void solver_conjugate_gradient_iterate(float * phi, float * phi0, float a, float c){
int i, j, k;
__m256 aScalar = _mm256_set1_ps(a);
__m256 cScalar = _mm256_set1_ps(c);
__m256 f, residual, stencil;
__m256 denominator, r0_squared, r1_squared, a_k, r_delta, beta, p0, p1, p_new, phi_old, phi1, r0, r1;
//iniitalize the r (residual) and p (search direction) arrays
for(k=1; k<DIM-1; k++){
for(j=1; j<DIM-1; j++){
int n = 0;
//solve as much as possible vectorized
for(i = 1; i < DIM-1; i=i+8){
//calculate the stencil applied to p
//get values from neighbors
stencil = _mm256_loadu_ps(&p[ode_index( i-1, j, k, DIM )]);
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i+1, j, k, DIM )]));
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i, j-1, k, DIM )]));
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i, j+1, k, DIM )]));
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i, j, k-1, DIM )]));
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index( i, j, k+1, DIM )]));
//multiply by a
stencil = _mm256_mul_ps(stencil,aScalar);
//add previous value
stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&p[ode_index(i,j,k,DIM)]));
//divide by 6
stencil = _mm256_div_ps(stencil,cScalar);
//load p0
p0 = _mm256_loadu_ps(&p[ode_index(i,j,k,DIM)]);
//calculate ak
denominator = _mm256_mul_ps(stencil,p0);
r0_squared = _mm256_mul_ps(_mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]),_mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]));
a_k = _mm256_div_ps(r0_squared,denominator);
//update phi
phi_old = _mm256_loadu_ps(&phi[ode_index(i,j,k,DIM)]);
phi1 = _mm256_mul_ps(p1,a_k);
phi1 = _mm256_add_ps(phi1,phi_old);
_mm256_storeu_ps(&phi[ode_index(i,j,k,DIM)],phi1);
//update residual
r_delta = _mm256_add_ps(stencil,a_k);
r0 = _mm256_loadu_ps(&r[ode_index(i,j,k,DIM)]);
r1 = _mm256_sub_ps(r0,r_delta);
_mm256_storeu_ps(&r[ode_index(i,j,k,DIM)],r1);
//calculate beta
r1_squared = _mm256_mul_ps(r1,r1);
beta = _mm256_div_ps(r1_squared,r0_squared);
//update p (search dir)
p_new = _mm256_mul_ps(beta,p0);
p1 = _mm256_add_ps(r1,p_new);
_mm256_storeu_ps(&p[ode_index(i,j,k,DIM)],p1);
}
}
}
}

105
src/main/c/src/math/ode/multigrid.c Normal file
View File

@ -0,0 +1,105 @@
#include <stdlib.h>
#include <stdio.h>
#include "fluid/queue/chunk.h"
#include "math/ode/gauss_seidel.h"
#include "math/ode/multigrid.h"
/**
* Dimension of the half resolution grid
*/
static int halfDim = ((DIM - 2) / 2) + 2;
/**
* Dimension of the quarter resolution grid
*/
static int quarterDim = ((DIM - 2) / 4) + 2;
/**
* The half resolution grids
*/
static float * halfGridPhi;
static float * halfGridPhi0;
/**
* The quarter resolution grids
*/
static float * quarterGridPhi;
static float * quarterGridPhi0;
/**
* Relaxes an ODE matrix by 1 iteration of multigrid method
* @param phi The phi array
* @param phi0 The phi array from the last frame
* @param a The a const
* @param c The c const
*/
void solver_multigrid_iterate(float * phi, float * phi0, float a, float c){
//
//half res
for(int x = 1; x < halfDim - 2; x++){
for(int y = 1; y < halfDim - 2; y++){
for(int z = 1; z < halfDim - 2; z++){
halfGridPhi[solver_gauss_seidel_get_index(x,y,z,halfDim)] = phi[solver_gauss_seidel_get_index(x*2,y*2,z*2,DIM)];
halfGridPhi0[solver_gauss_seidel_get_index(x,y,z,halfDim)] = phi0[solver_gauss_seidel_get_index(x*2,y*2,z*2,DIM)];
}
}
}
solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim);
//
//quarter res
//
//half res
for(int x = 1; x < quarterDim - 2; x++){
for(int y = 1; y < quarterDim - 2; y++){
for(int z = 1; z < quarterDim - 2; z++){
quarterGridPhi[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = halfGridPhi[solver_gauss_seidel_get_index(x*2,y*2,z*2,halfDim)];
quarterGridPhi0[solver_gauss_seidel_get_index(x,y,z,quarterDim)] = halfGridPhi0[solver_gauss_seidel_get_index(x*2,y*2,z*2,halfDim)];
}
}
}
solver_gauss_seidel_iterate_parallel(quarterGridPhi,quarterGridPhi0,a,c,quarterDim);
for(int x = 1; x < halfDim - 2; x++){
for(int y = 1; y < halfDim - 2; y++){
for(int z = 1; z < halfDim - 2; z++){
halfGridPhi[solver_gauss_seidel_get_index(x,y,z,halfDim)] = quarterGridPhi[solver_gauss_seidel_get_index(x/2,y/2,z/2,quarterDim)];
halfGridPhi0[solver_gauss_seidel_get_index(x,y,z,halfDim)] = quarterGridPhi0[solver_gauss_seidel_get_index(x/2,y/2,z/2,quarterDim)];
}
}
}
//
//half res
solver_gauss_seidel_iterate_parallel(halfGridPhi,halfGridPhi0,a,c,halfDim);
for(int x = 1; x < DIM - 2; x++){
for(int y = 1; y < DIM - 2; y++){
for(int z = 1; z < DIM - 2; z++){
phi[solver_gauss_seidel_get_index(x,y,z,DIM)] = halfGridPhi[solver_gauss_seidel_get_index(x/2,y/2,z/2,halfDim)];
phi0[solver_gauss_seidel_get_index(x,y,z,DIM)] = halfGridPhi0[solver_gauss_seidel_get_index(x/2,y/2,z/2,halfDim)];
}
}
}
}
/**
* Allocates the data for the multigrid solver
*/
LIBRARY_API void solver_multigrid_allocate(){
halfGridPhi = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
halfGridPhi0 = (float *)calloc(1,halfDim * halfDim * halfDim * sizeof(float));
quarterGridPhi = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
quarterGridPhi0 = (float *)calloc(1,quarterDim * quarterDim * quarterDim * sizeof(float));
}

23
src/main/c/src/math/ode/ode_utils.c Normal file
View File

@ -0,0 +1,23 @@
#include "immintrin.h"
#include "fluid/queue/chunk.h"
#include "math/ode/ode_utils.h"
// /**
// * Computes the stencil of a given source array at a given position
// */
// __m256 ode_poisson_stencil_parallel(float * source, float * sourcePrev, int x, int y, int z){
// __m256 POISSON_STENCIL_C_SCALAR = _mm256_set1_ps(6);
// //get values from neighbors
// __m256 stencil = _mm256_loadu_ps(&source[ode_index( x-1, y, z, DIM )]);
// stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x+1, y, z, DIM )]));
// stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x, y-1, z, DIM )]));
// stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x, y+1, z, DIM )]));
// stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x, y, z-1, DIM )]));
// stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&source[ode_index( x, y, z+1, DIM )]));
// //add previous value
// stencil = _mm256_add_ps(stencil,_mm256_loadu_ps(&sourcePrev[ode_index(x,y,z,DIM)]));
// //divide by 6
// stencil = _mm256_div_ps(stencil,POISSON_STENCIL_C_SCALAR);
// _mm256_storeu_ps(&source[ode_index(x,y,z,DIM)],stencil);
// }

75
src/test/c/fluid/sim/grid2/advect_projection_tests.c
View File

@ -11,6 +11,7 @@
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "math/ode/multigrid.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
@ -82,13 +83,87 @@ int fluid_sim_grid2_advect_projection_test1(){
return rVal;
}
/**
* Testing velocity advection
*/
int fluid_sim_grid2_advect_projection_compute_error_over_time(){
printf("fluid_sim_grid2_advect_projection_compute_error_over_time\n");
int rVal = 0;
Environment * env = fluid_environment_create();
Chunk ** queue = NULL;
queue = createChunkGrid(env,3,3,3);
for(int complexity = 2; complexity < DIM-1; complexity++){
//actually simulate
for(int frameCount = 1; frameCount < 50; frameCount++){
//setup chunk values
Chunk * currentChunk = queue[0];
chunk_fill(currentChunk,0);
for(int x = 1; x < complexity; x++){
for(int y = 1; y < complexity; y++){
for(int z = 1; z < complexity; z++){
currentChunk->d[CENTER_LOC][IX(x,y,z)] = MAX_FLUID_VALUE;
currentChunk->u[CENTER_LOC][IX(x,y,z)] = MAX_FLUID_VALUE;
}
}
}
float beforeSum = chunk_queue_sum_density(queue);
//solve
for(int frame = 0; frame < frameCount; frame++){
int chunkCount = arrlen(queue);
for(int chunkIndex = 0; chunkIndex < 1; chunkIndex++){
currentChunk = queue[chunkIndex];
//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);
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);
//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 afterSum = chunk_queue_sum_density(queue);
if(fabs(beforeSum - afterSum) > FLUID_GRID2_PROJECTION_ERROR_MARGIN){
printf("%f,",fabs(beforeSum - afterSum));
// rVal += assertEqualsFloat(beforeSum,afterSum,"Advection changed density! %f %f \n");
}
}
printf("\n");
}
rVal++;
return rVal;
}
/**
* Testing velocity advection
*/
int fluid_sim_grid2_advect_projection_tests(int argc, char **argv){
int rVal = 0;
solver_multigrid_allocate();
// rVal += fluid_sim_grid2_advect_projection_test1();
// rVal += fluid_sim_grid2_advect_projection_compute_error_over_time();
return rVal;
}

3
src/test/c/fluid/sim/grid2/density_advection_tests.c
View File

@ -10,6 +10,7 @@
#include "fluid/sim/grid2/density.h"
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "math/ode/multigrid.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
@ -227,6 +228,8 @@ int fluid_sim_grid2_density_advection_test5(){
int fluid_sim_grid2_density_advection_tests(int argc, char **argv){
int rVal = 0;
solver_multigrid_allocate();
rVal += fluid_sim_grid2_density_advection_test1();
rVal += fluid_sim_grid2_density_advection_test2();
rVal += fluid_sim_grid2_density_advection_test3();

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

@ -9,6 +9,7 @@
#include "fluid/sim/grid2/density.h"
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "math/ode/multigrid.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
@ -161,6 +162,8 @@ int fluid_sim_grid2_density_diffuse_test2(){
int fluid_sim_grid2_density_diffuse_tests(int argc, char **argv){
int rVal = 0;
solver_multigrid_allocate();
rVal += fluid_sim_grid2_density_diffuse_test1();
rVal += fluid_sim_grid2_density_diffuse_test2();

3
src/test/c/fluid/sim/grid2/finalize_projection_tests.c
View File

@ -11,6 +11,7 @@
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "math/ode/multigrid.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
@ -79,6 +80,8 @@ int fluid_sim_grid2_finalize_projection_test1(){
int fluid_sim_grid2_finalize_projection_tests(int argc, char **argv){
int rVal = 0;
solver_multigrid_allocate();
rVal += fluid_sim_grid2_finalize_projection_test1();
return rVal;

7
src/test/c/fluid/sim/grid2/setup_projection_tests.c
View File

@ -11,6 +11,7 @@
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "math/ode/multigrid.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
@ -81,8 +82,10 @@ int fluid_sim_grid2_setup_projection_test2(){
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();
solver_multigrid_allocate();
// rVal += fluid_sim_grid2_setup_projection_test1();
// rVal += fluid_sim_grid2_setup_projection_test2();
return rVal;
}

3
src/test/c/fluid/sim/grid2/solve_projection_tests.c
View File

@ -11,6 +11,7 @@
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "math/ode/multigrid.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
@ -107,6 +108,8 @@ int fluid_sim_grid2_solve_projection_test1(){
int fluid_sim_grid2_solve_projection_tests(int argc, char **argv){
int rVal = 0;
solver_multigrid_allocate();
rVal += fluid_sim_grid2_solve_projection_test1();
return rVal;

45
src/test/c/fluid/sim/grid2/speed_tests.c Normal file
View File

@ -0,0 +1,45 @@
#define TARGET_SIM_RADIUS 3.0
#define TARGET_SIM_DIAMETER (TARGET_SIM_RADIUS * 2) + 1
/**
* Ie probably shouldn't be simulating the sky
*/
#define PERCENT_ACTIVE_IN_RADIUS 0.5
#define CHUNKS_PER_PLAYER ((TARGET_SIM_DIAMETER * TARGET_SIM_DIAMETER * TARGET_SIM_DIAMETER) * PERCENT_ACTIVE_IN_RADIUS)
/**
* Say 8 players are logged in
*/
#define PLAYERS_PER_SERVER 8.0
//worst case scenario is they're all on different boats trying to reach one another
#define TARGET_CHUNKS (PLAYERS_PER_SERVER * CHUNKS_PER_PLAYER)
/**
* Number of threads on server to dedicate to fluid sim
*/
#define SIM_THREADS 4.0
/**
* Allow 8 milli seconds per thread to simulate
*/
#define TOTAL_ALLOWED_TIME_IN_MILLI_SECONDS 8.0
/**
* Allowed time per chunk
*/
#define TIME_PER_CHUNK (TOTAL_ALLOWED_TIME_IN_MILLI_SECONDS / (TARGET_CHUNKS / SIM_THREADS))
int fluid_sim_grid2_speed_tests(){
int rVal = 0;
int itmePErChunk = TIME_PER_CHUNK;
return rVal;
}

3
src/test/c/fluid/sim/grid2/velocity_advection_tests.c
View File

@ -11,6 +11,7 @@
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "math/ode/multigrid.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
@ -182,6 +183,8 @@ int fluid_sim_grid2_velocity_advection_test3(){
int fluid_sim_grid2_velocity_advection_tests(int argc, char **argv){
int rVal = 0;
solver_multigrid_allocate();
rVal += fluid_sim_grid2_velocity_advection_test1();
rVal += fluid_sim_grid2_velocity_advection_test2();
rVal += fluid_sim_grid2_velocity_advection_test3();

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

@ -9,6 +9,7 @@
#include "fluid/sim/grid2/solver_consts.h"
#include "fluid/sim/grid2/utilities.h"
#include "fluid/sim/grid2/velocity.h"
#include "math/ode/multigrid.h"
#include "../../../util/chunk_test_utils.h"
#include "../../../util/test.h"
@ -149,6 +150,8 @@ int fluid_sim_grid2_velocity_diffuse_test2(){
int fluid_sim_grid2_velocity_diffuse_tests(int argc, char **argv){
int rVal = 0;
solver_multigrid_allocate();
rVal += fluid_sim_grid2_velocity_diffuse_test1();
return rVal;