Open boundary conditions

doc/notes.txt

@@ -0,0 +1,56 @@
+
+Colocated multigrid
+
+    +---+---+---+---+---+---+---+---+---+
+ 8  | x | x | x | x | x | x | x | x | x |
+    +---+---+---+---+---+---+---+---+---+
+ 7  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 6  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 5  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 4  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 3  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 2  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 1  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 0  | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | x |
+    +---+---+---+---+---+---+---+---+---+
+      0   1   2   3   4   5   6   7   8
+
+- Pressure is defined on all cells:
+	- The 0 cells are stored but not used (we return p = 0 on these cells)
+	- The x cells are not stored (we implicitly return p = 0 out of bounds)
+
+- Velocity is only defined in the inner cells, we always return 0 on the borders
+
+The Coarse grid is as follows:
+
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   x   |   x   |   x   |   x   |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   0   |       |       |       |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   0   |       |       |       |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   0   |       |       |       |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   0   |   0   |   0   |   0   |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+
+So, the size of the total grid is (2^N+1)*(2^N+1), but we store only (2^N)*(2^N), and explicitly
+set the pressure to 0 on the lower and left borders.

notes.txt

@@ -1,11 +0,0 @@
-On macOS, homebrew version of llvm must be used instead of xcode's one:
-
-> brew install llvm
-
-Have to run a local http server to serve files:
-
-> python3 -m http.server
-
-
-brew llvm toolchain seems to be broken and does not find open gl header files, and trying to include them from the framework path in /applications/xcode/ etc creates a bunch of other errors, so the workaround was to download mesa openGL headers, and tweak glext.h so that it includes stdint.h instead of inttypes.h (because brew's llvm can't seem to find this one). This was annoying as fuck.
-Zig cc seem to have everything needed but doesn't recognize --export-all and such, so I could'nt find a way to use it.

src/main.c

@@ -323,7 +323,6 @@ GLuint create_texture(int width, int height, GLenum internalFormat, GLenum forma
 	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
 	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
 	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
-
 	return(texture);
 }
 
@@ -695,8 +694,8 @@ void WAFNDraw()
 		float splatPosX = mouseInput.x/canvasWidth;
 		float splatPosY = 1 - mouseInput.y/canvasHeight;
 
-		float splatVelX = 10000.*DELTA*mouseInput.deltaX/canvasWidth;
-		float splatVelY = -10000.*DELTA*mouseInput.deltaY/canvasWidth;
+		float splatVelX = 5000.*DELTA*mouseInput.deltaX/canvasWidth;
+		float splatVelY = -5000.*DELTA*mouseInput.deltaY/canvasWidth;
 
 		float intensity = 100*sqrtf(square(mouseInput.deltaX/canvasWidth) + square(mouseInput.deltaY/canvasHeight));
 
@@ -712,7 +711,7 @@ void WAFNDraw()
 
 	//DEBUG
 	//*
-	if((int)(frameT*3/2.) % 6 == 0)
+	if((int)(frameT*3/2.) % 9 == 0)
 	{
 		float dirX = 0.2*cosf(2*3.141516*frameT*7.26);
 		float dirY = 0.3+0.1*sinf(frameT);

src/shaders/advect.glsl

@@ -11,9 +11,54 @@ uniform sampler2D velocity;
 uniform float delta;
 uniform float dissipation;
 
+vec2 u(ivec2 coord)
+{
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(velocity, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(velocity, 0).y)
+	{
+		return(vec2(0.));
+	}
+
+	return(texelFetch(velocity, coord, 0).xy);
+}
+
+vec4 q(ivec2 coord)
+{
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(src, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(src, 0).y)
+	{
+		return(vec4(0.));
+	}
+	return(texelFetch(src, coord, 0));
+}
+
+vec4 bilerpSrc(vec2 pos)
+{
+	vec2 offset = fract(pos);
+
+	ivec2 bl = ivec2(floor(pos));
+
+	ivec2 br = bl + ivec2(1, 0);
+	ivec2 tl = bl + ivec2(0, 1);
+	ivec2 tr = bl + ivec2(1, 1);
+
+	vec4 lerpTop = (1.-offset.x)*q(tl) + offset.x*q(tr);
+	vec4 lerpBottom = (1.-offset.x)*q(bl) + offset.x*q(br);
+	vec4 result = (1.-offset.y)*lerpBottom + offset.y*lerpTop;
+
+	return(result);
+}
+
 void main()
 {
-	vec2 u = texture(velocity, texCoord).xy;
-	vec2 samplePos = texCoord - delta * u;
-	fragColor = texture(src, samplePos) / (1. + dissipation*delta);
+	float texWidth = float(textureSize(velocity, 0).x);
+
+	ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy));
+
+	vec2 samplePos = vec2(pixelCoord) - texWidth * delta * u(pixelCoord);
+	fragColor = bilerpSrc(samplePos) / (1. + dissipation*delta);
 }

