Separate discretization and solver. Port Multigrid solver.

BasicSolver/2D-mpi/Makefile

@@ -1,5 +1,5 @@
 #=======================================================================================
-# Copyright (C)  NHR@FAU, University Erlangen-Nuremberg.
+# Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
 # All rights reserved.
 # Use of this source code is governed by a MIT-style
 # license that can be found in the LICENSE file.

BasicSolver/2D-seq/Makefile

@@ -1,5 +1,5 @@
 #=======================================================================================
-# Copyright (C)  NHR@FAU, University Erlangen-Nuremberg.
+# Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
 # All rights reserved.
 # Use of this source code is governed by a MIT-style
 # license that can be found in the LICENSE file.
@@ -21,6 +21,7 @@ VPATH     = $(SRC_DIR)
 SRC       = $(filter-out $(wildcard $(SRC_DIR)/*-*.c),$(wildcard $(SRC_DIR)/*.c))
 ASM       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s, $(SRC))
 OBJ       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC))
+OBJ      += $(BUILD_DIR)/solver-$(SOLVER).o
 SOURCES   = $(SRC) $(wildcard $(SRC_DIR)/*.h)
 CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
 

BasicSolver/2D-seq/config.mk

@@ -1,12 +1,12 @@
 # Supported: GCC, CLANG, ICC
 TAG ?= CLANG
 ENABLE_OPENMP ?= false
+# Supported: sor, mg
+SOLVER ?= sor
+# Run in debug settings
+DEBUG ?= false
 
 #Feature options
 OPTIONS +=  -DARRAY_ALIGNMENT=64
-# OPTIONS +=  -DVERBOSE
+#OPTIONS +=  -DVERBOSE
 #OPTIONS +=  -DDEBUG
-#OPTIONS +=  -DBOUNDCHECK
-#OPTIONS +=  -DVERBOSE_AFFINITY
-#OPTIONS +=  -DVERBOSE_DATASIZE
-#OPTIONS +=  -DVERBOSE_TIMER

BasicSolver/2D-seq/dcavity.par

@@ -36,11 +36,13 @@ te      10.0		# final time
 dt      0.02	    # time stepsize
 tau     0.5	    	# safety factor for time stepsize control (<0 constant delt)
 
-# Pressure Iteration Data:
+# Solver Data:
 # -----------------------
 
 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
 gamma    0.9		# upwind differencing factor gamma
+levels   5         # Multigrid levels
 #===============================================================================

BasicSolver/2D-seq/include_CLANG.mk

@@ -2,16 +2,18 @@ CC   = clang
 GCC  = cc
 LINKER = $(CC)
 
-ifeq ($(ENABLE_OPENMP),true)
+ifeq ($(strip $(ENABLE_OPENMP)),true)
 OPENMP   = -fopenmp
 #OPENMP   = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp
 LIBS     = # -lomp
 endif
+ifeq ($(strip $(DEBUG)),true)
+CFLAGS   = -O0 -g -std=c17
+else
+CFLAGS   = -O3 -std=c17 $(OPENMP)
+endif
 
 VERSION  = --version
-# CFLAGS   = -O3 -std=c17 $(OPENMP)
-CFLAGS   = -Ofast -std=c17
-#CFLAGS   = -Ofast -fnt-store=aggressive  -std=c99 $(OPENMP) #AMD CLANG
 LFLAGS   = $(OPENMP) -lm
-DEFINES  = -D_GNU_SOURCE# -DDEBUG
+DEFINES  = -D_GNU_SOURCE
 INCLUDES =

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

@@ -0,0 +1,435 @@
+/*
+ * 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 "discretization.h"
+#include "parameter.h"
+#include "util.h"
+
+static void print(Discretization* d, double* grid)
+{
+    int imax = d->grid.imax;
+
+    for (int j = 0; j < d->grid.jmax + 2; j++) {
+        printf("%02d: ", j);
+        for (int i = 0; i < d->grid.imax + 2; i++) {
+            printf("%12.8f  ", grid[j * (imax + 2) + i]);
+        }
+        printf("\n");
+    }
+    fflush(stdout);
+}
+
+static void printConfig(Discretization* d)
+{
+    printf("Parameters for #%s#\n", d->problem);
+    printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n",
+        d->bcLeft,
+        d->bcRight,
+        d->bcBottom,
+        d->bcTop);
+    printf("\tReynolds number: %.2f\n", d->re);
+    printf("\tGx Gy: %.2f %.2f\n", d->gx, d->gy);
+    printf("Geometry data:\n");
+    printf("\tDomain box size (x, y): %.2f, %.2f\n", d->grid.xlength, d->grid.ylength);
+    printf("\tCells (x, y): %d, %d\n", d->grid.imax, d->grid.jmax);
+    printf("Timestep parameters:\n");
+    printf("\tDefault stepsize: %.2f, Final time %.2f\n", d->dt, d->te);
+    printf("\tdt bound: %.6f\n", d->dtBound);
+    printf("\tTau factor: %.2f\n", d->tau);
+    printf("Iterative d parameters:\n");
+    printf("\tgamma factor: %f\n", d->gamma);
+}
+
+void initDiscretization(Discretization* d, Parameter* p)
+{
+    d->problem      = p->name;
+    d->bcLeft       = p->bcLeft;
+    d->bcRight      = p->bcRight;
+    d->bcBottom     = p->bcBottom;
+    d->bcTop        = p->bcTop;
+    d->grid.imax    = p->imax;
+    d->grid.jmax    = p->jmax;
+    d->grid.xlength = p->xlength;
+    d->grid.ylength = p->ylength;
+    d->grid.dx      = p->xlength / p->imax;
+    d->grid.dy      = p->ylength / p->jmax;
+    d->re           = p->re;
+    d->gx           = p->gx;
+    d->gy           = p->gy;
+    d->dt           = p->dt;
+    d->te           = p->te;
+    d->tau          = p->tau;
+    d->gamma        = p->gamma;
+
+    int imax    = d->grid.imax;
+    int jmax    = d->grid.jmax;
+    size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
+    d->u        = allocate(64, size);
+    d->v        = allocate(64, size);
+    d->p        = allocate(64, size);
+    d->rhs      = allocate(64, size);
+    d->f        = allocate(64, size);
+    d->g        = allocate(64, size);
+
+    for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
+        d->u[i]   = p->u_init;
+        d->v[i]   = p->v_init;
+        d->p[i]   = p->p_init;
+        d->rhs[i] = 0.0;
+        d->f[i]   = 0.0;
+        d->g[i]   = 0.0;
+    }
+
+    double dx        = d->grid.dx;
+    double dy        = d->grid.dy;
+    double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
+    d->dtBound       = 0.5 * d->re * 1.0 / invSqrSum;
+#ifdef VERBOSE
+    printConfig(d);
+#endif
+}
+
+void computeRHS(Discretization* d)
+{
+    int imax    = d->grid.imax;
+    int jmax    = d->grid.jmax;
+    double idx  = 1.0 / d->grid.dx;
+    double idy  = 1.0 / d->grid.dy;
+    double idt  = 1.0 / d->dt;
+    double* rhs = d->rhs;
+    double* f   = d->f;
+    double* g   = d->g;
+
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 1; i++) {
+            RHS(i, j) = idt *
+                        ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy);
+        }
+    }
+}
+
+static double maxElement(Discretization* d, double* m)
+{
+    int size      = (d->grid.imax + 2) * (d->grid.jmax + 2);
+    double maxval = DBL_MIN;
+
+    for (int i = 0; i < size; i++) {
+        maxval = MAX(maxval, fabs(m[i]));
+    }
+
+    return maxval;
+}
+
+void normalizePressure(Discretization* d)
+{
+    int size    = (d->grid.imax + 2) * (d->grid.jmax + 2);
+    double* p   = d->p;
+    double avgP = 0.0;
+
+    for (int i = 0; i < size; i++) {
+        avgP += p[i];
+    }
+    avgP /= size;
+
+    for (int i = 0; i < size; i++) {
+        p[i] = p[i] - avgP;
+    }
+}
+
+void computeTimestep(Discretization* d)
+{
+    double dt   = d->dtBound;
+    double dx   = d->grid.dx;
+    double dy   = d->grid.dy;
+    double umax = maxElement(d, d->u);
+    double vmax = maxElement(d, d->v);
+
+    if (umax > 0) {
+        dt = (dt > dx / umax) ? dx / umax : dt;
+    }
+    if (vmax > 0) {
+        dt = (dt > dy / vmax) ? dy / vmax : dt;
+    }
+
+    d->dt = dt * d->tau;
+}
+
+void setBoundaryConditions(Discretization* d)
+{
+    int imax  = d->grid.imax;
+    int jmax  = d->grid.jmax;
+    double* u = d->u;
+    double* v = d->v;
+
+    // Left boundary
+    switch (d->bcLeft) {
+    case NOSLIP:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(0, j) = 0.0;
+            V(0, j) = -V(1, j);
+        }
+        break;
+    case SLIP:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(0, j) = 0.0;
+            V(0, j) = V(1, j);
+        }
+        break;
+    case OUTFLOW:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(0, j) = U(1, j);
+            V(0, j) = V(1, j);
+        }
+        break;
+    case PERIODIC:
+        break;
+    }
+
+    // Right boundary
+    switch (d->bcRight) {
+    case NOSLIP:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(imax, j)     = 0.0;
+            V(imax + 1, j) = -V(imax, j);
+        }
+        break;
+    case SLIP:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(imax, j)     = 0.0;
+            V(imax + 1, j) = V(imax, j);
+        }
+        break;
+    case OUTFLOW:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(imax, j)     = U(imax - 1, j);
+            V(imax + 1, j) = V(imax, j);
+        }
+        break;
+    case PERIODIC:
+        break;
+    }
+
+    // Bottom boundary
+    switch (d->bcBottom) {
+    case NOSLIP:
+        for (int i = 1; i < imax + 1; i++) {
+            V(i, 0) = 0.0;
+            U(i, 0) = -U(i, 1);
+        }
+        break;
+    case SLIP:
+        for (int i = 1; i < imax + 1; i++) {
+            V(i, 0) = 0.0;
+            U(i, 0) = U(i, 1);
+        }
+        break;
+    case OUTFLOW:
+        for (int i = 1; i < imax + 1; i++) {
+            U(i, 0) = U(i, 1);
+            V(i, 0) = V(i, 1);
+        }
+        break;
+    case PERIODIC:
+        break;
+    }
+
+    // Top boundary
+    switch (d->bcTop) {
+    case NOSLIP:
+        for (int i = 1; i < imax + 1; i++) {
+            V(i, jmax)     = 0.0;
+            U(i, jmax + 1) = -U(i, jmax);
+        }
+        break;
+    case SLIP:
+        for (int i = 1; i < imax + 1; i++) {
+            V(i, jmax)     = 0.0;
+            U(i, jmax + 1) = U(i, jmax);
+        }
+        break;
+    case OUTFLOW:
+        for (int i = 1; i < imax + 1; i++) {
+            U(i, jmax + 1) = U(i, jmax);
+            V(i, jmax)     = V(i, jmax - 1);
+        }
+        break;
+    case PERIODIC:
+        break;
+    }
+}
+
+void setSpecialBoundaryCondition(Discretization* d)
+{
+    int imax   = d->grid.imax;
+    int jmax   = d->grid.jmax;
+    double mDy = d->grid.dy;
+    double* u  = d->u;
+
+    if (strcmp(d->problem, "dcavity") == 0) {
+        for (int i = 1; i < imax; i++) {
+            U(i, jmax + 1) = 2.0 - U(i, jmax);
+        }
+    } else if (strcmp(d->problem, "canal") == 0) {
+        double ylength = d->grid.ylength;
+        double y;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            y       = mDy * (j - 0.5);
+            U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
+        }
+    }
+}
+
+void computeFG(Discretization* d)
+{
+    double* u        = d->u;
+    double* v        = d->v;
+    double* f        = d->f;
+    double* g        = d->g;
+    int imax         = d->grid.imax;
+    int jmax         = d->grid.jmax;
+    double gx        = d->gx;
+    double gy        = d->gy;
+    double gamma     = d->gamma;
+    double dt        = d->dt;
+    double inverseRe = 1.0 / d->re;
+    double inverseDx = 1.0 / d->grid.dx;
+    double inverseDy = 1.0 / d->grid.dy;
+    double du2dx, dv2dy, duvdx, duvdy;
+    double du2dx2, du2dy2, dv2dx2, dv2dy2;
+
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 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 --------------------------- */
+    for (int j = 1; j < jmax + 1; j++) {
+        F(0, j)    = U(0, j);
+        F(imax, j) = U(imax, j);
+    }
+
+    /* ---------------------- boundary of G --------------------------- */
+    for (int i = 1; i < imax + 1; i++) {
+        G(i, 0)    = V(i, 0);
+        G(i, jmax) = V(i, jmax);
+    }
+}
+
+void adaptUV(Discretization* d)
+{
+    int imax       = d->grid.imax;
+    int jmax       = d->grid.jmax;
+    double* p      = d->p;
+    double* u      = d->u;
+    double* v      = d->v;
+    double* f      = d->f;
+    double* g      = d->g;
+    double factorX = d->dt / d->grid.dx;
+    double factorY = d->dt / d->grid.dy;
+
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 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* d)
+{
+    int imax  = d->grid.imax;
+    int jmax  = d->grid.jmax;
+    double dx = d->grid.dx;
+    double dy = d->grid.dy;
+    double* p = d->p;
+    double* u = d->u;
+    double* v = d->v;
+    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 + 1; j++) {
+        y = (double)(j - 0.5) * dy;
+        for (int i = 1; i < imax + 1; i++) {
+            x = (double)(i - 0.5) * dx;
+            fprintf(fp, "%.2f %.2f %f\n", x, y, P(i, j));
+        }
+        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 + 1; j++) {
+        y = dy * (j - 0.5);
+        for (int i = 1; i < imax + 1; i++) {
+            x           = dx * (i - 0.5);
+            double velU = (U(i, j) + U(i - 1, j)) / 2.0;
+            double velV = (V(i, j) + V(i, j - 1)) / 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-seq/src/discretization.h

@@ -0,0 +1,40 @@
+/*
+ * 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 "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;
+} Discretization;
+
+extern void initDiscretization(Discretization*, Parameter*);
+extern void computeRHS(Discretization*);
+extern void normalizePressure(Discretization*);
+extern void computeTimestep(Discretization*);
+extern void setBoundaryConditions(Discretization*);
+extern void setSpecialBoundaryCondition(Discretization*);
+extern void computeFG(Discretization*);
+extern void adaptUV(Discretization*);
+extern void writeResult(Discretization*);
+#endif

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

@@ -4,12 +4,11 @@
  * Use of this source code is governed by a MIT-style
  * license that can be found in the LICENSE file.
  */
-#include <float.h>
-#include <limits.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
 
+#include "discretization.h"
 #include "parameter.h"
 #include "progress.h"
 #include "solver.h"
@@ -17,39 +16,41 @@
 
 int main(int argc, char** argv)
 {
-    double S, E;
+    double timeStart, timeStop;
-    Parameter params;
+    Parameter p;
-    Solver solver;
+    Discretization d;
-    initParameter(&params);
+    Solver s;
+    initParameter(&p);
 
     if (argc != 2) {
         printf("Usage: %s <configFile>\n", argv[0]);
         exit(EXIT_SUCCESS);
     }
 
-    readParameter(&params, argv[1]);
+    readParameter(&p, argv[1]);
-    printParameter(&params);
+    printParameter(&p);
-    initSolver(&solver, &params);
+    initDiscretization(&d, &p);
+    initSolver(&s, &d, &p);
 #ifndef VERBOSE
-    initProgress(solver.te);
+    initProgress(d.te);
 #endif
 
-    double tau = solver.tau;
+    double tau = d.tau;
-    double te  = solver.te;
+    double te  = d.te;
     double t   = 0.0;
     int nt     = 0;
 
-    S = getTimeStamp();
+    timeStart = getTimeStamp();
     while (t <= te) {
-        if (tau > 0.0) computeTimestep(&solver);
+        if (tau > 0.0) computeTimestep(&d);
-        setBoundaryConditions(&solver);
+        setBoundaryConditions(&d);
-        setSpecialBoundaryCondition(&solver);
+        setSpecialBoundaryCondition(&d);
-        computeFG(&solver);
+        computeFG(&d);
-        computeRHS(&solver);
+        computeRHS(&d);
-        if (nt % 100 == 0) normalizePressure(&solver);
+        if (nt % 100 == 0) normalizePressure(&d);
-        solveRB(&solver);
+        solve(&s, d.p, d.rhs);
-        adaptUV(&solver);
+        adaptUV(&d);
-        t += solver.dt;
+        t += d.dt;
         nt++;
 
 #ifdef VERBOSE
@@ -58,9 +59,9 @@ int main(int argc, char** argv)
         printProgress(t);
 #endif
     }
-    E = getTimeStamp();
+    timeStop = getTimeStamp();
     stopProgress();
-    printf("Solution took %.2fs\n", E - S);
+    printf("Solution took %.2fs\n", timeStop - timeStart);
-    writeResult(&solver);
+    writeResult(&d);
     return EXIT_SUCCESS;
 }

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

@@ -24,6 +24,8 @@ void initParameter(Parameter* param)
     param->re      = 100.0;
     param->gamma   = 0.9;
     param->tau     = 0.5;
+    param->rho     = 0.99;
+    param->levels  = 5;
 }
 
 void readParameter(Parameter* param, const char* filename)
@@ -61,6 +63,7 @@ void readParameter(Parameter* param, const char* filename)
             PARSE_INT(imax);
             PARSE_INT(jmax);
             PARSE_INT(itermax);
+            PARSE_INT(levels);
             PARSE_REAL(eps);
             PARSE_REAL(omg);
             PARSE_REAL(re);
@@ -78,6 +81,7 @@ void readParameter(Parameter* param, const char* filename)
             PARSE_REAL(u_init);
             PARSE_REAL(v_init);
             PARSE_REAL(p_init);
+            PARSE_REAL(rho);
         }
     }
 
