moebiusband/NuSiF-Solver
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Requested changes

Browse Source
This commit is contained in:
Aditya Ujeniya 2024-09-16 13:51:10 +02:00
parent be3db085af
commit 1a1420890b
8 changed files with 75 additions and 190 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
BasicSolver/2D-mpi/config.mk
View File

@ -1,5 +1,5 @@
# Supported: GCC, CLANG, ICC
TAG ?= ICC
# Supported: GCC, CLANG, ICX
TAG ?= ICX
# Supported: true, false
ENABLE_MPI ?= true
ENABLE_OPENMP ?= false

6
BasicSolver/2D-mpi/include_ICC.mk → BasicSolver/2D-mpi/include_ICX.mk
View File

@ -1,8 +1,8 @@
ifeq ($(ENABLE_MPI),true)
CC = mpiicc
CC = mpiicx
DEFINES = -D_MPI
else
CC = icc
CC = icx
endif
GCC = gcc
@ -13,7 +13,7 @@ OPENMP = -qopenmp
endif
VERSION = --version
CFLAGS = -O3 -xHost -qopt-zmm-usage=high -std=c99 $(OPENMP)
CFLAGS = -O3 -xHost -qopt-zmm-usage=high -std=c99 $(OPENMP) -Wno-unused-command-line-argument
LFLAGS = $(OPENMP)
DEFINES += -D_GNU_SOURCE# -DDEBUG
INCLUDES =

105
BasicSolver/2D-mpi/src/comm-v1.c
View File

