From 92d56af3b4aa306b75fd2de82b515e90944291f4 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 25 Sep 2020 16:18:54 -0400 Subject: [PATCH 01/39] template for freelee model --- cpu/FreeLee.cpp | 2820 +++++++++++++++++++++++++++++++++++++++ models/FreeLeeModel.cpp | 632 +++++++++ models/FreeLeeModel.h | 83 ++ 3 files changed, 3535 insertions(+) create mode 100644 cpu/FreeLee.cpp create mode 100644 models/FreeLeeModel.cpp create mode 100644 models/FreeLeeModel.h diff --git a/cpu/FreeLee.cpp b/cpu/FreeLee.cpp new file mode 100644 index 00000000..35cbd5fd --- /dev/null +++ b/cpu/FreeLee.cpp @@ -0,0 +1,2820 @@ +#include + +#define STOKES + +extern "C" void ScaLBL_Color_Init(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz) +{ + int n,N; + + N = Nx*Ny*Nz; + + for (n=0; n 0){ + + // Retrieve the color gradient + nx = ColorGrad[n]; + ny = ColorGrad[N+n]; + nz = ColorGrad[2*N+n]; + //...........Normalize the Color Gradient................................. + C = sqrt(nx*nx+ny*ny+nz*nz); + if (C==0.0) C=1.0; + nx = nx/C; + ny = ny/C; + nz = nz/C; + //......No color gradient at z-boundary if pressure BC are set............. + // if (pBC && k==0) nx = ny = nz = 0.f; + // if (pBC && k==Nz-1) nx = ny = nz = 0.f; + //........................................................................ + // READ THE DISTRIBUTIONS + // (read from opposite array due to previous swap operation) + //........................................................................ + f2 = distodd[n]; + f4 = distodd[N+n]; + f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; + //........................................................................ + f0 = disteven[n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + f17 = disteven[9*N+n]; + //........................................................................ + // PERFORM RELAXATION PROCESS + //........................................................................ + //....................compute the moments............................................... + rho = f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17; + m1 = -30*f0-11*(f2+f1+f4+f3+f6+f5)+8*(f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18 +f17); + m2 = 12*f0-4*(f2+f1 +f4+f3+f6 +f5)+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17; + jx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14; + m4 = 4*(-f1+f2)+f7-f8+f9-f10+f11-f12+f13-f14; + jy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18; + m6 = -4*(f3-f4)+f7-f8-f9+f10+f15-f16+f17-f18; + jz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18; + m8 = -4*(f5-f6)+f11-f12-f13+f14+f15-f16-f17+f18; + m9 = 2*(f1+f2)-f3-f4-f5-f6+f7+f8+f9+f10+f11+f12+f13+f14-2*(f15+f16+f17+f18); + m10 = -4*(f1+f2)+2*(f4+f3+f6+f5)+f8+f7+f10+f9+f12+f11+f14+f13-2*(f16+f15+f18+f17); + m11 = f4+f3-f6-f5+f8+f7+f10+f9-f12-f11-f14-f13; + m12 = -2*(f4+f3-f6-f5)+f8+f7+f10+f9-f12-f11-f14-f13; + m13 = f8+f7-f10-f9; + m14 = f16+f15-f18-f17; + m15 = f12+f11-f14-f13; + m16 = f7-f8+f9-f10-f11+f12-f13+f14; + m17 = -f7+f8+f9-f10+f15-f16+f17-f18; + m18 = f11-f12-f13+f14-f15+f16+f17-f18; + //..........Toelke, Fruediger et. al. 2006............... + if (C == 0.0) nx = ny = nz = 1.0; +#ifdef STOKES + m1 = m1 + rlx_setA*(- 11*rho -alpha*C - m1); + m2 = m2 + rlx_setA*(3*rho - m2); + m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); + m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); + m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); + m9 = m9 + rlx_setA*( 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); + m10 = m10 + rlx_setA*( - m10); + m11 = m11 + rlx_setA*( 0.5*alpha*C*(ny*ny-nz*nz)- m11); + m12 = m12 + rlx_setA*( - m12); + m13 = m13 + rlx_setA*( 0.5*alpha*C*nx*ny - m13); + m14 = m14 + rlx_setA*( 0.5*alpha*C*ny*nz - m14); + m15 = m15 + rlx_setA*( 0.5*alpha*C*nx*nz - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); +#else + m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) -alpha*C - m1); + m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho)- m2); + m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); + m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); + m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); + m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); + m10 = m10 + rlx_setA*( - m10); + m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); + m12 = m12 + rlx_setA*( - m12); + m13 = m13 + rlx_setA*( (jx*jy/rho) + 0.5*alpha*C*nx*ny - m13); + m14 = m14 + rlx_setA*( (jy*jz/rho) + 0.5*alpha*C*ny*nz - m14); + m15 = m15 + rlx_setA*( (jx*jz/rho) + 0.5*alpha*C*nx*nz - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); +#endif + //.................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); + f2 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +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); + f4 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +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); + f6 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +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); + 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); + 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); + 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); + 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); + 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); + 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); + 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); + 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); + 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); + 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); + 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); + //....................................................................................................... + // 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[n] = jx; + Velocity[N+n] = jy; + Velocity[2*N+n] = jz; + /* 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 ScaLBL_D3Q19_ColorCollide( char *ID, double *disteven, double *distodd, double *phi, 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) +{ + + 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+10)) delta=0; + a1 = na*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nb*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = na*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nb*(0.1111111111111111*(1-4.5*ux))+delta; + + A_odd[n] = a1; + A_even[N+n] = a2; + B_odd[n] = b1; + B_even[N+n] = b2; + //............................................... + // q = 2 + // Cq = {0,1,0} + delta = beta*na*nb*nab*0.1111111111111111*ny; + if (!(na*nb*nab>0)) delta=0; + a1 = na*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nb*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = na*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nb*(0.1111111111111111*(1-4.5*uy))+delta; + + A_odd[N+n] = a1; + A_even[2*N+n] = a2; + B_odd[N+n] = b1; + B_even[2*N+n] = b2; + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*na*nb*nab*0.1111111111111111*nz; + if (!(na*nb*nab>0)) delta=0; + a1 = na*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nb*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = na*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nb*(0.1111111111111111*(1-4.5*uz))+delta; + + A_odd[2*N+n] = a1; + A_even[3*N+n] = a2; + B_odd[2*N+n] = b1; + B_even[3*N+n] = b2; + //............................................... + + /* // Construction and streaming for the components + for (idx=0; idx<3; idx++){ + //............................................... + // Distribution index + q = 2*idx; + // Associated discrete velocity + Cqx = D3Q7[idx][0]; + Cqy = D3Q7[idx][1]; + Cqz = D3Q7[idx][2]; + // Generate the Equilibrium Distribution + a1 = na*feq[q]; + b1 = nb*feq[q]; + a2 = na*feq[q+1]; + b2 = nb*feq[q+1]; + // Recolor the distributions + if (C > 0.0){ + sp = nx*double(Cqx)+ny*double(Cqy)+nz*double(Cqz); + //if (idx > 2) sp = 0.7071067811865475*sp; + //delta = sp*min( min(a1,a2), min(b1,b2) ); + delta = na*nb/(na+nb)*0.1111111111111111*sp; + //if (a1>0 && b1>0){ + a1 += beta*delta; + a2 -= beta*delta; + b1 -= beta*delta; + b2 += beta*delta; + } + // Save the re-colored distributions + A_odd[N*idx+n] = a1; + A_even[N*(idx+1)+n] = a2; + B_odd[N*idx+n] = b1; + B_even[N*(idx+1)+n] = b2; + //............................................... + } + */ + } + } +} + +//************************************************************************* +extern "C" void 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) +{ + char id; + + int idx; + int in,jn,kn,n,nn,N; + int q,Cqx,Cqy,Cqz; + // int sendLoc; + + double na,nb; // density values + double ux,uy,uz; // flow velocity + double nx,ny,nz,C; // color gradient components + double a1,a2,b1,b2; + double sp,delta; + double feq[6]; // equilibrium distributions + // Set of Discrete velocities for the D3Q19 Model + int D3Q7[3][3]={{1,0,0},{0,1,0},{0,0,1}}; + N = Nx*Ny*Nz; + + for (n=0; n 0 && na+nb > 0.0){ + //.......Back out the 3-D indices for node n.............. + int k = n/(Nx*Ny); + int j = (n-Nx*Ny*k)/Nx; + int i = n-Nx*Ny*k-Nx*j; + //.....Load the Color gradient......... + nx = ColorGrad[n]; + ny = ColorGrad[N+n]; + nz = ColorGrad[2*N+n]; + C = sqrt(nx*nx+ny*ny+nz*nz); + nx = nx/C; + ny = ny/C; + nz = nz/C; + //....Load the flow velocity........... + ux = Velocity[n]; + uy = Velocity[N+n]; + uz = Velocity[2*N+n]; + //....Instantiate the density distributions + // Generate Equilibrium Distributions and stream + // Stationary value - distribution 0 + // Den[2*n] += 0.3333333333333333*na; + // Den[2*n+1] += 0.3333333333333333*nb; + Den[2*n] += 0.3333333333333333*na; + Den[2*n+1] += 0.3333333333333333*nb; + // Non-Stationary equilibrium distributions + feq[0] = 0.1111111111111111*(1+3*ux); + feq[1] = 0.1111111111111111*(1-3*ux); + feq[2] = 0.1111111111111111*(1+3*uy); + feq[3] = 0.1111111111111111*(1-3*uy); + feq[4] = 0.1111111111111111*(1+3*uz); + feq[5] = 0.1111111111111111*(1-3*uz); + // Construction and streaming for the components + for (idx=0; idx<3; idx++){ + // Distribution index + q = 2*idx; + // Associated discrete velocity + Cqx = D3Q7[idx][0]; + Cqy = D3Q7[idx][1]; + Cqz = D3Q7[idx][2]; + // Generate the Equilibrium Distribution + a1 = na*feq[q]; + b1 = nb*feq[q]; + a2 = na*feq[q+1]; + b2 = nb*feq[q+1]; + // Recolor the distributions + if (C > 0.0){ + sp = nx*double(Cqx)+ny*double(Cqy)+nz*double(Cqz); + //if (idx > 2) sp = 0.7071067811865475*sp; + //delta = sp*min( min(a1,a2), min(b1,b2) ); + delta = na*nb/(na+nb)*0.1111111111111111*sp; + //if (a1>0 && b1>0){ + a1 += beta*delta; + a2 -= beta*delta; + b1 -= beta*delta; + b2 += beta*delta; + } + + // .......Get the neighbor node.............. + //nn = n + Stride[idx]; + in = i+Cqx; + jn = j+Cqy; + kn = k+Cqz; + + // Adjust for periodic BC, if necessary + // if (in<0) in+= Nx; + // if (jn<0) jn+= Ny; + // if (kn<0) kn+= Nz; + // if (!(in 0 ){ + // Get the density value (Streaming already performed) + Na = Den[n]; + Nb = Den[N+n]; + Phi[n] = (Na-Nb)/(Na+Nb); + } + } + //................................................................... +} + +extern "C" void ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, int Nz, int Slice){ + int n; + for (n=Slice*Nx*Ny; n<(Slice+1)*Nx*Ny; n++){ + Phi[n] = value; + } +} + + +//extern "C" void ScaLBL_D3Q19_AAeven_Color(double *dist, double *Aq, double *Bq, double *Den, double *Velocity, +// double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, +// double Fx, double Fy, double Fz, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ + + int ijk,nn,n; + double fq; + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + double m3,m5,m7; + double nA,nB; // number density + double a1,b1,a2,b2,nAB,delta; + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + double phi,tau,rho0,rlx_setA,rlx_setB; + + const double mrt_V1=0.05263157894736842; + const double mrt_V2=0.012531328320802; + const double mrt_V3=0.04761904761904762; + const double mrt_V4=0.004594820384294068; + const double mrt_V5=0.01587301587301587; + const double mrt_V6=0.0555555555555555555555555; + const double mrt_V7=0.02777777777777778; + const double mrt_V8=0.08333333333333333; + const double mrt_V9=0.003341687552213868; + const double mrt_V10=0.003968253968253968; + const double mrt_V11=0.01388888888888889; + const double mrt_V12=0.04166666666666666; + + + for (int n=start; n0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + Aq[1*Np+n] = a1; + Bq[1*Np+n] = b1; + Aq[2*Np+n] = a2; + Bq[2*Np+n] = b2; + + //............................................... + // q = 2 + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + Aq[3*Np+n] = a1; + Bq[3*Np+n] = b1; + Aq[4*Np+n] = a2; + Bq[4*Np+n] = b2; + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + Aq[5*Np+n] = a1; + Bq[5*Np+n] = b1; + Aq[6*Np+n] = a2; + Bq[6*Np+n] = b2; + //............................................... + + } + +} + +//extern "C" void ScaLBL_D3Q19_AAodd_Color(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Velocity, +// double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, +// double Fx, double Fy, double Fz, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, + double *Phi, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ + + int n,nn,ijk,nread; + int nr1,nr2,nr3,nr4,nr5,nr6; + int nr7,nr8,nr9,nr10; + int nr11,nr12,nr13,nr14; + //int nr15,nr16,nr17,nr18; + double fq; + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + double m3,m5,m7; + double nA,nB; // number density + double a1,b1,a2,b2,nAB,delta; + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + double phi,tau,rho0,rlx_setA,rlx_setB; + + const double mrt_V1=0.05263157894736842; + const double mrt_V2=0.012531328320802; + const double mrt_V3=0.04761904761904762; + const double mrt_V4=0.004594820384294068; + const double mrt_V5=0.01587301587301587; + const double mrt_V6=0.0555555555555555555555555; + const double mrt_V7=0.02777777777777778; + const double mrt_V8=0.08333333333333333; + const double mrt_V9=0.003341687552213868; + const double mrt_V10=0.003968253968253968; + const double mrt_V11=0.01388888888888889; + const double mrt_V12=0.04166666666666666; + + for (int n=start; n even part of dist) + //fq = dist[nread]; // reading the f2 data into register fq + nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) + fq = dist[nr2]; // reading the f2 data into register fq + rho += fq; + m1 -= 11.0*(fq); + m2 -= 4.0*(fq); + jx -= fq; + m4 += 4.0*(fq); + m9 += 2.0*(fq); + m10 -= 4.0*(fq); + + // q=3 + //nread = neighborList[n+2*Np]; // neighbor 4 + //fq = dist[nread]; + nr3 = neighborList[n+2*Np]; // neighbor 4 + fq = dist[nr3]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy = fq; + m6 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 = fq; + m12 = -2.0*fq; + + // q = 4 + //nread = neighborList[n+3*Np]; // neighbor 3 + //fq = dist[nread]; + nr4 = neighborList[n+3*Np]; // neighbor 3 + fq = dist[nr4]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy -= fq; + m6 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 += fq; + m12 -= 2.0*fq; + + // q=5 + //nread = neighborList[n+4*Np]; + //fq = dist[nread]; + nr5 = neighborList[n+4*Np]; + fq = dist[nr5]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz = fq; + m8 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + + // q = 6 + //nread = neighborList[n+5*Np]; + //fq = dist[nread]; + nr6 = neighborList[n+5*Np]; + fq = dist[nr6]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz -= fq; + m8 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + // q=7 + //nread = neighborList[n+6*Np]; + //fq = dist[nread]; + nr7 = neighborList[n+6*Np]; + fq = dist[nr7]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 = fq; + m16 = fq; + m17 = -fq; + + // q = 8 + //nread = neighborList[n+7*Np]; + //fq = dist[nread]; + nr8 = neighborList[n+7*Np]; + fq = dist[nr8]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 += fq; + m16 -= fq; + m17 += fq; + + // q=9 + //nread = neighborList[n+8*Np]; + //fq = dist[nread]; + nr9 = neighborList[n+8*Np]; + fq = dist[nr9]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 += fq; + m17 += fq; + + // q = 10 + //nread = neighborList[n+9*Np]; + //fq = dist[nread]; + nr10 = neighborList[n+9*Np]; + fq = dist[nr10]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 -= fq; + m17 -= fq; + + // q=11 + //nread = neighborList[n+10*Np]; + //fq = dist[nread]; + nr11 = neighborList[n+10*Np]; + fq = dist[nr11]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 = fq; + m16 -= fq; + m18 = fq; + + // q=12 + //nread = neighborList[n+11*Np]; + //fq = dist[nread]; + nr12 = neighborList[n+11*Np]; + fq = dist[nr12]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 += fq; + m16 += fq; + m18 -= fq; + + // q=13 + //nread = neighborList[n+12*Np]; + //fq = dist[nread]; + nr13 = neighborList[n+12*Np]; + fq = dist[nr13]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 -= fq; + m18 -= fq; + + // q=14 + //nread = neighborList[n+13*Np]; + //fq = dist[nread]; + nr14 = neighborList[n+13*Np]; + fq = dist[nr14]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 += fq; + m18 += fq; + + // q=15 + nread = neighborList[n+14*Np]; + fq = dist[nread]; + //fq = dist[17*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 = fq; + m17 += fq; + m18 -= fq; + + // q=16 + nread = neighborList[n+15*Np]; + fq = dist[nread]; + //fq = dist[8*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 += fq; + m17 -= fq; + m18 += fq; + + // q=17 + //fq = dist[18*Np+n]; + nread = neighborList[n+16*Np]; + fq = dist[nread]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 += fq; + m18 += fq; + + // q=18 + nread = neighborList[n+17*Np]; + fq = dist[nread]; + //fq = dist[9*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 -= fq; + m18 -= fq; + + //........................................................................ + //..............carry out relaxation process.............................. + //..........Toelke, Fruediger et. al. 2006................................ + if (C == 0.0) nx = ny = nz = 0.0; + m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -19*alpha*C - m1); + m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); + m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); + m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); + m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); + m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); + m10 = m10 + rlx_setA*( - m10); + m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); + m12 = m12 + rlx_setA*( - m12); + m13 = m13 + rlx_setA*( (jx*jy/rho0) + 0.5*alpha*C*nx*ny - m13); + m14 = m14 + rlx_setA*( (jy*jz/rho0) + 0.5*alpha*C*ny*nz - m14); + m15 = m15 + rlx_setA*( (jx*jz/rho0) + 0.5*alpha*C*nx*nz - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); + //.................inverse transformation...................................................... + + // q=0 + fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; + dist[n] = fq; + + // q = 1 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10)+0.16666666*Fx; + //nread = neighborList[n+Np]; + dist[nr2] = fq; + + // q=2 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10) - 0.16666666*Fx; + //nread = neighborList[n]; + dist[nr1] = fq; + + // q = 3 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) + 0.16666666*Fy; + //nread = neighborList[n+3*Np]; + dist[nr4] = fq; + + // q = 4 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) - 0.16666666*Fy; + //nread = neighborList[n+2*Np]; + dist[nr3] = fq; + + // q = 5 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) + 0.16666666*Fz; + //nread = neighborList[n+5*Np]; + dist[nr6] = fq; + + // q = 6 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) - 0.16666666*Fz; + //nread = neighborList[n+4*Np]; + dist[nr5] = fq; + + // q = 7 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17) + 0.08333333333*(Fx+Fy); + //nread = neighborList[n+7*Np]; + dist[nr8] = fq; + + // q = 8 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11 + +mrt_V12*m12+0.25*m13+0.125*(m17-m16) - 0.08333333333*(Fx+Fy); + //nread = neighborList[n+6*Np]; + dist[nr7] = fq; + + // q = 9 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17) + 0.08333333333*(Fx-Fy); + //nread = neighborList[n+9*Np]; + dist[nr10] = fq; + + // q = 10 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17)- 0.08333333333*(Fx-Fy); + //nread = neighborList[n+8*Np]; + dist[nr9] = fq; + + // q = 11 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12+0.25*m15+0.125*(m18-m16) + 0.08333333333*(Fx+Fz); + //nread = neighborList[n+11*Np]; + dist[nr12] = fq; + + // q = 12 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+ + mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18) - 0.08333333333*(Fx+Fz); + //nread = neighborList[n+10*Np]; + dist[nr11]= fq; + + // q = 13 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15-0.125*(m16+m18) + 0.08333333333*(Fx-Fz); + //nread = neighborList[n+13*Np]; + dist[nr14] = fq; + + // q= 14 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15+0.125*(m16+m18) - 0.08333333333*(Fx-Fz); + //nread = neighborList[n+12*Np]; + dist[nr13] = fq; + + + // q = 15 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18) + 0.08333333333*(Fy+Fz); + nread = neighborList[n+15*Np]; + dist[nread] = fq; + + // q = 16 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17)- 0.08333333333*(Fy+Fz); + nread = neighborList[n+14*Np]; + dist[nread] = fq; + + + // q = 17 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) + -mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18) + 0.08333333333*(Fy-Fz); + nread = neighborList[n+17*Np]; + dist[nread] = fq; + + // q = 18 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6) + -mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18) - 0.08333333333*(Fy-Fz); + nread = neighborList[n+16*Np]; + dist[nread] = fq; + + // write the velocity + ux = jx / rho0; + uy = jy / rho0; + uz = jz / rho0; + Vel[n] = ux; + Vel[Np+n] = uy; + Vel[2*Np+n] = uz; + + // Instantiate mass transport distributions + // Stationary value - distribution 0 + nAB = 1.0/(nA+nB); + Aq[n] = 0.3333333333333333*nA; + Bq[n] = 0.3333333333333333*nB; + + //............................................... + // q = 0,2,4 + // Cq = {1,0,0}, {0,1,0}, {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nx; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + // q = 1 + //nread = neighborList[n+Np]; + Aq[nr2] = a1; + Bq[nr2] = b1; + // q=2 + //nread = neighborList[n]; + Aq[nr1] = a2; + Bq[nr1] = b2; + + //............................................... + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + // q = 3 + //nread = neighborList[n+3*Np]; + Aq[nr4] = a1; + Bq[nr4] = b1; + // q = 4 + //nread = neighborList[n+2*Np]; + Aq[nr3] = a2; + Bq[nr3] = b2; + + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + // q = 5 + //nread = neighborList[n+5*Np]; + Aq[nr6] = a1; + Bq[nr6] = b1; + // q = 6 + //nread = neighborList[n+4*Np]; + Aq[nr5] = a2; + Bq[nr5] = b2; + //............................................... + } +} + +extern "C" void ScaLBL_D3Q7_AAodd_PhaseField(int *neighborList, int *Map, double *Aq, double *Bq, + double *Den, double *Phi, int start, int finish, int Np){ + + int idx,n,nread; + double fq,nA,nB; + + for (int n=start; n 1.f){ + nA = 1.0; nB = 0.f; + } + else if (phi < -1.f){ + nB = 1.0; nA = 0.f; + } + else{ + nA=0.5*(phi+1.f); + nB=0.5*(1.f-phi); + } + Den[idx] = nA; + Den[Np+idx] = nB; + + Aq[idx]=0.3333333333333333*nA; + Aq[Np+idx]=0.1111111111111111*nA; + Aq[2*Np+idx]=0.1111111111111111*nA; + Aq[3*Np+idx]=0.1111111111111111*nA; + Aq[4*Np+idx]=0.1111111111111111*nA; + Aq[5*Np+idx]=0.1111111111111111*nA; + Aq[6*Np+idx]=0.1111111111111111*nA; + + Bq[idx]=0.3333333333333333*nB; + Bq[Np+idx]=0.1111111111111111*nB; + Bq[2*Np+idx]=0.1111111111111111*nB; + Bq[3*Np+idx]=0.1111111111111111*nB; + Bq[4*Np+idx]=0.1111111111111111*nB; + Bq[5*Np+idx]=0.1111111111111111*nB; + Bq[6*Np+idx]=0.1111111111111111*nB; + } +} + +extern "C" void ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int Source, int Dest){ + int n; double value; + for (n=0; n +#include + +ScaLBL_FreeLeeModel::ScaLBL_FreeLeeModel(int RANK, int NP, MPI_Comm COMM): +rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rhoA(0),rhoB(0),W(0),gamma(0), +Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0), +Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM) +{ + +ScaLBL_FreeLeeModel::~ScaLBL_FreeLeeModel(){ + +} +void ScaLBL_FreeLeeModel::ReadParams(string filename){ + // read the input database + db = std::make_shared( filename ); + domain_db = db->getDatabase( "Domain" ); + freelee_db = db->getDatabase( "FreeLee" ); + analysis_db = db->getDatabase( "Analysis" ); + vis_db = db->getDatabase( "Visualization" ); + + // set defaults + timestepMax = 100000; + tauA = tauB = 1.0; + rhoA = rhoB = 1.0; + Fx = Fy = Fz = 0.0; + gamma=1e-3; + W=5; + Restart=false; + din=dout=1.0; + flux=0.0; + + // Color Model parameters + if (freelee_db->keyExists( "timestepMax" )){ + timestepMax = freelee_db->getScalar( "timestepMax" ); + } + if (freelee_db->keyExists( "tauA" )){ + tauA = freelee_db->getScalar( "tauA" ); + } + if (freelee_db->keyExists( "tauB" )){ + tauB = freelee_db->getScalar( "tauB" ); + } + if (freelee_db->keyExists( "rhoA" )){ + rhoA = freelee_db->getScalar( "rhoA" ); + } + if (freelee_db->keyExists( "rhoB" )){ + rhoB = freelee_db->getScalar( "rhoB" ); + } + if (freelee_db->keyExists( "F" )){ + Fx = freelee_db->getVector( "F" )[0]; + Fy = freelee_db->getVector( "F" )[1]; + Fz = freelee_db->getVector( "F" )[2]; + } + if (freelee_db->keyExists( "gamma" )){ + gamma = freelee_db->getScalar( "gamma" ); + } + if (freelee_db->keyExists( "W" )){ + W = freelee_db->getScalar( "W" ); + } + if (freelee_db->keyExists( "Restart" )){ + Restart = freelee_db->getScalar( "Restart" ); + } + if (freelee_db->keyExists( "din" )){ + din = freelee_db->getScalar( "din" ); + } + if (freelee_db->keyExists( "dout" )){ + dout = freelee_db->getScalar( "dout" ); + } + if (freelee_db->keyExists( "flux" )){ + flux = freelee_db->getScalar( "flux" ); + } + inletA=1.f; + inletB=0.f; + outletA=0.f; + outletB=1.f; + //if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details) + + BoundaryCondition = 0; + if (domain_db->keyExists( "BC" )){ + BoundaryCondition = domain_db->getScalar( "BC" ); + } +} + +void ScaLBL_FreeLeeModel::SetDomain(){ + Dm = std::shared_ptr(new Domain(domain_db,comm)); // full domain for analysis + Mask = std::shared_ptr(new Domain(domain_db,comm)); // mask domain removes immobile phases + // domain parameters + Nx = Dm->Nx; + Ny = Dm->Ny; + Nz = Dm->Nz; + Lx = Dm->Lx; + Ly = Dm->Ly; + Lz = Dm->Lz; + N = Nx*Ny*Nz; + Nxh = Nx+2; + Nyh = Ny+2; + Nzh = Nz+2; + Nh = Nxh*Nyh*Nzh; + id = new signed char [N]; + for (int i=0; iid[i] = 1; // initialize this way + //Averages = std::shared_ptr ( new TwoPhase(Dm) ); // TwoPhase analysis object + Averages = std::shared_ptr ( new SubPhase(Dm) ); // TwoPhase analysis object + MPI_Barrier(comm); + Dm->CommInit(); + MPI_Barrier(comm); + // Read domain parameters + rank = Dm->rank(); + nprocx = Dm->nprocx(); + nprocy = Dm->nprocy(); + nprocz = Dm->nprocz(); +} + +void ScaLBL_FreeLeeModel::ReadInput(){ + + sprintf(LocalRankString,"%05d",rank); + sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); + sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString); + + if (freelee_db->keyExists( "image_sequence" )){ + auto ImageList = freelee_db->getVector( "image_sequence"); + int IMAGE_INDEX = freelee_db->getWithDefault( "image_index", 0 ); + std::string first_image = ImageList[IMAGE_INDEX]; + Mask->Decomp(first_image); + IMAGE_INDEX++; + } + else if (domain_db->keyExists( "GridFile" )){ + // Read the local domain data + auto input_id = readMicroCT( *domain_db, MPI_COMM_WORLD ); + // Fill the halo (assuming GCW of 1) + array size0 = { (int) input_id.size(0), (int) input_id.size(1), (int) input_id.size(2) }; + ArraySize size1 = { (size_t) Mask->Nx, (size_t) Mask->Ny, (size_t) Mask->Nz }; + ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 ); + fillHalo fill( MPI_COMM_WORLD, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 ); + Array id_view; + id_view.viewRaw( size1, Mask->id ); + fill.copy( input_id, id_view ); + fill.fill( id_view ); + } + else if (domain_db->keyExists( "Filename" )){ + auto Filename = domain_db->getScalar( "Filename" ); + Mask->Decomp(Filename); + } + else{ + Mask->ReadIDs(); + } + for (int i=0; iid[i]; // save what was read + + // Generate the signed distance map + // Initialize the domain and communication + Array id_solid(Nx,Ny,Nz); + // Solve for the position of the solid phase + for (int k=0;kid[n]; + if (label > 0) id_solid(i,j,k) = 1; + else id_solid(i,j,k) = 0; + } + } + } + // Initialize the signed distance function + for (int k=0;kSDs(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; + } + } + } +// MeanFilter(Averages->SDs); + if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); + CalcDist(Averages->SDs,id_solid,*Mask); + + if (rank == 0) cout << "Domain set." << endl; + + Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta); +} + +void ScaLBL_FreeLeeModel::AssignComponentLabels(double *phase) +{ + size_t NLABELS=0; + signed char VALUE=0; + double AFFINITY=0.f; + + auto LabelList = freelee_db->getVector( "ComponentLabels" ); + auto AffinityList = freelee_db->getVector( "ComponentAffinity" ); + + NLABELS=LabelList.size(); + if (NLABELS != AffinityList.size()){ + ERROR("Error: ComponentLabels and ComponentAffinity must be the same length! \n"); + } + + double label_count[NLABELS]; + double label_count_global[NLABELS]; + // Assign the labels + + for (size_t idx=0; idxid[n] = 0; // set mask to zero since this is an immobile component + } + } + // fluid labels are reserved + if (VALUE == 1) AFFINITY=1.0; + else if (VALUE == 2) AFFINITY=-1.0; + phase[n] = AFFINITY; + } + } + } + + // Set Dm to match Mask + for (int i=0; iid[i] = Mask->id[i]; + + for (size_t idx=0; idxComm, label_count[idx]); + + if (rank==0){ + printf("Component labels: %lu \n",NLABELS); + for (unsigned int idx=0; idxid[i] = Mask->id[i]; + Mask->CommInit(); + Np=Mask->PoreCount(); + //........................................................................... + if (rank==0) printf ("Create ScaLBL_Communicator \n"); + // Create a communicator for the device (will use optimized layout) + // ScaLBL_Communicator ScaLBL_Comm(Mask); // original + ScaLBL_Comm = std::shared_ptr(new ScaLBL_Communicator(Mask)); + ScaLBL_Comm_Regular = std::shared_ptr(new ScaLBL_Communicator(Mask)); + + // create wide halo for phase field + //ScaLBL_Comm_Regular->WideHalo + + // create the layout for the LBM + int Npad=(Np/16 + 2)*16; + if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); + Map.resize(Nx,Ny,Nz); Map.fill(-2); + auto neighborList= new int[18*Npad]; + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + MPI_Barrier(comm); + + //........................................................................... + // MAIN VARIABLES ALLOCATED HERE + //........................................................................... + // LBM variables + if (rank==0) printf ("Allocating distributions \n"); + //......................device distributions................................. + dist_mem_size = Np*sizeof(double); + neighborSize=18*(Np*sizeof(int)); + //........................................................................... + ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); + ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np); + ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &hq, 7*dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &mu_phi, dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &Den, dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nh); + ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np); + //........................................................................... + // Update GPU data structures + if (rank==0) printf ("Setting up device map and neighbor list \n"); + fflush(stdout); + int *TmpMap; + TmpMap=new int[Np]; + for (int k=1; kLastExterior(); idx++){ + auto n = TmpMap[idx]; + if (n > Nx*Ny*Nz){ + printf("Bad value! idx=%i \n", n); + TmpMap[idx] = Nx*Ny*Nz-1; + } + } + for (int idx=ScaLBL_Comm->FirstInterior(); idxLastInterior(); idx++){ + auto n = TmpMap[idx]; + if ( n > Nx*Ny*Nz ){ + printf("Bad value! idx=%i \n",n); + TmpMap[idx] = Nx*Ny*Nz-1; + } + } + ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np); + ScaLBL_DeviceBarrier(); + delete [] TmpMap; + + // copy the neighbor list + ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); + // initialize phi based on PhaseLabel (include solid component labels) + double *PhaseLabel; + PhaseLabel = new double[N]; + AssignComponentLabels(PhaseLabel); + ScaLBL_CopyToDevice(Phi, PhaseLabel, N*sizeof(double)); +} + +/******************************************************** + * AssignComponentLabels * + ********************************************************/ + +void ScaLBL_FreeLeeModel::Initialize(){ + + if (rank==0) printf ("Initializing distributions \n"); + ScaLBL_D3Q19_Init(fq, Np); + /* + * This function initializes model + */ + if (Restart == true){ + if (rank==0){ + printf("Reading restart file! \n"); + } + + // Read in the restart file to CPU buffers + int *TmpMap; + TmpMap = new int[Np]; + + double *cPhi, *cDist, *cDen; + cPhi = new double[N]; + cDen = new double[2*Np]; + cDist = new double[19*Np]; + ScaLBL_CopyToHost(TmpMap, dvcMap, Np*sizeof(int)); + ScaLBL_CopyToHost(cPhi, Phi, N*sizeof(double)); + + ifstream File(LocalRestartFile,ios::binary); + int idx; + double value,va,vb; + for (int n=0; nLastExterior(); n++){ + va = cDen[n]; + vb = cDen[Np + n]; + value = (va-vb)/(va+vb); + idx = TmpMap[n]; + if (!(idx < 0) && idxFirstInterior(); nLastInterior(); n++){ + va = cDen[n]; + vb = cDen[Np + n]; + value = (va-vb)/(va+vb); + idx = TmpMap[n]; + if (!(idx < 0) && idxLastExterior(), Np); + ScaLBL_PhaseField_Init(dvcMap, Phi, Den, hq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + + // establish reservoirs for external bC + if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4 ){ + if (Dm->kproc()==0){ + ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); + ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1); + ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2); + } + if (Dm->kproc() == nprocz-1){ + ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); + ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2); + ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3); + } + } + ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double)); +} + +void ScaLBL_FreeLeeModel::Run(){ + int nprocs=nprocx*nprocy*nprocz; + const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + + if (rank==0){ + printf("********************************************************\n"); + printf("No. of timesteps: %i \n", timestepMax); + fflush(stdout); + } + + //.......create and start timer............ + double starttime,stoptime,cputime; + ScaLBL_DeviceBarrier(); + MPI_Barrier(comm); + starttime = MPI_Wtime(); + //......................................... + + //************ MAIN ITERATION LOOP ***************************************/ + PROFILE_START("Loop"); + while (timestep < timestepMax ) { + //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); } + PROFILE_START("Update"); + // *************ODD TIMESTEP************* + timestep++; + // Compute the Phase indicator field + // Read for hq, Bq happens in this routine (requires communication) + ScaLBL_Comm->BiSendD3Q7AA(hq,Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, hq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(hq,Bq); //WRITE INTO OPPOSITE + ScaLBL_DeviceBarrier(); + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, hq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); + + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + if (BoundaryCondition > 0 && BoundaryCondition < 5){ + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + // Halo exchange for phase field + ScaLBL_Comm_Regular->SendHalo(Phi); + + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_DeviceBarrier(); + // Set BCs + if (BoundaryCondition == 3){ + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + if (BoundaryCondition == 4){ + din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + else if (BoundaryCondition == 5){ + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_DeviceBarrier(); + MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL); + + // *************EVEN TIMESTEP************* + timestep++; + // Compute the Phase indicator field + ScaLBL_Comm->BiSendD3Q7AA(hq,Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, hq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(hq,Bq); //WRITE INTO OPPOSITE + ScaLBL_DeviceBarrier(); + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, hq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); + + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL + // Halo exchange for phase field + if (BoundaryCondition > 0 && BoundaryCondition < 5){ + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + ScaLBL_Comm_Regular->SendHalo(Phi); + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_DeviceBarrier(); + // Set boundary conditions + if (BoundaryCondition == 3){ + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + else if (BoundaryCondition == 4){ + din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + else if (BoundaryCondition == 5){ + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_DeviceBarrier(); + MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL); + //************************************************************************ + PROFILE_STOP("Update"); + } + PROFILE_STOP("Loop"); + PROFILE_SAVE("lbpm_color_simulator",1); + //************************************************************************ + 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(Np)/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"); + + // ************************************************************************ +} + + +void ScaLBL_FreeLeeModel::WriteDebug(){ + // Copy back final phase indicator field and convert to regular layout + DoubleArray PhaseField(Nx,Ny,Nz); + //ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField); + ScaLBL_CopyToHost(PhaseField.data(), Phi, sizeof(double)*N); + + FILE *OUTFILE; + sprintf(LocalRankFilename,"Phase.%05i.raw",rank); + OUTFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,OUTFILE); + fclose(OUTFILE); + + ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); + FILE *AFILE; + sprintf(LocalRankFilename,"A.%05i.raw",rank); + AFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,AFILE); + fclose(AFILE); + + ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField); + FILE *BFILE; + sprintf(LocalRankFilename,"B.%05i.raw",rank); + BFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,BFILE); + fclose(BFILE); + + ScaLBL_Comm->RegularLayout(Map,Pressure,PhaseField); + FILE *PFILE; + sprintf(LocalRankFilename,"Pressure.%05i.raw",rank); + PFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,PFILE); + fclose(PFILE); + + ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField); + FILE *VELX_FILE; + sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank); + VELX_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,VELX_FILE); + fclose(VELX_FILE); + + ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField); + FILE *VELY_FILE; + sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank); + VELY_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,VELY_FILE); + fclose(VELY_FILE); + + ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField); + FILE *VELZ_FILE; + sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank); + VELZ_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,VELZ_FILE); + fclose(VELZ_FILE); + +/* ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField); + FILE *CGX_FILE; + sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank); + CGX_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,CGX_FILE); + fclose(CGX_FILE); + + ScaLBL_Comm->RegularLayout(Map,&ColorGrad[Np],PhaseField); + FILE *CGY_FILE; + sprintf(LocalRankFilename,"Gradient_Y.%05i.raw",rank); + CGY_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,CGY_FILE); + fclose(CGY_FILE); + + ScaLBL_Comm->RegularLayout(Map,&ColorGrad[2*Np],PhaseField); + FILE *CGZ_FILE; + sprintf(LocalRankFilename,"Gradient_Z.%05i.raw",rank); + CGZ_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,CGZ_FILE); + fclose(CGZ_FILE); +*/ +} diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h new file mode 100644 index 00000000..42419afa --- /dev/null +++ b/models/FreeLeeModel.h @@ -0,0 +1,83 @@ +/* +Implementation of color lattice boltzmann model + */ +#include +#include +#include +#include +#include +#include +#include + +#include "common/Communication.h" +#include "common/MPI_Helpers.h" +#include "ProfilerApp.h" +#include "threadpool/thread_pool.h" + +class ScaLBL_FreeLeeModel{ +public: + ScaLBL_FreeLeeModel(int RANK, int NP, MPI_Comm COMM); + ~ScaLBL_FreeLeeModel(); + + // functions in they should be run + void ReadParams(string filename); + void ReadParams(std::shared_ptr db0); + void SetDomain(); + void ReadInput(); + void Create(); + void Initialize(); + void Run(); + void WriteDebug(); + + bool Restart,pBC; + int timestep,timestepMax; + int BoundaryCondition; + double tauA,tauB,rhoA,rhoB; + double W,gamma; + double Fx,Fy,Fz,flux; + double din,dout,inletA,inletB,outletA,outletB; + + int Nx,Ny,Nz,N,Np; + int Nxh,Nyh,Nzh,Nh; // extra halo width + int rank,nprocx,nprocy,nprocz,nprocs; + double Lx,Ly,Lz; + + std::shared_ptr Dm; // this domain is for analysis + std::shared_ptr Mask; // this domain is for lbm + std::shared_ptr ScaLBL_Comm; + std::shared_ptr ScaLBL_Comm_Regular; + //std::shared_ptr Averages; + std::shared_ptr Averages; + + // input database + std::shared_ptr db; + std::shared_ptr domain_db; + std::shared_ptr freelee_db; + std::shared_ptr analysis_db; + std::shared_ptr vis_db; + + IntArray Map; + signed char *id; + int *NeighborList; + int *dvcMap; + double *fq, *hq; + double *mu_phi, *Den, *Phi; + double *ColorGrad; + double *Velocity; + double *Pressure; + +private: + MPI_Comm comm; + + int dist_mem_size; + int neighborSize; + // filenames + char LocalRankString[8]; + char LocalRankFilename[40]; + char LocalRestartFile[40]; + + //int rank,nprocs; + void LoadParams(std::shared_ptr db0); + +}; + From 0b94f29d58dd3ff0d720be5d7f47aa62a2fc9b72 Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 28 Sep 2020 15:44:22 -0400 Subject: [PATCH 02/39] fixed weird type in ScaLBL --- common/ScaLBL.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 2e8eef1a..f76a8de5 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -373,9 +373,9 @@ private: int *dvcRecvList_xy, *dvcRecvList_yz, *dvcRecvList_xz, *dvcRecvList_Xy, *dvcRecvList_Yz, *dvcRecvList_xZ; int *dvcRecvList_xY, *dvcRecvList_yZ, *dvcRecvList_Xz, *dvcRecvList_XY, *dvcRecvList_YZ, *dvcRecvList_XZ; // Recieve buffers for the distributions - int *dvcRecvDist_x, *dvcRecvDist_y, *dvcRecvDist_z, *dvcRecvDist_X, *dvcRecvDist_Y, *dvcRecvDist_Z; - int *dvcRecvDist_xy, *dvcRecvDist_yz, *dvcRecvDist_xz, *dvcRecvDist_Xy, *dvcRecvDist_Yz, *dvcRecvDist_xZ; - int *dvcRecvDist_xY, *dvcRecvDist_yZ, *dvcRecvDist_Xz, *dvcRecvDist_XY, *dvcRecvDist_YZ, *dvcRecvDist_XZ; + double *dvcRecvDist_x, *dvcRecvDist_y, *dvcRecvDist_z, *dvcRecvDist_X, *dvcRecvDist_Y, *dvcRecvDist_Z; + double *dvcRecvDist_xy, *dvcRecvDist_yz, *dvcRecvDist_xz, *dvcRecvDist_Xy, *dvcRecvDist_Yz, *dvcRecvDist_xZ; + double *dvcRecvDist_xY, *dvcRecvDist_yZ, *dvcRecvDist_Xz, *dvcRecvDist_XY, *dvcRecvDist_YZ, *dvcRecvDist_XZ; //...................................................................................... }; From a69504f96a3eff1ae31e32c44443d9a7c6f22e2a Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 09:43:38 -0400 Subject: [PATCH 03/39] fixed weird type in ScaLBL (revert) --- common/ScaLBL.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index f76a8de5..2e8eef1a 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -373,9 +373,9 @@ private: int *dvcRecvList_xy, *dvcRecvList_yz, *dvcRecvList_xz, *dvcRecvList_Xy, *dvcRecvList_Yz, *dvcRecvList_xZ; int *dvcRecvList_xY, *dvcRecvList_yZ, *dvcRecvList_Xz, *dvcRecvList_XY, *dvcRecvList_YZ, *dvcRecvList_XZ; // Recieve buffers for the distributions - double *dvcRecvDist_x, *dvcRecvDist_y, *dvcRecvDist_z, *dvcRecvDist_X, *dvcRecvDist_Y, *dvcRecvDist_Z; - double *dvcRecvDist_xy, *dvcRecvDist_yz, *dvcRecvDist_xz, *dvcRecvDist_Xy, *dvcRecvDist_Yz, *dvcRecvDist_xZ; - double *dvcRecvDist_xY, *dvcRecvDist_yZ, *dvcRecvDist_Xz, *dvcRecvDist_XY, *dvcRecvDist_YZ, *dvcRecvDist_XZ; + int *dvcRecvDist_x, *dvcRecvDist_y, *dvcRecvDist_z, *dvcRecvDist_X, *dvcRecvDist_Y, *dvcRecvDist_Z; + int *dvcRecvDist_xy, *dvcRecvDist_yz, *dvcRecvDist_xz, *dvcRecvDist_Xy, *dvcRecvDist_Yz, *dvcRecvDist_xZ; + int *dvcRecvDist_xY, *dvcRecvDist_yZ, *dvcRecvDist_Xz, *dvcRecvDist_XY, *dvcRecvDist_YZ, *dvcRecvDist_XZ; //...................................................................................... }; From 6b335eaf287c7b1a8a2b9503256edc28330fd023 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 10:39:05 -0400 Subject: [PATCH 04/39] add wide halo class --- common/WideHalo.cpp | 387 ++++++++++++++++++++++++++++++++++++++++++++ common/WideHalo.h | 96 +++++++++++ 2 files changed, 483 insertions(+) create mode 100644 common/WideHalo.cpp create mode 100644 common/WideHalo.h diff --git a/common/WideHalo.cpp b/common/WideHalo.cpp new file mode 100644 index 00000000..2e216c3c --- /dev/null +++ b/common/WideHalo.cpp @@ -0,0 +1,387 @@ +/* +This class implements support for halo widths larger than 1 + */ + +ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr Dm, int width) +{ + //...................................................................................... + Lock=false; // unlock the communicator + //...................................................................................... + // Create a separate copy of the communicator for the device + MPI_Comm_dup(Dm->Comm,&MPI_COMM_SCALBL); + //...................................................................................... + // Copy the domain size and communication information directly from Dm + Nx = Dm->Nx; + Ny = Dm->Ny; + Nz = Dm->Nz; + N = Nx*Ny*Nz; + Nxh = Nx + 2*(width - 1); + Nyh = Ny + 2*(width - 1); + Nzh = Nz + 2*(width - 1); + Nh = Nxh*Nyh*Nzh; + next=0; + + rank=Dm->rank(); + iproc = Dm->iproc(); + jproc = Dm->jproc(); + kproc = Dm->kproc(); + nprocx = Dm->nprocx(); + nprocy = Dm->nprocy(); + nprocz = Dm->nprocz(); + rank_info = RankInfoStruct(myrank,nprocx,nprocy,nprocz); + rank = rank_info.rank[1][1][1]; + rank_X = rank_info.rank[2][1][1]; + rank_x = rank_info.rank[0][1][1]; + rank_Y = rank_info.rank[1][2][1]; + rank_y = rank_info.rank[1][0][1]; + rank_Z = rank_info.rank[1][1][2]; + rank_z = rank_info.rank[1][1][0]; + rank_XY = rank_info.rank[2][2][1]; + rank_xy = rank_info.rank[0][0][1]; + rank_Xy = rank_info.rank[2][0][1]; + rank_xY = rank_info.rank[0][2][1]; + rank_XZ = rank_info.rank[2][1][2]; + rank_xz = rank_info.rank[0][1][0]; + rank_Xz = rank_info.rank[2][1][0]; + rank_xZ = rank_info.rank[0][1][2]; + rank_YZ = rank_info.rank[1][2][2]; + rank_yz = rank_info.rank[1][0][0]; + rank_Yz = rank_info.rank[1][2][0]; + rank_yZ = rank_info.rank[1][0][2]; + rank_XYz = rank_info.rank[2][2][0]; + rank_xyz = rank_info.rank[0][0][0]; + rank_Xyz = rank_info.rank[2][0][0]; + rank_xYz = rank_info.rank[0][2][0]; + rank_XYZ = rank_info.rank[2][2][2]; + rank_xyZ = rank_info.rank[0][0][2]; + rank_XyZ = rank_info.rank[2][0][2]; + rank_xYZ = rank_info.rank[0][2][2]; + + sendCount_x = (Ny-2)*(Nz-2)*width; + sendCount_y = (Nx-2)*(Nz-2)*width; + sendCount_z = (Nx-2)*(Ny-2)*width; + sendCount_X = (Ny-2)*(Nz-2)*width; + sendCount_Y = (Nx-2)*(Nz-2)*width; + sendCount_Z = (Nx-2)*(Ny-2)*width; + sendCount_xy = (Nz-2)*width*width; + sendCount_yz = (Nx-2)*width*width; + sendCount_xz = (Ny-2)*width*width; + sendCount_Xy = (Nz-2)*width*width; + sendCount_Yz = (Nx-2)*width*width; + sendCount_xZ = (Ny-2)*width*width; + sendCount_xY = (Nz-2)*width*width; + sendCount_yZ = (Nx-2)*width*width; + sendCount_Xz = (Ny-2)*width*width; + sendCount_XY = (Nz-2)*width*width; + sendCount_YZ = (Nx-2)*width*width; + sendCount_XZ = (Ny-2)*width*width; + sendCount_xyz = width*width*width; + sendCount_Xyz = width*width*width; + sendCount_xYz = width*width*width; + sendCount_XYz = width*width*width; + sendCount_xyZ = width*width*width; + sendCount_XyZ = width*width*width; + sendCount_xYZ = width*width*width; + sendCount_XYZ = width*width*width; + + RecvCount_x = (Ny-2)*(Nz-2)*width; + RecvCount_y = (Nx-2)*(Nz-2)*width; + RecvCount_z = (Nx-2)*(Ny-2)*width; + RecvCount_X = (Ny-2)*(Nz-2)*width; + RecvCount_Y = (Nx-2)*(Nz-2)*width; + RecvCount_Z = (Nx-2)*(Ny-2)*width; + RecvCount_xy = (Nz-2)*width*width; + RecvCount_yz = (Nx-2)*width*width; + RecvCount_xz = (Ny-2)*width*width; + RecvCount_Xy = (Nz-2)*width*width; + RecvCount_Yz = (Nx-2)*width*width; + RecvCount_xZ = (Ny-2)*width*width; + RecvCount_xY = (Nz-2)*width*width; + RecvCount_yZ = (Nx-2)*width*width; + RecvCount_Xz = (Ny-2)*width*width; + RecvCount_XY = (Nz-2)*width*width; + RecvCount_YZ = (Nx-2)*width*width; + RecvCount_XZ = (Ny-2)*width*width; + RecvCount_xyz = width*width*width; + RecvCount_Xyz = width*width*width; + RecvCount_xYz = width*width*width; + RecvCount_XYz = width*width*width; + RecvCount_xyZ = width*width*width; + RecvCount_XyZ = width*width*width; + RecvCount_xYZ = width*width*width; + RecvCount_XYZ = width*width*width; + //...................................................................................... + ScaLBL_AllocateZeroCopy((void **) &sendbuf_x, sendCount_x*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_X, sendCount_X*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_y, sendCount_y*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Y, sendCount_Y*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_z, sendCount_z*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Z, sendCount_Z*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xy, sendCount_xy*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xY, sendCount_xY*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Xy, sendCount_Xy*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_XY, sendCount_XY*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xz, sendCount_xz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xZ, sendCount_xZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Xz, sendCount_Xz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_XZ, sendCount_XZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_yz, sendCount_yz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_yZ, sendCount_yZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Yz, sendCount_Yz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_YZ, sendCount_YZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xyz, sendCount_xyz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xYz, sendCount_xYz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Xyz, sendCount_Xyz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_XYz, sendCount_XYz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xyZ, sendCount_xyZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xYZ, sendCount_xYZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_XyZ, sendCount_XyZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_XYZ, sendCount_XYZ*sizeof(double)); // Allocate device memory + //...................................................................................... + ScaLBL_AllocateZeroCopy((void **) &recvbuf_x, recvCount_x*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_X, recvCount_X*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_y, recvCount_y*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Y, recvCount_Y*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_z, recvCount_z*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Z, recvCount_Z*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xy, recvCount_xy*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xY, recvCount_xY*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Xy, recvCount_Xy*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_XY, recvCount_XY*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xz, recvCount_xz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xZ, recvCount_xZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Xz, recvCount_Xz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_XZ, recvCount_XZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_yz, recvCount_yz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_yZ, recvCount_yZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Yz, recvCount_Yz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_YZ, recvCount_YZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xyz, recvCount_xyz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xYz, recvCount_xYz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Xyz, recvCount_Xyz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_XYz, recvCount_XYz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xyZ, recvCount_xyZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xYZ, recvCount_xYZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_XyZ, recvCount_XyZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_XYZ, recvCount_XYZ*sizeof(double)); // Allocate device memory + //...................................................................................... + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_x, sendCount_x*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_X, sendCount_X*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_y, sendCount_y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Y, sendCount_Y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_z, sendCount_z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Z, sendCount_Z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xy, sendCount_xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xY, sendCount_xY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Xy, sendCount_Xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_XY, sendCount_XY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xz, sendCount_xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xZ, sendCount_xZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Xz, sendCount_Xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_XZ, sendCount_XZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_yz, sendCount_yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_yZ, sendCount_yZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Yz, sendCount_Yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_YZ, sendCount_YZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xyz, sendCount_xyz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xYz, sendCount_xYz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Xyz, sendCount_Xyz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_XYz, sendCount_XYz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xyZ, sendCount_xyZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xYZ, sendCount_xYZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_XyZ, sendCount_XyZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_XYZ, sendCount_XYZ*sizeof(int)); // Allocate device memory + //...................................................................................... + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_x, recvCount_x*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_X, recvCount_X*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_y, recvCount_y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Y, recvCount_Y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_z, recvCount_z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Z, recvCount_Z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xy, recvCount_xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xY, recvCount_xY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Xy, recvCount_Xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_XY, recvCount_XY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xz, recvCount_xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xZ, recvCount_xZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Xz, recvCount_Xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_XZ, recvCount_XZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_yz, recvCount_yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_yZ, recvCount_yZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Yz, recvCount_Yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_YZ, recvCount_YZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xyz, recvCount_xyZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xYz, recvCount_xYZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Xyz, recvCount_XyZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_XYz, recvCount_XYZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xyZ, recvCount_xyZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xYZ, recvCount_xYZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_XyZ, recvCount_XyZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_XYZ, recvCount_XYZ*sizeof(int)); // Allocate device memory + //...................................................................................... + + MPI_Barrier(MPI_COMM_SCALBL); + + /* Fill in communications patterns for the lists */ + int *sendList, *recvList; + sendList = new int [width*max(Nxh,Nyh)*max(Nxh,Nzh)]; + /* x face */ + int count = 0; + for (int k=0; k Dm, int width); + ~ScaLBLWideHalo_Communicator(); + //...................................................................................... + MPI_Comm MPI_COMM_SCALBL; // MPI Communicator + unsigned long int CommunicationCount,SendCount,RecvCount; + int Nx,Ny,Nz,N; // original domain structure + int Nxh,Nyh,Nzh,Nh; // with wide halo + RankInfoStruct rank_info; + + int first_interior,last_interior; + //...................................................................................... + // Set up for D3Q19 distributions -- all 27 neighbors are needed + //...................................................................................... + // Buffers to store data sent and recieved by this MPI process + double *sendbuf_x, *sendbuf_y, *sendbuf_z, *sendbuf_X, *sendbuf_Y, *sendbuf_Z; + double *sendbuf_xy, *sendbuf_yz, *sendbuf_xz, *sendbuf_Xy, *sendbuf_Yz, *sendbuf_xZ; + double *sendbuf_xY, *sendbuf_yZ, *sendbuf_Xz, *sendbuf_XY, *sendbuf_YZ, *sendbuf_XZ; + double *sendbuf_xyz, *sendbuf_Xyz, *sendbuf_xYz, *sendbuf_XYy; + double *sendbuf_xyZ, *sendbuf_XyZ, *sendbuf_xYZ, *sendbuf_XYZ; + double *recvbuf_x, *recvbuf_y, *recvbuf_z, *recvbuf_X, *recvbuf_Y, *recvbuf_Z; + double *recvbuf_xy, *recvbuf_yz, *recvbuf_xz, *recvbuf_Xy, *recvbuf_Yz, *recvbuf_xZ; + double *recvbuf_xY, *recvbuf_yZ, *recvbuf_Xz, *recvbuf_XY, *recvbuf_YZ, *recvbuf_XZ; + double *recvbuf_xyz, *recvbuf_Xyz, *recvbuf_xYz, *recvbuf_XYy; + double *recvbuf_xyZ, *recvbuf_XyZ, *recvbuf_xYZ, *recvbuf_XYZ; + //...................................................................................... + + int LastExterior(); + int FirstInterior(); + int LastInterior(); + + void Send(double *data); + void Recv(double *data); + + // Debugging and unit testing functions + void PrintDebug()); + +private: + bool Lock; // use Lock to make sure only one call at a time to protect data in transit + // only one set of Send requests can be active at any time (per instance) + int i,j,k,n; + int iproc,jproc,kproc; + int nprocx,nprocy,nprocz; + int sendtag,recvtag; + // Give the object it's own MPI communicator + RankInfoStruct rank_info; + MPI_Group Group; // Group of processors associated with this domain + MPI_Request req1[26],req2[26]; + MPI_Status stat1[26],stat2[26]; + //...................................................................................... + // MPI ranks for all 18 neighbors + //...................................................................................... + // These variables are all private to prevent external things from modifying them!! + //...................................................................................... + int rank; + int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z; + int rank_xy,rank_XY,rank_xY,rank_Xy; + int rank_xz,rank_XZ,rank_xZ,rank_Xz; + int rank_yz,rank_YZ,rank_yZ,rank_Yz; + int rank_xyz,rank_Xyz,rank_xYz,rank_XYz; + int rank_xyZ,rank_XyZ,rank_xYZ,rank_XYZ; + //...................................................................................... + //...................................................................................... + int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z; + int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ; + int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ; + int sendCount_xyz,sendCount_Xyz,sendCount_xYz,sendCount_XYz; + int sendCount_xyZ,sendCount_XyZ,sendCount_xYZ,sendCount_XYZ; + //...................................................................................... + int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z; + int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ; + int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ; + int recvCount_xyz,recvCount_Xyz,recvCount_xYz,recvCount_XYz; + int recvCount_xyZ,recvCount_XyZ,recvCount_xYZ,recvCount_XYZ; + //...................................................................................... + // Send buffers that reside on the compute device + int *dvcSendList_x, *dvcSendList_y, *dvcSendList_z, *dvcSendList_X, *dvcSendList_Y, *dvcSendList_Z; + int *dvcSendList_xy, *dvcSendList_yz, *dvcSendList_xz, *dvcSendList_Xy, *dvcSendList_Yz, *dvcSendList_xZ; + int *dvcSendList_xY, *dvcSendList_yZ, *dvcSendList_Xz, *dvcSendList_XY, *dvcSendList_YZ, *dvcSendList_XZ; + int *dvcSendList_xyz,*dvcSendList_Xyz,*dvcSendList_xYz,*dvcSendList_XYz; + int *dvcSendList_xyZ,*dvcSendList_XyZ,*dvcSendList_xYZ,*dvcSendList_XYZ; + // Recieve buffers that reside on the compute device + int *dvcRecvList_x, *dvcRecvList_y, *dvcRecvList_z, *dvcRecvList_X, *dvcRecvList_Y, *dvcRecvList_Z; + int *dvcRecvList_xy, *dvcRecvList_yz, *dvcRecvList_xz, *dvcRecvList_Xy, *dvcRecvList_Yz, *dvcRecvList_xZ; + int *dvcRecvList_xY, *dvcRecvList_yZ, *dvcRecvList_Xz, *dvcRecvList_XY, *dvcRecvList_YZ, *dvcRecvList_XZ; + int *dvcRecvList_xyz,*dvcRecvList_Xyz,*dvcRecvList_xYz,*dvcRecvList_XYz; + int *dvcRecvList_xyZ,*dvcRecvList_XyZ,*dvcRecvList_xYZ,*dvcRecvList_XYZ; + //...................................................................................... + +}; From 23f6c089f0621426fec30a07de19713cd78502dd Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 10:43:31 -0400 Subject: [PATCH 05/39] fix name of function --- common/WideHalo.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/common/WideHalo.cpp b/common/WideHalo.cpp index 2e216c3c..78753d33 100644 --- a/common/WideHalo.cpp +++ b/common/WideHalo.cpp @@ -277,7 +277,7 @@ ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr Date: Tue, 29 Sep 2020 10:48:54 -0400 Subject: [PATCH 06/39] fixing compile bugs --- common/WideHalo.cpp | 3 ++- common/WideHalo.h | 5 ++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/common/WideHalo.cpp b/common/WideHalo.cpp index 78753d33..0a52bbd8 100644 --- a/common/WideHalo.cpp +++ b/common/WideHalo.cpp @@ -1,6 +1,7 @@ /* This class implements support for halo widths larger than 1 */ +#include "common/WideHalo.h" ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr Dm, int width) { @@ -383,5 +384,5 @@ void ScaLBLWideHalo_Communicator::Recv(double *data){ } inline int getHaloBlock(){ -} + } diff --git a/common/WideHalo.h b/common/WideHalo.h index 11252426..ac9f8833 100644 --- a/common/WideHalo.h +++ b/common/WideHalo.h @@ -1,6 +1,9 @@ /* This class implements support for halo widths larger than 1 */ +#ifndef WideHalo_H +#define WideHalo_H +#include "common/Domain.h" class ScaLBLWideHalo_Communicator{ public: @@ -92,5 +95,5 @@ private: int *dvcRecvList_xyz,*dvcRecvList_Xyz,*dvcRecvList_xYz,*dvcRecvList_XYz; int *dvcRecvList_xyZ,*dvcRecvList_XyZ,*dvcRecvList_xYZ,*dvcRecvList_XYZ; //...................................................................................... - }; +#endif From afa9dac22adfffc634d95dae2cfdef637e1be4a5 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 10:56:01 -0400 Subject: [PATCH 07/39] fixing compile bugs --- common/WideHalo.cpp | 76 ++++++++++++++++++++++----------------------- common/WideHalo.h | 3 +- 2 files changed, 39 insertions(+), 40 deletions(-) diff --git a/common/WideHalo.cpp b/common/WideHalo.cpp index 0a52bbd8..32401ef9 100644 --- a/common/WideHalo.cpp +++ b/common/WideHalo.cpp @@ -20,7 +20,6 @@ ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr rank(); iproc = Dm->iproc(); @@ -29,7 +28,7 @@ ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr nprocx(); nprocy = Dm->nprocy(); nprocz = Dm->nprocz(); - rank_info = RankInfoStruct(myrank,nprocx,nprocy,nprocz); + rank_info = RankInfoStruct(rank,nprocx,nprocy,nprocz); rank = rank_info.rank[1][1][1]; rank_X = rank_info.rank[2][1][1]; rank_x = rank_info.rank[0][1][1]; @@ -85,32 +84,32 @@ ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr Date: Tue, 29 Sep 2020 10:58:33 -0400 Subject: [PATCH 08/39] fixing compile bugs --- common/WideHalo.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/WideHalo.h b/common/WideHalo.h index beeca65c..9ca88271 100644 --- a/common/WideHalo.h +++ b/common/WideHalo.h @@ -3,7 +3,7 @@ This class implements support for halo widths larger than 1 */ #ifndef WideHalo_H #define WideHalo_H -#include "common/Domain.h" +#include "common/ScaLBL.h" class ScaLBLWideHalo_Communicator{ public: From 229accee9c381a1bfd966c37064fdf4e449c732b Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 11:00:42 -0400 Subject: [PATCH 09/39] fixing compile bugs --- common/WideHalo.cpp | 3 ++- common/WideHalo.h | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/common/WideHalo.cpp b/common/WideHalo.cpp index 32401ef9..9582b601 100644 --- a/common/WideHalo.cpp +++ b/common/WideHalo.cpp @@ -384,5 +384,6 @@ void ScaLBLWideHalo_Communicator::Recv(double *data){ } inline int getHaloBlock(){ - + int count = 0; + return count; } diff --git a/common/WideHalo.h b/common/WideHalo.h index 9ca88271..2dbe8cae 100644 --- a/common/WideHalo.h +++ b/common/WideHalo.h @@ -24,12 +24,12 @@ public: double *sendbuf_x, *sendbuf_y, *sendbuf_z, *sendbuf_X, *sendbuf_Y, *sendbuf_Z; double *sendbuf_xy, *sendbuf_yz, *sendbuf_xz, *sendbuf_Xy, *sendbuf_Yz, *sendbuf_xZ; double *sendbuf_xY, *sendbuf_yZ, *sendbuf_Xz, *sendbuf_XY, *sendbuf_YZ, *sendbuf_XZ; - double *sendbuf_xyz, *sendbuf_Xyz, *sendbuf_xYz, *sendbuf_XYy; + double *sendbuf_xyz, *sendbuf_Xyz, *sendbuf_xYz, *sendbuf_XYz; double *sendbuf_xyZ, *sendbuf_XyZ, *sendbuf_xYZ, *sendbuf_XYZ; double *recvbuf_x, *recvbuf_y, *recvbuf_z, *recvbuf_X, *recvbuf_Y, *recvbuf_Z; double *recvbuf_xy, *recvbuf_yz, *recvbuf_xz, *recvbuf_Xy, *recvbuf_Yz, *recvbuf_xZ; double *recvbuf_xY, *recvbuf_yZ, *recvbuf_Xz, *recvbuf_XY, *recvbuf_YZ, *recvbuf_XZ; - double *recvbuf_xyz, *recvbuf_Xyz, *recvbuf_xYz, *recvbuf_XYy; + double *recvbuf_xyz, *recvbuf_Xyz, *recvbuf_xYz, *recvbuf_XYz; double *recvbuf_xyZ, *recvbuf_XyZ, *recvbuf_xYZ, *recvbuf_XYZ; //...................................................................................... From e4bdc864c6b7b05f02e440e0c7db761d09780242 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 11:05:09 -0400 Subject: [PATCH 10/39] fixing compile bugs --- models/FreeLeeModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp index 7fc005e3..f3e5b42c 100644 --- a/models/FreeLeeModel.cpp +++ b/models/FreeLeeModel.cpp @@ -1,7 +1,7 @@ /* color lattice boltzmann model */ -#include "models/ColorModel.h" +#include "models/FreeLeeModel.h" #include "analysis/distance.h" #include "analysis/morphology.h" #include "common/Communication.h" From a3aadd947efaa659f24c8e319b7c83a6ae811356 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 13:25:53 -0400 Subject: [PATCH 11/39] fixing compile bugs --- models/FreeLeeModel.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h index 42419afa..7409946d 100644 --- a/models/FreeLeeModel.h +++ b/models/FreeLeeModel.h @@ -1,5 +1,5 @@ /* -Implementation of color lattice boltzmann model +Implementation of Lee et al JCP 2016 lattice boltzmann model */ #include #include From cf1cf4eb53e863d13223e2c8da38dfb9430dcffc Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 13:29:14 -0400 Subject: [PATCH 12/39] fixing compile bugs --- models/FreeLeeModel.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h index 7409946d..0d0e63c3 100644 --- a/models/FreeLeeModel.h +++ b/models/FreeLeeModel.h @@ -13,6 +13,8 @@ Implementation of Lee et al JCP 2016 lattice boltzmann model #include "common/MPI_Helpers.h" #include "ProfilerApp.h" #include "threadpool/thread_pool.h" +#include "common/ScaLBL.h" +#include "common/WideHalo.h" class ScaLBL_FreeLeeModel{ public: From 917be9f6c4ebe75286c92877c98d04cf98147972 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 13:37:02 -0400 Subject: [PATCH 13/39] fixing compile bugs --- models/FreeLeeModel.cpp | 12 ++++++------ models/FreeLeeModel.h | 2 ++ 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp index f3e5b42c..e9bca7be 100644 --- a/models/FreeLeeModel.cpp +++ b/models/FreeLeeModel.cpp @@ -14,7 +14,8 @@ rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rh Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0), Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM) { - + +} ScaLBL_FreeLeeModel::~ScaLBL_FreeLeeModel(){ } @@ -105,8 +106,7 @@ void ScaLBL_FreeLeeModel::SetDomain(){ Nh = Nxh*Nyh*Nzh; id = new signed char [N]; for (int i=0; iid[i] = 1; // initialize this way - //Averages = std::shared_ptr ( new TwoPhase(Dm) ); // TwoPhase analysis object - Averages = std::shared_ptr ( new SubPhase(Dm) ); // TwoPhase analysis object + MPI_Barrier(comm); Dm->CommInit(); MPI_Barrier(comm); @@ -167,18 +167,18 @@ void ScaLBL_FreeLeeModel::ReadInput(){ } } } + SignDist.resize(Nx,Ny,Nz); // Initialize the signed distance function for (int k=0;kSDs(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; + SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; } } } -// MeanFilter(Averages->SDs); if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); - CalcDist(Averages->SDs,id_solid,*Mask); + CalcDist(SignDist,id_solid,*Mask); if (rank == 0) cout << "Domain set." << endl; diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h index 0d0e63c3..837ee6a4 100644 --- a/models/FreeLeeModel.h +++ b/models/FreeLeeModel.h @@ -67,6 +67,8 @@ public: double *ColorGrad; double *Velocity; double *Pressure; + + DoubleArray SignDistance; private: MPI_Comm comm; From 5c723742789eb323f64286bd6d1ecefc3445814e Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 13:40:41 -0400 Subject: [PATCH 14/39] fixing compile bugs --- models/FreeLeeModel.cpp | 3 +-- models/FreeLeeModel.h | 6 ++---- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp index e9bca7be..3265411e 100644 --- a/models/FreeLeeModel.cpp +++ b/models/FreeLeeModel.cpp @@ -181,8 +181,7 @@ void ScaLBL_FreeLeeModel::ReadInput(){ CalcDist(SignDist,id_solid,*Mask); if (rank == 0) cout << "Domain set." << endl; - - Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta); + } void ScaLBL_FreeLeeModel::AssignComponentLabels(double *phase) diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h index 837ee6a4..d89fc1c5 100644 --- a/models/FreeLeeModel.h +++ b/models/FreeLeeModel.h @@ -40,7 +40,7 @@ public: double din,dout,inletA,inletB,outletA,outletB; int Nx,Ny,Nz,N,Np; - int Nxh,Nyh,Nzh,Nh; // extra halo width + int Nxh,Nyh,Nzh,Nh; // extra halo width int rank,nprocx,nprocy,nprocz,nprocs; double Lx,Ly,Lz; @@ -48,8 +48,6 @@ public: std::shared_ptr Mask; // this domain is for lbm std::shared_ptr ScaLBL_Comm; std::shared_ptr ScaLBL_Comm_Regular; - //std::shared_ptr Averages; - std::shared_ptr Averages; // input database std::shared_ptr db; @@ -68,7 +66,7 @@ public: double *Velocity; double *Pressure; - DoubleArray SignDistance; + DoubleArray SignDist; private: MPI_Comm comm; From e6e2c4e27f9c5d3c1883164d179bb686e9dce961 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Tue, 29 Sep 2020 13:51:23 -0400 Subject: [PATCH 15/39] FreeLee model compiles --- models/FreeLeeModel.cpp | 75 +++-------------------------------------- 1 file changed, 5 insertions(+), 70 deletions(-) diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp index 3265411e..f52802d3 100644 --- a/models/FreeLeeModel.cpp +++ b/models/FreeLeeModel.cpp @@ -184,68 +184,6 @@ void ScaLBL_FreeLeeModel::ReadInput(){ } -void ScaLBL_FreeLeeModel::AssignComponentLabels(double *phase) -{ - size_t NLABELS=0; - signed char VALUE=0; - double AFFINITY=0.f; - - auto LabelList = freelee_db->getVector( "ComponentLabels" ); - auto AffinityList = freelee_db->getVector( "ComponentAffinity" ); - - NLABELS=LabelList.size(); - if (NLABELS != AffinityList.size()){ - ERROR("Error: ComponentLabels and ComponentAffinity must be the same length! \n"); - } - - double label_count[NLABELS]; - double label_count_global[NLABELS]; - // Assign the labels - - for (size_t idx=0; idxid[n] = 0; // set mask to zero since this is an immobile component - } - } - // fluid labels are reserved - if (VALUE == 1) AFFINITY=1.0; - else if (VALUE == 2) AFFINITY=-1.0; - phase[n] = AFFINITY; - } - } - } - - // Set Dm to match Mask - for (int i=0; iid[i] = Mask->id[i]; - - for (size_t idx=0; idxComm, label_count[idx]); - - if (rank==0){ - printf("Component labels: %lu \n",NLABELS); - for (unsigned int idx=0; idxLastExterior(), Np); - ScaLBL_PhaseField_Init(dvcMap, Phi, Den, hq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + //ScaLBL_PhaseField_Init(dvcMap, Phi, Den, hq, Bq, 0, ScaLBL_Comm->LastExterior(), Np); + //ScaLBL_PhaseField_Init(dvcMap, Phi, Den, hq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); // establish reservoirs for external bC if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4 ){ @@ -423,7 +357,7 @@ void ScaLBL_FreeLeeModel::Initialize(){ ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3); } } - ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double)); + //ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double)); } void ScaLBL_FreeLeeModel::Run(){ @@ -450,7 +384,7 @@ void ScaLBL_FreeLeeModel::Run(){ PROFILE_START("Update"); // *************ODD TIMESTEP************* timestep++; - // Compute the Phase indicator field + /* // Compute the Phase indicator field // Read for hq, Bq happens in this routine (requires communication) ScaLBL_Comm->BiSendD3Q7AA(hq,Bq); //READ FROM NORMAL ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, hq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -527,6 +461,7 @@ void ScaLBL_FreeLeeModel::Run(){ } ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); + */ ScaLBL_DeviceBarrier(); MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL); //************************************************************************ From 3126318900431a2ef658ef59855291aa92e6d128 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Sep 2020 16:35:09 -0400 Subject: [PATCH 16/39] building wide halo class --- common/WideHalo.cpp | 339 ++++++++++++++++++++++------------------ common/WideHalo.h | 17 ++ models/FreeLeeModel.cpp | 4 +- models/FreeLeeModel.h | 3 +- 4 files changed, 211 insertions(+), 152 deletions(-) diff --git a/common/WideHalo.cpp b/common/WideHalo.cpp index 9582b601..5d375b7f 100644 --- a/common/WideHalo.cpp +++ b/common/WideHalo.cpp @@ -56,7 +56,7 @@ ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr (new ScaLBL_Communicator(Mask)); ScaLBL_Comm_Regular = std::shared_ptr(new ScaLBL_Communicator(Mask)); - - // create wide halo for phase field - //ScaLBL_Comm_Regular->WideHalo + ScaLBL_Comm_WideHalo = std::shared_ptr(new ScaLBLWideHalo_Communicator(Mask,2)); // create the layout for the LBM int Npad=(Np/16 + 2)*16; diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h index d89fc1c5..01fb54c3 100644 --- a/models/FreeLeeModel.h +++ b/models/FreeLeeModel.h @@ -48,7 +48,8 @@ public: std::shared_ptr Mask; // this domain is for lbm std::shared_ptr ScaLBL_Comm; std::shared_ptr ScaLBL_Comm_Regular; - + std::shared_ptr ScaLBL_Comm_WideHalo; + // input database std::shared_ptr db; std::shared_ptr domain_db; From e69f2db7774f7c631c39376c02ee6fc88c35cdd7 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 30 Sep 2020 09:06:15 -0400 Subject: [PATCH 17/39] add map to widehalo --- common/Domain.cpp | 2 +- common/WideHalo.cpp | 128 ++++++-------------------------------------- common/WideHalo.h | 2 +- 3 files changed, 17 insertions(+), 115 deletions(-) diff --git a/common/Domain.cpp b/common/Domain.cpp index 7a1a6230..fb1e1c50 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -1,5 +1,5 @@ // Created by James McClure -// Copyright 2008-2013 +// Copyright 2008-2020 #include #include #include diff --git a/common/WideHalo.cpp b/common/WideHalo.cpp index 5d375b7f..48f83cd9 100644 --- a/common/WideHalo.cpp +++ b/common/WideHalo.cpp @@ -21,6 +21,8 @@ ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr rank(); iproc = Dm->iproc(); jproc = Dm->jproc(); @@ -56,118 +58,7 @@ ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr Date: Wed, 30 Sep 2020 11:45:17 -0400 Subject: [PATCH 18/39] add test color grad --- tests/CMakeLists.txt | 3 ++- tests/TestColorGrad.cpp | 15 +-------------- 2 files changed, 3 insertions(+), 15 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c536a7ec..01efe997 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -48,7 +48,8 @@ ADD_LBPM_TEST( TestTopo3D ) ADD_LBPM_TEST( TestFluxBC ) ADD_LBPM_TEST( TestMap ) #ADD_LBPM_TEST( TestMRT ) -#ADD_LBPM_TEST( TestColorGrad ) +ADD_LBPM_TEST( TestColorGrad ) +ADD_LBPM_TEST( TestWideHalo ) #ADD_LBPM_TEST( TestColorGradDFH ) ADD_LBPM_TEST( TestColorGradDFH ) ADD_LBPM_TEST( TestBubbleDFH ../example/Bubble/input.db) diff --git a/tests/TestColorGrad.cpp b/tests/TestColorGrad.cpp index 5cd6d924..33c9f3ce 100644 --- a/tests/TestColorGrad.cpp +++ b/tests/TestColorGrad.cpp @@ -72,20 +72,7 @@ int main(int argc, char **argv) //....................................................................... // Reading the domain information file //....................................................................... - ifstream domain("Domain.in"); - if (domain.good()){ - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - } - else if (nprocs==1){ + if (nprocs==1){ nprocx=nprocy=nprocz=1; Nx=Ny=Nz=3; nspheres=0; From 8d8a22a374cdd5ce72d69060aeb63bc3e3636102 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 30 Sep 2020 11:46:23 -0400 Subject: [PATCH 19/39] add test wide halo --- tests/TestWideHalo.cpp | 269 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 269 insertions(+) create mode 100644 tests/TestWideHalo.cpp diff --git a/tests/TestWideHalo.cpp b/tests/TestWideHalo.cpp new file mode 100644 index 00000000..5cd6d924 --- /dev/null +++ b/tests/TestWideHalo.cpp @@ -0,0 +1,269 @@ + +//************************************************************************* +// Lattice Boltzmann Simulator for Single Phase Flow in Porous Media +// James E. McCLure +//************************************************************************* +#include +#include +#include +#include "common/ScaLBL.h" +#include "common/MPI_Helpers.h" + +using namespace std; + + +//*************************************************************************************** +int main(int argc, char **argv) +{ + //***************************************** + // ***** MPI STUFF **************** + //***************************************** + // Initialize MPI + int rank,nprocs; + MPI_Init(&argc,&argv); + MPI_Comm comm = MPI_COMM_WORLD; + MPI_Comm_rank(comm,&rank); + MPI_Comm_size(comm,&nprocs); + int check; + { + // parallel domain size (# of sub-domains) + int nprocx,nprocy,nprocz; + int iproc,jproc,kproc; + + + if (rank == 0){ + printf("********************************************************\n"); + printf("Running Color Model: TestColor \n"); + printf("********************************************************\n"); + } + + // BGK Model parameters + string FILENAME; + unsigned int nBlocks, nthreads; + int timestepMax, interval; + double Fx,Fy,Fz,tol; + // Domain variables + double Lx,Ly,Lz; + int nspheres; + int Nx,Ny,Nz; + int i,j,k,n; + int dim = 3; + //if (rank == 0) printf("dim=%d\n",dim); + int timestep = 0; + int timesteps = 100; + int centralNode = 2; + + double tauA = 1.0; + double tauB = 1.0; + double rhoA = 1.0; + double rhoB = 1.0; + double alpha = 0.005; + double beta = 0.95; + + double tau = 1.0; + double mu=(tau-0.5)/3.0; + double rlx_setA=1.0/tau; + double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); + + Fx = Fy = 0.f; + Fz = 0.f; + + if (rank==0){ + //....................................................................... + // Reading the domain information file + //....................................................................... + ifstream domain("Domain.in"); + if (domain.good()){ + domain >> nprocx; + domain >> nprocy; + domain >> nprocz; + domain >> Nx; + domain >> Ny; + domain >> Nz; + domain >> nspheres; + domain >> Lx; + domain >> Ly; + domain >> Lz; + } + else if (nprocs==1){ + nprocx=nprocy=nprocz=1; + Nx=Ny=Nz=3; + nspheres=0; + Lx=Ly=Lz=1; + } + else if (nprocs==2){ + nprocx=2; nprocy=1; + nprocz=1; + Nx=Ny=Nz=dim; + Nx = dim; Ny = dim; Nz = dim; + nspheres=0; + Lx=Ly=Lz=1; + } + else if (nprocs==4){ + nprocx=nprocy=2; + nprocz=1; + Nx=Ny=Nz=dim; + nspheres=0; + Lx=Ly=Lz=1; + } + else if (nprocs==8){ + nprocx=nprocy=nprocz=2; + Nx=Ny=Nz=dim; + nspheres=0; + Lx=Ly=Lz=1; + } + //....................................................................... + } + // ************************************************************** + // Broadcast simulation parameters from rank 0 to all other procs + MPI_Barrier(comm); + //................................................. + MPI_Bcast(&Nx,1,MPI_INT,0,comm); + MPI_Bcast(&Ny,1,MPI_INT,0,comm); + MPI_Bcast(&Nz,1,MPI_INT,0,comm); + MPI_Bcast(&nprocx,1,MPI_INT,0,comm); + MPI_Bcast(&nprocy,1,MPI_INT,0,comm); + MPI_Bcast(&nprocz,1,MPI_INT,0,comm); + MPI_Bcast(&nspheres,1,MPI_INT,0,comm); + MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); + //................................................. + MPI_Barrier(comm); + // ************************************************************** + // ************************************************************** + + if (nprocs != nprocx*nprocy*nprocz){ + printf("nprocx = %i \n",nprocx); + printf("nprocy = %i \n",nprocy); + printf("nprocz = %i \n",nprocz); + INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); + } + + if (rank==0){ + printf("********************************************************\n"); + printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz); + printf("********************************************************\n"); + } + + MPI_Barrier(comm); + + double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; + int BoundaryCondition=0; + + Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + + Nx += 2; + Ny += 2; + Nz += 2; + int N = Nx*Ny*Nz; + + int Np=0; // number of local pore nodes + double *PhaseLabel; + PhaseLabel = new double[N]; + //....................................................................... + for (k=0;k 0){ + int idx = Map(i,j,k); + CX=COLORGRAD[idx]; + CY=COLORGRAD[Np+idx]; + CZ=COLORGRAD[2*Np+idx]; + double error=sqrt((CX-1.0)*(CX-1.0)+(CY-2.0)*(CY-2.0)+ (CZ-3.0)*(CZ-3.0)); + if (error > 1e-8) + printf("i,j,k=%i,%i,%i: Color gradient=%f,%f,%f \n",i,j,k,CX,CY,CZ); + } + } + } + } + + } + // **************************************************** + MPI_Barrier(comm); + MPI_Finalize(); + // **************************************************** + + return check; +} + From fed2fc43042ba3301cbe9fdddf8560c6af1d2d23 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 30 Sep 2020 11:53:28 -0400 Subject: [PATCH 20/39] update pointer in tests --- tests/TestColorGrad.cpp | 3 +-- tests/TestWideHalo.cpp | 18 ++---------------- 2 files changed, 3 insertions(+), 18 deletions(-) diff --git a/tests/TestColorGrad.cpp b/tests/TestColorGrad.cpp index 33c9f3ce..58531601 100644 --- a/tests/TestColorGrad.cpp +++ b/tests/TestColorGrad.cpp @@ -138,8 +138,7 @@ int main(int argc, char **argv) double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; int BoundaryCondition=0; - Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); - + Dm = std::shared_ptr(Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); Nx += 2; Ny += 2; Nz += 2; diff --git a/tests/TestWideHalo.cpp b/tests/TestWideHalo.cpp index 5cd6d924..58531601 100644 --- a/tests/TestWideHalo.cpp +++ b/tests/TestWideHalo.cpp @@ -72,20 +72,7 @@ int main(int argc, char **argv) //....................................................................... // Reading the domain information file //....................................................................... - ifstream domain("Domain.in"); - if (domain.good()){ - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - } - else if (nprocs==1){ + if (nprocs==1){ nprocx=nprocy=nprocz=1; Nx=Ny=Nz=3; nspheres=0; @@ -151,8 +138,7 @@ int main(int argc, char **argv) double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; int BoundaryCondition=0; - Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); - + Dm = std::shared_ptr(Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); Nx += 2; Ny += 2; Nz += 2; From e54ab8e195485c4339ce3c0e5bf618dc2d47e845 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 30 Sep 2020 11:54:44 -0400 Subject: [PATCH 21/39] update pointer in tests --- tests/TestColorGrad.cpp | 2 +- tests/TestWideHalo.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/TestColorGrad.cpp b/tests/TestColorGrad.cpp index 58531601..6e584bf1 100644 --- a/tests/TestColorGrad.cpp +++ b/tests/TestColorGrad.cpp @@ -138,7 +138,7 @@ int main(int argc, char **argv) double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; int BoundaryCondition=0; - Dm = std::shared_ptr(Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); Nx += 2; Ny += 2; Nz += 2; diff --git a/tests/TestWideHalo.cpp b/tests/TestWideHalo.cpp index 58531601..6e584bf1 100644 --- a/tests/TestWideHalo.cpp +++ b/tests/TestWideHalo.cpp @@ -138,7 +138,7 @@ int main(int argc, char **argv) double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; int BoundaryCondition=0; - Dm = std::shared_ptr(Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); Nx += 2; Ny += 2; Nz += 2; From 9f1e170a476c0a8a4f778c0fa32cae0138926e81 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 30 Sep 2020 11:57:25 -0400 Subject: [PATCH 22/39] update pointer in tests --- tests/TestColorGrad.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/TestColorGrad.cpp b/tests/TestColorGrad.cpp index 6e584bf1..aa02cc37 100644 --- a/tests/TestColorGrad.cpp +++ b/tests/TestColorGrad.cpp @@ -138,7 +138,7 @@ int main(int argc, char **argv) double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; int BoundaryCondition=0; - Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + std::shared_ptr Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); Nx += 2; Ny += 2; Nz += 2; From 7c16a76d1c10e0e342679208ed3bca557e5c2354 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 30 Sep 2020 12:00:58 -0400 Subject: [PATCH 23/39] update pointer in tests --- tests/TestColorGrad.cpp | 10 +++++----- tests/TestWideHalo.cpp | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/TestColorGrad.cpp b/tests/TestColorGrad.cpp index aa02cc37..d2f1d6de 100644 --- a/tests/TestColorGrad.cpp +++ b/tests/TestColorGrad.cpp @@ -138,7 +138,7 @@ int main(int argc, char **argv) double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; int BoundaryCondition=0; - std::shared_ptr Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + std::shared_ptr Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition)); Nx += 2; Ny += 2; Nz += 2; @@ -152,7 +152,7 @@ int main(int argc, char **argv) for (j=0;jid[n]=1; Np++; // Initialize gradient ColorGrad = (1,2,3) double value=double(3*k+2*j+i); @@ -160,7 +160,7 @@ int main(int argc, char **argv) } } } - Dm.CommInit(); + Dm->CommInit(); MPI_Barrier(comm); if (rank == 0) cout << "Domain set." << endl; if (rank==0) printf ("Create ScaLBL_Communicator \n"); @@ -177,7 +177,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Np]; - ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np); + ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); MPI_Barrier(comm); //......................device distributions................................. @@ -231,7 +231,7 @@ int main(int argc, char **argv) for (j=1;j 0){ + if (Dm->id[n] > 0){ int idx = Map(i,j,k); CX=COLORGRAD[idx]; CY=COLORGRAD[Np+idx]; diff --git a/tests/TestWideHalo.cpp b/tests/TestWideHalo.cpp index 6e584bf1..d2f1d6de 100644 --- a/tests/TestWideHalo.cpp +++ b/tests/TestWideHalo.cpp @@ -138,7 +138,7 @@ int main(int argc, char **argv) double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; int BoundaryCondition=0; - Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + std::shared_ptr Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition)); Nx += 2; Ny += 2; Nz += 2; @@ -152,7 +152,7 @@ int main(int argc, char **argv) for (j=0;jid[n]=1; Np++; // Initialize gradient ColorGrad = (1,2,3) double value=double(3*k+2*j+i); @@ -160,7 +160,7 @@ int main(int argc, char **argv) } } } - Dm.CommInit(); + Dm->CommInit(); MPI_Barrier(comm); if (rank == 0) cout << "Domain set." << endl; if (rank==0) printf ("Create ScaLBL_Communicator \n"); @@ -177,7 +177,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Np]; - ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np); + ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); MPI_Barrier(comm); //......................device distributions................................. @@ -231,7 +231,7 @@ int main(int argc, char **argv) for (j=1;j 0){ + if (Dm->id[n] > 0){ int idx = Map(i,j,k); CX=COLORGRAD[idx]; CY=COLORGRAD[Np+idx]; From 0c594fa33d7e872ea507ee8714e7f12e3055b6f4 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 30 Sep 2020 14:45:51 -0400 Subject: [PATCH 24/39] add halo to memory optimized layout --- common/ScaLBL.cpp | 24 +++++++++++------------- common/ScaLBL.h | 2 +- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index 13508dd4..a2bd1d56 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -352,7 +352,7 @@ void ScaLBL_Communicator::D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, int *list, i delete [] ReturnDist; } -int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np){ +int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np, int width){ /* * Generate a memory optimized layout * id[n] == 0 implies that site n should be ignored (treat as a mask) @@ -391,28 +391,26 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis n = k*Nx*Ny+j*Nx+i; if (id[n] > 0){ // Counts for the six faces - if (i==1) Map(n)=idx++; - else if (j==1) Map(n)=idx++; - else if (k==1) Map(n)=idx++; - else if (i==Nx-2) Map(n)=idx++; - else if (j==Ny-2) Map(n)=idx++; - else if (k==Nz-2) Map(n)=idx++; + if (i>0 && i<=width) Map(n)=idx++; + else if (j>0 && j<=width)) Map(n)=idx++; + else if (k>0 && k<=width)) Map(n)=idx++; + else if (i>Nx-width-1 && iNy-width-1 && jNz-width-1 && k 0 ){ diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 2e8eef1a..dc199639 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -295,7 +295,7 @@ public: int FirstInterior(); int LastInterior(); - int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np); + int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np, int width); void SendD3Q19AA(double *dist); void RecvD3Q19AA(double *dist); // void BiSendD3Q7(double *A_even, double *A_odd, double *B_even, double *B_odd); From 624ec5a14591edb6b486fc1a62517226d1a794f9 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 30 Sep 2020 15:10:18 -0400 Subject: [PATCH 25/39] memory optimized layout with halo width --- common/ScaLBL.cpp | 6 +++--- models/ColorModel.cpp | 2 +- models/DFHModel.cpp | 2 +- models/FreeLeeModel.cpp | 2 +- models/GreyscaleColorModel.cpp | 2 +- models/GreyscaleModel.cpp | 2 +- models/MRTModel.cpp | 2 +- tests/TestBubbleDFH.cpp | 2 +- tests/TestColorGrad.cpp | 3 +-- tests/TestColorGradDFH.cpp | 2 +- tests/TestColorMassBounceback.cpp | 2 +- tests/TestCommD3Q19.cpp | 2 +- tests/TestFluxBC.cpp | 2 +- tests/TestForceMoments.cpp | 2 +- tests/TestMRT.cpp | 2 +- tests/TestMap.cpp | 2 +- tests/TestPressVel.cpp | 2 +- tests/TestWideHalo.cpp | 2 +- 18 files changed, 20 insertions(+), 21 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index a2bd1d56..28264d17 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -391,9 +391,9 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis n = k*Nx*Ny+j*Nx+i; if (id[n] > 0){ // Counts for the six faces - if (i>0 && i<=width) Map(n)=idx++; - else if (j>0 && j<=width)) Map(n)=idx++; - else if (k>0 && k<=width)) Map(n)=idx++; + if (i>0 && i<=width) Map(n)=idx++; + else if (j>0 && j<=width) Map(n)=idx++; + else if (k>0 && k<=width) Map(n)=idx++; else if (i>Nx-width-1 && iNy-width-1 && jNz-width-1 && kMemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); MPI_Barrier(comm); //........................................................................... diff --git a/models/DFHModel.cpp b/models/DFHModel.cpp index 4eb03bea..1f796f54 100644 --- a/models/DFHModel.cpp +++ b/models/DFHModel.cpp @@ -205,7 +205,7 @@ void ScaLBL_DFHModel::Create(){ if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); MPI_Barrier(comm); //........................................................................... diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp index e6525bf4..547885b8 100644 --- a/models/FreeLeeModel.cpp +++ b/models/FreeLeeModel.cpp @@ -206,7 +206,7 @@ void ScaLBL_FreeLeeModel::Create(){ if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,2); MPI_Barrier(comm); //........................................................................... diff --git a/models/GreyscaleColorModel.cpp b/models/GreyscaleColorModel.cpp index d97c844d..736012ce 100644 --- a/models/GreyscaleColorModel.cpp +++ b/models/GreyscaleColorModel.cpp @@ -641,7 +641,7 @@ void ScaLBL_GreyscaleColorModel::Create(){ if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); MPI_Barrier(comm); //........................................................................... diff --git a/models/GreyscaleModel.cpp b/models/GreyscaleModel.cpp index 85832a4b..e61bdec8 100644 --- a/models/GreyscaleModel.cpp +++ b/models/GreyscaleModel.cpp @@ -313,7 +313,7 @@ void ScaLBL_GreyscaleModel::Create(){ if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); MPI_Barrier(comm); //........................................................................... diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index acfb8821..8b18b238 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -170,7 +170,7 @@ void ScaLBL_MRTModel::Create(){ if (rank==0) printf ("Set up memory efficient layout \n"); Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); MPI_Barrier(comm); //........................................................................... // MAIN VARIABLES ALLOCATED HERE diff --git a/tests/TestBubbleDFH.cpp b/tests/TestBubbleDFH.cpp index a8ba0cde..0841311e 100644 --- a/tests/TestBubbleDFH.cpp +++ b/tests/TestBubbleDFH.cpp @@ -249,7 +249,7 @@ int main(int argc, char **argv) if (rank==0) printf ("Set up memory efficient layout Npad=%i \n",Npad); IntArray Map(Nx,Ny,Nz); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); MPI_Barrier(comm); //........................................................................... diff --git a/tests/TestColorGrad.cpp b/tests/TestColorGrad.cpp index d2f1d6de..f7808f93 100644 --- a/tests/TestColorGrad.cpp +++ b/tests/TestColorGrad.cpp @@ -1,4 +1,3 @@ - //************************************************************************* // Lattice Boltzmann Simulator for Single Phase Flow in Porous Media // James E. McCLure @@ -177,7 +176,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Np]; - ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); + ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1); MPI_Barrier(comm); //......................device distributions................................. diff --git a/tests/TestColorGradDFH.cpp b/tests/TestColorGradDFH.cpp index d6376d82..cd18e401 100644 --- a/tests/TestColorGradDFH.cpp +++ b/tests/TestColorGradDFH.cpp @@ -104,7 +104,7 @@ int main(int argc, char **argv) int *neighborList; IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1); MPI_Barrier(comm); //......................device distributions................................. diff --git a/tests/TestColorMassBounceback.cpp b/tests/TestColorMassBounceback.cpp index c05c245e..678b4dfb 100644 --- a/tests/TestColorMassBounceback.cpp +++ b/tests/TestColorMassBounceback.cpp @@ -169,7 +169,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); Npad=Np+32; neighborList= new int[18*Npad]; - Np=ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); + Np=ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1); MPI_Barrier(comm); //......................device distributions................................. diff --git a/tests/TestCommD3Q19.cpp b/tests/TestCommD3Q19.cpp index e1fa821f..f74eac31 100644 --- a/tests/TestCommD3Q19.cpp +++ b/tests/TestCommD3Q19.cpp @@ -284,7 +284,7 @@ int main(int argc, char **argv) auto neighborList= new int[18*Npad]; IntArray Map(Nx,Ny,Nz); Map.fill(-2); - Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); + Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1); MPI_Barrier(comm); int neighborSize=18*Np*sizeof(int); //......................device distributions................................. diff --git a/tests/TestFluxBC.cpp b/tests/TestFluxBC.cpp index 020bbd89..6eac35d2 100644 --- a/tests/TestFluxBC.cpp +++ b/tests/TestFluxBC.cpp @@ -88,7 +88,7 @@ int main (int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1); MPI_Barrier(comm); //......................device distributions................................. diff --git a/tests/TestForceMoments.cpp b/tests/TestForceMoments.cpp index 1fb1e0a4..54324103 100644 --- a/tests/TestForceMoments.cpp +++ b/tests/TestForceMoments.cpp @@ -183,7 +183,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Np]; - ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); + ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1); if (rank == 0) PrintNeighborList(neighborList,Np, rank); diff --git a/tests/TestMRT.cpp b/tests/TestMRT.cpp index 30f46689..2b10411e 100644 --- a/tests/TestMRT.cpp +++ b/tests/TestMRT.cpp @@ -705,7 +705,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Np]; - ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np); + ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np,1); MPI_Barrier(comm); //......................device distributions................................. diff --git a/tests/TestMap.cpp b/tests/TestMap.cpp index a47c0d9e..ff0ba811 100644 --- a/tests/TestMap.cpp +++ b/tests/TestMap.cpp @@ -93,7 +93,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1); MPI_Barrier(comm); // Check the neighborlist diff --git a/tests/TestPressVel.cpp b/tests/TestPressVel.cpp index e655ced9..f66c1d2c 100644 --- a/tests/TestPressVel.cpp +++ b/tests/TestPressVel.cpp @@ -132,7 +132,7 @@ int main(int argc, char **argv) int *neighborList; IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1); MPI_Barrier(comm); //......................device distributions................................. diff --git a/tests/TestWideHalo.cpp b/tests/TestWideHalo.cpp index d2f1d6de..d5bb218c 100644 --- a/tests/TestWideHalo.cpp +++ b/tests/TestWideHalo.cpp @@ -177,7 +177,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Np]; - ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); + ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,2); MPI_Barrier(comm); //......................device distributions................................. From 5ed82a54f46afecb1db05a83ee0ab7425a878b38 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 30 Sep 2020 15:14:28 -0400 Subject: [PATCH 26/39] fix failing commm test --- common/ScaLBL.cpp | 6 +++--- tests/TestCommD3Q19.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index a2bd1d56..f0beb25b 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -394,9 +394,9 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis if (i>0 && i<=width) Map(n)=idx++; else if (j>0 && j<=width)) Map(n)=idx++; else if (k>0 && k<=width)) Map(n)=idx++; - else if (i>Nx-width-1 && iNy-width-1 && jNz-width-1 && kNx-width-2 && iNy-width-2 && jNz-width-2 && k Date: Wed, 30 Sep 2020 15:22:37 -0400 Subject: [PATCH 27/39] update Memory Optimized Layout -tests passing --- common/ScaLBL.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index 09fef5d4..ecc9647f 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -392,8 +392,8 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis if (id[n] > 0){ // Counts for the six faces if (i>0 && i<=width) Map(n)=idx++; - else if (j>0 && j<=width)) Map(n)=idx++; - else if (k>0 && k<=width)) Map(n)=idx++; + else if (j>0 && j<=width) Map(n)=idx++; + else if (k>0 && k<=width) Map(n)=idx++; else if (i>Nx-width-2 && iNy-width-2 && jNz-width-2 && k Date: Thu, 1 Oct 2020 16:36:16 -0400 Subject: [PATCH 28/39] adding wide halo gradient test --- common/ScaLBL.h | 100 +---------------------------------------- tests/TestWideHalo.cpp | 29 ++++++------ 2 files changed, 18 insertions(+), 111 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index dc199639..65211440 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -77,104 +77,6 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_MRT(double *dist, int start, int f extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz, double *Poros,double *Perm, double *Velocity,double Den,double *Pressure); -// GREYSCALE FREE-ENERGY MODEL (Two-component) - -//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleFE(double *dist, double *Aq, double *Bq, double *Den, -// double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np, -// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz, -// double *Poros,double *Perm, double *Velocity,double *Pressure); -// -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleFE(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, -// double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np, -// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz, -// double *Poros,double *Perm, double *Velocity,double *Pressure); -// -//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleFEChem(double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np, -// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, -// double Gx, double Gy, double Gz, -// double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap); -// -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleFEChem(int *neighborList, double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np, -// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, -// double Gx, double Gy, double Gz, -// double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap); -// -//extern "C" void ScaLBL_D3Q7_GreyscaleFE_Init(double *Den, double *Cq, double *PhiLap, double gamma, double kappaA, double kappaB, double lambdaA, double lambdaB, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_GreyscaleFE_IMRT_Init(double *dist, double *Den, double rhoA, double rhoB, int Np); -// -//extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleFEDensity(int *NeighborList, double *Aq, double *Bq, double *Den, double *Phi, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleFEDensity(double *Aq, double *Bq, double *Den, double *Phi, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleFEPhi(int *NeighborList, double *Cq, double *Phi, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleFEPhi(double *Cq, double *Phi, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_GreyscaleFE_Gradient(int *neighborList, double *Den, double *DenGrad, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_GreyscaleFE_Laplacian(int *neighborList, double *Den, double *DenLap, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_GreyscaleFE_Pressure(double *dist, double *Den, double *Porosity,double *Velocity, -// double *Pressure, double rhoA,double rhoB, int Np); -// -//extern "C" void ScaLBL_D3Q19_GreyscaleFE_PressureTensor(int *neighborList, double *Phi,double *Pressure, double *PressTensor, double *PhiLap, -// double kappaA,double kappaB,double lambdaA,double lambdaB, int start, int finish, int Np); - -// GREYSCALE SHAN-CHEN MODEL (Two-component) - -//extern "C" void ScaLBL_D3Q19_GreyscaleSC_Init(int *Map, double *distA, double *distB, double *DenA, double *DenB, int Np); -// -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(int *NeighborList, int *Map, double *distA, double *distB, double *DenA, double *DenB, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(int *Map, double *distA, double *distB, double *DenA, double *DenB, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, int *Mpa, double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB, -// double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, -// double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, -// int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(int *Map,double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB, -// double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, -// double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, -// int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK(int *neighborList, int *Map, double *distA, double *distB, double *DenA, double *DenB, double *DenGradA, double *DenGradB, -// double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, -// double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, -// int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK(int *Map, double *distA, double *distB, double *DenA, double *DenB, double *DenGradA, double *DenGradB, -// double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, -// double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, -// int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_GreyscaleSC_Gradient(int *neighborList, int *Map, double *Den, double *DenGrad, int strideY, int strideZ,int start, int finish, int Np); -// -//extern "C" void ScaLBL_GreyscaleSC_BC_z(int *list, int *Map, double *DenA, double *DenB, double vA, double vB, int count); -// -//extern "C" void ScaLBL_GreyscaleSC_BC_Z(int *list, int *Map, double *DenA, double *DenB, double vA, double vB, int count); -// -//extern "C" void ScaLBL_GreyscaleSC_AAeven_Pressure_BC_z(int *list, double *distA, double *distB, double dinA, double dinB, int count, int N); -// -//extern "C" void ScaLBL_GreyscaleSC_AAeven_Pressure_BC_Z(int *list, double *distA, double *distB, double doutA, double doutB, int count, int N); -// -//extern "C" void ScaLBL_GreyscaleSC_AAodd_Pressure_BC_z(int *neighborList, int *list, double *distA, double *distB, double dinA, double dinB, int count, int N); -// -//extern "C" void ScaLBL_GreyscaleSC_AAodd_Pressure_BC_Z(int *neighborList, int *list, double *distA, double *distB, double doutA, double doutB, int count, int N); - -// GREYSCALE COLOR MODEL (Two-component) -//extern "C" void ScaLBL_D3Q19_GreyscaleColor_Init(double *dist, double *Porosity, int Np); - -//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den, -// double *ColorGrad,double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel, -// double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, -// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); -// -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, -// double *ColorGrad,double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel, -// double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta, -// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel, @@ -211,6 +113,8 @@ extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq, extern "C" void ScaLBL_D3Q19_Gradient(int *Map, double *Phi, double *ColorGrad, int start, int finish, int Np, int Nx, int Ny, int Nz); +extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradient, int start, int finish, int Np, int Nx, int Ny, int Nz); + extern "C" void ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np); // Density functional hydrodynamics LBM diff --git a/tests/TestWideHalo.cpp b/tests/TestWideHalo.cpp index d5bb218c..b6dd970b 100644 --- a/tests/TestWideHalo.cpp +++ b/tests/TestWideHalo.cpp @@ -168,6 +168,7 @@ int main(int argc, char **argv) //Create a second communicator based on the regular data layout ScaLBL_Communicator ScaLBL_Comm_Regular(Dm); ScaLBL_Communicator ScaLBL_Comm(Dm); + ScaLBLWideHalo_Communicator WideHalo(Dm); // LBM variables if (rank==0) printf ("Set up the neighborlist \n"); @@ -197,20 +198,19 @@ int main(int argc, char **argv) //........................................................................... // Update GPU data structures if (rank==0) printf ("Setting up device map and neighbor list \n"); - int *TmpMap; - TmpMap=new int[Np*sizeof(int)]; - for (k=1; k Date: Thu, 1 Oct 2020 16:41:17 -0400 Subject: [PATCH 29/39] adding wide halo gradient test --- common/ScaLBL.h | 1 + tests/TestWideHalo.cpp | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 65211440..998ad74b 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -10,6 +10,7 @@ #ifndef ScalLBL_H #define ScalLBL_H #include "common/Domain.h" +#include "common/WideHalo.h" extern "C" int ScaLBL_SetDevice(int rank); diff --git a/tests/TestWideHalo.cpp b/tests/TestWideHalo.cpp index b6dd970b..b78426e0 100644 --- a/tests/TestWideHalo.cpp +++ b/tests/TestWideHalo.cpp @@ -218,9 +218,9 @@ int main(int argc, char **argv) ScaLBL_CopyToDevice(Phi, PhaseLabel, N*sizeof(double)); //........................................................................... - Nxh = Nx+2; - Nyh = Ny+2; - Nzh = Nz+2; + int Nxh = Nx+2; + int Nyh = Ny+2; + int Nzh = Nz+2; ScaLBL_D3Q19_MixedGradient(dvcMap, Phi, ColorGrad, 0, Np, Np, Nxh, Nyh, Nzh); double *COLORGRAD; From 8b095b24ada5970a3b051f7b7ed79c3f8213ec6b Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 5 Oct 2020 14:36:26 -0400 Subject: [PATCH 30/39] added mix gradient --- cpu/MixedGradient.cpp | 48 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 cpu/MixedGradient.cpp diff --git a/cpu/MixedGradient.cpp b/cpu/MixedGradient.cpp new file mode 100644 index 00000000..8aa7c53f --- /dev/null +++ b/cpu/MixedGradient.cpp @@ -0,0 +1,48 @@ +/* Implement Mixed Gradient (Lee et al. JCP 2016)*/ + +extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradient, int start, int finish, int Np, int Nx, int Ny, int Nz) +{ + static int D3Q19[18][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}, + {1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0}, + {1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1}, + {0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}}; + + int i,j,k,n,N; + int np,np2,nm; // neighbors + double v,vp,vp2,vm; // values at neighbors + double grad; + for (int idx=start; idx Date: Mon, 5 Oct 2020 14:36:50 -0400 Subject: [PATCH 31/39] fixing build --- common/ScaLBL.h | 1 - common/WideHalo.h | 5 +++-- tests/TestWideHalo.cpp | 3 ++- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 998ad74b..65211440 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -10,7 +10,6 @@ #ifndef ScalLBL_H #define ScalLBL_H #include "common/Domain.h" -#include "common/WideHalo.h" extern "C" int ScaLBL_SetDevice(int rank); diff --git a/common/WideHalo.h b/common/WideHalo.h index 24429e3f..601eda13 100644 --- a/common/WideHalo.h +++ b/common/WideHalo.h @@ -106,8 +106,9 @@ private: } } } - ScaLBL_AllocateZeroCopy((void **) &dvcList, count*sizeof(int)); // Allocate device memory - ScaLBL_CopyToZeroCopy(dvcList,List,count*sizeof(int)); + size_t numbytes=count*sizeof(int); + ScaLBL_AllocateZeroCopy((void **) &dvcList, numbytes); // Allocate device memory + ScaLBL_CopyToZeroCopy(dvcList,List,numbytes); return count; } diff --git a/tests/TestWideHalo.cpp b/tests/TestWideHalo.cpp index b78426e0..cc29a15d 100644 --- a/tests/TestWideHalo.cpp +++ b/tests/TestWideHalo.cpp @@ -7,6 +7,7 @@ #include #include #include "common/ScaLBL.h" +#include "common/WideHalo.h" #include "common/MPI_Helpers.h" using namespace std; @@ -168,7 +169,7 @@ int main(int argc, char **argv) //Create a second communicator based on the regular data layout ScaLBL_Communicator ScaLBL_Comm_Regular(Dm); ScaLBL_Communicator ScaLBL_Comm(Dm); - ScaLBLWideHalo_Communicator WideHalo(Dm); + ScaLBLWideHalo_Communicator WideHalo(Dm,2); // LBM variables if (rank==0) printf ("Set up the neighborlist \n"); From db5579530de88836f84a026689682ff29bb1051f Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Tue, 6 Oct 2020 14:05:26 -0400 Subject: [PATCH 32/39] mix gradient compiles --- cpu/MixedGradient.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpu/MixedGradient.cpp b/cpu/MixedGradient.cpp index 8aa7c53f..841dbdf1 100644 --- a/cpu/MixedGradient.cpp +++ b/cpu/MixedGradient.cpp @@ -29,7 +29,7 @@ extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradie vp = Phi[np]; vp2 = Phi[np2]; vm = Phi[nm]; - grad += 0.25*(5.0*vp1-vp2-3.0*v-vm); + grad += 0.25*(5.0*vp-vp2-3.0*v-vm); } for (int q=6; q<18; q++){ int iqx = D3Q19[q][0]; @@ -41,7 +41,7 @@ extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradie vp = Phi[np]; vp2 = Phi[np2]; vm = Phi[nm]; - grad += 0.125*(5.0*vp1-vp2-3.0*v-vm); + grad += 0.125*(5.0*vp-vp2-3.0*v-vm); } Gradient[n] = grad; } From 9b8dd50f78db9c1109f919471aaa59a6bf0a1c84 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Thu, 8 Oct 2020 13:17:15 -0400 Subject: [PATCH 33/39] wide halo compiles --- common/WideHalo.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/common/WideHalo.cpp b/common/WideHalo.cpp index 48f83cd9..ecc563f3 100644 --- a/common/WideHalo.cpp +++ b/common/WideHalo.cpp @@ -291,6 +291,11 @@ void ScaLBLWideHalo_Communicator::Send(double *data){ MPI_Irecv(recvbuf_XYz, recvCount_XYz,MPI_DOUBLE,rank_XYz,recvtag,MPI_COMM_SCALBL,&req2[25]); //................................................................................... } + + +ScaLBLWideHalo_Communicator::~ScaLBLWideHalo_Communicator() +{ +} void ScaLBLWideHalo_Communicator::Recv(double *data){ //................................................................................... From 35447181dfeb81e0b3e8de854f7e1fcbcc313190 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Fri, 9 Oct 2020 15:19:59 -0400 Subject: [PATCH 34/39] merging multi-halo with electrokinetic --- common/ScaLBL.h | 10 ---------- models/IonModel.cpp | 2 +- models/PoissonSolver.cpp | 2 +- models/StokesModel.cpp | 2 +- 4 files changed, 3 insertions(+), 13 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 9992f226..39905ee6 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -268,13 +268,8 @@ public: int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np, int width); void SendD3Q19AA(double *dist); void RecvD3Q19AA(double *dist); -<<<<<<< HEAD -// void BiSendD3Q7(double *A_even, double *A_odd, double *B_even, double *B_odd); -// void BiRecvD3Q7(double *A_even, double *A_odd, double *B_even, double *B_odd); -======= void SendD3Q7AA(double *fq, int Component); void RecvD3Q7AA(double *fq, int Component); ->>>>>>> electrokinetic void BiSendD3Q7AA(double *Aq, double *Bq); void BiRecvD3Q7AA(double *Aq, double *Bq); void TriSendD3Q7AA(double *Aq, double *Bq, double *Cq); @@ -295,21 +290,16 @@ public: void D3Q19_Reflection_BC_z(double *fq); void D3Q19_Reflection_BC_Z(double *fq); double D3Q19_Flux_BC_z(int *neighborList, double *fq, double flux, int time); -<<<<<<< HEAD void GreyscaleSC_BC_z(int *Map, double *DenA, double *DenB, double vA, double vB); void GreyscaleSC_BC_Z(int *Map, double *DenA, double *DenB, double vA, double vB); void GreyscaleSC_Pressure_BC_z(int *neighborList, double *fqA, double *fqB, double dinA, double dinB, int time); void GreyscaleSC_Pressure_BC_Z(int *neighborList, double *fqA, double *fqB, double doutA, double doutB, int time); -// void TestSendD3Q19(double *f_even, double *f_odd); -// void TestRecvD3Q19(double *f_even, double *f_odd); -======= void D3Q7_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time); void D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time); void Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin); void Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout); void D3Q7_Ion_Concentration_BC_z(int *neighborList, double *fq, double Cin, int time); void D3Q7_Ion_Concentration_BC_Z(int *neighborList, double *fq, double Cout, int time); ->>>>>>> electrokinetic // Debugging and unit testing functions void PrintD3Q19(); diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 24c4d8bd..eb9101fa 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -593,7 +593,7 @@ void ScaLBL_IonModel::Create(){ if (rank==0) printf ("LB Ion Solver: Set up memory efficient layout \n"); Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); MPI_Barrier(comm); //........................................................................... // MAIN VARIABLES ALLOCATED HERE diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 3e11de0a..636b88d1 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -271,7 +271,7 @@ void ScaLBL_Poisson::Create(){ if (rank==0) printf ("LB-Poisson Solver: Set up memory efficient layout \n"); Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); MPI_Barrier(comm); //........................................................................... // MAIN VARIABLES ALLOCATED HERE diff --git a/models/StokesModel.cpp b/models/StokesModel.cpp index 964baaae..891ea480 100644 --- a/models/StokesModel.cpp +++ b/models/StokesModel.cpp @@ -278,7 +278,7 @@ void ScaLBL_StokesModel::Create(){ if (rank==0) printf ("LB Single-Fluid Solver: Set up memory efficient layout \n"); Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); MPI_Barrier(comm); //........................................................................... // MAIN VARIABLES ALLOCATED HERE From 29e486a1a078bb54dcbca7e29b924712cc07390d Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Mon, 12 Oct 2020 09:08:53 -0400 Subject: [PATCH 35/39] moving electrochem model tests to executables --- tests/CMakeLists.txt | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index fd66d8e4..28a02749 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -35,8 +35,13 @@ ADD_LBPM_EXECUTABLE( GenerateSphereTest ) #ADD_LBPM_EXECUTABLE( BlobIdentifyParallel ) #ADD_LBPM_EXECUTABLE( convertIO ) #ADD_LBPM_EXECUTABLE( DataAggregator ) -#ADD_LBPM_EXECUTABLE( BlobAnalyzeParallel ) +#ADD_LBPM_EXECUTABLE( BlobAnalyzeParallel )( ADD_LBPM_EXECUTABLE( lbpm_minkowski_scalar ) +ADD_LBPM_EXECUTABLE( TestPoissonSolver ) +ADD_LBPM_EXECUTABLE( TestIonModel ) +ADD_LBPM_EXECUTABLE( TestNernstPlanck ) +ADD_LBPM_EXECUTABLE( TestPNP_Stokes ) + CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_DIR}/cylindertest COPYONLY ) @@ -48,10 +53,6 @@ ADD_LBPM_TEST( TestTorusEvolve ) ADD_LBPM_TEST( TestTopo3D ) ADD_LBPM_TEST( TestFluxBC ) ADD_LBPM_TEST( TestMap ) -ADD_LBPM_TEST( TestPoissonSolver ) -ADD_LBPM_TEST( TestIonModel ) -ADD_LBPM_TEST( TestNernstPlanck ) -ADD_LBPM_TEST( TestPNP_Stokes ) #ADD_LBPM_TEST( TestMRT ) ADD_LBPM_TEST( TestColorGrad ) ADD_LBPM_TEST( TestWideHalo ) From a066fa66062a30d2c119b76a3ac0b2a960d1adb5 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 13 Jan 2021 22:06:08 -0500 Subject: [PATCH 36/39] merge Lee model --- models/FreeLeeModel.cpp | 14 +++++++------- models/FreeLeeModel.h | 4 ++-- tests/CMakeLists.txt | 5 ----- tests/TestMRT.cpp | 5 ----- 4 files changed, 9 insertions(+), 19 deletions(-) diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp index 547885b8..5d0f64f6 100644 --- a/models/FreeLeeModel.cpp +++ b/models/FreeLeeModel.cpp @@ -9,7 +9,7 @@ color lattice boltzmann model #include #include -ScaLBL_FreeLeeModel::ScaLBL_FreeLeeModel(int RANK, int NP, MPI_Comm COMM): +ScaLBL_FreeLeeModel::ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM): rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rhoA(0),rhoB(0),W(0),gamma(0), Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0), Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM) @@ -107,9 +107,9 @@ void ScaLBL_FreeLeeModel::SetDomain(){ id = new signed char [N]; for (int i=0; iid[i] = 1; // initialize this way - MPI_Barrier(comm); + comm.barrier() Dm->CommInit(); - MPI_Barrier(comm); + comm.barrier() // Read domain parameters rank = Dm->rank(); nprocx = Dm->nprocx(); @@ -206,8 +206,8 @@ void ScaLBL_FreeLeeModel::Create(){ if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,2); - MPI_Barrier(comm); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,2); + comm.barrier() //........................................................................... // MAIN VARIABLES ALLOCATED HERE @@ -335,7 +335,7 @@ void ScaLBL_FreeLeeModel::Initialize(){ ScaLBL_CopyToDevice(Phi,cPhi,N*sizeof(double)); ScaLBL_DeviceBarrier(); - MPI_Barrier(comm); + comm.barrier() } if (rank==0) printf ("Initializing phase field \n"); @@ -371,7 +371,7 @@ void ScaLBL_FreeLeeModel::Run(){ //.......create and start timer............ double starttime,stoptime,cputime; ScaLBL_DeviceBarrier(); - MPI_Barrier(comm); + comm.barrier() starttime = MPI_Wtime(); //......................................... diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h index 01fb54c3..437f48b3 100644 --- a/models/FreeLeeModel.h +++ b/models/FreeLeeModel.h @@ -18,7 +18,7 @@ Implementation of Lee et al JCP 2016 lattice boltzmann model class ScaLBL_FreeLeeModel{ public: - ScaLBL_FreeLeeModel(int RANK, int NP, MPI_Comm COMM); + ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM); ~ScaLBL_FreeLeeModel(); // functions in they should be run @@ -70,7 +70,7 @@ public: DoubleArray SignDist; private: - MPI_Comm comm; + Utilities::MPI comm; int dist_mem_size; int neighborSize; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2bb810bd..414a13cf 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -4,13 +4,8 @@ ADD_LBPM_EXECUTABLE( lbpm_color_simulator ) ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator ) ADD_LBPM_EXECUTABLE( lbpm_greyscale_simulator ) -<<<<<<< HEAD -ADD_LBPM_EXECUTABLE( lbpm_greyscaleColor_simulator ) -ADD_LBPM_EXECUTABLE( lbpm_electrokinetic_SingleFluid_simulator ) -======= ADD_LBPM_EXECUTABLE( lbpm_electrokinetic_SingleFluid_simulator ) ADD_LBPM_EXECUTABLE( lbpm_greyscaleColor_simulator ) ->>>>>>> FOM #ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator ) #ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator ) ADD_LBPM_EXECUTABLE( lbpm_dfh_simulator ) diff --git a/tests/TestMRT.cpp b/tests/TestMRT.cpp index 5b58130e..172f6c55 100644 --- a/tests/TestMRT.cpp +++ b/tests/TestMRT.cpp @@ -804,13 +804,8 @@ int main(int argc, char **argv) } // **************************************************** -<<<<<<< HEAD comm.barrier(); Utilities::shutdown(); -======= - MPI_Barrier(comm); - MPI_Finalize(); ->>>>>>> electrokinetic // **************************************************** return check; From 1d85db88ed2c1c764924f3d6f69996e825d215ea Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Fri, 15 Jan 2021 15:38:15 -0500 Subject: [PATCH 37/39] adding Lee model to FOM --- common/WideHalo.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/common/WideHalo.h b/common/WideHalo.h index 601eda13..d5c97c2f 100644 --- a/common/WideHalo.h +++ b/common/WideHalo.h @@ -11,7 +11,8 @@ public: ScaLBLWideHalo_Communicator(std::shared_ptr Dm, int width); ~ScaLBLWideHalo_Communicator(); //...................................................................................... - MPI_Comm MPI_COMM_SCALBL; // MPI Communicator + //MPI_Comm MPI_COMM_SCALBL; // MPI Communicator + Utilities::MPI MPI_COMM_SCALBL; unsigned long int CommunicationCount,SendCount,RecvCount; int Nx,Ny,Nz,N; // original domain structure int Nxh,Nyh,Nzh,Nh; // with wide halo @@ -32,7 +33,6 @@ public: double *recvbuf_xyz, *recvbuf_Xyz, *recvbuf_xYz, *recvbuf_XYz; double *recvbuf_xyZ, *recvbuf_XyZ, *recvbuf_xYZ, *recvbuf_XYZ; //...................................................................................... - int LastExterior(); int FirstInterior(); int LastInterior(); From 9bc8100a1da945e43bee14bb961670df898f3a61 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 15 Jan 2021 16:33:57 -0500 Subject: [PATCH 38/39] debugging wide halo --- common/WideHalo.cpp | 109 ++++++++++++++++++++-------------------- models/FreeLeeModel.cpp | 15 +++--- models/FreeLeeModel.h | 2 +- tests/TestBubbleDFH.cpp | 7 +-- tests/TestColorGrad.cpp | 29 ++++++++--- 5 files changed, 85 insertions(+), 77 deletions(-) diff --git a/common/WideHalo.cpp b/common/WideHalo.cpp index ecc563f3..ed4a451d 100644 --- a/common/WideHalo.cpp +++ b/common/WideHalo.cpp @@ -9,7 +9,7 @@ ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr Comm,&MPI_COMM_SCALBL); + MPI_COMM_SCALBL = Dm->Comm.dup(); //...................................................................................... // Copy the domain size and communication information directly from Dm Nx = Dm->Nx; @@ -59,7 +59,7 @@ ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr id[i] = 1; // initialize this way - comm.barrier() + comm.barrier(); Dm->CommInit(); - comm.barrier() + comm.barrier(); // Read domain parameters rank = Dm->rank(); nprocx = Dm->nprocx(); @@ -139,7 +139,7 @@ void ScaLBL_FreeLeeModel::ReadInput(){ ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 ); fillHalo fill( MPI_COMM_WORLD, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 ); Array id_view; - id_view.viewRaw( size1, Mask->id ); + id_view.viewRaw( size1, Mask->id.data() ); fill.copy( input_id, id_view ); fill.fill( id_view ); } @@ -207,7 +207,7 @@ void ScaLBL_FreeLeeModel::Create(){ Map.resize(Nx,Ny,Nz); Map.fill(-2); auto neighborList= new int[18*Npad]; Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,2); - comm.barrier() + comm.barrier(); //........................................................................... // MAIN VARIABLES ALLOCATED HERE @@ -335,7 +335,7 @@ void ScaLBL_FreeLeeModel::Initialize(){ ScaLBL_CopyToDevice(Phi,cPhi,N*sizeof(double)); ScaLBL_DeviceBarrier(); - comm.barrier() + comm.barrier(); } if (rank==0) printf ("Initializing phase field \n"); @@ -371,7 +371,7 @@ void ScaLBL_FreeLeeModel::Run(){ //.......create and start timer............ double starttime,stoptime,cputime; ScaLBL_DeviceBarrier(); - comm.barrier() + comm.barrier(); starttime = MPI_Wtime(); //......................................... @@ -460,8 +460,7 @@ void ScaLBL_FreeLeeModel::Run(){ ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); */ - ScaLBL_DeviceBarrier(); - MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL); + ScaLBL_Comm->Barrier(); //************************************************************************ PROFILE_STOP("Update"); } diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h index 437f48b3..5aa2d30a 100644 --- a/models/FreeLeeModel.h +++ b/models/FreeLeeModel.h @@ -10,7 +10,7 @@ Implementation of Lee et al JCP 2016 lattice boltzmann model #include #include "common/Communication.h" -#include "common/MPI_Helpers.h" +#include "common/MPI.h" #include "ProfilerApp.h" #include "threadpool/thread_pool.h" #include "common/ScaLBL.h" diff --git a/tests/TestBubbleDFH.cpp b/tests/TestBubbleDFH.cpp index 77a8aadb..fe5c5a3d 100644 --- a/tests/TestBubbleDFH.cpp +++ b/tests/TestBubbleDFH.cpp @@ -247,13 +247,8 @@ int main(int argc, char **argv) if (rank==0) printf ("Set up memory efficient layout Npad=%i \n",Npad); IntArray Map(Nx,Ny,Nz); auto neighborList= new int[18*Npad]; -<<<<<<< HEAD - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np,1); - MPI_Barrier(comm); -======= - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np); + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1); comm.barrier(); ->>>>>>> FOM //........................................................................... // MAIN VARIABLES ALLOCATED HERE diff --git a/tests/TestColorGrad.cpp b/tests/TestColorGrad.cpp index 4fdedb7b..ac59fc38 100644 --- a/tests/TestColorGrad.cpp +++ b/tests/TestColorGrad.cpp @@ -1,3 +1,4 @@ + //************************************************************************* // Lattice Boltzmann Simulator for Single Phase Flow in Porous Media // James E. McCLure @@ -6,7 +7,7 @@ #include #include #include "common/ScaLBL.h" -#include "common/MPI_Helpers.h" +#include "common/MPI.h" using namespace std; @@ -70,7 +71,20 @@ int main(int argc, char **argv) //....................................................................... // Reading the domain information file //....................................................................... - if (nprocs==1){ + ifstream domain("Domain.in"); + if (domain.good()){ + domain >> nprocx; + domain >> nprocy; + domain >> nprocz; + domain >> Nx; + domain >> Ny; + domain >> Nz; + domain >> nspheres; + domain >> Lx; + domain >> Ly; + domain >> Lz; + } + else if (nprocs==1){ nprocx=nprocy=nprocz=1; Nx=Ny=Nz=3; nspheres=0; @@ -136,7 +150,8 @@ int main(int argc, char **argv) double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; int BoundaryCondition=0; - std::shared_ptr Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition)); + Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + Nx += 2; Ny += 2; Nz += 2; @@ -150,7 +165,7 @@ int main(int argc, char **argv) for (j=0;jid[n]=1; + Dm.id[n]=1; Np++; // Initialize gradient ColorGrad = (1,2,3) double value=double(3*k+2*j+i); @@ -158,7 +173,7 @@ int main(int argc, char **argv) } } } - Dm->CommInit(); + Dm.CommInit(); MPI_Barrier(comm); if (rank == 0) cout << "Domain set." << endl; if (rank==0) printf ("Create ScaLBL_Communicator \n"); @@ -175,7 +190,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Np]; - ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1); + ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np); MPI_Barrier(comm); //......................device distributions................................. @@ -229,7 +244,7 @@ int main(int argc, char **argv) for (j=1;jid[n] > 0){ + if (Dm.id[n] > 0){ int idx = Map(i,j,k); CX=COLORGRAD[idx]; CY=COLORGRAD[Np+idx]; From c911a4f3d2240f18e34aa680195fa57bc82a7c43 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 15 Jan 2021 17:12:38 -0500 Subject: [PATCH 39/39] add wide halo to fom --- tests/CMakeLists.txt | 3 +- tests/TestColorGrad.cpp | 109 +++-------------------------- tests/TestWideHalo.cpp | 147 +++++++++++++--------------------------- 3 files changed, 58 insertions(+), 201 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 414a13cf..0a8074a3 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -55,9 +55,8 @@ ADD_LBPM_TEST( TestTopo3D ) ADD_LBPM_TEST( TestFluxBC ) ADD_LBPM_TEST( TestMap ) #ADD_LBPM_TEST( TestMRT ) -ADD_LBPM_TEST( TestColorGrad ) +#ADD_LBPM_TEST( TestColorGrad ) ADD_LBPM_TEST( TestWideHalo ) -#ADD_LBPM_TEST( TestColorGradDFH ) ADD_LBPM_TEST( TestColorGradDFH ) ADD_LBPM_TEST( TestBubbleDFH ../example/Bubble/input.db) #ADD_LBPM_TEST( TestColorMassBounceback ../example/Bubble/input.db) diff --git a/tests/TestColorGrad.cpp b/tests/TestColorGrad.cpp index ac59fc38..9ea243ea 100644 --- a/tests/TestColorGrad.cpp +++ b/tests/TestColorGrad.cpp @@ -3,7 +3,12 @@ // Lattice Boltzmann Simulator for Single Phase Flow in Porous Media // James E. McCLure //************************************************************************* -#include +#include // Initialize MPI + Utilities::startup( argc, argv ); + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + int check; #include #include #include "common/ScaLBL.h" @@ -21,8 +26,8 @@ int main(int argc, char **argv) // Initialize MPI Utilities::startup( argc, argv ); Utilities::MPI comm( MPI_COMM_WORLD ); - int rank = comm.getRank(); - int nprocs = comm.getSize(); + int rank = comm.getRank(); + int nprocs = comm.getSize(); int check; { // parallel domain size (# of sub-domains) @@ -47,105 +52,13 @@ int main(int argc, char **argv) int Nx,Ny,Nz; int i,j,k,n; int dim = 3; + Nx = Ny = Nz = 32; + Lx = Ly = Lz = 1.0; //if (rank == 0) printf("dim=%d\n",dim); int timestep = 0; int timesteps = 100; int centralNode = 2; - double tauA = 1.0; - double tauB = 1.0; - double rhoA = 1.0; - double rhoB = 1.0; - double alpha = 0.005; - double beta = 0.95; - - double tau = 1.0; - double mu=(tau-0.5)/3.0; - double rlx_setA=1.0/tau; - double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); - - Fx = Fy = 0.f; - Fz = 0.f; - - if (rank==0){ - //....................................................................... - // Reading the domain information file - //....................................................................... - ifstream domain("Domain.in"); - if (domain.good()){ - domain >> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> Nx; - domain >> Ny; - domain >> Nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - } - else if (nprocs==1){ - nprocx=nprocy=nprocz=1; - Nx=Ny=Nz=3; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==2){ - nprocx=2; nprocy=1; - nprocz=1; - Nx=Ny=Nz=dim; - Nx = dim; Ny = dim; Nz = dim; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==4){ - nprocx=nprocy=2; - nprocz=1; - Nx=Ny=Nz=dim; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==8){ - nprocx=nprocy=nprocz=2; - Nx=Ny=Nz=dim; - nspheres=0; - Lx=Ly=Lz=1; - } - //....................................................................... - } - // ************************************************************** - // Broadcast simulation parameters from rank 0 to all other procs - MPI_Barrier(comm); - //................................................. - MPI_Bcast(&Nx,1,MPI_INT,0,comm); - MPI_Bcast(&Ny,1,MPI_INT,0,comm); - MPI_Bcast(&Nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); - // ************************************************************** - // ************************************************************** - - if (nprocs != nprocx*nprocy*nprocz){ - printf("nprocx = %i \n",nprocx); - printf("nprocy = %i \n",nprocy); - printf("nprocz = %i \n",nprocz); - INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!"); - } - - if (rank==0){ - printf("********************************************************\n"); - printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz); - printf("********************************************************\n"); - } - - MPI_Barrier(comm); double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; int BoundaryCondition=0; @@ -190,7 +103,7 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Np]; - ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np); + ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np,1); MPI_Barrier(comm); //......................device distributions................................. diff --git a/tests/TestWideHalo.cpp b/tests/TestWideHalo.cpp index cc29a15d..767aeaeb 100644 --- a/tests/TestWideHalo.cpp +++ b/tests/TestWideHalo.cpp @@ -8,7 +8,7 @@ #include #include "common/ScaLBL.h" #include "common/WideHalo.h" -#include "common/MPI_Helpers.h" +#include "common/MPI.h" using namespace std; @@ -20,105 +20,55 @@ int main(int argc, char **argv) // ***** MPI STUFF **************** //***************************************** // Initialize MPI - int rank,nprocs; - MPI_Init(&argc,&argv); - MPI_Comm comm = MPI_COMM_WORLD; - MPI_Comm_rank(comm,&rank); - MPI_Comm_size(comm,&nprocs); - int check; + Utilities::startup( argc, argv ); + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + int check=0; { - // parallel domain size (# of sub-domains) - int nprocx,nprocy,nprocz; - int iproc,jproc,kproc; - if (rank == 0){ printf("********************************************************\n"); printf("Running Color Model: TestColor \n"); printf("********************************************************\n"); } - - // BGK Model parameters - string FILENAME; - unsigned int nBlocks, nthreads; - int timestepMax, interval; - double Fx,Fy,Fz,tol; // Domain variables + int nprocx, nprocy, nprocz; double Lx,Ly,Lz; - int nspheres; int Nx,Ny,Nz; int i,j,k,n; - int dim = 3; - //if (rank == 0) printf("dim=%d\n",dim); - int timestep = 0; - int timesteps = 100; - int centralNode = 2; + int dim = 16; + Lx = Ly = Lz = 1.0; + int BoundaryCondition=0; - double tauA = 1.0; - double tauB = 1.0; - double rhoA = 1.0; - double rhoB = 1.0; - double alpha = 0.005; - double beta = 0.95; - - double tau = 1.0; - double mu=(tau-0.5)/3.0; - double rlx_setA=1.0/tau; - double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); - - Fx = Fy = 0.f; - Fz = 0.f; - - if (rank==0){ - //....................................................................... - // Reading the domain information file - //....................................................................... - if (nprocs==1){ - nprocx=nprocy=nprocz=1; - Nx=Ny=Nz=3; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==2){ - nprocx=2; nprocy=1; - nprocz=1; - Nx=Ny=Nz=dim; - Nx = dim; Ny = dim; Nz = dim; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==4){ - nprocx=nprocy=2; - nprocz=1; - Nx=Ny=Nz=dim; - nspheres=0; - Lx=Ly=Lz=1; - } - else if (nprocs==8){ - nprocx=nprocy=nprocz=2; - Nx=Ny=Nz=dim; - nspheres=0; - Lx=Ly=Lz=1; - } - //....................................................................... + //....................................................................... + // Reading the domain information file + //....................................................................... + nprocx=nprocy=nprocz=1; + if (nprocs==1){ + nprocx=nprocy=nprocz=1; + Nx=Ny=Nz=dim; + Lx=Ly=Lz=1; } - // ************************************************************** - // Broadcast simulation parameters from rank 0 to all other procs - MPI_Barrier(comm); - //................................................. - MPI_Bcast(&Nx,1,MPI_INT,0,comm); - MPI_Bcast(&Ny,1,MPI_INT,0,comm); - MPI_Bcast(&Nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - //................................................. - MPI_Barrier(comm); - // ************************************************************** + else if (nprocs==2){ + nprocx=2; nprocy=1; + nprocz=1; + Nx=Ny=Nz=dim; + Nx = dim; Ny = dim; Nz = dim; + Lx=Ly=Lz=1; + } + else if (nprocs==4){ + nprocx=nprocy=2; + nprocz=1; + Nx=Ny=Nz=dim; + Lx=Ly=Lz=1; + } + else if (nprocs==8){ + nprocx=nprocy=nprocz=2; + Nx=Ny=Nz=dim; + Lx=Ly=Lz=1; + } + //....................................................................... // ************************************************************** if (nprocs != nprocx*nprocy*nprocz){ @@ -134,10 +84,7 @@ int main(int argc, char **argv) printf("********************************************************\n"); } - MPI_Barrier(comm); - - double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; - int BoundaryCondition=0; + comm.barrier(); std::shared_ptr Dm = std::shared_ptr(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition)); Nx += 2; @@ -162,7 +109,7 @@ int main(int argc, char **argv) } } Dm->CommInit(); - MPI_Barrier(comm); + comm.barrier(); if (rank == 0) cout << "Domain set." << endl; if (rank==0) printf ("Create ScaLBL_Communicator \n"); @@ -179,12 +126,8 @@ int main(int argc, char **argv) IntArray Map(Nx,Ny,Nz); neighborList= new int[18*Np]; - ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,2); - MPI_Barrier(comm); - - //......................device distributions................................. - int dist_mem_size = Np*sizeof(double); - if (rank==0) printf ("Allocating distributions \n"); + ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,2); + comm.barrier(); int *NeighborList; int *dvcMap; @@ -241,8 +184,10 @@ int main(int argc, char **argv) CY=COLORGRAD[Np+idx]; CZ=COLORGRAD[2*Np+idx]; double error=sqrt((CX-1.0)*(CX-1.0)+(CY-2.0)*(CY-2.0)+ (CZ-3.0)*(CZ-3.0)); - if (error > 1e-8) + if (error > 1e-8){ + check++; printf("i,j,k=%i,%i,%i: Color gradient=%f,%f,%f \n",i,j,k,CX,CY,CZ); + } } } } @@ -250,8 +195,8 @@ int main(int argc, char **argv) } // **************************************************** - MPI_Barrier(comm); - MPI_Finalize(); + comm.barrier(); + Utilities::shutdown(); // **************************************************** return check;