Separate solver and discretization modules

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

@@ -0,0 +1,453 @@
+/*
+ * Copyright (C)  NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "allocate.h"
+#include "comm.h"
+#include "discretization.h"
+#include "parameter.h"
+#include "util.h"
+
+static void printConfig(Discretization* s)
+{
+    if (commIsMaster(&s->comm)) {
+        printf("Parameters for #%s#\n", s->problem);
+        printf("BC Left:%d Right:%d Bottom:%d Top:%d\n",
+            s->bcLeft,
+            s->bcRight,
+            s->bcBottom,
+            s->bcTop);
+        printf("\tReynolds number: %.2f\n", s->re);
+        printf("\tGx Gy: %.2f %.2f\n", s->gx, s->gy);
+        printf("Geometry data:\n");
+        printf("\tDomain box size (x, y): %.2f, %.2f\n",
+            s->grid.xlength,
+            s->grid.ylength);
+        printf("\tCells (x, y): %d, %d\n", s->grid.imax, s->grid.jmax);
+        printf("\tCell size (dx, dy): %f, %f\n", s->grid.dx, s->grid.dy);
+        printf("Timestep parameters:\n");
+        printf("\tDefault stepsize: %.2f, Final time %.2f\n", s->dt, s->te);
+        printf("\tdt bound: %.6f\n", s->dtBound);
+        printf("\tTau factor: %.2f\n", s->tau);
+        printf("Iterative s parameters:\n");
+        printf("\tgamma factor: %f\n", s->gamma);
+    }
+    commPrintConfig(&s->comm);
+}
+
+void initDiscretiztion(Discretization* s, Parameter* params)
+{
+    s->problem      = params->name;
+    s->bcLeft       = params->bcLeft;
+    s->bcRight      = params->bcRight;
+    s->bcBottom     = params->bcBottom;
+    s->bcTop        = params->bcTop;
+    s->grid.imax    = params->imax;
+    s->grid.jmax    = params->jmax;
+    s->grid.xlength = params->xlength;
+    s->grid.ylength = params->ylength;
+    s->grid.dx      = params->xlength / params->imax;
+    s->grid.dy      = params->ylength / params->jmax;
+    s->re           = params->re;
+    s->gx           = params->gx;
+    s->gy           = params->gy;
+    s->dt           = params->dt;
+    s->te           = params->te;
+    s->tau          = params->tau;
+    s->gamma        = params->gamma;
+
+    /* allocate arrays */
+    int imaxLocal = s->comm.imaxLocal;
+    int jmaxLocal = s->comm.jmaxLocal;
+    size_t size   = (imaxLocal + 2) * (jmaxLocal + 2);
+
+    s->u   = allocate(64, size * sizeof(double));
+    s->v   = allocate(64, size * sizeof(double));
+    s->p   = allocate(64, size * sizeof(double));
+    s->rhs = allocate(64, size * sizeof(double));
+    s->f   = allocate(64, size * sizeof(double));
+    s->g   = allocate(64, size * sizeof(double));
+
+    for (int i = 0; i < size; i++) {
+        s->u[i]   = params->u_init;
+        s->v[i]   = params->v_init;
+        s->p[i]   = params->p_init;
+        s->rhs[i] = 0.0;
+        s->f[i]   = 0.0;
+        s->g[i]   = 0.0;
+    }
+
+    double dx = s->grid.dx;
+    double dy = s->grid.dy;
+
+    double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
+    s->dtBound       = 0.5 * s->re * 1.0 / invSqrSum;
+#ifdef VERBOSE
+    printConfig(s);
+#endif
+}
+
+void computeRHS(Discretization* s)
+{
+    int imaxLocal = s->comm.imaxLocal;
+    int jmaxLocal = s->comm.jmaxLocal;
+    double idx    = 1.0 / s->grid.dx;
+    double idy    = 1.0 / s->grid.dy;
+    double idt    = 1.0 / s->dt;
+    double* rhs   = s->rhs;
+    double* f     = s->f;
+    double* g     = s->g;
+
+    commShift(&s->comm, f, g);
+
+    for (int j = 1; j < jmaxLocal + 1; j++) {
+        for (int i = 1; i < imaxLocal + 1; i++) {
+            RHS(i, j) = ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy) *
+                        idt;
+        }
+    }
+}
+
+static double maxElement(Discretization* s, double* m)
+{
+    int imaxLocal = s->comm.imaxLocal;
+    int jmaxLocal = s->comm.jmaxLocal;
+    int size      = (imaxLocal + 2) * (jmaxLocal + 2);
+    double maxval = DBL_MIN;
+
+    for (int i = 0; i < size; i++) {
+        maxval = MAX(maxval, fabs(m[i]));
+    }
+
+    commReduction(&maxval, MAX);
+    return maxval;
+}
+
+void computeTimestep(Discretization* s)
+{
+    double dt   = s->dtBound;
+    double dx   = s->grid.dx;
+    double dy   = s->grid.dy;
+    double umax = maxElement(s, s->u);
+    double vmax = maxElement(s, s->v);
+
+    if (umax > 0) {
+        dt = (dt > dx / umax) ? dx / umax : dt;
+    }
+    if (vmax > 0) {
+        dt = (dt > dy / vmax) ? dy / vmax : dt;
+    }
+
+    s->dt = dt * s->tau;
+}
+
+void setBoundaryConditions(Discretization* s)
+{
+    int imaxLocal = s->comm.imaxLocal;
+    int jmaxLocal = s->comm.jmaxLocal;
+    double* u     = s->u;
+    double* v     = s->v;
+
+    if (commIsBoundary(&s->comm, TOP)) {
+        switch (s->bcTop) {
+        case NOSLIP:
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                V(i, jmaxLocal)     = 0.0;
+                U(i, jmaxLocal + 1) = -U(i, jmaxLocal);
+            }
+            break;
+        case SLIP:
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                V(i, jmaxLocal)     = 0.0;
+                U(i, jmaxLocal + 1) = U(i, jmaxLocal);
+            }
+            break;
+        case OUTFLOW:
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                U(i, jmaxLocal + 1) = U(i, jmaxLocal);
+                V(i, jmaxLocal)     = V(i, jmaxLocal - 1);
+            }
+            break;
+        case PERIODIC:
+            break;
+        }
+    }
+
+    if (commIsBoundary(&s->comm, BOTTOM)) {
+        switch (s->bcBottom) {
+        case NOSLIP:
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                V(i, 0) = 0.0;
+                U(i, 0) = -U(i, 1);
+            }
+            break;
+        case SLIP:
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                V(i, 0) = 0.0;
+                U(i, 0) = U(i, 1);
+            }
+            break;
+        case OUTFLOW:
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                U(i, 0) = U(i, 1);
+                V(i, 0) = V(i, 1);
+            }
+            break;
+        case PERIODIC:
+            break;
+        }
+    }
+
+    if (commIsBoundary(&s->comm, RIGHT)) {
+        switch (s->bcRight) {
+        case NOSLIP:
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                U(imaxLocal, j)     = 0.0;
+                V(imaxLocal + 1, j) = -V(imaxLocal, j);
+            }
+            break;
+        case SLIP:
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                U(imaxLocal, j)     = 0.0;
+                V(imaxLocal + 1, j) = V(imaxLocal, j);
+            }
+            break;
+        case OUTFLOW:
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                U(imaxLocal, j)     = U(imaxLocal - 1, j);
+                V(imaxLocal + 1, j) = V(imaxLocal, j);
+            }
+            break;
+        case PERIODIC:
+            break;
+        }
+    }
+
+    if (commIsBoundary(&s->comm, LEFT)) {
+        switch (s->bcLeft) {
+        case NOSLIP:
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                U(0, j) = 0.0;
+                V(0, j) = -V(1, j);
+            }
+            break;
+        case SLIP:
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                U(0, j) = 0.0;
+                V(0, j) = V(1, j);
+            }
+            break;
+        case OUTFLOW:
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                U(0, j) = U(1, j);
+                V(0, j) = V(1, j);
+            }
+            break;
+        case PERIODIC:
+            break;
+        }
+    }
+}
+
+void setSpecialBoundaryCondition(Discretization* s)
+{
+    int imaxLocal = s->comm.imaxLocal;
+    int jmaxLocal = s->comm.jmaxLocal;
+    double* u     = s->u;
+
+    if (strcmp(s->problem, "dcavity") == 0) {
+        if (commIsBoundary(&s->comm, TOP)) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                U(i, jmaxLocal + 1) = 2.0 - U(i, jmaxLocal);
+            }
+        }
+    } else if (strcmp(s->problem, "canal") == 0) {
+        if (commIsBoundary(&s->comm, LEFT)) {
+            double ylength = s->grid.ylength;
+            double dy      = s->grid.dy;
+            int rest       = s->grid.jmax % s->comm.size;
+            int yc         = s->comm.rank * (s->grid.jmax / s->comm.size) +
+                     MIN(rest, s->comm.rank);
+            double ys = dy * (yc + 0.5);
+            double y;
+
+            /* printf("RANK %d yc: %d ys: %f\n", solver->rank, yc, ys); */
+
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                y       = ys + dy * (j - 0.5);
+                U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
+            }
+        }
+    }
+    /* print(solver, solver->u); */
+}
+
+void computeFG(Discretization* s)
+{
+    double* u = s->u;
+    double* v = s->v;
+    double* f = s->f;
+    double* g = s->g;
+
+    int imaxLocal = s->comm.imaxLocal;
+    int jmaxLocal = s->comm.jmaxLocal;
+
+    double gx        = s->gx;
+    double gy        = s->gy;
+    double gamma     = s->gamma;
+    double dt        = s->dt;
+    double inverseRe = 1.0 / s->re;
+    double inverseDx = 1.0 / s->grid.dx;
+    double inverseDy = 1.0 / s->grid.dy;
+    double du2dx, dv2dy, duvdx, duvdy;
+    double du2dx2, du2dy2, dv2dx2, dv2dy2;
+
+    commExchange(&s->comm, u);
+    commExchange(&s->comm, v);
+
+    for (int j = 1; j < jmaxLocal + 1; j++) {
+        for (int i = 1; i < imaxLocal + 1; i++) {
+            du2dx = inverseDx * 0.25 *
+                        ((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) -
+                            (U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) +
+                    gamma * inverseDx * 0.25 *
+                        (fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) +
+                            fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j)));
+
+            duvdy = inverseDy * 0.25 *
+                        ((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) -
+                            (V(i, j - 1) + V(i + 1, j - 1)) * (U(i, j) + U(i, j - 1))) +
+                    gamma * inverseDy * 0.25 *
+                        (fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) +
+                            fabs(V(i, j - 1) + V(i + 1, j - 1)) *
+                                (U(i, j) - U(i, j - 1)));
+
+            du2dx2  = inverseDx * inverseDx * (U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j));
+            du2dy2  = inverseDy * inverseDy * (U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1));
+            F(i, j) = U(i, j) + dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx);
+
+            duvdx = inverseDx * 0.25 *
+                        ((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) -
+                            (U(i - 1, j) + U(i - 1, j + 1)) * (V(i, j) + V(i - 1, j))) +
+                    gamma * inverseDx * 0.25 *
+                        (fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) +
+                            fabs(U(i - 1, j) + U(i - 1, j + 1)) *
+                                (V(i, j) - V(i - 1, j)));
+
+            dv2dy = inverseDy * 0.25 *
+                        ((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) -
+                            (V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) +
+                    gamma * inverseDy * 0.25 *
+                        (fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) +
+                            fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1)));
+
+            dv2dx2  = inverseDx * inverseDx * (V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j));
+            dv2dy2  = inverseDy * inverseDy * (V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1));
+            G(i, j) = V(i, j) + dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy);
+        }
+    }
+
+    /* ----------------------------- boundary of F --------------------------- */
+    if (commIsBoundary(&s->comm, LEFT)) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            F(0, j) = U(0, j);
+        }
+    }
+
+    if (commIsBoundary(&s->comm, RIGHT)) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            F(imaxLocal, j) = U(imaxLocal, j);
+        }
+    }
+
+    /* ----------------------------- boundary of G --------------------------- */
+    if (commIsBoundary(&s->comm, BOTTOM)) {
+        for (int i = 1; i < imaxLocal + 1; i++) {
+            G(i, 0) = V(i, 0);
+        }
+    }
+
+    if (commIsBoundary(&s->comm, TOP)) {
+        for (int i = 1; i < imaxLocal + 1; i++) {
+            G(i, jmaxLocal) = V(i, jmaxLocal);
+        }
+    }
+}
+
+void adaptUV(Discretization* s)
+{
+    int imaxLocal = s->comm.imaxLocal;
+    int jmaxLocal = s->comm.jmaxLocal;
+
+    double* p = s->p;
+    double* u = s->u;
+    double* v = s->v;
+    double* f = s->f;
+    double* g = s->g;
+
+    double factorX = s->dt / s->grid.dx;
+    double factorY = s->dt / s->grid.dy;
+
+    for (int j = 1; j < jmaxLocal + 1; j++) {
+        for (int i = 1; i < imaxLocal + 1; i++) {
+            U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX;
+            V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY;
+        }
+    }
+}
+
+void writeResult(Discretization* s, double* u, double* v, double* p)
+{
+    int imax  = s->grid.imax;
+    int jmax  = s->grid.jmax;
+    double dx = s->grid.dx;
+    double dy = s->grid.dy;
+    double x = 0.0, y = 0.0;
+
+    FILE* fp;
+    fp = fopen("pressure.dat", "w");
+
+    if (fp == NULL) {
+        printf("Error!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    for (int j = 1; j <= jmax; j++) {
+        y = (double)(j - 0.5) * dy;
+        for (int i = 1; i <= imax; i++) {
+            x = (double)(i - 0.5) * dx;
+            fprintf(fp, "%.2f %.2f %f\n", x, y, p[j * (imax + 2) + i]);
+        }
+        fprintf(fp, "\n");
+    }
+
+    fclose(fp);
+
+    fp = fopen("velocity.dat", "w");
+
+    if (fp == NULL) {
+        printf("Error!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    for (int j = 1; j <= jmax; j++) {
+        y = dy * (j - 0.5);
+        for (int i = 1; i <= imax; i++) {
+            x           = dx * (i - 0.5);
+            double velU = (u[j * (imax + 2) + i] + u[j * (imax + 2) + (i - 1)]) / 2.0;
+            double velV = (v[j * (imax + 2) + i] + v[(j - 1) * (imax + 2) + i]) / 2.0;
+            double len  = sqrt((velU * velU) + (velV * velV));
+            fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, velU, velV, len);
+        }
+    }
+
+    fclose(fp);
+}

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

@@ -0,0 +1,43 @@
+/*
+ * Copyright (C)  NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __DISCRETIZATION_H_
+#define __DISCRETIZATION_H_
+#include "comm.h"
+#include "grid.h"
+#include "parameter.h"
+
+enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
+
+typedef struct {
+    /* geometry and grid information */
+    Grid grid;
+    /* arrays */
+    double *p, *rhs;
+    double *f, *g;
+    double *u, *v;
+    /* parameters */
+    double re, tau, gamma;
+    double gx, gy;
+    /* time stepping */
+    double dt, te;
+    double dtBound;
+    char* problem;
+    int bcLeft, bcRight, bcBottom, bcTop;
+    /* communication */
+    Comm comm;
+} Discretization;
+
+void initDiscretiztion(Discretization*, Parameter*);
+void computeRHS(Discretization*);
+void normalizePressure(Discretization*);
+void computeTimestep(Discretization*);
+void setBoundaryConditions(Discretization*);
+void setSpecialBoundaryCondition(Discretization*);
+void computeFG(Discretization*);
+void adaptUV(Discretization*);
+void writeResult(Discretization* s, double* u, double* v, double* p);
+#endif

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

