Too Much Idle | OpenGL FDTD Tutorial: Episode 2 – Simulation Cycle and Shader Programs

Here we go, the second episode of my tutorial on how to implement real-time physical modelling synthesis on the GPU, using C++ and OpenGL shaders. The first episode can be found here.

We’ll finish with the main file, working out the simulation cycle. Then, we’ll work on the shaders finally! This means having a close look to the solver and its visual representation.
The list of covered files:

main.cpp [second–and last–part]
fbo_vs.glsl and fbo_fs.glsl [FBO shader program]
render_vs.glsl and render_fs.glsl [render shader program]

Main is still in C++, while the shader files are written in OpenGL Shading Language [GLSL to friends].

Straight to Business!

//—————————————————————————————————————————————————————————————————————
We have to start by explaining the STRUCTURE of the SIMULATION CYCLE
//—————————————————————————————————————————————————————————————————————

We left off from a perfectly initialized code, ready to flip audio samples in a continuous simulation cycle!

In the ideal situation of running our drum head in real-time [and pushing the computed samples directly to the audio output], we would syncronize the advancement of the simulation with the audio rate set by the audio interface. This could be easily done by implementing the simulation cycle inside the audio callback function likely present in the real-time audio programming environment we are using. In particular, per each audio call we would:

advance the simulation [FBO shader program calls] for as many steps as the size of the audio buffer we have to fill
retrieve the block of computed samples from the audio texture portion and copy them into the audio buffer
render the latest status of the simulation on screen [a single render shader program call]

Then return and chill out until the next callback happens.

Well, this can be theoretically done on any OS, by using a cross-platform real-time audio programming environment like PortAudio, or native solutions like Core Audio, Alsa and the Windows equivalent whose name I don’t know. However, all these solutions will have slightly different details [e.g., multi-period buffers, extra callback parameters]. To avoid confusion, for this tutorial I remained as generic as possible and instead of choosing a specific platform I simply made the code run offline, but maintaining a buffer-based structure. In other words, a number of total samples is defined a priori together with the size of a dummy audio buffer; then the steps needed for the simulation cycle to compute a buffer [the 3 bullet points] are looped, until the number of synthesized samples satisfies the chosen limit; per each iteration, instead of pushing the computed buffer to the audio card, the samples are appended to a global array where, in the end, the whole simulated audio is at our disposal.
Too many words, here is some pseudo code:

_Set number of total samples
_Set audio buffer size
_Compute number of needed buffers
_While buffer counter is less than needed buffers
___While sample counter is less than buffer size
_____Advance simulation
___Retrieve audio texture samples
___Append audio texture samples to global array
___Render to screen
_End of simulation

Not that complicated at this level of detail!
Furthermore, moving back from this approach to a real-time architecture should be a piece of cake.

Before analyzing the code, it’s good to write down some more words about the “Advance simulation” line of the pseudo code, since it’s the core of the algorithm.
As explained in this section of Episode 1, we deal with 2 copies of the domain, a snapshot of the model at time step N and a snapshot at time step N-1, at start respectively stored in quad0 and quad1. When we want to advance the simulation of one step, we use both the snapshots to calculate the new pressure values and store them in the N-1 domain; this is done in the shader, as we’ll soon see. Then we swap the indices, so that the N-1 domain becomes N and, vice versa, the N domain becomes N-1. This is done on the CPU side. According to this scheme, quad0 and quad1 alternatively contain the most recent pressure values.
This swapping also affects the way we fill the audio texture portion in the shader. When copying the audio from the listener’s position into the buffer [quad2], we have to alternatively pick from quad0 and quad1, so that all the synthesized samples are correctly lined up in the buffer. The safest way to do this, is to always sample the LEAST UPDATED quad. It sounds weird, I know, but there is a reason…
When we advance the simulation we have to execute a simulation step [update of one of the main quads] followed by an audio step [to fill the audio texture]. The 2 steps are done with 2 OpenGL calls that are non blocking, this means that if in the simulation step we update quad0 and then, in the audio step, we try to immediately read quad0 listener’s pressure value, we run the risk to pick an old value, cos the shader might not have finished to update all the fragments yet. This indeterminacy resolves into noise in the output audio stream [sometimes the picked sample is correct, sometimes it’s not…]. In this example, the solution is to pick from quad1 [which is done for sure] and, more in general, to read from the least updated quad.
This introduces a delay of one sample. Can you hear the difference?

To handle this swapping scheme into the code, we define 4 states on the CPU side, that drive the operations within the shader:

state0: draw quad0 [simulation step]
state1: read audio from quad1, cos quad0 might not be ready yet [ audio step]
state2: draw quad1 [simulation step]
state3: read audio from quad0, cos quad1 might not be ready yet [ audio step]

