/* * Copyright (C) 2022 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 #include "allocate.h" #include "parameter.h" #include "solver.h" #include "util.h" #define P(i, j) p[(j) * (imax + 2) + (i)] #define F(i, j) f[(j) * (imax + 2) + (i)] #define G(i, j) g[(j) * (imax + 2) + (i)] #define U(i, j) u[(j) * (imax + 2) + (i)] #define V(i, j) v[(j) * (imax + 2) + (i)] #define RHS(i, j) rhs[(j) * (imax + 2) + (i)] static int sizeOfRank(int rank, int size, int N) { return N / size + ((N % size > rank) ? 1 : 0); } static void print(Solver* solver, double* grid) { int imax = solver->imax; for (int i = 0; i < solver->size; i++) { if (i == solver->rank) { printf("### RANK %d " "#######################################################\n", solver->rank); for (int j = 0; j < solver->jmaxLocal + 2; j++) { printf("%02d: ", j); for (int i = 0; i < solver->imax + 2; i++) { printf("%12.8f ", grid[j * (imax + 2) + i]); } printf("\n"); } fflush(stdout); } MPI_Barrier(MPI_COMM_WORLD); } } static void exchange(Solver* solver, double* grid) { MPI_Request requests[4] = { MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL }; /* exchange ghost cells with top neighbor */ if (solver->rank + 1 < solver->size) { int top = solver->rank + 1; double* src = grid + (solver->jmaxLocal) * (solver->imax + 2) + 1; double* dst = grid + (solver->jmaxLocal + 1) * (solver->imax + 2) + 1; MPI_Isend(src, solver->imax, MPI_DOUBLE, top, 1, MPI_COMM_WORLD, &requests[0]); MPI_Irecv(dst, solver->imax, MPI_DOUBLE, top, 2, MPI_COMM_WORLD, &requests[1]); } /* exchange ghost cells with bottom neighbor */ if (solver->rank > 0) { int bottom = solver->rank - 1; double* src = grid + (solver->imax + 2) + 1; double* dst = grid + 1; MPI_Isend(src, solver->imax, MPI_DOUBLE, bottom, 2, MPI_COMM_WORLD, &requests[2]); MPI_Irecv(dst, solver->imax, MPI_DOUBLE, bottom, 1, MPI_COMM_WORLD, &requests[3]); } MPI_Waitall(4, requests, MPI_STATUSES_IGNORE); } static void shift(Solver* solver) { MPI_Request requests[2] = { MPI_REQUEST_NULL, MPI_REQUEST_NULL }; double* g = solver->g; /* shift G */ /* receive ghost cells from bottom neighbor */ if (solver->rank > 0) { int bottom = solver->rank - 1; MPI_Irecv(g + 1, solver->imax, MPI_DOUBLE, bottom, 0, MPI_COMM_WORLD, &requests[0]); } if (solver->rank + 1 < solver->size) { int top = solver->rank + 1; double* buf = g + (solver->jmaxLocal) * (solver->imax + 2) + 1; /* send ghost cells to top neighbor */ MPI_Isend(buf, solver->imax, MPI_DOUBLE, top, 0, MPI_COMM_WORLD, &requests[1]); } MPI_Waitall(2, requests, MPI_STATUSES_IGNORE); } void collectResult(Solver* solver) { double* Pall = NULL; double* Uall = NULL; double* Vall = NULL; int *rcvCounts, *displs; if (solver->rank == 0) { Pall = allocate(64, (solver->imax + 2) * (solver->jmax + 2) * sizeof(double)); Uall = allocate(64, (solver->imax + 2) * (solver->jmax + 2) * sizeof(double)); Vall = allocate(64, (solver->imax + 2) * (solver->jmax + 2) * sizeof(double)); rcvCounts = (int*)malloc(solver->size * sizeof(int)); displs = (int*)malloc(solver->size * sizeof(int)); rcvCounts[0] = solver->jmaxLocal * (solver->imax + 2); displs[0] = 0; int cursor = rcvCounts[0]; for (int i = 1; i < solver->size; i++) { rcvCounts[i] = sizeOfRank(i, solver->size, solver->jmax) * (solver->imax + 2); displs[i] = cursor; cursor += rcvCounts[i]; } } int cnt = solver->jmaxLocal * (solver->imax + 2); double* sendbuffer = solver->p + (solver->imax + 2); MPI_Gatherv(sendbuffer, cnt, MPI_DOUBLE, Pall, rcvCounts, displs, MPI_DOUBLE, 0, MPI_COMM_WORLD); sendbuffer = solver->u + (solver->imax + 2); MPI_Gatherv(sendbuffer, cnt, MPI_DOUBLE, Uall, rcvCounts, displs, MPI_DOUBLE, 0, MPI_COMM_WORLD); sendbuffer = solver->v + (solver->imax + 2); MPI_Gatherv(sendbuffer, cnt, MPI_DOUBLE, Vall, rcvCounts, displs, MPI_DOUBLE, 0, MPI_COMM_WORLD); if (solver->rank == 0) { writeResult(solver, Pall, Uall, Vall); } } static void printConfig(Solver* solver) { if (solver->rank == 0) { printf("Parameters for #%s#\n", solver->problem); printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n", solver->bcLeft, solver->bcRight, solver->bcBottom, solver->bcTop); printf("\tReynolds number: %.2f\n", solver->re); printf("\tGx Gy: %.2f %.2f\n", solver->gx, solver->gy); printf("Geometry data:\n"); printf("\tDomain box size (x, y): %.2f, %.2f\n", solver->xlength, solver->ylength); printf("\tCells (x, y): %d, %d\n", solver->imax, solver->jmax); printf("Timestep parameters:\n"); printf("\tDefault stepsize: %.2f, Final time %.2f\n", solver->dt, solver->te); printf("\tdt bound: %.6f\n", solver->dtBound); printf("\tTau factor: %.2f\n", solver->tau); printf("Iterative solver parameters:\n"); printf("\tMax iterations: %d\n", solver->itermax); printf("\tepsilon (stopping tolerance) : %f\n", solver->eps); printf("\tgamma factor: %f\n", solver->gamma); printf("\tomega (SOR relaxation): %f\n", solver->omega); printf("Communication parameters:\n"); } for (int i = 0; i < solver->size; i++) { if (i == solver->rank) { printf("\tRank %d of %d\n", solver->rank, solver->size); printf("\tLocal domain size: %dx%d\n", solver->imax, solver->jmaxLocal); fflush(stdout); } } } void initSolver(Solver* solver, Parameter* params) { MPI_Comm_rank(MPI_COMM_WORLD, &(solver->rank)); MPI_Comm_size(MPI_COMM_WORLD, &(solver->size)); solver->problem = params->name; solver->bcLeft = params->bcLeft; solver->bcRight = params->bcRight; solver->bcBottom = params->bcBottom; solver->bcTop = params->bcTop; solver->imax = params->imax; solver->jmax = params->jmax; solver->jmaxLocal = sizeOfRank(solver->rank, solver->size, solver->jmax); solver->xlength = params->xlength; solver->ylength = params->ylength; solver->dx = params->xlength / params->imax; solver->dy = params->ylength / params->jmax; solver->eps = params->eps; solver->omega = params->omg; solver->itermax = params->itermax; solver->re = params->re; solver->gx = params->gx; solver->gy = params->gy; solver->dt = params->dt; solver->te = params->te; solver->tau = params->tau; solver->gamma = params->gamma; solver->rho = params->rho; int imax = solver->imax; int jmaxLocal = solver->jmaxLocal; size_t bytesize = (imax + 2) * (jmaxLocal + 2) * sizeof(double); solver->u = allocate(64, bytesize); solver->v = allocate(64, bytesize); solver->p = allocate(64, bytesize); solver->rhs = allocate(64, bytesize); solver->f = allocate(64, bytesize); solver->g = allocate(64, bytesize); for (int i = 0; i < (imax + 2) * (jmaxLocal + 2); i++) { solver->u[i] = params->u_init; solver->v[i] = params->v_init; solver->p[i] = params->p_init; solver->rhs[i] = 0.0; solver->f[i] = 0.0; solver->g[i] = 0.0; } double dx = solver->dx; double dy = solver->dy; double inv_sqr_sum = 1.0 / (dx * dx) + 1.0 / (dy * dy); solver->dtBound = 0.5 * solver->re * 1.0 / inv_sqr_sum; #ifdef VERBOSE printConfig(solver); #endif } void computeRHS(Solver* solver) { int imax = solver->imax; int jmaxLocal = solver->jmaxLocal; double idx = 1.0 / solver->dx; double idy = 1.0 / solver->dy; double idt = 1.0 / solver->dt; double* rhs = solver->rhs; double* f = solver->f; double* g = solver->g; shift(solver); for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imax + 1; i++) { RHS(i, j) = ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy) * idt; } } } void solve(Solver* solver) { int imax = solver->imax; int jmax = solver->jmax; int jmaxLocal = solver->jmaxLocal; 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; while ((res >= epssq) && (it < itermax)) { res = 0.0; exchange(solver, p); for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imax + 1; i++) { 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); } } if (solver->rank == 0) { for (int i = 1; i < imax + 1; i++) { P(i, 0) = P(i, 1); } } if (solver->rank == (solver->size - 1)) { for (int i = 1; i < imax + 1; i++) { P(i, jmaxLocal + 1) = P(i, jmaxLocal); } } for (int j = 1; j < jmaxLocal + 1; j++) { P(0, j) = P(1, j); P(imax + 1, j) = P(imax, j); } MPI_Allreduce(MPI_IN_PLACE, &res, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); res = res / (double)(imax * jmax); #ifdef DEBUG if (solver->rank == 0) { printf("%d Residuum: %e\n", it, res); } #endif it++; } #ifdef VERBOSE if (solver->rank == 0) { printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); } #endif } void solveRB(Solver* solver) { int imax = solver->imax; int jmax = solver->jmax; int jmaxLocal = solver->jmaxLocal; 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)) { jsw = 1; for (pass = 0; pass < 2; pass++) { isw = jsw; exchange(solver, p); for (int j = 1; j < jmaxLocal + 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; } if (solver->rank == 0) { for (int i = 1; i < imax + 1; i++) { P(i, 0) = P(i, 1); } } if (solver->rank == (solver->size - 1)) { for (int i = 1; i < imax + 1; i++) { P(i, jmaxLocal + 1) = P(i, jmaxLocal); } } for (int j = 1; j < jmaxLocal + 1; j++) { P(0, j) = P(1, j); P(imax + 1, j) = P(imax, j); } MPI_Allreduce(MPI_IN_PLACE, &res, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); res = res / (double)(imax * jmax); #ifdef DEBUG if (solver->rank == 0) { printf("%d Residuum: %e\n", it, res); } #endif it++; } #ifdef VERBOSE if (solver->rank == 0) { 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; int jmaxLocal = solver->jmaxLocal; 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* p = solver->p; double* rhs = solver->rhs; double rho = solver->rho; double epssq = eps * eps; int it = 0; double res = 1.0; double omega = 1.0; int pass, jsw, isw; while ((res >= epssq) && (it < itermax)) { res = 0.0; jsw = 1; for (pass = 0; pass < 2; pass++) { isw = jsw; exchange(solver, p); for (int j = 1; j < jmaxLocal + 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) -= (omega * factor * r); res += (r * r); } isw = 3 - isw; } jsw = 3 - jsw; omega = (it == 0 && pass == 0 ? 1.0 / (1.0 - 0.5 * rho * rho) : 1.0 / (1.0 - 0.25 * rho * rho * omega)); } if (solver->rank == 0) { for (int i = 1; i < imax + 1; i++) { P(i, 0) = P(i, 1); } } if (solver->rank == (solver->size - 1)) { for (int i = 1; i < imax + 1; i++) { P(i, jmaxLocal + 1) = P(i, jmaxLocal); } } for (int j = 1; j < jmaxLocal + 1; j++) { P(0, j) = P(1, j); P(imax + 1, j) = P(imax, j); } MPI_Allreduce(MPI_IN_PLACE, &res, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); res = res / (double)(imax * jmax); #ifdef DEBUG if (solver->rank == 0) { printf("%d Residuum: %e\n", it, res); } #endif it++; } #ifdef VERBOSE if (solver->rank == 0) { printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); } #endif } static double maxElement(Solver* solver, double* m) { int size = (solver->imax + 2) * (solver->jmaxLocal + 2); double maxval = DBL_MIN; for (int i = 0; i < size; i++) { maxval = MAX(maxval, fabs(m[i])); } MPI_Allreduce(MPI_IN_PLACE, &maxval, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); return maxval; } void normalizePressure(Solver* solver) { int size = (solver->imax + 2) * (solver->jmaxLocal + 2); double* p = solver->p; double avgP = 0.0; for (int i = 0; i < size; i++) { avgP += p[i]; } MPI_Allreduce(MPI_IN_PLACE, &avgP, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); avgP /= (solver->imax + 2) * (solver->jmax + 2); for (int i = 0; i < size; i++) { p[i] = p[i] - avgP; } } void computeTimestep(Solver* solver) { double dt = solver->dtBound; double dx = solver->dx; double dy = solver->dy; double umax = maxElement(solver, solver->u); double vmax = maxElement(solver, solver->v); if (umax > 0) { dt = (dt > dx / umax) ? dx / umax : dt; } if (vmax > 0) { dt = (dt > dy / vmax) ? dy / vmax : dt; } solver->dt = dt * solver->tau; } void setBoundaryConditions(Solver* solver) { int imax = solver->imax; int jmaxLocal = solver->jmaxLocal; double* u = solver->u; double* v = solver->v; // Left boundary switch (solver->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; } // Right boundary switch (solver->bcRight) { case NOSLIP: for (int j = 1; j < jmaxLocal + 1; j++) { U(imax, j) = 0.0; V(imax + 1, j) = -V(imax, j); } break; case SLIP: for (int j = 1; j < jmaxLocal + 1; j++) { U(imax, j) = 0.0; V(imax + 1, j) = V(imax, j); } break; case OUTFLOW: for (int j = 1; j < jmaxLocal + 1; j++) { U(imax, j) = U(imax - 1, j); V(imax + 1, j) = V(imax, j); } break; case PERIODIC: break; } // Bottom boundary if (solver->rank == 0) { switch (solver->bcBottom) { case NOSLIP: for (int i = 1; i < imax + 1; i++) { V(i, 0) = 0.0; U(i, 0) = -U(i, 1); } break; case SLIP: for (int i = 1; i < imax + 1; i++) { V(i, 0) = 0.0; U(i, 0) = U(i, 1); } break; case OUTFLOW: for (int i = 1; i < imax + 1; i++) { U(i, 0) = U(i, 1); V(i, 0) = V(i, 1); } break; case PERIODIC: break; } } // Top boundary if (solver->rank == (solver->size - 1)) { switch (solver->bcTop) { case NOSLIP: for (int i = 1; i < imax + 1; i++) { V(i, jmaxLocal) = 0.0; U(i, jmaxLocal + 1) = -U(i, jmaxLocal); } break; case SLIP: for (int i = 1; i < imax + 1; i++) { V(i, jmaxLocal) = 0.0; U(i, jmaxLocal + 1) = U(i, jmaxLocal); } break; case OUTFLOW: for (int i = 1; i < imax + 1; i++) { U(i, jmaxLocal + 1) = U(i, jmaxLocal); V(i, jmaxLocal) = V(i, jmaxLocal - 1); } break; case PERIODIC: break; } } } void setSpecialBoundaryCondition(Solver* solver) { int imax = solver->imax; int jmaxLocal = solver->jmaxLocal; double* u = solver->u; if (strcmp(solver->problem, "dcavity") == 0) { if (solver->rank == (solver->size - 1)) { for (int i = 1; i < imax; i++) { U(i, jmaxLocal + 1) = 2.0 - U(i, jmaxLocal); } } } else if (strcmp(solver->problem, "canal") == 0) { double ylength = solver->ylength; double dy = solver->dy; int rest = solver->jmax % solver->size; int yc = solver->rank * (solver->jmax / solver->size) + MIN(rest, solver->rank); double ys = dy * (yc + 0.5); double y; /* printf("RANK %d yc: %d ys: %f\n", solver->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); } } /* print(solver, solver->u); */ } void computeFG(Solver* solver) { double* u = solver->u; double* v = solver->v; double* f = solver->f; double* g = solver->g; int imax = solver->imax; int jmaxLocal = solver->jmaxLocal; double gx = solver->gx; double gy = solver->gy; double gamma = solver->gamma; double dt = solver->dt; double inverseRe = 1.0 / solver->re; double inverseDx = 1.0 / solver->dx; double inverseDy = 1.0 / solver->dy; double du2dx, dv2dy, duvdx, duvdy; double du2dx2, du2dy2, dv2dx2, dv2dy2; exchange(solver, u); exchange(solver, v); for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imax + 1; i++) { 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); } } /* ----------------------------- boundary of F --------------------------- */ for (int j = 1; j < jmaxLocal + 1; j++) { F(0, j) = U(0, j); F(imax, j) = U(imax, j); } /* ----------------------------- boundary of G --------------------------- */ if (solver->rank == 0) { for (int i = 1; i < imax + 1; i++) { G(i, 0) = V(i, 0); } } if (solver->rank == (solver->size - 1)) { for (int i = 1; i < imax + 1; i++) { G(i, jmaxLocal) = V(i, jmaxLocal); } } } void adaptUV(Solver* solver) { int imax = solver->imax; int jmaxLocal = solver->jmaxLocal; double* p = solver->p; double* u = solver->u; double* v = solver->v; double* f = solver->f; double* g = solver->g; double factorX = solver->dt / solver->dx; double factorY = solver->dt / solver->dy; for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imax + 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(Solver* solver, double* p, double* u, double* v) { int imax = solver->imax; int jmax = solver->jmax; double dx = solver->dx; double dy = solver->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 + 1; j++) { y = (double)(j - 0.5) * dy; for (int i = 1; i < imax + 1; i++) { x = (double)(i - 0.5) * dx; fprintf(fp, "%.2f %.2f %f\n", x, y, P(i, j)); } 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 + 1; j++) { y = dy * (j - 0.5); for (int i = 1; i < imax + 1; i++) { x = dx * (i - 0.5); double vel_u = (U(i, j) + U(i - 1, j)) / 2.0; double vel_v = (V(i, j) + V(i, j - 1)) / 2.0; double len = sqrt((vel_u * vel_u) + (vel_v * vel_v)); fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, vel_u, vel_v, len); } } fclose(fp); }