Split red-balack and standard sor solvers

BasicSolver/2D-seq/config.mk

@@ -1,8 +1,8 @@
 # Supported: GCC, CLANG, ICC
 TAG ?= CLANG
 ENABLE_OPENMP ?= false
-# Supported: sor, mg
-SOLVER ?= sor
+# Supported: sor, rb, mg
+SOLVER ?= rb
 # Run in debug settings
 DEBUG ?= false
 

BasicSolver/2D-seq/dcavity.par

@@ -26,8 +26,8 @@ p_init    0.0		# initial value for pressure
 
 xlength    1.0		# domain size in x-direction
 ylength    1.0		# domain size in y-direction
-imax       100		# number of interior cells in x-direction
-jmax       100		# number of interior cells in y-direction
+imax       128		# number of interior cells in x-direction
+jmax       128		# number of interior cells in y-direction
 
 # Time Data:
 # ---------

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

@@ -0,0 +1,76 @@
+/*
+ * Copyright (C)  NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#include "solver.h"
+#include "util.h"
+
+void initSolver(Solver* s, Discretization* d, Parameter* p)
+{
+    s->grid    = &d->grid;
+    s->itermax = p->itermax;
+    s->eps     = p->eps;
+    s->omega   = p->omg;
+}
+
+void solve(Solver* solver, double* p, double* rhs)
+{
+    int imax      = solver->grid->imax;
+    int jmax      = solver->grid->jmax;
+    double eps    = solver->eps;
+    int itermax   = solver->itermax;
+    double dx2    = solver->grid->dx * solver->grid->dx;
+    double dy2    = solver->grid->dy * solver->grid->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    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++;
+    }
+
+#ifdef VERBOSE
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+#endif
+}

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

@@ -15,58 +15,6 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     s->omega   = p->omg;
 }
 
-void solveSOR(Solver* solver, double* p, double* rhs)
-{
-    int imax      = solver->grid->imax;
-    int jmax      = solver->grid->jmax;
-    double eps    = solver->eps;
-    int itermax   = solver->itermax;
-    double dx2    = solver->grid->dx * solver->grid->dx;
-    double dy2    = solver->grid->dy * solver->grid->dy;
-    double idx2   = 1.0 / dx2;
-    double idy2   = 1.0 / dy2;
-    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
-    double epssq  = eps * eps;
-    int it        = 0;
-    double res    = 1.0;
-
-    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, res);
-#endif
-        it++;
-    }
-
-#ifdef VERBOSE
-    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
-#endif
-}
-
 void solve(Solver* solver, double* p, double* rhs)
 {
     int imax      = solver->grid->imax;
@@ -81,17 +29,12 @@ void solve(Solver* solver, double* p, double* 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) {
+            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 +
@@ -100,9 +43,6 @@ void solve(Solver* solver, double* p, double* rhs)
                 P(i, j) -= (factor * r);
                 res += (r * r);
             }
-                isw = 3 - isw;
-            }
-            jsw = 3 - jsw;
         }
 
         for (int i = 1; i < imax + 1; i++) {

EnhancedSolver/2D-seq/config.mk

@@ -1,7 +1,7 @@
 # Supported: GCC, CLANG, ICC
 TAG ?= CLANG
 ENABLE_OPENMP ?= false
-# Supported: sor, mg
+# Supported: sor, rb, mg
 SOLVER ?= mg
 # Run in debug settings
 DEBUG ?= false

EnhancedSolver/2D-seq/src/solver-rb.c

@@ -0,0 +1,93 @@
+/*
+ * Copyright (C)  NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#include "grid.h"
+#include "solver.h"
+#include "util.h"
+
+void initSolver(Solver* s, Discretization* d, Parameter* p)
+{
+    s->grid    = &d->grid;
+    s->itermax = p->itermax;
+    s->eps     = p->eps;
+    s->omega   = p->omg;
+
+    Grid* g  = s->grid;
+    int imax = s->grid->imax;
+    int jmax = s->grid->jmax;
+
+    s->totalFluidCells = 0;
+    for (int j = 0; j < jmax + 2; j++) {
+        for (int i = 0; i < imax + 2; i++) {
+            if (gridIsFluid(g, i, j)) {
+                s->totalFluidCells++;
+            }
+        }
+    }
+}
+
+void solve(Solver* s, double* p, double* rhs)
+{
+    int imax      = s->grid->imax;
+    int jmax      = s->grid->jmax;
+    double eps    = s->eps;
+    int itermax   = s->itermax;
+    double dx2    = s->grid->dx * s->grid->dx;
+    double dy2    = s->grid->dy * s->grid->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double epssq  = eps * eps;
+    int it        = 0;
+    Grid* g       = s->grid;
+    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) {
+                    if (gridIsFluid(g, i, j)) {
+                        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)(s->totalFluidCells);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, res);
+#endif
+        it++;
+    }
+
+#ifdef VERBOSE
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+#endif
+}

EnhancedSolver/2D-seq/src/solver-sor.c

@@ -29,7 +29,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     }
 }
 
-void solveSOR(Solver* s, double* p, double* rhs)
+void solve(Solver* s, double* p, double* rhs)
 {
     int imax      = s->grid->imax;
     int jmax      = s->grid->jmax;
@@ -42,6 +42,7 @@ void solveSOR(Solver* s, double* p, double* rhs)
     double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
     double epssq  = eps * eps;
     int it        = 0;
+    Grid* g       = s->grid;
     double res    = 1.0;
 
     while ((res >= epssq) && (it < itermax)) {
@@ -49,7 +50,7 @@ void solveSOR(Solver* s, double* p, double* rhs)
 
         for (int j = 1; j < jmax + 1; j++) {
             for (int i = 1; i < imax + 1; i++) {
-
+                if (gridIsFluid(g, i, j)) {
                     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);
@@ -58,6 +59,7 @@ void solveSOR(Solver* s, double* p, double* rhs)
                     res += (r * r);
                 }
             }
+        }
 
         for (int i = 1; i < imax + 1; i++) {
             P(i, 0)        = P(i, 1);
@@ -80,66 +82,3 @@ void solveSOR(Solver* s, double* p, double* rhs)
     printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
 #endif
 }
-
-void solve(Solver* solver, double* p, double* rhs)
-{
-    int imax      = solver->grid->imax;
-    int jmax      = solver->grid->jmax;
-    double eps    = solver->eps;
-    int itermax   = solver->itermax;
-    double dx2    = solver->grid->dx * solver->grid->dx;
-    double dy2    = solver->grid->dy * solver->grid->dy;
-    double idx2   = 1.0 / dx2;
-    double idy2   = 1.0 / dy2;
-    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
-    double epssq  = eps * eps;
-    int it        = 0;
-    Grid* g       = solver->grid;
-    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) {
-                    if (gridIsFluid(g, i, j)) {
-                        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++;
-    }
-
-#ifdef VERBOSE
-    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
-#endif
-}