[たのしいスパコンシミュレーション]部屋のほこりのシミュレーション

# はじめに
部屋を掃除するのがめんどくさい人は多いと思います(私もそうです)。掃除しないでほうっておくとどうなるのでしょうか。

観察したところによれば、ほこりは部屋の隅っこに集まっているように見えます。時々部屋の隅っこに集まっているほこりをごみ箱に捨てればわざわざ掃除しなくてもそれなりに部屋はきれいになりそうな気もします。

この現象を数値計算でシミュレーションしてみましょう。

# シミュレーションのしくみ

実際の部屋にはいろいろなモノが置いてあったりしますが、計算の便宜上正方形の何もない部屋を想定します。

まず最初に、部屋にランダムにほこりの粒子をばらまきます。ある程度たくさんばらまいたほうがマスとしての挙動がわかりやすいです。

部屋の中の空気はランダムに動いているものとしますが、部屋の中心部分が一番空気の動きが早く、壁に近いほど遅くなります。

ループ1回ごとに各ほこりの速度と位置を求め、最後に部屋を升目に区切って集計をとればほこりの分布がわかります。

なお、観察によれば部屋の隅っこ近傍でほこり同士がくっつきあって「デカほこり」になっていることが多いですが、今回は「デカほこり」現象は無視します(ほこりの粒子同士の相互作用は無視する)。

# プログラム

プログラム(FORTRAN IV)は次の通りです。前回NanoPi-NEO4でスパコンを作ったときのプログラムと同じです。

各ほこり粒子の相互作用は考えないプログラムなので並列化しやすいこと、やっていることがしょーもない割に計算回数が多いことからMPI並列とします。

```hokori_mpi.f
C     HOKORI
C     SIMULATES MOVEMENT OF DUST IN ROOM
      INCLUDE 'mpif.h'
      INTEGER IERR,IRANK,ISIZE,LEN
      DOUBLE PRECISION TSTART,TEND,TD,MFLOPS
      DIMENSION ISTATUS(MPI_STATUS_SIZE)
      DIMENSION ROOM(20,20)
      DIMENSION IROOM(20,20)
      DIMENSION DUST(2,10000)
      DIMENSION JROOM(20,20)
C     SETUP FOR SIMULATION
      CALL MPI_INIT(IERR)
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR)
      IF(IRANK.NE.0) GO TO 5
      READ(5,100) NSTEP,NDUST
  100 FORMAT(I8,I5)
    5 CALL MPI_BCAST(NSTEP,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      CALL MPI_BCAST(NDUST,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
      IF(IRANK.NE.0) GO TO 90
C     STEP 1 GENERATES DUST IN ROOM
      DO 10 I=1,NDUST
      X = RAND(0)*2.0-1.0
      Y = RAND(0)*2.0-1.0
      DUST(1,I)=X
      DUST(2,I)=Y
   10 CONTINUE
C     STEP 2 GENERATES ROOM CONDX
      DO 20 I=1,20
      DO 20 J=1,20
      IF(I.LT.8.OR.I.GE.12) GO TO 21
      IF(J.LT.8.OR.J.GE.12) GO TO 21
      ROOM(I,J)=0.05
      GO TO 20
   21 IF(I.LT.5.OR.I.GE.15) GO TO 22
      IF(J.LT.5.OR.J.GE.15) GO TO 22
      ROOM(I,J)= 0.025
      GO TO 20
   22 IF(I.LT.2.OR.I.GE.18) GO TO 23
      IF(J.LT.2.OR.J.GE.18) GO TO 23
      ROOM(I,J)= 0.0125
      GO TO 20
   23 ROOM(I,J)= RAND(0)*0.01
   20 CONTINUE
   90 CALL MPI_BCAST(DUST,2*NDUST,MPI_REAL,0,MPI_COMM_WORLD,IERR)
      CALL MPI_BCAST(ROOM,20*20,MPI_REAL,0,MPI_COMM_WORLD,IERR)
C     START SIMULATION
      TSTART=MPI_WTIME()
      DO 30 I=1,NSTEP
      DO 30 N=IRANK+1,NDUST,ISIZE
      IX = IFIX(DUST(1,N)*10.0+10.0)+1
      IY = IFIX(DUST(2,N)*10.0+10.0)+1
      X = DUST(1,N)+(RAND(0)-0.5)*ROOM(IX,IY)
      Y = DUST(2,N)+(RAND(0)-0.5)*ROOM(IX,IY)
      IF(ABS(X).LE.0.95.AND.ABS(Y).LE.0.95) GO TO 35
      IF(X.LE.0.0) X= -0.99
      IF(X.GT.0.0) X= 0.99
      IF(Y.LE.0.0) Y= -0.99
      IF(Y.GT.0.0) Y = 0.99
   35 DUST(1,N)=X
      DUST(2,N)=Y
   30 CONTINUE
C     PRINT SIMULATION RESULT
      DO 40 I=1,20
      DO 40 J=1,20
      IROOM(I,J)=0
      JROOM(1,J)=0
   40 CONTINUE
      DO 50 N=IRANK+1,NDUST,ISIZE
      IX = IFIX(DUST(1,N)*10.0+10.0)+1
      IY = IFIX(DUST(2,N)*10.0+10.0)+1
      IROOM(IX,IY)=IROOM(IX,IY)+1
   50 CONTINUE
      CALL MPI_REDUCE(IROOM,JROOM,20*20,MPI_INTEGER,MPI_SUM,0,
     1MPI_COMM_WORLD,IERR)
      TEND=MPI_WTIME()
      TD=TEND-TSTART
      IF(IRANK.NE.0) GO TO 70
      MFLOPS=DBLE(FLOAT(NSTEP))*DBLE(FLOAT(NDUST))*14.0
      MFLOPS=(MFLOPS+DBLE(FLOAT(NDUST))*6.0)*1.0E-6/TD
      WRITE(6,200) NSTEP,NDUST
  200 FORMAT(1H ,8HNSTEP = ,I8,9H NDUST = ,I5)
      DO 60 I=1,20
      WRITE(6,210) (JROOM(I,J),J=1,20)
  210 FORMAT(1H ,20I3)
   60 CONTINUE
      WRITE(6,220) TD,MFLOPS
  220 FORMAT(1H ,12HTIME(SEC) : ,F8.3/1H ,12HMFLOPS    : ,F8.3)
   70 CALL MPI_FINALIZE(IERR)
      STOP
      END
```

