This is the thesis

[matches/honours.git] / thesis / simulation / egun / main.c~

diff --git a/thesis/simulation/egun/main.c~ b/thesis/simulation/egun/main.c~
new file mode 100644 (file)
index 0000000..ff92a23
--- /dev/null
+++ b/thesis/simulation/egun/main.c~
@@ -0,0 +1,115 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <stdbool.h>
+
+#include <CImg.h>
+
+#include <omp.h>
+
+#define DELTA_X 1.0
+#define DELTA_Y 1.0
+
+#define WIDTH 100
+#define HEIGHT 100
+#define ITERATIONS 2000
+#define EPSILON 0.00
+
+#define OMEGA 1.0
+
+
+float grid[2][WIDTH][HEIGHT];
+bool fixed[WIDTH][HEIGHT];
+
+
+int main(int argc, char ** argv)
+{
+       float beta2 = pow(DELTA_X / DELTA_Y, 2);
+       //printf("Beta2 is %f\n", beta2);
+       omp_set_nested(1);
+       
+       for (int i = 0; i < WIDTH; ++i)
+       {
+               for (int j = 0; j < HEIGHT; ++j)
+               {
+                       grid[0][i][j] = 0;
+                       fixed[i][j] = (i == 0 || j == 0 || i == WIDTH-1 || j == HEIGHT-1);
+               }
+       }
+
+       // Add the electrodes
+       
+
+       int n = 1; int p = 0;
+       int k = 0;
+       float maxChange = 0.0;
+       for (k = 0; k < ITERATIONS; ++k)
+       {
+               maxChange = 0.0;
+               #pragma omp parallel for
+               for (int i = 0; i < WIDTH; ++i)
+               {
+                       #pragma omp parallel for
+                       for (int j = 0; j < WIDTH; ++j)
+                       {
+                               if (fixed[i][j] == true)
+                               {
+                                       grid[n][i][j] = grid[p][i][j];
+                                       continue;
+                               }
+                               
+                               grid[n][i][j] = (1.0 - OMEGA) * grid[p][i][j];
+                       
+
+                               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]));
+                               
+                               grid[n][i][j] = grid[n][i][j] / (2.0 * (1.0 + beta2));
+
+                               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];
+                               if (grid[p][i][j] != 0.00)
+                                       change = change / grid[p][i][j];
+                               else
+                                       change = 1.0;
+                               if (change >= maxChange)
+                                       maxChange = change;
+                               
+                               
+
+                       }
+               }
+               
+               //fprintf(stdout, "# Complete iteration %d; max change was %f\n", k, maxChange);
+
+               
+       
+
+               char file[100];
+               sprintf(file, "%d.dat", k);
+               FILE * out = fopen(file, "w");
+
+               fprintf(out, "# Grid dimensions %d %d\n", WIDTH, HEIGHT);
+               fprintf(out, "# dx = %f, dy = %f, beta2 = %f\n", DELTA_X, DELTA_Y, beta2);
+               fprintf(out, "# Ran for %d iterations\n", k);
+               fprintf(out, "# Largest change is %f\n", maxChange);
+               for (int i = 0; i < WIDTH; ++i)
+               {
+                       for (int j = 0; j < HEIGHT; ++j)
+                       {
+                               fprintf(out, "%f\t%f\t%f\n", i*DELTA_X, j*DELTA_Y, grid[n][i][j]);
+                       }
+                       fprintf(out, "\n");
+               }
+               fprintf(out, "# DONE\n");
+               if (out != stdout && out != stderr) fclose(out);
+
+               if (maxChange < EPSILON)
+                       break;
+
+               n = (n == 1) ? 0 : 1;
+               p = (p == 1) ? 0 : 1;
+       }
+
+}
+
+//EOF