From aa26fcafdaac09beb2e572efd9c1d98bc20f6131 Mon Sep 17 00:00:00 2001
From: Rex Zhe Li <zhe.rex.li@gmail.com>
Date: Fri, 28 Aug 2020 11:15:55 -0400
Subject: [PATCH] fix miscellaneous bugs and update the data structure of
 electric potential

---
 common/ScaLBL.cpp                           |  33 ++
 common/ScaLBL.h                             |  40 +-
 cpu/D3Q7BC.cpp                              | 107 ++++++
 cpu/Poisson.cpp                             | 350 +++++++++++++++---
 cpu/Stokes.cpp                              | 162 ++++-----
 models/IonModel.cpp                         |  31 +-
 models/IonModel.h                           |   3 +-
 models/PoissonSolver.cpp                    | 383 ++++++++++++++++++--
 models/PoissonSolver.h                      |  11 +-
 models/StokesModel.cpp                      |  90 +++--
 models/StokesModel.h                        |   7 +-
 tests/lbpm_electrokinetic_dfh_simulator.cpp |  16 +-
 12 files changed, 1027 insertions(+), 206 deletions(-)

diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp
index 854c8b49..ed7c5e59 100644
--- a/common/ScaLBL.cpp
+++ b/common/ScaLBL.cpp
@@ -2053,3 +2053,36 @@ void ScaLBL_Communicator::PrintD3Q19(){
 	delete [] TempBuffer;
 }
 
+void ScaLBL_Communicator::D3Q7_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time){
+	if (kproc == 0) {
+		if (time%2==0){
+			ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(dvcSendList_z, fq, Vin, sendCount_z, N);
+		}
+		else{
+			ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_z(neighborList, dvcSendList_z, fq, Vin, sendCount_z, N);
+		}
+	}
+}
+
+void ScaLBL_Communicator::D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time){
+	if (kproc == nprocz-1){
+		if (time%2==0){
+			ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_Z(dvcSendList_Z, fq, Vout, sendCount_Z, N);
+		}
+		else{
+			ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(neighborList, dvcSendList_Z, fq, Vout, sendCount_Z, N);
+		}
+	}
+}
+
+void ScaLBL_Communicator::Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin){
+	if (kproc == 0) {
+		ScaLBL_Poisson_D3Q7_BC_z(dvcSendList_z, Map, Psi, Vin, sendCount_z);
+	}
+}
+
+void ScaLBL_Communicator::Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout){
+	if (kproc == nprocz-1){
+		ScaLBL_Poisson_D3Q7_BC_Z(dvcSendList_Z, Map, Psi, Vout, sendCount_Z);
+	}
+}
diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index 4d8a3dc3..414653a8 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -46,11 +46,8 @@ extern "C" void ScaLBL_UnpackDenD3Q7(int *list, int count, double *recvbuf, int
 
 extern "C" void ScaLBL_D3Q19_Init(double *Dist, int Np);
 
-
 extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np);
 
-extern "C" void ScaLBL_D3Q19_Momentum_Phys(double *dist, double *vel, double h, double time_conv, int Np);
-
 extern "C" void ScaLBL_D3Q19_Pressure(double *dist, double *press, int Np);
 
 // BGK MODEL
@@ -95,21 +92,30 @@ extern "C" void ScaLBL_IonConcentration_Phys(double *Den, double h, int ion_comp
 
 // LBM Poisson solver
 
-extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,
         int start, int finish, int Np);
 
-extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,
         int start, int finish, int Np);
 
-extern "C" void ScaLBL_D3Q7_Poisson_Init(double *dist, int Np);
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(int *neighborList,int *Map, double *dist, double *Psi, int start, int finish, int Np);
+
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *dist, double *Psi, int start, int finish, int Np);
+
+extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np);
+
+extern "C" void ScaLBL_D3Q7_Poisson_getElectricField(double *dist, double *ElectricField, double tau, int Np);
+
+extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
+        int strideY, int strideZ,int start, int finish, int Np);
 
 // LBM Stokes Model (adapted from MRT model)
 
 extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, 
-                                              double Gx, double Gy, double Gz, double Ex, double Ey, double Ez, int start, int finish, int Np);
+                double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np);
 
-extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,  
-                                             double Gx, double Gy, double Gz, double Ex, double Ey, double Ez, int start, int finish, int Np);
+extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, 
+                double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np);
 
 // MRT MODEL
 extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx,
@@ -190,6 +196,18 @@ extern "C" void ScaLBL_Solid_Dirichlet_D3Q7(double *dist,double *BoundaryValue,i
 
 extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int *BounceBackDist_list,int *BounceBackSolid_list,int N);
 
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np);
+
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np);
+
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np);
+
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np);
+
+extern "C" void ScaLBL_Poisson_D3Q7_BC_z(int *list, int *Map, double *Psi, double Vin, int count);
+
+extern "C" void ScaLBL_Poisson_D3Q7_BC_Z(int *list, int *Map, double *Psi, double Vout, int count);
+
 class ScaLBL_Communicator{
 public:
 	//......................................................................................
@@ -249,6 +267,10 @@ public:
 	void D3Q19_Reflection_BC_z(double *fq);
 	void D3Q19_Reflection_BC_Z(double *fq);
 	double D3Q19_Flux_BC_z(int *neighborList, double *fq, double flux, int time);
+	void D3Q7_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time);
+	void D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time);
+	void Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin);
+	void Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout);
 
 	// Debugging and unit testing functions
 	void PrintD3Q19();
diff --git a/cpu/D3Q7BC.cpp b/cpu/D3Q7BC.cpp
index e7bfd3a4..8c2588d8 100644
--- a/cpu/D3Q7BC.cpp
+++ b/cpu/D3Q7BC.cpp
@@ -30,3 +30,110 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int
 	}
 }
 
