31:T47a,{ const t0 = newTexture(seed); const fb0 = newFramebuffer(t0); const thalf = newTexture(); const fbhalf = newFramebuffer(thalf); const t1 = newTexture(); const fb1 = newFramebuffer(t1); for (let i = 0; true; ++i) { gl.useProgram(timestepProgram); const u_Dt = gl.getUniformLocation(timestepProgram, "u_Dt"); gl.viewport(0, 0, W, H); for (let j = 0; j < 8; ++j) { if (midpointIntegration) { // RK-2 predictor and corrector steps: step(t0, t0, fbhalf, u_Dt, timestep * 0.5); step(t0, thalf, fb1, u_Dt, timestep); // Same thing, framebuffers reversed: step(t1, t1, fbhalf, u_Dt, timestep * 0.5); step(t1, thalf, fb0, u_Dt, timestep); } else { // Euler step: step(t0, t0, fb1, u_Dt, timestep); // Same thing, framebuffers reversed: step(t1, t1, fb0, u_Dt, timestep); } } gl.useProgram(renderProgram); gl.viewport(0, 0, canvas.width, canvas.height); gl.bindTexture(gl.TEXTURE_2D, t0); gl.bindFramebuffer(gl.FRAMEBUFFER, null); gl.drawArrays(gl.TRIANGLES, 0, 3); yield i; } }32:T4be,timestepShader = newShader(gl.FRAGMENT_SHADER, ` precision mediump float; uniform sampler2D u_f, u_deriv; uniform vec2 u_size; uniform float u_F, u_K, u_Dt; const float D_a = 0.2; const float D_b = 0.1; void main() { vec2 position = gl_FragCoord.xy; vec2 value = texture2D(u_deriv, position / u_size).xy; vec2 laplacian = -3.0 * value + 0.5 * ( texture2D(u_deriv, (position + vec2(0.0, 1.0)) / u_size).xy + texture2D(u_deriv, (position + vec2(1.0, 0.0)) / u_size).xy + texture2D(u_deriv, (position + vec2(0.0, -1.0)) / u_size).xy + texture2D(u_deriv, (position + vec2(-1.0, 0.0)) / u_size).xy ) + 0.25 * ( texture2D(u_deriv, (position + vec2(1.0, 1.0)) / u_size).xy + texture2D(u_deriv, (position + vec2(-1.0, 1.0)) / u_size).xy + texture2D(u_deriv, (position + vec2(-1.0, -1.0)) / u_size).xy + texture2D(u_deriv, (position + vec2(1.0, -1.0)) / u_size).xy ); vec2 delta = vec2( D_a * laplacian.x - value.x * value.y * value.y + u_F * (1.0 - value.x), D_b * laplacian.y + value.x * value.y * value.y - (u_K + u_F) * value.y ); gl_FragColor = vec4(texture2D(u_f, position / u_size).xy + delta * u_Dt, 0.0, 1.0); }`)33:[] 2f:{"title":"Gray–Scott Reaction Diffusion","cells":[{"id":0,"value":"md`# Gray–Scott Reaction Diffusion`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1016,"value":"md`\nThis fork of [@mbostock's notebook](https://beta.observablehq.com/@mbostock/gray-scott-reaction-diffusion) adds two things. First, it adds a [nine-point finite difference stencil](https://en.wikipedia.org/wiki/Discrete_Laplace_operator#Implementation_via_operator_discretization), and second it adds optional second order midpoint time integration.\n\nI haven't run the stability analysis, but the nine-point stencil seems to increase the time step stability limit by about ~60%, meaning you can run it a fair amount faster.\n\nI thought perhaps second order time integration would similarly increase the stability limit, but the only conclusion is that second order time integration is twice as much more work and perhaps a bit more accurate, but actually makes the stability somewhat worse.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":80,"value":"canvas = DOM.canvas(width, 720)","pinned":false,"hidden":false,"mode":"ojs"},{"id":853,"value":"md`Ref. [@pmneila](https://pmneila.github.io/jsexp/grayscott/), [@robinhouston](https://bl.ocks.org/robinhouston/ed597847175cf692ecce).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":539,"value":"viewof FK = {\n const options = [\n {name: \"Default\", feed: 0.037, kill: 0.06},\n {name: \"Fast\", feed: 0.018, kill: 0.051},\n {name: \"First\", feed: 0.0545, kill: 0.062},\n {name: \"Solitons\", feed: 0.03, kill: 0.062},\n {name: \"Pulsating solitons\", feed: 0.025, kill: 0.06},\n {name: \"Worms\", feed: 0.078, kill: 0.061},\n {name: \"Mazes\", feed: 0.029, kill: 0.057},\n {name: \"Holes\", feed: 0.039, kill: 0.058},\n {name: \"Chaos\", feed: 0.026, kill: 0.051},\n {name: \"Chaos and holes\", feed: 0.034, kill: 0.056},\n {name: \"Moving spots\", feed: 0.014, kill: 0.054},\n {name: \"Spots and loops\", feed: 0.018, kill: 0.051},\n {name: \"Waves\", feed: 0.014, kill: 0.045},\n {name: \"The U-Skate World\", feed: 0.062, kill: 0.06093}\n ];\n const form = html`

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":957,"value":"function step (tex1, tex2, dest, u_Dt, dt) {\n gl.activeTexture(gl.TEXTURE0);\n gl.bindTexture(gl.TEXTURE_2D, tex1);\n\n gl.activeTexture(gl.TEXTURE1);\n gl.bindTexture(gl.TEXTURE_2D, tex2);\n\n gl.bindFramebuffer(gl.FRAMEBUFFER, dest);\n gl.uniform1f(u_Dt, dt);\n gl.drawArrays(gl.TRIANGLES, 0, 3);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":144,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":261,"value":"md`## Textures`","pinned":false,"hidden":false,"mode":"ojs"},{"id":102,"value":"W = 512","pinned":false,"hidden":false,"mode":"ojs"},{"id":653,"value":"H = W","pinned":false,"hidden":false,"mode":"ojs"},{"id":110,"value":"seed = {\n restart;\n const a = new Float32Array(4 * W * H);\n\n for (let y = 0; y < H; ++y) {\n for (let x = 0; x < W; ++x) {\n const i = W * y + x << 2;\n a[i + 0] = 0.9;\n }\n }\n\n function dot(cx, cy, r = 5) {\n const r2 = r ** 2;\n for (let y = cy - r; y < cy + r; ++y) {\n for (let x = cx - r; x < cx + r; ++x) {\n if ((x - cx) ** 2 + (y - cy) ** 2 < r2) {\n const i = W * y + x << 2;\n a[i + 1] = 0.9;\n }\n }\n }\n }\n\n for (let i = 0; i < 20; ++i) {\n dot(W * Math.random(), H * Math.random());\n }\n\n return a;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":117,"value":"function newTexture(source) {\n const texture = gl.createTexture();\n gl.bindTexture(gl.TEXTURE_2D, texture);\n gl.texImage2D(gl.TEXTURE_2D, 0, gl.RGBA, W, H, 0, gl.RGBA, gl.FLOAT, source);\n gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);\n gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);\n gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.REPEAT);\n gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.REPEAT);\n return texture;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":121,"value":"function newFramebuffer(texture) {\n const framebuffer = gl.createFramebuffer();\n gl.bindFramebuffer(gl.FRAMEBUFFER, framebuffer);\n gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, texture, 0);\n return framebuffer;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":257,"value":"md`## Shaders`","pinned":false,"hidden":false,"mode":"ojs"},{"id":93,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":97,"value":"renderShader = newShader(gl.FRAGMENT_SHADER, `\nprecision mediump float;\nuniform sampler2D u_f;\nuniform vec2 u_size;\n\nvoid main() {\n vec2 x = texture2D(u_f, gl_FragCoord.xy / u_size).xy;\n x.y = 1.0 - x.y;\n gl_FragColor = vec4(x, 1.0, 1.0);\n}`)","pinned":true,"hidden":false,"mode":"ojs"},{"id":89,"value":"vertexShader = newShader(gl.VERTEX_SHADER, `\nattribute vec2 a_position;\n\nvoid main() {\n gl_Position = vec4(a_position, 0.0, 1.0);\n}`)","pinned":true,"hidden":false,"mode":"ojs"},{"id":133,"value":"function newShader(type, source) {\n const shader = gl.createShader(type);\n gl.shaderSource(shader, source);\n gl.compileShader(shader);\n if (!gl.getShaderParameter(shader, gl.COMPILE_STATUS)) throw new Error(\"cannot compile shader\");\n return shader;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":254,"value":"md`## Programs`","pinned":false,"hidden":false,"mode":"ojs"},{"id":192,"value":"renderProgram = newProgram(vertexShader, renderShader)","pinned":true,"hidden":false,"mode":"ojs"},{"id":196,"value":"timestepProgram = newProgram(vertexShader, timestepShader)","pinned":true,"hidden":false,"mode":"ojs"},{"id":562,"value":"{\n gl.useProgram(timestepProgram);\n gl.uniform1f(gl.getUniformLocation(timestepProgram, \"u_F\"), FK.feed);\n gl.uniform1f(gl.getUniformLocation(timestepProgram, \"u_K\"), FK.kill);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":129,"value":"function newProgram(vertexShader, fragmentShader) {\n const program = gl.createProgram();\n gl.attachShader(program, vertexShader);\n gl.attachShader(program, fragmentShader);\n gl.linkProgram(program);\n if (!gl.getProgramParameter(program, gl.LINK_STATUS)) throw new Error(\"cannot link program\");\n const a_position = gl.getAttribLocation(program, \"a_position\");\n const u_size = gl.getUniformLocation(program, \"u_size\");\n const u_Dt = gl.getUniformLocation(program, \"u_Dt\");\n const u_deriv = gl.getUniformLocation(program, \"u_deriv\");\n gl.useProgram(program);\n gl.bindBuffer(gl.ARRAY_BUFFER, gl.createBuffer());\n gl.bufferData(gl.ARRAY_BUFFER, new Float32Array([-4, -4, 0, 4, 4, -4]), gl.STATIC_DRAW);\n gl.enableVertexAttribArray(a_position);\n gl.vertexAttribPointer(a_position, 2, gl.FLOAT, false, 0, 0);\n gl.uniform2f(u_size, W, H);\n gl.uniform1f(u_Dt, 1.0);\n gl.uniform1i(u_deriv, 1);\n return program;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":282,"value":"gl = {\n const gl = canvas.getContext(\"webgl\", {depth: false});\n if (!gl) throw new Error(\"WebGL unsupported\");\n if (!gl.getExtension(\"OES_texture_float\")) throw new Error(\"OES_texture_float unsupported\");\n if (gl.getParameter(gl.MAX_TEXTURE_SIZE) < W) throw new Error(\"big textures unsupported\");\n return gl;\n}","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}