/*
 * SOR_itit.cpp
 *
 *   Created on: Feb 10, 2014
 *  Modified on: Jul 05, 2020
 *
 *      Author: chuan
 */

#include "solver.h"

void CDelphiSolver::SOR_itit()
{
    delphi_real rmsch, rmsch2, rmxch, rmxch2;
    int itr, ires;
    delphi_real grden, grdn[5] = { 0.0, 0.0, 0.0, 0.0, 0.0 };
    delphi_integer ix, iy, iz;
    delphi_real maxres = (fRmsc > fMaxc) ? fRmsc : fMaxc;
    maxres = (maxres > fGridConverge) ? maxres : fGridConverge;
    vector<delphi_real> rmsl(nxran, 0.0), rmaxl(nxran, 0.0);
    string strLine60 = " ----------------------------------------------------------------";

    cout << strLine60 << endl;
    if (0.0 < fGridConverge)
        cout << "      " << " rms-change   max change    grid energy    #iterations" << endl;
    else
        cout << "      " << " rms-change   max change       #iterations" << endl;
    cout << strLine60 << endl;


    if (0 == iConvergeFract) 
    {
        iIterateInterval = 10;
        iConvergeFract   = 1;
    }

    if (iIterateInterval > iLinIterateNum) iIterateInterval = iLinIterateNum;

    initOddEvenItr(1); // forWhom = 1

    if (iGaussian != 0)
    {
        //for linear iteration, the pure non-linear part equals to zero 
        gaussianBoundaryNonlinear.assign(iDielecBndyOdd, 0.0);
        gaussianChargeNonlinear.assign(iCrgedGridSum, 0.0);

        //for linear iteration, there is no ion. Thus we won't need a density here
        //if(fabs(fIonStrength)<fZero) gaussianBoundaryDensity.assign(iDielecBndyOdd, 1.0);
    }

    if (debug_solver) 
    {
        cout << "gaussianBoundaryDielec.size= "  << gaussianBoundaryDielec.size()  << endl;
        cout << "gaussianBoundaryDensity.size= " << gaussianBoundaryDensity.size() << endl;
        cout << "gaussianChargeDielec.size= "    << gaussianChargeDielec.size()    << endl;
        cout << "gaussianChargeDensity.size= "   << gaussianChargeDensity.size()   << endl;
        cout << "iDielecBndyEven= " << iDielecBndyEven << endl;
        cout << "iDielecBndyOdd= " << iDielecBndyOdd << endl;
        cout << "iCrgedGridSum= " << iCrgedGridSum << endl;
    }

    itr = 1;
    ires = 0;
    do 
    {
        rmsch = 0.0;
        rmxch = 0.0;

        /*
         * iterate over odd points
         */
        #ifndef PARALLEL
        SOR_itrOddPoints(1, itr); // forWhom = 1, single CPU
        #else
        omp_SOR_itrOddPoints(1, itr); // forWhom = 1, multi CPU
        #endif
        
        if (bFixedRelaxParam) 
        {
            int itr2 = 2 * itr - 1;
            om3 = 1.0 / (1.0 - om2 * fSpec * 0.25);
            if (fZero > om1) om3 = 1.0 / (1.0 - om2 * fSpec * 0.5);
            om4 = om3 / om2;
            om2 = om3;
            om1 = 1.0 - om2;

            if (0.0 < fIonStrength) 
            {
                if (1 == itr2 % 2) 
                {
                    for (vector<delphi_real>::iterator it = prgfSaltMap1.begin(); it != prgfSaltMap1.end(); ++it)
                        *it = (*it) * om4;
                } 
                else 
                {
                    for (vector<delphi_real>::iterator it = prgfSaltMap2.begin(); it != prgfSaltMap2.end(); ++it)
                        *it = (*it) * om4;
                }
            }

            for (vector<delphi_real>::iterator it = prgfCrgValA.begin(); it != prgfCrgValA.end(); ++it)
                *it = (*it) * om4;

            for (delphi_integer iy = 0; iy < iDielecBndyOdd; iy++)
                for (delphi_integer ix = 0; ix < 6; ix++)
                    prgfBndyDielec[iy][ix] = prgfBndyDielec[iy][ix] * om4;

            sixth = sixth * om4;
        }
        
        /*
         * Next update phimap2 using the new phimap1
         */
        #ifndef PARALLEL
        SOR_itrEvenPoints(1, itr); // forWhom = 1, single CPU
        #else
        omp_SOR_itrEvenPoints(1, itr); // forWhom = 1, multi CPU
        #endif
        
        if (bFixedRelaxParam) 
        {
            int itr2 = 2 * itr;
            om3 = 1.0 / (1.0 - om2 * fSpec * 0.25);
            if (fZero > om1) om3 = 1.0 / (1.0 - om2 * fSpec * 0.5);
            om4 = om3 / om2;
            om2 = om3;
            om1 = 1.0 - om2;

            if (0.0 < fIonStrength) 
            {
                if (1 == itr2 % 2) 
                {
                    for (vector<delphi_real>::iterator it = prgfSaltMap1.begin(); it != prgfSaltMap1.end(); ++it)
                        *it = (*it) * om4;
                } 
                else 
                {
                    for (vector<delphi_real>::iterator it = prgfSaltMap2.begin(); it != prgfSaltMap2.end(); ++it)
                        *it = (*it) * om4;
                }
            }

            for (vector<delphi_real>::iterator it = prgfCrgValA.begin(); it != prgfCrgValA.end(); ++it)
                *it = (*it) * om4;

            for (delphi_integer iy = 0; iy < iDielecBndyOdd; iy++)
                for (delphi_integer ix = 0; ix < 6; ix++)
                    prgfBndyDielec[iy][ix] = prgfBndyDielec[iy][ix] * om4;

            sixth = sixth * om4;
        }

        #ifdef DEBUG_DELPHI_SOLVER_ITIT
        if (1 == itr)
        {
            string strTestFile = "test_itit.dat";
            ofstream ofTestStream(strTestFile.c_str());
            ofTestStream << boolalpha;
            ofTestStream << fixed << setprecision(7);

            ix = 1;
            for (vector<delphi_real>::iterator it = phimap1.begin(); it != phimap1.end(); ++it)
            {
                ofTestStream << "phimap1[" << setw(6) << right << ix << "] = " << setw(11) << right << *it << endl;
                ix++;
            }

            ix = 1;
            for (vector<delphi_real>::iterator it = phimap2.begin(); it != phimap2.end(); ++it)
            {
                ofTestStream << "phimap2[" << setw(6) << right << ix << "] = " << setw(11) << right << *it << endl;
                ix++;
            }

            ofTestStream.close();
        }
        #endif // DEBUG_DELPHI_SOLVER_ITIT
        
        /*
         * we also save time by only checking convergence every 10 iterations, rather than every single iteration.
         * store phi2 in phi3 to compare against next iteration
         */
        if (iIterateInterval - 1 == itr % iIterateInterval) // itr = 9,19,29,...
        {
            for (ix = 1; ix < iHalfGridNum; ix += iConvergeFract) 
                prgfPhiMap[ix] = phimap2[ix];
        }

        if (0.0 < fGridConverge) 
        {
            grden = 0.0;

            for (ix = 0; ix < iCrgedGridEven; ix++) 
            {
                iy     = prgiCrgPose[ix];
                grden += phimap1[iy - 1] * prgfCrgValG[ix];
            }

            for (ix = iCrgedGridEven; ix < iCrgedGridSum; ix++) 
            {
                iy     = prgiCrgPose[ix];
                grden += phimap2[iy - 1] * prgfCrgValG[ix];
            }

            grdn[itr % 5] = grden / 2.0;
            if (10 < itr) 
            {
                bool igt = true;
                
                for (int i = 0; i < 5; i++)
                    for (int j = 0; j < 5; j++)
                        if (abs(grdn[j] - grdn[i]) > fGridConverge) igt = false;
                
                if (igt) 
                {
                    cout << grdn[0] << " " << grdn[1] << " " << grdn[2] << " " << grdn[3] << " " << grdn[4] << endl;
                    ires = 1;
                }
            }
        }

        if (0 == itr % iIterateInterval || 1 == ires) //----- check to see if accuracy is sufficient
        {            
            delphi_real rnorm2 = 0.0, temp2;

            for (ix = 1; ix < iHalfGridNum; ix += iConvergeFract) 
            {
                temp2 = prgfPhiMap[ix] - phimap2[ix];
                rnorm2 += temp2 * temp2;

                //if(rmxch<abs(temp2)) flag_itr=ix; //Lin Li: to see which grid is the key grid
                rmxch = max(rmxch, abs(temp2));
            }

            rmsch = sqrt((delphi_real) iConvergeFract * rnorm2 / ((iGrid - 2) * (iGrid - 2) * (iGrid - 2)));
            //rnormch = sqrt(rnorm2);
            rmsch2 = rmsch;
            rmxch2 = rmxch;

            if (0.0 < fGridConverge)
                cout << "        " << scientific << rmsch2 << "  " << rmxch2 << "  " << grden << "  at  " 
                     << setw(5) << left << itr << " iterations\n";
            else
                cout << "        " << scientific << rmsch2 << "  " << rmxch2 << "  at  " << setw(5) 
                     << left << itr << " iterations\n";

            if (fRmsc > rmsch || fMaxc > rmxch) ires = 1;

            if (bLogGraph) 
            {
                int ibin;
                for (int j = itr - 9; j <= itr; j++) 
                {
                    ibin = (j - 1) * (60 - 1) / (iLinIterateNum - 1) + 1;
                    rmsl[ibin - 1] = rmsch;
                    rmaxl[ibin - 1] = rmxch;
                }
            }

//            //if (9 == itr)
//            {
//                string strTestFile = "rank1_solver_itit.dat";
//                ofstream ofTestStream(strTestFile.c_str());
//                ofTestStream << boolalpha;
//                ofTestStream << fixed << setprecision(7);
//
//                ix = 0;
//                for (vector<delphi_real>::iterator it = phimap1.begin(); it != phimap1.end(); ++it)
//                {
//                    ofTestStream << "mpi_phimap1[" << setw(6) << right << ix << "] = " << setw(11) << right << *it << endl;
//                    ix++;
//                }
//
//                ix = 0;
//                for (vector<delphi_real>::iterator it = phimap2.begin(); it != phimap2.end(); ++it)
//                {
//                    ofTestStream << "mpi_phimap2[" << setw(6) << right << ix << "] = " << setw(11) << right << *it << endl;
//                    ix++;
//                }
//
//                ofTestStream.close();
//                
//                return;
//            }
        } 
        
        itr++;
        
        /*
         * check to see if accuracy is sufficient
         */
        if (1.0e-7 > maxres) 
        {
            if (iLinIterateNum >= itr && 0 == ires)
                continue;
            else
                break;
        } 
        else 
        {
            if ((iLinIterateNum >= itr || bAutoConverge) && (0 == ires))
                continue;
            else
                break;
        }
        
    } while (true);  
    
    //Argo: Printing the last rms-change of the iteration
    //Needed for GAUSSIAN runs since many were found to diverge
    cout << strLine60 << endl;
    cout << left << " (" << strNumSolver << ")" << infoString << "Iteration ended with final rms-change of " << scientific << rmsch2 << endl;
    
    if (rmsch2 > fMaxc) 
    {
        //cout << " WARNING !!! Run probably DIVERGED" << endl;
        //cout << " WARNING !!! Final rms-change is " << rmsch2 << " greater than " << fMaxc << endl;
        //cout << " WARNING !!! Try increasing the scale or lower the maxc value or both" << endl;
        stringstream str_Diverged;
        str_Diverged << "Run probably DIVERGED! Final rms-change is " << rmsch2 << " greater than " << fMaxc << endl;
        str_Diverged << "             Try increasing the scale or lower the maxc value or both" << endl;
        CDiverged warn_div(str_Diverged);
    } 
    else 
    {
        cout << left << " (" << strNumSolver << ")" << infoString << "Run converged in the provided limits" << endl;
    }

    postItr(rmaxl, rmsl);

    /*
     * code phimap corner, for use in transference from irises to convex and via versa
     */
    {
        delphi_real ap1, ap2, ap3, ap4;
        ap1 = prgfPhiMap[0];
        ap2 = ap1 * 10000.0;
        ap3 = (int) ap2;
        if (0 < ap3)
            ap4 = (ap3 + 0.8) / 10000.0;
        else
            ap4 = (ap3 - 0.8) / 10000.0;
        prgfPhiMap[0] = ap4;
    }

    #ifdef DEBUG_DELPHI_SOLVER_ITIT
    {
        string strTestFile = "test_itit.dat";
        ofstream ofTestStream(strTestFile.c_str());
        ofTestStream << boolalpha;
        ofTestStream << fixed << setprecision(7);

        const delphi_real *** phimap = pdc->getKey_constPtr<delphi_real>("phimap",iGrid,iGrid,iGrid); // const pointer to 3D phimap
        for (iz = 0; iz < iGrid; iz++)
            for (iy = 0; iy < iGrid; iy++)
                for (ix = 0; ix < iGrid; ix++)
                    ofTestStream << "phimap[" << setw(6) << right << iz+1 << "," << setw(6) 
                                 << right << iy+1 << "," << setw(6) << right << ix+1 << "] = " 
                                 << setw(12) << right << phimap[iz][iy][ix] << endl;

        ofTestStream.close();
    }
    #endif // DEBUG_DELPHI_SOLVER_ITIT

}