+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];
+		double f0 = dist[n];
+		double f1 = dist[2*Np+n];
+		double f2 = dist[1*Np+n];
+		double f3 = dist[4*Np+n];
+		double f4 = dist[3*Np+n];
+		double f6 = dist[5*Np+n];
+		//...................................................
+		double f5 = Vin - (f0+f1+f2+f3+f4+f6);
+		dist[6*Np+n] = f5;
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np){
+	for (int idx=0; idx<count; idx++){
+		int n = list[idx];
+		double f0 = dist[n];
+		double f1 = dist[2*Np+n];
+		double f2 = dist[1*Np+n];
+		double f3 = dist[4*Np+n];
+		double f4 = dist[3*Np+n];
+		double f5 = dist[6*Np+n];
+		//...................................................
+		double f6 = Vout - (f0+f1+f2+f3+f4+f5);
+		dist[5*Np+n] = f6;
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np){
+    int nread,nr5;
+	for (int idx=0; idx<count; idx++){
+		int n = list[idx];
+		double f0 = dist[n];
+
+		nread = d_neighborList[n];
+		double f1 = dist[nread];
+
+		nread = d_neighborList[n+2*Np];
+		double f3 = dist[nread];
+
+		nread = d_neighborList[n+Np];
+		double f2 = dist[nread];
+
+		nread = d_neighborList[n+3*Np];
+		double f4 = dist[nread];
+
+		nread = d_neighborList[n+5*Np];
+		double f6 = dist[nread];
+
+		// Unknown distributions
+		nr5 = d_neighborList[n+4*Np];
+		double f5 = Vin - (f0+f1+f2+f3+f4+f6);
+		dist[nr5] = f5;
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np){
+    int nread,nr6;
+	for (int idx=0; idx<count; idx++){
+		int n = list[idx];
+		double f0 = dist[n];
+
+		nread = d_neighborList[n];
+		double f1 = dist[nread];
+
+		nread = d_neighborList[n+2*Np];
+		double f3 = dist[nread];
+
+		nread = d_neighborList[n+4*Np];
+		double f5 = dist[nread];
+
+		nread = d_neighborList[n+Np];
+		double f2 = dist[nread];
+
+		nread = d_neighborList[n+3*Np];
+		double f4 = dist[nread];
+
+		// unknown distributions
+		nr6 = d_neighborList[n+5*Np];
+		double f6 = Vout - (f0+f1+f2+f3+f4+f5);
+		dist[nr6] = f6;
+	}
+}
+
+extern "C" void ScaLBL_Poisson_D3Q7_BC_z(int *list, int *Map, double *Psi, double Vin, int count)
+{
+	int idx,n,nm;
+
+	for (idx=0; idx<count; idx++){
+		n = list[idx];
+		nm = Map[n];
+		Psi[nm] = Vin;
+	}
+}
+
+extern "C" void ScaLBL_Poisson_D3Q7_BC_Z(int *list, int *Map, double *Psi, double Vout, int count)
+{
+	int idx,n,nm;
+
+	for (idx=0; idx<count; idx++){
+		n = list[idx];
+		nm = Map[n];
+		Psi[nm] = Vout;
+	}
+}
diff --git a/cpu/Poisson.cpp b/cpu/Poisson.cpp
index d33ab658..8deaac8e 100644
--- a/cpu/Poisson.cpp
+++ b/cpu/Poisson.cpp
@@ -1,19 +1,110 @@
 
-extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,
-        int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(int *neighborList,int *Map, double *dist, double *Psi, int start, int finish, int Np){
 	int n;
 	double psi;//electric potential
-    double Ex,Ey,Ez;//electrical field
+	double fq;
+	int nread;
+    int idx;
+
+	for (n=start; n<finish; n++){
+
+		// q=0
+		fq = dist[n];
+        psi = fq;
+
+		// q=1
+		nread = neighborList[n]; 
+		fq = dist[nread];
+		psi += fq;
+		
+		// q=2
+		nread = neighborList[n+Np]; 
+		fq = dist[nread];  
+		psi += fq;
+
+		// q=3
+		nread = neighborList[n+2*Np]; 
+		fq = dist[nread];
+		psi += fq;
+
+		// q = 4
+		nread = neighborList[n+3*Np]; 
+		fq = dist[nread];
+		psi += fq;
+
+		// q=5
+		nread = neighborList[n+4*Np];
+		fq = dist[nread];
+		psi += fq;
+
+		// q = 6
+		nread = neighborList[n+5*Np];
+		fq = dist[nread];
+		psi += fq;
+        
+        idx=Map[n];
+        Psi[idx] = psi;
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *dist, double *Psi, int start, int finish, int Np){
+	int n;
+	double psi;//electric potential
+	double fq;
+    int idx;
+
+	for (n=start; n<finish; n++){
+
+		// q=0
+		fq = dist[n];
+		psi = fq;
+		
+		// q=1
+		fq = dist[2*Np+n];
+		psi += fq;
+
+		// q=2 
+		fq = dist[1*Np+n];
+		psi += fq;
+
+		// q=3
+		fq = dist[4*Np+n];
+		psi += fq;
+
+		// q=4
+		fq = dist[3*Np+n];
+		psi += fq;
+
+		// q=5
+		fq = dist[6*Np+n];
+		psi += fq;
+
+		// q=6
+		fq = dist[5*Np+n];
+		psi += fq;
+
+        idx=Map[n];
+        Psi[idx] = psi;
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,int start, int finish, int Np){
+	int n;
+	double psi;//electric potential
+    double Ex,Ey,Ez;//electric field
     double rho_e;//local charge density
 	double f0,f1,f2,f3,f4,f5,f6;
 	int nr1,nr2,nr3,nr4,nr5,nr6;
     double rlx=1.0/tau;
+    int idx;
 
 	for (n=start; n<finish; n++){
 
         //Load data
         rho_e = Den_charge[n];
         rho_e = gamma*rho_e/epsilon_LB;
+        idx=Map[n];
+        psi = Psi[idx];
 		
 		// q=0
 		f0 = dist[n];
@@ -40,53 +131,63 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, doubl
 		nr6 = neighborList[n+5*Np];
 		f6 = dist[nr6];
 		
-		psi = f0+f2+f1+f4+f3+f6+f5;
-		Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
-		Ey = (f3-f4)*rlx*4.5;
-		Ez = (f5-f6)*rlx*4.5;
-        ElectricField[n+0*Np] = Ex;
-        ElectricField[n+1*Np] = Ey;
-        ElectricField[n+2*Np] = Ez;
-        Psi[n] = psi;
+		//Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
+		//Ey = (f3-f4)*rlx*4.5;
+		//Ez = (f5-f6)*rlx*4.5;
+		//Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
+		//Ey = (f3-f4)*rlx*4.0;
+		//Ez = (f5-f6)*rlx*4.0;
+        //ElectricField[n+0*Np] = Ex;
+        //ElectricField[n+1*Np] = Ey;
+        //ElectricField[n+2*Np] = Ez;
 
 		// q = 0
-		dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e);
+		//dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e);
+		dist[n] = f0*(1.0-rlx) + 0.25*(rlx*psi+rho_e);
 
 		// q = 1
-		dist[nr2] = f1*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		//dist[nr2] = f1*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		dist[nr2] = f1*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
 
 		// q = 2
-		dist[nr1] = f2*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		//dist[nr1] = f2*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		dist[nr1] = f2*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
 
 		// q = 3
-		dist[nr4] = f3*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		//dist[nr4] = f3*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		dist[nr4] = f3*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
 
 		// q = 4
-		dist[nr3] = f4*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		//dist[nr3] = f4*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		dist[nr3] = f4*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
 
 		// q = 5
-		dist[nr6] = f5*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		//dist[nr6] = f5*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		dist[nr6] = f5*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
 
 		// q = 6
-		dist[nr5] = f6*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		//dist[nr5] = f6*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		dist[nr5] = f6*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
 		//........................................................................
 	}
 }
 
-extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,
-        int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,int start, int finish, int Np){
 	int n;
 	double psi;//electric potential
-    double Ex,Ey,Ez;//electrical field
+    double Ex,Ey,Ez;//electric field
     double rho_e;//local charge density
 	double f0,f1,f2,f3,f4,f5,f6;
     double rlx=1.0/tau;
+    int idx;
 
 	for (n=start; n<finish; n++){
 
         //Load data
         rho_e = Den_charge[n];
         rho_e = gamma*rho_e/epsilon_LB;
+        idx=Map[n];
+        psi = Psi[idx];
 
 		f0 = dist[n];
 		f1 = dist[2*Np+n];
@@ -97,49 +198,208 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, dou
 		f6 = dist[5*Np+n];
 
 
-		psi = f0+f2+f1+f4+f3+f6+f5;
-		Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
-		Ey = (f3-f4)*rlx*4.5;
-		Ez = (f5-f6)*rlx*4.5;
-        ElectricField[n+0*Np] = Ex;
-        ElectricField[n+1*Np] = Ey;
-        ElectricField[n+2*Np] = Ez;
-        Psi[n] = psi;
+		//Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
+		//Ey = (f3-f4)*rlx*4.5;
+		//Ez = (f5-f6)*rlx*4.5;
+		//Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
+		//Ey = (f3-f4)*rlx*4.0;
+		//Ez = (f5-f6)*rlx*4.0;
+        //ElectricField[n+0*Np] = Ex;
+        //ElectricField[n+1*Np] = Ey;
+        //ElectricField[n+2*Np] = Ez;
 
 		// q = 0
-		dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e);
+		//dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e);
+		dist[n] = f0*(1.0-rlx) + 0.25*(rlx*psi+rho_e);
 
 		// q = 1
-		dist[1*Np+n] = f1*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		//dist[1*Np+n] = f1*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		dist[1*Np+n] = f1*(1.0-rlx) + 0.125*(rlx*psi+rho_e); 
 
 		// q = 2
-		dist[2*Np+n] = f2*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		//dist[2*Np+n] = f2*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		dist[2*Np+n] = f2*(1.0-rlx) + 0.125*(rlx*psi+rho_e); 
 
 		// q = 3
-		dist[3*Np+n] = f3*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		//dist[3*Np+n] = f3*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
+		dist[3*Np+n] = f3*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
 
 		// q = 4
-		dist[4*Np+n] = f4*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		//dist[4*Np+n] = f4*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		dist[4*Np+n] = f4*(1.0-rlx) + 0.125*(rlx*psi+rho_e); 
 
 		// q = 5
-		dist[5*Np+n] = f5*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		//dist[5*Np+n] = f5*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		dist[5*Np+n] = f5*(1.0-rlx) + 0.125*(rlx*psi+rho_e); 
 
 		// q = 6
-		dist[6*Np+n] = f6*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		//dist[6*Np+n] = f6*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); 
+		dist[6*Np+n] = f6*(1.0-rlx) + 0.125*(rlx*psi+rho_e); 
 		//........................................................................
 	}
 }
 
-extern "C" void ScaLBL_D3Q7_Poisson_Init(double *dist, int Np)
+extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np)
 {
 	int n;
-	for (n=0; n<Np; n++){
-		dist[0*Np+n] = 0.3333333333333333;
-		dist[1*Np+n] = 0.1111111111111111;		
-		dist[2*Np+n] = 0.1111111111111111;	
-		dist[3*Np+n] = 0.1111111111111111;	
-		dist[4*Np+n] = 0.1111111111111111;	
-		dist[5*Np+n] = 0.1111111111111111;	
-		dist[6*Np+n] = 0.1111111111111111;	
+    int ijk;
+	for (n=start; n<finish; n++){
+        ijk = Map[n];
+		//dist[0*Np+n] = 0.3333333333333333*Psi[n];
+		//dist[1*Np+n] = 0.1111111111111111*Psi[n];		
+		//dist[2*Np+n] = 0.1111111111111111*Psi[n];	
+		//dist[3*Np+n] = 0.1111111111111111*Psi[n];	
+		//dist[4*Np+n] = 0.1111111111111111*Psi[n];	
+		//dist[5*Np+n] = 0.1111111111111111*Psi[n];	
+		//dist[6*Np+n] = 0.1111111111111111*Psi[n];	
+		//dist[0*Np+n] = 0.25*Psi[n];
+		//dist[1*Np+n] = 0.125*Psi[n];		
+		//dist[2*Np+n] = 0.125*Psi[n];	
+		//dist[3*Np+n] = 0.125*Psi[n];	
+		//dist[4*Np+n] = 0.125*Psi[n];	
+		//dist[5*Np+n] = 0.125*Psi[n];	
+		//dist[6*Np+n] = 0.125*Psi[n];	
+
+		dist[0*Np+n] = 0.25*Psi[ijk];
+		dist[1*Np+n] = 0.125*Psi[ijk];		
+		dist[2*Np+n] = 0.125*Psi[ijk];	
+		dist[3*Np+n] = 0.125*Psi[ijk];	
+		dist[4*Np+n] = 0.125*Psi[ijk];	
+		dist[5*Np+n] = 0.125*Psi[ijk];	
+		dist[6*Np+n] = 0.125*Psi[ijk];	
 	}
 }
+
+extern "C" void ScaLBL_D3Q7_Poisson_getElectricField(double *dist, double *ElectricField, double tau, int Np){
+	int n;
+	// distributions
+	double f1,f2,f3,f4,f5,f6;
+	double Ex,Ey,Ez;
+    double rlx=1.0/tau;
+
+	for (n=0; n<Np; n++){
+		//........................................................................
+		// Registers to store the distributions
+		//........................................................................
+		f1 = dist[Np+n];
+		f2 = dist[2*Np+n];
+		f3 = dist[3*Np+n];
+		f4 = dist[4*Np+n];
+		f5 = dist[5*Np+n];
+		f6 = dist[6*Np+n];
+		//.................Compute the Electric Field...................................
+		//Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
+		//Ey = (f3-f4)*rlx*4.5;
+		//Ez = (f5-f6)*rlx*4.5;
+		Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
+		Ey = (f3-f4)*rlx*4.0;
+		Ez = (f5-f6)*rlx*4.0;
+		//..................Write the Electric Field.....................................
+		ElectricField[0*Np+n] = Ex;
+		ElectricField[1*Np+n] = Ey;
+		ElectricField[2*Np+n] = Ez;
+		//........................................................................
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
+        int strideY, int strideZ,int start, int finish, int Np){
+
+	int n,nn;
+    int ijk;
+    int id;
+	// distributions
+	double m1,m2,m3,m4,m5,m6,m7,m8,m9;
+	double m10,m11,m12,m13,m14,m15,m16,m17,m18;
+	double nx,ny,nz;
+
+	for (n=start; n<finish; n++){
+
+		// Get the 1D index based on regular data layout
+		ijk = Map[n];
+		//					COMPUTE THE COLOR GRADIENT
+		//........................................................................
+		//.................Read Phase Indicator Values............................
+		//........................................................................
+		nn = ijk-1;							// neighbor index (get convention)
+        id = ID[nn];
+		m1 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 1
+		//........................................................................
+		nn = ijk+1;							// neighbor index (get convention)
+        id = ID[nn];
+		m2 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 2
+		//........................................................................
+		nn = ijk-strideY;							// neighbor index (get convention)
+        id = ID[nn];
+		m3 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 3
+		//........................................................................
+		nn = ijk+strideY;							// neighbor index (get convention)
+        id = ID[nn];
+		m4 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 4
+		//........................................................................
+		nn = ijk-strideZ;						// neighbor index (get convention)
+        id = ID[nn];
+		m5 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 5
+		//........................................................................
+		nn = ijk+strideZ;						// neighbor index (get convention)
+        id = ID[nn];
+		m6 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 6
+		//........................................................................
+		nn = ijk-strideY-1;						// neighbor index (get convention)
+        id = ID[nn];
+		m7 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 7
+		//........................................................................
+		nn = ijk+strideY+1;						// neighbor index (get convention)
+        id = ID[nn];
+		m8 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 8
+		//........................................................................
+		nn = ijk+strideY-1;						// neighbor index (get convention)
+        id = ID[nn];
+		m9 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 9
+		//........................................................................
+		nn = ijk-strideY+1;						// neighbor index (get convention)
+        id = ID[nn];
+		m10 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 10
+		//........................................................................
+		nn = ijk-strideZ-1;						// neighbor index (get convention)
+        id = ID[nn];
+		m11 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 11
+		//........................................................................
+		nn = ijk+strideZ+1;						// neighbor index (get convention)
+        id = ID[nn];
+		m12 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 12
+		//........................................................................
+		nn = ijk+strideZ-1;						// neighbor index (get convention)
+        id = ID[nn];
+		m13 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 13
+		//........................................................................
+		nn = ijk-strideZ+1;						// neighbor index (get convention)
+        id = ID[nn];
+		m14 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 14
+		//........................................................................
+		nn = ijk-strideZ-strideY;					// neighbor index (get convention)
+        id = ID[nn];
+		m15 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 15
+		//........................................................................
+		nn = ijk+strideZ+strideY;					// neighbor index (get convention)
+        id = ID[nn];
+		m16 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 16
+		//........................................................................
+		nn = ijk+strideZ-strideY;					// neighbor index (get convention)
+        id = ID[nn];
+		m17 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 17
+		//........................................................................
+		nn = ijk-strideZ+strideY;					// neighbor index (get convention)
+        id = ID[nn];
+		m18 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 18
+		//............Compute the Color Gradient...................................
+		nx = -1.f/18.f*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
+		ny = -1.f/18.f*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
+		nz = -1.f/18.f*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
+		
+		ElectricField[n] = nx;
+		ElectricField[Np+n] = ny;
+		ElectricField[2*Np+n] = nz;
+	}
+}
+
diff --git a/cpu/Stokes.cpp b/cpu/Stokes.cpp
index a31a8bed..a3842345 100644
--- a/cpu/Stokes.cpp
+++ b/cpu/Stokes.cpp
@@ -1,7 +1,6 @@
 #include <stdio.h>
 
-extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, 
-                                              double Gx, double Gy, double Gz, double Ex_const, double Ey_const, double Ez_const, int start, int finish, int Np)
+extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np)
 {
     double fq;
 	// conserved momemnts
@@ -32,13 +31,13 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
         
         //Load data
         rhoE = ChargeDensity[n];
-        Ex = ElectricField[n+0*Np]+Ex_const;
-        Ey = ElectricField[n+1*Np]+Ey_const;
-        Ez = ElectricField[n+2*Np]+Ez_const;
+        Ex = ElectricField[n+0*Np];
+        Ey = ElectricField[n+1*Np];
+        Ez = ElectricField[n+2*Np];
         //compute total body force, including input body force (Gx,Gy,Gz)
-        Fx = Gx + rhoE*Ex;
-        Fy = Gy + rhoE*Ey;
-        Fz = Gz + rhoE*Ez;
+        Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB
+        Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
+        Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
 
 		// q=0
 		fq = dist[n];
@@ -311,9 +310,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
 		m18 -= fq;
 
 		// write the velocity 
-		ux = jx / rho;
-		uy = jy / rho;
-		uz = jz / rho;
+		ux = jx / rho0;
+		uy = jy / rho0;
+		uz = jz / rho0;
 		Velocity[n] = ux;
 		Velocity[Np+n] = uy;
 		Velocity[2*Np+n] = uz;
@@ -326,18 +325,18 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
 
 		//..............incorporate external force................................................
 		//..............carry out relaxation process...............................................
-		m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) - m1);
-		m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho) - m2);
+		m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) - m1);
+		m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0) - m2);
 		m4 = m4 + rlx_setB*((-0.6666666666666666*jx) - m4);
 		m6 = m6 + rlx_setB*((-0.6666666666666666*jy) - m6);
 		m8 = m8 + rlx_setB*((-0.6666666666666666*jz) - m8);