The simulation starts from state0, then the states are incremented and looped [0 and 1, then 2 and 3, then 0 and 1, then 2 and 3…] until the application is running.
Hopefully the code will help understand this simple concept that sounds scary in words!

Note that we haven’t talked about how these states are actually implemented in the FBO shader yet! We are simply describing how we cycle through the commands in the simulation cycle.

//—————————————————————————————————————————————————————————————————————
Now it’s time for the CODE for the SIMULATION CYCLE [main.cpp]
//—————————————————————————————————————————————————————————————————————

Let’s see in detail how the central portion of the pseudo code—the lines “While sample counter is less than buffer size” “Advance simulation”—translates into code:

 
// w<em>hile sample counter is less than buffer size</em>
for(int n=0; n<buffer_size; n++) {
       
    // advance simulation:

	// pass next excitation value
	glUniform1f(excitation_loc, excitation);

	// simulation step [state0 or state2]
	state = quad*2;
	glUniform1i(state_loc, state);
	glDrawArrays(GL_TRIANGLE_STRIP, vertices[quad][0], vertices[quad][1]); // draw quad0 or quad1

	// audio step [state1 or state3]
	glUniform2fv(wrCoord_loc, 1, wrCoord); // fragment and channel
	glUniform1i(state_loc, state+1);
	glDrawArrays(GL_TRIANGLE_STRIP, vertices[quad2][0], vertices[quad2][1]); // draw quad2, appending audio from quad1 or from quad0

	// prepare next cycle
	quad = 1-quad;
	wrCoord[1] = int(wrCoord[1]+1)%4;
	if(wrCoord[1]==0)
		wrCoord[0]+=deltaCoordX;

	// next excitation
	if(++sampleCnt == excitationDt) {
		excitation = maxExcitation; // strike!
		sampleCnt = 0;
	}
	else
		excitation = 0;

	// re-sync all parallel GPU threads [this is also called implicitly when buffers are swapped]
	glMemoryBarrier(GL_TEXTURE_FETCH_BARRIER_BIT);
}

// w<em>hile sample counter is less than buffer size</em>

for(int n=0; n<buffer_size; n++) {

// advance simulation:

// pass next excitation value

glUniform1f(excitation_loc, excitation);

// simulation step [state0 or state2]

state = quad*2;

glUniform1i(state_loc, state);

glDrawArrays(GL_TRIANGLE_STRIP, vertices[quad][0], vertices[quad][1]); // draw quad0 or quad1

// audio step [state1 or state3]

glUniform2fv(wrCoord_loc, 1, wrCoord); // fragment and channel

glUniform1i(state_loc, state+1);

glDrawArrays(GL_TRIANGLE_STRIP, vertices[quad2][0], vertices[quad2][1]); // draw quad2, appending audio from quad1 or from quad0

// prepare next cycle

quad = 1-quad;

wrCoord[1] = int(wrCoord[1]+1)%4;

if(wrCoord[1]==0)

wrCoord[0]+=deltaCoordX;

// next excitation

if(++sampleCnt == excitationDt) {

excitation = maxExcitation; // strike!

sampleCnt = 0;

}

else

excitation = 0;

// re-sync all parallel GPU threads [this is also called implicitly when buffers are swapped]

glMemoryBarrier(GL_TEXTURE_FETCH_BARRIER_BIT);

}

This is the cycle that makes the simulation advance enough to fill a single audio buffer. It is a straightforward implementation of the pseudo code, we simply add some simple extra features.

First, we pass the current excitation value to the shader as a single floating point number. More info about the excitation in a few lines.

Then, it comes the simulation step. To advance the simulation we first have to tell the shader that we are either in state0 or in state2, according to which quad needs to be updated next [i.e., the N-1 portion of the domain]; once we have sent state as uniform to the shader, we can trigger the actual update of the chosen quad by passing its vertices to the function glDrawArrays(…). At this point the FBO shader program asynchronously starts to compute the new pressure values [more details in the shader program section].

Next is the audio step. The shader needs 2 uniforms for this step. The first one is the next available audio texture’s fragment, to incrementally fill the buffer; we pass wrCoord[2], which contains both the x coordinate of the fragment [y is fixed, on top of the texture] and the index of the next available RGBA channel, so that we line up 4 samples in each fragment [one in each channel]. The second uniform is again state, which before being passed is incremented to turn state0 in state1, or state2 in state3. The shader is then called by passing to glDrawArrays(…) the vertices of quad2; the call will trigger the copy of the pressure value of the listener point in the coordinates specified in wrCoord[2] [more details in the shader program section].

The next cycle is prepared by swapping the quads’ indices and properly incrementing wrCoord[2]. We also have to compute the next excitation sample. In this simple application the excitation is an impulse train, i.e., a long line of zeros that periodically turn into a spike, according to a fixed frequency. In the above code we use two user define values to achieve this signal: the frequency of the train [expressed as the number of zero-samples between spikes, excitationDt] and the amplitude of the spike [maxExcitation].