@@ -108,4 +112,6 @@ void printParameter(Parameter* param)
     printf("\tepsilon (stopping tolerance) : %f\n", param->eps);
     printf("\tgamma (stopping tolerance) : %f\n", param->gamma);
     printf("\tomega (SOR relaxation): %f\n", param->omg);
+    printf("\trho (SOR relaxation): %f\n", param->rho);
+    printf("\tMultiGrid levels : %d\n", param->levels);
 }

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

@@ -10,8 +10,8 @@
 typedef struct {
     double xlength, ylength;
     int imax, jmax;
-    int itermax;
+    int itermax, levels;
-    double eps, omg;
+    double eps, omg, rho;
     double re, tau, gamma;
     double te, dt;
     double gx, gy;

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

@@ -0,0 +1,192 @@
+/*
+ * 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 <stdio.h>
+#include <stdlib.h>
+
+#include "allocate.h"
+#include "discretization.h"
+#include "parameter.h"
+#include "solver.h"
+#include "util.h"
+
+#define FINEST_LEVEL 0
+#define S(i, j)      s[(j) * (imax + 2) + (i)]
+#define E(i, j)      e[(j) * (imax + 2) + (i)]
+#define R(i, j)      r[(j) * (imax + 2) + (i)]
+#define OLD(i, j)    old[(j) * (imax + 2) + (i)]
+
+static void restrictMG(Solver* s, int level, int imax, int jmax)
+{
+    double* r   = s->r[level + 1];
+    double* old = s->r[level];
+
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 1; i++) {
+            R(i, j) = (OLD(2 * i - 1, 2 * j - 1) + OLD(2 * i, 2 * j - 1) * 2 +
+                          OLD(2 * i + 1, 2 * j - 1) + OLD(2 * i - 1, 2 * j) * 2 +
+                          OLD(2 * i, 2 * j) * 4 + OLD(2 * i + 1, 2 * j) * 2 +
+                          OLD(2 * i - 1, 2 * j + 1) + OLD(2 * i, 2 * j + 1) * 2 +
+                          OLD(2 * i + 1, 2 * j + 1)) /
+                      16.0;
+        }
+    }
+}
+
+static void prolongate(Solver* s, int level, int imax, int jmax)
+{
+    double* old = s->r[level + 1];
+    double* e   = s->r[level];
+
+    for (int j = 2; j < jmax + 1; j += 2) {
+        for (int i = 2; i < imax + 1; i += 2) {
+            E(i, j) = OLD(i / 2, j / 2);
+        }
+    }
+}
+
+static void correct(Solver* s, double* p, int level, int imax, int jmax)
+{
+    double* e = s->e[level];
+    for (int j = 1; j < jmax + 1; ++j) {
+        for (int i = 1; i < imax + 1; ++i) {
+            P(i, j) += E(i, j);
+        }
+    }
+}
+
+static void setBoundaryCondition(double* p, int imax, int jmax)
+{
+    for (int i = 1; i < imax + 1; i++) {
+        P(i, 0)        = P(i, 1);
+        P(i, jmax + 1) = P(i, jmax);
+    }
+
+    for (int j = 1; j < jmax + 1; j++) {
+        P(0, j)        = P(1, j);
+        P(imax + 1, j) = P(imax, j);
+    }
+}
+
+static double smooth(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;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = isw; i < imax + 1; i += 2) {
+
+                R(i, j) = 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(i, j));
+                res += (R(i, j) * R(i, j));
+            }
+            isw = 3 - isw;
+        }
+        jsw = 3 - jsw;
+    }
+
+    res = res / (double)(imax * jmax);
+    return res;
+}
+
+void initSolver(Solver* s, Discretization* d, Parameter* p)
+{
+    s->eps     = p->eps;
+    s->omega   = p->omg;
+    s->itermax = p->itermax;
+    s->rho     = p->rho;
+    s->levels  = p->levels;
+    s->grid    = &d->grid;
+
+    int imax   = s->grid->imax;
+    int jmax   = s->grid->jmax;
+    int levels = s->levels;
+    printf("Using Multigrid solver with %d levels\n", levels);
+
+    s->r = malloc(levels * sizeof(double*));
+    s->e = malloc(levels * sizeof(double*));
+
+    size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
+
+    for (int j = 0; j < levels; j++) {
+        s->r[j] = allocate(64, size);
+        s->e[j] = allocate(64, size);
+
+        for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
+            s->r[j][i] = 0.0;
+            s->e[j][i] = 0.0;
+        }
+    }
+}
+
+double multiGrid(Solver* solver, double* p, double* rhs, int level, int imax, int jmax)
+{
+    double res = 0.0;
+
+    // coarsest level TODO: Use direct solver?
+    if (level == (solver->levels - 1)) {
+        for (int i = 0; i < 5; i++) {
+            smooth(solver, p, rhs, level, imax, jmax);
+        }
+        return res;
+    }
+
+    // pre-smoothing TODO: Make smoothing steps configurable?
+    for (int i = 0; i < 5; i++) {
+        smooth(solver, p, rhs, level, imax, jmax);
+        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
+    }
+
+    // restrict
+    restrictMG(solver, level, imax, jmax);
+
+    // MGSolver on residual and error.
+    // TODO: What if there is a rest?
+    multiGrid(solver,
+        solver->e[level + 1],
+        solver->r[level],
+        level + 1,
+        imax / 2,
+        jmax / 2);
+
+    // prolongate
+    prolongate(solver, level, imax, jmax);
+
+    // correct p on finer level using residual
+    correct(solver, p, level, imax, jmax);
+    if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
+
+    // post-smoothing
+    for (int i = 0; i < 5; i++) {
+        res = smooth(solver, p, rhs, level, imax, jmax);
+        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
+    }
+
+    return res;
+}
+
+void solve(Solver* s, double* p, double* rhs)
+{
+    double res = multiGrid(s, p, rhs, 0, s->grid->imax, s->grid->jmax);
+
+#ifdef VERBOSE
+    printf("Residuum: %.6f\n", res);
+#endif
+}

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

@@ -0,0 +1,129 @@
+/*
+ * 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 "discretization.h"
+#include "solver.h"
+#include "util.h"
+
+void initSolver(Solver* s, Discretization* d, Parameter* p)
+{
+    s->grid    = &d->grid;
+    s->itermax = p->itermax;
+    s->eps     = p->eps;
+    s->omega   = p->omg;
+}
+
+void solveSOR(Solver* solver, double* p, double* rhs)
+{
+    int imax      = solver->grid->imax;
+    int jmax      = solver->grid->jmax;
+    double eps    = solver->eps;
+    int itermax   = solver->itermax;
+    double dx2    = solver->grid->dx * solver->grid->dx;
+    double dy2    = solver->grid->dy * solver->grid->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double epssq  = eps * eps;
+    int it        = 0;
+    double res    = 1.0;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+
+                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);
+            }
+        }
+
+        for (int i = 1; i < imax + 1; i++) {
+            P(i, 0)        = P(i, 1);
+            P(i, jmax + 1) = P(i, jmax);
+        }
+
+        for (int j = 1; j < jmax + 1; j++) {
+            P(0, j)        = P(1, j);
+            P(imax + 1, j) = P(imax, j);
+        }
+
+        res = res / (double)(imax * jmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, res);
+#endif
+        it++;
+    }
+
+#ifdef VERBOSE
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+#endif
+}
+
+void solve(Solver* solver, double* p, double* rhs)
+{
+    int imax      = solver->grid->imax;
+    int jmax      = solver->grid->jmax;
+    double eps    = solver->eps;
+    int itermax   = solver->itermax;
+    double dx2    = solver->grid->dx * solver->grid->dx;
+    double dy2    = solver->grid->dy * solver->grid->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double epssq  = eps * eps;
+    int it        = 0;
+    double res    = 1.0;
+    int pass, jsw, isw;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+        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) {
+
+                    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;
+        }
+
+        for (int i = 1; i < imax + 1; i++) {
+            P(i, 0)        = P(i, 1);
+            P(i, jmax + 1) = P(i, jmax);
+        }
+
+        for (int j = 1; j < jmax + 1; j++) {
+            P(0, j)        = P(1, j);
+            P(imax + 1, j) = P(imax, j);
+        }
+
+        res = res / (double)(imax * jmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, res);
+#endif
+        it++;
+    }
+
+#ifdef VERBOSE
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+#endif
+}

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

@@ -1,564 +0,0 @@
-/*
- * 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 "parameter.h"
-#include "solver.h"
-#include "util.h"
-
-#define P(i, j)   p[(j) * (imax + 2) + (i)]
-#define F(i, j)   f[(j) * (imax + 2) + (i)]
-#define G(i, j)   g[(j) * (imax + 2) + (i)]
-#define U(i, j)   u[(j) * (imax + 2) + (i)]
-#define V(i, j)   v[(j) * (imax + 2) + (i)]
-#define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
-
-static void print(Solver* solver, double* grid)
-{
-    int imax = solver->imax;
-
-    for (int j = 0; j < solver->jmax + 2; j++) {
-        printf("%02d: ", j);
-        for (int i = 0; i < solver->imax + 2; i++) {
-            printf("%12.8f  ", grid[j * (imax + 2) + i]);
-        }
-        printf("\n");
-    }
-    fflush(stdout);
-}
-
-static void printConfig(Solver* solver)
-{
-    printf("Parameters for #%s#\n", solver->problem);
-    printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n",
-        solver->bcLeft,
-        solver->bcRight,
-        solver->bcBottom,
-        solver->bcTop);
-    printf("\tReynolds number: %.2f\n", solver->re);
-    printf("\tGx Gy: %.2f %.2f\n", solver->gx, solver->gy);
-    printf("Geometry data:\n");
-    printf("\tDomain box size (x, y): %.2f, %.2f\n", solver->xlength, solver->ylength);
-    printf("\tCells (x, y): %d, %d\n", solver->imax, solver->jmax);
-    printf("Timestep parameters:\n");
-    printf("\tDefault stepsize: %.2f, Final time %.2f\n", solver->dt, solver->te);
-    printf("\tdt bound: %.6f\n", solver->dtBound);
-    printf("\tTau factor: %.2f\n", solver->tau);
-    printf("Iterative solver parameters:\n");
-    printf("\tMax iterations: %d\n", solver->itermax);
-    printf("\tepsilon (stopping tolerance) : %f\n", solver->eps);
-    printf("\tgamma factor: %f\n", solver->gamma);
-    printf("\tomega (SOR relaxation): %f\n", solver->omega);
-}
-
-void initSolver(Solver* solver, Parameter* params)
-{
-    solver->problem  = params->name;
-    solver->bcLeft   = params->bcLeft;
-    solver->bcRight  = params->bcRight;
-    solver->bcBottom = params->bcBottom;
-    solver->bcTop    = params->bcTop;
-    solver->imax     = params->imax;
-    solver->jmax     = params->jmax;
-    solver->xlength  = params->xlength;
-    solver->ylength  = params->ylength;
-    solver->dx       = params->xlength / params->imax;
-    solver->dy       = params->ylength / params->jmax;
-    solver->eps      = params->eps;
-    solver->omega    = params->omg;
-    solver->itermax  = params->itermax;
-    solver->re       = params->re;
-    solver->gx       = params->gx;
-    solver->gy       = params->gy;
-    solver->dt       = params->dt;
-    solver->te       = params->te;
-    solver->tau      = params->tau;
-    solver->gamma    = params->gamma;
-
-    int imax    = solver->imax;
-    int jmax    = solver->jmax;
-    size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
-    solver->u   = allocate(64, size);
-    solver->v   = allocate(64, size);
-    solver->p   = allocate(64, size);
-    solver->rhs = allocate(64, size);
-    solver->f   = allocate(64, size);
-    solver->g   = allocate(64, size);
-
-    for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
-        solver->u[i]   = params->u_init;
-        solver->v[i]   = params->v_init;
-        solver->p[i]   = params->p_init;
-        solver->rhs[i] = 0.0;
-        solver->f[i]   = 0.0;
-        solver->g[i]   = 0.0;
-    }
-
-    double dx        = solver->dx;
-    double dy        = solver->dy;
-    double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
-    solver->dtBound  = 0.5 * solver->re * 1.0 / invSqrSum;
-#ifdef VERBOSE
-    printConfig(solver);
-#endif
-}
-
-void computeRHS(Solver* solver)
-{
-    int imax    = solver->imax;
-    int jmax    = solver->jmax;
-    double idx  = 1.0 / solver->dx;
-    double idy  = 1.0 / solver->dy;
-    double idt  = 1.0 / solver->dt;
-    double* rhs = solver->rhs;
-    double* f   = solver->f;
-    double* g   = solver->g;
-
-    for (int j = 1; j < jmax + 1; j++) {
-        for (int i = 1; i < imax + 1; i++) {
-            RHS(i, j) = idt *
-                        ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy);
-        }
-    }
-}
-
-void solve(Solver* solver)
-{
-    int imax      = solver->imax;
-    int jmax      = solver->jmax;
-    double eps    = solver->eps;
-    int itermax   = solver->itermax;
-    double dx2    = solver->dx * solver->dx;
-    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* p     = solver->p;
-    double* rhs   = solver->rhs;
-    double epssq  = eps * eps;
-    int it        = 0;
-    double res    = 1.0;
-
-    while ((res >= epssq) && (it < itermax)) {
-        res = 0.0;
-
-        for (int j = 1; j < jmax + 1; j++) {
-            for (int i = 1; i < imax + 1; i++) {
-
-                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);
-            }
-        }
-
-        for (int i = 1; i < imax + 1; i++) {
-            P(i, 0)        = P(i, 1);
-            P(i, jmax + 1) = P(i, jmax);
-        }
-
-        for (int j = 1; j < jmax + 1; j++) {
-            P(0, j)        = P(1, j);
-            P(imax + 1, j) = P(imax, j);
-        }
-
-        res = res / (double)(imax * jmax);
-#ifdef DEBUG
-        printf("%d Residuum: %e\n", it, res);
-#endif
-        it++;
-    }
-
-#ifdef VERBOSE
-    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
-#endif
-}
-
-void solveRB(Solver* solver)
-{
-    int imax      = solver->imax;
-    int jmax      = solver->jmax;
-    double eps    = solver->eps;
-    int itermax   = solver->itermax;
-    double dx2    = solver->dx * solver->dx;
-    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* p     = solver->p;
-    double* rhs   = solver->rhs;
-    double epssq  = eps * eps;
-    int it        = 0;
-    double res    = 1.0;
-    int pass, jsw, isw;
-
-    while ((res >= epssq) && (it < itermax)) {
-        res = 0.0;
-        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) {
-
-                    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;
-        }
-
-        for (int i = 1; i < imax + 1; i++) {
-            P(i, 0)        = P(i, 1);
-            P(i, jmax + 1) = P(i, jmax);
-        }
-
-        for (int j = 1; j < jmax + 1; j++) {
-            P(0, j)        = P(1, j);
-            P(imax + 1, j) = P(imax, j);
-        }
-
-        res = res / (double)(imax * jmax);
-#ifdef DEBUG
-        printf("%d Residuum: %e\n", it, res);
-#endif
-        it++;
-    }
-
-#ifdef VERBOSE
-    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
-#endif
-}
-
-static double maxElement(Solver* solver, double* m)
-{
-    int size      = (solver->imax + 2) * (solver->jmax + 2);
-    double maxval = DBL_MIN;
-
-    for (int i = 0; i < size; i++) {
-        maxval = MAX(maxval, fabs(m[i]));
-    }
-
-    return maxval;
-}
-
-void normalizePressure(Solver* solver)
-{
-    int size    = (solver->imax + 2) * (solver->jmax + 2);
-    double* p   = solver->p;
-    double avgP = 0.0;
-
-    for (int i = 0; i < size; i++) {
-        avgP += p[i];
-    }
-    avgP /= size;
-
-    for (int i = 0; i < size; i++) {
-        p[i] = p[i] - avgP;
-    }
-}
-
-void computeTimestep(Solver* solver)
-{
-    double dt   = solver->dtBound;
-    double dx   = solver->dx;
-    double dy   = solver->dy;
-    double umax = maxElement(solver, solver->u);
-    double vmax = maxElement(solver, solver->v);
-
-    if (umax > 0) {
-        dt = (dt > dx / umax) ? dx / umax : dt;
-    }
-    if (vmax > 0) {
-        dt = (dt > dy / vmax) ? dy / vmax : dt;
-    }
-
-    solver->dt = dt * solver->tau;
-}
-
-void setBoundaryConditions(Solver* solver)
-{
-    int imax  = solver->imax;
-    int jmax  = solver->jmax;
-    double* u = solver->u;
-    double* v = solver->v;
-
-    // Left boundary
-    switch (solver->bcLeft) {
-    case NOSLIP:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(0, j) = 0.0;
-            V(0, j) = -V(1, j);
-        }
-        break;
-    case SLIP:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(0, j) = 0.0;
-            V(0, j) = V(1, j);
-        }
-        break;
-    case OUTFLOW:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(0, j) = U(1, j);
-            V(0, j) = V(1, j);
-        }
-        break;
-    case PERIODIC:
-        break;
-    }
-
-    // Right boundary
-    switch (solver->bcRight) {
-    case NOSLIP:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(imax, j)     = 0.0;
-            V(imax + 1, j) = -V(imax, j);
-        }
-        break;
-    case SLIP:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(imax, j)     = 0.0;
-            V(imax + 1, j) = V(imax, j);
-        }
-        break;
-    case OUTFLOW:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(imax, j)     = U(imax - 1, j);
-            V(imax + 1, j) = V(imax, j);
-        }
-        break;
-    case PERIODIC:
-        break;
-    }
-
-    // Bottom boundary
-    switch (solver->bcBottom) {
-    case NOSLIP:
-        for (int i = 1; i < imax + 1; i++) {
-            V(i, 0) = 0.0;
-            U(i, 0) = -U(i, 1);
-        }
-        break;
-    case SLIP:
-        for (int i = 1; i < imax + 1; i++) {
-            V(i, 0) = 0.0;
-            U(i, 0) = U(i, 1);
-        }
-        break;
-    case OUTFLOW:
-        for (int i = 1; i < imax + 1; i++) {
-            U(i, 0) = U(i, 1);
-            V(i, 0) = V(i, 1);
-        }
-        break;
-    case PERIODIC:
-        break;
-    }
-
-    // Top boundary
-    switch (solver->bcTop) {
-    case NOSLIP:
-        for (int i = 1; i < imax + 1; i++) {
-            V(i, jmax)     = 0.0;
-            U(i, jmax + 1) = -U(i, jmax);
-        }
-        break;
-    case SLIP:
-        for (int i = 1; i < imax + 1; i++) {
-            V(i, jmax)     = 0.0;
-            U(i, jmax + 1) = U(i, jmax);
-        }
-        break;
-    case OUTFLOW:
-        for (int i = 1; i < imax + 1; i++) {
-            U(i, jmax + 1) = U(i, jmax);
-            V(i, jmax)     = V(i, jmax - 1);
-        }
-        break;
-    case PERIODIC:
-        break;
-    }
-}
-
-void setSpecialBoundaryCondition(Solver* solver)
-{
-    int imax   = solver->imax;
-    int jmax   = solver->jmax;
-    double mDy = solver->dy;
-    double* u  = solver->u;
-
-    if (strcmp(solver->problem, "dcavity") == 0) {
-        for (int i = 1; i < imax; i++) {
-            U(i, jmax + 1) = 2.0 - U(i, jmax);
-        }
-    } else if (strcmp(solver->problem, "canal") == 0) {
-        double ylength = solver->ylength;
-        double y;
-
-        for (int j = 1; j < jmax + 1; j++) {
-            y       = mDy * (j - 0.5);
-            U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
-        }
-    }
-}
-
-void computeFG(Solver* solver)
-{
-    double* u        = solver->u;
-    double* v        = solver->v;
-    double* f        = solver->f;
-    double* g        = solver->g;
-    int imax         = solver->imax;
-    int jmax         = solver->jmax;
-    double gx        = solver->gx;
-    double gy        = solver->gy;
-    double gamma     = solver->gamma;
-    double dt        = solver->dt;
-    double inverseRe = 1.0 / solver->re;
-    double inverseDx = 1.0 / solver->dx;
-    double inverseDy = 1.0 / solver->dy;
-    double du2dx, dv2dy, duvdx, duvdy;
-    double du2dx2, du2dy2, dv2dx2, dv2dy2;
-
-    for (int j = 1; j < jmax + 1; j++) {
-        for (int i = 1; i < imax + 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 --------------------------- */
-    for (int j = 1; j < jmax + 1; j++) {
-        F(0, j)    = U(0, j);
-        F(imax, j) = U(imax, j);
-    }
-
-    /* ---------------------- boundary of G --------------------------- */
-    for (int i = 1; i < imax + 1; i++) {
-        G(i, 0)    = V(i, 0);
-        G(i, jmax) = V(i, jmax);
-    }
-}
-
-void adaptUV(Solver* solver)
-{
-    int imax       = solver->imax;
-    int jmax       = solver->jmax;
-    double* p      = solver->p;
-    double* u      = solver->u;
-    double* v      = solver->v;
-    double* f      = solver->f;
-    double* g      = solver->g;
-    double factorX = solver->dt / solver->dx;
-    double factorY = solver->dt / solver->dy;
-
-    for (int j = 1; j < jmax + 1; j++) {
-        for (int i = 1; i < imax + 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* solver)
-{
-    int imax  = solver->imax;
-    int jmax  = solver->jmax;
-    double dx = solver->dx;
-    double dy = solver->dy;
-    double* p = solver->p;
-    double* u = solver->u;
-    double* v = solver->v;
-    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 + 1; j++) {
-        y = (double)(j - 0.5) * dy;
-        for (int i = 1; i < imax + 1; i++) {
-            x = (double)(i - 0.5) * dx;
-            fprintf(fp, "%.2f %.2f %f\n", x, y, P(i, j));
-        }
-        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 + 1; j++) {
-        y = dy * (j - 0.5);
-        for (int i = 1; i < imax + 1; i++) {
-            x            = dx * (i - 0.5);
-            double vel_u = (U(i, j) + U(i - 1, j)) / 2.0;
-            double vel_v = (V(i, j) + V(i, j - 1)) / 2.0;
-            double len   = sqrt((vel_u * vel_u) + (vel_v * vel_v));
-            fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, vel_u, vel_v, len);
-        }
-    }
-
-    fclose(fp);
-}

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

@@ -6,41 +6,21 @@
  */
 #ifndef __SOLVER_H_
 #define __SOLVER_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 eps, omega, rho;
-    double re, tau, gamma;
-    double gx, gy;
-    /* time stepping */
     int itermax;
-    double dt, te;
+    int levels;
-    double dtBound;
+    double **r, **e;
-    char* problem;
-    int bcLeft, bcRight, bcBottom, bcTop;
 } Solver;
 
-extern void initSolver(Solver*, Parameter*);
+extern void initSolver(Solver*, Discretization*, Parameter*);
-extern void computeRHS(Solver*);
+extern void solve(Solver*, double*, double*);
-extern void solve(Solver*);
+
-extern void solveRB(Solver*);
-extern void solveRBA(Solver*);
-extern void normalizePressure(Solver*);
-extern void computeTimestep(Solver*);
-extern void setBoundaryConditions(Solver*);
-extern void setSpecialBoundaryCondition(Solver*);
-extern void computeFG(Solver*);
-extern void adaptUV(Solver*);
-extern void writeResult(Solver*);
 #endif

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

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