31:Tf70, This notebook performs a more or less justified but not particularly accurate two-dimensional simulation of ${particleTextureRadius * particleTextureRadius} bodies interacting gravitationally on a periodic ${tex`${gridTextureRadius} \times ${gridTextureRadius}`} grid. It visualizes the potential computed under **two-dimensional** Newtonian gravity, which means two things. First, there's nothing fancy like relativity or spacetime involved, and second that [Poisson's equation](https://en.wikipedia.org/wiki/Poisson%27s_equation) which describes gravitational potential works differently in two dimensions than in three so that gravity won't work quite like it does in our three-dimensional universe. Typically when we simulate gravity in a flat plane we use the three dimensional inverse square law and restrict motion to a plane. This notebook is different and computes gravity that really only lives in two dimensions. In this flat world, energy must still be conserved, of course, but we don't get an inverse square law; things don't pull on each other with quite the right strength and don't fall into nice elliptical orbits. It's still a great and fun exercise though! Pairwise interactions of ${tex`n`} bodies scale with ${tex`O(n^2)`} so that adding up pairwise forces and integrating in time gets out of hand very quickly as ${tex`n`} grows. You can get around that with trees and the approximate [Barnes-Hut Method](https://en.wikipedia.org/wiki/Barnes%E2%80%93Hut_simulation) or the significantly more advanced [Fast Multipole Method](https://en.wikipedia.org/wiki/Fast_multipole_method). This notebook takes a different approach. An obsolete approach, to be clear. It uses the [Particle Mesh (PM) method](https://en.wikipedia.org/wiki/Particle_Mesh) to compute forces. The core principle in this simulation is the [Poisson equation for Newtonian gravity](https://en.wikipedia.org/wiki/Poisson%27s_equation#Newtonian_gravity), ${tex.block`\nabla^2 \phi = 4 \pi G \rho,`} where ${tex`\phi`} is the gravitational potential, ${tex`G`} is the gravitational constant, and ${tex`\rho`} is the local mass density. Instead of computing forces, we aggregate the total mass density in each grid cell, solve for the gravitational potential ${tex`\phi`}, and compute for each body the resulting gravitational force ${tex`\mathbf{g}`}, ${tex.block`\mathbf{g} = -\nabla \phi.`} *Edit: After some googling, I learned this method has a name. It's called the *[Particle Mesh Method](https://en.wikipedia.org/wiki/Particle_Mesh)*, appropriately enough. From the article,* > PM is considered an obsolete method as it does not model close interaction between particles well. It has been supplanted by the [Particle-Particle Particle-Mesh method](https://en.wikipedia.org/wiki/P3M), which uses a straight particle-particle sum between nearby particles in addition to the PM calculation. The rough steps of this method as implemented on the GPU are 1. Compute the density by drawing each body as a point on a texture. Use additive blending so that the output is a texture where each pixel counts the number of bodies in that cell. 2. Solve Poisson's equation for the gravitational potential. Solving this equation in the frequency domain means using the Fast Fourier Transform (see below for details). 3. Update particle positions and velocities where the gradient of the potential from step 2 is the acceleration felt by each particle. This method tends to look a little wild since hundreds of thousands of bodies fly in every direction, so we perform one additional step and simulate Galileian-invariant (but physically questionable) viscosity. To do this, we also compute a weighted average of the velocity as part of step 1, then in step 3 add to each body a force proportional to the difference between its velocity and the cell-average velocity. This pulls each mass toward the average velocity of that cell, though it does not conserve energy.32:T908,onFrame = { restart; camera.taint(); const onFrame = ({ tick }) => { if (run) { configureUniforms({ dt }, () => { // Step 1. Accumulate particles onto a texture densityFBO.use(() => { regl.clear({ color: [0, 0, 0, 1] }); accumulateParticles({ src: particleFBO[0] }); }); configureMap(() => { // Step 2. Solve Poisson's equation for the potential performFFT(forwardFFT); phiFBO[1].use(() => computePotential({ src: phiFBO[0] })); performFFT(inverseFFT); // Step 3. Integrate in time particleFBO[1].use(() => { integrate({ dt, viscosity, phi: phiFBO[0], src: particleFBO[0], density: densityFBO, globalDecay: decayTime }); }); swap(particleFBO, 0, 1); // Step 4. Rasterize the bodies onto a texture we'll loop and draw // to the screen loopFBO.use(() => { regl.clear({ color: view3d ? [0, 0, 0, 1] : backgroundColorRGBA }); drawParticlesToTexture({ src: particleFBO[0], density: densityFBO, viscosityColor: viscosityColorRGB, color: particleColorRGB, opacity, viscosity }); }); }); }); } // Step 5. either render a three dimensional view or just transfer // the loop framebuffer to the screen camera({ rotationCenter: camera.params.center }, state => { if (!run && !state.dirty) return; if (view3d) { regl.clear({ color: [0, 0, 0, 1], depth: 1 }); drawMesh( modelMatrices.map(m => ({ src: phiFBO[0], particleTexture: loopFBO, model: m, gamma: gamma })) ); } else { configureMap(() => { blit({ src: loopFBO, gamma }); }); } }); if (run) { if (saveFrames && tick % frameSaveIncrement === 0) { downloadURI( canvas.toDataURL(), `frame-${tick.toString().padStart(8, '0')}.png` ); } } }; const frame = regl.frame(onFrame); invalidation.then(frame.cancel); }33:T45a,accumulateParticles = regl({ vert: ` precision highp float; attribute vec2 aXY; varying vec2 vVelocity; uniform sampler2D uSrc; void main () { vec4 state = texture2D(uSrc, aXY); vVelocity = state.zw; gl_Position = vec4(state.xy * 2.0 - 1.0, 0, 1); gl_PointSize = 1.0; } `, frag: ` precision highp float; varying vec2 vVelocity; void main () { // We accumulate mass = 1 for each particle into the x channel, // but we also accumulate the velocity into the last two channels. // This will allow us to compute an average velocity in each cell // in order to apply a sort of viscosity which pulls each point in // the cell toward the average velocity of its neighbors. gl_FragColor = vec4(vVelocity, 1, 0); }`, attributes: { aXY: particleTextureLookup }, uniforms: { uSrc: regl.prop('src') }, blend: { enable: true, func: { srcRGB: 1, srcAlpha: 1, dstRGB: 1, dstAlpha: 1 } }, depth: { enable: false }, count: particleTextureRadius * particleTextureRadius, primitive: 'points' })34:T1aef,md`## Poisson's equation for the gravitational potential This is the step which perhaps requires the most explanation. Given the density ${tex`\rho`} computed when we accumulate particles above, the equation we seek to solve for the potential ${tex`\phi`} in ${tex.block`\nabla^2 \phi = 4 \pi G \rho.`} A commonly used strategy for solving [Poisson's equation](https://en.wikipedia.org/wiki/Poisson%27s_equation) is [Jacobi iteration](https://en.wikipedia.org/wiki/Jacobi_method), which you can think of roughly like repeatedly averaging grid points with their neighbors with addition of the source term ${tex`4 \pi G \rho`}. It's like a blurring or diffusion process, but since there's no time component to this equation, we have to iterate until we reach steady state. The process is excruciatingly slow! Since Newtonian gravity travels instantaneously, it's rather annoying to use a numerical method that travels one grid cell per iteration. Convergence typically takes hundreds or thousands of iterations. The _ideal_ technique is a [multigrid method](https://en.wikipedia.org/wiki/Multigrid_method) which (remarkably) requires ${tex`O(n)`} operations, where n is the number of grid cells. There's lots of bookkeeping involved though, so in this notebook we'll use the [Fast Fourier Transform](https://en.wikipedia.org/wiki/Fast_Fourier_transform) to solve Poisson's equation in the frequency domain. This method is almost as good, requiring ${tex`O(n \log n)`} iterations. Solving Poisson's equation in the frequency domain requires a bit of reorganization—in addition to the tricks it requires to solve it in WebGL! We start by representing the per-grid-cell density we've computed in terms of frequency components. You can think each frequency component as ripples in a given direction, similar to [JPEG components](https://en.wikipedia.org/wiki/JPEG#Encoding) if you're familiar with that. Each component is just a sine wave pointing some direction in two dimensions. We talk about them in terms of wavenumbers ${tex`k`} which are related to wavelength ${tex`\lambda`} as ${tex`k = \frac{2\pi}{\lambda}`}. Omitting some uninteresting details of integration, we consider only a single frequency component at a time and substitute for gravitational potential ${tex`\phi`} and density ${tex`\rho`} ${tex.block`\begin{aligned} \rho &= \hat{\rho} e^{i(k_x x + k_y y)} \\ \phi &= \hat{\phi} e^{i(k_x x + k_y y)} \end{aligned}`} where ${tex`k_x`} and ${tex`k_y`} are x and y wavenumbers. You can think of ${tex`\hat{\phi}`} and ${tex`\hat{\rho}`} as complex numbers representing how much of each frequency component we have. At last we can do some math! We write the [Laplacian operator](https://en.wikipedia.org/wiki/Laplace_operator) ${tex`\nabla^2`} as ${tex.block`\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}.`} Then we have for our equation ${tex.block`\begin{aligned} \nabla^2 \left(\hat{\phi} e^{i(k_x x + k_y y)}\right) &= -(k_x^2 + k_y^2)\hat{\phi}e^{i(k_x x + k_y y)} \\ &= 4 \pi G \hat{\rho} e^{i(k_x x + k_y y)} \end{aligned}`} Note that ${tex`\hat{\phi}`} and ${tex`\hat{\rho}`} are unaffected by the partial derivatives since they're constant complex numbers for a given wavenumber component. We cancel the complex exponential and write ${tex.block`-(k_x^2 + k_y^2)\hat{\phi} = 4 \pi G \hat{\rho},`} finally dividing through by the wavenumbers to write ${tex.block`\hat{\phi} = -\frac{4 \pi G \hat{\rho}}{k_x^2 + k_y^2}.`} If we started with the Fourier transform of the density, this is our answer as a Fourier transform of the potential. Once we've computed the the gravitational potential ${tex`\hat{\phi}`} in the frequency domain, we perform the inverse Fast Fourier Transform and obtain the potential ${tex`\phi`}. There's just one additional step. The gravitational potential is nice for visualizing potential wells, but for updating the particle positions we require its gradient ${tex.block`\mathbf{g} = -\nabla \phi,`}where ${tex`\mathbf{g}`} is the acceleration of a body due to gravity. We compute this similarly, writing the gradient as ${tex.block`\nabla \phi = \mathbf{i} \frac{\partial \phi}{\partial x} + \mathbf{j} \frac{\partial \phi}{\partial y}.`} Considering the components one at a time, we have in the frequency domain for the x-component of the gradient ${tex.block`\frac{\partial}{\partial x}\left(\hat{\phi} e^{i(k_x x + k_y y)}\right) = i k_x \hat{\phi}e^{i(k_x x + k_y y)}`} Cancelling the exponentials, we're left with ${tex.block`\frac{\partial \hat{\phi}}{\partial x} = i k_x \hat{\phi}.`} Similarly for the y-component, ${tex.block`\frac{\partial \hat{\phi}}{\partial y} = i k_y \hat{\phi}.`} To make a long story short, differentiation in the spatial domain is equivalent to multiplication by wavenumber and a phase shift in the frequency domain. Once we've computed ${tex`\hat{\phi}`}, we can directly compute the force on a given body! The actual computation of this gradient is where things get fun. Since we perform all of this computation using WebGL textures, we have four channels to work with, which for complex numbers ${tex`a`} and ${tex`b`} is encoded in RGBA as ${tex`(\mathrm{Re}(a), \mathrm{Im}(a), \mathrm{Re}(b), \mathrm{Im}(b))`}. Instead of four real numbers, we'll write this as a tuple of two complex numbers, ${tex`(a, b)`}. Our input is the density in the spatial comain, ${tex`(\rho, 0)`}. Once Fourier transformed, we have its frequency-domain representation, ${tex`(\hat{\rho}, 0)`}. We then compute the Fourier transform of the potential ${tex`\hat{\phi}`} via our equation ${tex`\hat{\phi} = -\frac{4 \pi G }{k_x^2 + k_y^2}\hat{\rho}`} and place this result in the first two components, yielding ${tex`(\hat{\phi}, 0)`}. The last two channels are left empty, so we use them during the inverse Fourier transform to compute the gradient of the potential! If we hope to retrieve them in the last two of four channels, then we seek the sum ${tex.block`\frac{\partial \phi}{\partial x} + i \frac{\partial \phi}{\partial y}`} which we write in the frequency domain as ${tex.block`i k_x \hat{\phi} + i (i k_y \hat{\phi}),`} or ${tex.block`i k_x \hat{\phi} - k_y \hat{\phi}.`} Thus our final value before inverse Fourier transforming is the tuple of two complex numbers, ${tex.block`(\hat{\phi}, i k_x \hat{\phi} - k_y \hat{\phi}).`} After inverse Fourier transforming, we recover our result, ${tex.block`\left(\phi, \frac{\partial \phi}{\partial x} + i\frac{\partial \phi}{\partial y}\right).`} If we break this tuple of two complex numbers into its four components as are stored in the texture, then we write this as ${tex.block`\left(\phi, 0, \frac{\partial \phi}{\partial x}, \frac{\partial \phi}{\partial y}\right)`} We plot the gravitational potential ${tex`\phi`} as a height field and sample the gradient with linear interpolation to compute the force acting on each body. `35:T6b8,integrate = regl({ frag: ` precision highp float; varying vec2 vUV; uniform sampler2D uSrc, uDensity, uPhi; uniform float dx, dxInv, dt, uGlobalDecay, uViscosity; vec4 deriv (vec4 y) { vec2 phiGradient = texture2D(uPhi, y.xy).zw; vec4 dydt = vec4(y.zw, phiGradient); // This step is a bit half-baked. I intended to compute the average // nearby velocity and move particles toward it, but it just doesn't // seem to have quite the dramatic effect I hoped vec3 density = texture2D(uDensity, y.xy).zxy; //vec3 densityN = texture2D(uDensity, y.xy + vec2(0, dx)).xzw; //vec3 densityS = texture2D(uDensity, y.xy - vec2(0, dx)).xzw; //vec3 densityE = texture2D(uDensity, y.xy + vec2(dx, 0)).xzw; //vec3 densityW = texture2D(uDensity, y.xy - vec2(dx, 0)).xzw; //vec3 avgDensity = 0.5 * density + 0.125 * (densityN + densityS + densityE + densityW); vec3 avgDensity = density; vec2 vMean = avgDensity.yz / max(avgDensity.x, 1.0); vec2 vdiff = vMean - y.zw; dydt.zw += vdiff * (1.0 - exp(-avgDensity.x * uViscosity)); dydt.zw -= uGlobalDecay * y.zw; return dydt; } void main () { vec4 y = texture2D(uSrc, vUV); // Midpoint integration // Predict: vec4 yHalf = y + deriv(y) * dt * 0.5; yHalf.xy = mod(yHalf.xy, vec2(1, 1)); // Correct: vec4 yNext = y + deriv(yHalf) * dt; yNext.xy = mod(yNext.xy, vec2(1, 1)); gl_FragColor = yNext; }`, uniforms: { uSrc: regl.prop('src'), uDensity: regl.prop('density'), uPhi: regl.prop('phi'), uViscosity: regl.prop('viscosity'), uGlobalDecay: regl.prop('globalDecay') } })36:T577,drawParticlesToTexture = regl({ vert: ` precision highp float; attribute vec2 aXY; varying vec2 vVelocity; varying float vFriction; uniform float uViscosity, uViscosityOpacity; uniform sampler2D uSrc, uDensity; void main () { vec4 state = texture2D(uSrc, aXY); vec3 density = texture2D(uDensity, state.xy).zxy; vec2 vavg = density.yz / max(density.x, 1.0); vec2 dv = state.zw - vavg; vFriction = pow(length(dv), 2.0) * 0.2 * uViscosityOpacity; gl_Position = vec4(state.xy * 2.0 - 1.0, 0, 1); gl_PointSize = 1.0; } `, frag: ` precision highp float; uniform float uAlpha; uniform vec3 uColor, uViscosityColor; varying float vFriction; void main () { gl_FragColor = vec4((uColor + vFriction * uViscosityColor) * uAlpha, 1.0); }`, attributes: { aXY: particleTextureLookup }, uniforms: { uViscosity: regl.prop('viscosity'), uViscosityColor: regl.prop('viscosityColor'), uColor: regl.prop('color'), uSrc: regl.prop('src'), uDensity: regl.prop('density'), uAlpha: regl.prop('opacity'), uViscosityOpacity: () => Math.pow(gridTextureRadius / 256, 2) }, blend: { enable: true, func: { srcRGB: 1, srcAlpha: 1, dstRGB: 1, dstAlpha: 1 } }, depth: { enable: false }, count: particleTextureRadius * particleTextureRadius, primitive: 'points' })37:Tc40,drawMesh = { function createTextureLookupTable(w, h, stride) { stride = stride || 2; var n = w * h * stride; var out = new Float32Array(n); for (var i = 0, iStride = 0; iStride < n; i++, iStride += stride) { out[iStride] = ((i % w) + 0.5) / (w - 1); out[iStride + 1] = (((i / w) | 0) + 0.5) / (h - 1); } return out; } var meshElements = []; var n = gridTextureRadius + 1; for (var j = 0; j < n - 1; j++) { for (var i = 0; i < n - 1; i++) { meshElements.push( [i + n * j, i + 1 + n * j, i + 1 + n * (j + 1)], [i + n * j, i + 1 + n * (j + 1), i + n * (j + 1)] ); } } return regl({ vert: ` precision highp float; attribute vec2 aXY; uniform mat4 uProjectionView, uModel; uniform float uVerticalScale; uniform sampler2D uSrc; varying vec3 vXYZ; varying vec2 vUV; varying float r, vY; void main () { vY = texture2D(uSrc, aXY).x; vUV = aXY; float yMean = 0.25 * ( texture2D(uSrc, vec2(0.25, 0.25)).x + texture2D(uSrc, vec2(0.25, 0.75)).x + texture2D(uSrc, vec2(0.75, 0.25)).x + texture2D(uSrc, vec2(0.75, 0.75)).x ); vY -= yMean; vY *= uVerticalScale; vXYZ = vec3(aXY.x, vY * 0.03, aXY.y); vec4 pos = uModel * vec4(vXYZ * 2.0 - 1.0, 1); r = length(pos.xz); gl_Position = uProjectionView * pos; }`, frag: ` #extension GL_OES_standard_derivatives : enable precision highp float; uniform sampler2D uParticleTex; float gridFactor (float parameter, float width, float feather) { float w1 = width - feather * 0.5; float d = fwidth(parameter); float looped = 0.5 - abs(mod(parameter, 1.0) - 0.5); return smoothstep(d * w1, d * (w1 + feather), looped); } float gridFactor (vec2 parameter, float width, float feather) { float w1 = width - feather * 0.5; vec2 d = fwidth(parameter); vec2 looped = 0.5 - abs(mod(parameter, 1.0) - 0.5); vec2 a2 = smoothstep(d * w1, d * (w1 + feather), looped); return min(a2.x, a2.y); } varying vec3 vXYZ; varying vec2 vUV; varying float r, vY; uniform float uGamma; void main () { vec3 particle = texture2D(uParticleTex, vUV).rgb; particle.r = pow(particle.r, uGamma); particle.g = pow(particle.g, uGamma); particle.b = pow(particle.b, uGamma); gl_FragColor = vec4(smoothstep(3.0, 1.5, r) * ( vec3( 0.10 + 0.10 * (1.0 - gridFactor(vXYZ.xz * 4.0, 0.5, 1.0)) + 0.07 * (1.0 - gridFactor(vXYZ.xz * 40.0, 0.5, 1.0)) + 0.85 * vec3(smoothstep(2.0, -40.0, vY)) ) + 0.6 * particle), 1); }`, attributes: { aXY: createTextureLookupTable(n, n) }, uniforms: { uVerticalScale: Math.pow(gridTextureRadius / 256, 2), uSrc: regl.prop('src'), uParticleTex: regl.prop('particleTexture'), uModel: regl.prop('model'), uGamma: (ctx, props) => 1 / props.gamma }, depth: { enable: true }, primitive: 'triangles', elements: meshElements }); }38:[] 2f:{"title":"2D (Non-physical) N-body Gravity with Poisson's Equation","cells":[{"id":3084,"value":"# 2D (Non-physical) N-body Gravity with Poisson's Equation","pinned":false,"hidden":false,"mode":"md"},{"id":3090,"value":"*[Read this notebook ➡️](https://rreusser.github.io/notebooks/particle-mesh-gravity/) on my personal website via [Observable Notebook Kit](https://observablehq.com/notebook-kit/kit).*","pinned":false,"hidden":false,"mode":"md"},{"id":930,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":2199,"value":"md`\n## Known issues\n- Viscosity is a bit half-baked and not quite scale-independent\n- Half floating point precision has some major precision issues which could be avoided with more careful arithmetic and scaling. ${\n !canWriteToFloatFBO\n ? `Warning: your device may support floating point textures but does not support writing to them. Half-floating point precision will be used instead, even though this computation has major precision issues with half-float precision.`\n : ''\n}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2691,"value":"md`Comments? Questions? [Let me know](https://twitter.com/rickyreusser/status/1189746847427907584).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":17,"value":"canvas = {\n // Hold onto the canvas and resize it\n var canvas = this || document.createElement('canvas');\n canvas.width = shape.width;\n canvas.height = shape.height;\n return canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":374,"value":"viewof restart = button({ value: 'Restart' })","pinned":false,"hidden":false,"mode":"ojs"},{"id":314,"value":"viewof run = checkbox({ options: ['Simulate'], value: 'Simulate' })","pinned":false,"hidden":false,"mode":"ojs"},{"id":1479,"value":"viewof view3d = checkbox({\n options: ['View potential in third dimension'],\n value: 'View potential in third dimension'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":625,"value":"viewof particleTextureRadius = {\n const form = html``;\n form.i.oninput = () => {\n var radius = Math.pow(2, form.i.valueAsNumber);\n var count = radius * radius;\n form.value = radius;\n form.o.value = `${count}`;\n };\n form.i.oninput();\n return form;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":658,"value":"viewof gridTextureRadius = {\n const form = html``;\n form.i.oninput = () => {\n var radius = Math.pow(2, form.i.valueAsNumber);\n form.value = radius;\n form.o.value = `${radius}`;\n };\n form.i.oninput();\n return form;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":671,"value":"viewof screenTextureRadius = {\n const form = html``;\n form.i.oninput = () => {\n var radius = Math.round(Math.pow(2, form.i.valueAsNumber));\n form.value = radius;\n form.o.value = `${radius}`;\n };\n form.i.oninput();\n return form;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":424,"value":"viewof viscosity = slider({\n min: 0.0,\n max: 0.5,\n step: 0.01,\n value: 0.0,\n description: 'Viscosity (not a very strong effect)'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":447,"value":"viewof decayTime = slider({\n min: 0.0,\n max: 0.5,\n step: 0.01,\n value: 0.0,\n description: 'Global velocity damping'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1031,"value":"viewof alpha = slider({\n min: 0.0,\n max: 1,\n step: 0.01,\n value: 0.5,\n description: 'alpha'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":849,"value":"viewof gamma = slider({\n min: 0.01,\n max: 3,\n step: 0.01,\n value: 2.2,\n description: 'gamma'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":995,"value":"viewof dt = {\n const form = html``;\n form.i.oninput = () => {\n form.value = form.i.valueAsNumber * 1e-4;\n form.o.value = `${form.value}`;\n };\n form.i.oninput();\n return form;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":813,"value":"viewof particleColor = color({\n value: '#574ffe',\n name: 'color',\n description: 'Particle color'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":871,"value":"viewof viscosityColor = color({\n value: '#ff6650',\n name: 'color',\n description: 'Viscosity color'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2027,"value":"md`## Frame output\n\nBe careful here. It'll slam your browser with downloads, but if you _do_ take the plunge, my suggested command is:\n\n\\`\\`\\`\nffmpeg -r 30 -pattern_type glob -i \"frame-*.png\" -vcodec libx264 -crf 25 -pix_fmt yuv420p -crf 23 output.mp4\n\\`\\`\\`\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1929,"value":"viewof saveFrames = checkbox({\n options: ['Save frames'],\n description:\n 'Warning! This starts downloading *lots* of files immediately. You should disable \"ask to save each file\" before checking this. Srsly. Be careful.'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1934,"value":"viewof frameSaveIncrement = slider({\n min: 1,\n max: 100,\n step: 1,\n value: 5,\n description: 'frame save skip'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2432,"value":"md`## Parameters`","pinned":false,"hidden":false,"mode":"ojs"},{"id":180,"value":"G = 1 * Math.pow(512 / particleTextureRadius, 2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":622,"value":"opacity = Math.min(\n 4,\n 0.2 *\n alpha *\n Math.pow(512 / particleTextureRadius, 2) *\n Math.pow(screenTextureRadius / 512, 2)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2350,"value":"particleColorRGB = hexRgbToFloat(particleColor).map(c =>\n Math.pow(c, gamma)\n )","pinned":false,"hidden":false,"mode":"ojs"},{"id":2352,"value":"viscosityColorRGB = hexRgbToFloat(viscosityColor).map(c =>\n Math.pow(c, gamma)\n )","pinned":false,"hidden":false,"mode":"ojs"},{"id":2346,"value":"backgroundColorRGBA = [0.04, 0.04, 0.04]\n .map(c => Math.pow(c, view3d ? 1.0 : gamma))\n .concat([1])","pinned":false,"hidden":false,"mode":"ojs"},{"id":1969,"value":"md`## WebGL setup\n\nWe perform all computation using the [regl](https://github.com/regl-project/regl) WebGL library.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2435,"value":"viewof dataType = select({\n options: ['float', 'half float'],\n value: canWriteToFloatFBO ? 'float' : 'half float',\n title: 'Data type',\n description:\n 'Framebuffer precision. Falls back to half float if floating point is not available, though the quality of results does suffer.'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"createREGL = require('regl')","pinned":false,"hidden":false,"mode":"ojs"},{"id":27,"value":"regl = {\n // Hold onto the instance and resize it. If we change extensions,\n // we'll have to comment out `this ||` in order to reinstantiate.\n // This just prevents excessive context creation when the width\n // changes.\n var regl =\n this ||\n createREGL({\n extensions: [\n 'OES_texture_float',\n 'OES_texture_float_linear',\n 'OES_element_index_uint',\n 'OES_standard_derivatives'\n ],\n optionalExtensions: [\n 'OES_texture_half_float',\n 'OES_texture_half_float_linear'\n ],\n canvas: canvas,\n pixelRatio: devicePixelRatio,\n attributes: { antialias: true, preserveDrawingBuffer: true }\n });\n return regl;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":66,"value":"configureUniforms = regl({\n uniforms: {\n resolution: ctx => [ctx.viewportWidth, ctx.viewportHeight],\n dx: ctx => 1 / gridTextureRadius,\n dxInv: ctx => gridTextureRadius,\n dt: regl.prop('dt')\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":58,"value":"configureMap = regl({\n vert: `\n precision highp float;\n attribute vec2 aXY;\n varying vec2 vUV;\n void main () {\n vUV = aXY * 0.5 + 0.5;\n gl_Position = vec4(aXY, 0, 1);\n }\n `,\n attributes: {\n aXY: [-4, -4, 4, -4, 0, 4]\n },\n depth: { enable: false },\n primitive: 'triangle',\n count: 3\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2169,"value":"configureCamera = regl({\n uniforms: {\n uProjectionView: regl.prop('projectionView')\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":750,"value":"blit = regl({\n frag: `\n precision highp float;\n varying vec2 vUV;\n uniform float uGamma, uAspect;\n uniform sampler2D uSrc;\n uniform vec2 uScale;\n void main () {\n vec3 color = texture2D(uSrc, mod(0.5 + (vUV - 0.5) * uScale, vec2(1))).xyz;\n gl_FragColor = vec4(\n pow(color.r, uGamma),\n pow(color.g, uGamma),\n pow(color.b, uGamma),\n 1);\n }`,\n uniforms: {\n uGamma: (ctx, props) => 1 / props.gamma,\n uSrc: regl.prop('src'),\n uAspect: ctx => ctx.framebufferWidth / ctx.framebufferHeight,\n uScale: ctx => [\n ctx.viewportWidth / screenTextureRadius,\n ctx.viewportHeight / screenTextureRadius\n ]\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2152,"value":"md`## Iteration\n\nWe define our main iteration loop below.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":114,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":2336,"value":"initialize = {\n // Trigger reinitialization when 'restart' is clicked\n restart;\n\n particleFBO.forEach(fbo =>\n transferAttrToTexture({\n src: initialParticleStateBuffer,\n dst: fbo\n })\n );\n\n camera.taint();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2000,"value":"md`## Initialization`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2637,"value":"function pdf(x, y) {\n // PDF used for sampling initial positions. x and y in [0, 1], and returned value is in [0, 1].\n return 1.0;\n\n // Lots of little bumps:22\n //return 0.9 + 0.1 * Math.cos(x * Math.PI * 16) * Math.sin(y * Math.PI * 16);\n\n // Circular blob:\n //var r2 = x * x + y * y;\n //if (r2 > 0.2) return 0;\n //return Math.exp(-r2 / Math.pow(0.3, 2));\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2701,"value":"function initialVelocity(out, x, y) {\n out[0] = 0;\n out[1] = 0;\n return;\n\n /*return [\n 0.0 * -Math.cos(x * Math.PI * 16) * Math.cos(y * Math.PI * 16),\n 0.0 * Math.sin(x * Math.PI * 16) * Math.sin(y * Math.PI * 16)\n ];*/\n\n // Roughly circular orbit speed:\n //var r = Math.sqrt(x * x + y * y);\n //var speed =\n // (4.5 - 3.5 * r + 0.5 * r * (1 - r)) * Math.pow(512 / gridTextureRadius, 1);\n //out[0] = y * speed;\n //out[1] = -x * speed;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":80,"value":"initialParticleState = {\n var v = [0, 0];\n var xn, yn;\n var n = particleTextureRadius * particleTextureRadius;\n var data = new Float32Array(n * 4);\n for (var i = 0, i4 = 0; i < n; i++, i4 += 4) {\n var acceptable = false;\n while (!acceptable) {\n var x = Math.random();\n var y = Math.random();\n xn = x * 2 - 1;\n yn = y * 2 - 1;\n acceptable = Math.random() < pdf(xn, yn);\n }\n data[i4] = x;\n data[i4 + 1] = y;\n\n // Initial velocity\n initialVelocity(v, xn, yn);\n data[i4 + 2] = v[0];\n data[i4 + 3] = v[1];\n }\n return data;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2413,"value":"initialParticleStateBuffer = regl.buffer(initialParticleState)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2325,"value":"transferAttrToTexture = regl({\n vert: `\n precision highp float;\n attribute vec2 aXY;\n attribute vec4 aSrc;\n varying vec4 vSrc;\n void main () {\n vSrc = aSrc;\n gl_Position = vec4(2.0 * aXY - 1.0, 0, 1);\n gl_PointSize = 1.0;\n }`,\n frag: `\n precision highp float;\n varying vec4 vSrc;\n void main () {\n gl_FragColor = vSrc;\n }`,\n attributes: {\n aXY: particleTextureLookup,\n aSrc: regl.prop('src')\n },\n depth: { enable: false },\n framebuffer: regl.prop('dst'),\n primitive: 'points',\n count: particleTextureRadius * particleTextureRadius\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1973,"value":"md`## Accumulate particles\n\nThe first step is to accumulate the particle density onto a texture. We rasterize all bodies as a single-pixel point onto a ${tex`${gridTextureRadius} \\times ${gridTextureRadius}`} texture.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":220,"value":"densityFBO = (this || regl.framebuffer)({\n color: regl.texture({\n radius: gridTextureRadius,\n wrapS: 'repeat',\n wrapT: 'repeat',\n mag: 'linear',\n min: 'linear',\n type: dataType\n })\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":131,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":42,"value":"particleTextureLookup = createTextureLookupTable(\n particleTextureRadius,\n particleTextureRadius\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1977,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":186,"value":"phiFBO = [0, 1, 2, 3].map(i =>\n regl.framebuffer({\n color: regl.texture({\n radius: gridTextureRadius,\n wrapS: 'repeat',\n wrapT: 'repeat',\n mag: i === 0 ? 'linear' : 'nearest',\n min: i === 0 ? 'linear' : 'nearest',\n type: dataType\n })\n })\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1638,"value":"performFFT = createFFTPassCommand(regl)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1626,"value":"forwardFFT = planFFT({\n width: gridTextureRadius,\n height: gridTextureRadius,\n input: densityFBO,\n ping: phiFBO[2],\n pong: phiFBO[3],\n output: phiFBO[0],\n forward: true\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1635,"value":"inverseFFT = planFFT({\n width: gridTextureRadius,\n height: gridTextureRadius,\n input: phiFBO[1],\n ping: phiFBO[2],\n pong: phiFBO[3],\n output: phiFBO[0],\n forward: false\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1641,"value":"computePotential = regl({\n frag: `\n precision highp float;\n varying vec2 vUV;\n uniform sampler2D uSrc;\n uniform vec2 uResolution;\n uniform float dx, uCoeff;\n\n ${glslWavenumber}\n \n void main () {\n vec2 yfft = texture2D(uSrc, vUV).zw;\n\n vec2 k = wavenumber(uResolution, vec2(dx));\n float kx2ky2 = dot(k, k);\n \n // Ensure we don't divide by zero:\n if (kx2ky2 < 1.0e-8) kx2ky2 = 1.0;\n\n vec2 phifft = uCoeff / kx2ky2 * yfft;\n\n vec2 dphidx = vec2(-k.x, k.x) * phifft.yx;\n vec2 dphidy = vec2(-k.y, k.y) * phifft.yx;\n vec2 gradient = dphidx + vec2(-dphidy.y, dphidy.x);\n\n gl_FragColor = vec4(phifft, gradient);\n }\n `,\n uniforms: {\n uSrc: regl.prop('src'),\n uResolution: forwardFFT[0].resolution,\n uCoeff: -4.0 * Math.PI * G\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2006,"value":"md`## Time integration\n\nIntegration in time is straightforward. We use second order midpoint integration, though we fudge it a little and don't fully recompute the potential between the two steps (because that'd be very expensive!).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":78,"value":"particleFBO = [0, 1].map(() => {\n return regl.framebuffer({\n color: regl.texture({\n radius: particleTextureRadius,\n mag: 'nearest',\n min: 'nearest',\n type: dataType\n })\n });\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":94,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":2010,"value":"md`## Output`","pinned":false,"hidden":false,"mode":"ojs"},{"id":745,"value":"loopFBO = (this || regl.framebuffer)({\n color: regl.texture({\n radius: screenTextureRadius,\n type: dataType,\n min: 'linear',\n mag: 'linear',\n wrapS: 'repeat',\n wrapT: 'repeat'\n })\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":395,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":1854,"value":"camera = {\n var camera = createReglCamera(regl, {\n phi: 0.4,\n distance: 4,\n center: [0, -1, 0]\n });\n createInteractions(camera);\n return camera;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1186,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":2263,"value":"md`*TO DO*: instanced rendering instead of multiple draw calls with a different model matrix.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2165,"value":"modelMatrices = [\n [-2, 0, -2],\n [-2, 0, 0],\n [-2, 0, 2],\n [0, 0, -2],\n [0, 0, 0],\n [0, 0, 2],\n [2, 0, -2],\n [2, 0, 0],\n [2, 0, 2]\n].map(v => mat4fromTranslation(mat4create(), v))","pinned":false,"hidden":false,"mode":"ojs"},{"id":2020,"value":"md`## Miscellaneous imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":790,"value":"function swap(obj, key1, key2) {\n key1 = key1 === undefined ? 0 : key1;\n key2 = key2 === undefined ? 1 : key2;\n var tmp = obj[key1];\n obj[key1] = obj[key2];\n obj[key2] = tmp;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1138,"value":"import {\n mat4perspective,\n mat4create,\n mat4lookAt,\n mat4multiply,\n mat4fromTranslation\n} from '@rreusser/gl-mat4'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1850,"value":"import { createReglCamera, createInteractions } from '@rreusser/regl-tools'","pinned":false,"hidden":false,"mode":"ojs"},{"id":834,"value":"import { hexRgbToFloat, downloadURI } from '@rreusser/utils'","pinned":false,"hidden":false,"mode":"ojs"},{"id":43,"value":"import {\n createTextureLookupTable,\n canWriteToFBOOfType\n} from '@rreusser/regl-tools'","pinned":false,"hidden":false,"mode":"ojs"},{"id":311,"value":"import { checkbox, button, slider, color, select } from '@jashkenas/inputs'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1623,"value":"import {\n planFFT,\n createFFTPassCommand,\n glslWavenumber\n} from '@rreusser/glsl-fft'","pinned":false,"hidden":false,"mode":"ojs"},{"id":2449,"value":"canWriteToFloatFBO = canWriteToFBOOfType(regl, 'float')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1168,"value":"shape = {\n //return { width: 800, height: 600 };\n var w = width;\n\n return {\n width: 2 * Math.floor(w / 2),\n height: Math.max(360, 2 * Math.floor(Math.floor((w * 3) / 4) / 2))\n };\n}","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W38","stdlib":"1","titling":"auto"}