diff --git a/common/Domain.h b/common/Domain.h index 5ed10d21..8f88a3c1 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -8,12 +8,492 @@ #include #include // std::exception #include -#include +#include #include "common/Utilities.h" +int MAX_BLOB_COUNT=500; 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+2; Ny = ny+2; Nz = nz+2; + 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); + BlobGraph.New(18,MAX_BLOB_COUNT); + } + + // 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; + //********************************** + //...................................................................................... + // 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; + //...................................................................................... + 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 *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; + //...................................................................................... + + // Solid indicator function + 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 + int idx,n; + for (idx=0; idx +#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; + +struct BlobInfo{ + + BlobInfo(int Nx, int Ny,int Nz){ + ID.New(Nx,Ny,Nz); + } + 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); + +}; + +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; idxvf|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"); + } + + 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, + 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 ){ + 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,Dm.Blobs,Phase,SignDist,vF,vS,i,j,k,temp); + nblobs++; + } + } + } + + } + } + } + printf("Number of blobs %i, Rank = %i \n",nblobs,rank); + MPI_Barrier(MPI_COMM_WORLD); + + /// Compute volume averages + for (k=kstart; k 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