Reaction-diffusion system

Reaction-diffusion systems are aptly named as they are primarily composed of two terms, reaction and diffusion. Here we will implement the the Gray-Scott model which models two substances, $U$ and $V$. The change of concentration of there variables is given by these differential equations.

This may look a bit dense so let's break it down. Starting with the first equation, here $u$ is the concentration of $U$, and the equation describes the change of concentration over time. There are three terms in the equation. The first, $D_u \Delta u$ is the diffusion term, the substance spreading naturally spreading out. Here $\Delta$ is the Laplacian and $D_u$ is the diffusion coefficient. Second is the reaction term, $-uv^2$, and the third term, $f(1 - u)$, represents how quickly $U$ is added to the system. $f$ is the feed rate of $u$. The second equation is similar to the first but the third term, $- (f + k)v$, differs. Here $V$ is instead removed from the system and $k$ is the kill rate.

In order to implement these equations on a compute we will need to discretize time and space.

The Laplacian and the discrete Laplacian

The Laplacian, or the Laplace operator, can be thought of as a second derivative for multivariable functions. It can also be interpreted as a deviation from the average of nearby values. I will am using $\Delta$ but $\nabla \cdot \nabla$ and $\nabla^2$ are also commonly used. The Laplacian is equal to the sum of the unmixed second partial derivatives. This is an example of the Laplacian of a function of two variables.

The Laplace operator can be discretized in different ways depending on the context and required stability. In our case we are working with a regular 2D grid. The idea is to approximate the Laplacian at a vertex in the grid by sampling the neighbouring vertices, approximating the second derivatives along each axis and summing them.

This is a formula approximating the second derivatives using central finite differences.

Here $h$ is the distance between adjacent vertices in the grid.
By adding the derivatives we get a formula for the discrete Laplacian in 2D.

Which can be written compactly as a convolution kernel.

Implementation

We will implement the equations in two dimensions as this will allow for complex behaviour while still keeping straightforward to implement. We wil discretize space into a regular grid with periodic boundary conditions. We will use a Webgl texture to store the concentrations, each pixel of the texture corresponding to a vertex in the grid.

As it is not possible for a shader to read from the same texture it is writing to we will need to textures and alternate between them when advancing the state of the system.

class ComputeVariable {
  constructor(textures, framebuffers) {
    this._length = textures.length;
    this._textures = textures;
    this._framebuffers = framebuffers;
    this._head = 0;
  }

  advance() {
    this._head = this._head + 1;
    if (this._head >= this._length) {
      this._head = 0;
    }
  }

  getFramebuffer() {
    return this._framebuffers[this._head];
  }

  getTexture() {
    return this._textures[this._head];
  }
}

We will need two main programs, one to update system state and one program to draw a visualization to the canvas.

Fragment Shader

use EXT_color_buffer_float extension use RG32f texture

#version 300 es

precision highp float;

uniform float u_feed;
uniform float u_kill;
uniform float u_deltaTime;
uniform vec2 u_diffusion;
uniform vec2 u_resolution;
uniform sampler2D u_concentrationTexture;

in vec2 v_textureCoord;
out vec2 fragColor;

vec2 laplacian5(sampler2D tex, vec2 coord, vec2 texResolution, vec2 centralValue) {
    vec2 result = -4.0 * centralValue;
    vec2 h = 1.0 / texResolution;
    vec2 c = coord;
    result += texture( tex, vec2(c.x + h.x, c.y)).xy;
    result += texture( tex, vec2(c.x - h.x, c.y)).xy;
    result += texture( tex, vec2(c.x, c.y + h.y)).xy;
    result += texture( tex, vec2(c.x, c.y - h.y)).xy;
    return result;
}

void main()	{
    vec2 concentration = texture( u_concentrationTexture, v_textureCoord ).xy;
    vec2 lap = laplacian5(u_concentrationTexture, v_textureCoord, u_resolution, concentration);
    float u = concentration.x;
    float v = concentration.y;
    u += (u_diffusion.x * lap.x - u * v * v + u_feed * (1.0 - u)) * u_deltaTime;
    v += (u_diffusion.y * lap.y + u * v * v - (u_kill + u_feed) * v) * u_deltaTime;
    concentration = clamp(vec2(u, v), 0.0, 1.0);
    fragColor = concentration;
}

Webgl basics

Lets start an overview of webgl. I will only cover the general concepts, for more details I recommend the excellent webglfundamentals and webgl2fundamentals. I will be using webgl2 as it is finally on its way to universal adoption with the release of safari 15.

We start with setting up attributes, two clip space coordinates and two texture coordinates. We create four points, one in each corner of the viewport.

const attributeData = new Float32Array([
  -1.0,  1.0,   0.0, 1.0, // x, y,   u, v
  -1.0, -1.0,   0.0, 0.0,
   1.0,  1.0,   1.0, 1.0,
   1.0, -1.0,   1.0, 0.0,
]);

We set up a vertex shader which passes sets the vertex position and passes on the texture coordiate as a varying.

#version 300 es

in vec3 a_position;
in vec2 a_textureCoord;

out vec2 v_textureCoord;

void main(){
  gl_Position = vec4(a_position, 1.0);
  v_textureCoord = a_textureCoord;
}

In effect this creates a quad covering the entire viewport as well as passing on convenient texture coordinates to our fragment shaders which will do our calculations.

This will draw two triangles covering the viewport, in effect calling the fragment shader once for every pixel of the viewport. The vertex shader also passes on texture coordinates as a varying to the fragment shader.

We will use a texture to represent the concentrations of the chemicals. As it is not possible for a fragment shader to read from the same texture it is writing to we'll actually need two textures, alternating which one we are reading from and which we are writing to.

render to a texture, float textures, sampling textures.