-		m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) - m9);
+		m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) - m9);
 		m10 = m10 + rlx_setA*(-0.5*((2*jx*jx-jy*jy-jz*jz)/rho) - m10);
-		m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) - m11);
-		m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho) - m12);
-		m13 = m13 + rlx_setA*((jx*jy/rho) - m13);
-		m14 = m14 + rlx_setA*((jy*jz/rho) - m14);
-		m15 = m15 + rlx_setA*((jx*jz/rho) - m15);
+		m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) - m11);
+		m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho0) - m12);
+		m13 = m13 + rlx_setA*((jx*jy/rho0) - m13);
+		m14 = m14 + rlx_setA*((jy*jz/rho0) - m14);
+		m15 = m15 + rlx_setA*((jx*jz/rho0) - m15);
 		m16 = m16 + rlx_setB*( - m16);
 		m17 = m17 + rlx_setB*( - m17);
 		m18 = m18 + rlx_setB*( - m18);
@@ -454,8 +453,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
 	}
 }
 
-extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,  
-                                             double Gx, double Gy, double Gz, double Ex_const, double Ey_const, double Ez_const, int start, int finish, int Np)
+extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np)
 {
     double fq;
 	// conserved momemnts
@@ -487,13 +485,13 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
 
         //Load data
         rhoE = ChargeDensity[n];
-        Ex = ElectricField[n+0*Np]+Ex_const;
-        Ey = ElectricField[n+1*Np]+Ey_const;
-        Ez = ElectricField[n+2*Np]+Ez_const;
+        Ex = ElectricField[n+0*Np];
+        Ey = ElectricField[n+1*Np];
+        Ez = ElectricField[n+2*Np];
         //compute total body force, including input body force (Gx,Gy,Gz)
-        Fx = Gx + rhoE*Ex;
-        Fy = Gy + rhoE*Ey;
-        Fz = Gz + rhoE*Ez;
+        Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
+        Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
+        Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
 
 		// q=0
 		fq = dist[n];
@@ -803,27 +801,27 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
 		m18 -= fq;
 
 		// write the velocity 
-		ux = jx / rho;
-		uy = jy / rho;
-		uz = jz / rho;
+		ux = jx / rho0;
+		uy = jy / rho0;
+		uz = jz / rho0;
 		Velocity[n] = ux;
 		Velocity[Np+n] = uy;
 		Velocity[2*Np+n] = uz;
 
 		//..............incorporate external force................................................
 		//..............carry out relaxation process...............................................
-		m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) - m1);
-		m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho) - m2);
+		m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) - m1);
+		m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0) - m2);
 		m4 = m4 + rlx_setB*((-0.6666666666666666*jx) - m4);
 		m6 = m6 + rlx_setB*((-0.6666666666666666*jy) - m6);
 		m8 = m8 + rlx_setB*((-0.6666666666666666*jz) - m8);
