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

#include "solver.h"
#include "parameter.h"
#include "allocate.h"

#define PI 3.14159265358979323846
#define P(i,j) p[(j)*(imax+2) + (i)]
#define RHS(i,j) rhs[(j)*(imax+2) + (i)]

void initSolver(Solver *solver, Parameter *params, int problem)
{
    solver->imax = params->imax;
    solver->jmax = params->jmax;
    solver->dx = params->xlength/params->imax;
    solver->dy = params->ylength/params->jmax;
    solver->eps = params->eps;
    solver->omega = params->omg;
    solver->itermax = params->itermax;

    int imax = solver->imax;
    int jmax = solver->jmax;
    size_t bytesize = (imax+2) * (jmax+2) * sizeof(double);
    solver->p = allocate(64, bytesize);
    solver->rhs = allocate(64, bytesize);

    double dx = solver->dx;
    double dy = solver->dy;
    double* p = solver->p;
    double* rhs = solver->rhs;

    for( int j=0; j<jmax+2; j++ ) {
        for( int i=0; i<imax+2; i++ ) {
            P(i,j) = sin(2.0*PI*i*dx*2.0)+sin(2.0*PI*j*dy*2.0);
        }
    }

    if(problem == 2) {
        for( int j=0; j<jmax+2; j++ ) {
            for( int i=0; i<imax+2; i++ ) {
                RHS(i,j) = sin(2.0*PI*i*dx);
            }
        }
    } else {
        for( int j=0; j<jmax+2; j++ ) {
            for( int i=0; i<imax+2; i++ ) {
                RHS(i,j) = 0.0;
            }
        }
    }
}

int solve(Solver *solver)
{
    double r;
    int it = 0;
    double res;

    int imax = solver->imax;
    int jmax = solver->jmax;
    double eps= solver->eps;
    double omega = solver->omega;
    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 = omega * 0.5 * (dx2*dy2) / (dx2+dy2);
    double* p = solver->p;
    double* rhs = solver->rhs;

    res = eps + 1.0;

    while((res >= eps) && (it < itermax)) {
        res = 0.0;

        for( int j=1; j<jmax+1; j++ ) {
            for( int i=1; i<imax+1; i++ ) {

                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);
            }
        }

        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);
        }

        it++;

        //Residuum
        res = sqrt(res / (imax*jmax));
#ifdef DEBUG
        printf("%d Residuum: %e\n",it, res);
#endif
    }

    printf("Solver took %d iterations\n",it);
    if( res < eps ){
        return 1;
    } else{
        return 0;
    }
}

void writeResult(Solver* solver)
{
    int imax = solver->imax;
    int jmax = solver->jmax;
    double* p = solver->p;

    FILE *fp;
    fp= fopen("p.dat", "w");

    if (fp== NULL) {
        printf("Error!\n");
        exit(EXIT_FAILURE);
    }

    for( int j=0; j<jmax+2; j++ ) {
        for( int i=0; i<imax+2; i++ ) {
            fprintf(fp, "%f ", P(i,j));
        }
        fprintf(fp, "\n");
    }

    fclose(fp);
}