Porting working 2D-seq and 2D-mpi with MG

BasicSolver/2D-seq/canal.par

@@ -36,6 +36,13 @@ te      100.0   # final time
 dt      0.02    # time stepsize
 tau     0.5     # safety factor for time stepsize control (<0 constant delt)
 
+# Multigrid data:
+# ---------
+
+levels        3         # Multigrid levels
+presmooth     5         # Pre-smoothning iterations
+postsmooth    5         # Post-smoothning iterations
+
 # Pressure Iteration Data:
 # -----------------------
 
@@ -43,5 +50,4 @@ itermax  500       # maximal number of pressure iteration in one time step
 eps      0.00001   # stopping tolerance for pressure iteration
 omg      1.8       # relaxation parameter for SOR iteration
 gamma    0.9       # upwind differencing factor gamma
-levels   5         # Multigrid levels
 #===============================================================================

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       800		# number of interior cells in x-direction
+imax       128		# number of interior cells in x-direction
-jmax       800		# number of interior cells in y-direction
+jmax       128		# number of interior cells in y-direction
 
 # Time Data:
 # ---------
@@ -36,6 +36,13 @@ te      10.0		# final time
 dt      0.02	    # time stepsize
 tau     0.5	    	# safety factor for time stepsize control (<0 constant delt)
 
+# Multigrid data:
+# ---------
+
+levels        2         # Multigrid levels
+presmooth     20         # Pre-smoothning iterations
+postsmooth    5         # Post-smoothning iterations
+
 # Solver Data:
 # -----------------------
 
@@ -44,5 +51,4 @@ eps      0.001		# stopping tolerance for pressure iteration
 rho      0.5
 omg      1.7		# relaxation parameter for SOR iteration
 gamma    0.9		# upwind differencing factor gamma
-levels   5         # Multigrid levels
 #===============================================================================

BasicSolver/2D-seq/src/parameter.c

@@ -24,6 +24,8 @@ void initParameter(Parameter* param)
     param->gamma   = 0.9;
     param->tau     = 0.5;
     param->levels  = 5;
+    param->presmooth = 5;
+    param->postsmooth = 5;
 }
 
 void readParameter(Parameter* param, const char* filename)
@@ -79,6 +81,8 @@ void readParameter(Parameter* param, const char* filename)
             PARSE_REAL(u_init);
             PARSE_REAL(v_init);
             PARSE_REAL(p_init);
+            PARSE_INT(presmooth);
+            PARSE_INT(postsmooth);
         }
     }
 

BasicSolver/2D-seq/src/parameter.h

@@ -18,6 +18,7 @@ typedef struct {
     char* name;
     int bcLeft, bcRight, bcBottom, bcTop;
     double u_init, v_init, p_init;
+    int presmooth, postsmooth;
 } Parameter;
 
 void initParameter(Parameter*);

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

@@ -84,6 +84,37 @@ static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int
 
     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) {
+
+                P(i, j) -= factor * (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));
+
+            }
+            isw = 3 - isw;
+        }
+        jsw = 3 - jsw;
+    }
+
+}
+
+static double calculateResidual(Solver* s, double* p, double* rhs, int level, int imax, int jmax)
+{
+    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* r     = s->r[level];
+    double res    = 1.0;
+    int pass, jsw, isw;
+
+    jsw = 1;
+
     for (pass = 0; pass < 2; pass++) {
         isw = jsw;
 
@@ -94,7 +125,6 @@ static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int
                           ((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(i, j));
                 res += (R(i, j) * R(i, j));
             }
             isw = 3 - isw;
@@ -110,7 +140,7 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax,
 {
     double res = 0.0;
 
-    // coarsest level TODO: Use direct solver?
+    // coarsest level
     if (level == COARSEST_LEVEL) {
         for (int i = 0; i < 5; i++) {
             smooth(s, p, rhs, level, imax, jmax);
@@ -118,17 +148,18 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax,
         return res;
     }
 
-    // pre-smoothing TODO: Make smoothing steps configurable?
+    // pre-smoothing
-    for (int i = 0; i < 5; i++) {
+    for (int i = 0; i < s->presmooth; i++) {
         smooth(s, p, rhs, level, imax, jmax);
         if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
     }
 
+    res = calculateResidual(s, p, rhs, level, imax, jmax);
+
     // restrict
     restrictMG(s, level, imax, jmax);
 
     // MGSolver on residual and error.
-    // TODO: What if there is a rest?
     multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, imax / 2, jmax / 2);
 
     // prolongate
@@ -139,8 +170,8 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax,
     if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
 
     // post-smoothing
-    for (int i = 0; i < 5; i++) {
+    for (int i = 0; i < s->postsmooth; i++) {
-        res = smooth(s, p, rhs, level, imax, jmax);
+        smooth(s, p, rhs, level, imax, jmax);
         if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
     }
 
@@ -154,6 +185,8 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     s->itermax = p->itermax;
     s->levels  = p->levels;
     s->grid    = &d->grid;
+    s->presmooth = p->presmooth;
+    s->postsmooth = p->postsmooth;
 
     int imax   = s->grid->imax;
     int jmax   = s->grid->jmax;

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

@@ -18,6 +18,7 @@ typedef struct {
     int itermax;
     int levels;
     double **r, **e;
+        int presmooth, postsmooth;
 } Solver;
 
 extern void initSolver(Solver*, Discretization*, Parameter*);