From fce45fece3052fee205f0fe183398f32f1bee0b4 Mon Sep 17 00:00:00 2001
From: Rex Zhe Li <zhe.rex.li@gmail.com>
Date: Thu, 20 May 2021 20:31:40 -0400
Subject: [PATCH 1/2] add slipping velocity BC for CPU; to be built

---
 cpu/D3Q7BC.cpp | 51 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 51 insertions(+)

diff --git a/cpu/D3Q7BC.cpp b/cpu/D3Q7BC.cpp
index 2917bab4..5bc7b589 100644
--- a/cpu/D3Q7BC.cpp
+++ b/cpu/D3Q7BC.cpp
@@ -30,6 +30,57 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int
 	}
 }
 
+extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
+                                                      double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
+                                                      int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
+                                                      double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
+                                                      int count, int Np){
+
+    int idx;
+    int iq,ib,ifluidBC;
+    double value_b,value_q;
+    double Ex,Ey,Ez;
+    double Etx,Ety,Etz;//tangential part of electric field
+    double E_mag_normal;
+    double nsx,nsy,nsz;//unit normal solid gradient
+    double ubx,uby,ubz;//slipping velocity at fluid boundary nodes
+    float cx,cy,cz;//lattice velocity (D3Q19)
+    double LB_weight;//lattice weighting coefficient (D3Q19)
+    double cs2_inv = 3.0;//inverse of cs^2 for D3Q19
+    double nu_LB = (tau-0.5)/cs2_inv;
+
+	for (idx=0; idx<Np; idx++){
+		iq       = BounceBackDist_list[idx];
+        ib       = BounceBackSolid_list[idx];
+        ifluidBC = FluidBoundary_list[idx];
+		value_b = zeta_potential[ib];//get zeta potential from a solid site
+        value_q = dist[iq];
+
+        //Load electric field and compute its tangential componet
+        Ex = ElectricField[ifluidBC+0*Np]; 
+        Ey = ElectricField[ifluidBC+1*Np];
+        Ez = ElectricField[ifluidBC+2*Np];
+        nsx = SolidGrad[ifluidBC+0*Np]; 
+        nsy = SolidGrad[ifluidBC+1*Np];
+        nsz = SolidGrad[ifluidBC+2*Np];
+        E_mag_normal = Ex*nsx+Ey*nsy+Ez*nsz;//magnitude of electric field in the direction normal to solid nodes
+        //compute tangential electric field
+        Etx = Ex - E_mag_normal*nsx;
+        Ety = Ey - E_mag_normal*nsy;
+        Etz = Ez - E_mag_normal*nsz;
+        ubx = -epsilon_LB*value_b*Etx/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;                                                                                                        
+        uby = -epsilon_LB*value_b*Ety/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;                                                                                                        
+        ubz = -epsilon_LB*value_b*Etz/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;                                                                                                        
+
+        //compute bounce-back distribution
+        LB_weight = lattice_weight[idx];
+        cx = lattice_cx[idx];
+        cy = lattice_cy[idx];
+        cz = lattice_cz[idx];
+		dist[iq] = value_q - 2.0*LB_weight*rho0*cs2_inv*(cx*ubx+cy*uby+cz*ubz);
+	}
+}
+
 extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){
 	for (int idx=0; idx<count; idx++){
 		int n = list[idx];

From 367e98be9bcf25d7bf095f2054ee39c73d1742df Mon Sep 17 00:00:00 2001
From: Rex Zhe Li <zhe.rex.li@gmail.com>
Date: Fri, 21 May 2021 00:38:06 -0400
Subject: [PATCH 2/2] fix a minor bug; CPU slipping velocity BC is now working

---
 cpu/D3Q7BC.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cpu/D3Q7BC.cpp b/cpu/D3Q7BC.cpp
index 5bc7b589..aad999f2 100644
--- a/cpu/D3Q7BC.cpp
+++ b/cpu/D3Q7BC.cpp
@@ -49,7 +49,7 @@ extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta
     double cs2_inv = 3.0;//inverse of cs^2 for D3Q19
     double nu_LB = (tau-0.5)/cs2_inv;
 
-	for (idx=0; idx<Np; idx++){
+	for (idx=0; idx<count; idx++){
 		iq       = BounceBackDist_list[idx];
         ib       = BounceBackSolid_list[idx];
         ifluidBC = FluidBoundary_list[idx];