-		m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) - m9);
+		m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) - m9);
 		m10 = m10 + rlx_setA*(-0.5*((2*jx*jx-jy*jy-jz*jz)/rho) - m10);
-		m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) - m11);
-		m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho) - m12);
-		m13 = m13 + rlx_setA*((jx*jy/rho) - m13);
-		m14 = m14 + rlx_setA*((jy*jz/rho) - m14);
-		m15 = m15 + rlx_setA*((jx*jz/rho) - m15);
+		m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) - m11);
+		m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho0) - m12);
+		m13 = m13 + rlx_setA*((jx*jy/rho0) - m13);
+		m14 = m14 + rlx_setA*((jy*jz/rho0) - m14);
+		m15 = m15 + rlx_setA*((jx*jz/rho0) - m15);
 		m16 = m16 + rlx_setB*( - m16);
 		m17 = m17 + rlx_setB*( - m17);
 		m18 = m18 + rlx_setB*( - m18);
@@ -955,47 +953,47 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
 	}
 }
 
-extern "C" void ScaLBL_D3Q19_Momentum_Phys(double *dist, double *vel, double h, double time_conv, int Np)
-{
-    //h: resolution [um/lu]
-    //time_conv: time conversion factor [sec/lt]
-	int n;
-	// distributions
-	double f1,f2,f3,f4,f5,f6,f7,f8,f9;
-	double f10,f11,f12,f13,f14,f15,f16,f17,f18;
-	double vx,vy,vz;
-
-	for (n=0; n<Np; n++){
-		//........................................................................
-		// Registers to store the distributions
-		//........................................................................
-		f2 = dist[2*Np+n];
-		f4 = dist[4*Np+n];
-		f6 = dist[6*Np+n];
-		f8 = dist[8*Np+n];
-		f10 = dist[10*Np+n];
-		f12 = dist[12*Np+n];
-		f14 = dist[14*Np+n];
-		f16 = dist[16*Np+n];
-		f18 = dist[18*Np+n];
-		//........................................................................
-		f1 = dist[Np+n];
-		f3 = dist[3*Np+n];
-		f5 = dist[5*Np+n];
-		f7 = dist[7*Np+n];
-		f9 = dist[9*Np+n];
-		f11 = dist[11*Np+n];
-		f13 = dist[13*Np+n];
-		f15 = dist[15*Np+n];
-		f17 = dist[17*Np+n];
-		//.................Compute the velocity...................................
-		vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
-		vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
-		vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
-		//..................Write the velocity.....................................
-		vel[0*Np+n] = vx*(h*1.0e-6)/time_conv;
-		vel[1*Np+n] = vy*(h*1.0e-6)/time_conv;
-		vel[2*Np+n] = vz*(h*1.0e-6)/time_conv;
-		//........................................................................
-	}
-}
+//extern "C" void ScaLBL_D3Q19_Momentum_Phys(double *dist, double *vel, double h, double time_conv, int Np)
+//{
+//    //h: resolution [um/lu]
+//    //time_conv: time conversion factor [sec/lt]
+//	int n;
+//	// distributions
+//	double f1,f2,f3,f4,f5,f6,f7,f8,f9;
+//	double f10,f11,f12,f13,f14,f15,f16,f17,f18;
+//	double vx,vy,vz;
+//
+//	for (n=0; n<Np; n++){
+//		//........................................................................
+//		// Registers to store the distributions
+//		//........................................................................
+//		f2 = dist[2*Np+n];
+//		f4 = dist[4*Np+n];
+//		f6 = dist[6*Np+n];
+//		f8 = dist[8*Np+n];
+//		f10 = dist[10*Np+n];
+//		f12 = dist[12*Np+n];
+//		f14 = dist[14*Np+n];
+//		f16 = dist[16*Np+n];
+//		f18 = dist[18*Np+n];
+//		//........................................................................
+//		f1 = dist[Np+n];
+//		f3 = dist[3*Np+n];
+//		f5 = dist[5*Np+n];
+//		f7 = dist[7*Np+n];
+//		f9 = dist[9*Np+n];
+//		f11 = dist[11*Np+n];
+//		f13 = dist[13*Np+n];
+//		f15 = dist[15*Np+n];
+//		f17 = dist[17*Np+n];
+//		//.................Compute the velocity...................................
+//		vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
+//		vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
+//		vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
+//		//..................Write the velocity.....................................
+//		vel[0*Np+n] = vx*(h*1.0e-6)/time_conv;
+//		vel[1*Np+n] = vy*(h*1.0e-6)/time_conv;
+//		vel[2*Np+n] = vz*(h*1.0e-6)/time_conv;
+//		//........................................................................
+//	}
+//}
diff --git a/models/IonModel.cpp b/models/IonModel.cpp
index 55c848b9..c8c1fba8 100644
--- a/models/IonModel.cpp
+++ b/models/IonModel.cpp
@@ -364,11 +364,16 @@ void ScaLBL_IonModel::Initialize(){
 
 void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
 
+    //Input parameter:
+    //1. Velocity is from StokesModel
+    //2. ElectricField is from Poisson model
+
     //LB-related parameter
     vector<double> rlx(tau.begin(),tau.end());
     for (double item : rlx){
         item = 1.0/item; 
     }
+    
 	//.......create and start timer............
 	//double starttime,stoptime,cputime;
 	//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
@@ -462,7 +467,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
 
 }
 
