Bird’s DSMC0.FOR in Java
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:
- 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..
- Middle Counter is set to 40, and controls the number of sampling calls before output of simulation results to the output file DSMC0.OUT
- 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.937604The 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.946756The 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) }
Download the Source Codes
Download link: DSMC0-java The archive includes the Java v1 and v3 codes as well as comparison output log files. To download the original Fortran code, please go to Prof. Bird’s website, http://www.gab.com.au/legacy.html