Solving the discretized problem using WebGL shaders

If we use a texture to represent the domain and a fragment shader to do the Jacobi relaxation steps, the shader will follow this general pseudocode:

Check if this fragment is a boundary point. If it’s one, return the previous value of this point.

Get the four nearest neighbors’ values.

Return the average of their values.

To flesh out this pseudocode, we need to define a specific representation for the discretized domain. Taking into account that the currently available WebGL versions don’t support floating point textures, we can use 32 bits RGBA fragments and do the following mapping:

Most of the code is straightforward, but doing the multiprecision arithmetic is tricky, as the quantities we are working with behave as floating point numbers in the shaders but are stored as integers. More specifically, the color numbers in the normal range, [0.0, 1.0], are multiplied by 255 and rounded to the nearest byte value when stored at the target texture.

My first idea was to start by reconstructing the floating point numbers for each input value, do the required operations with the floating numbers and convert the floating point numbers to color components that can be reliably stored (without losing precision). This gives us the following pseudocode for the iteration shader:

The reason why we multiply by 255 in place of 256 is that we need val_lo to keep track of the part of val that will be lost when we store it as a color component. As each byte value of a discrete color component will be associated with a range of size 1/255 in its continuous counterpart, we need to use the “low byte” to store the position of the continuous component within that range.

Solving the Laplace's equation using a 32x32 grid. Click the picture to see the live solving process (if your browser supports WebGL).

As can be seen, it has quite low resolution but converges fast. But if we just crank up the number of points, the convergence gets slower:

Incompletely converged solution in a 512x512 grid. Click the picture to see a live version.

How can we reconcile these approaches?

Multigrid

The basic idea behind multigrid methods is to apply the relaxation method on a hierarchy of increasingly finer discretizations of the problem, using in each step the coarse solution obtained in the previous grid as the “starting guess”. In this mode, the long wavelength parts of the solution (those that converge slowly in the finer grids) are obtained in the first coarse iterations, and the last iterations just add the finer parts of the solution (those that converge relatively easily in the finer grids).

Multigrid solution using grids from 8x8 to 512x512. Click the picture to see the live version.

Conclusions

It’s quite viable to use WebGL to do at least basic GPGPU tasks, though it is, in a certain sense, a step backward in time, as there is no CUDA, floating point textures or any feature that helps when working with non-graphic problems: you are on your own. But with the growing presence of WebGL support in modern browsers, it’s an interesting way of partially accessing the enormous computational power present in modern video cards from any JS application, without requiring the installation of a native application.

In the next posts we will explore other kinds of problem-solving where WebGL can provide a great performance boost.

This method is unarguably very simple and makes accessing strings very time efficient, as getting the string pointer only requires a multiplication and an addition. But the disadvantages are substantial:

the ids for each string are not positioned near the strings, making synchronization errors more likely;

the same storage space is used for all strings, leading to massive memory waste in most cases.

The rest of this post will be dedicated to explore an alternative method that solves these problems at the expense of requiring an initialization step and making string accesses a bit slower.

can be evaluated first as defining an enumeration value, and later as defining a string constant. This allows us to avoid the first problem, making the synchronization between the string ids an their associated string constants much easier.

To solve the second problem we can concatenate all the strings and make the accesses via an offset table. We can fill the sizes of each element using the sizeof operator over the string constants, but computing the offsets will require a runtime initialization step to do the necessary additions.

Testing

We cannot run a multiple file C program online (AFAIK), but we can simulate the inclusions manually and run it at Codepad, where it gives this output:

This is the first message.
This is the second message.

Edit 2/14 21:30

Reading this comment by Arseny Kapoulkine I realized that my previous solution is, in fact, badly overengineered and wasteful. 😀 In fact, for most purposes, we can avoid using enums at all and just use this well known solution (though it’s not really a string table…):

Source file – str_table.c

#define STR_TABLE_C
#include "str_table.h"

But if we want to be able to iterate through the string table, Arseny’s functional solution is a good solution (though it’s hard to follow for people like me that doesn’t know functional programming very well :-D).