diff --git a/common/Domain.cpp b/common/Domain.cpp index 2e808411..f0553011 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -1543,7 +1543,7 @@ void Domain::ReadFromFile(const std::string &Filename, } else { // Recieve the subdomain from rank = 0 //printf("Ready to recieve data %i at process %i \n", N,rank); - Comm.recv(id.data(), N, 0, 15); + Comm.recv(UserData, N, 0, 15); } Comm.barrier(); } diff --git a/common/Membrane.cpp b/common/Membrane.cpp new file mode 100644 index 00000000..999c9389 --- /dev/null +++ b/common/Membrane.cpp @@ -0,0 +1,1415 @@ +/* Flow adaptor class for multiphase flow methods */ + +#include "common/Membrane.h" +#include "analysis/distance.h" + +Membrane::Membrane(std::shared_ptr Dm, int *dvcNeighborList, int Nsites) { + + Np = Nsites; + initialNeighborList = new int[18*Np]; + ScaLBL_AllocateDeviceMemory((void **)&NeighborList, 18*Np*sizeof(int)); + Lock=false; // unlock the communicator + //...................................................................................... + // Create a separate copy of the communicator for the device + MPI_COMM_SCALBL = Dm->Comm.dup(); + + ScaLBL_CopyToHost(initialNeighborList, dvcNeighborList, 18*Np*sizeof(int)); + Dm->Comm.barrier(); + ScaLBL_CopyToDevice(NeighborList, initialNeighborList, 18*Np*sizeof(int)); + + /* Copy communication lists */ + //...................................................................................... + //Lock=false; // unlock the communicator + //...................................................................................... + // Create a separate copy of the communicator for the device + //MPI_COMM_SCALBL = Dm->Comm.dup(); + //...................................................................................... + // Copy the domain size and communication information directly from Dm + Nx = Dm->Nx; + Ny = Dm->Ny; + Nz = Dm->Nz; + N = Nx*Ny*Nz; + //next=0; + rank=Dm->rank(); + rank_x=Dm->rank_x(); + rank_y=Dm->rank_y(); + rank_z=Dm->rank_z(); + rank_X=Dm->rank_X(); + rank_Y=Dm->rank_Y(); + rank_Z=Dm->rank_Z(); + rank_xy=Dm->rank_xy(); + rank_XY=Dm->rank_XY(); + rank_xY=Dm->rank_xY(); + rank_Xy=Dm->rank_Xy(); + rank_xz=Dm->rank_xz(); + rank_XZ=Dm->rank_XZ(); + rank_xZ=Dm->rank_xZ(); + rank_Xz=Dm->rank_Xz(); + rank_yz=Dm->rank_yz(); + rank_YZ=Dm->rank_YZ(); + rank_yZ=Dm->rank_yZ(); + rank_Yz=Dm->rank_Yz(); + sendCount_x=Dm->sendCount("x"); + sendCount_y=Dm->sendCount("y"); + sendCount_z=Dm->sendCount("z"); + sendCount_X=Dm->sendCount("X"); + sendCount_Y=Dm->sendCount("Y"); + sendCount_Z=Dm->sendCount("Z"); + sendCount_xy=Dm->sendCount("xy"); + sendCount_yz=Dm->sendCount("yz"); + sendCount_xz=Dm->sendCount("xz"); + sendCount_Xy=Dm->sendCount("Xy"); + sendCount_Yz=Dm->sendCount("Yz"); + sendCount_xZ=Dm->sendCount("xZ"); + sendCount_xY=Dm->sendCount("xY"); + sendCount_yZ=Dm->sendCount("yZ"); + sendCount_Xz=Dm->sendCount("Xz"); + sendCount_XY=Dm->sendCount("XY"); + sendCount_YZ=Dm->sendCount("YZ"); + sendCount_XZ=Dm->sendCount("XZ"); + recvCount_x=Dm->recvCount("x"); + recvCount_y=Dm->recvCount("y"); + recvCount_z=Dm->recvCount("z"); + recvCount_X=Dm->recvCount("X"); + recvCount_Y=Dm->recvCount("Y"); + recvCount_Z=Dm->recvCount("Z"); + recvCount_xy=Dm->recvCount("xy"); + recvCount_yz=Dm->recvCount("yz"); + recvCount_xz=Dm->recvCount("xz"); + recvCount_Xy=Dm->recvCount("Xy"); + recvCount_Yz=Dm->recvCount("Yz"); + recvCount_xZ=Dm->recvCount("xZ"); + recvCount_xY=Dm->recvCount("xY"); + recvCount_yZ=Dm->recvCount("yZ"); + recvCount_Xz=Dm->recvCount("Xz"); + recvCount_XY=Dm->recvCount("XY"); + recvCount_YZ=Dm->recvCount("YZ"); + recvCount_XZ=Dm->recvCount("XZ"); + + if (rank == 0){ + printf("**** Creating membrane data structure ****** \n"); + } + printf(" Number of active lattice sites (rank = %i): %i \n",rank, Np); + + /* check symmetry for send / recv counts */ + if (sendCount_x != recvCount_X) printf("WARNING: rank %i send/recv mismatch (x/X)! \n",rank); + if (sendCount_y != recvCount_Y) printf("WARNING: rank %i send/recv mismatch (y/Y)! \n",rank); + if (sendCount_z != recvCount_Z) printf("WARNING: rank %i send/recv mismatch (z/Z)! \n",rank); + if (sendCount_X != recvCount_x) printf("WARNING: rank %i send/recv mismatch (X/x)! \n",rank); + if (sendCount_Y != recvCount_y) printf("WARNING: rank %i send/recv mismatch (Y/y)! \n",rank); + if (sendCount_x != recvCount_z) printf("WARNING: rank %i send/recv mismatch (Z/z)! \n",rank); + if (sendCount_xy != recvCount_XY) printf("WARNING: rank %i send/recv mismatch (xy/XY)! \n",rank); + if (sendCount_Xy != recvCount_xY) printf("WARNING: rank %i send/recv mismatch (Xy/xY)! \n",rank); + if (sendCount_xY != recvCount_Xy) printf("WARNING: rank %i send/recv mismatch (xY/Xy)! \n",rank); + if (sendCount_XY != recvCount_xy) printf("WARNING: rank %i send/recv mismatch (XY/xy)! \n",rank); + if (sendCount_xz != recvCount_XZ) printf("WARNING: rank %i send/recv mismatch (xz/XZ)! \n",rank); + if (sendCount_Xz != recvCount_xZ) printf("WARNING: rank %i send/recv mismatch (Xz/xZ)! \n",rank); + if (sendCount_xZ != recvCount_Xz) printf("WARNING: rank %i send/recv mismatch (xZ/Xz)! \n",rank); + if (sendCount_XZ != recvCount_xz) printf("WARNING: rank %i send/recv mismatch (XZ/xz)! \n",rank); + if (sendCount_yz != recvCount_YZ) printf("WARNING: rank %i send/recv mismatch (yz/YZ)! \n",rank); + if (sendCount_Yz != recvCount_yZ) printf("WARNING: rank %i send/recv mismatch (Yz/yZ)! \n",rank); + if (sendCount_yZ != recvCount_Yz) printf("WARNING: rank %i send/recv mismatch (yZ/Yz)! \n",rank); + if (sendCount_YZ != recvCount_yz) printf("WARNING: rank %i send/recv mismatch (YZ/yz)! \n",rank); + + iproc = Dm->iproc(); + jproc = Dm->jproc(); + kproc = Dm->kproc(); + nprocx = Dm->nprocx(); + nprocy = Dm->nprocy(); + nprocz = Dm->nprocz(); + //BoundaryCondition = Dm->BoundaryCondition; + //...................................................................................... + + ScaLBL_AllocateZeroCopy((void **) &sendbuf_x, 2*5*sendCount_x*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_X, 2*5*sendCount_X*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_y, 2*5*sendCount_y*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Y, 2*5*sendCount_Y*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_z, 2*5*sendCount_z*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Z, 2*5*sendCount_Z*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xy, 2*sendCount_xy*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xY, 2*sendCount_xY*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Xy, 2*sendCount_Xy*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_XY, 2*sendCount_XY*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xz, 2*sendCount_xz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_xZ, 2*sendCount_xZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Xz, 2*sendCount_Xz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_XZ, 2*sendCount_XZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_yz, 2*sendCount_yz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_yZ, 2*sendCount_yZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_Yz, 2*sendCount_Yz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &sendbuf_YZ, 2*sendCount_YZ*sizeof(double)); // Allocate device memory + //...................................................................................... + ScaLBL_AllocateZeroCopy((void **) &recvbuf_x, 2*5*recvCount_x*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_X, 2*5*recvCount_X*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_y, 2*5*recvCount_y*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Y, 2*5*recvCount_Y*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_z, 2*5*recvCount_z*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Z, 2*5*recvCount_Z*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xy, 2*recvCount_xy*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xY, 2*recvCount_xY*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Xy, 2*recvCount_Xy*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_XY, 2*recvCount_XY*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xz, 2*recvCount_xz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_xZ, 2*recvCount_xZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Xz, 2*recvCount_Xz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_XZ, 2*recvCount_XZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_yz, 2*recvCount_yz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_yZ, 2*recvCount_yZ*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_Yz, 2*recvCount_Yz*sizeof(double)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &recvbuf_YZ, 2*recvCount_YZ*sizeof(double)); // Allocate device memory + //...................................................................................... + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_x, sendCount_x*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_X, sendCount_X*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_y, sendCount_y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Y, sendCount_Y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_z, sendCount_z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Z, sendCount_Z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xy, sendCount_xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xY, sendCount_xY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Xy, sendCount_Xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_XY, sendCount_XY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xz, sendCount_xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_xZ, sendCount_xZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Xz, sendCount_Xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_XZ, sendCount_XZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_yz, sendCount_yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_yZ, sendCount_yZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_Yz, sendCount_Yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcSendList_YZ, sendCount_YZ*sizeof(int)); // Allocate device memory + //...................................................................................... + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_x, recvCount_x*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_X, recvCount_X*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_y, recvCount_y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Y, recvCount_Y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_z, recvCount_z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Z, recvCount_Z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xy, recvCount_xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xY, recvCount_xY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Xy, recvCount_Xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_XY, recvCount_XY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xz, recvCount_xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_xZ, recvCount_xZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Xz, recvCount_Xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_XZ, recvCount_XZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_yz, recvCount_yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_yZ, recvCount_yZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_Yz, recvCount_Yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvList_YZ, recvCount_YZ*sizeof(int)); // Allocate device memory + //...................................................................................... + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_x, 5*recvCount_x*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_X, 5*recvCount_X*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_y, 5*recvCount_y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_Y, 5*recvCount_Y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_z, 5*recvCount_z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_Z, 5*recvCount_Z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_xy, recvCount_xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_xY, recvCount_xY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_Xy, recvCount_Xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_XY, recvCount_XY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_xz, recvCount_xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_xZ, recvCount_xZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_Xz, recvCount_Xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_XZ, recvCount_XZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_yz, recvCount_yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_yZ, recvCount_yZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_Yz, recvCount_Yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvLinks_YZ, recvCount_YZ*sizeof(int)); // Allocate device memory + //...................................................................................... + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_x, 5*recvCount_x*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_X, 5*recvCount_X*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_y, 5*recvCount_y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_Y, 5*recvCount_Y*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_z, 5*recvCount_z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_Z, 5*recvCount_Z*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_xy, recvCount_xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_xY, recvCount_xY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_Xy, recvCount_Xy*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_XY, recvCount_XY*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_xz, recvCount_xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_xZ, recvCount_xZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_Xz, recvCount_Xz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_XZ, recvCount_XZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_yz, recvCount_yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_yZ, recvCount_yZ*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_Yz, recvCount_Yz*sizeof(int)); // Allocate device memory + ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_YZ, recvCount_YZ*sizeof(int)); // Allocate device memory + //...................................................................................... + + //...................................................................................... + ScaLBL_CopyToZeroCopy(dvcSendList_x,Dm->sendList("x"),sendCount_x*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_X,Dm->sendList("X"),sendCount_X*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_y,Dm->sendList("y"),sendCount_y*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_Y,Dm->sendList("Y"),sendCount_Y*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_z,Dm->sendList("z"),sendCount_z*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_Z,Dm->sendList("Z"),sendCount_Z*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_xy,Dm->sendList("xy"),sendCount_xy*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_XY,Dm->sendList("XY"),sendCount_XY*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_xY,Dm->sendList("xY"),sendCount_xY*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_Xy,Dm->sendList("Xy"),sendCount_Xy*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_xz,Dm->sendList("xz"),sendCount_xz*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_XZ,Dm->sendList("XZ"),sendCount_XZ*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_xZ,Dm->sendList("xZ"),sendCount_xZ*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_Xz,Dm->sendList("Xz"),sendCount_Xz*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_yz,Dm->sendList("yz"),sendCount_yz*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_YZ,Dm->sendList("YZ"),sendCount_YZ*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_yZ,Dm->sendList("yZ"),sendCount_yZ*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcSendList_Yz,Dm->sendList("Yz"),sendCount_Yz*sizeof(int)); + //...................................................................................... + ScaLBL_CopyToZeroCopy(dvcRecvList_x,Dm->recvList("x"),recvCount_x*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_X,Dm->recvList("X"),recvCount_X*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_y,Dm->recvList("y"),recvCount_y*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_Y,Dm->recvList("Y"),recvCount_Y*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_z,Dm->recvList("z"),recvCount_z*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_Z,Dm->recvList("Z"),recvCount_Z*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_xy,Dm->recvList("xy"),recvCount_xy*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_XY,Dm->recvList("XY"),recvCount_XY*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_xY,Dm->recvList("xY"),recvCount_xY*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_Xy,Dm->recvList("Xy"),recvCount_Xy*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_xz,Dm->recvList("xz"),recvCount_xz*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_XZ,Dm->recvList("XZ"),recvCount_XZ*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_xZ,Dm->recvList("xZ"),recvCount_xZ*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_Xz,Dm->recvList("Xz"),recvCount_Xz*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_yz,Dm->recvList("yz"),recvCount_yz*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_YZ,Dm->recvList("YZ"),recvCount_YZ*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_yZ,Dm->recvList("yZ"),recvCount_yZ*sizeof(int)); + ScaLBL_CopyToZeroCopy(dvcRecvList_Yz,Dm->recvList("Yz"),recvCount_Yz*sizeof(int)); + //...................................................................................... + +} + +Membrane::~Membrane() { + + delete [] initialNeighborList; + delete [] membraneLinks; + delete [] membraneTag; + delete [] membraneDist; + + ScaLBL_FreeDeviceMemory( coefficient_x ); + ScaLBL_FreeDeviceMemory( coefficient_X ); + ScaLBL_FreeDeviceMemory( coefficient_y ); + ScaLBL_FreeDeviceMemory( coefficient_Y ); + ScaLBL_FreeDeviceMemory( coefficient_z ); + ScaLBL_FreeDeviceMemory( coefficient_Z ); + + ScaLBL_FreeDeviceMemory( NeighborList ); + ScaLBL_FreeDeviceMemory( MembraneLinks ); + ScaLBL_FreeDeviceMemory( MembraneCoef ); + + ScaLBL_FreeDeviceMemory( sendbuf_x ); + ScaLBL_FreeDeviceMemory( sendbuf_X ); + ScaLBL_FreeDeviceMemory( sendbuf_y ); + ScaLBL_FreeDeviceMemory( sendbuf_Y ); + ScaLBL_FreeDeviceMemory( sendbuf_z ); + ScaLBL_FreeDeviceMemory( sendbuf_Z ); + ScaLBL_FreeDeviceMemory( sendbuf_xy ); + ScaLBL_FreeDeviceMemory( sendbuf_xY ); + ScaLBL_FreeDeviceMemory( sendbuf_Xy ); + ScaLBL_FreeDeviceMemory( sendbuf_XY ); + ScaLBL_FreeDeviceMemory( sendbuf_xz ); + ScaLBL_FreeDeviceMemory( sendbuf_xZ ); + ScaLBL_FreeDeviceMemory( sendbuf_Xz ); + ScaLBL_FreeDeviceMemory( sendbuf_XZ ); + ScaLBL_FreeDeviceMemory( sendbuf_yz ); + ScaLBL_FreeDeviceMemory( sendbuf_yZ ); + ScaLBL_FreeDeviceMemory( sendbuf_Yz ); + ScaLBL_FreeDeviceMemory( sendbuf_YZ ); + ScaLBL_FreeDeviceMemory( recvbuf_x ); + ScaLBL_FreeDeviceMemory( recvbuf_X ); + ScaLBL_FreeDeviceMemory( recvbuf_y ); + ScaLBL_FreeDeviceMemory( recvbuf_Y ); + ScaLBL_FreeDeviceMemory( recvbuf_z ); + ScaLBL_FreeDeviceMemory( recvbuf_Z ); + ScaLBL_FreeDeviceMemory( recvbuf_xy ); + ScaLBL_FreeDeviceMemory( recvbuf_xY ); + ScaLBL_FreeDeviceMemory( recvbuf_Xy ); + ScaLBL_FreeDeviceMemory( recvbuf_XY ); + ScaLBL_FreeDeviceMemory( recvbuf_xz ); + ScaLBL_FreeDeviceMemory( recvbuf_xZ ); + ScaLBL_FreeDeviceMemory( recvbuf_Xz ); + ScaLBL_FreeDeviceMemory( recvbuf_XZ ); + ScaLBL_FreeDeviceMemory( recvbuf_yz ); + ScaLBL_FreeDeviceMemory( recvbuf_yZ ); + ScaLBL_FreeDeviceMemory( recvbuf_Yz ); + ScaLBL_FreeDeviceMemory( recvbuf_YZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_x ); + ScaLBL_FreeDeviceMemory( dvcSendList_X ); + ScaLBL_FreeDeviceMemory( dvcSendList_y ); + ScaLBL_FreeDeviceMemory( dvcSendList_Y ); + ScaLBL_FreeDeviceMemory( dvcSendList_z ); + ScaLBL_FreeDeviceMemory( dvcSendList_Z ); + ScaLBL_FreeDeviceMemory( dvcSendList_xy ); + ScaLBL_FreeDeviceMemory( dvcSendList_xY ); + ScaLBL_FreeDeviceMemory( dvcSendList_Xy ); + ScaLBL_FreeDeviceMemory( dvcSendList_XY ); + ScaLBL_FreeDeviceMemory( dvcSendList_xz ); + ScaLBL_FreeDeviceMemory( dvcSendList_xZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_Xz ); + ScaLBL_FreeDeviceMemory( dvcSendList_XZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_yz ); + ScaLBL_FreeDeviceMemory( dvcSendList_yZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_Yz ); + ScaLBL_FreeDeviceMemory( dvcSendList_YZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_x ); + ScaLBL_FreeDeviceMemory( dvcRecvList_X ); + ScaLBL_FreeDeviceMemory( dvcRecvList_y ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Y ); + ScaLBL_FreeDeviceMemory( dvcRecvList_z ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Z ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xy ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xY ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Xy ); + ScaLBL_FreeDeviceMemory( dvcRecvList_XY ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Xz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_XZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_yz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_yZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Yz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_YZ ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_x ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_X ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_y ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_Y ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_z ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_Z ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_xy ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_xY ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_Xy ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_XY ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_xz ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_xZ ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_Xz ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_XZ ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_yz ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_yZ ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_Yz ); + ScaLBL_FreeDeviceMemory( dvcRecvLinks_YZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_x ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_X ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_y ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Y ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_z ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Z ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xy ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xY ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Xy ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_XY ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Xz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_XZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_yz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_yZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Yz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_YZ ); +} + +int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map){ + int mlink = 0; + int i,j,k,n; + int idx, neighbor; + double dist, locdist; + + if (rank == 0) printf(" Copy initial neighborlist... \n"); + int * neighborList = new int[18*Np]; + /* Copy neighborList */ + for (int idx=0; idxrecvList("X"),0,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); + linkCount_X[1] = D3Q19_MapRecv(-1,-1,0,Dm->recvList("X"),recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); + linkCount_X[2] = D3Q19_MapRecv(-1,1,0, Dm->recvList("X"),2*recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); + linkCount_X[3] = D3Q19_MapRecv(-1,0,-1,Dm->recvList("X"),3*recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); + linkCount_X[4] = D3Q19_MapRecv(-1,0,1, Dm->recvList("X"),4*recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); + //................................................................................... + //...Map recieve list for the x face: q=1,7,9,11,13.................................. + linkCount_x[0] = D3Q19_MapRecv(1,0,0, Dm->recvList("x"),0,recvCount_x,dvcRecvDist_x,dvcRecvLinks_x,Distance,Map); + linkCount_x[1] = D3Q19_MapRecv(1,1,0, Dm->recvList("x"),recvCount_x,recvCount_x,dvcRecvDist_x,dvcRecvLinks_x,Distance,Map); + linkCount_x[2] = D3Q19_MapRecv(1,-1,0,Dm->recvList("x"),2*recvCount_x,recvCount_x,dvcRecvDist_x,dvcRecvLinks_x,Distance,Map); + linkCount_x[3] = D3Q19_MapRecv(1,0,1, Dm->recvList("x"),3*recvCount_x,recvCount_x,dvcRecvDist_x,dvcRecvLinks_x,Distance,Map); + linkCount_x[4] = D3Q19_MapRecv(1,0,-1,Dm->recvList("x"),4*recvCount_x,recvCount_x,dvcRecvDist_x,dvcRecvLinks_x,Distance,Map); + //................................................................................... + //...Map recieve list for the y face: q=4,8,9,16,18 ................................... + linkCount_Y[0] = D3Q19_MapRecv(0,-1,0, Dm->recvList("Y"),0,recvCount_Y,dvcRecvDist_Y,dvcRecvLinks_Y,Distance,Map); + linkCount_Y[1] = D3Q19_MapRecv(-1,-1,0,Dm->recvList("Y"),recvCount_Y,recvCount_Y,dvcRecvDist_Y,dvcRecvLinks_Y,Distance,Map); + linkCount_Y[2] = D3Q19_MapRecv(1,-1,0, Dm->recvList("Y"),2*recvCount_Y,recvCount_Y,dvcRecvDist_Y,dvcRecvLinks_Y,Distance,Map); + linkCount_Y[3] = D3Q19_MapRecv(0,-1,-1,Dm->recvList("Y"),3*recvCount_Y,recvCount_Y,dvcRecvDist_Y,dvcRecvLinks_Y,Distance,Map); + linkCount_Y[4] = D3Q19_MapRecv(0,-1,1, Dm->recvList("Y"),4*recvCount_Y,recvCount_Y,dvcRecvDist_Y,dvcRecvLinks_Y,Distance,Map); + //................................................................................... + //...Map recieve list for the Y face: q=3,7,10,15,17 .................................. + linkCount_y[0] = D3Q19_MapRecv(0,1,0, Dm->recvList("y"),0,recvCount_y,dvcRecvDist_y,dvcRecvLinks_y,Distance,Map); + linkCount_y[1] = D3Q19_MapRecv(1,1,0, Dm->recvList("y"),recvCount_y,recvCount_y,dvcRecvDist_y,dvcRecvLinks_y,Distance,Map); + linkCount_y[2] = D3Q19_MapRecv(-1,1,0,Dm->recvList("y"),2*recvCount_y,recvCount_y,dvcRecvDist_y,dvcRecvLinks_y,Distance,Map); + linkCount_y[3] = D3Q19_MapRecv(0,1,1, Dm->recvList("y"),3*recvCount_y,recvCount_y,dvcRecvDist_y,dvcRecvLinks_y,Distance,Map); + linkCount_y[4] = D3Q19_MapRecv(0,1,-1,Dm->recvList("y"),4*recvCount_y,recvCount_y,dvcRecvDist_y,dvcRecvLinks_y,Distance,Map); + //................................................................................... + //...Map recieve list for the z face<<<6,12,13,16,17).............................................. + linkCount_Z[0] = D3Q19_MapRecv(0,0,-1, Dm->recvList("Z"),0,recvCount_Z,dvcRecvDist_Z,dvcRecvLinks_Z,Distance,Map); + linkCount_Z[1] = D3Q19_MapRecv(-1,0,-1,Dm->recvList("Z"),recvCount_Z,recvCount_Z,dvcRecvDist_Z,dvcRecvLinks_Z,Distance,Map); + linkCount_Z[2] = D3Q19_MapRecv(1,0,-1, Dm->recvList("Z"),2*recvCount_Z,recvCount_Z,dvcRecvDist_Z,dvcRecvLinks_Z,Distance,Map); + linkCount_Z[3] = D3Q19_MapRecv(0,-1,-1,Dm->recvList("Z"),3*recvCount_Z,recvCount_Z,dvcRecvDist_Z,dvcRecvLinks_Z,Distance,Map); + linkCount_Z[4] = D3Q19_MapRecv(0,1,-1, Dm->recvList("Z"),4*recvCount_Z,recvCount_Z,dvcRecvDist_Z,dvcRecvLinks_Z,Distance,Map); + //...Map recieve list for the Z face<<<5,11,14,15,18).............................................. + linkCount_z[0] = D3Q19_MapRecv(0,0,1, Dm->recvList("z"),0,recvCount_z,dvcRecvDist_z,dvcRecvLinks_z,Distance,Map); + linkCount_z[1] = D3Q19_MapRecv(1,0,1, Dm->recvList("z"),recvCount_z,recvCount_z,dvcRecvDist_z,dvcRecvLinks_z,Distance,Map); + linkCount_z[2] = D3Q19_MapRecv(-1,0,1,Dm->recvList("z"),2*recvCount_z,recvCount_z,dvcRecvDist_z,dvcRecvLinks_z,Distance,Map); + linkCount_z[3] = D3Q19_MapRecv(0,1,1, Dm->recvList("z"),3*recvCount_z,recvCount_z,dvcRecvDist_z,dvcRecvLinks_z,Distance,Map); + linkCount_z[4] = D3Q19_MapRecv(0,-1,1,Dm->recvList("z"),4*recvCount_z,recvCount_z,dvcRecvDist_z,dvcRecvLinks_z,Distance,Map); + //.................................................................................. + //...Map recieve list for the xy edge <<<8)................................ + linkCount_XY = D3Q19_MapRecv(-1,-1,0,Dm->recvList("XY"),0,recvCount_XY,dvcRecvDist_XY,dvcRecvLinks_XY,Distance,Map); + //...Map recieve list for the Xy edge <<<9)................................ + linkCount_xY = D3Q19_MapRecv(1,-1,0,Dm->recvList("xY"),0,recvCount_xY,dvcRecvDist_xY,dvcRecvLinks_xY,Distance,Map); + //...Map recieve list for the xY edge <<<10)................................ + linkCount_Xy = D3Q19_MapRecv(-1,1,0,Dm->recvList("Xy"),0,recvCount_Xy,dvcRecvDist_Xy,dvcRecvLinks_Xy,Distance,Map); + //...Map recieve list for the XY edge <<<7)................................ + linkCount_xy = D3Q19_MapRecv(1,1,0,Dm->recvList("xy"),0,recvCount_xy,dvcRecvDist_xy,dvcRecvLinks_xy,Distance,Map); + //...Map recieve list for the xz edge <<<12)................................ + linkCount_XZ = D3Q19_MapRecv(-1,0,-1,Dm->recvList("XZ"),0,recvCount_XZ,dvcRecvDist_XZ,dvcRecvLinks_XZ,Distance,Map); + //...Map recieve list for the xZ edge <<<14)................................ + linkCount_Xz = D3Q19_MapRecv(-1,0,1,Dm->recvList("Xz"),0,recvCount_Xz,dvcRecvDist_Xz,dvcRecvLinks_Xz,Distance,Map); + //...Map recieve list for the Xz edge <<<13)................................ + linkCount_xZ = D3Q19_MapRecv(1,0,-1,Dm->recvList("xZ"),0,recvCount_xZ,dvcRecvDist_xZ,dvcRecvLinks_xZ,Distance,Map); + //...Map recieve list for the XZ edge <<<11)................................ + linkCount_xz = D3Q19_MapRecv(1,0,1,Dm->recvList("xz"),0,recvCount_xz,dvcRecvDist_xz,dvcRecvLinks_xz,Distance,Map); + //...Map recieve list for the yz edge <<<16)................................ + linkCount_YZ = D3Q19_MapRecv(0,-1,-1,Dm->recvList("YZ"),0,recvCount_YZ,dvcRecvDist_YZ,dvcRecvLinks_YZ,Distance,Map); + //...Map recieve list for the yZ edge <<<18)................................ + linkCount_Yz = D3Q19_MapRecv(0,-1,1,Dm->recvList("Yz"),0,recvCount_Yz,dvcRecvDist_Yz,dvcRecvLinks_Yz,Distance,Map); + //...Map recieve list for the Yz edge <<<17)................................ + linkCount_yZ = D3Q19_MapRecv(0,1,-1,Dm->recvList("yZ"),0,recvCount_yZ,dvcRecvDist_yZ,dvcRecvLinks_yZ,Distance,Map); + //...Map recieve list for the YZ edge <<<15)................................ + linkCount_yz = D3Q19_MapRecv(0,1,1,Dm->recvList("yz"),0,recvCount_yz,dvcRecvDist_yz,dvcRecvLinks_yz,Distance,Map); + //................................................................................... + + //...................................................................................... + MPI_COMM_SCALBL.barrier(); + ScaLBL_DeviceBarrier(); + //....................................................................... + SendCount = sendCount_x+sendCount_X+sendCount_y+sendCount_Y+sendCount_z+sendCount_Z+ + sendCount_xy+sendCount_Xy+sendCount_xY+sendCount_XY+ + sendCount_xZ+sendCount_Xz+sendCount_xZ+sendCount_XZ+ + sendCount_yz+sendCount_Yz+sendCount_yZ+sendCount_YZ; + + RecvCount = recvCount_x+recvCount_X+recvCount_y+recvCount_Y+recvCount_z+recvCount_Z+ + recvCount_xy+recvCount_Xy+recvCount_xY+recvCount_XY+ + recvCount_xZ+recvCount_Xz+recvCount_xZ+recvCount_XZ+ + recvCount_yz+recvCount_Yz+recvCount_yZ+recvCount_YZ; + + CommunicationCount = SendCount+RecvCount; + //...................................................................................... + int *TempBuffer; + TempBuffer = new int [5*RecvCount]; + //....................................................................... + // Re-index the send lists + ScaLBL_CopyToHost(TempBuffer,dvcSendList_x,sendCount_x*sizeof(int)); + for (i=0; i db){ +void Membrane::AssignCoefficients(int *Map, double *Psi, string method){ + /* Assign mass transfer coefficients to the membrane data structure */ + + double Threshold; + double MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut; + + Threshold = -55.0; + MassFractionIn = 0.0; + MassFractionOut = 0.0; + ThresholdMassFractionOut = 0.0; + ThresholdMassFractionIn = 0.0; + + if (method == "Voltage Gated Potassium"){ + MassFractionIn = 0.0; + MassFractionOut = 0.0; + ThresholdMassFractionOut = 0.0; + ThresholdMassFractionIn = 1.0; + } + ScaLBL_D3Q7_Membrane_AssignLinkCoef(MembraneLinks, Map, MembraneDistance, Psi, MembraneCoef, + Threshold, MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut, + membraneLinkCount, Nx, Ny, Nz, Np); + + ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(-1,0,0,Map,MembraneDistance,Psi,Threshold, + MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut, + dvcRecvDist_X,dvcRecvLinks_X,coefficient_X,0,linkCount_X[0],recvCount_X, + Np,Nx,Ny,Nz); + + ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(1,0,0,Map,MembraneDistance,Psi,Threshold, + MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut, + dvcRecvDist_x,dvcRecvLinks_x,coefficient_x,0,linkCount_x[0],recvCount_x, + Np,Nx,Ny,Nz); + + ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,-1,0,Map,MembraneDistance,Psi,Threshold, + MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut, + dvcRecvDist_Y,dvcRecvLinks_Y,coefficient_Y,0,linkCount_Y[0],recvCount_Y, + Np,Nx,Ny,Nz); + + ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,1,0,Map,MembraneDistance,Psi,Threshold, + MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut, + dvcRecvDist_y,dvcRecvLinks_y,coefficient_y,0,linkCount_y[0],recvCount_y, + Np,Nx,Ny,Nz); + + ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,0,-1,Map,MembraneDistance,Psi,Threshold, + MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut, + dvcRecvDist_Z,dvcRecvLinks_Z,coefficient_Z,0,linkCount_Z[0],recvCount_Z, + Np,Nx,Ny,Nz); + + ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,0,1,Map,MembraneDistance,Psi,Threshold, + MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut, + dvcRecvDist_z,dvcRecvLinks_z,coefficient_z,0,linkCount_z[0],recvCount_z, + Np,Nx,Ny,Nz); + +} diff --git a/common/Membrane.h b/common/Membrane.h new file mode 100644 index 00000000..a153b8d6 --- /dev/null +++ b/common/Membrane.h @@ -0,0 +1,183 @@ +/* Flow adaptor class for multiphase flow methods */ + +#ifndef ScaLBL_Membrane_INC +#define ScaLBL_Membrane_INC +#include +#include +#include +#include +#include +#include +#include + +#include "common/ScaLBL.h" + +/** +* \brief Unpack D3Q19 distributions after communication using links determined based on membrane location +* @param q - index for distribution based on D3Q19 discrete velocity structure +* @param list - list of distributions to communicate +* @param links - list of active links based on the membrane location +* @param start - index to start parsing the list +* @param count - number of values to unppack +* @param recvbuf - memory buffer where recieved values have been stored +* @param dist - memory buffer to hold the distributions +* @param N - size of the distributions (derived from Domain structure) +*/ +extern "C" void Membrane_D3Q19_Unpack(int q, int *list, int *links, int start, int count, double *recvbuf, double *dist, int N); + + +/** +* \brief Set custom link rules for D3Q19 distribution based on membrane location +* @param q - index for distribution based on D3Q19 discrete velocity structure +* @param list - list of distributions to communicate +* @param links - list of active links based on the membrane location +* @param coef - coefficient to determine the local mass transport for each membrane link +* @param start - index to start parsing the list +* @param offset - offset to start reading membrane links +* @param count - number of values to unppack +* @param recvbuf - memory buffer where recieved values have been stored +* @param dist - memory buffer to hold the distributions +* @param N - size of the distributions (derived from Domain structure) +*/ +extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links, double *coef, int start, int offset, + int linkCount, double *recvbuf, double *dist, int N); + +/** + * \class Membrane + * @brief + * The Membrane class operates on ScaLBL data structures to insert membrane + * + */ + +class Membrane { +public: + int Np; + int Nx,Ny,Nz,N; + int membraneLinkCount; + + int *initialNeighborList; // original neighborlist + int *NeighborList; // modified neighborlist + + /* host data structures */ + int *membraneLinks; // D3Q19 links that cross membrane + int *membraneTag; // label each link in the membrane + double *membraneDist; // distance to membrane for each linked site + + /* + * Device data structures + */ + int *MembraneLinks; + double *MembraneCoef; // mass transport coefficient for the membrane + double *MembraneDistance; + + /** + * \brief Create a flow adaptor to operate on the LB model + * @param ScaLBL - originating data structures + * @param neighborList - list of neighbors for each site + */ + Membrane(std::shared_ptr Dm, int *initialNeighborList, int Nsites); + + /** + * \brief Destructor + */ + ~Membrane(); + + /** + * \brief Create membrane + * \details Create membrane structure from signed distance function + * @param Dm - domain structure + * @param Distance - signed distance to membrane + * @param Map - mapping between regular layout and compact layout + */ + int Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map); + + void SendD3Q19AA(double *dist); + void RecvD3Q19AA(double *dist); + void SendD3Q7AA(double *dist); + void RecvD3Q7AA(double *dist); + void AssignCoefficients(int *Map, double *Psi, std::string method); + //...................................................................................... + // Buffers to store data sent and recieved by this MPI process + double *sendbuf_x, *sendbuf_y, *sendbuf_z, *sendbuf_X, *sendbuf_Y, *sendbuf_Z; + double *sendbuf_xy, *sendbuf_yz, *sendbuf_xz, *sendbuf_Xy, *sendbuf_Yz, *sendbuf_xZ; + double *sendbuf_xY, *sendbuf_yZ, *sendbuf_Xz, *sendbuf_XY, *sendbuf_YZ, *sendbuf_XZ; + double *recvbuf_x, *recvbuf_y, *recvbuf_z, *recvbuf_X, *recvbuf_Y, *recvbuf_Z; + double *recvbuf_xy, *recvbuf_yz, *recvbuf_xz, *recvbuf_Xy, *recvbuf_Yz, *recvbuf_xZ; + double *recvbuf_xY, *recvbuf_yZ, *recvbuf_Xz, *recvbuf_XY, *recvbuf_YZ, *recvbuf_XZ; + //...................................................................................... + +private: + bool Lock; // use Lock to make sure only one call at a time to protect data in transit + int sendtag, recvtag; + int iproc,jproc,kproc; + int nprocx,nprocy,nprocz; + // Give the object it's own MPI communicator + RankInfoStruct rank_info; + Utilities::MPI MPI_COMM_SCALBL; // MPI Communicator for this domain + MPI_Request req1[18],req2[18]; + /** + * \brief Set up membrane communication + * \details associate p2p communication links to membrane where necessary + * returns the number of membrane links + * regular communications are stored in the first part of the list + * membrane communications are stored in the last part of the list + * @param Cqx - discrete velocity (x) + * @param Cqy - discrete velocity (y) + * @param Cqz - discrete velocity (z) + * @param list - list of recieved values + * @param count - number recieved values + * @param d3q19_recvlist - device array with the saved list + * @param d3q19_linkList - sorted list with regular and membrane links + * @param Distance - signed distance to membrane + * @param Map - data structure used to define mapping between dense and sparse representation + * */ + int D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list, int start, int count, + int *d3q19_recvlist, int *d3q19_linkList, DoubleArray &Distance, IntArray &Map); + //...................................................................................... + // MPI ranks for all 18 neighbors + //...................................................................................... + // These variables are all private to prevent external things from modifying them!! + //...................................................................................... + int rank; + int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z; + int rank_xy,rank_XY,rank_xY,rank_Xy; + int rank_xz,rank_XZ,rank_xZ,rank_Xz; + int rank_yz,rank_YZ,rank_yZ,rank_Yz; + //...................................................................................... + int SendCount, RecvCount, CommunicationCount; + //...................................................................................... + int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z; + int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ; + int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ; + //...................................................................................... + int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z; + int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ; + int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ; + //...................................................................................... + int linkCount_x[5], linkCount_y[5], linkCount_z[5], linkCount_X[5], linkCount_Y[5], linkCount_Z[5]; + int linkCount_xy, linkCount_yz, linkCount_xz, linkCount_Xy, linkCount_Yz, linkCount_xZ; + int linkCount_xY, linkCount_yZ, linkCount_Xz, linkCount_XY, linkCount_YZ, linkCount_XZ; + //...................................................................................... + // Send buffers that reside on the compute device + int *dvcSendList_x, *dvcSendList_y, *dvcSendList_z, *dvcSendList_X, *dvcSendList_Y, *dvcSendList_Z; + int *dvcSendList_xy, *dvcSendList_yz, *dvcSendList_xz, *dvcSendList_Xy, *dvcSendList_Yz, *dvcSendList_xZ; + int *dvcSendList_xY, *dvcSendList_yZ, *dvcSendList_Xz, *dvcSendList_XY, *dvcSendList_YZ, *dvcSendList_XZ; + // Recieve buffers that reside on the compute device + int *dvcRecvList_x, *dvcRecvList_y, *dvcRecvList_z, *dvcRecvList_X, *dvcRecvList_Y, *dvcRecvList_Z; + int *dvcRecvList_xy, *dvcRecvList_yz, *dvcRecvList_xz, *dvcRecvList_Xy, *dvcRecvList_Yz, *dvcRecvList_xZ; + int *dvcRecvList_xY, *dvcRecvList_yZ, *dvcRecvList_Xz, *dvcRecvList_XY, *dvcRecvList_YZ, *dvcRecvList_XZ; + // Link lists that reside on the compute device + int *dvcRecvLinks_x, *dvcRecvLinks_y, *dvcRecvLinks_z, *dvcRecvLinks_X, *dvcRecvLinks_Y, *dvcRecvLinks_Z; + int *dvcRecvLinks_xy, *dvcRecvLinks_yz, *dvcRecvLinks_xz, *dvcRecvLinks_Xy, *dvcRecvLinks_Yz, *dvcRecvLinks_xZ; + int *dvcRecvLinks_xY, *dvcRecvLinks_yZ, *dvcRecvLinks_Xz, *dvcRecvLinks_XY, *dvcRecvLinks_YZ, *dvcRecvLinks_XZ; + // Recieve buffers for the distributions + int *dvcRecvDist_x, *dvcRecvDist_y, *dvcRecvDist_z, *dvcRecvDist_X, *dvcRecvDist_Y, *dvcRecvDist_Z; + int *dvcRecvDist_xy, *dvcRecvDist_yz, *dvcRecvDist_xz, *dvcRecvDist_Xy, *dvcRecvDist_Yz, *dvcRecvDist_xZ; + int *dvcRecvDist_xY, *dvcRecvDist_yZ, *dvcRecvDist_Xz, *dvcRecvDist_XY, *dvcRecvDist_YZ, *dvcRecvDist_XZ; + //...................................................................................... + // mass transfer coefficient arrays + double *coefficient_x, *coefficient_X, *coefficient_y, *coefficient_Y, *coefficient_z, *coefficient_Z; + //...................................................................................... + +}; +#endif \ No newline at end of file diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index b1aec304..dc455c38 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -1,5 +1,5 @@ /* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University + Copyright 2013--2022 James E. McClure, Virginia Polytechnic & State University Copyright Equnior ASA This file is part of the Open Porous Media project (OPM). @@ -15,10 +15,8 @@ along with OPM. If not, see . */ #include "common/ScaLBL.h" - #include - ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr Dm){ //...................................................................................... Lock=false; // unlock the communicator @@ -366,7 +364,6 @@ ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr Dm){ } - ScaLBL_Communicator::~ScaLBL_Communicator() { @@ -888,8 +885,6 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis idx=Map(n); //if (rank == 0) printf("r: mapped n=%d\n",idx); TempBuffer[i]=idx; - - } ScaLBL_CopyToDevice(dvcRecvDist_x,TempBuffer,5*recvCount_x*sizeof(int)); @@ -1046,7 +1041,6 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis return(Np); } - void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC) { @@ -1448,7 +1442,6 @@ void ScaLBL_Communicator::SolidSlippingVelocityBCD3Q19(double *fq, double *zeta_ void ScaLBL_Communicator::SendD3Q19AA(double *dist){ - // 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?"); } diff --git a/common/ScaLBL.h b/common/ScaLBL.h index c2413dad..3fd0942b 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -217,6 +217,25 @@ extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int */ extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz); +// MEMBRANE MODEL + +extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, double *dist, double *Den, int memLinks, int Np); + +extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef, + double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, + int memLinks, int Nx, int Ny, int Nz, int Np); + +extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo( + const int Cqx, const int Cqy, int const Cqz, + int *Map, double *Distance, double *Psi, double Threshold, + double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, + int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count, + const int N, const int Nx, const int Ny, const int Nz); + +extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q, + int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count, + double *recvbuf, double *dist, int N, double *coef); + // GREYSCALE MODEL (Single-component) extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *Dist, int Np, double Den); @@ -702,6 +721,14 @@ public: double GetPerformance(int *NeighborList, double *fq, int Np); int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np, int width); + /** + * \brief Create membrane data structure + * - cut lattice links based on distance map + * @param Distance - signed distance to membrane + * @param neighborList - data structure that retains lattice links + * @param Np - number of lattice sites + * @param width - halo width for the model + */ void Barrier(){ ScaLBL_DeviceBarrier(); MPI_COMM_SCALBL.barrier(); @@ -782,7 +809,6 @@ private: int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ; int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ; //...................................................................................... - int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z; int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ; int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ; diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp index bc2dd70f..a10b631c 100644 --- a/cpu/D3Q19.cpp +++ b/cpu/D3Q19.cpp @@ -48,6 +48,7 @@ extern "C" void ScaLBL_D3Q19_Unpack(int q, int *list, int start, int count, } } + extern "C" void ScaLBL_D3Q19_AA_Init(double *f_even, double *f_odd, int Np) { int n; for (n = 0; n < Np; n++) { diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index eb5e3ef5..915edda3 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -1,4 +1,156 @@ #include +#include + +extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef, + double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, + int memLinks, int Nx, int Ny, int Nz, int Np){ + + int link,iq,ip,nq,np,nqm,npm; + double aq, ap, membranePotential; + /* Interior Links */ + for (link=0; link Threshold){ + aq = ThresholdMassFractionIn; ap = ThresholdMassFractionOut; + } + + /* Save the mass transfer coefficients */ + coef[2*link] = aq; coef[2*link+1] = ap; + } +} + +extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo( + const int Cqx, const int Cqy, int const Cqz, + int *Map, double *Distance, double *Psi, double Threshold, + double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, + int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count, + const int N, const int Nx, const int Ny, const int Nz) { + //.................................................................................... + // Unack distribution from the recv buffer + // Distribution q matche Cqx, Cqy, Cqz + // swap rule means that the distributions in recvbuf are OPPOSITE of q + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int n, idx, link, nqm, npm, i, j, k; + double distanceLocal, distanceNonlocal; + double psiLocal, psiNonlocal, membranePotential; + double ap,aq; // coefficient + + /* second enforce custom rule for membrane links */ + for (link = nlinks; link < count; link++) { + // get the index for the recv list (deal with reordering of links) + idx = d3q7_linkList[link]; // THINK start NEEDS TO BE HERE + // get the distribution index + n = d3q7_recvlist[start+idx]; + // get the index in strided layout + nqm = Map[n]; + distanceLocal = Distance[nqm]; + psiLocal = Psi[nqm]; + + // Get the 3-D indices from the send process + k = nqm/(Nx*Ny); j = (nqm-Nx*Ny*k)/Nx; i = nqm-Nx*Ny*k-Nx*j; + // Streaming link the non-local distribution + i -= Cqx; j -= Cqy; k -= Cqz; + npm = k*Nx*Ny + j*Nx + i; + distanceNonlocal = Distance[npm]; + psiNonlocal = Psi[npm]; + + membranePotential = psiLocal - psiNonlocal; + aq = MassFractionIn; + ap = MassFractionOut; + + /* link is inside membrane */ + if (distanceLocal > 0.0){ + if (membranePotential < Threshold*(-1.0)){ + ap = MassFractionIn; + aq = MassFractionOut; + } + else { + ap = ThresholdMassFractionIn; + aq = ThresholdMassFractionOut; + } + } + else if (membranePotential > Threshold){ + aq = ThresholdMassFractionIn; + ap = ThresholdMassFractionOut; + } + + // update link based on mass transfer coefficients + coef[2*(link-nlinks)] = aq; + coef[2*(link-nlinks)+1] = ap; + } +} + + +extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q, + int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count, + double *recvbuf, double *dist, int N, double *coef) { + //.................................................................................... + // Unack distribution from the recv buffer + // Distribution q matche Cqx, Cqy, Cqz + // swap rule means that the distributions in recvbuf are OPPOSITE of q + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int n, idx, link; + double fq,fp,fqq,ap,aq; // coefficient + /* First unpack the regular links */ + for (link = 0; link < nlinks; link++) { + // get the index for the recv list (deal with reordering of links) + idx = d3q7_linkList[link]; + // get the distribution index + n = d3q7_recvlist[start+idx]; + if (!(n < 0)){ + fp = recvbuf[start + idx]; + dist[q * N + n] = fp; + } + //printf(" site=%i, index=%i, value = %e \n",n,idx,fp); + } + /* second enforce custom rule for membrane links */ + for (link = nlinks; link < count; link++) { + // get the index for the recv list (deal with reordering of links) + idx = d3q7_linkList[link]; + // get the distribution index + n = d3q7_recvlist[start+idx]; + // update link based on mass transfer coefficients + if (!(n < 0)){ + aq = coef[2*(link-nlinks)]; + ap = coef[2*(link-nlinks)+1]; + fq = dist[q * N + n]; + fp = recvbuf[start + idx]; + fqq = (1-aq)*fq+ap*fp; + dist[q * N + n] = fqq; + } + //printf(" LINK: site=%i, index=%i \n", n, idx); + + } +} + +extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, + double *dist, double *Den, int memLinks, int Np){ + + int link,iq,ip,nq,np; + double aq, ap, fq, fp, fqq, fpp, Cq, Cp; + for (link=0; link Threshold){ + aq = ThresholdMassFractionIn; ap = ThresholdMassFractionOut; + } + + /* Save the mass transfer coefficients */ + coef[2*link] = aq; coef[2*link+1] = ap; + } + } +} + +__global__ void dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo( + const int Cqx, const int Cqy, int const Cqz, + int *Map, double *Distance, double *Psi, double Threshold, + double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, + int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count, + const int N, const int Nx, const int Ny, const int Nz) { + //.................................................................................... + // Unack distribution from the recv buffer + // Distribution q matche Cqx, Cqy, Cqz + // swap rule means that the distributions in recvbuf are OPPOSITE of q + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int n, idx, link, nqm, npm, i, j, k; + double distanceLocal, distanceNonlocal; + double psiLocal, psiNonlocal, membranePotential; + double ap,aq; // coefficient + + /* second enforce custom rule for membrane links */ + int S = (count-nlinks)/NBLOCKS/NTHREADS + 1; + for (int s=0; s 0.0){ + if (membranePotential < Threshold*(-1.0)){ + ap = MassFractionIn; + aq = MassFractionOut; + } + else { + ap = ThresholdMassFractionIn; + aq = ThresholdMassFractionOut; + } + } + else if (membranePotential > Threshold){ + aq = ThresholdMassFractionIn; + ap = ThresholdMassFractionOut; + } + + // update link based on mass transfer coefficients + coef[2*(link-nlinks)] = aq; + coef[2*(link-nlinks)+1] = ap; + } + } +} + +__global__ void dvc_ScaLBL_D3Q7_Membrane_Unpack(int q, + int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count, + double *recvbuf, double *dist, int N, double *coef) { + //.................................................................................... + // Unack distribution from the recv buffer + // Distribution q matche Cqx, Cqy, Cqz + // swap rule means that the distributions in recvbuf are OPPOSITE of q + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int n, idx, link; + double fq,fp,fqq,ap,aq; // coefficient + + /* second enforce custom rule for membrane links */ + int S = count/NBLOCKS/NTHREADS + 1; + for (int s=0; s>>(membrane, Map, Distance, Psi, coef, + Threshold, MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut, + memLinks, Nx, Ny, Nz, Np); + + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef: %s \n",cudaGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo( + const int Cqx, const int Cqy, int const Cqz, + int *Map, double *Distance, double *Psi, double Threshold, + double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, + int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count, + const int N, const int Nx, const int Ny, const int Nz) { + + dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo<<>>( + Cqx, Cqy, Cqz, Map, Distance, Psi, Threshold, + MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut, + d3q7_recvlist, d3q7_linkList, coef, start, nlinks, count, N, Nx, Ny, Nz); + + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo: %s \n",cudaGetErrorString(err)); + } +} + + +extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q, + int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count, + double *recvbuf, double *dist, int N, double *coef) { + + dvc_ScaLBL_D3Q7_Membrane_Unpack<<>>(q, d3q7_recvlist, d3q7_linkList, start, nlinks, count, + recvbuf, dist, N, coef) ; + + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_Unpack: %s \n",cudaGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, + double *dist, double *Den, int memLinks, int Np){ + + dvc_ScaLBL_D3Q7_Membrane_IonTransport<<>>(membrane, coef, dist, Den, memLinks, Np); + + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_IonTransport: %s \n",cudaGetErrorString(err)); + } +} + diff --git a/example/Bubble/CreateBubble.py b/example/Bubble/CreateBubble.py new file mode 100644 index 00000000..76fbefab --- /dev/null +++ b/example/Bubble/CreateBubble.py @@ -0,0 +1,18 @@ +import numpy as np +import matplotlib.pylab as plt + +D=np.ones((40,40,40),dtype="uint8") + +cx = 20 +cy = 20 +cz = 20 + +for i in range(0, 40): + for j in range (0, 40): + for k in range (0,40): + dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz)) + if (dist < 12.5 ) : + D[i,j,k] = 2 + +D.tofile("bubble_40x40x40.raw") + diff --git a/example/Bubble/CreateCell.py b/example/Bubble/CreateCell.py new file mode 100644 index 00000000..ac62863f --- /dev/null +++ b/example/Bubble/CreateCell.py @@ -0,0 +1,77 @@ +import numpy as np +import matplotlib.pylab as plt + +D=np.ones((40,40,40),dtype="uint8") + +cx = 20 +cy = 20 +cz = 20 + +for i in range(0, 40): + for j in range (0, 40): + for k in range (0,40): + dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz)) + if (dist < 15.5 ) : + D[i,j,k] = 2 + +D.tofile("cell_40x40x40.raw") + + +C1=np.zeros((40,40,40),dtype="double") +C2=np.zeros((40,40,40),dtype="double") +C3=np.zeros((40,40,40),dtype="double") +C4=np.zeros((40,40,40),dtype="double") +C5=np.zeros((40,40,40),dtype="double") +C6=np.zeros((40,40,40),dtype="double") + +for i in range(0, 40): + for j in range (0, 40): + for k in range (0,40): + #outside the cell + C1[i,j,k] = 4.0e-6 # K + C2[i,j,k] = 150.0e-6 # Na + C3[i,j,k] = 116.0e-6 # Cl + C4[i,j,k] = 29.0e-6 # HC03 + #C5[i,j,k] = 2.4e-6 # Ca + dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz)) + # inside the cell + if (dist < 15.5 ) : + C1[i,j,k] = 145.0e-6 + C2[i,j,k] = 12.0e-6 + C3[i,j,k] = 4.0e-6 + C4[i,j,k] = 12.0e-6 # 12 mmol / L + #C5[i,j,k] = 0.10e-6 # 100 nmol / L + + +# add up the total charge to make sure it is zero +TotalCharge = 0 +for i in range(0, 40): + for j in range (0, 40): + for k in range (0,40): + TotalCharge += C1[i,j,k] + C2[i,j,k] - C3[i,j,k] - C4[i,j,k] + +TotalCharge /= (40*40*40) + +print("Total charge " + str(TotalCharge)) + + +for i in range(0, 40): + for j in range (0, 40): + for k in range (0,40): + if TotalCharge < 0 : + # need more cation + C5[i,j,k] = abs(TotalCharge) + C6[i,j,k] = 0.0 + else : + # need more anion + C5[i,j,k] = 0.0 + C6[i,j,k] = abs(TotalCharge) + + +C1.tofile("cell_concentration_K_40x40x40.raw") +C2.tofile("cell_concentration_Na_40x40x40.raw") +C3.tofile("cell_concentration_Cl_40x40x40.raw") +C4.tofile("cell_concentration_HCO3_40x40x40.raw") +C5.tofile("cell_concentration_cation_40x40x40.raw") +C6.tofile("cell_concentration_anion_40x40x40.raw") + diff --git a/example/Bubble/cell.db b/example/Bubble/cell.db new file mode 100644 index 00000000..ccbf37e6 --- /dev/null +++ b/example/Bubble/cell.db @@ -0,0 +1,75 @@ +MultiphysController { + timestepMax = 60 + num_iter_Ion_List = 2 + analysis_interval = 50 + tolerance = 1.0e-9 + visualization_interval = 100 // Frequency to write visualization data + analysis_interval = 50 // Frequency to perform analysis +} +Stokes { + tau = 1.0 + F = 0, 0, 0 + ElectricField = 0, 0, 0 //body electric field; user-input unit: [V/m] + nu_phys = 0.889e-6 //fluid kinematic viscosity; user-input unit: [m^2/sec] +} +Ions { + IonConcentrationFile = "cell_concentration_K_40x40x40.raw", "double", "cell_concentration_Na_40x40x40.raw", "double", "cell_concentration_Cl_40x40x40.raw", "double", "cell_concentration_HCO3_40x40x40.raw", "double", "cell_concentration_anion_40x40x40.raw", "double", "cell_concentration_cation_40x40x40.raw", "double" + temperature = 293.15 //unit [K] + number_ion_species = 6 //number of ions + tauList = 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 + IonDiffusivityList = 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec] + IonValenceList = 1, 1, -1, -1, 1, -1 //valence charge of ions; dimensionless; positive/negative integer + IonConcentrationList = 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6 //user-input unit: [mol/m^3] + BC_Solid = 0 //solid boundary condition; 0=non-flux BC; 1=surface ion concentration + //SolidLabels = 0 //solid labels for assigning solid boundary condition; ONLY for BC_Solid=1 + //SolidValues = 1.0e-5 // user-input surface ion concentration unit: [mol/m^2]; ONLY for BC_Solid=1 + FluidVelDummy = 0.0, 0.0, 1.0e-2 // dummy fluid velocity for debugging +} +Poisson { + epsilonR = 78.5 //fluid dielectric constant [dimensionless] + BC_Inlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential + BC_Outlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential + //-------------------------------------------------------------------------- + //-------------------------------------------------------------------------- + BC_Solid = 2 //solid boundary condition; 1=surface potential; 2=surface charge density + SolidLabels = 0 //solid labels for assigning solid boundary condition + SolidValues = 0 //if surface potential, unit=[V]; if surface charge density, unit=[C/m^2] + WriteLog = true //write convergence log for LB-Poisson solver + // ------------------------------- Testing Utilities ---------------------------------------- + // ONLY for code debugging; the followings test sine/cosine voltage BCs; disabled by default + TestPeriodic = false + TestPeriodicTime = 1.0 //unit:[sec] + TestPeriodicTimeConv = 0.01 //unit:[sec] + TestPeriodicSaveInterval = 0.2 //unit:[sec] + //------------------------------ advanced setting ------------------------------------ + timestepMax = 100000 //max timestep for obtaining steady-state electrical potential + analysis_interval = 200 //timestep checking steady-state convergence + tolerance = 1.0e-6 //stopping criterion for steady-state solution +} +Domain { + Filename = "cell_40x40x40.raw" + nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz) + n = 40, 40, 40 // Size of local domain (Nx,Ny,Nz) + N = 40, 40, 40 // size of the input image + voxel_length = 1.0 //resolution; user-input unit: [um] + BC = 0 // Boundary condition type + ReadType = "8bit" + ReadValues = 0, 1, 2 + WriteValues = 0, 1, 2 +} +Analysis { + analysis_interval = 100 + subphase_analysis_interval = 50 // Frequency to perform analysis + restart_interval = 5000 // Frequency to write restart data + restart_file = "Restart" // Filename to use for restart file (will append rank) + N_threads = 4 // Number of threads to use + load_balance = "independent" // Load balance method to use: "none", "default", "independent" +} +Visualization { + save_electric_potential = true + save_concentration = true + save_velocity = true +} +Membrane { + MembraneLabels = 2 +} diff --git a/hip/Ion.hip b/hip/Ion.hip index b1d9636e..133e2da8 100644 --- a/hip/Ion.hip +++ b/hip/Ion.hip @@ -5,6 +5,179 @@ #define NBLOCKS 1024 #define NTHREADS 256 +__global__ void dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef, + double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, + int memLinks, int Nx, int Ny, int Nz, int Np){ + + int link,iq,ip,nq,np,nqm,npm; + double aq, ap, membranePotential; + /* Interior Links */ + + int S = memLinks/NBLOCKS/NTHREADS + 1; + for (int s=0; s Threshold){ + aq = ThresholdMassFractionIn; ap = ThresholdMassFractionOut; + } + + /* Save the mass transfer coefficients */ + coef[2*link] = aq; coef[2*link+1] = ap; + } + } +} + +__global__ void dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo( + const int Cqx, const int Cqy, int const Cqz, + int *Map, double *Distance, double *Psi, double Threshold, + double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, + int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count, + const int N, const int Nx, const int Ny, const int Nz) { + //.................................................................................... + // Unack distribution from the recv buffer + // Distribution q matche Cqx, Cqy, Cqz + // swap rule means that the distributions in recvbuf are OPPOSITE of q + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int n, idx, link, nqm, npm, i, j, k; + double distanceLocal, distanceNonlocal; + double psiLocal, psiNonlocal, membranePotential; + double ap,aq; // coefficient + + /* second enforce custom rule for membrane links */ + int S = (count-nlinks)/NBLOCKS/NTHREADS + 1; + for (int s=0; s 0.0){ + if (membranePotential < Threshold*(-1.0)){ + ap = MassFractionIn; + aq = MassFractionOut; + } + else { + ap = ThresholdMassFractionIn; + aq = ThresholdMassFractionOut; + } + } + else if (membranePotential > Threshold){ + aq = ThresholdMassFractionIn; + ap = ThresholdMassFractionOut; + } + + // update link based on mass transfer coefficients + coef[2*(link-nlinks)] = aq; + coef[2*(link-nlinks)+1] = ap; + } + } +} + +__global__ void dvc_ScaLBL_D3Q7_Membrane_Unpack(int q, + int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count, + double *recvbuf, double *dist, int N, double *coef) { + //.................................................................................... + // Unack distribution from the recv buffer + // Distribution q matche Cqx, Cqy, Cqz + // swap rule means that the distributions in recvbuf are OPPOSITE of q + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int n, idx, link; + double fq,fp,fqq,ap,aq; // coefficient + + /* second enforce custom rule for membrane links */ + int S = count/NBLOCKS/NTHREADS + 1; + for (int s=0; s>>(membrane, Map, Distance, Psi, coef, + Threshold, MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut, + memLinks, Nx, Ny, Nz, Np); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo( + const int Cqx, const int Cqy, int const Cqz, + int *Map, double *Distance, double *Psi, double Threshold, + double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, + int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count, + const int N, const int Nx, const int Ny, const int Nz) { + + dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo<<>>( + Cqx, Cqy, Cqz, Map, Distance, Psi, Threshold, + MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut, + d3q7_recvlist, d3q7_linkList, coef, start, nlinks, count, N, Nx, Ny, Nz); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo: %s \n",hipGetErrorString(err)); + } +} + + +extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q, + int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count, + double *recvbuf, double *dist, int N, double *coef) { + + dvc_ScaLBL_D3Q7_Membrane_Unpack<<>>(q, d3q7_recvlist, d3q7_linkList, start, nlinks, count, + recvbuf, dist, N, coef) ; + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_Unpack: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, + double *dist, double *Den, int memLinks, int Np){ + + dvc_ScaLBL_D3Q7_Membrane_IonTransport<<>>(membrane, coef, dist, Den, memLinks, Np); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_IonTransport: %s \n",hipGetErrorString(err)); + } +} diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index a09fcd6e..8cfed4bb 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -387,6 +387,10 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase) { AFFINITY, volume_fraction); } } + + // clean up + delete [] label_count; + delete [] label_count_global; } void ScaLBL_ColorModel::Create() { diff --git a/models/IonModel.cpp b/models/IonModel.cpp index dc08bd1d..db9ca039 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -583,6 +583,72 @@ void ScaLBL_IonModel::SetDomain() { nprocz = Dm->nprocz(); } +void ScaLBL_IonModel::SetMembrane() { + size_t NLABELS = 0; + + membrane_db = db->getDatabase("Membrane"); + + /* set distance based on labels inside the membrane (all other labels will be outside) */ + auto MembraneLabels = membrane_db->getVector("MembraneLabels"); + + IonMembrane = std::shared_ptr(new Membrane(Dm, NeighborList, Np)); + + signed char LABEL = 0; + double *label_count; + double *label_count_global; + Array membrane_id(Nx,Ny,Nz); + label_count = new double[NLABELS]; + label_count_global = new double[NLABELS]; + // Assign the labels + for (size_t idx = 0; idx < NLABELS; idx++) + label_count[idx] = 0; + /* set the distance to the membrane */ + MembraneDistance.resize(Nx, Ny, Nz); + MembraneDistance.fill(0); + for (int k = 0; k < Nz; k++) { + for (int j = 0; j < Ny; j++) { + for (int i = 0; i < Nx; i++) { + membrane_id(i,j,k) = 1; // default value + LABEL = Dm->id[k*Nx*Ny + j*Nx + i]; + for (size_t m=0; mComm.sumReduce(label_count[m]); + } + if (rank == 0) { + printf(" Membrane labels: %lu \n", MembraneLabels.size()); + for (size_t m=0; mCreate(Dm, MembraneDistance, Map); + + // clean up + delete [] label_count; + delete [] label_count_global; +} + void ScaLBL_IonModel::ReadInput() { sprintf(LocalRankString, "%05d", Dm->rank()); @@ -778,6 +844,7 @@ void ScaLBL_IonModel::Create() { int neighborSize = 18 * (Np * sizeof(int)); //........................................................................... ScaLBL_AllocateDeviceMemory((void **)&NeighborList, neighborSize); + ScaLBL_AllocateDeviceMemory((void **)&dvcMap, sizeof(int) * Np); ScaLBL_AllocateDeviceMemory((void **)&fq, number_ion_species * 7 * dist_mem_size); ScaLBL_AllocateDeviceMemory((void **)&Ci, @@ -794,6 +861,37 @@ void ScaLBL_IonModel::Create() { if (rank == 0) printf("LB Ion Solver: Setting up device map and neighbor list \n"); // copy the neighbor list + int *TmpMap; + TmpMap = new int[Np]; + for (int k = 1; k < Nz - 1; k++) { + for (int j = 1; j < Ny - 1; j++) { + for (int i = 1; i < Nx - 1; i++) { + int idx = Map(i, j, k); + if (!(idx < 0)) + TmpMap[idx] = k * Nx * Ny + j * Nx + i; + } + } + } + // check that TmpMap is valid + for (int idx = 0; idx < ScaLBL_Comm->LastExterior(); idx++) { + auto n = TmpMap[idx]; + if (n > Nx * Ny * Nz) { + printf("Bad value! idx=%i \n", n); + TmpMap[idx] = Nx * Ny * Nz - 1; + } + } + for (int idx = ScaLBL_Comm->FirstInterior(); + idx < ScaLBL_Comm->LastInterior(); idx++) { + auto n = TmpMap[idx]; + if (n > Nx * Ny * Nz) { + printf("Bad value! idx=%i \n", n); + TmpMap[idx] = Nx * Ny * Nz - 1; + } + } + ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int) * Np); + ScaLBL_Comm->Barrier(); + delete[] TmpMap; + ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); comm.barrier(); @@ -825,7 +923,7 @@ void ScaLBL_IonModel::Initialize() { if (ion_db->keyExists("IonConcentrationFile")) { //NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype" auto File_ion = ion_db->getVector("IonConcentrationFile"); - if (File_ion.size() == 2 * number_ion_species) { + if (File_ion.size() == 2*number_ion_species) { double *Ci_host; Ci_host = new double[number_ion_species * Np]; for (size_t ic = 0; ic < number_ion_species; ic++) { @@ -1181,6 +1279,225 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField) { //if (rank==0) printf("********************************************************\n"); } +void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, double *Psi) { + + //Input parameter: + //1. Velocity is from StokesModel + //2. ElectricField is from Poisson model + + //LB-related parameter + vector rlx; + for (size_t ic = 0; ic < tau.size(); ic++) { + rlx.push_back(1.0 / tau[ic]); + } + + //.......create and start timer............ + //double starttime,stoptime,cputime; + //ScaLBL_Comm->Barrier(); comm.barrier(); + //auto t1 = std::chrono::system_clock::now(); + /* set the mass transfer coefficients for the membrane */ + IonMembrane->AssignCoefficients(dvcMap, Psi, "default"); + + for (size_t ic = 0; ic < number_ion_species; ic++) { + timestep = 0; + while (timestep < timestepMax[ic]) { + //************************************************************************/ + // *************ODD TIMESTEP*************// + timestep++; + //Update ion concentration and charge density + IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL + ScaLBL_D3Q7_AAodd_IonConcentration( + IonMembrane->NeighborList, &fq[ic * Np * 7], &Ci[ic * Np], + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE + /* ScaLBL_Comm->Barrier(); + //--------------------------------------- Set boundary conditions -------------------------------------// + if (BoundaryConditionInlet[ic] > 0) { + switch (BoundaryConditionInlet[ic]) { + case 1: + ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], timestep); + break; + case 21: + ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 22: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 23: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], &ElectricField[2 * Np], + IonDiffusivity[ic], IonValence[ic], Vt, timestep); + break; + } + } + if (BoundaryConditionOutlet[ic] > 0) { + switch (BoundaryConditionOutlet[ic]) { + case 1: + ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], timestep); + break; + case 21: + ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 22: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 23: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], &ElectricField[2 * Np], + IonDiffusivity[ic], IonValence[ic], Vt, timestep); + break; + } + } + */ + //----------------------------------------------------------------------------------------------------// + ScaLBL_D3Q7_AAodd_IonConcentration(IonMembrane->NeighborList, &fq[ic * Np * 7], + &Ci[ic * Np], 0, + ScaLBL_Comm->LastExterior(), Np); + + //LB-Ion collison + ScaLBL_D3Q7_AAodd_Ion( + IonMembrane->NeighborList, &fq[ic * Np * 7], &Ci[ic * Np], + &FluxDiffusive[3 * ic * Np], &FluxAdvective[3 * ic * Np], + &FluxElectrical[3 * ic * Np], Velocity, ElectricField, + IonDiffusivity[ic], IonValence[ic], rlx[ic], Vt, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q7_AAodd_Ion( + IonMembrane->NeighborList, &fq[ic * Np * 7], &Ci[ic * Np], + &FluxDiffusive[3 * ic * Np], &FluxAdvective[3 * ic * Np], + &FluxElectrical[3 * ic * Np], Velocity, ElectricField, + IonDiffusivity[ic], IonValence[ic], rlx[ic], Vt, 0, + ScaLBL_Comm->LastExterior(), Np); + + if (BoundaryConditionSolid == 1) { + //TODO IonSolid may also be species-dependent + ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic * Np * 7], IonSolid); + } + ScaLBL_Comm->Barrier(); + comm.barrier(); + + // *************EVEN TIMESTEP*************// + timestep++; + //Update ion concentration and charge density + IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL + ScaLBL_D3Q7_AAeven_IonConcentration( + &fq[ic * Np * 7], &Ci[ic * Np], ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + //--------------------------------------- Set boundary conditions -------------------------------------// + /*if (BoundaryConditionInlet[ic] > 0) { + switch (BoundaryConditionInlet[ic]) { + case 1: + ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], timestep); + break; + case 21: + ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],q + &Velocity[2 * Np], timestep); + break; + case 22: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 23: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], &ElectricField[2 * Np], + IonDiffusivity[ic], IonValence[ic], Vt, timestep); + break; + } + } + if (BoundaryConditionOutlet[ic] > 0) { + switch (BoundaryConditionOutlet[ic]) { + case 1: + ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], timestep); + break; + case 21: + ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 22: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 23: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z( + IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], &ElectricField[2 * Np], + IonDiffusivity[ic], IonValence[ic], Vt, timestep); + break; + } + } + */ + //----------------------------------------------------------------------------------------------------// + ScaLBL_D3Q7_AAeven_IonConcentration(&fq[ic * Np * 7], &Ci[ic * Np], + 0, ScaLBL_Comm->LastExterior(), + Np); + + //LB-Ion collison + ScaLBL_D3Q7_AAeven_Ion( + &fq[ic * Np * 7], &Ci[ic * Np], &FluxDiffusive[3 * ic * Np], + &FluxAdvective[3 * ic * Np], &FluxElectrical[3 * ic * Np], + Velocity, ElectricField, IonDiffusivity[ic], IonValence[ic], + rlx[ic], Vt, ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q7_AAeven_Ion( + &fq[ic * Np * 7], &Ci[ic * Np], &FluxDiffusive[3 * ic * Np], + &FluxAdvective[3 * ic * Np], &FluxElectrical[3 * ic * Np], + Velocity, ElectricField, IonDiffusivity[ic], IonValence[ic], + rlx[ic], Vt, 0, ScaLBL_Comm->LastExterior(), Np); + + if (BoundaryConditionSolid == 1) { + //TODO IonSolid may also be species-dependent + ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic * Np * 7], IonSolid); + } + ScaLBL_Comm->Barrier(); + comm.barrier(); + } + } + + //Compute charge density for Poisson equation + for (size_t ic = 0; ic < number_ion_species; ic++) { + ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, + ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, + ScaLBL_Comm->LastExterior(), Np); + } + //************************************************************************/ + //if (rank==0) printf("-------------------------------------------------------------------\n"); + //// Compute the walltime per timestep + //auto t2 = std::chrono::system_clock::now(); + //double cputime = std::chrono::duration( t2 - t1 ).count() / timestep; + //// Performance obtained from each node + //double MLUPS = double(Np)/cputime/1000000; + + //if (rank==0) printf("********************************************************\n"); + //if (rank==0) printf("CPU time = %f \n", cputime); + //if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); + //MLUPS *= nprocs; + //if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); + //if (rank==0) printf("********************************************************\n"); +} + + void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const size_t ic) { //This function wirte out the data in a normal layout (by aggregating all decomposed domains) diff --git a/models/IonModel.h b/models/IonModel.h index 8930b1df..5b899c6c 100644 --- a/models/IonModel.h +++ b/models/IonModel.h @@ -1,7 +1,6 @@ /* * Ion transporte LB Model */ - #ifndef ScaLBL_IonModel_INC #define ScaLBL_IonModel_INC @@ -16,6 +15,7 @@ #include "common/ScaLBL.h" #include "common/Communication.h" +#include "common/Membrane.h" #include "common/MPI.h" #include "analysis/Minkowski.h" #include "ProfilerApp.h" @@ -30,10 +30,12 @@ public: void ReadParams(string filename); void ReadParams(std::shared_ptr db0); void SetDomain(); + void SetMembrane(); void ReadInput(); void Create(); void Initialize(); void Run(double *Velocity, double *ElectricField); + void RunMembrane(double *Velocity, double *ElectricField, double *Psi); void getIonConcentration(DoubleArray &IonConcentration, const size_t ic); void getIonConcentration_debug(int timestep); void getIonFluxDiffusive(DoubleArray &IonFlux_x, DoubleArray &IonFlux_y, @@ -88,6 +90,7 @@ public: IntArray Map; DoubleArray Distance; int *NeighborList; + int *dvcMap; double *fq; double *Ci; double *ChargeDensity; @@ -97,6 +100,12 @@ public: double *FluxDiffusive; double *FluxAdvective; double *FluxElectrical; + + /* these support membrane capabilities */ + std::shared_ptr membrane_db; + std::shared_ptr IonMembrane; + DoubleArray MembraneDistance; + int MembraneCount; // number of links the cross the membrane private: Utilities::MPI comm; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9e9dcf0a..f5d7943f 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -6,6 +6,7 @@ ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator ) ADD_LBPM_EXECUTABLE( lbpm_greyscale_simulator ) ADD_LBPM_EXECUTABLE( lbpm_greyscaleColor_simulator ) ADD_LBPM_EXECUTABLE( lbpm_electrokinetic_SingleFluid_simulator ) +ADD_LBPM_EXECUTABLE( lbpm_cell_simulator ) ADD_LBPM_EXECUTABLE( lbpm_freelee_simulator ) ADD_LBPM_EXECUTABLE( lbpm_freelee_SingleFluidBGK_simulator ) ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator ) @@ -58,6 +59,7 @@ ADD_LBPM_TEST( TestTopo3D ) ADD_LBPM_TEST( TestFluxBC ) ADD_LBPM_TEST( TestFlowAdaptor ) ADD_LBPM_TEST( TestMap ) +ADD_LBPM_TEST( TestMembrane ) #ADD_LBPM_TEST( TestMRT ) #ADD_LBPM_TEST( TestColorGrad ) ADD_LBPM_TEST( TestWideHalo ) diff --git a/tests/TestMembrane.cpp b/tests/TestMembrane.cpp new file mode 100644 index 00000000..b2abd78e --- /dev/null +++ b/tests/TestMembrane.cpp @@ -0,0 +1,256 @@ + +//************************************************************************* +// Lattice Boltzmann Simulator for Single Phase Flow in Porous Media +// James E. McCLure +//************************************************************************* +#include +#include +#include +#include "common/MPI.h" +#include "common/Membrane.h" +#include "common/ScaLBL.h" + +using namespace std; + +std::shared_ptr loadInputs( int nprocs ) +{ + //auto db = std::make_shared( "Domain.in" ); + auto db = std::make_shared(); + db->putScalar( "BC", 0 ); + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { 32, 32, 32 } ); + db->putScalar( "nspheres", 1 ); + db->putVector( "L", { 1, 1, 1 } ); + return db; +} + +//*************************************************************************************** +int main(int argc, char **argv) +{ + // Initialize MPI + Utilities::startup( argc, argv ); + Utilities::MPI comm( MPI_COMM_WORLD ); + int check=0; + { + + int i,j,k,n; + + int rank = comm.getRank(); + if (rank == 0){ + printf("********************************************************\n"); + printf("Running unit test: TestMembrane \n"); + printf("********************************************************\n"); + } + + // Load inputs + auto db = loadInputs( comm.getSize() ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + auto Dm = std::make_shared(db,comm); + + Nx += 2; + Ny += 2; + Nz += 2; + int N = Nx*Ny*Nz; + //....................................................................... + int Np = 0; + double distance,radius; + DoubleArray Distance(Nx,Ny,Nz); + for (k=0;kid[n] = 1; + radius = double(Nx)/4; + distance = sqrt(double((i-0.5*Nx)*(i-0.5*Nx)+ (j-0.5*Ny)*(j-0.5*Ny)+ (k-0.5*Nz)*(k-0.5*Nz)))-radius; + if (distance < 0.0 ){ + Dm->id[n] = 1; + } + Distance(i,j,k) = distance; + Np++; + } + } + } + Dm->CommInit(); + + // Create a communicator for the device (will use optimized layout) + std::shared_ptr ScaLBL_Comm(new ScaLBL_Communicator(Dm)); + //Create a second communicator based on the regular data layout + std::shared_ptr ScaLBL_Comm_Regular(new ScaLBL_Communicator(Dm)); + + if (rank==0){ + printf("Total domain size = %i \n",N); + printf("Reduced domain size = %i \n",Np); + } + + // LBM variables + if (rank==0) printf ("Set up the neighborlist \n"); + int Npad=Np+32; + int neighborSize=18*Npad*sizeof(int); + int *neighborList; + IntArray Map(Nx,Ny,Nz); + neighborList= new int[18*Npad]; + + //......................device distributions................................. + int *NeighborList; + int *dvcMap; + //........................................................................... + ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); + ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Npad); + ScaLBL_CopyToDevice(NeighborList, neighborList, 18*Np*sizeof(int)); + + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1); + comm.barrier(); + + double *dist; + dist = new double [19*Np]; + + // Check the neighborlist + printf("Check neighborlist: exterior %i, first interior %i last interior %i \n",ScaLBL_Comm->LastExterior(),ScaLBL_Comm->FirstInterior(),ScaLBL_Comm->LastInterior()); + for (int idx=0; idxLastExterior(); idx++){ + for (int q=0; q<18; q++){ + int nn = neighborList[q*Np+idx]%Np; + if (nn>Np) printf("neighborlist error (exterior) at q=%i, idx=%i \n",q,idx); + dist[q*Np + idx] = 0.0; + } + } + for (int idx=ScaLBL_Comm->FirstInterior(); idxLastInterior(); idx++){ + for (int q=0; q<18; q++){ + int nn = neighborList[q*Np+idx]%Np; + if (nn>Np) printf("neighborlist error (exterior) at q=%i, idx=%i \n",q,idx); + dist[q*Np + idx] = 0.0; + } + } + + /* create a membrane data structure */ + Membrane M(Dm, NeighborList, Np); + + int MembraneCount = M.Create(Dm, Distance, Map); + if (rank==0) printf (" Number of membrane links: %i \n", MembraneCount); + + /* create a tagged array to show where the mebrane is*/ + double *MembraneLinks; + MembraneLinks = new double [Nx*Ny*Nz]; + for (int n=0; n 0.f){ + Dm->id[n] = 127; + } + if (sum < 0.f){ + Dm->id[n] = 64; + } + } + } + } + if (argc > 1) + Dm->AggregateLabels("membrane.raw"); + + //........................................................................... + // Update GPU data structures + if (rank==0) printf ("Setting up device map and neighbor list \n"); + int *TmpMap; + TmpMap=new int[Np*sizeof(int)]; + for (k=1; kLastExterior(), + Np); + DoubleArray Result(Nx,Ny,Nz); + + ScaLBL_Comm->RegularLayout(Map, Ci, Result); + +/* for (k=1; k +#include +#include +#include +#include +#include +#include +#include + +#include "models/IonModel.h" +#include "models/StokesModel.h" +#include "models/PoissonSolver.h" +#include "models/MultiPhysController.h" +#include "common/Utilities.h" +#include "analysis/ElectroChemistry.h" + +using namespace std; + +//*************************************************************************** +// Test lattice-Boltzmann Ion Model coupled with Poisson equation +//*************************************************************************** + +int main(int argc, char **argv) +{ + // Initialize MPI and error handlers + Utilities::startup( argc, argv ); + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + + { // Limit scope so variables that contain communicators will free before MPI_Finialize + + if (rank == 0){ + printf("********************************************************\n"); + printf("Running LBPM electrokinetic single-fluid solver \n"); + printf("********************************************************\n"); + } + // Initialize compute device + int device=ScaLBL_SetDevice(rank); + NULL_USE( device ); + ScaLBL_DeviceBarrier(); + comm.barrier(); + + PROFILE_ENABLE(1); + //PROFILE_ENABLE_TRACE(); + //PROFILE_ENABLE_MEMORY(); + PROFILE_SYNCHRONIZE(); + PROFILE_START("Main"); + Utilities::setErrorHandlers(); + + auto filename = argv[1]; + ScaLBL_StokesModel StokesModel(rank,nprocs,comm); + ScaLBL_IonModel IonModel(rank,nprocs,comm); + ScaLBL_Poisson PoissonSolver(rank,nprocs,comm); + ScaLBL_Multiphys_Controller Study(rank,nprocs,comm);//multiphysics controller coordinating multi-model coupling + + // Load controller information + Study.ReadParams(filename); + + // Load user input database files for Navier-Stokes and Ion solvers + StokesModel.ReadParams(filename); + IonModel.ReadParams(filename); + + // Setup other model specific structures + StokesModel.SetDomain(); + StokesModel.ReadInput(); + StokesModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables + + IonModel.SetDomain(); + IonModel.ReadInput(); + IonModel.Create(); + IonModel.SetMembrane(); + + // Create analysis object + ElectroChemistryAnalyzer Analysis(IonModel.Dm); + + // Get internal iteration number + StokesModel.timestepMax = Study.getStokesNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv); + StokesModel.Initialize(); // initializing the model will set initial conditions for variables + + IonModel.timestepMax = Study.getIonNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv); + IonModel.Initialize(); + // Get maximal time converting factor based on Sotkes and Ion solvers + Study.getTimeConvMax_PNP_coupling(StokesModel.time_conv,IonModel.time_conv); + + // Initialize LB-Poisson model + PoissonSolver.ReadParams(filename); + PoissonSolver.SetDomain(); + PoissonSolver.ReadInput(); + PoissonSolver.Create(); + PoissonSolver.Initialize(Study.time_conv_max); + + int timestep=0; + while (timestep < Study.timestepMax){ + + timestep++; + PoissonSolver.Run(IonModel.ChargeDensity,timestep);//solve Poisson equtaion to get steady-state electrical potental + StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity + IonModel.RunMembrane(StokesModel.Velocity,PoissonSolver.ElectricField,PoissonSolver.Psi); //solve for ion transport with membrane + + timestep++;//AA operations + + if (timestep%Study.analysis_interval==0){ + Analysis.Basic(IonModel,PoissonSolver,StokesModel,timestep); + } + if (timestep%Study.visualization_interval==0){ + Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep); + /* PoissonSolver.getElectricPotential(timestep); + PoissonSolver.getElectricField(timestep); + IonModel.getIonConcentration(timestep); + StokesModel.getVelocity(timestep); + */ + } + } + + if (rank==0) printf("Save simulation raw data at maximum timestep\n"); + Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep); + + if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n"); + if (rank==0) printf("*************************************************************\n"); + + PROFILE_STOP("Main"); + PROFILE_SAVE("lbpm_electrokinetic_SingleFluid_simulator",1); + // **************************************************** + + } // Limit scope so variables that contain communicators will free before MPI_Finialize + + Utilities::shutdown(); +} + +