#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>

#include "analysis/pmmc.h"
#include "common/ScaLBL.h"
#include "common/MPI_Helpers.h"
#include "common/Communication.h"
#include "IO/Mesh.h"
#include "IO/Writer.h"
#include "ProfilerApp.h"

using namespace std;


std::shared_ptr<Database> readInput( const std::string& file )
{
    if ( exists( file )
        return std::make_shared<Database>( file );
    auto db = std::make_shared<Database>();
    db.putData( "Color", std::make_shared<Database>() );
    db.putData( "Domain", std::make_shared<Database>() );
    return db;
}


int main(int argc, char **argv)
{
  // Initialize MPI
  int provided_thread_support = -1;
  MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
  MPI_Comm comm;
  MPI_Comm_dup(MPI_COMM_WORLD,&comm);
  int rank = comm_rank(comm);
  int nprocs = comm_size(comm);
  if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE )
    std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
  { // Limit scope so variables that contain communicators will free before MPI_Finialize

    // parallel domain size (# of sub-domains)
    int nprocx,nprocy,nprocz;

    MPI_Request req1[18],req2[18];
    MPI_Status stat1[18],stat2[18];

    if (rank == 0){
        printf("********************************************************\n");
        printf("Running Hybrid Implementation of Color LBM    \n");
        printf("********************************************************\n");
    }

    PROFILE_ENABLE(1);
    //PROFILE_ENABLE_TRACE();
    //PROFILE_ENABLE_MEMORY();
    PROFILE_SYNCHRONIZE();
    PROFILE_START("Main");
    Utilities::setErrorHandlers();
    // Variables that specify the computational domain  
    string FILENAME;
    unsigned int nBlocks, nthreads;
    int Nx,Ny,Nz;
    int nspheres;
    double Lx,Ly,Lz;
    // Color Model parameters
    int i,j,k,n;

    // pmmc threshold values
    double fluid_isovalue,solid_isovalue;
    fluid_isovalue = 0.0;
    solid_isovalue = 0.0;
    nBlocks = 32;
    nthreads = 128;
    
    int RESTART_INTERVAL=1000;

    auto db = readInput( "input.in" );

    // Read variables from Color
    auto color_db = db->getDatabase( "Color" );
    auto tau   = color_db->getScalarWithDefault<double>( "tau", 1.0 );  
    auto alpha = color_db->getScalarWithDefault<double>( "alpha", 1e-3 );
    auto beta  = color_db->getScalarWithDefault<double>( "beta", 0.95 );
    auto phi_s = color_db->getScalarWithDefault<double>( "phi_s", 0.0 );
    auto Fx    = color_db->getVectorWithDefault<double>( "F", { 0, 0, 0 } )[0];
    auto Fy    = color_db->getScalarWithDefault<double>( "F", { 0, 0, 0 } )[1];
    auto Fz    = color_db->getScalarWithDefault<double>( "F", { 0, 0, 0 } )[2];
    auto Restart = color_db->getScalarWithDefault<bool>( "Restart", 0 );
    auto pBC   = color_db->getScalarWithDefault<double>( "pBC", 0 );
    auto din   = color_db->getScalarWithDefault<double>( "din", 1.0 );
    auto dout  = color_db->getScalarWithDefault<double>( "dout", 1.0 );
    auto timestepMax = color_db->getScalarWithDefault<int>( "timestepMax", 3 );
    auto interval = color_db->getScalarWithDefault<int>( "interval", 100 );
    auto tol   = color_db->getScalarWithDefault<double>( "tol", 1e-6 );
    auto das   = color_db->getScalarWithDefault<double>( "das", 0.1 );
    auto dbs   = color_db->getScalarWithDefault<double>( "dab", 0.9 );

    // Read variables from Domain
    auto domain_db = db->getDatabase( "Domain" );
    auto n = domain_db->getVectorWithDefault<int>( "n", { 50, 50, 50 } );
    auto L = domain_db->getVectorWithDefault<double>( "L", { 1.0, 1.0, 1.0 } );
    auto nproc = domain_db->getVectorWithDefault<int>( "nproc", { 1, 1, 1 } );
    auto nspheres = domain_db->getScalarWithDefault<int>( "nspheres", 1 );
    int Nx = n[0];
    int Ny = n[1];
    int Nz = n[2];
    int nprocx = nproc[0];
    int nprocy = nproc[1];
    int nprocz = nproc[2];
    double Lx = L[0];
    double Ly = L[1];
    double Lz = L[2];

    // Get the rank info
    const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
    int iproc = rank_info.ix;
    int jproc = rank_info.jy;
    int kproc = rank_info.kz;

    MPI_Barrier(comm);
    // **************************************************************
    // **************************************************************
    double Ps = -(das-dbs)/(das+dbs);
    double rlxA = 1.f/tau;
    double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
    double xIntPos;
    xIntPos = log((1.0+phi_s)/(1.0-phi_s))/(2.0*beta);     
    
    if (nprocs != nprocx*nprocy*nprocz){
        printf("Fatal error in processor number! \n");
        printf("   nprocx =  %i \n",nprocx);
        printf("   nprocy =  %i \n",nprocy);
        printf("   nprocz =  %i \n",nprocz);
        return 1;
    }

    if (rank==0){
        printf("********************************************************\n");
        printf("tau = %f \n", tau);
        printf("alpha = %f \n", alpha);        
        printf("beta = %f \n", beta);
        printf("das = %f \n", das);
        printf("dbs = %f \n", dbs);
        printf("Value of phi at solid surface = %f \n", phi_s);
        printf("Distance to phi = 0.0: %f \n", xIntPos);
        printf("gamma_{wn} = %f \n", 5.796*alpha);
        // printf("cos theta_c = %f \n", 1.05332*Ps);
        printf("Force(x) = %f \n", Fx);
        printf("Force(y) = %f \n", Fy);
        printf("Force(z) = %f \n", Fz);
        printf("Sub-domain size = %i x %i x %i\n",Nz,Nz,Nz);
        printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
        printf("********************************************************\n");
    }

    // Full domain used for averaging (do not use mask for analysis)
    Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,pBC);
    Dm.CommInit(comm);

    // Mask that excludes the solid phase
    Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,pBC);
     
     MPI_Barrier(comm);

    Nx+=2; Ny+=2; Nz += 2;

    int N = Nx*Ny*Nz;
    int dist_mem_size = N*sizeof(double);

    int S = N/nthreads/nBlocks+1;

    // unsigned int nBlocks = N/nthreads + (N%nthreads == 0?0:1);
    // dim3 grid(nBlocks,1,1);

    if (rank==0) printf("Number of blocks = %i \n", nBlocks);
    if (rank==0) printf("Threads per block = %i \n", nthreads);
    if (rank==0) printf("Sweeps per thread = %i \n", S);
    if (rank==0) printf("Number of nodes per side = %i \n", Nx);
    if (rank==0) printf("Total Number of nodes = %i \n", N);
    if (rank==0) printf("********************************************************\n");

    //.......................................................................
    if (rank == 0)    printf("Read input media... \n");
    //.......................................................................
    
    //.......................................................................
    // Filenames used
    char LocalRankString[8];
    char LocalRankFilename[40];
    char LocalRestartFile[40];
    sprintf(LocalRankString,"%05d",rank);
    sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
    sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
    
