Add Jacobi and Gauss-Seidel. Add target and infrastructure to generate Videos

.gitignore

@@ -53,3 +53,5 @@ Module.symvers
 Mkfile.old
 dkms.conf
 
+*.mp4
+

PoissonSolver/2D-seq/Makefile

@@ -10,6 +10,7 @@ TARGET	   = exe-$(TAG)
 BUILD_DIR  = ./$(TAG)
 SRC_DIR    = ./src
 MAKE_DIR   = ./
+VIDEO_DIR  = ./video
 Q         ?= @
 
 #DO NOT EDIT BELOW
@@ -22,7 +23,7 @@ ASM       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s,$(wildcard $(SRC_DIR)/*.
 OBJ       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o,$(wildcard $(SRC_DIR)/*.c))
 CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
 
-${TARGET}: $(BUILD_DIR) $(OBJ)
+${TARGET}: $(BUILD_DIR) $(VIDEO_DIR) $(OBJ)
 	$(info ===>  LINKING  $(TARGET))
 	$(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS)
 
@@ -35,17 +36,22 @@ $(BUILD_DIR)/%.s:  %.c
 	$(info ===>  GENERATE ASM  $@)
 	$(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@
 
-.PHONY: clean distclean tags info asm
+.PHONY: clean distclean mkvid tags info asm
 
 clean:
 	$(info ===>  CLEAN)
 	@rm -rf $(BUILD_DIR)
+	@rm -rf $(VIDEO_DIR)
 	@rm -f tags
 
 distclean: clean
 	$(info ===>  DIST CLEAN)
 	@rm -f $(TARGET)
 
+mkvid:
+	gnuplot animate.plot
+	ffmpeg -framerate 10 -i ./video/%04d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p out.mp4
+
 info:
 	$(info $(CFLAGS))
 	$(Q)$(CC) $(VERSION)
@@ -59,4 +65,7 @@ tags:
 $(BUILD_DIR):
 	@mkdir $(BUILD_DIR)
 
+$(VIDEO_DIR):
+	@mkdir $(VIDEO_DIR)
+
 -include $(OBJ:.o=.d)

PoissonSolver/2D-seq/animate.plot

@@ -4,12 +4,13 @@ set grid
 set hidden3d
 set xrange [0:40]
 set yrange [0:40]
-set zrange [-2:2]
 
-input(n) = sprintf("p-%d.dat", n)
+input(n) = sprintf("./video/%04d.dat", n)
-output(n) = sprintf("%03d.png", n)
+output(n) = sprintf("./video/%04d.png", n)
+title(n) = sprintf("Iteration %d", n)
 
-do for [i=1:50] {
+do for [i=0:400] {
-    set output output(i)
+  set output output(i)
+  set title title(i)
   splot input(i) matrix using 1:2:3 with lines
 }

PoissonSolver/2D-seq/src/main.c

@@ -17,7 +17,7 @@
     LIKWID_MARKER_STOP(#tag);                                                            \
     endTime = getTimeStamp();
 
-enum VARIANT { SOR = 1, RB, RBA };
+enum VARIANT { JACOBI = 1, GS, SOR, RB, RBA };
 
 int main(int argc, char** argv)
 {
@@ -43,12 +43,20 @@ int main(int argc, char** argv)
     }
 
     initSolver(&solver, &params, 2);
-    writeResult(&solver, "p-init.dat");
+    writeResult(&solver, "./video/0000.dat");
 
     switch (variant) {
+    case JACOBI:
+        printf("JACOBI\n");
+        LIKWID_PROFILE("JACOBI", solveJacobi);
+        break;
+    case GS:
+        printf("GAUSS-SEIDEL\n");
+        LIKWID_PROFILE("GS", solveGS);
+        break;
     case SOR:
         printf("Plain SOR\n");
-        LIKWID_PROFILE("SOR", solve);
+        LIKWID_PROFILE("SOR", solveSOR);
         break;
     case RB:
         printf("Red-black SOR\n");

PoissonSolver/2D-seq/src/solver.c

@@ -10,9 +10,11 @@
 #include "parameter.h"
 #include "solver.h"
 
-#define PI        3.14159265358979323846
+#define PI         3.14159265358979323846
-#define P(i, j)   p[(j) * (imax + 2) + (i)]
+#define P(i, j)    p[(j) * (imax + 2) + (i)]
-#define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
+#define Pnew(i, j) pNew[(j) * (imax + 2) + (i)]
+#define Pold(i, j) pOld[(j) * (imax + 2) + (i)]
+#define RHS(i, j)  rhs[(j) * (imax + 2) + (i)]
 
 void initSolver(Solver* solver, Parameter* params, int problem)
 {
@@ -57,7 +59,135 @@ void initSolver(Solver* solver, Parameter* params, int problem)
     }
 }
 
-void solve(Solver* solver)
+static inline void swapPtr(double** p1, double** p2)
+{
+    double* tmp = *p1;
+    *p1         = *p2;
+    *p2         = tmp;
+}
+
+void solveJacobi(Solver* solver)
+{
+    int imax        = solver->imax;
+    int jmax        = solver->jmax;
+    double eps      = solver->eps;
+    int itermax     = solver->itermax;
+    double dx2      = solver->dx * solver->dx;
+    double dy2      = solver->dy * solver->dy;
+    double idx2     = 1.0 / dx2;
+    double idy2     = 1.0 / dy2;
+    double factor   = 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double* p       = solver->p;
+    double* rhs     = solver->rhs;
+    double epssq    = eps * eps;
+    int it          = 0;
+    double res      = 1.0;
+    size_t bytesize = (imax + 2) * (jmax + 2) * sizeof(double);
+    double* pOld    = p;
+    double* pNew    = allocate(64, bytesize);
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+
+                double r = RHS(i, j) -
+                           ((Pold(i - 1, j) - 2.0 * Pold(i, j) + Pold(i + 1, j)) * idx2 +
+                               (Pold(i, j - 1) - 2.0 * Pold(i, j) + Pold(i, j + 1)) *
+                                   idy2);
+
+                Pnew(i, j) -= factor * r;
+                res += (r * r);
+            }
+        }
+
+        for (int i = 1; i < imax + 1; i++) {
+            Pnew(i, 0)        = Pnew(i, 1);
+            Pnew(i, jmax + 1) = Pnew(i, jmax);
+        }
+
+        for (int j = 1; j < jmax + 1; j++) {
+            Pnew(0, j)        = Pnew(1, j);
+            Pnew(imax + 1, j) = Pnew(imax, j);
+        }
+
+        swapPtr(&pNew, &pOld);
+        res = res / (double)(imax * jmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, res);
+#endif
+#ifdef ANIMATE
+        char filename[50];
+        sprintf(filename, "./video/%04d.dat", it);
+        writeResult(solver, filename);
+#endif
+        it++;
+    }
+
+    free(pNew);
+    printf("%d\n", it);
+}
+
+void solveGS(Solver* solver)
+{
+    int imax        = solver->imax;
+    int jmax        = solver->jmax;
+    double eps      = solver->eps;
+    int itermax     = solver->itermax;
+    double dx2      = solver->dx * solver->dx;
+    double dy2      = solver->dy * solver->dy;
+    double idx2     = 1.0 / dx2;
+    double idy2     = 1.0 / dy2;
+    double factor   = 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double* p       = solver->p;
+    double* rhs     = solver->rhs;
+    double epssq    = eps * eps;
+    int it          = 0;
+    double res      = 1.0;
+    size_t bytesize = (imax + 2) * (jmax + 2) * sizeof(double);
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+
+                double r = RHS(i, j) -
+                           ((P(i - 1, j) - 2.0 * P(i, j) + P(i + 1, j)) * idx2 +
+                               (P(i, j - 1) - 2.0 * P(i, j) + P(i, j + 1)) * idy2);
+
+                P(i, j) -= factor * r;
+                res += (r * r);
+            }
+        }
+
+        for (int i = 1; i < imax + 1; i++) {
+            P(i, 0)        = P(i, 1);
+            P(i, jmax + 1) = P(i, jmax);
+        }
+
+        for (int j = 1; j < jmax + 1; j++) {
+            P(0, j)        = P(1, j);
+            P(imax + 1, j) = P(imax, j);
+        }
+
+        res = res / (double)(imax * jmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, sqrt(res));
+#endif
+#ifdef ANIMATE
+        char filename[50];
+        sprintf(filename, "video/%04d.dat", it);
+        writeResult(solver, filename);
+#endif
+        it++;
+    }
+
+    printf("%d\n", it);
+}
+
+void solveSOR(Solver* solver)
 {
     int imax      = solver->imax;
     int jmax      = solver->jmax;
@@ -73,7 +203,6 @@ void solve(Solver* solver)
     double epssq  = eps * eps;
     int it        = 0;
     double res    = 1.0;
-    char filename[20];
 
     while ((res >= epssq) && (it < itermax)) {
         res = 0.0;
@@ -105,7 +234,8 @@ void solve(Solver* solver)
         printf("%d Residuum: %e\n", it, res);
 #endif
 #ifdef ANIMATE
-        sprintf(filename, "p-%d.dat", it);
+        char filename[50];
+        sprintf(filename, "./video/%04d.dat", it);
         writeResult(solver, filename);
 #endif
         it++;
@@ -169,7 +299,8 @@ void solveRB(Solver* solver)
         printf("%d Residuum: %e\n", it, res);
 #endif
 #ifdef ANIMATE
-        sprintf(filename, "p-%d.dat", it);
+        char filename[50];
+        sprintf(filename, "./video/%04d.dat", it);
         writeResult(solver, filename);
 #endif
         it++;
@@ -238,7 +369,8 @@ void solveRBA(Solver* solver)
         printf("%d Residuum: %e Omega: %e\n", it, res, omega);
 #endif
 #ifdef ANIMATE
-        sprintf(filename, "p-%d.dat", it);
+        char filename[50];
+        sprintf(filename, "video/%04d.dat", it);
         writeResult(solver, filename);
 #endif
         it++;
@@ -259,7 +391,7 @@ void writeResult(Solver* solver, char* filename)
     fp = fopen(filename, "w");
 
     if (fp == NULL) {
-        printf("Error!\n");
+        printf("Error opening file %s!\n", filename);
         exit(EXIT_FAILURE);
     }
 

PoissonSolver/2D-seq/src/solver.h

@@ -16,7 +16,9 @@ typedef struct {
 
 extern void initSolver(Solver*, Parameter*, int problem);
 extern void writeResult(Solver*, char*);
-extern void solve(Solver*);
+extern void solveJacobi(Solver*);
+extern void solveGS(Solver*);
+extern void solveSOR(Solver*);
 extern void solveRB(Solver*);
 extern void solveRBA(Solver*);
 #endif