# Bird’s DSMC0.FOR in Java

Posted on June 29th, 2013
Previous Article :: Next Article

## Introduction

G.A. Bird’s Molecular Gas Dynamics and the Direct Simulation of Gas Flows has become the de facto bible of particle-based neutral gas modeling. Prof. Bird maintains an active website dedicated to DSMC at www.gab.com.au. On this website you will find his latest solvers and source codes. You will also find the full archive of the legacy demonstration codes that are described in the book. These legacy codes can be obtained from http://www.gab.com.au/legacy.html.

## DSMC0.FOR

Among the legacy codes is DSMC0.FOR, which is the code listed in the book itself. Notwithstanding some of the other issues Prof. Bird discusses in www.gab.com.au/issues.html and DSMC07notes.pdf, my main difficulty with this program is that it is written in a pretty archaic Fortran. Those of you who had the “pleasure” of working with Fortran 77 surely remember the limits on variable names and the overall line length. Due to these, we get codes such as this SELECT subroutine:

SUBROUTINE SELECT
*--selects a potential collision pair and calculates the product of the
*--collision cross-section and relative speed
*
PARAMETER (MNM=1000,MNC=50,MNSC=400,MNSP=5,MNSG=1)
*
COMMON /MOLS  / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM)
COMMON /CELLS / CC(MNC),CG(3,MNC),IC(2,MNC,MNSG),ISC(MNSC),
&                CCG(2,MNC,MNSG,MNSG),ISCG(2,MNSC,MNSG),IG(2,MNSG)
COMMON /GAS   / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP)
COMMON /CONST / PI,SPI,BOLTZ
COMMON /ELAST / VRC(3),VRR,VR,L,M,LS,MS,CVR,MM,NN,N
*
K=INT(RF(0)*(IC(2,N,NN)-0.0001))+IC(1,N,NN)+1
L=IR(K)
*--the first molecule L has been chosen at random from group NN in cell
100   MSC=IPL(L)
IF ((NN.EQ.MM.AND.ISCG(2,MSC,MM).EQ.1).OR.
&    (NN.NE.MM.AND.ISCG(2,MSC,MM).EQ.0)) THEN
*--if MSC has no type MM molecule find the nearest sub-cell with one
NST=1
NSG=1
150     INC=NSG*NST
NSG=-NSG
NST=NST+1
MSC=MSC+INC
IF (MSC.LT.1.OR.MSC.GT.MNSC) GO TO 150
IF (ISC(MSC).NE.N.OR.ISCG(2,MSC,MM).LT.1) GO TO 150
END IF
*--the second molecule M is now chosen at random from the group MM
*--molecules that are in the sub-cell MSC
K=INT(RF(0)*(ISCG(2,MSC,MM)-0.0001))+ISCG(1,MSC,MM)+1
M=IR(K)
IF (L.EQ.M) GO TO 100
*--choose a new second molecule if the first is again chosen
DO 200 K=1,3
VRC(K)=PV(K,L)-PV(K,M)
200   CONTINUE
*--VRC(1 to 3) are the components of the relative velocity
VRR=VRC(1)**2+VRC(2)**2+VRC(3)**2
VR=SQRT(VRR)
*--VR is the relative speed
LS=IPS(L)
MS=IPS(M)
CVR=VR*SPM(1,LS,MS)*((2.*BOLTZ*SPM(2,LS,MS)/(SPM(5,LS,MS)*VRR))
&    **(SPM(3,LS,MS)-0.5))/SPM(6,LS,MS)
*--the collision cross-section is based on eqn (4.63)
RETURN
END

I wanted to make sure I have a good understanding of Bird’s DSMC algorith. However, as you can see above, it takes a bit of time and serious note taking to decipher the cryptic variable names. In addition, one peculiarity of Fortran is automatic variable declaration. Not only are the variables declared automatically, the variable name plays an important role as well. Fortran variables starting with I or N or automatically type cast as integers so we get some inherent data conversion that may not be apparent at first.