@@ -0,0 +1,16 @@
+/*
+ * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __GRID_H_
+#define __GRID_H_
+
+typedef struct {
+    double dx, dy;
+    int imax, jmax;
+    double xlength, ylength;
+} Grid;
+
+#endif // __GRID_H_

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

@@ -10,21 +10,22 @@
 
 #include "allocate.h"
 #include "comm.h"
+#include "discretization.h"
 #include "parameter.h"
 #include "progress.h"
 #include "solver.h"
 #include "timing.h"
 
-static void writeResults(Solver* s)
+static void writeResults(Discretization* s)
 {
 #ifdef _MPI
-    size_t bytesize = (s->imax + 2) * (s->jmax + 2) * sizeof(double);
+    size_t bytesize = (s->grid.imax + 2) * (s->grid.jmax + 2) * sizeof(double);
 
     double* ug = allocate(64, bytesize);
     double* vg = allocate(64, bytesize);
     double* pg = allocate(64, bytesize);
 
-    commCollectResult(&s->comm, ug, vg, pg, s->u, s->v, s->p, s->imax, s->jmax);
+    commCollectResult(&s->comm, ug, vg, pg, s->u, s->v, s->p, s->grid.imax, s->grid.jmax);
     writeResult(s, ug, vg, pg);
 
     free(ug);
@@ -40,9 +41,10 @@ int main(int argc, char** argv)
     int rank;
     double timeStart, timeStop;
     Parameter p;
+    Discretization d;
     Solver s;
 
-    commInit(&s.comm, argc, argv);
+    commInit(&d.comm, argc, argv);
     initParameter(&p);
 
     if (argc != 2) {
@@ -51,44 +53,45 @@ int main(int argc, char** argv)
     }
 
     readParameter(&p, argv[1]);
-    commPartition(&s.comm, p.jmax, p.imax);
+    commPartition(&d.comm, p.jmax, p.imax);
-    if (commIsMaster(&s.comm)) {
+    if (commIsMaster(&d.comm)) {
         printParameter(&p);
     }
 
-    initSolver(&s, &p);
+    initDiscretiztion(&d, &p);
+    initSolver(&s, &d, &p);
 #ifdef TEST
-    commPrintConfig(&s.comm);
+    commPrintConfig(&d.comm);
-    commTestInit(&s.comm, s.p, s.f, s.g);
+    commTestInit(&d.comm, d.p, d.f, d.g);
-    commExchange(&s.comm, s.p);
+    commExchange(&d.comm, d.p);
-    commShift(&s.comm, s.f, s.g);
+    commShift(&d.comm, d.f, d.g);
-    commTestWrite(&s.comm, s.p, s.f, s.g);
+    commTestWrite(&d.comm, d.p, d.f, d.g);
-    writeResults(&s);
+    writeResults(&d);
-    commFinalize(&s.comm);
+    commFinalize(&d.comm);
     exit(EXIT_SUCCESS);
 #endif
 #ifndef VERBOSE
-    initProgress(s.te);
+    initProgress(d.te);
 #endif
 
-    double tau = s.tau;
+    double tau = d.tau;
-    double te  = s.te;
+    double te  = d.te;
     double t   = 0.0;
 
     timeStart = getTimeStamp();
     while (t <= te) {
-        if (tau > 0.0) computeTimestep(&s);
+        if (tau > 0.0) computeTimestep(&d);
-        setBoundaryConditions(&s);
+        setBoundaryConditions(&d);
-        setSpecialBoundaryCondition(&s);
+        setSpecialBoundaryCondition(&d);
-        computeFG(&s);
+        computeFG(&d);
-        computeRHS(&s);
+        computeRHS(&d);
-        solve(&s);
+        solve(&s, d.p, d.rhs);
-        adaptUV(&s);
+        adaptUV(&d);
-        t += s.dt;
+        t += d.dt;
 
 #ifdef VERBOSE
-        if (commIsMaster(&s.comm)) {
+        if (commIsMaster(s.comm)) {
-            printf("TIME %f , TIMESTEP %f\n", t, s.dt);
+            printf("TIME %f , TIMESTEP %f\n", t, d.dt);
         }
 #else
         printProgress(t);
@@ -98,11 +101,11 @@ int main(int argc, char** argv)
 #ifndef VERBOSE
     stopProgress();
 #endif
-    if (commIsMaster(&s.comm)) {
+    if (commIsMaster(s.comm)) {
         printf("Solution took %.2fs\n", timeStop - timeStart);
     }
 
-    writeResults(&s);
+    writeResults(&d);
-    commFinalize(&s.comm);
+    commFinalize(s.comm);
     return EXIT_SUCCESS;
 }

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

@@ -4,144 +4,38 @@
  * Use of this source code is governed by a MIT style
  * license that can be found in the LICENSE file.
  */
