/* * 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 #include #include #include #include #include "allocate.h" #include "comm.h" #include "discretization.h" #include "parameter.h" #include "util.h" static double distance(double i, double j, double iCenter, double jCenter) { return sqrt(pow(iCenter - i, 2) + pow(jCenter - j, 2) * 1.0); } double sumOffset(double* sizes, int init, int offset, int coord) { double sum = 0; for (int i = init - offset; coord > 0; i -= offset, --coord) { sum += sizes[i]; } return sum; } void print(Discretization* d, double* grid) { int imaxLocal = d->comm.imaxLocal; for (int i = 0; i < d->comm.size; i++) { if (i == d->comm.rank) { sleep(1 * d->comm.rank); printf("### RANK %d LVL " "###################################################### #\n ", d->comm.rank); for (int j = 0; j < d->comm.jmaxLocal + 2; j++) { printf("%02d: ", j); for (int i = 0; i < d->comm.imaxLocal + 2; i++) { printf("%2.2f ", grid[j * (imaxLocal + 2) + i]); } printf("\n"); } fflush(stdout); } } } static void printConfig(Discretization* d) { if (commIsMaster(&d->comm)) { printf("Parameters for #%s#\n", d->problem); printf("BC Left:%d Right:%d Bottom:%d Top:%d\n", d->bcLeft, d->bcRight, d->bcBottom, d->bcTop); printf("\tReynolds number: %.2f\n", d->re); printf("\tGx Gy: %.2f %.2f\n", d->gx, d->gy); printf("Geometry data:\n"); printf("\tDomain box size (x, y): %.2f, %.2f\n", d->grid.xlength, d->grid.ylength); printf("\tCells (x, y): %d, %d\n", d->grid.imax, d->grid.jmax); printf("\tCell size (dx, dy): %f, %f\n", d->grid.dx, d->grid.dy); printf("Timestep parameters:\n"); printf("\tDefault stepsize: %.2f, Final time %.2f\n", d->dt, d->te); printf("\tdt bound: %.6f\n", d->dtBound); printf("\tTau factor: %.2f\n", d->tau); printf("Iterative s parameters:\n"); printf("\tgamma factor: %f\n", d->gamma); } commPrintConfig(&d->comm); } void initDiscretiztion(Discretization* d, Parameter* params) { d->problem = params->name; d->bcLeft = params->bcLeft; d->bcRight = params->bcRight; d->bcBottom = params->bcBottom; d->bcTop = params->bcTop; d->grid.imax = params->imax; d->grid.jmax = params->jmax; d->grid.xlength = params->xlength; d->grid.ylength = params->ylength; d->grid.dx = params->xlength / params->imax; d->grid.dy = params->ylength / params->jmax; d->re = params->re; d->gx = params->gx; d->gy = params->gy; d->dt = params->dt; d->te = params->te; d->tau = params->tau; d->gamma = params->gamma; /* allocate arrays */ int imaxLocal = d->comm.imaxLocal; int jmaxLocal = d->comm.jmaxLocal; size_t size = (imaxLocal + 2) * (jmaxLocal + 2); d->u = allocate(64, size * sizeof(double)); d->v = allocate(64, size * sizeof(double)); d->p = allocate(64, size * sizeof(double)); d->rhs = allocate(64, size * sizeof(double)); d->f = allocate(64, size * sizeof(double)); d->g = allocate(64, size * sizeof(double)); d->grid.s = allocate(64, size * sizeof(double)); for (int i = 0; i < size; i++) { d->u[i] = params->u_init; d->v[i] = params->v_init; d->p[i] = params->p_init; d->rhs[i] = 0.0; d->f[i] = 0.0; d->g[i] = 0.0; d->grid.s[i] = FLUID; } double dx = d->grid.dx; double dy = d->grid.dy; double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy); d->dtBound = 0.5 * d->re * 1.0 / invSqrSum; d->xLocal = d->comm.imaxLocal * d->grid.dx; d->yLocal = d->comm.jmaxLocal * d->grid.dy; double xLocal[d->comm.size]; double yLocal[d->comm.size]; #ifdef _MPI MPI_Allgather(&d->xLocal, 1, MPI_DOUBLE, xLocal, 1, MPI_DOUBLE, d->comm.comm); MPI_Allgather(&d->yLocal, 1, MPI_DOUBLE, yLocal, 1, MPI_DOUBLE, d->comm.comm); d->xOffset = sumOffset(xLocal, d->comm.rank, d->comm.dims[JDIM], d->comm.coords[IDIM]); d->yOffset = sumOffset(yLocal, d->comm.rank, 1, d->comm.coords[JDIM]); d->xOffsetEnd = d->xOffset + d->xLocal; d->yOffsetEnd = d->yOffset + d->yLocal; #else d->xOffset = 0; d->yOffset = 0; d->xOffsetEnd = d->xOffset + d->xLocal; d->yOffsetEnd = d->yOffset + d->yLocal; #endif printf("Rank : %d, xOffset : %.2f, yOffset : %.2f, xOffsetEnd : %.2f, yOffsetEnd : " "%.2f\n", d->comm.rank, d->xOffset, d->yOffset, d->xOffsetEnd, d->yOffsetEnd); double* s = d->grid.s; int iOffset = 0, jOffset = 0; double xCenter = 0, yCenter = 0, radius = 0; double x1 = 0, x2 = 0, y1 = 0, y2 = 0; switch (params->shape) { case NOSHAPE: break; case RECT: x1 = params->xCenter - params->xRectLength / 2; x2 = params->xCenter + params->xRectLength / 2; y1 = params->yCenter - params->yRectLength / 2; y2 = params->yCenter + params->yRectLength / 2; iOffset = d->xOffset / dx; jOffset = d->yOffset / dy; for (int j = 1; j < jmaxLocal + 1; ++j) { for (int i = 1; i < imaxLocal + 1; ++i) { if ((x1 <= ((i + iOffset) * dx)) && (((i + iOffset) * dx) <= x2) && (y1 <= ((j + jOffset) * dy)) && (((j + jOffset) * dy) <= y2)) { S(i, j) = OBSTACLE; } } } break; case CIRCLE: xCenter = params->xCenter; yCenter = params->yCenter; radius = params->circleRadius; iOffset = d->xOffset / dx; jOffset = d->yOffset / dy; for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { if (distance(((i + iOffset) * dx), ((j + jOffset) * dy), xCenter, yCenter) <= radius) { S(i, j) = OBSTACLE; } } } break; default: break; } #ifdef _MPI commExchange(&d->comm, s); #endif for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { if (S(i, j - 1) == FLUID && S(i, j + 1) == OBSTACLE && S(i, j) == OBSTACLE) S(i, j) = BOTTOM; // BOTTOM if (S(i - 1, j) == FLUID && S(i + 1, j) == OBSTACLE && S(i, j) == OBSTACLE) S(i, j) = LEFT; // LEFT if (S(i + 1, j) == FLUID && S(i - 1, j) == OBSTACLE && S(i, j) == OBSTACLE) S(i, j) = RIGHT; // RIGHT if (S(i, j + 1) == FLUID && S(i, j - 1) == OBSTACLE && S(i, j) == OBSTACLE) S(i, j) = TOP; // TOP if (S(i - 1, j - 1) == FLUID && S(i, j - 1) == FLUID && S(i - 1, j) == FLUID && S(i + 1, j + 1) == OBSTACLE && (S(i, j) == OBSTACLE || S(i, j) == LEFT || S(i, j) == BOTTOM)) S(i, j) = BOTTOMLEFT; // BOTTOMLEFT if (S(i + 1, j - 1) == FLUID && S(i, j - 1) == FLUID && S(i + 1, j) == FLUID && S(i - 1, j + 1) == OBSTACLE && (S(i, j) == OBSTACLE || S(i, j) == RIGHT || S(i, j) == BOTTOM)) S(i, j) = BOTTOMRIGHT; // BOTTOMRIGHT if (S(i - 1, j + 1) == FLUID && S(i - 1, j) == FLUID && S(i, j + 1) == FLUID && S(i + 1, j - 1) == OBSTACLE && (S(i, j) == OBSTACLE || S(i, j) == LEFT || S(i, j) == TOP)) S(i, j) = TOPLEFT; // TOPLEFT if (S(i + 1, j + 1) == FLUID && S(i + 1, j) == FLUID && S(i, j + 1) == FLUID && S(i - 1, j - 1) == OBSTACLE && (S(i, j) == OBSTACLE || S(i, j) == RIGHT || S(i, j) == TOP)) S(i, j) = TOPRIGHT; // TOPRIGHT } } #ifdef VERBOSE printConfig(d); #endif } void computeRHS(Discretization* d) { int imaxLocal = d->comm.imaxLocal; int jmaxLocal = d->comm.jmaxLocal; double idx = 1.0 / d->grid.dx; double idy = 1.0 / d->grid.dy; double idt = 1.0 / d->dt; double* rhs = d->rhs; double* f = d->f; double* g = d->g; commShift(&d->comm, f, g); for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { RHS(i, j) = ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy) * idt; } } } static double maxElement(Discretization* d, double* m) { int imaxLocal = d->comm.imaxLocal; int jmaxLocal = d->comm.jmaxLocal; int size = (imaxLocal + 2) * (jmaxLocal + 2); double maxval = DBL_MIN; for (int i = 0; i < size; i++) { maxval = MAX(maxval, fabs(m[i])); } commReduction(&maxval, MAX); return maxval; } void computeTimestep(Discretization* d) { double dt = d->dtBound; double dx = d->grid.dx; double dy = d->grid.dy; double umax = maxElement(d, d->u); double vmax = maxElement(d, d->v); if (umax > 0) { dt = (dt > dx / umax) ? dx / umax : dt; } if (vmax > 0) { dt = (dt > dy / vmax) ? dy / vmax : dt; } d->dt = dt * d->tau; } void setBoundaryConditions(Discretization* d) { int imaxLocal = d->comm.imaxLocal; int jmaxLocal = d->comm.jmaxLocal; double* u = d->u; double* v = d->v; if (commIsBoundary(&d->comm, T)) { switch (d->bcTop) { case NOSLIP: for (int i = 1; i < imaxLocal + 1; i++) { V(i, jmaxLocal) = 0.0; U(i, jmaxLocal + 1) = -U(i, jmaxLocal); } break; case SLIP: for (int i = 1; i < imaxLocal + 1; i++) { V(i, jmaxLocal) = 0.0; U(i, jmaxLocal + 1) = U(i, jmaxLocal); } break; case OUTFLOW: for (int i = 1; i < imaxLocal + 1; i++) { U(i, jmaxLocal + 1) = U(i, jmaxLocal); V(i, jmaxLocal) = V(i, jmaxLocal - 1); } break; case PERIODIC: break; } } if (commIsBoundary(&d->comm, B)) { switch (d->bcBottom) { case NOSLIP: for (int i = 1; i < imaxLocal + 1; i++) { V(i, 0) = 0.0; U(i, 0) = -U(i, 1); } break; case SLIP: for (int i = 1; i < imaxLocal + 1; i++) { V(i, 0) = 0.0; U(i, 0) = U(i, 1); } break; case OUTFLOW: for (int i = 1; i < imaxLocal + 1; i++) { U(i, 0) = U(i, 1); V(i, 0) = V(i, 1); } break; case PERIODIC: break; } } if (commIsBoundary(&d->comm, R)) { switch (d->bcRight) { case NOSLIP: for (int j = 1; j < jmaxLocal + 1; j++) { U(imaxLocal, j) = 0.0; V(imaxLocal + 1, j) = -V(imaxLocal, j); } break; case SLIP: for (int j = 1; j < jmaxLocal + 1; j++) { U(imaxLocal, j) = 0.0; V(imaxLocal + 1, j) = V(imaxLocal, j); } break; case OUTFLOW: for (int j = 1; j < jmaxLocal + 1; j++) { U(imaxLocal, j) = U(imaxLocal - 1, j); V(imaxLocal + 1, j) = V(imaxLocal, j); } break; case PERIODIC: break; } } if (commIsBoundary(&d->comm, L)) { switch (d->bcLeft) { case NOSLIP: for (int j = 1; j < jmaxLocal + 1; j++) { U(0, j) = 0.0; V(0, j) = -V(1, j); } break; case SLIP: for (int j = 1; j < jmaxLocal + 1; j++) { U(0, j) = 0.0; V(0, j) = V(1, j); } break; case OUTFLOW: for (int j = 1; j < jmaxLocal + 1; j++) { U(0, j) = U(1, j); V(0, j) = V(1, j); } break; case PERIODIC: break; } } } void setSpecialBoundaryCondition(Discretization* d) { int imaxLocal = d->comm.imaxLocal; int jmaxLocal = d->comm.jmaxLocal; double* u = d->u; double* s = d->grid.s; if (strcmp(d->problem, "dcavity") == 0) { if (commIsBoundary(&d->comm, T)) { for (int i = 1; i < imaxLocal + 1; i++) { U(i, jmaxLocal + 1) = 2.0 - U(i, jmaxLocal); } } } else if (strcmp(d->problem, "canal") == 0) { if (commIsBoundary(&d->comm, L)) { double ylength = d->grid.ylength; double dy = d->grid.dy; int rest = d->grid.jmax % d->comm.dims[JDIM]; int yc = d->comm.rank * (d->grid.jmax / d->comm.dims[JDIM]) + MIN(rest, d->comm.rank); double ys = dy * (yc + 0.5); double y; // printf("RANK %d yc: %d ys: %f\n",d->comm.rank, yc, ys); for (int j = 1; j < jmaxLocal + 1; j++) { y = ys + dy * (j - 0.5); U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength); } } } else if (strcmp(d->problem, "backstep") == 0) { for (int j = 1; j < jmaxLocal + 1; j++) { if (S(0, j) == FLUID) U(0, j) = 1.0; } } else if (strcmp(d->problem, "karman") == 0) { for (int j = 1; j < jmaxLocal + 1; j++) { U(0, j) = 1.0; } } /* print(solver, solver->u); */ } void setObjectBoundaryCondition(Discretization* d) { int imaxLocal = d->comm.imaxLocal; int jmaxLocal = d->comm.jmaxLocal; double* u = d->u; double* v = d->v; double* s = d->grid.s; for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { switch ((int)S(i, j)) { case TOP: U(i, j) = -U(i, j + 1); U(i - 1, j) = -U(i - 1, j + 1); V(i, j) = 0.0; break; case BOTTOM: U(i, j) = -U(i, j - 1); U(i - 1, j) = -U(i - 1, j - 1); V(i, j) = 0.0; break; case LEFT: U(i - 1, j) = 0.0; V(i, j) = -V(i - 1, j); V(i, j - 1) = -V(i - 1, j - 1); break; case RIGHT: U(i, j) = 0.0; V(i, j) = -V(i + 1, j); V(i, j - 1) = -V(i + 1, j - 1); break; case TOPLEFT: U(i, j) = -U(i, j + 1); U(i - 1, j) = 0.0; V(i, j) = 0.0; V(i, j - 1) = -V(i - 1, j - 1); break; case TOPRIGHT: U(i, j) = 0.0; U(i - 1, j) = -U(i - 1, j + 1); V(i, j) = 0.0; V(i, j - 1) = -V(i + 1, j - 1); break; case BOTTOMLEFT: U(i, j) = -U(i, j - 1); U(i - 1, j) = 0.0; V(i, j) = -V(i - 1, j); V(i, j - 1) = 0.0; break; case BOTTOMRIGHT: U(i, j) = 0.0; U(i - 1, j) = -U(i - 1, j - 1); V(i, j) = -V(i, j + 1); V(i, j - 1) = 0.0; break; } } } } void computeFG(Discretization* d) { double* u = d->u; double* v = d->v; double* f = d->f; double* g = d->g; double* s = d->grid.s; int imaxLocal = d->comm.imaxLocal; int jmaxLocal = d->comm.jmaxLocal; double gx = d->gx; double gy = d->gy; double gamma = d->gamma; double dt = d->dt; double inverseRe = 1.0 / d->re; double inverseDx = 1.0 / d->grid.dx; double inverseDy = 1.0 / d->grid.dy; double du2dx, dv2dy, duvdx, duvdy; double du2dx2, du2dy2, dv2dx2, dv2dy2; commExchange(&d->comm, u); commExchange(&d->comm, v); for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { if (S(i, j) == FLUID) { du2dx = inverseDx * 0.25 * ((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) - (U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) + gamma * inverseDx * 0.25 * (fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) + fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j))); duvdy = inverseDy * 0.25 * ((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) - (V(i, j - 1) + V(i + 1, j - 1)) * (U(i, j) + U(i, j - 1))) + gamma * inverseDy * 0.25 * (fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) + fabs(V(i, j - 1) + V(i + 1, j - 1)) * (U(i, j) - U(i, j - 1))); du2dx2 = inverseDx * inverseDx * (U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j)); du2dy2 = inverseDy * inverseDy * (U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1)); F(i, j) = U(i, j) + dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx); duvdx = inverseDx * 0.25 * ((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) - (U(i - 1, j) + U(i - 1, j + 1)) * (V(i, j) + V(i - 1, j))) + gamma * inverseDx * 0.25 * (fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) + fabs(U(i - 1, j) + U(i - 1, j + 1)) * (V(i, j) - V(i - 1, j))); dv2dy = inverseDy * 0.25 * ((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) - (V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) + gamma * inverseDy * 0.25 * (fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) + fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1))); dv2dx2 = inverseDx * inverseDx * (V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j)); dv2dy2 = inverseDy * inverseDy * (V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1)); G(i, j) = V(i, j) + dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy); } else { switch ((int)S(i, j)) { case TOP: G(i, j) = V(i, j); break; case BOTTOM: G(i, j - 1) = V(i, j - 1); break; case LEFT: F(i - 1, j) = U(i - 1, j); break; case RIGHT: F(i, j) = U(i, j); break; case TOPLEFT: F(i - 1, j) = U(i - 1, j); G(i, j) = V(i, j); break; case TOPRIGHT: F(i, j) = U(i, j); G(i, j) = V(i, j); break; case BOTTOMLEFT: F(i - 1, j) = U(i - 1, j); G(i, j - 1) = V(i, j - 1); break; case BOTTOMRIGHT: F(i, j) = U(i, j); G(i, j - 1) = V(i, j - 1); break; } } } } /* ----------------------------- boundary of F --------------------------- */ if (commIsBoundary(&d->comm, L)) { for (int j = 1; j < jmaxLocal + 1; j++) { F(0, j) = U(0, j); } } if (commIsBoundary(&d->comm, R)) { for (int j = 1; j < jmaxLocal + 1; j++) { F(imaxLocal, j) = U(imaxLocal, j); } } /* ----------------------------- boundary of G --------------------------- */ if (commIsBoundary(&d->comm, B)) { for (int i = 1; i < imaxLocal + 1; i++) { G(i, 0) = V(i, 0); } } if (commIsBoundary(&d->comm, T)) { for (int i = 1; i < imaxLocal + 1; i++) { G(i, jmaxLocal) = V(i, jmaxLocal); } } } void adaptUV(Discretization* d) { int imaxLocal = d->comm.imaxLocal; int jmaxLocal = d->comm.jmaxLocal; double* p = d->p; double* u = d->u; double* v = d->v; double* f = d->f; double* g = d->g; double factorX = d->dt / d->grid.dx; double factorY = d->dt / d->grid.dy; for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX; V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY; } } } void writeResult(Discretization* d, double* u, double* v, double* p) { int imax = d->grid.imax; int jmax = d->grid.jmax; double dx = d->grid.dx; double dy = d->grid.dy; double x = 0.0, y = 0.0; FILE* fp; fp = fopen("pressure.dat", "w"); if (fp == NULL) { printf("Error!\n"); exit(EXIT_FAILURE); } for (int j = 1; j <= jmax; j++) { y = (double)(j - 0.5) * dy; for (int i = 1; i <= imax; i++) { x = (double)(i - 0.5) * dx; fprintf(fp, "%.2f %.2f %f\n", x, y, p[j * (imax + 2) + i]); } fprintf(fp, "\n"); } fclose(fp); fp = fopen("velocity.dat", "w"); if (fp == NULL) { printf("Error!\n"); exit(EXIT_FAILURE); } for (int j = 1; j <= jmax; j++) { y = dy * (j - 0.5); for (int i = 1; i <= imax; i++) { x = dx * (i - 0.5); double velU = (u[j * (imax + 2) + i] + u[j * (imax + 2) + (i - 1)]) / 2.0; double velV = (v[j * (imax + 2) + i] + v[(j - 1) * (imax + 2) + i]) / 2.0; double len = sqrt((velU * velU) + (velV * velV)); fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, velU, velV, len); } } fclose(fp); }