//
//  slbmpi.h
//  MHDA
//
//  Created by Dr. Jacques C. Richard on 1/31/12.
//

#ifndef _slbmpi_h
#define _slbmpi_h

#include <boost/iterator/indirect_iterator.hpp>
#include <boost/iterator/transform_iterator.hpp>

const double c1_3  = 1.e0/3.e0;
const double c1_18 = 1.e0/18.e0;
const double c1_36 = 1.e0/36.e0;

const double w_epsi=3.0e0;
const double w_epsi_j=-11.0e0/2.0e0;
const double w_xx=-1.e0/2.0e0;    

//const double c4=1.e0/4.e0; //JCR: un-used
const double c8=1.e0/8.e0;
const double c10=1.e0/10.0e0;
const double c12=1.e0/12.0;
const double c18=1.e0/18.0;
const double c19=1.e0/19.0;
const double c21=1.e0/21.0;
const double c23=2.0e0/3.0;
const double c24=1.e0/24.0;
const double c36=1.e0/36.0;
const double c40=1.e0/40.0;
const double c63=1.e0/63.0;
const double c252=1.e0/252.0;
const double c399=5.e0/399.0;
const double c1197=4.e0/1197.0;
const double c2394=11.e0/2394.0;

const double rh0=1.0e0; 

//KGK Function to set lattice ID for a lattice
struct lattice_ID
{
        int mm,nn,rr;
        const int MM,NN,RR;
        lattice_ID (int& mm1, int& nn1, int& rr1, int& MM1, int& NN1, int& RR1) : mm(mm1),nn(nn1),rr(rr1),MM(MM1),NN(NN1),RR(RR1) {}

    void operator()(lattice_type& lattice)
    {
        if (rr > RR) 
        {
            rr = 1; nn++;
            if (nn > NN) 
            {
                nn = 1; mm++;
            }
        }
        lattice.rr = rr++; lattice.nn = nn; lattice.mm = mm;
        if (mm > MM || nn > NN) std::cout << " Error in lattice_ID in math_functions.h " << std::endl;
    }
};

//KGK Interface function to set lattice ID for all lattices in the grid
void set_lattice_ID(const xyz_type<int> PhysicalSpaceSize, grid_type& grid)
{
    int MM = PhysicalSpaceSize.x+2;
    int NN = PhysicalSpaceSize.y+2;
    int RR = PhysicalSpaceSize.z+2;
    int mm=1;
    int nn=1;
    int rr=1;
    
    std::for_each(grid.begin(), grid.end(), lattice_ID(mm, nn, rr, MM, NN, RR));
};