-#include <float.h>
 #include <math.h>
 #include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
 
 #include "allocate.h"
 #include "comm.h"
+#include "discretization.h"
 #include "parameter.h"
 #include "solver.h"
 #include "util.h"
 
-#define P(i, j)   p[(j) * (imaxLocal + 2) + (i)]
+void initSolver(Solver* s, Discretization* d, Parameter* p)
-#define F(i, j)   f[(j) * (imaxLocal + 2) + (i)]
-#define G(i, j)   g[(j) * (imaxLocal + 2) + (i)]
-#define U(i, j)   u[(j) * (imaxLocal + 2) + (i)]
-#define V(i, j)   v[(j) * (imaxLocal + 2) + (i)]
-#define RHS(i, j) rhs[(j) * (imaxLocal + 2) + (i)]
-
-static void printConfig(Solver* s)
 {
-    if (commIsMaster(&s->comm)) {
+    s->grid    = &d->grid;
-        printf("Parameters for #%s#\n", s->problem);
+    s->eps     = p->eps;
-        printf("BC Left:%d Right:%d Bottom:%d Top:%d\n",
+    s->omega   = p->omg;
-            s->bcLeft,
+    s->itermax = p->itermax;
-            s->bcRight,
+    s->comm    = &d->comm;
-            s->bcBottom,
-            s->bcTop);
-        printf("\tReynolds number: %.2f\n", s->re);
-        printf("\tGx Gy: %.2f %.2f\n", s->gx, s->gy);
-        printf("Geometry data:\n");
-        printf("\tDomain box size (x, y): %.2f, %.2f\n", s->xlength, s->ylength);
-        printf("\tCells (x, y): %d, %d\n", s->imax, s->jmax);
-        printf("\tCell size (dx, dy): %f, %f\n", s->dx, s->dy);
-        printf("Timestep parameters:\n");
-        printf("\tDefault stepsize: %.2f, Final time %.2f\n", s->dt, s->te);
-        printf("\tdt bound: %.6f\n", s->dtBound);
-        printf("\tTau factor: %.2f\n", s->tau);
-        printf("Iterative s parameters:\n");
-        printf("\tMax iterations: %d\n", s->itermax);
-        printf("\tepsilon (stopping tolerance) : %f\n", s->eps);
-        printf("\tgamma factor: %f\n", s->gamma);
-        printf("\tomega (SOR relaxation): %f\n", s->omega);
-    }
-    commPrintConfig(&s->comm);
 }
 
