/*
  slbmpi.cpp
  MHDA

  Created by Dr. Jacques C. Richard on 1/31/12.
*/

#include <math.h>
#include <float.h>
#include <algorithm>
#include <numeric>
#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <ostream>
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
#include <boost/lambda/algorithm.hpp>
#include <boost/lambda/if.hpp>
#include <boost/math/tools/roots.hpp>
#include <string>

/*  * /
//#undef _OPENMP
/ * * /
#ifdef _OPENMP
#include <omp.h> //http://openmp.org/ g++: -fopenmp; IBM: -qsmp=omp (conflict w/Array.h); icc: -openmp (conflict w/multi_array.h); flags #define _OPENMP
#endif
/ * */
#include "xyz_type.h"

/* to comment out boost::mpi */
#include <mpi.h>                                                                //http://www.open-mpi.org/ or http://www.mpi-forum.org/ include 1st for system MPI?
#include <boost/mpi.hpp>                                                        //http://www.boost.org/ #defines BOOST_MPI_HPP!=BOOST_MPI_ENVIRONMENT_HPP!=..

typedef struct {int xs; int xe; int ys; int ye; int zs; int ze;} Internal_Coords;

void mpe_decomp1d(int global_num_y, int num_procs, int myid, int *start_y, int *end_y) 
{
    int local_num_y, deficit;
    local_num_y  = global_num_y / num_procs;
    *start_y     = myid * local_num_y + 1;
    deficit      = global_num_y % num_procs;
    if(myid < deficit) 
    {
        *start_y += myid;
        local_num_y++;
    }
    else 
    {
        *start_y += deficit;
    }
    *end_y = *start_y + local_num_y - 1;
    if (*end_y > global_num_y || myid == num_procs - 1) { *end_y = global_num_y; }
    *start_y = *start_y - 1 ; //JCR: to ensure start @ 0 instead of 1!
} // mpe_decomp1d  from http://www.cct.lsu.edu/~scheinin/MPI_Tutorial/Exercises/Jacobi/jacobi.c 
 /* no boost::mpi if commented to here */

using namespace boost;
using namespace std;
using namespace boost::lambda;
/*
#include "fourier.h" //http://vision.eng.shu.ac.uk/jan/mimas/docs/group__fourierTransforms.html but needs work; #includes "boost/multi_array.hpp"
#include "multi_array_op.h"
typedef boost::array<double , 1 > array1d;
typedef boost::multi_array<double , 2 > array2d;	//http://www.cs.brown.edu/~jwicks/boost/libs/multi_array/doc/user.html
typedef boost::multi_array<double , 3 > array3d; 
typedef boost::multi_array< std::complex< double > , 3 > array3dc;

using namespace mimas;
*/
#ifndef LINE_MAX
#define LINE_MAX 256
#endif

const std::string CodeBG = "\n*** \n\
...\n" ; //https://parasol.tamu.edu/stapl/

const std::string help = "               HELP = Sample code input:\n\
1000000             = BIG # as marker of last line of input\n" ;

static inline int IDX(const int i, const int j, const int k, const int J, const int K){return((i*J+j)*K+k);}

#ifdef DEBUG
#include <stdexcept>
#include <exception>
#define ErP(a,b) continue; // static inline void ErP( const std::string a, const std::string b ) {std::cout << a << b << "\n";}
static inline void Pause( const std::string at )
{
    std::cout<<"\nPause @ \""<<at<<"\" if 1 process: <Enter> or <Return> continues; ^C aborts:";
    //cin.get( ); // comment this line to actually not pause but show where would have paused & continue
}
#endif

#ifndef DEBUG
#define ErP(a,b) continue;
static inline void Pause( const std::string at ){ return (void)( at ) ; }
#endif

#define MAX max //#define MAX(x,y) (x>y?x:y)  //quickly defines MAX and MIN functions; JCR: switch to C++ versions/templates
#define MAX3(x,y,z) MAX(MAX(x,y),z)
#define MIN min //#define MIN(x,y) (x>y?y:x)
#define MIN3(x,y,z) MIN(MIN(x,y),z)

//KGK Global constants-Need in math_functions.h // JCR: further globalized for spectral Boltzmann solver
const double Near0 = 1.0e-17 ;                                                  // JCR: small enough!?
const double pi = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651 ; //http://boost-sandbox.sourceforge.net/libs/math_constants/doc/html/index.html
//#include <boost/math/constants/constants.hpp>
//const double pi = boost::math::constants::pi<double>(); //http://stackoverflow.com/questions/1727881/how-to-use-the-pi-constant-in-c
const double Lambda = 2.0/(3.0+sqrt(double(2.0))), C0=1.0/(4.0*pi), Cnst = 1.0/(C0*16.0*pi*pi*pow(2.0*pi*Lambda,3.0));//in front of collision operator/kernel;smooth

//Global sizes for arrays & structs for chosen algorithm, i.e., less directly connected to physics 
xyz_type<int> Mode( 2 , 2 , 2 ) ;             //JCR: For Spectral Boltzmann solver; 2 1 1 closest 8MxMyMz to en=19; half-modes (-Mode..Mode)
int Mode3 = 8 * Mode.x * Mode.y * Mode.z + 2 ;// velocity space volume size = # of modes (product of)

xyz_type<int> problem_size( 4, 4, 4) ;        //32^3,10^3; 64,32,32; 160,121,81//KGK Define min(4,4,4)Refer:init_boundary_neighbors_wf boundary.h

