Add solver variants

BasicSolver/2D-mpi-v3/src/solver.c

@@ -428,6 +428,7 @@ int solve(Solver* solver)
     double dy2    = solver->dy * solver->dy;
     double idx2   = 1.0 / dx2;
     double idy2   = 1.0 / dy2;
+    // identical to 1/((2/dx2)+(2/dy2))
     double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
     double* p     = solver->p;
     double* rhs   = solver->rhs;

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

@@ -247,68 +247,6 @@ void solveRB(Solver* solver)
 #endif
 }
 
-void solveRBA(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 omega  = solver->omega;
-    double* p     = solver->p;
-    double* rhs   = solver->rhs;
-    double epssq  = eps * eps;
-    int it        = 0;
-    double anorm = 1.0, anormf = 0.0;
-    int pass, jsw, isw;
-
-    while ((anorm >= eps * anormf) && (it < itermax)) {
-        anorm = 0.0;
-        jsw   = 1;
-
-        for (pass = 0; pass < 2; pass++) {
-            isw = jsw;
-
-            for (int j = 1; j < jmax + 1; j++) {
-                for (int i = isw; i < imax + 1; i += 2) {
-
-                    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);
-
-                    anorm += fabs(r);
-                    P(i, j) -= (omega * factor * r);
-                }
-                isw = 3 - isw;
-            }
-            jsw = 3 - jsw;
-        }
-
-        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);
-        }
-
-#ifdef DEBUG
-        printf("%d Residuum: %e\n", it, res);
-#endif
-        it++;
-    }
-
-#ifdef VERBOSE
-    printf("Solver took %d iterations to reach %f\n", it, anorm);
-#endif
-}
-
 static double maxElement(Solver* solver, double* m)
 {
     int size      = (solver->imax + 2) * (solver->jmax + 2);

PoissonSolver/2D-seq/poisson.par

@@ -16,6 +16,7 @@ jmax       200		# number of interior cells in y-direction
 
 itermax  10000		# maximal number of pressure iteration in one time step
 eps      0.000001	# stopping tolerance for pressure iteration
-omg      1.9		# relaxation parameter for SOR iteration
+# omg      0.9999		# relaxation parameter for SOR iteration
+omg      1.93		# relaxation parameter for SOR iteration
 
 #===============================================================================

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

@@ -4,18 +4,19 @@
  * Use of this source code is governed by a MIT-style
  * license that can be found in the LICENSE file.
  */
-#include <stdlib.h>
-#include <stdio.h>
-#include <unistd.h>
-#include <limits.h>
 #include <float.h>
+#include <limits.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
 
 #include "parameter.h"
 #include "solver.h"
-
+#include "timing.h"
 
 int main(int argc, char** argv)
 {
+    double startTime, endTime;
     Parameter params;
     Solver solver;
     initParameter(&params);
@@ -29,7 +30,10 @@ int main (int argc, char** argv)
     printParameter(&params);
 
     initSolver(&solver, &params, 2);
+    startTime = getTimeStamp();
     solve(&solver);
+    endTime = getTimeStamp();
+    printf("Solution took %.2fs\n", endTime - startTime);
     writeResult(&solver);
 
     return EXIT_SUCCESS;

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

@@ -4,13 +4,13 @@
  * Use of this source code is governed by a MIT style
  * license that can be found in the LICENSE file.
  */
-#include "stdlib.h"
-#include "stdio.h"
 #include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
 
-#include "solver.h"
-#include "parameter.h"
 #include "allocate.h"
+#include "parameter.h"
+#include "solver.h"
 
 #define PI        3.14159265358979323846
 #define P(i, j)   p[(j) * (imax + 2) + (i)]
@@ -58,35 +58,31 @@ void initSolver(Solver *solver, Parameter *params, int problem)
     }
 }
 
-int solve(Solver *solver)
+void solve(Solver* solver)
 {
-    double r;
-    int it = 0;
-    double res;
-
     int imax      = solver->imax;
     int jmax      = solver->jmax;
     double eps    = solver->eps;
-    double omega = solver->omega;
     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 = omega * 0.5 * (dx2*dy2) / (dx2+dy2);
+    double factor = solver->omega * 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;
 
-    res = eps + 1.0;
+    while ((res >= epssq) && (it < itermax)) {
-
-    while((res >= eps) && (it < itermax)) {
         res = 0.0;
 
         for (int j = 1; j < jmax + 1; j++) {
             for (int i = 1; i < imax + 1; i++) {
 
-                r = RHS(i,j) - ((P(i-1,j) - 2.0 * P(i,j) + P(i+1,j)) * idx2 +
+                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);
@@ -104,21 +100,139 @@ int solve(Solver *solver)
             P(imax + 1, j) = P(imax, j);
         }
 
-        it++;
+        res = res / (double)(imax * jmax);
-
-        //Residuum
-        res = sqrt(res / (imax*jmax));
 #ifdef DEBUG
         printf("%d Residuum: %e\n", it, res);
 #endif
+        it++;
     }
 
-    printf("Solver took %d iterations\n",it);
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
-    if( res < eps ){
-        return 1;
-    } else{
-        return 0;
 }
+
+void solveRB(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 = solver->omega * 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;
+    int pass, jsw, isw;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+        jsw = 1;
+
+        for (pass = 0; pass < 2; pass++) {
+            isw = jsw;
+
+            for (int j = 1; j < jmax + 1; j++) {
+                for (int i = isw; i < imax + 1; i += 2) {
+
+                    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);
+                }
+                isw = 3 - isw;
+            }
+            jsw = 3 - jsw;
+        }
+
+        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, res);
+#endif
+        it++;
+    }
+
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+}
+
+void solveRBA(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 rho    = solver->omega;
+    double* p     = solver->p;
+    double* rhs   = solver->rhs;
+    double epssq  = eps * eps;
+    int it        = 0;
+    double res    = 1.0;
+    int pass, jsw, isw;
+    double omega = 1.0;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+        jsw = 1;
+
+        for (pass = 0; pass < 2; pass++) {
+            isw = jsw;
+
+            for (int j = 1; j < jmax + 1; j++) {
+                for (int i = isw; i < imax + 1; i += 2) {
+
+                    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) -= (omega * factor * r);
+                    res += (r * r);
+                }
+                isw = 3 - isw;
+            }
+            jsw   = 3 - jsw;
+            omega = (it == 0 && pass == 0 ? 1.0 / (1.0 - 0.5 * rho * rho)
+                                          : 1.0 / (1.0 - 0.25 * rho * rho * omega));
+        }
+
+        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 Omega: %e\n", it, res, omega);
+#endif
+        it++;
+    }
+
+    printf("Final omega: %f\n", omega);
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
 }
 
 void writeResult(Solver* solver)

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

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