diff --git a/BasicSolver/2D-seq/src/affinity.c b/BasicSolver/2D-seq/src/affinity.c deleted file mode 100644 index ef8355a..0000000 --- a/BasicSolver/2D-seq/src/affinity.c +++ /dev/null @@ -1,61 +0,0 @@ -/* - * 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. - */ -#ifdef __linux__ -#ifdef _OPENMP -#include -#include -#include -#include -#include -#include -#include - -#define MAX_NUM_THREADS 128 -#define gettid() syscall(SYS_gettid) - -static int getProcessorID(cpu_set_t* cpu_set) -{ - int processorId; - - for (processorId = 0; processorId < MAX_NUM_THREADS; processorId++) { - if (CPU_ISSET(processorId, cpu_set)) { - break; - } - } - return processorId; -} - -int affinity_getProcessorId() -{ - cpu_set_t cpu_set; - CPU_ZERO(&cpu_set); - sched_getaffinity(gettid(), sizeof(cpu_set_t), &cpu_set); - - return getProcessorID(&cpu_set); -} - -void affinity_pinThread(int processorId) -{ - cpu_set_t cpuset; - pthread_t thread; - - thread = pthread_self(); - CPU_ZERO(&cpuset); - CPU_SET(processorId, &cpuset); - pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset); -} - -void affinity_pinProcess(int processorId) -{ - cpu_set_t cpuset; - - CPU_ZERO(&cpuset); - CPU_SET(processorId, &cpuset); - sched_setaffinity(0, sizeof(cpu_set_t), &cpuset); -} -#endif /*_OPENMP*/ -#endif /*__linux__*/ diff --git a/BasicSolver/2D-seq/src/affinity.h b/BasicSolver/2D-seq/src/affinity.h deleted file mode 100644 index 84bf733..0000000 --- a/BasicSolver/2D-seq/src/affinity.h +++ /dev/null @@ -1,14 +0,0 @@ -/* - * 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. - */ -#ifndef AFFINITY_H -#define AFFINITY_H - -extern int affinity_getProcessorId(); -extern void affinity_pinProcess(int); -extern void affinity_pinThread(int); - -#endif /*AFFINITY_H*/ diff --git a/BasicSolver/2D-seq/src/parameter.c b/BasicSolver/2D-seq/src/parameter.c index 777a274..704f24e 100644 --- a/BasicSolver/2D-seq/src/parameter.c +++ b/BasicSolver/2D-seq/src/parameter.c @@ -24,7 +24,6 @@ void initParameter(Parameter* param) param->re = 100.0; param->gamma = 0.9; param->tau = 0.5; - param->rho = 0.99; param->levels = 5; } @@ -81,7 +80,6 @@ void readParameter(Parameter* param, const char* filename) PARSE_REAL(u_init); PARSE_REAL(v_init); PARSE_REAL(p_init); - PARSE_REAL(rho); } } @@ -112,6 +110,5 @@ 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); } diff --git a/BasicSolver/2D-seq/src/progress.h b/BasicSolver/2D-seq/src/progress.h index 1cdb9a3..314e921 100644 --- a/BasicSolver/2D-seq/src/progress.h +++ b/BasicSolver/2D-seq/src/progress.h @@ -9,6 +9,6 @@ extern void initProgress(double); extern void printProgress(double); -extern void stopProgress(); +extern void stopProgress(void); #endif diff --git a/BasicSolver/2D-seq/src/solver-mg.c b/BasicSolver/2D-seq/src/solver-mg.c index 721be56..ccd884c 100644 --- a/BasicSolver/2D-seq/src/solver-mg.c +++ b/BasicSolver/2D-seq/src/solver-mg.c @@ -8,16 +8,15 @@ #include #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)] +#define FINEST_LEVEL 0 +#define COARSEST_LEVEL (solver->levels - 1) +#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) { @@ -51,6 +50,7 @@ static void prolongate(Solver* s, int level, int imax, int jmax) 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); @@ -106,37 +106,8 @@ static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int 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) +static double multiGrid( + Solver* solver, double* p, double* rhs, int level, int imax, int jmax) { double res = 0.0; @@ -161,7 +132,7 @@ double multiGrid(Solver* solver, double* p, double* rhs, int level, int imax, in // TODO: What if there is a rest? multiGrid(solver, solver->e[level + 1], - solver->r[level], + solver->r[level + 1], level + 1, imax / 2, jmax / 2); @@ -182,6 +153,35 @@ double multiGrid(Solver* solver, double* p, double* rhs, int level, int imax, in return res; } +void initSolver(Solver* s, Discretization* d, Parameter* p) +{ + s->eps = p->eps; + s->omega = p->omg; + s->itermax = p->itermax; + 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; + } + } +} + void solve(Solver* s, double* p, double* rhs) { double res = multiGrid(s, p, rhs, 0, s->grid->imax, s->grid->jmax); diff --git a/BasicSolver/2D-seq/src/solver-sor.c b/BasicSolver/2D-seq/src/solver-sor.c index b516e5b..f817b5f 100644 --- a/BasicSolver/2D-seq/src/solver-sor.c +++ b/BasicSolver/2D-seq/src/solver-sor.c @@ -4,7 +4,6 @@ * 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" diff --git a/BasicSolver/3D-seq/Makefile b/BasicSolver/3D-seq/Makefile index 84c5eff..64735f6 100644 --- a/BasicSolver/3D-seq/Makefile +++ b/BasicSolver/3D-seq/Makefile @@ -18,9 +18,10 @@ include $(MAKE_DIR)/include_$(TAG).mk INCLUDES += -I$(SRC_DIR) -I$(BUILD_DIR) VPATH = $(SRC_DIR) -SRC = $(wildcard $(SRC_DIR)/*.c) +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) @@ -47,6 +48,8 @@ clean: distclean: clean $(info ===> DIST CLEAN) @rm -f $(TARGET) + @rm -f *.dat + @rm -f *.png info: $(info $(CFLAGS)) diff --git a/BasicSolver/3D-seq/config.mk b/BasicSolver/3D-seq/config.mk index af5f1a0..ccdfa2c 100644 --- a/BasicSolver/3D-seq/config.mk +++ b/BasicSolver/3D-seq/config.mk @@ -1,12 +1,12 @@ # Supported: GCC, CLANG, ICC TAG ?= CLANG ENABLE_OPENMP ?= false +# Supported: sor, mg +SOLVER ?= mg +# 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 diff --git a/BasicSolver/3D-seq/dcavity.par b/BasicSolver/3D-seq/dcavity.par index 6d1ef6d..8224887 100644 --- a/BasicSolver/3D-seq/dcavity.par +++ b/BasicSolver/3D-seq/dcavity.par @@ -47,6 +47,8 @@ tau 0.5 # safety factor for time stepsize control (<0 constant delt) itermax 1000 # maximal number of pressure iteration in one time step eps 0.001 # stopping tolerance for pressure iteration +rho 0.5 omg 1.7 # relaxation parameter for SOR iteration gamma 0.9 # upwind differencing factor gamma +levels 5 # Multigrid levels #=============================================================================== diff --git a/BasicSolver/3D-seq/include_CLANG.mk b/BasicSolver/3D-seq/include_CLANG.mk index 5abaf06..a124053 100644 --- a/BasicSolver/3D-seq/include_CLANG.mk +++ b/BasicSolver/3D-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 #-Weverything -#CFLAGS = -Ofast -fnt-store=aggressive -std=c99 $(OPENMP) #AMD CLANG LFLAGS = $(OPENMP) -lm -DEFINES = -D_GNU_SOURCE# -DDEBUG +DEFINES = -D_GNU_SOURCE INCLUDES = diff --git a/BasicSolver/3D-seq/src/solver.c b/BasicSolver/3D-seq/src/discretization.c similarity index 66% rename from BasicSolver/3D-seq/src/solver.c rename to BasicSolver/3D-seq/src/discretization.c index 3f9cb5e..eacac17 100644 --- a/BasicSolver/3D-seq/src/solver.c +++ b/BasicSolver/3D-seq/src/discretization.c @@ -1,5 +1,5 @@ /* - * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. + * 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. @@ -10,131 +10,120 @@ #include #include "allocate.h" +#include "discretization.h" #include "parameter.h" -#include "solver.h" #include "util.h" -#define P(i, j, k) p[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] -#define F(i, j, k) f[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] -#define G(i, j, k) g[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] -#define H(i, j, k) h[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] -#define U(i, j, k) u[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] -#define V(i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] -#define W(i, j, k) w[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] -#define RHS(i, j, k) rhs[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] - -static void printConfig(Solver* s) +static void printConfig(Discretization* d) { - printf("Parameters for #%s#\n", s->problem); + printf("Parameters for #%s#\n", d->problem); printf("BC Left:%d Right:%d Bottom:%d Top:%d Front:%d Back:%d\n", - s->bcLeft, - s->bcRight, - s->bcBottom, - s->bcTop, - s->bcFront, - s->bcBack); - printf("\tReynolds number: %.2f\n", s->re); - printf("\tGx Gy: %.2f %.2f %.2f\n", s->gx, s->gy, s->gz); + d->bcLeft, + d->bcRight, + d->bcBottom, + d->bcTop, + d->bcFront, + d->bcBack); + printf("\tReynolds number: %.2f\n", d->re); + printf("\tGx Gy: %.2f %.2f %.2f\n", d->gx, d->gy, d->gz); printf("Geometry data:\n"); printf("\tDomain box size (x, y, z): %.2f, %.2f, %.2f\n", - s->grid.xlength, - s->grid.ylength, - s->grid.zlength); - printf("\tCells (x, y, z): %d, %d, %d\n", s->grid.imax, s->grid.jmax, s->grid.kmax); - printf("\tCell size (dx, dy, dz): %f, %f, %f\n", s->grid.dx, s->grid.dy, s->grid.dz); + d->grid.xlength, + d->grid.ylength, + d->grid.zlength); + printf("\tCells (x, y, z): %d, %d, %d\n", d->grid.imax, d->grid.jmax, d->grid.kmax); + printf("\tCell size (dx, dy, dz): %f, %f, %f\n", d->grid.dx, d->grid.dy, d->grid.dz); 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("\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 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); + printf("\tepsilon (stopping tolerance) : %f\n", d->eps); + printf("\tgamma factor: %f\n", d->gamma); + printf("\tomega (SOR relaxation): %f\n", d->omega); } -void initSolver(Solver* s, Parameter* params) +void initDiscretization(Discretization* d, Parameter* p) { - s->problem = params->name; - s->bcLeft = params->bcLeft; - s->bcRight = params->bcRight; - s->bcBottom = params->bcBottom; - s->bcTop = params->bcTop; - s->bcFront = params->bcFront; - s->bcBack = params->bcBack; + d->problem = p->name; + d->bcLeft = p->bcLeft; + d->bcRight = p->bcRight; + d->bcBottom = p->bcBottom; + d->bcTop = p->bcTop; + d->bcFront = p->bcFront; + d->bcBack = p->bcBack; - s->grid.imax = params->imax; - s->grid.jmax = params->jmax; - s->grid.kmax = params->kmax; - s->grid.xlength = params->xlength; - s->grid.ylength = params->ylength; - s->grid.zlength = params->zlength; - s->grid.dx = params->xlength / params->imax; - s->grid.dy = params->ylength / params->jmax; - s->grid.dz = params->zlength / params->kmax; + d->grid.imax = p->imax; + d->grid.jmax = p->jmax; + d->grid.kmax = p->kmax; + d->grid.xlength = p->xlength; + d->grid.ylength = p->ylength; + d->grid.zlength = p->zlength; + d->grid.dx = p->xlength / p->imax; + d->grid.dy = p->ylength / p->jmax; + d->grid.dz = p->zlength / p->kmax; - 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->gz = params->gz; - s->dt = params->dt; - s->te = params->te; - s->tau = params->tau; - s->gamma = params->gamma; + d->eps = p->eps; + d->omega = p->omg; + d->re = p->re; + d->gx = p->gx; + d->gy = p->gy; + d->gz = p->gz; + d->dt = p->dt; + d->te = p->te; + d->tau = p->tau; + d->gamma = p->gamma; - int imax = s->grid.imax; - int jmax = s->grid.jmax; - int kmax = s->grid.kmax; + int imax = d->grid.imax; + int jmax = d->grid.jmax; + int kmax = d->grid.kmax; size_t bytesize = (imax + 2) * (jmax + 2) * (kmax + 2) * sizeof(double); - s->u = allocate(64, bytesize); - s->v = allocate(64, bytesize); - s->w = allocate(64, bytesize); - s->p = allocate(64, bytesize); - s->rhs = allocate(64, bytesize); - s->f = allocate(64, bytesize); - s->g = allocate(64, bytesize); - s->h = allocate(64, bytesize); + d->u = allocate(64, bytesize); + d->v = allocate(64, bytesize); + d->w = allocate(64, bytesize); + d->p = allocate(64, bytesize); + d->rhs = allocate(64, bytesize); + d->f = allocate(64, bytesize); + d->g = allocate(64, bytesize); + d->h = allocate(64, bytesize); for (int i = 0; i < (imax + 2) * (jmax + 2) * (kmax + 2); i++) { - s->u[i] = params->u_init; - s->v[i] = params->v_init; - s->w[i] = params->w_init; - s->p[i] = params->p_init; - s->rhs[i] = 0.0; - s->f[i] = 0.0; - s->g[i] = 0.0; - s->h[i] = 0.0; + d->u[i] = p->u_init; + d->v[i] = p->v_init; + d->w[i] = p->w_init; + d->p[i] = p->p_init; + d->rhs[i] = 0.0; + d->f[i] = 0.0; + d->g[i] = 0.0; + d->h[i] = 0.0; } - double dx = s->grid.dx; - double dy = s->grid.dy; - double dz = s->grid.dz; + double dx = d->grid.dx; + double dy = d->grid.dy; + double dz = d->grid.dz; double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy) + 1.0 / (dz * dz); - s->dtBound = 0.5 * s->re * 1.0 / invSqrSum; + d->dtBound = 0.5 * d->re * 1.0 / invSqrSum; #ifdef VERBOSE printConfig(s); #endif /* VERBOSE */ } -void computeRHS(Solver* s) +void computeRHS(Discretization* d) { - int imax = s->grid.imax; - int jmax = s->grid.jmax; - int kmax = s->grid.kmax; - double idx = 1.0 / s->grid.dx; - double idy = 1.0 / s->grid.dy; - double idz = 1.0 / s->grid.dz; - double idt = 1.0 / s->dt; + int imax = d->grid.imax; + int jmax = d->grid.jmax; + int kmax = d->grid.kmax; + double idx = 1.0 / d->grid.dx; + double idy = 1.0 / d->grid.dy; + double idz = 1.0 / d->grid.dz; + double idt = 1.0 / d->dt; - double* rhs = s->rhs; - double* f = s->f; - double* g = s->g; - double* h = s->h; + double* rhs = d->rhs; + double* f = d->f; + double* g = d->g; + double* h = d->h; for (int k = 1; k < kmax + 1; k++) { for (int j = 1; j < jmax + 1; j++) { @@ -148,94 +137,9 @@ void computeRHS(Solver* s) } } -void solve(Solver* s) +static double maxElement(Discretization* d, double* m) { - int imax = s->grid.imax; - int jmax = s->grid.jmax; - int kmax = s->grid.kmax; - double eps = s->eps; - int itermax = s->itermax; - double dx2 = s->grid.dx * s->grid.dx; - double dy2 = s->grid.dy * s->grid.dy; - double dz2 = s->grid.dz * s->grid.dz; - double idx2 = 1.0 / dx2; - double idy2 = 1.0 / dy2; - double idz2 = 1.0 / dz2; - double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) / - (dy2 * dz2 + dx2 * dz2 + dx2 * dy2); - double* p = s->p; - double* rhs = s->rhs; - double epssq = eps * eps; - int it = 0; - double res = 1.0; - int pass, ksw, jsw, isw; - - while ((res >= epssq) && (it < itermax)) { - res = 0.0; - ksw = 1; - - for (pass = 0; pass < 2; pass++) { - jsw = ksw; - - for (int k = 1; k < kmax + 1; k++) { - isw = jsw; - for (int j = 1; j < jmax + 1; j++) { - for (int i = isw; i < imax + 1; i += 2) { - - double r = - RHS(i, j, k) - - ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + - (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * - idy2 + - (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * - idz2); - - P(i, j, k) -= (factor * r); - res += (r * r); - } - isw = 3 - isw; - } - jsw = 3 - jsw; - } - ksw = 3 - ksw; - } - - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; i++) { - P(i, j, 0) = P(i, j, 1); - P(i, j, kmax + 1) = P(i, j, kmax); - } - } - - for (int k = 1; k < kmax + 1; k++) { - for (int i = 1; i < imax + 1; i++) { - P(i, 0, k) = P(i, 1, k); - P(i, jmax + 1, k) = P(i, jmax, k); - } - } - - for (int k = 1; k < kmax + 1; k++) { - for (int j = 1; j < jmax + 1; j++) { - P(0, j, k) = P(1, j, k); - P(imax + 1, j, k) = P(imax, j, k); - } - } - - res = res / (double)(imax * jmax * kmax); -#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* s, double* m) -{ - int size = (s->grid.imax + 2) * (s->grid.jmax + 2) * (s->grid.kmax + 2); + int size = (d->grid.imax + 2) * (d->grid.jmax + 2) * (d->grid.kmax + 2); double maxval = DBL_MIN; for (int i = 0; i < size; i++) { @@ -245,10 +149,10 @@ static double maxElement(Solver* s, double* m) return maxval; } -void normalizePressure(Solver* s) +void normalizePressure(Discretization* d) { - int size = (s->grid.imax + 2) * (s->grid.jmax + 2) * (s->grid.kmax + 2); - double* p = s->p; + int size = (d->grid.imax + 2) * (d->grid.jmax + 2) * (d->grid.kmax + 2); + double* p = d->p; double avgP = 0.0; for (int i = 0; i < size; i++) { @@ -261,16 +165,16 @@ void normalizePressure(Solver* s) } } -void computeTimestep(Solver* s) +void computeTimestep(Discretization* d) { - double dt = s->dtBound; - double dx = s->grid.dx; - double dy = s->grid.dy; - double dz = s->grid.dz; + double dt = d->dtBound; + double dx = d->grid.dx; + double dy = d->grid.dy; + double dz = d->grid.dz; - double umax = maxElement(s, s->u); - double vmax = maxElement(s, s->v); - double wmax = maxElement(s, s->w); + double umax = maxElement(d, d->u); + double vmax = maxElement(d, d->v); + double wmax = maxElement(d, d->w); if (umax > 0) { dt = (dt > dx / umax) ? dx / umax : dt; @@ -282,20 +186,20 @@ void computeTimestep(Solver* s) dt = (dt > dz / wmax) ? dz / wmax : dt; } - s->dt = dt * s->tau; + d->dt = dt * d->tau; } -void setBoundaryConditions(Solver* s) +void setBoundaryConditions(Discretization* d) { - int imax = s->grid.imax; - int jmax = s->grid.jmax; - int kmax = s->grid.kmax; + int imax = d->grid.imax; + int jmax = d->grid.jmax; + int kmax = d->grid.kmax; - double* u = s->u; - double* v = s->v; - double* w = s->w; + double* u = d->u; + double* v = d->v; + double* w = d->w; - switch (s->bcTop) { + switch (d->bcTop) { case NOSLIP: for (int k = 1; k < kmax + 1; k++) { for (int i = 1; i < imax + 1; i++) { @@ -327,7 +231,7 @@ void setBoundaryConditions(Solver* s) break; } - switch (s->bcBottom) { + switch (d->bcBottom) { case NOSLIP: for (int k = 1; k < kmax + 1; k++) { for (int i = 1; i < imax + 1; i++) { @@ -359,7 +263,7 @@ void setBoundaryConditions(Solver* s) break; } - switch (s->bcLeft) { + switch (d->bcLeft) { case NOSLIP: for (int k = 1; k < kmax + 1; k++) { for (int j = 1; j < jmax + 1; j++) { @@ -391,7 +295,7 @@ void setBoundaryConditions(Solver* s) break; } - switch (s->bcRight) { + switch (d->bcRight) { case NOSLIP: for (int k = 1; k < kmax + 1; k++) { for (int j = 1; j < jmax + 1; j++) { @@ -423,7 +327,7 @@ void setBoundaryConditions(Solver* s) break; } - switch (s->bcFront) { + switch (d->bcFront) { case NOSLIP: for (int j = 1; j < jmax + 1; j++) { for (int i = 1; i < imax + 1; i++) { @@ -455,7 +359,7 @@ void setBoundaryConditions(Solver* s) break; } - switch (s->bcBack) { + switch (d->bcBack) { case NOSLIP: for (int j = 1; j < jmax + 1; j++) { for (int i = 1; i < imax + 1; i++) { @@ -488,23 +392,23 @@ void setBoundaryConditions(Solver* s) } } -void setSpecialBoundaryCondition(Solver* s) +void setSpecialBoundaryCondition(Discretization* d) { - int imax = s->grid.imax; - int jmax = s->grid.jmax; - int kmax = s->grid.kmax; + int imax = d->grid.imax; + int jmax = d->grid.jmax; + int kmax = d->grid.kmax; - double mDy = s->grid.dy; - double* u = s->u; + double mDy = d->grid.dy; + double* u = d->u; - if (strcmp(s->problem, "dcavity") == 0) { + if (strcmp(d->problem, "dcavity") == 0) { for (int k = 1; k < kmax; k++) { for (int i = 1; i < imax; i++) { U(i, jmax + 1, k) = 2.0 - U(i, jmax, k); } } - } else if (strcmp(s->problem, "canal") == 0) { - double ylength = s->grid.ylength; + } else if (strcmp(d->problem, "canal") == 0) { + double ylength = d->grid.ylength; double y; for (int k = 1; k < kmax + 1; k++) { @@ -516,29 +420,29 @@ void setSpecialBoundaryCondition(Solver* s) } } -void computeFG(Solver* s) +void computeFG(Discretization* d) { - int imax = s->grid.imax; - int jmax = s->grid.jmax; - int kmax = s->grid.kmax; + int imax = d->grid.imax; + int jmax = d->grid.jmax; + int kmax = d->grid.kmax; - double* u = s->u; - double* v = s->v; - double* w = s->w; - double* f = s->f; - double* g = s->g; - double* h = s->h; + double* u = d->u; + double* v = d->v; + double* w = d->w; + double* f = d->f; + double* g = d->g; + double* h = d->h; - double gx = s->gx; - double gy = s->gy; - double gz = s->gz; - double dt = s->dt; + double gx = d->gx; + double gy = d->gy; + double gz = d->gz; + double dt = d->dt; - double gamma = s->gamma; - double inverseRe = 1.0 / s->re; - double inverseDx = 1.0 / s->grid.dx; - double inverseDy = 1.0 / s->grid.dy; - double inverseDz = 1.0 / s->grid.dz; + double gamma = d->gamma; + double inverseRe = 1.0 / d->re; + double inverseDx = 1.0 / d->grid.dx; + double inverseDy = 1.0 / d->grid.dy; + double inverseDz = 1.0 / d->grid.dz; double du2dx, dv2dy, dw2dz; double duvdx, duwdx, duvdy, dvwdy, duwdz, dvwdz; double du2dx2, du2dy2, du2dz2; @@ -705,23 +609,23 @@ void computeFG(Solver* s) } } -void adaptUV(Solver* s) +void adaptUV(Discretization* d) { - int imax = s->grid.imax; - int jmax = s->grid.jmax; - int kmax = s->grid.kmax; + int imax = d->grid.imax; + int jmax = d->grid.jmax; + int kmax = d->grid.kmax; - double* p = s->p; - double* u = s->u; - double* v = s->v; - double* w = s->w; - double* f = s->f; - double* g = s->g; - double* h = s->h; + double* p = d->p; + double* u = d->u; + double* v = d->v; + double* w = d->w; + double* f = d->f; + double* g = d->g; + double* h = d->h; - double factorX = s->dt / s->grid.dx; - double factorY = s->dt / s->grid.dy; - double factorZ = s->dt / s->grid.dz; + double factorX = d->dt / d->grid.dx; + double factorY = d->dt / d->grid.dy; + double factorZ = d->dt / d->grid.dz; for (int k = 1; k < kmax + 1; k++) { for (int j = 1; j < jmax + 1; j++) { diff --git a/BasicSolver/3D-seq/src/discretization.h b/BasicSolver/3D-seq/src/discretization.h new file mode 100644 index 0000000..17f97fd --- /dev/null +++ b/BasicSolver/3D-seq/src/discretization.h @@ -0,0 +1,41 @@ +/* + * 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, *h; + double *u, *v, *w; + /* parameters */ + double eps, omega; + double re, tau, gamma; + double gx, gy, gz; + /* time stepping */ + double dt, te; + double dtBound; + char* problem; + int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack; +} 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*); +#endif diff --git a/BasicSolver/3D-seq/src/main.c b/BasicSolver/3D-seq/src/main.c index 789e69d..c289631 100644 --- a/BasicSolver/3D-seq/src/main.c +++ b/BasicSolver/3D-seq/src/main.c @@ -9,6 +9,7 @@ #include #include "allocate.h" +#include "discretization.h" #include "parameter.h" #include "progress.h" #include "solver.h" @@ -17,7 +18,8 @@ #define G(v, i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] -static void createBulkArrays(Solver* s, double* pg, double* ug, double* vg, double* wg) +static void createBulkArrays( + Discretization* s, double* pg, double* ug, double* vg, double* wg) { int imax = s->grid.imax; int jmax = s->grid.jmax; @@ -67,6 +69,7 @@ int main(int argc, char** argv) { double timeStart, timeStop; Parameter p; + Discretization d; Solver s; initParameter(&p); @@ -77,51 +80,53 @@ int main(int argc, char** argv) readParameter(&p, argv[1]); printParameter(&p); - initSolver(&s, &p); + initDiscretization(&d, &p); + initSolver(&s, &d, &p); #ifndef VERBOSE - initProgress(s.te); + initProgress(d.te); #endif - double tau = s.tau; - double te = s.te; + double tau = d.tau; + double te = d.te; double t = 0.0; + int nt = 0; timeStart = getTimeStamp(); while (t <= te) { - if (tau > 0.0) computeTimestep(&s); - setBoundaryConditions(&s); - setSpecialBoundaryCondition(&s); - computeFG(&s); - computeRHS(&s); - solve(&s); - adaptUV(&s); - t += s.dt; + if (tau > 0.0) computeTimestep(&d); + setBoundaryConditions(&d); + setSpecialBoundaryCondition(&d); + computeFG(&d); + computeRHS(&d); + if (nt % 100 == 0) normalizePressure(&d); + solve(&s, d.p, d.rhs); + adaptUV(&d); + t += d.dt; + nt++; #ifdef VERBOSE - printf("TIME %f , TIMESTEP %f\n", t, s.dt); + printf("TIME %f , TIMESTEP %f\n", t, solver.dt); #else printProgress(t); #endif } timeStop = getTimeStamp(); -#ifndef VERBOSE stopProgress(); -#endif printf("Solution took %.2fs\n", timeStop - timeStart); timeStart = getTimeStamp(); double *pg, *ug, *vg, *wg; - size_t bytesize = (size_t)(s.grid.imax * s.grid.jmax * s.grid.kmax) * sizeof(double); + size_t bytesize = (size_t)(d.grid.imax * d.grid.jmax * d.grid.kmax) * sizeof(double); pg = allocate(64, bytesize); ug = allocate(64, bytesize); vg = allocate(64, bytesize); wg = allocate(64, bytesize); - createBulkArrays(&s, pg, ug, vg, wg); - VtkOptions opts = { .grid = s.grid }; - vtkOpen(&opts, s.problem); + createBulkArrays(&d, pg, ug, vg, wg); + VtkOptions opts = { .grid = d.grid }; + vtkOpen(&opts, d.problem); vtkScalar(&opts, "pressure", pg); vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg }); vtkClose(&opts); diff --git a/BasicSolver/3D-seq/src/parameter.c b/BasicSolver/3D-seq/src/parameter.c index 2128c2a..005b69a 100644 --- a/BasicSolver/3D-seq/src/parameter.c +++ b/BasicSolver/3D-seq/src/parameter.c @@ -26,6 +26,7 @@ void initParameter(Parameter* param) param->re = 100.0; param->gamma = 0.9; param->tau = 0.5; + param->levels = 5; } void readParameter(Parameter* param, const char* filename) @@ -65,6 +66,7 @@ void readParameter(Parameter* param, const char* filename) PARSE_INT(jmax); PARSE_INT(kmax); PARSE_INT(itermax); + PARSE_INT(levels); PARSE_REAL(eps); PARSE_REAL(omg); PARSE_REAL(re); @@ -123,4 +125,5 @@ 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("\tMultiGrid levels : %d\n", param->levels); } diff --git a/BasicSolver/3D-seq/src/parameter.h b/BasicSolver/3D-seq/src/parameter.h index d946d5c..83917c0 100644 --- a/BasicSolver/3D-seq/src/parameter.h +++ b/BasicSolver/3D-seq/src/parameter.h @@ -10,8 +10,8 @@ typedef struct { int imax, jmax, kmax; double xlength, ylength, zlength; - int itermax; - double eps, omg; + int itermax, levels; + double eps, omg, rho; double re, tau, gamma; double te, dt; double gx, gy, gz; diff --git a/BasicSolver/3D-seq/src/solver-mg.c b/BasicSolver/3D-seq/src/solver-mg.c new file mode 100644 index 0000000..9644f59 --- /dev/null +++ b/BasicSolver/3D-seq/src/solver-mg.c @@ -0,0 +1,250 @@ +/* + * 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 +#include + +#include "allocate.h" +#include "solver.h" +#include "util.h" + +#define FINEST_LEVEL 0 +#define COARSEST_LEVEL (s->levels - 1) +#define S(i, j, k) s[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define E(i, j, k) e[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define R(i, j, k) r[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define OLD(i, j, k) old[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] + +static void restrictMG(Solver* s, int level, int imax, int jmax, int kmax) +{ + double* r = s->r[level + 1]; + double* old = s->r[level]; + + for (int k = 1; k < kmax + 1; k++) { + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; ++i) { + R(i, j, k) = (OLD(2 * i - 1, 2 * j - 1, 2 * k) + + OLD(2 * i, 2 * j - 1, 2 * k) * 2 + + OLD(2 * i + 1, 2 * j - 1, 2 * k) + + OLD(2 * i - 1, 2 * j, 2 * k) * 2 + + OLD(2 * i, 2 * j, 2 * k) * 8 + + OLD(2 * i + 1, 2 * j, 2 * k) * 2 + + OLD(2 * i - 1, 2 * j + 1, 2 * k) + + OLD(2 * i, 2 * j + 1, 2 * k) * 2 + + OLD(2 * i + 1, 2 * j + 1, 2 * k) + + + OLD(2 * i - 1, 2 * j - 1, 2 * k - 1) + + OLD(2 * i, 2 * j - 1, 2 * k - 1) * 2 + + OLD(2 * i + 1, 2 * j - 1, 2 * k - 1) + + OLD(2 * i - 1, 2 * j, 2 * k - 1) * 2 + + OLD(2 * i, 2 * j, 2 * k - 1) * 4 + + OLD(2 * i + 1, 2 * j, 2 * k - 1) * 2 + + OLD(2 * i - 1, 2 * j + 1, 2 * k - 1) + + OLD(2 * i, 2 * j + 1, 2 * k - 1) * 2 + + OLD(2 * i + 1, 2 * j + 1, 2 * k - 1) + + + OLD(2 * i - 1, 2 * j - 1, 2 * k + 1) + + OLD(2 * i, 2 * j - 1, 2 * k + 1) * 2 + + OLD(2 * i + 1, 2 * j - 1, 2 * k + 1) + + OLD(2 * i - 1, 2 * j, 2 * k + 1) * 2 + + OLD(2 * i, 2 * j, 2 * k + 1) * 4 + + OLD(2 * i + 1, 2 * j, 2 * k + 1) * 2 + + OLD(2 * i - 1, 2 * j + 1, 2 * k + 1) + + OLD(2 * i, 2 * j + 1, 2 * k + 1) * 2 + + OLD(2 * i + 1, 2 * j + 1, 2 * k + 1)) / + 64.0; + } + } + } +} + +static void prolongate(Solver* s, int level, int imax, int jmax, int kmax) +{ + double* old = s->r[level + 1]; + double* e = s->r[level]; + + for (int k = 2; k < kmax + 1; k += 2) { + for (int j = 2; j < jmax + 1; j += 2) { + for (int i = 2; i < imax + 1; i += 2) { + E(i, j, k) = OLD(i / 2, j / 2, k / 2); + } + } + } +} + +static void correct(Solver* s, double* p, int level, int imax, int jmax, int kmax) +{ + double* e = s->e[level]; + + for (int k = 1; k < kmax + 1; ++k) { + for (int j = 1; j < jmax + 1; ++j) { + for (int i = 1; i < imax + 1; ++i) { + P(i, j, k) += E(i, j, k); + } + } + } +} + +static void setBoundaryCondition(double* p, int imax, int jmax, int kmax) +{ + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + P(i, j, 0) = P(i, j, 1); + P(i, j, kmax + 1) = P(i, j, kmax); + } + } + + for (int k = 1; k < kmax + 1; k++) { + for (int i = 1; i < imax + 1; i++) { + P(i, 0, k) = P(i, 1, k); + P(i, jmax + 1, k) = P(i, jmax, k); + } + } + + for (int k = 1; k < kmax + 1; k++) { + for (int j = 1; j < jmax + 1; j++) { + P(0, j, k) = P(1, j, k); + P(imax + 1, j, k) = P(imax, j, k); + } + } +} + +static double smooth( + Solver* s, double* p, double* rhs, int level, int imax, int jmax, int kmax) +{ + double eps = s->eps; + int itermax = s->itermax; + double dx2 = s->grid->dx * s->grid->dx; + double dy2 = s->grid->dy * s->grid->dy; + double dz2 = s->grid->dz * s->grid->dz; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double idz2 = 1.0 / dz2; + double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) / + (dy2 * dz2 + dx2 * dz2 + dx2 * dy2); + double* r = s->r[level]; + double epssq = eps * eps; + int it = 0; + int pass, ksw, jsw, isw; + double res = 1.0; + + ksw = 1; + + for (pass = 0; pass < 2; pass++) { + jsw = ksw; + + for (int k = 1; k < kmax + 1; k++) { + isw = jsw; + for (int j = 1; j < jmax + 1; j++) { + for (int i = isw; i < imax + 1; i += 2) { + + R(i, j, k) = + RHS(i, j, k) - + ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + + (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * idy2 + + (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * idz2); + + P(i, j, k) -= (factor * R(i, j, k)); + res += (R(i, j, k) * R(i, j, k)); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + ksw = 3 - ksw; + } + + res = res / (double)(imax * jmax * kmax); + + return res; +} + +static double multiGrid( + Solver* s, double* p, double* rhs, int level, int imax, int jmax, int kmax) +{ + double res = 0.0; + + // coarsest level TODO: Use direct solver? + if (level == COARSEST_LEVEL) { + for (int i = 0; i < 5; i++) { + smooth(s, p, rhs, level, imax, jmax, kmax); + } + return res; + } + + // pre-smoothing TODO: Make smoothing steps configurable? + for (int i = 0; i < 5; i++) { + smooth(s, p, rhs, level, imax, jmax, kmax); + if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax); + } + + // restrict + restrictMG(s, level, imax, jmax, kmax); + + // 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, + kmax / 2); + + // prolongate + prolongate(s, level, imax, jmax, kmax); + + // correct p on finer level using residual + correct(s, p, level, imax, jmax, kmax); + if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax); + + // post-smoothing + for (int i = 0; i < 5; i++) { + res = smooth(s, p, rhs, level, imax, jmax, kmax); + if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax); + } + + return res; +} + +void initSolver(Solver* s, Discretization* d, Parameter* p) +{ + s->eps = p->eps; + s->omega = p->omg; + s->itermax = p->itermax; + s->levels = p->levels; + s->grid = &d->grid; + + int imax = s->grid->imax; + int jmax = s->grid->jmax; + int kmax = s->grid->kmax; + 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) * (kmax + 2); + + for (int j = 0; j < levels; j++) { + s->r[j] = allocate(64, size * sizeof(double)); + s->e[j] = allocate(64, size * sizeof(double)); + + for (size_t i = 0; i < size; i++) { + s->r[j][i] = 0.0; + s->e[j][i] = 0.0; + } + } +} + +void solve(Solver* s, double* p, double* rhs) +{ + double res = multiGrid(s, p, rhs, 0, s->grid->imax, s->grid->jmax, s->grid->kmax); + +#ifdef VERBOSE + printf("Residuum: %.6f\n", res); +#endif +} diff --git a/BasicSolver/3D-seq/src/solver-sor.c b/BasicSolver/3D-seq/src/solver-sor.c new file mode 100644 index 0000000..6f954f5 --- /dev/null +++ b/BasicSolver/3D-seq/src/solver-sor.c @@ -0,0 +1,99 @@ +/* + * 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 "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 solve(Solver* s, double* p, double* rhs) +{ + int imax = s->grid->imax; + int jmax = s->grid->jmax; + int kmax = s->grid->kmax; + double eps = s->eps; + int itermax = s->itermax; + double dx2 = s->grid->dx * s->grid->dx; + double dy2 = s->grid->dy * s->grid->dy; + double dz2 = s->grid->dz * s->grid->dz; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double idz2 = 1.0 / dz2; + double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) / + (dy2 * dz2 + dx2 * dz2 + dx2 * dy2); + double epssq = eps * eps; + int it = 0; + double res = 1.0; + int pass, ksw, jsw, isw; + + while ((res >= epssq) && (it < itermax)) { + res = 0.0; + ksw = 1; + + for (pass = 0; pass < 2; pass++) { + jsw = ksw; + + for (int k = 1; k < kmax + 1; k++) { + isw = jsw; + for (int j = 1; j < jmax + 1; j++) { + for (int i = isw; i < imax + 1; i += 2) { + + double r = + RHS(i, j, k) - + ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + + (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * + idy2 + + (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * + idz2); + + P(i, j, k) -= (factor * r); + res += (r * r); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + ksw = 3 - ksw; + } + + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + P(i, j, 0) = P(i, j, 1); + P(i, j, kmax + 1) = P(i, j, kmax); + } + } + + for (int k = 1; k < kmax + 1; k++) { + for (int i = 1; i < imax + 1; i++) { + P(i, 0, k) = P(i, 1, k); + P(i, jmax + 1, k) = P(i, jmax, k); + } + } + + for (int k = 1; k < kmax + 1; k++) { + for (int j = 1; j < jmax + 1; j++) { + P(0, j, k) = P(1, j, k); + P(imax + 1, j, k) = P(imax, j, k); + } + } + + res = res / (double)(imax * jmax * kmax); +#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 +} diff --git a/BasicSolver/3D-seq/src/solver.h b/BasicSolver/3D-seq/src/solver.h index c68ae9f..d7ce513 100644 --- a/BasicSolver/3D-seq/src/solver.h +++ b/BasicSolver/3D-seq/src/solver.h @@ -6,38 +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 */ - Grid grid; - /* arrays */ - double *p, *rhs; - double *f, *g, *h; - double *u, *v, *w; + Grid* grid; /* parameters */ - double eps, omega; - double re, tau, gamma; - double gx, gy, gz; - /* time stepping */ + double eps, omega, rho; int itermax; - double dt, te; - double dtBound; - char* problem; - int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack; + int levels; + double **r, **e; } Solver; -extern void initSolver(Solver*, Parameter*); -extern void computeRHS(Solver*); -extern void solve(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 initSolver(Solver*, Discretization*, Parameter*); +extern void solve(Solver*, double*, double*); + #endif diff --git a/BasicSolver/3D-seq/src/util.h b/BasicSolver/3D-seq/src/util.h index c27b2ba..1131237 100644 --- a/BasicSolver/3D-seq/src/util.h +++ b/BasicSolver/3D-seq/src/util.h @@ -19,4 +19,13 @@ #define ABS(a) ((a) >= 0 ? (a) : -(a)) #endif +#define P(i, j, k) p[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define F(i, j, k) f[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define G(i, j, k) g[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define H(i, j, k) h[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define U(i, j, k) u[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define V(i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define W(i, j, k) w[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] +#define RHS(i, j, k) rhs[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)] + #endif // __UTIL_H_