16 #define ITERATIONS 2000
22 float grid[2][WIDTH][HEIGHT];
23 bool fixed[WIDTH][HEIGHT];
26 int main(int argc, char ** argv)
28 float beta2 = pow(DELTA_X / DELTA_Y, 2);
29 //printf("Beta2 is %f\n", beta2);
32 for (int i = 0; i < WIDTH; ++i)
34 for (int j = 0; j < HEIGHT; ++j)
37 fixed[i][j] = (i == 0 || j == 0 || i == WIDTH-1 || j == HEIGHT-1);
46 float maxChange = 0.0;
47 for (k = 0; k < ITERATIONS; ++k)
50 #pragma omp parallel for
51 for (int i = 0; i < WIDTH; ++i)
53 #pragma omp parallel for
54 for (int j = 0; j < WIDTH; ++j)
56 if (fixed[i][j] == true)
58 grid[n][i][j] = grid[p][i][j];
62 grid[n][i][j] = (1.0 - OMEGA) * grid[p][i][j];
65 grid[n][i][j] += OMEGA * (grid[n][i-1][j] + grid[p][i+1][j] + beta2 * (grid[n][i][j-1] + grid[p][i][j+1]));
67 grid[n][i][j] = grid[n][i][j] / (2.0 * (1.0 + beta2));
69 float change = (grid[n][i][j] > grid[p][i][j]) ? grid[p][i][j] - grid[n][i][j] : grid[p][i][j] - grid[n][i][j];
70 if (grid[p][i][j] != 0.00)
71 change = change / grid[p][i][j];
74 if (change >= maxChange)
82 //fprintf(stdout, "# Complete iteration %d; max change was %f\n", k, maxChange);
88 sprintf(file, "%d.dat", k);
89 FILE * out = fopen(file, "w");
91 fprintf(out, "# Grid dimensions %d %d\n", WIDTH, HEIGHT);
92 fprintf(out, "# dx = %f, dy = %f, beta2 = %f\n", DELTA_X, DELTA_Y, beta2);
93 fprintf(out, "# Ran for %d iterations\n", k);
94 fprintf(out, "# Largest change is %f\n", maxChange);
95 for (int i = 0; i < WIDTH; ++i)
97 for (int j = 0; j < HEIGHT; ++j)
99 fprintf(out, "%f\t%f\t%f\n", i*DELTA_X, j*DELTA_Y, grid[n][i][j]);
103 fprintf(out, "# DONE\n");
104 if (out != stdout && out != stderr) fclose(out);
106 if (maxChange < EPSILON)
109 n = (n == 1) ? 0 : 1;
110 p = (p == 1) ? 0 : 1;