## Study of the Original Fortran Code

So I decided to convert the code to a more modern programming language – Java. However, the first thing needed was to get the original code running. I used the Fortran 77 compiler available under Cygwin (you may need to download the gfortran package):

$f77 DSMC0.for -o DSMC0.exe$ ./DSMC0.exe

The code starts by asking for a numeric input to indicate new or restart run. After you press “1” for new simulation, you should see the following output:

INPUT 0,1 FOR CONTINUING,NEW CALCULATION:-
1000 MOLECULES
DSMC0:- Move     1    1   of     4   40            0. Collisions
DSMC0:- Move     2    1   of     4   40          146. Collisions
DSMC0:- Move     3    1   of     4   40          477. Collisions
DSMC0:- Move     4    1   of     4   40          819. Collisions
DSMC0:- Move     1    2   of     4   40         1148. Collisions
DSMC0:- Move     2    2   of     4   40         1495. Collisions
DSMC0:- Move     3    2   of     4   40         1820. Collisions
DSMC0:- Move     4    2   of     4   40         2156. Collisions
DSMC0:- Move     1    3   of     4   40         2491. Collisions
DSMC0:- Move     2    3   of     4   40         2844. Collisions
DSMC0:- Move     3    3   of     4   40         3188. Collisions
DSMC0:- Move     4    3   of     4   40         3548. Collisions
DSMC0:- Move     1    4   of     4   40         3871. Collisions
...

The code will crunch for few seconds. Bird uses an interesting loop control. Instead of having just one counter, there are actually three:

1. Inner Counter is set to 4, and controls the number of time steps between diagnostic sampling of per-cell gas parameters such as mean density, velocity, momentum, etc..
2. Middle Counter is set to 40, and controls the number of sampling calls before output of simulation results to the output file DSMC0.OUT
3. Outer Counter is set to 500, and controls the number of file outputs before the simulation ends.

Hence, by default, the code will perform 4*40*500 time steps with 1000 particles, for a total of 80,000,000 molecular moves.After the simulation ends, the output file DSMC0.OUT will contain on top

FROM ZERO TIME TO TIME 2.0008986
COLLISIONS:-
17632800.    6939509.    3109089.    2685787.     716351.
6939509.    2692564.    1225935.    1052972.     278493.
3109089.    1225935.     503056.     432324.     124803.
2685787.    1052972.     432324.     360758.     107785.
716351.     278493.     124803.     107785.      26846.
TOTAL NUMBER OF SAMPLES  20000
1000  MOLECULES
8.E+7  TOTAL MOLECULAR MOVES
60804196  SELECTIONS  27281060  COLLISION EVENTS, RATIO   0.4486707
MEAN COLLISION SEPARATION  0.00097568764
...
AVERAGE TEMPERATURE  304.1602
...
RATIO OF COLLISION NUMBER TO THEORETICAL VALUE

0.999211  1.000274  1.001023  1.000735  1.000186
1.000274  0.997429  1.002613  0.999148  1.000843
1.001023  1.002613  0.987719  1.000369  0.998635
1.000735  0.999148  1.000369  0.989916  1.003636
1.000186  1.000843  0.998635  1.003636  0.937604

The DSMC0 program collides particles of 5 different species with individual partial pressures of 0.6, 0.2, 0.1, 0.08, and 0.02, respectively. The table above shows the number of collisions between particles of different species. Notice the table is symmetric – the number of collisions between species 3 and species 2 is the same as the number of collisions between species 3 and species 2. The “selections” counter indicates the number of collision pairs that were selected. Of these, the “collision events” number actually collided. The ratio should be about 0.5 – this is the logic behind the NTC scheme. Bird also uses sub-cells to select particles located closest to each other. The mean molecular separation is given by the last line.

This file then continues with listing of averaged velocities and temperatures in each cell. I am not listing this here, however, the average temperature is listed. Finally, Bird uses an analytical expression to come up with the theoretical collision rate for each species. As we can see, the ratios are close to 1.