The final sync call to glMemoryBarrier(…) is fundamental. As I wrote before, whenever we trigger the shaders with glDrawArrays(…) the pressure updates happen asynchronously on the GPU, while the CPU program keeps going forward. With the memory barrier, the CPU code stops and waits for all the GPU threads to be done before resuming its execution. This allows us to be sure that the computation of the new pressure values all over the domain is done and the current listener sample has been copied in the audio texture; then we can safely start the computation of the next iteration [simulation and audio step].

The previous code block is contained in the following lines, again translated from the pseudo code [plus some extra operations]:

 
// user defined settings
int samplerate = 44100;
float maxExcitation = 10; // dimensionless
int excitationFreq = 1; // we strike every second

// set number of total samples
int totalSampleNum = samplerate*5; // let's run for ~5 seconds

// set audio buffer size
int buffer_size = 128;

// compute number of needed buffers
int bufferNum = totalSampleNum/buffer_size;

// misc init
int state = -1;
int sampleCnt = 0;

// allocate a global audio buffer big enough to contain all the samples synthesized during the simulation
totalAudioBuffer = new float[bufferNum*buffer_size];

// excitation
float excitation = maxExcitation; // we start with a strike!
int excitationDt = samplerate/excitationFreq; // in samples

float wrCoord[2] = {0, 0}; // to store audio [sampled simulation output] in next available texture fragment and channel


// while buffer counter is less than needed buffers
for(int i=0; i<bufferNum; i++) {
        
	// switch to FBO shader
	glUseProgram(shader_program_fbo); // activate FBO shader program
	glBindFramebuffer(GL_FRAMEBUFFER, fbo); // render to our framebuffer!
	glViewport(0, 0, textureWidth, textureHeight); // full viewport, to give access to all texture




    //------------------------------------------------------
    // previous code block
    //------------------------------------------------------
    // while sample counter is less than buffer size
    for(int n=0; n<buffer_size; n++) {
            // advance simulation
    }

    //------------------------------------------------------
    //------------------------------------------------------



	// reset write coordinates at the end of each buffer
	wrCoord[0] = 0;
	wrCoord[1] = 0;


	// retrieve audio texture samples, stored in quad2 ---> intrinsically uses pbo!
	glReadPixels(0, textureHeight-1, buffer_size/4, 1, GL_RGBA, GL_FLOAT, 0); // quad2 is a single row on top of the texture, with 4 samples in each pixel
	float *sampleBuffer = (float*)glMapBuffer(GL_PIXEL_PACK_BUFFER, GL_READ_ONLY); // we need to do this...
	glUnmapBuffer(GL_PIXEL_PACK_BUFFER); //...and this, to have a physical copy of the samples in CPU memory
	
        // append audio texture samples to global array
	memcpy(totalAudioBuffer+i*(buffer_size), sampleBuffer, buffer_size*sizeof(float)); // buffer by buffer


        // switch to render shader
	glUseProgram(shader_program_render);  // activate render shader program
	glBindFramebuffer(GL_FRAMEBUFFER, 0); // disable FBO to render to screen
	glViewport(0, 0, textureWidth*magnifier, textureHeight*magnifier); // magnified viewport, so that window matches rightmost quad

    // render to screen 
    glDrawArrays(GL_TRIANGLE_STRIP, vertices[quad0][0], vertices[quad0][1]); // we pass quad0 vertices indeed
	glfwSwapBuffers(window); // put the stuff we've been drawing onto the buffer that will be displayed
}

// end of simulation!

// user defined settings

int samplerate = 44100;

float maxExcitation = 10; // dimensionless

int excitationFreq = 1; // we strike every second

// set number of total samples

int totalSampleNum = samplerate*5; // let's run for ~5 seconds

// set audio buffer size

int buffer_size = 128;

// compute number of needed buffers

int bufferNum = totalSampleNum/buffer_size;

// misc init

int state = -1;

int sampleCnt = 0;

// allocate a global audio buffer big enough to contain all the samples synthesized during the simulation

totalAudioBuffer = new float[bufferNum*buffer_size];

// excitation

float excitation = maxExcitation; // we start with a strike!

int excitationDt = samplerate/excitationFreq; // in samples

float wrCoord[2] = {0, 0}; // to store audio [sampled simulation output] in next available texture fragment and channel

// while buffer counter is less than needed buffers

