This commit is contained in:
Jan Eitzinger 2024-01-31 18:02:55 +01:00
parent cc2aef66b4
commit 57193243f3
18 changed files with 57 additions and 237 deletions

View File

@ -1,90 +0,0 @@
/*
* Copyright (C) 2024 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 <float.h>
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include "allocate.h"
#include "comm.h"
#include "parameter.h"
#include "progress.h"
#include "solver.h"
#include "timing.h"
#include "vtkWriter.h"
int main(int argc, char** argv)
{
double timeStart, timeStop;
Parameter p;
Solver s;
commInit(&s.comm, argc, argv);
initParameter(&p);
if (argc != 2) {
printf("Usage: %s <configFile>\n", argv[0]);
exit(EXIT_SUCCESS);
}
readParameter(&p, argv[1]);
commPartition(&s.comm, p.kmax, p.jmax, p.imax);
if (commIsMaster(&s.comm)) {
printParameter(&p);
}
initSolver(&s, &p);
#ifndef VERBOSE
initProgress(s.te);
#endif
double tau = s.tau;
double te = s.te;
double t = 0.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;
#ifdef VERBOSE
if (commIsMaster(&s.comm)) {
printf("TIME %f , TIMESTEP %f\n", t, s.dt);
}
#else
printProgress(t);
#endif
}
timeStop = getTimeStamp();
#ifndef VERBOSE
stopProgress();
#endif
if (commIsMaster(&s.comm)) {
printf("Solution took %.2fs\n", timeStop - timeStart);
}
timeStart = getTimeStamp();
VtkOptions opts = { .grid = s.grid, .comm = s.comm };
vtkOpen(&opts, s.problem);
vtkScalar(&opts, "pressure", s.p);
vtkVector(&opts, "velocity", (VtkVector) { s.u, s.v, s.w });
vtkClose(&opts);
timeStop = getTimeStamp();
if (commIsMaster(&s.comm)) {
printf("Result output took %.2fs\n", timeStop - timeStart);
}
commFinalize(&s.comm);
return EXIT_SUCCESS;
}

View File

@ -1,5 +1,5 @@
#======================================================================================= #=======================================================================================
# Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. # Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
# All rights reserved. # All rights reserved.
# Use of this source code is governed by a MIT-style # Use of this source code is governed by a MIT-style
# license that can be found in the LICENSE file. # license that can be found in the LICENSE file.

View File

@ -10,7 +10,7 @@ endif
VERSION = --version VERSION = --version
# CFLAGS = -O3 -std=c17 $(OPENMP) # CFLAGS = -O3 -std=c17 $(OPENMP)
CFLAGS = -Ofast -std=c17 -Weverything CFLAGS = -Ofast -std=c17 #-Weverything
#CFLAGS = -Ofast -fnt-store=aggressive -std=c99 $(OPENMP) #AMD CLANG #CFLAGS = -Ofast -fnt-store=aggressive -std=c99 $(OPENMP) #AMD CLANG
LFLAGS = $(OPENMP) -lm LFLAGS = $(OPENMP) -lm
DEFINES = -D_GNU_SOURCE# -DDEBUG DEFINES = -D_GNU_SOURCE# -DDEBUG

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. * All rights reserved.
* Use of this source code is governed by a MIT-style * Use of this source code is governed by a MIT-style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. * All rights reserved.
* Use of this source code is governed by a MIT-style * Use of this source code is governed by a MIT-style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver. * All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style * Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,11 +1,9 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. * All rights reserved.
* Use of this source code is governed by a MIT-style * Use of this source code is governed by a MIT-style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.
*/ */
#include <float.h>
#include <limits.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <unistd.h> #include <unistd.h>
@ -95,7 +93,7 @@ int main(int argc, char** argv)
setSpecialBoundaryCondition(&s); setSpecialBoundaryCondition(&s);
computeFG(&s); computeFG(&s);
computeRHS(&s); computeRHS(&s);
solveRB(&s); solve(&s);
adaptUV(&s); adaptUV(&s);
t += s.dt; t += s.dt;
@ -111,7 +109,7 @@ int main(int argc, char** argv)
#endif #endif
printf("Solution took %.2fs\n", timeStop - timeStart); printf("Solution took %.2fs\n", timeStop - timeStart);
timeStart = getTimeStamp(); timeStart = getTimeStamp();
double *pg, *ug, *vg, *wg; 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)(s.grid.imax * s.grid.jmax * s.grid.kmax) * sizeof(double);
@ -127,11 +125,9 @@ int main(int argc, char** argv)
vtkScalar(&opts, "pressure", pg); vtkScalar(&opts, "pressure", pg);
vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg }); vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
vtkClose(&opts); vtkClose(&opts);
timeStop = getTimeStamp();
if (commIsMaster(&s.comm)) { timeStop = getTimeStamp();
printf("Result output took %.2fs\n", timeStop - timeStart); printf("Result output took %.2fs\n", timeStop - timeStart);
}
return EXIT_SUCCESS; return EXIT_SUCCESS;
} }

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver. * All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style * Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver. * All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style * Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver. * All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style * Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. * All rights reserved.
* Use of this source code is governed by a MIT-style * Use of this source code is governed by a MIT-style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver. * All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style * Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.
@ -7,7 +7,6 @@
#include <float.h> #include <float.h>
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h>
#include <string.h> #include <string.h>
#include "allocate.h" #include "allocate.h"
@ -41,19 +40,13 @@ static void printConfig(Solver* s)
s->grid.xlength, s->grid.xlength,
s->grid.ylength, s->grid.ylength,
s->grid.zlength); s->grid.zlength);
printf("\tCells (x, y, z): %d, %d, %d\n", printf("\tCells (x, y, z): %d, %d, %d\n", s->grid.imax, s->grid.jmax, s->grid.kmax);
s->grid.imax, printf("\tCell size (dx, dy, dz): %f, %f, %f\n", s->grid.dx, s->grid.dy, s->grid.dz);
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);
printf("Timestep parameters:\n"); printf("Timestep parameters:\n");
printf("\tDefault stepsize: %.2f, Final time %.2f\n", s->dt, s->te); printf("\tDefault stepsize: %.2f, Final time %.2f\n", s->dt, s->te);
printf("\tdt bound: %.6f\n", s->dtBound); printf("\tdt bound: %.6f\n", s->dtBound);
printf("\tTau factor: %.2f\n", s->tau); printf("\tTau factor: %.2f\n", s->tau);
printf("Iterative parameters:\n"); printf("Iterative parameters:\n");
printf("\tMax iterations: %d\n", s->itermax); printf("\tMax iterations: %d\n", s->itermax);
printf("\tepsilon (stopping tolerance) : %f\n", s->eps); printf("\tepsilon (stopping tolerance) : %f\n", s->eps);
printf("\tgamma factor: %f\n", s->gamma); printf("\tgamma factor: %f\n", s->gamma);
@ -62,13 +55,13 @@ static void printConfig(Solver* s)
void initSolver(Solver* s, Parameter* params) void initSolver(Solver* s, Parameter* params)
{ {
s->problem = params->name; s->problem = params->name;
s->bcLeft = params->bcLeft; s->bcLeft = params->bcLeft;
s->bcRight = params->bcRight; s->bcRight = params->bcRight;
s->bcBottom = params->bcBottom; s->bcBottom = params->bcBottom;
s->bcTop = params->bcTop; s->bcTop = params->bcTop;
s->bcFront = params->bcFront; s->bcFront = params->bcFront;
s->bcBack = params->bcBack; s->bcBack = params->bcBack;
s->grid.imax = params->imax; s->grid.imax = params->imax;
s->grid.jmax = params->jmax; s->grid.jmax = params->jmax;
@ -80,17 +73,17 @@ void initSolver(Solver* s, Parameter* params)
s->grid.dy = params->ylength / params->jmax; s->grid.dy = params->ylength / params->jmax;
s->grid.dz = params->zlength / params->kmax; s->grid.dz = params->zlength / params->kmax;
s->eps = params->eps; s->eps = params->eps;
s->omega = params->omg; s->omega = params->omg;
s->itermax = params->itermax; s->itermax = params->itermax;
s->re = params->re; s->re = params->re;
s->gx = params->gx; s->gx = params->gx;
s->gy = params->gy; s->gy = params->gy;
s->gz = params->gz; s->gz = params->gz;
s->dt = params->dt; s->dt = params->dt;
s->te = params->te; s->te = params->te;
s->tau = params->tau; s->tau = params->tau;
s->gamma = params->gamma; s->gamma = params->gamma;
int imax = s->grid.imax; int imax = s->grid.imax;
int jmax = s->grid.jmax; int jmax = s->grid.jmax;
@ -130,13 +123,13 @@ void initSolver(Solver* s, Parameter* params)
void computeRHS(Solver* s) void computeRHS(Solver* s)
{ {
int imax = s->grid.imax; int imax = s->grid.imax;
int jmax = s->grid.jmax; int jmax = s->grid.jmax;
int kmax = s->grid.kmax; int kmax = s->grid.kmax;
double idx = 1.0 / s->grid.dx; double idx = 1.0 / s->grid.dx;
double idy = 1.0 / s->grid.dy; double idy = 1.0 / s->grid.dy;
double idz = 1.0 / s->grid.dz; double idz = 1.0 / s->grid.dz;
double idt = 1.0 / s->dt; double idt = 1.0 / s->dt;
double* rhs = s->rhs; double* rhs = s->rhs;
double* f = s->f; double* f = s->f;
@ -156,83 +149,6 @@ void computeRHS(Solver* s)
} }
void solve(Solver* s) void solve(Solver* s)
{
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;
while ((res >= epssq) && (it < itermax)) {
res = 0.0;
for (int k = 1; k < kmax + 1; k++) {
for (int j = 1; j < jmax + 1; j++) {
for (int i = 1; i < imax + 1; i++) {
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);
}
}
}
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
}
void solveRB(Solver* s)
{ {
int imax = s->grid.imax; int imax = s->grid.imax;
int jmax = s->grid.jmax; int jmax = s->grid.jmax;
@ -347,10 +263,10 @@ void normalizePressure(Solver* s)
void computeTimestep(Solver* s) void computeTimestep(Solver* s)
{ {
double dt = s->dtBound; double dt = s->dtBound;
double dx = s->grid.dx; double dx = s->grid.dx;
double dy = s->grid.dy; double dy = s->grid.dy;
double dz = s->grid.dz; double dz = s->grid.dz;
double umax = maxElement(s, s->u); double umax = maxElement(s, s->u);
double vmax = maxElement(s, s->v); double vmax = maxElement(s, s->v);
@ -613,10 +529,10 @@ void computeFG(Solver* s)
double* g = s->g; double* g = s->g;
double* h = s->h; double* h = s->h;
double gx = s->gx; double gx = s->gx;
double gy = s->gy; double gy = s->gy;
double gz = s->gz; double gz = s->gz;
double dt = s->dt; double dt = s->dt;
double gamma = s->gamma; double gamma = s->gamma;
double inverseRe = 1.0 / s->re; double inverseRe = 1.0 / s->re;

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver. * All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style * Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.
@ -34,12 +34,10 @@ typedef struct {
extern void initSolver(Solver*, Parameter*); extern void initSolver(Solver*, Parameter*);
extern void computeRHS(Solver*); extern void computeRHS(Solver*);
extern void solve(Solver*); extern void solve(Solver*);
extern void solveRB(Solver*);
extern void normalizePressure(Solver*); extern void normalizePressure(Solver*);
extern void computeTimestep(Solver*); extern void computeTimestep(Solver*);
extern void setBoundaryConditions(Solver*); extern void setBoundaryConditions(Solver*);
extern void setSpecialBoundaryCondition(Solver*); extern void setSpecialBoundaryCondition(Solver*);
extern void computeFG(Solver*); extern void computeFG(Solver*);
extern void adaptUV(Solver*); extern void adaptUV(Solver*);
extern void writeResult(Solver*);
#endif #endif

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. * All rights reserved.
* Use of this source code is governed by a MIT-style * Use of this source code is governed by a MIT-style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. * All rights reserved.
* Use of this source code is governed by a MIT-style * Use of this source code is governed by a MIT-style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. * All rights reserved.
* Use of this source code is governed by a MIT-style * Use of this source code is governed by a MIT-style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver. * All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style * Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.

View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver. * All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style * Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file. * license that can be found in the LICENSE file.