diff --git a/.gitignore b/.gitignore index 71520c4..ebea632 100644 --- a/.gitignore +++ b/.gitignore @@ -53,3 +53,5 @@ Module.symvers Mkfile.old dkms.conf +*.mp4 + diff --git a/PoissonSolver/2D-seq/Makefile b/PoissonSolver/2D-seq/Makefile index 53fcb50..640561f 100644 --- a/PoissonSolver/2D-seq/Makefile +++ b/PoissonSolver/2D-seq/Makefile @@ -10,6 +10,7 @@ TARGET = exe-$(TAG) BUILD_DIR = ./$(TAG) SRC_DIR = ./src MAKE_DIR = ./ +VIDEO_DIR = ./video Q ?= @ #DO NOT EDIT BELOW @@ -22,7 +23,7 @@ ASM = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s,$(wildcard $(SRC_DIR)/*. OBJ = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o,$(wildcard $(SRC_DIR)/*.c)) CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES) -${TARGET}: $(BUILD_DIR) $(OBJ) +${TARGET}: $(BUILD_DIR) $(VIDEO_DIR) $(OBJ) $(info ===> LINKING $(TARGET)) $(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS) @@ -35,17 +36,22 @@ $(BUILD_DIR)/%.s: %.c $(info ===> GENERATE ASM $@) $(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@ -.PHONY: clean distclean tags info asm +.PHONY: clean distclean mkvid tags info asm clean: $(info ===> CLEAN) @rm -rf $(BUILD_DIR) + @rm -rf $(VIDEO_DIR) @rm -f tags distclean: clean $(info ===> DIST CLEAN) @rm -f $(TARGET) +mkvid: + gnuplot animate.plot + ffmpeg -framerate 10 -i ./video/%04d.png -vcodec libx264 -crf 25 -pix_fmt yuv420p out.mp4 + info: $(info $(CFLAGS)) $(Q)$(CC) $(VERSION) @@ -59,4 +65,7 @@ tags: $(BUILD_DIR): @mkdir $(BUILD_DIR) +$(VIDEO_DIR): + @mkdir $(VIDEO_DIR) + -include $(OBJ:.o=.d) diff --git a/PoissonSolver/2D-seq/animate.plot b/PoissonSolver/2D-seq/animate.plot index 4082e54..80f3830 100644 --- a/PoissonSolver/2D-seq/animate.plot +++ b/PoissonSolver/2D-seq/animate.plot @@ -4,12 +4,13 @@ set grid set hidden3d set xrange [0:40] set yrange [0:40] -set zrange [-2:2] -input(n) = sprintf("p-%d.dat", n) -output(n) = sprintf("%03d.png", n) +input(n) = sprintf("./video/%04d.dat", n) +output(n) = sprintf("./video/%04d.png", n) +title(n) = sprintf("Iteration %d", n) -do for [i=1:50] { - set output output(i) +do for [i=0:400] { + set output output(i) + set title title(i) splot input(i) matrix using 1:2:3 with lines } diff --git a/PoissonSolver/2D-seq/src/main.c b/PoissonSolver/2D-seq/src/main.c index a5fe991..945f5bc 100644 --- a/PoissonSolver/2D-seq/src/main.c +++ b/PoissonSolver/2D-seq/src/main.c @@ -17,7 +17,7 @@ LIKWID_MARKER_STOP(#tag); \ endTime = getTimeStamp(); -enum VARIANT { SOR = 1, RB, RBA }; +enum VARIANT { JACOBI = 1, GS, SOR, RB, RBA }; int main(int argc, char** argv) { @@ -43,12 +43,20 @@ int main(int argc, char** argv) } initSolver(&solver, ¶ms, 2); - writeResult(&solver, "p-init.dat"); + writeResult(&solver, "./video/0000.dat"); switch (variant) { + case JACOBI: + printf("JACOBI\n"); + LIKWID_PROFILE("JACOBI", solveJacobi); + break; + case GS: + printf("GAUSS-SEIDEL\n"); + LIKWID_PROFILE("GS", solveGS); + break; case SOR: printf("Plain SOR\n"); - LIKWID_PROFILE("SOR", solve); + LIKWID_PROFILE("SOR", solveSOR); break; case RB: printf("Red-black SOR\n"); diff --git a/PoissonSolver/2D-seq/src/solver.c b/PoissonSolver/2D-seq/src/solver.c index 06dda9b..01875e9 100644 --- a/PoissonSolver/2D-seq/src/solver.c +++ b/PoissonSolver/2D-seq/src/solver.c @@ -10,9 +10,11 @@ #include "parameter.h" #include "solver.h" -#define PI 3.14159265358979323846 -#define P(i, j) p[(j) * (imax + 2) + (i)] -#define RHS(i, j) rhs[(j) * (imax + 2) + (i)] +#define PI 3.14159265358979323846 +#define P(i, j) p[(j) * (imax + 2) + (i)] +#define Pnew(i, j) pNew[(j) * (imax + 2) + (i)] +#define Pold(i, j) pOld[(j) * (imax + 2) + (i)] +#define RHS(i, j) rhs[(j) * (imax + 2) + (i)] void initSolver(Solver* solver, Parameter* params, int problem) { @@ -57,7 +59,135 @@ void initSolver(Solver* solver, Parameter* params, int problem) } } -void solve(Solver* solver) +static inline void swapPtr(double** p1, double** p2) +{ + double* tmp = *p1; + *p1 = *p2; + *p2 = tmp; +} + +void solveJacobi(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 = 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; + size_t bytesize = (imax + 2) * (jmax + 2) * sizeof(double); + double* pOld = p; + double* pNew = allocate(64, bytesize); + + 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) - + ((Pold(i - 1, j) - 2.0 * Pold(i, j) + Pold(i + 1, j)) * idx2 + + (Pold(i, j - 1) - 2.0 * Pold(i, j) + Pold(i, j + 1)) * + idy2); + + Pnew(i, j) -= factor * r; + res += (r * r); + } + } + + for (int i = 1; i < imax + 1; i++) { + Pnew(i, 0) = Pnew(i, 1); + Pnew(i, jmax + 1) = Pnew(i, jmax); + } + + for (int j = 1; j < jmax + 1; j++) { + Pnew(0, j) = Pnew(1, j); + Pnew(imax + 1, j) = Pnew(imax, j); + } + + swapPtr(&pNew, &pOld); + res = res / (double)(imax * jmax); +#ifdef DEBUG + printf("%d Residuum: %e\n", it, res); +#endif +#ifdef ANIMATE + char filename[50]; + sprintf(filename, "./video/%04d.dat", it); + writeResult(solver, filename); +#endif + it++; + } + + free(pNew); + printf("%d\n", it); +} + +void solveGS(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 = 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; + size_t bytesize = (imax + 2) * (jmax + 2) * sizeof(double); + + 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, sqrt(res)); +#endif +#ifdef ANIMATE + char filename[50]; + sprintf(filename, "video/%04d.dat", it); + writeResult(solver, filename); +#endif + it++; + } + + printf("%d\n", it); +} + +void solveSOR(Solver* solver) { int imax = solver->imax; int jmax = solver->jmax; @@ -73,7 +203,6 @@ void solve(Solver* solver) double epssq = eps * eps; int it = 0; double res = 1.0; - char filename[20]; while ((res >= epssq) && (it < itermax)) { res = 0.0; @@ -105,7 +234,8 @@ void solve(Solver* solver) printf("%d Residuum: %e\n", it, res); #endif #ifdef ANIMATE - sprintf(filename, "p-%d.dat", it); + char filename[50]; + sprintf(filename, "./video/%04d.dat", it); writeResult(solver, filename); #endif it++; @@ -169,7 +299,8 @@ void solveRB(Solver* solver) printf("%d Residuum: %e\n", it, res); #endif #ifdef ANIMATE - sprintf(filename, "p-%d.dat", it); + char filename[50]; + sprintf(filename, "./video/%04d.dat", it); writeResult(solver, filename); #endif it++; @@ -238,7 +369,8 @@ void solveRBA(Solver* solver) printf("%d Residuum: %e Omega: %e\n", it, res, omega); #endif #ifdef ANIMATE - sprintf(filename, "p-%d.dat", it); + char filename[50]; + sprintf(filename, "video/%04d.dat", it); writeResult(solver, filename); #endif it++; @@ -259,7 +391,7 @@ void writeResult(Solver* solver, char* filename) fp = fopen(filename, "w"); if (fp == NULL) { - printf("Error!\n"); + printf("Error opening file %s!\n", filename); exit(EXIT_FAILURE); } diff --git a/PoissonSolver/2D-seq/src/solver.h b/PoissonSolver/2D-seq/src/solver.h index c84099b..1bda6d8 100644 --- a/PoissonSolver/2D-seq/src/solver.h +++ b/PoissonSolver/2D-seq/src/solver.h @@ -16,7 +16,9 @@ typedef struct { extern void initSolver(Solver*, Parameter*, int problem); extern void writeResult(Solver*, char*); -extern void solve(Solver*); +extern void solveJacobi(Solver*); +extern void solveGS(Solver*); +extern void solveSOR(Solver*); extern void solveRB(Solver*); extern void solveRBA(Solver*); #endif