diff --git a/.gitignore b/.gitignore index cd531cf..71520c4 100644 --- a/.gitignore +++ b/.gitignore @@ -30,6 +30,7 @@ *.dylib # Executables +exe-* *.exe *.out *.app diff --git a/BasicSolver/2D-seq/src/main.c b/BasicSolver/2D-seq/src/main.c index 0a28cde..65b42ab 100644 --- a/BasicSolver/2D-seq/src/main.c +++ b/BasicSolver/2D-seq/src/main.c @@ -47,7 +47,7 @@ int main(int argc, char** argv) computeFG(&solver); computeRHS(&solver); if (nt % 100 == 0) normalizePressure(&solver); - solve(&solver); + solveRB(&solver); adaptUV(&solver); t += solver.dt; nt++; diff --git a/BasicSolver/2D-seq/src/solver.c b/BasicSolver/2D-seq/src/solver.c index 6102821..1005a93 100644 --- a/BasicSolver/2D-seq/src/solver.c +++ b/BasicSolver/2D-seq/src/solver.c @@ -185,6 +185,130 @@ void solve(Solver* solver) #endif } +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++; + } + +#ifdef VERBOSE + printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); +#endif +} + +void solveRBA(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 = 0.5 * (dx2 * dy2) / (dx2 + dy2); + double omega = solver->omega; + double* p = solver->p; + double* rhs = solver->rhs; + double epssq = eps * eps; + int it = 0; + double anorm = 1.0, anormf = 0.0; + int pass, jsw, isw; + + while ((anorm >= eps * anormf) && (it < itermax)) { + anorm = 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); + + anorm += fabs(r); + P(i, j) -= (omega * factor * 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); + } + +#ifdef DEBUG + printf("%d Residuum: %e\n", it, res); +#endif + it++; + } + +#ifdef VERBOSE + printf("Solver took %d iterations to reach %f\n", it, anorm); +#endif +} + static double maxElement(Solver* solver, double* m) { int size = (solver->imax + 2) * (solver->jmax + 2); diff --git a/BasicSolver/2D-seq/src/solver.h b/BasicSolver/2D-seq/src/solver.h index 1869d6d..56f0297 100644 --- a/BasicSolver/2D-seq/src/solver.h +++ b/BasicSolver/2D-seq/src/solver.h @@ -31,14 +31,16 @@ typedef struct { int bcLeft, bcRight, bcBottom, bcTop; } Solver; -void initSolver(Solver*, Parameter*); -void computeRHS(Solver*); -void solve(Solver*); -void normalizePressure(Solver*); -void computeTimestep(Solver*); -void setBoundaryConditions(Solver*); -void setSpecialBoundaryCondition(Solver*); -void computeFG(Solver*); -void adaptUV(Solver*); -void writeResult(Solver*); +extern void initSolver(Solver*, Parameter*); +extern void computeRHS(Solver*); +extern void solve(Solver*); +extern void solveRB(Solver*); +extern void solveRBA(Solver*); +extern void normalizePressure(Solver*); +extern void computeTimestep(Solver*); +extern void setBoundaryConditions(Solver*); +extern void setSpecialBoundaryCondition(Solver*); +extern void computeFG(Solver*); +extern void adaptUV(Solver*); +extern void writeResult(Solver*); #endif