From a5d9837443de35587072efd0b339081bbe47c59c Mon Sep 17 00:00:00 2001 From: Aditya ujeniya Date: Wed, 10 Jul 2024 09:46:17 +0200 Subject: [PATCH] Porting working 2D-seq and 2D-mpi with MG --- BasicSolver/2D-seq/canal.par | 8 ++++- BasicSolver/2D-seq/dcavity.par | 12 ++++++-- BasicSolver/2D-seq/src/parameter.c | 4 +++ BasicSolver/2D-seq/src/parameter.h | 1 + BasicSolver/2D-seq/src/solver-mg.c | 47 +++++++++++++++++++++++++----- BasicSolver/2D-seq/src/solver.h | 1 + 6 files changed, 62 insertions(+), 11 deletions(-) diff --git a/BasicSolver/2D-seq/canal.par b/BasicSolver/2D-seq/canal.par index f113e4a..6ed28e6 100644 --- a/BasicSolver/2D-seq/canal.par +++ b/BasicSolver/2D-seq/canal.par @@ -36,6 +36,13 @@ te 100.0 # final time dt 0.02 # time stepsize tau 0.5 # safety factor for time stepsize control (<0 constant delt) +# Multigrid data: +# --------- + +levels 3 # Multigrid levels +presmooth 5 # Pre-smoothning iterations +postsmooth 5 # Post-smoothning iterations + # Pressure Iteration Data: # ----------------------- @@ -43,5 +50,4 @@ itermax 500 # maximal number of pressure iteration in one time step eps 0.00001 # stopping tolerance for pressure iteration omg 1.8 # relaxation parameter for SOR iteration gamma 0.9 # upwind differencing factor gamma -levels 5 # Multigrid levels #=============================================================================== diff --git a/BasicSolver/2D-seq/dcavity.par b/BasicSolver/2D-seq/dcavity.par index d8e311e..c4366a1 100644 --- a/BasicSolver/2D-seq/dcavity.par +++ b/BasicSolver/2D-seq/dcavity.par @@ -26,8 +26,8 @@ p_init 0.0 # initial value for pressure xlength 1.0 # domain size in x-direction ylength 1.0 # domain size in y-direction -imax 800 # number of interior cells in x-direction -jmax 800 # number of interior cells in y-direction +imax 128 # number of interior cells in x-direction +jmax 128 # number of interior cells in y-direction # Time Data: # --------- @@ -36,6 +36,13 @@ te 10.0 # final time dt 0.02 # time stepsize tau 0.5 # safety factor for time stepsize control (<0 constant delt) +# Multigrid data: +# --------- + +levels 2 # Multigrid levels +presmooth 20 # Pre-smoothning iterations +postsmooth 5 # Post-smoothning iterations + # Solver Data: # ----------------------- @@ -44,5 +51,4 @@ eps 0.001 # stopping tolerance for pressure iteration rho 0.5 omg 1.7 # relaxation parameter for SOR iteration gamma 0.9 # upwind differencing factor gamma -levels 5 # Multigrid levels #=============================================================================== diff --git a/BasicSolver/2D-seq/src/parameter.c b/BasicSolver/2D-seq/src/parameter.c index e559545..b7c32ca 100644 --- a/BasicSolver/2D-seq/src/parameter.c +++ b/BasicSolver/2D-seq/src/parameter.c @@ -24,6 +24,8 @@ void initParameter(Parameter* param) param->gamma = 0.9; param->tau = 0.5; param->levels = 5; + param->presmooth = 5; + param->postsmooth = 5; } void readParameter(Parameter* param, const char* filename) @@ -79,6 +81,8 @@ void readParameter(Parameter* param, const char* filename) PARSE_REAL(u_init); PARSE_REAL(v_init); PARSE_REAL(p_init); + PARSE_INT(presmooth); + PARSE_INT(postsmooth); } } diff --git a/BasicSolver/2D-seq/src/parameter.h b/BasicSolver/2D-seq/src/parameter.h index 5098794..3714f96 100644 --- a/BasicSolver/2D-seq/src/parameter.h +++ b/BasicSolver/2D-seq/src/parameter.h @@ -18,6 +18,7 @@ typedef struct { char* name; int bcLeft, bcRight, bcBottom, bcTop; double u_init, v_init, p_init; + int presmooth, postsmooth; } Parameter; void initParameter(Parameter*); diff --git a/BasicSolver/2D-seq/src/solver-mg.c b/BasicSolver/2D-seq/src/solver-mg.c index d1e56b4..92537b8 100644 --- a/BasicSolver/2D-seq/src/solver-mg.c +++ b/BasicSolver/2D-seq/src/solver-mg.c @@ -84,6 +84,37 @@ static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int 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) { + + P(i, j) -= factor * (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)); + + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + +} + +static double calculateResidual(Solver* s, double* p, double* rhs, int level, int imax, int jmax) +{ + 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* r = s->r[level]; + double res = 1.0; + int pass, jsw, isw; + + jsw = 1; + for (pass = 0; pass < 2; pass++) { isw = jsw; @@ -94,7 +125,6 @@ static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int ((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(i, j)); res += (R(i, j) * R(i, j)); } isw = 3 - isw; @@ -110,7 +140,7 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax, { double res = 0.0; - // coarsest level TODO: Use direct solver? + // coarsest level if (level == COARSEST_LEVEL) { for (int i = 0; i < 5; i++) { smooth(s, p, rhs, level, imax, jmax); @@ -118,17 +148,18 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax, return res; } - // pre-smoothing TODO: Make smoothing steps configurable? - for (int i = 0; i < 5; i++) { + // pre-smoothing + for (int i = 0; i < s->presmooth; i++) { smooth(s, p, rhs, level, imax, jmax); if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax); } + res = calculateResidual(s, p, rhs, level, imax, jmax); + // restrict restrictMG(s, level, imax, jmax); // MGSolver on residual and error. - // TODO: What if there is a rest? multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, imax / 2, jmax / 2); // prolongate @@ -139,8 +170,8 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax, if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax); // post-smoothing - for (int i = 0; i < 5; i++) { - res = smooth(s, p, rhs, level, imax, jmax); + for (int i = 0; i < s->postsmooth; i++) { + smooth(s, p, rhs, level, imax, jmax); if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax); } @@ -154,6 +185,8 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) s->itermax = p->itermax; s->levels = p->levels; s->grid = &d->grid; + s->presmooth = p->presmooth; + s->postsmooth = p->postsmooth; int imax = s->grid->imax; int jmax = s->grid->jmax; diff --git a/BasicSolver/2D-seq/src/solver.h b/BasicSolver/2D-seq/src/solver.h index d7ce513..955b88d 100644 --- a/BasicSolver/2D-seq/src/solver.h +++ b/BasicSolver/2D-seq/src/solver.h @@ -18,6 +18,7 @@ typedef struct { int itermax; int levels; double **r, **e; + int presmooth, postsmooth; } Solver; extern void initSolver(Solver*, Discretization*, Parameter*);