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