diff --git a/PoissonSolver/2D-mpi-sq/Makefile b/PoissonSolver/2D-mpi-sq/Makefile new file mode 100644 index 0000000..53fcb50 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/Makefile @@ -0,0 +1,62 @@ +#======================================================================================= +# Copyright (C) 2022 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. +#======================================================================================= + +#CONFIGURE BUILD SYSTEM +TARGET = exe-$(TAG) +BUILD_DIR = ./$(TAG) +SRC_DIR = ./src +MAKE_DIR = ./ +Q ?= @ + +#DO NOT EDIT BELOW +include $(MAKE_DIR)/config.mk +include $(MAKE_DIR)/include_$(TAG).mk +INCLUDES += -I$(SRC_DIR)/includes -I$(BUILD_DIR) + +VPATH = $(SRC_DIR) +ASM = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s,$(wildcard $(SRC_DIR)/*.c)) +OBJ = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o,$(wildcard $(SRC_DIR)/*.c)) +CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES) + +${TARGET}: $(BUILD_DIR) $(OBJ) + $(info ===> LINKING $(TARGET)) + $(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS) + +$(BUILD_DIR)/%.o: %.c $(MAKE_DIR)/include_$(TAG).mk + $(info ===> COMPILE $@) + $(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@ + $(Q)$(GCC) $(CPPFLAGS) -MT $(@:.d=.o) -MM $< > $(BUILD_DIR)/$*.d + +$(BUILD_DIR)/%.s: %.c + $(info ===> GENERATE ASM $@) + $(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@ + +.PHONY: clean distclean tags info asm + +clean: + $(info ===> CLEAN) + @rm -rf $(BUILD_DIR) + @rm -f tags + +distclean: clean + $(info ===> DIST CLEAN) + @rm -f $(TARGET) + +info: + $(info $(CFLAGS)) + $(Q)$(CC) $(VERSION) + +asm: $(BUILD_DIR) $(ASM) + +tags: + $(info ===> GENERATE TAGS) + $(Q)ctags -R + +$(BUILD_DIR): + @mkdir $(BUILD_DIR) + +-include $(OBJ:.o=.d) diff --git a/PoissonSolver/2D-mpi-sq/README.md b/PoissonSolver/2D-mpi-sq/README.md new file mode 100644 index 0000000..b0a80a6 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/README.md @@ -0,0 +1,48 @@ +# C source skeleton + +## Build + +1. Configure the toolchain and additional options in `config.mk`: +``` +# Supported: GCC, CLANG, ICC +TAG ?= GCC +ENABLE_OPENMP ?= false + +OPTIONS += -DARRAY_ALIGNMENT=64 +#OPTIONS += -DVERBOSE_AFFINITY +#OPTIONS += -DVERBOSE_DATASIZE +#OPTIONS += -DVERBOSE_TIMER +``` + +The verbosity options enable detailed output about affinity settings, allocation sizes and timer resolution. + + +2. Build with: +``` +make +``` + +You can build multiple toolchains in the same directory, but notice that the Makefile is only acting on the one currently set. +Intermediate build results are located in the `` directory. + +To output the executed commands use: +``` +make Q= +``` + +3. Clean up with: +``` +make clean +``` +to clean intermediate build results. + +``` +make distclean +``` +to clean intermediate build results and binary. + +4. (Optional) Generate assembler: +``` +make asm +``` +The assembler files will also be located in the `` directory. diff --git a/PoissonSolver/2D-mpi-sq/animate.plot b/PoissonSolver/2D-mpi-sq/animate.plot new file mode 100644 index 0000000..4082e54 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/animate.plot @@ -0,0 +1,15 @@ +set term png size 1024,768 enhanced font ,12 +set datafile separator whitespace +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) + +do for [i=1:50] { + set output output(i) + splot input(i) matrix using 1:2:3 with lines +} diff --git a/PoissonSolver/2D-mpi-sq/config.mk b/PoissonSolver/2D-mpi-sq/config.mk new file mode 100644 index 0000000..17f876a --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/config.mk @@ -0,0 +1,8 @@ +# Supported: GCC, CLANG, ICC +TAG ?= CLANG + +#Feature options +OPTIONS += -DARRAY_ALIGNMENT=64 +#OPTIONS += -DVERBOSE_AFFINITY +#OPTIONS += -DVERBOSE_DATASIZE +#OPTIONS += -DVERBOSE_TIMER diff --git a/PoissonSolver/2D-mpi-sq/include_CLANG.mk b/PoissonSolver/2D-mpi-sq/include_CLANG.mk new file mode 100644 index 0000000..c950c8c --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/include_CLANG.mk @@ -0,0 +1,18 @@ +CC = mpicc +GCC = cc +LINKER = $(CC) + +ifeq ($(ENABLE_OPENMP),true) +OPENMP = -fopenmp +#OPENMP = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp +LIBS = # -lomp +endif + +VERSION = --version +CFLAGS = -Ofast -std=c99 $(OPENMP) +#CFLAGS = -Ofast -fnt-store=aggressive -std=c99 $(OPENMP) #AMD CLANG +LFLAGS = $(OPENMP) -lm +DEFINES = -D_GNU_SOURCE +DEFINES += -DANIMATE +# DEFINES += -DDEBUG +INCLUDES = diff --git a/PoissonSolver/2D-mpi-sq/include_GCC.mk b/PoissonSolver/2D-mpi-sq/include_GCC.mk new file mode 100644 index 0000000..427e798 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/include_GCC.mk @@ -0,0 +1,14 @@ +CC = gcc +GCC = gcc +LINKER = $(CC) + +ifeq ($(ENABLE_OPENMP),true) +OPENMP = -fopenmp +endif + +VERSION = --version +CFLAGS = -Ofast -ffreestanding -std=c99 $(OPENMP) +LFLAGS = $(OPENMP) +DEFINES = -D_GNU_SOURCE +INCLUDES = +LIBS = diff --git a/PoissonSolver/2D-mpi-sq/include_ICC.mk b/PoissonSolver/2D-mpi-sq/include_ICC.mk new file mode 100644 index 0000000..94b8e20 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/include_ICC.mk @@ -0,0 +1,14 @@ +CC = icc +GCC = gcc +LINKER = $(CC) + +ifeq ($(ENABLE_OPENMP),true) +OPENMP = -qopenmp +endif + +VERSION = --version +CFLAGS = -O3 -xHost -qopt-zmm-usage=high -std=c99 $(OPENMP) +LFLAGS = $(OPENMP) +DEFINES = -D_GNU_SOURCE +INCLUDES = +LIBS = diff --git a/PoissonSolver/2D-mpi-sq/poisson.par b/PoissonSolver/2D-mpi-sq/poisson.par new file mode 100644 index 0000000..7f2301a --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/poisson.par @@ -0,0 +1,22 @@ +# Problem specific Data: +# --------------------- + +name poisson + +# Geometry Data: +# ------------- + +xlength 1.0 # domain size in x-direction +ylength 1.0 # domain size in y-direction +imax 100 # number of interior cells in x-direction +jmax 100 # number of interior cells in y-direction + +# Pressure Iteration Data: +# ----------------------- + +itermax 1000000 # maximal number of pressure iteration in one time step +eps 0.000001 # stopping tolerance for pressure iteration +rho 0.99999 # relaxation parameter for SOR iteration +omg 1.2 # relaxation parameter for SOR iteration + +#=============================================================================== diff --git a/PoissonSolver/2D-mpi-sq/src/allocate.c b/PoissonSolver/2D-mpi-sq/src/allocate.c new file mode 100644 index 0000000..513a693 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/allocate.c @@ -0,0 +1,37 @@ +/* + * Copyright (C) 2022 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. + */ +#include +#include +#include + +void* allocate (int alignment, size_t bytesize) +{ + int errorCode; + void* ptr; + + errorCode = posix_memalign(&ptr, alignment, bytesize); + + if (errorCode) { + if (errorCode == EINVAL) { + fprintf(stderr, + "Error: Alignment parameter is not a power of two\n"); + exit(EXIT_FAILURE); + } + if (errorCode == ENOMEM) { + fprintf(stderr, + "Error: Insufficient memory to fulfill the request\n"); + exit(EXIT_FAILURE); + } + } + + if (ptr == NULL) { + fprintf(stderr, "Error: posix_memalign failed!\n"); + exit(EXIT_FAILURE); + } + + return ptr; +} diff --git a/PoissonSolver/2D-mpi-sq/src/allocate.h b/PoissonSolver/2D-mpi-sq/src/allocate.h new file mode 100644 index 0000000..8e996bb --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/allocate.h @@ -0,0 +1,11 @@ +/* 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 __ALLOCATE_H_ +#define __ALLOCATE_H_ +#include + +extern void* allocate(int alignment, size_t bytesize); + +#endif diff --git a/PoissonSolver/2D-mpi-sq/src/likwid-marker.h b/PoissonSolver/2D-mpi-sq/src/likwid-marker.h new file mode 100644 index 0000000..a35b495 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/likwid-marker.h @@ -0,0 +1,53 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#ifndef LIKWID_MARKERS_H +#define LIKWID_MARKERS_H + +#ifdef LIKWID_PERFMON +#include +#define LIKWID_MARKER_INIT likwid_markerInit() +#define LIKWID_MARKER_THREADINIT likwid_markerThreadInit() +#define LIKWID_MARKER_SWITCH likwid_markerNextGroup() +#define LIKWID_MARKER_REGISTER(regionTag) likwid_markerRegisterRegion(regionTag) +#define LIKWID_MARKER_START(regionTag) likwid_markerStartRegion(regionTag) +#define LIKWID_MARKER_STOP(regionTag) likwid_markerStopRegion(regionTag) +#define LIKWID_MARKER_CLOSE likwid_markerClose() +#define LIKWID_MARKER_RESET(regionTag) likwid_markerResetRegion(regionTag) +#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count) likwid_markerGetRegion(regionTag, nevents, events, time, count) +#else /* LIKWID_PERFMON */ +#define LIKWID_MARKER_INIT +#define LIKWID_MARKER_THREADINIT +#define LIKWID_MARKER_SWITCH +#define LIKWID_MARKER_REGISTER(regionTag) +#define LIKWID_MARKER_START(regionTag) +#define LIKWID_MARKER_STOP(regionTag) +#define LIKWID_MARKER_CLOSE +#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count) +#define LIKWID_MARKER_RESET(regionTag) +#endif /* LIKWID_PERFMON */ + +#endif /*LIKWID_MARKERS_H*/ diff --git a/PoissonSolver/2D-mpi-sq/src/main.c b/PoissonSolver/2D-mpi-sq/src/main.c new file mode 100644 index 0000000..af09d50 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/main.c @@ -0,0 +1,46 @@ +/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.ke + * All rights reserved. + * Use of this source code is governed by a MIT-style + * license that can be found in the LICENSE file. */ +#include +#include +#include + +#include "parameter.h" +#include "solver.h" +#include "timing.h" + +int main(int argc, char** argv) +{ + + MPI_Init(&argc, &argv); + /******* Debug Params *********/ + int rank, size; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + printf("Process %d of %d \n", rank, size); + + double startTime, endTime; + Parameter params; + Solver solver; + initParameter(¶ms); + + if (argc < 2) { + printf("Usage: %s \n", argv[0]); + exit(EXIT_SUCCESS); + } + readParameter(¶ms, argv[1]); + initSolver(&solver, ¶ms, 2); + // writeResult(&solver, "p-0.dat"); + // startTime = getTimeStamp(); + // solveRB(&solver); + // endTime = getTimeStamp(); + // printf(" %.2fs\n", endTime - startTime); + // writeResult(&solver, "p.dat"); + + if (rank == 0) { + } + + MPI_Finalize(); + return EXIT_SUCCESS; +} diff --git a/PoissonSolver/2D-mpi-sq/src/parameter.c b/PoissonSolver/2D-mpi-sq/src/parameter.c new file mode 100644 index 0000000..0dec654 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/parameter.c @@ -0,0 +1,79 @@ +/* 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. */ +#include +#include +#include +//--- +#include "parameter.h" +#include "util.h" +#define MAXLINE 4096 + +void initParameter(Parameter* param) +{ + param->xlength = 1.0; + param->ylength = 1.0; + param->imax = 100; + param->jmax = 100; + param->itermax = 1000; + param->eps = 0.0001; + param->omg = 1.8; + param->rho = 0.99; +} + +void readParameter(Parameter* param, const char* filename) +{ + FILE* fp = fopen(filename, "r"); + char line[MAXLINE]; + int i; + + if (!fp) { + fprintf(stderr, "Could not open parameter file: %s\n", filename); + exit(EXIT_FAILURE); + } + + while (!feof(fp)) { + line[0] = '\0'; + fgets(line, MAXLINE, fp); + for (i = 0; line[i] != '\0' && line[i] != '#'; i++) + ; + line[i] = '\0'; + + char* tok = strtok(line, " "); + char* val = strtok(NULL, " "); + +#define PARSE_PARAM(p, f) \ + if (strncmp(tok, #p, sizeof(#p) / sizeof(#p[0]) - 1) == 0) { \ + param->p = f(val); \ + } +#define PARSE_STRING(p) PARSE_PARAM(p, strdup) +#define PARSE_INT(p) PARSE_PARAM(p, atoi) +#define PARSE_REAL(p) PARSE_PARAM(p, atof) + + if (tok != NULL && val != NULL) { + PARSE_REAL(xlength); + PARSE_REAL(ylength); + PARSE_INT(imax); + PARSE_INT(jmax); + PARSE_INT(itermax); + PARSE_REAL(eps); + PARSE_REAL(omg); + PARSE_REAL(rho); + } + } + + fclose(fp); +} + +void printParameter(Parameter* param) +{ + printf("Parameters:\n"); + printf("Geometry data:\n"); + printf("\tDomain box size (x, y): %e, %e\n", param->xlength, param->ylength); + printf("\tCells (x, y): %d, %d\n", param->imax, param->jmax); + printf("Iterative solver parameters:\n"); + printf("\tMax iterations: %d\n", param->itermax); + printf("\tepsilon (stopping tolerance) : %e\n", param->eps); + printf("\tomega (SOR relaxation): %e\n", param->omg); +} diff --git a/PoissonSolver/2D-mpi-sq/src/parameter.h b/PoissonSolver/2D-mpi-sq/src/parameter.h new file mode 100644 index 0000000..1f8cdcc --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/parameter.h @@ -0,0 +1,18 @@ +/* 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 __PARAMETER_H_ +#define __PARAMETER_H_ + +typedef struct { + double xlength, ylength; + int imax, jmax; + int itermax; + double eps, omg, rho, gamma; +} Parameter; + +void initParameter(Parameter*); +void readParameter(Parameter*, const char*); +void printParameter(Parameter*); +#endif diff --git a/PoissonSolver/2D-mpi-sq/src/solver.c b/PoissonSolver/2D-mpi-sq/src/solver.c new file mode 100644 index 0000000..3346091 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/solver.c @@ -0,0 +1,222 @@ +/* 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. */ +#include "math.h" +#include "stdio.h" +#include "stdlib.h" + +#include "allocate.h" +#include "parameter.h" +#include "solver.h" +#include + +#define PI 3.14159265358979323846 +#define P(i, j) p[(j) * (solver->limax + 2) + (i)] +#define RHS(i, j) rhs[(j) * (solver->limax + 2) + (i)] + +int distribute_dim(int dim_rank, int dim_comm_size, int dim_size) +{ + int base_size = dim_size / dim_comm_size; + return dim_rank < (dim_size % dim_comm_size) ? base_size + 1 : base_size; +} + +int get_dim_start(int dim_rank, int dim_comm_size, int dim_size) +{ + int dim_start = 0; + for (int i = 0; i < dim_rank; i++) + dim_start += distribute_dim(i, dim_comm_size, dim_size); + return dim_start; +} + +void initSolver(Solver* solver, Parameter* params, int problem) +{ + solver->imax = params->imax; + solver->jmax = params->jmax; + solver->dx = params->xlength / params->imax; + solver->dy = params->ylength / params->jmax; + solver->eps = params->eps; + solver->omega = params->omg; + solver->rho = params->rho; + solver->itermax = params->itermax; + + int imax = solver->imax; + int jmax = solver->jmax; + + int rank, size; + solver->dims[0] = 0; + solver->dims[1] = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + MPI_Dims_create(size, 2, solver->dims); + MPI_Cart_create(MPI_COMM_WORLD, + 2, + solver->dims, + (int[2]) { 0, 0 }, + 0, + &solver->cart_comm); + MPI_Cart_coords(solver->cart_comm, rank, 2, solver->coords); + solver->limax = distribute_dim(solver->coords[0], solver->dims[0], solver->imax); + solver->ljmax = distribute_dim(solver->coords[1], solver->dims[1], solver->jmax); + + solver->snd_displacements[0] = ((solver->limax + 2) + 1) * sizeof(double); + solver->snd_displacements[1] = ((solver->limax + 2) * (solver->ljmax) + 1) * + sizeof(double); + solver->snd_displacements[2] = ((solver->limax + 2) + 1) * sizeof(double); + solver->snd_displacements[3] = ((solver->limax + 2) + solver->limax) * sizeof(double); + + solver->rcv_displacements[0] = (1) * sizeof(double); + solver->rcv_displacements[1] = ((solver->limax + 2) * (solver->ljmax + 1) + 1) * + sizeof(double); + solver->rcv_displacements[2] = ((solver->limax + 2)) * sizeof(double); + solver->rcv_displacements[3] = ((solver->limax + 2) + solver->limax + 1) * + sizeof(double); + + MPI_Datatype contiguos_boundary; + MPI_Datatype vector_boundary; + MPI_Datatype internal; + + MPI_Type_contiguous(solver->limax, MPI_DOUBLE, &contiguos_boundary); + MPI_Type_vector(solver->ljmax, 1, solver->limax + 2, MPI_DOUBLE, &vector_boundary); + MPI_Type_create_subarray(2, + (int[2]){solver->limax+2,solver->ljmax+2}, + (int[2]){solver->limax,solver->ljmax}, + (int[2]){1,1}, + MPI_ORDER_C, + MPI_DOUBLE, + &internal); + MPI_Type_commit(&internal); + + MPI_Type_commit(&contiguos_boundary); + MPI_Type_commit(&vector_boundary); + solver->types[0] = contiguos_boundary; + solver->types[1] = contiguos_boundary; + solver->types[2] = vector_boundary; + solver->types[3] = vector_boundary; + solver->types[4] = internal; + + size_t bytesize = (solver->limax + 2) * (solver->ljmax + 2) * sizeof(double); + solver->p = allocate(64, bytesize); + solver->rhs = allocate(64, bytesize); + double dx = solver->dx; + double dy = solver->dy; + double* p = solver->p; + double* rhs = solver->rhs; + int i_start = get_dim_start(solver->coords[0], solver->dims[0], solver->imax); + int j_start = get_dim_start(solver->coords[1], solver->dims[1], solver->jmax); + printf("Prc %d; {%d,%d}, {%d,%d}, {%d,%d}\n", + rank, + solver->coords[0], + solver->coords[1], + solver->limax, + solver->ljmax, + i_start, + j_start); + fflush(stdout); + + for (int j = 0; j < solver->ljmax + 2; j++) { + for (int i = 0; i < solver->limax + 2; i++) { + P(i, j) = sin(2.0 * PI * (i_start + i) * dx * 2.0) + + sin(2.0 * PI * (j_start + j) * dy * 2.0); + } + } + + if (problem == 2) { + for (int j = 0; j < solver->ljmax + 2; j++) { + for (int i = 0; i < solver->limax + 2; i++) { + RHS(i, j) = sin(2.0 * PI * (i_start + i) * dx); + } + } + } else { + for (int j = 0; j < solver->ljmax + 2; j++) { + for (int i = 0; i < solver->limax + 2; i++) { + RHS(i, j) = 0.0; + } + } + } +} + +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++; + } + + // printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); + printf("%d, %f\n", it, solver->omega); +} + +void writeResult(Solver* solver, char* filename) +{ + int imax = solver->imax; + int jmax = solver->jmax; + double* p = solver->p; + + FILE* fp; + fp = fopen(filename, "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + for (int j = 0; j < jmax + 2; j++) { + for (int i = 0; i < imax + 2; i++) { + fprintf(fp, "%f ", P(i, j)); + } + fprintf(fp, "\n"); + } + + fclose(fp); +} diff --git a/PoissonSolver/2D-mpi-sq/src/solver.h b/PoissonSolver/2D-mpi-sq/src/solver.h new file mode 100644 index 0000000..01961b3 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/solver.h @@ -0,0 +1,28 @@ +/* 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 __SOLVER_H_ +#define __SOLVER_H_ +#include "parameter.h" +#include + +typedef struct { + double dx, dy; + int imax, jmax; + int limax,ljmax; + double *p, *rhs; + double eps, omega, rho; + int itermax; + MPI_Comm cart_comm; + MPI_Datatype types[5]; + MPI_Aint rcv_displacements[4]; + MPI_Aint snd_displacements[4]; + int coords[2]; + int dims[2]; +} Solver; + +extern void initSolver(Solver*, Parameter*, int problem); +extern void writeResult(Solver*, char*); +extern void solveRB(Solver*); +#endif diff --git a/PoissonSolver/2D-mpi-sq/src/timing.c b/PoissonSolver/2D-mpi-sq/src/timing.c new file mode 100644 index 0000000..09620e1 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/timing.c @@ -0,0 +1,22 @@ +/* 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. */ +#include +#include + +double getTimeStamp() +{ + struct timespec ts; + clock_gettime(CLOCK_MONOTONIC, &ts); + return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9; +} + +double getTimeResolution() +{ + struct timespec ts; + clock_getres(CLOCK_MONOTONIC, &ts); + return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9; +} + +double getTimeStamp_() { return getTimeStamp(); } diff --git a/PoissonSolver/2D-mpi-sq/src/timing.h b/PoissonSolver/2D-mpi-sq/src/timing.h new file mode 100644 index 0000000..ecfe8a9 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/timing.h @@ -0,0 +1,11 @@ +/* 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 __TIMING_H_ +#define __TIMING_H_ + +extern double getTimeStamp(); +extern double getTimeResolution(); + +#endif // __TIMING_H_ diff --git a/PoissonSolver/2D-mpi-sq/src/util.h b/PoissonSolver/2D-mpi-sq/src/util.h new file mode 100644 index 0000000..3fa11a6 --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/src/util.h @@ -0,0 +1,20 @@ +/* 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 __UTIL_H_ +#define __UTIL_H_ +#define HLINE \ + "----------------------------------------------------------------------------\n" + +#ifndef MIN +#define MIN(x, y) ((x) < (y) ? (x) : (y)) +#endif +#ifndef MAX +#define MAX(x, y) ((x) > (y) ? (x) : (y)) +#endif +#ifndef ABS +#define ABS(a) ((a) >= 0 ? (a) : -(a)) +#endif + +#endif // __UTIL_H_ diff --git a/PoissonSolver/2D-mpi-sq/surface.plot b/PoissonSolver/2D-mpi-sq/surface.plot new file mode 100644 index 0000000..7d56c7d --- /dev/null +++ b/PoissonSolver/2D-mpi-sq/surface.plot @@ -0,0 +1,7 @@ +set terminal png size 1024,768 enhanced font ,12 +set output 'p.png' +set datafile separator whitespace + +set grid +set hidden3d +splot 'p.dat' matrix using 1:2:3 with lines