Compare commits

..

No commits in common. "8091c714e2d78863cb6e5a7f15496c98e345efa4" and "a834a57bf8dc404ce4c988a1b1f93db7c60d80c2" have entirely different histories.

7 changed files with 136 additions and 184 deletions

View File

@ -1,8 +1,8 @@
# Supported: GCC, CLANG, ICC # Supported: GCC, CLANG, ICC
TAG ?= CLANG TAG ?= CLANG
ENABLE_OPENMP ?= false ENABLE_OPENMP ?= false
# Supported: sor, rb, mg # Supported: sor, mg
SOLVER ?= rb SOLVER ?= sor
# Run in debug settings # Run in debug settings
DEBUG ?= false DEBUG ?= false

View File

@ -26,8 +26,8 @@ p_init 0.0 # initial value for pressure
xlength 1.0 # domain size in x-direction xlength 1.0 # domain size in x-direction
ylength 1.0 # domain size in y-direction ylength 1.0 # domain size in y-direction
imax 128 # number of interior cells in x-direction imax 100 # number of interior cells in x-direction
jmax 128 # number of interior cells in y-direction jmax 100 # number of interior cells in y-direction
# Time Data: # Time Data:
# --------- # ---------

View File

@ -1,76 +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 "solver.h"
#include "util.h"
void initSolver(Solver* s, Discretization* d, Parameter* p)
{
s->grid = &d->grid;
s->itermax = p->itermax;
s->eps = p->eps;
s->omega = p->omg;
}
void solve(Solver* solver, double* p, double* rhs)
{
int imax = solver->grid->imax;
int jmax = solver->grid->jmax;
double eps = solver->eps;
int itermax = solver->itermax;
double dx2 = solver->grid->dx * solver->grid->dx;
double dy2 = solver->grid->dy * solver->grid->dy;
double idx2 = 1.0 / dx2;
double idy2 = 1.0 / dy2;
double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
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
}

View File

@ -15,7 +15,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
s->omega = p->omg; s->omega = p->omg;
} }
void solve(Solver* solver, double* p, double* rhs) void solveSOR(Solver* solver, double* p, double* rhs)
{ {
int imax = solver->grid->imax; int imax = solver->grid->imax;
int jmax = solver->grid->jmax; int jmax = solver->grid->jmax;
@ -66,3 +66,63 @@ void solve(Solver* solver, double* p, double* rhs)
printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
#endif #endif
} }
void solve(Solver* solver, double* p, double* rhs)
{
int imax = solver->grid->imax;
int jmax = solver->grid->jmax;
double eps = solver->eps;
int itermax = solver->itermax;
double dx2 = solver->grid->dx * solver->grid->dx;
double dy2 = solver->grid->dy * solver->grid->dy;
double idx2 = 1.0 / dx2;
double idy2 = 1.0 / dy2;
double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
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
}

View File

@ -1,7 +1,7 @@
# Supported: GCC, CLANG, ICC # Supported: GCC, CLANG, ICC
TAG ?= CLANG TAG ?= CLANG
ENABLE_OPENMP ?= false ENABLE_OPENMP ?= false
# Supported: sor, rb, mg # Supported: sor, mg
SOLVER ?= mg SOLVER ?= mg
# Run in debug settings # Run in debug settings
DEBUG ?= false DEBUG ?= false

View File

@ -1,93 +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 "grid.h"
#include "solver.h"
#include "util.h"
void initSolver(Solver* s, Discretization* d, Parameter* p)
{
s->grid = &d->grid;
s->itermax = p->itermax;
s->eps = p->eps;
s->omega = p->omg;
Grid* g = s->grid;
int imax = s->grid->imax;
int jmax = s->grid->jmax;
s->totalFluidCells = 0;
for (int j = 0; j < jmax + 2; j++) {
for (int i = 0; i < imax + 2; i++) {
if (gridIsFluid(g, i, j)) {
s->totalFluidCells++;
}
}
}
}
void solve(Solver* s, double* p, double* rhs)
{
int imax = s->grid->imax;
int jmax = s->grid->jmax;
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 it = 0;
Grid* g = s->grid;
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) {
if (gridIsFluid(g, i, j)) {
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)(s->totalFluidCells);
#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
}

View File

@ -29,7 +29,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
} }
} }
void solve(Solver* s, double* p, double* rhs) void solveSOR(Solver* s, double* p, double* rhs)
{ {
int imax = s->grid->imax; int imax = s->grid->imax;
int jmax = s->grid->jmax; int jmax = s->grid->jmax;
@ -42,7 +42,6 @@ void solve(Solver* s, double* p, double* rhs)
double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
double epssq = eps * eps; double epssq = eps * eps;
int it = 0; int it = 0;
Grid* g = s->grid;
double res = 1.0; double res = 1.0;
while ((res >= epssq) && (it < itermax)) { while ((res >= epssq) && (it < itermax)) {
@ -50,7 +49,7 @@ void solve(Solver* s, double* p, double* rhs)
for (int j = 1; j < jmax + 1; j++) { for (int j = 1; j < jmax + 1; j++) {
for (int i = 1; i < imax + 1; i++) { for (int i = 1; i < imax + 1; i++) {
if (gridIsFluid(g, i, j)) {
double r = RHS(i, j) - double r = RHS(i, j) -
((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + ((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 + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
@ -59,7 +58,6 @@ void solve(Solver* s, double* p, double* rhs)
res += (r * r); res += (r * r);
} }
} }
}
for (int i = 1; i < imax + 1; i++) { for (int i = 1; i < imax + 1; i++) {
P(i, 0) = P(i, 1); P(i, 0) = P(i, 1);
@ -82,3 +80,66 @@ void solve(Solver* s, double* p, double* rhs)
printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
#endif #endif
} }
void solve(Solver* solver, double* p, double* rhs)
{
int imax = solver->grid->imax;
int jmax = solver->grid->jmax;
double eps = solver->eps;
int itermax = solver->itermax;
double dx2 = solver->grid->dx * solver->grid->dx;
double dy2 = solver->grid->dy * solver->grid->dy;
double idx2 = 1.0 / dx2;
double idy2 = 1.0 / dy2;
double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
double epssq = eps * eps;
int it = 0;
Grid* g = solver->grid;
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) {
if (gridIsFluid(g, i, j)) {
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
}