enum face_type {Xn, Xp, Yn, Yp, Zn, Zp, END};//makes a list of constants; each item is assigned an ID #:0..6; JCR: lattice/cube faces for BC

// For a struct that has choices of solver, IC/BC, # of velocities, lattice, algorithm, mesh, etc.: 
#define D3Q19
#ifdef D3Q19
const int en = 19 ;							 //size of velocity set as per dflt LBM -- BUT THIS SHOULD NOT BE HARD-CODED!!!!
const int cix[en]={0, 1, -1,  0,  0,  0,  0,  1, -1,  1, -1,  1, -1,  1, -1,  0,  0,  0,  0};
const int ciy[en]={0, 0,  0,  1, -1,  0,  0,  1,  1, -1, -1,  0,  0,  0,  0,  1, -1,  1, -1};
const int ciz[en]={0, 0,  0,  0,  0,  1, -1,  0,  0,  0,  0,  1,  1, -1, -1,  1,  1, -1, -1};
#endif 

#ifndef D3Q19
const int en = 513 ; // (8*Mode.x*Mode.y*Mode.z+1)	//size of velocity set >= # of spectral modes
std::vector<int> cix ( en , 0 ) ; 
std::vector<int> ciy ( en , 0 ) ; 
std::vector<int> ciz ( en , 0 ) ; 
//typedef std::vector<double> PDF_type ;
#endif 
typedef std::vector<double> v_Space ;       //Global containers (SLOW if local) vs malloc 2 avoid malloc errors/conflict; fixed w/fftw++1.06!?
typedef std::vector<double> r_Space ;       //http://www.cplusplus.com/reference/stl/vector/ Compared to arrays, ..almost the same performance
typedef std::vector<double> tVec ;
typedef std::vector<double> WaveVec ;

const int gn = 7 ;                          //size of velocity set as per dflt LBM analogous vector PDF g solver for macro B

typedef xyz_type<double> xyz_dt;            //xyz_dt is an alias of type xyz_type using double numbers (for x, y, z)
typedef double           PDF_type[ en ] ;	//declares PDF_type as an alias for an array w/en indices; f!=force, f:=PDF
typedef xyz_dt       vec_PDF_type[ gn ] ;	//declares vec_PDF_type as an array of vectors in physical space (xyz) w/gn indices.

static inline double mag(const xyz_dt v){return(sqrt(v.x*v.x+v.y*v.y+v.z*v.z));}//Magnitude of a 3D vector ;
static inline std::vector<double> Mag( const std::vector<xyz_dt>& a )
{
    std::vector<double> b( a.size() , 0.0 ) ;
    for( std::vector<double>::iterator i = b.begin() ; i < b.end() ; i++ ) b[ (int)(*i) ] = mag( a[ (int)(*i) ] ) ;
    return b ;
    
}  //Magnitudes of 3D vectors in r_Space ;

const std::complex<double> ImUnit( 0.0 , 1.0 );	// the imaginary unit complex ImUnit = sqrt(-1.0);

//"Universal", hence, "global" constants: starting w/physical constants
const double unit_charge = 1.602176487e-19 ;//1.602176487e-19 Coulomb or Amp*seconds
const double k_B         = 1.380650424e-23 ;//Boltzmann's constant in units: 1.380650424e-23 J/K
const double epsilon0    = 8.854187817e-12 ;//Permitivity of free space = 8.854187817e-12 A*A s*s*s*s /kg m/^3 = C*C/N/m^2 or F/m
const double mu0         = 4.0e-7 * pi ;    //magnetic permeability of free space = 4.0e-7 * pi N/A^2 = H/m ≈ 1.2566×10−6 H/m (or T·m/A)
const double c_L         = 1.0 / sqrt( mu0 * epsilon0 ) ; // m/s, speed of light
const double Avogadro	 = 6.0221417930e+23 ;// Avogadro constant is the # of "elementary entities" (atoms or molecules) in 1 mole
const double mass_e      = 9.1094e-28 ;      // kg, mass of electron
const double mass_p      = 1.6726e-24 ;      // kg, mass of proton
/*
 //Dflt case non-dimensionalization/normalization parameters: VASIMR data
 double n0_V				 = 1.0e+20 ;                // /m^3; 
 double b0_V				 = 1.0 ;                    // T or Tesla
 double en0V				 = 100.0 ;                  // eV ion energy where 1 eV = 1.6021765314e-19 J = 11604.50520 K = e / k_B
 double j0_V				 = 10.0 ;                   // A/m^2
 double P0_V				 = 200.0e+3 ;               // W for VX-200 VASIMR
 double ndot0_V			 = 3.1e+20 ;              // ions/sec
 double Vdot0_V			 = 100.0e-6 / 60 ;        // VASIMR range: 30 - 150 sccm
 */
//  Fluid-specific kinematic, thermodynamic & electromagnetic properties; some could b global but some should really b part of class/struct
double mass				 = 2.163517e-25 ;     //2.163517e-25 kg = FluidCpt mass; 
double e_m				 = 1.0;//unit_charge / mass ; // shows up often in plasmas; 1 if unitless for unit_charge/mass;
double MolarMass		 = 131.30 ;           // Xenon Molar Mass = Moleculat Wt = 131.30 g/mol = 131.30 kg/kmol 
double AtomicNum		 = 54.0 ;             // Xenon Atomic Number = Z
double rhoXe			 = 5.894e+3 ;         // Density: 5.9g/L @ 273K & 1atm = 5.894e+3 kg/m^3
double GasConst = 1.0;//e3 * Avogadro*k_B/MolarMass;// Specific gas constant (mult by 1e3 since N_A is per mol & need kmol)
double Normalizer		 = 1.0 ;              // to normalize 0th moment of f, (rho), to yield rho=1 incomp. cases, using 0.019 for a while
const int Max_species	 = 5; 
int num_species	 = 1;

