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

Implemented Function pointer for SOR, RB and RBA variants

Browse Source
This commit is contained in:
Aditya Ujeniya
2023-07-14 21:50:38 +02:00
parent 9f55413efb
commit b50dbe7d4b
30 changed files with 61658 additions and 61142 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20050
BasicSolver/2D-mpi-v1/pressure.dat
View File

File diff suppressed because it is too large Load Diff

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

@@ -22,7 +22,7 @@ enum VARIANT { SOR = 1, RB, RBA };
int main (int argc, char** argv)
{
int rank;
int variant = SOR;
int variant = RB;
double S, E;
Parameter params;
Solver solver;
@@ -47,15 +47,17 @@ int main (int argc, char** argv)
initSolver(&solver, &params);
initProgress(solver.te);
void (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
double tau = solver.tau;
double te = solver.te;
double t = 0.0;
S = getTimeStamp();
switch (variant) {
case SOR:
printf("Plain SOR\n");
while (t <= te) {
while (t <= te)
{
if (tau > 0.0) {
computeTimestep(&solver);
}
@@ -64,49 +66,9 @@ int main (int argc, char** argv)
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
solve(&solver);
(*solver_generic[variant - 1])(&solver);
adaptUV(&solver);
t += solver.dt;
}
break;
case RB:
printf("Red-black SOR\n");
while (t <= te) {
if (tau > 0.0) {
computeTimestep(&solver);
}
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
solveRB(&solver);
adaptUV(&solver);
t += solver.dt;
}
break;
case RBA:
printf("Red-black SOR with acceleration\n");
while (t <= te) {
if (tau > 0.0) {
computeTimestep(&solver);
}
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
solveRBA(&solver);
adaptUV(&solver);
t += solver.dt;
}
break;
}
#ifdef VERBOSE
if (rank == 0) {
printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
@@ -114,6 +76,7 @@ int main (int argc, char** argv)
#else
printProgress(t);
#endif
}
E = getTimeStamp();

20000
BasicSolver/2D-mpi-v1/velocity.dat
View File

File diff suppressed because it is too large Load Diff

BIN
BasicSolver/2D-mpi-v1/velocity.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 19 KiB

After

Width:  |  Height:  |  Size: 206 KiB

19602
BasicSolver/2D-mpi-v2/pressure.dat
View File

File diff suppressed because it is too large Load Diff

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

@@ -21,7 +21,7 @@ enum VARIANT { SOR = 1, RB, RBA };
int main (int argc, char** argv)
{
int rank;
int variant = SOR;
int variant = RB;
double S, E;
Parameter params;
Solver solver;
@@ -52,12 +52,12 @@ int main (int argc, char** argv)
double tau = solver.tau;
double te = solver.te;
double t = 0.0;
int (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
S = getTimeStamp();
switch (variant) {
case SOR:
printf("Plain SOR\n");
while (t <= te) {
while (t <= te) {
if (tau > 0.0) {
computeTimestep(&solver);
}
@@ -66,48 +66,9 @@ int main (int argc, char** argv)
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
solve(&solver);
(*solver_generic[variant - 1])(&solver);
adaptUV(&solver);
t += solver.dt;
}
break;
case RB:
printf("Red-black SOR\n");
while (t <= te) {
if (tau > 0.0) {
computeTimestep(&solver);
}
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
solveRB(&solver);
adaptUV(&solver);
t += solver.dt;
}
break;
case RBA:
printf("Red-black SOR with acceleration\n");
while (t <= te) {
if (tau > 0.0) {
computeTimestep(&solver);
}
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
solveRBA(&solver);
adaptUV(&solver);
t += solver.dt;
}
break;
}
#ifdef VERBOSE
if (rank == 0) {
printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
@@ -115,6 +76,8 @@ int main (int argc, char** argv)
#else
printProgress(t);
#endif
}
E = getTimeStamp();
stopProgress();

19552
BasicSolver/2D-mpi-v2/velocity.dat
View File

File diff suppressed because it is too large Load Diff

BIN
BasicSolver/2D-mpi-v2/velocity.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 19 KiB

After

Width:  |  Height:  |  Size: 184 KiB

19602
BasicSolver/2D-mpi-v3/pressure.dat
View File

File diff suppressed because it is too large Load Diff

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

@@ -22,7 +22,7 @@ enum VARIANT { SOR = 1, RB, RBA };
int main (int argc, char** argv)
{
int rank;
int variant = SOR;
int variant = RB;
double S, E;
Parameter params;
@@ -52,11 +52,14 @@ int main (int argc, char** argv)
double te = solver.te;
double t = 0.0;
int (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
S = getTimeStamp();
switch (variant) {
case SOR:
printf("Plain SOR\n");
while (t <= te) {
while (t <= te)
{
if (tau > 0.0) {
computeTimestep(&solver);
}
@@ -65,47 +68,10 @@ int main (int argc, char** argv)
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
solve(&solver);
(*solver_generic[variant - 1])(&solver);
adaptUV(&solver);
t += solver.dt;
}
break;
case RB:
printf("Red-black SOR\n");
while (t <= te) {
if (tau > 0.0) {
computeTimestep(&solver);
}
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
solveRB(&solver);
adaptUV(&solver);
t += solver.dt;
}
break;
case RBA:
printf("Red-black SOR with acceleration\n");
while (t <= te) {
if (tau > 0.0) {
computeTimestep(&solver);
}
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
solveRBA(&solver);
adaptUV(&solver);
t += solver.dt;
}
break;
}
#ifdef VERBOSE
if (rank == 0) {
@@ -114,6 +80,7 @@ int main (int argc, char** argv)
#else
printProgress(t);
#endif
}
E = getTimeStamp();
stopProgress();

19552
BasicSolver/2D-mpi-v3/velocity.dat
View File

File diff suppressed because it is too large Load Diff

3200
BasicSolver/2D-seq/pressure.dat
View File

File diff suppressed because it is too large Load Diff

91
BasicSolver/2D-seq/src/main.c
View File

@@ -20,7 +20,7 @@ enum VARIANT { SOR = 1, RB, RBA };
int main (int argc, char** argv)
{
int rank;
int variant = SOR;
int variant = RB;
double S, E;
Parameter params;
@@ -49,76 +49,29 @@ int main (int argc, char** argv)
double t = 0.0;
int nt = 0;
void (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
S = getTimeStamp();
switch (variant) {
case SOR:
printf("Plain SOR\n");
while (t <= te)
{
if (tau > 0.0) computeTimestep(&solver);
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
if (nt % 100 == 0) normalizePressure(&solver);
solve(&solver);
adaptUV(&solver);
t += solver.dt;
nt++;
#ifdef VERBOSE
printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
#else
printProgress(t);
#endif
}
break;
case RB:
printf("Red-black SOR\n");
while (t <= te)
{
if (tau > 0.0) computeTimestep(&solver);
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
if (nt % 100 == 0) normalizePressure(&solver);
solveRB(&solver);
adaptUV(&solver);
t += solver.dt;
nt++;
#ifdef VERBOSE
printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
#else
printProgress(t);
#endif
}
break;
case RBA:
printf("Red-black SOR with acceleration\n");
while (t <= te)
{
if (tau > 0.0) computeTimestep(&solver);
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
if (nt % 100 == 0) normalizePressure(&solver);
solveRBA(&solver);
adaptUV(&solver);
t += solver.dt;
nt++;
#ifdef VERBOSE
printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
#else
printProgress(t);
#endif
}
break;
while (t <= te)
{
if (tau > 0.0) computeTimestep(&solver);
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
if (nt % 100 == 0) normalizePressure(&solver);
(*solver_generic[variant - 1])(&solver);
adaptUV(&solver);
t += solver.dt;
nt++;
#ifdef VERBOSE
printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
#else
printProgress(t);
#endif
}
E = getTimeStamp();

36
BasicSolver/2D-seq/velocity.dat
View File

@@ -24,7 +24,7 @@
0.59 0.01 -0.008094 -0.000200 0.008096
0.61 0.01 -0.007642 -0.000251 0.007647
0.64 0.01 -0.007095 -0.000297 0.007101
0.66 0.01 -0.006463 -0.000335 0.006471
0.66 0.01 -0.006463 -0.000335 0.006472
0.69 0.01 -0.005762 -0.000365 0.005774
0.71 0.01 -0.005011 -0.000386 0.005026
0.74 0.01 -0.004229 -0.000396 0.004248
@@ -76,7 +76,7 @@
0.89 0.04 -0.001666 -0.000909 0.001898
0.91 0.04 -0.000758 -0.000580 0.000954
0.94 0.04 -0.000202 -0.000248 0.000320
0.96 0.04 0.000029 0.000027 0.000039
0.96 0.04 0.000029 0.000026 0.000039
0.99 0.04 0.000035 0.000132 0.000136
0.01 0.06 -0.000034 -0.000123 0.000127
0.04 0.06 -0.000273 0.000256 0.000374
@@ -268,7 +268,7 @@
0.69 0.16 -0.061532 -0.022524 0.065525
0.71 0.16 -0.055295 -0.023898 0.060238
0.74 0.16 -0.048728 -0.024670 0.054617
0.76 0.16 -0.041987 -0.024793 0.048760
0.76 0.16 -0.041986 -0.024793 0.048760
0.79 0.16 -0.035235 -0.024238 0.042766
0.81 0.16 -0.028644 -0.022992 0.036731
0.84 0.16 -0.022386 -0.021072 0.030744
@@ -294,7 +294,7 @@
0.34 0.19 -0.075708 0.026077 0.080073
0.36 0.19 -0.081370 0.023110 0.084588
0.39 0.19 -0.086279 0.019648 0.088488
0.41 0.19 -0.090340 0.015772 0.091707
0.41 0.19 -0.090340 0.015773 0.091707
0.44 0.19 -0.093475 0.011565 0.094187
0.46 0.19 -0.095621 0.007112 0.095885
0.49 0.19 -0.096736 0.002500 0.096768
@@ -389,7 +389,7 @@
0.71 0.24 -0.077691 -0.046320 0.090451
0.74 0.24 -0.068991 -0.048114 0.084111
0.76 0.24 -0.059958 -0.048721 0.077257
0.79 0.24 -0.050796 -0.048075 0.069939
0.79 0.24 -0.050796 -0.048075 0.069938
0.81 0.24 -0.041727 -0.046141 0.062210
0.84 0.24 -0.032981 -0.042925 0.054132
0.86 0.24 -0.024797 -0.038478 0.045776
@@ -408,7 +408,7 @@
0.19 0.26 -0.045674 0.054903 0.071418
0.21 0.26 -0.055333 0.056832 0.079320
0.24 0.26 -0.065054 0.057332 0.086712
0.26 0.26 -0.074620 0.056461 0.093574
0.26 0.26 -0.074621 0.056461 0.093574
0.29 0.26 -0.083835 0.054302 0.099885
0.31 0.26 -0.092517 0.050954 0.105620
0.34 0.26 -0.100504 0.046530 0.110752
@@ -606,7 +606,7 @@
0.14 0.39 -0.040517 0.098443 0.106455
0.16 0.39 -0.053081 0.106833 0.119293
0.19 0.39 -0.066194 0.112215 0.130284
0.21 0.39 -0.079501 0.114706 0.139563
0.21 0.39 -0.079502 0.114706 0.139563
0.24 0.39 -0.092687 0.114466 0.147287
0.26 0.39 -0.105476 0.111690 0.153623
0.29 0.39 -0.117633 0.106594 0.158744
@@ -630,7 +630,7 @@
0.74 0.39 -0.111319 -0.114027 0.159355
0.76 0.39 -0.098068 -0.117318 0.152908
0.79 0.39 -0.084249 -0.117805 0.144830
0.81 0.39 -0.070186 -0.115275 0.134960
0.81 0.39 -0.070185 -0.115275 0.134960
0.84 0.39 -0.056251 -0.109576 0.123171
0.86 0.39 -0.042866 -0.100636 0.109385
0.89 0.39 -0.030488 -0.088477 0.093583
@@ -696,7 +696,7 @@
0.39 0.44 -0.170816 0.082592 0.189735
0.41 0.44 -0.177729 0.066561 0.189784
0.44 0.44 -0.183191 0.049483 0.189756
0.46 0.44 -0.187155 0.031606 0.189805
0.46 0.44 -0.187155 0.031605 0.189805
0.49 0.44 -0.189584 0.013174 0.190041
0.51 0.44 -0.190442 -0.005565 0.190523
0.54 0.44 -0.189701 -0.024360 0.191258
@@ -763,7 +763,7 @@
0.06 0.49 -0.013496 0.088857 0.089876
0.09 0.49 -0.024333 0.115469 0.118005
0.11 0.49 -0.037314 0.136818 0.141815
0.14 0.49 -0.051761 0.152957 0.161477
0.14 0.49 -0.051761 0.152956 0.161477
0.16 0.49 -0.067064 0.164072 0.177249
0.19 0.49 -0.082692 0.170453 0.189452
0.21 0.49 -0.098201 0.172458 0.198457
@@ -797,7 +797,7 @@
0.91 0.49 -0.027263 -0.120196 0.123249
0.94 0.49 -0.015090 -0.091937 0.093167
0.96 0.49 -0.006161 -0.058588 0.058911
0.99 0.49 -0.001306 -0.020753 0.020794
0.99 0.49 -0.001307 -0.020753 0.020794
0.01 0.51 -0.001267 0.022927 0.022962
0.04 0.51 -0.005931 0.063665 0.063941
0.06 0.51 -0.014434 0.098757 0.099806
@@ -818,7 +818,7 @@
0.44 0.51 -0.197267 0.063695 0.207295
0.46 0.51 -0.201225 0.041450 0.205450
0.49 0.51 -0.203822 0.018572 0.204666
0.51 0.51 -0.205038 -0.004702 0.205092
0.51 0.51 -0.205038 -0.004703 0.205092
0.54 0.51 -0.204841 -0.028133 0.206764
0.56 0.51 -0.203186 -0.051464 0.209603
0.59 0.51 -0.200018 -0.074415 0.213412
@@ -930,7 +930,7 @@
0.24 0.59 -0.127806 0.228465 0.261784
0.26 0.59 -0.140687 0.218257 0.259671
0.29 0.59 -0.152059 0.204463 0.254808
0.31 0.59 -0.161946 0.187726 0.247927
0.31 0.59 -0.161947 0.187726 0.247927
0.34 0.59 -0.170423 0.168596 0.239726
0.36 0.59 -0.177595 0.147530 0.230879
0.39 0.59 -0.183579 0.124907 0.222043
@@ -1022,7 +1022,7 @@
0.54 0.64 -0.176772 -0.028584 0.179068
0.56 0.64 -0.177447 -0.056962 0.186366
0.59 0.64 -0.177740 -0.085477 0.197226
0.61 0.64 -0.177485 -0.113961 0.210922
0.61 0.64 -0.177485 -0.113961 0.210921
0.64 0.64 -0.176458 -0.142170 0.226604
0.66 0.64 -0.174383 -0.169764 0.243370
0.69 0.64 -0.170939 -0.196274 0.260276
@@ -1058,7 +1058,7 @@
0.44 0.66 -0.153866 0.082763 0.174713
0.46 0.66 -0.154549 0.055784 0.164308
0.49 0.66 -0.155334 0.028385 0.157906
0.51 0.66 -0.156288 0.000597 0.156290
0.51 0.66 -0.156288 0.000597 0.156289
0.54 0.66 -0.157423 -0.027557 0.159817
0.56 0.66 -0.158690 -0.056046 0.168296
0.59 0.66 -0.159981 -0.084810 0.181071
@@ -1125,7 +1125,7 @@
0.11 0.71 -0.064931 0.286424 0.293692
0.14 0.71 -0.081730 0.303029 0.313857
0.16 0.71 -0.095581 0.307708 0.322211
0.19 0.71 -0.106003 0.303115 0.321115
0.19 0.71 -0.106002 0.303115 0.321115
0.21 0.71 -0.113003 0.291579 0.312711
0.24 0.71 -0.116908 0.275028 0.298844
0.26 0.71 -0.118215 0.254982 0.281053
@@ -1186,7 +1186,7 @@
0.64 0.74 -0.077679 -0.130864 0.152182
0.66 0.74 -0.085621 -0.160050 0.181513
0.69 0.74 -0.093957 -0.189964 0.211930
0.71 0.74 -0.102094 -0.220307 0.242814
0.71 0.74 -0.102095 -0.220307 0.242814
0.74 0.74 -0.109290 -0.250517 0.273319
0.76 0.74 -0.114650 -0.279666 0.302255
0.79 0.74 -0.117166 -0.306352 0.327993
@@ -1486,7 +1486,7 @@
0.14 0.94 0.236018 0.149848 0.279569
0.16 0.94 0.308503 0.117843 0.330244
0.19 0.94 0.370556 0.093841 0.382253
0.21 0.94 0.422981 0.075644 0.429691
0.21 0.94 0.422980 0.075644 0.429691
0.24 0.94 0.467064 0.061599 0.471109
0.26 0.94 0.504083 0.050535 0.506610
0.29 0.94 0.535145 0.041634 0.536762

7
BasicSolver/3D-mpi-io/canal.par
View File

@@ -31,9 +31,9 @@ p_init 0.0 # initial value for pressure
xlength 30.0 # domain size in x-direction
ylength 4.0 # domain size in y-direction
zlength 4.0 # domain size in z-direction
imax 200 # number of interior cells in x-direction
jmax 50 # number of interior cells in y-direction
kmax 50 # number of interior cells in z-direction
imax 40 # number of interior cells in x-direction
jmax 10 # number of interior cells in y-direction
kmax 10 # number of interior cells in z-direction
# Time Data:
# ---------
@@ -47,6 +47,7 @@ tau 0.5 # safety factor for time stepsize control (<0 constant delt)
itermax 500 # maximal number of pressure iteration in one time step
eps 0.0001 # stopping tolerance for pressure iteration
rho 0.52
omg 1.3 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
#===============================================================================

2
BasicSolver/3D-mpi-io/config.mk
View File

@@ -1,5 +1,5 @@
# Supported: GCC, CLANG, ICC
TAG ?= ICC
TAG ?= CLANG
ENABLE_OPENMP ?= false
#Feature options

BIN
BasicSolver/3D-mpi-io/dcavity-p4.vtk Normal file
View File

Binary file not shown.

7
BasicSolver/3D-mpi-io/dcavity.par
View File

@@ -31,9 +31,9 @@ p_init 0.0 # initial value for pressure
xlength 1.0 # domain size in x-direction
ylength 1.0 # domain size in y-direction
zlength 1.0 # domain size in z-direction
imax 128 # number of interior cells in x-direction
jmax 128 # number of interior cells in y-direction
kmax 128 # number of interior cells in z-direction
imax 40 # number of interior cells in x-direction
jmax 40 # number of interior cells in y-direction
kmax 40 # number of interior cells in z-direction
# Time Data:
# ---------
@@ -47,6 +47,7 @@ tau 0.5 # safety factor for time stepsize control (<0 constant delt)
itermax 1000 # maximal number of pressure iteration in one time step
eps 0.001 # stopping tolerance for pressure iteration
rho 0.5
omg 1.8 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
#===============================================================================

19
BasicSolver/3D-mpi-io/src/main.c
View File

@@ -19,9 +19,15 @@
#include "timing.h"
#include "vtkWriter.h"
enum VARIANT { SOR = 1, RB, RBA };
int main(int argc, char** argv)
{
int variant = RB;
int rank;
double timeStart, timeStop;
Parameter params;
Solver solver;
@@ -30,12 +36,21 @@ int main(int argc, char** argv)
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
initParameter(&params);
if (argc != 2) {
if (argc < 2) {
printf("Usage: %s <configFile>\n", argv[0]);
exit(EXIT_SUCCESS);
}
readParameter(&params, argv[1]);
if (argc == 3)
{
variant = atoi(argv[2]);
}
void (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
if (commIsMaster(&solver.comm)) {
printParameter(&params);
}
@@ -57,7 +72,7 @@ int main(int argc, char** argv)
computeFG(&solver);
computeRHS(&solver);
// if (nt % 100 == 0) normalizePressure(&solver);
solve(&solver);
(*solver_generic[variant - 1])(&solver);
adaptUV(&solver);
t += solver.dt;
nt++;

3
BasicSolver/3D-mpi-io/src/parameter.c
View File

@@ -26,6 +26,7 @@ void initParameter(Parameter* param)
param->re = 100.0;
param->gamma = 0.9;
param->tau = 0.5;
param->rho = 0.99;
}
void readParameter(Parameter* param, const char* filename)
@@ -86,6 +87,8 @@ void readParameter(Parameter* param, const char* filename)
PARSE_REAL(v_init);
PARSE_REAL(w_init);
PARSE_REAL(p_init);
PARSE_REAL(rho);
}
}

2
BasicSolver/3D-mpi-io/src/parameter.h
View File

@@ -11,7 +11,7 @@ typedef struct {
int imax, jmax, kmax;
double xlength, ylength, zlength;
int itermax;
double eps, omg;
double eps, omg, rho;
double re, tau, gamma;
double te, dt;
double gx, gy, gz;

263
BasicSolver/3D-mpi-io/src/solver.c
View File

@@ -103,6 +103,8 @@ void initSolver(Solver* s, Parameter* params)
s->te = params->te;
s->tau = params->tau;
s->gamma = params->gamma;
s->rho = params->rho;
commInit(&s->comm, s->grid.kmax, s->grid.jmax, s->grid.imax);
/* allocate arrays */
@@ -288,6 +290,267 @@ void solve(Solver* s)
#endif
}
void solveRB(Solver* s)
{
int imaxLocal = s->comm.imaxLocal;
int jmaxLocal = s->comm.jmaxLocal;
int kmaxLocal = s->comm.kmaxLocal;
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;
int pass, ksw, jsw, isw;
while ((res >= epssq) && (it < itermax)) {
res = 0.0;
ksw = 1;
for (pass = 0; pass < 2; pass++) {
jsw = ksw;
commExchange(&s->comm, p);
for (int k = 1; k < kmaxLocal + 1; k++) {
isw = jsw;
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = isw; i < imaxLocal + 1; i += 2) {
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);
}
isw = 3 - isw;
}
jsw = 3 - jsw;
}
ksw = 3 - ksw;
}
if (commIsBoundary(&s->comm, FRONT)) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, j, 0) = P(i, j, 1);
}
}
}
if (commIsBoundary(&s->comm, BACK)) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal);
}
}
}
if (commIsBoundary(&s->comm, BOTTOM)) {
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, 0, k) = P(i, 1, k);
}
}
}
if (commIsBoundary(&s->comm, TOP)) {
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k);
}
}
}
if (commIsBoundary(&s->comm, LEFT)) {
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
P(0, j, k) = P(1, j, k);
}
}
}
if (commIsBoundary(&s->comm, RIGHT)) {
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
P(imaxLocal + 1, j, k) = P(imaxLocal, j, k);
}
}
}
commReduction(&res, SUM);
res = res / (double)(imax * jmax * kmax);
#ifdef DEBUG
if (commIsMaster(&s->comm)) {
printf("%d Residuum: %e\n", it, res);
}
#endif
//commExchange(&s->comm, p);
it++;
}
#ifdef VERBOSE
if (commIsMaster(&s->comm)) {
printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
}
#endif
}
void solveRBA(Solver* s)
{
int imaxLocal = s->comm.imaxLocal;
int jmaxLocal = s->comm.jmaxLocal;
int kmaxLocal = s->comm.kmaxLocal;
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 = 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;
double rho = s->rho;
double omega = 1.0;
int pass, ksw, jsw, isw;
while ((res >= epssq) && (it < itermax)) {
res = 0.0;
ksw = 1;
for (pass = 0; pass < 2; pass++) {
jsw = ksw;
commExchange(&s->comm, p);
for (int k = 1; k < kmaxLocal + 1; k++) {
isw = jsw;
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = isw; i < imaxLocal + 1; i += 2) {
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) -= (omega * factor * r);
res += (r * r);
}
isw = 3 - isw;
}
jsw = 3 - jsw;
}
ksw = 3 - ksw;
omega = (it == 0 && pass == 0 ? 1.0 / (1.0 - 0.5 * rho * rho)
: 1.0 / (1.0 - 0.25 * rho * rho * omega));
}
if (commIsBoundary(&s->comm, FRONT)) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, j, 0) = P(i, j, 1);
}
}
}
if (commIsBoundary(&s->comm, BACK)) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal);
}
}
}
if (commIsBoundary(&s->comm, BOTTOM)) {
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, 0, k) = P(i, 1, k);
}
}
}
if (commIsBoundary(&s->comm, TOP)) {
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int i = 1; i < imaxLocal + 1; i++) {
P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k);
}
}
}
if (commIsBoundary(&s->comm, LEFT)) {
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
P(0, j, k) = P(1, j, k);
}
}
}
if (commIsBoundary(&s->comm, RIGHT)) {
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
P(imaxLocal + 1, j, k) = P(imaxLocal, j, k);
}
}
}
commReduction(&res, SUM);
res = res / (double)(imax * jmax * kmax);
#ifdef DEBUG
if (commIsMaster(&s->comm)) {
printf("%d Residuum: %e\n", it, res);
}
#endif
//commExchange(&s->comm, p);
it++;
}
#ifdef VERBOSE
if (commIsMaster(&s->comm)) {
printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
}
#endif
}
static double maxElement(Solver* s, double* m)
{
int size = (s->comm.imaxLocal + 2) * (s->comm.jmaxLocal + 2) *

4
BasicSolver/3D-mpi-io/src/solver.h
View File

@@ -20,7 +20,7 @@ typedef struct {
double *f, *g, *h;
double *u, *v, *w;
/* parameters */
double eps, omega;
double eps, omega, rho;
double re, tau, gamma;
double gx, gy, gz;
/* time stepping */
@@ -36,6 +36,8 @@ typedef struct {
void initSolver(Solver*, Parameter*);
void computeRHS(Solver*);
void solve(Solver*);
void solveRB(Solver*);
void solveRBA(Solver*);
void normalizePressure(Solver*);
void computeTimestep(Solver*);
void setBoundaryConditions(Solver*);

53
BasicSolver/3D-mpi/src/main.c
View File

@@ -24,7 +24,7 @@ enum VARIANT { SOR = 1, RB, RBA };
int main(int argc, char** argv)
{
int rank;
int variant = SOR;
int variant = RB;
double timeStart, timeStop;
Parameter params;
@@ -57,60 +57,22 @@ int main(int argc, char** argv)
double t = 0.0;
int nt = 0;
void (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
timeStart = getTimeStamp();
switch (variant) {
case SOR:
printf("Plain SOR\n");
while (t <= te) {
while (t <= te) {
if (tau > 0.0) computeTimestep(&solver);
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
// if (nt % 100 == 0) normalizePressure(&solver);
solve(&solver);
(*solver_generic[variant - 1])(&solver);
adaptUV(&solver);
t += solver.dt;
nt++;
}
break;
case RB:
printf("Red-black SOR\n");
while (t <= te) {
if (tau > 0.0) computeTimestep(&solver);
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
// if (nt % 100 == 0) normalizePressure(&solver);
solveRB(&solver);
adaptUV(&solver);
t += solver.dt;
nt++;
}
break;
case RBA:
printf("Red-black SOR with acceleration\n");
while (t <= te) {
if (tau > 0.0) computeTimestep(&solver);
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
// if (nt % 100 == 0) normalizePressure(&solver);
solveRBA(&solver);
adaptUV(&solver);
t += solver.dt;
nt++;
}
break;
}
#ifdef VERBOSE
if (rank == 0) {
@@ -119,6 +81,9 @@ int main(int argc, char** argv)
#else
printProgress(t);
#endif
}
timeStop = getTimeStamp();
stopProgress();

1
BasicSolver/3D-mpi/src/solver.c
View File

@@ -449,7 +449,6 @@ void solveRBA(Solver* s)
int it = 0;
double res = 1.0;
int pass, ksw, jsw, isw;
commExchange(&s->comm, p);
while ((res >= epssq) && (it < itermax)) {
res = 0.0;

49
BasicSolver/3D-seq/src/main.c
View File

@@ -71,7 +71,7 @@ static void createBulkArrays(Solver* s, double* pg, double* ug, double* vg, doub
int main(int argc, char** argv)
{
int rank;
int variant = SOR;
int variant = RB;
double timeStart, timeStop;
Parameter params;
@@ -100,53 +100,18 @@ int main(int argc, char** argv)
timeStart = getTimeStamp();
switch (variant) {
case SOR:
printf("Plain SOR\n");
while (t <= te) {
void (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
while (t <= te) {
if (tau > 0.0) computeTimestep(&s);
setBoundaryConditions(&s);
setSpecialBoundaryCondition(&s);
computeFG(&s);
computeRHS(&s);
solve(&s);
(*solver_generic[variant - 1])(&s);
adaptUV(&s);
t += s.dt;
}
break;
case RB:
printf("Red-black SOR\n");
while (t <= te) {
if (tau > 0.0) computeTimestep(&s);
setBoundaryConditions(&s);
setSpecialBoundaryCondition(&s);
computeFG(&s);
computeRHS(&s);
solveRB(&s);
adaptUV(&s);
t += s.dt;
}
break;
case RBA:
printf("Red-black SOR with acceleration\n");
while (t <= te) {
if (tau > 0.0) computeTimestep(&s);
setBoundaryConditions(&s);
setSpecialBoundaryCondition(&s);
computeFG(&s);
computeRHS(&s);
solveRBA(&s);
adaptUV(&s);
t += s.dt;
}
break;
}
#ifdef VERBOSE
if (rank == 0) {
printf("TIME %f , TIMESTEP %f\n", t, s.dt);
@@ -154,7 +119,7 @@ int main(int argc, char** argv)
#else
printProgress(t);
#endif
}
timeStop = getTimeStamp();
#ifndef VERBOSE
stopProgress();