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

save the work; to be continued to clean up the code

Browse Source
This commit is contained in:
Rex Zhe Li 2020-08-19 11:31:32 -04:00
parent 2180f7b2bf
commit d8a1837ba0
16 changed files with 373 additions and 17889 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

501
common/ScaLBL.cpp
View File

@ -1094,368 +1094,6 @@ void ScaLBL_Communicator::RecvD3Q19AA(double *dist){
}
void ScaLBL_Communicator::BiSendD3Q19AA(double *Aq, double *Bq){
// NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2
if (Lock==true){
ERROR("ScaLBL Error (SendD3Q19): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?");
}
else{
Lock=true;
}
// assign tag of 19 to D3Q19 communication
sendtag = recvtag = 38;
ScaLBL_DeviceBarrier();
// Pack the distributions
//...Packing for x face(2,8,10,12,14)................................
ScaLBL_D3Q19_Pack(2 ,dvcSendList_x,0*sendCount_x,sendCount_x,sendbuf_x,Aq,N);
ScaLBL_D3Q19_Pack(8 ,dvcSendList_x,1*sendCount_x,sendCount_x,sendbuf_x,Aq,N);
ScaLBL_D3Q19_Pack(10,dvcSendList_x,2*sendCount_x,sendCount_x,sendbuf_x,Aq,N);
ScaLBL_D3Q19_Pack(12,dvcSendList_x,3*sendCount_x,sendCount_x,sendbuf_x,Aq,N);
ScaLBL_D3Q19_Pack(14,dvcSendList_x,4*sendCount_x,sendCount_x,sendbuf_x,Aq,N);
ScaLBL_D3Q19_Pack(2 ,dvcSendList_x,5*sendCount_x,sendCount_x,sendbuf_x,Bq,N);
ScaLBL_D3Q19_Pack(8 ,dvcSendList_x,6*sendCount_x,sendCount_x,sendbuf_x,Bq,N);
ScaLBL_D3Q19_Pack(10,dvcSendList_x,7*sendCount_x,sendCount_x,sendbuf_x,Bq,N);
ScaLBL_D3Q19_Pack(12,dvcSendList_x,8*sendCount_x,sendCount_x,sendbuf_x,Bq,N);
ScaLBL_D3Q19_Pack(14,dvcSendList_x,9*sendCount_x,sendCount_x,sendbuf_x,Bq,N);
req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, 10*sendCount_x,rank_x,sendtag);
req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X, 10*recvCount_X,rank_X,recvtag);
//...Packing for X face(1,7,9,11,13)................................
ScaLBL_D3Q19_Pack(1 ,dvcSendList_X,0*sendCount_X,sendCount_X,sendbuf_X,Aq,N);
ScaLBL_D3Q19_Pack(7 ,dvcSendList_X,1*sendCount_X,sendCount_X,sendbuf_X,Aq,N);
ScaLBL_D3Q19_Pack(9 ,dvcSendList_X,2*sendCount_X,sendCount_X,sendbuf_X,Aq,N);
ScaLBL_D3Q19_Pack(11,dvcSendList_X,3*sendCount_X,sendCount_X,sendbuf_X,Aq,N);
ScaLBL_D3Q19_Pack(13,dvcSendList_X,4*sendCount_X,sendCount_X,sendbuf_X,Aq,N);
ScaLBL_D3Q19_Pack(1 ,dvcSendList_X,5*sendCount_X,sendCount_X,sendbuf_X,Bq,N);
ScaLBL_D3Q19_Pack(7 ,dvcSendList_X,6*sendCount_X,sendCount_X,sendbuf_X,Bq,N);
ScaLBL_D3Q19_Pack(9 ,dvcSendList_X,7*sendCount_X,sendCount_X,sendbuf_X,Bq,N);
ScaLBL_D3Q19_Pack(11,dvcSendList_X,8*sendCount_X,sendCount_X,sendbuf_X,Bq,N);
ScaLBL_D3Q19_Pack(13,dvcSendList_X,9*sendCount_X,sendCount_X,sendbuf_X,Bq,N);
req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, 10*sendCount_X,rank_X,sendtag);
req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, 10*recvCount_x,rank_x,recvtag);
//...Packing for y face(4,8,9,16,18).................................
ScaLBL_D3Q19_Pack(4 ,dvcSendList_y,0*sendCount_y,sendCount_y,sendbuf_y,Aq,N);
ScaLBL_D3Q19_Pack(8 ,dvcSendList_y,1*sendCount_y,sendCount_y,sendbuf_y,Aq,N);
ScaLBL_D3Q19_Pack(9 ,dvcSendList_y,2*sendCount_y,sendCount_y,sendbuf_y,Aq,N);
ScaLBL_D3Q19_Pack(16,dvcSendList_y,3*sendCount_y,sendCount_y,sendbuf_y,Aq,N);
ScaLBL_D3Q19_Pack(18,dvcSendList_y,4*sendCount_y,sendCount_y,sendbuf_y,Aq,N);
ScaLBL_D3Q19_Pack(4 ,dvcSendList_y,5*sendCount_y,sendCount_y,sendbuf_y,Bq,N);
ScaLBL_D3Q19_Pack(8 ,dvcSendList_y,6*sendCount_y,sendCount_y,sendbuf_y,Bq,N);
ScaLBL_D3Q19_Pack(9 ,dvcSendList_y,7*sendCount_y,sendCount_y,sendbuf_y,Bq,N);
ScaLBL_D3Q19_Pack(16,dvcSendList_y,8*sendCount_y,sendCount_y,sendbuf_y,Bq,N);
ScaLBL_D3Q19_Pack(18,dvcSendList_y,9*sendCount_y,sendCount_y,sendbuf_y,Bq,N);
req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y, 10*sendCount_y,rank_y,sendtag);
req2[2] = MPI_COMM_SCALBL.Irecv(recvbuf_Y, 10*recvCount_Y,rank_Y,recvtag);
//...Packing for Y face(3,7,10,15,17).................................
ScaLBL_D3Q19_Pack(3 ,dvcSendList_Y,0*sendCount_Y,sendCount_Y,sendbuf_Y,Aq,N);
ScaLBL_D3Q19_Pack(7 ,dvcSendList_Y,1*sendCount_Y,sendCount_Y,sendbuf_Y,Aq,N);
ScaLBL_D3Q19_Pack(10,dvcSendList_Y,2*sendCount_Y,sendCount_Y,sendbuf_Y,Aq,N);
ScaLBL_D3Q19_Pack(15,dvcSendList_Y,3*sendCount_Y,sendCount_Y,sendbuf_Y,Aq,N);
ScaLBL_D3Q19_Pack(17,dvcSendList_Y,4*sendCount_Y,sendCount_Y,sendbuf_Y,Aq,N);
ScaLBL_D3Q19_Pack(3 ,dvcSendList_Y,5*sendCount_Y,sendCount_Y,sendbuf_Y,Bq,N);
ScaLBL_D3Q19_Pack(7 ,dvcSendList_Y,6*sendCount_Y,sendCount_Y,sendbuf_Y,Bq,N);
ScaLBL_D3Q19_Pack(10,dvcSendList_Y,7*sendCount_Y,sendCount_Y,sendbuf_Y,Bq,N);
ScaLBL_D3Q19_Pack(15,dvcSendList_Y,8*sendCount_Y,sendCount_Y,sendbuf_Y,Bq,N);
ScaLBL_D3Q19_Pack(17,dvcSendList_Y,9*sendCount_Y,sendCount_Y,sendbuf_Y,Bq,N);
req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, 10*sendCount_Y,rank_Y,sendtag);
req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y, 10*recvCount_y,rank_y,recvtag);
//...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q19_Pack(6 ,dvcSendList_z,0*sendCount_z,sendCount_z,sendbuf_z,Aq,N);
ScaLBL_D3Q19_Pack(12,dvcSendList_z,1*sendCount_z,sendCount_z,sendbuf_z,Aq,N);
ScaLBL_D3Q19_Pack(13,dvcSendList_z,2*sendCount_z,sendCount_z,sendbuf_z,Aq,N);
ScaLBL_D3Q19_Pack(16,dvcSendList_z,3*sendCount_z,sendCount_z,sendbuf_z,Aq,N);
ScaLBL_D3Q19_Pack(17,dvcSendList_z,4*sendCount_z,sendCount_z,sendbuf_z,Aq,N);
ScaLBL_D3Q19_Pack(6 ,dvcSendList_z,5*sendCount_z,sendCount_z,sendbuf_z,Bq,N);
ScaLBL_D3Q19_Pack(12,dvcSendList_z,6*sendCount_z,sendCount_z,sendbuf_z,Bq,N);
ScaLBL_D3Q19_Pack(13,dvcSendList_z,7*sendCount_z,sendCount_z,sendbuf_z,Bq,N);
ScaLBL_D3Q19_Pack(16,dvcSendList_z,8*sendCount_z,sendCount_z,sendbuf_z,Bq,N);
ScaLBL_D3Q19_Pack(17,dvcSendList_z,9*sendCount_z,sendCount_z,sendbuf_z,Bq,N);
req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z, 10*sendCount_z,rank_z,sendtag);
req2[4] = MPI_COMM_SCALBL.Irecv(recvbuf_Z, 10*recvCount_Z,rank_Z,recvtag);
//...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q19_Pack(5 ,dvcSendList_Z,0*sendCount_Z,sendCount_Z,sendbuf_Z,Aq,N);
ScaLBL_D3Q19_Pack(11,dvcSendList_Z,1*sendCount_Z,sendCount_Z,sendbuf_Z,Aq,N);
ScaLBL_D3Q19_Pack(14,dvcSendList_Z,2*sendCount_Z,sendCount_Z,sendbuf_Z,Aq,N);
ScaLBL_D3Q19_Pack(15,dvcSendList_Z,3*sendCount_Z,sendCount_Z,sendbuf_Z,Aq,N);
ScaLBL_D3Q19_Pack(18,dvcSendList_Z,4*sendCount_Z,sendCount_Z,sendbuf_Z,Aq,N);
ScaLBL_D3Q19_Pack(5 ,dvcSendList_Z,5*sendCount_Z,sendCount_Z,sendbuf_Z,Bq,N);
ScaLBL_D3Q19_Pack(11,dvcSendList_Z,6*sendCount_Z,sendCount_Z,sendbuf_Z,Bq,N);
ScaLBL_D3Q19_Pack(14,dvcSendList_Z,7*sendCount_Z,sendCount_Z,sendbuf_Z,Bq,N);
ScaLBL_D3Q19_Pack(15,dvcSendList_Z,8*sendCount_Z,sendCount_Z,sendbuf_Z,Bq,N);
ScaLBL_D3Q19_Pack(18,dvcSendList_Z,9*sendCount_Z,sendCount_Z,sendbuf_Z,Bq,N);
req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, 10*sendCount_Z,rank_Z,sendtag);
req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z, 10*recvCount_z,rank_z,recvtag);
//...Pack the xy edge (8)................................
ScaLBL_D3Q19_Pack(8,dvcSendList_xy,0*sendCount_xy,sendCount_xy,sendbuf_xy,Aq,N);
ScaLBL_D3Q19_Pack(8,dvcSendList_xy,1*sendCount_xy,sendCount_xy,sendbuf_xy,Bq,N);
req1[6] = MPI_COMM_SCALBL.Isend(sendbuf_xy, 2*sendCount_xy,rank_xy,sendtag);
req2[6] = MPI_COMM_SCALBL.Irecv(recvbuf_XY, 2*recvCount_XY,rank_XY,recvtag);
//...Pack the Xy edge (9)................................
ScaLBL_D3Q19_Pack(9,dvcSendList_Xy,0*sendCount_Xy,sendCount_Xy,sendbuf_Xy,Aq,N);
ScaLBL_D3Q19_Pack(9,dvcSendList_Xy,1*sendCount_Xy,sendCount_Xy,sendbuf_Xy,Bq,N);
req1[8] = MPI_COMM_SCALBL.Isend(sendbuf_Xy, 2*sendCount_Xy,rank_Xy,sendtag);
req2[8] = MPI_COMM_SCALBL.Irecv(recvbuf_xY, 2*recvCount_xY,rank_xY,recvtag);
//...Pack the xY edge (10)................................
ScaLBL_D3Q19_Pack(10,dvcSendList_xY,0*sendCount_xY,sendCount_xY,sendbuf_xY,Aq,N);
ScaLBL_D3Q19_Pack(10,dvcSendList_xY,1*sendCount_xY,sendCount_xY,sendbuf_xY,Bq,N);
req1[9] = MPI_COMM_SCALBL.Isend(sendbuf_xY, 2*sendCount_xY,rank_xY,sendtag);
req2[9] = MPI_COMM_SCALBL.Irecv(recvbuf_Xy, 2*recvCount_Xy,rank_Xy,recvtag);
//...Pack the XY edge (7)................................
ScaLBL_D3Q19_Pack(7,dvcSendList_XY,0*sendCount_XY,sendCount_XY,sendbuf_XY,Aq,N);
ScaLBL_D3Q19_Pack(7,dvcSendList_XY,1*sendCount_XY,sendCount_XY,sendbuf_XY,Bq,N);
req1[7] = MPI_COMM_SCALBL.Isend(sendbuf_XY, 2*sendCount_XY,rank_XY,sendtag);
req2[7] = MPI_COMM_SCALBL.Irecv(recvbuf_xy, 2*recvCount_xy,rank_xy,recvtag);
//...Pack the xz edge (12)................................
ScaLBL_D3Q19_Pack(12,dvcSendList_xz,0*sendCount_xz,sendCount_xz,sendbuf_xz,Aq,N);
ScaLBL_D3Q19_Pack(12,dvcSendList_xz,1*sendCount_xz,sendCount_xz,sendbuf_xz,Bq,N);
req1[10] = MPI_COMM_SCALBL.Isend(sendbuf_xz, 2*sendCount_xz,rank_xz,sendtag);
req2[10] = MPI_COMM_SCALBL.Irecv(recvbuf_XZ, 2*recvCount_XZ,rank_XZ,recvtag);
//...Pack the xZ edge (14)................................
ScaLBL_D3Q19_Pack(14,dvcSendList_xZ,0*sendCount_xZ,sendCount_xZ,sendbuf_xZ,Aq,N);
ScaLBL_D3Q19_Pack(14,dvcSendList_xZ,1*sendCount_xZ,sendCount_xZ,sendbuf_xZ,Bq,N);
req1[13] = MPI_COMM_SCALBL.Isend(sendbuf_xZ, 2*sendCount_xZ,rank_xZ,sendtag);
req2[13] = MPI_COMM_SCALBL.Irecv(recvbuf_Xz, 2*recvCount_Xz,rank_Xz,recvtag);
//...Pack the Xz edge (13)................................
ScaLBL_D3Q19_Pack(13,dvcSendList_Xz,0*sendCount_Xz,sendCount_Xz,sendbuf_Xz,Aq,N);
ScaLBL_D3Q19_Pack(13,dvcSendList_Xz,1*sendCount_Xz,sendCount_Xz,sendbuf_Xz,Bq,N);
req1[12] = MPI_COMM_SCALBL.Isend(sendbuf_Xz, 2*sendCount_Xz,rank_Xz,sendtag);
req2[12] = MPI_COMM_SCALBL.Irecv(recvbuf_xZ, 2*recvCount_xZ,rank_xZ,recvtag);
//...Pack the XZ edge (11)................................
ScaLBL_D3Q19_Pack(11,dvcSendList_XZ,0*sendCount_XZ,sendCount_XZ,sendbuf_XZ,Aq,N);
ScaLBL_D3Q19_Pack(11,dvcSendList_XZ,1*sendCount_XZ,sendCount_XZ,sendbuf_XZ,Bq,N);
req1[11] = MPI_COMM_SCALBL.Isend(sendbuf_XZ, 2*sendCount_XZ,rank_XZ,sendtag);
req2[11] = MPI_COMM_SCALBL.Irecv(recvbuf_xz, 2*recvCount_xz,rank_xz,recvtag);
//...Pack the yz edge (16)................................
ScaLBL_D3Q19_Pack(16,dvcSendList_yz,0*sendCount_yz,sendCount_yz,sendbuf_yz,Aq,N);
ScaLBL_D3Q19_Pack(16,dvcSendList_yz,1*sendCount_yz,sendCount_yz,sendbuf_yz,Bq,N);
req1[14] = MPI_COMM_SCALBL.Isend(sendbuf_yz, 2*sendCount_yz,rank_yz,sendtag);
req2[14] = MPI_COMM_SCALBL.Irecv(recvbuf_YZ, 2*recvCount_YZ,rank_YZ,recvtag);
//...Pack the yZ edge (18)................................
ScaLBL_D3Q19_Pack(18,dvcSendList_yZ,0*sendCount_yZ,sendCount_yZ,sendbuf_yZ,Aq,N);
ScaLBL_D3Q19_Pack(18,dvcSendList_yZ,1*sendCount_yZ,sendCount_yZ,sendbuf_yZ,Bq,N);
req1[17] = MPI_COMM_SCALBL.Isend(sendbuf_yZ, 2*sendCount_yZ,rank_yZ,sendtag);
req2[17] = MPI_COMM_SCALBL.Irecv(recvbuf_Yz, 2*recvCount_Yz,rank_Yz,recvtag);
//...Pack the Yz edge (17)................................
ScaLBL_D3Q19_Pack(17,dvcSendList_Yz,0*sendCount_Yz,sendCount_Yz,sendbuf_Yz,Aq,N);
ScaLBL_D3Q19_Pack(17,dvcSendList_Yz,1*sendCount_Yz,sendCount_Yz,sendbuf_Yz,Bq,N);
req1[16] = MPI_COMM_SCALBL.Isend(sendbuf_Yz, 2*sendCount_Yz,rank_Yz,sendtag);
req2[16] = MPI_COMM_SCALBL.Irecv(recvbuf_yZ, 2*recvCount_yZ,rank_yZ,recvtag);
//...Pack the YZ edge (15)................................
ScaLBL_D3Q19_Pack(15,dvcSendList_YZ,0*sendCount_YZ,sendCount_YZ,sendbuf_YZ,Aq,N);
ScaLBL_D3Q19_Pack(15,dvcSendList_YZ,1*sendCount_YZ,sendCount_YZ,sendbuf_YZ,Bq,N);
req1[15] = MPI_COMM_SCALBL.Isend(sendbuf_YZ, 2*sendCount_YZ,rank_YZ,sendtag);
req2[15] = MPI_COMM_SCALBL.Irecv(recvbuf_yz, 2*recvCount_yz,rank_yz,recvtag);
//...................................................................................
}
void ScaLBL_Communicator::BiRecvD3Q19AA(double *Aq, double *Bq){
// NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2
//...................................................................................
// Wait for completion of D3Q19 communication
MPI_COMM_SCALBL.waitAll(18,req1);
MPI_COMM_SCALBL.waitAll(18,req2);
ScaLBL_DeviceBarrier();
//...................................................................................
// NOTE: AA Routine writes to opposite
// Unpack the distributions on the device
//...................................................................................
//...Unpacking for x face(2,8,10,12,14)................................
ScaLBL_D3Q19_Unpack(2, dvcRecvDist_x,0*recvCount_x,recvCount_x,recvbuf_x,Aq,N);
ScaLBL_D3Q19_Unpack(8, dvcRecvDist_x,1*recvCount_x,recvCount_x,recvbuf_x,Aq,N);
ScaLBL_D3Q19_Unpack(10,dvcRecvDist_x,2*recvCount_x,recvCount_x,recvbuf_x,Aq,N);
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_x,3*recvCount_x,recvCount_x,recvbuf_x,Aq,N);
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_x,4*recvCount_x,recvCount_x,recvbuf_x,Aq,N);
ScaLBL_D3Q19_Unpack(2, dvcRecvDist_x,0*recvCount_x,recvCount_x,&recvbuf_x[5*recvCount_x],Bq,N);
ScaLBL_D3Q19_Unpack(8, dvcRecvDist_x,1*recvCount_x,recvCount_x,&recvbuf_x[5*recvCount_x],Bq,N);
ScaLBL_D3Q19_Unpack(10,dvcRecvDist_x,2*recvCount_x,recvCount_x,&recvbuf_x[5*recvCount_x],Bq,N);
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_x,3*recvCount_x,recvCount_x,&recvbuf_x[5*recvCount_x],Bq,N);
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_x,4*recvCount_x,recvCount_x,&recvbuf_x[5*recvCount_x],Bq,N);
//...................................................................................
//...Packing for X face(1,7,9,11,13)................................
ScaLBL_D3Q19_Unpack(1, dvcRecvDist_X,0*recvCount_X,recvCount_X,recvbuf_X,Aq,N);
ScaLBL_D3Q19_Unpack(7, dvcRecvDist_X,1*recvCount_X,recvCount_X,recvbuf_X,Aq,N);
ScaLBL_D3Q19_Unpack(9, dvcRecvDist_X,2*recvCount_X,recvCount_X,recvbuf_X,Aq,N);
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_X,3*recvCount_X,recvCount_X,recvbuf_X,Aq,N);
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_X,4*recvCount_X,recvCount_X,recvbuf_X,Aq,N);
ScaLBL_D3Q19_Unpack(1, dvcRecvDist_X,0*recvCount_X,recvCount_X,&recvbuf_X[5*recvCount_X],Bq,N);
ScaLBL_D3Q19_Unpack(7, dvcRecvDist_X,1*recvCount_X,recvCount_X,&recvbuf_X[5*recvCount_X],Bq,N);
ScaLBL_D3Q19_Unpack(9, dvcRecvDist_X,2*recvCount_X,recvCount_X,&recvbuf_X[5*recvCount_X],Bq,N);
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_X,3*recvCount_X,recvCount_X,&recvbuf_X[5*recvCount_X],Bq,N);
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_X,4*recvCount_X,recvCount_X,&recvbuf_X[5*recvCount_X],Bq,N);
//...................................................................................
//...Packing for y face(4,8,9,16,18).................................
ScaLBL_D3Q19_Unpack(4, dvcRecvDist_y,0*recvCount_y,recvCount_y,recvbuf_y,Aq,N);
ScaLBL_D3Q19_Unpack(8, dvcRecvDist_y,1*recvCount_y,recvCount_y,recvbuf_y,Aq,N);
ScaLBL_D3Q19_Unpack(9, dvcRecvDist_y,2*recvCount_y,recvCount_y,recvbuf_y,Aq,N);
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_y,3*recvCount_y,recvCount_y,recvbuf_y,Aq,N);
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_y,4*recvCount_y,recvCount_y,recvbuf_y,Aq,N);
ScaLBL_D3Q19_Unpack(4, dvcRecvDist_y,0*recvCount_y,recvCount_y,&recvbuf_y[5*recvCount_y],Bq,N);
ScaLBL_D3Q19_Unpack(8, dvcRecvDist_y,1*recvCount_y,recvCount_y,&recvbuf_y[5*recvCount_y],Bq,N);
ScaLBL_D3Q19_Unpack(9, dvcRecvDist_y,2*recvCount_y,recvCount_y,&recvbuf_y[5*recvCount_y],Bq,N);
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_y,3*recvCount_y,recvCount_y,&recvbuf_y[5*recvCount_y],Bq,N);
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_y,4*recvCount_y,recvCount_y,&recvbuf_y[5*recvCount_y],Bq,N);
//...................................................................................
//...Packing for Y face(3,7,10,15,17).................................
ScaLBL_D3Q19_Unpack(3, dvcRecvDist_Y,0*recvCount_Y,recvCount_Y,recvbuf_Y,Aq,N);
ScaLBL_D3Q19_Unpack(7, dvcRecvDist_Y,1*recvCount_Y,recvCount_Y,recvbuf_Y,Aq,N);
ScaLBL_D3Q19_Unpack(10,dvcRecvDist_Y,2*recvCount_Y,recvCount_Y,recvbuf_Y,Aq,N);
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Y,3*recvCount_Y,recvCount_Y,recvbuf_Y,Aq,N);
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Y,4*recvCount_Y,recvCount_Y,recvbuf_Y,Aq,N);
ScaLBL_D3Q19_Unpack(3, dvcRecvDist_Y,0*recvCount_Y,recvCount_Y,&recvbuf_Y[5*recvCount_Y],Bq,N);
ScaLBL_D3Q19_Unpack(7, dvcRecvDist_Y,1*recvCount_Y,recvCount_Y,&recvbuf_Y[5*recvCount_Y],Bq,N);
ScaLBL_D3Q19_Unpack(10,dvcRecvDist_Y,2*recvCount_Y,recvCount_Y,&recvbuf_Y[5*recvCount_Y],Bq,N);
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Y,3*recvCount_Y,recvCount_Y,&recvbuf_Y[5*recvCount_Y],Bq,N);
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Y,4*recvCount_Y,recvCount_Y,&recvbuf_Y[5*recvCount_Y],Bq,N);
//...................................................................................
//...Pack the xy edge (8)................................
ScaLBL_D3Q19_Unpack(8,dvcRecvDist_xy,0,recvCount_xy,recvbuf_xy,Aq,N);
ScaLBL_D3Q19_Unpack(8,dvcRecvDist_xy,0,recvCount_xy,&recvbuf_xy[recvCount_xy],Bq,N);
//...Pack the Xy edge (9)................................
ScaLBL_D3Q19_Unpack(9,dvcRecvDist_Xy,0,recvCount_Xy,recvbuf_Xy,Aq,N);
ScaLBL_D3Q19_Unpack(9,dvcRecvDist_Xy,0,recvCount_Xy,&recvbuf_Xy[recvCount_Xy],Bq,N);
//...Pack the xY edge (10)................................
ScaLBL_D3Q19_Unpack(10,dvcRecvDist_xY,0,recvCount_xY,recvbuf_xY,Aq,N);
ScaLBL_D3Q19_Unpack(10,dvcRecvDist_xY,0,recvCount_xY,&recvbuf_xY[recvCount_xY],Bq,N);
//...Pack the XY edge (7)................................
ScaLBL_D3Q19_Unpack(7,dvcRecvDist_XY,0,recvCount_XY,recvbuf_XY,Aq,N);
ScaLBL_D3Q19_Unpack(7,dvcRecvDist_XY,0,recvCount_XY,&recvbuf_XY[recvCount_XY],Bq,N);
if (BoundaryCondition > 0 && kproc == 0){
// don't unpack little z
//...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q19_Unpack(5, dvcRecvDist_Z,0*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_Z,1*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_Z,4*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(5, dvcRecvDist_Z,0*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_Z,1*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_Z,2*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Z,3*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_Z,4*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
//...Pack the xZ edge (14)................................
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_xZ,0,recvCount_xZ,recvbuf_xZ,Aq,N);
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_xZ,0,recvCount_xZ,&recvbuf_xZ[recvCount_xZ],Bq,N);
//...Pack the XZ edge (11)................................
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_XZ,0,recvCount_XZ,recvbuf_XZ,Aq,N);
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_XZ,0,recvCount_XZ,&recvbuf_XZ[recvCount_XZ],Bq,N);
//...Pack the yZ edge (18)................................
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_yZ,0,recvCount_yZ,recvbuf_yZ,Aq,N);
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_yZ,0,recvCount_yZ,&recvbuf_yZ[recvCount_yZ],Bq,N);
//...Pack the YZ edge (15)................................
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_YZ,0,recvCount_YZ,recvbuf_YZ,Aq,N);
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_YZ,0,recvCount_YZ,&recvbuf_YZ[recvCount_YZ],Bq,N);
}
else if (BoundaryCondition > 0 && kproc == nprocz-1){
// don't unpack big Z
//...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q19_Unpack(6, dvcRecvDist_z,0*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_z,1*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_z,2*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_z,3*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_z,4*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(6, dvcRecvDist_z,0*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_z,1*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_z,2*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_z,3*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_z,4*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
//...Pack the xz edge (12)................................
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_xz,0,recvCount_xz,recvbuf_xz,Aq,N);
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_xz,0,recvCount_xz,&recvbuf_xz[recvCount_xz],Bq,N);
//...Pack the Xz edge (13)................................
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_Xz,0,recvCount_Xz,recvbuf_Xz,Aq,N);
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_Xz,0,recvCount_Xz,&recvbuf_Xz[recvCount_Xz],Bq,N);
//...Pack the yz edge (16)................................
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_yz,0,recvCount_yz,recvbuf_yz,Aq,N);
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_yz,0,recvCount_yz,&recvbuf_yz[recvCount_yz],Bq,N);
//...Pack the Yz edge (17)................................
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Yz,0,recvCount_Yz,recvbuf_Yz,Aq,N);
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Yz,0,recvCount_Yz,&recvbuf_Yz[recvCount_Yz],Bq,N);
}
else {
//...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q19_Unpack(6, dvcRecvDist_z,0*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_z,1*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_z,2*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_z,3*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_z,4*recvCount_z,recvCount_z,recvbuf_z,Aq,N);
ScaLBL_D3Q19_Unpack(6, dvcRecvDist_z,0*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_z,1*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_z,2*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_z,3*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_z,4*recvCount_z,recvCount_z,&recvbuf_z[5*recvCount_z],Bq,N);
//...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q19_Unpack(5, dvcRecvDist_Z,0*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_Z,1*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_Z,4*recvCount_Z,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q19_Unpack(5, dvcRecvDist_Z,0*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_Z,1*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_Z,2*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Z,3*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_Z,4*recvCount_Z,recvCount_Z,&recvbuf_Z[5*recvCount_Z],Bq,N);
//...Pack the xZ edge (14)................................
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_xZ,0,recvCount_xZ,recvbuf_xZ,Aq,N);
ScaLBL_D3Q19_Unpack(14,dvcRecvDist_xZ,0,recvCount_xZ,&recvbuf_xZ[recvCount_xZ],Bq,N);
//...Pack the XZ edge (11)................................
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_XZ,0,recvCount_XZ,recvbuf_XZ,Aq,N);
ScaLBL_D3Q19_Unpack(11,dvcRecvDist_XZ,0,recvCount_XZ,&recvbuf_XZ[recvCount_XZ],Bq,N);
//...Pack the yZ edge (18)................................
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_yZ,0,recvCount_yZ,recvbuf_yZ,Aq,N);
ScaLBL_D3Q19_Unpack(18,dvcRecvDist_yZ,0,recvCount_yZ,&recvbuf_yZ[recvCount_yZ],Bq,N);
//...Pack the YZ edge (15)................................
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_YZ,0,recvCount_YZ,recvbuf_YZ,Aq,N);
ScaLBL_D3Q19_Unpack(15,dvcRecvDist_YZ,0,recvCount_YZ,&recvbuf_YZ[recvCount_YZ],Bq,N);
//...Pack the xz edge (12)................................
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_xz,0,recvCount_xz,recvbuf_xz,Aq,N);
ScaLBL_D3Q19_Unpack(12,dvcRecvDist_xz,0,recvCount_xz,&recvbuf_xz[recvCount_xz],Bq,N);
//...Pack the Xz edge (13)................................
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_Xz,0,recvCount_Xz,recvbuf_Xz,Aq,N);
ScaLBL_D3Q19_Unpack(13,dvcRecvDist_Xz,0,recvCount_Xz,&recvbuf_Xz[recvCount_Xz],Bq,N);
//...Pack the yz edge (16)................................
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_yz,0,recvCount_yz,recvbuf_yz,Aq,N);
ScaLBL_D3Q19_Unpack(16,dvcRecvDist_yz,0,recvCount_yz,&recvbuf_yz[recvCount_yz],Bq,N);
//...Pack the Yz edge (17)................................
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Yz,0,recvCount_Yz,recvbuf_Yz,Aq,N);
ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Yz,0,recvCount_Yz,&recvbuf_Yz[recvCount_Yz],Bq,N);
}
//...................................................................................
Lock=false; // unlock the communicator after communications complete
//...................................................................................
}
void ScaLBL_Communicator::RecvGrad(double *phi, double *grad){
// Recieves halo and incorporates into D3Q19 based stencil gradient computation
@ -1539,106 +1177,6 @@ void ScaLBL_Communicator::RecvGrad(double *phi, double *grad){
}
void ScaLBL_Communicator::SendD3Q7AA(double *Aq){
// NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2
if (Lock==true){
ERROR("ScaLBL Error (SendD3Q7): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?");
}
else{
Lock=true;
}
// assign tag of 19 to D3Q19 communication
sendtag = recvtag = 7;
ScaLBL_DeviceBarrier();
// Pack the distributions
//...Packing for x face(2,8,10,12,14)................................
ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,Aq,N);
req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x,rank_x,sendtag);
req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X, recvCount_X,rank_X,recvtag);
//...Packing for X face(1,7,9,11,13)................................
ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,Aq,N);
req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, sendCount_X,rank_X,sendtag);
req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, recvCount_x,rank_x,recvtag);
//...Packing for y face(4,8,9,16,18).................................
ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,Aq,N);
req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y, sendCount_y,rank_y,sendtag);
req2[2] = MPI_COMM_SCALBL.Irecv(recvbuf_Y, recvCount_Y,rank_Y,recvtag);
//...Packing for Y face(3,7,10,15,17).................................
ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,Aq,N);
req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, sendCount_Y,rank_Y,sendtag);
req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y, recvCount_y,rank_y,recvtag);
//...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q19_Pack(6,dvcSendList_z,0,sendCount_z,sendbuf_z,Aq,N);
req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z, sendCount_z,rank_z,sendtag);
req2[4] = MPI_COMM_SCALBL.Irecv(recvbuf_Z, recvCount_Z,rank_Z,recvtag);
//...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q19_Pack(5,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,Aq,N);
req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, sendCount_Z,rank_Z,sendtag);
req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z, recvCount_z,rank_z,recvtag);
//...................................................................................
}
void ScaLBL_Communicator::RecvD3Q7AA(double *Aq){
// NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2
//...................................................................................
// Wait for completion of D3Q19 communication
MPI_COMM_SCALBL.waitAll(6,req1);
MPI_COMM_SCALBL.waitAll(6,req2);
ScaLBL_DeviceBarrier();
//...................................................................................
// NOTE: AA Routine writes to opposite
// Unpack the distributions on the device
//...................................................................................
//...Unpacking for x face(2,8,10,12,14)................................
ScaLBL_D3Q7_Unpack(2,dvcRecvDist_x,0,recvCount_x,recvbuf_x,Aq,N);
//...................................................................................
//...Packing for X face(1,7,9,11,13)................................
ScaLBL_D3Q7_Unpack(1,dvcRecvDist_X,0,recvCount_X,recvbuf_X,Aq,N);
//...................................................................................
//...Packing for y face(4,8,9,16,18).................................
ScaLBL_D3Q7_Unpack(4,dvcRecvDist_y,0,recvCount_y,recvbuf_y,Aq,N);
//...................................................................................
//...Packing for Y face(3,7,10,15,17).................................
ScaLBL_D3Q7_Unpack(3,dvcRecvDist_Y,0,recvCount_Y,recvbuf_Y,Aq,N);
//...................................................................................
if (BoundaryCondition > 0 && kproc == 0){
// don't unpack little z
//...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,Aq,N);
}
else if (BoundaryCondition > 0 && kproc == nprocz-1){
// don't unpack big z
//...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,Aq,N);
}
else {
//...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,Aq,N);
//...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,Aq,N);
}
//...................................................................................
Lock=false; // unlock the communicator after communications complete
//...................................................................................
}
void ScaLBL_Communicator::BiSendD3Q7AA(double *Aq, double *Bq){
// NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2
@ -2171,42 +1709,3 @@ void ScaLBL_Communicator::PrintD3Q19(){
delete [] TempBuffer;
}
void ScaLBL_Communicator::GreyscaleSC_BC_z(int *Map, double *DenA, double *DenB, double vA, double vB)
{
if (kproc == 0) {
// Set the density field on the z inlet
ScaLBL_GreyscaleSC_BC_z(dvcSendList_z, Map, DenA, DenB, vA, vB, sendCount_z);
}
}
void ScaLBL_Communicator::GreyscaleSC_BC_Z(int *Map, double *DenA, double *DenB, double vA, double vB)
{
if (kproc == nprocz-1){
// Set the density field on the Z outlet
ScaLBL_GreyscaleSC_BC_Z(dvcSendList_Z, Map, DenA, DenB, vA, vB, sendCount_Z);
}
}
void ScaLBL_Communicator::GreyscaleSC_Pressure_BC_z(int *neighborList, double *fqA, double *fqB, double dinA, double dinB, int time)
{
if (kproc == 0) {
if (time%2==0){
ScaLBL_GreyscaleSC_AAeven_Pressure_BC_z(dvcSendList_z, fqA, fqB, dinA, dinB, sendCount_z, N);
}
else{
ScaLBL_GreyscaleSC_AAodd_Pressure_BC_z(neighborList, dvcSendList_z, fqA, fqB, dinA, dinB, sendCount_z, N);
}
}
}
void ScaLBL_Communicator::GreyscaleSC_Pressure_BC_Z(int *neighborList, double *fqA, double *fqB, double doutA, double doutB, int time)
{
if (kproc == nprocz-1){
if (time%2==0){
ScaLBL_GreyscaleSC_AAeven_Pressure_BC_Z(dvcSendList_Z, fqA, fqB, doutA, doutB, sendCount_Z, N);
}
else{
ScaLBL_GreyscaleSC_AAodd_Pressure_BC_Z(neighborList, dvcSendList_Z, fqA, fqB, doutA, doutB, sendCount_Z, N);
}
}
}

164
common/ScaLBL.h
View File

@ -79,89 +79,89 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
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);
//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);
//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);
@ -298,12 +298,8 @@ public:
int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np);
void SendD3Q19AA(double *dist);
void RecvD3Q19AA(double *dist);
void BiSendD3Q19AA(double *Aq, double *Bq);
void BiRecvD3Q19AA(double *Aq, double *Bq);
// 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 *Aq);
void RecvD3Q7AA(double *Aq);
void BiSendD3Q7AA(double *Aq, double *Bq);
void BiRecvD3Q7AA(double *Aq, double *Bq);
void TriSendD3Q7AA(double *Aq, double *Bq, double *Cq);

