diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp
index 1a9a986a..b46a7473 100644
--- a/analysis/ElectroChemistry.cpp
+++ b/analysis/ElectroChemistry.cpp
@@ -176,8 +176,16 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
     auto VxVar = std::make_shared<IO::Variable>();
     auto VyVar = std::make_shared<IO::Variable>();
     auto VzVar = std::make_shared<IO::Variable>();
+    std::vector<shared_ptr<IO::Variable>> IonFluxDiffusive;
+    for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+        //push in x-,y-, and z-component for each ion species
+        IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
+        IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
+        IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
+    }
+    //--------------------------------------------------------------------------------------------------------------------
 
-    
+    //-------------------------------------Create Names for Variables------------------------------------------------------
     if (vis_db->getWithDefault<bool>( "save_electric_potential", true )){
     	ElectricPotential->name = "ElectricPotential";
     	ElectricPotential->type = IO::VariableType::VolumeVariable;
@@ -214,7 +222,35 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
         VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
         visData[0].vars.push_back(VzVar);
     }
+
+    if (vis_db->getWithDefault<bool>( "save_ion_flux_diffusive", false )){
+    	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+            // x-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxDiffusive_x",ion+1);
+    		IonFluxDiffusive[3*ion+0]->name = VisName;
+    		IonFluxDiffusive[3*ion+0]->type = IO::VariableType::VolumeVariable;
+    		IonFluxDiffusive[3*ion+0]->dim = 1;
+    		IonFluxDiffusive[3*ion+0]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    		visData[0].vars.push_back(IonFluxDiffusive[3*ion+0]);
+            // y-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxDiffusive_y",ion+1);
+    		IonFluxDiffusive[3*ion+1]->name = VisName;
+    		IonFluxDiffusive[3*ion+1]->type = IO::VariableType::VolumeVariable;
+    		IonFluxDiffusive[3*ion+1]->dim = 1;
+    		IonFluxDiffusive[3*ion+1]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    		visData[0].vars.push_back(IonFluxDiffusive[3*ion+1]);
+            // z-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxDiffusive_z",ion+1);
+    		IonFluxDiffusive[3*ion+2]->name = VisName;
+    		IonFluxDiffusive[3*ion+2]->type = IO::VariableType::VolumeVariable;
+    		IonFluxDiffusive[3*ion+2]->dim = 1;
+    		IonFluxDiffusive[3*ion+2]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    		visData[0].vars.push_back(IonFluxDiffusive[3*ion+2]);
+    	}
+    }
+    //--------------------------------------------------------------------------------------------------------------------
     
+    //------------------------------------Save All Variables--------------------------------------------------------------
     if (vis_db->getWithDefault<bool>( "save_electric_potential", true )){
     	ASSERT(visData[0].vars[0]->name=="ElectricPotential");
     	Poisson.getElectricPotential(ElectricalPotential);
@@ -225,7 +261,7 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
     if (vis_db->getWithDefault<bool>( "save_concentration", true )){
     	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
     		sprintf(VisName,"IonConcentration_%zu",ion+1);
-    		IonConcentration[ion]->name = VisName;
+    		//IonConcentration[ion]->name = VisName;
     		ASSERT(visData[0].vars[1+ion]->name==VisName);
     		Array<double>& IonConcentrationData = visData[0].vars[1+ion]->data;
     		Ion.getIonConcentration(Rho,ion);
@@ -245,10 +281,36 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
     	fillData.copy(Vel_y,VelyData);
     	fillData.copy(Vel_z,VelzData);
     }
+
+    if (vis_db->getWithDefault<bool>( "save_ion_flux_diffusive", false )){
+    	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+            
+            // x-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxDiffusive_x",ion+1);
+    		//IonFluxDiffusive[3*ion+0]->name = VisName;
+    		ASSERT(visData[0].vars[4+Ion.numer_ion_species+3*ion+0]->name==VisName);
+            // y-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxDiffusive_y",ion+1);
+    		//IonFluxDiffusive[3*ion+1]->name = VisName;
+    		ASSERT(visData[0].vars[4+Ion.numer_ion_species+3*ion+1]->name==VisName);
+            // z-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxDiffusive_z",ion+1);
+    		//IonFluxDiffusive[3*ion+2]->name = VisName;
+    		ASSERT(visData[0].vars[4+Ion.numer_ion_species+3*ion+2]->name==VisName);
+
+    		Array<double>& IonFluxData_x = visData[0].vars[4+Ion.numer_ion_species+3*ion+0]->data;
+    		Array<double>& IonFluxData_y = visData[0].vars[4+Ion.numer_ion_species+3*ion+1]->data;
+    		Array<double>& IonFluxData_z = visData[0].vars[4+Ion.numer_ion_species+3*ion+2]->data;
+    		Ion.getIonFluxDiffusive(IonFluxDiffusive_x,IonFluxDiffusive_y,IonFluxDiffusive_z,ion);
+    		fillData.copy(IonFluxDiffusive_x,IonFluxData_x);
+    		fillData.copy(IonFluxDiffusive_y,IonFluxData_y);
+    		fillData.copy(IonFluxDiffusive_z,IonFluxData_z);
+    	}
+    }
     
     if (vis_db->getWithDefault<bool>( "write_silo", true ))
     	IO::writeData( timestep, visData, Dm->Comm );
