How does this code implement the gauss-siedel method?

Question

I have been trying to port and understand the code used in the paper Real-Time Fluid Dynamics for Games. One function in the paper, the method to diffuse one of the vector fields is said to use the Gauss-Siedel relaxation to approximate the solutions.

Here is the function:

void diffuse ( int N, int b, float * x, float * x0, float diff, float dt ) 
{ 
    int i, j, k; 
    float a=dt*diff*N*N; 
 
    for ( k=0 ; k<20 ; k++ ) { 
        for ( i=1 ; i<=N ; i++ ) { 
            for ( j=1 ; j<=N ; j++ ) { 
                x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+ 
  x[IX(i,j-1)]+x[IX(i,j+1)]))/(1+4*a); 
            } 
        } 
        set_bnd ( N, b, x ); 
    } 
}

I understand what the set_bnd(N, b, x) call is doing (handling interactions with the edges of the simulation), and I know that the outermost for-loop is doing the iteration 20 times, but how does this relate to the formula for the gauss-siedel method? Also, what is the term a doing here, and why is the sum of the surrounding cells in the matrix multiplied by $\frac{a}{1+4a}$?

Answer 1

This is exactly diffusion process there's a change of matter between every two adjacent cells, thus each element tends towards the average of it's neighbors in the given rate $a$, which is proportional to the timestep, the diffusion rate, and inversely proportional to the cell area.

$ x_{n+1} (i,j) = a * \sum_{neighbors} ({x_n(neighbor) - x_n(i,j)})$

Summing it up gives

$ x_{n+1} (i,j) = a * (x_n(i+1, j) + x_n(i-1, j) + x_n(i, j+1) + x_n(i, j-1) -4x_n(i,j)) $

Substituting this exactly gives an unstable solution (what was called "diffuse bad" in the article), since this method is explicit. In order to make the method implicit, we will replace $4x_n(i,j)$ on the RHS by $4x_{n+1}(i,j)$, move the $4x_{n+1}(i,j)$ to the LHS obtaining $(1+4a)x_{n+1}(i,j)$, and dividing by $1+4a$ to get the formula you wrote. This implicitization procedure is common to improve stability of iterative methods.

The notation "Gauss-Seidel" here is correct, but misleading in my opinion, since the Gauss-Seidel algorithm works in an entirely different "framework" (when you're solving a system of linear equations $Ax=b$), though they are equivalent (in GS you decompose $A$ into $U+D+L$ and then replace $Dx_{n}$ by $Dx_{n+1}$, which exactly what we did). Please think of it purely as a method of making the iterative algorithm more implicit to improve stability.