3113
cpu/GreyscaleFE.cpp
View File

File diff suppressed because it is too large Load Diff

3625
cpu/GreyscaleSC.cpp
View File

File diff suppressed because it is too large Load Diff

3338
gpu/GreyscaleFE.cu
View File

File diff suppressed because it is too large Load Diff

3819
gpu/GreyscaleSC.cu
View File

File diff suppressed because it is too large Load Diff

495
models/GreyscaleColorModel.cpp
View File

@ -1,5 +1,5 @@
/*
color lattice boltzmann model
Two-fluid greyscale color lattice boltzmann model
*/
#include "models/GreyscaleColorModel.h"
#include "analysis/distance.h"
@ -9,7 +9,7 @@ color lattice boltzmann model
#include <stdlib.h>
#include <time.h>
ScaLBL_GreyscaleColorModel::ScaLBL_GreyscaleColorModel(int RANK, int NP, const Utilities::MPI& COMM):
ScaLBL_GreyscaleColorModel::ScaLBL_GreyscaleColorModel(int RANK, int NP, MPI_Comm COMM):
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),tauA_eff(0),tauB_eff(0),rhoA(0),rhoB(0),alpha(0),beta(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0),GreyPorosity(0),
Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),greyMode(0),Lx(0),Ly(0),Lz(0),comm(COMM)
@ -56,8 +56,6 @@ void ScaLBL_GreyscaleColorModel::ReadCheckpoint(char *FILENAME, double *cPhi, do
File.close();
}
*/
void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
// read the input database
db = std::make_shared<Database>( filename );
@ -139,21 +137,21 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
// Override user-specified boundary condition for specific protocols
auto protocol = greyscaleColor_db->getWithDefault<std::string>( "protocol", "none" );
if (protocol == "seed water"){
if (BoundaryCondition != 0 ){
if (BoundaryCondition != 0 && BoundaryCondition != 5){
BoundaryCondition = 0;
if (rank==0) printf("WARNING: protocol (seed water) supports only full periodic boundary condition \n");
}
domain_db->putScalar<int>( "BC", BoundaryCondition );
}
else if (protocol == "open connected oil"){
if (BoundaryCondition != 0 ){
if (BoundaryCondition != 0 && BoundaryCondition != 5){
BoundaryCondition = 0;
if (rank==0) printf("WARNING: protocol (open connected oil) supports only full periodic boundary condition \n");
}
domain_db->putScalar<int>( "BC", BoundaryCondition );
}
else if (protocol == "shell aggregation"){
if (BoundaryCondition != 0 ){
if (BoundaryCondition != 0 && BoundaryCondition != 5){
BoundaryCondition = 0;
if (rank==0) printf("WARNING: protocol (shell aggregation) supports only full periodic boundary condition \n");
}
@ -176,9 +174,9 @@ void ScaLBL_GreyscaleColorModel::SetDomain(){
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
Averages = std::shared_ptr<SubPhase> ( new SubPhase(Dm) ); // TwoPhase analysis object
comm.barrier();
MPI_Barrier(comm);
Dm->CommInit();
comm.barrier();
MPI_Barrier(comm);
// Read domain parameters
rank = Dm->rank();
nprocx = Dm->nprocx();
@ -201,12 +199,12 @@ void ScaLBL_GreyscaleColorModel::ReadInput(){
}
else if (domain_db->keyExists( "GridFile" )){
// Read the local domain data
auto input_id = readMicroCT( *domain_db, comm );
auto input_id = readMicroCT( *domain_db, MPI_COMM_WORLD );
// Fill the halo (assuming GCW of 1)
array<int,3> 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<signed char> fill( comm, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
fillHalo<signed char> fill( MPI_COMM_WORLD, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
Array<signed char> id_view;
id_view.viewRaw( size1, Mask->id );
fill.copy( input_id, id_view );
@ -315,7 +313,7 @@ void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
for (size_t idx=0; idx<NLABELS; idx++)
label_count_global[idx] = Dm->Comm.sumReduce( label_count[idx] );
label_count_global[idx] = sumReduce( Dm->Comm, label_count[idx]);
if (rank==0){
printf("Number of component labels: %lu \n",NLABELS);
@ -329,6 +327,7 @@ void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
ScaLBL_CopyToDevice(Phi, phase, N*sizeof(double));
ScaLBL_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
delete [] phase;
}
@ -586,7 +585,7 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
// Set Dm to match Mask
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
for (int idx=0; idx<NLABELS; idx++) label_count_global[idx]=Dm->Comm.sumReduce(label_count[idx]);
for (int idx=0; idx<NLABELS; idx++) label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
//Initialize a weighted porosity after considering grey voxels
GreyPorosity=0.0;
@ -597,7 +596,7 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
if (rank==0){
printf("Image resolution: %.5g [um/voxel]\n",Dm->voxel_length);
printf("Component labels: %lu \n",NLABELS);
printf("Number of component labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
POROSITY=PorosityList[idx];
@ -617,7 +616,6 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
delete [] Permeability;
}
void ScaLBL_GreyscaleColorModel::Create(){
/*
* This function creates the variables needed to run a LBM
@ -644,7 +642,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
Map.resize(Nx,Ny,Nz); Map.fill(-2);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
comm.barrier();
MPI_Barrier(comm);
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
@ -716,7 +714,6 @@ void ScaLBL_GreyscaleColorModel::Create(){
}
}
void ScaLBL_GreyscaleColorModel::Initialize(){
if (rank==0) printf ("Initializing distributions \n");
@ -782,14 +779,15 @@ void ScaLBL_GreyscaleColorModel::Initialize(){
ScaLBL_CopyToDevice(Phi,cPhi,N*sizeof(double));
ScaLBL_DeviceBarrier();
comm.barrier();
MPI_Barrier(comm);
}
if (rank==0) printf ("Initializing phase field \n");
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition >0 ){
// 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);
@ -812,6 +810,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
int IMAGE_COUNT = 0;
std::vector<std::string> ImageList;
bool SET_CAPILLARY_NUMBER = false;
bool RESCALE_FORCE = false;
bool MORPH_ADAPT = false;
bool USE_MORPH = false;
bool USE_SEED = false;
@ -820,6 +819,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
int MAX_MORPH_TIMESTEPS = 50000; // maximum number of LBM timesteps to spend in morphological adaptation routine
int MIN_STEADY_TIMESTEPS = 100000;
int MAX_STEADY_TIMESTEPS = 200000;
int RESCALE_FORCE_AFTER_TIMESTEP = 0;
int RAMP_TIMESTEPS = 0;//50000; // number of timesteps to run initially (to get a reasonable velocity field before other pieces kick in)
int CURRENT_MORPH_TIMESTEPS=0; // counter for number of timesteps spent in morphological adaptation routine (reset each time)
int CURRENT_STEADY_TIMESTEPS=0; // counter for number of timesteps spent in morphological adaptation routine (reset each time)
@ -834,17 +834,13 @@ void ScaLBL_GreyscaleColorModel::Run(){
double initial_volume = 0.0;
double delta_volume = 0.0;
double delta_volume_target = 0.0;
double RESIDUAL_ENDPOINT_THRESHOLD = 0.04;
double NOISE_THRESHOLD = 0.0;
double BUMP_RATE = 2.0;
bool USE_BUMP_RATE = false;
/* history for morphological algoirthm */
double KRA_MORPH_FACTOR=0.8;
double KRA_MORPH_FACTOR=0.5;
double volA_prev = 0.0;
double log_krA_prev = 1.0;
double log_krA_target = 1.0;
double log_krA = 0.0;
double log_krA = 1.0;
double slope_krA_volume = 0.0;
if (greyscaleColor_db->keyExists( "vol_A_previous" )){
volA_prev = greyscaleColor_db->getScalar<double>( "vol_A_previous" );
@ -868,46 +864,33 @@ void ScaLBL_GreyscaleColorModel::Run(){
USE_MORPH = true;
}
else if (protocol == "seed water"){
morph_delta = 0.05;
morph_delta = -0.05;
seed_water = 0.01;
USE_SEED = true;
USE_MORPH = true;
}
else if (protocol == "open connected oil"){
morph_delta = 0.05;
morph_delta = -0.05;
USE_MORPH = true;
USE_MORPHOPEN_OIL = true;
}
else if (protocol == "shell aggregation"){
morph_delta = 0.05;
morph_delta = -0.05;
USE_MORPH = true;
}
if (greyscaleColor_db->keyExists( "residual_endpoint_threshold" )){
RESIDUAL_ENDPOINT_THRESHOLD = greyscaleColor_db->getScalar<double>( "residual_endpoint_threshold" );
}
NULL_USE( RESIDUAL_ENDPOINT_THRESHOLD );
if (greyscaleColor_db->keyExists( "noise_threshold" )){
NOISE_THRESHOLD = greyscaleColor_db->getScalar<double>( "noise_threshold" );
USE_BUMP_RATE = true;
}
if (greyscaleColor_db->keyExists( "bump_rate" )){
BUMP_RATE = greyscaleColor_db->getScalar<double>( "bump_rate" );
USE_BUMP_RATE = true;
}
if (greyscaleColor_db->keyExists( "capillary_number" )){
capillary_number = greyscaleColor_db->getScalar<double>( "capillary_number" );
SET_CAPILLARY_NUMBER=true;
//RESCALE_FORCE_MAX = 1;
}
// if (analysis_db->keyExists( "rescale_force_count" )){
// RESCALE_FORCE_MAX = analysis_db->getScalar<int>( "rescale_force_count" );
// }
if (greyscaleColor_db->keyExists( "rescale_force_after_timestep" )){
RESCALE_FORCE_AFTER_TIMESTEP = greyscaleColor_db->getScalar<int>( "rescale_force_after_timestep" );
RESCALE_FORCE = true;
}
if (greyscaleColor_db->keyExists( "timestep" )){
timestep = greyscaleColor_db->getScalar<int>( "timestep" );
}
if (BoundaryCondition != 0 && SET_CAPILLARY_NUMBER==true){
if (rank == 0) printf("WARINING: capillary number target only supported for BC = 0 \n");
if (BoundaryCondition != 0 && BoundaryCondition != 5 && SET_CAPILLARY_NUMBER==true){
if (rank == 0) printf("WARINING: capillary number target only supported for BC = 0 or 5 \n");
SET_CAPILLARY_NUMBER=false;
}
if (analysis_db->keyExists( "seed_water" )){
@ -984,8 +967,8 @@ void ScaLBL_GreyscaleColorModel::Run(){
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier();
comm.barrier();
starttime = Utilities::MPI::time();
MPI_Barrier(comm);
starttime = MPI_Wtime();
//.........................................
//************ MAIN ITERATION LOOP ***************************************/
@ -1010,7 +993,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
if (BoundaryCondition > 0){
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);
}
@ -1042,6 +1025,10 @@ void ScaLBL_GreyscaleColorModel::Run(){
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);
}
if (greyMode==true){
//Model-1&4
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,
@ -1057,7 +1044,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
}
ScaLBL_DeviceBarrier();
comm.barrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
// *************EVEN TIMESTEP*************
timestep++;
@ -1071,7 +1058,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
// Halo exchange for phase field
if (BoundaryCondition > 0){
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);
}
@ -1102,6 +1089,10 @@ void ScaLBL_GreyscaleColorModel::Run(){
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);
}
if (greyMode==true){
//Model-1&4
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,
@ -1117,18 +1108,16 @@ void ScaLBL_GreyscaleColorModel::Run(){
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
}
ScaLBL_DeviceBarrier();
comm.barrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
//************************************************************************
PROFILE_STOP("Update");
if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition > 0){
if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition == 4){
printf("%i %f \n",timestep,din);
}
// Run the analysis
analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
// allow initial ramp-up to get closer to steady state
if (timestep > RAMP_TIMESTEPS && timestep%analysis_interval == 0 && USE_MORPH){
analysis.finish();
@ -1138,7 +1127,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
double volA = Averages->gnb.V;
volA /= Dm->Volume;
volB /= Dm->Volume;;
initial_volume = volA*Dm->Volume;
//initial_volume = volA*Dm->Volume;
double vA_x = Averages->gnb.Px/Averages->gnb.M;
double vA_y = Averages->gnb.Py/Averages->gnb.M;
double vA_z = Averages->gnb.Pz/Averages->gnb.M;
@ -1163,10 +1152,22 @@ void ScaLBL_GreyscaleColorModel::Run(){
double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha);
if (SET_CAPILLARY_NUMBER && CURRENT_STEADY_TIMESTEPS%MIN_STEADY_TIMESTEPS < analysis_interval ){
Fx *= capillary_number / Ca;
Fy *= capillary_number / Ca;
Fz *= capillary_number / Ca;
if ( morph_timesteps > morph_interval ){
bool isSteady = false;
if ( (fabs((Ca - Ca_previous)/Ca) < tolerance && CURRENT_STEADY_TIMESTEPS > MIN_STEADY_TIMESTEPS))
isSteady = true;
if (CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS)
isSteady = true;
if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && CURRENT_STEADY_TIMESTEPS > RESCALE_FORCE_AFTER_TIMESTEP){
RESCALE_FORCE = false;
double RESCALE_FORCE_FACTOR = capillary_number / Ca;
if (RESCALE_FORCE_FACTOR > 2.0) RESCALE_FORCE_FACTOR = 2.0;
if (RESCALE_FORCE_FACTOR < 0.5) RESCALE_FORCE_FACTOR = 0.5;
Fx *= RESCALE_FORCE_FACTOR;
Fy *= RESCALE_FORCE_FACTOR;
Fz *= RESCALE_FORCE_FACTOR;
force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
if (force_mag > 1e-3){
Fx *= 1e-3/force_mag; // impose ceiling for stability
Fy *= 1e-3/force_mag;
@ -1176,28 +1177,31 @@ void ScaLBL_GreyscaleColorModel::Run(){
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
greyscaleColor_db->putVector<double>("F",{Fx,Fy,Fz});
}
if ( morph_timesteps > morph_interval ){
bool isSteady = false;
if ( (fabs((Ca - Ca_previous)/Ca) < tolerance && CURRENT_STEADY_TIMESTEPS > MIN_STEADY_TIMESTEPS))
isSteady = true;
if (CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS)
isSteady = true;
if ( isSteady ){
MORPH_ADAPT = true;
CURRENT_MORPH_TIMESTEPS=0;
//delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
/** morphological target based on relative permeability for A **/
delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
//****** ENDPOINT ADAPTATION ********/
double krA_TMP= fabs(muA*flow_rate_A / force_mag);
double krB_TMP= fabs(muB*flow_rate_B / force_mag);
log_krA = log(krA_TMP);
if (krA_TMP < 0.0){
// cannot do endpoint adaptation if kr is negative
log_krA = log_krA_prev;
}
else if (krA_TMP < krB_TMP && morph_delta > 0.0){
/** morphological target based on relative permeability for A **/
log_krA_target = log(KRA_MORPH_FACTOR*(krA_TMP));
slope_krA_volume = (log_krA - log_krA_prev)/(Dm->Volume*(volA - volA_prev));
delta_volume_target=Dm->Volume*(volA+(log_krA_target - log_krA)/slope_krA_volume);
delta_volume_target=min(delta_volume_target,Dm->Volume*(volA+(log_krA_target - log_krA)/slope_krA_volume));
if (rank==0){
printf(" Enabling endpoint adaptation: krA = %f, krB = %f \n",krA_TMP,krB_TMP);
printf(" log(kr)=%f, volume=%f, TARGET log(kr)=%f, volume change=%f \n",log_krA, volA, log_krA_target, delta_volume_target/(volA*Dm->Volume));
}
}
log_krA_prev = log_krA;
volA_prev = volA;
printf(" log(kr)=%f, volume=%f, TARGET log(kr)=%f, volume change=%f \n",log_krA, volA, log_krA_target, delta_volume_target/(volA*Dm->Volume));
//******************************** **/
/** compute averages & write data **/
Averages->Full();
Averages->Write(timestep);
@ -1213,8 +1217,8 @@ void ScaLBL_GreyscaleColorModel::Run(){
double pB = Averages->gwb.p;
double pAc = Averages->gnc.p;
double pBc = Averages->gwc.p;
double pAB = (pA-pB)/(h*5.796*alpha);
double pAB_connected = (pAc-pBc)/(h*5.796*alpha);
double pAB = (pA-pB)/(h*6.0*alpha);
double pAB_connected = (pAc-pBc)/(h*6.0*alpha);
// connected contribution
double Vol_nc = Averages->gnc.V/Dm->Volume;
double Vol_wc = Averages->gwc.V/Dm->Volume;
@ -1277,16 +1281,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
Fy *= 1e-3/force_mag;
Fz *= 1e-3/force_mag;
}
if (flow_rate_A < NOISE_THRESHOLD && USE_BUMP_RATE){
if (rank==0) printf("Hit noise threshold (%f): bumping capillary number by %f X \n",NOISE_THRESHOLD,BUMP_RATE);
Fx *= BUMP_RATE; // impose bump condition
Fy *= BUMP_RATE;
Fz *= BUMP_RATE;
capillary_number *= BUMP_RATE;
greyscaleColor_db->putScalar<int>("capillary_number",capillary_number);
current_db->putDatabase("Color", greyscaleColor_db);
MORPH_ADAPT = false; // re-run current point if below noise threshold
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
greyscaleColor_db->putVector<double>("F",{Fx,Fy,Fz});
@ -1314,7 +1308,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
std::string next_image = ImageList[IMAGE_INDEX];
if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX);
greyscaleColor_db->putScalar<int>("image_index",IMAGE_INDEX);
//ImageInit(next_image);
ImageInit(next_image);
}
else{
if (rank==0) printf("Finished simulating image sequence \n");
@ -1343,42 +1337,30 @@ void ScaLBL_GreyscaleColorModel::Run(){
CURRENT_STEADY_TIMESTEPS=0;
initial_volume = volA*Dm->Volume;
delta_volume = 0.0;
if (USE_DIRECT){
//BoundaryCondition = 0;
//ScaLBL_Comm->BoundaryCondition = 0;
//ScaLBL_Comm_Regular->BoundaryCondition = 0;
//Fx = capillary_number*dir_x*force_mag / Ca;
//Fy = capillary_number*dir_y*force_mag / Ca;
//Fz = capillary_number*dir_z*force_mag / Ca;
}
if (RESCALE_FORCE_AFTER_TIMESTEP > 0)
RESCALE_FORCE = true;
}
else if (!(USE_DIRECT) && CURRENT_MORPH_TIMESTEPS > MAX_MORPH_TIMESTEPS) {
MORPH_ADAPT = false;
CURRENT_STEADY_TIMESTEPS=0;
initial_volume = volA*Dm->Volume;
delta_volume = 0.0;
RESCALE_FORCE = true;
if (RESCALE_FORCE_AFTER_TIMESTEP > 0)
RESCALE_FORCE = true;
}
if ( REVERSE_FLOW_DIRECTION ){
//if (rank==0) printf("*****REVERSE FLOW DIRECTION***** \n");
delta_volume = 0.0;
// flow direction will reverse after next steady point
MORPH_ADAPT = false;
CURRENT_STEADY_TIMESTEPS=0;
//morph_delta *= (-1.0);
REVERSE_FLOW_DIRECTION = false;
}
comm.barrier();
}
morph_timesteps += analysis_interval;
}
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
}
analysis.finish();
PROFILE_STOP("Loop");
PROFILE_SAVE("lbpm_color_simulator",1);
//************************************************************************
ScaLBL_DeviceBarrier();
comm.barrier();
stoptime = Utilities::MPI::time();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
stoptime = MPI_Wtime();
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
@ -1395,52 +1377,47 @@ void ScaLBL_GreyscaleColorModel::Run(){
// ************************************************************************
}
//double ScaLBL_GreyscaleColorModel::ImageInit(std::string Filename){
//
// if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str());
// Mask->Decomp(Filename);
// for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
// for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; // save what was read
//
// //double *PhaseLabel;
// //PhaseLabel = new double[Nx*Ny*Nz];
// //AssignComponentLabels(PhaseLabel);
// AssignComponentLabels();
//
// double Count = 0.0;
// double PoreCount = 0.0;
// for (int k=1; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
// if (id[Nx*Ny*k+Nx*j+i] == 2){
// PoreCount++;
// Count++;
// }
// else if (id[Nx*Ny*k+Nx*j+i] == 1){
// PoreCount++;
// }
// }
// }
// }
//
// Count = Dm->Comm.sumReduce( Count );
// PoreCount = Dm->Comm.sumReduce( PoreCount );
//
// if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
// ScaLBL_CopyToDevice(Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double));
// comm.barrier();
//
// ScaLBL_D3Q19_Init(fq, Np);
// ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
// ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
// comm.barrier();
//
// ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double));
//
// double saturation = Count/PoreCount;
// return saturation;
//
//}
double ScaLBL_GreyscaleColorModel::ImageInit(std::string Filename){
if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str());
Mask->Decomp(Filename);
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; // save what was read
AssignComponentLabels();
double Count = 0.0;
double PoreCount = 0.0;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
if (id[Nx*Ny*k+Nx*j+i] == 2){
PoreCount++;
Count++;
}
else if (id[Nx*Ny*k+Nx*j+i] == 1){
PoreCount++;
}
}
}
}
Count=sumReduce( Dm->Comm, Count);
PoreCount=sumReduce( Dm->Comm, PoreCount);
if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
ScaLBL_D3Q19_Init(fq, Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double));
double saturation = Count/PoreCount;
return saturation;
}
double ScaLBL_GreyscaleColorModel::MorphOpenConnected(double target_volume_change){
@ -1465,7 +1442,7 @@ double ScaLBL_GreyscaleColorModel::MorphOpenConnected(double target_volume_chang
BlobIDstruct new_index;
double vF=0.0; double vS=0.0;
ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm);
Dm->Comm.barrier();
MPI_Barrier(Dm->Comm);
long long count_connected=0;
long long count_porespace=0;
@ -1487,9 +1464,9 @@ double ScaLBL_GreyscaleColorModel::MorphOpenConnected(double target_volume_chang
}
}
}
count_connected = Dm->Comm.sumReduce( count_connected);
count_porespace = Dm->Comm.sumReduce( count_porespace);
count_water = Dm->Comm.sumReduce( count_water);
count_connected=sumReduce( Dm->Comm, count_connected);
count_porespace=sumReduce( Dm->Comm, count_porespace);
count_water=sumReduce( Dm->Comm, count_water);
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
@ -1561,7 +1538,7 @@ double ScaLBL_GreyscaleColorModel::MorphOpenConnected(double target_volume_chang
}
}
}
count_morphopen = Dm->Comm.sumReduce( count_morphopen);
count_morphopen=sumReduce( Dm->Comm, count_morphopen);
volume_change = double(count_morphopen - count_connected);
if (rank==0) printf(" opening of connected oil %f \n",volume_change/count_connected);
@ -1569,7 +1546,7 @@ double ScaLBL_GreyscaleColorModel::MorphOpenConnected(double target_volume_chang
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition >0 ){
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);
@ -1584,32 +1561,21 @@ double ScaLBL_GreyscaleColorModel::MorphOpenConnected(double target_volume_chang
}
return(volume_change);
}
double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil){
srand(time(NULL));
double mass_loss =0.f;
double count =0.f;
double *Aq_tmp, *Bq_tmp;
double *Vel_tmp;
Aq_tmp = new double [7*Np];
Bq_tmp = new double [7*Np];
Vel_tmp = new double [3*Np];
ScaLBL_CopyToHost(Aq_tmp, Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, Bq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Vel_tmp, Velocity, 7*Np*sizeof(double));
//Extract averged velocity
double vx_glb = (Averages->gnb.Px+Averages->gwb.Px)/(Averages->gnb.M+Averages->gwb.M);
double vy_glb = (Averages->gnb.Py+Averages->gwb.Py)/(Averages->gnb.M+Averages->gwb.M);
double vz_glb = (Averages->gnb.Pz+Averages->gwb.Pz)/(Averages->gnb.M+Averages->gwb.M);
double v_mag_glb = sqrt(vx_glb*vx_glb+vy_glb*vy_glb+vz_glb*vz_glb);
for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
double v_mag_local = sqrt(Vel_tmp[n]*Vel_tmp[n]+Vel_tmp[n+1*Np]*Vel_tmp[n+1*Np]+Vel_tmp[n+2*Np]*Vel_tmp[n+2*Np]);
double weight = (v_mag_local<v_mag_glb) ? v_mag_local/v_mag_glb : 1.0;
double random_value = weight*seed_water_in_oil*double(rand())/ RAND_MAX;
double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
double phase_id = (dA - dB) / (dA + dB);
@ -1634,9 +1600,7 @@ double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil
}
for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){
double v_mag_local = sqrt(Vel_tmp[n]*Vel_tmp[n]+Vel_tmp[n+1*Np]*Vel_tmp[n+1*Np]+Vel_tmp[n+2*Np]*Vel_tmp[n+2*Np]);
double weight = (v_mag_local<v_mag_glb) ? v_mag_local/v_mag_glb : 1.0;
double random_value = weight*seed_water_in_oil*double(rand())/ RAND_MAX;
double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
double phase_id = (dA - dB) / (dA + dB);
@ -1660,8 +1624,8 @@ double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil
mass_loss += random_value*seed_water_in_oil;
}
count = Dm->Comm.sumReduce( count );
mass_loss = Dm->Comm.sumReduce( mass_loss );
count= sumReduce( Dm->Comm, count);
mass_loss= sumReduce( Dm->Comm, mass_loss);
if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count);
// Need to initialize Aq, Bq, Den, Phi directly
@ -1672,106 +1636,14 @@ double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil
return(mass_loss);
}
//double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil){
// srand(time(NULL));
// double mass_loss =0.f;
// double count =0.f;
// double *Aq_tmp, *Bq_tmp;
//
// Aq_tmp = new double [7*Np];
// Bq_tmp = new double [7*Np];
//
// ScaLBL_CopyToHost(Aq_tmp, Aq, 7*Np*sizeof(double));
// ScaLBL_CopyToHost(Bq_tmp, Bq, 7*Np*sizeof(double));
//
///* for (int k=1; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
// double random_value = double(rand())/ RAND_MAX;
//
// if (Averages->SDs(i,j,k) < 0.f){
// // skip
// }
// else if (phase(i,j,k) > 0.f ){
// phase(i,j,k) -= random_value*seed_water_in_oil;
// mass_loss += random_value*seed_water_in_oil;
// count++;
// }
// else {
//
// }
// }
// }
// }
// */
// for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
// double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
// double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
// double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
// double phase_id = (dA - dB) / (dA + dB);
// if (phase_id > 0.0){
// Aq_tmp[n] -= 0.3333333333333333*random_value;
// Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
//
// Bq_tmp[n] += 0.3333333333333333*random_value;
// Bq_tmp[n+Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
// }
// mass_loss += random_value*seed_water_in_oil;
// }
//
// for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){
// double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
// double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
// double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
// double phase_id = (dA - dB) / (dA + dB);
// if (phase_id > 0.0){
// Aq_tmp[n] -= 0.3333333333333333*random_value;
// Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
// Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
//
// Bq_tmp[n] += 0.3333333333333333*random_value;
// Bq_tmp[n+Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
// Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
// }
// mass_loss += random_value*seed_water_in_oil;
// }
//
// count = Dm->Comm.sumReduce( count );
// mass_loss = Dm->Comm.sumReduce( mass_loss );
// if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count);
//
// // Need to initialize Aq, Bq, Den, Phi directly
// //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
// ScaLBL_CopyToDevice(Aq, Aq_tmp, 7*Np*sizeof(double));
// ScaLBL_CopyToDevice(Bq, Bq_tmp, 7*Np*sizeof(double));
//
// return(mass_loss);
//}
double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double target_delta_volume){
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
double vF = 0.f;
double vS = 0.f;
double delta_volume;
double WallFactor = 0.0;
bool USE_CONNECTED_NWP = false;
DoubleArray phase(Nx,Ny,Nz);
IntArray phase_label(Nx,Ny,Nz);;
@ -1792,7 +1664,7 @@ double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double tar
}
}
}
double volume_initial = Dm->Comm.sumReduce( count);
double volume_initial = sumReduce( Dm->Comm, count);
/*
sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank);
FILE *INPUT = fopen(LocalRankFilename,"wb");
@ -1800,13 +1672,16 @@ double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double tar
fclose(INPUT);
*/
// 2. Identify connected components of phase field -> phase_label
double volume_connected = 0.0;
double second_biggest = 0.0;
if (USE_CONNECTED_NWP){
BlobIDstruct new_index;
ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm);
comm.barrier();
MPI_Barrier(Dm->Comm);
// only operate on component "0"
count = 0.0;
double second_biggest = 0.0;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
@ -1824,8 +1699,29 @@ double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double tar
}
}
}
double volume_connected = Dm->Comm.sumReduce( count );
second_biggest = Dm->Comm.sumReduce( second_biggest );
volume_connected = sumReduce( Dm->Comm, count);
second_biggest = sumReduce( Dm->Comm, second_biggest);
}
else {
// use the whole NWP
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (Averages->SDs(i,j,k) > 0.f){
if (phase(i,j,k) > 0.f ){
phase_id(i,j,k) = 0;
}
else {
phase_id(i,j,k) = 1;
}
}
else {
phase_id(i,j,k) = 1;
}
}
}
}
}
/*int reach_x, reach_y, reach_z;
for (int k=0; k<Nz; k++){
@ -1860,18 +1756,21 @@ double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double tar
}
}
if (USE_CONNECTED_NWP){
if (volume_connected - second_biggest < 2.0*fabs(target_delta_volume) && target_delta_volume < 0.0){
// if connected volume is less than 2% just delete the whole thing
if (rank==0) printf("Connected region has shrunk! \n");
REVERSE_FLOW_DIRECTION = true;
}
/* else{*/
if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest );
}
if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
double target_delta_volume_incremental = target_delta_volume;
if (fabs(target_delta_volume) > 0.01*volume_initial)
target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental);
delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor);
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
@ -1912,7 +1811,7 @@ double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double tar
}
}
}
double volume_final = Dm->Comm.sumReduce( count );
double volume_final= sumReduce( Dm->Comm, count);
delta_volume = (volume_final-volume_initial);
if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial);
@ -1935,7 +1834,7 @@ double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double tar
// 7. Re-initialize phase field and density
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition >0 ){
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);
@ -2039,27 +1938,27 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
fwrite(PhaseField.data(),8,N,GreySG_Z_FILE);
fclose(GreySG_Z_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);
/* 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);
*/
}
//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-1

