thesis/simulation/egun/main.c~

   1 #include <stdlib.h>
   2 #include <stdio.h>
   3 #include <math.h>
   4 #include <string.h>
   5 #include <stdbool.h>
   6
   7 #include <CImg.h>
   8
   9 #include <omp.h>
  10
  11 #define DELTA_X 1.0
  12 #define DELTA_Y 1.0
  13
  14 #define WIDTH 100
  15 #define HEIGHT 100
  16 #define ITERATIONS 2000
  17 #define EPSILON 0.00
  18
  19 #define OMEGA 1.0
  20
  21
  22 float grid[2][WIDTH][HEIGHT];
  23 bool fixed[WIDTH][HEIGHT];
  24
  25
  26 int main(int argc, char ** argv)
  27 {
  28         float beta2 = pow(DELTA_X / DELTA_Y, 2);
  29         //printf("Beta2 is %f\n", beta2);
  30         omp_set_nested(1);
  31
  32         for (int i = 0; i < WIDTH; ++i)
  33         {
  34                 for (int j = 0; j < HEIGHT; ++j)
  35                 {
  36                         grid[0][i][j] = 0;
  37                         fixed[i][j] = (i == 0 || j == 0 || i == WIDTH-1 || j == HEIGHT-1);
  38                 }
  39         }
  40
  41         // Add the electrodes
  42
  43
  44         int n = 1; int p = 0;
  45         int k = 0;
  46         float maxChange = 0.0;
  47         for (k = 0; k < ITERATIONS; ++k)
  48         {
  49                 maxChange = 0.0;
  50                 #pragma omp parallel for
  51                 for (int i = 0; i < WIDTH; ++i)
  52                 {
  53                         #pragma omp parallel for
  54                         for (int j = 0; j < WIDTH; ++j)
  55                         {
  56                                 if (fixed[i][j] == true)
  57                                 {
  58                                         grid[n][i][j] = grid[p][i][j];
  59                                         continue;
  60                                 }
  61
  62                                 grid[n][i][j] = (1.0 - OMEGA) * grid[p][i][j];
  63
  64
  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]));
  66
  67                                 grid[n][i][j] = grid[n][i][j] / (2.0 * (1.0 + beta2));
  68
  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];
  72                                 else
  73                                         change = 1.0;
  74                                 if (change >= maxChange)
  75                                         maxChange = change;
  76
  77
  78
  79                         }
  80                 }
  81
  82                 //fprintf(stdout, "# Complete iteration %d; max change was %f\n", k, maxChange);
  83
  84
  85
  86
  87                 char file[100];
  88                 sprintf(file, "%d.dat", k);
  89                 FILE * out = fopen(file, "w");
  90
  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)
  96                 {
  97                         for (int j = 0; j < HEIGHT; ++j)
  98                         {
  99                                 fprintf(out, "%f\t%f\t%f\n", i*DELTA_X, j*DELTA_Y, grid[n][i][j]);
 100                         }
 101                         fprintf(out, "\n");
 102                 }
 103                 fprintf(out, "# DONE\n");
 104                 if (out != stdout && out != stderr) fclose(out);
 105
 106                 if (maxChange < EPSILON)
 107                         break;
 108
 109                 n = (n == 1) ? 0 : 1;
 110                 p = (p == 1) ? 0 : 1;
 111         }
 112
 113 }
 114
 115 //EOF