class Fluid_Species 
{
private:
    PDF_type fi, old_fi ; 
public:
    static inline double r_EqState(const double P, const double t){ return ( ( P /(GasConst*t + Near0) )/Normalizer ) ; }//density normalized
    static inline double p_EqState(const double r, const double t){ return ( ( r * GasConst*t )			*Normalizer ) ; }//density normalized?
    static inline double f_Eq(double c_x,double c_y,double c_z,double cm_x,double cm_y,double cm_z,const double sigma,const double Konst)
    {	//3D Maxwell
        double Exp_Limit = 7.0 ;											// JCR: limit range of xp for exp( xp )
        double u2        = ( c_x - cm_x ) * ( c_x - cm_x ) + ( c_y - cm_y ) * ( c_y - cm_y ) + ( c_z - cm_z ) * ( c_z - cm_z ) ;
        double xp        = - 0.5 * ( u2 ) / ( sigma ) ;						//std::cout<<"xp="<<xp<<" u^2="<<u2<<" "; 
        xp        = ((xp > - Exp_Limit ) ? (xp) : ( - Exp_Limit )) ;		//std::cout<<"xp="<<xp<<" ";
        xp        = ((xp <   Exp_Limit ) ? (xp) : (   Exp_Limit )) ;		//std::cout<<"xp="<<xp<<" ";
        return ( Konst * exp( xp ) ) ; 
    } ;
    double rho , T , p ;
    double sigma( double Tmp = 1.0 ){return(Tmp);}//http://research.nianet.org/~luo/Reprints-luo/Notes/TUBraunschweig_03/TUBraunschweig_03.pdf
    double Konst( double Tmp = 1.0 , double Rho = 1.0 ){ return( Rho * pow( 2.0 * pi * GasConst * Tmp , -1.5 ) ) ;} //Dflt form as @ Luo
    //double Konst( double Temp = 1.0 , double Rho = 1.0 ){ return( Rho * pow( mass / ( 2.0 * pi * k_B * Temp ) , 1.5 ) ) ;} //Luo
    xyz_dt u, old_u, grad_rho, J, curlJ ;
    inline void Initialize_f( PDF_type f )
    {
        //Normalizer = 0.0 ;
        for (int i = 0; i <= en - 1 ; i++)
        {
            f[ i ] = f_Eq( double(cix[i]), double(ciy[i]), double(ciz[i]), 0.,0.,0., sigma(1.), Konst(1.,1.) ) ;
            //Normalizer += f[ i ] ; std::cout<<" Normalizer="<<Normalizer ;
        }
        Normalizer = *max_element( &f[ 0 ] , &f[ en - 1 ] ) ;
        //Normalizer = std::accumulate( f, f , 0.0 ) ;			//std::cout<<" Normalizer="<<Normalizer ;
    } ;
    Fluid_Species( const int Cpt = num_species )
    { 
        Initialize_f(     fi ) ; 
        Initialize_f( old_fi ) ; 
        rho = std::accumulate(fi, fi + en , 0.0); 
        //Normalizer = rho ; 	std::cout<<" Normalizer="<<Normalizer ;
        rho = rho / Normalizer ;
        T = 0.0 ;
        if ( Cpt == 1 ) for(int i=0; i<=en-1;i++) T += ( cix[i] * cix[i] + ciy[i] * ciy[i] + ciz[i] * ciz[i] ) * fi[i] / ( 3.*en*rho ) ;
        if ( ( Cpt == num_species ) && ( Cpt < Max_species ) ) p = p_EqState( rho , T ) ;
        u		 = ( 0.0 ) ;
        old_u	 = ( 0.0 ) ;
        J		 = rho * u ; //std::cout<<"J="<<J;
        grad_rho = ( 0.0 ) ;
        curlJ	 = ( 0.0 ) ;
    } ;// Constructor
    ~Fluid_Species() 
    {
        //delete [ ] fi; delete [ ] old_fi ; // 4-1-2010: IBM's xlC_r: delete operator not allowed for type "double [19]"
    } ; // Destructor
}; 

struct fluid // single fluid macroscopic equivalent (e.g., density-averaged)
{
public:
	static inline double gamma(const double P, const double t){ return ( ((P<100)&&(P>0.1)&&(t<1000)&&(t>0.1))?1.4:1.33 ) ; } //dflt Std "Atm"
	double rho ;// (double R[ Max_species ] ){ return( std::accumulate( R, R + Max_species, 0.0 ) / num_species ) ; } ; // average
	double GasCon ;
	double T ; 
	double p ;
    static inline double A_over_A_star( const double M, const double g ) { return( pow( 2.0*(1.0+0.5*(g-1.0)*M*M)/(g+1.0) , 0.5*(g+1.0)/(g-1.0) )/M ) ; }
    static inline double M_of_A( const double A_ratio = 1.0, double lo_guess_M = 1.001e-0, double hi_guess_M = 3.001e+1, const double g = 1.4 ) ;
	xyz_dt u, old_u, grad_rho, J, curlJ ;
};