for(int i=0; i<bufferNum; i++) {

// switch to FBO shader

glUseProgram(shader_program_fbo); // activate FBO shader program

glBindFramebuffer(GL_FRAMEBUFFER, fbo); // render to our framebuffer!

glViewport(0, 0, textureWidth, textureHeight); // full viewport, to give access to all texture

//------------------------------------------------------

// previous code block

//------------------------------------------------------

// while sample counter is less than buffer size

for(int n=0; n<buffer_size; n++) {

// advance simulation

}

//------------------------------------------------------

// reset write coordinates at the end of each buffer

wrCoord[0] = 0;

wrCoord[1] = 0;

// retrieve audio texture samples, stored in quad2 ---> intrinsically uses pbo!

glReadPixels(0, textureHeight-1, buffer_size/4, 1, GL_RGBA, GL_FLOAT, 0); // quad2 is a single row on top of the texture, with 4 samples in each pixel

float *sampleBuffer = (float*)glMapBuffer(GL_PIXEL_PACK_BUFFER, GL_READ_ONLY); // we need to do this...

glUnmapBuffer(GL_PIXEL_PACK_BUFFER); //...and this, to have a physical copy of the samples in CPU memory

// append audio texture samples to global array

memcpy(totalAudioBuffer+i*(buffer_size), sampleBuffer, buffer_size*sizeof(float)); // buffer by buffer

// switch to render shader

glUseProgram(shader_program_render); // activate render shader program

glBindFramebuffer(GL_FRAMEBUFFER, 0); // disable FBO to render to screen

glViewport(0, 0, textureWidth*magnifier, textureHeight*magnifier); // magnified viewport, so that window matches rightmost quad

// render to screen

glDrawArrays(GL_TRIANGLE_STRIP, vertices[quad0][0], vertices[quad0][1]); // we pass quad0 vertices indeed

glfwSwapBuffers(window); // put the stuff we've been drawing onto the buffer that will be displayed

}

// end of simulation!

Some remarks.
To compute audio samples we make use of the FBO shader, while the render shader is needed to display the result on screen. This means that we have to constantly switch back and forth between the 2 shader programs and their settings [e.g., viewport size]. In this example, we decide to render only at the end of each full buffer; going faster would be useless, since the frame rate is physically bounded to the refresh rate of the monitor we use [VSync is disabled, wait for Episode 3 to see how].

The 2 shaders have slightly different viewport settings. For the FBO shader we pass to glViewport(…) the whole texture, to grant full access to it; indeed the simulation step updates one quad and reads from the other one, while the audio step has to fill the audio texture—so we need it all. When rendering to screen we don’t need the full texture, since what we need to show is the domain only [unless we are debugging this weird structure]; however, to make the window match either quad0 or quad1 [the domain], we have to go for the full viewport again, also including the magnifier factor we applied to the window when created. Since we put no offset in the viewport, the quad that will overlap with the window will be quad0…so there’s no point in drawing quad1. Since the buffer size is generally a power of 2, at the end of the inner for loop the most updated is quad1…so why we pass quad0’s vertices? There is a trick that explains this…the most attentive of you might have guessed it already, but I will disclose it at the end of this post.
Last comment, at the end of the simulation we have all of our freshly made samples stored in totalAudioBuffer, ready to be used as we please. This marks the end of our big, inelegant main function.