@ -10,18 +10,6 @@
#include "comm.h"
#ifdef _MPI
// subroutines local to this module
static int sum(int* sizes, int position)
{
int sum = 0;
for (int i = 0; i < position; i += position) {
sum += sizes[i];
}
return sum;
}
static void gatherArray(
Comm* c, int cnt, int* rcvCounts, int* displs, double* src, double* dst)
{
@ -42,6 +30,52 @@ static void gatherArray(
MPI_COMM_WORLD);
}
#endif // defined _MPI
void commCollectResult(Comm* c,
double* ug,
double* vg,
double* pg,
double* u,
double* v,
double* p,
int imax,
int jmax)
{
#ifdef _MPI
int *rcvCounts, *displs;
int cnt = c->jmaxLocal * (imax + 2);
if (c->rank == 0) {
rcvCounts = (int*)malloc(c->size * sizeof(int));
displs = (int*)malloc(c->size * sizeof(int));
for (int i = 0; i < c->size; ++i) {
rcvCounts[i] = 0;
displs[i] = 0;
}
}
if (c->rank == 0 && c->size == 1) {
cnt = (c->jmaxLocal + 2) * (imax + 2);
} else if (c->rank == 0 || c->rank == (c->size - 1)) {
cnt = (c->jmaxLocal + 1) * (imax + 2);
}
MPI_Gather(&cnt, 1, MPI_INTEGER, rcvCounts, 1, MPI_INTEGER, 0, MPI_COMM_WORLD);
if (c->rank == 0) {
displs[0] = 0;
int cursor = rcvCounts[0];
for (int i = 1; i < c->size; i++) {
displs[i] = cursor;
cursor += rcvCounts[i];
}
}
gatherArray(c, cnt, rcvCounts, displs, p, pg);
gatherArray(c, cnt, rcvCounts, displs, u, ug);
gatherArray(c, cnt, rcvCounts, displs, v, vg);
#endif
}
// exported subroutines
int commIsBoundary(Comm* c, int direction)
@ -128,49 +162,6 @@ void commShift(Comm* c, double* f, double* g)
#endif
}
void commCollectResult(Comm* c,
double* ug,
double* vg,
double* pg,
double* u,
double* v,
double* p,
int jmax,
int imax)
{
#ifdef _MPI
int *rcvCounts, *displs;
int cnt = c->jmaxLocal * (imax + 2);
if (c->rank == 0) {
rcvCounts = (int*)malloc(c->size * sizeof(int));
displs = (int*)malloc(c->size * sizeof(int));
}
if (c->rank == 0 && c->size == 1) {
cnt = (c->jmaxLocal + 2) * (imax + 2);
} else if (c->rank == 0 || c->rank == (c->size - 1)) {
cnt = (c->jmaxLocal + 1) * (imax + 2);
}
MPI_Gather(&cnt, 1, MPI_INTEGER, rcvCounts, 1, MPI_INTEGER, 0, MPI_COMM_WORLD);
if (c->rank == 0) {
displs[0] = 0;
int cursor = rcvCounts[0];
for (int i = 1; i < c->size; i++) {
displs[i] = cursor;
cursor += rcvCounts[i];
}
}
gatherArray(c, cnt, rcvCounts, displs, p, pg);
gatherArray(c, cnt, rcvCounts, displs, u, ug);
gatherArray(c, cnt, rcvCounts, displs, v, vg);
#endif
}
void commPartition(Comm* c, int jmax, int imax)
{
#ifdef _MPI
@ -211,7 +202,8 @@ void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLo
newcomm->jmaxLocal = jmaxLocal / 2;
newcomm->neighbours[BOTTOM] = newcomm->rank == 0 ? -1 : newcomm->rank - 1;
newcomm->neighbours[TOP] = newcomm->rank == (newcomm->size - 1) ? -1 : newcomm->rank + 1;
newcomm->neighbours[TOP] = newcomm->rank == (newcomm->size - 1) ? -1
: newcomm->rank + 1;
newcomm->neighbours[LEFT] = -1;
newcomm->neighbours[RIGHT] = -1;
@ -221,7 +213,6 @@ void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLo
newcomm->dims[IDIM] = 1;
newcomm->dims[JDIM] = newcomm->size;
#endif
newcomm->imaxLocal = imaxLocal;
newcomm->jmaxLocal = jmaxLocal;

18
BasicSolver/2D-mpi/src/main.c
View File

@ -16,6 +16,24 @@
#include "solver.h"
#include "timing.h"
static FILE* initResidualWriter()
{
FILE* fp;
fp = fopen("residual.dat", "w");
if (fp == NULL) {
printf("Error!\n");
exit(EXIT_FAILURE);
}
return fp;
}
static void writeResidual(FILE* fp, double ts, double res)
{
fprintf(fp, "%f, %f\n", ts, res);
}
static void writeResults(Discretization* s)
{
#ifdef _MPI

19
BasicSolver/2D-mpi/src/progress.c
View File

@ -49,22 +49,3 @@ void stopProgress()
printf("\n");
fflush(stdout);
}
FILE* initResidualWriter()
{
FILE* fp;
fp = fopen("residual.dat", "w");
if (fp == NULL) {
printf("Error!\n");
exit(EXIT_FAILURE);
}
return fp;
}
void writeResidual(FILE* fp, double ts, double res)
{
fprintf(fp, "%f, %f\n", ts, res);
}

2
BasicSolver/2D-mpi/src/progress.h
View File

@ -10,6 +10,4 @@
extern void initProgress(double);
extern void printProgress(double);
extern void stopProgress();
extern FILE* initResidualWriter(void);
extern void writeResidual(FILE*, double, double);
#endif

5
BasicSolver/2D-mpi/src/solver-mg.c
View File

@ -109,7 +109,7 @@ static void setBoundaryCondition(Solver* s, double* p, int imaxLocal, int jmaxLo
#endif
}
static double smooth(Solver* s, double* p, double* rhs, int level, Comm* comm)
static void smooth(Solver* s, double* p, double* rhs, int level, Comm* comm)
{
int imaxLocal = comm->imaxLocal;
int jmaxLocal = comm->jmaxLocal;
@ -228,10 +228,11 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, Comm* comm
// restrict
restrictMG(s, level, comm);
// Create a new comm object withupdated imaxLocal and jmaxLocal
// along with their updated bufferTypes, sdispls, rdispls
Comm newcomm;
commUpdateDatatypes(s->comm, &newcomm, comm->imaxLocal, comm->jmaxLocal);
// MGSolver on residual and error.
multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, &newcomm);

104
BasicSolver/2D-mpi/src/solver.c
View File

@ -1,104 +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 <math.h>
#include <stdio.h>
#include "allocate.h"
#include "comm.h"
#include "discretization.h"
#include "parameter.h"
#include "solver.h"
#include "util.h"
void initSolver(Solver* s, Discretization* d, Parameter* p)
{
s->grid = &d->grid;
s->eps = p->eps;
s->omega = p->omg;
s->itermax = p->itermax;
s->comm = &d->comm;
}
void solve(Solver* s, double* p, double* rhs)
{
int imax = s->grid->imax;
int jmax = s->grid->jmax;
int imaxLocal = s->comm->imaxLocal;
int jmaxLocal = s->comm->jmaxLocal;
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 idx2 = 1.0 / dx2;
double idy2 = 1.0 / dy2;
double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
double epssq = eps * eps;
int pass, jsw, isw;
int it = 0;
double res = 1.0;
while ((res >= epssq) && (it < itermax)) {
jsw = 1;
for (pass = 0; pass < 2; pass++) {
isw = jsw;
commExchange(s->comm, p);
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = isw; i < imaxLocal + 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;
}
if (commIsBoundary(s->comm, BOTTOM)) { // set bottom bc
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, 0) = P(i, 1);
}
}
if (commIsBoundary(s->comm, TOP)) { // set top bc
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, jmaxLocal + 1) = P(i, jmaxLocal);
}
}
if (commIsBoundary(s->comm, LEFT)) { // set left bc
for (int j = 1; j < jmaxLocal + 1; j++) {
P(0, j) = P(1, j);
}
}
if (commIsBoundary(s->comm, RIGHT)) { // set right bc
for (int j = 1; j < jmaxLocal + 1; j++) {
P(imaxLocal + 1, j) = P(imaxLocal, j);
}
}
commReduction(&res, SUM);
res = res / (double)(imax * jmax);
#ifdef DEBUG
if (commIsMaster(s->comm)) {
printf("%d Residuum: %e\n", it, res);
}
#endif
it++;
}
#ifdef VERBOSE
if (commIsMaster(s->comm)) {
printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
}
#endif
}