template <class T> // area ratio eq. to get Mach # for any fluid (i.e., any gamma = g)
struct A_over_A_star_functor
{
    A_over_A_star_functor( const T& gam, const T& Aratio) : g(gam), A_ratio(Aratio) {}
    T operator()( const T& M)
    {
        return ( fluid::A_over_A_star( M , g ) - A_ratio ) ;
    }
private:
    T g, A_ratio, Result;
};

inline double fluid::M_of_A( const double A_ratio , double lo_guess_M , double hi_guess_M, const double g ) // find M(A/A*) by bi-section/interval halving 
{
    boost::math::tools::eps_tolerance<double> tol_preci( std::numeric_limits<double>::digits ) ;
    int iter  = 999 ;
    std::pair< double , double > M_root ( 0.0 , 0.0 ) ;
    M_root = boost::math::tools::bisect( A_over_A_star_functor<double>( g , A_ratio ) , lo_guess_M , hi_guess_M , tol_preci , (boost::uintmax_t&)(iter) ) ;
    return M_root.first ;                         //Subsonic or supersonic based on input, so generally M_root.first = M_root.second
} 

#include "lattice_type.h"

typedef std::vector<lattice_type> grid_type;	//grid_type is a vector of lattice_type;

//Case parameters
const int MaxLoops		 = 10 ;                 //Max Number of Loops
struct Interpolate{ int pre ; int post ; } ; //int Interpolate::pre = 2 ; int Interpolate::post = 2 ; //interpolate pre/post spectral
//struct Interpolate{ static int pre ; static int post ; } ; int Interpolate::pre = 2 ; int Interpolate::post = 2 ; //interpolate pre/post spectral

struct IBC_data 
{	
	int	   IC_BC ;
	xyz_dt u0,B0,E0;                            //initial vector fields for IC's
	double T_u, T_d, p_u , p_d  ;               //initial scalar variables for IC's
	int		when, when2begin, when2end ;        //when to apply, start & stop applying shock &/or perturbation in inlet velocity
    double  PerturbAmp ;                        // the perturbation strength
    double  PerturbFreq ;                       // the perturbation frequency
    double  PerturbPhase ;                      //
    int     PerturbType ;                       //Const=1,Lin=2,Sin=3; 
    int     PerturbVar ;                        //u=1,Dp=2,B=3,E=4
	int     periodicx, periodicy, periodicz ;   // 1 = periodic
	xyz_type<int> where ;                       //coordinates (in grid index) for where to apply shock or perturbation
	xyz_type<int> Physical_Grid ;               //physical space size
	double rinlet ;					
	double rexit ;
	double nlength ;
	double xshift ;
} ;	

static inline double Perturb( const int t , const int t_s, const int t_e, const IBC_data &u_p)
{return( ((t>=t_s) && (t<=t_e))? 1.0+u_p.PerturbAmp*sin((u_p.PerturbFreq*(double)(t)+u_p.PerturbPhase)*pi/180.0) : 0. );}

struct Collision_Integrator
{
	double uxx_LBM ;                            // c_s_LBM > u_LBE = v_rms
	double tau_f ;                              // LBM relaxation time; 1/tau_f = fluid frequency. 
	double tau_b ;                              // 1/tau_b = magnetic field frequency w/
	double dpdx ;                               // Pressure gradient in x direction
	unsigned int solver ;                       // choice of solver: MRT=1|SRT=2|L'B',SB=3|FD E&M,SB=4|Spectral=5|potentials,spectral=6;
	int Compres ;                               // Compressible = 1 for yes else 0=no
	double Dt ;                                 // Time step for non-LBM temporal integration or scaling MHDT energy/spectra evolution
	Interpolate PolyInterpolate ;               // whether or not to interpolate to/from LBM to SB & back
	unsigned int Mode3 ;                        // 8*Mode.x*y*z
	unsigned int Smooth ;                       // whether to use a smoother or not & which type
	xyz_type<int> HMode ;                       // the SB half-mode sizes in x, y, z
} ;

#include "no_boundary_iterator.h"
#include "slbmpi.h"
#include "boundary.h"