template<typename T>
struct init_internal_neighbors_wf
{
#ifdef BOOST_MPI_HPP
    const boost::mpi::communicator &world ;
    init_internal_neighbors_wf( const boost::mpi::communicator &W ) : world( W ) { }
#endif
    void operator()(T& lattice) const
    { //JCR: NT, for some reason, used (x,y,z)_offset[1..en] = - (cx, cy, cz)[1..en] here?! Is that to use "+" in idx._ + _offset?
        const int x_offset[ en - 1 ] = {-1,  1,  0,  0,  0,  0, -1,  1, -1, 1, -1,  1, -1,  1,  0,  0,  0,  0};
        const int y_offset[ en - 1 ] = { 0,  0, -1,  1,  0,  0, -1, -1,  1, 1,  0,  0,  0,  0, -1,  1, -1,  1};
        const int z_offset[ en - 1 ] = { 0,  0,  0,  0, -1,  1,  0,  0,  0, 0, -1, -1,  1,  1, -1, -1,  1,  1};
        
        const xyz_type<int>& idx  = lattice.idx;
        const xyz_type<int>& dims = lattice.dims;
        
        for (int i=0; i < en - 1; i++)
        {
#ifdef BOOST_MPI_HPP
            if ( world.size() > 1 ) //( ( idx > T::dims / world.size() ) || ( idx.x < 0 ) || ( idx.y < 0 ) || ( idx.z < 0 ) )
            {
                std::vector< boost::mpi::request > reqs ; 
                std::vector< T > Neighbor2Proc ;
                std::vector< T > Neighbor4Proc( en - 1 ) ;
                
                int Sender   = world.rank( ) ; // also = tag & dest & src
                int Receiver = Sender + 1 ; if(Receiver >= world.size())Receiver = Sender - 1 ;if(( Receiver<0)||(Receiver==Sender)) return ;
                int SendTag  = 0 ;
                int RecvTag  = 1 ;
                int NewOffset_x = 0 ; // ( idx.x > T::dims.x / world.size() ? abs( x_offset[i] ) : x_offset[i] ) ;
                int NewOffset_y = 0 ; // ( idx.y > T::dims.y / world.size() ? abs( y_offset[i] ) : y_offset[i] ) ;
                int NewOffset_z = 0 ; // ( idx.z > T::dims.z / world.size() ? abs( z_offset[i] ) : z_offset[i] ) ;
                Neighbor2Proc.push_back (*(T::begin_iterator + T::offset(idx.x+x_offset[i]-NewOffset_x, idx.y+y_offset[i]-NewOffset_y, idx.z+z_offset[i]-NewOffset_z)));
                //Neighbor2Proc[ Sender ] = (*(T::begin_iterator + T::offset(idx.x+x_offset[i]-NewOffset_x, idx.y+y_offset[i]-NewOffset_y, idx.z+z_offset[i]-NewOffset_z)));
                long int Msg2Send_size = Neighbor2Proc.size() ; // T::begin_iterator + ;
                Msg2Send_size += (idx.x+x_offset[i]-NewOffset_x)*(idx.y+y_offset[i]-NewOffset_y)*(idx.z+z_offset[i]-NewOffset_z) ;
                std::copy( Neighbor2Proc.begin() , Neighbor2Proc.end() , Neighbor4Proc.begin() ) ; // a dflt !?
#ifdef DEBUG
                std::cout<< "\ninit_internal_neighbors_wf: Process "<<world.rank()<<" of "<<world.size()<<" about to exchange (if necessary) w/+/-1!" 
                         << " Sender="<< Sender 
                         <<", Receiver="<< Receiver 
                         <<", Neighbor2Proc.size()="<<Neighbor2Proc.size()
                         <<", Msg2Send_size=" << Msg2Send_size
                         <<", i="<<i<<" @ idx="<<idx<<"\n" ;
                Pause( "init_internal_neighbors_wf: _slbmpi_h: 108: pre-exchange" ) ;
#endif
                if ( Msg2Send_size > 0 ) // if @ a location w/something 2 exchange?!
                {
                  std::cerr << "DEBUG: slbmpi.h:" << __LINE__  << std::endl << std::flush;
                  // ORIG: reqs.push_back( world.isend( Sender   , SendTag , &Neighbor2Proc[ i ] , Msg2Send_size ) ) ;
                  reqs.push_back( world.isend( Receiver   , SendTag , Neighbor2Proc[i] ) );
#ifdef DEBUG
                    std::cout<< "\ninit_internal_neighbors_wf: Process "<<world.rank()<<" of "<<world.size()<<" after msg sent!" 
                    << " Receiver="<< Receiver <<", Sender="<< Sender <<", SendTag="<< SendTag
                    <<", Msg2Send_size="<<Msg2Send_size<<", Neighbor2Proc.size()="<<Neighbor2Proc.size()<<", i="<<i<<" @ idx="<<idx<<"\n" ;
                    Pause( "init_internal_neighbors_wf: _slbmpi_h: 167: exchanged?" ) ;
#endif
                  std::cerr << "DEBUG: slbmpi.h:" << __LINE__  << std::endl << std::flush;
                    boost::mpi::status msg = world.probe() ;
                  std::cerr << "DEBUG: slbmpi.h:" << __LINE__  << std::endl << std::flush;
#ifdef DEBUG
                    if ( msg.error() ) {
                      std::cerr << "DEBUG: MPI error detected!" << std::endl << std::flush;
                    }
                    else {
                      std::cerr << "DEBUG: received msg of length " << msg.count<T>()
                                << " from " << msg.source() 
                                << " with tag " << msg.tag()
                                << std::endl << std::flush;
                    };
#endif
                    SendTag = ( SendTag == msg.source() ? msg.source() : SendTag ) ;
                    //SendTag = msg.source() ;
                    RecvTag = msg.tag() ;
                    
                    reqs.push_back( world.irecv( msg.source() , msg.tag() ,  Neighbor4Proc[ i ] ) ) ; //  , Msg2Send_size: not compiling:  request irecv(int source, int tag, T * values, int n) const; @ http://www.boost.org/doc/libs/1_48_0/doc/html/boost/mpi/communicator.html#id444949-bb
                    
#ifdef DEBUG
                    std::cout<< "\ninit_internal_neighbors_wf: Process "<<world.rank()<<" of "<<world.size()<<" after msg received!" 
                    << " Receiver="<< Receiver <<", Sender="<< Sender <<", RecvTag="<< RecvTag <<", Neighbor4Proc.size()="<<Neighbor4Proc.size()<<", i="<<i<<" @ idx="<<idx<<"\n" ;
                    Pause( "init_internal_neighbors_wf: _slbmpi_h: 178: exchanged!?" ) ;
#endif
                    for( unsigned int ir = 0 ; ir < reqs.size() ; ir++ ) if ( reqs[ ir ].test() ) reqs[ ir ].wait() ;
                    //      boost::mpi::wait_any(reqs.begin(), reqs.end()); // wait_all hangs/seems to wait forever; wait_any, or w/o anything here, segflt's
                    //     }
#ifdef DEBUG
                    std::cout<< "\ninit_internal_neighbors_wf: Process "<<world.rank()<<" of "<<world.size()<<" exchanging w/+1!" 
                    << " Receiver="<< Receiver <<", Sender="<< Sender <<"! Next: lattice.neighbors[i] = &(Neighbor4Proc[ Receiver ])\n" ;
                    Pause( "init_internal_neighbors_wf: _slbmpi_h: 191: exchanged!?!?" ) ;
#endif
                }
#ifdef DEBUG
                else {
                  std::cerr << "DEBUG: Rank " << world.rank() << " skipping send because Msg2Send_size==0" << std::endl << std::flush;
                };
#endif
                lattice.neighbors[i] = &( Neighbor4Proc[ i ] ) ;
#ifdef DEBUG
                std::cout<< "\ninit_internal_neighbors_wf: Process "<<world.rank()<<" of "<<world.size()<<" exchanging w/+1!" 
                << " Receiver="<< Receiver <<", Sender="<< Sender <<"! After: lattice.neighbors[i] = &(Neighbor4Proc[ Receiver ])\n" ;
                Pause( "init_internal_neighbors_wf: _slbmpi_h: 198: exchanged?" ) ;
#endif
            }
            else
#endif
                //KGK: For boundary lattices, definition of neighbors depends on the boundary conditions
                //KGK: This is taken care of in the boundary functions. Hence, set dummies (self) for boundary lattices
                //KGK: gcc is forgiving (not crashing) but this change is required for visual C++
                if ((idx.x+x_offset[i])>=0     && (idx.y+y_offset[i])>=0     && (idx.z+z_offset[i])>=0 && 
                    (idx.x+x_offset[i])<dims.x && (idx.y+y_offset[i])<dims.y && (idx.z+z_offset[i])<dims.z)
                    lattice.neighbors[i] = &(*(T::begin_iterator + T::offset(idx.x+x_offset[i], idx.y+y_offset[i], idx.z+z_offset[i])));
                else 
                    lattice.neighbors[i] = &(*(T::begin_iterator));
            //KGK Store the indices for the neighbors for positive identification/checking
            lattice.neighborID[i][0] = lattice.neighbors[i]->mm; lattice.neighborID[i][1] = lattice.neighbors[i]->nn; 
            lattice.neighborID[i][2] = lattice.neighbors[i]->rr;
            
        }
    }   
};