## DSMC0.FOR Details

So far so good. Before I could start with the conversion, I wanted to make sure I have at least the minimal understanding of the code organization. The code consists of the following functions:

• INIT0: calls DATA, computes inter-species collision properties such as cross-section parameters. The values are computed by averaging parameters for species i and j. Also sets up the simulation mesh. Finally it loads particles at the prescribed temperature and density.
• SAMPIO: initializes data sampling – basically sets various arrays to zero
• MOVE0: moves particles on 1-D mesh and performs reflection at the boundaries.
• INDEXM: sorts particles by molecular groups, then cells, and finally sub-cells
• COLLM: the heart of the NTC algorithm. Computes the number of pairs to sample and then calls SELECT for each pair to select two particles. Random number is then used to determine if collision occurred and if yes, calls ELASTIC to perform the collision.
• SELECT: picks two particles of possibly different species (but belonging to the same collision group) and also computes their relative velocity. The particles are picked within subcell if suffucient number of particles, otherwise the algorithm fans out in alternating directions looking for candidates.
• ELASTIC: sets the post-collision velocities by using either the VHS or VSS algorithm. The algorithm is switched based on a user-specified parameter for each species.
• SAMPLE0: samples energy and momentum in each cell
• OUT0: uses sampled data to output progressive simulation results. This function writes the DSMC0.OUT file.
• RVELC: Samples two velocity components for a gas at a prescribed temperature. Note, although this function generates both U and V components, Bird uses just one of them, and instead calls this function three times to set all three components of velocity.
• GAM: Evaluates the Gamma function. Bird uses it to compute the cross-section per equation 4.63 in his book.
• RF: a pseudo-random number generator. I converted this code to Java to allow for direct comparison of results. However, functionally this is equivalent to your favorite RNG.
• DATA: sets user inputs such as species properties and number of time steps.

The main loop basically consists of the following

DO 200 JJJ=1,NSP
DO 150 III=1,NIS
TIME=TIME+DTM
*
WRITE (*,99001) III,JJJ,NIS,NSP,NCOL
99001     FORMAT (' DSMC0:- Move ',2I5,'   of ',2I5,F14.0,' Collisions')
*
CALL MOVE0
CALL INDEXM
CALL COLLM
*
150     CONTINUE
*
CALL SAMPLE0
*
200   CONTINUE

## First Steps: DSMC0v1.java

My first step in the conversion was to do a direct conversion of the Fortran code to Java. This involved rewriting the syntax to be compatible with Java, however, all function and variable names were left intact. The array indexing was also left in the Fortran syntax of [1:N] vs the Java/C [0:N-1]. This version also heavily relies on global variables.

For comparison, here is again the Select subroutine. Still quite ugly, but at least it runs.

public void SELECT(int N)
{
//selects a potential collision pair and calculates the product of the
//collision cross-section and relative speed
int K=(int)(RF(0)*(IC[2][N][NN]-0.0001))+IC[1][N][NN]+1;
L=IR[K];
M = L;

//the first molecule L has been chosen at random from group NN in cell
do {
MSC=IPL[L];
if ((NN==MM && ISCG[2][MSC][MM]==1) || (NN!=MM&&ISCG[2][MSC][MM]==0))
{
//if MSC has no type MM molecule find the nearest sub-cell with one
int NST=1;
int NSG=1;
do
{
int INC=NSG*NST;
NSG=-NSG;
NST=NST+1;
MSC=MSC+INC;
}	while (MSC<1 || MSC>MNSC ||
ISC[MSC] != N || ISCG[2][MSC][MM]<1);
}

//the second molecule M is now chosen at random from the group MM
//molecules that are in the sub-cell MSC
K=(int)(RF(0)*(ISCG[2][MSC][MM]-0.0001))+ISCG[1][MSC][MM]+1;
M=IR[K];
//choose a new second molecule if the first is again chosen
} while (L==M);

for (K=1;K<=3;K++)
{
VRC[K] = PV[K][L]-PV[K][M];
//VRC(1 to 3) are the components of the relative velocity
}

VRR=VRC[1]*VRC[1]+VRC[2]*VRC[2]+VRC[3]*VRC[3];
VR=Math.sqrt(VRR);//VR is the relative speed
LS=IPS[L];
MS=IPS[M];
CVR=VR*SPM[1][LS][MS]*Math.pow(2.*BOLTZ*SPM[2][LS][MS]/(SPM[5][LS][MS]*VRR),SPM[3][LS][MS]-0.5)/SPM[6][LS][MS];
//the collision cross-section is based on eqn (4.63)
}

