If any man will draw up his case, and put his name at the foot of the first page, I will give him an immediate reply. Where he compels me to turn over the sheet, he must wait my leisure.
John Montagu, 4th Earl of Sandwich

Nick Trefethen has suggested that numerical programming can often best be accomplished through Ten Digit Algorithms, for which ten digits of accuracy should be achieved within one page of code to run in under five seconds. I strongly believe that such code can aid clarity, communication and experimentation. It is also fun to write. While software often seems to be designed with the philosophy that it’s perfected when there’s nothing left to add, TDAs are perfected when there’s nothing left to take away.

I have written a lattice Boltzmann code for fluid flow, whose implementation fits on a single page and runs in five seconds, producing useful results. To my knowledge, this is the most concise implementation of LBM, and could serve as an introduction to the method.

I chose to write it in Matlab, since its support for array operations and visualisation made the code clearer and concise.

2-D Model of Porous Medium Flows

% 2D Lattice Boltzmann (BGK) model of a fluid.
%  c4  c3   c2  D2Q9 model. At each timestep, particle densities propagate
%    \  |  /    outwards in the directions indicated in the figure. An
%  c5 -c9 - c1  equivalent 'equilibrium' density is found, and the densities
%    /  |  \    relax towards that state, in a proportion governed by omega.
%  c6  c7   c8      Iain Haslam, March 2006.
omega=1.0; density=1.0; t1=4/9; t2=1/9; t3=1/36; c_squ=1/3; nx=31; ny=31;
F=repmat(density/9,[nx ny 9]); FEQ=F; msize=nx*ny; CI=[0:msize:msize*7];
%a=repmat(-15:15,[31,1]);BOUND=(a.^2+a'.^2)<16;BOUND(1:nx,[1 ny])=1;
%BOUND=zeros(nx,ny);BOUND(1:nx,1)=1;%open channel
BOUND=rand(nx,ny)>0.7; %extremely porous random domain
ON=find(BOUND); %matrix offset of each Occupied Node
TO_REFLECT=[ON+CI(1) ON+CI(2) ON+CI(3) ON+CI(4) ...
            ON+CI(5) ON+CI(6) ON+CI(7) ON+CI(8)];
REFLECTED= [ON+CI(5) ON+CI(6) ON+CI(7) ON+CI(8) ...
            ON+CI(1) ON+CI(2) ON+CI(3) ON+CI(4)];
avu=1; prevavu=1; ts=0; deltaU=1e-7; numactivenodes=sum(sum(1-BOUND));
while (ts<4000 & 1e-10<abs((prevavu-avu)/avu)) | ts<100
    % Propagate
    F(:,:,4)=F([2:nx 1],[ny 1:ny-1],4);F(:,:,3)=F(:,[ny 1:ny-1],3);
    F(:,:,2)=F([nx 1:nx-1],[ny 1:ny-1],2);F(:,:,5)=F([2:nx 1],:,5);
    F(:,:,1)=F([nx 1:nx-1],:,1);F(:,:,6)=F([2:nx 1],[2:ny 1],6);
    F(:,:,7)=F(:,[2:ny 1],7); F(:,:,8)=F([nx 1:nx-1],[2:ny 1],8);
    BOUNCEDBACK=F(TO_REFLECT); %Densities bouncing back at next timestep
    DENSITY=sum(F,3);
    UX=(sum(F(:,:,[1 2 8]),3)-sum(F(:,:,[4 5 6]),3))./DENSITY;
    UY=(sum(F(:,:,[2 3 4]),3)-sum(F(:,:,[6 7 8]),3))./DENSITY;
    UX(1,1:ny)=UX(1,1:ny)+deltaU; %Increase inlet pressure
    UX(ON)=0; UY(ON)=0; DENSITY(ON)=0;
    U_SQU=UX.^2+UY.^2; U_C2=UX+UY; U_C4=-UX+UY; U_C6=-U_C2; U_C8=-U_C4;
    % Calculate equilibrium distribution: stationary
    FEQ(:,:,9)=t1*DENSITY.*(1-U_SQU/(2*c_squ));
    % nearest-neighbours
    FEQ(:,:,1)=t2*DENSITY.*(1+UX/c_squ+0.5*(UX/c_squ).^2-U_SQU/(2*c_squ));
    FEQ(:,:,3)=t2*DENSITY.*(1+UY/c_squ+0.5*(UY/c_squ).^2-U_SQU/(2*c_squ));
    FEQ(:,:,5)=t2*DENSITY.*(1-UX/c_squ+0.5*(UX/c_squ).^2-U_SQU/(2*c_squ));
    FEQ(:,:,7)=t2*DENSITY.*(1-UY/c_squ+0.5*(UY/c_squ).^2-U_SQU/(2*c_squ));
    % next-nearest neighbours
    FEQ(:,:,2)=t3*DENSITY.*(1+U_C2/c_squ+0.5*(U_C2/c_squ).^2-U_SQU/(2*c_squ));
    FEQ(:,:,4)=t3*DENSITY.*(1+U_C4/c_squ+0.5*(U_C4/c_squ).^2-U_SQU/(2*c_squ));
    FEQ(:,:,6)=t3*DENSITY.*(1+U_C6/c_squ+0.5*(U_C6/c_squ).^2-U_SQU/(2*c_squ));
    FEQ(:,:,8)=t3*DENSITY.*(1+U_C8/c_squ+0.5*(U_C8/c_squ).^2-U_SQU/(2*c_squ));
    F=omega*FEQ+(1-omega)*F;
    F(REFLECTED)=BOUNCEDBACK;
    prevavu=avu;avu=sum(sum(UX))/numactivenodes; ts=ts+1;
end
figure;colormap(gray(2));image(2-BOUND');hold on;
quiver(2:nx,1:ny,UX(2:nx,:)',UY(2:nx,:)');
title(['Flow field after ',num2str(ts),'\deltat']);xlabel('x');ylabel('y');
Flow through random grid of cells
Flow through (not quite) random grid of cells - click for details Typical Flow Field Calculated from the 2-D LBM Code.

3-D Model of Porous Medium Flows

An analagous LBM in 3-D uses the same governing equations, but requires 19 velocities. This code won’t quite fit on a single page, but I think the similarities with the 2-D code, and deliberate repetition during calculation of the equilibrium particle densities redeem the second page. Repetition is clearer and more efficient than creating a single array of weightings and concatenating the velocity arrays in order to calculate FEQ in a single line.

19 velocity vectors pointing to adjacent cubes One stationary vector, six nearest neighbours, and twelve next-nearest in 3-D.

Download the 3-D code here.

Validation

The flow profile in a fully saturated channel is calculated using the LBM and compared with the analytical solution, to demonstrate that the method works as advertised.

Note that providing ten digits of accuracy when solving engineering problems, rather than the mathematics for which TDAs were proposed, is generally neither achievable nor required.

Download the validation code here.

Comments

This page has been around since late 2005, originally at my University of Durham page, and is relatively popular. I am grateful to receive messages with corrections, simplifications or other suggestions, directed to iain@exolete.com.

Requests for answers to homework questions are ignored. I recommend that you read and play with the code, perhaps altering the parameters until you understand what it is doing. When in doubt, remember that now is a good time to put your work on a firm theoretical foundation. The fifth chapter of Sauro Succi’s “The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond” covers the Lattice BGK model used here.

And Finally…

It is almost a shame not to have been able to use the venerable Pascal for this, since Blaise emphasised the value of brevity when he noted in a letter, I would have written a shorter one, but I did not have the time.