From a49a31f281241c84923e73dd0a4a7ae2ef94b3c4 Mon Sep 17 00:00:00 2001
From: Jan Eitzinger <jan@moebiusband.org>
Date: Sat, 17 Jun 2023 09:19:04 +0200
Subject: [PATCH] Refactor and add red-black solver
---
BasicSolver/3D-seq/src/main.c | 37 +++++++-------
BasicSolver/3D-seq/src/progress.h | 2 +-
BasicSolver/3D-seq/src/solver.c | 85 +++++++++++++++++++++++++++++++
BasicSolver/3D-seq/src/solver.h | 1 +
4 files changed, 104 insertions(+), 21 deletions(-)
diff --git a/BasicSolver/3D-seq/src/main.c b/BasicSolver/3D-seq/src/main.c
index 4aed6bd..73f9df2 100644
--- a/BasicSolver/3D-seq/src/main.c
+++ b/BasicSolver/3D-seq/src/main.c
@@ -69,7 +69,7 @@ int main(int argc, char** argv)
{
double timeStart, timeStop;
Parameter params;
- Solver solver;
+ Solver s;
initParameter(¶ms);
if (argc != 2) {
@@ -79,30 +79,28 @@ int main(int argc, char** argv)
readParameter(¶ms, argv[1]);
printParameter(¶ms);
- initSolver(&solver, ¶ms);
+ initSolver(&s, ¶ms);
#ifndef VERBOSE
- initProgress(solver.te);
+ initProgress(s.te);
#endif
- double tau = solver.tau;
- double te = solver.te;
+ double tau = s.tau;
+ double te = s.te;
double t = 0.0;
- int nt = 0;
timeStart = getTimeStamp();
while (t <= te) {
- if (tau > 0.0) computeTimestep(&solver);
- setBoundaryConditions(&solver);
- setSpecialBoundaryCondition(&solver);
- computeFG(&solver);
- computeRHS(&solver);
- solve(&solver);
- adaptUV(&solver);
- t += solver.dt;
- nt++;
+ if (tau > 0.0) computeTimestep(&s);
+ setBoundaryConditions(&s);
+ setSpecialBoundaryCondition(&s);
+ computeFG(&s);
+ computeRHS(&s);
+ solveRB(&s);
+ adaptUV(&s);
+ t += s.dt;
#ifdef VERBOSE
- printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
+ printf("TIME %f , TIMESTEP %f\n", t, s.dt);
#else
printProgress(t);
#endif
@@ -115,8 +113,7 @@ int main(int argc, char** argv)
double *pg, *ug, *vg, *wg;
- size_t bytesize = solver.grid.imax * solver.grid.jmax * solver.grid.kmax *
- sizeof(double);
+ size_t bytesize = (size_t)(s.grid.imax * s.grid.jmax * s.grid.kmax) * sizeof(double);
pg = allocate(64, bytesize);
ug = allocate(64, bytesize);
@@ -124,8 +121,8 @@ int main(int argc, char** argv)
wg = allocate(64, bytesize);
createBulkArrays(&s, pg, ug, vg, wg);
- VtkOptions opts = { .grid = solver.grid };
- vtkOpen(&opts, solver.problem);
+ VtkOptions opts = { .grid = s.grid };
+ vtkOpen(&opts, s.problem);
vtkScalar(&opts, "pressure", pg);
vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
vtkClose(&opts);
diff --git a/BasicSolver/3D-seq/src/progress.h b/BasicSolver/3D-seq/src/progress.h
index 9ef2d96..4d02cdb 100644
--- a/BasicSolver/3D-seq/src/progress.h
+++ b/BasicSolver/3D-seq/src/progress.h
@@ -9,6 +9,6 @@
extern void initProgress(double);
extern void printProgress(double);
-extern void stopProgress();
+extern void stopProgress(void);
#endif
diff --git a/BasicSolver/3D-seq/src/solver.c b/BasicSolver/3D-seq/src/solver.c
index bc3e39f..c741ece 100644
--- a/BasicSolver/3D-seq/src/solver.c
+++ b/BasicSolver/3D-seq/src/solver.c
@@ -221,6 +221,91 @@ void solve(Solver* s)
#endif
}
+void solveRB(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;
+ int pass, ksw, jsw, isw;
+
+ while ((res >= epssq) && (it < itermax)) {
+ res = 0.0;
+ ksw = 1;
+
+ for (pass = 0; pass < 2; pass++) {
+ jsw = ksw;
+
+ for (int k = 1; k < kmax + 1; k++) {
+ isw = jsw;
+ for (int j = 1; j < jmax + 1; j++) {
+ for (int i = isw; i < imax + 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;
+ }
+
+ 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
+}
+
static double maxElement(Solver* s, double* m)
{
int size = (s->grid.imax + 2) * (s->grid.jmax + 2) * (s->grid.kmax + 2);
diff --git a/BasicSolver/3D-seq/src/solver.h b/BasicSolver/3D-seq/src/solver.h
index 63505d0..4cec259 100644
--- a/BasicSolver/3D-seq/src/solver.h
+++ b/BasicSolver/3D-seq/src/solver.h
@@ -34,6 +34,7 @@ typedef struct {
extern void initSolver(Solver*, Parameter*);
extern void computeRHS(Solver*);
extern void solve(Solver*);
+extern void solveRB(Solver*);
extern void normalizePressure(Solver*);
extern void computeTimestep(Solver*);
extern void setBoundaryConditions(Solver*);