Newer
Older
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<link rel="stylesheet" type="text/css" href="dependencies/bootstrap.min.css">
<link rel="stylesheet" type="text/css" href="dependencies/flat-ui.min.css">
<link rel="stylesheet" type="text/css" href="main.css">
<script id="2d-vertex-shader" type="x-shader/x-vertex">
attribute vec2 a_position;
void main() {
gl_Position = vec4(a_position, 0, 1);
}
</script>
<script id="boundaryConditionsShader" type="x-shader/x-fragment">
precision mediump float;
uniform sampler2D u_texture;
uniform float u_scale;
uniform vec2 u_textureSize;
uniform vec2 u_obstaclePosition;//scaled to texture size
uniform float u_obstacleRad;
void main() {
vec2 fragCoord = gl_FragCoord.xy;
vec2 dir = fragCoord - 3.0*vec2(0.5, 0.5) - u_obstaclePosition;//not sure where this fac of 3 came from?
gl_FragColor = u_scale*texture2D(u_texture, (fragCoord + dir/dist*2.0)/u_textureSize);
return;
}
gl_FragColor = texture2D(u_texture, fragCoord/u_textureSize);
}
</script>
<script id="2d-render-shader" type="x-shader/x-fragment">
precision mediump float;
if (texture2D(u_obstacle, fragCoord/u_textureSize).x == 1.0){
float mat1 = texture2D(u_material, fragCoord/u_textureSize).x;
vec3 color = vec3(250.0/255.0, 235.0/255.0, 215.0/255.0);
gl_FragColor = vec4(mat1*color, 1);
<script id="gradientSubtractionShader" type="x-shader/x-fragment">
precision mediump float;
uniform sampler2D u_velocity;
uniform sampler2D u_pressure;
uniform vec2 u_textureSize;
uniform float u_const;
vec2 repeatBoundary(vec2 coord){
coord -= vec2(0.5, 0.5);
//if (coord.x < 0.0) coord.x = u_textureSize.x-1.0;
if (coord.x >= u_textureSize.x-1.0) coord.x = 0.0;
if (coord.y < 0.0) coord.y = u_textureSize.y-1.0;
else if (coord.y >= u_textureSize.y-1.0) coord.y = 0.0;
coord += vec2(0.5, 0.5);
return coord;
}
void main() {
vec2 fragCoord = gl_FragCoord.xy;
vec2 currentVelocity = texture2D(u_velocity, fragCoord/u_textureSize).xy;
float n = texture2D(u_pressure, repeatBoundary(fragCoord+vec2(0.0, 1.0))/u_textureSize).x;
float s = texture2D(u_pressure, repeatBoundary(fragCoord+vec2(0.0, -1.0))/u_textureSize).x;
float e = texture2D(u_pressure, (fragCoord+vec2(1.0, 0.0))/u_textureSize).x;
float w = texture2D(u_pressure, (fragCoord+vec2(-1.0, 0.0))/u_textureSize).x;
gl_FragColor = vec4(currentVelocity-u_const*vec2(e-w, n-s), 0, 0);
}
</script>
<script id="divergenceShader" type="x-shader/x-fragment">
precision mediump float;
uniform sampler2D u_velocity;
uniform vec2 u_textureSize;
vec2 repeatBoundary(vec2 coord){
coord -= vec2(0.5, 0.5);
if (coord.x < 0.0) coord.x = u_textureSize.x-1.0;
else if (coord.x >= u_textureSize.x-1.0) coord.x = 0.0;
if (coord.y < 0.0) coord.y = u_textureSize.y-1.0;
else if (coord.y >= u_textureSize.y-1.0) coord.y = 0.0;
coord += vec2(0.5, 0.5);
return coord;
}
void main() {
vec2 fragCoord = gl_FragCoord.xy;
//finite difference formulation of divergence
//periodic boundary
float n = texture2D(u_velocity, repeatBoundary(fragCoord+vec2(0.0, 1.0))/u_textureSize).y;
float s = texture2D(u_velocity, repeatBoundary(fragCoord+vec2(0.0, -1.0))/u_textureSize).y;
float e = texture2D(u_velocity, repeatBoundary(fragCoord+vec2(1.0, 0.0))/u_textureSize).x;
float w = texture2D(u_velocity, (fragCoord+vec2(-1.0, 0.0))/u_textureSize).x;
<script id="forceShader" type="x-shader/x-fragment">
precision mediump float;
uniform sampler2D u_velocity;
uniform vec2 u_textureSize;
uniform vec2 u_mouseCoord;
uniform vec2 u_mouseDir;
uniform float u_mouseEnable;
void main() {
vec2 fragCoord = gl_FragCoord.xy;
vec2 currentVelocity = texture2D(u_velocity, fragCoord/u_textureSize).xy;
if (u_mouseEnable == 1.0){
vec2 pxDist = fragCoord - u_mouseCoord;
currentVelocity += u_mouseDir*u_dt*exp(-(pxDist.x*pxDist.x+pxDist.y*pxDist.y)*u_reciprocalRadius);
}
gl_FragColor = vec4(currentVelocity, 0, 0);
}
</script>
uniform vec2 u_textureSize;
uniform float u_alpha;
uniform float u_reciprocalBeta;
void main() {
vec2 fragCoord = gl_FragCoord.xy;
vec2 currentState = texture2D(u_b, fragCoord/u_textureSize).xy;
vec2 n = texture2D(u_x, (fragCoord+vec2(0.0, 1.0))/u_textureSize).xy;
vec2 s = texture2D(u_x, (fragCoord+vec2(0.0, -1.0))/u_textureSize).xy;
vec2 e = texture2D(u_x, (fragCoord+vec2(1.0, 0.0))/u_textureSize).xy;
vec2 w = texture2D(u_x, (fragCoord+vec2(-1.0, 0.0))/u_textureSize).xy;
vec2 nextState = (n + s + e + w + u_alpha * currentState) * u_reciprocalBeta;
gl_FragColor = vec4(nextState, 0, 0);
}
<script id="advectShaderMat" type="x-shader/x-fragment">
uniform sampler2D u_velocity;
uniform sampler2D u_material;
vec2 bilinearInterp(vec2 pos, sampler2D texture, vec2 size){
//bilinear interp between nearest cells
vec2 pxCenter = vec2(0.5, 0.5);
vec2 ceiled = ceil(pos);
vec2 floored = floor(pos);
vec2 n = texture2D(texture, (ceiled+pxCenter)/size).xy;//actually ne
vec2 s = texture2D(texture, (floored+pxCenter)/size).xy;//actually sw
if (ceiled.x != floored.x){
vec2 se = texture2D(texture, (vec2(ceiled.x, floored.y)+pxCenter)/size).xy;
vec2 nw = texture2D(texture, (vec2(floored.x, ceiled.y)+pxCenter)/size).xy;
n = n*(pos.x-floored.x) + nw*(ceiled.x-pos.x);
s = se*(pos.x-floored.x) + s*(ceiled.x-pos.x);
}
vec2 materialVal = n;
if (ceiled.y != floored.y){
materialVal = n*(pos.y-floored.y) + s*(ceiled.y-pos.y);
}
return materialVal;
}
vec2 pxCenter = vec2(0.5, 0.5);
//bilinear interp
//vec2 currentVelocity = 1.0/u_scale*texture2D(u_velocity, fragCoord/u_textureSize).xy;
vec2 currentVelocity = 1.0/u_scale*bilinearInterp((fragCoord-pxCenter)*u_scale + pxCenter, u_velocity, u_textureSize*u_scale);
if (fragCoord.x < 1.0 || length(currentVelocity) == 0.0) {//boundary or no velocity
gl_FragColor = vec4(texture2D(u_material, fragCoord/u_textureSize).xy, 0, 0);
if (pos.x >= u_textureSize.x-1.0) {
gl_FragColor = vec4(0, 0, 0, 0);
return;
}
//periodic boundary in y
if (pos.y < 0.0) pos.y += u_textureSize.y-1.0;
if (pos.y >= u_textureSize.y-1.0) pos.y -= u_textureSize.y-1.0;
gl_FragColor = vec4(bilinearInterp(pos, u_material, u_textureSize), 0, 0);
}
</script>
<script id="advectShaderVel" type="x-shader/x-fragment">
precision mediump float;
uniform sampler2D u_velocity;
uniform sampler2D u_material;
uniform vec2 u_textureSize;
uniform float u_scale;
uniform float u_dt;
void main() {
vec2 fragCoord = gl_FragCoord.xy;
vec2 currentVelocity = u_scale*texture2D(u_velocity, fragCoord/u_textureSize).xy;
//implicitly solve advection
if (length(currentVelocity) == 0.0) {//boundary or no velocity
gl_FragColor = vec4(texture2D(u_material, fragCoord/u_textureSize).xy, 0, 0);
return;
}
vec2 pxCenter = vec2(0.5, 0.5);
vec2 pos = fragCoord - pxCenter - u_dt*currentVelocity;
if (pos.x < 1.0) {
gl_FragColor = vec4(1.0, 0, 0, 0);
return;
}
if (pos.x >= u_textureSize.x-1.0) {
gl_FragColor = vec4(0, 0, 0, 0);
if (pos.x >= u_textureSize.x-1.0) pos.x -= u_textureSize.x-1.0;
if (pos.y < 0.0) {
pos.y += u_textureSize.y-1.0;
}
if (pos.y >= u_textureSize.y-1.0) {
pos.y -= u_textureSize.y-1.0;
}
//bilinear interp between nearest cells
vec2 ceiled = ceil(pos);
vec2 floored = floor(pos);
vec2 n = texture2D(u_material, (ceiled+pxCenter)/u_textureSize).xy;//actually ne
vec2 s = texture2D(u_material, (floored+pxCenter)/u_textureSize).xy;//actually sw
vec2 se = texture2D(u_material, (vec2(ceiled.x, floored.y)+pxCenter)/u_textureSize).xy;
vec2 nw = texture2D(u_material, (vec2(floored.x, ceiled.y)+pxCenter)/u_textureSize).xy;
n = n*(pos.x-floored.x) + nw*(ceiled.x-pos.x);
s = se*(pos.x-floored.x) + s*(ceiled.x-pos.x);
}
if (ceiled.y != floored.y){
materialVal = n*(pos.y-floored.y) + s*(ceiled.y-pos.y);
}
}
</script>
<script type="text/javascript" src="dependencies/jquery-3.1.0.min.js"></script>
<script type="text/javascript" src="dependencies/flat-ui.min.js"></script>
<script type="text/javascript" src="GlBoilerplate.js"></script>
<script type="text/javascript" src="GPUMath.js"></script>
<script type="text/javascript" src="main.js"></script>
</head>
<body>
<canvas id="glcanvas"></canvas>
<a href="#" id="about">?</a>
<div class="modal fade" id="aboutModal" tabindex="-1" role="dialog" aria-labelledby="basicModal" aria-hidden="true">
<div class="modal-dialog modal-lg">
<div class="modal-content">
<div class="modal-body">
This simulation solves the <a href="https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations" target="_blank">Navier-Stokes equations</a> for incompressible fluid flow past an obstacle in a GPU fragment shader.
It exhibits a phenomenon called <a href="https://en.wikipedia.org/wiki/Vortex_shedding" target="_blank">vortex shedding</a>,
where vortices of alternating spin spontaneously emerge behind the obstacle.
To increase performance, I solved for the velocity vector field of the fluid at a lower resolution than I used to compute the distribution of material moving through the fluid (shown in black and white).
I used bilinear interpolation to smooth out artifacts caused by this speedup.
I ignored the viscous diffusion term from the Navier-Stokes formula to encourage better vortex formation (the implicit advection solving I'm using creates enough diffusion on its own for this system).
<br/><br/>
Click and drag to apply a force to the fluid.
<br/><br/>
To learn more about the math involved, check out the following sources:<br/>
<a href="http://http.developer.nvidia.com/GPUGems/gpugems_ch38.html" target="_blank">Fast Fluid Dynamics Simulation on the GPU</a> - a very well written tutorial about programming the Navier-Stokes equations on a GPU.
Though not WebGL specific, it was still very useful.<br/>
<a href="http://jamie-wong.com/2016/08/05/webgl-fluid-simulation/" target="_blank">Fluid Simulation (with WebGL demo)</a> - this article has some nice, interactive graphics that helped me debug my code.<br/>
<a href="http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/ns.pdf" target="_blank">Stable Fluids</a> - a paper about stable numerical methods for evaluating Navier-Stokes on a discrete grid.<br/>
Written by <a href="http://www.amandaghassaei.com/" target="_blank">Amanda Ghassaei</a> as a homework assingment for <a href="http://fab.cba.mit.edu/classes/MAS.864/" target="_blank">The Nature of Mathematical Modeling</a>, code on <a href="https://github.com/amandaghassaei/VortexShedding" target="_blank">Github</a>.