From 9383ffb106ef7b0afdcd005456c05c0270768729 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 24 May 2014 13:29:32 -0400 Subject: [PATCH 01/43] created files to contain thermal LBM --- cpu/thermal.cpp | 2 ++ cpu/thermal.h | 2 ++ 2 files changed, 4 insertions(+) create mode 100644 cpu/thermal.cpp create mode 100644 cpu/thermal.h diff --git a/cpu/thermal.cpp b/cpu/thermal.cpp new file mode 100644 index 00000000..29e7d2e3 --- /dev/null +++ b/cpu/thermal.cpp @@ -0,0 +1,2 @@ +// cpu implementation for thermal lattice boltzmann methods +// copyright James McClure, 2014 diff --git a/cpu/thermal.h b/cpu/thermal.h new file mode 100644 index 00000000..15ecde06 --- /dev/null +++ b/cpu/thermal.h @@ -0,0 +1,2 @@ +// header for thermal lattice boltzmann methods +// copyright James McClure, 2014 From 306b31e0a2c6a2ddffc7eadef40eb929c1449d58 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 11 Jun 2014 13:52:51 -0400 Subject: [PATCH 02/43] minor change to cpu/thermal.cpp --- cpu/thermal.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/cpu/thermal.cpp b/cpu/thermal.cpp index 29e7d2e3..24a7de81 100644 --- a/cpu/thermal.cpp +++ b/cpu/thermal.cpp @@ -1,2 +1,5 @@ // cpu implementation for thermal lattice boltzmann methods // copyright James McClure, 2014 + + + From a3bdb6efc18ec055247672e2f1ac631ed672bba3 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 10:22:01 -0400 Subject: [PATCH 03/43] Basic things to tests/lb2_Blob_wia_mpi.cpp --- tests/lb2_Blob_wia_mpi.cpp | 2832 ++++++++++++++++++++++++++++++++++++ 1 file changed, 2832 insertions(+) create mode 100644 tests/lb2_Blob_wia_mpi.cpp diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp new file mode 100644 index 00000000..12bfed69 --- /dev/null +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -0,0 +1,2832 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include "pmmc.h" +#include "Domain.h" +#include "Extras.h" +#include "D3Q19.h" +#include "D3Q7.h" +#include "Color.h" +#include "Communication.h" + +//#define CBUB +#define USE_EXP_CONTACT_ANGLE +#define WRITE_SURFACES +#define MAX_LOCAL_BLOB_COUNT 500 + +using namespace std; + + +inline int ComputeLocalBlob(IntArray blobs, int &nblobs, int &ncubes, IntArray indicator, + DoubleArray F, DoubleArray S, double vf, double vs, int startx, int starty, + int startz, IntArray temp) +{ + // Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob + // F>vf => oil phase S>vs => in porespace + // update the list of blobs, indicator mesh + int m = F.m; // maxima for the meshes + int n = F.n; + int o = F.o; + + int cubes_in_blob=0; + int nrecent = 1; // number of nodes added at most recent sweep + temp(0,0) = startx; // Set the initial point as a "seed" for the sweeps + temp(1,0) = starty; + temp(2,0) = startz; + int ntotal = 1; // total number of nodes in blob + indicator(startx,starty,startz) = nblobs; + + int p,s,x,y,z,start,finish,nodx,nody,nodz; + int imin=startx,imax=startx,jmin=starty,jmax=starty; // initialize maxima / minima + int kmin=startz,kmax=startz; + int d[26][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}, + {1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0},{1,0,1},{-1,0,1}, + {1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}, + {1,1,1},{1,1,-1},{1,-1,1},{1,-1,-1},{-1,1,1},{-1,1,-1}, + {-1,-1,1},{-1,-1,-1}}; // directions to neighbors + 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 + bool status = 1; // status == true => continue to look for points + while (status == 1){ + start = ntotal - nrecent; + finish = ntotal; + nrecent = 0; // set recent points back to zero for next sweep through + for (s=start;s m-1 ){ nodx = 0; } + nody=y+d[p][1]; + if (nody < 0 ){ nody = n-1; } // Periodic BC for y + if (nody > n-1 ){ nody = 0; } + nodz=z+d[p][2]; + if (nodz < 0 ){ nodz = 0; } // No periodic BC for z + if (nodz > o-1 ){ nodz = o-1; } + if ( F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs + && indicator(nodx,nody,nodz) == -1 ){ + // Node is a part of the blob - add it to the list + temp(0,ntotal) = nodx; + temp(1,ntotal) = nody; + temp(2,ntotal) = nodz; + ntotal++; + nrecent++; + // Update the indicator map + indicator(nodx,nody,nodz) = nblobs; + // Update the min / max for the cube loop + if ( nodx < imin ){ imin = nodx; } + if ( nodx > imax ){ imax = nodx; } + if ( nody < jmin ){ jmin = nody; } + if ( nody > jmax ){ jmax = nody; } + if ( nodz < kmin ){ kmin = nodz; } + if ( nodz > kmax ){ kmax = nodz; } + } + else if (F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs + && indicator(nodx,nody,nodz) > -1 && indicator(nodx,nody,nodz) != nblobs){ + // Some kind of error in algorithm + printf("Error in blob search algorithm!"); + } + } + + } + if ( nrecent == 0){ + status = 0; + } + } + // Use points in temporary storage array to add cubes to the list of blobs + if ( imin > 0) { imin = imin-1; } +// if ( imax < m-1) { imax = imax+1; } + if ( jmin > 0) { jmin = jmin-1; } +// if ( jmax < n-1) { jmax = jmax+1; } + if ( kmin > 0) { kmin = kmin-1; } +// if ( kmax < o-1) { kmax = kmax+1; } + int i,j,k; + bool add; + for (k=kmin;k -1 && indicator(nodx,nody,nodz) != nblobs){ + printf("Overlapping cube!"); + cout << i << ", " << j << ", " << k << endl; + } + } + } + } + } + } + + return cubes_in_blob; +} + +//************************************************************************* +// 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> FILENAME; + // Line 2: domain size (Nx, Ny, Nz) +// input >> Nz; // number of nodes (x,y,z) +// input >> nBlocks; +// input >> nthreads; + // Line 3: model parameters (tau, alpha, beta, das, dbs) + 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 >> das; +// input >> dbs; + // Line 4: wetting phase saturation to initialize + input >> wp_saturation; + // Line 5: External force components (Fx,Fy, Fz) + input >> Fx; + input >> Fy; + input >> Fz; + // Line 6: Pressure Boundary conditions + input >> Restart; + input >> pBC; + input >> din; + input >> dout; + // Line 7: time-stepping criteria + input >> timestepMax; // max no. of timesteps + input >> interval; // restart interval + input >> tol; // error tolerance + //............................................................. + + //....................................................................... + // Reading the domain information file + //....................................................................... + ifstream domain("Domain.in"); + domain >> nprocx; + domain >> nprocy; + domain >> nprocz; + domain >> Nx; + domain >> Ny; + domain >> Nz; + domain >> nspheres; + domain >> Lx; + domain >> Ly; + domain >> Lz; + //....................................................................... + + } + // ************************************************************** + // Broadcast simulation parameters from rank 0 to all other procs + MPI_Barrier(MPI_COMM_WORLD); + //................................................. + MPI_Bcast(&tau,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&alpha,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&beta,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&das,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&dbs,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&phi_s,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&wp_saturation,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&pBC,1,MPI_LOGICAL,0,MPI_COMM_WORLD); + MPI_Bcast(&Restart,1,MPI_LOGICAL,0,MPI_COMM_WORLD); + MPI_Bcast(&din,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&dout,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&Fx,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&Fy,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&Fz,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(×tepMax,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(&interval,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(&tol,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + // Computational domain + MPI_Bcast(&Nx,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(&Ny,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(&Nz,1,MPI_INT,0,MPI_COMM_WORLD); +// MPI_Bcast(&nBlocks,1,MPI_INT,0,MPI_COMM_WORLD); +// MPI_Bcast(&nthreads,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(&nprocx,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(&nprocy,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(&nprocz,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(&nspheres,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(&Lx,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&Ly,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + MPI_Bcast(&Lz,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + //................................................. + MPI_Barrier(MPI_COMM_WORLD); + + RESTART_INTERVAL=interval; + // ************************************************************** + // ************************************************************** + 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); + + // Set the density values inside the solid based on the input value phi_s + das = (phi_s+1.0)*0.5; + dbs = 1.0 - das; + + if (nprocs != nprocx*nprocy*nprocz){ + printf("nprocx = %i \n",nprocx); + printf("nprocy = %i \n",nprocy); + printf("nprocz = %i \n",nprocz); + INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); + } + + 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(MPI_COMM_WORLD); + + Nz += 2; + Nx = Ny = Nz; // Cubic domain + + int N = Nx*Ny*Nz; + int dist_mem_size = N*sizeof(double); + +// unsigned int nBlocks = 32; +// int nthreads = 128; + 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]; + char tmpstr[10]; + 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 sum_local; + double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs); + double porosity; +/* //....................................................................... + ifstream PM(LocalRankFilename,ios::binary); + for (k=0;k 0) sum++; + } + } + } + PM.close(); +// printf("File porosity = %f\n", double(sum)/N); +*/ //........................................................................... +// double *SignDist; +// SignDist = new double[N]; + DoubleArray SignDist(Nx,Ny,Nz); + //....................................................................... +#ifdef CBUB + // Initializes a constrained bubble test + double BubbleBot = 20.0; // How big to make the NWP bubble + double BubbleTop = 60.0; // How big to make the NWP bubble + double TubeRadius = 15.5; // Radius of the capillary tube + sum=0; + for (k=0;k 0.0){ + id[n] = 2; + } + // compute the porosity (actual interface location used) + if (SignDist.data[n] > 1.0){ + sum++; + } + } + } + } + //......................................................... + // If pressure boundary conditions are applied remove solid + if (pBC && kproc == 0){ + for (k=0; k<3; k++){ + for (j=0;j 0.0 && !pBC){ + phi_s = -phi_s; + das = (phi_s+1.0)*0.5; + dbs = 1.0 - das; + if (rank == 0) printf("Resetting phi_s = %f, das = %f, dbs = %f \n", phi_s, das, dbs); + FlipID(id,Nx*Ny*Nz); + } +#endif + + // Set up MPI communication structurese + if (rank==0) printf ("Setting up communication control structures \n"); + //...................................................................................... + // Get the actual D3Q19 communication counts (based on location of solid phase) + // Discrete velocity set symmetry implies the sendcount = recvcount + int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z; + int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ; + int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ; + sendCount_x = sendCount_y = sendCount_z = sendCount_X = sendCount_Y = sendCount_Z = 0; + sendCount_xy = sendCount_yz = sendCount_xz = sendCount_Xy = sendCount_Yz = sendCount_xZ = 0; + sendCount_xY = sendCount_yZ = sendCount_Xz = sendCount_XY = sendCount_YZ = sendCount_XZ = 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 DistValues(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) + // Set up kstart, kfinish so that the reservoirs are excluded from averaging + int kstart,kfinish; + kstart = 1; + kfinish = Nz-1; +// if (pBC && kproc==0) kstart = 4; +// if (pBC && kproc==nprocz-1) kfinish = Nz-4; + for (k=kstart; k CPU + //........................................................................... + DeviceBarrier(); + ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz); + CopyToHost(Phase.data,Phi,N*sizeof(double)); + CopyToHost(Press.data,Pressure,N*sizeof(double)); + CopyToHost(Vel_x.data,&Velocity[0],N*sizeof(double)); + CopyToHost(Vel_y.data,&Velocity[N],N*sizeof(double)); + CopyToHost(Vel_z.data,&Velocity[2*N],N*sizeof(double)); + MPI_Barrier(MPI_COMM_WORLD); + //........................................................................... + + 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(MPI_COMM_WORLD); + starttime = MPI_Wtime(); + //......................................... + + sendtag = recvtag = 5; + FILE *TIMELOG; + if (rank==0){ + TIMELOG= fopen("timelog.tcat","a+"); + if (fseek(TIMELOG,0,SEEK_SET) == fseek(TIMELOG,0,SEEK_CUR)){ + // If timelog is empty, write a short header to list the averages + fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); + fprintf(TIMELOG,"timestep dEs "); // Timestep, Change in Surface Energy + fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns efawns "); // Scalar averages + fprintf(TIMELOG,"vw[x, y, z] vn[x, y, z] vwn[x, y, z]"); // Velocity averages + fprintf(TIMELOG,"Gwn [xx, yy, zz, xy, xz, yz] "); // Orientation tensors + fprintf(TIMELOG,"Gws [xx, yy, zz, xy, xz, yz] "); + fprintf(TIMELOG,"Gns [xx, yy, zz, xy, xz, yz] "); + fprintf(TIMELOG,"trJwn trawn trRwn\n"); // trimmed curvature for wn surface + fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); + } + } + + err = 1.0; + sat_w_previous = 1.01; // slightly impossible value! + if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol); + //************ MAIN ITERATION LOOP ***************************************/ + while (timestep < timestepMax && err > tol){ + + //************************************************************************* + // Fused Color Gradient and Collision + //************************************************************************* + ColorCollideOpt( ID,f_even,f_odd,Phi,ColorGrad, + Velocity,Nx,Ny,Nz,rlxA,rlxB,alpha,beta,Fx,Fy,Fz); + //************************************************************************* + + //................................................................................... + PackDist(1,dvcSendList_x,0,sendCount_x,sendbuf_x,f_even,N); + PackDist(4,dvcSendList_x,sendCount_x,sendCount_x,sendbuf_x,f_even,N); + PackDist(5,dvcSendList_x,2*sendCount_x,sendCount_x,sendbuf_x,f_even,N); + PackDist(6,dvcSendList_x,3*sendCount_x,sendCount_x,sendbuf_x,f_even,N); + PackDist(7,dvcSendList_x,4*sendCount_x,sendCount_x,sendbuf_x,f_even,N); + //...Packing for X face(1,7,9,11,13)................................ + PackDist(0,dvcSendList_X,0,sendCount_X,sendbuf_X,f_odd,N); + PackDist(3,dvcSendList_X,sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + PackDist(4,dvcSendList_X,2*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + PackDist(5,dvcSendList_X,3*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + PackDist(6,dvcSendList_X,4*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + //...Packing for y face(4,8,9,16,18)................................. + PackDist(2,dvcSendList_y,0,sendCount_y,sendbuf_y,f_even,N); + PackDist(4,dvcSendList_y,sendCount_y,sendCount_y,sendbuf_y,f_even,N); + PackDist(4,dvcSendList_y,2*sendCount_y,sendCount_y,sendbuf_y,f_odd,N); + PackDist(8,dvcSendList_y,3*sendCount_y,sendCount_y,sendbuf_y,f_even,N); + PackDist(9,dvcSendList_y,4*sendCount_y,sendCount_y,sendbuf_y,f_even,N); + //...Packing for Y face(3,7,10,15,17)................................. + PackDist(1,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,f_odd,N); + PackDist(3,dvcSendList_Y,sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); + PackDist(5,dvcSendList_Y,2*sendCount_Y,sendCount_Y,sendbuf_Y,f_even,N); + PackDist(7,dvcSendList_Y,3*sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); + PackDist(8,dvcSendList_Y,4*sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); + //...Packing for z face(6,12,13,16,17)................................ + PackDist(3,dvcSendList_z,0,sendCount_z,sendbuf_z,f_even,N); + PackDist(6,dvcSendList_z,sendCount_z,sendCount_z,sendbuf_z,f_even,N); + PackDist(6,dvcSendList_z,2*sendCount_z,sendCount_z,sendbuf_z,f_odd,N); + PackDist(8,dvcSendList_z,3*sendCount_z,sendCount_z,sendbuf_z,f_even,N); + PackDist(8,dvcSendList_z,4*sendCount_z,sendCount_z,sendbuf_z,f_odd,N); + //...Packing for Z face(5,11,14,15,18)................................ + PackDist(2,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,f_odd,N); + PackDist(5,dvcSendList_Z,sendCount_Z,sendCount_Z,sendbuf_Z,f_odd,N); + PackDist(7,dvcSendList_Z,2*sendCount_Z,sendCount_Z,sendbuf_Z,f_even,N); + PackDist(7,dvcSendList_Z,3*sendCount_Z,sendCount_Z,sendbuf_Z,f_odd,N); + PackDist(9,dvcSendList_Z,4*sendCount_Z,sendCount_Z,sendbuf_Z,f_even,N); + //...Pack the xy edge (8)................................ + PackDist(4,dvcSendList_xy,0,sendCount_xy,sendbuf_xy,f_even,N); + //...Pack the Xy edge (9)................................ + PackDist(4,dvcSendList_Xy,0,sendCount_Xy,sendbuf_Xy,f_odd,N); + //...Pack the xY edge (10)................................ + PackDist(5,dvcSendList_xY,0,sendCount_xY,sendbuf_xY,f_even,N); + //...Pack the XY edge (7)................................ + PackDist(3,dvcSendList_XY,0,sendCount_XY,sendbuf_XY,f_odd,N); + //...Pack the xz edge (12)................................ + PackDist(6,dvcSendList_xz,0,sendCount_xz,sendbuf_xz,f_even,N); + //...Pack the xZ edge (14)................................ + PackDist(7,dvcSendList_xZ,0,sendCount_xZ,sendbuf_xZ,f_even,N); + //...Pack the Xz edge (13)................................ + PackDist(6,dvcSendList_Xz,0,sendCount_Xz,sendbuf_Xz,f_odd,N); + //...Pack the XZ edge (11)................................ + PackDist(5,dvcSendList_XZ,0,sendCount_XZ,sendbuf_XZ,f_odd,N); + //...Pack the xz edge (12)................................ + //...Pack the yz edge (16)................................ + PackDist(8,dvcSendList_yz,0,sendCount_yz,sendbuf_yz,f_even,N); + //...Pack the yZ edge (18)................................ + PackDist(9,dvcSendList_yZ,0,sendCount_yZ,sendbuf_yZ,f_even,N); + //...Pack the Yz edge (17)................................ + PackDist(8,dvcSendList_Yz,0,sendCount_Yz,sendbuf_Yz,f_odd,N); + //...Pack the YZ edge (15)................................ + PackDist(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,MPI_COMM_WORLD,&req1[0]); + MPI_Irecv(recvbuf_X, 5*recvCount_X,MPI_DOUBLE,rank_X,recvtag,MPI_COMM_WORLD,&req2[0]); + MPI_Isend(sendbuf_X, 5*sendCount_X,MPI_DOUBLE,rank_X,sendtag,MPI_COMM_WORLD,&req1[1]); + MPI_Irecv(recvbuf_x, 5*recvCount_x,MPI_DOUBLE,rank_x,recvtag,MPI_COMM_WORLD,&req2[1]); + MPI_Isend(sendbuf_y, 5*sendCount_y,MPI_DOUBLE,rank_y,sendtag,MPI_COMM_WORLD,&req1[2]); + MPI_Irecv(recvbuf_Y, 5*recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,MPI_COMM_WORLD,&req2[2]); + MPI_Isend(sendbuf_Y, 5*sendCount_Y,MPI_DOUBLE,rank_Y,sendtag,MPI_COMM_WORLD,&req1[3]); + MPI_Irecv(recvbuf_y, 5*recvCount_y,MPI_DOUBLE,rank_y,recvtag,MPI_COMM_WORLD,&req2[3]); + MPI_Isend(sendbuf_z, 5*sendCount_z,MPI_DOUBLE,rank_z,sendtag,MPI_COMM_WORLD,&req1[4]); + MPI_Irecv(recvbuf_Z, 5*recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,MPI_COMM_WORLD,&req2[4]); + MPI_Isend(sendbuf_Z, 5*sendCount_Z,MPI_DOUBLE,rank_Z,sendtag,MPI_COMM_WORLD,&req1[5]); + MPI_Irecv(recvbuf_z, 5*recvCount_z,MPI_DOUBLE,rank_z,recvtag,MPI_COMM_WORLD,&req2[5]); + MPI_Isend(sendbuf_xy, sendCount_xy,MPI_DOUBLE,rank_xy,sendtag,MPI_COMM_WORLD,&req1[6]); + MPI_Irecv(recvbuf_XY, recvCount_XY,MPI_DOUBLE,rank_XY,recvtag,MPI_COMM_WORLD,&req2[6]); + MPI_Isend(sendbuf_XY, sendCount_XY,MPI_DOUBLE,rank_XY,sendtag,MPI_COMM_WORLD,&req1[7]); + MPI_Irecv(recvbuf_xy, recvCount_xy,MPI_DOUBLE,rank_xy,recvtag,MPI_COMM_WORLD,&req2[7]); + MPI_Isend(sendbuf_Xy, sendCount_Xy,MPI_DOUBLE,rank_Xy,sendtag,MPI_COMM_WORLD,&req1[8]); + MPI_Irecv(recvbuf_xY, recvCount_xY,MPI_DOUBLE,rank_xY,recvtag,MPI_COMM_WORLD,&req2[8]); + MPI_Isend(sendbuf_xY, sendCount_xY,MPI_DOUBLE,rank_xY,sendtag,MPI_COMM_WORLD,&req1[9]); + MPI_Irecv(recvbuf_Xy, recvCount_Xy,MPI_DOUBLE,rank_Xy,recvtag,MPI_COMM_WORLD,&req2[9]); + MPI_Isend(sendbuf_xz, sendCount_xz,MPI_DOUBLE,rank_xz,sendtag,MPI_COMM_WORLD,&req1[10]); + MPI_Irecv(recvbuf_XZ, recvCount_XZ,MPI_DOUBLE,rank_XZ,recvtag,MPI_COMM_WORLD,&req2[10]); + MPI_Isend(sendbuf_XZ, sendCount_XZ,MPI_DOUBLE,rank_XZ,sendtag,MPI_COMM_WORLD,&req1[11]); + MPI_Irecv(recvbuf_xz, recvCount_xz,MPI_DOUBLE,rank_xz,recvtag,MPI_COMM_WORLD,&req2[11]); + MPI_Isend(sendbuf_Xz, sendCount_Xz,MPI_DOUBLE,rank_Xz,sendtag,MPI_COMM_WORLD,&req1[12]); + MPI_Irecv(recvbuf_xZ, recvCount_xZ,MPI_DOUBLE,rank_xZ,recvtag,MPI_COMM_WORLD,&req2[12]); + MPI_Isend(sendbuf_xZ, sendCount_xZ,MPI_DOUBLE,rank_xZ,sendtag,MPI_COMM_WORLD,&req1[13]); + MPI_Irecv(recvbuf_Xz, recvCount_Xz,MPI_DOUBLE,rank_Xz,recvtag,MPI_COMM_WORLD,&req2[13]); + MPI_Isend(sendbuf_yz, sendCount_yz,MPI_DOUBLE,rank_yz,sendtag,MPI_COMM_WORLD,&req1[14]); + MPI_Irecv(recvbuf_YZ, recvCount_YZ,MPI_DOUBLE,rank_YZ,recvtag,MPI_COMM_WORLD,&req2[14]); + MPI_Isend(sendbuf_YZ, sendCount_YZ,MPI_DOUBLE,rank_YZ,sendtag,MPI_COMM_WORLD,&req1[15]); + MPI_Irecv(recvbuf_yz, recvCount_yz,MPI_DOUBLE,rank_yz,recvtag,MPI_COMM_WORLD,&req2[15]); + MPI_Isend(sendbuf_Yz, sendCount_Yz,MPI_DOUBLE,rank_Yz,sendtag,MPI_COMM_WORLD,&req1[16]); + MPI_Irecv(recvbuf_yZ, recvCount_yZ,MPI_DOUBLE,rank_yZ,recvtag,MPI_COMM_WORLD,&req2[16]); + MPI_Isend(sendbuf_yZ, sendCount_yZ,MPI_DOUBLE,rank_yZ,sendtag,MPI_COMM_WORLD,&req1[17]); + MPI_Irecv(recvbuf_Yz, recvCount_Yz,MPI_DOUBLE,rank_Yz,recvtag,MPI_COMM_WORLD,&req2[17]); + //................................................................................... + + //************************************************************************* + // Carry out the density streaming step for mass transport + //************************************************************************* +// DensityStreamD3Q7(ID, Den, Copy, Phi, ColorGrad, Velocity, beta, Nx, Ny, Nz, pBC, S); + //************************************************************************* + MassColorCollideD3Q7(ID, A_even, A_odd, B_even, B_odd, Den, Phi, + ColorGrad, Velocity, beta, N, pBC); + + + //************************************************************************* + // Swap the distributions for momentum transport + //************************************************************************* + SwapD3Q19(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 ................................. + UnpackDist(0,-1,0,0,dvcRecvList_X,0,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + UnpackDist(3,-1,-1,0,dvcRecvList_X,recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + UnpackDist(4,-1,1,0,dvcRecvList_X,2*recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + UnpackDist(5,-1,0,-1,dvcRecvList_X,3*recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + UnpackDist(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.................................. + UnpackDist(1,1,0,0,dvcRecvList_x,0,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + UnpackDist(4,1,1,0,dvcRecvList_x,recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + UnpackDist(5,1,-1,0,dvcRecvList_x,2*recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + UnpackDist(6,1,0,1,dvcRecvList_x,3*recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + UnpackDist(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 ................................... + UnpackDist(1,0,-1,0,dvcRecvList_Y,0,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); + UnpackDist(3,-1,-1,0,dvcRecvList_Y,recvCount_Y,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); + UnpackDist(5,1,-1,0,dvcRecvList_Y,2*recvCount_Y,recvCount_Y,recvbuf_Y,f_even,Nx,Ny,Nz); + UnpackDist(7,0,-1,-1,dvcRecvList_Y,3*recvCount_Y,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); + UnpackDist(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 .................................. + UnpackDist(2,0,1,0,dvcRecvList_y,0,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); + UnpackDist(4,1,1,0,dvcRecvList_y,recvCount_y,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); + UnpackDist(4,-1,1,0,dvcRecvList_y,2*recvCount_y,recvCount_y,recvbuf_y,f_odd,Nx,Ny,Nz); + UnpackDist(8,0,1,1,dvcRecvList_y,3*recvCount_y,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); + UnpackDist(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).............................................. + UnpackDist(2,0,0,-1,dvcRecvList_Z,0,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); + UnpackDist(5,-1,0,-1,dvcRecvList_Z,recvCount_Z,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); + UnpackDist(7,1,0,-1,dvcRecvList_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,f_even,Nx,Ny,Nz); + UnpackDist(7,0,-1,-1,dvcRecvList_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); + UnpackDist(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).............................................. + UnpackDist(3,0,0,1,dvcRecvList_z,0,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); + UnpackDist(6,1,0,1,dvcRecvList_z,recvCount_z,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); + UnpackDist(6,-1,0,1,dvcRecvList_z,2*recvCount_z,recvCount_z,recvbuf_z,f_odd,Nx,Ny,Nz); + UnpackDist(8,0,1,1,dvcRecvList_z,3*recvCount_z,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); + UnpackDist(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)................................ + UnpackDist(3,-1,-1,0,dvcRecvList_XY,0,recvCount_XY,recvbuf_XY,f_odd,Nx,Ny,Nz); + //...Map recieve list for the Xy edge <<<9)................................ + UnpackDist(5,1,-1,0,dvcRecvList_xY,0,recvCount_xY,recvbuf_xY,f_even,Nx,Ny,Nz); + //...Map recieve list for the xY edge <<<10)................................ + UnpackDist(4,-1,1,0,dvcRecvList_Xy,0,recvCount_Xy,recvbuf_Xy,f_odd,Nx,Ny,Nz); + //...Map recieve list for the XY edge <<<7)................................ + UnpackDist(4,1,1,0,dvcRecvList_xy,0,recvCount_xy,recvbuf_xy,f_even,Nx,Ny,Nz); + //...Map recieve list for the xz edge <<<12)................................ + UnpackDist(5,-1,0,-1,dvcRecvList_XZ,0,recvCount_XZ,recvbuf_XZ,f_odd,Nx,Ny,Nz); + //...Map recieve list for the xZ edge <<<14)................................ + UnpackDist(6,-1,0,1,dvcRecvList_Xz,0,recvCount_Xz,recvbuf_Xz,f_odd,Nx,Ny,Nz); + //...Map recieve list for the Xz edge <<<13)................................ + UnpackDist(7,1,0,-1,dvcRecvList_xZ,0,recvCount_xZ,recvbuf_xZ,f_even,Nx,Ny,Nz); + //...Map recieve list for the XZ edge <<<11)................................ + UnpackDist(6,1,0,1,dvcRecvList_xz,0,recvCount_xz,recvbuf_xz,f_even,Nx,Ny,Nz); + //...Map recieve list for the yz edge <<<16)................................ + UnpackDist(7,0,-1,-1,dvcRecvList_YZ,0,recvCount_YZ,recvbuf_YZ,f_odd,Nx,Ny,Nz); + //...Map recieve list for the yZ edge <<<18)................................ + UnpackDist(8,0,-1,1,dvcRecvList_Yz,0,recvCount_Yz,recvbuf_Yz,f_odd,Nx,Ny,Nz); + //...Map recieve list for the Yz edge <<<17)................................ + UnpackDist(9,0,1,-1,dvcRecvList_yZ,0,recvCount_yZ,recvbuf_yZ,f_even,Nx,Ny,Nz); + //...Map recieve list for the YZ edge <<<15)................................ + UnpackDist(8,0,1,1,dvcRecvList_yz,0,recvCount_yz,recvbuf_yz,f_even,Nx,Ny,Nz); + //................................................................................... + + //................................................................................... + PackDist(1,dvcSendList_x,0,sendCount_x,sendbuf_x,A_even,N); + PackDist(1,dvcSendList_x,sendCount_x,sendCount_x,sendbuf_x,B_even,N); + //...Packing for X face(1,7,9,11,13)................................ + PackDist(0,dvcSendList_X,0,sendCount_X,sendbuf_X,A_odd,N); + PackDist(0,dvcSendList_X,sendCount_X,sendCount_X,sendbuf_X,B_odd,N); + //...Packing for y face(4,8,9,16,18)................................. + PackDist(2,dvcSendList_y,0,sendCount_y,sendbuf_y,A_even,N); + PackDist(2,dvcSendList_y,sendCount_y,sendCount_y,sendbuf_y,B_even,N); + //...Packing for Y face(3,7,10,15,17)................................. + PackDist(1,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,A_odd,N); + PackDist(1,dvcSendList_Y,sendCount_Y,sendCount_Y,sendbuf_Y,B_odd,N); + //...Packing for z face(6,12,13,16,17)................................ + PackDist(3,dvcSendList_z,0,sendCount_z,sendbuf_z,A_even,N); + PackDist(3,dvcSendList_z,sendCount_z,sendCount_z,sendbuf_z,B_even,N); + //...Packing for Z face(5,11,14,15,18)................................ + PackDist(2,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,A_odd,N); + PackDist(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,MPI_COMM_WORLD,&req1[0]); + MPI_Irecv(recvbuf_X, 2*recvCount_X,MPI_DOUBLE,rank_X,recvtag,MPI_COMM_WORLD,&req2[0]); + MPI_Isend(sendbuf_X, 2*sendCount_X,MPI_DOUBLE,rank_X,sendtag,MPI_COMM_WORLD,&req1[1]); + MPI_Irecv(recvbuf_x, 2*recvCount_x,MPI_DOUBLE,rank_x,recvtag,MPI_COMM_WORLD,&req2[1]); + MPI_Isend(sendbuf_y, 2*sendCount_y,MPI_DOUBLE,rank_y,sendtag,MPI_COMM_WORLD,&req1[2]); + MPI_Irecv(recvbuf_Y, 2*recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,MPI_COMM_WORLD,&req2[2]); + MPI_Isend(sendbuf_Y, 2*sendCount_Y,MPI_DOUBLE,rank_Y,sendtag,MPI_COMM_WORLD,&req1[3]); + MPI_Irecv(recvbuf_y, 2*recvCount_y,MPI_DOUBLE,rank_y,recvtag,MPI_COMM_WORLD,&req2[3]); + MPI_Isend(sendbuf_z, 2*sendCount_z,MPI_DOUBLE,rank_z,sendtag,MPI_COMM_WORLD,&req1[4]); + MPI_Irecv(recvbuf_Z, 2*recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,MPI_COMM_WORLD,&req2[4]); + MPI_Isend(sendbuf_Z, 2*sendCount_Z,MPI_DOUBLE,rank_Z,sendtag,MPI_COMM_WORLD,&req1[5]); + MPI_Irecv(recvbuf_z, 2*recvCount_z,MPI_DOUBLE,rank_z,recvtag,MPI_COMM_WORLD,&req2[5]); + //................................................................................... + + SwapD3Q7(ID, A_even, A_odd, Nx, Ny, Nz); + SwapD3Q7(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 ................................. + UnpackDist(0,-1,0,0,dvcRecvList_X,0,recvCount_X,recvbuf_X,A_odd,Nx,Ny,Nz); + UnpackDist(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.................................. + UnpackDist(1,1,0,0,dvcRecvList_x,0,recvCount_x,recvbuf_x,A_even,Nx,Ny,Nz); + UnpackDist(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 ................................... + UnpackDist(1,0,-1,0,dvcRecvList_Y,0,recvCount_Y,recvbuf_Y,A_odd,Nx,Ny,Nz); + UnpackDist(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 .................................. + UnpackDist(2,0,1,0,dvcRecvList_y,0,recvCount_y,recvbuf_y,A_even,Nx,Ny,Nz); + UnpackDist(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).............................................. + UnpackDist(2,0,0,-1,dvcRecvList_Z,0,recvCount_Z,recvbuf_Z,A_odd,Nx,Ny,Nz); + UnpackDist(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).............................................. + UnpackDist(3,0,0,1,dvcRecvList_z,0,recvCount_z,recvbuf_z,A_even,Nx,Ny,Nz); + UnpackDist(3,0,0,1,dvcRecvList_z,recvCount_z,recvCount_z,recvbuf_z,B_even,Nx,Ny,Nz); + //.................................................................................. + + //.................................................................................. + ComputeDensityD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz); + ComputeDensityD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz); + + //************************************************************************* + // Compute the phase indicator field + //************************************************************************* +// ComputePhi(ID, Phi, Copy, Den, N); + ComputePhi(ID, Phi, Den, N); + //************************************************************************* + + //................................................................................... + PackValues(dvcSendList_x, sendCount_x,sendbuf_x, Phi, N); + PackValues(dvcSendList_y, sendCount_y,sendbuf_y, Phi, N); + PackValues(dvcSendList_z, sendCount_z,sendbuf_z, Phi, N); + PackValues(dvcSendList_X, sendCount_X,sendbuf_X, Phi, N); + PackValues(dvcSendList_Y, sendCount_Y,sendbuf_Y, Phi, N); + PackValues(dvcSendList_Z, sendCount_Z,sendbuf_Z, Phi, N); + PackValues(dvcSendList_xy, sendCount_xy,sendbuf_xy, Phi, N); + PackValues(dvcSendList_xY, sendCount_xY,sendbuf_xY, Phi, N); + PackValues(dvcSendList_Xy, sendCount_Xy,sendbuf_Xy, Phi, N); + PackValues(dvcSendList_XY, sendCount_XY,sendbuf_XY, Phi, N); + PackValues(dvcSendList_xz, sendCount_xz,sendbuf_xz, Phi, N); + PackValues(dvcSendList_xZ, sendCount_xZ,sendbuf_xZ, Phi, N); + PackValues(dvcSendList_Xz, sendCount_Xz,sendbuf_Xz, Phi, N); + PackValues(dvcSendList_XZ, sendCount_XZ,sendbuf_XZ, Phi, N); + PackValues(dvcSendList_yz, sendCount_yz,sendbuf_yz, Phi, N); + PackValues(dvcSendList_yZ, sendCount_yZ,sendbuf_yZ, Phi, N); + PackValues(dvcSendList_Yz, sendCount_Yz,sendbuf_Yz, Phi, N); + PackValues(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,MPI_COMM_WORLD,&req1[0]); + MPI_Irecv(recvbuf_X, recvCount_X,MPI_DOUBLE,rank_X,recvtag,MPI_COMM_WORLD,&req2[0]); + MPI_Isend(sendbuf_X, sendCount_X,MPI_DOUBLE,rank_X,sendtag,MPI_COMM_WORLD,&req1[1]); + MPI_Irecv(recvbuf_x, recvCount_x,MPI_DOUBLE,rank_x,recvtag,MPI_COMM_WORLD,&req2[1]); + MPI_Isend(sendbuf_y, sendCount_y,MPI_DOUBLE,rank_y,sendtag,MPI_COMM_WORLD,&req1[2]); + MPI_Irecv(recvbuf_Y, recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,MPI_COMM_WORLD,&req2[2]); + MPI_Isend(sendbuf_Y, sendCount_Y,MPI_DOUBLE,rank_Y,sendtag,MPI_COMM_WORLD,&req1[3]); + MPI_Irecv(recvbuf_y, recvCount_y,MPI_DOUBLE,rank_y,recvtag,MPI_COMM_WORLD,&req2[3]); + MPI_Isend(sendbuf_z, sendCount_z,MPI_DOUBLE,rank_z,sendtag,MPI_COMM_WORLD,&req1[4]); + MPI_Irecv(recvbuf_Z, recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,MPI_COMM_WORLD,&req2[4]); + MPI_Isend(sendbuf_Z, sendCount_Z,MPI_DOUBLE,rank_Z,sendtag,MPI_COMM_WORLD,&req1[5]); + MPI_Irecv(recvbuf_z, recvCount_z,MPI_DOUBLE,rank_z,recvtag,MPI_COMM_WORLD,&req2[5]); + MPI_Isend(sendbuf_xy, sendCount_xy,MPI_DOUBLE,rank_xy,sendtag,MPI_COMM_WORLD,&req1[6]); + MPI_Irecv(recvbuf_XY, recvCount_XY,MPI_DOUBLE,rank_XY,recvtag,MPI_COMM_WORLD,&req2[6]); + MPI_Isend(sendbuf_XY, sendCount_XY,MPI_DOUBLE,rank_XY,sendtag,MPI_COMM_WORLD,&req1[7]); + MPI_Irecv(recvbuf_xy, recvCount_xy,MPI_DOUBLE,rank_xy,recvtag,MPI_COMM_WORLD,&req2[7]); + MPI_Isend(sendbuf_Xy, sendCount_Xy,MPI_DOUBLE,rank_Xy,sendtag,MPI_COMM_WORLD,&req1[8]); + MPI_Irecv(recvbuf_xY, recvCount_xY,MPI_DOUBLE,rank_xY,recvtag,MPI_COMM_WORLD,&req2[8]); + MPI_Isend(sendbuf_xY, sendCount_xY,MPI_DOUBLE,rank_xY,sendtag,MPI_COMM_WORLD,&req1[9]); + MPI_Irecv(recvbuf_Xy, recvCount_Xy,MPI_DOUBLE,rank_Xy,recvtag,MPI_COMM_WORLD,&req2[9]); + MPI_Isend(sendbuf_xz, sendCount_xz,MPI_DOUBLE,rank_xz,sendtag,MPI_COMM_WORLD,&req1[10]); + MPI_Irecv(recvbuf_XZ, recvCount_XZ,MPI_DOUBLE,rank_XZ,recvtag,MPI_COMM_WORLD,&req2[10]); + MPI_Isend(sendbuf_XZ, sendCount_XZ,MPI_DOUBLE,rank_XZ,sendtag,MPI_COMM_WORLD,&req1[11]); + MPI_Irecv(recvbuf_xz, recvCount_xz,MPI_DOUBLE,rank_xz,recvtag,MPI_COMM_WORLD,&req2[11]); + MPI_Isend(sendbuf_Xz, sendCount_Xz,MPI_DOUBLE,rank_Xz,sendtag,MPI_COMM_WORLD,&req1[12]); + MPI_Irecv(recvbuf_xZ, recvCount_xZ,MPI_DOUBLE,rank_xZ,recvtag,MPI_COMM_WORLD,&req2[12]); + MPI_Isend(sendbuf_xZ, sendCount_xZ,MPI_DOUBLE,rank_xZ,sendtag,MPI_COMM_WORLD,&req1[13]); + MPI_Irecv(recvbuf_Xz, recvCount_Xz,MPI_DOUBLE,rank_Xz,recvtag,MPI_COMM_WORLD,&req2[13]); + MPI_Isend(sendbuf_yz, sendCount_yz,MPI_DOUBLE,rank_yz,sendtag,MPI_COMM_WORLD,&req1[14]); + MPI_Irecv(recvbuf_YZ, recvCount_YZ,MPI_DOUBLE,rank_YZ,recvtag,MPI_COMM_WORLD,&req2[14]); + MPI_Isend(sendbuf_YZ, sendCount_YZ,MPI_DOUBLE,rank_YZ,sendtag,MPI_COMM_WORLD,&req1[15]); + MPI_Irecv(recvbuf_yz, recvCount_yz,MPI_DOUBLE,rank_yz,recvtag,MPI_COMM_WORLD,&req2[15]); + MPI_Isend(sendbuf_Yz, sendCount_Yz,MPI_DOUBLE,rank_Yz,sendtag,MPI_COMM_WORLD,&req1[16]); + MPI_Irecv(recvbuf_yZ, recvCount_yZ,MPI_DOUBLE,rank_yZ,recvtag,MPI_COMM_WORLD,&req2[16]); + MPI_Isend(sendbuf_yZ, sendCount_yZ,MPI_DOUBLE,rank_yZ,sendtag,MPI_COMM_WORLD,&req1[17]); + MPI_Irecv(recvbuf_Yz, recvCount_Yz,MPI_DOUBLE,rank_Yz,recvtag,MPI_COMM_WORLD,&req2[17]); + //................................................................................... + //................................................................................... + // Wait for completion of Indicator Field communication + //................................................................................... + MPI_Waitall(18,req1,stat1); + MPI_Waitall(18,req2,stat2); + DeviceBarrier(); + //................................................................................... + //................................................................................... +/* UnpackValues(faceGrid, packThreads, dvcSendList_x, sendCount_x,sendbuf_x, Phi, N); + UnpackValues(faceGrid, packThreads, dvcSendList_y, sendCount_y,sendbuf_y, Phi, N); + UnpackValues(faceGrid, packThreads, dvcSendList_z, sendCount_z,sendbuf_z, Phi, N); + UnpackValues(faceGrid, packThreads, dvcSendList_X, sendCount_X,sendbuf_X, Phi, N); + UnpackValues(faceGrid, packThreads, dvcSendList_Y, sendCount_Y,sendbuf_Y, Phi, N); + UnpackValues(faceGrid, packThreads, dvcSendList_Z, sendCount_Z,sendbuf_Z, Phi, N); +*/ + UnpackValues(dvcRecvList_x, recvCount_x,recvbuf_x, Phi, N); + UnpackValues(dvcRecvList_y, recvCount_y,recvbuf_y, Phi, N); + UnpackValues(dvcRecvList_z, recvCount_z,recvbuf_z, Phi, N); + UnpackValues(dvcRecvList_X, recvCount_X,recvbuf_X, Phi, N); + UnpackValues(dvcRecvList_Y, recvCount_Y,recvbuf_Y, Phi, N); + UnpackValues(dvcRecvList_Z, recvCount_Z,recvbuf_Z, Phi, N); + UnpackValues(dvcRecvList_xy, recvCount_xy,recvbuf_xy, Phi, N); + UnpackValues(dvcRecvList_xY, recvCount_xY,recvbuf_xY, Phi, N); + UnpackValues(dvcRecvList_Xy, recvCount_Xy,recvbuf_Xy, Phi, N); + UnpackValues(dvcRecvList_XY, recvCount_XY,recvbuf_XY, Phi, N); + UnpackValues(dvcRecvList_xz, recvCount_xz,recvbuf_xz, Phi, N); + UnpackValues(dvcRecvList_xZ, recvCount_xZ,recvbuf_xZ, Phi, N); + UnpackValues(dvcRecvList_Xz, recvCount_Xz,recvbuf_Xz, Phi, N); + UnpackValues(dvcRecvList_XZ, recvCount_XZ,recvbuf_XZ, Phi, N); + UnpackValues(dvcRecvList_yz, recvCount_yz,recvbuf_yz, Phi, N); + UnpackValues(dvcRecvList_yZ, recvCount_yZ,recvbuf_yZ, Phi, N); + UnpackValues(dvcRecvList_Yz, recvCount_Yz,recvbuf_Yz, Phi, N); + UnpackValues(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, Phi, N); + //................................................................................... + + + if (pBC && kproc == 0) { + PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz); + ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + } + + if (pBC && kproc == nprocz-1){ + PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); + ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + } + + //................................................................................... + + MPI_Barrier(MPI_COMM_WORLD); + + // Timestep completed! + timestep++; + //................................................................... + if (timestep%1000 == 995){ + //........................................................................... + // Copy the phase indicator field for the earlier timestep + DeviceBarrier(); + CopyToHost(Phase_tplus.data,Phi,N*sizeof(double)); + //........................................................................... + } + if (timestep%1000 == 0){ + //........................................................................... + // Copy the data for for the analysis timestep + //........................................................................... + // Copy the phase from the GPU -> CPU + //........................................................................... + DeviceBarrier(); + ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz); + CopyToHost(Phase.data,Phi,N*sizeof(double)); + CopyToHost(Press.data,Pressure,N*sizeof(double)); + CopyToHost(Vel_x.data,&Velocity[0],N*sizeof(double)); + CopyToHost(Vel_y.data,&Velocity[N],N*sizeof(double)); + CopyToHost(Vel_z.data,&Velocity[2*N],N*sizeof(double)); + MPI_Barrier(MPI_COMM_WORLD); + } + if (timestep%1000 == 5){ + //........................................................................... + // Copy the phase indicator field for the later timestep + DeviceBarrier(); + CopyToHost(Phase_tminus.data,Phi,N*sizeof(double)); + //........................................................................... + // Calculate the time derivative of the phase indicator field + for (n=0; n 0.0 ){ + + // 1-D index for this cube corner + n = i + j*Nx + k*Nx*Ny; + // Compute the non-wetting phase volume contribution + if ( Phase(i,j,k) > 0.0 ) + nwp_volume += 1.0; + + // volume averages over the non-wetting phase + if ( Phase(i,j,k) > 0.999 ){ + // volume the excludes the interfacial region + vol_n += 1.0; + // pressure + pan += Press(i,j,k); + // velocity + van(0) += Vel_x(i,j,k); + van(1) += Vel_y(i,j,k); + van(2) += Vel_z(i,j,k); + } + + // volume averages over the wetting phase + if ( Phase(i,j,k) < -0.999 ){ + // volume the excludes the interfacial region + vol_w += 1.0; + // pressure + paw += Press(i,j,k); + // velocity + vaw(0) += Vel_x(i,j,k); + vaw(1) += Vel_y(i,j,k); + vaw(2) += Vel_z(i,j,k); + } + } + } + } + } + for (c=0;c 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 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 averages over the non-wetting phase + if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0.99 ){ + // volume the excludes the interfacial region + vol_n += 0.125; + // pressure + pan += 0.125*Press.data[n]; + // velocity + van(0) += 0.125*Vel_x.data[n]; + van(1) += 0.125*Vel_y.data[n]; + van(2) += 0.125*Vel_z.data[n]; + } + + // volume averages over the wetting phase + if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < -0.99 ){ + // volume the excludes the interfacial region + vol_w += 0.125; + // pressure + paw += 0.125*Press.data[n]; + // velocity + vaw(0) += 0.125*Vel_x.data[n]; + vaw(1) += 0.125*Vel_y.data[n]; + vaw(2) += 0.125*Vel_z.data[n]; + } + } + } + +*/ //........................................................................... + // Construct the interfaces and common curve + pmmc_ConstructLocalCube(SignDist, 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); + + // Integrate the contact angle + efawns += pmmc_CubeContactAngle(CubeValues,Values,Phase_x,Phase_y,Phase_z,SignDist_x,SignDist_y,SignDist_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); + Kwn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); + + // Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels) + pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SignDist,nw_pts,nw_tris,Values,DistValues, + i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn); + + pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SignDist,nw_pts,nw_tris,Values,DistValues, + i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn); + + // Compute the normal speed of the interface + 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); + + 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(MPI_COMM_WORLD); + MPI_Allreduce(&nwp_volume,&nwp_volume_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&awn,&awn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&ans,&ans_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&aws,&aws_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&lwns,&lwns_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&As,&As_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&Jwn,&Jwn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&Kwn,&Kwn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + // Phase averages + MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&paw,&paw_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&pan,&pan_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&vaw(0),&vaw_global(0),3,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&van(0),&van_global(0),3,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&vawn(0),&vawn_global(0),3,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&Gwn(0),&Gwn_global(0),6,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&Gns(0),&Gns_global(0),6,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&Gws(0),&Gws_global(0),6,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&trawn,&trawn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&trJwn,&trJwn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Allreduce(&trRwn,&trRwn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + MPI_Barrier(MPI_COMM_WORLD); + //......................................................................... + // Compute the change in the total surface energy based on the defined interval + // See McClure, Prins and Miller (2014) + //......................................................................... + 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; + Kwn_global /= awn_global; + efawns_global /= lwns_global; + + if (trawn_global > 0.0) trJwn_global /= trawn_global; + if (trawn_global > 0.0) trRwn_global /= trawn_global; + trRwn_global = 2.0*fabs(trRwn_global); + trJwn_global = fabs(trJwn_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 (dimensionless per unit volume) + awn_global = awn_global*iVol_global*D; + ans_global = ans_global*iVol_global*D; + aws_global = aws_global*iVol_global*D; + dEs = dEs*iVol_global*D; + lwns_global = lwns_global*iVol_global*D*D; + + //......................................................................... + if (rank==0){ +/* printf("%i %.5g ",timestep-5,dEs); // change in surface energy + printf("%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure + printf("%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas + printf("%.5g %5g ",Jwn_global, Kwn_global); // curvature of wn interface + printf("%.5g ",lwns_global); // common curve length + printf("%.5g ",efawns_global); // average contact angle + printf("%.5g %.5g %.5g ",vaw_global(0),vaw_global(1),vaw_global(2)); // average velocity of w phase + printf("%.5g %.5g %.5g ",van_global(0),van_global(1),van_global(2)); // average velocity of n phase + printf("%.5g %.5g %.5g ",vawn_global(0),vawn_global(1),vawn_global(2)); // velocity of wn interface + printf("%.5g %.5g %.5g %.5g %.5g %.5g ", + Gwn_global(0),Gwn_global(1),Gwn_global(2),Gwn_global(3),Gwn_global(4),Gwn_global(5)); // orientation of wn interface + printf("%.5g %.5g %.5g %.5g %.5g %.5g ", + Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5)); // orientation of ns interface + printf("%.5g %.5g %.5g %.5g %.5g %.5g ", + Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface + printf("%.5g %5g \n",trawn_global, trJwn_global); // Trimmed curvature + +*/ + fprintf(TIMELOG,"%i %.5g ",timestep-5,dEs); // change in surface energy + fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure + fprintf(TIMELOG,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas + fprintf(TIMELOG,"%.5g %5g ",Jwn_global, Kwn_global); // curvature of wn interface + fprintf(TIMELOG,"%.5g ",lwns_global); // common curve length + fprintf(TIMELOG,"%.5g ",efawns_global); // average contact angle + fprintf(TIMELOG,"%.5g %.5g %.5g ",vaw_global(0),vaw_global(1),vaw_global(2)); // average velocity of w phase + fprintf(TIMELOG,"%.5g %.5g %.5g ",van_global(0),van_global(1),van_global(2)); // average velocity of n phase + fprintf(TIMELOG,"%.5g %.5g %.5g ",vawn_global(0),vawn_global(1),vawn_global(2)); // velocity of wn interface + fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", + Gwn_global(0),Gwn_global(1),Gwn_global(2),Gwn_global(3),Gwn_global(4),Gwn_global(5)); // orientation of wn interface + fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", + Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5)); // orientation of ns interface + fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", + Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface + fprintf(TIMELOG,"%.5g %5g %5g\n",trawn_global, trJwn_global, trRwn_global); // Trimmed curvature + fflush(TIMELOG); + } + } + + if (timestep%RESTART_INTERVAL == 0){ + if (pBC){ + err = fabs(sat_w - sat_w_previous); + sat_w_previous = sat_w; + if (rank==0) printf("Timestep %i: change in saturation since last checkpoint is %f \n", timestep, err); + } + else{ + // Not clear yet + } + // Copy the data to the CPU + CopyToHost(cDistEven,f_even,10*N*sizeof(double)); + CopyToHost(cDistOdd,f_odd,9*N*sizeof(double)); + CopyToHost(cDen,Den,2*N*sizeof(double)); + // Read in the restart file to CPU buffers + WriteCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N); + +#ifdef WRITE_SURFACES + + sprintf(tmpstr,"vis%03d",logcount); + if (rank==0){ + mkdir(tmpstr,0777); + } + MPI_Barrier(MPI_COMM_WORLD); + + FILE *WN_TRIS; + sprintf(LocalRankFilename,"%s/%s%s",tmpstr,"wn-tris.",LocalRankString); + WN_TRIS = fopen(LocalRankFilename,"wb"); + + FILE *NS_TRIS; + sprintf(LocalRankFilename,"%s/%s%s",tmpstr,"ns-tris.",LocalRankString); + NS_TRIS = fopen(LocalRankFilename,"wb"); + + FILE *WS_TRIS; + sprintf(LocalRankFilename,"%s/%s%s",tmpstr,"ws-tris.",LocalRankString); + WS_TRIS = fopen(LocalRankFilename,"wb"); + + FILE *WNS_PTS; + sprintf(LocalRankFilename,"%s/%s%s",tmpstr,"wns-crv.",LocalRankString); + WNS_PTS = fopen(LocalRankFilename,"wb"); + + for (c=0;c Date: Wed, 25 Jun 2014 11:04:44 -0400 Subject: [PATCH 04/43] Basic changes to blob code --- tests/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6c84e729..7bea272a 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,5 +1,6 @@ # Copy files for the tests INSTALL_LBPM_EXE( lb2_Color_wia_mpi ) +INSTALL_LBPM_EXE( lb2_Blob_wia_mpi ) INSTALL_LBPM_EXE( TestBubble ) CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_DIR}/cylindertest COPYONLY ) From 3927c83c1d4aa72f0217146fd9dd51c8e253abff Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 12:45:12 -0400 Subject: [PATCH 05/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 12bfed69..6895abf5 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -1254,6 +1254,7 @@ int main(int argc, char **argv) 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); @@ -1321,7 +1322,9 @@ int main(int argc, char **argv) } } } - + IntArray BlobList(3,ncubes); // store indices for blobs (cubes) + IntArray temp(3,ncubes); // temporary storage array + IntArray BlobSizes(MAX_LOCAL_BLOB_COUNT); // number of nodes in each blob int logcount = 0; // number of surface write-outs @@ -2305,6 +2308,26 @@ int main(int argc, char **argv) Jwn = Kwn = efawns = 0.0; trJwn = trawn = trRwn = 0.0; + nblobs=0; + double vF,vS; + vF = vS = 0.0; + for (k=0;k 0.0 ){ + if ( SignDist(i,j,k) > 0.0 ){ + // node i,j,k is in the porespace + BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,ncubes,LocalBlobID,Phase,SignDist,vF,vS,i,j,k,temp); + nblobs++; + } + } + } + + } + } + } + /// Compute volume averages for (k=kstart; k Date: Wed, 25 Jun 2014 12:58:48 -0400 Subject: [PATCH 06/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 6895abf5..abfd08f3 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -1322,6 +1322,7 @@ int main(int argc, char **argv) } } } + int nblobs = 0; IntArray BlobList(3,ncubes); // store indices for blobs (cubes) IntArray temp(3,ncubes); // temporary storage array IntArray BlobSizes(MAX_LOCAL_BLOB_COUNT); // number of nodes in each blob @@ -2308,13 +2309,21 @@ int main(int argc, char **argv) Jwn = Kwn = efawns = 0.0; trJwn = trawn = trRwn = 0.0; + for (k=0;k 0.0 ){ if ( SignDist(i,j,k) > 0.0 ){ // node i,j,k is in the porespace From 5ce9a9e30ef0751a7360c5b4fc8e541ef7123039 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 13:00:32 -0400 Subject: [PATCH 07/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index abfd08f3..747587ab 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -2817,6 +2817,12 @@ int main(int argc, char **argv) fwrite(Press.data,8,N,PRESS); fclose(PRESS); + sprintf(LocalRankFilename,"%s%s","BlobID.",LocalRankString); + FILE *BLOB; + BLOB = fopen(LocalRankFilename,"wb"); + fwrite(LocalBlobID.data,4,N,BLOB); + fclose(BLOB); + /* sprintf(LocalRankFilename,"%s%s","dPdt.",LocalRankString); FILE *SPEED; SPEED = fopen(LocalRankFilename,"wb"); From 86c74310378da48ec174a8aa54689ea2fb6caf82 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 14:44:54 -0400 Subject: [PATCH 08/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 747587ab..b8dec517 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -96,7 +96,6 @@ inline int ComputeLocalBlob(IntArray blobs, int &nblobs, int &ncubes, IntArray i printf("Error in blob search algorithm!"); } } - } if ( nrecent == 0){ status = 0; @@ -2320,13 +2319,14 @@ int main(int argc, char **argv) nblobs=0; double vF,vS; vF = vS = 0.0; - for (k=0;k 0.0 ){ if ( SignDist(i,j,k) > 0.0 ){ // node i,j,k is in the porespace + printf(); BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,ncubes,LocalBlobID,Phase,SignDist,vF,vS,i,j,k,temp); nblobs++; } From 820c773ffadbc07b6e257106e0d1eb2f83609264 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 14:45:55 -0400 Subject: [PATCH 09/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index b8dec517..93ba31ea 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -2326,7 +2326,7 @@ int main(int argc, char **argv) if ( Phase(i,j,k) > 0.0 ){ if ( SignDist(i,j,k) > 0.0 ){ // node i,j,k is in the porespace - printf(); + printf("Compute blob %i, \n",nblobs); BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,ncubes,LocalBlobID,Phase,SignDist,vF,vS,i,j,k,temp); nblobs++; } From a459367d130af34d64fadbde1e3b6fb9015719ab Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 14:48:33 -0400 Subject: [PATCH 10/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 93ba31ea..a6d7fc5e 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -2326,7 +2326,7 @@ int main(int argc, char **argv) if ( Phase(i,j,k) > 0.0 ){ if ( SignDist(i,j,k) > 0.0 ){ // node i,j,k is in the porespace - printf("Compute blob %i, \n",nblobs); + // printf("Compute blob %i, \n",nblobs); BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,ncubes,LocalBlobID,Phase,SignDist,vF,vS,i,j,k,temp); nblobs++; } @@ -2336,7 +2336,9 @@ int main(int argc, char **argv) } } } - + printf("Number of blobs %i, Rank = %i \n",nblobs,rank); + MPI_Barrier(MPI_COMM_WORLD); + /// Compute volume averages for (k=kstart; k Date: Wed, 25 Jun 2014 14:51:02 -0400 Subject: [PATCH 11/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index a6d7fc5e..1234cbc3 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -2316,7 +2316,7 @@ int main(int argc, char **argv) } } - nblobs=0; + nblobs=0, nc=0; double vF,vS; vF = vS = 0.0; for (k=0;k 0.0 ){ // node i,j,k is in the porespace // printf("Compute blob %i, \n",nblobs); - BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,ncubes,LocalBlobID,Phase,SignDist,vF,vS,i,j,k,temp); + BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,nc,LocalBlobID,Phase,SignDist,vF,vS,i,j,k,temp); nblobs++; } } From 27f52898adf14be5093f54799b347cd4c2e0c1bb Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 15:56:23 -0400 Subject: [PATCH 12/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 142 ++++++++++++++++++++++++++++++++++++- 1 file changed, 140 insertions(+), 2 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 1234cbc3..debb927d 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -22,6 +22,142 @@ using namespace std; +struct BlobInfo{ + + BlobInfo(int Nx,int Ny,int Nz){ + IntArray.New(Nx,Ny,Nz); + } + IntArray ID; + + +}; + +inline void PackBlobData(int *list, int count, int *sendbuf, int *data){ + // 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 0.0 ){ // node i,j,k is in the porespace // printf("Compute blob %i, \n",nblobs); - BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,nc,LocalBlobID,Phase,SignDist,vF,vS,i,j,k,temp); + BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,nc,LocalBlobs.ID,Phase,SignDist,vF,vS,i,j,k,temp); nblobs++; } } From 996ad21bda1bdc995be43d52e0189517873d8652 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 15:56:50 -0400 Subject: [PATCH 13/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index debb927d..ef0a49ba 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -2960,7 +2960,7 @@ int main(int argc, char **argv) sprintf(LocalRankFilename,"%s%s","BlobID.",LocalRankString); FILE *BLOB; BLOB = fopen(LocalRankFilename,"wb"); - fwrite(LocalBlobID.data,4,N,BLOB); + fwrite(LocalBlobs.ID.data,4,N,BLOB); fclose(BLOB); /* sprintf(LocalRankFilename,"%s%s","dPdt.",LocalRankString); From 00cd19092d56d05d684f58b5b6fd7f581d6c2023 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 15:57:28 -0400 Subject: [PATCH 14/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index ef0a49ba..dc76f9a1 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -24,7 +24,7 @@ using namespace std; struct BlobInfo{ - BlobInfo(int Nx,int Ny,int Nz){ + void BlobInfo(int Nx, int Ny,int Nz){ IntArray.New(Nx,Ny,Nz); } IntArray ID; From 38031e6ceb0bbe27fcb329d4cfcc263477f0b4fa Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 15:58:05 -0400 Subject: [PATCH 15/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index dc76f9a1..95d69a02 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -25,7 +25,7 @@ using namespace std; struct BlobInfo{ void BlobInfo(int Nx, int Ny,int Nz){ - IntArray.New(Nx,Ny,Nz); + ID.New(Nx,Ny,Nz); } IntArray ID; From 32d0bc3a126c19ee2a762b8d0d77b57b08654eb8 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 15:58:20 -0400 Subject: [PATCH 16/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 95d69a02..1ae526f7 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -24,7 +24,7 @@ using namespace std; struct BlobInfo{ - void BlobInfo(int Nx, int Ny,int Nz){ + BlobInfo(int Nx, int Ny,int Nz){ ID.New(Nx,Ny,Nz); } IntArray ID; From d6b6772303fdd1fdbdac45d22083cd926d0cd41e Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 16:22:24 -0400 Subject: [PATCH 17/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 45 +++++++++++++++++++++++++++++++++----- 1 file changed, 40 insertions(+), 5 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 1ae526f7..bb931251 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -22,6 +22,35 @@ using namespace std; +struct Domain{ + + Domain(int nx, int ny, int nz, int rnk, int npx, int npy, int npz, + double lx, double ly, double lz){ + Nx = nx; Ny = ny; Nz = nz; + Lx = lx, Ly = ly, Lz = lz; + rank = rnk; + nprocx=npx; nprocy=npy; nprocz=npz; + N = Nx*Ny*Nz; + ID = new char [N]; + Blobs.New(Nx,Ny,Nz); + } + + // Basic domain information + int Nx,Ny,Nz; + int rank; + int iproc,jproc,kproc; + int nprocx,nprocy,nprocz; + double Lx,Ly,Lz; + char *ID; + // Neighbors and MPI information + int rank; + int rank_x,rank_X; + + // Blob information + IntArray Blobs; + +}; + struct BlobInfo{ BlobInfo(int Nx, int Ny,int Nz){ @@ -29,6 +58,11 @@ struct BlobInfo{ } IntArray ID; + void GenerateTable(); + void CommunicateHalo(); + + void Pack(int *list, int count, int *sendbuf, int *data); + void UnpackBlobData(int *list, int count, int *recvbuf, int *data); }; @@ -544,6 +578,8 @@ int main(int argc, char **argv) printf("********************************************************\n"); } + Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz); + 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, @@ -1330,7 +1366,6 @@ int main(int argc, char **argv) // Vel = new double[3*N]; // fluid velocity // Press = new double[N]; // fluid pressure - BlobInfo LocalBlobs(Nx,Ny,Nz); IntArray LocalBlobID(Nx,Ny,Nz); DoubleArray Press(Nx,Ny,Nz); @@ -2449,7 +2484,7 @@ int main(int argc, char **argv) for (k=0;k 0.0 ){ if ( SignDist(i,j,k) > 0.0 ){ // node i,j,k is in the porespace // printf("Compute blob %i, \n",nblobs); - BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,nc,LocalBlobs.ID,Phase,SignDist,vF,vS,i,j,k,temp); + BlobSizes(nblobs) = ComputeBlob(BlobList,nblobs,nc,Dm.Blobs,Phase,SignDist,vF,vS,i,j,k,temp); nblobs++; } } @@ -2960,7 +2995,7 @@ int main(int argc, char **argv) sprintf(LocalRankFilename,"%s%s","BlobID.",LocalRankString); FILE *BLOB; BLOB = fopen(LocalRankFilename,"wb"); - fwrite(LocalBlobs.ID.data,4,N,BLOB); + fwrite(Dm.Blobs.data,4,N,BLOB); fclose(BLOB); /* sprintf(LocalRankFilename,"%s%s","dPdt.",LocalRankString); From d40e2eab5fc3e3ccb5250d445411658246b69479 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 16:23:44 -0400 Subject: [PATCH 18/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index bb931251..267998fd 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -36,7 +36,7 @@ struct Domain{ } // Basic domain information - int Nx,Ny,Nz; + int Nx,Ny,Nz,N; int rank; int iproc,jproc,kproc; int nprocx,nprocy,nprocz; @@ -2495,7 +2495,7 @@ int main(int argc, char **argv) for (k=0;k 0.0 ){ if ( SignDist(i,j,k) > 0.0 ){ // node i,j,k is in the porespace From 059218db600a56f0d2b69f04bd57e20f8c4ba59c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 16:24:27 -0400 Subject: [PATCH 19/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 267998fd..f6e92756 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -37,12 +37,10 @@ struct Domain{ // Basic domain information int Nx,Ny,Nz,N; - int rank; int iproc,jproc,kproc; int nprocx,nprocy,nprocz; double Lx,Ly,Lz; char *ID; - // Neighbors and MPI information int rank; int rank_x,rank_X; From 33fead0904746f1e1414d7f312941a6804cd4afc Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 16:31:13 -0400 Subject: [PATCH 20/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 53 +++++++++++++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 3 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index f6e92756..51e7bb73 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -34,19 +34,66 @@ struct Domain{ ID = new char [N]; Blobs.New(Nx,Ny,Nz); } - + + void InitializeRanks() + { + // map the rank to the block index + iproc = rank%nprocx; + jproc = (rank/nprocx)%nprocy; + kproc = rank/(nprocx*nprocy); + + // set up the neighbor ranks + int i = iproc; + int j = jproc; + int k = kproc; + rank_X = getRankForBlock(nprocx,nprocy,nprocz,i+1,j,k); + rank_x = getRankForBlock(nprocx,nprocy,nprocz,i-1,j,k); + rank_Y = getRankForBlock(nprocx,nprocy,nprocz,i,j+1,k); + rank_y = getRankForBlock(nprocx,nprocy,nprocz,i,j-1,k); + rank_Z = getRankForBlock(nprocx,nprocy,nprocz,i,j,k+1); + rank_z = getRankForBlock(nprocx,nprocy,nprocz,i,j,k-1); + rank_XY = getRankForBlock(nprocx,nprocy,nprocz,i+1,j+1,k); + rank_xy = getRankForBlock(nprocx,nprocy,nprocz,i-1,j-1,k); + rank_Xy = getRankForBlock(nprocx,nprocy,nprocz,i+1,j-1,k); + rank_xY = getRankForBlock(nprocx,nprocy,nprocz,i-1,j+1,k); + rank_XZ = getRankForBlock(nprocx,nprocy,nprocz,i+1,j,k+1); + rank_xz = getRankForBlock(nprocx,nprocy,nprocz,i-1,j,k-1); + rank_Xz = getRankForBlock(nprocx,nprocy,nprocz,i+1,j,k-1); + rank_xZ = getRankForBlock(nprocx,nprocy,nprocz,i-1,j,k+1); + rank_YZ = getRankForBlock(nprocx,nprocy,nprocz,i,j+1,k+1); + rank_yz = getRankForBlock(nprocx,nprocy,nprocz,i,j-1,k-1); + rank_Yz = getRankForBlock(nprocx,nprocy,nprocz,i,j+1,k-1); + rank_yZ = getRankForBlock(nprocx,nprocy,nprocz,i,j-1,k+1); + } // Basic domain information int Nx,Ny,Nz,N; int iproc,jproc,kproc; int nprocx,nprocy,nprocz; double Lx,Ly,Lz; - char *ID; int rank; - int rank_x,rank_X; + //********************************** + // MPI ranks for all 18 neighbors + //********************************** + int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z; + int rank_xy,rank_XY,rank_xY,rank_Xy; + int rank_xz,rank_XZ,rank_xZ,rank_Xz; + int rank_yz,rank_YZ,rank_yZ,rank_Yz; + //********************************** + // Solid indicator function + char *ID; // Blob information IntArray Blobs; +private: + int getRankForBlock( int i, int j, int k ) + { + int i2 = (i+nprocx)%nprocx; + int j2 = (j+nprocy)%nprocy; + int k2 = (k+nprocz)%nprocz; + return i2 + j2*nprocx + k2*nprocx*nprocy; + } + }; struct BlobInfo{ From da9b84837dbb67aaff585d9a3afd572c1fbff95c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 16:32:11 -0400 Subject: [PATCH 21/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 51e7bb73..f9afe992 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -46,24 +46,24 @@ struct Domain{ int i = iproc; int j = jproc; int k = kproc; - rank_X = getRankForBlock(nprocx,nprocy,nprocz,i+1,j,k); - rank_x = getRankForBlock(nprocx,nprocy,nprocz,i-1,j,k); - rank_Y = getRankForBlock(nprocx,nprocy,nprocz,i,j+1,k); - rank_y = getRankForBlock(nprocx,nprocy,nprocz,i,j-1,k); - rank_Z = getRankForBlock(nprocx,nprocy,nprocz,i,j,k+1); - rank_z = getRankForBlock(nprocx,nprocy,nprocz,i,j,k-1); - rank_XY = getRankForBlock(nprocx,nprocy,nprocz,i+1,j+1,k); - rank_xy = getRankForBlock(nprocx,nprocy,nprocz,i-1,j-1,k); - rank_Xy = getRankForBlock(nprocx,nprocy,nprocz,i+1,j-1,k); - rank_xY = getRankForBlock(nprocx,nprocy,nprocz,i-1,j+1,k); - rank_XZ = getRankForBlock(nprocx,nprocy,nprocz,i+1,j,k+1); - rank_xz = getRankForBlock(nprocx,nprocy,nprocz,i-1,j,k-1); - rank_Xz = getRankForBlock(nprocx,nprocy,nprocz,i+1,j,k-1); - rank_xZ = getRankForBlock(nprocx,nprocy,nprocz,i-1,j,k+1); - rank_YZ = getRankForBlock(nprocx,nprocy,nprocz,i,j+1,k+1); - rank_yz = getRankForBlock(nprocx,nprocy,nprocz,i,j-1,k-1); - rank_Yz = getRankForBlock(nprocx,nprocy,nprocz,i,j+1,k-1); - rank_yZ = getRankForBlock(nprocx,nprocy,nprocz,i,j-1,k+1); + rank_X = getRankForBlock(i+1,j,k); + rank_x = getRankForBlock(i-1,j,k); + rank_Y = getRankForBlock(i,j+1,k); + rank_y = getRankForBlock(i,j-1,k); + rank_Z = getRankForBlock(i,j,k+1); + rank_z = getRankForBlock(i,j,k-1); + rank_XY = getRankForBlock(i+1,j+1,k); + rank_xy = getRankForBlock(i-1,j-1,k); + rank_Xy = getRankForBlock(i+1,j-1,k); + rank_xY = getRankForBlock(i-1,j+1,k); + rank_XZ = getRankForBlock(i+1,j,k+1); + rank_xz = getRankForBlock(i-1,j,k-1); + rank_Xz = getRankForBlock(i+1,j,k-1); + rank_xZ = getRankForBlock(i-1,j,k+1); + rank_YZ = getRankForBlock(i,j+1,k+1); + rank_yz = getRankForBlock(i,j-1,k-1); + rank_Yz = getRankForBlock(i,j+1,k-1); + rank_yZ = getRankForBlock(i,j-1,k+1); } // Basic domain information int Nx,Ny,Nz,N; From 2becaf215e2bd8660a065c338e39fa7c5e845715 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 16:33:01 -0400 Subject: [PATCH 22/43] updated lb2_Blob_wia_mpi --- tests/lb2_Blob_wia_mpi.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index f9afe992..1212ccd3 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -624,8 +624,9 @@ int main(int argc, char **argv) } Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz); + Dm.InitializeRanks(); - InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc, + 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 ); From db967b000f6f39671202eba23ed1135969a373c2 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 16:35:53 -0400 Subject: [PATCH 23/43] Added Domain class to common/Domain.h --- common/Domain.h | 76 ++++++++++++++++++++++++++++++++++++++ tests/lb2_Blob_wia_mpi.cpp | 72 ------------------------------------ 2 files changed, 76 insertions(+), 72 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index 5ed10d21..241fb2ad 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -14,6 +14,82 @@ using namespace std; +struct Domain{ + + Domain(int nx, int ny, int nz, int rnk, int npx, int npy, int npz, + double lx, double ly, double lz){ + Nx = nx; Ny = ny; Nz = nz; + Lx = lx, Ly = ly, Lz = lz; + rank = rnk; + nprocx=npx; nprocy=npy; nprocz=npz; + N = Nx*Ny*Nz; + ID = new char [N]; + Blobs.New(Nx,Ny,Nz); + } + + + // Basic domain information + int Nx,Ny,Nz,N; + int iproc,jproc,kproc; + int nprocx,nprocy,nprocz; + double Lx,Ly,Lz; + int rank; + //********************************** + // MPI ranks for all 18 neighbors + //********************************** + int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z; + int rank_xy,rank_XY,rank_xY,rank_Xy; + int rank_xz,rank_XZ,rank_xZ,rank_Xz; + int rank_yz,rank_YZ,rank_yZ,rank_Yz; + //********************************** + + // Solid indicator function + char *ID; + // Blob information + IntArray Blobs; + + void InitializeRanks() + { + // map the rank to the block index + iproc = rank%nprocx; + jproc = (rank/nprocx)%nprocy; + kproc = rank/(nprocx*nprocy); + + // set up the neighbor ranks + int i = iproc; + int j = jproc; + int k = kproc; + rank_X = getRankForBlock(i+1,j,k); + rank_x = getRankForBlock(i-1,j,k); + rank_Y = getRankForBlock(i,j+1,k); + rank_y = getRankForBlock(i,j-1,k); + rank_Z = getRankForBlock(i,j,k+1); + rank_z = getRankForBlock(i,j,k-1); + rank_XY = getRankForBlock(i+1,j+1,k); + rank_xy = getRankForBlock(i-1,j-1,k); + rank_Xy = getRankForBlock(i+1,j-1,k); + rank_xY = getRankForBlock(i-1,j+1,k); + rank_XZ = getRankForBlock(i+1,j,k+1); + rank_xz = getRankForBlock(i-1,j,k-1); + rank_Xz = getRankForBlock(i+1,j,k-1); + rank_xZ = getRankForBlock(i-1,j,k+1); + rank_YZ = getRankForBlock(i,j+1,k+1); + rank_yz = getRankForBlock(i,j-1,k-1); + rank_Yz = getRankForBlock(i,j+1,k-1); + rank_yZ = getRankForBlock(i,j-1,k+1); + } + +private: + int getRankForBlock( int i, int j, int k ) + { + int i2 = (i+nprocx)%nprocx; + int j2 = (j+nprocy)%nprocy; + int k2 = (k+nprocz)%nprocz; + return i2 + j2*nprocx + k2*nprocx*nprocy; + } + +}; + inline void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad) { // Read in the full sphere pack diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 1212ccd3..f2cc69e9 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -22,79 +22,7 @@ using namespace std; -struct Domain{ - - Domain(int nx, int ny, int nz, int rnk, int npx, int npy, int npz, - double lx, double ly, double lz){ - Nx = nx; Ny = ny; Nz = nz; - Lx = lx, Ly = ly, Lz = lz; - rank = rnk; - nprocx=npx; nprocy=npy; nprocz=npz; - N = Nx*Ny*Nz; - ID = new char [N]; - Blobs.New(Nx,Ny,Nz); - } - void InitializeRanks() - { - // map the rank to the block index - iproc = rank%nprocx; - jproc = (rank/nprocx)%nprocy; - kproc = rank/(nprocx*nprocy); - - // set up the neighbor ranks - int i = iproc; - int j = jproc; - int k = kproc; - rank_X = getRankForBlock(i+1,j,k); - rank_x = getRankForBlock(i-1,j,k); - rank_Y = getRankForBlock(i,j+1,k); - rank_y = getRankForBlock(i,j-1,k); - rank_Z = getRankForBlock(i,j,k+1); - rank_z = getRankForBlock(i,j,k-1); - rank_XY = getRankForBlock(i+1,j+1,k); - rank_xy = getRankForBlock(i-1,j-1,k); - rank_Xy = getRankForBlock(i+1,j-1,k); - rank_xY = getRankForBlock(i-1,j+1,k); - rank_XZ = getRankForBlock(i+1,j,k+1); - rank_xz = getRankForBlock(i-1,j,k-1); - rank_Xz = getRankForBlock(i+1,j,k-1); - rank_xZ = getRankForBlock(i-1,j,k+1); - rank_YZ = getRankForBlock(i,j+1,k+1); - rank_yz = getRankForBlock(i,j-1,k-1); - rank_Yz = getRankForBlock(i,j+1,k-1); - rank_yZ = getRankForBlock(i,j-1,k+1); - } - // Basic domain information - int Nx,Ny,Nz,N; - int iproc,jproc,kproc; - int nprocx,nprocy,nprocz; - double Lx,Ly,Lz; - int rank; - //********************************** - // MPI ranks for all 18 neighbors - //********************************** - int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z; - int rank_xy,rank_XY,rank_xY,rank_Xy; - int rank_xz,rank_XZ,rank_xZ,rank_Xz; - int rank_yz,rank_YZ,rank_yZ,rank_Yz; - //********************************** - - // Solid indicator function - char *ID; - // Blob information - IntArray Blobs; - -private: - int getRankForBlock( int i, int j, int k ) - { - int i2 = (i+nprocx)%nprocx; - int j2 = (j+nprocy)%nprocy; - int k2 = (k+nprocz)%nprocz; - return i2 + j2*nprocx + k2*nprocx*nprocy; - } - -}; struct BlobInfo{ From 54d441241262b632b9d67474a5c0f1ae85248f92 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 25 Jun 2014 16:45:42 -0400 Subject: [PATCH 24/43] updated lb2_Blob_wia_mpi --- common/Domain.h | 3 +-- tests/lb2_Blob_wia_mpi.cpp | 2 -- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index 241fb2ad..b268f47c 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -18,7 +18,7 @@ struct Domain{ Domain(int nx, int ny, int nz, int rnk, int npx, int npy, int npz, double lx, double ly, double lz){ - Nx = nx; Ny = ny; Nz = nz; + Nx = nx+2; Ny = ny+2; Nz = nz+2; Lx = lx, Ly = ly, Lz = lz; rank = rnk; nprocx=npx; nprocy=npy; nprocz=npz; @@ -27,7 +27,6 @@ struct Domain{ Blobs.New(Nx,Ny,Nz); } - // Basic domain information int Nx,Ny,Nz,N; int iproc,jproc,kproc; diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index f2cc69e9..141425b0 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -22,8 +22,6 @@ using namespace std; - - struct BlobInfo{ BlobInfo(int Nx, int Ny,int Nz){ From 4cdbea6f3189e4cf656ee18b399ceb573cce182e Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 18:16:12 -0400 Subject: [PATCH 25/43] Basic things to tests/lb2_Blob_wia_mpi.cpp --- tests/lb2_Blob_wia_mpi.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index 141425b0..ade3389f 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -740,6 +740,7 @@ int main(int argc, char **argv) for (i=0;i Date: Wed, 25 Jun 2014 18:33:01 -0400 Subject: [PATCH 26/43] Modified signed distance in the reservoirs when pBC is true --- tests/lb2_Color_wia_mpi.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/lb2_Color_wia_mpi.cpp b/tests/lb2_Color_wia_mpi.cpp index 7502b98a..7f521b43 100644 --- a/tests/lb2_Color_wia_mpi.cpp +++ b/tests/lb2_Color_wia_mpi.cpp @@ -467,6 +467,7 @@ int main(int argc, char **argv) for (i=0;i Date: Wed, 25 Jun 2014 20:36:35 -0400 Subject: [PATCH 27/43] Added InitalizeRanks() and CommInit() to common/Domain.h --- common/Domain.h | 171 +++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 140 insertions(+), 31 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index b268f47c..0ee27ea0 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -41,42 +41,25 @@ struct Domain{ int rank_xz,rank_XZ,rank_xZ,rank_Xz; int rank_yz,rank_YZ,rank_yZ,rank_Yz; //********************************** + //...................................................................................... + // Get the actual D3Q19 communication counts (based on location of solid phase) + // Discrete velocity set symmetry implies the sendcount = recvcount + int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z; + int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ; + int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ; + //...................................................................................... + int *sendList_x, *sendList_y, *sendList_z, *sendList_X, *sendList_Y, *sendList_Z; + int *sendList_xy, *sendList_yz, *sendList_xz, *sendList_Xy, *sendList_Yz, *sendList_xZ; + int *sendList_xY, *sendList_yZ, *sendList_Xz, *sendList_XY, *sendList_YZ, *sendList_XZ; + //...................................................................................... // Solid indicator function - char *ID; + char *id; // Blob information IntArray Blobs; - void InitializeRanks() - { - // map the rank to the block index - iproc = rank%nprocx; - jproc = (rank/nprocx)%nprocy; - kproc = rank/(nprocx*nprocy); - - // set up the neighbor ranks - int i = iproc; - int j = jproc; - int k = kproc; - rank_X = getRankForBlock(i+1,j,k); - rank_x = getRankForBlock(i-1,j,k); - rank_Y = getRankForBlock(i,j+1,k); - rank_y = getRankForBlock(i,j-1,k); - rank_Z = getRankForBlock(i,j,k+1); - rank_z = getRankForBlock(i,j,k-1); - rank_XY = getRankForBlock(i+1,j+1,k); - rank_xy = getRankForBlock(i-1,j-1,k); - rank_Xy = getRankForBlock(i+1,j-1,k); - rank_xY = getRankForBlock(i-1,j+1,k); - rank_XZ = getRankForBlock(i+1,j,k+1); - rank_xz = getRankForBlock(i-1,j,k-1); - rank_Xz = getRankForBlock(i+1,j,k-1); - rank_xZ = getRankForBlock(i-1,j,k+1); - rank_YZ = getRankForBlock(i,j+1,k+1); - rank_yz = getRankForBlock(i,j-1,k-1); - rank_Yz = getRankForBlock(i,j+1,k-1); - rank_yZ = getRankForBlock(i,j-1,k+1); - } + void InitializeRanks(); + void CommInit(); private: int getRankForBlock( int i, int j, int k ) @@ -89,6 +72,132 @@ private: }; +void Domain::InitializeRanks() +{ + // map the rank to the block index + iproc = rank%nprocx; + jproc = (rank/nprocx)%nprocy; + kproc = rank/(nprocx*nprocy); + + // set up the neighbor ranks + int i = iproc; + int j = jproc; + int k = kproc; + rank_X = getRankForBlock(i+1,j,k); + rank_x = getRankForBlock(i-1,j,k); + rank_Y = getRankForBlock(i,j+1,k); + rank_y = getRankForBlock(i,j-1,k); + rank_Z = getRankForBlock(i,j,k+1); + rank_z = getRankForBlock(i,j,k-1); + rank_XY = getRankForBlock(i+1,j+1,k); + rank_xy = getRankForBlock(i-1,j-1,k); + rank_Xy = getRankForBlock(i+1,j-1,k); + rank_xY = getRankForBlock(i-1,j+1,k); + rank_XZ = getRankForBlock(i+1,j,k+1); + rank_xz = getRankForBlock(i-1,j,k-1); + rank_Xz = getRankForBlock(i+1,j,k-1); + rank_xZ = getRankForBlock(i-1,j,k+1); + rank_YZ = getRankForBlock(i,j+1,k+1); + rank_yz = getRankForBlock(i,j-1,k-1); + rank_Yz = getRankForBlock(i,j+1,k-1); + rank_yZ = getRankForBlock(i,j-1,k+1); +} + + +void Domain::CommInit(){ + int i,j,k,n; + sendCount_x = sendCount_y = sendCount_z = sendCount_X = sendCount_Y = sendCount_Z = 0; + sendCount_xy = sendCount_yz = sendCount_xz = sendCount_Xy = sendCount_Yz = sendCount_xZ = 0; + sendCount_xY = sendCount_yZ = sendCount_Xz = sendCount_XY = sendCount_YZ = sendCount_XZ = 0; + //...................................................................................... + for (k=0; k Date: Wed, 25 Jun 2014 20:38:43 -0400 Subject: [PATCH 28/43] fixed little bug in common/Domain,h --- common/Domain.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index 0ee27ea0..45df2617 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -23,7 +23,7 @@ struct Domain{ rank = rnk; nprocx=npx; nprocy=npy; nprocz=npz; N = Nx*Ny*Nz; - ID = new char [N]; + id = new char [N]; Blobs.New(Nx,Ny,Nz); } From 580d98b5c8a537a982409f42b86a486f2839e748 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 20:44:04 -0400 Subject: [PATCH 29/43] Added integer send buffer allocation to Domain::CommInit common/Domain.h --- common/Domain.h | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index 45df2617..766de552 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -44,6 +44,7 @@ struct Domain{ //...................................................................................... // Get the actual D3Q19 communication counts (based on location of solid phase) // Discrete velocity set symmetry implies the sendcount = recvcount + //...................................................................................... int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z; int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ; int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ; @@ -52,6 +53,10 @@ struct Domain{ int *sendList_xy, *sendList_yz, *sendList_xz, *sendList_Xy, *sendList_Yz, *sendList_xZ; int *sendList_xY, *sendList_yZ, *sendList_Xz, *sendList_XY, *sendList_YZ, *sendList_XZ; //...................................................................................... + int *sendBuf_x, *sendBuf_y, *sendBuf_z, *sendBuf_X, *sendBuf_Y, *sendBuf_Z; + int *sendBuf_xy, *sendBuf_yz, *sendBuf_xz, *sendBuf_Xy, *sendBuf_Yz, *sendBuf_xZ; + int *sendBuf_xY, *sendBuf_yZ, *sendBuf_Xz, *sendBuf_XY, *sendBuf_YZ, *sendBuf_XZ; + //...................................................................................... // Solid indicator function char *id; @@ -141,7 +146,7 @@ void Domain::CommInit(){ } } } - // allocate send buffers + // allocate send lists sendList_x = new int [sendCount_x]; sendList_y = new int [sendCount_y]; sendList_z = new int [sendCount_z]; @@ -196,6 +201,25 @@ void Domain::CommInit(){ } } } + // allocate send buffers + sendBuf_x = new int [sendCount_x]; + sendBuf_y = new int [sendCount_y]; + sendBuf_z = new int [sendCount_z]; + sendBuf_X = new int [sendCount_X]; + sendBuf_Y = new int [sendCount_Y]; + sendBuf_Z = new int [sendCount_Z]; + sendBuf_xy = new int [sendCount_xy]; + sendBuf_yz = new int [sendCount_yz]; + sendBuf_xz = new int [sendCount_xz]; + sendBuf_Xy = new int [sendCount_Xy]; + sendBuf_Yz = new int [sendCount_Yz]; + sendBuf_xZ = new int [sendCount_xZ]; + sendBuf_xY = new int [sendCount_xY]; + sendBuf_yZ = new int [sendCount_yZ]; + sendBuf_Xz = new int [sendCount_Xz]; + sendBuf_XY = new int [sendCount_XY]; + sendBuf_YZ = new int [sendCount_YZ]; + sendBuf_XZ = new int [sendCount_XZ]; } inline void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad) From 3599a568dc4ceba937c22ba6f7b43b380a6eeea1 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 20:49:37 -0400 Subject: [PATCH 30/43] Adding mpi to Domain.h --- common/Domain.h | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index 766de552..06dbfdf0 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -8,7 +8,7 @@ #include #include // std::exception #include -#include +#include #include "common/Utilities.h" @@ -56,6 +56,14 @@ struct Domain{ int *sendBuf_x, *sendBuf_y, *sendBuf_z, *sendBuf_X, *sendBuf_Y, *sendBuf_Z; int *sendBuf_xy, *sendBuf_yz, *sendBuf_xz, *sendBuf_Xy, *sendBuf_Yz, *sendBuf_xZ; int *sendBuf_xY, *sendBuf_yZ, *sendBuf_Xz, *sendBuf_XY, *sendBuf_YZ, *sendBuf_XZ; + //...................................................................................... + int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z; + int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ; + int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ; + //...................................................................................... + int *recvBuf_x, *recvBuf_y, *recvBuf_z, *recvBuf_X, *recvBuf_Y, *recvBuf_Z; + int *recvBuf_xy, *recvBuf_yz, *recvBuf_xz, *recvBuf_Xy, *recvBuf_Yz, *recvBuf_xZ; + int *recvBuf_xY, *recvBuf_yZ, *recvBuf_Xz, *recvBuf_XY, *recvBuf_YZ, *recvBuf_XZ; //...................................................................................... // Solid indicator function From 15c3ae7cd66f45782849b58cd44ca3d13e85356a Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 20:58:29 -0400 Subject: [PATCH 31/43] MPI Initialization of send and recieve countounts in class Domain --- common/Domain.h | 58 ++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 55 insertions(+), 3 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index 06dbfdf0..a3850260 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -26,13 +26,15 @@ struct Domain{ id = new char [N]; Blobs.New(Nx,Ny,Nz); } - + // Basic domain information int Nx,Ny,Nz,N; int iproc,jproc,kproc; int nprocx,nprocy,nprocz; double Lx,Ly,Lz; int rank; + MPI_Comm Communicator; + //********************************** // MPI ranks for all 18 neighbors //********************************** @@ -72,7 +74,7 @@ struct Domain{ IntArray Blobs; void InitializeRanks(); - void CommInit(); + void CommInit(MPI_Comm comm); private: int getRankForBlock( int i, int j, int k ) @@ -117,8 +119,15 @@ void Domain::InitializeRanks() } -void Domain::CommInit(){ +void Domain::CommInit(MPI_Comm Communicator){ int i,j,k,n; + int sendtag = 21; + int recvtag = 21; + + //...................................................................................... + MPI_Request req1[18], req2[18]; + MPI_Status stat1[18],stat2[18]; + //...................................................................................... sendCount_x = sendCount_y = sendCount_z = sendCount_X = sendCount_Y = sendCount_Z = 0; sendCount_xy = sendCount_yz = sendCount_xz = sendCount_Xy = sendCount_Yz = sendCount_xZ = 0; sendCount_xY = sendCount_yZ = sendCount_Xz = sendCount_XY = sendCount_YZ = sendCount_XZ = 0; @@ -228,6 +237,49 @@ void Domain::CommInit(){ sendBuf_XY = new int [sendCount_XY]; sendBuf_YZ = new int [sendCount_YZ]; sendBuf_XZ = new int [sendCount_XZ]; + + MPI_Isend(&sendCount_x, 1,MPI_INT,rank_x,sendtag+0,Communicator,&req1[0]); + MPI_Irecv(&recvCount_X, 1,MPI_INT,rank_X,recvtag+0,Communicator,&req2[0]); + MPI_Isend(&sendCount_X, 1,MPI_INT,rank_X,sendtag+1,Communicator,&req1[1]); + MPI_Irecv(&recvCount_x, 1,MPI_INT,rank_x,recvtag+1,Communicator,&req2[1]); + MPI_Isend(&sendCount_y, 1,MPI_INT,rank_y,sendtag+2,Communicator,&req1[2]); + MPI_Irecv(&recvCount_Y, 1,MPI_INT,rank_Y,recvtag+2,Communicator,&req2[2]); + MPI_Isend(&sendCount_Y, 1,MPI_INT,rank_Y,sendtag+3,Communicator,&req1[3]); + MPI_Irecv(&recvCount_y, 1,MPI_INT,rank_y,recvtag+3,Communicator,&req2[3]); + MPI_Isend(&sendCount_z, 1,MPI_INT,rank_z,sendtag+4,Communicator,&req1[4]); + MPI_Irecv(&recvCount_Z, 1,MPI_INT,rank_Z,recvtag+4,Communicator,&req2[4]); + MPI_Isend(&sendCount_Z, 1,MPI_INT,rank_Z,sendtag+5,Communicator,&req1[5]); + MPI_Irecv(&recvCount_z, 1,MPI_INT,rank_z,recvtag+5,Communicator,&req2[5]); + + MPI_Isend(&sendCount_xy, 1,MPI_INT,rank_xy,sendtag+6,Communicator,&req1[6]); + MPI_Irecv(&recvCount_XY, 1,MPI_INT,rank_XY,recvtag+6,Communicator,&req2[6]); + MPI_Isend(&sendCount_XY, 1,MPI_INT,rank_XY,sendtag+7,Communicator,&req1[7]); + MPI_Irecv(&recvCount_xy, 1,MPI_INT,rank_xy,recvtag+7,Communicator,&req2[7]); + MPI_Isend(&sendCount_Xy, 1,MPI_INT,rank_Xy,sendtag+8,Communicator,&req1[8]); + MPI_Irecv(&recvCount_xY, 1,MPI_INT,rank_xY,recvtag+8,Communicator,&req2[8]); + MPI_Isend(&sendCount_xY, 1,MPI_INT,rank_xY,sendtag+9,Communicator,&req1[9]); + MPI_Irecv(&recvCount_Xy, 1,MPI_INT,rank_Xy,recvtag+9,Communicator,&req2[9]); + + MPI_Isend(&sendCount_xz, 1,MPI_INT,rank_xz,sendtag+10,Communicator,&req1[10]); + MPI_Irecv(&recvCount_XZ, 1,MPI_INT,rank_XZ,recvtag+10,Communicator,&req2[10]); + MPI_Isend(&sendCount_XZ, 1,MPI_INT,rank_XZ,sendtag+11,Communicator,&req1[11]); + MPI_Irecv(&recvCount_xz, 1,MPI_INT,rank_xz,recvtag+11,Communicator,&req2[11]); + MPI_Isend(&sendCount_Xz, 1,MPI_INT,rank_Xz,sendtag+12,Communicator,&req1[12]); + MPI_Irecv(&recvCount_xZ, 1,MPI_INT,rank_xZ,recvtag+12,Communicator,&req2[12]); + MPI_Isend(&sendCount_xZ, 1,MPI_INT,rank_xZ,sendtag+13,Communicator,&req1[13]); + MPI_Irecv(&recvCount_Xz, 1,MPI_INT,rank_Xz,recvtag+13,Communicator,&req2[13]); + + MPI_Isend(&sendCount_yz, 1,MPI_INT,rank_yz,sendtag+14,Communicator,&req1[14]); + MPI_Irecv(&recvCount_YZ, 1,MPI_INT,rank_YZ,recvtag+14,Communicator,&req2[14]); + MPI_Isend(&sendCount_YZ, 1,MPI_INT,rank_YZ,sendtag+15,Communicator,&req1[15]); + MPI_Irecv(&recvCount_yz, 1,MPI_INT,rank_yz,recvtag+15,Communicator,&req2[15]); + MPI_Isend(&sendCount_Yz, 1,MPI_INT,rank_Yz,sendtag+16,Communicator,&req1[16]); + MPI_Irecv(&recvCount_yZ, 1,MPI_INT,rank_yZ,recvtag+16,Communicator,&req2[16]); + MPI_Isend(&sendCount_yZ, 1,MPI_INT,rank_yZ,sendtag+17,Communicator,&req1[17]); + MPI_Irecv(&recvCount_Yz, 1,MPI_INT,rank_Yz,recvtag+17,Communicator,&req2[17]); + MPI_Waitall(18,req1,stat1); + MPI_Waitall(18,req2,stat2); + MPI_Barrier(Communicator); } inline void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad) From a664a6304c7f552771f6b56034e24fa7faa2a41e Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 20:59:22 -0400 Subject: [PATCH 32/43] fixed bug in Dm.CommInit() --- tests/lb2_Blob_wia_mpi.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/lb2_Blob_wia_mpi.cpp b/tests/lb2_Blob_wia_mpi.cpp index ade3389f..c86e0777 100644 --- a/tests/lb2_Blob_wia_mpi.cpp +++ b/tests/lb2_Blob_wia_mpi.cpp @@ -551,6 +551,7 @@ int main(int argc, char **argv) Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz); Dm.InitializeRanks(); + Dm.CommInit(MPI_COMM_WORLD); InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc, rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z, From 1264a5c063509c53d08971c55888fcfd57eae983 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 21:04:43 -0400 Subject: [PATCH 33/43] Mapped recieve counts in Domain::CommInit() --- common/Domain.h | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/common/Domain.h b/common/Domain.h index a3850260..d691128c 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -280,6 +280,30 @@ void Domain::CommInit(MPI_Comm Communicator){ MPI_Waitall(18,req1,stat1); MPI_Waitall(18,req2,stat2); MPI_Barrier(Communicator); + + //...................................................................................... + for (int idx=0; idx Date: Wed, 25 Jun 2014 21:09:35 -0400 Subject: [PATCH 34/43] Updated Domain::CommInit() --- common/Domain.h | 73 +++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 65 insertions(+), 8 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index d691128c..c460e654 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -63,6 +63,10 @@ struct Domain{ int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ; int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ; //...................................................................................... + int *recvList_x, *recvList_y, *recvList_z, *recvList_X, *recvList_Y, *recvList_Z; + int *recvList_xy, *recvList_yz, *recvList_xz, *recvList_Xy, *recvList_Yz, *recvList_xZ; + int *recvList_xY, *recvList_yZ, *recvList_Xz, *recvList_XY, *recvList_YZ, *recvList_XZ; + //...................................................................................... int *recvBuf_x, *recvBuf_y, *recvBuf_z, *recvBuf_X, *recvBuf_Y, *recvBuf_Z; int *recvBuf_xy, *recvBuf_yz, *recvBuf_xz, *recvBuf_Xy, *recvBuf_Yz, *recvBuf_xZ; int *recvBuf_xY, *recvBuf_yZ, *recvBuf_Xz, *recvBuf_XY, *recvBuf_YZ, *recvBuf_XZ; @@ -250,7 +254,6 @@ void Domain::CommInit(MPI_Comm Communicator){ MPI_Irecv(&recvCount_Z, 1,MPI_INT,rank_Z,recvtag+4,Communicator,&req2[4]); MPI_Isend(&sendCount_Z, 1,MPI_INT,rank_Z,sendtag+5,Communicator,&req1[5]); MPI_Irecv(&recvCount_z, 1,MPI_INT,rank_z,recvtag+5,Communicator,&req2[5]); - MPI_Isend(&sendCount_xy, 1,MPI_INT,rank_xy,sendtag+6,Communicator,&req1[6]); MPI_Irecv(&recvCount_XY, 1,MPI_INT,rank_XY,recvtag+6,Communicator,&req2[6]); MPI_Isend(&sendCount_XY, 1,MPI_INT,rank_XY,sendtag+7,Communicator,&req1[7]); @@ -259,7 +262,6 @@ void Domain::CommInit(MPI_Comm Communicator){ MPI_Irecv(&recvCount_xY, 1,MPI_INT,rank_xY,recvtag+8,Communicator,&req2[8]); MPI_Isend(&sendCount_xY, 1,MPI_INT,rank_xY,sendtag+9,Communicator,&req1[9]); MPI_Irecv(&recvCount_Xy, 1,MPI_INT,rank_Xy,recvtag+9,Communicator,&req2[9]); - MPI_Isend(&sendCount_xz, 1,MPI_INT,rank_xz,sendtag+10,Communicator,&req1[10]); MPI_Irecv(&recvCount_XZ, 1,MPI_INT,rank_XZ,recvtag+10,Communicator,&req2[10]); MPI_Isend(&sendCount_XZ, 1,MPI_INT,rank_XZ,sendtag+11,Communicator,&req1[11]); @@ -268,7 +270,6 @@ void Domain::CommInit(MPI_Comm Communicator){ MPI_Irecv(&recvCount_xZ, 1,MPI_INT,rank_xZ,recvtag+12,Communicator,&req2[12]); MPI_Isend(&sendCount_xZ, 1,MPI_INT,rank_xZ,sendtag+13,Communicator,&req1[13]); MPI_Irecv(&recvCount_Xz, 1,MPI_INT,rank_Xz,recvtag+13,Communicator,&req2[13]); - MPI_Isend(&sendCount_yz, 1,MPI_INT,rank_yz,sendtag+14,Communicator,&req1[14]); MPI_Irecv(&recvCount_YZ, 1,MPI_INT,rank_YZ,recvtag+14,Communicator,&req2[14]); MPI_Isend(&sendCount_YZ, 1,MPI_INT,rank_YZ,sendtag+15,Communicator,&req1[15]); @@ -280,30 +281,86 @@ void Domain::CommInit(MPI_Comm Communicator){ MPI_Waitall(18,req1,stat1); MPI_Waitall(18,req2,stat2); MPI_Barrier(Communicator); - + //...................................................................................... + // recv buffers + recvList_x = new int [recvCount_x]; + recvList_y = new int [recvCount_y]; + recvList_z = new int [recvCount_z]; + recvList_X = new int [recvCount_X]; + recvList_Y = new int [recvCount_Y]; + recvList_Z = new int [recvCount_Z]; + recvList_xy = new int [recvCount_xy]; + recvList_yz = new int [recvCount_yz]; + recvList_xz = new int [recvCount_xz]; + recvList_Xy = new int [recvCount_Xy]; + recvList_Yz = new int [recvCount_Yz]; + recvList_xZ = new int [recvCount_xZ]; + recvList_xY = new int [recvCount_xY]; + recvList_yZ = new int [recvCount_yZ]; + recvList_Xz = new int [recvCount_Xz]; + recvList_XY = new int [recvCount_XY]; + recvList_YZ = new int [recvCount_YZ]; + recvList_XZ = new int [recvCount_XZ]; + //...................................................................................... + MPI_Isend(sendList_x, sendCount_x,MPI_INT,rank_x,sendtag,Communicator,&req1[0]); + MPI_Irecv(recvList_X, recvCount_X,MPI_INT,rank_X,recvtag,Communicator,&req2[0]); + MPI_Isend(sendList_X, sendCount_X,MPI_INT,rank_X,sendtag,Communicator,&req1[1]); + MPI_Irecv(recvList_x, recvCount_x,MPI_INT,rank_x,recvtag,Communicator,&req2[1]); + MPI_Isend(sendList_y, sendCount_y,MPI_INT,rank_y,sendtag,Communicator,&req1[2]); + MPI_Irecv(recvList_Y, recvCount_Y,MPI_INT,rank_Y,recvtag,Communicator,&req2[2]); + MPI_Isend(sendList_Y, sendCount_Y,MPI_INT,rank_Y,sendtag,Communicator,&req1[3]); + MPI_Irecv(recvList_y, recvCount_y,MPI_INT,rank_y,recvtag,Communicator,&req2[3]); + MPI_Isend(sendList_z, sendCount_z,MPI_INT,rank_z,sendtag,Communicator,&req1[4]); + MPI_Irecv(recvList_Z, recvCount_Z,MPI_INT,rank_Z,recvtag,Communicator,&req2[4]); + MPI_Isend(sendList_Z, sendCount_Z,MPI_INT,rank_Z,sendtag,Communicator,&req1[5]); + MPI_Irecv(recvList_z, recvCount_z,MPI_INT,rank_z,recvtag,Communicator,&req2[5]); + MPI_Isend(sendList_xy, sendCount_xy,MPI_INT,rank_xy,sendtag,Communicator,&req1[6]); + MPI_Irecv(recvList_XY, recvCount_XY,MPI_INT,rank_XY,recvtag,Communicator,&req2[6]); + MPI_Isend(sendList_XY, sendCount_XY,MPI_INT,rank_XY,sendtag,Communicator,&req1[7]); + MPI_Irecv(recvList_xy, recvCount_xy,MPI_INT,rank_xy,recvtag,Communicator,&req2[7]); + MPI_Isend(sendList_Xy, sendCount_Xy,MPI_INT,rank_Xy,sendtag,Communicator,&req1[8]); + MPI_Irecv(recvList_xY, recvCount_xY,MPI_INT,rank_xY,recvtag,Communicator,&req2[8]); + MPI_Isend(sendList_xY, sendCount_xY,MPI_INT,rank_xY,sendtag,Communicator,&req1[9]); + MPI_Irecv(recvList_Xy, recvCount_Xy,MPI_INT,rank_Xy,recvtag,Communicator,&req2[9]); + MPI_Isend(sendList_xz, sendCount_xz,MPI_INT,rank_xz,sendtag,Communicator,&req1[10]); + MPI_Irecv(recvList_XZ, recvCount_XZ,MPI_INT,rank_XZ,recvtag,Communicator,&req2[10]); + MPI_Isend(sendList_XZ, sendCount_XZ,MPI_INT,rank_XZ,sendtag,Communicator,&req1[11]); + MPI_Irecv(recvList_xz, recvCount_xz,MPI_INT,rank_xz,recvtag,Communicator,&req2[11]); + MPI_Isend(sendList_Xz, sendCount_Xz,MPI_INT,rank_Xz,sendtag,Communicator,&req1[12]); + MPI_Irecv(recvList_xZ, recvCount_xZ,MPI_INT,rank_xZ,recvtag,Communicator,&req2[12]); + MPI_Isend(sendList_xZ, sendCount_xZ,MPI_INT,rank_xZ,sendtag,Communicator,&req1[13]); + MPI_Irecv(recvList_Xz, recvCount_Xz,MPI_INT,rank_Xz,recvtag,Communicator,&req2[13]); + MPI_Isend(sendList_yz, sendCount_yz,MPI_INT,rank_yz,sendtag,Communicator,&req1[14]); + MPI_Irecv(recvList_YZ, recvCount_YZ,MPI_INT,rank_YZ,recvtag,Communicator,&req2[14]); + MPI_Isend(sendList_YZ, sendCount_YZ,MPI_INT,rank_YZ,sendtag,Communicator,&req1[15]); + MPI_Irecv(recvList_yz, recvCount_yz,MPI_INT,rank_yz,recvtag,Communicator,&req2[15]); + MPI_Isend(sendList_Yz, sendCount_Yz,MPI_INT,rank_Yz,sendtag,Communicator,&req1[16]); + MPI_Irecv(recvList_yZ, recvCount_yZ,MPI_INT,rank_yZ,recvtag,Communicator,&req2[16]); + MPI_Isend(sendList_yZ, sendCount_yZ,MPI_INT,rank_yZ,sendtag,Communicator,&req1[17]); + MPI_Irecv(recvList_Yz, recvCount_Yz,MPI_INT,rank_Yz,recvtag,Communicator,&req2[17]); + MPI_Waitall(18,req1,stat1); + MPI_Waitall(18,req2,stat2); //...................................................................................... for (int idx=0; idx Date: Wed, 25 Jun 2014 21:11:53 -0400 Subject: [PATCH 35/43] Updated Domain::CommInit() --- common/Domain.h | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index c460e654..f29eef2e 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -240,8 +240,7 @@ void Domain::CommInit(MPI_Comm Communicator){ sendBuf_Xz = new int [sendCount_Xz]; sendBuf_XY = new int [sendCount_XY]; sendBuf_YZ = new int [sendCount_YZ]; - sendBuf_XZ = new int [sendCount_XZ]; - + sendBuf_XZ = new int [sendCount_XZ]; MPI_Isend(&sendCount_x, 1,MPI_INT,rank_x,sendtag+0,Communicator,&req1[0]); MPI_Irecv(&recvCount_X, 1,MPI_INT,rank_X,recvtag+0,Communicator,&req2[0]); MPI_Isend(&sendCount_X, 1,MPI_INT,rank_X,sendtag+1,Communicator,&req1[1]); @@ -360,6 +359,26 @@ void Domain::CommInit(MPI_Comm Communicator){ for (int idx=0; idx Date: Wed, 25 Jun 2014 21:29:18 -0400 Subject: [PATCH 36/43] Updated Domain::CommInit() --- common/Domain.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index f29eef2e..04e8f2f0 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -127,7 +127,6 @@ void Domain::CommInit(MPI_Comm Communicator){ int i,j,k,n; int sendtag = 21; int recvtag = 21; - //...................................................................................... MPI_Request req1[18], req2[18]; MPI_Status stat1[18],stat2[18]; @@ -240,7 +239,8 @@ void Domain::CommInit(MPI_Comm Communicator){ sendBuf_Xz = new int [sendCount_Xz]; sendBuf_XY = new int [sendCount_XY]; sendBuf_YZ = new int [sendCount_YZ]; - sendBuf_XZ = new int [sendCount_XZ]; + sendBuf_XZ = new int [sendCount_XZ]; + //...................................................................................... MPI_Isend(&sendCount_x, 1,MPI_INT,rank_x,sendtag+0,Communicator,&req1[0]); MPI_Irecv(&recvCount_X, 1,MPI_INT,rank_X,recvtag+0,Communicator,&req2[0]); MPI_Isend(&sendCount_X, 1,MPI_INT,rank_X,sendtag+1,Communicator,&req1[1]); From a42034ec1c7740d2de29ca243e1d3e245f13d19f Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 21:40:24 -0400 Subject: [PATCH 37/43] Updated Domain::CommInit() --- common/Domain.h | 101 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) diff --git a/common/Domain.h b/common/Domain.h index 04e8f2f0..e2696444 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -79,8 +79,27 @@ struct Domain{ void InitializeRanks(); void CommInit(MPI_Comm comm); + void BlobComm(MPI_Comm comm); private: + void PackBlobData(int *list, int count, int *sendbuf, int *data){ + // 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 Date: Wed, 25 Jun 2014 21:41:44 -0400 Subject: [PATCH 38/43] Updated Domain::CommInit() --- common/Domain.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index e2696444..f0102f6e 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -401,7 +401,7 @@ void Domain::CommInit(MPI_Comm Communicator){ } -Domain::BlobComm(MPI_Comm Communicator){ +void Domain::BlobComm(MPI_Comm Communicator){ //...................................................................................... int sendtag, recvtag; sendtag = recvtag = 51; From 1d2a392fe9015d84e805516a684cd2b19eba6026 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 21:42:38 -0400 Subject: [PATCH 39/43] Updated Domain::CommInit() --- common/Domain.h | 72 ++++++++++++++++++++++++------------------------- 1 file changed, 36 insertions(+), 36 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index f0102f6e..70790799 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -406,24 +406,24 @@ void Domain::BlobComm(MPI_Comm Communicator){ int sendtag, recvtag; sendtag = recvtag = 51; //...................................................................................... - PackBlobsData(sendList_x, sendCount_x ,sendBuf_x, Blobs.data); - PackBlobsData(sendList_X, sendCount_X ,sendBuf_X, Blobs.data); - PackBlobsData(sendList_y, sendCount_y ,sendBuf_y, Blobs.data); - PackBlobsData(sendList_Y, sendCount_Y ,sendBuf_Y, Blobs.data); - PackBlobsData(sendList_z, sendCount_z ,sendBuf_z, Blobs.data); - PackBlobsData(sendList_Z, sendCount_Z ,sendBuf_Z, Blobs.data); - PackBlobsData(sendList_xy, sendCount_xy ,sendBuf_xy, Blobs.data); - PackBlobsData(sendList_Xy, sendCount_Xy ,sendBuf_Xy, Blobs.data); - PackBlobsData(sendList_xY, sendCount_xY ,sendBuf_xY, Blobs.data); - PackBlobsData(sendList_XY, sendCount_XY ,sendBuf_XY, Blobs.data); - PackBlobsData(sendList_xz, sendCount_xz ,sendBuf_xz, Blobs.data); - PackBlobsData(sendList_Xz, sendCount_Xz ,sendBuf_Xz, Blobs.data); - PackBlobsData(sendList_xZ, sendCount_xZ ,sendBuf_xZ, Blobs.data); - PackBlobsData(sendList_XZ, sendCount_XZ ,sendBuf_XZ, Blobs.data); - PackBlobsData(sendList_yz, sendCount_yz ,sendBuf_yz, Blobs.data); - PackBlobsData(sendList_Yz, sendCount_Yz ,sendBuf_Yz, Blobs.data); - PackBlobsData(sendList_yZ, sendCount_yZ ,sendBuf_yZ, Blobs.data); - PackBlobsData(sendList_YZ, sendCount_YZ ,sendBuf_YZ, Blobs.data); + PackBlobData(sendList_x, sendCount_x ,sendBuf_x, Blobs.data); + PackBlobData(sendList_X, sendCount_X ,sendBuf_X, Blobs.data); + PackBlobData(sendList_y, sendCount_y ,sendBuf_y, Blobs.data); + PackBlobData(sendList_Y, sendCount_Y ,sendBuf_Y, Blobs.data); + PackBlobData(sendList_z, sendCount_z ,sendBuf_z, Blobs.data); + PackBlobData(sendList_Z, sendCount_Z ,sendBuf_Z, Blobs.data); + PackBlobData(sendList_xy, sendCount_xy ,sendBuf_xy, Blobs.data); + PackBlobData(sendList_Xy, sendCount_Xy ,sendBuf_Xy, Blobs.data); + PackBlobData(sendList_xY, sendCount_xY ,sendBuf_xY, Blobs.data); + PackBlobData(sendList_XY, sendCount_XY ,sendBuf_XY, Blobs.data); + PackBlobData(sendList_xz, sendCount_xz ,sendBuf_xz, Blobs.data); + PackBlobData(sendList_Xz, sendCount_Xz ,sendBuf_Xz, Blobs.data); + PackBlobData(sendList_xZ, sendCount_xZ ,sendBuf_xZ, Blobs.data); + PackBlobData(sendList_XZ, sendCount_XZ ,sendBuf_XZ, Blobs.data); + PackBlobData(sendList_yz, sendCount_yz ,sendBuf_yz, Blobs.data); + PackBlobData(sendList_Yz, sendCount_Yz ,sendBuf_Yz, Blobs.data); + PackBlobData(sendList_yZ, sendCount_yZ ,sendBuf_yZ, Blobs.data); + PackBlobData(sendList_YZ, sendCount_YZ ,sendBuf_YZ, Blobs.data); //...................................................................................... MPI_Sendrecv(sendBuf_x,sendCount_x,MPI_INT,rank_x,sendtag, recvBuf_X,recvCount_X,MPI_INT,rank_X,recvtag,Communicator,MPI_STATUS_IGNORE); @@ -462,24 +462,24 @@ void Domain::BlobComm(MPI_Comm Communicator){ MPI_Sendrecv(sendBuf_yZ,sendCount_yZ,MPI_INT,rank_yZ,sendtag, recvBuf_Yz,recvCount_Yz,MPI_INT,rank_Yz,recvtag,Communicator,MPI_STATUS_IGNORE); //........................................................................................ - UnpackBlobsData(recvList_x, recvCount_x ,recvBuf_x, Blobs.data); - UnpackBlobsData(recvList_X, recvCount_X ,recvBuf_X, Blobs.data); - UnpackBlobsData(recvList_y, recvCount_y ,recvBuf_y, Blobs.data); - UnpackBlobsData(recvList_Y, recvCount_Y ,recvBuf_Y, Blobs.data); - UnpackBlobsData(recvList_z, recvCount_z ,recvBuf_z, Blobs.data); - UnpackBlobsData(recvList_Z, recvCount_Z ,recvBuf_Z, Blobs.data); - UnpackBlobsData(recvList_xy, recvCount_xy ,recvBuf_xy, Blobs.data); - UnpackBlobsData(recvList_Xy, recvCount_Xy ,recvBuf_Xy, Blobs.data); - UnpackBlobsData(recvList_xY, recvCount_xY ,recvBuf_xY, Blobs.data); - UnpackBlobsData(recvList_XY, recvCount_XY ,recvBuf_XY, Blobs.data); - UnpackBlobsData(recvList_xz, recvCount_xz ,recvBuf_xz, Blobs.data); - UnpackBlobsData(recvList_Xz, recvCount_Xz ,recvBuf_Xz, Blobs.data); - UnpackBlobsData(recvList_xZ, recvCount_xZ ,recvBuf_xZ, Blobs.data); - UnpackBlobsData(recvList_XZ, recvCount_XZ ,recvBuf_XZ, Blobs.data); - UnpackBlobsData(recvList_yz, recvCount_yz ,recvBuf_yz, Blobs.data); - UnpackBlobsData(recvList_Yz, recvCount_Yz ,recvBuf_Yz, Blobs.data); - UnpackBlobsData(recvList_yZ, recvCount_yZ ,recvBuf_yZ, Blobs.data); - UnpackBlobsData(recvList_YZ, recvCount_YZ ,recvBuf_YZ, Blobs.data); + UnpackBlobData(recvList_x, recvCount_x ,recvBuf_x, Blobs.data); + UnpackBlobData(recvList_X, recvCount_X ,recvBuf_X, Blobs.data); + UnpackBlobData(recvList_y, recvCount_y ,recvBuf_y, Blobs.data); + UnpackBlobData(recvList_Y, recvCount_Y ,recvBuf_Y, Blobs.data); + UnpackBlobData(recvList_z, recvCount_z ,recvBuf_z, Blobs.data); + UnpackBlobData(recvList_Z, recvCount_Z ,recvBuf_Z, Blobs.data); + UnpackBlobData(recvList_xy, recvCount_xy ,recvBuf_xy, Blobs.data); + UnpackBlobData(recvList_Xy, recvCount_Xy ,recvBuf_Xy, Blobs.data); + UnpackBlobData(recvList_xY, recvCount_xY ,recvBuf_xY, Blobs.data); + UnpackBlobData(recvList_XY, recvCount_XY ,recvBuf_XY, Blobs.data); + UnpackBlobData(recvList_xz, recvCount_xz ,recvBuf_xz, Blobs.data); + UnpackBlobData(recvList_Xz, recvCount_Xz ,recvBuf_Xz, Blobs.data); + UnpackBlobData(recvList_xZ, recvCount_xZ ,recvBuf_xZ, Blobs.data); + UnpackBlobData(recvList_XZ, recvCount_XZ ,recvBuf_XZ, Blobs.data); + UnpackBlobData(recvList_yz, recvCount_yz ,recvBuf_yz, Blobs.data); + UnpackBlobData(recvList_Yz, recvCount_Yz ,recvBuf_Yz, Blobs.data); + UnpackBlobData(recvList_yZ, recvCount_yZ ,recvBuf_yZ, Blobs.data); + UnpackBlobData(recvList_YZ, recvCount_YZ ,recvBuf_YZ, Blobs.data); //...................................................................................... } From 0bfbeb9b9d34f3997af4450427d8e5aef9890813 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 22:30:49 -0400 Subject: [PATCH 40/43] Updated Domain::CommInit() --- common/Domain.h | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index 70790799..444fa200 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -11,6 +11,7 @@ #include #include "common/Utilities.h" +int MAX_BLOB_COUNT=500; using namespace std; @@ -25,6 +26,7 @@ struct Domain{ N = Nx*Ny*Nz; id = new char [N]; Blobs.New(Nx,Ny,Nz); + BlobGraph.New(18,MAX_BLOB_COUNT); } // Basic domain information @@ -33,7 +35,6 @@ struct Domain{ int nprocx,nprocy,nprocz; double Lx,Ly,Lz; int rank; - MPI_Comm Communicator; //********************************** // MPI ranks for all 18 neighbors @@ -76,12 +77,23 @@ struct Domain{ char *id; // Blob information IntArray Blobs; + IntArray BlobGraph; void InitializeRanks(); void CommInit(MPI_Comm comm); void BlobComm(MPI_Comm comm); + + void getBlobConnections(); +} private: + int getRankForBlock( int i, int j, int k ) + { + int i2 = (i+nprocx)%nprocx; + int j2 = (j+nprocy)%nprocy; + int k2 = (k+nprocz)%nprocz; + return i2 + j2*nprocx + k2*nprocx*nprocy; + } void PackBlobData(int *list, int count, int *sendbuf, int *data){ // Fill in the phase ID values from neighboring processors // This packs up the values that need to be sent from one processor to another @@ -100,13 +112,6 @@ private: data[n] = recvbuf[idx]; } } - int getRankForBlock( int i, int j, int k ) - { - int i2 = (i+nprocx)%nprocx; - int j2 = (j+nprocy)%nprocy; - int k2 = (k+nprocz)%nprocz; - return i2 + j2*nprocx + k2*nprocx*nprocy; - } }; @@ -483,6 +488,13 @@ void Domain::BlobComm(MPI_Comm Communicator){ //...................................................................................... } +Domain::getBlobConnections(){ + + BlobGraph(0,nblob) = rank of the connecting blob; + BlobGraph(1,nblob) = ID of the connecting blob; + +} + inline void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad) { // Read in the full sphere pack From 1aa533f8584691a81a7b494d479dda1f68954a3c Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 22:31:12 -0400 Subject: [PATCH 41/43] Updated Domain::CommInit() --- common/Domain.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/common/Domain.h b/common/Domain.h index 444fa200..3a180cf5 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -490,8 +490,8 @@ void Domain::BlobComm(MPI_Comm Communicator){ Domain::getBlobConnections(){ - BlobGraph(0,nblob) = rank of the connecting blob; - BlobGraph(1,nblob) = ID of the connecting blob; +// BlobGraph(0,nblob) = rank of the connecting blob; +// BlobGraph(1,nblob) = ID of the connecting blob; } From c790232151c1563a19b5868e46df365e67116f21 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 22:31:50 -0400 Subject: [PATCH 42/43] Updated Domain::CommInit() --- common/Domain.h | 1 - 1 file changed, 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index 3a180cf5..b0db6299 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -84,7 +84,6 @@ struct Domain{ void BlobComm(MPI_Comm comm); void getBlobConnections(); -} private: int getRankForBlock( int i, int j, int k ) From 414d18f52737f5701dcf467be853b2b789abd487 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 25 Jun 2014 22:32:24 -0400 Subject: [PATCH 43/43] Updated Domain::CommInit() --- common/Domain.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/Domain.h b/common/Domain.h index b0db6299..8f88a3c1 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -487,7 +487,7 @@ void Domain::BlobComm(MPI_Comm Communicator){ //...................................................................................... } -Domain::getBlobConnections(){ +void Domain::getBlobConnections(){ // BlobGraph(0,nblob) = rank of the connecting blob; // BlobGraph(1,nblob) = ID of the connecting blob;