/*
 * NWT_nitit.cpp
 *
 *  Created on: Aug 09, 2019
 *      Author: chuan
 *
 *  Remark: The nature of NWT iteration is more similar to that of the linear SOR so that it is constructed
 *  based on SOR_itit.
 */

#include "solver.h"

void CDelphiSolver::NWT()
{
    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 = " ----------------------------------------------------------------";

    iIterateInterval = 10;
    iConvergeFract   = 1;

    /*
     * NWT method uses the same initialization as that used in the delphi method
     */
    initOddEvenItr(1); // forWhom = 1

    cout << strLine60 << endl;
    cout << "      " << " rms-change   max change       #iterations" << endl;
    cout << strLine60 << endl;

    itr  = 1;
    ires = 0;

    do
    {
        rmsch = 0.0;
        rmxch = 0.0;

        /*
         * iterate over odd points
         */
        NWT_itrOddPoints();

        /*
         * Next update phimap2 using the new phimap1
         */
        NWT_itrEvenPoints();

        /*
         * test code
         */
//        if (1 == itr)
//        {
//            string strTestFile = "test_NWT.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();
//        }

        /*
         * 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];
        }

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

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

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

            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) / (_ntIterateNum - 1) + 1;
                    rmsl[ibin - 1]  = rmsch;
                    rmaxl[ibin - 1] = rmxch;
                }
            }

        } //----- end of convergence check

        itr++;

        /*
         * check to see if accuracy is sufficient
         */
        if (1.0e-7 > maxres)
        {
            if (_ntIterateNum >= itr && 0 == ires)
                continue;
            else
                break;
        }
        else
        {
            if ((_ntIterateNum >= itr || bAutoConverge) && (0 == ires))
                continue;
            else
                break;
        }

    } while (true);

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