/*
 * 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 <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#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));
    #ifdef _OPENMP
    #pragma omp parallel for  
    #endif
    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);
}