10
models/GreyscaleColorModel.h
View File

@ -1,5 +1,5 @@
/*
Implementation of color lattice boltzmann model
Implementation of two-fluid greyscale color lattice boltzmann model
*/
#include <stdio.h>
#include <stdlib.h>
@ -12,13 +12,13 @@ Implementation of color lattice boltzmann model
#include "common/Communication.h"
#include "analysis/TwoPhase.h"
#include "analysis/runAnalysis.h"
#include "common/MPI.h"
#include "common/MPI_Helpers.h"
#include "ProfilerApp.h"
#include "threadpool/thread_pool.h"
class ScaLBL_GreyscaleColorModel{
public:
ScaLBL_GreyscaleColorModel(int RANK, int NP, const Utilities::MPI& COMM);
ScaLBL_GreyscaleColorModel(int RANK, int NP, MPI_Comm COMM);
~ScaLBL_GreyscaleColorModel();
// functions in they should be run
@ -67,7 +67,7 @@ public:
double *fq, *Aq, *Bq;
double *Den, *Phi;
//double *GreySolidPhi; //Model 2 & 3
double *GreySolidGrad;//Model 1
double *GreySolidGrad;//Model 1 & 4
//double *ColorGrad;
double *Velocity;
double *Pressure;
@ -75,7 +75,7 @@ public:
double *Permeability_dvc;
private:
Utilities::MPI comm;
MPI_Comm comm;
int dist_mem_size;
int neighborSize;

1247
models/GreyscaleFEModel.cpp
View File

File diff suppressed because it is too large Load Diff

108
models/GreyscaleFEModel.h
View File

@ -1,108 +0,0 @@
/*
Implementation of multicomponent greyscale free-energy based lattice boltzmann model
*/
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/Communication.h"
#include "common/MPI.h"
#include "common/Database.h"
#include "common/ScaLBL.h"
#include "ProfilerApp.h"
#include "threadpool/thread_pool.h"
class ScaLBL_GreyscaleFEModel{
public:
ScaLBL_GreyscaleFEModel(int RANK, int NP, MPI_Comm COMM);
~ScaLBL_GreyscaleFEModel();
// functions in they should be run
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
void SetDomain();
void ReadInput();
void Create();
void Initialize();
void Run();
void WriteDebug();
void VelocityField();
bool Restart,pBC;
int timestep,timestepMax;
int BoundaryCondition;
double tauA,tauB;
double tauA_eff,tauB_eff;
double rhoA,rhoB;
double tolerance;
double Fx,Fy,Fz,flux;
double din,dout;
double dp;//solid particle diameter, unit in voxel
double GreyPorosity;
//double Gsc;
double gamma;
double kappaA,kappaB;
double lambdaA,lambdaB;
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
double Lx,Ly,Lz;
std::shared_ptr<Domain> Dm; // this domain is for analysis
std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
// input database
std::shared_ptr<Database> db;
std::shared_ptr<Database> domain_db;
std::shared_ptr<Database> greyscaleFE_db;
std::shared_ptr<Database> analysis_db;
std::shared_ptr<Database> vis_db;
signed char *id;
int *NeighborList;
//double *fq,*Aq,*Bq;
double *fq,*Cq;
double *Den;
double *Permeability;//grey voxel permeability
double *Porosity;
double *Velocity;
double *SolidForce;
double *Pressure_dvc;
double *PressureGrad;// gradiant of pressure
double *PressTensor;//pressure tensor
double *PressTensorGrad;// gradient of pressure tensor
double *Phi;
double *PhiLap;//laplacian of phase field phi
// double *DenGradA;
// double *DenGradB;
// double *DenLapA;
// double *DenLapB;
IntArray Map;
DoubleArray SignDist;
DoubleArray Velocity_x;
DoubleArray Velocity_y;
DoubleArray Velocity_z;
DoubleArray PorosityMap;
DoubleArray Pressure;
private:
MPI_Comm comm;
int dist_mem_size;
int neighborSize;
// filenames
char LocalRankString[8];
char LocalRankFilename[40];
char LocalRestartFile[40];
void AssignComponentLabels(double *Porosity, double *Permeablity, double *SolidPotential);
void AssignSolidForce(double *SolidPotential, double *SolidForce);
void Density_and_Phase_Init();
};

13
models/GreyscaleModel.cpp
View File

@ -584,15 +584,10 @@ void ScaLBL_GreyscaleModel::Run(){
}
}
}
// Disable the this, if the one below does not work, then enable this
//vax = sumReduce( Mask->Comm, vax_loc);
//vay = sumReduce( Mask->Comm, vay_loc);
//vaz = sumReduce( Mask->Comm, vaz_loc);
//count = sumReduce( Mask->Comm, count_loc);
vax = Mask->Comm.sumReduce( vax_loc );
vay = Mask->Comm.sumReduce( vay_loc );
vaz = Mask->Comm.sumReduce( vaz_loc );
count = Mask->Comm.sumReduce( count_loc );
vax = sumReduce( Mask->Comm, vax_loc);
vay = sumReduce( Mask->Comm, vay_loc);
vaz = sumReduce( Mask->Comm, vaz_loc);
count = sumReduce( Mask->Comm, count_loc);
vax /= count;
vay /= count;