-void initSolver(Solver* s, Parameter* params)
+void solve(Solver* s, double* p, double* rhs)
 {
-    s->problem  = params->name;
+    int imax      = s->grid->imax;
-    s->bcLeft   = params->bcLeft;
+    int jmax      = s->grid->jmax;
-    s->bcRight  = params->bcRight;
+    int imaxLocal = s->comm->imaxLocal;
-    s->bcBottom = params->bcBottom;
+    int jmaxLocal = s->comm->jmaxLocal;
-    s->bcTop    = params->bcTop;
-    s->imax     = params->imax;
-    s->jmax     = params->jmax;
-    s->xlength  = params->xlength;
-    s->ylength  = params->ylength;
-    s->dx       = params->xlength / params->imax;
-    s->dy       = params->ylength / params->jmax;
-    s->eps      = params->eps;
-    s->omega    = params->omg;
-    s->itermax  = params->itermax;
-    s->re       = params->re;
-    s->gx       = params->gx;
-    s->gy       = params->gy;
-    s->dt       = params->dt;
-    s->te       = params->te;
-    s->tau      = params->tau;
-    s->gamma    = params->gamma;
-
-    /* allocate arrays */
-    int imaxLocal = s->comm.imaxLocal;
-    int jmaxLocal = s->comm.jmaxLocal;
-    size_t size   = (imaxLocal + 2) * (jmaxLocal + 2);
-
-    s->u   = allocate(64, size * sizeof(double));
-    s->v   = allocate(64, size * sizeof(double));
-    s->p   = allocate(64, size * sizeof(double));
-    s->rhs = allocate(64, size * sizeof(double));
-    s->f   = allocate(64, size * sizeof(double));
-    s->g   = allocate(64, size * sizeof(double));
-
-    for (int i = 0; i < size; i++) {
-        s->u[i]   = params->u_init;
-        s->v[i]   = params->v_init;
-        s->p[i]   = params->p_init;
-        s->rhs[i] = 0.0;
-        s->f[i]   = 0.0;
-        s->g[i]   = 0.0;
-    }
-
-    double dx = s->dx;
-    double dy = s->dy;
-
-    double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
-    s->dtBound       = 0.5 * s->re * 1.0 / invSqrSum;
-#ifdef VERBOSE
-    printConfig(s);
-#endif
-}
-
-void computeRHS(Solver* s)
-{
-    int imaxLocal = s->comm.imaxLocal;
-    int jmaxLocal = s->comm.jmaxLocal;
-    double idx    = 1.0 / s->dx;
-    double idy    = 1.0 / s->dy;
-    double idt    = 1.0 / s->dt;
-    double* rhs   = s->rhs;
-    double* f     = s->f;
-    double* g     = s->g;
-
-    commShift(&s->comm, f, g);
-
-    for (int j = 1; j < jmaxLocal + 1; j++) {
-        for (int i = 1; i < imaxLocal + 1; i++) {
-            RHS(i, j) = ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy) *
-                        idt;
-        }
-    }
-}
-
-void solve(Solver* s)
-{
-    int imax      = s->imax;
-    int jmax      = s->jmax;
-    int imaxLocal = s->comm.imaxLocal;
-    int jmaxLocal = s->comm.jmaxLocal;
     double eps    = s->eps;
     int itermax   = s->itermax;
-    double dx2    = s->dx * s->dx;
+    double dx2    = s->grid->dx * s->grid->dx;
-    double dy2    = s->dy * s->dy;
+    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* p     = s->p;
-    double* rhs   = s->rhs;
     double epssq  = eps * eps;
     int pass, jsw, isw;
     int it     = 0;
@@ -151,7 +45,7 @@ void solve(Solver* s)
         jsw = 1;
         for (pass = 0; pass < 2; pass++) {
             isw = jsw;
-            commExchange(&s->comm, p);
+            commExchange(s->comm, p);
 
             for (int j = 1; j < jmaxLocal + 1; j++) {
                 for (int i = isw; i < imaxLocal + 1; i += 2) {
@@ -168,25 +62,25 @@ void solve(Solver* s)
             jsw = 3 - jsw;
         }
 
-        if (commIsBoundary(&s->comm, BOTTOM)) { // set bottom bc
+        if (commIsBoundary(s->comm, BOTTOM)) { // set bottom bc
             for (int i = 1; i < imaxLocal + 1; i++) {
                 P(i, 0) = P(i, 1);
             }
         }
 
-        if (commIsBoundary(&s->comm, TOP)) { // set top bc
+        if (commIsBoundary(s->comm, TOP)) { // set top bc
             for (int i = 1; i < imaxLocal + 1; i++) {
                 P(i, jmaxLocal + 1) = P(i, jmaxLocal);
             }
         }
 
-        if (commIsBoundary(&s->comm, LEFT)) { // set left bc
+        if (commIsBoundary(s->comm, LEFT)) { // set left bc
             for (int j = 1; j < jmaxLocal + 1; j++) {
                 P(0, j) = P(1, j);
             }
         }
 
-        if (commIsBoundary(&s->comm, RIGHT)) { // set right bc
+        if (commIsBoundary(s->comm, RIGHT)) { // set right bc
             for (int j = 1; j < jmaxLocal + 1; j++) {
                 P(imaxLocal + 1, j) = P(imaxLocal, j);
             }
@@ -195,7 +89,7 @@ void solve(Solver* s)
         commReduction(&res, SUM);
         res = res / (double)(imax * jmax);
 #ifdef DEBUG
-        if (commIsMaster(&s->comm)) {
+        if (commIsMaster(s->comm)) {
             printf("%d Residuum: %e\n", it, res);
         }
 #endif
@@ -203,343 +97,8 @@ void solve(Solver* s)
     }
 
 #ifdef VERBOSE
-    if (commIsMaster(&s->comm)) {
+    if (commIsMaster(s->comm)) {
         printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
     }
 #endif
 }
-
-static double maxElement(Solver* s, double* m)
-{
-    int imaxLocal = s->comm.imaxLocal;
-    int jmaxLocal = s->comm.jmaxLocal;
-    int size      = (imaxLocal + 2) * (jmaxLocal + 2);
-    double maxval = DBL_MIN;
-
-    for (int i = 0; i < size; i++) {
-        maxval = MAX(maxval, fabs(m[i]));
-    }
-
-    commReduction(&maxval, MAX);
-    return maxval;
-}
-
-void computeTimestep(Solver* s)
-{
-    double dt   = s->dtBound;
-    double dx   = s->dx;
-    double dy   = s->dy;
-    double umax = maxElement(s, s->u);
-    double vmax = maxElement(s, s->v);
-
-    if (umax > 0) {
-        dt = (dt > dx / umax) ? dx / umax : dt;
-    }
-    if (vmax > 0) {
-        dt = (dt > dy / vmax) ? dy / vmax : dt;
-    }
-
-    s->dt = dt * s->tau;
-}
-
-void setBoundaryConditions(Solver* s)
-{
-    int imaxLocal = s->comm.imaxLocal;
-    int jmaxLocal = s->comm.jmaxLocal;
-    double* u     = s->u;
-    double* v     = s->v;
-
-    if (commIsBoundary(&s->comm, TOP)) {
-        switch (s->bcTop) {
-        case NOSLIP:
-            for (int i = 1; i < imaxLocal + 1; i++) {
-                V(i, jmaxLocal)     = 0.0;
-                U(i, jmaxLocal + 1) = -U(i, jmaxLocal);
-            }
-            break;
-        case SLIP:
-            for (int i = 1; i < imaxLocal + 1; i++) {
-                V(i, jmaxLocal)     = 0.0;
-                U(i, jmaxLocal + 1) = U(i, jmaxLocal);
-            }
-            break;
-        case OUTFLOW:
-            for (int i = 1; i < imaxLocal + 1; i++) {
-                U(i, jmaxLocal + 1) = U(i, jmaxLocal);
-                V(i, jmaxLocal)     = V(i, jmaxLocal - 1);
-            }
-            break;
-        case PERIODIC:
-            break;
-        }
-    }
-
-    if (commIsBoundary(&s->comm, BOTTOM)) {
-        switch (s->bcBottom) {
-        case NOSLIP:
-            for (int i = 1; i < imaxLocal + 1; i++) {
-                V(i, 0) = 0.0;
-                U(i, 0) = -U(i, 1);
-            }
-            break;
-        case SLIP:
-            for (int i = 1; i < imaxLocal + 1; i++) {
-                V(i, 0) = 0.0;
-                U(i, 0) = U(i, 1);
-            }
-            break;
-        case OUTFLOW:
-            for (int i = 1; i < imaxLocal + 1; i++) {
-                U(i, 0) = U(i, 1);
-                V(i, 0) = V(i, 1);
-            }
-            break;
-        case PERIODIC:
-            break;
-        }
-    }
-
-    if (commIsBoundary(&s->comm, RIGHT)) {
-        switch (s->bcRight) {
-        case NOSLIP:
-            for (int j = 1; j < jmaxLocal + 1; j++) {
-                U(imaxLocal, j)     = 0.0;
-                V(imaxLocal + 1, j) = -V(imaxLocal, j);
-            }
-            break;
-        case SLIP:
-            for (int j = 1; j < jmaxLocal + 1; j++) {
-                U(imaxLocal, j)     = 0.0;
-                V(imaxLocal + 1, j) = V(imaxLocal, j);
-            }
-            break;
-        case OUTFLOW:
-            for (int j = 1; j < jmaxLocal + 1; j++) {
-                U(imaxLocal, j)     = U(imaxLocal - 1, j);
-                V(imaxLocal + 1, j) = V(imaxLocal, j);
-            }
-            break;
-        case PERIODIC:
-            break;
-        }
-    }
-
-    if (commIsBoundary(&s->comm, LEFT)) {
-        switch (s->bcLeft) {
-        case NOSLIP:
-            for (int j = 1; j < jmaxLocal + 1; j++) {
-                U(0, j) = 0.0;
-                V(0, j) = -V(1, j);
-            }
-            break;
-        case SLIP:
-            for (int j = 1; j < jmaxLocal + 1; j++) {
-                U(0, j) = 0.0;
-                V(0, j) = V(1, j);
-            }
-            break;
-        case OUTFLOW:
-            for (int j = 1; j < jmaxLocal + 1; j++) {
-                U(0, j) = U(1, j);
-                V(0, j) = V(1, j);
-            }
-            break;
-        case PERIODIC:
-            break;
-        }
-    }
-}
-
-void setSpecialBoundaryCondition(Solver* s)
-{
-    int imaxLocal = s->comm.imaxLocal;
-    int jmaxLocal = s->comm.jmaxLocal;
-    double* u     = s->u;
-
-    if (strcmp(s->problem, "dcavity") == 0) {
-        if (commIsBoundary(&s->comm, TOP)) {
-            for (int i = 1; i < imaxLocal + 1; i++) {
-                U(i, jmaxLocal + 1) = 2.0 - U(i, jmaxLocal);
-            }
-        }
-    } else if (strcmp(s->problem, "canal") == 0) {
-        if (commIsBoundary(&s->comm, LEFT)) {
-            double ylength = s->ylength;
-            double dy      = s->dy;
-            int rest       = s->jmax % s->comm.size;
-            int yc    = s->comm.rank * (s->jmax / s->comm.size) + MIN(rest, s->comm.rank);
-            double ys = dy * (yc + 0.5);
-            double y;
-
-            /* printf("RANK %d yc: %d ys: %f\n", solver->rank, yc, ys); */
-
-            for (int j = 1; j < jmaxLocal + 1; j++) {
-                y       = ys + dy * (j - 0.5);
-                U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
-            }
-        }
-    }
-    /* print(solver, solver->u); */
-}
-
-void computeFG(Solver* s)
-{
-    double* u = s->u;
-    double* v = s->v;
-    double* f = s->f;
-    double* g = s->g;
-
-    int imaxLocal = s->comm.imaxLocal;
-    int jmaxLocal = s->comm.jmaxLocal;
-
-    double gx        = s->gx;
-    double gy        = s->gy;
-    double gamma     = s->gamma;
-    double dt        = s->dt;
-    double inverseRe = 1.0 / s->re;
-    double inverseDx = 1.0 / s->dx;
-    double inverseDy = 1.0 / s->dy;
-    double du2dx, dv2dy, duvdx, duvdy;
-    double du2dx2, du2dy2, dv2dx2, dv2dy2;
-
-    commExchange(&s->comm, u);
-    commExchange(&s->comm, v);
-
-    for (int j = 1; j < jmaxLocal + 1; j++) {
-        for (int i = 1; i < imaxLocal + 1; i++) {
-            du2dx = inverseDx * 0.25 *
-                        ((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) -
-                            (U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) +
-                    gamma * inverseDx * 0.25 *
-                        (fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) +
-                            fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j)));
-
-            duvdy = inverseDy * 0.25 *
-                        ((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) -
-                            (V(i, j - 1) + V(i + 1, j - 1)) * (U(i, j) + U(i, j - 1))) +
-                    gamma * inverseDy * 0.25 *
-                        (fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) +
-                            fabs(V(i, j - 1) + V(i + 1, j - 1)) *
-                                (U(i, j) - U(i, j - 1)));
-
-            du2dx2  = inverseDx * inverseDx * (U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j));
-            du2dy2  = inverseDy * inverseDy * (U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1));
-            F(i, j) = U(i, j) + dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx);
-
-            duvdx = inverseDx * 0.25 *
-                        ((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) -
-                            (U(i - 1, j) + U(i - 1, j + 1)) * (V(i, j) + V(i - 1, j))) +
-                    gamma * inverseDx * 0.25 *
-                        (fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) +
-                            fabs(U(i - 1, j) + U(i - 1, j + 1)) *
-                                (V(i, j) - V(i - 1, j)));
-
-            dv2dy = inverseDy * 0.25 *
-                        ((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) -
-                            (V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) +
-                    gamma * inverseDy * 0.25 *
-                        (fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) +
-                            fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1)));
-
-            dv2dx2  = inverseDx * inverseDx * (V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j));
-            dv2dy2  = inverseDy * inverseDy * (V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1));
-            G(i, j) = V(i, j) + dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy);
-        }
-    }
-
-    /* ----------------------------- boundary of F --------------------------- */
-    if (commIsBoundary(&s->comm, LEFT)) {
-        for (int j = 1; j < jmaxLocal + 1; j++) {
-            F(0, j) = U(0, j);
-        }
-    }
-
-    if (commIsBoundary(&s->comm, RIGHT)) {
-        for (int j = 1; j < jmaxLocal + 1; j++) {
-            F(imaxLocal, j) = U(imaxLocal, j);
-        }
-    }
-
-    /* ----------------------------- boundary of G --------------------------- */
-    if (commIsBoundary(&s->comm, BOTTOM)) {
-        for (int i = 1; i < imaxLocal + 1; i++) {
-            G(i, 0) = V(i, 0);
-        }
-    }
-
-    if (commIsBoundary(&s->comm, TOP)) {
-        for (int i = 1; i < imaxLocal + 1; i++) {
-            G(i, jmaxLocal) = V(i, jmaxLocal);
-        }
-    }
-}
-
-void adaptUV(Solver* s)
-{
-    int imaxLocal = s->comm.imaxLocal;
-    int jmaxLocal = s->comm.jmaxLocal;
-
-    double* p = s->p;
-    double* u = s->u;
-    double* v = s->v;
-    double* f = s->f;
-    double* g = s->g;
-
-    double factorX = s->dt / s->dx;
-    double factorY = s->dt / s->dy;
-
-    for (int j = 1; j < jmaxLocal + 1; j++) {
-        for (int i = 1; i < imaxLocal + 1; i++) {
-            U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX;
-            V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY;
-        }
-    }
-}
-
-void writeResult(Solver* s, double* u, double* v, double* p)
-{
-    int imax  = s->imax;
-    int jmax  = s->jmax;
-    double dx = s->dx;
-    double dy = s->dy;
-    double x = 0.0, y = 0.0;
-
-    FILE* fp;
-    fp = fopen("pressure.dat", "w");
-
-    if (fp == NULL) {
-        printf("Error!\n");
-        exit(EXIT_FAILURE);
-    }
-
-    for (int j = 1; j <= jmax; j++) {
-        y = (double)(j - 0.5) * dy;
-        for (int i = 1; i <= imax; i++) {
-            x = (double)(i - 0.5) * dx;
-            fprintf(fp, "%.2f %.2f %f\n", x, y, p[j * (imax + 2) + i]);
-        }
-        fprintf(fp, "\n");
-    }
-
-    fclose(fp);
-
-    fp = fopen("velocity.dat", "w");
-
-    if (fp == NULL) {
-        printf("Error!\n");
-        exit(EXIT_FAILURE);
-    }
-
-    for (int j = 1; j <= jmax; j++) {
-        y = dy * (j - 0.5);
-        for (int i = 1; i <= imax; i++) {
-            x           = dx * (i - 0.5);
-            double velU = (u[j * (imax + 2) + i] + u[j * (imax + 2) + (i - 1)]) / 2.0;
-            double velV = (v[j * (imax + 2) + i] + v[(j - 1) * (imax + 2) + i]) / 2.0;
-            double len  = sqrt((velU * velU) + (velV * velV));
-            fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, velU, velV, len);
-        }
-    }
-
-    fclose(fp);
-}

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

