Completed porting, fixing bugs and testing

BasicSolver/2D-mpi-v2/canal.par

@@ -41,7 +41,7 @@ tau      0.5     # safety factor for time stepsize control (<0 constant delt)
 
 itermax       500       # maximal number of pressure iteration in one time step
 eps           0.00001   # stopping tolerance for pressure iteration
-rho           0.9999
+rho           0.52
 omg           1.8       # relaxation parameter for SOR iteration
 gamma         0.9       # upwind differencing factor gamma
 #===============================================================================

BasicSolver/2D-mpi-v2/dcavity.par

@@ -41,7 +41,7 @@ tau     0.5		# safety factor for time stepsize control (<0 constant delt)
 
 itermax  1000		# maximal number of pressure iteration in one time step
 eps      0.001		# stopping tolerance for pressure iteration
+rho      0.5		
 omg      1.7		# relaxation parameter for SOR iteration
-rho      0.9999		# 
 gamma    0.9		# upwind differencing factor gamma
 #===============================================================================

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

@@ -69,15 +69,8 @@ int main (int argc, char** argv)
         solve(&solver);
         adaptUV(&solver);
         t += solver.dt;
-
-#ifdef VERBOSE
-        if (rank == 0) {
-            printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
         }
-#else
-        printProgress(t);
-#endif
-    }
+
         break;
     case RB:
         printf("Red-black SOR\n");
@@ -93,15 +86,8 @@ int main (int argc, char** argv)
         solveRB(&solver);
         adaptUV(&solver);
         t += solver.dt;
-
-#ifdef VERBOSE
-        if (rank == 0) {
-            printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
         }
-#else
-        printProgress(t);
-#endif
-    }
+
         break;
     case RBA:
         printf("Red-black SOR with acceleration\n");
@@ -117,6 +103,10 @@ int main (int argc, char** argv)
         solveRBA(&solver);
         adaptUV(&solver);
         t += solver.dt;
+        }
+
+        break;
+    }
 
 #ifdef VERBOSE
         if (rank == 0) {
@@ -125,9 +115,7 @@ int main (int argc, char** argv)
 #else
         printProgress(t);
 #endif
-    }
-        break;
-    }
+
     E = getTimeStamp();
     stopProgress();
     if (rank == 0) {

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

@@ -423,8 +423,8 @@ void initSolver(Solver* solver, Parameter* params)
         &solver->jNeighbours[1]);
     MPI_Cart_get(solver->comm, NDIMS, solver->dims, periods, solver->coords);
 
-    solver->imaxLocal = sizeOfRank(solver->rank, dims[IDIM], solver->imax);
-    solver->jmaxLocal = sizeOfRank(solver->rank, dims[JDIM], solver->jmax);
+    solver->imaxLocal = sizeOfRank(solver->coords[IDIM], dims[IDIM], solver->imax);
+    solver->jmaxLocal = sizeOfRank(solver->coords[JDIM], dims[JDIM], solver->jmax);
 
     MPI_Type_contiguous(solver->imaxLocal, MPI_DOUBLE, &solver->jBufferType);
     MPI_Type_commit(&solver->jBufferType);
@@ -585,21 +585,28 @@ int solveRB(Solver* solver)
     double epssq  = eps * eps;
     int it        = 0;
     double res    = 1.0;
+    int pass, jsw, isw;
 
     while ((res >= epssq) && (it < itermax)) {
-        res = 0.0;
-        exchange(solver, p);
+        jsw = 1;
 
-        for (int j = 1; j < jmaxLocal + 1; j++) {
-            for (int i = 1; i < imaxLocal + 1; i++) {
+        for (pass = 0; pass < 2; pass++) {
+            isw = jsw;
+            exchange(solver, p);
 
-                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);
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                for (int i = isw; i < imaxLocal + 1; i += 2) {
 
-                P(i, j) -= (factor * r);
-                res += (r * r);
+                    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;
         }
 
         if (solver->coords[JDIM] == 0) { // set bottom bc
@@ -660,27 +667,39 @@ int solveRBA(Solver* solver)
     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 factor = 0.5 * (dx2 * dy2) / (dx2 + dy2);
     double* p     = solver->p;
     double* rhs   = solver->rhs;
+    double rho    = solver->rho;
     double epssq  = eps * eps;
     int it        = 0;
     double res    = 1.0;
+    double omega  = 1.0;
+    int pass, jsw, isw;
 
     while ((res >= epssq) && (it < itermax)) {
         res = 0.0;
-        exchange(solver, p);
+        jsw = 1;
 
-        for (int j = 1; j < jmaxLocal + 1; j++) {
-            for (int i = 1; i < imaxLocal + 1; i++) {
+        for (pass = 0; pass < 2; pass++) {
+            isw = jsw;
+            exchange(solver, p);
 
-                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);
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                for (int i = isw; i < imaxLocal + 1; i += 2) {
 
-                P(i, j) -= (factor * r);
-                res += (r * r);
+                    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));
         }
 
         if (solver->coords[JDIM] == 0) { // set bottom bc