C のコードでの周期境界条件の剰余系による実装

#はじめに

周期境界条件の２次元の拡散方程式を陽解法で実装する場合、境界部分の処理(n ->0, -1 ->n-1)が意外と面倒である。この部分を剰余系 (mod, %) を使って実装してみる。

#include <stdio.h>
#include <stdlib.h>
# define n 100

int main()
{
  printf("n mod n : %d\n", n%n);
  printf("n+1 mod n : %d\n", (n+1)%n);
  printf("-1 mod n : %d\n", (-1)%n);
}

この出力は以下のようになる。

n mod n : 0
n+1 mod n : 1
-1 mod n : -1

インデックスの値が負になってしまうと都合が悪い。そこで、元々の値にnを加える。

#include <stdio.h>
#include <stdlib.h>
# define n 100

int main()
{
  printf("(n+n) mod n : %d\n", (n+n)%n);
  printf("(n+n+1) mod n : %d\n", (n+n+1)%n);
  printf("(n-1) mod n : %d\n", (n-1)%n);
}

この出力は以下のようになる。

(n+n) mod n : 0
(n+n+1) mod n : 1
(n-1) mod n : 99

２次元拡散方程式によるベンチマーク

実際に計算を行う際、条件分岐を用いるのと剰余を用いるのでは、剰余の方が割り算を多数行うので遅くなる可能性がある。ベンチマークを行ってみる。

//
// 2D diffusion equation
//  Comparison of periodic boundary condition implementation
//

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define _USE_MATH_DEFINES
#include <math.h>


#define n 100
#define totalLoop 100000

long i, j, x, y, t, start,end;
double du,dx,dt;

double u[n][n];
double uNew[n][n];

int main(void) {
  //
  //  Parameters
  //

  dx=1;
  du=1;
  dt=0.1;


  //
  // initial distribution
  //

  for(x=0;x<n;x++){
    for(y=0;y<n;y++){
      u[y][x] = 0;
    }
  }
  u[(int)(n/2)][(int)(n/2)]=1;




  //
  // calculation
  //


  start = clock();

  for(t=0;t<totalLoop;t++)
  {
    for(x=1;x<n-1;x++){
      for(y=1;y<n-1;y++){
        uNew[y][x] = u[y][x]+dt/dx/dx*
        (
          u[y+1][x]+
          u[y-1][x]+
          u[y][x+1]+
          u[y][x-1]-
          4*u[y][x]
        );
      }
    };

    //Calculating y boundary
    for(x=1;x<n-1;x++){
      uNew[0][x] = u[0][x]+dt/dx/dx*
      (
        u[1][x]+
        u[n-1][x]+
        u[0][x+1]+
        u[0][x-1]
        -4*u[0][x]
      );
      uNew[n-1][x] = u[n-1][x]+dt/dx/dx*
      (
        u[0][x]+
        u[n-2][x]+
        u[n-1][x+1]+
        u[n-1][x-1]
        -4*u[n-1][x]
      );
    };

    //calculating x boundary
    for(y=1;y<n-1;y++){
      uNew[y][0] = u[y][0]+dt/dx/dx*
      (
        u[y+1][0]+
        u[y-1][0]+
        u[y][1]+
        u[y][n-1]-
        4*u[y][0]
      );
      uNew[y][n-1] = u[y][n-1]+dt/dx/dx*
      (
        u[y+1][n-1]+
        u[y-1][n-1]+
        u[y][0]+
        u[y][n-2]-
        4*u[y][n-1]
      );
    };

    //Calculating four courners
    uNew[0][0] = u[0][0]+dt/dx/dx*
    (
      u[0][n-1]+
      u[0][1]+
      u[n-1][0]+
      u[1][0]-
      4*u[0][0]
    );

    uNew[0][n-1] = u[0][n-1]+dt/dx/dx*
    (
      u[0][n-2]+
      u[0][0]+
      u[n-1][n-1]+
      u[1][n-1]-
      4*u[0][n-1]
    );

    uNew[n-1][0] = u[n-1][0]+dt/dx/dx*
    (
      u[0][0]+
      u[n-2][0]+
      u[n-1][1]+
      u[n-1][n-1]-
      4*u[n-1][0]
    );
    uNew[n-1][n-1] = u[n-1][n-1]+dt/dx/dx*
    (
      u[n-1][0]+
      u[n-1][n-2]+
      u[0][n-1]+
      u[n-2][n-1]-
      4*u[n-1][n-1]
    );
  }

  end = clock();

  printf("n=%d\n",n);
  printf("Without MOD:　%lf s\n", (double)(end-start)/ CLOCKS_PER_SEC);


  start = clock();

  for(t=0;t<totalLoop;t++)
  {
    for(x=0;x<n;x++){
      for(y=0;y<n;y++){
        uNew[y][x] = u[y][x]+dt/dx/dx*
        (
          u[(y+1+n)%n][x]+
          u[(y-1+n)%n][x]+
          u[y][(x+1+n)%n]+
          u[y][(x-1+n)%n]-
          4*u[y][x]
        );
      }
    }
  }

  end = clock();

  printf("MOD system :　%lf s\n", (double)(end-start)/ CLOCKS_PER_SEC);

  return 0;


}

n=100のとき、出力は以下のようになる。

n=100
Without MOD: 1.820518 s
MOD system : 2.445765 s

2-30％のスピードダウンになっている。しかし、n=2^mとなるように選ぶと、

n=128
Without MOD: 4.366525 s
MOD system : 4.545593 s

数%しか差がなくなる。これは、コンパイラ内部で2^mの割り算はビットシフトを用いて高速に行っているからと推定される。コードのコンパクトさを考えると実用的な手法と思われる。