clan/tests/swim.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4
   5 #include <assert.h>
   6
   7 #include "decls.h"
   8 #include "util.h"
   9
  10 double t_start, t_end;
  11
  12 main ()
  13 {
  14     int t, i, j;
  15     int N3 = P_N3 ;
  16     int M = P_M ;
  17     int N = P_N ;
  18     int ALPHA = P_ALPHA ;
  19     int DX = P_DX ;
  20     int DY = P_DY ;
  21     int TDT = P_TDT ;
  22     int DT = P_DT  ;
  23
  24     init();
  25
  26 #ifdef PERFCTR
  27     PERF_INIT;
  28 #endif
  29
  30     IF_TIME(t_start = rtclock());
  31
  32
  33     /* pluto start (M,N,N3) */
  34 #pragma scop
  35
  36     for (t=0; t<N3; t++) {
  37
  38         FSDX = 4/DX;
  39         FSDY = 4/DY;
  40
  41         for (i=0; i<M; i++) {
  42             for (j=0; j<N; j++) {
  43                 CU[i+1][j] = 0.5*(P[i+1][j]+P[i][j])*U[i+1][j];
  44                 CV[i][j+1] = 0.5*(P[i][j+1]+P[i][j])*V[i][j+1];
  45                 Z[i+1][j+1] = (FSDX*(V[i+1][j+1]-V[i][j+1])-FSDY*(U[i+1][j+1]
  46                             -U[i+1][j]))/(P[i][j]+P[i+1][j]+P[i+1][j+1]+P[i][j+1]);
  47                 H[i][j] = P[i][j]+0.25*(U[i+1][j]*U[i+1][j]+U[i][j]*U[i][j]
  48                         +V[i][j+1]*V[i][j+1]+V[i][j]*V[i][j]);
  49             }
  50         }
  51
  52         for (j=0; j<N; j++)  {
  53             CU[0][j] = CU[M+1][j];
  54             CV[M][j+1] = CV[0][j+1];
  55             Z[0][j+1] = Z[M][j+1];
  56             H[M][j] = H[0][j];
  57         }
  58
  59         for (i=0; i<M; i++) {
  60             CU[i+1][N] = CU[i+1][0];
  61             CV[i][0] = CV[i][N];
  62             Z[i+1][0] = Z[i+1][N];
  63             H[i][N] = H[i][0];
  64         }
  65
  66
  67         CU[0][N] = CU[M][0];
  68         CV[M][0] = CV[0][N];
  69         Z[0][0] = Z[M][N];
  70         H[M][N] = H[0][0];
  71
  72         TDTS8 = TDT/8;
  73         TDTSDX = TDT/DX;
  74         TDTSDY = TDT/DY;
  75
  76         for (i=0; i<M; i++) {
  77             for (j=0; j<N; j++) {
  78                 UNEW[i+1][j] = UOLD[i+1][j]+
  79                     TDTS8*(Z[i+1][j+1]+Z[i+1][j])*(CV[i+1][j+1]+CV[i][j+1]+CV[i][j]
  80                             +CV[i+1][j])-TDTSDX*(H[i+1][j]-H[i][j]);
  81                 VNEW[i][j+1] = VOLD[i][j+1]-TDTS8*(Z[i+1][j+1]+Z[i][j+1])
  82                     *(CU[i+1][j+1]+CU[i][j+1]+CU[i][j]+CU[i+1][j])
  83                     -TDTSDY*(H[i][j+1]-H[i][j]);
  84                 PNEW[i][j] = POLD[i][j]-TDTSDX*(CU[i+1][j]-CU[i][j])
  85                     -TDTSDY*(CV[i][j+1]-CV[i][j]);
  86             }
  87         }
  88
  89         for (j=0; j<N; j++) {
  90             UNEW[0][j] = UNEW[M][j];
  91             VNEW[M][j+1] = VNEW[0][j+1];
  92             PNEW[M][j] = PNEW[0][j];
  93         }
  94
  95         for (i=0; i<M; i++) {
  96             UNEW[i+1][N] = UNEW[i+1][0];
  97             VNEW[i][0] = VNEW[i][N];
  98             PNEW[i][N] = PNEW[i][0];
  99         }
 100
 101         UNEW[0][N] = UNEW[M][0];
 102         VNEW[M][0] = VNEW[0][N];
 103         PNEW[M][N] = PNEW[0][0];
 104         time = time + DT;
 105
 106
 107         for (i=0; i<M; i++) {
 108             for (j=0; j<N; j++) {
 109                 UOLD[i][j] = U[i][j]+ALPHA*(UNEW[i][j]-2*U[i][j]+UOLD[i][j]);
 110                 VOLD[i][j] = V[i][j]+ALPHA*(VNEW[i][j]-2*V[i][j]+VOLD[i][j]);
 111                 POLD[i][j] = P[i][j]+ALPHA*(PNEW[i][j]-2*P[i][j]+POLD[i][j]);
 112                 U[i][j] = UNEW[i][j];
 113                 V[i][j] = VNEW[i][j];
 114                 P[i][j] = PNEW[i][j];
 115             }
 116         }
 117
 118         for (j=0; j<N; j++) {
 119             UOLD[M][j] = UOLD[0][j];
 120             VOLD[M][j] = VOLD[0][j];
 121             POLD[M][j] = POLD[0][j];
 122             U[M][j] = U[0][j];
 123             V[M][j] = V[0][j];
 124             P[M][j] = P[0][j];
 125         }
 126
 127         for (i=0; i<M; i++) {
 128             UOLD[i][N] = UOLD[i][0];
 129             VOLD[i][N] = VOLD[i][0];
 130             POLD[i][N] = POLD[i][0];
 131             U[i][N] = U[i][0];
 132             V[i][N] = V[i][0];
 133             P[i][N] = P[i][0];
 134         }
 135
 136         UOLD[M][N] = UOLD[0][0];
 137         VOLD[M][N] = VOLD[0][0];
 138         POLD[M][N] = POLD[0][0];
 139         U[M][N] = U[0][0];
 140         V[M][N] = V[0][0];
 141         P[M][N] = P[0][0];
 142     }
 143 #pragma endscop
 144     /* pluto end */
 145
 146     IF_TIME(t_end = rtclock());
 147     IF_TIME(printf("%0.6lfs\n", t_end - t_start));
 148
 149 #ifdef TEST
 150     print();
 151 #endif
 152 }
 153