-void ScaLBL_IonModel::getIonConcentration(){
+//TODO this ruin the ion concentration on device
+//need to do something similar to electric field
+void ScaLBL_IonModel::getIonConcentration(int timestep){
 	for (int ic=0; ic<number_ion_species; ic++){
         ScaLBL_IonConcentration_Phys(Ci, h, ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         ScaLBL_IonConcentration_Phys(Ci, h, ic, 0, ScaLBL_Comm->LastExterior(), Np);
@@ -474,7 +481,7 @@ void ScaLBL_IonModel::getIonConcentration(){
         ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 
         FILE *OUTFILE;
-        sprintf(LocalRankFilename,"Ion%02i.%05i.raw",ic+1,rank);
+        sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank);
         OUTFILE = fopen(LocalRankFilename,"wb");
         fwrite(PhaseField.data(),8,N,OUTFILE);
         fclose(OUTFILE);
@@ -482,3 +489,23 @@ void ScaLBL_IonModel::getIonConcentration(){
 
 }
 
+//void ScaLBL_IonModel::getIonConcentration(){
+//	for (int ic=0; ic<number_ion_species; ic++){
+//        ScaLBL_IonConcentration_Phys(Ci, h, ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+//        ScaLBL_IonConcentration_Phys(Ci, h, ic, 0, ScaLBL_Comm->LastExterior(), Np);
+//    }
+//
+//    DoubleArray PhaseField(Nx,Ny,Nz);
+//	for (int ic=0; ic<number_ion_species; ic++){
+//	    ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
+//        ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+//
+//        FILE *OUTFILE;
+//        sprintf(LocalRankFilename,"Ion%02i.%05i.raw",ic+1,rank);
+//        OUTFILE = fopen(LocalRankFilename,"wb");
+//        fwrite(PhaseField.data(),8,N,OUTFILE);
+//        fclose(OUTFILE);
+//    }
+//
+//}
+
diff --git a/models/IonModel.h b/models/IonModel.h
index 6232d105..59e5b6e6 100644
--- a/models/IonModel.h
+++ b/models/IonModel.h
@@ -29,7 +29,7 @@ public:
 	void Create();
 	void Initialize();
 	void Run(double *Velocity, double *ElectricField);
-    void getIonConcentration();
+    void getIonConcentration(int timestep);
 	
 	//bool Restart,pBC;
 	int timestep,timestepMax;
@@ -40,6 +40,7 @@ public:
     double kb,electron_charge,T,Vt;
     double k2_inv;
     double tolerance;
+    double Ex,Ey,Ez;
 	
 	int number_ion_species;
     vector<double> IonDiffusivity;//User input unit [m^2/sec]
diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp
index 701c009c..c3e5c019 100644
--- a/models/PoissonSolver.cpp
+++ b/models/PoissonSolver.cpp
@@ -7,7 +7,7 @@
 
 ScaLBL_Poisson::ScaLBL_Poisson(int RANK, int NP, MPI_Comm COMM):
 rank(RANK), nprocs(NP),timestep(0),timestepMax(0),tau(0),k2_inv(0),gamma(0),tolerance(0),h(0),
-epsilon0(0),epsilon0_LB(0),epsilonR(0),epsilon_LB(0),Nx(0),Ny(0),Nz(0),N(0),Np(0),analysis_interval(0),
+epsilon0(0),epsilon0_LB(0),epsilonR(0),epsilon_LB(0),Vin(0),Vout(0),Nx(0),Ny(0),Nz(0),N(0),Np(0),analysis_interval(0),
 nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),BoundaryConditionSolid(0),Lx(0),Ly(0),Lz(0),comm(COMM)
 {
 
@@ -22,17 +22,20 @@ void ScaLBL_Poisson::ReadParams(string filename){
 	domain_db = db->getDatabase( "Domain" );
 	electric_db = db->getDatabase( "Poisson" );
 	
-    k2_inv = 4.5;//speed of sound for D3Q7 lattice 
-    gamma = 0.3;//time step of LB-Poisson equation
+    //k2_inv = 4.5;//speed of sound for D3Q7 lattice 
+    k2_inv = 4.0;//speed of sound for D3Q7 lattice 
+    gamma = 1.0;//time step of LB-Poisson equation
 	tau = 0.5+k2_inv*gamma;
 	timestepMax = 100000;
 	tolerance = 1.0e-6;//stopping criterion for obtaining steady-state electricla potential
     h = 1.0;//resolution; unit: um/lu
-    epsilon0 = 8.85e-12;//electrical permittivity of vaccum; unit:[C/(V*m)]
+    epsilon0 = 8.85e-12;//electric permittivity of vaccum; unit:[C/(V*m)]
     epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
     epsilonR = 78.4;//default dielectric constant of water
-    epsilon_LB = epsilon0_LB*epsilonR;//electrical permittivity 
+    epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity 
     analysis_interval = 1000; 
+    Vin  = 1.0; //Boundary-z (inlet)  electric potential
+    Vout = 1.0; //Boundary-Z (outlet) electric potential
 
 	// LB-Poisson Model parameters
 	if (electric_db->keyExists( "timestepMax" )){
@@ -56,20 +59,23 @@ void ScaLBL_Poisson::ReadParams(string filename){
 	if (electric_db->keyExists( "BC_Solid" )){
 		BoundaryConditionSolid = electric_db->getScalar<int>( "BC_Solid" );
 	}
+    // Read boundary condition for electric potentiona
+    // BC = 0: normal periodic BC
+    // BC = 1: fixed inlet and outlet potential
+    BoundaryCondition = 0;
+	if (electric_db->keyExists( "BC" )){
+		BoundaryCondition = electric_db->getScalar<int>( "BC" );
+	}
 
 	// Read domain parameters
 	if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
 		h = domain_db->getScalar<double>( "voxel_length" );
 	}
 
-    BoundaryCondition = 0;
-	if (domain_db->keyExists( "BC" )){
-		BoundaryCondition = domain_db->getScalar<int>( "BC" );
-	}
 
     //Re-calcualte model parameters if user updates input
     epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
-    epsilon_LB = epsilon0_LB*epsilonR;//electrical permittivity 
+    epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity 
 	tau = 0.5+k2_inv*gamma;
 
 	if (rank==0) printf("***********************************************************************************\n");
@@ -202,13 +208,13 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid)
                 AFFINITY=0.f;
 				// Assign the affinity from the paired list
 				for (unsigned int idx=0; idx < NLABELS; idx++){
-				      //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
 					if (VALUE == LabelList[idx]){
 						AFFINITY=AffinityList[idx];
                         //NOTE need to convert the user input phys unit to LB unit
                         if (BoundaryConditionSolid==2){
                             //for BCS=1, i.e. Dirichlet-type, no need for unit conversion
-                            AFFINITY = AFFINITY*(h*h*1.0e-12); 
+                            //TODO maybe there is a factor of gamm missing here ?
+                            AFFINITY = AFFINITY*(h*h*1.0e-12)/epsilon_LB; 
                         }
 						label_count[idx] += 1.0;
 						idx = NLABELS;
@@ -244,7 +250,6 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid)
 	}
 }
 
-
 void ScaLBL_Poisson::Create(){
 	/*
 	 *  This function creates the variables needed to run a LBM 
@@ -260,6 +265,7 @@ void ScaLBL_Poisson::Create(){
 	// Create a communicator for the device (will use optimized layout)
 	// ScaLBL_Communicator ScaLBL_Comm(Mask); // original
 	ScaLBL_Comm  = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
+	ScaLBL_Comm_Regular  = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
 
 	int Npad=(Np/16 + 2)*16;
 	if (rank==0)    printf ("LB-Poisson Solver: Set up memory efficient layout \n");
@@ -277,37 +283,125 @@ void ScaLBL_Poisson::Create(){
 	int neighborSize=18*(Np*sizeof(int));
 	//...........................................................................
 	ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
+	ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &dvcID, sizeof(signed char)*Nx*Ny*Nz);
 	ScaLBL_AllocateDeviceMemory((void **) &fq, 7*dist_mem_size);  
-	ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz);
 	ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np);
-	ScaLBL_AllocateDeviceMemory((void **) &PoissonSolid, sizeof(double)*Nx*Ny*Nz);  
+	//ScaLBL_AllocateDeviceMemory((void **) &PoissonSolid, sizeof(double)*Nx*Ny*Nz);  
 	//...........................................................................
 	
 	// Update GPU data structures
 	if (rank==0)    printf ("LB-Poisson Solver: Setting up device map and neighbor list \n");
+	fflush(stdout);
+	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_DeviceBarrier();
+	delete [] TmpMap;
 	// copy the neighbor list 
 	ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
 	ScaLBL_DeviceBarrier();
 	MPI_Barrier(comm);
+	delete [] neighborList;
+    // copy node ID
+	ScaLBL_CopyToDevice(dvcID, Mask->id, sizeof(signed char)*Nx*Ny*Nz);
+	ScaLBL_DeviceBarrier();
 	
-    //Initialize solid boundary for electrical potential
+    //Initialize solid boundary for electric potential
     ScaLBL_Comm->SetupBounceBackList(Map, Mask->id, Np);
 	MPI_Barrier(comm);
 
-    double *PoissonSolid_host;
-    PoissonSolid_host = new double[Nx*Ny*Nz];
-    AssignSolidBoundary(PoissonSolid_host);
-	ScaLBL_CopyToDevice(PoissonSolid, PoissonSolid_host, Nx*Ny*Nz*sizeof(double));
-	ScaLBL_DeviceBarrier();
-    delete [] PoissonSolid_host;
+    //double *PoissonSolid_host;
+    //PoissonSolid_host = new double[Nx*Ny*Nz];
+    //AssignSolidBoundary(PoissonSolid_host);
+	//ScaLBL_CopyToDevice(PoissonSolid, PoissonSolid_host, Nx*Ny*Nz*sizeof(double));
+	//ScaLBL_DeviceBarrier();
+    //delete [] PoissonSolid_host;
 }        
+// Method 1
+// Psi - size N
+// ID_dvc - size N
+// Method 2
+// Psi - size Np
+// PoissonSolid size N
+
+void ScaLBL_Poisson::Potential_Init(double *psi_init){
+
+    if (BoundaryCondition==1){
+	    if (electric_db->keyExists( "Vin" )){
+	    	Vin = electric_db->getScalar<double>( "Vin" );
+	    }
+	    if (electric_db->keyExists( "Vout" )){
+	    	Vout = electric_db->getScalar<double>( "Vout" );
+	    }
+    }
+    //By default only periodic BC is applied and Vin=Vout=1.0, i.e. there is no potential gradient along Z-axis
+    double slope = (Vout-Vin)/(Nz-2);
+    double psi_linearized;
+	for (int k=0;k<Nz;k++){
+        if (k==0 || k==1){
+            psi_linearized = Vin;
+        }
+        else if (k==Nz-1 || k==Nz-2){
+            psi_linearized = Vout;
+        }
+        else{
+            psi_linearized = slope*(k-1)+Vin;
+        }
+		for (int j=0;j<Ny;j++){
+			for (int i=0;i<Nx;i++){
+				int n = k*Nx*Ny+j*Nx+i;
+				if (Mask->id[n]>0){
+                    psi_init[n] = psi_linearized;
+                }
+            }
+        }
+    }
+}
+
 
 void ScaLBL_Poisson::Initialize(){
 	/*
 	 * This function initializes model
 	 */
     if (rank==0)    printf ("LB-Poisson Solver: initializing D3Q7 distributions\n");
-    ScaLBL_D3Q7_Poisson_Init(fq, Np);
+    //NOTE the initialization involves two steps:
+    //1. assign solid boundary value (surface potential or surface change density)
+    //2. Initialize electric potential for pore nodes
+    double *psi_host;
+    psi_host = new double [Nx*Ny*Nz];
+    AssignSolidBoundary(psi_host);//step1
+    Potential_Init(psi_host);//step2
+	ScaLBL_CopyToDevice(Psi, psi_host, Nx*Ny*Nz*sizeof(double));
+	ScaLBL_DeviceBarrier();
+    ScaLBL_D3Q7_Poisson_Init(dvcMap, fq, Psi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+    ScaLBL_D3Q7_Poisson_Init(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np);
+    delete [] psi_host;
 }
 
 void ScaLBL_Poisson::Run(double *ChargeDensity){
@@ -325,30 +419,83 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
 		// *************ODD TIMESTEP*************//
         timestep++;
 		ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FROM NORMAL
-		ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(NeighborList, dvcMap, fq, Psi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm->RecvD3Q7AA(fq, 0); //WRITE INTO OPPOSITE
+        ScaLBL_DeviceBarrier();
 		// Set boundary conditions
-		/* ... */
-		ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, 0, ScaLBL_Comm->LastExterior(), Np);
-		ScaLBL_Comm->SolidDirichletD3Q7(fq, PoissonSolid);
+		if (BoundaryCondition == 1){
+			ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq,  Vin, timestep);
+			ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep);
+		}
+        //-------------------------//
+		ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(NeighborList, dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np);
+
+        //compute electric field
+        ScaLBL_Comm_Regular->SendHalo(Psi);
+        ScaLBL_D3Q7_Poisson_ElectricField(NeighborList, dvcMap, dvcID, Psi, ElectricField, BoundaryConditionSolid,
+                                          Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+        ScaLBL_Comm_Regular->RecvHalo(Psi);
+        ScaLBL_DeviceBarrier();
+        if (BoundaryCondition == 1){
+            ScaLBL_Comm->Poisson_D3Q7_BC_z(dvcMap,Psi,Vin);
+            ScaLBL_Comm->Poisson_D3Q7_BC_Z(dvcMap,Psi,Vout);
+        }
+        ScaLBL_D3Q7_Poisson_ElectricField(NeighborList, dvcMap, dvcID, Psi, ElectricField, BoundaryConditionSolid, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
+
+        //perform collision
+		ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, 0, ScaLBL_Comm->LastExterior(), Np);
+        if (BoundaryConditionSolid==1){
+		    ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi);
+        }
+        else if (BoundaryConditionSolid==2){
+		    ScaLBL_Comm->SolidNeumannD3Q7(fq, Psi);
+        }
 		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 
 		// *************EVEN TIMESTEP*************//
 		timestep++;
 		ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FORM NORMAL
-		ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(dvcMap, fq, Psi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm->RecvD3Q7AA(fq, 0); //WRITE INTO OPPOSITE
+        ScaLBL_DeviceBarrier();
 		// Set boundary conditions
-		/* ... */
-		ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, 0, ScaLBL_Comm->LastExterior(), Np);
-		ScaLBL_Comm->SolidDirichletD3Q7(fq, PoissonSolid);
+		if (BoundaryCondition == 1){
+			ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq,  Vin, timestep);
+			ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep);
+		}
+        //-------------------------//
+		ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np);
+
+        //compute electric field
+        ScaLBL_Comm_Regular->SendHalo(Psi);
+        ScaLBL_D3Q7_Poisson_ElectricField(NeighborList, dvcMap, dvcID, Psi, ElectricField, BoundaryConditionSolid,
+                                          Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+        ScaLBL_Comm_Regular->RecvHalo(Psi);
+        ScaLBL_DeviceBarrier();
+        if (BoundaryCondition == 1){
+            ScaLBL_Comm->Poisson_D3Q7_BC_z(dvcMap,Psi,Vin);
+            ScaLBL_Comm->Poisson_D3Q7_BC_Z(dvcMap,Psi,Vout);
+        }
+        ScaLBL_D3Q7_Poisson_ElectricField(NeighborList, dvcMap, dvcID, Psi, ElectricField, BoundaryConditionSolid, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
+
+        //perform collision
+		ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, 0, ScaLBL_Comm->LastExterior(), Np);
+        if (BoundaryConditionSolid==1){
+		    ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi);
+        }
+        else if (BoundaryConditionSolid==2){
+		    ScaLBL_Comm->SolidNeumannD3Q7(fq, Psi);
+        }
 		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 		//************************************************************************/
 
         // Check convergence of steady-state solution
         if (timestep%analysis_interval==0){
         
-			ScaLBL_Comm->RegularLayout(Map,Psi,Psi_host);
+			//ScaLBL_Comm->RegularLayout(Map,Psi,Psi_host);
+            ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
 			double count_loc=0;
 			double count;
             double psi_avg;
@@ -393,14 +540,180 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
 
 }
 
-void ScaLBL_Poisson::getElectricalPotential(){
+void ScaLBL_Poisson::getElectricPotential(int timestep){
 
         DoubleArray PhaseField(Nx,Ny,Nz);
-	    ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
+	    //ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
+        ScaLBL_CopyToHost(PhaseField.data(),Psi,sizeof(double)*Nx*Ny*Nz);
         //ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
         FILE *OUTFILE;
-        sprintf(LocalRankFilename,"Electrical_Potential.%05i.raw",rank);
+        sprintf(LocalRankFilename,"Electric_Potential_Time_%i.%05i.raw",timestep,rank);
         OUTFILE = fopen(LocalRankFilename,"wb");
         fwrite(PhaseField.data(),8,N,OUTFILE);
         fclose(OUTFILE);
 }
+
+void ScaLBL_Poisson::getElectricField(int timestep){
+
+        //ScaLBL_D3Q7_Poisson_getElectricField(fq,ElectricField,tau,Np);
+        //ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+
+        DoubleArray PhaseField(Nx,Ny,Nz);
+	    ScaLBL_Comm->RegularLayout(Map,&ElectricField[0*Np],PhaseField);
+        ElectricField_LB_to_Phys(PhaseField);
+        FILE *EX;
+        sprintf(LocalRankFilename,"ElectricField_X_Time_%i.%05i.raw",timestep,rank);
+        EX = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,EX);
+        fclose(EX);
+
+	    ScaLBL_Comm->RegularLayout(Map,&ElectricField[1*Np],PhaseField);
+        ElectricField_LB_to_Phys(PhaseField);
+        FILE *EY;
+        sprintf(LocalRankFilename,"ElectricField_Y_Time_%i.%05i.raw",timestep,rank);
+        EY = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,EY);
+        fclose(EY);
+
+	    ScaLBL_Comm->RegularLayout(Map,&ElectricField[2*Np],PhaseField);
+        ElectricField_LB_to_Phys(PhaseField);
+        FILE *EZ;
+        sprintf(LocalRankFilename,"ElectricField_Z_Time_%i.%05i.raw",timestep,rank);
+        EZ = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,EZ);
+        fclose(EZ);
+}
+
+
+void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){
+	for (int k=0;k<Nz;k++){
+		for (int j=0;j<Ny;j++){
+			for (int i=0;i<Nx;i++){
+                int idx=Map(i,j,k);
+				if (!(idx < 0)){
+                    Efield_reg(i,j,k) = Efield_reg(i,j,k)/(h*1.0e-6); 
+                }
+            }
+        }
+    }
+}
+
+
+//void ScaLBL_Poisson::getElectricPotential(){
+//
+//        DoubleArray PhaseField(Nx,Ny,Nz);
+//	    ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
+//        //ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+//        FILE *OUTFILE;
+//        sprintf(LocalRankFilename,"Electric_Potential.%05i.raw",rank);
+//        OUTFILE = fopen(LocalRankFilename,"wb");
+//        fwrite(PhaseField.data(),8,N,OUTFILE);
+//        fclose(OUTFILE);
+//}
+
+//old version where Psi is of size Np
+//void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid)
+//{
+//	size_t NLABELS=0;
+//	signed char VALUE=0;
+//	double AFFINITY=0.f;
+//
+//	auto LabelList = electric_db->getVector<int>( "SolidLabels" );
+//	auto AffinityList = electric_db->getVector<double>( "SolidValues" );
+//
+//	NLABELS=LabelList.size();
+//	if (NLABELS != AffinityList.size()){
+//		ERROR("Error: LB-Poisson Solver: SolidLabels and SolidValues must be the same length! \n");
+//	}
+//
+//	double label_count[NLABELS];
+//	double label_count_global[NLABELS];
+//	// Assign the labels
+//
+//	for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
+//
+//	for (int k=0;k<Nz;k++){
+//		for (int j=0;j<Ny;j++){
+//			for (int i=0;i<Nx;i++){
+//				int n = k*Nx*Ny+j*Nx+i;
+//				VALUE=Mask->id[n];
+//                AFFINITY=0.f;
+//				// Assign the affinity from the paired list
+//				for (unsigned int idx=0; idx < NLABELS; idx++){
+//				      //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
+//					if (VALUE == LabelList[idx]){
+//						AFFINITY=AffinityList[idx];
+//                        //NOTE need to convert the user input phys unit to LB unit
+//                        if (BoundaryConditionSolid==2){
+//                            //for BCS=1, i.e. Dirichlet-type, no need for unit conversion
+//                            //TODO maybe there is a factor of gamm missing here ?
+//                            AFFINITY = AFFINITY*(h*h*1.0e-12)/epsilon_LB; 
+//                        }
+//						label_count[idx] += 1.0;
+//						idx = NLABELS;
+//						//Mask->id[n] = 0; // set mask to zero since this is an immobile component
+//					}
+//				}
+//				poisson_solid[n] = AFFINITY;
+//			}
+//		}
+//	}
+//
+//	for (size_t idx=0; idx<NLABELS; idx++)
+//		label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
+//
+//	if (rank==0){
+//		printf("LB-Poisson Solver: number of Poisson solid labels: %lu \n",NLABELS);
+//		for (unsigned int idx=0; idx<NLABELS; idx++){
+//			VALUE=LabelList[idx];
+//			AFFINITY=AffinityList[idx];
+//			double volume_fraction  = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
+//            switch (BoundaryConditionSolid){
+//                case 1:
+//			        printf("   label=%d, surface potential=%.3g [V], volume fraction=%.2g\n",VALUE,AFFINITY,volume_fraction); 
+//                    break;
+//                case 2: 
+//			        printf("   label=%d, surface charge density=%.3g [C/m^2], volume fraction=%.2g\n",VALUE,AFFINITY,volume_fraction); 
+//                    break;
+//                default:
+//			        printf("   label=%d, surface potential=%.3g [V], volume fraction=%.2g\n",VALUE,AFFINITY,volume_fraction); 
+//                    break;
+//            }
+//		}
+//	}
+//}
+
+// old version where Psi is of size Np
+//void ScaLBL_Poisson::Potential_Init(double *psi_init){
+//
+//    if (BoundaryCondition==1){
+//	    if (electric_db->keyExists( "Vin" )){
+//	    	Vin = electric_db->getScalar<double>( "Vin" );
+//	    }
+//	    if (electric_db->keyExists( "Vout" )){
+//	    	Vout = electric_db->getScalar<double>( "Vout" );
+//	    }
+//    }
+//    //By default only periodic BC is applied and Vin=Vout=1.0, i.e. there is no potential gradient along Z-axis
+//    double slope = (Vout-Vin)/(Nz-2);
+//    double psi_linearized;
+//	for (int k=0;k<Nz;k++){
+//        if (k==0 || k==1){
+//            psi_linearized = Vin;
+//        }
+//        else if (k==Nz-1 || k==Nz-2){
+//            psi_linearized = Vout;
+//        }
+//        else{
+//            psi_linearized = slope*(k-1)+Vin;
+//        }
+//		for (int j=0;j<Ny;j++){
+//			for (int i=0;i<Nx;i++){
+//				int idx = Map(i,j,k);
+//				if (!(idx < 0)){
+//                    psi_init[idx] = psi_linearized;
+//                }
+//            }
+//        }
+//    }
+//}
diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h
index 9f66c09c..eec17868 100644
--- a/models/PoissonSolver.h
+++ b/models/PoissonSolver.h
@@ -28,7 +28,8 @@ public:
 	void Create();
 	void Initialize();
 	void Run(double *ChargeDensity);