初期化はMPIで並列化したプロセスのランク0で行います。

ソース中のStep1 (DO 10 .... 10 CONTINUE)で部屋にランダムにほこりをばらまいています。実際の部屋は3次元ですがここでは2次元の部屋を考えているので、ほこりの数だけX座標とY座標を求めています。

Step2 (DO 20 .... 20 CONTINUE)は部屋の各部分の空気の動きの重みを決めています。

以上の下準備ができたらMPI_BCAST()で各ランクに部屋の状態データを放送して共有させます。

DO 30 .... 30 CONTINUEの部分では乱数を求め、その乱数と部屋の各部分の空気の動きの重みを掛け算して次ステップ時のほこりの位置を求めます。

あらかじめ指定したステップ数だけループしたら、各ランクでほこりの位置の集計をとり、最後にMPI_REDUCE()を呼び出してランク0に集計をとって出力します。

(なお、このプログラムはFORTRAN IVで書いてありますが、一部変数名がFORTRAN IVの制限である6文字を超えています)

# 実行結果

プログラムは次のコマンドでコンパイルできます。

```
$ mpif77 -o hokori_mpi -Os hokori_mpi.f
```

実行は次のようにします。このプログラムは実行開始時にほこりの数とステップ数を読み込みますので、echoコマンドの出力を食わせています。商用のスパコンで実行するときはそのシステムのやり方に合わせてください。

なお、実行は自作スパコン「SUGIBAC S-820」(NanoPi-NEO4 x3台並列 18コア)で行いました。

```
$ echo " 1000000 1000" | mpiexec -np 18 --hostfile ./nodelist ./hokorii_mpi
```

1000個のほこりを用意して1千万ステップのシミュレーションを実行しました。実行結果はこちら。