1433
models/GreyscaleSCModel.cpp
View File

File diff suppressed because it is too large Load Diff

103
models/GreyscaleSCModel.h
View File

@ -1,103 +0,0 @@
/*
Implementation of color lattice boltzmann model
*/
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/Communication.h"
#include "common/MPI.h"
#include "common/Database.h"
#include "common/ScaLBL.h"
#include "ProfilerApp.h"
#include "threadpool/thread_pool.h"
class ScaLBL_GreyscaleSCModel{
public:
ScaLBL_GreyscaleSCModel(int RANK, int NP, MPI_Comm COMM);
~ScaLBL_GreyscaleSCModel();
// functions in they should be run
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
void SetDomain();
void ReadInput();
void Create();
void Initialize();
void Run();
void WriteDebug();
void WriteOutput();
bool Restart,pBC;
int timestep,timestepMax;
int BoundaryCondition;
int CollisionType;
double tauA,tauB;
double tauA_eff,tauB_eff;
double Gsc;
double rhoA,rhoB;
double rhoA_minor,rhoB_minor;//dissolved density
double tolerance;
double Fx,Fy,Fz;
double fluxA,fluxB;
double dinA,doutA;
double dinB,doutB;
double GreyPorosity;
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
double Lx,Ly,Lz;
std::shared_ptr<Domain> Dm; // this domain is for analysis
std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
// input database
std::shared_ptr<Database> db;
std::shared_ptr<Database> domain_db;
std::shared_ptr<Database> greyscaleSC_db;
std::shared_ptr<Database> analysis_db;
std::shared_ptr<Database> vis_db;
signed char *id;
int *NeighborList;
int *dvcMap;
double *fqA, *fqB;
double *Permeability;//grey voxel permeability
//double relPermA,relPermB;//grey voxel relperm
double *Porosity;
double *Velocity;
double *Pressure_dvc;
double *DenA, *DenB;
double *DenGradA,*DenGradB;
double *SolidForceA,*SolidForceB;
IntArray Map;
DoubleArray SignDist;
DoubleArray Velocity_x;
DoubleArray Velocity_y;
DoubleArray Velocity_z;
DoubleArray PorosityMap;
DoubleArray Pressure;
DoubleArray DenA_data;
DoubleArray DenB_data;
private:
MPI_Comm comm;
int dist_mem_size;
int neighborSize;
// filenames
char LocalRankString[8];
char LocalRankFilename[40];
char LocalRestartFile[40];
void AssignGreyscaleAndSolidLabels();
void Density_Init();
};

