diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 372cb73a..36b49063 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -481,8 +481,100 @@ ScaLBL_Communicator::ScaLBL_Communicator(Domain &Dm){ } ScaLBL_Communicator::~ScaLBL_Communicator(){ - // destrutor does nothing (bad idea) - // -- note that there needs to be a way to free memory allocated on the device!!! + // Free communicator + MPI_Comm_free(&MPI_COMM_SCALBL); + MPI_Group_free(&Group); + // Free internal memory + ScaLBL_FreeDeviceMemory( sendbuf_x ); + ScaLBL_FreeDeviceMemory( sendbuf_X ); + ScaLBL_FreeDeviceMemory( sendbuf_y ); + ScaLBL_FreeDeviceMemory( sendbuf_Y ); + ScaLBL_FreeDeviceMemory( sendbuf_z ); + ScaLBL_FreeDeviceMemory( sendbuf_Z ); + ScaLBL_FreeDeviceMemory( sendbuf_xy ); + ScaLBL_FreeDeviceMemory( sendbuf_xY ); + ScaLBL_FreeDeviceMemory( sendbuf_Xy ); + ScaLBL_FreeDeviceMemory( sendbuf_XY ); + ScaLBL_FreeDeviceMemory( sendbuf_xz ); + ScaLBL_FreeDeviceMemory( sendbuf_xZ ); + ScaLBL_FreeDeviceMemory( sendbuf_Xz ); + ScaLBL_FreeDeviceMemory( sendbuf_XZ ); + ScaLBL_FreeDeviceMemory( sendbuf_yz ); + ScaLBL_FreeDeviceMemory( sendbuf_yZ ); + ScaLBL_FreeDeviceMemory( sendbuf_Yz ); + ScaLBL_FreeDeviceMemory( sendbuf_YZ ); + ScaLBL_FreeDeviceMemory( recvbuf_x ); + ScaLBL_FreeDeviceMemory( recvbuf_X ); + ScaLBL_FreeDeviceMemory( recvbuf_y ); + ScaLBL_FreeDeviceMemory( recvbuf_Y ); + ScaLBL_FreeDeviceMemory( recvbuf_z ); + ScaLBL_FreeDeviceMemory( recvbuf_Z ); + ScaLBL_FreeDeviceMemory( recvbuf_xy ); + ScaLBL_FreeDeviceMemory( recvbuf_xY ); + ScaLBL_FreeDeviceMemory( recvbuf_Xy ); + ScaLBL_FreeDeviceMemory( recvbuf_XY ); + ScaLBL_FreeDeviceMemory( recvbuf_xz ); + ScaLBL_FreeDeviceMemory( recvbuf_xZ ); + ScaLBL_FreeDeviceMemory( recvbuf_Xz ); + ScaLBL_FreeDeviceMemory( recvbuf_XZ ); + ScaLBL_FreeDeviceMemory( recvbuf_yz ); + ScaLBL_FreeDeviceMemory( recvbuf_yZ ); + ScaLBL_FreeDeviceMemory( recvbuf_Yz ); + ScaLBL_FreeDeviceMemory( recvbuf_YZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_x ); + ScaLBL_FreeDeviceMemory( dvcSendList_X ); + ScaLBL_FreeDeviceMemory( dvcSendList_y ); + ScaLBL_FreeDeviceMemory( dvcSendList_Y ); + ScaLBL_FreeDeviceMemory( dvcSendList_z ); + ScaLBL_FreeDeviceMemory( dvcSendList_Z ); + ScaLBL_FreeDeviceMemory( dvcSendList_xy ); + ScaLBL_FreeDeviceMemory( dvcSendList_xY ); + ScaLBL_FreeDeviceMemory( dvcSendList_Xy ); + ScaLBL_FreeDeviceMemory( dvcSendList_XY ); + ScaLBL_FreeDeviceMemory( dvcSendList_xz ); + ScaLBL_FreeDeviceMemory( dvcSendList_xZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_Xz ); + ScaLBL_FreeDeviceMemory( dvcSendList_XZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_yz ); + ScaLBL_FreeDeviceMemory( dvcSendList_yZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_Yz ); + ScaLBL_FreeDeviceMemory( dvcSendList_YZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_x ); + ScaLBL_FreeDeviceMemory( dvcRecvList_X ); + ScaLBL_FreeDeviceMemory( dvcRecvList_y ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Y ); + ScaLBL_FreeDeviceMemory( dvcRecvList_z ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Z ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xy ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xY ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Xy ); + ScaLBL_FreeDeviceMemory( dvcRecvList_XY ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Xz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_XZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_yz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_yZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Yz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_YZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_x ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_X ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_y ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Y ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_z ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Z ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xy ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xY ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Xy ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_XY ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Xz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_XZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_yz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_yZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Yz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_YZ ); } void ScaLBL_Communicator::D3Q19_MapRecv(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count, diff --git a/cpu/Extras.cpp b/cpu/Extras.cpp index c2b17cac..05c1175f 100644 --- a/cpu/Extras.cpp +++ b/cpu/Extras.cpp @@ -13,6 +13,13 @@ extern "C" void ScaLBL_AllocateDeviceMemory(void** address, size_t size){ } } + +extern "C" void ScaLBL_FreeDeviceMemory(void* address){ + if ( address != NULL ) + free( address ); +} + + extern "C" void ScaLBL_CopyToDevice(void* dest, const void* source, size_t size){ // cudaMemcpy(dest,source,size,cudaMemcpyHostToDevice); memcpy(dest, source, size); diff --git a/gpu/Extras.cu b/gpu/Extras.cu index 08429d0c..2d27a689 100644 --- a/gpu/Extras.cu +++ b/gpu/Extras.cu @@ -10,6 +10,11 @@ extern "C" void ScaLBL_AllocateDeviceMemory(void** address, size_t size){ } } +extern "C" void ScaLBL_FreeDeviceMemory(void* address){ + if ( address != NULL ) + cudaFree( address ); +} + extern "C" void ScaLBL_CopyToDevice(void* dest, const void* source, size_t size){ cudaMemcpy(dest,source,size,cudaMemcpyHostToDevice); cudaError_t err = cudaGetLastError(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index de7cb14b..77e08021 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -20,7 +20,7 @@ ADD_LBPM_EXECUTABLE( lbpm_inkbottle_pp ) ADD_LBPM_EXECUTABLE( lbpm_plates_pp ) ADD_LBPM_EXECUTABLE( lbpm_squaretube_pp ) ADD_LBPM_EXECUTABLE( lbpm_BlobAnalysis ) -#ADD_LBPM_EXECUTABLE( TestBubble ) +ADD_LBPM_EXECUTABLE( TestBubble ) ADD_LBPM_EXECUTABLE( ComponentLabel ) ADD_LBPM_EXECUTABLE( ColorToBinary ) ADD_LBPM_EXECUTABLE( BlobAnalysis ) @@ -33,7 +33,7 @@ CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_ # Add the tests ADD_LBPM_TEST( pmmc_cylinder ) -#ADD_LBPM_TEST( TestBubble ) +ADD_LBPM_TEST( TestBubble ) ADD_LBPM_TEST( TestTorus ) ADD_LBPM_TEST( TestFluxBC ) ADD_LBPM_TEST( TestInterfaceSpeed ) diff --git a/tests/TestBubble.cpp b/tests/TestBubble.cpp index 727e2c67..e5c5aa4e 100644 --- a/tests/TestBubble.cpp +++ b/tests/TestBubble.cpp @@ -13,2232 +13,1046 @@ #include "IO/Writer.h" #include "ProfilerApp.h" -#define USE_NEW_WRITER - using namespace std; -//************************************************************************* -// Implementation of Two-Phase Immiscible LBM using CUDA -//************************************************************************* -inline void PackID(int *list, int count, char *sendbuf, char *ID){ - // Fill in the phase ID values from neighboring processors - // This packs up the values that need to be sent from one processor to another - int idx,n; - - for (idx=0; idx> tau; // Viscosity parameter - input >> alpha; // Surface Tension parameter - input >> beta; // Width of the interface - input >> phi_s; // value of phi at the solid surface - input >> wp_saturation; - input >> Fx; - input >> Fy; - input >> Fz; - input >> Restart; - input >> pBC; - input >> din; - input >> dout; - input >> timestepMax; // max no. of timesteps - input >> interval; // restart interval - input >> tol; // error tolerance - } - else { - printf("WARNING: No valid Color.in file, using hard-coded values \n"); - tau=1.0; alpha=1e-3; beta=0.95, phi_s=0.0; - wp_saturation=Fx=Fy=Fz=0; - Restart=0; pBC=0; - din=dout=1.0; - timestepMax=100; - interval=100; - tol=1e-6; - das = 0.1; dbs = 0.9; // hard coded for density initialization - // should be OK to remove these parameters - // they should have no impact with the - // current boundary condition - } - //....................................................................... - // Reading the domain information file - //....................................................................... - ifstream domain("Domain.in"); - if (domain.good()){ - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - } - else { - //....................................................................... - // Set the domain for single processor test - printf("WARNING: No valid Domain.in file, using hard-coded values \n"); - nprocx=nprocy=nprocz=1; - Nx=Ny=Nz=50; - nspheres=1; - Lx=Ly=Lz=1; - } + if (rank == 0){ + printf("********************************************************\n"); + printf("Running Hybrid Implementation of Color LBM \n"); + printf("********************************************************\n"); + } - } - // ************************************************************** - // Broadcast simulation parameters from rank 0 to all other procs - MPI_Barrier(comm); - //................................................. - MPI_Bcast(&tau,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&das,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&dbs,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&phi_s,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&wp_saturation,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&pBC,1,MPI_LOGICAL,0,comm); - MPI_Bcast(&Restart,1,MPI_LOGICAL,0,comm); - MPI_Bcast(&din,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&dout,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Fx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Fy,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Fz,1,MPI_DOUBLE,0,comm); - MPI_Bcast(×tepMax,1,MPI_INT,0,comm); - MPI_Bcast(&interval,1,MPI_INT,0,comm); - MPI_Bcast(&tol,1,MPI_DOUBLE,0,comm); - // Computational domain - MPI_Bcast(&Nz,1,MPI_INT,0,comm); -// MPI_Bcast(&nBlocks,1,MPI_INT,0,comm); -// MPI_Bcast(&nthreads,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - 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); + 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 timestepMax, interval; + double tau,Fx,Fy,Fz,tol; + double alpha, beta; + double das, dbs, phi_s; + double din,dout; + double wp_saturation; + bool pBC,Restart; + 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; + + if (rank==0){ + //............................................................. + // READ SIMULATION PARMAETERS FROM INPUT FILE + //............................................................. + ifstream input("Color.in"); + if (input.good()){ + input >> tau; // Viscosity parameter + input >> alpha; // Surface Tension parameter + input >> beta; // Width of the interface + input >> phi_s; // value of phi at the solid surface + input >> wp_saturation; + input >> Fx; + input >> Fy; + input >> Fz; + input >> Restart; + input >> pBC; + input >> din; + input >> dout; + input >> timestepMax; // max no. of timesteps + input >> interval; // restart interval + input >> tol; // error tolerance + } + else { + printf("WARNING: No valid Color.in file, using hard-coded values \n"); + tau = 1.0; + alpha = 1e-3; + beta = 0.95; + phi_s = 0.0; + wp_saturation=Fx=Fy=Fz=0; + Restart = 0; + pBC = 0; + din=dout = 1.0; + timestepMax = 3; + interval = 100; + tol = 1e-6; + das = 0.1; dbs = 0.9; // hard coded for density initialization + // should be OK to remove these parameters + // they should have no impact with the + // current boundary condition + } + //....................................................................... + // Reading the domain information file + //....................................................................... + ifstream domain("Domain.in"); + if (domain.good()){ + domain >> nprocx; + domain >> nprocy; + domain >> nprocz; + domain >> Nx; + domain >> Ny; + domain >> Nz; + domain >> nspheres; + domain >> Lx; + domain >> Ly; + domain >> Lz; + } + else { + //....................................................................... + // Set the domain for single processor test + printf("WARNING: No valid Domain.in file, using hard-coded values \n"); + nprocx=nprocy=nprocz=1; + Nx=Ny=Nz=50; + nspheres=1; + Lx=Ly=Lz=1; + } + + } + // ************************************************************** + // Broadcast simulation parameters from rank 0 to all other procs + MPI_Barrier(comm); + //................................................. + MPI_Bcast(&tau,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&das,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&dbs,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&phi_s,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&wp_saturation,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&pBC,1,MPI_LOGICAL,0,comm); + MPI_Bcast(&Restart,1,MPI_LOGICAL,0,comm); + MPI_Bcast(&din,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&dout,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Fx,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Fy,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Fz,1,MPI_DOUBLE,0,comm); + MPI_Bcast(×tepMax,1,MPI_INT,0,comm); + MPI_Bcast(&interval,1,MPI_INT,0,comm); + MPI_Bcast(&tol,1,MPI_DOUBLE,0,comm); + // Computational domain + MPI_Bcast(&Nx,1,MPI_INT,0,comm); + MPI_Bcast(&Ny,1,MPI_INT,0,comm); + MPI_Bcast(&Nz,1,MPI_INT,0,comm); + // MPI_Bcast(&nBlocks,1,MPI_INT,0,comm); + // MPI_Bcast(&nthreads,1,MPI_INT,0,comm); + MPI_Bcast(&nprocx,1,MPI_INT,0,comm); + MPI_Bcast(&nprocy,1,MPI_INT,0,comm); + MPI_Bcast(&nprocz,1,MPI_INT,0,comm); + MPI_Bcast(&nspheres,1,MPI_INT,0,comm); + MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); + //................................................. + + // 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"); - } + 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"); + } - InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc, - rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z, - rank_xy, rank_XY, rank_xY, rank_Xy, rank_xz, rank_XZ, rank_xZ, rank_Xz, - rank_yz, rank_YZ, rank_yZ, rank_Yz ); - - MPI_Barrier(comm); + // 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); - Nz += 2; - Nx = Ny = Nz; // Cubic domain + // Mask that excludes the solid phase + Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,pBC); + + MPI_Barrier(comm); - int N = Nx*Ny*Nz; - int dist_mem_size = N*sizeof(double); + Nx+=2; Ny+=2; Nz += 2; + //Nx = Ny = Nz; // Cubic domain -// unsigned int nBlocks = 32; -// int nthreads = 128; - int S = N/nthreads/nBlocks+1; + int N = Nx*Ny*Nz; + int dist_mem_size = N*sizeof(double); -// unsigned int nBlocks = N/nthreads + (N%nthreads == 0?0:1); -// dim3 grid(nBlocks,1,1); + int S = N/nthreads/nBlocks+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"); + // unsigned int nBlocks = N/nthreads + (N%nthreads == 0?0:1); + // dim3 grid(nBlocks,1,1); - //....................................................................... - 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; + 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"); - 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 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 nw_pts(20); - DTMutableList ns_pts(20); - DTMutableList ws_pts(20); - DTMutableList 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 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 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 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 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 nw_pts(20); + DTMutableList ns_pts(20); + DTMutableList ws_pts(20); + DTMutableList 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 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 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 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 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; - sendtag = recvtag = 5; - - //************ 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_D3Q19_Pack(1,dvcSendList_x,0,sendCount_x,sendbuf_x,f_even,N); - ScaLBL_D3Q19_Pack(4,dvcSendList_x,sendCount_x,sendCount_x,sendbuf_x,f_even,N); - ScaLBL_D3Q19_Pack(5,dvcSendList_x,2*sendCount_x,sendCount_x,sendbuf_x,f_even,N); - ScaLBL_D3Q19_Pack(6,dvcSendList_x,3*sendCount_x,sendCount_x,sendbuf_x,f_even,N); - ScaLBL_D3Q19_Pack(7,dvcSendList_x,4*sendCount_x,sendCount_x,sendbuf_x,f_even,N); - //...Packing for X face(1,7,9,11,13)................................ - ScaLBL_D3Q19_Pack(0,dvcSendList_X,0,sendCount_X,sendbuf_X,f_odd,N); - ScaLBL_D3Q19_Pack(3,dvcSendList_X,sendCount_X,sendCount_X,sendbuf_X,f_odd,N); - ScaLBL_D3Q19_Pack(4,dvcSendList_X,2*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); - ScaLBL_D3Q19_Pack(5,dvcSendList_X,3*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); - ScaLBL_D3Q19_Pack(6,dvcSendList_X,4*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); - //...Packing for y face(4,8,9,16,18)................................. - ScaLBL_D3Q19_Pack(2,dvcSendList_y,0,sendCount_y,sendbuf_y,f_even,N); - ScaLBL_D3Q19_Pack(4,dvcSendList_y,sendCount_y,sendCount_y,sendbuf_y,f_even,N); - ScaLBL_D3Q19_Pack(4,dvcSendList_y,2*sendCount_y,sendCount_y,sendbuf_y,f_odd,N); - ScaLBL_D3Q19_Pack(8,dvcSendList_y,3*sendCount_y,sendCount_y,sendbuf_y,f_even,N); - ScaLBL_D3Q19_Pack(9,dvcSendList_y,4*sendCount_y,sendCount_y,sendbuf_y,f_even,N); - //...Packing for Y face(3,7,10,15,17)................................. - ScaLBL_D3Q19_Pack(1,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,f_odd,N); - ScaLBL_D3Q19_Pack(3,dvcSendList_Y,sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); - ScaLBL_D3Q19_Pack(5,dvcSendList_Y,2*sendCount_Y,sendCount_Y,sendbuf_Y,f_even,N); - ScaLBL_D3Q19_Pack(7,dvcSendList_Y,3*sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); - ScaLBL_D3Q19_Pack(8,dvcSendList_Y,4*sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); - //...Packing for z face(6,12,13,16,17)................................ - ScaLBL_D3Q19_Pack(3,dvcSendList_z,0,sendCount_z,sendbuf_z,f_even,N); - ScaLBL_D3Q19_Pack(6,dvcSendList_z,sendCount_z,sendCount_z,sendbuf_z,f_even,N); - ScaLBL_D3Q19_Pack(6,dvcSendList_z,2*sendCount_z,sendCount_z,sendbuf_z,f_odd,N); - ScaLBL_D3Q19_Pack(8,dvcSendList_z,3*sendCount_z,sendCount_z,sendbuf_z,f_even,N); - ScaLBL_D3Q19_Pack(8,dvcSendList_z,4*sendCount_z,sendCount_z,sendbuf_z,f_odd,N); - //...Packing for Z face(5,11,14,15,18)................................ - ScaLBL_D3Q19_Pack(2,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,f_odd,N); - ScaLBL_D3Q19_Pack(5,dvcSendList_Z,sendCount_Z,sendCount_Z,sendbuf_Z,f_odd,N); - ScaLBL_D3Q19_Pack(7,dvcSendList_Z,2*sendCount_Z,sendCount_Z,sendbuf_Z,f_even,N); - ScaLBL_D3Q19_Pack(7,dvcSendList_Z,3*sendCount_Z,sendCount_Z,sendbuf_Z,f_odd,N); - ScaLBL_D3Q19_Pack(9,dvcSendList_Z,4*sendCount_Z,sendCount_Z,sendbuf_Z,f_even,N); - //...Pack the xy edge (8)................................ - ScaLBL_D3Q19_Pack(4,dvcSendList_xy,0,sendCount_xy,sendbuf_xy,f_even,N); - //...Pack the Xy edge (9)................................ - ScaLBL_D3Q19_Pack(4,dvcSendList_Xy,0,sendCount_Xy,sendbuf_Xy,f_odd,N); - //...Pack the xY edge (10)................................ - ScaLBL_D3Q19_Pack(5,dvcSendList_xY,0,sendCount_xY,sendbuf_xY,f_even,N); - //...Pack the XY edge (7)................................ - ScaLBL_D3Q19_Pack(3,dvcSendList_XY,0,sendCount_XY,sendbuf_XY,f_odd,N); - //...Pack the xz edge (12)................................ - ScaLBL_D3Q19_Pack(6,dvcSendList_xz,0,sendCount_xz,sendbuf_xz,f_even,N); - //...Pack the xZ edge (14)................................ - ScaLBL_D3Q19_Pack(7,dvcSendList_xZ,0,sendCount_xZ,sendbuf_xZ,f_even,N); - //...Pack the Xz edge (13)................................ - ScaLBL_D3Q19_Pack(6,dvcSendList_Xz,0,sendCount_Xz,sendbuf_Xz,f_odd,N); - //...Pack the XZ edge (11)................................ - ScaLBL_D3Q19_Pack(5,dvcSendList_XZ,0,sendCount_XZ,sendbuf_XZ,f_odd,N); - //...Pack the xz edge (12)................................ - //...Pack the yz edge (16)................................ - ScaLBL_D3Q19_Pack(8,dvcSendList_yz,0,sendCount_yz,sendbuf_yz,f_even,N); - //...Pack the yZ edge (18)................................ - ScaLBL_D3Q19_Pack(9,dvcSendList_yZ,0,sendCount_yZ,sendbuf_yZ,f_even,N); - //...Pack the Yz edge (17)................................ - ScaLBL_D3Q19_Pack(8,dvcSendList_Yz,0,sendCount_Yz,sendbuf_Yz,f_odd,N); - //...Pack the YZ edge (15)................................ - ScaLBL_D3Q19_Pack(7,dvcSendList_YZ,0,sendCount_YZ,sendbuf_YZ,f_odd,N); - //................................................................................... - - //................................................................................... - // Send all the distributions - MPI_Isend(sendbuf_x, 5*sendCount_x,MPI_DOUBLE,rank_x,sendtag,comm,&req1[0]); - MPI_Irecv(recvbuf_X, 5*recvCount_X,MPI_DOUBLE,rank_X,recvtag,comm,&req2[0]); - MPI_Isend(sendbuf_X, 5*sendCount_X,MPI_DOUBLE,rank_X,sendtag,comm,&req1[1]); - MPI_Irecv(recvbuf_x, 5*recvCount_x,MPI_DOUBLE,rank_x,recvtag,comm,&req2[1]); - MPI_Isend(sendbuf_y, 5*sendCount_y,MPI_DOUBLE,rank_y,sendtag,comm,&req1[2]); - MPI_Irecv(recvbuf_Y, 5*recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,comm,&req2[2]); - MPI_Isend(sendbuf_Y, 5*sendCount_Y,MPI_DOUBLE,rank_Y,sendtag,comm,&req1[3]); - MPI_Irecv(recvbuf_y, 5*recvCount_y,MPI_DOUBLE,rank_y,recvtag,comm,&req2[3]); - MPI_Isend(sendbuf_z, 5*sendCount_z,MPI_DOUBLE,rank_z,sendtag,comm,&req1[4]); - MPI_Irecv(recvbuf_Z, 5*recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,comm,&req2[4]); - MPI_Isend(sendbuf_Z, 5*sendCount_Z,MPI_DOUBLE,rank_Z,sendtag,comm,&req1[5]); - MPI_Irecv(recvbuf_z, 5*recvCount_z,MPI_DOUBLE,rank_z,recvtag,comm,&req2[5]); - MPI_Isend(sendbuf_xy, sendCount_xy,MPI_DOUBLE,rank_xy,sendtag,comm,&req1[6]); - MPI_Irecv(recvbuf_XY, recvCount_XY,MPI_DOUBLE,rank_XY,recvtag,comm,&req2[6]); - MPI_Isend(sendbuf_XY, sendCount_XY,MPI_DOUBLE,rank_XY,sendtag,comm,&req1[7]); - MPI_Irecv(recvbuf_xy, recvCount_xy,MPI_DOUBLE,rank_xy,recvtag,comm,&req2[7]); - MPI_Isend(sendbuf_Xy, sendCount_Xy,MPI_DOUBLE,rank_Xy,sendtag,comm,&req1[8]); - MPI_Irecv(recvbuf_xY, recvCount_xY,MPI_DOUBLE,rank_xY,recvtag,comm,&req2[8]); - MPI_Isend(sendbuf_xY, sendCount_xY,MPI_DOUBLE,rank_xY,sendtag,comm,&req1[9]); - MPI_Irecv(recvbuf_Xy, recvCount_Xy,MPI_DOUBLE,rank_Xy,recvtag,comm,&req2[9]); - MPI_Isend(sendbuf_xz, sendCount_xz,MPI_DOUBLE,rank_xz,sendtag,comm,&req1[10]); - MPI_Irecv(recvbuf_XZ, recvCount_XZ,MPI_DOUBLE,rank_XZ,recvtag,comm,&req2[10]); - MPI_Isend(sendbuf_XZ, sendCount_XZ,MPI_DOUBLE,rank_XZ,sendtag,comm,&req1[11]); - MPI_Irecv(recvbuf_xz, recvCount_xz,MPI_DOUBLE,rank_xz,recvtag,comm,&req2[11]); - MPI_Isend(sendbuf_Xz, sendCount_Xz,MPI_DOUBLE,rank_Xz,sendtag,comm,&req1[12]); - MPI_Irecv(recvbuf_xZ, recvCount_xZ,MPI_DOUBLE,rank_xZ,recvtag,comm,&req2[12]); - MPI_Isend(sendbuf_xZ, sendCount_xZ,MPI_DOUBLE,rank_xZ,sendtag,comm,&req1[13]); - MPI_Irecv(recvbuf_Xz, recvCount_Xz,MPI_DOUBLE,rank_Xz,recvtag,comm,&req2[13]); - MPI_Isend(sendbuf_yz, sendCount_yz,MPI_DOUBLE,rank_yz,sendtag,comm,&req1[14]); - MPI_Irecv(recvbuf_YZ, recvCount_YZ,MPI_DOUBLE,rank_YZ,recvtag,comm,&req2[14]); - MPI_Isend(sendbuf_YZ, sendCount_YZ,MPI_DOUBLE,rank_YZ,sendtag,comm,&req1[15]); - MPI_Irecv(recvbuf_yz, recvCount_yz,MPI_DOUBLE,rank_yz,recvtag,comm,&req2[15]); - MPI_Isend(sendbuf_Yz, sendCount_Yz,MPI_DOUBLE,rank_Yz,sendtag,comm,&req1[16]); - MPI_Irecv(recvbuf_yZ, recvCount_yZ,MPI_DOUBLE,rank_yZ,recvtag,comm,&req2[16]); - MPI_Isend(sendbuf_yZ, sendCount_yZ,MPI_DOUBLE,rank_yZ,sendtag,comm,&req1[17]); - MPI_Irecv(recvbuf_Yz, recvCount_Yz,MPI_DOUBLE,rank_Yz,recvtag,comm,&req2[17]); - //................................................................................... - - //************************************************************************* - // Carry out the density streaming step for mass transport - //************************************************************************* - // DensityStreamD3Q7(ID, Den, Copy, Phi, ColorGrad, Velocity, beta, Nx, Ny, Nz, pBC); - //************************************************************************* - ScaLBL_D3Q7_ColorCollideMass(ID, A_even, A_odd, B_even, B_odd, Den, Phi, - ColorGrad, Velocity, beta, N, pBC); - - - //************************************************************************* - // Swap the distributions for momentum transport - //************************************************************************* - ScaLBL_D3Q19_Swap(ID, f_even, f_odd, Nx, Ny, Nz); - //************************************************************************* - - //................................................................................... - // Wait for completion of D3Q19 communication - MPI_Waitall(18,req1,stat1); - MPI_Waitall(18,req2,stat2); - - //................................................................................... - // Unpack the distributions on the device - //................................................................................... - //...Map recieve list for the X face: q=2,8,10,12,13 ................................. - ScaLBL_D3Q19_Unpack(0,-1,0,0,dvcRecvList_X,0,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(3,-1,-1,0,dvcRecvList_X,recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(4,-1,1,0,dvcRecvList_X,2*recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(5,-1,0,-1,dvcRecvList_X,3*recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(6,-1,0,1,dvcRecvList_X,4*recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); - //................................................................................... - //...Map recieve list for the x face: q=1,7,9,11,13.................................. - ScaLBL_D3Q19_Unpack(1,1,0,0,dvcRecvList_x,0,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(4,1,1,0,dvcRecvList_x,recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(5,1,-1,0,dvcRecvList_x,2*recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(6,1,0,1,dvcRecvList_x,3*recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(7,1,0,-1,dvcRecvList_x,4*recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); - //................................................................................... - //...Map recieve list for the y face: q=4,8,9,16,18 ................................... - ScaLBL_D3Q19_Unpack(1,0,-1,0,dvcRecvList_Y,0,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(3,-1,-1,0,dvcRecvList_Y,recvCount_Y,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(5,1,-1,0,dvcRecvList_Y,2*recvCount_Y,recvCount_Y,recvbuf_Y,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(7,0,-1,-1,dvcRecvList_Y,3*recvCount_Y,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(8,0,-1,1,dvcRecvList_Y,4*recvCount_Y,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); - //................................................................................... - //...Map recieve list for the Y face: q=3,7,10,15,17 .................................. - ScaLBL_D3Q19_Unpack(2,0,1,0,dvcRecvList_y,0,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(4,1,1,0,dvcRecvList_y,recvCount_y,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(4,-1,1,0,dvcRecvList_y,2*recvCount_y,recvCount_y,recvbuf_y,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(8,0,1,1,dvcRecvList_y,3*recvCount_y,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(9,0,1,-1,dvcRecvList_y,4*recvCount_y,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); - //................................................................................... - //...Map recieve list for the z face<<<6,12,13,16,17).............................................. - ScaLBL_D3Q19_Unpack(2,0,0,-1,dvcRecvList_Z,0,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(5,-1,0,-1,dvcRecvList_Z,recvCount_Z,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(7,1,0,-1,dvcRecvList_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(7,0,-1,-1,dvcRecvList_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(9,0,1,-1,dvcRecvList_Z,4*recvCount_Z,recvCount_Z,recvbuf_Z,f_even,Nx,Ny,Nz); - //...Map recieve list for the Z face<<<5,11,14,15,18).............................................. - ScaLBL_D3Q19_Unpack(3,0,0,1,dvcRecvList_z,0,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(6,1,0,1,dvcRecvList_z,recvCount_z,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(6,-1,0,1,dvcRecvList_z,2*recvCount_z,recvCount_z,recvbuf_z,f_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(8,0,1,1,dvcRecvList_z,3*recvCount_z,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(8,0,-1,1,dvcRecvList_z,4*recvCount_z,recvCount_z,recvbuf_z,f_odd,Nx,Ny,Nz); - //.................................................................................. - //...Map recieve list for the xy edge <<<8)................................ - ScaLBL_D3Q19_Unpack(3,-1,-1,0,dvcRecvList_XY,0,recvCount_XY,recvbuf_XY,f_odd,Nx,Ny,Nz); - //...Map recieve list for the Xy edge <<<9)................................ - ScaLBL_D3Q19_Unpack(5,1,-1,0,dvcRecvList_xY,0,recvCount_xY,recvbuf_xY,f_even,Nx,Ny,Nz); - //...Map recieve list for the xY edge <<<10)................................ - ScaLBL_D3Q19_Unpack(4,-1,1,0,dvcRecvList_Xy,0,recvCount_Xy,recvbuf_Xy,f_odd,Nx,Ny,Nz); - //...Map recieve list for the XY edge <<<7)................................ - ScaLBL_D3Q19_Unpack(4,1,1,0,dvcRecvList_xy,0,recvCount_xy,recvbuf_xy,f_even,Nx,Ny,Nz); - //...Map recieve list for the xz edge <<<12)................................ - ScaLBL_D3Q19_Unpack(5,-1,0,-1,dvcRecvList_XZ,0,recvCount_XZ,recvbuf_XZ,f_odd,Nx,Ny,Nz); - //...Map recieve list for the xZ edge <<<14)................................ - ScaLBL_D3Q19_Unpack(6,-1,0,1,dvcRecvList_Xz,0,recvCount_Xz,recvbuf_Xz,f_odd,Nx,Ny,Nz); - //...Map recieve list for the Xz edge <<<13)................................ - ScaLBL_D3Q19_Unpack(7,1,0,-1,dvcRecvList_xZ,0,recvCount_xZ,recvbuf_xZ,f_even,Nx,Ny,Nz); - //...Map recieve list for the XZ edge <<<11)................................ - ScaLBL_D3Q19_Unpack(6,1,0,1,dvcRecvList_xz,0,recvCount_xz,recvbuf_xz,f_even,Nx,Ny,Nz); - //...Map recieve list for the yz edge <<<16)................................ - ScaLBL_D3Q19_Unpack(7,0,-1,-1,dvcRecvList_YZ,0,recvCount_YZ,recvbuf_YZ,f_odd,Nx,Ny,Nz); - //...Map recieve list for the yZ edge <<<18)................................ - ScaLBL_D3Q19_Unpack(8,0,-1,1,dvcRecvList_Yz,0,recvCount_Yz,recvbuf_Yz,f_odd,Nx,Ny,Nz); - //...Map recieve list for the Yz edge <<<17)................................ - ScaLBL_D3Q19_Unpack(9,0,1,-1,dvcRecvList_yZ,0,recvCount_yZ,recvbuf_yZ,f_even,Nx,Ny,Nz); - //...Map recieve list for the YZ edge <<<15)................................ - ScaLBL_D3Q19_Unpack(8,0,1,1,dvcRecvList_yz,0,recvCount_yz,recvbuf_yz,f_even,Nx,Ny,Nz); - //................................................................................... - - //................................................................................... - ScaLBL_D3Q19_Pack(1,dvcSendList_x,0,sendCount_x,sendbuf_x,A_even,N); - ScaLBL_D3Q19_Pack(1,dvcSendList_x,sendCount_x,sendCount_x,sendbuf_x,B_even,N); - //...Packing for X face(1,7,9,11,13)................................ - ScaLBL_D3Q19_Pack(0,dvcSendList_X,0,sendCount_X,sendbuf_X,A_odd,N); - ScaLBL_D3Q19_Pack(0,dvcSendList_X,sendCount_X,sendCount_X,sendbuf_X,B_odd,N); - //...Packing for y face(4,8,9,16,18)................................. - ScaLBL_D3Q19_Pack(2,dvcSendList_y,0,sendCount_y,sendbuf_y,A_even,N); - ScaLBL_D3Q19_Pack(2,dvcSendList_y,sendCount_y,sendCount_y,sendbuf_y,B_even,N); - //...Packing for Y face(3,7,10,15,17)................................. - ScaLBL_D3Q19_Pack(1,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,A_odd,N); - ScaLBL_D3Q19_Pack(1,dvcSendList_Y,sendCount_Y,sendCount_Y,sendbuf_Y,B_odd,N); - //...Packing for z face(6,12,13,16,17)................................ - ScaLBL_D3Q19_Pack(3,dvcSendList_z,0,sendCount_z,sendbuf_z,A_even,N); - ScaLBL_D3Q19_Pack(3,dvcSendList_z,sendCount_z,sendCount_z,sendbuf_z,B_even,N); - //...Packing for Z face(5,11,14,15,18)................................ - ScaLBL_D3Q19_Pack(2,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,A_odd,N); - ScaLBL_D3Q19_Pack(2,dvcSendList_Z,sendCount_Z,sendCount_Z,sendbuf_Z,B_odd,N); - //................................................................................... - - //................................................................................... - // Send all the distributions - MPI_Isend(sendbuf_x, 2*sendCount_x,MPI_DOUBLE,rank_x,sendtag,comm,&req1[0]); - MPI_Irecv(recvbuf_X, 2*recvCount_X,MPI_DOUBLE,rank_X,recvtag,comm,&req2[0]); - MPI_Isend(sendbuf_X, 2*sendCount_X,MPI_DOUBLE,rank_X,sendtag,comm,&req1[1]); - MPI_Irecv(recvbuf_x, 2*recvCount_x,MPI_DOUBLE,rank_x,recvtag,comm,&req2[1]); - MPI_Isend(sendbuf_y, 2*sendCount_y,MPI_DOUBLE,rank_y,sendtag,comm,&req1[2]); - MPI_Irecv(recvbuf_Y, 2*recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,comm,&req2[2]); - MPI_Isend(sendbuf_Y, 2*sendCount_Y,MPI_DOUBLE,rank_Y,sendtag,comm,&req1[3]); - MPI_Irecv(recvbuf_y, 2*recvCount_y,MPI_DOUBLE,rank_y,recvtag,comm,&req2[3]); - MPI_Isend(sendbuf_z, 2*sendCount_z,MPI_DOUBLE,rank_z,sendtag,comm,&req1[4]); - MPI_Irecv(recvbuf_Z, 2*recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,comm,&req2[4]); - MPI_Isend(sendbuf_Z, 2*sendCount_Z,MPI_DOUBLE,rank_Z,sendtag,comm,&req1[5]); - MPI_Irecv(recvbuf_z, 2*recvCount_z,MPI_DOUBLE,rank_z,recvtag,comm,&req2[5]); - //................................................................................... - - ScaLBL_D3Q7_Swap(ID, A_even, A_odd, Nx, Ny, Nz); - ScaLBL_D3Q7_Swap(ID, B_even, B_odd, Nx, Ny, Nz); - - //................................................................................... - // Wait for completion of D3Q19 communication - MPI_Waitall(6,req1,stat1); - MPI_Waitall(6,req2,stat2); - //................................................................................... - // Unpack the distributions on the device - //................................................................................... - //...Map recieve list for the X face: q=2,8,10,12,13 ................................. - ScaLBL_D3Q19_Unpack(0,-1,0,0,dvcRecvList_X,0,recvCount_X,recvbuf_X,A_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(0,-1,0,0,dvcRecvList_X,recvCount_X,recvCount_X,recvbuf_X,B_odd,Nx,Ny,Nz); - //................................................................................... - //...Map recieve list for the x face: q=1,7,9,11,13.................................. - ScaLBL_D3Q19_Unpack(1,1,0,0,dvcRecvList_x,0,recvCount_x,recvbuf_x,A_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(1,1,0,0,dvcRecvList_x,recvCount_x,recvCount_x,recvbuf_x,B_even,Nx,Ny,Nz); - //................................................................................... - //...Map recieve list for the y face: q=4,8,9,16,18 ................................... - ScaLBL_D3Q19_Unpack(1,0,-1,0,dvcRecvList_Y,0,recvCount_Y,recvbuf_Y,A_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(1,0,-1,0,dvcRecvList_Y,recvCount_Y,recvCount_Y,recvbuf_Y,B_odd,Nx,Ny,Nz); - //................................................................................... - //...Map recieve list for the Y face: q=3,7,10,15,17 .................................. - ScaLBL_D3Q19_Unpack(2,0,1,0,dvcRecvList_y,0,recvCount_y,recvbuf_y,A_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(2,0,1,0,dvcRecvList_y,recvCount_y,recvCount_y,recvbuf_y,B_even,Nx,Ny,Nz); - //................................................................................... - //...Map recieve list for the z face<<<6,12,13,16,17).............................................. - ScaLBL_D3Q19_Unpack(2,0,0,-1,dvcRecvList_Z,0,recvCount_Z,recvbuf_Z,A_odd,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(2,0,0,-1,dvcRecvList_Z,recvCount_Z,recvCount_Z,recvbuf_Z,B_odd,Nx,Ny,Nz); - //...Map recieve list for the Z face<<<5,11,14,15,18).............................................. - ScaLBL_D3Q19_Unpack(3,0,0,1,dvcRecvList_z,0,recvCount_z,recvbuf_z,A_even,Nx,Ny,Nz); - ScaLBL_D3Q19_Unpack(3,0,0,1,dvcRecvList_z,recvCount_z,recvCount_z,recvbuf_z,B_even,Nx,Ny,Nz); - //.................................................................................. - - //.................................................................................. - 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_ComputePhaseField(ID, Phi, Den, N); - //************************************************************************* - - //................................................................................... - ScaLBL_Scalar_Pack(dvcSendList_x, sendCount_x,sendbuf_x, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_y, sendCount_y,sendbuf_y, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_z, sendCount_z,sendbuf_z, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_X, sendCount_X,sendbuf_X, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_Y, sendCount_Y,sendbuf_Y, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_Z, sendCount_Z,sendbuf_Z, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_xy, sendCount_xy,sendbuf_xy, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_xY, sendCount_xY,sendbuf_xY, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_Xy, sendCount_Xy,sendbuf_Xy, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_XY, sendCount_XY,sendbuf_XY, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_xz, sendCount_xz,sendbuf_xz, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_xZ, sendCount_xZ,sendbuf_xZ, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_Xz, sendCount_Xz,sendbuf_Xz, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_XZ, sendCount_XZ,sendbuf_XZ, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_yz, sendCount_yz,sendbuf_yz, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_yZ, sendCount_yZ,sendbuf_yZ, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_Yz, sendCount_Yz,sendbuf_Yz, Phi, N); - ScaLBL_Scalar_Pack(dvcSendList_YZ, sendCount_YZ,sendbuf_YZ, Phi, N); - //................................................................................... - // Send / Recv all the phase indcator field values - //................................................................................... - MPI_Isend(sendbuf_x, sendCount_x,MPI_DOUBLE,rank_x,sendtag,comm,&req1[0]); - MPI_Irecv(recvbuf_X, recvCount_X,MPI_DOUBLE,rank_X,recvtag,comm,&req2[0]); - MPI_Isend(sendbuf_X, sendCount_X,MPI_DOUBLE,rank_X,sendtag,comm,&req1[1]); - MPI_Irecv(recvbuf_x, recvCount_x,MPI_DOUBLE,rank_x,recvtag,comm,&req2[1]); - MPI_Isend(sendbuf_y, sendCount_y,MPI_DOUBLE,rank_y,sendtag,comm,&req1[2]); - MPI_Irecv(recvbuf_Y, recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,comm,&req2[2]); - MPI_Isend(sendbuf_Y, sendCount_Y,MPI_DOUBLE,rank_Y,sendtag,comm,&req1[3]); - MPI_Irecv(recvbuf_y, recvCount_y,MPI_DOUBLE,rank_y,recvtag,comm,&req2[3]); - MPI_Isend(sendbuf_z, sendCount_z,MPI_DOUBLE,rank_z,sendtag,comm,&req1[4]); - MPI_Irecv(recvbuf_Z, recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,comm,&req2[4]); - MPI_Isend(sendbuf_Z, sendCount_Z,MPI_DOUBLE,rank_Z,sendtag,comm,&req1[5]); - MPI_Irecv(recvbuf_z, recvCount_z,MPI_DOUBLE,rank_z,recvtag,comm,&req2[5]); - MPI_Isend(sendbuf_xy, sendCount_xy,MPI_DOUBLE,rank_xy,sendtag,comm,&req1[6]); - MPI_Irecv(recvbuf_XY, recvCount_XY,MPI_DOUBLE,rank_XY,recvtag,comm,&req2[6]); - MPI_Isend(sendbuf_XY, sendCount_XY,MPI_DOUBLE,rank_XY,sendtag,comm,&req1[7]); - MPI_Irecv(recvbuf_xy, recvCount_xy,MPI_DOUBLE,rank_xy,recvtag,comm,&req2[7]); - MPI_Isend(sendbuf_Xy, sendCount_Xy,MPI_DOUBLE,rank_Xy,sendtag,comm,&req1[8]); - MPI_Irecv(recvbuf_xY, recvCount_xY,MPI_DOUBLE,rank_xY,recvtag,comm,&req2[8]); - MPI_Isend(sendbuf_xY, sendCount_xY,MPI_DOUBLE,rank_xY,sendtag,comm,&req1[9]); - MPI_Irecv(recvbuf_Xy, recvCount_Xy,MPI_DOUBLE,rank_Xy,recvtag,comm,&req2[9]); - MPI_Isend(sendbuf_xz, sendCount_xz,MPI_DOUBLE,rank_xz,sendtag,comm,&req1[10]); - MPI_Irecv(recvbuf_XZ, recvCount_XZ,MPI_DOUBLE,rank_XZ,recvtag,comm,&req2[10]); - MPI_Isend(sendbuf_XZ, sendCount_XZ,MPI_DOUBLE,rank_XZ,sendtag,comm,&req1[11]); - MPI_Irecv(recvbuf_xz, recvCount_xz,MPI_DOUBLE,rank_xz,recvtag,comm,&req2[11]); - MPI_Isend(sendbuf_Xz, sendCount_Xz,MPI_DOUBLE,rank_Xz,sendtag,comm,&req1[12]); - MPI_Irecv(recvbuf_xZ, recvCount_xZ,MPI_DOUBLE,rank_xZ,recvtag,comm,&req2[12]); - MPI_Isend(sendbuf_xZ, sendCount_xZ,MPI_DOUBLE,rank_xZ,sendtag,comm,&req1[13]); - MPI_Irecv(recvbuf_Xz, recvCount_Xz,MPI_DOUBLE,rank_Xz,recvtag,comm,&req2[13]); - MPI_Isend(sendbuf_yz, sendCount_yz,MPI_DOUBLE,rank_yz,sendtag,comm,&req1[14]); - MPI_Irecv(recvbuf_YZ, recvCount_YZ,MPI_DOUBLE,rank_YZ,recvtag,comm,&req2[14]); - MPI_Isend(sendbuf_YZ, sendCount_YZ,MPI_DOUBLE,rank_YZ,sendtag,comm,&req1[15]); - MPI_Irecv(recvbuf_yz, recvCount_yz,MPI_DOUBLE,rank_yz,recvtag,comm,&req2[15]); - MPI_Isend(sendbuf_Yz, sendCount_Yz,MPI_DOUBLE,rank_Yz,sendtag,comm,&req1[16]); - MPI_Irecv(recvbuf_yZ, recvCount_yZ,MPI_DOUBLE,rank_yZ,recvtag,comm,&req2[16]); - MPI_Isend(sendbuf_yZ, sendCount_yZ,MPI_DOUBLE,rank_yZ,sendtag,comm,&req1[17]); - MPI_Irecv(recvbuf_Yz, recvCount_Yz,MPI_DOUBLE,rank_Yz,recvtag,comm,&req2[17]); - //................................................................................... - //................................................................................... - // Wait for completion of Indicator Field communication - //................................................................................... - MPI_Waitall(18,req1,stat1); - MPI_Waitall(18,req2,stat2); - ScaLBL_DeviceBarrier(); - //................................................................................... - //................................................................................... - /* ScaLBL_Scalar_Unpack(faceGrid, packThreads, dvcSendList_x, sendCount_x,sendbuf_x, Phi, N); - ScaLBL_Scalar_Unpack(faceGrid, packThreads, dvcSendList_y, sendCount_y,sendbuf_y, Phi, N); - ScaLBL_Scalar_Unpack(faceGrid, packThreads, dvcSendList_z, sendCount_z,sendbuf_z, Phi, N); - ScaLBL_Scalar_Unpack(faceGrid, packThreads, dvcSendList_X, sendCount_X,sendbuf_X, Phi, N); - ScaLBL_Scalar_Unpack(faceGrid, packThreads, dvcSendList_Y, sendCount_Y,sendbuf_Y, Phi, N); - ScaLBL_Scalar_Unpack(faceGrid, packThreads, dvcSendList_Z, sendCount_Z,sendbuf_Z, Phi, N); - */ - ScaLBL_Scalar_Unpack(dvcRecvList_x, recvCount_x,recvbuf_x, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_y, recvCount_y,recvbuf_y, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_z, recvCount_z,recvbuf_z, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_X, recvCount_X,recvbuf_X, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_Y, recvCount_Y,recvbuf_Y, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_Z, recvCount_Z,recvbuf_Z, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_xy, recvCount_xy,recvbuf_xy, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_xY, recvCount_xY,recvbuf_xY, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_Xy, recvCount_Xy,recvbuf_Xy, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_XY, recvCount_XY,recvbuf_XY, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_xz, recvCount_xz,recvbuf_xz, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_xZ, recvCount_xZ,recvbuf_xZ, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_Xz, recvCount_Xz,recvbuf_Xz, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_XZ, recvCount_XZ,recvbuf_XZ, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_yz, recvCount_yz,recvbuf_yz, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_yZ, recvCount_yZ,recvbuf_yZ, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_Yz, recvCount_Yz,recvbuf_Yz, Phi, N); - ScaLBL_Scalar_Unpack(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, Phi, N); - //................................................................................... - - if (pBC && kproc == 0) { - ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); - ScaLBL_Color_BC_z(Phi,Den,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,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 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 - } - - - #ifdef USE_NEW_WRITER - shared_ptr mesh( new IO::TriList() ); - mesh->A.reserve(8*ncubes); - mesh->B.reserve(8*ncubes); - mesh->C.reserve(8*ncubes); - #else - FILE *WN_TRIS; - sprintf(LocalRankFilename,"%s%s","wn-tris.",LocalRankString); - WN_TRIS = fopen(LocalRankFilename,"wb"); - #endif - for (c=0;cA.push_back(A); - mesh->B.push_back(B); - mesh->C.push_back(C); - #else - //fprintf(WN_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z); - fwrite(&A.x,sizeof(A.x),1,WN_TRIS); - fwrite(&A.y,sizeof(A.y),1,WN_TRIS); - fwrite(&A.z,sizeof(A.z),1,WN_TRIS); - fwrite(&B.x,sizeof(B.x),1,WN_TRIS); - fwrite(&B.y,sizeof(B.y),1,WN_TRIS); - fwrite(&B.z,sizeof(B.z),1,WN_TRIS); - fwrite(&C.x,sizeof(C.x),1,WN_TRIS); - fwrite(&C.y,sizeof(C.y),1,WN_TRIS); - fwrite(&C.z,sizeof(C.z),1,WN_TRIS); - #endif - } - } - #ifdef USE_NEW_WRITER - std::vector meshData(1); - meshData[0].meshName = "wn-tris"; - meshData[0].mesh = mesh; - for (size_t k=0; k 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; iA.size(); i++) { - const Point& a = mesh->A[i]; - const Point& b = mesh->B[i]; - const Point& c = mesh->C[i]; - dist->data(0,i) = sqrt(a.x*a.x+a.y*a.y+a.z*a.z); - dist->data(1,i) = sqrt(b.x*b.x+b.y*b.y+b.z*b.z); - dist->data(2,i) = sqrt(c.x*c.x+c.y*c.y+c.z*c.z); + // Initialize the bubble + for (k=0;k 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_Color_BC_z(Phi,Den,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,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 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 mesh( new IO::TriList() ); + mesh->A.reserve(8*ncubes); + mesh->B.reserve(8*ncubes); + mesh->C.reserve(8*ncubes); + for (c=0;cA.push_back(A); + mesh->B.push_back(B); + mesh->C.push_back(C); + } + } + std::vector meshData(1); + meshData[0].meshName = "wn-tris"; + meshData[0].mesh = mesh; + for (size_t k=0; k 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; iA.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)); + //************************************************************************/ + 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); -*/ //************************************************************************/ + 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); - MPI_Finalize(); - // **************************************************** - return 0; + + // **************************************************** + MPI_Barrier(comm); + } // Limit scope so variables that contain communicators will free before MPI_Finialize + MPI_Comm_free(&comm); + MPI_Finalize(); + return 0; } diff --git a/tests/TestFluxBC.cpp b/tests/TestFluxBC.cpp index e0302a58..6d9b1ec9 100644 --- a/tests/TestFluxBC.cpp +++ b/tests/TestFluxBC.cpp @@ -117,7 +117,7 @@ int main (int argc, char **argv) printf("Copying velocity to host \n"); double * vel; - vel = new double [N]; + vel = new double [3*N]; ScaLBL_CopyToHost(vel,dvc_vel,3*N*sizeof(double)); // Check the first layer diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index 409debb7..f47b0b73 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -113,7 +113,6 @@ int main(int argc, char **argv) // parallel domain size (# of sub-domains) int nprocx,nprocy,nprocz; - int iproc,jproc,kproc; MPI_Request req1[18],req2[18]; MPI_Status stat1[18],stat2[18]; @@ -268,7 +267,6 @@ int main(int argc, char **argv) // Get the rank info const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); - MPI_Barrier(comm); // **************************************************************