From 42338802c523c65d2ee6746dd0ff9bb9a840c302 Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 27 Jan 2014 11:43:24 -0500 Subject: [PATCH] Updated CPU branch --- cpu/Color.cpp | 833 ++++++++-- cpu/Color.h | 41 +- cpu/D3Q19.cpp | 46 +- cpu/D3Q19.h | 8 +- cpu/D3Q7.cpp | 10 +- cpu/D3Q7.h | 8 +- cpu/Extras.cpp | 28 + cpu/Extras.h | 7 + cpu/lb2_Color_wia_mpi.cpp | 2447 +++++++++++++++++++++++++++++ cpu/lb2_Color_wia_mpi_bubble.cpp | 2492 ++++++++++++++++++++++++++++++ include/Array.h | 12 +- include/Domain.h | 15 + include/pmmc.h | 75 +- pmmc/Analysis.cpp | 127 +- tests/TestContactAngle.cpp | 2 - tests/TestInterfaceSpeed.cpp | 14 +- 16 files changed, 5921 insertions(+), 244 deletions(-) create mode 100644 cpu/Extras.cpp create mode 100644 cpu/Extras.h create mode 100644 cpu/lb2_Color_wia_mpi.cpp create mode 100644 cpu/lb2_Color_wia_mpi_bubble.cpp diff --git a/cpu/Color.cpp b/cpu/Color.cpp index 67228f6f..63010746 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -1,10 +1,18 @@ #include -extern void InitDenColor(char *ID, double *Den, double *Phi, double das, double dbs, int N) +extern "C" void dvc_InitDenColor(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz, int S) { - int n; + int i,j,k,n,N; + + N = Nx*Ny*Nz; + for (n=0; n 0){ - f_even[n] = 0.3333333333333333; - f_odd[n] = 0.055555555555555555; //double(100*n)+1.f; - f_even[N+n] = 0.055555555555555555; //double(100*n)+2.f; - f_odd[N+n] = 0.055555555555555555; //double(100*n)+3.f; - f_even[2*N+n] = 0.055555555555555555; //double(100*n)+4.f; - f_odd[2*N+n] = 0.055555555555555555; //double(100*n)+5.f; - f_even[3*N+n] = 0.055555555555555555; //double(100*n)+6.f; - f_odd[3*N+n] = 0.0277777777777778; //double(100*n)+7.f; - f_even[4*N+n] = 0.0277777777777778; //double(100*n)+8.f; - f_odd[4*N+n] = 0.0277777777777778; //double(100*n)+9.f; - f_even[5*N+n] = 0.0277777777777778; //double(100*n)+10.f; - f_odd[5*N+n] = 0.0277777777777778; //double(100*n)+11.f; - f_even[6*N+n] = 0.0277777777777778; //double(100*n)+12.f; - f_odd[6*N+n] = 0.0277777777777778; //double(100*n)+13.f; - f_even[7*N+n] = 0.0277777777777778; //double(100*n)+14.f; - f_odd[7*N+n] = 0.0277777777777778; //double(100*n)+15.f; - f_even[8*N+n] = 0.0277777777777778; //double(100*n)+16.f; - f_odd[8*N+n] = 0.0277777777777778; //double(100*n)+17.f; - f_even[9*N+n] = 0.0277777777777778; //double(100*n)+18.f; + + //.......Back out the 3-D indices for node n.............. + k = n/(Nx*Ny); + j = (n-Nx*Ny*k)/Nx; + i = n-Nx*Ny*k-Nx*j; + + if ( ID[n] == 1){ + Den[2*n] = 1.0; + Den[2*n+1] = 0.0; + Phi[n] = 1.0; + } + else if ( ID[n] == 2){ + Den[2*n] = 0.0; + Den[2*n+1] = 1.0; + Phi[n] = -1.0; } else{ - for(int q=0; q<9; q++){ - f_even[q*N+n] = -1.0; - f_odd[q*N+n] = -1.0; - } - f_even[9*N+n] = -1.0; + Den[2*n] = das; + Den[2*n+1] = dbs; + Phi[n] = (das-dbs)/(das+dbs); + d = fabs(Distance[n]); + Phi[n] = (2.f*(exp(-2.f*beta*(d+xp)))/(1.f+exp(-2.f*beta*(d+xp))) - 1.f); + } + if (i == 0 || j == 0 || k == 0 || i == Nx-1 || j == Ny-1 || k == Nz-1){ + Den[2*n] = 0.0; + Den[2*n+1] = 0.0; } } } -extern void Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz) + +extern "C" void dvc_Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz, int S) { int n,N; // distributions double f1,f2,f3,f4,f5,f6,f7,f8,f9; double f10,f11,f12,f13,f14,f15,f16,f17,f18; double vx,vy,vz; - + N = Nx*Ny*Nz; - + for (n=0; n 0){ //........................................................................ // Registers to store the distributions @@ -103,26 +117,70 @@ extern void Compute_VELOCITY(char *ID, double *disteven, double *distodd, double vel[N+n] = vy; vel[2*N+n] = vz; //........................................................................ - + + } + } +} + +extern "C" void dvc_ComputePressureD3Q19(char *ID, double *disteven, double *distodd, double *Pressure, + int Nx, int Ny, int Nz, int S) +{ + int n,N; + // distributions + double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9; + double f10,f11,f12,f13,f14,f15,f16,f17,f18; + + N = Nx*Ny*Nz; + + for (n=0; n 0){ + //........................................................................ + // Registers to store the distributions + //........................................................................ + f0 = disteven[n]; + f2 = disteven[N+n]; + f4 = disteven[2*N+n]; + f6 = disteven[3*N+n]; + f8 = disteven[4*N+n]; + f10 = disteven[5*N+n]; + f12 = disteven[6*N+n]; + f14 = disteven[7*N+n]; + f16 = disteven[8*N+n]; + f18 = disteven[9*N+n]; + //........................................................................ + f1 = distodd[n]; + f3 = distodd[1*N+n]; + f5 = distodd[2*N+n]; + f7 = distodd[3*N+n]; + f9 = distodd[4*N+n]; + f11 = distodd[5*N+n]; + f13 = distodd[6*N+n]; + f15 = distodd[7*N+n]; + f17 = distodd[8*N+n]; + //.................Compute the velocity................................... + Pressure[n] = 0.3333333333333333*(f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+ + f9+f12+f11+f14+f13+f16+f15+f18+f17); + } } } //************************************************************************* //************************************************************************* -extern void PressureBC_inlet(double *disteven, double *distodd, double din, - int Nx, int Ny, int Nz) +extern "C" void dvc_PressureBC_inlet(double *disteven, double *distodd, double din, + int Nx, int Ny, int Nz, int S) { int n,N; // distributions double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9; double f10,f11,f12,f13,f14,f15,f16,f17,f18; double uz; - + N = Nx*Ny*Nz; - + for (n=0; n 0){ - + // Retrieve the color gradient - nx = ColorGrad[3*n]; - ny = ColorGrad[3*n+1]; - nz = ColorGrad[3*n+2]; + nx = ColorGrad[n]; + ny = ColorGrad[N+n]; + nz = ColorGrad[2*N+n]; //...........Normalize the Color Gradient................................. C = sqrt(nx*nx+ny*ny+nz*nz); nx = nx/C; @@ -474,57 +534,464 @@ extern void ColorCollide( char *ID, double *disteven, double *distodd, double *C //.................inverse transformation...................................................... f0 = 0.05263157894736842*rho-0.012531328320802*m1+0.04761904761904762*m2; f1 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 - +0.1*(jx-m4)+0.0555555555555555555555555*(m9-m10); + +0.1*(jx-m4)+0.0555555555555555555555555*(m9-m10); f2 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 - +0.1*(m4-jx)+0.0555555555555555555555555*(m9-m10); + +0.1*(m4-jx)+0.0555555555555555555555555*(m9-m10); f3 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 - +0.1*(jy-m6)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m11-m12); + +0.1*(jy-m6)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m11-m12); f4 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 - +0.1*(m6-jy)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m11-m12); + +0.1*(m6-jy)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m11-m12); f5 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 - +0.1*(jz-m8)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m12-m11); + +0.1*(jz-m8)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m12-m11); f6 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 - +0.1*(m8-jz)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m12-m11); + +0.1*(m8-jz)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m12-m11); f7 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2+0.1*(jx+jy)+0.025*(m4+m6) - +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 - +0.04166666666666666*m12+0.25*m13+0.125*(m16-m17); + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12+0.25*m13+0.125*(m16-m17); f8 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2-0.1*(jx+jy)-0.025*(m4+m6) - +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 - +0.04166666666666666*m12+0.25*m13+0.125*(m17-m16); + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12+0.25*m13+0.125*(m17-m16); f9 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2+0.1*(jx-jy)+0.025*(m4-m6) - +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 - +0.04166666666666666*m12-0.25*m13+0.125*(m16+m17); + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12-0.25*m13+0.125*(m16+m17); f10 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2+0.1*(jy-jx)+0.025*(m6-m4) - +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 - +0.04166666666666666*m12-0.25*m13-0.125*(m16+m17); + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12-0.25*m13-0.125*(m16+m17); f11 = 0.05263157894736842*rho+0.003341687552213868*m1 - +0.003968253968253968*m2+0.1*(jx+jz)+0.025*(m4+m8) - +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 - -0.04166666666666666*m12+0.25*m15+0.125*(m18-m16); + +0.003968253968253968*m2+0.1*(jx+jz)+0.025*(m4+m8) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12+0.25*m15+0.125*(m18-m16); f12 = 0.05263157894736842*rho+0.003341687552213868*m1 - +0.003968253968253968*m2-0.1*(jx+jz)-0.025*(m4+m8) - +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 - -0.04166666666666666*m12+0.25*m15+0.125*(m16-m18); + +0.003968253968253968*m2-0.1*(jx+jz)-0.025*(m4+m8) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12+0.25*m15+0.125*(m16-m18); f13 = 0.05263157894736842*rho+0.003341687552213868*m1 - +0.003968253968253968*m2+0.1*(jx-jz)+0.025*(m4-m8) - +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 - -0.04166666666666666*m12-0.25*m15-0.125*(m16+m18); + +0.003968253968253968*m2+0.1*(jx-jz)+0.025*(m4-m8) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12-0.25*m15-0.125*(m16+m18); f14 = 0.05263157894736842*rho+0.003341687552213868*m1 - +0.003968253968253968*m2+0.1*(jz-jx)+0.025*(m8-m4) - +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 - -0.04166666666666666*m12-0.25*m15+0.125*(m16+m18); + +0.003968253968253968*m2+0.1*(jz-jx)+0.025*(m8-m4) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12-0.25*m15+0.125*(m16+m18); f15 = 0.05263157894736842*rho+0.003341687552213868*m1 - +0.003968253968253968*m2+0.1*(jy+jz)+0.025*(m6+m8) - -0.0555555555555555555555555*m9-0.02777777777777778*m10+0.25*m14+0.125*(m17-m18); + +0.003968253968253968*m2+0.1*(jy+jz)+0.025*(m6+m8) + -0.0555555555555555555555555*m9-0.02777777777777778*m10+0.25*m14+0.125*(m17-m18); f16 = 0.05263157894736842*rho+0.003341687552213868*m1 - +0.003968253968253968*m2-0.1*(jy+jz)-0.025*(m6+m8) - -0.0555555555555555555555555*m9-0.02777777777777778*m10+0.25*m14+0.125*(m18-m17); + +0.003968253968253968*m2-0.1*(jy+jz)-0.025*(m6+m8) + -0.0555555555555555555555555*m9-0.02777777777777778*m10+0.25*m14+0.125*(m18-m17); f17 = 0.05263157894736842*rho+0.003341687552213868*m1 - +0.003968253968253968*m2+0.1*(jy-jz)+0.025*(m6-m8) - -0.0555555555555555555555555*m9-0.02777777777777778*m10-0.25*m14+0.125*(m17+m18); + +0.003968253968253968*m2+0.1*(jy-jz)+0.025*(m6-m8) + -0.0555555555555555555555555*m9-0.02777777777777778*m10-0.25*m14+0.125*(m17+m18); f18 = 0.05263157894736842*rho+0.003341687552213868*m1 - +0.003968253968253968*m2+0.1*(jz-jy)+0.025*(m8-m6) - -0.0555555555555555555555555*m9-0.02777777777777778*m10-0.25*m14-0.125*(m17+m18); + +0.003968253968253968*m2+0.1*(jz-jy)+0.025*(m8-m6) + -0.0555555555555555555555555*m9-0.02777777777777778*m10-0.25*m14-0.125*(m17+m18); + //....................................................................................................... + // incorporate external force + f1 += 0.16666666*Fx; + f2 -= 0.16666666*Fx; + f3 += 0.16666666*Fy; + f4 -= 0.16666666*Fy; + f5 += 0.16666666*Fz; + f6 -= 0.16666666*Fz; + f7 += 0.08333333333*(Fx+Fy); + f8 -= 0.08333333333*(Fx+Fy); + f9 += 0.08333333333*(Fx-Fy); + f10 -= 0.08333333333*(Fx-Fy); + f11 += 0.08333333333*(Fx+Fz); + f12 -= 0.08333333333*(Fx+Fz); + f13 += 0.08333333333*(Fx-Fz); + f14 -= 0.08333333333*(Fx-Fz); + f15 += 0.08333333333*(Fy+Fz); + f16 -= 0.08333333333*(Fy+Fz); + f17 += 0.08333333333*(Fy-Fz); + f18 -= 0.08333333333*(Fy-Fz); + //*********** WRITE UPDATED VALUES TO MEMORY ****************** + // Write the updated distributions + //....EVEN..................................... + disteven[n] = f0; + disteven[N+n] = f2; + disteven[2*N+n] = f4; + disteven[3*N+n] = f6; + disteven[4*N+n] = f8; + disteven[5*N+n] = f10; + disteven[6*N+n] = f12; + disteven[7*N+n] = f14; + disteven[8*N+n] = f16; + disteven[9*N+n] = f18; + //....ODD...................................... + distodd[n] = f1; + distodd[N+n] = f3; + distodd[2*N+n] = f5; + distodd[3*N+n] = f7; + distodd[4*N+n] = f9; + distodd[5*N+n] = f11; + distodd[6*N+n] = f13; + distodd[7*N+n] = f15; + distodd[8*N+n] = f17; + + //...Store the Velocity.......................... + Velocity[3*n] = jx; + Velocity[3*n+1] = jy; + Velocity[3*n+2] = jz; + //...Store the Color Gradient.................... + // ColorGrad[3*n] = nx*C; + // ColorGrad[3*n+1] = ny*C; + // ColorGrad[3*n+2] = nz*C; + //............................................... + //*************************************************************** + } // check if n is in the solid + } // loop over n +} + +extern "C" void dvc_ColorCollideOpt( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad, + double *Velocity, int Nx, int Ny, int Nz, int S,double rlx_setA, double rlx_setB, + double alpha, double beta, double Fx, double Fy, double Fz) +{ + + int i,j,k,n,nn,N; + // distributions + double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9; + double f10,f11,f12,f13,f14,f15,f16,f17,f18; + + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + // additional variables needed for computations + double rho,jx,jy,jz,C,nx,ny,nz; + + N = Nx*Ny*Nz; + char id; + + for (n=0; n 0){ + + //.......Back out the 3-D indices for node n.............. + k = n/(Nx*Ny); + j = (n-Nx*Ny*k)/Nx; + i = n-Nx*Ny*k-Nx*j; + //........................................................................ + //........Get 1-D index for this thread.................... + // n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x; + //........................................................................ + // COMPUTE THE COLOR GRADIENT + //........................................................................ + //.................Read Phase Indicator Values............................ + //........................................................................ + nn = n-1; // neighbor index (get convention) + if (i-1<0) nn += Nx; // periodic BC along the x-boundary + f1 = phi[nn]; // get neighbor for phi - 1 + //........................................................................ + nn = n+1; // neighbor index (get convention) + if (!(i+1 0){ + for (n=0; n 0 && n>> (ID, Den, Phi, das, dbs, Nx, Ny, Nz, S); +} +//************************************************************************* +extern "C" void dvc_ComputeColorGradient(int nBlocks, int nthreads, int S, + char *ID, double *Phi, double *ColorGrad, int Nx, int Ny, int Nz) +{ + ComputeColorGradient<<>>(ID, Phi, ColorGrad, Nx, Ny, Nz, S); +} +//************************************************************************* +extern "C" void dvc_ColorCollide(int nBlocks, int nthreads, int S, + char *ID, double *f_even, double *f_odd, double *ColorGrad, double *Velocity, + double rlxA, double rlxB,double alpha, double beta, double Fx, double Fy, double Fz, + int Nx, int Ny, int Nz, bool pBC) +{ + ColorCollide<<>>(ID, f_even, f_odd, ColorGrad, Velocity, Nx, Ny, Nz, S, + rlxA, rlxB, alpha, beta, Fx, Fy, Fz, pBC); +} +//************************************************************************* +extern "C" void dvc_ColorCollideOpt(int nBlocks, int nthreads, int S, + char *ID, double *f_even, double *f_odd, double *Phi, double *ColorGrad, + double *Velocity, int Nx, int Ny, int Nz,double rlxA, double rlxB, + double alpha, double beta, double Fx, double Fy, double Fz) +{ + ColorCollideOpt<<>>(ID, f_even, f_odd, Phi, ColorGrad, Velocity, Nx, Ny, Nz, S, + rlxA, rlxB, alpha, beta, Fx, Fy, Fz); +// bool pBC = false; +// ColorCollide<<>>(ID, f_even, f_odd, ColorGrad, Velocity, Nx, Ny, Nz, S, +// rlxA, rlxB, alpha, beta, Fx, Fy, Fz, pBC); +} + +//************************************************************************* +extern "C" void dvc_DensityStreamD3Q7(int nBlocks, int nthreads, int S, + char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity, + double beta, int Nx, int Ny, int Nz, bool pBC) +{ + DensityStreamD3Q7<<>>(ID,Den,Copy,Phi,ColorGrad,Velocity,beta,Nx,Ny,Nz,pBC,S); +} +//************************************************************************* +extern "C" void dvc_ComputePhi(int nBlocks, int nthreads, int S, + char *ID, double *Phi, double *Copy, double *Den, int N) +{ + ComputePhi<<>>(ID,Phi,Copy,Den,N,S); +} +//************************************************************************* +extern "C" void dvc_ComputePressure(int nBlocks, int nthreads, int S, + char *ID, double *disteven, double *distodd, + double *Pressure, int Nx, int Ny, int Nz) +{ + + ComputePressureD3Q19<<>>(ID,disteven,distodd,Pressure,Nx,Ny,Nz,S); + +} +*/ diff --git a/cpu/Color.h b/cpu/Color.h index a63128c0..29fc5561 100644 --- a/cpu/Color.h +++ b/cpu/Color.h @@ -1,21 +1,20 @@ -extern void InitDenColor(char *ID, double *Den, double *Phi, double das, double dbs, int N); -extern void InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz); - -extern void Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz); - -//************************************************************************* -//************************************************************************* -extern void PressureBC_inlet(double *disteven, double *distodd, double din, - int Nx, int Ny, int Nz); -extern void PressureBC_outlet(double *disteven, double *distodd, double dout, - int Nx, int Ny, int Nz, int S, int outlet); -//************************************************************************* -extern void ComputeColorGradient(char *ID, double *phi, double *ColorGrad, int Nx, int Ny, int Nz); -//************************************************************************* -extern void ColorCollide( char *ID, double *disteven, double *distodd, double *ColorGrad, - double *Velocity, int Nx, int Ny, int Nz, double rlx_setA, double rlx_setB, - double alpha, double beta, double Fx, double Fy, double Fz, bool pBC); -//************************************************************************* -extern void DensityStreamD3Q7(char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity, - double beta, int Nx, int Ny, int Nz, bool pBC); -extern void ComputePhi(char *ID, double *Phi, double *Copy, double *Den, int N); +extern "C" void dvc_InitDenColor(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz, int S); +extern "C" void dvc_InitDenColorDistance(char *ID, double *Den, double *Phi, double *Distance, + double das, double dbs, double beta, double xp, int Nx, int Ny, int Nz, int S); +extern "C" void dvc_Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz, int S); +extern "C" void dvc_ComputePressureD3Q19(char *ID, double *disteven, double *distodd, double *Pressure, + int Nx, int Ny, int Nz, int S); +extern "C" void dvc_PressureBC_inlet(double *disteven, double *distodd, double din, + int Nx, int Ny, int Nz, int S); +extern "C" void dvc_PressureBC_outlet(double *disteven, double *distodd, double dout, + int Nx, int Ny, int Nz, int S, int outlet); +extern "C" void dvc_ComputeColorGradient(char *ID, double *phi, double *ColorGrad, int Nx, int Ny, int Nz, int S); +extern "C" void dvc_ColorCollide( char *ID, double *disteven, double *distodd, double *ColorGrad, + double *Velocity, int Nx, int Ny, int Nz, int S,double rlx_setA, double rlx_setB, + double alpha, double beta, double Fx, double Fy, double Fz, bool pBC); +extern "C" void dvc_ColorCollideOpt( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad, + double *Velocity, int Nx, int Ny, int Nz, int S,double rlx_setA, double rlx_setB, + double alpha, double beta, double Fx, double Fy, double Fz); +extern "C" void dvc_DensityStreamD3Q7(char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity, + double beta, int Nx, int Ny, int Nz, bool pBC, int S); +extern "C" void dvc_ComputePhi(char *ID, double *Phi, double *Copy, double *Den, int N, int S); diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp index b1872135..c4267fbc 100644 --- a/cpu/D3Q19.cpp +++ b/cpu/D3Q19.cpp @@ -1,4 +1,4 @@ -extern void PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N){ +extern "C" void dvc_PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N){ //.................................................................................... // Pack distribution q into the send buffer for the listed lattice sites // dist may be even or odd distributions stored by stream layout @@ -10,9 +10,7 @@ extern void PackDist(int q, int *list, int start, int count, double *sendbuf, do } } - - -extern void MapRecvDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count, +extern "C" void dvc_UnpackDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count, double *recvbuf, double *dist, int Nx, int Ny, int Nz){ //.................................................................................... // Unack distribution from the recv buffer @@ -54,8 +52,46 @@ extern void MapRecvDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, } } +extern "C" void dvc_InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz, int S) +{ + int n,N; + N = Nx*Ny*Nz; + + for (n=0; n 0){ + f_even[n] = 0.3333333333333333; + f_odd[n] = 0.055555555555555555; //double(100*n)+1.f; + f_even[N+n] = 0.055555555555555555; //double(100*n)+2.f; + f_odd[N+n] = 0.055555555555555555; //double(100*n)+3.f; + f_even[2*N+n] = 0.055555555555555555; //double(100*n)+4.f; + f_odd[2*N+n] = 0.055555555555555555; //double(100*n)+5.f; + f_even[3*N+n] = 0.055555555555555555; //double(100*n)+6.f; + f_odd[3*N+n] = 0.0277777777777778; //double(100*n)+7.f; + f_even[4*N+n] = 0.0277777777777778; //double(100*n)+8.f; + f_odd[4*N+n] = 0.0277777777777778; //double(100*n)+9.f; + f_even[5*N+n] = 0.0277777777777778; //double(100*n)+10.f; + f_odd[5*N+n] = 0.0277777777777778; //double(100*n)+11.f; + f_even[6*N+n] = 0.0277777777777778; //double(100*n)+12.f; + f_odd[6*N+n] = 0.0277777777777778; //double(100*n)+13.f; + f_even[7*N+n] = 0.0277777777777778; //double(100*n)+14.f; + f_odd[7*N+n] = 0.0277777777777778; //double(100*n)+15.f; + f_even[8*N+n] = 0.0277777777777778; //double(100*n)+16.f; + f_odd[8*N+n] = 0.0277777777777778; //double(100*n)+17.f; + f_even[9*N+n] = 0.0277777777777778; //double(100*n)+18.f; + } + else{ + for(int q=0; q<9; q++){ + f_even[q*N+n] = -1.0; + f_odd[q*N+n] = -1.0; + } + f_even[9*N+n] = -1.0; + } + } +} + //************************************************************************* -extern void SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz) +extern "C" void dvc_SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz, int S) { int n,nn,N; // distributions diff --git a/cpu/D3Q19.h b/cpu/D3Q19.h index 7cb75617..3a888553 100644 --- a/cpu/D3Q19.h +++ b/cpu/D3Q19.h @@ -1,6 +1,6 @@ - -extern void PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N); -extern void MapRecvDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count, +extern "C" void dvc_PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N); +extern "C" void dvc_UnpackDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count, double *recvbuf, double *dist, int Nx, int Ny, int Nz); //************************************************************************* -extern void SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz); +extern "C" void dvc_InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz, int S); +extern "C" void dvc_SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz, int S); diff --git a/cpu/D3Q7.cpp b/cpu/D3Q7.cpp index f9cba98c..8718fbaa 100644 --- a/cpu/D3Q7.cpp +++ b/cpu/D3Q7.cpp @@ -1,6 +1,6 @@ -// GPU Functions for D3Q7 Lattice Boltzmann Methods +// CPU Functions for D3Q7 Lattice Boltzmann Methods -extern void PackValues(int *list, int count, double *sendbuf, double *Data, int N){ +extern "C" void dvc_PackValues(int *list, int count, double *sendbuf, double *Data, int N){ //.................................................................................... // Pack distribution q into the send buffer for the listed lattice sites // dist may be even or odd distributions stored by stream layout @@ -11,7 +11,7 @@ extern void PackValues(int *list, int count, double *sendbuf, double *Data, int sendbuf[idx] = Data[n]; } } -extern void UnpackValues(int *list, int count, double *recvbuf, double *Data, int N){ +extern "C" void dvc_UnpackValues(int *list, int count, double *recvbuf, double *Data, int N){ //.................................................................................... // Pack distribution q into the send buffer for the listed lattice sites // dist may be even or odd distributions stored by stream layout @@ -23,7 +23,7 @@ extern void UnpackValues(int *list, int count, double *recvbuf, double *Data, in } } -extern void PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N){ +extern "C" void dvc_PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N){ //.................................................................................... // Pack distribution into the send buffer for the listed lattice sites //.................................................................................... @@ -38,7 +38,7 @@ extern void PackDenD3Q7(int *list, int count, double *sendbuf, int number, doubl } -extern void UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N){ +extern "C" void dvc_UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N){ //.................................................................................... // Unack distribution from the recv buffer // Sum to the existing density value diff --git a/cpu/D3Q7.h b/cpu/D3Q7.h index 83309bf8..bce793cc 100644 --- a/cpu/D3Q7.h +++ b/cpu/D3Q7.h @@ -1,9 +1,9 @@ // CPU Functions for D3Q7 Lattice Boltzmann Methods -extern void PackValues(int *list, int count, double *sendbuf, double *Data, int N); +extern "C" void dvc_PackValues(int *list, int count, double *sendbuf, double *Data, int N); -extern void UnpackValues(int *list, int count, double *recvbuf, double *Data, int N); +extern "C" void dvc_UnpackValues(int *list, int count, double *recvbuf, double *Data, int N); -extern void PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N); +extern "C" void dvc_PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N); -extern void UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N); +extern "C" void dvc_UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N); diff --git a/cpu/Extras.cpp b/cpu/Extras.cpp new file mode 100644 index 00000000..1d925495 --- /dev/null +++ b/cpu/Extras.cpp @@ -0,0 +1,28 @@ +// Basic cuda functions callable from C/C++ code +#include +#include +#include + +extern "C" void dvc_AllocateDeviceMemory(void** address, size_t size){ + //cudaMalloc(address,size); + (*address) = malloc(size); + + if (*address==NULL){ + printf("Memory allocation failed! \n"); + } +} + +extern "C" void dvc_CopyToDevice(void* dest, void* source, size_t size){ +// cudaMemcpy(dest,source,size,cudaMemcpyHostToDevice); + memcpy(dest, source, size); +} + + +extern "C" void dvc_CopyToHost(void* dest, void* source, size_t size){ +// cudaMemcpy(dest,source,size,cudaMemcpyDeviceToHost); + memcpy(dest, source, size); +} + +extern "C" void dvc_Barrier(){ +// cudaDeviceSynchronize(); +} diff --git a/cpu/Extras.h b/cpu/Extras.h new file mode 100644 index 00000000..c2f3363f --- /dev/null +++ b/cpu/Extras.h @@ -0,0 +1,7 @@ +extern "C" void dvc_AllocateDeviceMemory(void** address, size_t size); + +extern "C" void dvc_CopyToDevice(void* dest, void* source, size_t size); + +extern "C" void dvc_CopyToHost(void* dest, void* source, size_t size); + +extern "C" void dvc_Barrier(); diff --git a/cpu/lb2_Color_wia_mpi.cpp b/cpu/lb2_Color_wia_mpi.cpp new file mode 100644 index 00000000..986b1042 --- /dev/null +++ b/cpu/lb2_Color_wia_mpi.cpp @@ -0,0 +1,2447 @@ +#include +#include +#include +#include +#include + +#include "pmmc.h" +#include "Domain.h" +#include "Extras.h" +#include "D3Q19.h" +#include "D3Q7.h" +#include "Color.h" + +using namespace std; + +//************************************************************************* +// Implementation of Two-Phase Immiscible LBM using CUDA +//************************************************************************* +inline void PackID(int *list, int count, char *sendbuf, char *ID){ + // Fill in the phase ID values from neighboring processors + // This packs up the values that need to be sent from one processor to another + int idx,n; + + for (idx=0; idx> 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 >> xIntPos; // Contact angle parameter +// 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; // error interval + input >> tol; // error tolerance + //............................................................. + das = 0.1; dbs = 0.9; // hard coded for density initialization + // should be OK to remove these parameters + // they should have no impact with the + // current boundary condition + //....................................................................... + // Reading the domain information file + //....................................................................... + ifstream domain("Domain.in"); + 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(&xIntPos,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(&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); + // ************************************************************** + // ************************************************************** + double Ps = -(das-dbs)/(das+dbs); + double rlxA = 1.f/tau; + double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA); + + if (nprocs != nprocx*nprocy*nprocz){ + printf("Fatal error in processor number! \n"); + printf("nprocx = %i \n",nprocx); + printf("nprocy = %i \n",nprocy); + printf("nprocz = %i \n",nprocz); + } + + 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("phi_s = %f \n", Ps); + printf("gamma_{wn} = %f \n", 6.01603*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"); + } + + MPI_Barrier(MPI_COMM_WORLD); + kproc = rank/(nprocx*nprocy); + jproc = (rank-nprocx*nprocy*kproc)/nprocx; + iproc = rank-nprocx*nprocy*kproc-nprocz*jproc; + + //.......................................... + // set up the neighbor ranks + //.......................................... + i=iproc; j=jproc; k =kproc; + i+=1; + j+=0; + k+=0; + if (i<0) i+=nprocx; + if (j<0) j+=nprocy; + if (k<0) k+=nprocz; + if (!(i 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] = 1; + sum++; + } + } + } + } + sum_local = 1.0*sum; + MPI_Allreduce(&sum_local,&porosity,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + porosity = porosity*iVol_global; + if (rank==0) printf("Media porosity = %f \n",porosity); + + // Generate the residual NWP + if (rank==0) printf("Initializing with NWP saturation = %f \n",wp_saturation); + GenerateResidual(id,Nx,Ny,Nz,wp_saturation); +#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 InterfaceSpeed(20); + DoubleArray NormalVector(60); + + // IntArray store; + + int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0; + int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0; + +// double s,s1,s2,s3; // Triangle sides (lengths) + Point A,B,C,P; +// double area; + + // Initialize arrays for local solid surface + DTMutableList local_sol_pts(20); + int n_local_sol_pts = 0; + IntArray local_sol_tris(3,18); + int n_local_sol_tris; + DoubleArray values(20); + DTMutableList local_nws_pts(20); + int n_local_nws_pts; + + //int n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg; + int c; + //int newton_steps = 0; + //........................................................................... + int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo + IntArray cubeList(3,ncubes); + int nc=0; + //........................................................................... + // Set up the cube list (very regular in this case due to lack of blob-ID) + for (k=0; k CPU + //........................................................................... + dvc_Barrier(); + dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S); + dvc_CopyToHost(Phase.data,Phi,N*sizeof(double)); + dvc_CopyToHost(Press,Pressure,N*sizeof(double)); + dvc_CopyToHost(Vel,Velocity,3*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; + if (rank==0){ + printf("--------------------------------------------------------------------------------------\n"); + printf("timestep dEs "); // Timestep, Change in Surface Energy + printf("sw pw pn awn ans aws Jwn lwns efawns "); // Scalar averages + printf("vw[x, y, z] vn[x, y, z] vwn[x, y, z]"); // Velocity averages + printf("Gwn [xx, yy, zz, xy, xz, yz] "); // Orientation tensors + printf("Gws [xx, yy, zz, xy, xz, yz] "); + printf("Gns [xx, yy, zz, xy, xz, yz] \n"); + printf("--------------------------------------------------------------------------------------\n"); + } + + //************ MAIN ITERATION LOOP ***************************************/ + while (timestep < timestepMax){ + + //************************************************************************* + // Compute the color gradient + //************************************************************************* + //dvc_ComputeColorGradient(nBlocks, nthreads, S, + // ID, Phi, ColorGrad, Nx, Ny, Nz); + //************************************************************************* + + //************************************************************************* + // Perform collision step for the momentum transport + //************************************************************************* +// dvc_ColorCollide(nBlocks, nthreads, S, ID, f_even, f_odd, ColorGrad, Velocity, +// rlxA, rlxB,alpha, beta, Fx, Fy, Fz, Nx, Ny, Nz, pBC); + //************************************************************************* + + //************************************************************************* + // Fused Color Gradient and Collision + //************************************************************************* + dvc_ColorCollideOpt( ID,f_even,f_odd,Phi,ColorGrad, + Velocity,Nx,Ny,Nz,S,rlxA,rlxB,alpha,beta,Fx,Fy,Fz); + //************************************************************************* + + //................................................................................... + dvc_PackDist(1,dvcSendList_x,0,sendCount_x,sendbuf_x,f_even,N); + dvc_PackDist(4,dvcSendList_x,sendCount_x,sendCount_x,sendbuf_x,f_even,N); + dvc_PackDist(5,dvcSendList_x,2*sendCount_x,sendCount_x,sendbuf_x,f_even,N); + dvc_PackDist(6,dvcSendList_x,3*sendCount_x,sendCount_x,sendbuf_x,f_even,N); + dvc_PackDist(7,dvcSendList_x,4*sendCount_x,sendCount_x,sendbuf_x,f_even,N); + //...Packing for X face(1,7,9,11,13)................................ + dvc_PackDist(0,dvcSendList_X,0,sendCount_X,sendbuf_X,f_odd,N); + dvc_PackDist(3,dvcSendList_X,sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + dvc_PackDist(4,dvcSendList_X,2*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + dvc_PackDist(5,dvcSendList_X,3*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + dvc_PackDist(6,dvcSendList_X,4*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + //...Packing for y face(4,8,9,16,18)................................. + dvc_PackDist(2,dvcSendList_y,0,sendCount_y,sendbuf_y,f_even,N); + dvc_PackDist(4,dvcSendList_y,sendCount_y,sendCount_y,sendbuf_y,f_even,N); + dvc_PackDist(4,dvcSendList_y,2*sendCount_y,sendCount_y,sendbuf_y,f_odd,N); + dvc_PackDist(8,dvcSendList_y,3*sendCount_y,sendCount_y,sendbuf_y,f_even,N); + dvc_PackDist(9,dvcSendList_y,4*sendCount_y,sendCount_y,sendbuf_y,f_even,N); + //...Packing for Y face(3,7,10,15,17)................................. + dvc_PackDist(1,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,f_odd,N); + dvc_PackDist(3,dvcSendList_Y,sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); + dvc_PackDist(5,dvcSendList_Y,2*sendCount_Y,sendCount_Y,sendbuf_Y,f_even,N); + dvc_PackDist(7,dvcSendList_Y,3*sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); + dvc_PackDist(8,dvcSendList_Y,4*sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); + //...Packing for z face(6,12,13,16,17)................................ + dvc_PackDist(3,dvcSendList_z,0,sendCount_z,sendbuf_z,f_even,N); + dvc_PackDist(6,dvcSendList_z,sendCount_z,sendCount_z,sendbuf_z,f_even,N); + dvc_PackDist(6,dvcSendList_z,2*sendCount_z,sendCount_z,sendbuf_z,f_odd,N); + dvc_PackDist(8,dvcSendList_z,3*sendCount_z,sendCount_z,sendbuf_z,f_even,N); + dvc_PackDist(8,dvcSendList_z,4*sendCount_z,sendCount_z,sendbuf_z,f_odd,N); + //...Packing for Z face(5,11,14,15,18)................................ + dvc_PackDist(2,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,f_odd,N); + dvc_PackDist(5,dvcSendList_Z,sendCount_Z,sendCount_Z,sendbuf_Z,f_odd,N); + dvc_PackDist(7,dvcSendList_Z,2*sendCount_Z,sendCount_Z,sendbuf_Z,f_even,N); + dvc_PackDist(7,dvcSendList_Z,3*sendCount_Z,sendCount_Z,sendbuf_Z,f_odd,N); + dvc_PackDist(9,dvcSendList_Z,4*sendCount_Z,sendCount_Z,sendbuf_Z,f_even,N); + //...Pack the xy edge (8)................................ + dvc_PackDist(4,dvcSendList_xy,0,sendCount_xy,sendbuf_xy,f_even,N); + //...Pack the Xy edge (9)................................ + dvc_PackDist(4,dvcSendList_Xy,0,sendCount_Xy,sendbuf_Xy,f_odd,N); + //...Pack the xY edge (10)................................ + dvc_PackDist(5,dvcSendList_xY,0,sendCount_xY,sendbuf_xY,f_even,N); + //...Pack the XY edge (7)................................ + dvc_PackDist(3,dvcSendList_XY,0,sendCount_XY,sendbuf_XY,f_odd,N); + //...Pack the xz edge (12)................................ + dvc_PackDist(6,dvcSendList_xz,0,sendCount_xz,sendbuf_xz,f_even,N); + //...Pack the xZ edge (14)................................ + dvc_PackDist(7,dvcSendList_xZ,0,sendCount_xZ,sendbuf_xZ,f_even,N); + //...Pack the Xz edge (13)................................ + dvc_PackDist(6,dvcSendList_Xz,0,sendCount_Xz,sendbuf_Xz,f_odd,N); + //...Pack the XZ edge (11)................................ + dvc_PackDist(5,dvcSendList_XZ,0,sendCount_XZ,sendbuf_XZ,f_odd,N); + //...Pack the xz edge (12)................................ + //...Pack the yz edge (16)................................ + dvc_PackDist(8,dvcSendList_yz,0,sendCount_yz,sendbuf_yz,f_even,N); + //...Pack the yZ edge (18)................................ + dvc_PackDist(9,dvcSendList_yZ,0,sendCount_yZ,sendbuf_yZ,f_even,N); + //...Pack the Yz edge (17)................................ + dvc_PackDist(8,dvcSendList_Yz,0,sendCount_Yz,sendbuf_Yz,f_odd,N); + //...Pack the YZ edge (15)................................ + dvc_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 + //************************************************************************* + dvc_DensityStreamD3Q7(ID, Den, Copy, Phi, ColorGrad, Velocity, beta, Nx, Ny, Nz, pBC, S); + //************************************************************************* + + //************************************************************************* + // Swap the distributions for momentum transport + //************************************************************************* + dvc_SwapD3Q19(ID, f_even, f_odd, Nx, Ny, Nz, S); + //************************************************************************* + + //................................................................................... + // 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 ................................. + dvc_UnpackDist(0,-1,0,0,dvcRecvList_X,0,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(3,-1,-1,0,dvcRecvList_X,recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(4,-1,1,0,dvcRecvList_X,2*recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(5,-1,0,-1,dvcRecvList_X,3*recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + dvc_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.................................. + dvc_UnpackDist(1,1,0,0,dvcRecvList_x,0,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + dvc_UnpackDist(4,1,1,0,dvcRecvList_x,recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + dvc_UnpackDist(5,1,-1,0,dvcRecvList_x,2*recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + dvc_UnpackDist(6,1,0,1,dvcRecvList_x,3*recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + dvc_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 ................................... + dvc_UnpackDist(1,0,-1,0,dvcRecvList_Y,0,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(3,-1,-1,0,dvcRecvList_Y,recvCount_Y,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(5,1,-1,0,dvcRecvList_Y,2*recvCount_Y,recvCount_Y,recvbuf_Y,f_even,Nx,Ny,Nz); + dvc_UnpackDist(7,0,-1,-1,dvcRecvList_Y,3*recvCount_Y,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); + dvc_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 .................................. + dvc_UnpackDist(2,0,1,0,dvcRecvList_y,0,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); + dvc_UnpackDist(4,1,1,0,dvcRecvList_y,recvCount_y,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); + dvc_UnpackDist(4,-1,1,0,dvcRecvList_y,2*recvCount_y,recvCount_y,recvbuf_y,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(8,0,1,1,dvcRecvList_y,3*recvCount_y,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); + dvc_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).............................................. + dvc_UnpackDist(2,0,0,-1,dvcRecvList_Z,0,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(5,-1,0,-1,dvcRecvList_Z,recvCount_Z,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(7,1,0,-1,dvcRecvList_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,f_even,Nx,Ny,Nz); + dvc_UnpackDist(7,0,-1,-1,dvcRecvList_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); + dvc_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).............................................. + dvc_UnpackDist(3,0,0,1,dvcRecvList_z,0,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); + dvc_UnpackDist(6,1,0,1,dvcRecvList_z,recvCount_z,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); + dvc_UnpackDist(6,-1,0,1,dvcRecvList_z,2*recvCount_z,recvCount_z,recvbuf_z,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(8,0,1,1,dvcRecvList_z,3*recvCount_z,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_UnpackDist(8,0,1,1,dvcRecvList_yz,0,recvCount_yz,recvbuf_yz,f_even,Nx,Ny,Nz); + //................................................................................... + + //................................................................................... + dvc_PackDenD3Q7(dvcRecvList_x,recvCount_x,recvbuf_x,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_y,recvCount_y,recvbuf_y,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_z,recvCount_z,recvbuf_z,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_X,recvCount_X,recvbuf_X,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_Y,recvCount_Y,recvbuf_Y,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_Z,recvCount_Z,recvbuf_Z,2,Den,N); + //................................................................................... + + //................................................................................... + // Send all the D3Q7 distributions + MPI_Isend(recvbuf_x, 2*recvCount_x,MPI_DOUBLE,rank_x,sendtag,MPI_COMM_WORLD,&req1[0]); + MPI_Irecv(sendbuf_X, 2*sendCount_X,MPI_DOUBLE,rank_X,recvtag,MPI_COMM_WORLD,&req2[0]); + MPI_Isend(recvbuf_X, 2*recvCount_X,MPI_DOUBLE,rank_X,sendtag,MPI_COMM_WORLD,&req1[1]); + MPI_Irecv(sendbuf_x, 2*sendCount_x,MPI_DOUBLE,rank_x,recvtag,MPI_COMM_WORLD,&req2[1]); + MPI_Isend(recvbuf_y, 2*recvCount_y,MPI_DOUBLE,rank_y,sendtag,MPI_COMM_WORLD,&req1[2]); + MPI_Irecv(sendbuf_Y, 2*sendCount_Y,MPI_DOUBLE,rank_Y,recvtag,MPI_COMM_WORLD,&req2[2]); + MPI_Isend(recvbuf_Y, 2*recvCount_Y,MPI_DOUBLE,rank_Y,sendtag,MPI_COMM_WORLD,&req1[3]); + MPI_Irecv(sendbuf_y, 2*sendCount_y,MPI_DOUBLE,rank_y,recvtag,MPI_COMM_WORLD,&req2[3]); + MPI_Isend(recvbuf_z, 2*recvCount_z,MPI_DOUBLE,rank_z,sendtag,MPI_COMM_WORLD,&req1[4]); + MPI_Irecv(sendbuf_Z, 2*sendCount_Z,MPI_DOUBLE,rank_Z,recvtag,MPI_COMM_WORLD,&req2[4]); + MPI_Isend(recvbuf_Z, 2*recvCount_Z,MPI_DOUBLE,rank_Z,sendtag,MPI_COMM_WORLD,&req1[5]); + MPI_Irecv(sendbuf_z, 2*sendCount_z,MPI_DOUBLE,rank_z,recvtag,MPI_COMM_WORLD,&req2[5]); + //................................................................................... + //................................................................................... + // Wait for completion of D3Q7 communication + MPI_Waitall(6,req1,stat1); + MPI_Waitall(6,req2,stat2); + //................................................................................... + //................................................................................... + dvc_UnpackDenD3Q7(dvcSendList_x,sendCount_x,sendbuf_x,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_y,sendCount_y,sendbuf_y,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_z,sendCount_z,sendbuf_z,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_X,sendCount_X,sendbuf_X,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_Y,sendCount_Y,sendbuf_Y,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_Z,sendCount_Z,sendbuf_Z,2,Den,N); + //................................................................................... + + //************************************************************************* + // Compute the phase indicator field and reset Copy, Den + //************************************************************************* + dvc_ComputePhi(ID, Phi, Copy, Den, N, S); + //************************************************************************* + + //................................................................................... + dvc_PackValues(dvcSendList_x, sendCount_x,sendbuf_x, Phi, N); + dvc_PackValues(dvcSendList_y, sendCount_y,sendbuf_y, Phi, N); + dvc_PackValues(dvcSendList_z, sendCount_z,sendbuf_z, Phi, N); + dvc_PackValues(dvcSendList_X, sendCount_X,sendbuf_X, Phi, N); + dvc_PackValues(dvcSendList_Y, sendCount_Y,sendbuf_Y, Phi, N); + dvc_PackValues(dvcSendList_Z, sendCount_Z,sendbuf_Z, Phi, N); + dvc_PackValues(dvcSendList_xy, sendCount_xy,sendbuf_xy, Phi, N); + dvc_PackValues(dvcSendList_xY, sendCount_xY,sendbuf_xY, Phi, N); + dvc_PackValues(dvcSendList_Xy, sendCount_Xy,sendbuf_Xy, Phi, N); + dvc_PackValues(dvcSendList_XY, sendCount_XY,sendbuf_XY, Phi, N); + dvc_PackValues(dvcSendList_xz, sendCount_xz,sendbuf_xz, Phi, N); + dvc_PackValues(dvcSendList_xZ, sendCount_xZ,sendbuf_xZ, Phi, N); + dvc_PackValues(dvcSendList_Xz, sendCount_Xz,sendbuf_Xz, Phi, N); + dvc_PackValues(dvcSendList_XZ, sendCount_XZ,sendbuf_XZ, Phi, N); + dvc_PackValues(dvcSendList_yz, sendCount_yz,sendbuf_yz, Phi, N); + dvc_PackValues(dvcSendList_yZ, sendCount_yZ,sendbuf_yZ, Phi, N); + dvc_PackValues(dvcSendList_Yz, sendCount_Yz,sendbuf_Yz, Phi, N); + dvc_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); + dvc_Barrier(); + //................................................................................... + //................................................................................... +/* dvc_UnpackValues(faceGrid, packThreads, dvcSendList_x, sendCount_x,sendbuf_x, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_y, sendCount_y,sendbuf_y, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_z, sendCount_z,sendbuf_z, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_X, sendCount_X,sendbuf_X, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_Y, sendCount_Y,sendbuf_Y, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_Z, sendCount_Z,sendbuf_Z, Phi, N); +*/ + dvc_UnpackValues(dvcRecvList_x, recvCount_x,recvbuf_x, Phi, N); + dvc_UnpackValues(dvcRecvList_y, recvCount_y,recvbuf_y, Phi, N); + dvc_UnpackValues(dvcRecvList_z, recvCount_z,recvbuf_z, Phi, N); + dvc_UnpackValues(dvcRecvList_X, recvCount_X,recvbuf_X, Phi, N); + dvc_UnpackValues(dvcRecvList_Y, recvCount_Y,recvbuf_Y, Phi, N); + dvc_UnpackValues(dvcRecvList_Z, recvCount_Z,recvbuf_Z, Phi, N); + dvc_UnpackValues(dvcRecvList_xy, recvCount_xy,recvbuf_xy, Phi, N); + dvc_UnpackValues(dvcRecvList_xY, recvCount_xY,recvbuf_xY, Phi, N); + dvc_UnpackValues(dvcRecvList_Xy, recvCount_Xy,recvbuf_Xy, Phi, N); + dvc_UnpackValues(dvcRecvList_XY, recvCount_XY,recvbuf_XY, Phi, N); + dvc_UnpackValues(dvcRecvList_xz, recvCount_xz,recvbuf_xz, Phi, N); + dvc_UnpackValues(dvcRecvList_xZ, recvCount_xZ,recvbuf_xZ, Phi, N); + dvc_UnpackValues(dvcRecvList_Xz, recvCount_Xz,recvbuf_Xz, Phi, N); + dvc_UnpackValues(dvcRecvList_XZ, recvCount_XZ,recvbuf_XZ, Phi, N); + dvc_UnpackValues(dvcRecvList_yz, recvCount_yz,recvbuf_yz, Phi, N); + dvc_UnpackValues(dvcRecvList_yZ, recvCount_yZ,recvbuf_yZ, Phi, N); + dvc_UnpackValues(dvcRecvList_Yz, recvCount_Yz,recvbuf_Yz, Phi, N); + dvc_UnpackValues(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, Phi, N); + //................................................................................... + MPI_Barrier(MPI_COMM_WORLD); + + // Iteration completed! + timestep++; + //................................................................... + if (timestep%1000 == 995){ + //........................................................................... + // Copy the phase indicator field for the earlier timestep + dvc_Barrier(); + dvc_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 + //........................................................................... + dvc_Barrier(); + dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S); + dvc_CopyToHost(Phase.data,Phi,N*sizeof(double)); + dvc_CopyToHost(Press,Pressure,N*sizeof(double)); + dvc_CopyToHost(Vel,Velocity,3*N*sizeof(double)); + MPI_Barrier(MPI_COMM_WORLD); + } + if (timestep%1000 == 5){ + //........................................................................... + // Copy the phase indicator field for the later timestep + dvc_Barrier(); + dvc_CopyToHost(Phase_tminus.data,Phi,N*sizeof(double)); + //........................................................................... + // Calculate the time derivative of the phase indicator field + for (n=0; n 0 ){ + + // 1-D index for this cube corner + n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny; + + // Compute the 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.9999 ){ + // volume the excludes the interfacial region + vol_n += 0.125; + // pressure + pan += 0.125*Press[n]; + // velocity + van(0) += 0.125*Vel[3*n]; + van(1) += 0.125*Vel[3*n+1]; + van(2) += 0.125*Vel[3*n+2]; + } + + // volume averages over the wetting phase + if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < -0.9999 ){ + // volume the excludes the interfacial region + vol_w += 0.125; + // pressure + paw += 0.125*Press[n]; + // velocity + vaw(0) += 0.125*Vel[3*n]; + vaw(1) += 0.125*Vel[3*n+1]; + vaw(2) += 0.125*Vel[3*n+2]; + } + } + } + + //........................................................................... + // 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); + + pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris, + NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris); + + //........................................................................... + // Compute the Interfacial Areas, Common Line length + /* awn += pmmc_CubeSurfaceArea(nw_pts,nw_tris,n_nw_tris); + ans += pmmc_CubeSurfaceArea(ns_pts,ns_tris,n_ns_tris); + aws += pmmc_CubeSurfaceArea(ws_pts,ws_tris,n_ws_tris); + */ + As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris); + + // Compute the surface orientation and the interfacial area + awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris); + ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris); + aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris); + lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts); + //........................................................................... + } + //........................................................................... + MPI_Barrier(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(&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_Barrier(MPI_COMM_WORLD); + //......................................................................... + // Compute the change in the total surface energy based on the defined interval + // See McClure, Prins and Miller (2013) + //......................................................................... + dAwn += awn_global; + dAns += ans_global; + dEs = 6.01603*alpha*(dAwn + 1.05332*Ps*dAns); + dAwn = -awn_global; // Get ready for the next analysis interval + dAns = -ans_global; + + // Normalize the phase averages + // (density of both components = 1.0) + paw_global = paw_global / vol_w_global; + vaw_global(0) = vaw_global(0) / vol_w_global; + vaw_global(1) = vaw_global(1) / vol_w_global; + vaw_global(2) = vaw_global(2) / vol_w_global; + pan_global = pan_global / vol_n_global; + van_global(0) = van_global(0) / vol_n_global; + van_global(1) = van_global(1) / vol_n_global; + van_global(2) = van_global(2) / vol_n_global; + + // Normalize surface averages by the interfacial area + Jwn_global /= awn_global; + efawns_global /= lwns_global; + + if (awn_global > 0.0) for (i=0; i<3; i++) vawn_global(i) /= awn_global; + if (awn_global > 0.0) for (i=0; i<6; i++) Gwn_global(i) /= awn_global; + if (ans_global > 0.0) for (i=0; i<6; i++) Gns_global(i) /= ans_global; + if (aws_global > 0.0) for (i=0; i<6; i++) Gws_global(i) /= aws_global; + + sat_w = 1.0 - nwp_volume_global*iVol_global/porosity; + // Compute the specific interfacial areas and common line length (per unit volume) + awn_global = awn_global*iVol_global; + ans_global = ans_global*iVol_global; + aws_global = aws_global*iVol_global; + lwns_global = lwns_global*iVol_global; + dEs = dEs*iVol_global; + + //......................................................................... + if (rank==0){ +/* printf("-------------------------------- \n"); + printf("Timestep = %i \n", timestep); + printf("NWP volume = %f \n", nwp_volume_global); + printf("Area wn = %f \n", awn_global); + printf("Area ns = %f \n", ans_global); + printf("Area ws = %f \n", aws_global); + printf("Change in surface energy = %f \n", dEs); + printf("-------------------------------- \n"); + + printf("%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g \n", + timestep,dEs,sat_w, + awn_global,ans_global,aws_global, lwns_global, p_w_global, p_n_global, + vx_w_global, vy_w_global, vz_w_global, + vx_n_global, vy_n_global, vz_n_global); +*/ + 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 ",Jwn_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 \n", + Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface + } + } + + if (timestep%RESTART_INTERVAL == 0){ + // Copy the data to the CPU + dvc_CopyToHost(cDistEven,f_even,10*N*sizeof(double)); + dvc_CopyToHost(cDistOdd,f_odd,9*N*sizeof(double)); + dvc_CopyToHost(cDen,Copy,2*N*sizeof(double)); + // Read in the restart file to CPU buffers + WriteCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N); + } + } + //************************************************************************/ + dvc_Barrier(); + MPI_Barrier(MPI_COMM_WORLD); + stoptime = MPI_Wtime(); + if (rank==0) printf("-------------------------------------------------------------------\n"); + // Compute the walltime per timestep + cputime = (stoptime - starttime)/timestep; + // Performance obtained from each node + double MLUPS = double(Nx*Ny*Nz)/cputime/1000000; + + if (rank==0) printf("********************************************************\n"); + if (rank==0) printf("CPU time = %f \n", cputime); + if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); + MLUPS *= nprocs; + if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); + if (rank==0) printf("********************************************************\n"); + + //************************************************************************/ + // Write out the phase indicator field + //************************************************************************/ + sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString); + // printf("Local File Name = %s \n",LocalRankFilename); +// dvc_CopyToHost(Phase.data,Phi,N*sizeof(double)); + +//#ifdef WriteOutput + FILE *PHASE; + PHASE = fopen(LocalRankFilename,"wb"); + fwrite(Phase.data,8,N,PHASE); +// fwrite(MeanCurvature.data,8,N,PHASE); + fclose(PHASE); +//#endif + +/* double *DensityValues; + DensityValues = new double [2*N]; + dvc_CopyToHost(DensityValues,Copy,2*N*sizeof(double)); + FILE *PHASE; + PHASE = fopen(LocalRankFilename,"wb"); + fwrite(DensityValues,8,2*N,PHASE); + fclose(PHASE); +*/ //************************************************************************/ + + // **************************************************** + MPI_Barrier(MPI_COMM_WORLD); + MPI_Finalize(); + // **************************************************** +} diff --git a/cpu/lb2_Color_wia_mpi_bubble.cpp b/cpu/lb2_Color_wia_mpi_bubble.cpp new file mode 100644 index 00000000..2be9fbb4 --- /dev/null +++ b/cpu/lb2_Color_wia_mpi_bubble.cpp @@ -0,0 +1,2492 @@ +#include +#include +#include +#include +#include + +#include "pmmc.h" +#include "Domain.h" +#include "Extras.h" +#include "D3Q19.h" +#include "D3Q7.h" +#include "Color.h" + +using namespace std; + +//************************************************************************* +// Implementation of Two-Phase Immiscible LBM using CUDA +//************************************************************************* +inline void PackID(int *list, int count, char *sendbuf, char *ID){ + // Fill in the phase ID values from neighboring processors + // This packs up the values that need to be sent from one processor to another + int idx,n; + + for (idx=0; idx> 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 >> xIntPos; // Contact angle parameter +// 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; // error interval + input >> tol; // error tolerance + //............................................................. + das = 0.1; dbs = 0.9; // hard coded for density initialization + // should be OK to remove these parameters + // they should have no impact with the + // current boundary condition + //....................................................................... + // Reading the domain information file + //....................................................................... + ifstream domain("Domain.in"); + 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(&xIntPos,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(&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); + // ************************************************************** + // ************************************************************** + double Ps = -(das-dbs)/(das+dbs); + double rlxA = 1.f/tau; + double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA); + + if (nprocs != nprocx*nprocy*nprocz){ + printf("Fatal error in processor number! \n"); + printf("nprocx = %i \n",nprocx); + printf("nprocy = %i \n",nprocy); + printf("nprocz = %i \n",nprocz); + } + + 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("phi_s = %f \n", Ps); + printf("gamma_{wn} = %f \n", 6.01603*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"); + } + + MPI_Barrier(MPI_COMM_WORLD); + kproc = rank/(nprocx*nprocy); + jproc = (rank-nprocx*nprocy*kproc)/nprocx; + iproc = rank-nprocx*nprocy*kproc-nprocz*jproc; + + //.......................................... + // set up the neighbor ranks + //.......................................... + i=iproc; j=jproc; k =kproc; + i+=1; + j+=0; + k+=0; + if (i<0) i+=nprocx; + if (j<0) j+=nprocy; + if (k<0) k+=nprocz; + if (!(i 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] = 1; + sum++; + } + } + } + } + sum_local = 1.0*sum; + MPI_Allreduce(&sum_local,&porosity,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + porosity = porosity*iVol_global; + if (rank==0) printf("Media porosity = %f \n",porosity); + + // Generate the residual NWP + if (rank==0) printf("Initializing with NWP saturation = %f \n",wp_saturation); + GenerateResidual(id,Nx,Ny,Nz,wp_saturation); +#endif + + double BubbleRadius = 15.5; // Radius of the capillary tube + sum=0; + for (k=0;k fluid_isovalue) + int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners + DoubleArray CubeValues(2,2,2); +// int count_in=0,count_out=0; +// int nodx,nody,nodz; + // initialize lists for vertices for surfaces, common line + DTMutableList nw_pts(20); + DTMutableList ns_pts(20); + DTMutableList ws_pts(20); + DTMutableList nws_pts(20); + // initialize triangle lists for surfaces + IntArray nw_tris(3,20); + IntArray ns_tris(3,20); + IntArray ws_tris(3,20); + // initialize list for line segments + IntArray nws_seg(2,20); + DTMutableList tmp(20); + DoubleArray Values(20); + DoubleArray ContactAngle(20); + DoubleArray Curvature(20); + DoubleArray InterfaceSpeed(20); + DoubleArray NormalVector(60); + + // IntArray store; + + int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0; + int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0; + +// double s,s1,s2,s3; // Triangle sides (lengths) + Point A,B,C,P; +// double area; + + // Initialize arrays for local solid surface + DTMutableList local_sol_pts(20); + int n_local_sol_pts = 0; + IntArray local_sol_tris(3,18); + int n_local_sol_tris; + DoubleArray values(20); + DTMutableList local_nws_pts(20); + int n_local_nws_pts; + + //int n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg; + int c; + //int newton_steps = 0; + //........................................................................... + int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo + IntArray cubeList(3,ncubes); + int nc=0; + //........................................................................... + // Set up the cube list (very regular in this case due to lack of blob-ID) + for (k=0; k CPU + //........................................................................... + dvc_Barrier(); + dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S); + dvc_CopyToHost(Phase.data,Phi,N*sizeof(double)); + dvc_CopyToHost(Press.data,Pressure,N*sizeof(double)); + dvc_CopyToHost(Vel,Velocity,3*N*sizeof(double)); + MPI_Barrier(MPI_COMM_WORLD); + //........................................................................... + + timestep=0; + + sendtag = recvtag = 5; + + //************ MAIN ITERATION LOOP ***************************************/ + while (timestep < timestepMax){ + + //************************************************************************* + // Compute the color gradient + //************************************************************************* + //dvc_ComputeColorGradient(nBlocks, nthreads, S, + // ID, Phi, ColorGrad, Nx, Ny, Nz); + //************************************************************************* + + //************************************************************************* + // Perform collision step for the momentum transport + //************************************************************************* + // dvc_ColorCollide(nBlocks, nthreads, S, ID, f_even, f_odd, ColorGrad, Velocity, + // rlxA, rlxB,alpha, beta, Fx, Fy, Fz, Nx, Ny, Nz, pBC); + //************************************************************************* + + //************************************************************************* + // Fused Color Gradient and Collision + //************************************************************************* + dvc_ColorCollideOpt( ID,f_even,f_odd,Phi,ColorGrad, + Velocity,Nx,Ny,Nz,S,rlxA,rlxB,alpha,beta,Fx,Fy,Fz); + //************************************************************************* + + //................................................................................... + dvc_PackDist(1,dvcSendList_x,0,sendCount_x,sendbuf_x,f_even,N); + dvc_PackDist(4,dvcSendList_x,sendCount_x,sendCount_x,sendbuf_x,f_even,N); + dvc_PackDist(5,dvcSendList_x,2*sendCount_x,sendCount_x,sendbuf_x,f_even,N); + dvc_PackDist(6,dvcSendList_x,3*sendCount_x,sendCount_x,sendbuf_x,f_even,N); + dvc_PackDist(7,dvcSendList_x,4*sendCount_x,sendCount_x,sendbuf_x,f_even,N); + //...Packing for X face(1,7,9,11,13)................................ + dvc_PackDist(0,dvcSendList_X,0,sendCount_X,sendbuf_X,f_odd,N); + dvc_PackDist(3,dvcSendList_X,sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + dvc_PackDist(4,dvcSendList_X,2*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + dvc_PackDist(5,dvcSendList_X,3*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + dvc_PackDist(6,dvcSendList_X,4*sendCount_X,sendCount_X,sendbuf_X,f_odd,N); + //...Packing for y face(4,8,9,16,18)................................. + dvc_PackDist(2,dvcSendList_y,0,sendCount_y,sendbuf_y,f_even,N); + dvc_PackDist(4,dvcSendList_y,sendCount_y,sendCount_y,sendbuf_y,f_even,N); + dvc_PackDist(4,dvcSendList_y,2*sendCount_y,sendCount_y,sendbuf_y,f_odd,N); + dvc_PackDist(8,dvcSendList_y,3*sendCount_y,sendCount_y,sendbuf_y,f_even,N); + dvc_PackDist(9,dvcSendList_y,4*sendCount_y,sendCount_y,sendbuf_y,f_even,N); + //...Packing for Y face(3,7,10,15,17)................................. + dvc_PackDist(1,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,f_odd,N); + dvc_PackDist(3,dvcSendList_Y,sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); + dvc_PackDist(5,dvcSendList_Y,2*sendCount_Y,sendCount_Y,sendbuf_Y,f_even,N); + dvc_PackDist(7,dvcSendList_Y,3*sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); + dvc_PackDist(8,dvcSendList_Y,4*sendCount_Y,sendCount_Y,sendbuf_Y,f_odd,N); + //...Packing for z face(6,12,13,16,17)................................ + dvc_PackDist(3,dvcSendList_z,0,sendCount_z,sendbuf_z,f_even,N); + dvc_PackDist(6,dvcSendList_z,sendCount_z,sendCount_z,sendbuf_z,f_even,N); + dvc_PackDist(6,dvcSendList_z,2*sendCount_z,sendCount_z,sendbuf_z,f_odd,N); + dvc_PackDist(8,dvcSendList_z,3*sendCount_z,sendCount_z,sendbuf_z,f_even,N); + dvc_PackDist(8,dvcSendList_z,4*sendCount_z,sendCount_z,sendbuf_z,f_odd,N); + //...Packing for Z face(5,11,14,15,18)................................ + dvc_PackDist(2,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,f_odd,N); + dvc_PackDist(5,dvcSendList_Z,sendCount_Z,sendCount_Z,sendbuf_Z,f_odd,N); + dvc_PackDist(7,dvcSendList_Z,2*sendCount_Z,sendCount_Z,sendbuf_Z,f_even,N); + dvc_PackDist(7,dvcSendList_Z,3*sendCount_Z,sendCount_Z,sendbuf_Z,f_odd,N); + dvc_PackDist(9,dvcSendList_Z,4*sendCount_Z,sendCount_Z,sendbuf_Z,f_even,N); + //...Pack the xy edge (8)................................ + dvc_PackDist(4,dvcSendList_xy,0,sendCount_xy,sendbuf_xy,f_even,N); + //...Pack the Xy edge (9)................................ + dvc_PackDist(4,dvcSendList_Xy,0,sendCount_Xy,sendbuf_Xy,f_odd,N); + //...Pack the xY edge (10)................................ + dvc_PackDist(5,dvcSendList_xY,0,sendCount_xY,sendbuf_xY,f_even,N); + //...Pack the XY edge (7)................................ + dvc_PackDist(3,dvcSendList_XY,0,sendCount_XY,sendbuf_XY,f_odd,N); + //...Pack the xz edge (12)................................ + dvc_PackDist(6,dvcSendList_xz,0,sendCount_xz,sendbuf_xz,f_even,N); + //...Pack the xZ edge (14)................................ + dvc_PackDist(7,dvcSendList_xZ,0,sendCount_xZ,sendbuf_xZ,f_even,N); + //...Pack the Xz edge (13)................................ + dvc_PackDist(6,dvcSendList_Xz,0,sendCount_Xz,sendbuf_Xz,f_odd,N); + //...Pack the XZ edge (11)................................ + dvc_PackDist(5,dvcSendList_XZ,0,sendCount_XZ,sendbuf_XZ,f_odd,N); + //...Pack the xz edge (12)................................ + //...Pack the yz edge (16)................................ + dvc_PackDist(8,dvcSendList_yz,0,sendCount_yz,sendbuf_yz,f_even,N); + //...Pack the yZ edge (18)................................ + dvc_PackDist(9,dvcSendList_yZ,0,sendCount_yZ,sendbuf_yZ,f_even,N); + //...Pack the Yz edge (17)................................ + dvc_PackDist(8,dvcSendList_Yz,0,sendCount_Yz,sendbuf_Yz,f_odd,N); + //...Pack the YZ edge (15)................................ + dvc_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 + //************************************************************************* + dvc_DensityStreamD3Q7(ID, Den, Copy, Phi, ColorGrad, Velocity, beta, Nx, Ny, Nz, pBC, S); + //************************************************************************* + + //************************************************************************* + // Swap the distributions for momentum transport + //************************************************************************* + dvc_SwapD3Q19(ID, f_even, f_odd, Nx, Ny, Nz, S); + //************************************************************************* + + //................................................................................... + // 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 ................................. + dvc_UnpackDist(0,-1,0,0,dvcRecvList_X,0,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(3,-1,-1,0,dvcRecvList_X,recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(4,-1,1,0,dvcRecvList_X,2*recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(5,-1,0,-1,dvcRecvList_X,3*recvCount_X,recvCount_X,recvbuf_X,f_odd,Nx,Ny,Nz); + dvc_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.................................. + dvc_UnpackDist(1,1,0,0,dvcRecvList_x,0,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + dvc_UnpackDist(4,1,1,0,dvcRecvList_x,recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + dvc_UnpackDist(5,1,-1,0,dvcRecvList_x,2*recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + dvc_UnpackDist(6,1,0,1,dvcRecvList_x,3*recvCount_x,recvCount_x,recvbuf_x,f_even,Nx,Ny,Nz); + dvc_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 ................................... + dvc_UnpackDist(1,0,-1,0,dvcRecvList_Y,0,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(3,-1,-1,0,dvcRecvList_Y,recvCount_Y,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(5,1,-1,0,dvcRecvList_Y,2*recvCount_Y,recvCount_Y,recvbuf_Y,f_even,Nx,Ny,Nz); + dvc_UnpackDist(7,0,-1,-1,dvcRecvList_Y,3*recvCount_Y,recvCount_Y,recvbuf_Y,f_odd,Nx,Ny,Nz); + dvc_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 .................................. + dvc_UnpackDist(2,0,1,0,dvcRecvList_y,0,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); + dvc_UnpackDist(4,1,1,0,dvcRecvList_y,recvCount_y,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); + dvc_UnpackDist(4,-1,1,0,dvcRecvList_y,2*recvCount_y,recvCount_y,recvbuf_y,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(8,0,1,1,dvcRecvList_y,3*recvCount_y,recvCount_y,recvbuf_y,f_even,Nx,Ny,Nz); + dvc_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).............................................. + dvc_UnpackDist(2,0,0,-1,dvcRecvList_Z,0,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(5,-1,0,-1,dvcRecvList_Z,recvCount_Z,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(7,1,0,-1,dvcRecvList_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,f_even,Nx,Ny,Nz); + dvc_UnpackDist(7,0,-1,-1,dvcRecvList_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,f_odd,Nx,Ny,Nz); + dvc_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).............................................. + dvc_UnpackDist(3,0,0,1,dvcRecvList_z,0,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); + dvc_UnpackDist(6,1,0,1,dvcRecvList_z,recvCount_z,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); + dvc_UnpackDist(6,-1,0,1,dvcRecvList_z,2*recvCount_z,recvCount_z,recvbuf_z,f_odd,Nx,Ny,Nz); + dvc_UnpackDist(8,0,1,1,dvcRecvList_z,3*recvCount_z,recvCount_z,recvbuf_z,f_even,Nx,Ny,Nz); + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_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)................................ + dvc_UnpackDist(8,0,1,1,dvcRecvList_yz,0,recvCount_yz,recvbuf_yz,f_even,Nx,Ny,Nz); + //................................................................................... + + //................................................................................... + dvc_PackDenD3Q7(dvcRecvList_x,recvCount_x,recvbuf_x,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_y,recvCount_y,recvbuf_y,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_z,recvCount_z,recvbuf_z,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_X,recvCount_X,recvbuf_X,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_Y,recvCount_Y,recvbuf_Y,2,Den,N); + dvc_PackDenD3Q7(dvcRecvList_Z,recvCount_Z,recvbuf_Z,2,Den,N); + //................................................................................... + + //................................................................................... + // Send all the D3Q7 distributions + MPI_Isend(recvbuf_x, 2*recvCount_x,MPI_DOUBLE,rank_x,sendtag,MPI_COMM_WORLD,&req1[0]); + MPI_Irecv(sendbuf_X, 2*sendCount_X,MPI_DOUBLE,rank_X,recvtag,MPI_COMM_WORLD,&req2[0]); + MPI_Isend(recvbuf_X, 2*recvCount_X,MPI_DOUBLE,rank_X,sendtag,MPI_COMM_WORLD,&req1[1]); + MPI_Irecv(sendbuf_x, 2*sendCount_x,MPI_DOUBLE,rank_x,recvtag,MPI_COMM_WORLD,&req2[1]); + MPI_Isend(recvbuf_y, 2*recvCount_y,MPI_DOUBLE,rank_y,sendtag,MPI_COMM_WORLD,&req1[2]); + MPI_Irecv(sendbuf_Y, 2*sendCount_Y,MPI_DOUBLE,rank_Y,recvtag,MPI_COMM_WORLD,&req2[2]); + MPI_Isend(recvbuf_Y, 2*recvCount_Y,MPI_DOUBLE,rank_Y,sendtag,MPI_COMM_WORLD,&req1[3]); + MPI_Irecv(sendbuf_y, 2*sendCount_y,MPI_DOUBLE,rank_y,recvtag,MPI_COMM_WORLD,&req2[3]); + MPI_Isend(recvbuf_z, 2*recvCount_z,MPI_DOUBLE,rank_z,sendtag,MPI_COMM_WORLD,&req1[4]); + MPI_Irecv(sendbuf_Z, 2*sendCount_Z,MPI_DOUBLE,rank_Z,recvtag,MPI_COMM_WORLD,&req2[4]); + MPI_Isend(recvbuf_Z, 2*recvCount_Z,MPI_DOUBLE,rank_Z,sendtag,MPI_COMM_WORLD,&req1[5]); + MPI_Irecv(sendbuf_z, 2*sendCount_z,MPI_DOUBLE,rank_z,recvtag,MPI_COMM_WORLD,&req2[5]); + //................................................................................... + //................................................................................... + // Wait for completion of D3Q7 communication + MPI_Waitall(6,req1,stat1); + MPI_Waitall(6,req2,stat2); + //................................................................................... + //................................................................................... + dvc_UnpackDenD3Q7(dvcSendList_x,sendCount_x,sendbuf_x,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_y,sendCount_y,sendbuf_y,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_z,sendCount_z,sendbuf_z,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_X,sendCount_X,sendbuf_X,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_Y,sendCount_Y,sendbuf_Y,2,Den,N); + dvc_UnpackDenD3Q7(dvcSendList_Z,sendCount_Z,sendbuf_Z,2,Den,N); + //................................................................................... + + //************************************************************************* + // Compute the phase indicator field and reset Copy, Den + //************************************************************************* + dvc_ComputePhi(ID, Phi, Copy, Den, N, S); + //************************************************************************* + + //................................................................................... + dvc_PackValues(dvcSendList_x, sendCount_x,sendbuf_x, Phi, N); + dvc_PackValues(dvcSendList_y, sendCount_y,sendbuf_y, Phi, N); + dvc_PackValues(dvcSendList_z, sendCount_z,sendbuf_z, Phi, N); + dvc_PackValues(dvcSendList_X, sendCount_X,sendbuf_X, Phi, N); + dvc_PackValues(dvcSendList_Y, sendCount_Y,sendbuf_Y, Phi, N); + dvc_PackValues(dvcSendList_Z, sendCount_Z,sendbuf_Z, Phi, N); + dvc_PackValues(dvcSendList_xy, sendCount_xy,sendbuf_xy, Phi, N); + dvc_PackValues(dvcSendList_xY, sendCount_xY,sendbuf_xY, Phi, N); + dvc_PackValues(dvcSendList_Xy, sendCount_Xy,sendbuf_Xy, Phi, N); + dvc_PackValues(dvcSendList_XY, sendCount_XY,sendbuf_XY, Phi, N); + dvc_PackValues(dvcSendList_xz, sendCount_xz,sendbuf_xz, Phi, N); + dvc_PackValues(dvcSendList_xZ, sendCount_xZ,sendbuf_xZ, Phi, N); + dvc_PackValues(dvcSendList_Xz, sendCount_Xz,sendbuf_Xz, Phi, N); + dvc_PackValues(dvcSendList_XZ, sendCount_XZ,sendbuf_XZ, Phi, N); + dvc_PackValues(dvcSendList_yz, sendCount_yz,sendbuf_yz, Phi, N); + dvc_PackValues(dvcSendList_yZ, sendCount_yZ,sendbuf_yZ, Phi, N); + dvc_PackValues(dvcSendList_Yz, sendCount_Yz,sendbuf_Yz, Phi, N); + dvc_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); + dvc_Barrier(); + //................................................................................... + //................................................................................... + /* dvc_UnpackValues(faceGrid, packThreads, dvcSendList_x, sendCount_x,sendbuf_x, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_y, sendCount_y,sendbuf_y, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_z, sendCount_z,sendbuf_z, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_X, sendCount_X,sendbuf_X, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_Y, sendCount_Y,sendbuf_Y, Phi, N); + dvc_UnpackValues(faceGrid, packThreads, dvcSendList_Z, sendCount_Z,sendbuf_Z, Phi, N); + */ + dvc_UnpackValues(dvcRecvList_x, recvCount_x,recvbuf_x, Phi, N); + dvc_UnpackValues(dvcRecvList_y, recvCount_y,recvbuf_y, Phi, N); + dvc_UnpackValues(dvcRecvList_z, recvCount_z,recvbuf_z, Phi, N); + dvc_UnpackValues(dvcRecvList_X, recvCount_X,recvbuf_X, Phi, N); + dvc_UnpackValues(dvcRecvList_Y, recvCount_Y,recvbuf_Y, Phi, N); + dvc_UnpackValues(dvcRecvList_Z, recvCount_Z,recvbuf_Z, Phi, N); + dvc_UnpackValues(dvcRecvList_xy, recvCount_xy,recvbuf_xy, Phi, N); + dvc_UnpackValues(dvcRecvList_xY, recvCount_xY,recvbuf_xY, Phi, N); + dvc_UnpackValues(dvcRecvList_Xy, recvCount_Xy,recvbuf_Xy, Phi, N); + dvc_UnpackValues(dvcRecvList_XY, recvCount_XY,recvbuf_XY, Phi, N); + dvc_UnpackValues(dvcRecvList_xz, recvCount_xz,recvbuf_xz, Phi, N); + dvc_UnpackValues(dvcRecvList_xZ, recvCount_xZ,recvbuf_xZ, Phi, N); + dvc_UnpackValues(dvcRecvList_Xz, recvCount_Xz,recvbuf_Xz, Phi, N); + dvc_UnpackValues(dvcRecvList_XZ, recvCount_XZ,recvbuf_XZ, Phi, N); + dvc_UnpackValues(dvcRecvList_yz, recvCount_yz,recvbuf_yz, Phi, N); + dvc_UnpackValues(dvcRecvList_yZ, recvCount_yZ,recvbuf_yZ, Phi, N); + dvc_UnpackValues(dvcRecvList_Yz, recvCount_Yz,recvbuf_Yz, Phi, N); + dvc_UnpackValues(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, Phi, N); + //................................................................................... + MPI_Barrier(MPI_COMM_WORLD); + + // Iteration completed! + timestep++; + + if (timestep%RESTART_INTERVAL == 0){ + // Copy the data to the CPU + dvc_CopyToHost(cDistEven,f_even,10*N*sizeof(double)); + dvc_CopyToHost(cDistOdd,f_odd,9*N*sizeof(double)); + dvc_CopyToHost(cDen,Copy,2*N*sizeof(double)); + // Read in the restart file to CPU buffers + WriteCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N); + } + } + // End the bubble loop + //........................................................................... + // Copy the phase indicator field for the later timestep + dvc_Barrier(); + dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S); + dvc_CopyToHost(Phase_tminus.data,Phi,N*sizeof(double)); + dvc_CopyToHost(Phase_tplus.data,Phi,N*sizeof(double)); + dvc_CopyToHost(Phase.data,Phi,N*sizeof(double)); + dvc_CopyToHost(Press.data,Pressure,N*sizeof(double)); + + //........................................................................... + // Calculate the time derivative of the phase indicator field + for (n=0; n 0 ){ + + // 1-D index for this cube corner + n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny; + + // Compute the 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.99999 ){ + // volume the excludes the interfacial region + vol_n += 0.125; + // pressure + pan += 0.125*Press.data[n]; + // velocity + van(0) += 0.125*Vel[3*n]; + van(1) += 0.125*Vel[3*n+1]; + van(2) += 0.125*Vel[3*n+2]; + } + + // volume averages over the wetting phase + if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < -0.99999 ){ + // volume the excludes the interfacial region + vol_w += 0.125; + // pressure + paw += 0.125*Press.data[n]; + // velocity + vaw(0) += 0.125*Vel[3*n]; + vaw(1) += 0.125*Vel[3*n+1]; + vaw(2) += 0.125*Vel[3*n+2]; + } + } + } + + //........................................................................... + // 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); + + pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris, + NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris); + + //........................................................................... + // Compute the Interfacial Areas, Common Line length + /* awn += pmmc_CubeSurfaceArea(nw_pts,nw_tris,n_nw_tris); + ans += pmmc_CubeSurfaceArea(ns_pts,ns_tris,n_ns_tris); + aws += pmmc_CubeSurfaceArea(ws_pts,ws_tris,n_ws_tris); + */ + As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris); + + // Compute the surface orientation and the interfacial area + awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris); + ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris); + aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris); + lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts); + //........................................................................... + } + //........................................................................... + MPI_Barrier(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(&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_Barrier(MPI_COMM_WORLD); + //......................................................................... + // Compute the change in the total surface energy based on the defined interval + // See McClure, Prins and Miller (2013) + //......................................................................... + dAwn += awn_global; + dAns += ans_global; + dEs = 6.01603*alpha*(dAwn + 1.05332*Ps*dAns); + dAwn = -awn_global; // Get ready for the next analysis interval + dAns = -ans_global; + + // Normalize the phase averages + // (density of both components = 1.0) + paw_global = paw_global / vol_w_global; + vaw_global(0) = vaw_global(0) / vol_w_global; + vaw_global(1) = vaw_global(1) / vol_w_global; + vaw_global(2) = vaw_global(2) / vol_w_global; + pan_global = pan_global / vol_n_global; + van_global(0) = van_global(0) / vol_n_global; + van_global(1) = van_global(1) / vol_n_global; + van_global(2) = van_global(2) / vol_n_global; + + // Normalize surface averages by the interfacial area + Jwn_global /= awn_global; + efawns_global /= lwns_global; + + if (awn_global > 0.0) for (i=0; i<3; i++) vawn_global(i) /= awn_global; + if (awn_global > 0.0) for (i=0; i<6; i++) Gwn_global(i) /= awn_global; + if (ans_global > 0.0) for (i=0; i<6; i++) Gns_global(i) /= ans_global; + if (aws_global > 0.0) for (i=0; i<6; i++) Gws_global(i) /= aws_global; + + sat_w = 1.0 - nwp_volume_global*iVol_global/porosity; + // Compute the specific interfacial areas and common line length (per unit volume) + awn_global = awn_global*iVol_global; + ans_global = ans_global*iVol_global; + aws_global = aws_global*iVol_global; + lwns_global = lwns_global*iVol_global; + dEs = dEs*iVol_global; + + //......................................................................... + if (rank==0){ + /* printf("-------------------------------- \n"); + printf("Timestep = %i \n", timestep); + printf("NWP volume = %f \n", nwp_volume_global); + printf("Area wn = %f \n", awn_global); + printf("Area ns = %f \n", ans_global); + printf("Area ws = %f \n", aws_global); + printf("Change in surface energy = %f \n", dEs); + printf("-------------------------------- \n"); + + printf("%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g \n", + timestep,dEs,sat_w, + awn_global,ans_global,aws_global, lwns_global, p_w_global, p_n_global, + vx_w_global, vy_w_global, vz_w_global, + vx_n_global, vy_n_global, vz_n_global); + */ + printf("%.5g ",BubbleRadius); // bubble radius + printf("%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure + printf("%.5g ",awn_global); // interfacial area + printf("%.5g ",Jwn_global); // curvature of wn interface + printf("%.5g %.5g %.5g %.5g %.5g %.5g \n", + Gwn_global(0),Gwn_global(1),Gwn_global(2),Gwn_global(3),Gwn_global(4),Gwn_global(5)); // orientation of wn interface + } + + } + //************************************************************************/ + dvc_Barrier(); + MPI_Barrier(MPI_COMM_WORLD); + stoptime = MPI_Wtime(); + if (rank==0) printf("-------------------------------------------------------------------\n"); + // Compute the walltime per timestep + cputime = (stoptime - starttime)/timestep; + // Performance obtained from each node + double MLUPS = double(Nx*Ny*Nz)/cputime/1000000; + + if (rank==0) printf("********************************************************\n"); + if (rank==0) printf("CPU time = %f \n", cputime); + if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); + MLUPS *= nprocs; + if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); + if (rank==0) printf("********************************************************\n"); + + //************************************************************************/ + // Write out the phase indicator field + //************************************************************************/ + sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString); + // printf("Local File Name = %s \n",LocalRankFilename); +// dvc_CopyToHost(Phase.data,Phi,N*sizeof(double)); + + FILE *PHASE; + PHASE = fopen(LocalRankFilename,"wb"); + fwrite(Press.data,8,N,PHASE); +// fwrite(MeanCurvature.data,8,N,PHASE); + fclose(PHASE); + +/* double *DensityValues; + DensityValues = new double [2*N]; + dvc_CopyToHost(DensityValues,Copy,2*N*sizeof(double)); + FILE *PHASE; + PHASE = fopen(LocalRankFilename,"wb"); + fwrite(DensityValues,8,2*N,PHASE); + fclose(PHASE); +*/ //************************************************************************/ + + // **************************************************** + MPI_Barrier(MPI_COMM_WORLD); + MPI_Finalize(); + // **************************************************** +} diff --git a/include/Array.h b/include/Array.h index 75fd70c6..2059c205 100644 --- a/include/Array.h +++ b/include/Array.h @@ -312,14 +312,14 @@ double DoubleArray::e(int i, int j, int k) extern DoubleArray IncreaseSize(DoubleArray &A, int addLength) { if (addLength<0) { - printf("IncreaseSize(Array,Length)","Length needs to be >0."); + printf("IncreaseSize(Array,Length) Length needs to be >0."); return DoubleArray(); } int newM,newN,newO; if (A.o>1) { if (addLength%(A.m*A.n)!=0) { - printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n"); + printf("IncreaseSize(Array,Length) Length needs to be a multiple of m*n"); return DoubleArray(); } newM = A.m; @@ -328,7 +328,7 @@ extern DoubleArray IncreaseSize(DoubleArray &A, int addLength) } else if (A.n>1) { if (addLength%(A.m)!=0) { - printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m"); + printf("IncreaseSize(Array,Length) Length needs to be a multiple of m"); return DoubleArray(); } newM = A.m; @@ -348,14 +348,14 @@ extern DoubleArray IncreaseSize(DoubleArray &A, int addLength) extern IntArray IncreaseSize(IntArray &A, int addLength) { if (addLength<0) { - printf("IncreaseSize(Array,Length)","Length needs to be >0."); + printf("IncreaseSize(Array,Length) Length needs to be >0."); return IntArray(); } int newM,newN,newO; if (A.o>1) { if (addLength%(A.m*A.n)!=0) { - printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n"); + printf("IncreaseSize(Array,Length) Length needs to be a multiple of m*n"); return IntArray(); } newM = A.m; @@ -364,7 +364,7 @@ extern IntArray IncreaseSize(IntArray &A, int addLength) } else if (A.n>1) { if (addLength%(A.m)!=0) { - printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m"); + printf("IncreaseSize(Array,Length) Length needs to be a multiple of m"); return IntArray(); } newM = A.m; diff --git a/include/Domain.h b/include/Domain.h index 19cf74c7..12fdc279 100755 --- a/include/Domain.h +++ b/include/Domain.h @@ -399,3 +399,18 @@ inline void ReadCheckpoint(char *FILENAME, double *cDen, double *cDistEven, doub } File.close(); } + +inline void ReadBinaryFile(char *FILENAME, double *Data, int N) +{ + int n; + double value; + ifstream File(FILENAME,ios::binary); + for (n=0; n 0.0){ + // Surface value (speed) + vA = SurfaceValues(Triangles(0,r)); + vB = SurfaceValues(Triangles(1,r)); + vC = SurfaceValues(Triangles(2,r)); // Increment the averaged values // x component - vA = SurfaceVector(Triangles(0,r))*SurfaceValues(Triangles(0,r)); - vB = SurfaceVector(Triangles(1,r))*SurfaceValues(Triangles(1,r)); - vC = SurfaceVector(Triangles(2,r))*SurfaceValues(Triangles(2,r)); - AvgVel(0) += sqrt(temp)*0.33333333333333333*(vA+vB+vC); + vAx = SurfaceVector(Triangles(0,r))*vA; + vBx = SurfaceVector(Triangles(1,r))*vB; + vCx = SurfaceVector(Triangles(2,r))*vC; // y component - vA = SurfaceVector(npts+Triangles(0,r))*SurfaceValues(Triangles(0,r)); - vB = SurfaceVector(npts+Triangles(1,r))*SurfaceValues(Triangles(1,r)); - vC = SurfaceVector(npts+Triangles(2,r))*SurfaceValues(Triangles(2,r)); - AvgVel(1) += sqrt(temp)*0.33333333333333333*(vA+vB+vC); + vAy = SurfaceVector(npts+Triangles(0,r))*vA; + vBy = SurfaceVector(npts+Triangles(1,r))*vB; + vCy = SurfaceVector(npts+Triangles(2,r))*vC; // z component - vA = SurfaceVector(2*npts+Triangles(0,r))*SurfaceValues(Triangles(0,r)); - vB = SurfaceVector(2*npts+Triangles(1,r))*SurfaceValues(Triangles(1,r)); - vC = SurfaceVector(2*npts+Triangles(2,r))*SurfaceValues(Triangles(2,r)); - AvgVel(2) += sqrt(temp)*0.33333333333333333*(vA+vB+vC); + vAz = SurfaceVector(2*npts+Triangles(0,r))*vA; + vBz = SurfaceVector(2*npts+Triangles(1,r))*vB; + vCz = SurfaceVector(2*npts+Triangles(2,r))*vC; + + AvgVel(0) += sqrt(temp)*0.33333333333333333*(vAx+vBx+vCx); + AvgVel(1) += sqrt(temp)*0.33333333333333333*(vAy+vBy+vCy); + AvgVel(2) += sqrt(temp)*0.33333333333333333*(vAz+vBz+vCz); + + // Update the Averages. Differentiate between advancing (0,1,2) and receding (3,4,5) interfaces + // All points on a triangle have the same orientation in the color gradient +/* if (vA > 0.0){ + // Advancing interface + AvgVel(0) += sqrt(temp)*0.33333333333333333*(vAx+vBx+vCx); + AvgVel(1) += sqrt(temp)*0.33333333333333333*(vAy+vBy+vCy); + AvgVel(2) += sqrt(temp)*0.33333333333333333*(vAz+vBz+vCz); + } + else{ + // Receding interface + AvgVel(3) += sqrt(temp)*0.33333333333333333*(vAx+vBx+vCx); + AvgVel(4) += sqrt(temp)*0.33333333333333333*(vAy+vBy+vCy); + AvgVel(5) += sqrt(temp)*0.33333333333333333*(vAz+vBz+vCz); + } +*/ } } //............................................................................. diff --git a/pmmc/Analysis.cpp b/pmmc/Analysis.cpp index 17e6d636..0c2def7e 100644 --- a/pmmc/Analysis.cpp +++ b/pmmc/Analysis.cpp @@ -15,12 +15,17 @@ int main(int argc, char **argv) { //....................................................................... // printf("Radius = %s \n,"RADIUS); - int Nx,Ny,Nz; + int Nx,Ny,Nz,N; int i,j,k,p,q,r,n; int nspheres; double Lx,Ly,Lz; //....................................................................... Nx = Ny = Nz = 60; + cout << "Enter Domain size " << endl; + cout << "Nx = " << endl; + cin >> Nx; + Ny = Nz = Nx; + N = Nx*Ny*Nz; //....................................................................... // Reading the domain information file /* //....................................................................... @@ -129,6 +134,8 @@ int main(int argc, char **argv) DoubleArray Gns(6); DoubleArray Gws(6); + double iVol = 1.0/Nx/Ny/Nz; + int c; //........................................................................... int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo @@ -159,8 +166,8 @@ int main(int argc, char **argv) dist1 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)) - RADIUS; dist2 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)+(k-Cz)*(k-Cz)) - CAPRAD; - SignDist(i,j,k) = -dist1; - Phase(i,j,k) = dist2; + //SignDist(i,j,k) = -dist1; + //Phase(i,j,k) = dist2; } } } @@ -175,7 +182,41 @@ int main(int argc, char **argv) } } - //........................................................................... + awn = aws = ans = lwns = 0.0; + nwp_volume = 0.0; + As = 0.0; + Jwn = 0.0; + efawns = 0.0; + // Compute phase averages + pan = paw = 0.0; + vaw(0) = vaw(1) = vaw(2) = 0.0; + van(0) = van(1) = van(2) = 0.0; + vawn(0) = vawn(1) = vawn(2) = 0.0; + Gwn(0) = Gwn(1) = Gwn(2) = 0.0; + Gwn(3) = Gwn(4) = Gwn(5) = 0.0; + Gws(0) = Gws(1) = Gws(2) = 0.0; + Gws(3) = Gws(4) = Gws(5) = 0.0; + Gns(0) = Gns(1) = Gns(2) = 0.0; + Gns(3) = Gns(4) = Gns(5) = 0.0; + vol_w = vol_n =0.0; + + // Read the input files for the phase, distance and pressure field + char PHASEFILE[16]; + sprintf(PHASEFILE,"Phase.in"); + ReadBinaryFile(PHASEFILE,Phase.data,Nx*Ny*Nz); + char DISTFILE[16]; + sprintf(DISTFILE,"SignDist.in"); + ReadBinaryFile(DISTFILE,SignDist.data,Nx*Ny*Nz); + /* FILE *PRESS + PRESS = fopen("Pressure.in","wb"); + fread(Phase.data,8,N,PRESS); + fclose(PRESS); + + FILE *VEL; + VEL = fopen("Pressure.in","wb"); + fread(Phase.data,8,3*N,VEL); + fclose(VEL); +*/ //........................................................................... // Calculate the time derivative of the phase indicator field for (int n=0; n