//    printf("Local File Name =  %s \n",LocalRankFilename);
    // .......... READ THE INPUT FILE .......................................
//    char value;
    char *id;
    id = new char[N];
    int sum = 0;
    double iVol_global = 1.0/(1.0*Nx*Ny*Nz*nprocs);
    double porosity = 0;

    DoubleArray SDs(Nx,Ny,Nz);
    DoubleArray SDn(Nx,Ny,Nz);
    //.......................................................................
    
    double BubbleRadius = 15.5; // Radius of the capillary tube
    sum=0;
    for (k=0;k<Nz;k++){
        for (j=0;j<Ny;j++){
            for (i=0;i<Nx;i++){
                n = k*Nx*Ny + j*Nz + i;
                // Cylindrical capillary tube aligned with the z direction
                SDs(i,j,k) = 100;
                // Initialize phase positions field
                if (SDs(i,j,k) < 0.0){
                    id[n] = 0;
                }
                else {
                    id[n] = 1;
                    sum++;
                }
            }
        }
    }
    
    // Set up MPI communication structurese
    if (rank==0)    printf ("Setting up communication control structures \n");

	// Initialize communication structures in averaging domain
	for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
	Mask.CommInit(comm);

    //......................................................................................
    // Create a communicator for the device
    ScaLBL_Communicator ScaLBL_Comm(Mask);

    //...........device phase ID.................................................
    if (rank==0)    printf ("Copying phase ID to device \n");
    char *ID;
    ScaLBL_AllocateDeviceMemory((void **) &ID, N);                        // Allocate device memory
    // Copy to the device
    ScaLBL_CopyToDevice(ID, id, N);
    ScaLBL_DeviceBarrier();
    //...........................................................................

    //...........................................................................
    //                MAIN  VARIABLES ALLOCATED HERE
    //...........................................................................
    // LBM variables
    if (rank==0)    printf ("Allocating distributions \n");
    //......................device distributions.................................
    double *f_even,*f_odd;
    double *A_even,*A_odd,*B_even,*B_odd;
    //...........................................................................
    ScaLBL_AllocateDeviceMemory((void **) &f_even, 10*dist_mem_size);    // Allocate device memory
    ScaLBL_AllocateDeviceMemory((void **) &f_odd, 9*dist_mem_size);    // Allocate device memory
    ScaLBL_AllocateDeviceMemory((void **) &A_even, 4*dist_mem_size);    // Allocate device memory
    ScaLBL_AllocateDeviceMemory((void **) &A_odd, 3*dist_mem_size);    // Allocate device memory
    ScaLBL_AllocateDeviceMemory((void **) &B_even, 4*dist_mem_size);    // Allocate device memory
    ScaLBL_AllocateDeviceMemory((void **) &B_odd, 3*dist_mem_size);    // Allocate device memory
    //...........................................................................
    double *Phi,*Den;
    double *ColorGrad, *Velocity, *Pressure;
    //...........................................................................
    ScaLBL_AllocateDeviceMemory((void **) &Phi, dist_mem_size);
    ScaLBL_AllocateDeviceMemory((void **) &Pressure, dist_mem_size);
    ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
    ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
    ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
    //...........................................................................
    //...........................................................................
    // Phase indicator (in array form as needed by PMMC algorithm)
    DoubleArray Phase(Nx,Ny,Nz);

    // Extra copies of phi needed to compute time derivatives on CPU
    DoubleArray Phase_tminus(Nx,Ny,Nz);
    DoubleArray Phase_tplus(Nx,Ny,Nz);
    DoubleArray dPdt(Nx,Ny,Nz);

    //copies of data needed to perform checkpointing from cpu
    double *cDen, *cDistEven, *cDistOdd;
    cDen = new double[2*N];
    cDistEven = new double[10*N];
    cDistOdd = new double[9*N];

    // data needed to perform CPU-based averaging
    //double *Vel = new double[3*N];        // fluid velocity

    DoubleArray MeanCurvature(Nx,Ny,Nz);
    DoubleArray GaussCurvature(Nx,Ny,Nz);
    DoubleArray SDs_x(Nx,Ny,Nz);        // Gradient of the signed distance
    DoubleArray SDs_y(Nx,Ny,Nz);
    DoubleArray SDs_z(Nx,Ny,Nz);
    DoubleArray Phase_x(Nx,Ny,Nz);            // Gradient of the phase indicator field
    DoubleArray Phase_y(Nx,Ny,Nz);
    DoubleArray Phase_z(Nx,Ny,Nz);
    DoubleArray Press(Nx,Ny,Nz);
    MeanCurvature.fill(0);
    GaussCurvature.fill(0);
    SDs_x.fill(0);
    SDs_y.fill(0);
    SDs_z.fill(0);
    Phase_x.fill(0);
    Phase_y.fill(0);
    Phase_z.fill(0);
    Press.fill(0);
    
    /*****************************************************************
     VARIABLES FOR THE PMMC ALGORITHM
     ****************************************************************** */
    //...........................................................................
    // Averaging variables
    //...........................................................................
    // local averages (to each MPI process)
    double awn,ans,aws,lwns,nwp_volume;
    double As;
    double vol_w, vol_n;                        // volumes the exclude the interfacial region
    double sat_w;
    double pan,paw;                                // local phase averaged pressure
    // double vx_w,vy_w,vz_w,vx_n,vy_n,vz_n;          // local phase averaged velocity
    // Global averages (all processes)
    double vol_w_global, vol_n_global;            // volumes the exclude the interfacial region
    double awn_global,ans_global,aws_global;
    double lwns_global;
    double efawns,efawns_global;                // averaged contact angle
    double Jwn,Jwn_global;                        // average mean curavture - wn interface
    DoubleArray van(3);
    DoubleArray vaw(3);
    DoubleArray vawn(3);
    DoubleArray Gwn(6);
    DoubleArray Gns(6);
    DoubleArray Gws(6);
    
    double nwp_volume_global;                    // volume for the wetting phase (for saturation)
    // double p_n_global,p_w_global;                // global phase averaged pressure
    // double vx_w_global,vy_w_global,vz_w_global;    // global phase averaged velocity
    // double vx_n_global,vy_n_global,vz_n_global;    // global phase averaged velocity
    double As_global;
    double dEs,dAwn,dAns;                        // Global surface energy (calculated by rank=0)
    double pan_global,paw_global;                                // local phase averaged pressure
    DoubleArray van_global(3);
    DoubleArray vaw_global(3);
    DoubleArray vawn_global(3);
    DoubleArray Gwn_global(6);
    DoubleArray Gns_global(6);
    DoubleArray Gws_global(6);
    //...........................................................................

    // bool add=1;            // Set to false if any corners contain nw-phase ( F > fluid_isovalue)
    int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};  // cube corners
    DoubleArray CubeValues(2,2,2);
    // int count_in=0,count_out=0;
    // int nodx,nody,nodz;
    // initialize lists for vertices for surfaces, common line
    DTMutableList<Point> nw_pts(20);
    DTMutableList<Point> ns_pts(20);
    DTMutableList<Point> ws_pts(20);
    DTMutableList<Point> nws_pts(20);
    // initialize triangle lists for surfaces
    IntArray nw_tris(3,20);
    IntArray ns_tris(3,20);
    IntArray ws_tris(3,20);
    // initialize list for line segments
    IntArray nws_seg(2,20);
    DTMutableList<Point> tmp(20);
    DoubleArray Values(20);
    DoubleArray ContactAngle(20);
    DoubleArray Curvature(20);
    DoubleArray InterfaceSpeed(20);
    DoubleArray NormalVector(60);
    
    // IntArray store;
    
    int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0;
    int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
    
    // double s,s1,s2,s3;        // Triangle sides (lengths)
    Point A,B,C,P;
    // double area;
    
    // Initialize arrays for local solid surface
    DTMutableList<Point> local_sol_pts(20);
    int n_local_sol_pts = 0;
    IntArray local_sol_tris(3,18);
    int n_local_sol_tris;
    DoubleArray values(20);
    DTMutableList<Point> local_nws_pts(20);
    int n_local_nws_pts;
    
    //int n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg;
    int c;
    //int newton_steps = 0;
    //...........................................................................
    int ncubes = (Nx-2)*(Ny-2)*(Nz-2);    // Exclude the "upper" halo
    IntArray cubeList(3,ncubes);
    int nc=0;
    //...........................................................................
    // Set up the cube list (very regular in this case due to lack of blob-ID)
    for (k=0; k<Nz-2; k++){
        for (j=0; j<Ny-2; j++){
            for (i=0; i<Nx-2; i++){
                cubeList(0,nc) = i;
                cubeList(1,nc) = j;
                cubeList(2,nc) = k;
                nc++;
            }
        }
    }
    //...........................................................................
    // Grids used to pack faces on the GPU for MPI
    int faceGrid,edgeGrid,packThreads;
    packThreads=512;
    edgeGrid=1;
    faceGrid=Nx*Ny/packThreads;
    NULL_USE(faceGrid);
    NULL_USE(edgeGrid);
    //...........................................................................    

    
    //...........................................................................
    int timestep = 0;
    if (rank==0) printf("********************************************************\n");
    if (rank==0)    printf("No. of timesteps: %i \n", timestepMax);

    //.......create a stream for the LB calculation.......
    //    cudaStream_t stream;
    //    cudaStreamCreate(&stream);

    //.......create and start timer............
    double starttime,stoptime,cputime;
    MPI_Barrier(comm);
    starttime = MPI_Wtime();
    //.........................................
    //...........................................................................
    //                MAIN  VARIABLES INITIALIZED HERE
    //...........................................................................
    double BubRad[5];
    BubRad[0] = 8.0;
    BubRad[1] = 10.0;
    BubRad[2] = 12.0;
    BubRad[3] = 15.0;
    BubRad[4] = 20.0;
    //...........................................................................
    if (rank==0){
        printf("--------------------------------------------------------------------------------------\n");
        printf("radius ");                                // Timestep, Change in Surface Energy
        printf("sw pw pn awn Jwn ");                    // Scalar averages
        printf("Gwn [xx, yy, zz, xy, xz, yz] ");        // Orientation tensors
        printf("--------------------------------------------------------------------------------------\n");
    }

    for (int bubbleCount=0; bubbleCount<4; bubbleCount++){
        char bubbleCountName[20];
        sprintf(bubbleCountName,"bubbleCount-%i",bubbleCount);
        PROFILE_START(bubbleCountName);
        BubbleRadius = BubRad[bubbleCount]; // Radius of the current bubble

        // Initialize the bubble
        for (k=0;k<Nz;k++){
            for (j=0;j<Ny;j++){
                for (i=0;i<Nx;i++){
                    n = k*Nx*Ny + j*Nz + i;
                    int iglobal= i+(Nx-2)*iproc;
                    int jglobal= j+(Ny-2)*jproc;
                    int kglobal= k+(Nz-2)*kproc;
                    // Cylindrical capillary tube aligned with the z direction
                    SDs(i,j,k) = 100;
                    // Initialize phase positions field
                    // if ((i-0.5*Nx)*(i-0.5*Nx)+(j-0.5*Ny)*(j-0.5*Ny)+(k-0.5*Nz)*(k-0.5*Nz) < BubbleRadius*BubbleRadius){
                    //     id[n] = 2;
                    // }
                    // Initialize phase position field for parallel bubble test
                    if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
                            +(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)
                            +(kglobal-0.5*(Nz-2)*nprocz)*(kglobal-0.5*(Nz-2)*nprocz) < BubbleRadius*BubbleRadius){
                        id[n] = 2;
                    }
                    else{
                        id[n]=1;
                    }
                }
            }
        }
        // Copy the bubble to the device and initialize
        ScaLBL_CopyToDevice(ID, id, N);

        //...........................................................................
        //...........................................................................
        ScaLBL_D3Q19_Init(ID, f_even, f_odd, Nx, Ny, Nz);
        //......................................................................
        //ScaLBL_Color_InitDistance(ID, Den, Phi, SDs.data(), das, dbs, beta, xIntPos, Nx, Ny, Nz);
        ScaLBL_Color_Init(ID, Den, Phi, das, dbs, Nx, Ny, Nz);
        ScaLBL_D3Q7_Init(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
        ScaLBL_D3Q7_Init(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
        //......................................................................
        // Once phase has been initialized, map solid to account for 'smeared' interface
        //......................................................................
        for (i=0; i<N; i++)    SDs(i) -= (1.0); //
        //......................................................................
        //.......................................................................
        sprintf(LocalRankString,"%05d",rank);
        sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
        WriteLocalSolidID(LocalRankFilename, id, N);
        sprintf(LocalRankFilename,"%s%s","SDs.",LocalRankString);
        WriteLocalSolidDistance(LocalRankFilename, SDs.data(), N);
        //.......................................................................
        if (Restart == true){
            if (rank==0) printf("Reading restart file! \n");
            // Read in the restart file to CPU buffers
            ReadCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
            // Copy the restart data to the GPU
            ScaLBL_CopyToDevice(f_even,cDistEven,10*N*sizeof(double));
            ScaLBL_CopyToDevice(f_odd,cDistOdd,9*N*sizeof(double));
            ScaLBL_CopyToDevice(Den,cDen,2*N*sizeof(double));
            ScaLBL_DeviceBarrier();
            MPI_Barrier(comm);
        }

        //*************************************************************************
        //         Compute the phase indicator field and reset Copy, Den
        //*************************************************************************
        ScaLBL_ComputePhaseField(ID, Phi, Den, N);
        //*************************************************************************    
        ScaLBL_DeviceBarrier();
        ScaLBL_Comm.SendHalo(Phi);
        ScaLBL_Comm.RecvHalo(Phi);
        ScaLBL_DeviceBarrier();
        MPI_Barrier(comm);
        //*************************************************************************

        if (rank==0 && pBC){
            printf("Setting inlet pressure = %f \n", din);
            printf("Setting outlet pressure = %f \n", dout);
        }
        if (pBC && kproc == 0)    {
            ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);            
            ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
        }
            
        if (pBC && kproc == nprocz-1){
            ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
            ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
        }

        //...........................................................................
        // Copy the phase indicator field for the earlier timestep
        ScaLBL_DeviceBarrier();
        ScaLBL_CopyToHost(Phase_tplus.data(),Phi,N*sizeof(double));
        //...........................................................................
        //...........................................................................
        // Copy the data for for the analysis timestep
        //...........................................................................
        // Copy the phase from the GPU -> CPU
        //...........................................................................
        ScaLBL_DeviceBarrier();
        ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
        ScaLBL_CopyToHost(Phase.data(),Phi,N*sizeof(double));
        ScaLBL_CopyToHost(Press.data(),Pressure,N*sizeof(double));
        MPI_Barrier(comm);
        //...........................................................................
        
        int timestep = 0;

        //************ MAIN ITERATION LOOP ***************************************/
        PROFILE_START("Time-loop");
        while (timestep < timestepMax){


            //*************************************************************************
            // Fused Color Gradient and Collision 
            //*************************************************************************
            ScaLBL_D3Q19_ColorCollide( ID,f_even,f_odd,Phi,ColorGrad,
                                 Velocity,Nx,Ny,Nz,rlxA,rlxB,alpha,beta,Fx,Fy,Fz);
            //*************************************************************************

            ScaLBL_DeviceBarrier();
            //*************************************************************************
            // Pack and send the D3Q19 distributions
            ScaLBL_Comm.SendD3Q19(f_even, f_odd);
            //*************************************************************************
            
            //*************************************************************************
            //         Carry out the density streaming step for mass transport
            //*************************************************************************
            ScaLBL_D3Q7_ColorCollideMass(ID, A_even, A_odd, B_even, B_odd, Den, Phi, 
                                    ColorGrad, Velocity, beta, N, pBC);
            //*************************************************************************

            ScaLBL_DeviceBarrier();
            MPI_Barrier(comm);
            //*************************************************************************
            //         Swap the distributions for momentum transport
            //*************************************************************************
            ScaLBL_D3Q19_Swap(ID, f_even, f_odd, Nx, Ny, Nz);
            //*************************************************************************

            ScaLBL_DeviceBarrier();
            MPI_Barrier(comm);
            //*************************************************************************
            // Wait for communications to complete and unpack the distributions
            ScaLBL_Comm.RecvD3Q19(f_even, f_odd);
            //*************************************************************************

            ScaLBL_DeviceBarrier();
            //*************************************************************************
            // Pack and send the D3Q7 distributions
            ScaLBL_Comm.BiSendD3Q7(A_even, A_odd, B_even, B_odd);
            //*************************************************************************

            ScaLBL_DeviceBarrier();
            ScaLBL_D3Q7_Swap(ID, A_even, A_odd, Nx, Ny, Nz);
            ScaLBL_D3Q7_Swap(ID, B_even, B_odd, Nx, Ny, Nz);

            ScaLBL_DeviceBarrier();
            MPI_Barrier(comm);

            //*************************************************************************
            // Wait for communication and unpack the D3Q7 distributions
            ScaLBL_Comm.BiRecvD3Q7(A_even, A_odd, B_even, B_odd);
            //*************************************************************************

            ScaLBL_DeviceBarrier();
            //..................................................................................
            ScaLBL_D3Q7_Density(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
            ScaLBL_D3Q7_Density(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
            
            //*************************************************************************
            //         Compute the phase indicator field 
            //*************************************************************************
            // ScaLBL_ComputePhaseField(ID, Phi, Copy, Den, N);
            ScaLBL_DeviceBarrier();
            MPI_Barrier(comm);
        
            ScaLBL_ComputePhaseField(ID, Phi, Den, N);
            //*************************************************************************
            ScaLBL_Comm.SendHalo(Phi);
            ScaLBL_DeviceBarrier();
            ScaLBL_Comm.RecvHalo(Phi);
            //*************************************************************************

            ScaLBL_DeviceBarrier();

            // Pressure boundary conditions
            if (pBC && kproc == 0)    {
                ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);            
	            ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
                ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
            }
                
            if (pBC && kproc == nprocz-1){
                ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
	            ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
                ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
            }
            
            //...................................................................................

            MPI_Barrier(comm);

            // Timestep completed!
            timestep++;

            if (timestep%RESTART_INTERVAL == 0){
                // Copy the data to the CPU
                ScaLBL_CopyToHost(cDistEven,f_even,10*N*sizeof(double));
                ScaLBL_CopyToHost(cDistOdd,f_odd,9*N*sizeof(double));
                ScaLBL_CopyToHost(cDen,Den,2*N*sizeof(double));
                // Read in the restart file to CPU buffers
                WriteCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
            }
        }
        PROFILE_STOP("Time-loop");
        // End the bubble loop
        //...........................................................................
        // Copy the phase indicator field for the later timestep
        ScaLBL_DeviceBarrier();
        ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
        ScaLBL_CopyToHost(Phase_tminus.data(),Phi,N*sizeof(double));
        ScaLBL_CopyToHost(Phase_tplus.data(),Phi,N*sizeof(double));
        ScaLBL_CopyToHost(Phase.data(),Phi,N*sizeof(double));
        ScaLBL_CopyToHost(Press.data(),Pressure,N*sizeof(double));

        double temp=0.5/beta;
        for (n=0; n<N; n++){
          double value = Phase.data()[n];
            SDn(n) = temp*log((1.0+value)/(1.0-value));
        }

        //...........................................................................
        // Calculate the time derivative of the phase indicator field
        for (n=0; n<N; n++)    dPdt(n) = 0.1*(Phase_tplus(n) - Phase_tminus(n));
        //...........................................................................

        // Compute the gradients of the phase indicator and signed distance fields
        pmmc_MeshGradient(Phase,Phase_x,Phase_y,Phase_z,Nx,Ny,Nz);
        pmmc_MeshGradient(SDs,SDs_x,SDs_y,SDs_z,Nx,Ny,Nz);
        //...........................................................................
        // Compute the mesh curvature of the phase indicator field
        pmmc_MeshCurvature(SDn, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
        //...........................................................................
        // Fill in the halo region for the mesh gradients and curvature
        //...........................................................................
        // Pressure
        Dm.CommunicateMeshHalo(Press);
        // Mean Curvature
        Dm.CommunicateMeshHalo(MeanCurvature);
        // Gaussian Curvature
        Dm.CommunicateMeshHalo(GaussCurvature);
        // Gradient of the phase indicator field
        Dm.CommunicateMeshHalo(Phase_x);
        Dm.CommunicateMeshHalo(Phase_y);
        Dm.CommunicateMeshHalo(Phase_z);

        //...........................................................................
        // Compute areas using porous medium marching cubes algorithm
        // McClure, Adalsteinsson, et al. (2007)
        //...........................................................................
        awn = aws = ans = lwns = 0.0;
        nwp_volume = 0.0;
        As = 0.0;

        // Compute phase averages
        pan = paw = 0.0;
        vaw(0) = vaw(1) = vaw(2) = 0.0;
        van(0) = van(1) = van(2) = 0.0;
        vawn(0) = vawn(1) = vawn(2) = 0.0;
        Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
        Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
        Gws(0) = Gws(1) = Gws(2) = 0.0;
        Gws(3) = Gws(4) = Gws(5) = 0.0;
        Gns(0) = Gns(1) = Gns(2) = 0.0;
        Gns(3) = Gns(4) = Gns(5) = 0.0;
        vol_w = vol_n =0.0;
        Jwn = efawns = 0.0;
        
        for (c=0;c<ncubes;c++){
            // Get cube from the list
            i = cubeList(0,c);
            j = cubeList(1,c);
            k = cubeList(2,c);

            //...........................................................................
            // Compute volume averages
            for (int p=0;p<8;p++){

                double delphi;
                if ( SDs(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
                    // 1-D index for this cube corner
                    n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
                    // compute the norm of the gradient of the phase indicator field
                    delphi = sqrt(Phase_x(n)*Phase_x(n)+Phase_y(n)*Phase_y(n)+Phase_z(n)*Phase_z(n));
                    // Compute the non-wetting phase volume contribution
                    if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
                        nwp_volume += 0.125;
                        // volume the excludes the interfacial region
                        if (delphi < 1e-4){
                            vol_n += 0.125;
                            // pressure
                            pan += 0.125*Press(n);
                        }
                    }
                    else if (delphi < 1e-4){
                        // volume the excludes the interfacial region
                        vol_w += 0.125;
                        // pressure
                        paw += 0.125*Press(n);
                    }
                }
            }

            //...........................................................................
            // Construct the interfaces and common curve
            pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue,
                    nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
                    local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
                    n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
                    n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
                    i, j, k, Nx, Ny, Nz);

            // Integrate the contact angle
            efawns += pmmc_CubeContactAngle(CubeValues,Values,Phase_x,Phase_y,Phase_z,SDs_x,SDs_y,SDs_z,
                    local_nws_pts,i,j,k,n_local_nws_pts);

            // Integrate the mean curvature
            Jwn    += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);

            pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
                    NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);

            //...........................................................................
            // Compute the Interfacial Areas, Common Line length
            /*        awn += pmmc_CubeSurfaceArea(nw_pts,nw_tris,n_nw_tris);
            ans += pmmc_CubeSurfaceArea(ns_pts,ns_tris,n_ns_tris);
            aws += pmmc_CubeSurfaceArea(ws_pts,ws_tris,n_ws_tris);
             */
            As  += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);

            // Compute the surface orientation and the interfacial area
            awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
            ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris);
            aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris);
            lwns +=  pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
            //...........................................................................
        }
        //...........................................................................
        MPI_Barrier(comm);
        MPI_Allreduce(&nwp_volume,&nwp_volume_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&awn,&awn_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&ans,&ans_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&aws,&aws_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&lwns,&lwns_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&As,&As_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&Jwn,&Jwn_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,comm);
        // Phase averages
        MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&paw,&paw_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&pan,&pan_global,1,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&vaw(0),&vaw_global(0),3,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&van(0),&van_global(0),3,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&vawn(0),&vawn_global(0),3,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&Gwn(0),&Gwn_global(0),6,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&Gns(0),&Gns_global(0),6,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Allreduce(&Gws(0),&Gws_global(0),6,MPI_DOUBLE,MPI_SUM,comm);
        MPI_Barrier(comm);
        //.........................................................................
        // Compute the change in the total surface energy based on the defined interval
        // See McClure, Prins and Miller (2013) 
        //.........................................................................
        dAwn += awn_global;
        dAns += ans_global;
        dEs = 6.01603*alpha*(dAwn + 1.05332*Ps*dAns);
        dAwn = -awn_global;        // Get ready for the next analysis interval
        dAns = -ans_global;

        // Normalize the phase averages 
        // (density of both components = 1.0)
        paw_global = paw_global / vol_w_global;
        vaw_global(0) = vaw_global(0) / vol_w_global;
        vaw_global(1) = vaw_global(1) / vol_w_global;
        vaw_global(2) = vaw_global(2) / vol_w_global;
        pan_global = pan_global / vol_n_global;
        van_global(0) = van_global(0) / vol_n_global;
        van_global(1) = van_global(1) / vol_n_global;
        van_global(2) = van_global(2) / vol_n_global;

        // Normalize surface averages by the interfacial area
        Jwn_global /= awn_global;
        efawns_global /= lwns_global;

        if (awn_global > 0.0)    for (i=0; i<3; i++)        vawn_global(i) /= awn_global;
        if (awn_global > 0.0)    for (i=0; i<6; i++)        Gwn_global(i) /= awn_global;
        if (ans_global > 0.0)    for (i=0; i<6; i++)        Gns_global(i) /= ans_global;
        if (aws_global > 0.0)    for (i=0; i<6; i++)        Gws_global(i) /= aws_global;

        sat_w = 1.0 - nwp_volume_global*iVol_global/porosity;
        // Compute the specific interfacial areas and common line length (per unit volume)
        awn_global = awn_global*iVol_global;
        ans_global = ans_global*iVol_global;
        aws_global = aws_global*iVol_global;
        lwns_global = lwns_global*iVol_global;
        dEs = dEs*iVol_global;

        //.........................................................................
        if (rank==0){
            printf("%.5g ",BubbleRadius);                                    // bubble radius
            printf("%.5g %.5g ",paw_global,pan_global);            // saturation and pressure
            printf("%.5g ",awn_global);                                        // interfacial area
            printf("%.5g ",Jwn_global);                                        // curvature of wn interface
            printf("%.5g %.5g %.5g %.5g %.5g %.5g \n",
                    Gwn_global(0),Gwn_global(1),Gwn_global(2),Gwn_global(3),Gwn_global(4),Gwn_global(5));        // orientation of wn interface
        }


        shared_ptr<IO::TriList> mesh( new IO::TriList() );
        mesh->A.reserve(8*ncubes);
        mesh->B.reserve(8*ncubes);
        mesh->C.reserve(8*ncubes);
        for (c=0;c<ncubes;c++){
            // Get cube from the list
            i = cubeList(0,c);
            j = cubeList(1,c);
            k = cubeList(2,c);
            //...........................................................................
            // Construct the interfaces and common curve
            pmmc_ConstructLocalCube(SDs, Phase, solid_isovalue, fluid_isovalue,
                    nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
                    local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
                    n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
                    n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
                    i, j, k, Nx, Ny, Nz);

            //.......................................................................................
            // Write the triangle lists to text file
            for (int r=0;r<n_nw_tris;r++){
                A = nw_pts(nw_tris(0,r));
                B = nw_pts(nw_tris(1,r));
                C = nw_pts(nw_tris(2,r));
                // compare the trianlge orientation against the color gradient
                // Orientation of the triangle
                double tri_normal_x = (A.y-B.y)*(B.z-C.z) - (A.z-B.z)*(B.y-C.y);
                double tri_normal_y = (A.z-B.z)*(B.x-C.x) - (A.x-B.x)*(B.z-C.z);
                double tri_normal_z = (A.x-B.x)*(B.y-C.y) - (A.y-B.y)*(B.x-C.x);

                double normal_x = Phase_x(i,j,k);
                double normal_y = Phase_y(i,j,k);
                double normal_z = Phase_z(i,j,k);

                // If the normals don't point in the same direction, flip the orientation of the triangle
                // Right hand rule for triangle orientation is used to determine rendering for most software
                if (normal_x*tri_normal_x + normal_y*tri_normal_y + normal_z*tri_normal_z < 0.0){
                    P = A;
                    A = C;
                    C = P;
                }
                mesh->A.push_back(A);
                mesh->B.push_back(B);
                mesh->C.push_back(C);
            }        
        }
        std::vector<IO::MeshDataStruct> meshData(1);
        meshData[0].meshName = "wn-tris";
        meshData[0].mesh = mesh;
        for (size_t k=0; k<meshData.size(); k++) {
            shared_ptr<IO::Variable> dist( new IO::Variable() );
            dist->name = "distance";
            dist->dim = 1;
            dist->type = IO::VariableType::NodeVariable;
            dist->data.resize(3*mesh->A.size());
            for (size_t i=0; i<mesh->A.size(); i++) {
                const Point& a = mesh->A[i];
                const Point& b = mesh->B[i];
                const Point& c = mesh->C[i];
                dist->data(3*i+0) = sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
                dist->data(3*i+1) = sqrt(b.x*b.x+b.y*b.y+b.z*b.z);
                dist->data(3*i+2) = sqrt(c.x*c.x+c.y*c.y+c.z*c.z);
            }
            meshData[k].vars.push_back(dist);
        }
        IO::writeData( bubbleCount, meshData, comm );

        PROFILE_STOP(bubbleCountName);
    }
    NULL_USE(sat_w);

    //************************************************************************/
    ScaLBL_DeviceBarrier();
    MPI_Barrier(comm);
    stoptime = MPI_Wtime();
    if (rank==0) printf("-------------------------------------------------------------------\n");
    // Compute the walltime per timestep
    cputime = (stoptime - starttime)/timestep;
    // Performance obtained from each node
    double MLUPS = double(Nx*Ny*Nz)/cputime/1000000;
    
    if (rank==0) printf("********************************************************\n");
    if (rank==0) printf("CPU time = %f \n", cputime);
    if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
    MLUPS *= nprocs;
    if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
    if (rank==0) printf("********************************************************\n");
    
    //************************************************************************/
    // Write out the phase indicator field 
    //************************************************************************/
    sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
    //    printf("Local File Name =  %s \n",LocalRankFilename);
    // ScaLBL_CopyToHost(Phase.data(),Phi,N*sizeof(double));

    FILE *PHASE;
    PHASE = fopen(LocalRankFilename,"wb");
    fwrite(Press.data(),8,N,PHASE);
    // fwrite(MeanCurvature.data(),8,N,PHASE);
    fclose(PHASE);
    
/*  double *DensityValues;
    DensityValues = new double [2*N];
    ScaLBL_CopyToHost(DensityValues,Copy,2*N*sizeof(double));
    FILE *PHASE;
    PHASE = fopen(LocalRankFilename,"wb");
    fwrite(DensityValues,8,2*N,PHASE);
    fclose(PHASE);
*/    //************************************************************************/
    PROFILE_SAVE("TestBubble");
    
    // ****************************************************
    MPI_Barrier(comm);
  } // Limit scope so variables that contain communicators will free before MPI_Finialize
  MPI_Comm_free(&comm);
  MPI_Finalize();
  return 0;
}