## Fortran vs. Java discrepancies

If we were to run this and the original Fortran code for 4 time steps, we get the following:

JAVA
DSMC0:- Move 1 4   of 4 4 3871.0 Collisions
DSMC0:- Move 2 4   of 4 4 4218.0 Collisions
DSMC0:- Move 3 4   of 4 4 4557.0 Collisions
DSMC0:- Move 4 4   of 4 4 4910.0 Collisions
5260.0 Collisions

FORTRAN
DSMC0:- Move     1    4   of     4    4         3871. Collisions
DSMC0:- Move     2    4   of     4    4         4218. Collisions
DSMC0:- Move     3    4   of     4    4         4557. Collisions
DSMC0:- Move     4    4   of     4    4         4910. Collisions
5260.0 Collisions

Looks pretty good so far. However, the collision counts start to diverge as the number of iterations is increased. Specifically, I found that the results are identical for the first 124 time steps. The Fortran code then however samples one more particle than Java. Since the additional sample is associated with an additional call to the random number generator to check for collision, the Java and Fortran pseudo random sequences go out of sync.

So where is this extra particle coming from? My best guess – since I could not find any other culprit, is that it’s due to differences in floating point mathematics between Java and Fortran. Bird uses the variable CCG[2][][][] to store fractional number of selection pairs, i.e. from the Java version:

double ASEL=0.5*IC[2][N][NN]*AVN*FNUM*CCG[1][N][NN][MM]*DTM/CC[N]+CCG[2][N][NN][MM]
int NSEL=(int)ASEL;
CCG[2][N][NN][MM]=ASEL-NSEL;

The values are initially set to random numbers and are continuously updated as particle selection pairs are computed. The Table below compares initial (time step 1) values
for selected cells:

ASEL_JAVA       CCG_JAVA        ASEL_FORTRAN    CCG_FORTRAN
29.78241059	0.919028079	29.7824078	0.919028044
29.4727082	0.609325688	29.4727058	0.609325647
29.65552519	0.792142683	29.6555233	0.792142689
29.78654117	0.923158662	29.7865391	0.923158646
29.45737197	0.593989464	29.4573708	0.593989432

Although the Java and Fortran values are similar, they are not identical. Then, this is what happens at the time step where Fortran creates the additional particle:

ASEL_JAVA	NSEL_J  CCG_JAVA	ASEL_FORTRAN    NSEL_F  CCG_FOR
143.2518608	143	0.463151892	143.252274	143	0.463546753
134.2325445	134	0.877979811	134.232513	134	0.877944946
137.8954049	137	0.052230695	137.893661	137	0.050521851
137.5810259	137	0.754412055	137.583359	137	0.756713867
146.9983001	146	0.499474699	147.001602	147	0.502731323

The fractional component is approximately 0.5 in both Java and Fortran, however, in Fortran it happens to be just slightly larger than this value. This small difference is just large enough to bring the real-valued ASEL to the next integer. But given that DSMC is a statistical method, small issues such as these should not matter in the long run. Comparing the Java and Fortran versions after 80,000 time steps, we get the following

--JAVA--
AVERAGE TEMPERATURE 304.065