Here we’re using standard OpenGL functions only (:

This concludes the main file, as promised here you can find the full C++ source file:

main.cpp

It is in an a slightly rearranged form [better structured hopefully] and it also includes some EXTRA features, useful for debugging (;
More files to come!

//—————————————————————————————————————————————————————————————————————
SHADER PROGRAMS! In GLSL
//—————————————————————————————————————————————————————————————————————

Alright, time to see how things are dealt with on the GPU side.
We have 2 shader programs, each composed of a vertex and a fragment shader. As you probably noticed, we work with VERY FEW vertices and a quite complicated texture layout, so it’s easy to guess that the vertex shaders will be quite simple, while most of the magic will happen in the fragment ones.
Let’s start from the FBO shader program and in particular by it’s vertex component:

 
in vec4 pos_and_texc;
in vec4 texl_and_texu;
in vec4 texr_and_texd;

out vec2 tex_c;
out vec2 tex_l;
out vec2 tex_u;
out vec2 tex_r;
out vec2 tex_d;

void main () {	
	gl_Position = vec4 (pos_and_texc.rg, 0.0, 1.0); // set vertex position
	// unpack and pass texture coords to fragment shader
	tex_c = pos_and_texc.ba;
	tex_l = texl_and_texu.rg;
	tex_u = texl_and_texu.ba;
	tex_r = texr_and_texd.rg;
	tex_d = texr_and_texd.ba;
}

in vec4 pos_and_texc;

in vec4 texl_and_texu;

in vec4 texr_and_texd;

out vec2 tex_c;

out vec2 tex_l;

out vec2 tex_u;

out vec2 tex_r;

out vec2 tex_d;

void main () {

gl_Position = vec4 (pos_and_texc.rg, 0.0, 1.0); // set vertex position

// unpack and pass texture coords to fragment shader

tex_c = pos_and_texc.ba;

tex_l = texl_and_texu.rg;

tex_u = texl_and_texu.ba;

tex_r = texr_and_texd.rg;

tex_d = texr_and_texd.ba;

}

This is quite self-explanatory, I just want to remark how we place all the vertices on the same plane and how the in and out variables are conceived to conveniently packed the data [vertex attributes] coming from the CPU and pass it along to the fragment shader.
Yes, we DON’T DO ANYTHING here!

Moving on to the fragment shader of the FBO program, here is the list of variables [in, out, uniforms, const and parameters] that we use:

 
in vec2 tex_c;
in vec2 tex_l;
in vec2 tex_u;
in vec2 tex_r;
in vec2 tex_d;

out vec4 frag_color;

//----------------------------------------------------------------------------------------------------------
// quads to draw
const int state0 = 0; // draw right quad 
const int state1 = 1; // read audio from left quad [cos right might not be ready yet]
const int state2 = 2; // draw left quad 
const int state3 = 3; // read audio from right quad [cos left might not be ready yet]

//----------------------------------------------------------------------------------------------------------
uniform sampler2D inOutTexture;
uniform int state;
uniform float excitation;
uniform vec2 wrCoord;   // [x coord of audioWrite pixel, RBGA index]
uniform vec2 listenerFragCoord[4];
uniform float deltaCoordX; // width of each fragment

//----------------------------------------------------------------------------------------------------------
// can modify this simulate different materials and types of boundaries
//----------------------------------------------------------------------------------------------------------
float mu = 0.001; // damping factor, the higher the quicker the damping. generally way below 1
float rho = 0.5;  // propagation factor, that combines spatial scale and speed in the medium. must be <= 0.5
float gamma = 0;  // 0 means clamped boundary [wall], 1 means free boundary 
//----------------------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------------------

in vec2 tex_c;

in vec2 tex_l;

in vec2 tex_u;

in vec2 tex_r;

in vec2 tex_d;

out vec4 frag_color;

//----------------------------------------------------------------------------------------------------------

// quads to draw

const int state0 = 0; // draw right quad

const int state1 = 1; // read audio from left quad [cos right might not be ready yet]

const int state2 = 2; // draw left quad

const int state3 = 3; // read audio from right quad [cos left might not be ready yet]

//----------------------------------------------------------------------------------------------------------

uniform sampler2D inOutTexture;

uniform int state;

uniform float excitation;

uniform vec2 wrCoord; // [x coord of audioWrite pixel, RBGA index]

uniform vec2 listenerFragCoord[4];

uniform float deltaCoordX; // width of each fragment

//----------------------------------------------------------------------------------------------------------

// can modify this simulate different materials and types of boundaries

//----------------------------------------------------------------------------------------------------------

float mu = 0.001; // damping factor, the higher the quicker the damping. generally way below 1

float rho = 0.5; // propagation factor, that combines spatial scale and speed in the medium. must be <= 0.5

float gamma = 0; // 0 means clamped boundary [wall], 1 means free boundary

//----------------------------------------------------------------------------------------------------------

The texture is a sampler2d variable, named inOutTexture to underline that we both read its pixels and update a portion of them.
The only out variable is frag_color, cos this shader simply updates the RGBA channels [a color…] of the fragments belonging to the passed quad.
The uniform deltaCoordX contains the width of each fragment and it will come super handy to sample points starting from the coordinated of the current one.
Finally, have a close look at the 3 float simulation parameters mu, rho and gamma cos they will be part of the FDTD implementation that I will show shortly!

The main function of this guy:

 
void main() {

	if( (state == state0) || (state == state2) ) {	
		// regular FDTD calculation
		frag_color = computeFDTD();
	}
	else {
		// save audio routine
		frag_color = saveAudio();
	}
};

void main() {

if( (state == state0) || (state == state2) ) {

// regular FDTD calculation

frag_color = computeFDTD();

}

else {

// save audio routine

frag_color = saveAudio();

}

};

Nice! —> at this point you should immediately get why we have this structure, right?!

Let’s focus on computeFDTD() first, a.k.a. the simulation step, a.k.a the FDTD solver.
Here is the mathematical implementation, which shows how we actually use all the data we made available between CPU and GPU:

Equation number 2 describes the boundary conditions.
As mentioned before, both the equations contain the 3 simulation parameters that are defined as user variables in the shader!

And here is the code implementation (:

 
vec4 computeFDTD() {
	// current point needs current pressure, previous pressure, type [checked later on]
	vec4 frag_c  = texture(inOutTexture, tex_c);
	vec4 p       = frag_c.rrrr; // p_n, this is turned into a vector cos used in parallel neighbor computation
	float p_prev = frag_c.; // p_n-1
	
	// neighbours [pl_n, pr_n, pu_n, pd_n] need current pressure and type [only b...is boundary?]
	vec4 p_neigh;
	vec4 b_neigh;
	
	// left
	vec4 frag_l = texture(inOutTexture, tex_l);
	p_neigh.r   = frag_l.r;
	b_neigh.r   = frag_l.b;
	
	// up
	vec4 frag_u = texture(inOutTexture, tex_u);
	p_neigh.g   = frag_u.r;
	b_neigh.g   = frag_u.b;
	
	// right
	vec4 frag_r = texture(inOutTexture, tex_r);
	p_neigh.b   = frag_r.r;
	b_neigh.b   = frag_r.b;

	// down
	vec4 frag_d = texture(inOutTexture, tex_d);
	p_neigh.a   = frag_d.r;
	b_neigh.a   = frag_d.b;
	
	// parallel computation of pLRUD
	vec4 pLRUD = p_neigh*b_neigh + p*(1-b_neigh)*gamma; // gamma to simulate clamped or free edge
	
	// assemble equation
	float p_next = 2*p.r + (mu-1) * p_prev;
	p_next += (pLRUD.x+pLRUD.y+pLRUD.z+pLRUD.w - 4*p.r) * rho;
	p_next /= mu+1;
			
	// add excitation if this is excitation point [piece of cake]
	int is_excitation = int(frag_c.a);
	p_next += excitation * is_excitation;
		
	//          p_n+1    p_n  boundary? excitation?
	return vec4(p_next,  p.r, frag_c.b, frag_c.a); // pack and return (: [we maintain channels B and A intact]
}

vec4 computeFDTD() {

// current point needs current pressure, previous pressure, type [checked later on]

vec4 frag_c = texture(inOutTexture, tex_c);

vec4 p = frag_c.rrrr; // p_n, this is turned into a vector cos used in parallel neighbor computation

float p_prev = frag_c.; // p_n-1

// neighbours [pl_n, pr_n, pu_n, pd_n] need current pressure and type [only b...is boundary?]

vec4 p_neigh;

vec4 b_neigh;

// left

vec4 frag_l = texture(inOutTexture, tex_l);

p_neigh.r = frag_l.r;

b_neigh.r = frag_l.b;

// up

vec4 frag_u = texture(inOutTexture, tex_u);

p_neigh.g = frag_u.r;

b_neigh.g = frag_u.b;

// right

vec4 frag_r = texture(inOutTexture, tex_r);

p_neigh.b = frag_r.r;

b_neigh.b = frag_r.b;

// down

vec4 frag_d = texture(inOutTexture, tex_d);

p_neigh.a = frag_d.r;

b_neigh.a = frag_d.b;

// parallel computation of pLRUD

vec4 pLRUD = p_neigh*b_neigh + p*(1-b_neigh)*gamma; // gamma to simulate clamped or free edge

// assemble equation

float p_next = 2*p.r + (mu-1) * p_prev;

p_next += (pLRUD.x+pLRUD.y+pLRUD.z+pLRUD.w - 4*p.r) * rho;

p_next /= mu+1;

// add excitation if this is excitation point [piece of cake]

int is_excitation = int(frag_c.a);

p_next += excitation * is_excitation;

// p_n+1 p_n boundary? excitation?

return vec4(p_next, p.r, frag_c.b, frag_c.a); // pack and return (: [we maintain channels B and A intact]

}

Nothing to highlight here, except for the parallel computation of the neighbors, which makes the code slightly more complicated.
Let’s conclude with saveAudio(), which implements the audio step, i.e., lines up the audio samples into the audio texture [quad2] sampling from the listener fragment:

 
vec4 saveAudio() {
	// write location
	float writeXCoord  = wrCoord[0]; // this helps chose what pixel to write audio to
	float writeChannel = wrCoord[1]; // this determines what pixel's channel to write to
		 
	// first of all copy all the 4 values stored in previous step 
	vec4 color = texture(inOutTexture, tex_c);
	
	// then see if this fragment is the one supposed to save the new audio sample...	
	float diffx = tex_c.x-writeXCoord;
	
	// is this the next available fragment?
	if( (diffx  < deltaCoordX) && (diffx >= 0) ) {
		// sample chosen grid point
		int readState = 1-(state/2); // index of the state we are reading from
		vec2 audioCoord = listenerFragCoord[readState]; 
		vec4 audioFrag = texture(inOutTexture, audioCoord); // get the audio info from the listener
		
		// silence boundaries, using b
		float audio = audioFrag.r * audioFrag.b;
		
		// put in the correct channel, according to the audio write command sent from cpu	
		color[int(writeChannel)] = audio;
	}	 
	return color;
}

vec4 saveAudio() {

// write location

float writeXCoord = wrCoord[0]; // this helps chose what pixel to write audio to

float writeChannel = wrCoord[1]; // this determines what pixel's channel to write to

// first of all copy all the 4 values stored in previous step

vec4 color = texture(inOutTexture, tex_c);

// then see if this fragment is the one supposed to save the new audio sample...

float diffx = tex_c.x-writeXCoord;

// is this the next available fragment?

if( (diffx < deltaCoordX) && (diffx >= 0) ) {

// sample chosen grid point

int readState = 1-(state/2); // index of the state we are reading from

vec2 audioCoord = listenerFragCoord[readState];

vec4 audioFrag = texture(inOutTexture, audioCoord); // get the audio info from the listener

// silence boundaries, using b

float audio = audioFrag.r * audioFrag.b;

// put in the correct channel, according to the audio write command sent from cpu

color[int(writeChannel)] = audio;

}

return color;

}

Again, nothing to declare here, other than the fact that we have to retrieve the index of the quad we read from [readState = 0 or 1] from the state uniform [whose values here are 1 or 3].

It’s time to conclude this long post with the render shader program!
The render shader program does nothing but sticking the main texture where we do all of our calculations on the passed quad and display it on the window…sort of…since each fragment will be characterized by a color that depends on the RGBA channels of the texture, according to the type of grid point it is related to, but it is not directly copied from the texture.
Here its the super simple vertex shader:

 
in vec4 pos_and_texc;

out vec2 tex_c;


void main () {
	gl_Position = vec4(pos_and_texc.rg, 0.0, 1.0); // set position
	tex_c       = pos_and_texc.ba;                 // pass texture coord along to fragment shader
}

in vec4 pos_and_texc;

out vec2 tex_c;

void main () {

gl_Position = vec4(pos_and_texc.rg, 0.0, 1.0); // set position

tex_c = pos_and_texc.ba; // pass texture coord along to fragment shader

}

Here we use only a subset of the attributes we pass from the main file, namely the position of the vertex and the coordinates of the current fragment. The attributes are shared between the 2 shader programs [as well as between the 2 shaders in each program], so we can only choose from the list of attributes we prepared for the FBO shaders.

Finally, here we have the fragment shader.
This guy’s main tasks consist of detecting what type of grid point the current texture coordinate is associated to and of coloring it accordingly, putting particular care in a proper representation of the pressure values contained in each regular point.
First are the used variables [in and out, uniforms, const and parameters] and a function that indeed turns the a pressure into a colored field:

 
in vec2 tex_c;

out vec4 frag_color;

//----------------------------------------------------------------------------------------------------------
uniform sampler2D inputTexture;
uniform vec2 listenerFragCoord;
uniform vec2 deltaCoord;
//----------------------------------------------------------------------------------------------------------



//----------------------------------------------------------------------------------------------------------
const vec4 color_white  = vec4(1);
const vec4 color_black  = vec4(0);
const vec4 color_red    = vec4(0.9, 0, 0, 1);
const vec4 color_green  = vec4(0, 1, 0, 1);
const vec4 color_blue   = vec4(0, 0, 1, 1);
const vec4 color_yellow = vec4(1, 1, 0, 1);
const vec4 color_pink   = vec4(1, 0.0, 1, 1);
const vec4 color_purple = vec4(0.6, 0.0, 1, 1);, 1, 1);

float maxVol = 20;
const float log_min = 10.0f;
//----------------------------------------------------------------------------------------------------------


vec4 get_color_from_pressure(float p) {
	if(p > maxVol)
		return color_red;
	else if (p<-maxVol)
		return color_purple;
		
	float slider = abs(p/maxVol);
	float log_scale = (log(slider + 1e-8) + log_min) / log_min;
	float clamp_scale = clamp(log_scale, 0, 1);

	if (p > 0)
		return mix(vec4(0), color_white, clamp_scale);
	return mix(vec4(0), color_blue, clamp_scale);
}

in vec2 tex_c;

out vec4 frag_color;

//----------------------------------------------------------------------------------------------------------

uniform sampler2D inputTexture;

uniform vec2 listenerFragCoord;

uniform vec2 deltaCoord;

//----------------------------------------------------------------------------------------------------------

const vec4 color_white = vec4(1);

const vec4 color_black = vec4(0);

const vec4 color_red = vec4(0.9, 0, 0, 1);

const vec4 color_green = vec4(0, 1, 0, 1);

const vec4 color_blue = vec4(0, 0, 1, 1);

const vec4 color_yellow = vec4(1, 1, 0, 1);

const vec4 color_pink = vec4(1, 0.0, 1, 1);

const vec4 color_purple = vec4(0.6, 0.0, 1, 1);, 1, 1);

float maxVol = 20;

const float log_min = 10.0f;

//----------------------------------------------------------------------------------------------------------

vec4 get_color_from_pressure(float p) {

if(p > maxVol)

return color_red;

else if (p<-maxVol)

return color_purple;

float slider = abs(p/maxVol);

float log_scale = (log(slider + 1e-8) + log_min) / log_min;

float clamp_scale = clamp(log_scale, 0, 1);

if (p > 0)

return mix(vec4(0), color_white, clamp_scale);

return mix(vec4(0), color_blue, clamp_scale);

}

Again, the only out var is frag_color, the new color of the current fragment. We computed it from the data carried in inputTexture. We also declare a bunch of useful colors and some settings that are used in get_color_from_pressure() to turn the passed pressure value into a color field that fluctuates between dark blue [smallest pressure value] to bright white [maximum pressure value], passing through pitch black.

Awesome, now let’s see how the main uses the RGBA channels of the texture to color the fragments:

 
void main () {
	frag_color = texture(inputTexture, tex_c); // read current RGBA values
	
	// decode type
	float boundary = 1-frag_color.b;
	float excitation = frag_color.a;
	
	// color cells according to type
	
	// regular transmissive point, the most likely case
	if(boundary == 0) {
		float p = frag_color.r;
		frafrag_colorg_c = get_color_from_pressure(p);
	}
	// boundary
	else if(boundary == 1) { 
		frag_color = color_green;
	}
	
	// excitation
	if(excitation == 1) {
		float p = frag_color.r;
		vec4 color = get_color_from_pressure(p);
		
		frag_color = mix(color, color_yellow, 0.5); // a bit transparent, to see propagation
	}
	
	
	// listener
	vec2 posDiff = vec2(tex_c.x - listenerFragCoord.x, tex_c.y - listenerFragCoord.y);
	vec2 absDiff = vec2(abs(posDiff.x), abs(posDiff.y));
	if(absDiff.x<deltaCoord.x/2 && absDiff.y<deltaCoord.y/2)
		frag_color = mix(frag_color, color_pink, 0.5); // a bit transparent, to see what's underneath
}

void main () {

frag_color = texture(inputTexture, tex_c); // read current RGBA values

// decode type

float boundary = 1-frag_color.b;

float excitation = frag_color.a;

// color cells according to type

// regular transmissive point, the most likely case

if(boundary == 0) {

float p = frag_color.r;

frafrag_colorg_c = get_color_from_pressure(p);

}

// boundary

else if(boundary == 1) {

frag_color = color_green;

}

// excitation

if(excitation == 1) {

float p = frag_color.r;

vec4 color = get_color_from_pressure(p);

frag_color = mix(color, color_yellow, 0.5); // a bit transparent, to see propagation

}

// listener

vec2 posDiff = vec2(tex_c.x - listenerFragCoord.x, tex_c.y - listenerFragCoord.y);

vec2 absDiff = vec2(abs(posDiff.x), abs(posDiff.y));

if(absDiff.x<deltaCoord.x/2 && absDiff.y<deltaCoord.y/2)

frag_color = mix(frag_color, color_pink, 0.5); // a bit transparent, to see what's underneath

}

Everything’s straightforward, the type of the texture fragment is decoded from the BA channels and we assign either a fixed color [to boundaries and excitations], or a dynamic color obtained by passing to get_color_from_pressure() the current pressure saved in the R channel.

There is only one important comment left!
If you paid a lot of attention to these first 2 episodes, you’ll probably remember that tex_c does NOT exactly contain the texture coordinates of the current fragment, rather it contains the coordinates fragment AT THE PREVIOUS TIME STEP [i.e., on the opposite side of the texture, in the “dual” domain]. This is done to optimize the number of vertex attributes passed through the VBO, as explained in Episode 1. Sooo, whenever, we check the type of the local grid point [to compute the fragment’s color] we actually access the corresponding fragment on the other side of the domain. Luckily, the type do not change over time, they are always the same at every time step, so checking either side of the texture is fine.
We also access the previous listener fragment with listenerFragCoord; this is a uniform, which means it may change dynamically in a real-time application, making concrete the risk of displaying a non-updated position for a frame…
beside the fact that this would not be that big of a deal [one frame delay, come on!], the listener fragments for quad0 and quad1 are crossed [see the initialization code in Episode 1], to make sure that every sample placed in the audio texture has been correctly updated before sampling. So, in the end we’re so lucky to get to display the latest listener position without bothering at all!

Holy cow, we’re done with the second episode too!
Before jumping into the third and last episode, here are the 4 GLSL files for the 2 shader programs.
FBO shader program:

Render shader program:

In the next episode you’re gonna get the source code for the full project.

c u later_