```
 NSTEP =        1 NDUST =  1000
  21  3  1  1  2  0  1  0  0  1  0  1  0  3  1  0  0  0  2 26
   0  2  2  4  1  1  2  3  2  4  4  4  1  3  0  1  1  4  0  1
   0  3  2  4  6  4  2  3  1  1  4  1  2  3  2  2  1  1  2  3
   1  0  2  3  3  0  2  2  2  3  1  2  1  3  1  4  1  1  4  0
   1  3  2  2  6  1  3  1  2  0  2  3  1  3  3  3  5  2  1  1
   3  2  3  6  1  1  2  7  4  4  3  3  2  3  2  2  2  1  3  3
   1  3  2  5  1  2  1  2  2  2  2  3  3  2  1  1  0  0  1  1
   0  3  2  4  4  2  2  6  1  3  6  6  3  2  1  2  7  4  2  0
   4  1  3  4  2  3  3  2  2  2  2  4  4  7  3  1  5  2  2  0
   0  4  4  2  1  4  1  6  2  2  3  1  2  1  2  5  3  3  4  1
   1  6  5  2  2  3  2  3  1  7  5  2  0  2  3  1  2  0  3  0
   2  1  3  1  2  2  2  3  4  1  6  2  3  2  3  1  7  3  2  3
   2  0  3  1  0  1  2  2  3  2  0  2  4  4  3  2  2  2  2  0
   0  3  3  1  5  1  6  0  2  3  2  4  2  5  4  3  5  1  3  2
   2  2  4  2  1  2  3  2  3  1  3  4  3  3  4  1  3  3  2  0
   2  3  1  2  4  2  0  4  0  6  2  1  3  3  6  4  1  1  4  1
   2  1  2  3  4  4  2  2  2  4  3  1  2  4  3  2  3  3  5  2
   2  2  2  2  3  5  4  0  7  2  2  2  1  4  1  5  1  5  0  1
   1  2  2  3  6  2  0  1  3  1  2  2  2  3  0  2  1  1  3  0
  28  1  3  2  3  2  1  0  1  1  1  2  1  0  2  0  1  0  1 25
 TIME(SEC) :    0.003
 MFLOPS    :    5.783
 NSTEP =     8000 NDUST =  1000
 114  7  2  1  0  2  3  0  1  1  3  9  0  0  1  0  0  0  0 84
   0  0  1  0  1  1  1  0  1  4  0  1  1  1  0  2  0 10  0  0
   0  0  1  2  1  2  3  4  4  1  4  1  0  1  2  2  1  0  1  0
   1  2  3  2  3  1  3  1  3  3  0  0  5  2  2  1  3  0 13  1
   1  0  2  0  1  0  0  0  0  1  1  0  1  0  1  0  1  2  1  0
   8  1  1  2  1  1  0  2  1  0  0  1  0  0  5  3  0  2  0  0
   0  1  0  0  0  0  2  0  0  0  0  0  1  3  2  1  0 34  0  0
   1  1  0  3  1  1  2  0  0  1  1  0  0  0  0  0  0  0  0  0
   0  3  3  1  0  2  0  1  1  0  1  0  0  1  1  0  0  0  0  0
   0  0  5  4  3  0  1  2  0  0  0  2  1  0  0  0  1  0  0  0
   8  0  3  5  0  2  0  1  1  0  1  0  0  1  0  0  0  0  0  0
   1  1  2  2  0  1  1  0  1  0  0  1  2  0  1  0  0 42  0  1
   2  1  2  2  0  0  1  3  1  1  1  0  0  0  4  2  1  2  0  0
   3  1  1  2  0  0  1  0  1  0  0  2  1  0  3  3  2  9  1  0
   1  0  0  1  6  0  4  4  2  1  4  2  0  1  3  3  1  2  1  5
   0  0  1  2  1  1  0  4  2  2  2  1  1  3  3  1  1  5  4  0
   0  0  0  1  0  0  1  0  2  1  3  0  3  2  2  2  0  2  0  1
   0  4  2  1 29  0 23  4  9  2  5  2  8  1  1  1  3  4  2  1
   0  0  1  0  0  0  1  5  0  1 14  1  0  8  4  0  6  0  2  0
  89  3  0  3  0  0  0  3  0  0  0  4  0  0  6  0  0  2  0 60
 TIME(SEC) :    0.179
 MFLOPS    :  624.223
 NSTEP =    16000 NDUST =  1000
 162  7  2  1  0  0  1  0  0  0  1  7  0  1  0  0  0  0  0107
   0  0  0  0  0  0  0  0  0  0  1  0  0  1  0  0  1 15  0  0
   0  0  0  1  0  0  0  2  0  2  1  1  1  0  0  0  0  0  1  0
   0  1  0  0  3  1  2  2  3  1  1  1  2  0  0  1  0  0 16  0
   0  1  0  1  0  1  1  0  0  0  0  0  0  0  2  3  0  2  0  0
   6  0  0  2  0  1  1  0  1  1  0  0  1  1  2  2  1  0  1  0
   0  0  0  0  0  0  0  0  0  0  0  1  1  0  0  1  1 48  0  0
   0  0  0  1  1  1  0  1  0  0  0  0  0  0  0  0  0  0  0  0
   0  1  1  1  0  0  0  0  0  1  1  0  0  0  0  0  0  0  0  0
   1  1  1  1  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0
   7  1  1  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0
   0  0  2  1  0  0  0  1  1  0  1  1  0  1  0  0  0 53  0  1
   2  1  0  3  0  0  0  1  0  0  1  1  0  0  2  0  0  1  0  0
   2  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  7  2  0
   0  1  0  0  2  0  0  0  2  1  2  1  1  3  0  0  0  3  3  5
   3  0  0  2  1  0  0  0  0  4  2  2  0  1  3  1  1  2  2  3
   0  1  0  0  0  1  1  0  1  1  0  2  2  1  1  1  1  2  2  1
   0  5  3  0 32  1 23  2  2  5  4  3  4  5  2  2  1  5  1  4
   0  1  1  0  0  0  0  9  1  5 23  3  1  7  2  1  7  0  1  0
 113  3  0  2  0  0  0  4  0  1  0  3  0  0 16  0  0  5  0 72
 TIME(SEC) :    0.328
 MFLOPS    :  682.674
 NSTEP =    24000 NDUST =  1000
 189  7  1  0  0  0  1  0  0  0  0  4  0  0  0  0  0  0  0124
   0  0  0  0  0  0  1  0  0  0  0  0  0  1  0  0  0 16  1  0
   0  0  1  1  1  0  0  0  0  0  1  0  0  0  1  0  0  1  0  0
   0  0  1  0  1  2  0  1  1  0  0  0  2  0  0  1  0  1 14  0
   1  1  0  0  0  0  0  1  0  0  0  0  0  0  1  0  0  2  1  0
   0  0  0  0  0  0  0  0  1  0  0  1  0  0  0  0  0  0  1  0
   0  0  0  0  0  0  1  0  0  0  0  0  0  0  1  0  1 53  0  0
   1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   1  0  0  1  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0
   2  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   1  0  0  0  0  0  0  1  0  0  0  0  0  1  0  0  1 57  0  1
   1  0  1  0  0  0  0  0  0  0  0  0  0  0  1  1  0  0  0  0
   0  0  2  0  0  0  0  1  1  0  0  0  0  0  0  0  0  1  2  0
   0  0  0  0  0  0  1  0  1  3  0  1  0  3  0  0  0  1  3  9
   0  0  1  0  0  1  0  1  0  1  0  0  1  0  0  0  1  1  1  0
   0  0  0  0  0  0  1  1  0  0  1  0  1  1  0  1  0  2  1  0
   0  3  1  0 34  0 24  0  0  1  1  6  4  2  0  1  1  4  1  3
   0  1  1  0  0  1  0  6  0  9 34  4  1  6  1  0 11  1  2  0
 135  3  0  2  0  0  0  6  0  0  0  3  0  0 22  0  0  6  0 98
 TIME(SEC) :    0.420
 MFLOPS    :  799.129
 NSTEP =    32000 NDUST =  1000
 203  5  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0134
   0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0 17  0  0
   0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0
   0  0  0  0  1  0  1  0  0  0  0  1  0  1  0  1  0  1 11  0
   0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  2  0  1
   2  0  0  1  0  0  0  0  0  0  0  0  0  0  0  1  0  1  0  0
   0  0  0  0  0  2  0  0  0  0  0  0  0  0  0  0  0 55  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0 59  1  0
   0  1  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0
   0  0  1  0  0  0  0  0  0  0  0  1  0  0  2  2  0  2  1 11
   0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  1  0  1  2
   0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  1  0  3  0  1
   1  2  1  0 38  0 26  1  0  0  1  1  0  2  1  0  0  3  0  0
   0  1  0  0  0  0  0  4  3  3 47  0  3  1  5  1  7  3  1  0
 142  3  0  1  0  0  0  7  0  0  0  5  0  0 28  0  0  8  0108
 TIME(SEC) :    0.559
 MFLOPS    :  802.105
 NSTEP =    40000 NDUST =  1000
 206  5  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0138
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 19  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0 10  0
   0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  1  1  0
   2  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  1 56  0  0
   0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
   0  0  0  0  0  0  0  0  0  1  0  0  0  0  1  0  0 60  1  0
   1  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0
   0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0
   0  0  0  0  0  1  0  0  1  0  1  0  0  0  0  0  0  0  0 11
   0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1
   0  0  0  0  0  0  0  0  0  1  0  3  2  0  0  0  0  0  1  1
   0  3  0  0 38  0 27  0  1  0  2  0  0  0  0  0  0  2  0  0
   0  0  0  0  1  0  0  0  0  2 49  0  0  0  1  1  6  1  1  0
 148  4  0  0  0  0  0 10  0  0  0  8  0  0 29  1  0 10  0117
 TIME(SEC) :    0.686
 MFLOPS    :  816.361
```

初期状態ではランダムに分布していたほこりが見事に部屋の隅っこに集まっています。一部変なところに固まっているやつもいますが、部屋の中央付近は十分きれいになっています。

# まとめ

上記結果より、掃除をしないでほうっておくだけで部屋のほこりは隅に固まることがわかります。

ときどき部屋の隅っこをチェックしてたまっているほこりを捨てるだけで部屋はそれなりにきれいになります。

でも部屋はこまめに掃除したほうがよいでしょう。これは生活態度の問題(コンピュータでは計算できない)です。(苦笑)

# 参照文献

+ パソコンで遊ぶ物理シミュレーション 量子力学からカオスまでが目で見える ブルーバックスＢ‐９２４ / 神原武志，樫森与志喜，内藤正美
+ FACOM FORTRAN 解説編 / 富士通
+ GNU Fortranマニュアル