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; + } }