-    void getElectricalPotential();
+    void getElectricPotential(int timestep);
+    void getElectricField(int timestep);
 
 	//bool Restart,pBC;
 	int timestep,timestepMax;
@@ -39,6 +40,7 @@ public:
 	double tolerance;
     double k2_inv,gamma;
     double epsilon0,epsilon0_LB,epsilonR,epsilon_LB;
+    double Vin, Vout;
 	
 	int Nx,Ny,Nz,N,Np;
 	int rank,nprocx,nprocy,nprocz,nprocs;
@@ -48,6 +50,7 @@ public:
 	std::shared_ptr<Domain> Dm;   // this domain is for analysis
 	std::shared_ptr<Domain> Mask; // this domain is for lbm
 	std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
+	std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
     // input database
     std::shared_ptr<Database> db;
     std::shared_ptr<Database> domain_db;
@@ -57,10 +60,12 @@ public:
     DoubleArray Distance;
     DoubleArray Psi_host;
     int *NeighborList;
+    int *dvcMap;
+    signed char *dvcID;
     double *fq;
     double *Psi; 
     double *ElectricField;
-    double *PoissonSolid;
+    //double *PoissonSolid;
 
 private:
 	MPI_Comm comm;
@@ -73,4 +78,6 @@ private:
     //int rank,nprocs;
     void LoadParams(std::shared_ptr<Database> db0);    	
     void AssignSolidBoundary(double *poisson_solid);
