ErikFabrizzi/NuSiF-Solver
1
0
Fork 0
forked from moebiusband/NuSiF-Solver
Code Pull Requests Activity

Branch Checkpoint

Browse Source
This commit is contained in:
indiano
2025-05-01 21:48:47 +02:00
parent 765883ab88
commit 5148069114
20 changed files with 755 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

37
PoissonSolver/2D-mpi-sq/src/allocate.c Normal file
View File

@@ -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 <stdlib.h>
#include <stdio.h>
#include <errno.h>
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;
}

11
PoissonSolver/2D-mpi-sq/src/allocate.h Normal file
View File

@@ -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 <stdlib.h>
extern void* allocate(int alignment, size_t bytesize);
#endif

53
PoissonSolver/2D-mpi-sq/src/likwid-marker.h Normal file
View File

@@ -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 <likwid.h>
#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*/

46
PoissonSolver/2D-mpi-sq/src/main.c Normal file
View File

@@ -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 <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#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(&params);
if (argc < 2) {
printf("Usage: %s <configFile>\n", argv[0]);
exit(EXIT_SUCCESS);
}
readParameter(&params, argv[1]);
initSolver(&solver, &params, 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;
}

79
PoissonSolver/2D-mpi-sq/src/parameter.c Normal file
View File

@@ -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 <stdio.h>
#include <stdlib.h>
#include <string.h>
//---
#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);
}

18
PoissonSolver/2D-mpi-sq/src/parameter.h Normal file
View File

@@ -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

222
PoissonSolver/2D-mpi-sq/src/solver.c Normal file
View File

@@ -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 <mpi.h>
#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);
}

28
PoissonSolver/2D-mpi-sq/src/solver.h Normal file
View File

@@ -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 <mpi.h>
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

22
PoissonSolver/2D-mpi-sq/src/timing.c Normal file
View File

@@ -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 <stdlib.h>
#include <time.h>
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(); }

11
PoissonSolver/2D-mpi-sq/src/timing.h Normal file
View File

@@ -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_

20
PoissonSolver/2D-mpi-sq/src/util.h Normal file
View File

@@ -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_