/*
 * solver_run.cpp
 *
 *  Created on: Feb 10, 2014
 *  Updated on: Aug 09, 2019 to include NWT method
 *      Author: chuan
 */

#include "solver.h"

void CDelphiSolver::run()
{
    debug_solver = false;
    const char * infoString = "Info> ";
    const char * timeString = "Time> ";
    const char * bndcString = "Bndc> ";
    size_t MAXWIDTH = 45;

    int noit;

    /*
     * select appropriate solver to use
     */
    if ( 0 == strNumSolver.compare("SOR") ) // SOR solver
    {
        bIsSOR  = true;
        bIsNWT  = false;

        if (0 == iNonIterateNum || fZero > fIonStrength || inhomo == 1) // LPBE
        {
            bIsLPBE = true;
            bIsNPBE = false;

            cout << endl;
            cout << left << " DelPhi solves the (LPBE) model with the (" << strNumSolver << ") method now ..." << endl << endl;
        }
        else
        {
            bIsLPBE = false;
            bIsNPBE = true;

            cout << endl;
            cout << left << " DelPhi solves the (NLPBE) model with the (" << strNumSolver << ") method now ..." << endl << endl;
        }
    }

    if ( 0 == strNumSolver.compare("NWT") ) // NWT solver
    {
        if (0 == iNonIterateNum || fZero > fIonStrength || inhomo == 1) // LPBE
        {

            bIsLPBE    = true;
            bIsNPBE    = false;
            bIsSOR     = true;
            bIsNWT     = false;

            cout << endl;
            cout << left << " DelPhi solves the (LPBE) model with the (" << strNumSolver << ") method, which is identical to the (SOR) method, now ..." << endl << endl;
        }
        else  // NWT solver and NPBE
        {
            bIsLPBE = false;
            bIsNPBE = true;
            bIsSOR  = false;
            bIsNWT  = true;

            cout << endl;
            cout << left << " DelPhi solves the (NLPBE) model with the (" << strNumSolver << ") method now ..." << endl << endl;
        }
    }

    validateInput();

    if (inhomo == 1) fIonStrength = 0;

    if (iGaussian == 1 && inhomo == 0) //reset some vectors for 2nd run of Gaussian
    {
        if (prgfGridCrg.size() > 0) 
        {
            if (debug_solver) cout << "prgfGridCrg.size(): " << prgfGridCrg.size() << endl;
            vector<delphi_real>().swap(prgfGridCrg); //gchrg
        }
        if (prgigGridCrgPose.size() > 0) 
        {
            if (debug_solver) cout << "prgigGridCrgPose.size(): " << prgigGridCrgPose.size() << endl;
            vector<SGrid<delphi_integer> >().swap(prgigGridCrgPose); //gchrgp
        }

        //vector<delphi_integer> ibndx,ibndy,ibndz;
        if (ibndx.size() > 0) 
        {
            if (debug_solver) cout << "ibndx.size(): " << ibndx.size() << endl;
            vector<delphi_integer>().swap(ibndx); //ibndx
        }
        if (ibndy.size() > 0) 
        {
            if (debug_solver) cout << "ibndy.size(): " << ibndy.size() << endl;
            vector<delphi_integer>().swap(ibndy); //ibndy
        }
        if (ibndz.size() > 0) 
        {
            if (debug_solver) cout << "ibndz.size(): " << ibndz.size() << endl;
            vector<delphi_integer>().swap(ibndz); //ibndz
        }
        if (phimap1.size() > 0) 
        {
            if (debug_solver) cout << "phimap1.size(): " << phimap1.size() << endl;
            vector<delphi_real>().swap(phimap1); //phimap1
        }
        if (phimap2.size() > 0) 
        {
            if (debug_solver) cout << "phimap2.size(): " << phimap2.size() << endl;
            vector<delphi_real>().swap(phimap2); //phimap2
        }

        if (prgiBndyDielecIndex.size() > 0) 
        {
            if (debug_solver) cout << "prgiBndyDielecIndex.size(): " << prgiBndyDielecIndex.size() << endl;
            vector<delphi_integer>().swap(prgiBndyDielecIndex); //prgiBndyDielecIndex
        }
        if (prgfBndyDielec.size() > 0) 
        {
            if (debug_solver) cout << "prgfBndyDielec.size(): " << prgfBndyDielec.size() << endl;
            vector<vector<delphi_real> >().swap(prgfBndyDielec); //prgfBndyDielec
        }
        //vector< SDoubleGridValue >& prgdgvCrgBndyGrid;   // cgbp(ibc)
        if (prgdgvCrgBndyGrid.size() > 0) 
        {
            if (debug_solver) cout << "prgdgvCrgBndyGrid.size(): " << prgdgvCrgBndyGrid.size() << endl;
            vector<SDoubleGridValue>().swap(prgdgvCrgBndyGrid); //prgdgvCrgBndyGrid
        }
        if (debug_solver) cout << "prgfPhiMap.size(): " << prgfPhiMap.size() << endl;

        if (gaussianBoundaryDensity.size() > 0) 
        {
            if (debug_solver) cout << "gaussianBoundaryDensity.size(): " << gaussianBoundaryDensity.size() << endl;
            vector<delphi_real>().swap(gaussianBoundaryDensity); //prgdgvCrgBndyGrid
        }

        iDielecBndyOdd = 0;
    }
    
    if (debug_solver) cout << "iDielecBndyOdd: " << iDielecBndyOdd << endl;

	/*
	 * Both SOR and NWT methods use the same routines to set boundary/salt maps and charges in order
	 * to minimize the impact of existing code. However, the constant coefficients will be alternated
	 * in the NWT method.
	 */
	setDielecBndySaltMap();
	setCrg();

	cout << left << " (" << strNumSolver << ") " << timeString << "iepsmp to db, and charging done on ";
	pTimer->showTime();
	cout << endl;
	cout << left << " (" << strNumSolver << ") Number of grid points assigned charge" << " : " << iCrgedGridSum << endl;

	if (bEpsOut && iGaussian == 0) //----- write dielectric map
	{
		unique_ptr < CIO > pio(new CIO()); // smart unique_ptr
		// pio->writeEpsMap(iAtomNum,iObjectNum,iGrid,fScale,fgBoxCenter,prgigEpsMap,prgbDielecMap,strEpsFile);
		pio->writeHomoEpsMap(iGrid, repsout, repsin, fScale, fgBoxCenter, prgbDielecMap, strEpsFile);
		pio.reset();
	}

	if (bEpsOut && iGaussian == 1) //----- write dielectric map
	{
		unique_ptr < CIO > pio(new CIO()); // smart unique_ptr
		// pio->writeEpsMap(iAtomNum,iObjectNum,iGrid,fScale,fgBoxCenter,prgigEpsMap,prgbDielecMap,strEpsFile);
		pio->writeGaussEpsMap(iGrid, repsout, fScale, fgBoxCenter, gepsmp2, strEpsFile);
		pio.reset();
	}

	//----- the order of SOR_calculateRelaxFactor and setBndy cannnot be reverted! prgfPhiMap is used
	//      as temporary array in SOR_calculateRelaxFactor...

	phimap1.assign((iGrid*iGrid*iGrid + 1) / 2, 0.0); //phimap1.assign(icgrid/2,0.0);
	phimap2.assign((iGrid*iGrid*iGrid + 1) / 2, 0.0); //phimap2.assign(icgrid/2,0.0);

    if ( bIsSOR ) // NWT solver and LPBE
    {
		if (bSpectralRadius)
		{
			fSpec = fSpectralRadius;
			cout << infoString << left << " (" << strNumSolver << ") Using entered value for relaxation of" << " : " << fSpec << endl;
		}
		else
		{
			fSpec = SOR_calculateRelaxFactor();
		}

		noit = (int) (7.8 / log(1.0 + sqrt(1 - fSpec)));

		cout << left << " (" << strNumSolver << ") Estimated iterations to convergence" << " : " << noit << endl;

		if (bAutoConverge) iLinIterateNum = noit;
    }
    else if ( bIsNWT ) // NWT solver and NPBE
    {
    	_ntIterateNum = iNonIterateNum;

        if (bManualRelaxParam)
        {
            om2 = fRelaxParam;
        }
        else
        {
        	//fSpec = 0.9982;
        	//om2   = 2.0 / (1.0 + sqrt(1.0 - fSpec));

        	om2 = 1.9186758;
        }
        cout << left << " (" << strNumSolver << ") Input fixed relaxation parameter" << " : " << om2 << endl;

    	om1   = 1.0 - om2;

    	cout << left << " (" << strNumSolver << ") Estimated iterations to convergence" << " : " << _ntIterateNum << endl;
    }

	/*
	 * setBndy() requires using phimap. move allocating phimap here to save memory usage
	 */
	prgfPhiMap.assign(iGrid * iGrid * iGrid, 0.0);

	setBndy();

	cout << endl;
	cout << left << " (" << strNumSolver << ") " << timeString << "Setup time done on ";
	pTimer->showTime();
	cout << endl;

	if (0.3 < (delphi_real) iCrgBndyGridNum / iBndyGridNum && bAutoConverge)
	{
		iLinIterateNum = iLinIterateNum * iCrgBndyGridNum / (0.3 * iBndyGridNum);
		cout << left << setw(MAXWIDTH) << " Re-estimated iterations now " << " : " << iLinIterateNum << endl;
	}

	cout << left << " (" << strNumSolver << ") " << timeString << "Now iterating on ";
	pTimer->showTime();
	cout << endl;

    if (bIsSOR) // SOR solver
    {
		if (bIsLPBE)
		{
		    if (0 >= iLinIterateNum) throw CZeorLinIterateNum(bAutoConverge, iLinIterateNum);
			SOR_itit();
		}
		else
		{
		    if (50 < noit || 0 >= iLinIterateNum)
			{
				iLinIterateNum = noit / 2;
				cout << left << left << " (" << strNumSolver << ") Re-estimated iterations now " << " : " << iLinIterateNum << endl;
			}

			delphi_real qfact = abs(fNetCrg) * (delphi_real) iCrgBndyGridNum / iBndyGridNum;

			SOR_nitit(qfact);
		}
    }

    if (bIsNWT) // NWT solver
    {
        NWT();
    }
}