-
+    //--------------------------------------------------------------------------------------------------------------------
 /*    if (vis_db->getWithDefault<bool>( "save_8bit_raw", true )){
     	char CurrentIDFilename[40];
     	sprintf(CurrentIDFilename,"id_t%d.raw",timestep);
diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index d37723c2..6e15b9f1 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -253,10 +253,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *di
 extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, int start, int finish, int Np);
 
 
-extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField, 
+extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, 
                                       double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
 
-extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField, 
+extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, 
                                        double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
 
 extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np);
diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp
index 5d7ef87a..672ad71c 100644
--- a/cpu/Ion.cpp
+++ b/cpu/Ion.cpp
@@ -80,13 +80,14 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i
     }
 }
 
-extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField, 
+extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, 
                                       double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
 	int n;
 	double Ci;
     double ux,uy,uz;
     double uEPx,uEPy,uEPz;//electrochemical induced velocity
     double Ex,Ey,Ez;//electrical field
+    double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
 	double f0,f1,f2,f3,f4,f5,f6;
 	int nr1,nr2,nr3,nr4,nr5,nr6;
 
@@ -124,6 +125,14 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
 		// q=6
 		nr6 = neighborList[n+5*Np];
 		f6 = dist[nr6];
+
+        // compute diffusive flux
+		flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
+		flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
+		flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
+        FluxDiffusive[n+0*Np] = flux_diffusive_x;
+        FluxDiffusive[n+1*Np] = flux_diffusive_y;
+        FluxDiffusive[n+2*Np] = flux_diffusive_z;
 		
 		// q=0
 		dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
@@ -149,13 +158,14 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
 	}
 }
 
-extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField, 
+extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, 
                                        double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
 	int n;
 	double Ci;
     double ux,uy,uz;
     double uEPx,uEPy,uEPz;//electrochemical induced velocity
     double Ex,Ey,Ez;//electrical field
+    double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
 	double f0,f1,f2,f3,f4,f5,f6;
 
 	for (n=start; n<finish; n++){
@@ -180,6 +190,14 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Veloci
 		f5 = dist[6*Np+n];
 		f6 = dist[5*Np+n];
 		
+        // compute diffusive flux
+		flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
+		flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
+		flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
+        FluxDiffusive[n+0*Np] = flux_diffusive_x;
+        FluxDiffusive[n+1*Np] = flux_diffusive_y;
+        FluxDiffusive[n+2*Np] = flux_diffusive_z;
+
 		// q=0
 		dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
 
diff --git a/cpu/Poisson.cpp b/cpu/Poisson.cpp
index 503c02ad..b508072f 100644
--- a/cpu/Poisson.cpp
+++ b/cpu/Poisson.cpp
@@ -132,7 +132,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
 		f6 = dist[nr6];
 		
 		Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
-		Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice speed of sound
+		Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice squared speed of sound
 		Ez = (f5-f6)*rlx*4.0;
         ElectricField[n+0*Np] = Ex;
         ElectricField[n+1*Np] = Ey;
@@ -189,7 +189,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_c
 
 
 		Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
-		Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice speed of sound
+		Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice squared speed of sound
 		Ez = (f5-f6)*rlx*4.0;
         ElectricField[n+0*Np] = Ex;
         ElectricField[n+1*Np] = Ey;
diff --git a/models/IonModel.cpp b/models/IonModel.cpp
index b14279a3..1f08186a 100644
--- a/models/IonModel.cpp
+++ b/models/IonModel.cpp
@@ -692,6 +692,7 @@ void ScaLBL_IonModel::Create(){
 	ScaLBL_AllocateDeviceMemory((void **) &fq, number_ion_species*7*dist_mem_size);  
 	ScaLBL_AllocateDeviceMemory((void **) &Ci, number_ion_species*sizeof(double)*Np);
 	ScaLBL_AllocateDeviceMemory((void **) &ChargeDensity, sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &FluxDiffusive, number_ion_species*3*sizeof(double)*Np);
 	//...........................................................................
 	// Update GPU data structures
 	if (rank==0)    printf ("LB Ion Solver: Setting up device map and neighbor list \n");
@@ -877,9 +878,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
             
 
             //LB-Ion collison
-            ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
+            ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
                                   rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
-            ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
+            ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
                                   rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
             
             if (BoundaryConditionSolid==1){
@@ -933,9 +934,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
             
 
             //LB-Ion collison
-            ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
+            ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
                                   rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
-            ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
+            ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
                                   rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
             
             if (BoundaryConditionSolid==1){
@@ -973,7 +974,22 @@ void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const s
 	ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration);
 	ScaLBL_Comm->Barrier(); comm.barrier();
 	IonConcentration_LB_to_Phys(IonConcentration);
+}
 
+void ScaLBL_IonModel::getIonFluxDiffusive(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic){
+	//This function wirte out the data in a normal layout (by aggregating all decomposed domains)
+
+	ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+0*Np],IonFlux_x);
+	IonFlux_LB_to_Phys(IonFlux_x,ic);
+	ScaLBL_Comm->Barrier(); comm.barrier();
+
+	ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+1*Np],IonFlux_y);
+	IonFlux_LB_to_Phys(IonFlux_y,ic);
+	ScaLBL_Comm->Barrier(); comm.barrier();
+
+	ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+2*Np],IonFlux_z);
+	IonFlux_LB_to_Phys(IonFlux_z,ic);
+	ScaLBL_Comm->Barrier(); comm.barrier();
 }
 
 void ScaLBL_IonModel::getIonConcentration_debug(int timestep){
@@ -992,13 +1008,67 @@ void ScaLBL_IonModel::getIonConcentration_debug(int timestep){
     }
 }
 
+void ScaLBL_IonModel::getIonFluxDiffusive_debug(int timestep){
+    //This function write out decomposed data
+
+    DoubleArray PhaseField(Nx,Ny,Nz);
+	for (size_t ic=0; ic<number_ion_species; ic++){
+        //x-component
+	    ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+0*Np],PhaseField);
+        ScaLBL_Comm->Barrier(); comm.barrier();
+        IonFlux_LB_to_Phys(PhaseField);
+
+        FILE *OUTFILE_X;
+        sprintf(LocalRankFilename,"IonFluxDiffusive_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE_X);
+        fclose(OUTFILE_X);
+
+        //y-component
+	    ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+1*Np],PhaseField);
+        ScaLBL_Comm->Barrier(); comm.barrier();
+        IonFlux_LB_to_Phys(PhaseField);
+
+        FILE *OUTFILE_Y;
+        sprintf(LocalRankFilename,"IonFluxDiffusive_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE_Y);
+        fclose(OUTFILE_Y);
+
+        //z-component
+	    ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+2*Np],PhaseField);
+        ScaLBL_Comm->Barrier(); comm.barrier();
+        IonFlux_LB_to_Phys(PhaseField);
+
+        FILE *OUTFILE_Z;
+        sprintf(LocalRankFilename,"IonFluxDiffusive_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE_Z);
+        fclose(OUTFILE_Z);
+    }
+}
+
+
 void ScaLBL_IonModel::IonConcentration_LB_to_Phys(DoubleArray &Den_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)){
-                    Den_reg(i,j,k) = Den_reg(i,j,k)/(h*h*h*1.0e-18); 
+                    Den_reg(i,j,k) = Den_reg(i,j,k)/(h*h*h*1.0e-18);//this converts the unit to [mol/m^3] 
+                }
+            }
+        }
+    }
+}
+
+void ScaLBL_IonModel::IonFlux_LB_to_Phys(DoubleArray &Den_reg, const size_t ic){
+	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)){
+                    Den_reg(i,j,k) = Den_reg(i,j,k)/(h*h*1.0e-12*time_conv[ic]);//this converts the unit to [mol/m^2/s] 
                 }
             }
         }
diff --git a/models/IonModel.h b/models/IonModel.h
index 6e2a71aa..33e3eb4b 100644
--- a/models/IonModel.h
+++ b/models/IonModel.h
@@ -83,6 +83,7 @@ public:
     double *IonSolid;
     double *FluidVelocityDummy;
     double *ElectricFieldDummy;
+    double *FluxDiffusive;
 
 private:
 	Utilities::MPI comm;
diff --git a/tests/TestIonModel.cpp b/tests/TestIonModel.cpp
index 58c051dc..622a0e8f 100644
--- a/tests/TestIonModel.cpp
+++ b/tests/TestIonModel.cpp
@@ -77,6 +77,7 @@ int main(int argc, char **argv)
             }
         }
         IonModel.getIonConcentration_debug(timestep);
+        IonModel.getIonFluxDiffusive_debug(timestep);
 
         if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
         if (rank==0) printf("*************************************************************\n");