COLLISIONS:-
17641762	 6941561	 3105652	 2683553	  717500
6941561	 2692212	 1224665	 1054237	  277452
3105652	 1224665	  505002	  432828	  125664
2683553	 1054237	  432828	  358914	  107301
717500	  277452	  125664	  107301	   27150

RATIO OF COLLISION NUMBER TO THEORETICAL VALUE

1.0002	  1.0011	  1.0004	  1.0004	  1.0023
1.0011	  0.9978	  1.0021	  1.0009	  0.9976
1.0004	  1.0021	  0.9921	  1.0021	  1.0061
1.0004	  1.0009	  1.0021	  0.9854	  0.9997
1.0023	  0.9976	  1.0061	  0.9997	  0.9487

--FORTRAN--
AVERAGE TEMPERATURE   304.158112

COLLISIONS:-
17633570.    6942742.    3109295.    2690434.     716643.
6942742.    2685244.    1220862.    1053511.     279321.
3109295.    1220862.     505602.     430802.     124565.
2690434.    1053511.     430802.     360802.     107069.
716643.     279321.     124565.     107069.      27108.

RATIO OF COLLISION NUMBER TO THEORETICAL VALUE

0.999257  1.000742  1.001091  1.002469  1.000595
1.000742  0.994719  0.998466  0.999661  1.003820
1.001091  0.998466  0.992720  0.996849  0.996733
1.002469  0.999661  0.996849  0.990038  0.996971
1.000595  1.003820  0.996733  0.996971  0.946756

The two sets are quite close, although it seems that my Java version slightly overestimates the theoretical values while Bird’s Fortran version slightly underestimates them.

## Final Java Version: DSMC0v3.java

My next step in the conversion was to make the code more Java like – the indexing was changed to [0:N-1] and the data arrays were renamed to more descriptive terms. However, this still resulted in code such as

CellGeom[2][M] = CellWidth;

Hence, the final step was to replace the data arrays with classes. Here again is the Select subroutine:

public void SELECT(int N)
{
//selects a potential collision pair and calculates the product of the
//collision cross-section and relative speed
int K=(int)(RF(0)*(cell_info[N][NN].count-0.0001))+cell_info[N][NN].start;
L=CrossRef[K];
M = L;

//the first molecule L has been chosen at random from group NN in cell
do {
MSC=SubCellNumber[L];
if ((NN==MM && spg_cell_info[MSC][MM].count==1) ||
(NN!=MM && spg_cell_info[MSC][MM].count==0))
{
//if MSC has no type MM molecule find the nearest sub-cell with one
int NST=1;	/*these were 1, keep 1 or 0?*/
int NSG=1;
do
{
int INC=NSG*NST;
NSG=-NSG;
NST++;
MSC=MSC+INC;
}	while (MSC<0 || MSC>=NumSubCells ||
Sub2Cell[MSC] != N || spg_cell_info[MSC][MM].count<1);
}

//the second molecule M is now chosen at random from the group MM
//molecules that are in the sub-cell MSC
K=(int)(RF(0)*(spg_cell_info[MSC][MM].count-0.0001))+spg_cell_info[MSC][MM].start;
M=CrossRef[K];
//choose a new second molecule if the first is again chosen
} while (L==M);

for (int i=0;i<3;i++)
{
VRC[i] = part[L].v[i]-part[M].v[i];
//VRC(1 to 3) are the components of the relative velocity
}

VRR=VRC[0]*VRC[0]+VRC[1]*VRC[1]+VRC[2]*VRC[2];
VR=Math.sqrt(VRR);//VR is the relative speed
LS=SpeciesCodeNumber[L];
MS=SpeciesCodeNumber[M];
CVR=VR*species_interaction[LS][MS].sigma*
Math.pow(2.*BOLTZ*species_interaction[LS][MS].ref_temp/(species_interaction[LS][MS].reduced_mass*VRR),
species_interaction[LS][MS].visc_temp_index-0.5)/species_interaction[LS][MS].gamma;
//the collision cross-section is based on eqn (4.63)
}