/*
 * mpi_solver_run.cpp
 *
 *  Created on: Aug 26, 2019
 *      Author: chuanli
 */

#include "solver.h"

/*
 * implementation of pure abstract function mpi_run() in class IAbstractModule
 */
void CDelphiSolver::mpi_run()
{
    int mpi_modflag = 0, mpi_slvflag = 0;
    /*
     * master process carries out necessary preparations for iterations
     */
    if (0 == mpi_rank)
    {
       /*
        * 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;
           }
       }

       if (bIsNPBE) mpi_modflag = 1;
       if ( bIsNWT) mpi_slvflag = 1;
    }

    MPI_Bcast(&mpi_modflag, 1, MPI_INT, 0, MPI_COMM_WORLD); // mpi solver flag
    MPI_Bcast(&mpi_slvflag, 1, MPI_INT, 0, MPI_COMM_WORLD); // mpi solver flag

    MPI_Barrier(MPI_COMM_WORLD);

    if (0 == mpi_modflag)
    {
        bIsLPBE = true;
        bIsNPBE = false;
    }

    if (1 == mpi_modflag)
    {
        bIsLPBE = false;
        bIsNPBE = true;
    }

    if (0 == mpi_slvflag)
    {
        bIsSOR = true;
        bIsNWT = false;
    }

    if (1 == mpi_slvflag)
    {
        bIsSOR = false;
        bIsNWT = true;
    }

    if (bIsNPBE && bIsNWT)
    {
        if (0 == mpi_rank)
            cout << " NO MPI-PARALLELIZED SOLVER = (" << strNumSolver <<"). CHOOSE ANOTHER ONE TO TRY... \n\n" << endl;
        MPI_Finalize();
        exit(1);
    }

    //-------------------------------------------------------------------------------------------------------------------------------

    int noit;

    if (inhomo == 1) fIonStrength = 0;

   /*
    * master process carries out necessary preparations for iterations
    */
   if (0 == mpi_rank)
   {
      debug_solver = false;
      const char * infoString = " Info> ";
      const char * timeString = " Time> ";
      const char * bndcString = " Bndc> ";
      size_t MAXWIDTH = 45;

      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;

      validateInput();

      /*
       * master syncs prgbDielecMap (a.k.a idebmap(igrid,igrid,igrid)) from global data container
       */
      pdc->readGlobalVector1D<char>("idebmap", 0, iGrid * iGrid * iGrid, prgbDielecMap);

      mpi_setDielecBndySaltMap();

      mpi_setCrg();

      /*
       * master frees prgbDielecMap (a.k.a idebmap(igrid,igrid,igrid)) in local data container
       * when running linear solver
       */
      if (0 == iNonIterateNum || fZero > fIonStrength || inhomo == 1)
          if (prgbDielecMap.size() > 0) vector<char>().swap(prgbDielecMap);

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

      /*
       * MPI implementation does not output dielectric maps
       */
//      if (bEpsOut && iGaussian == 0) //----- write dielectric map
//      {
//         unique_ptr<CIO> pio(new CIO()); // smart unique_ptr
//         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->writeGaussEpsMap(iGrid, repsout, fScale, fgBoxCenter, gepsmp2, strEpsFile);
//         pio.reset();
//      }

      //----- the order of SOR_calculateRelaxFactor and setBndy cannot 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 (bSpectralRadius)
      {
         fSpec = fSpectralRadius;
         cout << infoString << left << setw(MAXWIDTH) << "Using entered value for relaxation of" << " : " << fSpec << endl;
      }
      else
      {
         fSpec = SOR_calculateRelaxFactor();
      }

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

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

      if (bAutoConverge) iLinIterateNum = noit;

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

      /*
       * master syncs prggvCrgedAtom (a.k.a atmcrg(nqass)) from global data container if iBndyType = 4
       */
      if ( (7 == iBndyType && !(bIonsEng && 0 == iNonIterateNum && fZero < abs(fNetCrg)) ) || 4 == iBndyType )
          pdc->readGlobalVector1D< SGridValue<delphi_real> >("atmcrg", 0, iCrgGridNum, prggvCrgedAtom);

      setBndy();

      /*
       * master frees prggvCrgedAtom (a.k.a atmcrg(nqass)) in local data container to save memory usage
       */
      if (prggvCrgedAtom.size() > 0) vector< SGridValue<delphi_real> >().swap(prggvCrgedAtom);

      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 (0 == iNonIterateNum || fZero > fIonStrength || inhomo == 1)
      {
    	  	  if (0 >= iLinIterateNum) throw CZeorLinIterateNum(bAutoConverge,iLinIterateNum);
      }
      else
      {
          if (50 < noit || 0 >= iLinIterateNum)
          {
             iLinIterateNum = noit / 2;
             cout << left << " (" << strNumSolver << ") Re-estimated iterations now " << " : " << iLinIterateNum << endl;
          }
      }
   } // ----- end of if (0 == mpi_rank())

   /*
    * ---------- the following lines are performed on both master and salve processes ----------
    */

   if (iGaussian == 1) // Gaussian based runs
   {
       /*
        * Gaussian is not parallelized yet!
        */
       if (0 == mpi_rank)
       {
           if (0 == iNonIterateNum || fZero > fIonStrength || inhomo == 1)
               SOR_itit();
           else
           {
               delphi_real qfact = abs(fNetCrg) * (delphi_real) iCrgBndyGridNum / iBndyGridNum;
               SOR_nitit(qfact);
           }
       }
   }
   else // regular runs
   {
	   /*
	    * data container variables not to be changed
	    */
	   MPI_Bcast(&fIonStrength,     1, mpi_delphi_real, 0, MPI_COMM_WORLD); // rionst    (value from data container)

	   /*
        * reference to read-and-write variables from data container
        */
	   MPI_Bcast(&iGrid,            1, MPI_INT,         0, MPI_COMM_WORLD); // igrid     (reference to data container)
	   MPI_Bcast(&iBndyType,        1, MPI_INT,         0, MPI_COMM_WORLD); // ibctyp    (reference to data container)
	   MPI_Bcast(&fNetCrg,          1, mpi_delphi_real, 0, MPI_COMM_WORLD); // qnet      (reference to data container)
	   MPI_Bcast(&iLinIterateNum,   1, MPI_INT,         0, MPI_COMM_WORLD); // nlit      (reference to data container)
	   MPI_Bcast(&iIterateInterval, 1, MPI_INT,         0, MPI_COMM_WORLD); // icon1     (reference to data container)
	   MPI_Bcast(&iConvergeFract,   1, MPI_INT,         0, MPI_COMM_WORLD); // icon2     (reference to data container)
	   MPI_Bcast(&fRelaxParam,      1, mpi_delphi_real, 0, MPI_COMM_WORLD); // relpar    (reference to data container)

	   /*
	    * out into data container
	    */
	   MPI_Bcast(&iDielecBndyOdd,   1, MPI_INT,         0, MPI_COMM_WORLD); // icount2b  (reference to data container)
	   MPI_Bcast(&iCrgedGridSum,    1, MPI_INT,         0, MPI_COMM_WORLD); // icount1b  (reference to data container)
	   MPI_Bcast(&iCrgBndyGridNum,  1, MPI_INT,         0, MPI_COMM_WORLD); // ibc       (reference to data container)

	   /*
	    * class variables
	    */
	   MPI_Bcast(&fSpec,            1, mpi_delphi_real, 0, MPI_COMM_WORLD); // fSpec     (local value)
	   MPI_Bcast(&iDielecBndyEven,  1, MPI_INT,         0, MPI_COMM_WORLD); // icount2a  (local value)
	   MPI_Bcast(&iCrgedGridEven,   1, MPI_INT,         0, MPI_COMM_WORLD); // icount1a  (local value)

	   /*
	    * local function variables
	    */
	   MPI_Bcast(&noit,            1, MPI_INT,         0, MPI_COMM_WORLD); // noit      (local value)

       if (0 == iNonIterateNum || fZero > fIonStrength || inhomo == 1)
       {
          mpi_SOR_itit();
       }
       else
       {
          delphi_real qfact = abs(fNetCrg)*(delphi_real)iCrgBndyGridNum/iBndyGridNum;
          mpi_SOR_nitit(qfact);
       }
   }
}

