Refactor and add red-black solver

BasicSolver/3D-seq/src/main.c

@@ -69,7 +69,7 @@ int main(int argc, char** argv)
 {
     double timeStart, timeStop;
     Parameter params;
-    Solver solver;
+    Solver s;
     initParameter(&params);
 
     if (argc != 2) {
@@ -79,30 +79,28 @@ int main(int argc, char** argv)
 
     readParameter(&params, argv[1]);
     printParameter(&params);
-    initSolver(&solver, &params);
+    initSolver(&s, &params);
 #ifndef VERBOSE
-    initProgress(solver.te);
+    initProgress(s.te);
 #endif
 
-    double tau = solver.tau;
-    double te  = solver.te;
+    double tau = s.tau;
+    double te  = s.te;
     double t   = 0.0;
-    int nt     = 0;
 
     timeStart = getTimeStamp();
     while (t <= te) {
-        if (tau > 0.0) computeTimestep(&solver);
-        setBoundaryConditions(&solver);
-        setSpecialBoundaryCondition(&solver);
-        computeFG(&solver);
-        computeRHS(&solver);
-        solve(&solver);
-        adaptUV(&solver);
-        t += solver.dt;
-        nt++;
+        if (tau > 0.0) computeTimestep(&s);
+        setBoundaryConditions(&s);
+        setSpecialBoundaryCondition(&s);
+        computeFG(&s);
+        computeRHS(&s);
+        solveRB(&s);
+        adaptUV(&s);
+        t += s.dt;
 
 #ifdef VERBOSE
-        printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
+        printf("TIME %f , TIMESTEP %f\n", t, s.dt);
 #else
         printProgress(t);
 #endif
@@ -115,8 +113,7 @@ int main(int argc, char** argv)
 
     double *pg, *ug, *vg, *wg;
 
-    size_t bytesize = solver.grid.imax * solver.grid.jmax * solver.grid.kmax *
-                      sizeof(double);
+    size_t bytesize = (size_t)(s.grid.imax * s.grid.jmax * s.grid.kmax) * sizeof(double);
 
     pg = allocate(64, bytesize);
     ug = allocate(64, bytesize);
@@ -124,8 +121,8 @@ int main(int argc, char** argv)
     wg = allocate(64, bytesize);
 
     createBulkArrays(&s, pg, ug, vg, wg);
-    VtkOptions opts = { .grid = solver.grid };
-    vtkOpen(&opts, solver.problem);
+    VtkOptions opts = { .grid = s.grid };
+    vtkOpen(&opts, s.problem);
     vtkScalar(&opts, "pressure", pg);
     vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
     vtkClose(&opts);

BasicSolver/3D-seq/src/progress.h

@@ -9,6 +9,6 @@
 
 extern void initProgress(double);
 extern void printProgress(double);
-extern void stopProgress();
+extern void stopProgress(void);
 
 #endif

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

@@ -221,6 +221,91 @@ void solve(Solver* s)
 #endif
 }
 
+void solveRB(Solver* s)
+{
+    int imax      = s->grid.imax;
+    int jmax      = s->grid.jmax;
+    int kmax      = s->grid.kmax;
+    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 dz2    = s->grid.dz * s->grid.dz;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double idz2   = 1.0 / dz2;
+    double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) /
+                    (dy2 * dz2 + dx2 * dz2 + dx2 * dy2);
+    double* p    = s->p;
+    double* rhs  = s->rhs;
+    double epssq = eps * eps;
+    int it       = 0;
+    double res   = 1.0;
+    int pass, ksw, jsw, isw;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+        ksw = 1;
+
+        for (pass = 0; pass < 2; pass++) {
+            jsw = ksw;
+
+            for (int k = 1; k < kmax + 1; k++) {
+                isw = jsw;
+                for (int j = 1; j < jmax + 1; j++) {
+                    for (int i = isw; i < imax + 1; i += 2) {
+
+                        double r =
+                            RHS(i, j, k) -
+                            ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 +
+                                (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
+                                    idy2 +
+                                (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
+                                    idz2);
+
+                        P(i, j, k) -= (factor * r);
+                        res += (r * r);
+                    }
+                    isw = 3 - isw;
+                }
+                jsw = 3 - jsw;
+            }
+            ksw = 3 - ksw;
+        }
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                P(i, j, 0)        = P(i, j, 1);
+                P(i, j, kmax + 1) = P(i, j, kmax);
+            }
+        }
+
+        for (int k = 1; k < kmax + 1; k++) {
+            for (int i = 1; i < imax + 1; i++) {
+                P(i, 0, k)        = P(i, 1, k);
+                P(i, jmax + 1, k) = P(i, jmax, k);
+            }
+        }
+
+        for (int k = 1; k < kmax + 1; k++) {
+            for (int j = 1; j < jmax + 1; j++) {
+                P(0, j, k)        = P(1, j, k);
+                P(imax + 1, j, k) = P(imax, j, k);
+            }
+        }
+
+        res = res / (double)(imax * jmax * kmax);
+#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
+}
+
 static double maxElement(Solver* s, double* m)
 {
     int size      = (s->grid.imax + 2) * (s->grid.jmax + 2) * (s->grid.kmax + 2);

BasicSolver/3D-seq/src/solver.h

@@ -34,6 +34,7 @@ typedef struct {
 extern void initSolver(Solver*, Parameter*);
 extern void computeRHS(Solver*);
 extern void solve(Solver*);
+extern void solveRB(Solver*);
 extern void normalizePressure(Solver*);
 extern void computeTimestep(Solver*);
 extern void setBoundaryConditions(Solver*);