+    void Potential_Init(double *psi_init);
+    void ElectricField_LB_to_Phys(DoubleArray &Efield_reg);
 };
diff --git a/models/StokesModel.cpp b/models/StokesModel.cpp
index 41520253..caaf2877 100644
--- a/models/StokesModel.cpp
+++ b/models/StokesModel.cpp
@@ -7,7 +7,7 @@
 
 ScaLBL_StokesModel::ScaLBL_StokesModel(int RANK, int NP, MPI_Comm COMM):
 rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0),
-Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),mu(0),h(0),nu_phys(0),time_conv(0),tolerance(0),
+Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),mu(0),h(0),nu_phys(0),rho_phys(0),rho0(0),den_scale(0),time_conv(0),tolerance(0),
 Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
 {
 
@@ -27,17 +27,17 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
     //-------------------------------------------------------------------//
 
 	//---------------------- Default model parameters --------------------------//		
+    rho_phys = 1000.0; //by default use water density; unit [kg/m^3]
     nu_phys = 1.004e-6;//by default use water kinematic viscosity at 20C; unit [m^2/sec]    
     h = 1.0;//image resolution;[um]
 	tau = 1.0;
     mu = (tau-0.5)/3.0;//LB kinematic viscosity;unit [lu^2/lt]
     time_conv = h*h*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
+    rho0 = 1.0;//LB density
+    den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
 	tolerance = 1.0e-8;
 	Fx = Fy = 0.0;
 	Fz = 1.0e-5;
-    //Body electric field [V/lu]
-    Ex = Ey = 0.0;
-    Ez = 1.0e-3;
     //--------------------------------------------------------------------------//
 
 	// Read domain parameters
@@ -59,19 +59,20 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
 	if (stokes_db->keyExists( "tau" )){
 		tau = stokes_db->getScalar<double>( "tau" );
 	}
+	if (stokes_db->keyExists( "rho0" )){
+		rho0 = stokes_db->getScalar<double>( "rho0" );
+	}
 	if (stokes_db->keyExists( "nu_phys" )){
 		nu_phys = stokes_db->getScalar<double>( "nu_phys" );
 	}
+	if (stokes_db->keyExists( "rho_phys" )){
+		rho_phys = stokes_db->getScalar<double>( "rho_phys" );
+	}
 	if (stokes_db->keyExists( "F" )){
 		Fx = stokes_db->getVector<double>( "F" )[0];
 		Fy = stokes_db->getVector<double>( "F" )[1];
 		Fz = stokes_db->getVector<double>( "F" )[2];
 	}
-	if (stokes_db->keyExists( "ElectricField" )){//NOTE user-input has physical unit [V/m]
-		Ex = stokes_db->getVector<double>( "ElectricField" )[0];
-		Ey = stokes_db->getVector<double>( "ElectricField" )[1];
-		Ez = stokes_db->getVector<double>( "ElectricField" )[2];
-	}
 	if (stokes_db->keyExists( "Restart" )){
 		Restart = stokes_db->getScalar<bool>( "Restart" );
 	}
@@ -88,10 +89,7 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
     // Re-calculate model parameters due to parameter read
 	mu=(tau-0.5)/3.0;
     time_conv = (h*h*1.0e-12)*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
-    // convert user-input electric field ([V/m]) from physical unit to LB unit 
-    Ex = Ex*(h*1.0e-6);//LB electric field: V/lu 
-    Ey = Ey*(h*1.0e-6);
-    Ez = Ez*(h*1.0e-6);
+    den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
 
 	if (rank==0) printf("*****************************************************\n");
 	if (rank==0) printf("LB Single-Fluid Navier-Stokes Solver: \n");
@@ -246,7 +244,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
 	while (timestep < timestepMax) {
         //************************************************************************/
         ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
-        ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez, 
+        ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
                                      ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
         // Set boundary conditions
@@ -262,12 +260,14 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
             ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
             ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
         }
-        ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez, 0, ScaLBL_Comm->LastExterior(), Np);
+        ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv, 
+                                     0, ScaLBL_Comm->LastExterior(), Np);
         ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 
         timestep++;
         ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
