From d59c78b0ce4d459815a0552e881af674e78920da Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 25 Dec 2021 17:03:35 -0500 Subject: [PATCH 01/37] add membrane class --- analysis/Membrane.cpp | 339 ++++++++++++++++++++++++++++++++++++++++++ analysis/Membrane.h | 54 +++++++ 2 files changed, 393 insertions(+) create mode 100644 analysis/Membrane.cpp create mode 100644 analysis/Membrane.h diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp new file mode 100644 index 00000000..b387dc38 --- /dev/null +++ b/analysis/Membrane.cpp @@ -0,0 +1,339 @@ +/* Flow adaptor class for multiphase flow methods */ + +#include "analysis/FlowAdaptor.h" +#include "analysis/distance.h" +#include "analysis/morphology.h" + +Membrane::Membrane(ScaLBL_Communicator &ScaLBL, int *neighborList) { + +} + +Membrane::~Membrane() { + +} + +int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int *membrane, int Np){ + int mlink = 0; + int n, idx, neighbor; + double dist, locdist; + /* go through the neighborlist structure */ + for (idx=0; idx +#include +#include +#include +#include +#include +#include + +#include "common/ScaLBL.h" + +/** + * \class Membrane + * @brief + * The Membrane class operates on ScaLBL data structures to insert membrane + * + */ + +class Membrane { +public: + /** + * \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(ScaLBL_Communicator &ScaLBL, int *neighborList); + + /** + * \brief Destructor + */ + ~Membrane(); + + /** + * \brief Fast-forward interface motion + * \details Accelerate the movement of interfaces based on the time derivative + * Optional keys to control behavior can be specified in the input database: + * move_interface_cutoff -- identifies the diffuse interface region + * move_interface_factor -- determines how much to ``fast forward" + * @param Distance - signed distance to membrane + * @param Map - mapping between regular layout and compact layout + * @param neighborLIst - list of neighbors for each site + * @param membrane - links that form the membrane + * @param Np - number of sites in compact layout + */ + int Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int *membrane, int Np); + + +private: + +}; +#endif \ No newline at end of file From 26766ef69a8e8b16d8add25d7a463a443d9f1084 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 25 Dec 2021 17:58:53 -0500 Subject: [PATCH 02/37] update membrane structure --- analysis/Membrane.cpp | 151 +++++++++++------------------------------- analysis/Membrane.h | 7 +- common/ScaLBL.cpp | 6 +- common/ScaLBL.h | 8 +++ 4 files changed, 52 insertions(+), 120 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index b387dc38..a901e5e9 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -12,7 +12,7 @@ Membrane::~Membrane() { } -int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int *membrane, int Np){ +int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int Np){ int mlink = 0; int n, idx, neighbor; double dist, locdist; @@ -32,7 +32,6 @@ int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, in else if (!(idx<0)){ - int neighbor; // cycle through the neighbors of lattice site idx neighbor=Map(i-1,j,k); dist=Distance(i-1,j,k); if (dist*locdist < 0.0){ @@ -155,7 +154,8 @@ int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, in } /* allocate memory */ - membrane = new int [2*mlink]; + membraneLinks = new int [2*mlink]; + membraneDist = new double [2*mlink]; /* construct the membrane*/ mlink = 0; @@ -167,167 +167,94 @@ int Membrane::Create(DoubleArray &Distance, IntArray &Map, int *neighborList, in locdist=Distance(i,j,k); else if (!(idx<0)){ - - int neighbor; // cycle through the neighbors of lattice site idx - neighbor=Map(i-1,j,k); - dist=Distance(i-1,j,k); - if (dist*locdist < 0.0){ - neighborList[idx]=idx + 2*Np; - //membrane[2*mlink] = idx + 2*Np; - //membrane[2*mlink+1] = neighbor + 1*Np; - //mlink++; - } neighbor=Map(i+1,j,k); dist=Distance(i+1,j,k); if (dist*locdist < 0.0){ - neighborList[Np+idx] = idx + 1*Np; - membrane[2*mlink] = idx + 1*Np; - membrane[2*mlink+1] = neighbor + 2*Np; + membraneLinks[2*mlink] = idx + 1*Np; + membraneLinks[2*mlink+1] = neighbor + 2*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i,j-1,k); - dist=Distance(i,j-1,k); - if (dist*locdist < 0.0){ - neighborList[2*Np+idx]=idx + 4*Np; - //membrane[2*mlink] = idx + 4*Np; - //membrane[2*mlink+1] = neighbor + 3*Np; - //mlink++; - } - neighbor=Map(i,j+1,k); dist=Distance(i,j+1,k); if (dist*locdist < 0.0){ - neighborList[3*Np+idx]=idx + 3*Np; - membrane[2*mlink] = idx + 3*Np; - membrane[2*mlink+1] = neighbor + 4*Np; + membraneLinks[2*mlink] = idx + 3*Np; + membraneLinks[2*mlink+1] = neighbor + 4*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i,j,k-1); - dist=Distance(i,j,k-1); - if (dist*locdist < 0.0){ - neighborList[4*Np+idx]=idx + 6*Np; - //membrane[2*mlink] = idx + 6*Np; - //membrane[2*mlink+1] = neighbor + 5*Np; - //mlink++; - } - neighbor=Map(i,j,k+1); dist=Distance(i,j,k+1); if (dist*locdist < 0.0){ - neighborList[5*Np+idx]=idx + 5*Np; - membrane[2*mlink] = idx + 5*Np; - membrane[2*mlink+1] = neighbor + 6*Np; + membraneLinks[2*mlink] = idx + 5*Np; + membraneLinks[2*mlink+1] = neighbor + 6*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i-1,j-1,k); - dist=Distance(i-1,j-1,k); - if (dist*locdist < 0.0){ - neighborList[6*Np+idx]=idx + 8*Np; - //membrane[2*mlink] = idx + 8*Np; - //membrane[2*mlink+1] = neighbor + 7*Np - //mlink++; - } - neighbor=Map(i+1,j+1,k); dist=Distance(i+1,j+1,k); if (dist*locdist < 0.0){ - neighborList[7*Np+idx]=idx + 7*Np; - membrane[2*mlink] = idx + 7*Np; - membrane[2*mlink+1] = neighbor+8*Np; + membraneLinks[2*mlink] = idx + 7*Np; + membraneLinks[2*mlink+1] = neighbor+8*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i-1,j+1,k); - dist=Distance(i-1,j+1,k); - if (dist*locdist < 0.0){ - neighborList[8*Np+idx]=idx + 10*Np; - //membrane[2*mlink] = idx + 10*Np; - //membrane[2*mlink+1] = neighbor + 9*Np; - //mlink++; - } - neighbor=Map(i+1,j-1,k); dist=Distance(i+1,j-1,k); if (dist*locdist < 0.0){ - neighborList[9*Np+idx]=idx + 9*Np; - membrane[2*mlink] = idx + 9*Np; - membrane[2*mlink+1] = neighbor + 10*Np; + membraneLinks[2*mlink] = idx + 9*Np; + membraneLinks[2*mlink+1] = neighbor + 10*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i-1,j,k-1); - dist=Distance(i-1,j,k-1); - if (dist*locdist < 0.0){ - neighborList[10*Np+idx]=idx + 12*Np; - //membrane[2*mlink] = idx + 12*Np; - //membrane[2*mlink+1] = neighbor + 11*Np;; - //mlink++; - } - neighbor=Map(i+1,j,k+1); dist=Distance(i+1,j,k+1); if (dist*locdist < 0.0){ - neighborList[11*Np+idx]=idx + 11*Np; - membrane[2*mlink] = idx + 11*Np; - membrane[2*mlink+1] = neighbor + 12*Np; + membraneLinks[2*mlink] = idx + 11*Np; + membraneLinks[2*mlink+1] = neighbor + 12*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i-1,j,k+1); - dist=Distance(i-1,j,k+1); - if (dist*locdist < 0.0){ - neighborList[12*Np+idx]=idx + 14*Np; - //membrane[2*mlink] = idx + 14*Np; - //membrane[2*mlink+1] = neighbor + 13*Np; - //mlink++; - } - neighbor=Map(i+1,j,k-1); dist=Distance(i+1,j,k-1); if (dist*locdist < 0.0){ - neighborList[13*Np+idx]=idx + 13*Np; - membrane[2*mlink] = idx + 13*Np; - membrane[2*mlink+1] = neighbor + 14*Np; + membraneLinks[2*mlink] = idx + 13*Np; + membraneLinks[2*mlink+1] = neighbor + 14*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i,j-1,k-1); - dist=Distance(i,j-1,k-1); - if (dist*locdist < 0.0){ - neighborList[14*Np+idx]=idx + 16*Np; - //membrane[2*mlink] = idx + 16*Np; - //membrane[2*mlink+1] = neighbor + 15*Np; - //mlink++; - } - neighbor=Map(i,j+1,k+1); dist=Distance(i,j+1,k+1); if (dist*locdist < 0.0){ - neighborList[15*Np+idx]=idx + 15*Np; - membrane[2*mlink] = idx + 15*Np; - membrane[2*mlink+1] =neighbor + 16*Np; + membraneLinks[2*mlink] = idx + 15*Np; + membraneLinks[2*mlink+1] =neighbor + 16*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } - neighbor=Map(i,j-1,k+1); - dist=Distance(i,j-1,k+1); - if (dist*locdist < 0.0){ - neighborList[16*Np+idx]=idx + 18*Np; - //membrane[2*mlink] = idx + 18*Np; - //membrane[2*mlink+1] = neighbor + 17*Np; - //mlink++; - } - neighbor=Map(i,j+1,k-1); dist=Distance(i,j+1,k-1); if (dist*locdist < 0.0){ - neighborList[17*Np+idx]=idx + 17*Np; - membrane[2*mlink] = idx + 17*Np; - membrane[2*mlink+1] = neighbor + 18*Np; + membraneLinks[2*mlink] = idx + 17*Np; + membraneLinks[2*mlink+1] = neighbor + 18*Np; + membraneDist[2*mlink] = locdist; + membraneDist[2*mlink+1] = dist; mlink++; } } diff --git a/analysis/Membrane.h b/analysis/Membrane.h index 3d64899a..98c64715 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -21,6 +21,9 @@ class Membrane { public: + int *membraneLinks; // D3Q19 links that cross membrane + double *membraneDist; // distance to membrane for each linked site + /** * \brief Create a flow adaptor to operate on the LB model * @param ScaLBL - originating data structures @@ -45,9 +48,7 @@ public: * @param membrane - links that form the membrane * @param Np - number of sites in compact layout */ - int Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int *membrane, int Np); - - + int Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int Np); private: }; diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index e3c3baf7..80f39033 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 @@ -324,7 +322,6 @@ ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr Dm){ } - ScaLBL_Communicator::~ScaLBL_Communicator() { @@ -988,7 +985,6 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis return(Np); } - void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC) { diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 2be0121e..cb376eb9 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -699,6 +699,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(); From efd29363b6256bde2cdd01961404ad4afd070c7c Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 27 Dec 2021 19:47:49 -0500 Subject: [PATCH 03/37] membrane communications --- analysis/Membrane.cpp | 482 +++++++++++++++++++++++++++++++++++++++++- analysis/Membrane.h | 71 +++++-- common/ScaLBL.h | 1 - 3 files changed, 532 insertions(+), 22 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index a901e5e9..c91e6e0c 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -1,27 +1,399 @@ /* Flow adaptor class for multiphase flow methods */ -#include "analysis/FlowAdaptor.h" +#include "analysis/Membrane.h" #include "analysis/distance.h" -#include "analysis/morphology.h" -Membrane::Membrane(ScaLBL_Communicator &ScaLBL, int *neighborList) { +Membrane::Membrane(std::shared_ptr Dm, int *initialNeighborList) { + Np = Dm->Np; + neighborList = new int[18*Np]; + /* Copy neighborList */ + for (int idx=0; idxComm.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");\ + + /* 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 **) &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() { + 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( 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(DoubleArray &Distance, IntArray &Map, int *neighborList, int Np){ +int Membrane::D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list, int start, int count, + DoubleArray &Distance, double *recvDistance, int *d3q19_recvlist){ + int i,j,k,n,nn,idx; + double dist,locdist; + int * ReturnDist; + ReturnDist=new int [count]; + + int regularCount = 0; + int membraneCount = 0; + for (idx=0; idx Dm, DoubleArray &Distance, IntArray &Map){ int mlink = 0; int n, idx, neighbor; double dist, locdist; /* go through the neighborlist structure */ - for (idx=0; idxrecvList("X"),0,recvCount_X,dvcRecvDist_X); + D3Q19_MapRecv(-1,-1,0,Dm->recvList("X"),recvCount_X,recvCount_X,dvcRecvDist_X); + D3Q19_MapRecv(-1,1,0, Dm->recvList("X"),2*recvCount_X,recvCount_X,dvcRecvDist_X); + D3Q19_MapRecv(-1,0,-1,Dm->recvList("X"),3*recvCount_X,recvCount_X,dvcRecvDist_X); + D3Q19_MapRecv(-1,0,1, Dm->recvList("X"),4*recvCount_X,recvCount_X,dvcRecvDist_X); + //................................................................................... + //...Map recieve list for the x face: q=1,7,9,11,13.................................. + D3Q19_MapRecv(1,0,0, Dm->recvList("x"),0,recvCount_x,dvcRecvDist_x); + D3Q19_MapRecv(1,1,0, Dm->recvList("x"),recvCount_x,recvCount_x,dvcRecvDist_x); + D3Q19_MapRecv(1,-1,0,Dm->recvList("x"),2*recvCount_x,recvCount_x,dvcRecvDist_x); + D3Q19_MapRecv(1,0,1, Dm->recvList("x"),3*recvCount_x,recvCount_x,dvcRecvDist_x); + D3Q19_MapRecv(1,0,-1,Dm->recvList("x"),4*recvCount_x,recvCount_x,dvcRecvDist_x); + //................................................................................... + //...Map recieve list for the y face: q=4,8,9,16,18 ................................... + D3Q19_MapRecv(0,-1,0, Dm->recvList("Y"),0,recvCount_Y,dvcRecvDist_Y); + D3Q19_MapRecv(-1,-1,0,Dm->recvList("Y"),recvCount_Y,recvCount_Y,dvcRecvDist_Y); + D3Q19_MapRecv(1,-1,0, Dm->recvList("Y"),2*recvCount_Y,recvCount_Y,dvcRecvDist_Y); + D3Q19_MapRecv(0,-1,-1,Dm->recvList("Y"),3*recvCount_Y,recvCount_Y,dvcRecvDist_Y); + D3Q19_MapRecv(0,-1,1, Dm->recvList("Y"),4*recvCount_Y,recvCount_Y,dvcRecvDist_Y); + //................................................................................... + //...Map recieve list for the Y face: q=3,7,10,15,17 .................................. + D3Q19_MapRecv(0,1,0, Dm->recvList("y"),0,recvCount_y,dvcRecvDist_y); + D3Q19_MapRecv(1,1,0, Dm->recvList("y"),recvCount_y,recvCount_y,dvcRecvDist_y); + D3Q19_MapRecv(-1,1,0,Dm->recvList("y"),2*recvCount_y,recvCount_y,dvcRecvDist_y); + D3Q19_MapRecv(0,1,1, Dm->recvList("y"),3*recvCount_y,recvCount_y,dvcRecvDist_y); + D3Q19_MapRecv(0,1,-1,Dm->recvList("y"),4*recvCount_y,recvCount_y,dvcRecvDist_y); + //................................................................................... + //...Map recieve list for the z face<<<6,12,13,16,17).............................................. + D3Q19_MapRecv(0,0,-1, Dm->recvList("Z"),0,recvCount_Z,dvcRecvDist_Z); + D3Q19_MapRecv(-1,0,-1,Dm->recvList("Z"),recvCount_Z,recvCount_Z,dvcRecvDist_Z); + D3Q19_MapRecv(1,0,-1, Dm->recvList("Z"),2*recvCount_Z,recvCount_Z,dvcRecvDist_Z); + D3Q19_MapRecv(0,-1,-1,Dm->recvList("Z"),3*recvCount_Z,recvCount_Z,dvcRecvDist_Z); + D3Q19_MapRecv(0,1,-1, Dm->recvList("Z"),4*recvCount_Z,recvCount_Z,dvcRecvDist_Z); + //...Map recieve list for the Z face<<<5,11,14,15,18).............................................. + D3Q19_MapRecv(0,0,1, Dm->recvList("z"),0,recvCount_z,dvcRecvDist_z); + D3Q19_MapRecv(1,0,1, Dm->recvList("z"),recvCount_z,recvCount_z,dvcRecvDist_z); + D3Q19_MapRecv(-1,0,1,Dm->recvList("z"),2*recvCount_z,recvCount_z,dvcRecvDist_z); + D3Q19_MapRecv(0,1,1, Dm->recvList("z"),3*recvCount_z,recvCount_z,dvcRecvDist_z); + D3Q19_MapRecv(0,-1,1,Dm->recvList("z"),4*recvCount_z,recvCount_z,dvcRecvDist_z); + //.................................................................................. + //...Map recieve list for the xy edge <<<8)................................ + D3Q19_MapRecv(-1,-1,0,Dm->recvList("XY"),0,recvCount_XY,dvcRecvDist_XY); + //...Map recieve list for the Xy edge <<<9)................................ + D3Q19_MapRecv(1,-1,0,Dm->recvList("xY"),0,recvCount_xY,dvcRecvDist_xY); + //...Map recieve list for the xY edge <<<10)................................ + D3Q19_MapRecv(-1,1,0,Dm->recvList("Xy"),0,recvCount_Xy,dvcRecvDist_Xy); + //...Map recieve list for the XY edge <<<7)................................ + D3Q19_MapRecv(1,1,0,Dm->recvList("xy"),0,recvCount_xy,dvcRecvDist_xy); + //...Map recieve list for the xz edge <<<12)................................ + D3Q19_MapRecv(-1,0,-1,Dm->recvList("XZ"),0,recvCount_XZ,dvcRecvDist_XZ); + //...Map recieve list for the xZ edge <<<14)................................ + D3Q19_MapRecv(-1,0,1,Dm->recvList("Xz"),0,recvCount_Xz,dvcRecvDist_Xz); + //...Map recieve list for the Xz edge <<<13)................................ + D3Q19_MapRecv(1,0,-1,Dm->recvList("xZ"),0,recvCount_xZ,dvcRecvDist_xZ); + //...Map recieve list for the XZ edge <<<11)................................ + D3Q19_MapRecv(1,0,1,Dm->recvList("xz"),0,recvCount_xz,dvcRecvDist_xz); + //...Map recieve list for the yz edge <<<16)................................ + D3Q19_MapRecv(0,-1,-1,Dm->recvList("YZ"),0,recvCount_YZ,dvcRecvDist_YZ); + //...Map recieve list for the yZ edge <<<18)................................ + D3Q19_MapRecv(0,-1,1,Dm->recvList("Yz"),0,recvCount_Yz,dvcRecvDist_Yz); + //...Map recieve list for the Yz edge <<<17)................................ + D3Q19_MapRecv(0,1,-1,Dm->recvList("yZ"),0,recvCount_yZ,dvcRecvDist_yZ); + //...Map recieve list for the YZ edge <<<15)................................ + D3Q19_MapRecv(0,1,1,Dm->recvList("yz"),0,recvCount_yz,dvcRecvDist_yz); + //................................................................................... + + //...................................................................................... + 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; + //...................................................................................... + return mlink; } diff --git a/analysis/Membrane.h b/analysis/Membrane.h index 98c64715..fe207614 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -21,6 +21,7 @@ class Membrane { public: + int *neighborList; // modified neighborlist int *membraneLinks; // D3Q19 links that cross membrane double *membraneDist; // distance to membrane for each linked site @@ -29,7 +30,7 @@ public: * @param ScaLBL - originating data structures * @param neighborList - list of neighbors for each site */ - Membrane(ScaLBL_Communicator &ScaLBL, int *neighborList); + Membrane(ScaLBL_Communicator &ScaLBL, int *initialNeighborList); /** * \brief Destructor @@ -37,19 +38,65 @@ public: ~Membrane(); /** - * \brief Fast-forward interface motion - * \details Accelerate the movement of interfaces based on the time derivative - * Optional keys to control behavior can be specified in the input database: - * move_interface_cutoff -- identifies the diffuse interface region - * move_interface_factor -- determines how much to ``fast forward" + * \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 - * @param neighborLIst - list of neighbors for each site - * @param membrane - links that form the membrane - * @param Np - number of sites in compact layout + * @param Map - mapping between regular layout and compact layout */ - int Create(DoubleArray &Distance, IntArray &Map, int *neighborList, int Np); + int Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map); + private: - + int Np; + /** + * \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 Distance - signed distance to membrane + * @param recvDistance - distance values from neighboring processor + * @param d3q19_recvlist - device array with the saved list + * */ + int D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list, int start, int count, + DoubleArray &Distance, double *recvDistance, int *d3q19_recvlist); + //...................................................................................... + // 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_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; + //...................................................................................... + // 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; + // 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; + //...................................................................................... }; #endif \ No newline at end of file diff --git a/common/ScaLBL.h b/common/ScaLBL.h index cb376eb9..f10b150d 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -787,7 +787,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; From e491327c89b0be1a66d5bce92f44ab5e0c822be8 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 28 Dec 2021 11:09:00 -0500 Subject: [PATCH 04/37] working on new comm data structures --- analysis/Membrane.cpp | 41 +++++++++++++++++++++-------------------- 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index c91e6e0c..b63bab20 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -343,49 +343,50 @@ Membrane::~Membrane() { ScaLBL_FreeDeviceMemory( dvcRecvDist_YZ ); } -int Membrane::D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list, int start, int count, - DoubleArray &Distance, double *recvDistance, int *d3q19_recvlist){ +int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_p, int rank_q + const int *originalSendList, const int shift_x, const int shift_y, const int shift_z, + const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList){ + int i,j,k,n,nn,idx; double dist,locdist; - int * ReturnDist; - ReturnDist=new int [count]; - + int regularCount = 0; int membraneCount = 0; - for (idx=0; idx Dm, DoubleArray &Distance, IntArra } /* Re-organize communication based on membrane structure*/ - + membraneCount_x = D3Q19_MapSendRecv(Cx, Cy, Cz, Dm->sendList("x"), Distance, ); // MPI_COMM_SCALBL.barrier(); From 6d82dc83a2655cbc7b7c53a9315d6702b0e2d115 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 28 Dec 2021 16:12:49 -0500 Subject: [PATCH 05/37] add membrane unit test --- analysis/Membrane.cpp | 33 +++++-- analysis/Membrane.h | 6 +- tests/CMakeLists.txt | 1 + tests/TestMembrane.cpp | 198 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 230 insertions(+), 8 deletions(-) create mode 100644 tests/TestMembrane.cpp diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index b63bab20..6b9a99d8 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -343,12 +343,22 @@ Membrane::~Membrane() { ScaLBL_FreeDeviceMemory( dvcRecvDist_YZ ); } -int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_p, int rank_q - const int *originalSendList, const int shift_x, const int shift_y, const int shift_z, +int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_q, int rank_Q, + const int shift_x, const int shift_y, const int shift_z, + const int *originalSendList, const int sendCount, const int recvCount, const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList){ + int sendtag = 2389; + int recvtag = 2389; int i,j,k,n,nn,idx; double dist,locdist; + + int *SendList; + int *RecvList; + int *RecvList_q; + SendList = new int [sendCount]; // for this process + RecvList = new int [recvCount]; // filled in by the neighbor process + RecvList_q = new int [sendCount]; // goes to the neighbor process int regularCount = 0; int membraneCount = 0; @@ -377,16 +387,27 @@ int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int if (dist*locdist < 0.0){ /* store membrane values at the end */ membraneCount++; - membraneRecvList[sendCount-membraneCount] = nn; - membraneSendList[sendCount-membraneCount] = n; + RecvList_q[sendCount-membraneCount] = nn; + SendList[sendCount-membraneCount] = n; } else { - membraneRecvList[regularCount] = nn; - membraneSendList[regularCount++] = n; + RecvList_q[regularCount] = nn; + SendList[regularCount++] = n; } } + // send the recv list to the neighbor process + req1[0] = Dm->Comm.Isend(RecvList_q, sendCount, rank_q, sendtag); + req2[0] = Dm->Comm.Irecv(RecvList, recvCount, rank_Q, recvtag); + Dm->Barrier(); + // Return updated version to the device + ScaLBL_CopyToDevice(&membraneSendList[start], SendList, sendCount*sizeof(int)); + ScaLBL_CopyToDevice(&membraneRecvList[start], RecvList, recvCount*sizeof(int)); + // clean up the work arrays + delete [] SendList; + delete [] RecvList; + delete [] RecvList_q; return membraneCount; } diff --git a/analysis/Membrane.h b/analysis/Membrane.h index fe207614..92432c86 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -63,8 +63,10 @@ private: * @param recvDistance - distance values from neighboring processor * @param d3q19_recvlist - device array with the saved list * */ - int D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list, int start, int count, - DoubleArray &Distance, double *recvDistance, int *d3q19_recvlist); + int D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_q, int rank_Q, + const int shift_x, const int shift_y, const int shift_z, + const int *originalSendList, const int sendCount, const int recvCount, + const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList); //...................................................................................... // MPI ranks for all 18 neighbors //...................................................................................... diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5e4441e3..c391613b 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -59,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..bd7e8234 --- /dev/null +++ b/tests/TestMembrane.cpp @@ -0,0 +1,198 @@ + +//************************************************************************* +// Lattice Boltzmann Simulator for Single Phase Flow in Porous Media +// James E. McCLure +//************************************************************************* +#include +#include +#include +#include "common/ScaLBL.h" +#include "common/MPI.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", { 5, 5, 5 } ); + 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; + + static int D3Q19[18][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}, + {1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0}, + {1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1}, + {0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}}; + + int rank = comm.getRank(); + if (rank == 0){ + printf("********************************************************\n"); + printf("Running unit test: TestMap \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; + for (k=1;kid[n] = 1; + 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]; + + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1); + comm.barrier(); + + // 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); + + } + } + 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); + + } + } + + //......................device distributions................................. + int *NeighborList; + int *dvcMap; + //........................................................................... + ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); + ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Npad); + + //........................................................................... + // 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; kfirst_interior; idxlast_interior; idx++){ + n = TmpMap[idx]; + k = n/(Nx*Ny); + j = (n-Nx*Ny*k)/Nx; + i = n-Nx*Ny*k-Nx*j; + for (int q=1; q<19; q++){ + int nn = neighborList[(q-1)*Np+idx]; + double value=fq[nn]; + // 3D index of neighbor + int iq=i-D3Q19[q-1][0]; + int jq=j-D3Q19[q-1][1]; + int kq=k-D3Q19[q-1][2]; + if (iq==0) iq=1; + if (jq==0) jq=1; + if (kq==0) kq=1; + if (iq==Nx-1) iq=Nx-2; + if (jq==Ny-1) jq=Ny-2; + if (kq==Nz-1) kq=Nz-2; + double check = kq*100.f+jq*10.f+iq*1.f+q*0.01; + if (value != check) + printf("Neighbor q=%i, i=%i,j=%i,k=%i: %f \n",q,iq,jq,kq,value); + } + } + + delete [] TmpMap; + + } + Utilities::shutdown(); + + return check; +} + From 74bf7dbec45d6bb25072e56b08711df91f9c54d7 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 28 Dec 2021 16:59:44 -0500 Subject: [PATCH 06/37] membrane compiles --- analysis/Membrane.cpp | 28 +++++++++++++++------------- analysis/Membrane.h | 29 ++++++++++++++++++++++++----- tests/TestMembrane.cpp | 5 ++++- 3 files changed, 43 insertions(+), 19 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index 6b9a99d8..b8d05b01 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -3,9 +3,9 @@ #include "analysis/Membrane.h" #include "analysis/distance.h" -Membrane::Membrane(std::shared_ptr Dm, int *initialNeighborList) { +Membrane::Membrane(std::shared_ptr Dm, int *initialNeighborList, int Nsites) { - Np = Dm->Np; + Np = Nsites; neighborList = new int[18*Np]; /* Copy neighborList */ for (int idx=0; idx Dm, const int *originalSendList, const int sendCount, const int recvCount, + int startSend, int startRecv, const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList){ int sendtag = 2389; @@ -398,11 +399,11 @@ int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int // send the recv list to the neighbor process req1[0] = Dm->Comm.Isend(RecvList_q, sendCount, rank_q, sendtag); req2[0] = Dm->Comm.Irecv(RecvList, recvCount, rank_Q, recvtag); - Dm->Barrier(); + Dm->Comm.barrier(); // Return updated version to the device - ScaLBL_CopyToDevice(&membraneSendList[start], SendList, sendCount*sizeof(int)); - ScaLBL_CopyToDevice(&membraneRecvList[start], RecvList, recvCount*sizeof(int)); + ScaLBL_CopyToDevice(&membraneSendList[startSend], SendList, sendCount*sizeof(int)); + ScaLBL_CopyToDevice(&membraneRecvList[startRecv], RecvList, recvCount*sizeof(int)); // clean up the work arrays delete [] SendList; @@ -413,6 +414,7 @@ int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map){ int mlink = 0; + int i,j,k; int n, idx, neighbor; double dist, locdist; /* go through the neighborlist structure */ @@ -424,7 +426,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra idx=Map(i,j,k); locdist=Distance(i,j,k); - else if (!(idx<0)){ + if (!(idx<0)){ neighbor=Map(i-1,j,k); dist=Distance(i-1,j,k); @@ -560,7 +562,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra idx=Map(i,j,k); locdist=Distance(i,j,k); - else if (!(idx<0)){ + if (!(idx<0)){ neighbor=Map(i+1,j,k); dist=Distance(i+1,j,k); @@ -657,7 +659,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra } /* Re-organize communication based on membrane structure*/ - membraneCount_x = D3Q19_MapSendRecv(Cx, Cy, Cz, Dm->sendList("x"), Distance, ); + //membraneCount_x = D3Q19_MapSendRecv(Cx, Cy, Cz, Dm->sendList("x"), Distance, ); // MPI_COMM_SCALBL.barrier(); @@ -665,7 +667,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra // Set up the recieve distribution lists //................................................................................... //...Map recieve list for the X face: q=2,8,10,12,14 ................................. - D3Q19_MapRecv(-1,0,0, Dm->recvList("X"),0,recvCount_X,dvcRecvDist_X); +/* D3Q19_MapRecv(-1,0,0, Dm->recvList("X"),0,recvCount_X,dvcRecvDist_X); D3Q19_MapRecv(-1,-1,0,Dm->recvList("X"),recvCount_X,recvCount_X,dvcRecvDist_X); D3Q19_MapRecv(-1,1,0, Dm->recvList("X"),2*recvCount_X,recvCount_X,dvcRecvDist_X); D3Q19_MapRecv(-1,0,-1,Dm->recvList("X"),3*recvCount_X,recvCount_X,dvcRecvDist_X); @@ -747,6 +749,6 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra CommunicationCount = SendCount+RecvCount; //...................................................................................... - + */ return mlink; } diff --git a/analysis/Membrane.h b/analysis/Membrane.h index 92432c86..da655afe 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -21,6 +21,9 @@ class Membrane { public: + int Np; + int Nx,Ny,Nz,N; + int *neighborList; // modified neighborlist int *membraneLinks; // D3Q19 links that cross membrane double *membraneDist; // distance to membrane for each linked site @@ -30,7 +33,7 @@ public: * @param ScaLBL - originating data structures * @param neighborList - list of neighbors for each site */ - Membrane(ScaLBL_Communicator &ScaLBL, int *initialNeighborList); + Membrane(std::shared_ptr Dm, int *initialNeighborList, int Nsites); /** * \brief Destructor @@ -46,8 +49,23 @@ public: */ int Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map); + //...................................................................................... + // 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: - int Np; + 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 @@ -63,9 +81,10 @@ private: * @param recvDistance - distance values from neighboring processor * @param d3q19_recvlist - device array with the saved list * */ - int D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int rank_q, int rank_Q, - const int shift_x, const int shift_y, const int shift_z, - const int *originalSendList, const int sendCount, const int recvCount, + int D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, + int rank_q, int rank_Q, const int shift_x, const int shift_y, const int shift_z, + std::shared_ptr Dm, const int *originalSendList, const int sendCount, const int recvCount, + int startSend, int startRecv, const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList); //...................................................................................... // MPI ranks for all 18 neighbors diff --git a/tests/TestMembrane.cpp b/tests/TestMembrane.cpp index bd7e8234..eb3752c2 100644 --- a/tests/TestMembrane.cpp +++ b/tests/TestMembrane.cpp @@ -6,8 +6,8 @@ #include #include #include -#include "common/ScaLBL.h" #include "common/MPI.h" +#include "analysis/Membrane.h" using namespace std; @@ -107,6 +107,9 @@ int main(int argc, char **argv) } } + + /* create a membrane data structure */ + Membrane M(Dm, neighborList, Np); //......................device distributions................................. int *NeighborList; From ec76dc9a3b10ad276e52c0721d8974e81b22bc6d Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 28 Dec 2021 17:21:11 -0500 Subject: [PATCH 07/37] membrane test --- tests/TestMembrane.cpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/tests/TestMembrane.cpp b/tests/TestMembrane.cpp index eb3752c2..f7cdc50b 100644 --- a/tests/TestMembrane.cpp +++ b/tests/TestMembrane.cpp @@ -59,11 +59,19 @@ int main(int argc, char **argv) int N = Nx*Ny*Nz; //....................................................................... int Np = 0; - for (k=1;kid[n] = 1; + radius = double(Nx)/6; + 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++; } } From 0fb6fbfd7a41d0c37faeac37ad507dc921cbeb6e Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 12 Jan 2022 09:48:04 -0500 Subject: [PATCH 08/37] update membrane test --- analysis/Membrane.cpp | 11 ++++++----- analysis/Membrane.h | 8 +++++--- sample_scripts/configure_arden | 2 +- tests/TestMembrane.cpp | 17 ++++------------- 4 files changed, 16 insertions(+), 22 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index b8d05b01..4d355b2a 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -344,9 +344,10 @@ Membrane::~Membrane() { } int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, - int rank_q, int rank_Q, const int shift_x, const int shift_y, const int shift_z, - std::shared_ptr Dm, const int *originalSendList, const int sendCount, const int recvCount, - int startSend, int startRecv, + const int shift_x, const int shift_y, const int shift_z, + const int sendCount, const int recvCount, + int rank_q, int rank_Q, int startSend, int startRecv, + std::shared_ptr Dm, const int *originalSendList, const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList){ int sendtag = 2389; @@ -414,8 +415,8 @@ int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map){ int mlink = 0; - int i,j,k; - int n, idx, neighbor; + int i,j,k,n; + int idx, neighbor; double dist, locdist; /* go through the neighborlist structure */ /* count & cut the links */ diff --git a/analysis/Membrane.h b/analysis/Membrane.h index da655afe..3a842b75 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -49,6 +49,7 @@ public: */ int Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map); + void SendRecv(double *dist); //...................................................................................... // Buffers to store data sent and recieved by this MPI process double *sendbuf_x, *sendbuf_y, *sendbuf_z, *sendbuf_X, *sendbuf_Y, *sendbuf_Z; @@ -82,9 +83,10 @@ private: * @param d3q19_recvlist - device array with the saved list * */ int D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, - int rank_q, int rank_Q, const int shift_x, const int shift_y, const int shift_z, - std::shared_ptr Dm, const int *originalSendList, const int sendCount, const int recvCount, - int startSend, int startRecv, + const int shift_x, const int shift_y, const int shift_z, + const int sendCount, const int recvCount, + int rank_q, int rank_Q, int startSend, int startRecv, + std::shared_ptr Dm, const int *originalSendList, const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList); //...................................................................................... // MPI ranks for all 18 neighbors diff --git a/sample_scripts/configure_arden b/sample_scripts/configure_arden index 89ffb766..58d303f9 100755 --- a/sample_scripts/configure_arden +++ b/sample_scripts/configure_arden @@ -15,6 +15,6 @@ cmake -D CMAKE_C_COMPILER:PATH=/opt/arden/openmpi/3.1.2/bin/mpicc \ -D HDF5_DIRECTORY="/opt/arden/hdf5/1.8.12" \ -D USE_CUDA=0 \ -D USE_TIMER=0 \ - ~/Programs/LBPM + ~/Programs/LBPM-WIA # -D HDF5_LIB="/opt/arden/hdf5/1.8.12/lib/libhdf5.a"\ diff --git a/tests/TestMembrane.cpp b/tests/TestMembrane.cpp index f7cdc50b..ecdbd4d6 100644 --- a/tests/TestMembrane.cpp +++ b/tests/TestMembrane.cpp @@ -119,6 +119,9 @@ int main(int argc, char **argv) /* create a membrane data structure */ Membrane M(Dm, neighborList, Np); + int MembraneCount = M.Create(Dm, Distance, Map); + + //......................device distributions................................. int *NeighborList; int *dvcMap; @@ -159,19 +162,7 @@ int main(int argc, char **argv) } } } - /* for (int idx=0; idx Date: Fri, 14 Jan 2022 10:59:04 -0500 Subject: [PATCH 09/37] update membrane test --- analysis/Membrane.cpp | 6 ++---- tests/TestMembrane.cpp | 49 ++++++++++++++++++++++++++++++++++++------ 2 files changed, 44 insertions(+), 11 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index 4d355b2a..4618380f 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -372,7 +372,7 @@ int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, // if (rank ==0) printf("@ Get 3D indices from the send process: i=%d, j=%d, k=%d\n",i,j,k); /* distance to membrane at the send site */ - locdist = Distance(i,j,k); + dist = Distance(i,j,k); // Streaming for the non-local distribution i += Cqx; j += Cqy; k += Cqz; @@ -415,7 +415,7 @@ int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map){ int mlink = 0; - int i,j,k,n; + int i,j,k; int idx, neighbor; double dist, locdist; /* go through the neighborlist structure */ @@ -423,7 +423,6 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra for (k=1;k Dm, DoubleArray &Distance, IntArra for (k=1;k loadInputs( int nprocs ) auto db = std::make_shared(); db->putScalar( "BC", 0 ); db->putVector( "nproc", { 1, 1, 1 } ); - db->putVector( "n", { 5, 5, 5 } ); + db->putVector( "n", { 32, 32, 32 } ); db->putScalar( "nspheres", 1 ); db->putVector( "L", { 1, 1, 1 } ); return db; @@ -42,7 +42,7 @@ int main(int argc, char **argv) int rank = comm.getRank(); if (rank == 0){ printf("********************************************************\n"); - printf("Running unit test: TestMap \n"); + printf("Running unit test: TestMembrane \n"); printf("********************************************************\n"); } @@ -66,7 +66,7 @@ int main(int argc, char **argv) for (i=0;iid[n] = 1; - radius = double(Nx)/6; + 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; @@ -99,29 +99,64 @@ int main(int argc, char **argv) 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 (argc > 1) + Dm->AggregateLabels("membrane.raw"); + //......................device distributions................................. int *NeighborList; int *dvcMap; From 66fc8556088e68dda83bd4260ccf435c3379a2fe Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 19 Jan 2022 16:18:29 -0500 Subject: [PATCH 10/37] convention for inside / outside membrane link direction --- analysis/Membrane.cpp | 148 ++++++++++++++++++++++++++++++----------- analysis/Membrane.h | 1 + tests/TestMembrane.cpp | 6 +- 3 files changed, 116 insertions(+), 39 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index 4618380f..99e05234 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -552,9 +552,11 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra /* allocate memory */ membraneLinks = new int [2*mlink]; membraneDist = new double [2*mlink]; - + membraneTag = new int [mlink]; + /* construct the membrane*/ mlink = 0; + int insideMem = 0; int outsideMem = 0; for (k=1;k Dm, DoubleArray &Distance, IntArra neighbor=Map(i+1,j,k); dist=Distance(i+1,j,k); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 1*Np; - membraneLinks[2*mlink+1] = neighbor + 2*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 1*Np; + membraneLinks[outsideMem] = neighbor + 2*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i,j+1,k); dist=Distance(i,j+1,k); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 3*Np; - membraneLinks[2*mlink+1] = neighbor + 4*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 3*Np; + membraneLinks[outsideMem] = neighbor + 4*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i,j,k+1); dist=Distance(i,j,k+1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 5*Np; - membraneLinks[2*mlink+1] = neighbor + 6*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 5*Np; + membraneLinks[outsideMem] = neighbor + 6*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i+1,j+1,k); dist=Distance(i+1,j+1,k); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 7*Np; - membraneLinks[2*mlink+1] = neighbor+8*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 7*Np; + membraneLinks[outsideMem] = neighbor+8*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i+1,j-1,k); dist=Distance(i+1,j-1,k); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 9*Np; - membraneLinks[2*mlink+1] = neighbor + 10*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 9*Np; + membraneLinks[outsideMem] = neighbor + 10*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i+1,j,k+1); dist=Distance(i+1,j,k+1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 11*Np; - membraneLinks[2*mlink+1] = neighbor + 12*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 11*Np; + membraneLinks[outsideMem] = neighbor + 12*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i+1,j,k-1); dist=Distance(i+1,j,k-1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 13*Np; - membraneLinks[2*mlink+1] = neighbor + 14*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 13*Np; + membraneLinks[outsideMem] = neighbor + 14*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i,j+1,k+1); dist=Distance(i,j+1,k+1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 15*Np; - membraneLinks[2*mlink+1] =neighbor + 16*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 15*Np; + membraneLinks[outsideMem] =neighbor + 16*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } neighbor=Map(i,j+1,k-1); dist=Distance(i,j+1,k-1); if (dist*locdist < 0.0){ - membraneLinks[2*mlink] = idx + 17*Np; - membraneLinks[2*mlink+1] = neighbor + 18*Np; - membraneDist[2*mlink] = locdist; - membraneDist[2*mlink+1] = dist; + if (locdist < 0.0){ + insideMem = 2*mlink; + outsideMem = 2*mlink+1; + } + else{ + insideMem = 2*mlink+1; + outsideMem = 2*mlink; + } + membraneLinks[insideMem] = idx + 17*Np; + membraneLinks[outsideMem] = neighbor + 18*Np; + membraneDist[insideMem] = locdist; + membraneDist[outsideMem] = dist; mlink++; } } diff --git a/analysis/Membrane.h b/analysis/Membrane.h index 3a842b75..563af5b0 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -26,6 +26,7 @@ public: int *neighborList; // modified neighborlist int *membraneLinks; // D3Q19 links that cross membrane + int *membraneTag; // label each link in the membrane double *membraneDist; // distance to membrane for each linked site /** diff --git a/tests/TestMembrane.cpp b/tests/TestMembrane.cpp index 319992b0..048be117 100644 --- a/tests/TestMembrane.cpp +++ b/tests/TestMembrane.cpp @@ -118,7 +118,6 @@ int main(int argc, char **argv) dist[q*Np + idx] = 0.0; } } - /* create a membrane data structure */ Membrane M(Dm, neighborList, Np); @@ -135,7 +134,7 @@ int main(int argc, char **argv) for (int mlink=0; mlink 0.f){ Dm->id[n] = 127; } + if (sum < 0.f){ + Dm->id[n] = 64; + } } } } From 7d679571f19fc3288558bc36adf2a2117c681cd5 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 28 Jan 2022 05:16:37 -0500 Subject: [PATCH 11/37] working on membrane comm --- analysis/Membrane.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index 99e05234..2dd6e3b5 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -735,7 +735,13 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra //membraneCount_x = D3Q19_MapSendRecv(Cx, Cy, Cz, Dm->sendList("x"), Distance, ); // MPI_COMM_SCALBL.barrier(); - + /* int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, + const int shift_x, const int shift_y, const int shift_z, + const int sendCount, const int recvCount, + int rank_q, int rank_Q, int startSend, int startRecv, + std::shared_ptr Dm, const int *originalSendList, + const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList) + */ //................................................................................... // Set up the recieve distribution lists //................................................................................... From 5bbda39fdbbac7925b5f72d16d229c4151d38d89 Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 21 Feb 2022 16:39:24 -0500 Subject: [PATCH 12/37] working on membrane communication structures --- analysis/Membrane.cpp | 477 ++++++++++++++++++++++++++++++----------- analysis/Membrane.h | 59 ++++- common/ScaLBL.cpp | 1 - cpu/D3Q19.cpp | 1 + cpu/MembraneHelper.cpp | 42 ++++ 5 files changed, 443 insertions(+), 137 deletions(-) create mode 100644 cpu/MembraneHelper.cpp diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index 2dd6e3b5..e90709ec 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -5,6 +5,11 @@ Membrane::Membrane(std::shared_ptr Dm, int *initialNeighborList, int Nsites) { + Lock=false; // unlock the communicator + //...................................................................................... + // Create a separate copy of the communicator for the device + MPI_COMM_SCALBL = Dm->Comm.dup(); + Np = Nsites; neighborList = new int[18*Np]; /* Copy neighborList */ @@ -187,6 +192,25 @@ Membrane::Membrane(std::shared_ptr Dm, int *initialNeighborList, int Ns 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 @@ -323,6 +347,24 @@ Membrane::~Membrane() { 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 ); @@ -343,76 +385,6 @@ Membrane::~Membrane() { ScaLBL_FreeDeviceMemory( dvcRecvDist_YZ ); } -int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, - const int shift_x, const int shift_y, const int shift_z, - const int sendCount, const int recvCount, - int rank_q, int rank_Q, int startSend, int startRecv, - std::shared_ptr Dm, const int *originalSendList, - const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList){ - - int sendtag = 2389; - int recvtag = 2389; - int i,j,k,n,nn,idx; - double dist,locdist; - - int *SendList; - int *RecvList; - int *RecvList_q; - SendList = new int [sendCount]; // for this process - RecvList = new int [recvCount]; // filled in by the neighbor process - RecvList_q = new int [sendCount]; // goes to the neighbor process - - int regularCount = 0; - int membraneCount = 0; - for (idx=0; idxComm.Isend(RecvList_q, sendCount, rank_q, sendtag); - req2[0] = Dm->Comm.Irecv(RecvList, recvCount, rank_Q, recvtag); - Dm->Comm.barrier(); - - // Return updated version to the device - ScaLBL_CopyToDevice(&membraneSendList[startSend], SendList, sendCount*sizeof(int)); - ScaLBL_CopyToDevice(&membraneRecvList[startRecv], RecvList, recvCount*sizeof(int)); - - // clean up the work arrays - delete [] SendList; - delete [] RecvList; - delete [] RecvList_q; - return membraneCount; -} - int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map){ int mlink = 0; int i,j,k; @@ -718,6 +690,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra } else{ insideMem = 2*mlink+1; + outsideMem = 2*mlink; } membraneLinks[insideMem] = idx + 17*Np; @@ -732,84 +705,71 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra } /* Re-organize communication based on membrane structure*/ - //membraneCount_x = D3Q19_MapSendRecv(Cx, Cy, Cz, Dm->sendList("x"), Distance, ); - // MPI_COMM_SCALBL.barrier(); - - /* int Membrane::D3Q19_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, - const int shift_x, const int shift_y, const int shift_z, - const int sendCount, const int recvCount, - int rank_q, int rank_Q, int startSend, int startRecv, - std::shared_ptr Dm, const int *originalSendList, - const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList) - */ - //................................................................................... - // Set up the recieve distribution lists - //................................................................................... //...Map recieve list for the X face: q=2,8,10,12,14 ................................. -/* D3Q19_MapRecv(-1,0,0, Dm->recvList("X"),0,recvCount_X,dvcRecvDist_X); - D3Q19_MapRecv(-1,-1,0,Dm->recvList("X"),recvCount_X,recvCount_X,dvcRecvDist_X); - D3Q19_MapRecv(-1,1,0, Dm->recvList("X"),2*recvCount_X,recvCount_X,dvcRecvDist_X); - D3Q19_MapRecv(-1,0,-1,Dm->recvList("X"),3*recvCount_X,recvCount_X,dvcRecvDist_X); - D3Q19_MapRecv(-1,0,1, Dm->recvList("X"),4*recvCount_X,recvCount_X,dvcRecvDist_X); + 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 x face: q=1,7,9,11,13.................................. - D3Q19_MapRecv(1,0,0, Dm->recvList("x"),0,recvCount_x,dvcRecvDist_x); - D3Q19_MapRecv(1,1,0, Dm->recvList("x"),recvCount_x,recvCount_x,dvcRecvDist_x); - D3Q19_MapRecv(1,-1,0,Dm->recvList("x"),2*recvCount_x,recvCount_x,dvcRecvDist_x); - D3Q19_MapRecv(1,0,1, Dm->recvList("x"),3*recvCount_x,recvCount_x,dvcRecvDist_x); - D3Q19_MapRecv(1,0,-1,Dm->recvList("x"),4*recvCount_x,recvCount_x,dvcRecvDist_x); + 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 ................................... - D3Q19_MapRecv(0,-1,0, Dm->recvList("Y"),0,recvCount_Y,dvcRecvDist_Y); - D3Q19_MapRecv(-1,-1,0,Dm->recvList("Y"),recvCount_Y,recvCount_Y,dvcRecvDist_Y); - D3Q19_MapRecv(1,-1,0, Dm->recvList("Y"),2*recvCount_Y,recvCount_Y,dvcRecvDist_Y); - D3Q19_MapRecv(0,-1,-1,Dm->recvList("Y"),3*recvCount_Y,recvCount_Y,dvcRecvDist_Y); - D3Q19_MapRecv(0,-1,1, Dm->recvList("Y"),4*recvCount_Y,recvCount_Y,dvcRecvDist_Y); + 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 .................................. - D3Q19_MapRecv(0,1,0, Dm->recvList("y"),0,recvCount_y,dvcRecvDist_y); - D3Q19_MapRecv(1,1,0, Dm->recvList("y"),recvCount_y,recvCount_y,dvcRecvDist_y); - D3Q19_MapRecv(-1,1,0,Dm->recvList("y"),2*recvCount_y,recvCount_y,dvcRecvDist_y); - D3Q19_MapRecv(0,1,1, Dm->recvList("y"),3*recvCount_y,recvCount_y,dvcRecvDist_y); - D3Q19_MapRecv(0,1,-1,Dm->recvList("y"),4*recvCount_y,recvCount_y,dvcRecvDist_y); + 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).............................................. - D3Q19_MapRecv(0,0,-1, Dm->recvList("Z"),0,recvCount_Z,dvcRecvDist_Z); - D3Q19_MapRecv(-1,0,-1,Dm->recvList("Z"),recvCount_Z,recvCount_Z,dvcRecvDist_Z); - D3Q19_MapRecv(1,0,-1, Dm->recvList("Z"),2*recvCount_Z,recvCount_Z,dvcRecvDist_Z); - D3Q19_MapRecv(0,-1,-1,Dm->recvList("Z"),3*recvCount_Z,recvCount_Z,dvcRecvDist_Z); - D3Q19_MapRecv(0,1,-1, Dm->recvList("Z"),4*recvCount_Z,recvCount_Z,dvcRecvDist_Z); + 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).............................................. - D3Q19_MapRecv(0,0,1, Dm->recvList("z"),0,recvCount_z,dvcRecvDist_z); - D3Q19_MapRecv(1,0,1, Dm->recvList("z"),recvCount_z,recvCount_z,dvcRecvDist_z); - D3Q19_MapRecv(-1,0,1,Dm->recvList("z"),2*recvCount_z,recvCount_z,dvcRecvDist_z); - D3Q19_MapRecv(0,1,1, Dm->recvList("z"),3*recvCount_z,recvCount_z,dvcRecvDist_z); - D3Q19_MapRecv(0,-1,1,Dm->recvList("z"),4*recvCount_z,recvCount_z,dvcRecvDist_z); + 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)................................ - D3Q19_MapRecv(-1,-1,0,Dm->recvList("XY"),0,recvCount_XY,dvcRecvDist_XY); + 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)................................ - D3Q19_MapRecv(1,-1,0,Dm->recvList("xY"),0,recvCount_xY,dvcRecvDist_xY); + 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)................................ - D3Q19_MapRecv(-1,1,0,Dm->recvList("Xy"),0,recvCount_Xy,dvcRecvDist_Xy); + 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)................................ - D3Q19_MapRecv(1,1,0,Dm->recvList("xy"),0,recvCount_xy,dvcRecvDist_xy); + 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)................................ - D3Q19_MapRecv(-1,0,-1,Dm->recvList("XZ"),0,recvCount_XZ,dvcRecvDist_XZ); + 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)................................ - D3Q19_MapRecv(-1,0,1,Dm->recvList("Xz"),0,recvCount_Xz,dvcRecvDist_Xz); + 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)................................ - D3Q19_MapRecv(1,0,-1,Dm->recvList("xZ"),0,recvCount_xZ,dvcRecvDist_xZ); + 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)................................ - D3Q19_MapRecv(1,0,1,Dm->recvList("xz"),0,recvCount_xz,dvcRecvDist_xz); + 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)................................ - D3Q19_MapRecv(0,-1,-1,Dm->recvList("YZ"),0,recvCount_YZ,dvcRecvDist_YZ); + 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)................................ - D3Q19_MapRecv(0,-1,1,Dm->recvList("Yz"),0,recvCount_Yz,dvcRecvDist_Yz); + 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)................................ - D3Q19_MapRecv(0,1,-1,Dm->recvList("yZ"),0,recvCount_yZ,dvcRecvDist_yZ); + 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)................................ - D3Q19_MapRecv(0,1,1,Dm->recvList("yz"),0,recvCount_yz,dvcRecvDist_yz); + linkCount_yz = D3Q19_MapRecv(0,1,1,Dm->recvList("yz"),0,recvCount_yz,dvcRecvDist_yz,dvcRecvLinks_yz,Distance,Map); //................................................................................... //...................................................................................... @@ -828,6 +788,271 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra CommunicationCount = SendCount+RecvCount; //...................................................................................... - */ + return mlink; } + +int Membrane::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){ + + int linkCount = 0; + int memLinkCount=0; + int i,j,k,n,nn,idx; + int * ReturnDist; + double distanceNonLocal,distanceLocal; + ReturnDist=new int [count]; + int *LinkList=new int [count]; + int *memLinkList=new int [count]; + + for (idx=0; idx Dm, DoubleArray &Distance, IntArray &Map); - - void SendRecv(double *dist); + + void SendD3Q19AA(double *dist); + void RecvD3Q19AA(double *dist); //...................................................................................... // Buffers to store data sent and recieved by this MPI process double *sendbuf_x, *sendbuf_y, *sendbuf_z, *sendbuf_X, *sendbuf_Y, *sendbuf_Z; @@ -62,6 +93,8 @@ public: //...................................................................................... 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 @@ -79,16 +112,13 @@ private: * @param Cqz - discrete velocity (z) * @param list - list of recieved values * @param count - number recieved values - * @param Distance - signed distance to membrane - * @param recvDistance - distance values from neighboring processor * @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_MapSendRecv(const int Cqx, const int Cqy, const int Cqz, - const int shift_x, const int shift_y, const int shift_z, - const int sendCount, const int recvCount, - int rank_q, int rank_Q, int startSend, int startRecv, - std::shared_ptr Dm, const int *originalSendList, - const DoubleArray &Distance, int *membraneSendList, int *membraneRecvList); + 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 //...................................................................................... @@ -100,6 +130,7 @@ private: 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; @@ -109,6 +140,10 @@ private: 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; @@ -117,6 +152,10 @@ private: 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; diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index 80f39033..cefb0f85 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -1386,7 +1386,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/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/MembraneHelper.cpp b/cpu/MembraneHelper.cpp new file mode 100644 index 00000000..7c604d8b --- /dev/null +++ b/cpu/MembraneHelper.cpp @@ -0,0 +1,42 @@ + +extern "C" void Membrane_D3Q19_Unpack(int q, int *list, int *links, int start, int linkCount, + double *recvbuf, double *dist, int N) { + //.................................................................................... + // 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; + for (link=0; link Date: Mon, 21 Feb 2022 16:49:42 -0500 Subject: [PATCH 13/37] add cell simulator --- tests/CMakeLists.txt | 1 + tests/lbpm_cell_simulator.cpp | 131 ++++++++++++++++++++++++++++++++++ 2 files changed, 132 insertions(+) create mode 100644 tests/lbpm_cell_simulator.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c391613b..f17c3f6e 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 ) diff --git a/tests/lbpm_cell_simulator.cpp b/tests/lbpm_cell_simulator.cpp new file mode 100644 index 00000000..caeef89a --- /dev/null +++ b/tests/lbpm_cell_simulator.cpp @@ -0,0 +1,131 @@ +#include +#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(); + + // 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.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential + + 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(); +} + + From dd9da41d4b84426ee2b12745cc2396e830f478f2 Mon Sep 17 00:00:00 2001 From: James McClure Date: Thu, 24 Feb 2022 21:43:18 -0500 Subject: [PATCH 14/37] added cell simulator --- cpu/MembraneHelper.cpp | 2 +- tests/lbpm_cell_simulator.cpp | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/cpu/MembraneHelper.cpp b/cpu/MembraneHelper.cpp index 7c604d8b..960e2f20 100644 --- a/cpu/MembraneHelper.cpp +++ b/cpu/MembraneHelper.cpp @@ -39,4 +39,4 @@ extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links, double *c if (!(n < 0)) dist[q * N + n] = alpha*recvbuf[start + idx]; } -} \ No newline at end of file +} diff --git a/tests/lbpm_cell_simulator.cpp b/tests/lbpm_cell_simulator.cpp index caeef89a..304c9cad 100644 --- a/tests/lbpm_cell_simulator.cpp +++ b/tests/lbpm_cell_simulator.cpp @@ -89,7 +89,6 @@ int main(int argc, char **argv) PoissonSolver.Create(); PoissonSolver.Initialize(Study.time_conv_max); - int timestep=0; while (timestep < Study.timestepMax){ From f270604e6b2fa163c8bfe1edd85a9dd94ba2b3c3 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 7 Mar 2022 15:36:16 -0500 Subject: [PATCH 15/37] add membrane transport function for d3q7 --- analysis/Membrane.cpp | 158 ++++++++++++++++++++++-------------------- analysis/Membrane.h | 1 + common/ScaLBL.h | 5 ++ cpu/Ion.cpp | 19 +++++ 4 files changed, 106 insertions(+), 77 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index e90709ec..859c35bf 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -522,13 +522,18 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra } /* allocate memory */ + membraneTag = new int [mlink]; membraneLinks = new int [2*mlink]; membraneDist = new double [2*mlink]; - membraneTag = new int [mlink]; + membraneCoef = new double [2*mlink]; /* construct the membrane*/ + /* * + * Sites inside the membrane (negative distance) -- store at 2*mlink + * Sites outside the membrane (positive distance) -- store at 2*mlink+1 + */ mlink = 0; - int insideMem = 0; int outsideMem = 0; + int localSite = 0; int neighborSite = 0; for (k=1;k Dm, DoubleArray &Distance, IntArra dist=Distance(i+1,j,k); if (dist*locdist < 0.0){ if (locdist < 0.0){ - insideMem = 2*mlink; - outsideMem = 2*mlink+1; + localSite = 2*mlink; + neighborSite = 2*mlink+1; } else{ - insideMem = 2*mlink+1; - outsideMem = 2*mlink; + localSite = 2*mlink+1; + neighborSite = 2*mlink; } - membraneLinks[insideMem] = idx + 1*Np; - membraneLinks[outsideMem] = neighbor + 2*Np; - membraneDist[insideMem] = locdist; - membraneDist[outsideMem] = dist; + membraneLinks[localSite] = idx + 1*Np; + membraneLinks[neighborSite] = neighbor + 2*Np; + membraneDist[localSite] = locdist; + membraneDist[neighborSite] = dist; mlink++; } @@ -559,17 +564,17 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra dist=Distance(i,j+1,k); if (dist*locdist < 0.0){ if (locdist < 0.0){ - insideMem = 2*mlink; - outsideMem = 2*mlink+1; + localSite = 2*mlink; + neighborSite = 2*mlink+1; } else{ - insideMem = 2*mlink+1; - outsideMem = 2*mlink; + localSite = 2*mlink+1; + neighborSite = 2*mlink; } - membraneLinks[insideMem] = idx + 3*Np; - membraneLinks[outsideMem] = neighbor + 4*Np; - membraneDist[insideMem] = locdist; - membraneDist[outsideMem] = dist; + membraneLinks[localSite] = idx + 3*Np; + membraneLinks[neighborSite] = neighbor + 4*Np; + membraneDist[localSite] = locdist; + membraneDist[neighborSite] = dist; mlink++; } @@ -577,17 +582,17 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra dist=Distance(i,j,k+1); if (dist*locdist < 0.0){ if (locdist < 0.0){ - insideMem = 2*mlink; - outsideMem = 2*mlink+1; + localSite = 2*mlink; + neighborSite = 2*mlink+1; } else{ - insideMem = 2*mlink+1; - outsideMem = 2*mlink; + localSite = 2*mlink+1; + neighborSite = 2*mlink; } - membraneLinks[insideMem] = idx + 5*Np; - membraneLinks[outsideMem] = neighbor + 6*Np; - membraneDist[insideMem] = locdist; - membraneDist[outsideMem] = dist; + membraneLinks[localSite] = idx + 5*Np; + membraneLinks[neighborSite] = neighbor + 6*Np; + membraneDist[localSite] = locdist; + membraneDist[neighborSite] = dist; mlink++; } @@ -595,17 +600,17 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra dist=Distance(i+1,j+1,k); if (dist*locdist < 0.0){ if (locdist < 0.0){ - insideMem = 2*mlink; - outsideMem = 2*mlink+1; + localSite = 2*mlink; + neighborSite = 2*mlink+1; } else{ - insideMem = 2*mlink+1; - outsideMem = 2*mlink; + localSite = 2*mlink+1; + neighborSite = 2*mlink; } - membraneLinks[insideMem] = idx + 7*Np; - membraneLinks[outsideMem] = neighbor+8*Np; - membraneDist[insideMem] = locdist; - membraneDist[outsideMem] = dist; + membraneLinks[localSite] = idx + 7*Np; + membraneLinks[neighborSite] = neighbor+8*Np; + membraneDist[localSite] = locdist; + membraneDist[neighborSite] = dist; mlink++; } @@ -613,17 +618,17 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra dist=Distance(i+1,j-1,k); if (dist*locdist < 0.0){ if (locdist < 0.0){ - insideMem = 2*mlink; - outsideMem = 2*mlink+1; + localSite = 2*mlink; + neighborSite = 2*mlink+1; } else{ - insideMem = 2*mlink+1; - outsideMem = 2*mlink; + localSite = 2*mlink+1; + neighborSite = 2*mlink; } - membraneLinks[insideMem] = idx + 9*Np; - membraneLinks[outsideMem] = neighbor + 10*Np; - membraneDist[insideMem] = locdist; - membraneDist[outsideMem] = dist; + membraneLinks[localSite] = idx + 9*Np; + membraneLinks[neighborSite] = neighbor + 10*Np; + membraneDist[localSite] = locdist; + membraneDist[neighborSite] = dist; mlink++; } @@ -631,17 +636,17 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra dist=Distance(i+1,j,k+1); if (dist*locdist < 0.0){ if (locdist < 0.0){ - insideMem = 2*mlink; - outsideMem = 2*mlink+1; + localSite = 2*mlink; + neighborSite = 2*mlink+1; } else{ - insideMem = 2*mlink+1; - outsideMem = 2*mlink; + localSite = 2*mlink+1; + neighborSite = 2*mlink; } - membraneLinks[insideMem] = idx + 11*Np; - membraneLinks[outsideMem] = neighbor + 12*Np; - membraneDist[insideMem] = locdist; - membraneDist[outsideMem] = dist; + membraneLinks[localSite] = idx + 11*Np; + membraneLinks[neighborSite] = neighbor + 12*Np; + membraneDist[localSite] = locdist; + membraneDist[neighborSite] = dist; mlink++; } @@ -649,17 +654,17 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra dist=Distance(i+1,j,k-1); if (dist*locdist < 0.0){ if (locdist < 0.0){ - insideMem = 2*mlink; - outsideMem = 2*mlink+1; + localSite = 2*mlink; + neighborSite = 2*mlink+1; } else{ - insideMem = 2*mlink+1; - outsideMem = 2*mlink; + localSite = 2*mlink+1; + neighborSite = 2*mlink; } - membraneLinks[insideMem] = idx + 13*Np; - membraneLinks[outsideMem] = neighbor + 14*Np; - membraneDist[insideMem] = locdist; - membraneDist[outsideMem] = dist; + membraneLinks[localSite] = idx + 13*Np; + membraneLinks[neighborSite] = neighbor + 14*Np; + membraneDist[localSite] = locdist; + membraneDist[neighborSite] = dist; mlink++; } @@ -667,17 +672,17 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra dist=Distance(i,j+1,k+1); if (dist*locdist < 0.0){ if (locdist < 0.0){ - insideMem = 2*mlink; - outsideMem = 2*mlink+1; + localSite = 2*mlink; + neighborSite = 2*mlink+1; } else{ - insideMem = 2*mlink+1; - outsideMem = 2*mlink; + localSite = 2*mlink+1; + neighborSite = 2*mlink; } - membraneLinks[insideMem] = idx + 15*Np; - membraneLinks[outsideMem] =neighbor + 16*Np; - membraneDist[insideMem] = locdist; - membraneDist[outsideMem] = dist; + membraneLinks[localSite] = idx + 15*Np; + membraneLinks[neighborSite] = neighbor + 16*Np; + membraneDist[localSite] = locdist; + membraneDist[neighborSite] = dist; mlink++; } @@ -685,18 +690,17 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra dist=Distance(i,j+1,k-1); if (dist*locdist < 0.0){ if (locdist < 0.0){ - insideMem = 2*mlink; - outsideMem = 2*mlink+1; + localSite = 2*mlink; + neighborSite = 2*mlink+1; } else{ - insideMem = 2*mlink+1; - - outsideMem = 2*mlink; + localSite = 2*mlink+1; + neighborSite = 2*mlink; } - membraneLinks[insideMem] = idx + 17*Np; - membraneLinks[outsideMem] = neighbor + 18*Np; - membraneDist[insideMem] = locdist; - membraneDist[outsideMem] = dist; + membraneLinks[localSite] = idx + 17*Np; + membraneLinks[neighborSite] = neighbor + 18*Np; + membraneDist[localSite] = locdist; + membraneDist[neighborSite] = dist; mlink++; } } @@ -706,10 +710,10 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra /* Re-organize communication based on membrane structure*/ //...Map recieve list for the X face: q=2,8,10,12,14 ................................. - linkCount_X[0]= D3Q19_MapRecv(-1,0,0, Dm->recvList("X"),0,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); + 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[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.................................. diff --git a/analysis/Membrane.h b/analysis/Membrane.h index 1fc9edf5..aa9caa95 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -58,6 +58,7 @@ public: int *membraneLinks; // D3Q19 links that cross membrane int *membraneTag; // label each link in the membrane double *membraneDist; // distance to membrane for each linked site + double *membraneCoef; // mass transport coefficient for the membrane /** * \brief Create a flow adaptor to operate on the LB model diff --git a/common/ScaLBL.h b/common/ScaLBL.h index f10b150d..0e3320aa 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -217,6 +217,11 @@ 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); + + // GREYSCALE MODEL (Single-component) extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *Dist, int Np, double Den); diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index eb5e3ef5..4fe7e593 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -1,5 +1,24 @@ #include +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 Date: Mon, 7 Mar 2022 16:47:19 -0500 Subject: [PATCH 16/37] add membrane unpack function --- common/ScaLBL.h | 1 + cpu/Ion.cpp | 27 +++++++++++++++++++++++++++ 2 files changed, 28 insertions(+) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 0e3320aa..f13a91d3 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -221,6 +221,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int star 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_Unpack(int q, int *list, int start, int count, double *recvbuf, double *dist, int N, int nlinks, double *coef); // GREYSCALE MODEL (Single-component) diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index 4fe7e593..1d2c9e13 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -1,5 +1,32 @@ #include +extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q, int *list, int start, int count, + double *recvbuf, double *dist, int N, int nlinks, 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 + for (link = 0; link < nlinks; link++) { + // pick out links from the end of the list + idx = count - nlinks + link; + // Get the value from the list -- note that n is the index is from the send (non-local) process + n = list[idx]; + // update link based on mass transfer coefficients + if (!(n < 0)){ + aq = coef[2*link]; + ap = coef[2*link+1]; + fq = dist[q * N + n]; + fp = recvbuf[start + idx]; + fqq = (1-aq)*fq+ap*fp; + dist[q * N + n] = fqq; + } + } +} + extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, double *dist, double *Den, int memLinks, int Np){ From abf5823de6adf8ff5d47e3dbea47ef2390424116 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 11 Mar 2022 05:39:34 -0500 Subject: [PATCH 17/37] membrane data structures compiling --- analysis/Membrane.cpp | 112 ++++++++++++++++++++++++++++++++++++++++-- analysis/Membrane.h | 16 +++++- common/ScaLBL.h | 4 +- cpu/Ion.cpp | 36 ++++++++++---- 4 files changed, 154 insertions(+), 14 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index 859c35bf..27878d8b 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -275,6 +275,17 @@ Membrane::Membrane(std::shared_ptr Dm, int *initialNeighborList, int Ns Membrane::~Membrane() { + 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( dvcNeighborList ); + ScaLBL_FreeDeviceMemory( dvcMembraneLinks ); + ScaLBL_FreeDeviceMemory( dvcMembraneCoef ); + ScaLBL_FreeDeviceMemory( sendbuf_x ); ScaLBL_FreeDeviceMemory( sendbuf_X ); ScaLBL_FreeDeviceMemory( sendbuf_y ); @@ -525,7 +536,6 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra membraneTag = new int [mlink]; membraneLinks = new int [2*mlink]; membraneDist = new double [2*mlink]; - membraneCoef = new double [2*mlink]; /* construct the membrane*/ /* * @@ -707,6 +717,17 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra } } } + + /* Create device copies of data structures */ + ScaLBL_AllocateDeviceMemory((void **)&dvcNeighborList, 18*Np*sizeof(int)); + ScaLBL_AllocateDeviceMemory((void **)&dvcMembraneLinks, 2*mlink*sizeof(int)); + ScaLBL_AllocateDeviceMemory((void **)&dvcMembraneCoef, 2*mlink*sizeof(double)); + ScaLBL_AllocateDeviceMemory((void **)&dvcMembraneDistance, 2*mlink*sizeof(double)); + + ScaLBL_CopyToDevice(dvcNeighborList, neighborList, 18*Np*sizeof(int)); + ScaLBL_CopyToDevice(dvcMembraneLinks, membraneLinks, 2*mlink*sizeof(int)); + ScaLBL_CopyToDevice(dvcMembraneDistance, membraneDist, 2*mlink*sizeof(double)); + /* Re-organize communication based on membrane structure*/ //...Map recieve list for the X face: q=2,8,10,12,14 ................................. @@ -792,7 +813,15 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra CommunicationCount = SendCount+RecvCount; //...................................................................................... - + // Allocate membrane coefficient buffers (for d3q7 recv) + ScaLBL_AllocateZeroCopy((void **) &coefficient_x, (recvCount_x - linkCount_x[0])*sizeof(double)); + ScaLBL_AllocateZeroCopy((void **) &coefficient_X, (recvCount_X - linkCount_X[0])*sizeof(double)); + ScaLBL_AllocateZeroCopy((void **) &coefficient_y, (recvCount_y - linkCount_y[0])*sizeof(double)); + ScaLBL_AllocateZeroCopy((void **) &coefficient_Y, (recvCount_Y - linkCount_Y[0])*sizeof(double)); + ScaLBL_AllocateZeroCopy((void **) &coefficient_z, (recvCount_z - linkCount_z[0])*sizeof(double)); + ScaLBL_AllocateZeroCopy((void **) &coefficient_Z, (recvCount_Z - linkCount_Z[0])*sizeof(double)); + //...................................................................................... + return mlink; } @@ -1030,7 +1059,7 @@ void Membrane::RecvD3Q19AA(double *dist){ Membrane_D3Q19_Unpack(15,dvcRecvDist_Z, dvcRecvLinks_Z,3*recvCount_Z,linkCount_Z[3],recvbuf_Z,dist,N); Membrane_D3Q19_Unpack(18,dvcRecvDist_Z, dvcRecvLinks_Z,4*recvCount_Z,linkCount_Z[4],recvbuf_Z,dist,N); //.................................................................................. - //...Pack the xy edge (8)................................ + //...Pack the xy edge (8)............................... Membrane_D3Q19_Unpack(8,dvcRecvDist_xy, dvcRecvLinks_xy,0,recvCount_xy,recvbuf_xy,dist,N); //...Pack the Xy edge (9)................................ Membrane_D3Q19_Unpack(9,dvcRecvDist_Xy, dvcRecvLinks_Xy,0,recvCount_Xy,recvbuf_Xy,dist,N); @@ -1057,6 +1086,83 @@ void Membrane::RecvD3Q19AA(double *dist){ //................................................................................... Lock=false; // unlock the communicator after communications complete //................................................................................... +} + +void Membrane::SendD3Q7AA(double *dist){ + + if (Lock==true){ + ERROR("Membrane Error (SendD3Q7): Membrane communicator is locked -- did you forget to match Send/Recv calls?"); + } + else{ + Lock=true; + } + // assign tag of 37 to D3Q7 communication + sendtag = recvtag = 37; + ScaLBL_DeviceBarrier(); + // Pack the distributions + //...Packing for x face(q=2)................................ + ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,dist,N); + req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x,rank_x,sendtag); + req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X, recvCount_X,rank_X,recvtag); + //...Packing for X face(q=1)................................ + ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,dist,N); + req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, sendCount_X,rank_X,sendtag); + req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, recvCount_x,rank_x,recvtag); + //...Packing for y face(q=4)................................. + ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,dist,N); + req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y, sendCount_y,rank_y,sendtag); + req2[2] = MPI_COMM_SCALBL.Irecv(recvbuf_Y, recvCount_Y,rank_Y,recvtag); + //...Packing for Y face(q=3)................................. + ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,dist,N); + req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, sendCount_Y,rank_Y,sendtag); + req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y, recvCount_y,rank_y,recvtag); + //...Packing for z face(q=6)................................ + ScaLBL_D3Q19_Pack(6,dvcSendList_z,0,sendCount_z,sendbuf_z,dist,N); + req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z, sendCount_z,rank_z,sendtag); + req2[4] = MPI_COMM_SCALBL.Irecv(recvbuf_Z, recvCount_Z,rank_Z,recvtag); + //...Packing for Z face(q=5)................................ + ScaLBL_D3Q19_Pack(5,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,dist,N); + req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, sendCount_Z,rank_Z,sendtag); + req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z, recvCount_z,rank_z,recvtag); +} + +void Membrane::RecvD37AA(double *dist){ + + //................................................................................... + // Wait for completion of D3Q19 communication + MPI_COMM_SCALBL.waitAll(6,req1); + MPI_COMM_SCALBL.waitAll(6,req2); + ScaLBL_DeviceBarrier(); + //................................................................................... + // NOTE: AA Routine writes to opposite + // Unpack the distributions on the device + //................................................................................... + //...Unpacking for x face(q=2)................................ + ScaLBL_D3Q7_Membrane_Unpack(2,dvcRecvDist_x, dvcRecvLinks_x,0,linkCount_x[0],recvCount_x,recvbuf_x,dist,N,coefficient_x); + //................................................................................... + //...Packing for X face(q=1)................................ + ScaLBL_D3Q7_Membrane_Unpack(1,dvcRecvDist_X, dvcRecvLinks_X,0,linkCount_X[0],recvCount_X,recvbuf_X,dist,N,coefficient_X); + //................................................................................... + //...Packing for y face(q=4)................................. + ScaLBL_D3Q7_Membrane_Unpack(4,dvcRecvDist_y, dvcRecvLinks_y,0,linkCount_y[0],recvCount_y,recvbuf_y,dist,N,coefficient_y); + //................................................................................... + //...Packing for Y face(q=3)................................. + ScaLBL_D3Q7_Membrane_Unpack(3,dvcRecvDist_Y, dvcRecvLinks_Y,0,linkCount_Y[0],recvCount_Y,recvbuf_Y,dist,N,coefficient_Y); + //................................................................................... + //...Packing for z face(q=6)................................ + ScaLBL_D3Q7_Membrane_Unpack(6,dvcRecvDist_z, dvcRecvLinks_z,0,linkCount_z[0],recvCount_z,recvbuf_z,dist,N,coefficient_z); + //...Packing for Z face(q=5)................................ + ScaLBL_D3Q7_Membrane_Unpack(5,dvcRecvDist_Z, dvcRecvLinks_Z,0,linkCount_Z[0],recvCount_Z,recvbuf_Z,dist,N,coefficient_Z); + //.................................................................................. + + //................................................................................... + Lock=false; // unlock the communicator after communications complete + //................................................................................... } +void Membrane::AssignCoefficients(double *dvcPsi, double *dvcDistance, double *dvcMap, std::string method){ + /* Assign mass transfer coefficients to the membrane data structure */ + + +} diff --git a/analysis/Membrane.h b/analysis/Membrane.h index aa9caa95..3074dd46 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -58,8 +58,15 @@ public: int *membraneLinks; // D3Q19 links that cross membrane int *membraneTag; // label each link in the membrane double *membraneDist; // distance to membrane for each linked site - double *membraneCoef; // mass transport coefficient for the membrane + /* + * Device data structures + */ + int *dvcNeighborList; + double *dvcMembraneLinks; + double *dvcMembraneCoef; // mass transport coefficient for the membrane + double *dvcMembraneDistance; + /** * \brief Create a flow adaptor to operate on the LB model * @param ScaLBL - originating data structures @@ -83,6 +90,9 @@ public: void SendD3Q19AA(double *dist); void RecvD3Q19AA(double *dist); + void SendD3Q7AA(double *dist); + void RecvD37AA(double *dist); + void AssignCoefficients(double *dvcPsi, double *dvcDistance, double *dvcMap, 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; @@ -162,5 +172,9 @@ private: 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.h b/common/ScaLBL.h index f13a91d3..ebd5b3ad 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -221,7 +221,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int star 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_Unpack(int q, int *list, int start, int count, double *recvbuf, double *dist, int N, int nlinks, double *coef); +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) diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index 1d2c9e13..a07361bd 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -1,7 +1,13 @@ #include -extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q, int *list, int start, int count, - double *recvbuf, double *dist, int N, int nlinks, double *coef) { +extern "C" void ScaLBL_D3Q7_AssignLinkCoef(){ + +} + + +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 @@ -9,16 +15,28 @@ extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q, int *list, int start, int cou // dist may be even or odd distributions stored by stream layout //.................................................................................... int n, idx, link; - double fq fp,fqq,ap,aq; // coefficient + double fq,fp,fqq,ap,aq; // coefficient + /* First unpack the regular links */ for (link = 0; link < nlinks; link++) { - // pick out links from the end of the list - idx = count - nlinks + link; - // Get the value from the list -- note that n is the index is from the send (non-local) process - n = list[idx]; + // 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; + } + } + /* 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]; - ap = coef[2*link+1]; + 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; From 12026f54d437c21adbee2117adc3d8f239b3a2a3 Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 14 Mar 2022 16:16:47 -0400 Subject: [PATCH 18/37] update to membrane capability --- analysis/Membrane.cpp | 34 ++++++++++++++++- analysis/Membrane.h | 2 +- common/ScaLBL.h | 11 ++++++ cpu/Ion.cpp | 86 ++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 130 insertions(+), 3 deletions(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index 27878d8b..b67f191c 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -1161,8 +1161,40 @@ void Membrane::RecvD37AA(double *dist){ } -void Membrane::AssignCoefficients(double *dvcPsi, double *dvcDistance, double *dvcMap, std::string method){ +void Membrane::AssignCoefficients(int *Map, double *Psi, double *Distance, std::string method){ /* Assign mass transfer coefficients to the membrane data structure */ + double Threshold; + double MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut; + + ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(-1,0,0,Map,Distance,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,Distance,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,Distance,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,Distance,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,Distance,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,Distance,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/analysis/Membrane.h b/analysis/Membrane.h index 3074dd46..558122eb 100644 --- a/analysis/Membrane.h +++ b/analysis/Membrane.h @@ -92,7 +92,7 @@ public: void RecvD3Q19AA(double *dist); void SendD3Q7AA(double *dist); void RecvD37AA(double *dist); - void AssignCoefficients(double *dvcPsi, double *dvcDistance, double *dvcMap, std::string method); + void AssignCoefficients(int *Map, double *Psi, double *Distance, 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; diff --git a/common/ScaLBL.h b/common/ScaLBL.h index ebd5b3ad..8883054f 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -221,6 +221,17 @@ extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int star 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); diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index a07361bd..f8a9f79d 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -1,7 +1,91 @@ #include -extern "C" void ScaLBL_D3Q7_AssignLinkCoef(){ +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]; + // 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; + } } From 1abb9adea68cc7731ab5aaf41ebb2e7518713e70 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 18 Mar 2022 11:05:50 -0400 Subject: [PATCH 19/37] clean up relabel --- analysis/Membrane.cpp | 14 +++++++++++++- models/ColorModel.cpp | 4 ++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/analysis/Membrane.cpp b/analysis/Membrane.cpp index b67f191c..b0125ec7 100644 --- a/analysis/Membrane.cpp +++ b/analysis/Membrane.cpp @@ -1161,12 +1161,24 @@ void Membrane::RecvD37AA(double *dist){ } -void Membrane::AssignCoefficients(int *Map, double *Psi, double *Distance, std::string method){ +// std::shared_ptr db){ +void Membrane::AssignCoefficients(int *Map, double *Psi, double *Distance, string method){ +} /* Assign mass transfer coefficients to the membrane data structure */ double Threshold; double MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut; + membrane_db = db->getDatabase("Membrane"); + + if (method == "Voltage Gated Potassium"){ + MassFractionIn = 0.0; + MassFractionOut = 0.0; + MassFractionIn = 0.0; + ThresholdMassFractionOut = 0.0; + ThresholdMassFractionIn = -55.0; + } + ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(-1,0,0,Map,Distance,Psi,Threshold, MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut, dvcRecvDist_X,dvcRecvLinks_X,coefficient_X,0,linkCount_X[0],recvCount_X, diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 21315063..9e253f09 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -388,6 +388,10 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase) { AFFINITY, volume_fraction); } } + + // clean up + delete [] label_count; + delete [] label_count_global; } void ScaLBL_ColorModel::Create() { From e8d0b0b48a962f0e63214c6686d7cb96a8959167 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 18 Mar 2022 11:06:16 -0400 Subject: [PATCH 20/37] adding membrane functions --- models/IonModel.cpp | 64 +++++++++++++++++++++++++++++++++++++++++++++ models/IonModel.h | 11 +++++++- 2 files changed, 74 insertions(+), 1 deletion(-) diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 7213cfce..9ab02d1e 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -583,6 +583,70 @@ void ScaLBL_IonModel::SetDomain() { nprocz = Dm->nprocz(); } +void ScaLBL_IonModel::SetMembrane() { + membrane_db = db->getDatabase("Membrane"); + + /* set distance based on labels inside the membrane (all other labels will be outside) */ + auto InsideMembraneLabels = membrane_db->getVector("InsideMembraneLabels"); + 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; mrank()); diff --git a/models/IonModel.h b/models/IonModel.h index 8930b1df..5fcaa982 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, @@ -98,6 +100,13 @@ public: 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; From bfa6f5e5b81c4277d7b90b092dfdd5e81ec1dda0 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 18 Mar 2022 11:09:34 -0400 Subject: [PATCH 21/37] move membrane to common folder --- {analysis => common}/Membrane.cpp | 0 {analysis => common}/Membrane.h | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename {analysis => common}/Membrane.cpp (100%) rename {analysis => common}/Membrane.h (100%) diff --git a/analysis/Membrane.cpp b/common/Membrane.cpp similarity index 100% rename from analysis/Membrane.cpp rename to common/Membrane.cpp diff --git a/analysis/Membrane.h b/common/Membrane.h similarity index 100% rename from analysis/Membrane.h rename to common/Membrane.h From bfe1de6be2cb9a9b25921befbab3de6b466934a2 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 18 Mar 2022 11:20:09 -0400 Subject: [PATCH 22/37] membrane structure in IonModel --- common/Membrane.cpp | 5 +---- models/IonModel.cpp | 11 ++++++----- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/common/Membrane.cpp b/common/Membrane.cpp index b0125ec7..7959b73f 100644 --- a/common/Membrane.cpp +++ b/common/Membrane.cpp @@ -1,6 +1,6 @@ /* Flow adaptor class for multiphase flow methods */ -#include "analysis/Membrane.h" +#include "common/Membrane.h" #include "analysis/distance.h" Membrane::Membrane(std::shared_ptr Dm, int *initialNeighborList, int Nsites) { @@ -1163,14 +1163,11 @@ void Membrane::RecvD37AA(double *dist){ // std::shared_ptr db){ void Membrane::AssignCoefficients(int *Map, double *Psi, double *Distance, string method){ -} /* Assign mass transfer coefficients to the membrane data structure */ double Threshold; double MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut; - membrane_db = db->getDatabase("Membrane"); - if (method == "Voltage Gated Potassium"){ MassFractionIn = 0.0; MassFractionOut = 0.0; diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 9ab02d1e..4484e8f3 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -584,10 +584,12 @@ void ScaLBL_IonModel::SetDomain() { } 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 InsideMembraneLabels = membrane_db->getVector("InsideMembraneLabels"); + auto MembraneLabels = membrane_db->getVector("MembraneLabels"); IonMembrane = std::shared_ptr(new Membrane(Dm, NeighborList, Np)); signed char LABEL = 0; @@ -608,7 +610,7 @@ void ScaLBL_IonModel::SetMembrane() { membrane_id(i,j,k) = 1; // default value LABEL = Dm->id[k*Nx*Ny + j*Nx + i]; for (size_t m=0; mCreate(Dm, MembraneDistance, Map); // clean up delete [] label_count; From c94d4e51942331d7b984198109b4b12a3525d464 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 18 Mar 2022 11:23:01 -0400 Subject: [PATCH 23/37] membrane structure in IonModel --- tests/TestMembrane.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/TestMembrane.cpp b/tests/TestMembrane.cpp index 048be117..6b3953b8 100644 --- a/tests/TestMembrane.cpp +++ b/tests/TestMembrane.cpp @@ -7,7 +7,7 @@ #include #include #include "common/MPI.h" -#include "analysis/Membrane.h" +#include "common/Membrane.h" using namespace std; From 702eaae1c1ab81d40dee79610ddb103edf1a5d6d Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 18 Mar 2022 18:08:44 -0400 Subject: [PATCH 24/37] try at membrane simulator --- common/Membrane.cpp | 55 +++++---- common/Membrane.h | 14 ++- models/IonModel.cpp | 217 ++++++++++++++++++++++++++++++++++ tests/TestMembrane.cpp | 17 +-- tests/lbpm_cell_simulator.cpp | 3 +- 5 files changed, 270 insertions(+), 36 deletions(-) diff --git a/common/Membrane.cpp b/common/Membrane.cpp index 7959b73f..be6da668 100644 --- a/common/Membrane.cpp +++ b/common/Membrane.cpp @@ -3,21 +3,21 @@ #include "common/Membrane.h" #include "analysis/distance.h" -Membrane::Membrane(std::shared_ptr Dm, int *initialNeighborList, int Nsites) { +Membrane::Membrane(std::shared_ptr Dm, int *dvcNeighborList, int Nsites) { + + Np = Nsites; + int * 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(); - Np = Nsites; - neighborList = new int[18*Np]; - /* Copy neighborList */ - for (int idx=0; idxComm.barrier(); + ScaLBL_CopyToDevice(NeighborList, initialNeighborList, 18*Np*sizeof(int)); + /* Copy communication lists */ //...................................................................................... //Lock=false; // unlock the communicator @@ -275,6 +275,11 @@ Membrane::Membrane(std::shared_ptr Dm, int *initialNeighborList, int Ns Membrane::~Membrane() { + delete [] initialNeighborList; + delete [] membraneLinks; + delete [] membraneTag; + delete [] membraneDist; + ScaLBL_FreeDeviceMemory( coefficient_x ); ScaLBL_FreeDeviceMemory( coefficient_X ); ScaLBL_FreeDeviceMemory( coefficient_y ); @@ -282,9 +287,9 @@ Membrane::~Membrane() { ScaLBL_FreeDeviceMemory( coefficient_z ); ScaLBL_FreeDeviceMemory( coefficient_Z ); - ScaLBL_FreeDeviceMemory( dvcNeighborList ); - ScaLBL_FreeDeviceMemory( dvcMembraneLinks ); - ScaLBL_FreeDeviceMemory( dvcMembraneCoef ); + ScaLBL_FreeDeviceMemory( NeighborList ); + ScaLBL_FreeDeviceMemory( MembraneLinks ); + ScaLBL_FreeDeviceMemory( MembraneCoef ); ScaLBL_FreeDeviceMemory( sendbuf_x ); ScaLBL_FreeDeviceMemory( sendbuf_X ); @@ -401,6 +406,15 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra int i,j,k; int idx, neighbor; double dist, locdist; + + int * neighborList = new int[18*Np]; + /* Copy neighborList */ + for (int idx=0; idx Dm, DoubleArray &Distance, IntArra } } } - + /* allocate memory */ membraneTag = new int [mlink]; membraneLinks = new int [2*mlink]; @@ -719,14 +733,13 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra } /* Create device copies of data structures */ - ScaLBL_AllocateDeviceMemory((void **)&dvcNeighborList, 18*Np*sizeof(int)); - ScaLBL_AllocateDeviceMemory((void **)&dvcMembraneLinks, 2*mlink*sizeof(int)); - ScaLBL_AllocateDeviceMemory((void **)&dvcMembraneCoef, 2*mlink*sizeof(double)); - ScaLBL_AllocateDeviceMemory((void **)&dvcMembraneDistance, 2*mlink*sizeof(double)); + ScaLBL_AllocateDeviceMemory((void **)&MembraneLinks, 2*mlink*sizeof(int)); + ScaLBL_AllocateDeviceMemory((void **)&MembraneCoef, 2*mlink*sizeof(double)); + ScaLBL_AllocateDeviceMemory((void **)&MembraneDistance, 2*mlink*sizeof(double)); - ScaLBL_CopyToDevice(dvcNeighborList, neighborList, 18*Np*sizeof(int)); - ScaLBL_CopyToDevice(dvcMembraneLinks, membraneLinks, 2*mlink*sizeof(int)); - ScaLBL_CopyToDevice(dvcMembraneDistance, membraneDist, 2*mlink*sizeof(double)); + ScaLBL_CopyToDevice(NeighborList, neighborList, 18*Np*sizeof(int)); + ScaLBL_CopyToDevice(MembraneLinks, membraneLinks, 2*mlink*sizeof(int)); + ScaLBL_CopyToDevice(MembraneDistance, membraneDist, 2*mlink*sizeof(double)); /* Re-organize communication based on membrane structure*/ @@ -1126,7 +1139,7 @@ void Membrane::SendD3Q7AA(double *dist){ req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z, recvCount_z,rank_z,recvtag); } -void Membrane::RecvD37AA(double *dist){ +void Membrane::RecvD3Q7AA(double *dist){ //................................................................................... // Wait for completion of D3Q19 communication diff --git a/common/Membrane.h b/common/Membrane.h index 558122eb..81d03bf7 100644 --- a/common/Membrane.h +++ b/common/Membrane.h @@ -54,7 +54,10 @@ public: int Np; int Nx,Ny,Nz,N; - int *neighborList; // modified neighborlist + 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 @@ -62,10 +65,9 @@ public: /* * Device data structures */ - int *dvcNeighborList; - double *dvcMembraneLinks; - double *dvcMembraneCoef; // mass transport coefficient for the membrane - double *dvcMembraneDistance; + double *MembraneLinks; + double *MembraneCoef; // mass transport coefficient for the membrane + double *MembraneDistance; /** * \brief Create a flow adaptor to operate on the LB model @@ -91,7 +93,7 @@ public: void SendD3Q19AA(double *dist); void RecvD3Q19AA(double *dist); void SendD3Q7AA(double *dist); - void RecvD37AA(double *dist); + void RecvD3Q7AA(double *dist); void AssignCoefficients(int *Map, double *Psi, double *Distance, std::string method); //...................................................................................... // Buffers to store data sent and recieved by this MPI process diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 4484e8f3..a73db663 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -1246,6 +1246,223 @@ 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(); + + 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], + &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/tests/TestMembrane.cpp b/tests/TestMembrane.cpp index 6b3953b8..fc87b48f 100644 --- a/tests/TestMembrane.cpp +++ b/tests/TestMembrane.cpp @@ -95,6 +95,14 @@ int main(int argc, char **argv) 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(); @@ -120,7 +128,7 @@ int main(int argc, char **argv) } /* create a membrane data structure */ - Membrane M(Dm, neighborList, Np); + Membrane M(Dm, NeighborList, Np); int MembraneCount = M.Create(Dm, Distance, Map); if (rank==0) printf (" Number of membrane links: %i \n", MembraneCount); @@ -159,13 +167,6 @@ int main(int argc, char **argv) if (argc > 1) Dm->AggregateLabels("membrane.raw"); - //......................device distributions................................. - int *NeighborList; - int *dvcMap; - //........................................................................... - ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); - ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Npad); - //........................................................................... // Update GPU data structures if (rank==0) printf ("Setting up device map and neighbor list \n"); diff --git a/tests/lbpm_cell_simulator.cpp b/tests/lbpm_cell_simulator.cpp index 304c9cad..2ac647d5 100644 --- a/tests/lbpm_cell_simulator.cpp +++ b/tests/lbpm_cell_simulator.cpp @@ -69,6 +69,7 @@ int main(int argc, char **argv) IonModel.SetDomain(); IonModel.ReadInput(); IonModel.Create(); + IonModel.SetMembrane(); // Create analysis object ElectroChemistryAnalyzer Analysis(IonModel.Dm); @@ -95,7 +96,7 @@ int main(int argc, char **argv) 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.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential + IonModel.RunMembrane(StokesModel.Velocity,PoissonSolver.ElectricField,PoissonSolver.Psi); //solve for ion transport with membrane timestep++;//AA operations From a50d9e9aa64f83422fb7c58d51d71cca7ab98093 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sun, 20 Mar 2022 09:29:42 -0400 Subject: [PATCH 25/37] add python script to generate bubble --- example/Bubble/CreateBubble.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 example/Bubble/CreateBubble.py diff --git a/example/Bubble/CreateBubble.py b/example/Bubble/CreateBubble.py new file mode 100644 index 00000000..c58547c2 --- /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") + From 9e3a07d4192d2fa3184de2c6325e445937798ea2 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sun, 20 Mar 2022 09:30:46 -0400 Subject: [PATCH 26/37] add python script to generate bubble --- example/Bubble/CreateBubble.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/Bubble/CreateBubble.py b/example/Bubble/CreateBubble.py index c58547c2..76fbefab 100644 --- a/example/Bubble/CreateBubble.py +++ b/example/Bubble/CreateBubble.py @@ -11,7 +11,7 @@ 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 : + if (dist < 12.5 ) : D[i,j,k] = 2 D.tofile("bubble_40x40x40.raw") From 51c88f0055d93aa1258170bfbca6edebdd65dd57 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sun, 20 Mar 2022 11:22:46 -0400 Subject: [PATCH 27/37] cell simulator runs --- common/Membrane.cpp | 20 +++++++++++++++----- models/IonModel.cpp | 4 ++-- 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/common/Membrane.cpp b/common/Membrane.cpp index be6da668..4435662c 100644 --- a/common/Membrane.cpp +++ b/common/Membrane.cpp @@ -6,10 +6,9 @@ Membrane::Membrane(std::shared_ptr Dm, int *dvcNeighborList, int Nsites) { Np = Nsites; - int * initialNeighborList = new int[18*Np]; - ScaLBL_AllocateDeviceMemory((void **)&NeighborList, 18*Np*sizeof(int)); - - Lock=false; // unlock the communicator + 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(); @@ -85,7 +84,12 @@ Membrane::Membrane(std::shared_ptr Dm, int *dvcNeighborList, int Nsites recvCount_Xz=Dm->recvCount("Xz"); recvCount_XY=Dm->recvCount("XY"); recvCount_YZ=Dm->recvCount("YZ"); - recvCount_XZ=Dm->recvCount("XZ");\ + 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); @@ -407,6 +411,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra 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; idx Dm, DoubleArray &Distance, IntArra /* go through the neighborlist structure */ /* count & cut the links */ + if (rank == 0) printf(" Cut membrane links... \n"); for (k=1;k Dm, DoubleArray &Distance, IntArra * Sites inside the membrane (negative distance) -- store at 2*mlink * Sites outside the membrane (positive distance) -- store at 2*mlink+1 */ + if (rank == 0) printf(" Construct membrane data structures... \n"); mlink = 0; int localSite = 0; int neighborSite = 0; for (k=1;k Dm, DoubleArray &Distance, IntArra } } } + + if (rank == 0) printf(" Create device data structures... \n"); /* Create device copies of data structures */ ScaLBL_AllocateDeviceMemory((void **)&MembraneLinks, 2*mlink*sizeof(int)); @@ -742,6 +751,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra ScaLBL_CopyToDevice(MembraneDistance, membraneDist, 2*mlink*sizeof(double)); + if (rank == 0) printf(" Construct communication data structures... \n"); /* Re-organize communication based on membrane structure*/ //...Map recieve list for the X face: q=2,8,10,12,14 ................................. linkCount_X[0] = D3Q19_MapRecv(-1,0,0, Dm->recvList("X"),0,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); diff --git a/models/IonModel.cpp b/models/IonModel.cpp index a73db663..cc716a5a 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -623,12 +623,12 @@ void ScaLBL_IonModel::SetMembrane() { label_count_global[m] = Dm->Comm.sumReduce(label_count[m]); } if (rank == 0) { - printf("Membrane labels: %lu \n", MembraneLabels.size()); + printf(" Membrane labels: %lu \n", MembraneLabels.size()); for (size_t m=0; m Date: Sun, 20 Mar 2022 13:32:24 -0400 Subject: [PATCH 28/37] read input files --- models/IonModel.cpp | 3 ++- models/IonModel.h | 3 +-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/models/IonModel.cpp b/models/IonModel.cpp index cc716a5a..6f3a3d18 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -590,6 +590,7 @@ void ScaLBL_IonModel::SetMembrane() { /* 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; @@ -890,7 +891,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++) { diff --git a/models/IonModel.h b/models/IonModel.h index 5fcaa982..824dd000 100644 --- a/models/IonModel.h +++ b/models/IonModel.h @@ -99,14 +99,13 @@ 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; From 4927e5470781e82edd03efcaa783b21d97907a8c Mon Sep 17 00:00:00 2001 From: James McClure Date: Sun, 20 Mar 2022 15:13:53 -0400 Subject: [PATCH 29/37] add single cell example --- common/Membrane.cpp | 36 +++++++++++----- common/Membrane.h | 5 ++- example/Bubble/CreateCell.py | 36 ++++++++++++++++ example/Bubble/cell.db | 79 ++++++++++++++++++++++++++++++++++++ models/IonModel.cpp | 38 ++++++++++++++++- models/IonModel.h | 1 + 6 files changed, 180 insertions(+), 15 deletions(-) create mode 100644 example/Bubble/CreateCell.py create mode 100644 example/Bubble/cell.db diff --git a/common/Membrane.cpp b/common/Membrane.cpp index 4435662c..ab8afd73 100644 --- a/common/Membrane.cpp +++ b/common/Membrane.cpp @@ -556,6 +556,9 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra membraneTag = new int [mlink]; membraneLinks = new int [2*mlink]; membraneDist = new double [2*mlink]; + membraneLinkCount = mlink; + + if (rank == 0) printf(" (cut %i links crossing membrane) \n",mlink); /* construct the membrane*/ /* * @@ -744,11 +747,13 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra /* Create device copies of data structures */ ScaLBL_AllocateDeviceMemory((void **)&MembraneLinks, 2*mlink*sizeof(int)); ScaLBL_AllocateDeviceMemory((void **)&MembraneCoef, 2*mlink*sizeof(double)); - ScaLBL_AllocateDeviceMemory((void **)&MembraneDistance, 2*mlink*sizeof(double)); + //ScaLBL_AllocateDeviceMemory((void **)&MembraneDistance, 2*mlink*sizeof(double)); + ScaLBL_AllocateDeviceMemory((void **)&MembraneDistance, Nx*Ny*Nz*sizeof(double)); ScaLBL_CopyToDevice(NeighborList, neighborList, 18*Np*sizeof(int)); ScaLBL_CopyToDevice(MembraneLinks, membraneLinks, 2*mlink*sizeof(int)); - ScaLBL_CopyToDevice(MembraneDistance, membraneDist, 2*mlink*sizeof(double)); + //ScaLBL_CopyToDevice(MembraneDistance, membraneDist, 2*mlink*sizeof(double)); + ScaLBL_CopyToDevice(MembraneDistance, Distance.data(), Nx*Ny*Nz*sizeof(double)); if (rank == 0) printf(" Construct communication data structures... \n"); @@ -1185,46 +1190,55 @@ void Membrane::RecvD3Q7AA(double *dist){ } // std::shared_ptr db){ -void Membrane::AssignCoefficients(int *Map, double *Psi, double *Distance, string method){ +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; - MassFractionIn = 0.0; ThresholdMassFractionOut = 0.0; - ThresholdMassFractionIn = -55.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,Distance,Psi,Threshold, + 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,Distance,Psi,Threshold, + 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,Distance,Psi,Threshold, + 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,Distance,Psi,Threshold, + 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,Distance,Psi,Threshold, + 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,Distance,Psi,Threshold, + 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 index 81d03bf7..a153b8d6 100644 --- a/common/Membrane.h +++ b/common/Membrane.h @@ -53,6 +53,7 @@ class Membrane { public: int Np; int Nx,Ny,Nz,N; + int membraneLinkCount; int *initialNeighborList; // original neighborlist int *NeighborList; // modified neighborlist @@ -65,7 +66,7 @@ public: /* * Device data structures */ - double *MembraneLinks; + int *MembraneLinks; double *MembraneCoef; // mass transport coefficient for the membrane double *MembraneDistance; @@ -94,7 +95,7 @@ public: void RecvD3Q19AA(double *dist); void SendD3Q7AA(double *dist); void RecvD3Q7AA(double *dist); - void AssignCoefficients(int *Map, double *Psi, double *Distance, std::string method); + 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; diff --git a/example/Bubble/CreateCell.py b/example/Bubble/CreateCell.py new file mode 100644 index 00000000..c267f65b --- /dev/null +++ b/example/Bubble/CreateCell.py @@ -0,0 +1,36 @@ +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") + +for i in range(0, 40): + for j in range (0, 40): + for k in range (0,40): + C1[i,j,k] = 4.0e-6 + C2[i,j,k] = 150.0e-6 + dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz)) + if (dist < 15.5 ) : + C1[i,j,k] = 140.0e-6 + C2[i,j,k] = 10.0e-6 + + +C1.tofile("cell_concentration_K_40x40x40.raw") +C2.tofile("cell_concentration_Na_40x40x40.raw") + diff --git a/example/Bubble/cell.db b/example/Bubble/cell.db new file mode 100644 index 00000000..6581f682 --- /dev/null +++ b/example/Bubble/cell.db @@ -0,0 +1,79 @@ +MultiphysController { + timestepMax = 20 + 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" + temperature = 293.15 //unit [K] + number_ion_species = 2 //number of ions + tauList = 1.0, 1.0 + IonDiffusivityList = 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec] + IonValenceList = 1, 1 //valence charge of ions; dimensionless; positive/negative integer + IonConcentrationList = 1.0e-6, 1.0e-6 //user-input unit: [mol/m^3] + BC_InletList = 0, 0 //boundary condition for inlet; 0=periodic; 1=ion concentration; 2=ion flux + BC_OutletList = 0, 0 //boundary condition for outlet; 0=periodic; 1=ion concentration; 2=ion flux + InletValueList = 0, 0 //if ion concentration unit=[mol/m^3]; if flux (inward) unit=[mol/m^2/sec] + OutletValueList = 0, 0 //if ion concentration unit=[mol/m^3]; if flux (inward) unit=[mol/m^2/sec] + 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/models/IonModel.cpp b/models/IonModel.cpp index 6f3a3d18..9fdabe61 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -629,7 +629,7 @@ void ScaLBL_IonModel::SetMembrane() { LABEL = MembraneLabels[m]; double volume_fraction = double(label_count_global[m]) / double((Nx - 2) * (Ny - 2) * (Nz - 2) * nprocs); - printf(" label=%d, volume fraction==%f\n", LABEL, volume_fraction); + printf(" label=%d, volume fraction = %f\n", LABEL, volume_fraction); } } /* signed distance to the membrane ( - inside / + outside) */ @@ -844,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, @@ -860,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(); @@ -1263,7 +1295,9 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl //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]) { diff --git a/models/IonModel.h b/models/IonModel.h index 824dd000..5b899c6c 100644 --- a/models/IonModel.h +++ b/models/IonModel.h @@ -90,6 +90,7 @@ public: IntArray Map; DoubleArray Distance; int *NeighborList; + int *dvcMap; double *fq; double *Ci; double *ChargeDensity; From 286e779459538a47a961c9c7af45ca5c72ce19d7 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sun, 20 Mar 2022 18:01:19 -0400 Subject: [PATCH 30/37] refining cell example --- example/Bubble/CreateCell.py | 51 ++++++++++++++++++++++++++++++++---- example/Bubble/cell.db | 18 +++++-------- 2 files changed, 53 insertions(+), 16 deletions(-) diff --git a/example/Bubble/CreateCell.py b/example/Bubble/CreateCell.py index c267f65b..ac62863f 100644 --- a/example/Bubble/CreateCell.py +++ b/example/Bubble/CreateCell.py @@ -19,18 +19,59 @@ 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): - C1[i,j,k] = 4.0e-6 - C2[i,j,k] = 150.0e-6 + #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] = 140.0e-6 - C2[i,j,k] = 10.0e-6 + 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 index 6581f682..dbefcb2d 100644 --- a/example/Bubble/cell.db +++ b/example/Bubble/cell.db @@ -1,5 +1,5 @@ MultiphysController { - timestepMax = 20 + timestepMax = 40 num_iter_Ion_List = 2 analysis_interval = 50 tolerance = 1.0e-9 @@ -13,17 +13,13 @@ Stokes { 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" + 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 = 2 //number of ions - tauList = 1.0, 1.0 - IonDiffusivityList = 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec] - IonValenceList = 1, 1 //valence charge of ions; dimensionless; positive/negative integer - IonConcentrationList = 1.0e-6, 1.0e-6 //user-input unit: [mol/m^3] - BC_InletList = 0, 0 //boundary condition for inlet; 0=periodic; 1=ion concentration; 2=ion flux - BC_OutletList = 0, 0 //boundary condition for outlet; 0=periodic; 1=ion concentration; 2=ion flux - InletValueList = 0, 0 //if ion concentration unit=[mol/m^3]; if flux (inward) unit=[mol/m^2/sec] - OutletValueList = 0, 0 //if ion concentration unit=[mol/m^3]; if flux (inward) unit=[mol/m^2/sec] + 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 From 0a1057926d534acb5f4ca4b875f5442d1650549b Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 21 Mar 2022 19:44:04 -0400 Subject: [PATCH 31/37] start on cuda function --- cuda/Ion.cu | 201 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 201 insertions(+) diff --git a/cuda/Ion.cu b/cuda/Ion.cu index 49d0b80a..8db0ae05 100644 --- a/cuda/Ion.cu +++ b/cuda/Ion.cu @@ -5,6 +5,207 @@ #define NBLOCKS 1024 #define NTHREADS 256 + +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) { +} + +extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, + double *dist, double *Den, int memLinks, int Np){ +} + +__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 Date: Mon, 21 Mar 2022 19:44:21 -0400 Subject: [PATCH 32/37] werkin --- common/Membrane.cpp | 100 ++++++++++++++++++++++++----------------- cpu/Ion.cpp | 7 ++- tests/TestMembrane.cpp | 75 ++++++++++++++++++------------- 3 files changed, 107 insertions(+), 75 deletions(-) diff --git a/common/Membrane.cpp b/common/Membrane.cpp index ab8afd73..39171944 100644 --- a/common/Membrane.cpp +++ b/common/Membrane.cpp @@ -433,46 +433,46 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i-1,j,k); dist=Distance(i-1,j,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[idx]=idx + 2*Np; } neighbor=Map(i+1,j,k); dist=Distance(i+1,j,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[Np+idx] = idx + 1*Np; mlink++; } neighbor=Map(i,j-1,k); dist=Distance(i,j-1,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[2*Np+idx]=idx + 4*Np; } neighbor=Map(i,j+1,k); dist=Distance(i,j+1,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[3*Np+idx]=idx + 3*Np; mlink++; } neighbor=Map(i,j,k-1); dist=Distance(i,j,k-1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[4*Np+idx]=idx + 6*Np; } neighbor=Map(i,j,k+1); dist=Distance(i,j,k+1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[5*Np+idx]=idx + 5*Np; mlink++; } neighbor=Map(i-1,j-1,k); dist=Distance(i-1,j-1,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[6*Np+idx]=idx + 8*Np; } @@ -485,65 +485,65 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i-1,j+1,k); dist=Distance(i-1,j+1,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[8*Np+idx]=idx + 10*Np; } neighbor=Map(i+1,j-1,k); dist=Distance(i+1,j-1,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[9*Np+idx]=idx + 9*Np; mlink++; } neighbor=Map(i-1,j,k-1); dist=Distance(i-1,j,k-1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[10*Np+idx]=idx + 12*Np; } neighbor=Map(i+1,j,k+1); dist=Distance(i+1,j,k+1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[11*Np+idx]=idx + 11*Np; mlink++; } neighbor=Map(i-1,j,k+1); dist=Distance(i-1,j,k+1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[12*Np+idx]=idx + 14*Np; } neighbor=Map(i+1,j,k-1); dist=Distance(i+1,j,k-1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[13*Np+idx]=idx + 13*Np; mlink++; } neighbor=Map(i,j-1,k-1); dist=Distance(i,j-1,k-1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[14*Np+idx]=idx + 16*Np; } neighbor=Map(i,j+1,k+1); dist=Distance(i,j+1,k+1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[15*Np+idx]=idx + 15*Np; mlink++; } neighbor=Map(i,j-1,k+1); dist=Distance(i,j-1,k+1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[16*Np+idx]=idx + 18*Np; } neighbor=Map(i,j+1,k-1); dist=Distance(i,j+1,k-1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ neighborList[17*Np+idx]=idx + 17*Np; mlink++; } @@ -579,7 +579,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i+1,j,k); dist=Distance(i+1,j,k); if (dist*locdist < 0.0){ - if (locdist < 0.0){ + if (locdist < 0.0 && !(neighbor<0)){ localSite = 2*mlink; neighborSite = 2*mlink+1; } @@ -596,7 +596,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i,j+1,k); dist=Distance(i,j+1,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ if (locdist < 0.0){ localSite = 2*mlink; neighborSite = 2*mlink+1; @@ -614,7 +614,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i,j,k+1); dist=Distance(i,j,k+1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ if (locdist < 0.0){ localSite = 2*mlink; neighborSite = 2*mlink+1; @@ -632,7 +632,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i+1,j+1,k); dist=Distance(i+1,j+1,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ if (locdist < 0.0){ localSite = 2*mlink; neighborSite = 2*mlink+1; @@ -650,7 +650,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i+1,j-1,k); dist=Distance(i+1,j-1,k); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ if (locdist < 0.0){ localSite = 2*mlink; neighborSite = 2*mlink+1; @@ -668,7 +668,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i+1,j,k+1); dist=Distance(i+1,j,k+1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ if (locdist < 0.0){ localSite = 2*mlink; neighborSite = 2*mlink+1; @@ -686,7 +686,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i+1,j,k-1); dist=Distance(i+1,j,k-1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ if (locdist < 0.0){ localSite = 2*mlink; neighborSite = 2*mlink+1; @@ -704,7 +704,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i,j+1,k+1); dist=Distance(i,j+1,k+1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ if (locdist < 0.0){ localSite = 2*mlink; neighborSite = 2*mlink+1; @@ -722,7 +722,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra neighbor=Map(i,j+1,k-1); dist=Distance(i,j+1,k-1); - if (dist*locdist < 0.0){ + if (dist*locdist < 0.0 && !(neighbor<0)){ if (locdist < 0.0){ localSite = 2*mlink; neighborSite = 2*mlink+1; @@ -736,7 +736,7 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra membraneDist[localSite] = locdist; membraneDist[neighborSite] = dist; mlink++; - } + } } } } @@ -824,6 +824,12 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra //...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); //................................................................................... + if (rank == 0) printf(" x count = %i \n",linkCount_x[0]); + if (rank == 0) printf(" X count = %i \n",linkCount_X[0]); + if (rank == 0) printf(" y count = %i \n",linkCount_y[0]); + if (rank == 0) printf(" Y count = %i \n",linkCount_Y[0]); + if (rank == 0) printf(" z count = %i \n",linkCount_z[0]); + if (rank == 0) printf(" Z count = %i \n",linkCount_Z[0]); //...................................................................................... MPI_COMM_SCALBL.barrier(); @@ -881,8 +887,14 @@ int Membrane::D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list, int star nn = Map(i,j,k); distanceNonLocal = Distance(i,j,k); + + printf("CHECK: idx=%i, n=%i, (%i, %i, %i) shift {%i, %i, %i}, stored nn=%i \n",idx,n,i,j,k,Cqx,Cqy,Cqz,nn); + ReturnDist[idx] = nn; - + //if (nn < 0){ + // printf(" Check map for site (%i, %i, %i) based on Cq=(%i, %i, %i) \n", i,j,k,Cqx,Cqy,Cqz); + //} + /* tag the links to swap out later*/ if (distanceLocal*distanceNonLocal < 0.0){ memLinkList[memLinkCount++] = idx; @@ -944,6 +956,7 @@ void Membrane::SendD3Q19AA(double *dist){ req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, 5*sendCount_X,rank_X,sendtag); req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, 5*recvCount_x,rank_x,recvtag); + //...Packing for y face(4,8,9,16,18)................................. ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,dist,N); ScaLBL_D3Q19_Pack(8,dvcSendList_y,sendCount_y,sendCount_y,sendbuf_y,dist,N); @@ -962,6 +975,7 @@ void Membrane::SendD3Q19AA(double *dist){ req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, 5*sendCount_Y,rank_Y,sendtag); req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y, 5*recvCount_y,rank_y,recvtag); + //...Packing for z face(6,12,13,16,17)................................ ScaLBL_D3Q19_Pack(6,dvcSendList_z,0,sendCount_z,sendbuf_z,dist,N); ScaLBL_D3Q19_Pack(12,dvcSendList_z,sendCount_z,sendCount_z,sendbuf_z,dist,N); @@ -1129,29 +1143,32 @@ void Membrane::SendD3Q7AA(double *dist){ ScaLBL_DeviceBarrier(); // Pack the distributions //...Packing for x face(q=2)................................ - ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,dist,N); + ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,dist,Np); req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x,rank_x,sendtag); req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X, recvCount_X,rank_X,recvtag); //...Packing for X face(q=1)................................ - ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,dist,N); + ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,dist,Np); req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, sendCount_X,rank_X,sendtag); req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, recvCount_x,rank_x,recvtag); + //for (int idx=0; idx #include "common/MPI.h" #include "common/Membrane.h" +#include "common/ScaLBL.h" using namespace std; @@ -33,11 +34,6 @@ int main(int argc, char **argv) { int i,j,k,n; - - static int D3Q19[18][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}, - {1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0}, - {1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1}, - {0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}}; int rank = comm.getRank(); if (rank == 0){ @@ -185,8 +181,8 @@ int main(int argc, char **argv) ScaLBL_DeviceBarrier(); // Create a dummy distribution data structure - double *fq; - fq = new double[19*Np]; + double *fq_host; + fq_host = new double[19*Np]; if (rank==0) printf ("Setting up distributions \n"); for (k=1; kLastExterior(), + Np); + DoubleArray Result(Nx,Ny,Nz); +\ + ScaLBL_Comm->RegularLayout(Map, Ci, Result); - for (int idx=ScaLBL_Comm->first_interior; idxlast_interior; idx++){ - n = TmpMap[idx]; - k = n/(Nx*Ny); - j = (n-Nx*Ny*k)/Nx; - i = n-Nx*Ny*k-Nx*j; - for (int q=1; q<19; q++){ - int nn = neighborList[(q-1)*Np+idx]; - double value=fq[nn]; - // 3D index of neighbor - int iq=i-D3Q19[q-1][0]; - int jq=j-D3Q19[q-1][1]; - int kq=k-D3Q19[q-1][2]; - if (iq==0) iq=1; - if (jq==0) jq=1; - if (kq==0) kq=1; - if (iq==Nx-1) iq=Nx-2; - if (jq==Ny-1) jq=Ny-2; - if (kq==Nz-1) kq=Nz-2; - double check = kq*100.f+jq*10.f+iq*1.f+q*0.01; - if (value != check) - printf("Neighbor q=%i, i=%i,j=%i,k=%i: %f \n",q,iq,jq,kq,value); +/* for (k=1; k Date: Mon, 21 Mar 2022 19:51:09 -0400 Subject: [PATCH 33/37] start on cuda function --- cuda/Ion.cu | 83 +++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 58 insertions(+), 25 deletions(-) diff --git a/cuda/Ion.cu b/cuda/Ion.cu index 8db0ae05..f95f3ee4 100644 --- a/cuda/Ion.cu +++ b/cuda/Ion.cu @@ -6,31 +6,6 @@ #define NTHREADS 256 -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) { -} - -extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, - double *dist, double *Den, int memLinks, int Np){ -} - __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){ @@ -620,3 +595,61 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity } //cudaProfilerStop(); } + +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){ + + dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef<<>>(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)); + } +} + From 8ae69e4c1ee147c6f9ace97b234f67671645b34f Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 21 Mar 2022 19:53:05 -0400 Subject: [PATCH 34/37] start on hip function --- hip/Ion.cu | 230 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 230 insertions(+) diff --git a/hip/Ion.cu b/hip/Ion.cu index b1d9636e..ce1fb26d 100644 --- a/hip/Ion.cu +++ b/hip/Ion.cu @@ -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)); + } +} From 2c3272e4231ff55a853905da0f3aa2c4b55897ea Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 22 Mar 2022 17:31:12 -0400 Subject: [PATCH 35/37] updates and fix for user input reader --- common/Domain.cpp | 2 +- common/Membrane.cpp | 181 +++++++++++++++++++++++++++++++++++++---- common/ScaLBL.cpp | 2 - cpu/Ion.cpp | 4 +- tests/TestMembrane.cpp | 9 +- 5 files changed, 177 insertions(+), 21 deletions(-) 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 index 39171944..999c9389 100644 --- a/common/Membrane.cpp +++ b/common/Membrane.cpp @@ -407,7 +407,7 @@ Membrane::~Membrane() { int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArray &Map){ int mlink = 0; - int i,j,k; + int i,j,k,n; int idx, neighbor; double dist, locdist; @@ -824,17 +824,11 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra //...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); //................................................................................... - if (rank == 0) printf(" x count = %i \n",linkCount_x[0]); - if (rank == 0) printf(" X count = %i \n",linkCount_X[0]); - if (rank == 0) printf(" y count = %i \n",linkCount_y[0]); - if (rank == 0) printf(" Y count = %i \n",linkCount_Y[0]); - if (rank == 0) printf(" z count = %i \n",linkCount_z[0]); - if (rank == 0) printf(" Z count = %i \n",linkCount_Z[0]); //...................................................................................... 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+ @@ -846,6 +840,164 @@ int Membrane::Create(std::shared_ptr Dm, DoubleArray &Distance, IntArra 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 Dm, DoubleArray &Distance, IntArra ScaLBL_AllocateZeroCopy((void **) &coefficient_Z, (recvCount_Z - linkCount_Z[0])*sizeof(double)); //...................................................................................... + delete [] TempBuffer; return mlink; } @@ -864,13 +1017,14 @@ int Membrane::D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list, int star int linkCount = 0; int memLinkCount=0; - int i,j,k,n,nn,idx; + int i,j,k,n,nn,idx,link; int * ReturnDist; double distanceNonLocal,distanceLocal; ReturnDist=new int [count]; int *LinkList=new int [count]; int *memLinkList=new int [count]; + //printf("MAP %i RECV VALUES: Nx=%i, Ny=%i, Nz=%i \n",count,Nx,Ny,Nz); for (idx=0; idxLastExterior(), Np); DoubleArray Result(Nx,Ny,Nz); -\ + ScaLBL_Comm->RegularLayout(Map, Ci, Result); /* for (k=1; k Date: Tue, 22 Mar 2022 20:03:34 -0400 Subject: [PATCH 36/37] update cell example --- example/Bubble/cell.db | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/Bubble/cell.db b/example/Bubble/cell.db index dbefcb2d..ccbf37e6 100644 --- a/example/Bubble/cell.db +++ b/example/Bubble/cell.db @@ -1,5 +1,5 @@ MultiphysController { - timestepMax = 40 + timestepMax = 60 num_iter_Ion_List = 2 analysis_interval = 50 tolerance = 1.0e-9 From 8a0937e11170fc8b9d0e7e5f5af505950e8ba6f8 Mon Sep 17 00:00:00 2001 From: James McClure Date: Thu, 24 Mar 2022 07:38:06 -0400 Subject: [PATCH 37/37] add sigmoid to ion equilibrium dist --- cpu/Ion.cpp | 56 +++++++++++++++++++++++------- cuda/Ion.cu | 84 ++++++++++++++++++++++++++++++--------------- hip/Ion.cu | 84 ++++++++++++++++++++++++++++++--------------- models/IonModel.cpp | 4 +-- 4 files changed, 158 insertions(+), 70 deletions(-) diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index 2463bacc..915edda3 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -1,4 +1,5 @@ #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, @@ -250,6 +251,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double Ex, Ey, Ez; //electrical field double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z; double f0, f1, f2, f3, f4, f5, f6; + double X,Y,Z,factor_x, factor_y, factor_z; int nr1, nr2, nr3, nr4, nr5, nr6; for (n = start; n < finish; n++) { @@ -300,33 +302,48 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, FluxElectrical[n + 0 * Np] = uEPx * Ci; FluxElectrical[n + 1 * Np] = uEPy * Ci; FluxElectrical[n + 2 * Np] = uEPz * Ci; + + /* use logistic function to prevent negative distributions*/ + X = 4.0 * (ux + uEPx); + Y = 4.0 * (uy + uEPy); + Z = 4.0 * (uz + uEPz); + factor_x = X / sqrt(1 + X*X); + factor_y = Y / sqrt(1 + Y*Y); + factor_z = Z / sqrt(1 + Z*Z); // q=0 dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci; // q = 1 dist[nr2] = - f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx)); + f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x); + //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx)); + // q=2 dist[nr1] = - f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx)); + f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x); + //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx)); // q = 3 dist[nr4] = - f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy)); + f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y ); + //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy)); // q = 4 dist[nr3] = - f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy)); + f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y); + //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy)); // q = 5 dist[nr6] = - f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz)); + f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z); + //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz)); // q = 6 dist[nr5] = - f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz)); + f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z); + } } @@ -341,6 +358,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion( double Ex, Ey, Ez; //electrical field double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z; double f0, f1, f2, f3, f4, f5, f6; + double X,Y,Z, factor_x, factor_y, factor_z; for (n = start; n < finish; n++) { @@ -377,33 +395,47 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion( FluxElectrical[n + 0 * Np] = uEPx * Ci; FluxElectrical[n + 1 * Np] = uEPy * Ci; FluxElectrical[n + 2 * Np] = uEPz * Ci; + + /* use logistic function to prevent negative distributions*/ + X = 4.0 * (ux + uEPx); + Y = 4.0 * (uy + uEPy); + Z = 4.0 * (uz + uEPz); + factor_x = X / sqrt(1 + X*X); + factor_y = Y / sqrt(1 + Y*Y); + factor_z = Z / sqrt(1 + Z*Z); // q=0 dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci; // q = 1 dist[1 * Np + n] = - f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx)); + f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x); + //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx)); // q=2 dist[2 * Np + n] = - f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx)); + f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x); + //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx)); // q = 3 dist[3 * Np + n] = - f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy)); + f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y); + //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy)); // q = 4 dist[4 * Np + n] = - f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy)); + f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y); + //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy)); // q = 5 dist[5 * Np + n] = - f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz)); + f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z); + //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz)); // q = 6 dist[6 * Np + n] = - f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz)); + f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z); + //f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz)); } } diff --git a/cuda/Ion.cu b/cuda/Ion.cu index f95f3ee4..29362cd1 100644 --- a/cuda/Ion.cu +++ b/cuda/Ion.cu @@ -282,6 +282,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub double Ex,Ey,Ez;//electrical field double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z; double f0,f1,f2,f3,f4,f5,f6; + double X,Y,Z,factor_x,factor_y,factor_z; int nr1,nr2,nr3,nr4,nr5,nr6; int S = Np/NBLOCKS/NTHREADS + 1; @@ -336,34 +337,47 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub FluxElectrical[n+0*Np] = uEPx*Ci; FluxElectrical[n+1*Np] = uEPy*Ci; FluxElectrical[n+2*Np] = uEPz*Ci; + + /* use logistic function to prevent negative distributions*/ + X = 4.0 * (ux + uEPx); + Y = 4.0 * (uy + uEPy); + Z = 4.0 * (uz + uEPz); + factor_x = X / sqrt(1 + X*X); + factor_y = Y / sqrt(1 + Y*Y); + factor_z = Z / sqrt(1 + Z*Z); // q=0 - dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci; - //dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); + dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci; // q = 1 - dist[nr2] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx)); - //dist[nr2] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); + dist[nr2] = + f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x); + //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx)); + // q=2 - dist[nr1] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx)); - //dist[nr1] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); + dist[nr1] = + f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x); + //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx)); // q = 3 - dist[nr4] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy)); - //dist[nr4] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); + dist[nr4] = + f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y ); + //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy)); // q = 4 - dist[nr3] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy)); - //dist[nr3] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); + dist[nr3] = + f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y); + //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy)); // q = 5 - dist[nr6] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz)); - //dist[nr6] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); + dist[nr6] = + f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z); + //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz)); // q = 6 - dist[nr5] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz)); - //dist[nr5] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); + dist[nr5] = + f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z); } } } @@ -377,6 +391,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F double Ex,Ey,Ez;//electrical field double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z; double f0,f1,f2,f3,f4,f5,f6; + double X,Y,Z,factor_x,factor_y,factor_z; int S = Np/NBLOCKS/NTHREADS + 1; for (int s=0; sBarrier(); comm.barrier(); - // *************EVEN TIMESTEP*************// + // *************EVEN TIMESTEP*************// timestep++; //Update ion concentration and charge density IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL @@ -1405,7 +1405,7 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl break; case 21: ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z( - IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],q &Velocity[2 * Np], timestep); break; case 22: