OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Validated Bubble Test

Browse Source
This commit is contained in:
James McClure 2014-03-06 12:50:13 -05:00
parent 384edd7f58
commit 81764dfd3f
6 changed files with 726 additions and 1013 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
cpu/Color.cpp
View File

@ -756,10 +756,13 @@ extern "C" void dvc_ColorCollide( char *ID, double *disteven, double *distodd, d
distodd[8*N+n] = f17;
//...Store the Velocity..........................
Velocity[3*n] = jx;
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....................
*/ //...Store the Color Gradient....................
// ColorGrad[3*n] = nx*C;
// ColorGrad[3*n+1] = ny*C;
// ColorGrad[3*n+2] = nz*C;
@ -1162,9 +1165,9 @@ extern "C" void dvc_ColorCollideOpt( char *ID, double *disteven, double *distodd
distodd[7*N+n] = f15;
distodd[8*N+n] = f17;
//...Store the Velocity..........................
Velocity[3*n] = jx;
Velocity[3*n+1] = jy;
Velocity[3*n+2] = jz;
Velocity[n] = jx;
Velocity[N+n] = jy;
Velocity[2*N+n] = jz;
//***************************************************************
} // check if n is in the solid
} // loop over n
@ -1201,9 +1204,9 @@ extern "C" void dvc_MassColorCollideD3Q7(char *ID, double *A_even, double *A_odd
ny = ny/C;
nz = nz/C;
//....Load the flow velocity...........
ux = Velocity[3*n];
uy = Velocity[3*n+1];
uz = Velocity[3*n+2];
ux = Velocity[n];
uy = Velocity[N+n];
uz = Velocity[2*N+n];
//........................................................................
// READ THE DISTRIBUTIONS
// (read from opposite array due to previous swap operation)
@ -1315,9 +1318,9 @@ extern "C" void dvc_DensityStreamD3Q7(char *ID, double *Den, double *Copy, doubl
ny = ny/C;
nz = nz/C;
//....Load the flow velocity...........
ux = Velocity[3*n];
uy = Velocity[3*n+1];
uz = Velocity[3*n+2];
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

322
cpu/lb2_Color_wia_mpi.cpp
View File

@ -1016,10 +1016,11 @@ int main(int argc, char **argv)
cDistOdd = new double[9*N];
// data needed to perform CPU-based averaging
double *Vel,*Press;
Vel = new double[3*N]; // fluid velocity
Press = new double[N]; // fluid pressure
// double *Vel;
// Vel = new double[3*N]; // fluid velocity
// Press = new double[N]; // fluid pressure
DoubleArray Press(Nx,Ny,Nz);
DoubleArray MeanCurvature(Nx,Ny,Nz);
DoubleArray GaussCurvature(Nx,Ny,Nz);
DoubleArray SignDist_x(Nx,Ny,Nz); // Gradient of the signed distance
@ -1028,7 +1029,9 @@ int main(int argc, char **argv)
DoubleArray Phase_x(Nx,Ny,Nz); // Gradient of the phase indicator field
DoubleArray Phase_y(Nx,Ny,Nz);
DoubleArray Phase_z(Nx,Ny,Nz);
DoubleArray Vel_x(Nx,Ny,Nz); // Gradient of the phase indicator field
DoubleArray Vel_y(Nx,Ny,Nz);
DoubleArray Vel_z(Nx,Ny,Nz);
/*****************************************************************
VARIABLES FOR THE PMMC ALGORITHM
****************************************************************** */
@ -1276,6 +1279,10 @@ int main(int argc, char **argv)
dvc_UnpackValues(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, Phi, N);
//...................................................................................
if (rank==0 && pBC){
printf("Setting inlet pressure = %f \n", din);
printf("Setting outlet pressure = %f \n", dout);
}
if (pBC && kproc == 0) {
dvc_PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz,S);
dvc_ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz,S);
@ -1299,8 +1306,10 @@ int main(int argc, char **argv)
dvc_Barrier();
dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S);
dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
dvc_CopyToHost(Press,Pressure,N*sizeof(double));
dvc_CopyToHost(Vel,Velocity,3*N*sizeof(double));
dvc_CopyToHost(Press.data,Pressure,N*sizeof(double));
dvc_CopyToHost(Vel_x.data,&Velocity[0],N*sizeof(double));
dvc_CopyToHost(Vel_y.data,&Velocity[N],N*sizeof(double));
dvc_CopyToHost(Vel_z.data,&Velocity[2*N],N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
//...........................................................................
@ -1723,94 +1732,14 @@ int main(int argc, char **argv)
//...................................................................................
// ZeroHalo(Den,Nx,Ny,Nz);
// ZeroHalo(Copy,Nx,Ny,Nz);
if (pBC && kproc == 0) {
dvc_PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz,S);
dvc_ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz,S);
// Fill the inlet with component a
/* for (k=0; k<1; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = 1.0;
}
}
}
for (k=1; k<3; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = 1.0;
Den[n] = 1.0;
Den[N+n] = 0.0;
A_even[n] = 0.3333333333333333;
A_odd[n] = 0.1111111111111111;
A_even[N+n] = 0.1111111111111111;
A_odd[N+n] = 0.1111111111111111;
A_even[2*N+n] = 0.1111111111111111;
A_odd[2*N+n] = 0.1111111111111111;
A_even[3*N+n] = 0.1111111111111111;
B_even[n] = 0.0;
B_odd[n] = 0.0;
B_even[N+n] = 0.0;
B_odd[N+n] = 0.0;
B_even[2*N+n] = 0.0;
B_odd[2*N+n] = 0.0;
B_even[3*N+n] = 0.0;
}
}
}
*/
}
if (pBC && kproc == nprocz-1){
dvc_PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,S,Nx*Ny*(Nz-2));
dvc_ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz,S);
/* // Fill the outlet with component b
for (k=Nz-3; k<Nz-1; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = -1.0;
Den[n] = 0.0;
Den[N+n] = 1.0;
A_even[n] = 0.0;
A_odd[n] = 0.0;
A_even[N+n] = 0.0;
A_odd[N+n] = 0.0;
A_even[2*N+n] = 0.0;
A_odd[2*N+n] = 0.0;
A_even[3*N+n] = 0.0;
B_even[n] = 0.3333333333333333;
B_odd[n] = 0.1111111111111111;
B_even[N+n] = 0.1111111111111111;
B_odd[N+n] = 0.1111111111111111;
B_even[2*N+n] = 0.1111111111111111;
B_odd[2*N+n] = 0.1111111111111111;
B_even[3*N+n] = 0.1111111111111111;
}
}
}
for (k=Nz-1; k<Nz; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = -1.0;
}
}
}
*/
}
//...................................................................................
@ -1836,8 +1765,10 @@ int main(int argc, char **argv)
dvc_Barrier();
dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S);
dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
dvc_CopyToHost(Press,Pressure,N*sizeof(double));
dvc_CopyToHost(Vel,Velocity,3*N*sizeof(double));
dvc_CopyToHost(Press.data,Pressure,N*sizeof(double));
dvc_CopyToHost(Vel_x.data,&Velocity[0],N*sizeof(double));
dvc_CopyToHost(Vel_y.data,&Velocity[N],N*sizeof(double));
dvc_CopyToHost(Vel_z.data,&Velocity[2*N],N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
}
if (timestep%1000 == 5){
@ -1859,6 +1790,100 @@ int main(int argc, char **argv)
//...........................................................................
// Fill in the halo region for the mesh gradients and curvature
//...........................................................................
// Pressure
//...........................................................................
CommunicateMeshHalo(Vel_x, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
// Velocity
//...........................................................................
CommunicateMeshHalo(Press, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
//...........................................................................
CommunicateMeshHalo(Vel_y, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
//...........................................................................
CommunicateMeshHalo(Vel_z, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
// Mean Curvature
//...........................................................................
CommunicateMeshHalo(MeanCurvature, MPI_COMM_WORLD,
@ -2022,11 +2047,11 @@ int main(int argc, char **argv)
// volume the excludes the interfacial region
vol_n += 0.125;
// pressure
pan += 0.125*Press[n];
pan += 0.125*Press.data[n];
// velocity
van(0) += 0.125*Vel[3*n];
van(1) += 0.125*Vel[3*n+1];
van(2) += 0.125*Vel[3*n+2];
van(0) += 0.125*Vel_x.data[n];
van(1) += 0.125*Vel_y.data[n];
van(2) += 0.125*Vel_z.data[n];
}
// volume averages over the wetting phase
@ -2034,11 +2059,11 @@ int main(int argc, char **argv)
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
paw += 0.125*Press[n];
paw += 0.125*Press.data[n];
// velocity
vaw(0) += 0.125*Vel[3*n];
vaw(1) += 0.125*Vel[3*n+1];
vaw(2) += 0.125*Vel[3*n+2];
vaw(0) += 0.125*Vel_x.data[n];
vaw(1) += 0.125*Vel_y.data[n];
vaw(2) += 0.125*Vel_z.data[n];
}
}
}
@ -2202,104 +2227,6 @@ int main(int argc, char **argv)
// Write out the phase indicator field
//************************************************************************/
// printf("Local File Name = %s \n",LocalRankFilename);
// dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
//!!!!!!!!DEBUG HERE!!!!!!!!!!!!!
/* dvc_InitDenColorDistance(ID, Den, Phi, SignDist.data, das, dbs, beta, xIntPos, Nx, Ny, Nz, S);
dvc_InitD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz, S);
dvc_InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz, S);
dvc_ComputeDensityD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz, S);
dvc_ComputeDensityD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz, S);
dvc_ComputePhi(ID, Phi, Den, N, S);
if (pBC && kproc == 0) {
dvc_PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz,S);
// Fill the inlet with component a
for (k=0; k<1; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = 1.0;
}
}
}
for (k=1; k<4; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = 1.0;
Den[n] = 1.0;
Den[N+n] = 0.0;
A_even[n] = 0.3333333333333333;
A_odd[n] = 0.1111111111111111;
A_even[N+n] = 0.1111111111111111;
A_odd[N+n] = 0.1111111111111111;
A_even[2*N+n] = 0.1111111111111111;
A_odd[2*N+n] = 0.1111111111111111;
A_even[3*N+n] = 0.1111111111111111;
B_even[n] = 0.0;
B_odd[n] = 0.0;
B_even[N+n] = 0.0;
B_odd[N+n] = 0.0;
B_even[2*N+n] = 0.0;
B_odd[2*N+n] = 0.0;
B_even[3*N+n] = 0.0;
}
}
}
}
if (pBC && kproc == nprocz-1){
dvc_PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,S,Nx*Ny*(Nz-2));
// Fill the outlet with component b
for (k=Nz-4; k<Nz-1; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = -1.0;
Den[n] = 0.0;
Den[N+n] = 1.0;
A_even[n] = 0.0;
A_odd[n] = 0.0;
A_even[N+n] = 0.0;
A_odd[N+n] = 0.0;
A_even[2*N+n] = 0.0;
A_odd[2*N+n] = 0.0;
A_even[3*N+n] = 0.0;
B_even[n] = 0.3333333333333333;
B_odd[n] = 0.1111111111111111;
B_even[N+n] = 0.1111111111111111;
B_odd[N+n] = 0.1111111111111111;
B_even[2*N+n] = 0.1111111111111111;
B_odd[2*N+n] = 0.1111111111111111;
B_even[3*N+n] = 0.1111111111111111;
}
}
}
for (k=Nz-1; k<Nz; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = -1.0;
}
}
}
}
dvc_SwapD3Q7(ID, A_even, A_odd, Nx, Ny, Nz, S);
dvc_SwapD3Q7(ID, B_even, B_odd, Nx, Ny, Nz, S);
dvc_ComputeDensityD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz, S);
dvc_ComputeDensityD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz, S);
*/
//#ifdef WriteOutput
@ -2312,11 +2239,11 @@ int main(int argc, char **argv)
fclose(PHASE);
//#endif
dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S);
dvc_CopyToHost(Press,Pressure,N*sizeof(double));
dvc_CopyToHost(Press.data,Pressure,N*sizeof(double));
sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
FILE *PRESS;
PRESS = fopen(LocalRankFilename,"wb");
fwrite(Press,8,N,PRESS);
fwrite(Press.data,8,N,PRESS);
fclose(PRESS);
/* sprintf(LocalRankFilename,"%s%s","dPdt.",LocalRankString);
@ -2324,7 +2251,6 @@ int main(int argc, char **argv)
SPEED = fopen(LocalRankFilename,"wb");
fwrite(dPdt.data,8,N,SPEED);
fclose(SPEED);
*/
sprintf(LocalRankFilename,"%s%s","DenA.",LocalRankString);
FILE *DENA;
@ -2338,7 +2264,7 @@ int main(int argc, char **argv)
fwrite(&Den[N],8,N,DENB);
fclose(DENB);
/* sprintf(LocalRankFilename,"%s%s","GradMag.",LocalRankString);
sprintf(LocalRankFilename,"%s%s","GradMag.",LocalRankString);
FILE *GRAD;
GRAD = fopen(LocalRankFilename,"wb");
for (k=0; k<Nz; k++){

4
example/Bubble/Color.in
View File

@ -1,6 +1,6 @@
1.0
1.0e-2 0.95 0.1 0.9
1.0e-2 0.95 0.8
0.7
0.0 0.0 0.0
0 0 1.0 1.0
20000 1000 1e-5
200 1000 1e-5

2
example/Bubble/Domain.in
View File

@ -1,4 +1,4 @@
1 1 1
80 80 80
229
0
1.0 1.0 1.0

46
example/Bubble/ExampleOutput.txt Normal file
View File

@ -0,0 +1,46 @@
********************************************************
Running Hybrid Implementation of Color LBM
********************************************************
********************************************************
tau = 1.000000
alpha = 0.010000
beta = 0.950000
das = 0.100000
dbs = 0.900000
Value of phi at solid surface = 0.800000
Distance to phi = 0.0: 1.156434
gamma_{wn} = 0.057960
Force(x) = 0.000000
Force(y) = 0.000000
Force(z) = 0.000000
Sub-domain size = 80 x 80 x 80
Parallel domain size = 1 x 1 x 1
********************************************************
Number of blocks = 32
Threads per block = 128
Sweeps per thread = 135
Number of nodes per side = 82
Total Number of nodes = 551368
********************************************************
Read input media...
Setting up communication control structures
Preparing the sendlists
SendLists are ready on host
Prepare to copy send/recv Lists to device
Devices are ready to communicate.
Copying phase ID to device
Allocating distributions
********************************************************
No. of timesteps: 1000
--------------------------------------------------------------------------------------
radius sw pw pn awn Jwn Gwn [xx, yy, zz, xy, xz, yz] --------------------------------------------------------------------------------------
8 0.34906 0.33333 0.0013898 0.26946 0.33334 0.33333 0.33333 -5.3309e-06 1.1201e-05 -1.3209e-05
10 0.34568 0.33333 0.0022109 0.21364 0.33334 0.33333 0.33333 -4.1967e-06 1.0265e-05 -1.2941e-05
12 0.34355 0.33333 0.0031972 0.1779 0.33334 0.33333 0.33333 -5.6086e-06 1.1842e-05 -1.4445e-05
15 0.34143 0.33332 0.005043 0.14164 0.33334 0.33333 0.33333 -5.5737e-06 1.0229e-05 -1.2561e-05
-------------------------------------------------------------------
********************************************************
CPU time = inf
Lattice update rate (per core)= 0.000000 MLUPS
Lattice update rate (total)= 0.000000 MLUPS
********************************************************

764
tests/TestBubble.cpp
View File

@ -318,157 +318,17 @@ int main(int argc, char **argv)
double sum_local;
double iVol_global = 1.0/(1.0*Nx*Ny*Nz*nprocs);
double porosity;
/* //.......................................................................
ifstream PM(LocalRankFilename,ios::binary);
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 0;
}
}
}
for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
PM.read((char *) (&value), sizeof(value));
n = k*Nx*Ny+j*Nx+i;
id[n] = value;
if (value > 0) sum++;
}
}
}
PM.close();
// printf("File porosity = %f\n", double(sum)/N);
*/ //...........................................................................
// double *SignDist;
// SignDist = new double[N];
DoubleArray SignDist(Nx,Ny,Nz);
/* //.......................................................................
#ifdef CBUB
// Initializes a constrained bubble test
double BubbleBot = 20.0; // How big to make the NWP bubble
double BubbleTop = 60.0; // How big to make the NWP bubble
double TubeRadius = 15.5; // Radius of the capillary tube
//.......................................................................
double BubbleRadius = 15.5; // Radius of the capillary tube
sum=0;
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
// Cylindrical capillary tube aligned with the z direction
SignDist(i,j,k) = TubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)
+ (j-Ny/2)*(j-Ny/2)));
// Initialize phase positions field
if (SignDist(i,j,k) < 0.0){
id[n] = 0;
}
else if (k<BubbleBot){
id[n] = 2;
sum++;
}
else if (k<BubbleTop && rank == 0 && pBC == 0){
id[n] = 1;
sum++;
}
else{
id[n] = 2;
sum++;
}
}
}
}
porosity = double(sum)/double(1.0*N);
#else
// Read in sphere pack
if (rank==1) printf("nspheres =%i \n",nspheres);
//.......................................................................
double *cx,*cy,*cz,*rad;
cx = new double[nspheres];
cy = new double[nspheres];
cz = new double[nspheres];
rad = new double[nspheres];
//.......................................................................
if (rank == 0) printf("Reading the sphere packing \n");
if (rank == 0) ReadSpherePacking(nspheres,cx,cy,cz,rad);
MPI_Barrier(MPI_COMM_WORLD);
// Broadcast the sphere packing to all processes
MPI_Bcast(cx,nspheres,MPI_DOUBLE,0,MPI_COMM_WORLD);
MPI_Bcast(cy,nspheres,MPI_DOUBLE,0,MPI_COMM_WORLD);
MPI_Bcast(cz,nspheres,MPI_DOUBLE,0,MPI_COMM_WORLD);
MPI_Bcast(rad,nspheres,MPI_DOUBLE,0,MPI_COMM_WORLD);
//...........................................................................
MPI_Barrier(MPI_COMM_WORLD);
if (rank == 0) cout << "Domain set." << endl;
//.......................................................................
// sprintf(LocalRankString,"%05d",rank);
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
//.......................................................................
SignedDistance(SignDist.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
iproc,jproc,kproc,nprocx,nprocy,nprocz);
//.......................................................................
// Assign the phase ID field based on the signed distance
//.......................................................................
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 0;
}
}
}
sum=0;
for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (SignDist.data[n] > 0.0){
id[n] = 2;
sum++;
}
}
}
}
//.........................................................
// If pressure boundary conditions are applied remove solid
if (pBC && kproc == 0){
for (k=0; k<3; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 1;
}
}
}
}
if (pBC && kproc == nprocz-1){
for (k=Nz-3; k<Nz; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 2;
}
}
}
}
//.........................................................
sum_local = 1.0*sum;
MPI_Allreduce(&sum_local,&porosity,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
porosity = porosity*iVol_global;
if (rank==0) printf("Media porosity = %f \n",porosity);
// Generate the residual NWP
if (rank==0) printf("Initializing with NWP saturation = %f \n",wp_saturation);
// if (!pBC) GenerateResidual(id,Nx,Ny,Nz,wp_saturation);
#endif
*/
double BubbleRadius = 15.5;
sum=0;
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
// No Solid phase
SignDist(i,j,k) = 100;
// Initialize phase positions field
if (SignDist(i,j,k) < 0.0){
@ -482,7 +342,6 @@ int main(int argc, char **argv)
}
}
// Set up MPI communication structurese
if (rank==0) printf ("Setting up communication control structures \n");
//......................................................................................
@ -599,20 +458,89 @@ int main(int argc, char **argv)
//**********************************************************************************
// Fill in the recieve counts using MPI
sendtag = recvtag = 3;
CommunicateSendRecvCounts( MPI_COMM_WORLD, sendtag, recvtag,
rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z,
rank_xy, rank_XY, rank_xY, rank_Xy,
rank_xz, rank_XZ, rank_xZ, rank_Xz,
rank_yz, rank_YZ, rank_yZ, rank_Yz,
sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z,
sendCount_xy, sendCount_XY, sendCount_xY, sendCount_Xy,
sendCount_xz, sendCount_XZ, sendCount_xZ, sendCount_Xz,
sendCount_yz, sendCount_YZ, sendCount_yZ, sendCount_Yz,
recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z,
recvCount_xy, recvCount_XY, recvCount_xY, recvCount_Xy,
recvCount_xz, recvCount_XZ, recvCount_xZ, recvCount_Xz,
recvCount_yz, recvCount_YZ, recvCount_yZ, recvCount_Yz );
//**********************************************************************************
MPI_Isend(&sendCount_x, 1,MPI_INT,rank_x,sendtag,MPI_COMM_WORLD,&req1[0]);
MPI_Irecv(&recvCount_X, 1,MPI_INT,rank_X,recvtag,MPI_COMM_WORLD,&req2[0]);
MPI_Isend(&sendCount_X, 1,MPI_INT,rank_X,sendtag,MPI_COMM_WORLD,&req1[1]);
MPI_Irecv(&recvCount_x, 1,MPI_INT,rank_x,recvtag,MPI_COMM_WORLD,&req2[1]);
MPI_Isend(&sendCount_y, 1,MPI_INT,rank_y,sendtag,MPI_COMM_WORLD,&req1[2]);
MPI_Irecv(&recvCount_Y, 1,MPI_INT,rank_Y,recvtag,MPI_COMM_WORLD,&req2[2]);
MPI_Isend(&sendCount_Y, 1,MPI_INT,rank_Y,sendtag,MPI_COMM_WORLD,&req1[3]);
MPI_Irecv(&recvCount_y, 1,MPI_INT,rank_y,recvtag,MPI_COMM_WORLD,&req2[3]);
MPI_Isend(&sendCount_z, 1,MPI_INT,rank_z,sendtag,MPI_COMM_WORLD,&req1[4]);
MPI_Irecv(&recvCount_Z, 1,MPI_INT,rank_Z,recvtag,MPI_COMM_WORLD,&req2[4]);
MPI_Isend(&sendCount_Z, 1,MPI_INT,rank_Z,sendtag,MPI_COMM_WORLD,&req1[5]);
MPI_Irecv(&recvCount_z, 1,MPI_INT,rank_z,recvtag,MPI_COMM_WORLD,&req2[5]);
MPI_Isend(&sendCount_xy, 1,MPI_INT,rank_xy,sendtag,MPI_COMM_WORLD,&req1[6]);
MPI_Irecv(&recvCount_XY, 1,MPI_INT,rank_XY,recvtag,MPI_COMM_WORLD,&req2[6]);
MPI_Isend(&sendCount_XY, 1,MPI_INT,rank_XY,sendtag,MPI_COMM_WORLD,&req1[7]);
MPI_Irecv(&recvCount_xy, 1,MPI_INT,rank_xy,recvtag,MPI_COMM_WORLD,&req2[7]);
MPI_Isend(&sendCount_Xy, 1,MPI_INT,rank_Xy,sendtag,MPI_COMM_WORLD,&req1[8]);
MPI_Irecv(&recvCount_xY, 1,MPI_INT,rank_xY,recvtag,MPI_COMM_WORLD,&req2[8]);
MPI_Isend(&sendCount_xY, 1,MPI_INT,rank_xY,sendtag,MPI_COMM_WORLD,&req1[9]);
MPI_Irecv(&recvCount_Xy, 1,MPI_INT,rank_Xy,recvtag,MPI_COMM_WORLD,&req2[9]);
MPI_Isend(&sendCount_xz, 1,MPI_INT,rank_xz,sendtag,MPI_COMM_WORLD,&req1[10]);
MPI_Irecv(&recvCount_XZ, 1,MPI_INT,rank_XZ,recvtag,MPI_COMM_WORLD,&req2[10]);
MPI_Isend(&sendCount_XZ, 1,MPI_INT,rank_XZ,sendtag,MPI_COMM_WORLD,&req1[11]);
MPI_Irecv(&recvCount_xz, 1,MPI_INT,rank_xz,recvtag,MPI_COMM_WORLD,&req2[11]);
MPI_Isend(&sendCount_Xz, 1,MPI_INT,rank_Xz,sendtag,MPI_COMM_WORLD,&req1[12]);
MPI_Irecv(&recvCount_xZ, 1,MPI_INT,rank_xZ,recvtag,MPI_COMM_WORLD,&req2[12]);
MPI_Isend(&sendCount_xZ, 1,MPI_INT,rank_xZ,sendtag,MPI_COMM_WORLD,&req1[13]);
MPI_Irecv(&recvCount_Xz, 1,MPI_INT,rank_Xz,recvtag,MPI_COMM_WORLD,&req2[13]);
MPI_Isend(&sendCount_yz, 1,MPI_INT,rank_yz,sendtag,MPI_COMM_WORLD,&req1[14]);
MPI_Irecv(&recvCount_YZ, 1,MPI_INT,rank_YZ,recvtag,MPI_COMM_WORLD,&req2[14]);
MPI_Isend(&sendCount_YZ, 1,MPI_INT,rank_YZ,sendtag,MPI_COMM_WORLD,&req1[15]);
MPI_Irecv(&recvCount_yz, 1,MPI_INT,rank_yz,recvtag,MPI_COMM_WORLD,&req2[15]);
MPI_Isend(&sendCount_Yz, 1,MPI_INT,rank_Yz,sendtag,MPI_COMM_WORLD,&req1[16]);
MPI_Irecv(&recvCount_yZ, 1,MPI_INT,rank_yZ,recvtag,MPI_COMM_WORLD,&req2[16]);
MPI_Isend(&sendCount_yZ, 1,MPI_INT,rank_yZ,sendtag,MPI_COMM_WORLD,&req1[17]);
MPI_Irecv(&recvCount_Yz, 1,MPI_INT,rank_Yz,recvtag,MPI_COMM_WORLD,&req2[17]);
MPI_Waitall(18,req1,stat1);
MPI_Waitall(18,req2,stat2);
MPI_Barrier(MPI_COMM_WORLD);
/* MPI_Send(&sendCount_x,1,MPI_INT,rank_X,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_X,1,MPI_INT,rank_x,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_X,1,MPI_INT,rank_x,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_x,1,MPI_INT,rank_X,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_y,1,MPI_INT,rank_Y,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_Y,1,MPI_INT,rank_y,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_Y,1,MPI_INT,rank_y,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_y,1,MPI_INT,rank_Y,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_z,1,MPI_INT,rank_Z,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_Z,1,MPI_INT,rank_z,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_Z,1,MPI_INT,rank_z,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_z,1,MPI_INT,rank_Z,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_xy,1,MPI_INT,rank_XY,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_XY,1,MPI_INT,rank_xy,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_XY,1,MPI_INT,rank_xy,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_xy,1,MPI_INT,rank_XY,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_Xy,1,MPI_INT,rank_xY,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_xY,1,MPI_INT,rank_Xy,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_xY,1,MPI_INT,rank_Xy,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_Xy,1,MPI_INT,rank_xY,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_xz,1,MPI_INT,rank_XZ,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_XZ,1,MPI_INT,rank_xz,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_XZ,1,MPI_INT,rank_xz,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_xz,1,MPI_INT,rank_XZ,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_Xz,1,MPI_INT,rank_xZ,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_xZ,1,MPI_INT,rank_Xz,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_xZ,1,MPI_INT,rank_Xz,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_Xz,1,MPI_INT,rank_xZ,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_yz,1,MPI_INT,rank_YZ,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_YZ,1,MPI_INT,rank_yz,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_YZ,1,MPI_INT,rank_yz,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_yz,1,MPI_INT,rank_YZ,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_Yz,1,MPI_INT,rank_yZ,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_yZ,1,MPI_INT,rank_Yz,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Send(&sendCount_yZ,1,MPI_INT,rank_Yz,sendtag,MPI_COMM_WORLD);
MPI_Recv(&recvCount_Yz,1,MPI_INT,rank_yZ,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Barrier(MPI_COMM_WORLD);
*/ //**********************************************************************************
//......................................................................................
int *recvList_x, *recvList_y, *recvList_z, *recvList_X, *recvList_Y, *recvList_Z;
int *recvList_xy, *recvList_yz, *recvList_xz, *recvList_Xy, *recvList_Yz, *recvList_xZ;
@ -642,25 +570,48 @@ int main(int argc, char **argv)
// Use MPI to fill in the appropriate values for recvList
// Fill in the recieve lists using MPI
sendtag = recvtag = 4;
CommunicateRecvLists( MPI_COMM_WORLD, sendtag, recvtag,
sendList_x, sendList_y, sendList_z, sendList_X, sendList_Y, sendList_Z,
sendList_xy, sendList_XY, sendList_xY, sendList_Xy,
sendList_xz, sendList_XZ, sendList_xZ, sendList_Xz,
sendList_yz, sendList_YZ, sendList_yZ, sendList_Yz,
sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z,
sendCount_xy, sendCount_XY, sendCount_xY, sendCount_Xy,
sendCount_xz, sendCount_XZ, sendCount_xZ, sendCount_Xz,
sendCount_yz, sendCount_YZ, sendCount_yZ, sendCount_Yz,
recvList_x, recvList_y, recvList_z, recvList_X, recvList_Y, recvList_Z,
recvList_xy, recvList_XY, recvList_xY, recvList_Xy,
recvList_xz, recvList_XZ, recvList_xZ, recvList_Xz,
recvList_yz, recvList_YZ, recvList_yZ, recvList_Yz,
recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z,
recvCount_xy, recvCount_XY, recvCount_xY, recvCount_Xy,
recvCount_xz, recvCount_XZ, recvCount_xZ, recvCount_Xz,
recvCount_yz, recvCount_YZ, recvCount_yZ, recvCount_Yz,
rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z, rank_xy, rank_XY, rank_xY,
rank_Xy, rank_xz, rank_XZ, rank_xZ, rank_Xz, rank_yz, rank_YZ, rank_yZ, rank_Yz );
MPI_Isend(sendList_x, sendCount_x,MPI_INT,rank_x,sendtag,MPI_COMM_WORLD,&req1[0]);
MPI_Irecv(recvList_X, recvCount_X,MPI_INT,rank_X,recvtag,MPI_COMM_WORLD,&req2[0]);
MPI_Isend(sendList_X, sendCount_X,MPI_INT,rank_X,sendtag,MPI_COMM_WORLD,&req1[1]);
MPI_Irecv(recvList_x, recvCount_x,MPI_INT,rank_x,recvtag,MPI_COMM_WORLD,&req2[1]);
MPI_Isend(sendList_y, sendCount_y,MPI_INT,rank_y,sendtag,MPI_COMM_WORLD,&req1[2]);
MPI_Irecv(recvList_Y, recvCount_Y,MPI_INT,rank_Y,recvtag,MPI_COMM_WORLD,&req2[2]);
MPI_Isend(sendList_Y, sendCount_Y,MPI_INT,rank_Y,sendtag,MPI_COMM_WORLD,&req1[3]);
MPI_Irecv(recvList_y, recvCount_y,MPI_INT,rank_y,recvtag,MPI_COMM_WORLD,&req2[3]);
MPI_Isend(sendList_z, sendCount_z,MPI_INT,rank_z,sendtag,MPI_COMM_WORLD,&req1[4]);
MPI_Irecv(recvList_Z, recvCount_Z,MPI_INT,rank_Z,recvtag,MPI_COMM_WORLD,&req2[4]);
MPI_Isend(sendList_Z, sendCount_Z,MPI_INT,rank_Z,sendtag,MPI_COMM_WORLD,&req1[5]);
MPI_Irecv(recvList_z, recvCount_z,MPI_INT,rank_z,recvtag,MPI_COMM_WORLD,&req2[5]);
MPI_Isend(sendList_xy, sendCount_xy,MPI_INT,rank_xy,sendtag,MPI_COMM_WORLD,&req1[6]);
MPI_Irecv(recvList_XY, recvCount_XY,MPI_INT,rank_XY,recvtag,MPI_COMM_WORLD,&req2[6]);
MPI_Isend(sendList_XY, sendCount_XY,MPI_INT,rank_XY,sendtag,MPI_COMM_WORLD,&req1[7]);
MPI_Irecv(recvList_xy, recvCount_xy,MPI_INT,rank_xy,recvtag,MPI_COMM_WORLD,&req2[7]);
MPI_Isend(sendList_Xy, sendCount_Xy,MPI_INT,rank_Xy,sendtag,MPI_COMM_WORLD,&req1[8]);
MPI_Irecv(recvList_xY, recvCount_xY,MPI_INT,rank_xY,recvtag,MPI_COMM_WORLD,&req2[8]);
MPI_Isend(sendList_xY, sendCount_xY,MPI_INT,rank_xY,sendtag,MPI_COMM_WORLD,&req1[9]);
MPI_Irecv(recvList_Xy, recvCount_Xy,MPI_INT,rank_Xy,recvtag,MPI_COMM_WORLD,&req2[9]);
MPI_Isend(sendList_xz, sendCount_xz,MPI_INT,rank_xz,sendtag,MPI_COMM_WORLD,&req1[10]);
MPI_Irecv(recvList_XZ, recvCount_XZ,MPI_INT,rank_XZ,recvtag,MPI_COMM_WORLD,&req2[10]);
MPI_Isend(sendList_XZ, sendCount_XZ,MPI_INT,rank_XZ,sendtag,MPI_COMM_WORLD,&req1[11]);
MPI_Irecv(recvList_xz, recvCount_xz,MPI_INT,rank_xz,recvtag,MPI_COMM_WORLD,&req2[11]);
MPI_Isend(sendList_Xz, sendCount_Xz,MPI_INT,rank_Xz,sendtag,MPI_COMM_WORLD,&req1[12]);
MPI_Irecv(recvList_xZ, recvCount_xZ,MPI_INT,rank_xZ,recvtag,MPI_COMM_WORLD,&req2[12]);
MPI_Isend(sendList_xZ, sendCount_xZ,MPI_INT,rank_xZ,sendtag,MPI_COMM_WORLD,&req1[13]);
MPI_Irecv(recvList_Xz, recvCount_Xz,MPI_INT,rank_Xz,recvtag,MPI_COMM_WORLD,&req2[13]);
MPI_Isend(sendList_yz, sendCount_yz,MPI_INT,rank_yz,sendtag,MPI_COMM_WORLD,&req1[14]);
MPI_Irecv(recvList_YZ, recvCount_YZ,MPI_INT,rank_YZ,recvtag,MPI_COMM_WORLD,&req2[14]);
MPI_Isend(sendList_YZ, sendCount_YZ,MPI_INT,rank_YZ,sendtag,MPI_COMM_WORLD,&req1[15]);
MPI_Irecv(recvList_yz, recvCount_yz,MPI_INT,rank_yz,recvtag,MPI_COMM_WORLD,&req2[15]);
MPI_Isend(sendList_Yz, sendCount_Yz,MPI_INT,rank_Yz,sendtag,MPI_COMM_WORLD,&req1[16]);
MPI_Irecv(recvList_yZ, recvCount_yZ,MPI_INT,rank_yZ,recvtag,MPI_COMM_WORLD,&req2[16]);
MPI_Isend(sendList_yZ, sendCount_yZ,MPI_INT,rank_yZ,sendtag,MPI_COMM_WORLD,&req1[17]);
MPI_Irecv(recvList_Yz, recvCount_Yz,MPI_INT,rank_Yz,recvtag,MPI_COMM_WORLD,&req2[17]);
MPI_Waitall(18,req1,stat1);
MPI_Waitall(18,req2,stat2);
MPI_Barrier(MPI_COMM_WORLD);
//......................................................................................
for (int idx=0; idx<recvCount_x; idx++) recvList_x[idx] -= (Nx-2);
for (int idx=0; idx<recvCount_X; idx++) recvList_X[idx] += (Nx-2);
@ -995,6 +946,8 @@ int main(int argc, char **argv)
dvc_CopyToDevice(ID, id, N);
//...........................................................................
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
@ -1022,6 +975,7 @@ int main(int argc, char **argv)
dvc_AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
dvc_AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
//...........................................................................
//...........................................................................
// Phase indicator (in array form as needed by PMMC algorithm)
DoubleArray Phase(Nx,Ny,Nz);
@ -1037,9 +991,8 @@ int main(int argc, char **argv)
cDistOdd = new double[9*N];
// data needed to perform CPU-based averaging
double *Vel,*Press;
double *Vel;
Vel = new double[3*N]; // fluid velocity
Press = new double[N]; // fluid pressure
DoubleArray MeanCurvature(Nx,Ny,Nz);
DoubleArray GaussCurvature(Nx,Ny,Nz);
@ -1049,6 +1002,7 @@ int main(int argc, char **argv)
DoubleArray Phase_x(Nx,Ny,Nz); // Gradient of the phase indicator field
DoubleArray Phase_y(Nx,Ny,Nz);
DoubleArray Phase_z(Nx,Ny,Nz);
DoubleArray Press(Nx,Ny,Nz);
/*****************************************************************
VARIABLES FOR THE PMMC ALGORITHM
@ -1141,9 +1095,9 @@ int main(int argc, char **argv)
int nc=0;
//...........................................................................
// Set up the cube list (very regular in this case due to lack of blob-ID)
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
for (k=0; k<Nz-2; k++){
for (j=0; j<Ny-2; j++){
for (i=0; i<Nx-2; i++){
cubeList(0,nc) = i;
cubeList(1,nc) = j;
cubeList(2,nc) = k;
@ -1159,6 +1113,21 @@ int main(int argc, char **argv)
faceGrid=Nx*Ny/packThreads;
//...........................................................................
//...........................................................................
int timestep = 0;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
//.......create a stream for the LB calculation.......
// cudaStream_t stream;
// cudaStreamCreate(&stream);
//.......create and start timer............
double starttime,stoptime,cputime;
MPI_Barrier(MPI_COMM_WORLD);
starttime = MPI_Wtime();
//.........................................
//...........................................................................
// MAIN VARIABLES INITIALIZED HERE
//...........................................................................
@ -1200,8 +1169,8 @@ int main(int argc, char **argv)
}
// Copy the bubble to the device and initialize
dvc_CopyToDevice(ID, id, N);
//...........................................................................
if (rank==0) printf("Setting the distributions, size = %i\n", N);
//...........................................................................
dvc_InitD3Q19(ID, f_even, f_odd, Nx, Ny, Nz, S);
//......................................................................
@ -1335,6 +1304,10 @@ int main(int argc, char **argv)
dvc_UnpackValues(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, Phi, N);
//...................................................................................
if (rank==0 && pBC){
printf("Setting inlet pressure = %f \n", din);
printf("Setting outlet pressure = %f \n", dout);
}
if (pBC && kproc == 0) {
dvc_PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz,S);
dvc_ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz,S);
@ -1358,45 +1331,16 @@ int main(int argc, char **argv)
dvc_Barrier();
dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S);
dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
dvc_CopyToHost(Press,Pressure,N*sizeof(double));
dvc_CopyToHost(Vel,Velocity,3*N*sizeof(double));
dvc_CopyToHost(Press.data,Pressure,N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
//...........................................................................
int timestep = 0;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
//.......create a stream for the LB calculation.......
// cudaStream_t stream;
// cudaStreamCreate(&stream);
//.......create and start timer............
double starttime,stoptime,cputime;
MPI_Barrier(MPI_COMM_WORLD);
starttime = MPI_Wtime();
//.........................................
sendtag = recvtag = 5;
timestepMax =10005;
printf("Hard-coding timestepMax = 10,000 for bubble test \n");
//************ MAIN ITERATION LOOP ***************************************/
while (timestep < timestepMax){
//*************************************************************************
// Compute the color gradient
//*************************************************************************
//dvc_ComputeColorGradient(nBlocks, nthreads, S,
// ID, Phi, ColorGrad, Nx, Ny, Nz);
//*************************************************************************
//*************************************************************************
// Perform collision step for the momentum transport
//*************************************************************************
// dvc_ColorCollide(nBlocks, nthreads, S, ID, f_even, f_odd, ColorGrad, Velocity,
// rlxA, rlxB,alpha, beta, Fx, Fy, Fz, Nx, Ny, Nz, pBC);
//*************************************************************************
//*************************************************************************
// Fused Color Gradient and Collision
@ -1773,95 +1717,14 @@ int main(int argc, char **argv)
dvc_UnpackValues(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, Phi, N);
//...................................................................................
// ZeroHalo(Den,Nx,Ny,Nz);
// ZeroHalo(Copy,Nx,Ny,Nz);
if (pBC && kproc == 0) {
dvc_PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz,S);
dvc_ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz,S);
// Fill the inlet with component a
/* for (k=0; k<1; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = 1.0;
}
}
}
for (k=1; k<3; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = 1.0;
Den[n] = 1.0;
Den[N+n] = 0.0;
A_even[n] = 0.3333333333333333;
A_odd[n] = 0.1111111111111111;
A_even[N+n] = 0.1111111111111111;
A_odd[N+n] = 0.1111111111111111;
A_even[2*N+n] = 0.1111111111111111;
A_odd[2*N+n] = 0.1111111111111111;
A_even[3*N+n] = 0.1111111111111111;
B_even[n] = 0.0;
B_odd[n] = 0.0;
B_even[N+n] = 0.0;
B_odd[N+n] = 0.0;
B_even[2*N+n] = 0.0;
B_odd[2*N+n] = 0.0;
B_even[3*N+n] = 0.0;
}
}
}
*/
}
if (pBC && kproc == nprocz-1){
dvc_PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,S,Nx*Ny*(Nz-2));
dvc_ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz,S);
/* // Fill the outlet with component b
for (k=Nz-3; k<Nz-1; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = -1.0;
Den[n] = 0.0;
Den[N+n] = 1.0;
A_even[n] = 0.0;
A_odd[n] = 0.0;
A_even[N+n] = 0.0;
A_odd[N+n] = 0.0;
A_even[2*N+n] = 0.0;
A_odd[2*N+n] = 0.0;
A_even[3*N+n] = 0.0;
B_even[n] = 0.3333333333333333;
B_odd[n] = 0.1111111111111111;
B_even[N+n] = 0.1111111111111111;
B_odd[N+n] = 0.1111111111111111;
B_even[2*N+n] = 0.1111111111111111;
B_odd[2*N+n] = 0.1111111111111111;
B_even[3*N+n] = 0.1111111111111111;
}
}
}
for (k=Nz-1; k<Nz; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = -1.0;
}
}
}
*/
}
//...................................................................................
@ -1870,32 +1733,26 @@ int main(int argc, char **argv)
// Timestep completed!
timestep++;
//...................................................................
if (timestep%10000 == 995){
//...........................................................................
// Copy the phase indicator field for the earlier timestep
dvc_Barrier();
dvc_CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
//...........................................................................
if (timestep%RESTART_INTERVAL == 0){
// Copy the data to the CPU
dvc_CopyToHost(cDistEven,f_even,10*N*sizeof(double));
dvc_CopyToHost(cDistOdd,f_odd,9*N*sizeof(double));
dvc_CopyToHost(cDen,Den,2*N*sizeof(double));
// Read in the restart file to CPU buffers
WriteCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
}
if (timestep%10000 == 0){
//...........................................................................
// Copy the data for for the analysis timestep
//...........................................................................
// Copy the phase from the GPU -> CPU
//...........................................................................
dvc_Barrier();
dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S);
dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
dvc_CopyToHost(Press,Pressure,N*sizeof(double));
dvc_CopyToHost(Vel,Velocity,3*N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
}
if (timestep%10000 == 5){
// End the bubble loop
//...........................................................................
// Copy the phase indicator field for the later timestep
dvc_Barrier();
dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S);
dvc_CopyToHost(Phase_tminus.data,Phi,N*sizeof(double));
dvc_CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
dvc_CopyToHost(Press.data,Pressure,N*sizeof(double));
//...........................................................................
// Calculate the time derivative of the phase indicator field
for (n=0; n<N; n++) dPdt(n) = 0.1*(Phase_tplus(n) - Phase_tminus(n));
@ -1910,6 +1767,29 @@ int main(int argc, char **argv)
//...........................................................................
// Fill in the halo region for the mesh gradients and curvature
//...........................................................................
// Pressure
CommunicateMeshHalo(Press, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
// Mean Curvature
//...........................................................................
CommunicateMeshHalo(MeanCurvature, MPI_COMM_WORLD,
@ -2050,50 +1930,51 @@ int main(int argc, char **argv)
vol_w = vol_n =0.0;
Jwn = efawns = 0.0;
// Phase averages
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
// 1-D index for this cube corner
n = i + j*Nx + k*Nx*Ny;
// Compute the non-wetting phase volume contribution
if ( Phase(i,j,k) > 0 )
nwp_volume += 1.0;
// volume averages over the non-wetting phase
if ( Phase(i,j,k) > 0.99999 ){
// volume the excludes the interfacial region
vol_n += 1.0;
// pressure
pan += Press.data[n];
// velocity
van(0) += Vel[3*n];
van(1) += Vel[3*n+1];
van(2) += Vel[3*n+2];
}
// volume averages over the wetting phase
if ( Phase(i,j,k) < -0.99999 ){
// volume the excludes the interfacial region
vol_w += 1.0;
// pressure
paw += Press.data[n];
// velocity
vaw(0) += Vel[3*n];
vaw(1) += Vel[3*n+1];
vaw(2) += Vel[3*n+2];
}
}
}
}
for (c=0;c<ncubes;c++){
// Get cube from the list
i = cubeList(0,c);
j = cubeList(1,c);
k = cubeList(2,c);
//...........................................................................
// Compute volume averages
for (p=0;p<8;p++){
if ( SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
// 1-D index for this cube corner
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
// Compute the non-wetting phase volume contribution
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 )
nwp_volume += 0.125;
// volume averages over the non-wetting phase
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0.9999 ){
// volume the excludes the interfacial region
vol_n += 0.125;
// pressure
pan += 0.125*Press[n];
// velocity
van(0) += 0.125*Vel[3*n];
van(1) += 0.125*Vel[3*n+1];
van(2) += 0.125*Vel[3*n+2];
}
// volume averages over the wetting phase
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < -0.9999 ){
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
paw += 0.125*Press[n];
// velocity
vaw(0) += 0.125*Vel[3*n];
vaw(1) += 0.125*Vel[3*n+1];
vaw(2) += 0.125*Vel[3*n+2];
}
}
}
//...........................................................................
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(SignDist, Phase, solid_isovalue, fluid_isovalue,
@ -2190,188 +2071,45 @@ int main(int argc, char **argv)
//.........................................................................
if (rank==0){
/* printf("-------------------------------- \n");
printf("Timestep = %i \n", timestep);
printf("NWP volume = %f \n", nwp_volume_global);
printf("Area wn = %f \n", awn_global);
printf("Area ns = %f \n", ans_global);
printf("Area ws = %f \n", aws_global);
printf("Change in surface energy = %f \n", dEs);
printf("-------------------------------- \n");
printf("%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g \n",
timestep,dEs,sat_w,
awn_global,ans_global,aws_global, lwns_global, p_w_global, p_n_global,
vx_w_global, vy_w_global, vz_w_global,
vx_n_global, vy_n_global, vz_n_global);
*/
printf("%.5g ",BubbleRadius); // bubble radius
printf("%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure
printf("%.5g %.5g ",paw_global,pan_global); // saturation and pressure
printf("%.5g ",awn_global); // interfacial area
printf("%.5g ",Jwn_global); // curvature of wn interface
printf("%.5g %.5g %.5g %.5g %.5g %.5g \n",
Gwn_global(0),Gwn_global(1),Gwn_global(2),Gwn_global(3),Gwn_global(4),Gwn_global(5)); // orientation of wn interface
}
}
}
}
//************************************************************************/
dvc_Barrier();
MPI_Barrier(MPI_COMM_WORLD);
stoptime = MPI_Wtime();
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
// Performance obtained from each node
double MLUPS = double(Nx*Ny*Nz)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
//************************************************************************/
// Write out the phase indicator field
//************************************************************************/
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
// printf("Local File Name = %s \n",LocalRankFilename);
// dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
//!!!!!!!!DEBUG HERE!!!!!!!!!!!!!
/* dvc_InitDenColorDistance(ID, Den, Phi, SignDist.data, das, dbs, beta, xIntPos, Nx, Ny, Nz, S);
dvc_InitD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz, S);
dvc_InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz, S);
dvc_ComputeDensityD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz, S);
dvc_ComputeDensityD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz, S);
dvc_ComputePhi(ID, Phi, Den, N, S);
if (pBC && kproc == 0) {
dvc_PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz,S);
// Fill the inlet with component a
for (k=0; k<1; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = 1.0;
}
}
}
for (k=1; k<4; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = 1.0;
Den[n] = 1.0;
Den[N+n] = 0.0;
A_even[n] = 0.3333333333333333;
A_odd[n] = 0.1111111111111111;
A_even[N+n] = 0.1111111111111111;
A_odd[N+n] = 0.1111111111111111;
A_even[2*N+n] = 0.1111111111111111;
A_odd[2*N+n] = 0.1111111111111111;
A_even[3*N+n] = 0.1111111111111111;
B_even[n] = 0.0;
B_odd[n] = 0.0;
B_even[N+n] = 0.0;
B_odd[N+n] = 0.0;
B_even[2*N+n] = 0.0;
B_odd[2*N+n] = 0.0;
B_even[3*N+n] = 0.0;
}
}
}
}
if (pBC && kproc == nprocz-1){
dvc_PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,S,Nx*Ny*(Nz-2));
// Fill the outlet with component b
for (k=Nz-4; k<Nz-1; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = -1.0;
Den[n] = 0.0;
Den[N+n] = 1.0;
A_even[n] = 0.0;
A_odd[n] = 0.0;
A_even[N+n] = 0.0;
A_odd[N+n] = 0.0;
A_even[2*N+n] = 0.0;
A_odd[2*N+n] = 0.0;
A_even[3*N+n] = 0.0;
B_even[n] = 0.3333333333333333;
B_odd[n] = 0.1111111111111111;
B_even[N+n] = 0.1111111111111111;
B_odd[N+n] = 0.1111111111111111;
B_even[2*N+n] = 0.1111111111111111;
B_odd[2*N+n] = 0.1111111111111111;
B_even[3*N+n] = 0.1111111111111111;
}
}
}
for (k=Nz-1; k<Nz; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Phi[n] = -1.0;
}
}
}
}
dvc_SwapD3Q7(ID, A_even, A_odd, Nx, Ny, Nz, S);
dvc_SwapD3Q7(ID, B_even, B_odd, Nx, Ny, Nz, S);
dvc_ComputeDensityD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz, S);
dvc_ComputeDensityD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz, S);
*/
//#ifdef WriteOutput
dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
FILE *PHASE;
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Phase.data,8,N,PHASE);
fwrite(Press.data,8,N,PHASE);
// fwrite(MeanCurvature.data,8,N,PHASE);
fclose(PHASE);
//#endif
dvc_ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz,S);
dvc_CopyToHost(Press,Pressure,N*sizeof(double));
sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
FILE *PRESS;
PRESS = fopen(LocalRankFilename,"wb");
fwrite(Press,8,N,PRESS);
fclose(PRESS);
/* sprintf(LocalRankFilename,"%s%s","dPdt.",LocalRankString);
FILE *SPEED;
SPEED = fopen(LocalRankFilename,"wb");
fwrite(dPdt.data,8,N,SPEED);
fclose(SPEED);
*/
sprintf(LocalRankFilename,"%s%s","DenA.",LocalRankString);
FILE *DENA;
DENA = fopen(LocalRankFilename,"wb");
fwrite(&Den[0],8,N,DENA);
fclose(DENA);
sprintf(LocalRankFilename,"%s%s","DenB.",LocalRankString);
FILE *DENB;
DENB = fopen(LocalRankFilename,"wb");
fwrite(&Den[N],8,N,DENB);
fclose(DENB);
/* sprintf(LocalRankFilename,"%s%s","GradMag.",LocalRankString);
FILE *GRAD;
GRAD = fopen(LocalRankFilename,"wb");
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
Phase(i,j,k) = sqrt(Phase_x(i,j,k)*Phase_x(i,j,k) + Phase_y(i,j,k)*Phase_y(i,j,k) + Phase_z(i,j,k)*Phase_z(i,j,k));
}
}
}
fwrite(Phase.data,8,N,GRAD);
fclose(GRAD);
*/
/* double *DensityValues;
DensityValues = new double [2*N];
dvc_CopyToHost(DensityValues,Copy,2*N*sizeof(double));