More pressurecell work

src/main/c/includes/fluid/sim/pressurecell/solver_consts.h

@@ -104,7 +104,7 @@
 /**
  * Enables clamping small density values to 0
  */
-#define FLUID_PRESSURECELL_ENABLE_CLAMP_MIN_DENSITY 0
+#define FLUID_PRESSURECELL_ENABLE_CLAMP_MIN_DENSITY 1
 
 /**
  * Enables renormalizing the velocity field to a max value of 1 during projection

src/main/c/src/fluid/sim/pressurecell/pressure.c

@@ -4,7 +4,10 @@
 #include "fluid/sim/pressurecell/solver_consts.h"
 #include "math/ode/multigrid_parallel.h"
 
-
+/**
+ * Divergence level we warn at
+ */
+#define DIVERGENCE_WARNING_VALUE 100.0f
 
 /**
  * Approximates the pressure for this chunk
@@ -85,7 +88,7 @@ LIBRARY_API void pressurecell_approximate_pressure(Environment * environment, Ch
                 phi0[IX(x,y,z)] = divCache[IX(x,y,z)];
                 pressureTemp[IX(x,y,z)] = pressureCache[IX(x,y,z)] * FLUID_PRESSURECELL_PRESSURE_BACKDOWN_FACTOR;
                 // pressureTemp[IX(x,y,z)] = 0;
-                if(divCache[IX(x,y,z)] > 3){
+                if(divCache[IX(x,y,z)] > DIVERGENCE_WARNING_VALUE){
                     printf("invalid divergence!\n");
                     printf("%f \n", divCache[IX(x,y,z)]);
                     printf("\n");
@@ -294,7 +297,7 @@ LIBRARY_API void pressurecell_approximate_divergence(Environment * environment,
                 // }
                 newDivergence = (du+dv+dw) * (-0.5f * FLUID_PRESSURECELL_SPACING);
                 // divArr[IX(x,y,z)] = divArr[IX(x,y,z)] + newDivergence - FLUID_PRESSURECELL_RESIDUAL_MULTIPLIER * divArr[IX(x,y,z)] + outflowDiv;
-                if(newDivergence > FLUID_PRESSURECELL_MAX_DIVERGENCE || newDivergence < -FLUID_PRESSURECELL_MAX_DIVERGENCE){
+                if(newDivergence > DIVERGENCE_WARNING_VALUE || newDivergence < -DIVERGENCE_WARNING_VALUE){
                     printf("Invalid divergence! \n");
                     printf("%f \n",newDivergence);
                     printf("%f %f \n", uArr[IX(x+1,y,z)], uArr[IX(x-1,y,z)]);

src/main/c/src/fluid/sim/pressurecell/pressurecell.c

@@ -143,6 +143,7 @@ LIBRARY_API void fluid_pressurecell_simulate(
         Chunk * currentChunk = chunks[i];
         pressurecell_advect_density(environment,currentChunk);
         if(FLUID_PRESSURECELL_ENABLE_RECAPTURE){
+            //d->dTemp->d
             fluid_pressurecell_recapture_density(environment,currentChunk);
         }
     }
@@ -153,6 +154,7 @@ LIBRARY_API void fluid_pressurecell_simulate(
     //
     for(int i = 0; i < numChunks; i++){
         Chunk * currentChunk = chunks[i];
+        //d->d
         fluid_pressurecell_normalize_chunk(environment,currentChunk);
     }
 
@@ -171,6 +173,7 @@ LIBRARY_API void fluid_pressurecell_simulate(
     //
     for(int i = 0; i < numChunks; i++){
         Chunk * currentChunk = chunks[i];
+        //uTemp->u
         pressurecell_copy_for_next_frame(environment,currentChunk);
         fluid_pressurecell_clearArr(currentChunk->d0[CENTER_LOC]);
         fluid_pressurecell_clearArr(currentChunk->u0[CENTER_LOC]);

src/main/c/src/fluid/sim/pressurecell/velocity.c

@@ -6,6 +6,11 @@
 #include "math/ode/multigrid_parallel.h"
 
 
+/**
+ * Velocity magnitude at which we warn
+ */
+#define WARNING_VELOCITY_MAGNITUDE 1000
+
 /**
  * Adds velocity from the delta buffer to this chunk
 */
