update to membrane capability

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

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;

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);

cpu/Ion.cpp

@@ -1,7 +1,91 @@
 #include <stdio.h>
 
-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<memLinks; link++){
+		
+		// inside             	//outside
+		aq = MassFractionIn;	ap = MassFractionOut;  
+		iq = membrane[2*link]; 	ip = membrane[2*link+1];
+		nq = iq%Np;				np = ip%Np;
+		nqm = Map[nq];			npm = Map[np]; // strided layout
+		
+		/* membrane potential for this link */
+		membranePotential = Psi[nqm] - Psi[npm];
+		if (membranePotential > 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;
+    }
 }