src/shaders/divergence.glsl

@@ -10,11 +10,25 @@ uniform sampler2D src;
 
 float ux(ivec2 coord)
 {
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(src, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(src, 0).y)
+	{
+		return(0.);
+	}
 	return(texelFetch(src, coord, 0).x);
 }
 
 float uy(ivec2 coord)
 {
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(src, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(src, 0).y)
+	{
+		return(0.);
+	}
 	return(texelFetch(src, coord, 0).y);
 }
 
@@ -32,12 +46,5 @@ void main()
 	float t = uy(vt);
 	float b = uy(vb);
 
-/*
-	vec2 c = texelFetch(src, pixelCoord, 0).xy;
-	if(vr.x >= textureSize(src, 0).x) { r = -c.x; }
-	if(vl.x < 0) { l = -c.x; }
-	if(vt.y >= textureSize(src, 0).y) { t = -c.y; }
-	if(vb.y < 0) { b = -c.y; }
-*/
 	fragColor = vec4(-(r - l + t - b)/2., 0, 0, 1);
 }

src/shaders/jacobi_step.glsl

@@ -12,9 +12,9 @@ uniform sampler2D bTex;
 float x(ivec2 coord)
 {
 	if(  coord.x <= 0
-	  || coord.x > textureSize(xTex, 0).x
+	  || coord.x >= textureSize(xTex, 0).x
 	  || coord.y <= 0
-	  || coord.y > textureSize(xTex, 0).y)
+	  || coord.y >= textureSize(xTex, 0).y)
 	{
 		return(0.);
 	}
@@ -23,6 +23,13 @@ float x(ivec2 coord)
 
 float b(ivec2 coord)
 {
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(bTex, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(bTex, 0).y)
+	{
+		return(0.);
+	}
 	return(texelFetch(bTex, coord, 0).x);
 }
 

src/shaders/multigrid_correct.glsl

@@ -12,9 +12,28 @@ uniform float invGridSize;
 
 float e(ivec2 coord)
 {
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(error, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(error, 0).y)
+	{
+		return(0.);
+	}
 	return(texelFetch(error, coord, 0).x);
 }
 
+float p(ivec2 coord)
+{
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(src, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(src, 0).y)
+	{
+		return(0.);
+	}
+	return(texelFetch(src, coord, 0).x);
+}
+
 void main()
 {
 	ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy));
@@ -30,5 +49,5 @@ void main()
 	float bottomLerp = (1.-offset.x)*e(bl) + offset.x*e(br);
 	float bilerpError = (1.-offset.y)*bottomLerp + offset.y*topLerp;
 
-	fragColor = vec4(texture(src, texCoord).x + bilerpError, 0, 0, 1);
+	fragColor = vec4(p(pixelCoord) + bilerpError, 0, 0, 1);
 }

src/shaders/multigrid_restrict_residual.glsl

@@ -12,9 +12,9 @@ uniform sampler2D bTex;
 float x(ivec2 coord)
 {
 	if(  coord.x <= 0
-	  || coord.x > textureSize(xTex, 0).x
+	  || coord.x >= textureSize(xTex, 0).x
 	  || coord.y <= 0
-	  || coord.y > textureSize(xTex, 0).y)
+	  || coord.y >= textureSize(xTex, 0).y)
 	{
 		return(0.);
 	}
@@ -24,9 +24,9 @@ float x(ivec2 coord)
 float b(ivec2 coord)
 {
 	if(  coord.x <= 0
-	  || coord.x > textureSize(bTex, 0).x
+	  || coord.x >= textureSize(bTex, 0).x
 	  || coord.y <= 0
-	  || coord.y > textureSize(bTex, 0).y)
+	  || coord.y >= textureSize(bTex, 0).y)
 	{
 		return(0.);
 	}

src/shaders/subtract_pressure.glsl

@@ -12,16 +12,25 @@ uniform float invGridSize;
 
 vec2 u(ivec2 coord)
 {
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(src, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(src, 0).y)
+	{
+		return(vec2(0.));
+	}
 	return(texelFetch(src, coord, 0).xy);
 }
 
 float p(ivec2 coord)
 {
-	if(coord.x < 0) { coord.x = 0; }
-	if(coord.x >= textureSize(pressure, 0).x) { coord.x = textureSize(pressure, 0).x-1; }
-	if(coord.y < 0) { coord.y = 0; }
-	if(coord.y >= textureSize(pressure, 0).y) { coord.y = textureSize(pressure, 0).y-1; }
-
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(pressure, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(pressure, 0).y)
+	{
+		return(0.);
+	}
 	return(texelFetch(pressure, coord, 0).x);
 }