-        ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+        ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv, 
+                                      ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
         // Set boundary conditions
         if (BoundaryCondition == 3){
@@ -282,41 +282,87 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
             ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
             ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
         }
-        ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez, 0, ScaLBL_Comm->LastExterior(), Np);
+        ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
+                                      0, ScaLBL_Comm->LastExterior(), Np);
         ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
         //************************************************************************/
     }
 }
 
-void ScaLBL_StokesModel::getVelocity(){
+void ScaLBL_StokesModel::getVelocity(int timestep){
     //get velocity in physical unit [m/sec]
-	ScaLBL_D3Q19_Momentum_Phys(fq, Velocity, h, time_conv, Np);
+	ScaLBL_D3Q19_Momentum(fq, Velocity, Np);
 	ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 
     DoubleArray PhaseField(Nx,Ny,Nz);
 	ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
+    Velocity_LB_to_Phys(PhaseField);
 	FILE *VELX_FILE;
-	sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);
+	sprintf(LocalRankFilename,"Velocity_X_Time_%i.%05i.raw",timestep,rank);
 	VELX_FILE = fopen(LocalRankFilename,"wb");
 	fwrite(PhaseField.data(),8,N,VELX_FILE);
 	fclose(VELX_FILE);
 
 	ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
+    Velocity_LB_to_Phys(PhaseField);
 	FILE *VELY_FILE;
-	sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank);
+	sprintf(LocalRankFilename,"Velocity_Y_Time_%i.%05i.raw",timestep,rank);
 	VELY_FILE = fopen(LocalRankFilename,"wb");
 	fwrite(PhaseField.data(),8,N,VELY_FILE);
 	fclose(VELY_FILE);
 
 	ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
+    Velocity_LB_to_Phys(PhaseField);
 	FILE *VELZ_FILE;
-	sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank);
+	sprintf(LocalRankFilename,"Velocity_Z_Time_%i.%05i.raw",timestep,rank);
 	VELZ_FILE = fopen(LocalRankFilename,"wb");
 	fwrite(PhaseField.data(),8,N,VELZ_FILE);
 	fclose(VELZ_FILE);
 
 }
 
+void ScaLBL_StokesModel::Velocity_LB_to_Phys(DoubleArray &Vel_reg){
+	for (int k=0;k<Nz;k++){
+		for (int j=0;j<Ny;j++){
+			for (int i=0;i<Nx;i++){
+                int idx=Map(i,j,k);
+				if (!(idx < 0)){
+                    Vel_reg(i,j,k) = Vel_reg(i,j,k)*(h*1.0e-6)/time_conv; 
+                }
+            }
+        }
+    }
+}
+
+//void ScaLBL_StokesModel::getVelocity(){
+//    //get velocity in physical unit [m/sec]
+//	ScaLBL_D3Q19_Momentum_Phys(fq, Velocity, h, time_conv, Np);
+//	ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+//
+//    DoubleArray PhaseField(Nx,Ny,Nz);
+//	ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
+//	FILE *VELX_FILE;
+//	sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);
+//	VELX_FILE = fopen(LocalRankFilename,"wb");
+//	fwrite(PhaseField.data(),8,N,VELX_FILE);
+//	fclose(VELX_FILE);
+//
+//	ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
+//	FILE *VELY_FILE;
+//	sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank);
+//	VELY_FILE = fopen(LocalRankFilename,"wb");
+//	fwrite(PhaseField.data(),8,N,VELY_FILE);
+//	fclose(VELY_FILE);
+//
+//	ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
+//	FILE *VELZ_FILE;
+//	sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank);
+//	VELZ_FILE = fopen(LocalRankFilename,"wb");
+//	fwrite(PhaseField.data(),8,N,VELZ_FILE);
+//	fclose(VELZ_FILE);
+//
+//}
+
 void ScaLBL_StokesModel::Run(){
 	double rlx_setA=1.0/tau;
 	double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
diff --git a/models/StokesModel.h b/models/StokesModel.h
index f0a4de6a..d40df415 100644
--- a/models/StokesModel.h
+++ b/models/StokesModel.h
@@ -30,19 +30,21 @@ public:
 	void Run();
 	void Run_Lite(double *ChargeDensity, double *ElectricField);
 	void VelocityField();
-    void getVelocity();
+    void getVelocity(int timestep);
 	
 	bool Restart,pBC;
 	int timestep,timestepMax;
 	int BoundaryCondition;
 	double tau,mu;
+    double rho0;
 	double Fx,Fy,Fz,flux;
-    double Ex,Ey,Ez;
 	double din,dout;
 	double tolerance;
     double nu_phys;
+    double rho_phys;
     double time_conv;
     double h;//image resolution
+    double den_scale;//scale factor for density
 	
 	int Nx,Ny,Nz,N,Np;
 	int rank,nprocx,nprocy,nprocz,nprocs;
@@ -78,4 +80,5 @@ private:
    
     //int rank,nprocs;
     void LoadParams(std::shared_ptr<Database> db0);    	
+    void Velocity_LB_to_Phys(DoubleArray &Vel_reg);
 };
diff --git a/tests/lbpm_electrokinetic_dfh_simulator.cpp b/tests/lbpm_electrokinetic_dfh_simulator.cpp
index 156fbc18..1df5c5e1 100644
--- a/tests/lbpm_electrokinetic_dfh_simulator.cpp
+++ b/tests/lbpm_electrokinetic_dfh_simulator.cpp
@@ -78,14 +78,17 @@ int main(int argc, char **argv)
         while (timestep < Study.timestepMax){
             
             timestep++;
-            if (rank==0) printf("timestep=%i; running Poisson solver\n",timestep);    
+            //if (rank==0) printf("timestep=%i; running Poisson solver\n",timestep);    
             PoissonSolver.Run(IonModel.ChargeDensity);//solve Poisson equtaion to get steady-state electrical potental
+            //PoissonSolver.getElectricPotential(timestep);
 
-            if (rank==0) printf("timestep=%i; running StokesModel\n",timestep);    
+            //if (rank==0) printf("timestep=%i; running StokesModel\n",timestep);    
             StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
+            //StokesModel.getVelocity(timestep);
 
-            if (rank==0) printf("timestep=%i; running Ion model\n",timestep);    
+            //if (rank==0) printf("timestep=%i; running Ion model\n",timestep);    
             IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
+            //IonModel.getIonConcentration(timestep);
             
             
             timestep++;//AA operations
@@ -94,9 +97,10 @@ int main(int argc, char **argv)
             //--------------------------------------------
         }
 
-        StokesModel.getVelocity();
-        PoissonSolver.getElectricalPotential();
-        IonModel.getIonConcentration();
+        StokesModel.getVelocity(timestep);
+        PoissonSolver.getElectricPotential(timestep);
+        PoissonSolver.getElectricField(timestep);
+        IonModel.getIonConcentration(timestep);
 
         if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
         if (rank==0) printf("*************************************************************\n");