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 writing TDAs can be aids to clarity, communication and experimentation. They’re 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.
Pascal would also be an appropriate choice, since Blaise Pascal famously highlighted the importance of brevity when apologising for the length of a letter – I would have written a shorter one, but I did not have the time .
% 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');
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.
One stationary vector, six nearest neighbours, and twelve next-nearest in 3-D.
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.
1 From the penultimate paragraph of his 16th Provincial letter. Thanks to Dr. Craig Meskell for suggesting this quote.
2 Repetition here 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.
This page has been around since 2005 (originally at my University of Durham page), and has had thousands of visitors. I receive plenty of email about it, and I'm happy to answer questions, but I won't do your homework for you. I strongly recommend Sauro Succi's The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond - Chapter 5 covers the Lattice BGK model used here.