1d:{"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":"$1f","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":"$20","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":"$21","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":"$22","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":"$23","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":"$24","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":"$25","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":"$W26","stdlib":"1"}