2
tests/CMakeLists.txt
View File

@ -4,8 +4,6 @@
ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
ADD_LBPM_EXECUTABLE( lbpm_greyscale_simulator )
ADD_LBPM_EXECUTABLE( lbpm_greyscaleFE_simulator )
ADD_LBPM_EXECUTABLE( lbpm_greyscaleSC_simulator )
ADD_LBPM_EXECUTABLE( lbpm_greyscaleColor_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )

59
tests/lbpm_greyscaleFE_simulator.cpp
View File

@ -1,59 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI.h"
#include "models/GreyscaleFEModel.h"
//#define WRITE_SURFACES
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);
{
// parallel domain size (# of sub-domains)
int nprocx,nprocy,nprocz;
int iproc,jproc,kproc;
if (rank == 0){
printf("****************************************\n");
printf("Running Greyscale Two-Phase Calculation \n");
printf("****************************************\n");
}
// Initialize compute device
int device=ScaLBL_SetDevice(rank);
ScaLBL_DeviceBarrier();
MPI_Barrier(comm);
ScaLBL_GreyscaleFEModel GreyscaleFE(rank,nprocs,comm);
auto filename = argv[1];
GreyscaleFE.ReadParams(filename);
GreyscaleFE.SetDomain(); // this reads in the domain
GreyscaleFE.ReadInput();
GreyscaleFE.Create(); // creating the model will create data structure to match the pore structure and allocate variables
GreyscaleFE.Initialize(); // initializing the model will set initial conditions for variables
GreyscaleFE.Run();
//GreyscaleFE.VelocityField();
GreyscaleFE.WriteDebug();
}
// ****************************************************
MPI_Barrier(comm);
MPI_Finalize();
// ****************************************************
}

58
tests/lbpm_greyscaleSC_simulator.cpp
View File

@ -1,58 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI.h"
#include "models/GreyscaleSCModel.h"
//#define WRITE_SURFACES
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);
{
// parallel domain size (# of sub-domains)
int nprocx,nprocy,nprocz;
int iproc,jproc,kproc;
if (rank == 0){
printf("****************************************\n");
printf("Running Greyscale Two-Phase Calculation \n");
printf("****************************************\n");
}
// Initialize compute device
int device=ScaLBL_SetDevice(rank);
ScaLBL_DeviceBarrier();
MPI_Barrier(comm);
ScaLBL_GreyscaleSCModel GreyscaleSC(rank,nprocs,comm);
auto filename = argv[1];
GreyscaleSC.ReadParams(filename);
GreyscaleSC.SetDomain(); // this reads in the domain
GreyscaleSC.ReadInput();
GreyscaleSC.Create(); // creating the model will create data structure to match the pore structure and allocate variables
GreyscaleSC.Initialize(); // initializing the model will set initial conditions for variables
GreyscaleSC.Run();
GreyscaleSC.WriteDebug();
}
// ****************************************************
MPI_Barrier(comm);
MPI_Finalize();
// ****************************************************
}