Fix solver and cleanup

BasicSolver/2D-mpi/config.mk

@@ -2,7 +2,7 @@
 TAG ?= CLANG
 ENABLE_MPI ?= true
 ENABLE_OPENMP ?= false
-COMM_TYPE ?= v3
+COMM_TYPE ?= v2
 
 #Feature options
 OPTIONS +=  -DARRAY_ALIGNMENT=64

BasicSolver/2D-mpi/src/comm-v2.c

@@ -138,43 +138,27 @@ int commIsBoundary(Comm* c, int direction)
 void commExchange(Comm* c, double* grid)
 {
 #ifdef _MPI
-    double* sbuf[NDIRS];
-    double* rbuf[NDIRS];
     MPI_Request requests[8];
     for (int i = 0; i < 8; i++)
         requests[i] = MPI_REQUEST_NULL;
 
-    rbuf[LEFT]   = grid + (c->imaxLocal + 2);
-    rbuf[RIGHT]  = grid + (c->imaxLocal + 2) + (c->imaxLocal + 1);
-    rbuf[BOTTOM] = grid + 1;
-    rbuf[TOP]    = grid + (c->jmaxLocal + 1) * (c->imaxLocal + 2) + 1;
-    sbuf[LEFT]   = grid + (c->imaxLocal + 2) + 1;
-    sbuf[RIGHT]  = grid + (c->imaxLocal + 2) + (c->imaxLocal);
-    sbuf[BOTTOM] = grid + (c->imaxLocal + 2) + 1;
-    sbuf[TOP]    = grid + (c->jmaxLocal) * (c->imaxLocal + 2) + 1;
-    // buf[0] = grid + (c->imaxLocal + 2);                          // recv left
-    // buf[1] = grid + (c->imaxLocal + 2) + 1;                      // send left
-    // buf[2] = grid + (c->imaxLocal + 2) + (c->imaxLocal + 1);     // recv right
-    // buf[3] = grid + (c->imaxLocal + 2) + (c->imaxLocal);         // send right
-    // buf[4] = grid + 1;                                           // recv bottom
-    // buf[5] = grid + (c->imaxLocal + 2) + 1;                      // send bottom
-    // buf[6] = grid + (c->jmaxLocal + 1) * (c->imaxLocal + 2) + 1; // recv top
-    // buf[7] = grid + (c->jmaxLocal) * (c->imaxLocal + 2) + 1;     // send top
-
     for (int i = 0; i < NDIRS; i++) {
+        double* sbuf = grid + c->sdispls[i];
+        double* rbuf = grid + c->rdispls[i];
+
         int tag = 0;
         if (c->neighbours[i] != MPI_PROC_NULL) {
             // printf("DEBUG: Rank %d - SendRecv with %d\n", c->rank, c->neighbours[i]);
             tag = c->neighbours[i];
         }
-        MPI_Irecv(rbuf[i],
+        MPI_Irecv(rbuf,
             1,
             c->bufferTypes[i],
             c->neighbours[i],
             tag,
             c->comm,
             &requests[i * 2]);
-        MPI_Isend(sbuf[i],
+        MPI_Isend(sbuf,
             1,
             c->bufferTypes[i],
             c->neighbours[i],
@@ -283,6 +267,15 @@ void commPartition(Comm* c, int jmax, int imax)
     c->bufferTypes[RIGHT]  = iBufferType;
     c->bufferTypes[BOTTOM] = jBufferType;
     c->bufferTypes[TOP]    = jBufferType;
+
+    c->rdispls[LEFT]   = (c->imaxLocal + 2);
+    c->rdispls[RIGHT]  = (c->imaxLocal + 2) + (c->imaxLocal + 1);
+    c->rdispls[BOTTOM] = 1;
+    c->rdispls[TOP]    = (c->jmaxLocal + 1) * (c->imaxLocal + 2) + 1;
+    c->sdispls[LEFT]   = (c->imaxLocal + 2) + 1;
+    c->sdispls[RIGHT]  = (c->imaxLocal + 2) + (c->imaxLocal);
+    c->sdispls[BOTTOM] = (c->imaxLocal + 2) + 1;
+    c->sdispls[TOP]    = (c->jmaxLocal) * (c->imaxLocal + 2) + 1;
 #else
     c->imaxLocal = imax;
     c->jmaxLocal = jmax;

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

@@ -127,7 +127,7 @@ void computeRHS(Solver* s)
     }
 }
 
-int solve(Solver* s)
+void solve(Solver* s)
 {
     int imax      = s->imax;
     int jmax      = s->jmax;
@@ -143,23 +143,29 @@ int solve(Solver* s)
     double* p     = s->p;
     double* rhs   = s->rhs;
     double epssq  = eps * eps;
-    int it        = 0;
+    int pass, jsw, isw;
-    double res    = 1.0;
+    int it     = 0;
-    commExchange(&s->comm, p);
+    double res = 1.0;
 
     while ((res >= epssq) && (it < itermax)) {
-        res = 0.0;
+        jsw = 1;
+        for (pass = 0; pass < 2; pass++) {
+            isw = jsw;
+            commExchange(&s->comm, p);
 
-        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int j = 1; j < jmaxLocal + 1; j++) {
-            for (int i = 1; i < imaxLocal + 1; i++) {
+                for (int i = isw; i < imaxLocal + 1; i += 2) {
 
-                double r = RHS(i, j) -
+                    double r = RHS(i, j) -
-                           ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
+                               ((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);
+                    P(i, j) -= (factor * r);
-                res += (r * r);
+                    res += (r * r);
+                }
+                isw = 3 - isw;
             }
+            jsw = 3 - jsw;
         }
 
         if (commIsBoundary(&s->comm, BOTTOM)) { // set bottom bc
@@ -201,11 +207,6 @@ int solve(Solver* s)
         printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
     }
 #endif
-    if (res < eps) {
-        return 0;
-    } else {
-        return 1;
-    }
 }
 
 static double maxElement(Solver* s, double* m)

BasicSolver/2D-mpi/src/solver.h

@@ -36,7 +36,7 @@ typedef struct {
 
 void initSolver(Solver*, Parameter*);
 void computeRHS(Solver*);
-int solve(Solver*);
+void solve(Solver*);
 void normalizePressure(Solver*);
 void computeTimestep(Solver*);
 void setBoundaryConditions(Solver*);