Only do sor sweeps on fluid cells

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

@@ -4,6 +4,7 @@
  * 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"
 
@@ -13,19 +14,32 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     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 solveSOR(Solver* solver, double* p, double* rhs)
+void solveSOR(Solver* s, double* p, double* rhs)
 {
-    int imax      = solver->grid->imax;
+    int imax      = s->grid->imax;
-    int jmax      = solver->grid->jmax;
+    int jmax      = s->grid->jmax;
-    double eps    = solver->eps;
+    double eps    = s->eps;
-    int itermax   = solver->itermax;
+    int itermax   = s->itermax;
-    double dx2    = solver->grid->dx * solver->grid->dx;
+    double dx2    = s->grid->dx * s->grid->dx;
-    double dy2    = solver->grid->dy * solver->grid->dy;
+    double dy2    = s->grid->dy * s->grid->dy;
     double idx2   = 1.0 / dx2;
     double idy2   = 1.0 / dy2;
-    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
     double epssq  = eps * eps;
     int it        = 0;
     double res    = 1.0;
@@ -55,7 +69,7 @@ void solveSOR(Solver* solver, double* p, double* rhs)
             P(imax + 1, j) = P(imax, j);
         }
 
-        res = res / (double)(imax * jmax);
+        res = res / (double)(s->totalFluidCells);
 #ifdef DEBUG
         printf("%d Residuum: %e\n", it, res);
 #endif
@@ -80,6 +94,7 @@ void solve(Solver* solver, double* p, double* rhs)
     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;
 
@@ -92,14 +107,16 @@ void solve(Solver* solver, double* p, double* rhs)
 
             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 + 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;

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

@@ -17,6 +17,7 @@ typedef struct {
     double eps, omega, rho;
     int itermax;
     int levels;
+    int totalFluidCells;
     double **r, **e;
 } Solver;