struct init_f
{
        const IBC_data &IBc_Data ;
        const int solver ;
        const int compres ;
        const std::vector<xyz_dt> &v, &b, &e ;
        const r_Space &t_0, &r_0 ;
        
        init_f(const IBC_data &IBc_Dat, const int Solver, const int Compres , 
           const std::vector<xyz_dt> &V, const std::vector<xyz_dt> &B, const std::vector<xyz_dt> &E, 
           const r_Space &T0, const r_Space &R0 
           ) : IBc_Data(IBc_Dat), solver(Solver), compres(Compres), v(V), b(B), e(E), t_0(T0), r_0(R0) 
        { }
    
    void operator()(lattice_type& lattice) const 
    {
        
        const int n = lattice.offset(lattice.idx.x, lattice.idx.y, lattice.idx.z);//cout<<" Checking: idx="<<lattice.idx", n="<<n;
#ifdef DEBUG
        if((mag(v[n])<Near0)||isnan(mag(v[n]))||isinf(mag(v[n]))||isnan(mag(e[n]))||isnan(mag(b[n]))||isinf(mag(b[n])))
        {
            cout<<"\nENTERING init_f: v[n="<<n<<"]="<<v[n]<<", b="<<b[n]<<", e="<<e[n]<<", u0="<<IBc_Data.u0
            <<", B0="<<IBc_Data.B0<<", E0="<<IBc_Data.E0<<" @ "<<lattice.idx;
            //       Pause( "ENTERING init_f: _slbmpi_h: 243: entered w/arrays==...?" ) ;
        }
#endif
        xyz_dt w        = v[   n ] ;                                              // passing what's initialized @ main
        if ( lattice.is_velinlet() ) lattice.u = IBc_Data.u0; else lattice.u=v[n];//initialized @ main; redundant? Not if v_init varies..
        if ( IBc_Data.when > 0 ) lattice.u = v[ n ] ;                             //if when=t4ReStart b4 init_f; OVERRIDES is_velinlet ck
        lattice.B       = b[   n ] ;
        lattice.E       = e[   n ] ;
        lattice.T       = t_0[ n ] ;
        lattice.rho     = r_0[ n ] ;
        lattice.A       = 0.0 ;                                                   //  MUST initialize
        lattice.J       = lattice.rho * lattice.u ;
        lattice.curlJ   = 0.0 ;
        lattice.gradrho = 0.0 ;
        lattice.Phi     = 0.0 ;
        {
            Normalizer   = 1.0 ;
            lattice.rho    = rh0;
            lattice.fi[0]  = c1_3*rh0;
            std::fill(lattice.fi + 1, lattice.fi + 7,  c1_18 * rh0);
            std::fill(lattice.fi + 7, lattice.fi + en, c1_36 * rh0);
        }

        lattice.T   = 0.0 ;                                                     // to get from f instead of pre-set
        for(int i=0; i<=en-1;i++)
        { 
            lattice.T += ( cix[i] * cix[i] + ciy[i] * ciy[i] + ciz[i] * ciz[i] ) * lattice.fi[i] / lattice.rho ; //really m(v_molec_av)^2/2 =(3/2)k_B T but try this
        } 
        
        lattice.p   = Fluid_Species::p_EqState(lattice.rho,lattice.T) ;         //actually p_d or p_u could be used but don't want to over-specify p(r,T)
        
        std::copy(lattice.fi, lattice.fi + en, lattice.old_fi);
#ifdef DEBUG
        //  if( n == 0 ) std::cout<<"\n @ all n for 7>Solver>2 else @ n=0=(i,j,k), setting Normalizer="<<Normalizer << endl ;
        //std::cout<<" ..init_f.. ";
        for (int i=0; i<en; i++) if( (fabs(lattice.fi[i])>1.0e10) || isinf(lattice.fi[i]) || isnan(lattice.fi[i]) || (fabs(lattice.fi[i])<Near0) || (lattice.rho>10.0) || (lattice.rho<Near0) )
        {
            cout<<"init_f: lattice.fi[i]= "<<lattice.fi[i]<<" w="<<w<<" & lat.rh="<<lattice.rho<<" @ "<<lattice.idx<<" @ i="<<i<<endl;
            //       Pause( "init_f: _slbmpi_h: 281: lattice.f" ) ;
        } 
#endif
    }
};