int main(int argc, char **argv)
{
#ifdef BOOST_MPI_HPP
    boost::mpi::environment  env( argc , argv ) ;
    const int root = 0 ;
    int IO_Rank = root ; 
    
    cout <<"\nboost::mpi environment'ed w/env.host_rank()="<<env.host_rank()
    << " & env.io_rank()="<<env.io_rank()
    << " & IO_Rank=" << IO_Rank ;//still = 1!? Maybe all processes can IO!?
    boost::mpi::communicator world ; 
    cout << "\nProcess=" << world.rank( ) 
    <<"'s boost::mpi::communicator 'world' is set up w/"<<world.size()<<" processes" ;//error here @ .size() @ Mac Pro but not mini
    if ( env.initialized() ) 
    {
        cout <<"\nProcess="<<world.rank()<<" reports env.initialized!\n";    
    }
    else 
    {//} but Mac Pro still needs:
        /* valgrind MPI_Init errors like w/env( argc, argv ) if "argc, argv" included but MPI_Init...(arg.. below needed to run */
#if defined(_OPENMP)
        cout<<"\nMPI_Init_thread'ed w/code "<<MPI_Init_thread(&argc,&argv,MPI_THREAD_SERIALIZED,&num_species )<<" (0=MPI_SUCCESS)!\n";
#endif
#if !defined(_OPENMP) && !defined(__PETSC_H)
        cout<<"\nMPI_Init'ed w/code " << MPI_Init(&argc, &argv ) <<" (0=MPI_SUCCESS)!\n" ; 
#endif
#if defined(__PETSC_H)
        cout<<"\nPetscInitialize'd w/code " << PetscInitialize(&argc,&argv,(char*)0,PETSC_NULL) <<" (0=MPI_SUCCESS)!\n" ;
#endif
    } //non-Mac Pro may simply report calling MPI_Init twice as erroneous & may continue!
#endif
    
	Fluid_Species FluidCpt( 1 ) ;
	
	Interpolate   PlInt ;
	
	int t4ReStart=0 ;//JCR: Control variables for re-start files saved in case of a crash: SLOW!!!!
	
	double  tau_f             = 1.0;      //0.74; 1/tau_f = fluid frequency. 1/tau_b = magnetic field frequency w/LBM relaxation time
	double  tau_b             = 1.0;      //1.14; .66; h=8 Rm=20    //3.7; h=8 Rm=1  //1.14; h=8 Rm=5  //.58; h=8 Rm=40
	int     Smooth            = 3 ;       //JCR: Whether or not to smooth f; 1=none, 2=Lanczos, 3=Raised cosine, 4=Fejer, 5=WENO
	int		solver            = 5 ;       //0=dflt (do nothing); MRT=1|SRT=2|L'B',SB=3|FD E&M,SB=4|Spectral=5|potentials,spectral=6;
	int		Compress          = 0 ;       //0=dflt compressible flow(yes=1,no=0)
	double	Dt                = 1.00e-0 ; //JCR: step size in SB integrator
	int		t                 = t4ReStart;// loop start/initial time
    double Smagorinsky        = 0.1 ;
    double tMaxEd             = 2.0 ;
    
	IBC_data IBC_dat ;                    //set initial data
	IBC_dat.IC_BC             = 1 ;       //1=jet @x=0, extrap @NX, extrap; 2=DIT, periodic; 3=TG vortex, periodic; 4:uniform field(s);5:shk
	IBC_dat.u0                = 0.0;      //dflt initial velocity (inlet or shock)
	IBC_dat.B0                = 0.0;      //dflt initial magnetic field
	IBC_dat.E0                = 0.0;      //dflt initial electric field
	IBC_dat.where             = xyz_type<int>(0,0,0);                           //dflt location of moving shock front or perturbation
	IBC_dat.when              = t ;	
    IBC_dat.when2begin        = 20 ;
    IBC_dat.when2end          = 30 ;
    IBC_dat.PerturbAmp       = 0.01 ;
	IBC_dat.periodicx         = 1 ;       //1=periodic
	IBC_dat.periodicy         = 1 ;
	IBC_dat.periodicz         = 1 ;
	IBC_dat.T_u               = 16.0 ;
	IBC_dat.T_d               = 1.0 ;
	IBC_dat.p_u               = 8.0 ;
	IBC_dat.p_d               = 1.0 ;
	IBC_dat.rinlet            = 10;
	IBC_dat.rexit             = 6;
	IBC_dat.nlength           = 20;
	IBC_dat.xshift            = 0.0;
    
	int			 numLoops     = 0 ;                                             //FHE Number of Loops
	double loopX[MaxLoops]    = { 10.0 } ;                                      //KSB&FHE X-Location of the loops
	double loopR[MaxLoops]    = { 5.0 } ;                                       //KSB&FHE Radius of the loops
	double loopI[MaxLoops]    = { 00.0 };                                       //KSB&FHE Current of the loops
	double dpdx               = 1.0e-4;                                         //dpdx for pressure-driven flows
    
	lattice_type::width		  = 24;                                             //24;  // choose an even number, width of the inlet
	lattice_type::height	  = 16;                                             //16;  // choose an even number, lattice_type::height of the inlet
	lattice_type::inletradius = 5 ;                                             //end of dflts (e.g., RJF3_10.in).
    
    std::cout<< "\nChecking for input!\n" ;
#ifdef BOOST_MPI_HPP
    if ( world.rank() == IO_Rank ) 
    { 
        std::cout<< "\nProcess "<<world.rank()<<" is checking/reading input file!" ;
#endif
        
        /*   */ // comment i/o starting here if all input data w/in source code as dflt data
        std::cout << CodeBG.c_str() << endl;
         /*   */ // comment i/o ending here if not using!
        
#ifdef BOOST_MPI_HPP
        // then share inputed data w/the rest of the 'world'
        boost::mpi::broadcast( world , problem_size.x       , IO_Rank ) ;
        boost::mpi::broadcast( world , problem_size.y       , IO_Rank ) ;
        boost::mpi::broadcast( world , problem_size.z       , IO_Rank ) ;
        boost::mpi::broadcast( world , Mode.x               , IO_Rank ) ;
        boost::mpi::broadcast( world , Mode.y               , IO_Rank ) ;
        boost::mpi::broadcast( world , Mode.z               , IO_Rank ) ;
        boost::mpi::broadcast( world , Dt                   , IO_Rank ) ;
        boost::mpi::broadcast( world , solver               , IO_Rank ) ;
        boost::mpi::broadcast( world , Compress             , IO_Rank ) ;
        boost::mpi::broadcast( world , tau_f                , IO_Rank ) ;
        boost::mpi::broadcast( world , tau_b                , IO_Rank ) ;
        boost::mpi::broadcast( world , PlInt.pre            , IO_Rank ) ;
        boost::mpi::broadcast( world , PlInt.post           , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.IC_BC        , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.periodicx    , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.periodicy    , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.periodicz    , IO_Rank ) ;
        boost::mpi::broadcast( world , Smagorinsky          , IO_Rank ) ;
        boost::mpi::broadcast( world , tMaxEd               , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.when2begin   , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.when2end     , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.PerturbType  , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.PerturbVar   , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.PerturbAmp  , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.PerturbFreq  , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.PerturbPhase , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.u0.x         , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.u0.y         , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.u0.z         , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.B0.x         , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.B0.y         , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.B0.z         , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.E0.x         , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.E0.y         , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.E0.z         , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.T_u          , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.T_d          , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.p_u          , IO_Rank ) ;
        boost::mpi::broadcast( world , IBC_dat.p_d          , IO_Rank ) ;
        boost::mpi::broadcast( world , numLoops             , IO_Rank ) ;
        boost::mpi::broadcast( world , loopX , numLoops     , IO_Rank ) ;
        boost::mpi::broadcast( world , loopR , numLoops     , IO_Rank ) ;
        boost::mpi::broadcast( world , loopI , numLoops     , IO_Rank ) ;
        boost::mpi::broadcast( world , lattice_type::inlettype   , IO_Rank ) ;
        boost::mpi::broadcast( world , lattice_type::width       , IO_Rank ) ;
        boost::mpi::broadcast( world , lattice_type::height      , IO_Rank ) ;
        boost::mpi::broadcast( world , lattice_type::inletradius , IO_Rank ) ;
    }
#endif
	
    xyz_dt u = IBC_dat.u0 ; //jet velocity of incoming plasma; JCR: to set IC v field instead of u(0.01, 0.0, 0.0)
	
	Mode3 = 8 * Mode.x * Mode.y * Mode.z + 2 ; // ( 2 * Mode.x + 2 ) * ( 2 * Mode.y + 2 ) * ( 2 * Mode.z + 2 ) + 2 ; // + 2 as globally set above
    int ProblemSize = (problem_size.x+2)*(problem_size.y+2)*(problem_size.z+2)+2;// physical space volume size
	
	//KGK Grid is of type-grid_type (which is a Vector of type-lattice_type and of size 34*34*34 elements)
	grid_type grid((problem_size.x+2)*(problem_size.y+2)*(problem_size.z+2));   //JCR: phase space = vector=physical r_space of lattice=v_space
    /*
     #ifdef BOOST_MPI_HPP //decompose grid http://boost.2283326.n4.nabble.com/mpi-why-this-code-hangs-on-gather-td2572658.html
     boost::mpi::broadcast( world , grid , 0 ) ;
     unsigned long chunkSize = grid.size() / world.size() ; cout << "\nchunkSize=" << chunkSize ;
     unsigned long extraBits = grid.size() % world.size() ; cout << ", extraBits=" << extraBits ;
     for( int iProc = 1 ; iProc < world.size() ; ++iProc ){ cout << ", iProc=" << iProc ;
     unsigned long Pos = extraBits + iProc * chunkSize ; cout << ", Pos=" << Pos ;
     for(unsigned int jProc = 0 ; jProc < chunkSize; ++jProc ){ cout << ", jProc=" << jProc ;
     grid.push_back( grid[ Pos ]); 
     Pos++ ; 
     } 
     } 
     #endif
     */
	//KGK Set the "static" fields (width, height, diameter ..); unique - only one set for all lattices
	//KGK Should be done through an initializer/constructor function call to honor information hiding
    if(lattice_type::inlettype == 0)
    {
		lattice_type::diameter = 2.0*sqrt(lattice_type::width*lattice_type::height/pi); // Equivalent diameter of inlet
    }
    if(lattice_type::inlettype ==1)
    {
        lattice_type::diameter = 2.0*(lattice_type::inletradius);
    }
    
    if ( max( lattice_type::width , lattice_type::height ) / min( lattice_type::width , lattice_type::height ) > 10 )
    {
        lattice_type::diameter = max( lattice_type::width , lattice_type::height ) ;
        cout << "\nNOTE: The hydraulic [equivalent area] diameter = max( width , height ) since " 
        << "max( width , height ) / min( width , height ) > 10, i.e., assuming 2D, so effectively, \n" ;
    }
    
	lattice_type::center_x        = (problem_size.x-1) / 2 + 1;
	lattice_type::center_y        = (problem_size.y-1) / 2 + 1;
	lattice_type::center_z        = (problem_size.z-1) / 2 + 1;
	lattice_type::dims            =  problem_size + 2;
    /*
     #ifdef BOOST_MPI_HPP
     lattice_type::dims            =  problem_size * world.size( ) + 2 ;         //So init_indices generates for whole region?!
     #endif
     */
	lattice_type::begin_iterator  = grid.begin();                               //RAD: object is grid, function is begin()
	lattice_type::hh.x            = 1./problem_size.x;                          //KGK Reference grid (w/o boundary cells)size=1; 
	lattice_type::hh.y            = 1./problem_size.y;                          // and
	lattice_type::hh.z            = 1./problem_size.z;                          // hh will make solutions consistent w/grid refinement
	
	double L_ph                   = ( IBC_dat.IC_BC == 1 ? lattice_type::diameter : 2.0 * pi ) ;
	
	IBC_dat.when                  = t4ReStart ;
	IBC_dat.Physical_Grid         = problem_size ;
	
	Collision_Integrator IntegrateVars ;
	IntegrateVars.uxx_LBM         = L_ph / (double)( problem_size.x - 1 ) ;     //JCR: /1 not / time_ph L/NX for CFL? 0.061=2pi/(nnf+3)? Ben scaled by this for turb. physical non-dim. quantities!?
	IntegrateVars.tau_f           = tau_f ;
	IntegrateVars.tau_b           = tau_b;
	IntegrateVars.dpdx            = dpdx ;
	IntegrateVars.solver          = solver ;
	IntegrateVars.Compres         = Compress ;
	IntegrateVars.Dt              = Dt ;
	IntegrateVars.Mode3           = Mode3 ;
	IntegrateVars.HMode           = Mode ;
	IntegrateVars.PolyInterpolate = PlInt ; 
	IntegrateVars.Smooth          = Smooth ; 
	
	//6 vars below currently not used in computation, kept for output consistency; JCR: for checking run case parameters
	double viscosity              = (2.0 * tau_f-1.)/6.;
	double mag_resist             = (2.0 * tau_b-1.)/8.;
	double Re                     = IBC_dat.u0.x*lattice_type::diameter/viscosity;
	double Rem                    = IBC_dat.u0.x*lattice_type::diameter/mag_resist;
    cout     << "the Reynolds Number is " << Re << "\n";
    cout     << "the Magnetic Reynolds Number is " << Rem << "\n";
	
	std::vector<xyz_dt> vv( ProblemSize , 0.0 ) ;//3D vector fields; ~1%SLOWER local but for parallel;IBM xlC error:xyz_type<double> v[ProblemSize]
	std::vector<xyz_dt> Bb( ProblemSize , 0.0 ) ;
	std::vector<xyz_dt> ee( ProblemSize , 0.0 ) ;
	
	r_Space t_0( ProblemSize , Lambda * pi / 6.0 ) , r_0( ProblemSize , 1.0 ) ; // 3D scalar fields & Pareschi et al. "sigma"
    
    r_Space nu_t( ProblemSize , 0.0 ) ;                                         // Turbulent viscosity for LES
	
	v_Space b1( 64* Mode.x * Mode.y * Mode.z * Mode.x * Mode.y * Mode.z + 1 , 0.0 );
	v_Space b2( 8 * Mode.x * Mode.y * Mode.z + 1 , 0.0 );
	
	v_Space vx( Mode3 , 0.0 ); 
	v_Space vy( Mode3 , 0.0 ); 
	v_Space vz( Mode3 , 0.0 ); 
	v_Space cx( MAX( Mode3 , en ) , 0.0 ) ;
	v_Space cy( MAX( Mode3 , en ) , 0.0 ) ;
	v_Space cz( MAX( Mode3 , en ) , 0.0 ) ;
    
	std::vector< std::complex<double> > S1( 8 * Mode.x * Mode.y * Mode.z + 1 , ImUnit );//v/freq. space=transformation/translation index shift @ 
	std::vector< std::complex<double> > S2( 8 * Mode.x * Mode.y * Mode.z + 1 , ImUnit );//MatLab:sum -HalfMode to HalfMode but implementation-dep.?
    
	v_Space ff(		Mode3 , 0.0 ) ;				//whole collision integrated f product: \int{[f(v1)f(v2)-[f(v1')f(v2')]..} in v_Space freq. domain
	v_Space ff_S2L( Mode3 , 0.0 ) ;				//mapped SB 2 LBM 2 b integrated over t as RHS out of freq. domain; both r 0 if !global
    
	
#if defined(DEBUG)
	if( ( (solver>2) && (solver<7) ) && ( (ff.size()==0) || (ff_S2L.size()==0) || ((int)ff_S2L.size()!=Mode3) || ((int)ff.size()!=Mode3) ) )
		cout<<"\n WARNING: main: SIZES after resize(solver="<<solver<<"): ff:"<<ff.size()<<" ff_S2L:"<<ff_S2L.size()<<" & Mode3="<<Mode3<<"\n"; // checks OK!
#endif
    
#ifndef BOOST_MPI_HPP
    set_lattice_ID( problem_size , grid ) ; //KGK set lattice ID for each lattice, Required for SPECTRAL and FD methods; + 2 inside it
    
    std::for_each(grid.begin(),grid.end(),init_indices(problem_size+2));			//RAD: applies init_indices to each element from begin to end
    std::for_each(grid.begin(),grid.end(),init_internal_neighbors_wf<lattice_type>());//KGK RJ //RAD: offsets from previously dealt neighbor
#endif
#ifdef BOOST_MPI_HPP 
    xyz_type<int> MiniSize = lattice_type::dims / world.size() ; 
    cout << "\nProcess=" << world.rank( ) << "'s MiniGridSize=" << MiniSize ;
    Internal_Coords I_C ;
    /* a la http://boost.2283326.n4.nabble.com/mpi-why-this-code-hangs-on-gather-td2572658.html * /
     if ( MiniSize.x < 4 ) MiniSize.x = lattice_type::dims.x ;
     if ( MiniSize.y < 4 ) MiniSize.y = lattice_type::dims.y ;
     if ( MiniSize.z < 4 ) MiniSize.z = lattice_type::dims.z ;
     if ( MiniSize != lattice_type::dims ){
     xyz_type<int> xtraBits = xyz_type<int>( problem_size.x % world.size() , problem_size.y % world.size() , problem_size.z % world.size() ) ; cout << ", xtraBits=" << xtraBits ;
     for( int iProc = 1 ; iProc < world.size() ; ++iProc ){ cout << ", iProc=" << iProc ;
     xyz_type<int> Loc = xtraBits + iProc * MiniSize ; cout << ", Loc=" << Loc ;
     for(xyz_type<int> jProc = 0 ; jProc < MiniSize; jProc = jProc + 1 ){ cout << ", jProc=" << jProc ;
     grid_type MiniGrid( grid.size() / world.size() ) ;
     set_lattice_ID( MiniSize , MiniGrid ) ; 
     Loc = Loc + 1 ; 
     } 
     }
     Pause( "main: 868: Domain decomposition" ) ;
     } else
     set_lattice_ID( problem_size , grid ) ; //KGK set lattice ID for each lattice, Required for SPECTRAL and FD methods; + 2 inside it
     / * */
    mpe_decomp1d( lattice_type::dims.x , world.size( ), world.rank( ), &I_C.xs , &I_C.xe ) ;
    mpe_decomp1d( lattice_type::dims.y , world.size( ), world.rank( ), &I_C.ys , &I_C.ye ) ;
    mpe_decomp1d( lattice_type::dims.z , world.size( ), world.rank( ), &I_C.zs , &I_C.ze ) ;
    MiniSize = xyz_type<int>( I_C.xe - I_C.xs , I_C.ye - I_C.ys , I_C.ze - I_C.zs ) ; 
    cout << "\nProcess=" << world.rank( ) 
    << "'s mpe_decomp1d's Internal_Coords I_C="<<I_C.xs<<" "<<I_C.xe<<" "<<I_C.ys<<" "<<I_C.ye<<" "<<I_C.zs<<" "<<I_C.ze ;
    cout << "\nProcess=" << world.rank( ) << "'s mpe_decomp1d MiniGridSize=" << MiniSize << std::endl ;
    if ( MiniSize != lattice_type::dims )
    {
        grid_type MiniGrid( grid.size() / world.size() ) ;
        set_lattice_ID( MiniSize , MiniGrid ) ; 
        Pause( "main: 637: Domain decomposition" ) ;
    } else
        set_lattice_ID( lattice_type::dims , grid ) ; //KGK set lattice ID for each lattice, Required for SPECTRAL and FD methods; + 2 inside it
    
    Pause( "main: 641: pre-init_indices" ) ;
    std::for_each(grid.begin(),grid.end(),init_indices(lattice_type::dims));		//JCR: applies init_indices to each element of all region?!
    Pause( "main: 643: post-init_indices & pre-init_internal_neighbors_wf" ) ;
    std::for_each(grid.begin(),grid.end(),init_internal_neighbors_wf<lattice_type>( world ));//KGK RJ //RAD: offsets from previously dealt neighbor
    Pause( "main: 645: post-init_indices & -init_internal_neighbors_wf" ) ;
#endif
    
    //#ifdef D3Q19
    // if ( ( IBC_dat.IC_BC != 2 ) || ( IBC_dat.IC_BC != 3 ) || ( IBC_dat.IC_BC != 9 ) )   // these if's don't have much fx on RJ
    std::for_each(grid.begin(),grid.end(),init_boundary_neighbors_wf<lattice_type>());//JCR:Original RJ BC;see only far x BC w/o!Conflict periodic?NaN w/o
    //#endif
#ifndef D3Q19
    std::for_each(grid.begin(),grid.end(),init_internal_neighbors_SB<lattice_type>(Mode,vx,vy,vz));//JCR: offsets, in v_SB directions, from neighbor
    std::for_each(grid.begin(),grid.end(),Set_SB_stream_directions(Mode,vx,vy,vz));//set initial streaming directions
#endif
    //std::for_each(grid.begin(),grid.end(),init_boundary_neighbors_wall4BC<lattice_type>(IBC_dat));//JCR: Set wall 4 BC;breaks rho...!
    
    //JCR: dflt initializations or by IC_BC case @ switch below
    // if ( ( IBC_dat.IC_BC == 2 ) || ( IBC_dat.IC_BC == 3 ) || ( IBC_dat.IC_BC == 9 ) || ( IBC_dat.periodicx == 1 ) || ( IBC_dat.periodicy == 1 ) || ( IBC_dat.periodicz == 1 ) ) //these if's don't have much fx on RJ & need i_b_n_p_BC to set onBoundaryTag
    if ( ( solver != 5 ) || ( solver != 6 ) )
        std::for_each(grid.begin(),grid.end(),init_boundary_neighbors_periodic_BC<lattice_type>(IBC_dat));//KGK Set neighbors for Periodic BC;JCR:fixes rho & sets lattice.onBoundaryTag; needed by E&M spectral solvers; may NaN w/o; still reflect inlet @ other faces w/o it
	
    std::cout<<"Output files may exist based on case!\n";
    
#ifdef BOOST_MPI_HPP
    if ( env.finalized() ) //boost usually invokes the related destructor @ program end so may not be finalized here
    {
        cout     <<"\nProcess="<<world.rank()<<" reports env.finalized!\n"; 
    }
    else 
    {//} but Mac Pro still needs:
#if !defined(__PETSC_H)
        cout<<"\nMPI_Finalize'd w/code " << MPI_Finalize( ) <<" (0=MPI_SUCCESS, presumably threads too)!\n" ; 
#endif
#if defined(__PETSC_H)
        cout<<"\nPetscFinalize'd w/code "<< PetscFinalize() <<" (0=MPI_SUCCESS)!\n";//has MPI_Finalize(); 
#endif
	}
#endif
    
    return EXIT_SUCCESS ;
}