@@ -45,7 +50,7 @@ LIBRARY_API void pressurecell_add_velocity(Environment * environment, Chunk * ch
                 uTemp[IX(x,y,z)] = uArr[IX(x,y,z)] + uSourceArr[IX(x,y,z)];
                 vTemp[IX(x,y,z)] = vArr[IX(x,y,z)] + vSourceArr[IX(x,y,z)];
                 wTemp[IX(x,y,z)] = wArr[IX(x,y,z)] + wSourceArr[IX(x,y,z)];
-                if(vTemp[IX(x,y,z)] > 5){
+                if(vTemp[IX(x,y,z)] > WARNING_VELOCITY_MAGNITUDE){
                     printf("Invalid add velocity!\n");
                     printf("%f %f %f  \n",  vTemp[IX(x,y,z)], vArr[IX(x,y,z)], vSourceArr[IX(x,y,z)] );
                     printf("\n");
@@ -179,7 +184,7 @@ LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk
                     )
                 ) * a
                 ;
-                if(vArr[IX(x,y,z)] > 5){
+                if(vArr[IX(x,y,z)] > WARNING_VELOCITY_MAGNITUDE){
                     printf("Invalid diffuse!\n");
                     printf("%f \n", vArr[IX(x,y,z)]);
                     printf("%f\n", vTemp[IX(x,y,z)]);
@@ -242,6 +247,7 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk *
     float xp, yp, zp;
     float s0, s1, t0, t1, u0, u1;
     float interpolatedU, interpolatedV, interpolatedW;
+    float vecU, vecV, vecW;
     float magnitude, maxMagnitude;
     float interpConst = environment->consts.dt / (FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING);
     for(y = 1; y < DIM-1; y++){
@@ -258,10 +264,31 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk *
     for(y = 1; y < DIM-1; y++){
         for(z = 1; z < DIM-1; z++){
             for(x = 1; x < DIM-1; x++){
+                //figure how far we're advecting
+                vecU = uArr[IX(x,y,z)] * interpConst;
+                vecV = vArr[IX(x,y,z)] * interpConst;
+                vecW = wArr[IX(x,y,z)] * interpConst;
+                if(vecU > 0.999f){
+                    vecU = 0.999f;
+                } else if(vecU < -0.999f){
+                    vecU = -0.999f;
+                }
+                if(vecV > 0.999f){
+                    vecV = 0.999f;
+                } else if(vecV < -0.999f){
+                    vecV = -0.999f;
+                }
+                if(vecW > 0.999f){
+                    vecW = 0.999f;
+                } else if(vecW < -0.999f){
+                    vecW = -0.999f;
+                }
+
+
                 //calculate the real (float) position we are at
-                xp = x - uArr[IX(x,y,z)] * interpConst;
+                xp = x - vecU;
-                yp = y - vArr[IX(x,y,z)] * interpConst;
+                yp = y - vecV;
-                zp = z - wArr[IX(x,y,z)] * interpConst;
+                zp = z - vecW;
 
                 //clamp to border
                 x0 = xp;
@@ -380,7 +407,7 @@ LIBRARY_API void pressurecell_advect_velocity(Environment * environment, Chunk *
                 //     interpolatedW = interpolatedW / magnitude;
                 // }
 
-                if(magnitude > 10){
+                if(magnitude > WARNING_VELOCITY_MAGNITUDE){
                     printf("advect invalid set: %f  %f  %f  \n",  uArr[IX(x,y,z)],  vArr[IX(x,y,z)],  wArr[IX(x,y,z)]);
                     int a = 1;
                     int b = 0;
@@ -498,7 +525,7 @@ LIBRARY_API double pressurecell_project_velocity(Environment * environment, Chun
                     maxMagnitude = magnitude;
                 }
 
-                if(magnitude != magnitude || magnitude > 10){
+                if(magnitude != magnitude || magnitude > WARNING_VELOCITY_MAGNITUDE){
                     printf("invalid magnitude! %f\n", magnitude);
                     printf("%f %f %f\n",  pressureDifferenceX, pressureDifferenceY, pressureDifferenceZ);
                     printf("%f %f %f\n",  uArr[IX(x,y,z)], vArr[IX(x,y,z)], wArr[IX(x,y,z)]);