struct set_old_fi_wf
{
    void operator()(lattice_type& lattice) const 
    {
        std::copy(lattice.fi, lattice.fi + en, lattice.old_fi);
    }
};


struct physical_wf
{
    static double ux_max;
    
    const xyz_dt u_i ; 
    const Collision_Integrator Coll_Integrate ;
    const IBC_data IBC_D ;
    
    physical_wf(const xyz_dt& U_i, const Collision_Integrator Collide_Integrate , const IBC_data IBC_d) : 
    u_i(U_i), Coll_Integrate(Collide_Integrate), IBC_D(IBC_d) { }
    
    void operator()(lattice_type& lattice) const
    {
        xyz_dt& u          = lattice.u;
        const PDF_type& fi = lattice.fi;
        
#ifdef DEBUG
        if ((isnan(lattice.rho))||(lattice.rho>10.0)||(lattice.rho<Near0))
        {
            cout<<"\nINTO physical_wf: u_i="<<u_i<<", u= "<<u<<", rh="<<lattice.rho<<", T="<<lattice.T<<", p="<<lattice.p<<" @ "<<lattice.idx;
            for (int i=0; i<en; i++) if( (fabs(lattice.fi[i])>1.0e10) || isinf(lattice.fi[i]) || isnan(lattice.fi[i]) || (fabs(lattice.fi[i])>1.0e1) )
            {
                cout<<" fi[i="<<i<<"]= "<<lattice.fi[i];
                Pause( "INTO physical_wf: _slbmpi_h: 320: lattice.f" ) ;
            }
            cout<<"\nINTO physical_wf: To sum &/or plot in Excel or Matlab:\n lattice.fi        fi   lattice.old_fi";
            for (int i=0; i<en; i++) cout<<"\n"<<lattice.fi[i]<<" "<<fi[i]<<" "<<lattice.old_fi[i]; //to sum &/or plot in Excel or Matlab
            cout  << "\nfi[1]-fi[2] + fi[7]-fi[8] + fi[9]-fi[10] + fi[11]-fi[12] + fi[13]-fi[14]="
            <<    fi[1]-fi[2] + fi[7]-fi[8] + fi[9]-fi[10] + fi[11]-fi[12] + fi[13]-fi[14] ;
            cout<<endl;
            Pause( "physical_wf: _slbmpi_h: 327: lattice.f for u.x==u.x?" ) ;
        }
#endif
        
        lattice.rho = std::accumulate(fi, fi + en , 0.0) ; // / Normalizer; // allow small turb. rh variations

        if (lattice.is_velinlet()) 
        {   
            u = u_i ;
            if ( IBC_D.PerturbVar == 1 ) u = u_i * Perturb( IBC_D.when , IBC_D.when2begin, IBC_D.when2end, IBC_D ) ; //cout<<"\nu_i="<<u;
        } else {
            if (lattice.is_on_face(Xn)) 
            {
                if ( ( IBC_D.periodicx != 1 ) || ( IBC_D.periodicy != 1 ) || ( IBC_D.periodicz != 1 ) )
                    u = xyz_dt(0.0, 0.0, 0.0);
            } else {
                if ( ( IBC_D.periodicx != 1 ) || ( IBC_D.periodicy != 1 ) || ( IBC_D.periodicz != 1 ) )
#ifdef D3Q19
                    u = xyz_dt( inner_product(fi, fi + en, cix, 0.0 ) , inner_product(fi, fi + en, ciy, 0.0 ) , inner_product(fi, fi + en, ciz, 0.0 ) ) ;
#endif
#ifndef D3Q19
                u = xyz_dt( inner_product(fi, fi + en, cix.begin(), 0.0 ) , inner_product(fi, fi + en, ciy.begin(), 0.0 ) , inner_product(fi, fi + en, ciz.begin(), 0.0 ) ) ;
#endif
            }
        }
        if ( IBC_D.IC_BC == 10 ) for ( int a_face = Yn ; a_face < END ; a_face++ ) if ( lattice.is_on_face( a_face ) ) u = xyz_dt(0.0, 0.0, 0.0);
        u = u / lattice.rho; // u = (u + 0.5 * f) / lattice.rho; //JCR: 0.5f not EEF but for BFF from Breyannis & Valougeorgis 
        lattice.T   = 0.0 ; 
        for(int i=0; i<=en-1;i++)
        { 
            lattice.T += ( cix[i] * cix[i] + ciy[i] * ciy[i] + ciz[i] * ciz[i] ) * fi[i] / lattice.rho ;
        } 
        lattice.p = Fluid_Species::p_EqState(lattice.rho,lattice.T) ;//consistent w/init_f;//lattice.p = lattice.rho * lattice.T * k_B / mass ;
        
        if(lattice.is_on_centerline()) 
        {
            ux_max = std::max(ux_max, u.x);
        }
        
#ifdef DEBUG
        if((mag(u)>1.0e10)||isnan(mag(u))||isinf(mag(u))||isnan(lattice.rho)||(lattice.rho>10.0)||(lattice.rho<Near0))
        {
            cout<<"\nLEAVING physical_wf: u= "<<u<<" & rh="<<lattice.rho<<" & T="<<lattice.T<<" & p="<<lattice.p<<" @ "<<lattice.idx<<" fi=";
            for (int i=0; i<en; i++) cout<<" fi[i="<<i<<"]= "<<lattice.fi[i];
            cout<<"\nLEAVING physical_wf: To sum &/or plot in Excel or Matlab:\n lattice.fi         fi   lattice.old_fi";
            for (int i=0; i<en; i++) cout<<"\n"<<lattice.fi[i]<<" "<<fi[i]<<" "<<lattice.old_fi[i]; //to sum &/or plot in Excel or Matlab
            cout  << "\nfi[1]-fi[2] + fi[7]-fi[8] + fi[9]-fi[10] + fi[11]-fi[12] + fi[13]-fi[14]="
            <<    fi[1]-fi[2] + fi[7]-fi[8] + fi[9]-fi[10] + fi[11]-fi[12] + fi[13]-fi[14] 
            << " == u.x? = " << u.x ;
            Pause( "LEAVING physical_wf: _slbmpi_h: 376: lattice.f for u.x==u.x" ) ;
            cout<<endl;
        }
#endif
        
    }
};

double physical_wf::ux_max = 0.0;

struct collision_mrt_ben_wf
{
    const IBC_data &IBc ;
    const Collision_Integrator &CollIntegrate ;
    v_Space &ff ;                       //whole collision integrated f product: \int{[f(v1)f(v2)-[f(v1')f(v2')]..} in v_Space freq. domain
    const v_Space &ff_S2L ; //mapped SB 2 LBM 2 b integrated over t as RHS out of freq. domain; both r 0 if !global
    const r_Space &nu_t ;   // if LES
        
        collision_mrt_ben_wf(const IBC_data &ibC,const Collision_Integrator &CollIntegrateur,v_Space &FF,const v_Space &FF_S2L,const r_Space &Nu_t ) : 
        IBc(ibC), CollIntegrate(CollIntegrateur), ff(FF), ff_S2L(FF_S2L), nu_t( Nu_t ) { }
    
    void operator()(lattice_type& lattice) const
    { 
        (void)(lattice) ;
    }
}  ;

#endif