@@ -7,41 +7,20 @@
 #ifndef __SOLVER_H_
 #define __SOLVER_H_
 #include "comm.h"
+#include "discretization.h"
+#include "grid.h"
 #include "parameter.h"
 
-enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
-
 typedef struct {
     /* geometry and grid information */
-    double dx, dy;
+    Grid* grid;
-    int imax, jmax;
-    double xlength, ylength;
-    /* arrays */
-    double *p, *rhs;
-    double *f, *g;
-    double *u, *v;
     /* parameters */
     double eps, omega;
-    double re, tau, gamma;
-    double gx, gy;
-    /* time stepping */
     int itermax;
-    double dt, te;
-    double dtBound;
-    char* problem;
-    int bcLeft, bcRight, bcBottom, bcTop;
     /* communication */
-    Comm comm;
+    Comm* comm;
 } Solver;
 
-void initSolver(Solver*, Parameter*);
+void initSolver(Solver*, Discretization*, Parameter*);
-void computeRHS(Solver*);
+void solve(Solver*, double*, double*);
-void solve(Solver*);
-void normalizePressure(Solver*);
-void computeTimestep(Solver*);
-void setBoundaryConditions(Solver*);
-void setSpecialBoundaryCondition(Solver*);
-void computeFG(Solver*);
-void adaptUV(Solver*);
-void writeResult(Solver* s, double* u, double* v, double* p);
 #endif

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

@@ -19,4 +19,11 @@
 #define ABS(a) ((a) >= 0 ? (a) : -(a))
 #endif
 
+#define P(i, j)   p[(j) * (imaxLocal + 2) + (i)]
+#define F(i, j)   f[(j) * (imaxLocal + 2) + (i)]
+#define G(i, j)   g[(j) * (imaxLocal + 2) + (i)]
+#define U(i, j)   u[(j) * (imaxLocal + 2) + (i)]
+#define V(i, j)   v[(j) * (imaxLocal + 2) + (i)]
+#define RHS(i, j) rhs[(j) * (imaxLocal + 2) + (i)]
+
 #endif // __UTIL_H_

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

@@ -0,0 +1,16 @@
+/*
+ * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __GRID_H_
+#define __GRID_H_
+
+typedef struct {
+    double dx, dy;
+    int imax, jmax;
+    double xlength, ylength;
+} Grid;
+
+#endif // __GRID_H_