add routines to save advective and electromigrational flux;CPU only;to be verified

analysis/ElectroChemistry.cpp

@@ -10,6 +10,9 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr <Domain> dm):
 	
 	ChemicalPotential.resize(Nx,Ny,Nz);       ChemicalPotential.fill(0);
 	ElectricalPotential.resize(Nx,Ny,Nz);      ElectricalPotential.fill(0);
+	ElectricalField_x.resize(Nx,Ny,Nz);      ElectricalField_x.fill(0);
+	ElectricalField_y.resize(Nx,Ny,Nz);      ElectricalField_y.fill(0);
+	ElectricalField_z.resize(Nx,Ny,Nz);      ElectricalField_z.fill(0);
 	Pressure.resize(Nx,Ny,Nz);     	Pressure.fill(0);
 	Rho.resize(Nx,Ny,Nz);       	Rho.fill(0);
 	Vel_x.resize(Nx,Ny,Nz);         Vel_x.fill(0);	    // Gradient of the phase indicator field
@@ -169,13 +172,17 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
     visData[0].meshName = "domain";
     visData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
     auto ElectricPotential = std::make_shared<IO::Variable>();
+
+    //ion concentration
     std::vector<shared_ptr<IO::Variable>> IonConcentration;
     for (size_t ion=0; ion<Ion.number_ion_species; ion++){
         IonConcentration.push_back(std::make_shared<IO::Variable>());
     }
+    //fluid velocity
     auto VxVar = std::make_shared<IO::Variable>();
     auto VyVar = std::make_shared<IO::Variable>();
     auto VzVar = std::make_shared<IO::Variable>();
+    // diffusive ion flux
     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
@@ -183,6 +190,22 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
         IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
         IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
     }
+    // advective ion flux
+    std::vector<shared_ptr<IO::Variable>> IonFluxAdvective;
+    for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+        //push in x-,y-, and z-component for each ion species
+        IonFluxAdvective.push_back(std::make_shared<IO::Variable>());
+        IonFluxAdvective.push_back(std::make_shared<IO::Variable>());
+        IonFluxAdvective.push_back(std::make_shared<IO::Variable>());
+    }
+    // electro-migrational ion flux
+    std::vector<shared_ptr<IO::Variable>> IonFluxElectrical;
+    for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+        //push in x-,y-, and z-component for each ion species
+        IonFluxElectrical.push_back(std::make_shared<IO::Variable>());
+        IonFluxElectrical.push_back(std::make_shared<IO::Variable>());
+        IonFluxElectrical.push_back(std::make_shared<IO::Variable>());
+    }
     //--------------------------------------------------------------------------------------------------------------------
 
     //-------------------------------------Create Names for Variables------------------------------------------------------
@@ -248,6 +271,58 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
     		visData[0].vars.push_back(IonFluxDiffusive[3*ion+2]);
     	}
     }
+
+    if (vis_db->getWithDefault<bool>( "save_ion_flux_advective", false )){
+    	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+            // x-component of advective flux
+    		sprintf(VisName,"Ion%zu_FluxAdvective_x",ion+1);
+    		IonFluxAdvective[3*ion+0]->name = VisName;
+    		IonFluxAdvective[3*ion+0]->type = IO::VariableType::VolumeVariable;
+    		IonFluxAdvective[3*ion+0]->dim = 1;
+    		IonFluxAdvective[3*ion+0]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    		visData[0].vars.push_back(IonFluxAdvective[3*ion+0]);
+            // y-component of advective flux
+    		sprintf(VisName,"Ion%zu_FluxAdvective_y",ion+1);
+    		IonFluxAdvective[3*ion+1]->name = VisName;
+    		IonFluxAdvective[3*ion+1]->type = IO::VariableType::VolumeVariable;
+    		IonFluxAdvective[3*ion+1]->dim = 1;
+    		IonFluxAdvective[3*ion+1]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    		visData[0].vars.push_back(IonFluxAdvective[3*ion+1]);
+            // z-component of advective flux
+    		sprintf(VisName,"Ion%zu_FluxAdvective_z",ion+1);
+    		IonFluxAdvective[3*ion+2]->name = VisName;
+    		IonFluxAdvective[3*ion+2]->type = IO::VariableType::VolumeVariable;
+    		IonFluxAdvective[3*ion+2]->dim = 1;
+    		IonFluxAdvective[3*ion+2]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    		visData[0].vars.push_back(IonFluxAdvective[3*ion+2]);
+    	}
+    }
+
+    if (vis_db->getWithDefault<bool>( "save_ion_flux_electrical", false )){
+    	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+            // x-component of electro-migrational flux
+    		sprintf(VisName,"Ion%zu_FluxElectrical_x",ion+1);
+    		IonFluxElectrical[3*ion+0]->name = VisName;
+    		IonFluxElectrical[3*ion+0]->type = IO::VariableType::VolumeVariable;
+    		IonFluxElectrical[3*ion+0]->dim = 1;
+    		IonFluxElectrical[3*ion+0]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    		visData[0].vars.push_back(IonFluxElectrical[3*ion+0]);
+            // y-component of electro-migrational flux
+    		sprintf(VisName,"Ion%zu_FluxElectrical_y",ion+1);
+    		IonFluxElectrical[3*ion+1]->name = VisName;
+    		IonFluxElectrical[3*ion+1]->type = IO::VariableType::VolumeVariable;
+    		IonFluxElectrical[3*ion+1]->dim = 1;
+    		IonFluxElectrical[3*ion+1]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    		visData[0].vars.push_back(IonFluxElectrical[3*ion+1]);
+            // z-component of electro-migrational flux
+    		sprintf(VisName,"Ion%zu_FluxElectrical_z",ion+1);
+    		IonFluxElectrical[3*ion+2]->name = VisName;
+    		IonFluxElectrical[3*ion+2]->type = IO::VariableType::VolumeVariable;
+    		IonFluxElectrical[3*ion+2]->dim = 1;
+    		IonFluxElectrical[3*ion+2]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    		visData[0].vars.push_back(IonFluxElectrical[3*ion+2]);
+    	}
+    }
     //--------------------------------------------------------------------------------------------------------------------
     
     //------------------------------------Save All Variables--------------------------------------------------------------
@@ -307,7 +382,59 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
     		fillData.copy(IonFluxDiffusive_z,IonFluxData_z);
     	}
     }
+
+    if (vis_db->getWithDefault<bool>( "save_ion_flux_advective", false )){
+    	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+            
+            // x-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxAdvective_x",ion+1);
+    		//IonFluxDiffusive[3*ion+0]->name = VisName;
+    		ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+0]->name==VisName);
+            // y-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxAdvective_y",ion+1);
+    		//IonFluxDiffusive[3*ion+1]->name = VisName;
+    		ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+1]->name==VisName);
+            // z-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxAdvective_z",ion+1);
+    		//IonFluxDiffusive[3*ion+2]->name = VisName;
+    		ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+2]->name==VisName);
+
+    		Array<double>& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+0]->data;
+    		Array<double>& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+1]->data;
+    		Array<double>& IonFluxData_z = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+2]->data;
+    		Ion.getIonFluxAdvective(IonFluxAdvective_x,IonFluxAdvective_y,IonFluxAdvective_z,ion);
+    		fillData.copy(IonFluxAdvective_x,IonFluxData_x);
+    		fillData.copy(IonFluxAdvective_y,IonFluxData_y);
+    		fillData.copy(IonFluxAdvective_z,IonFluxData_z);
+    	}
+    }
     
+    if (vis_db->getWithDefault<bool>( "save_ion_flux_electrical", false )){
+    	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+            
+            // x-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxElectrical_x",ion+1);
+    		//IonFluxDiffusive[3*ion+0]->name = VisName;
+    		ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+0]->name==VisName);
+            // y-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxElectrical_y",ion+1);
+    		//IonFluxDiffusive[3*ion+1]->name = VisName;
+    		ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+1]->name==VisName);
+            // z-component of diffusive flux
+    		sprintf(VisName,"Ion%zu_FluxElectrical_z",ion+1);
+    		//IonFluxDiffusive[3*ion+2]->name = VisName;
+    		ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+2]->name==VisName);
+
+    		Array<double>& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+0]->data;
+    		Array<double>& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+1]->data;
+    		Array<double>& IonFluxData_z = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+2]->data;
+    		Ion.getIonFluxElectrical(IonFluxElectrical_x,IonFluxElectrical_y,IonFluxElectrical_z,ion);
+    		fillData.copy(IonFluxElectrical_x,IonFluxData_x);
+    		fillData.copy(IonFluxElectrical_y,IonFluxData_y);
+    		fillData.copy(IonFluxElectrical_z,IonFluxData_z);
+    	}
+    }
+
     if (vis_db->getWithDefault<bool>( "write_silo", true ))
     	IO::writeData( timestep, visData, Dm->Comm );
     //--------------------------------------------------------------------------------------------------------------------

analysis/ElectroChemistry.h

@@ -35,6 +35,9 @@ public:
 	DoubleArray Rho;	           // density field 
 	DoubleArray ChemicalPotential;	   // density field 
 	DoubleArray ElectricalPotential;	// density field 
+	DoubleArray ElectricalField_x;	// density field 
+	DoubleArray ElectricalField_y;	// density field 
+	DoubleArray ElectricalField_z;	// density field 
 	DoubleArray Pressure; 	// pressure field
 	DoubleArray Vel_x;		// velocity field
 	DoubleArray Vel_y;

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 *FluxDiffusive, double *Velocity, double *ElectricField, 
+extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, 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 *FluxDiffusive, double *Velocity, double *ElectricField, 
+extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, 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);

cpu/Ion.cpp

@@ -80,7 +80,7 @@ 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 *FluxDiffusive, double *Velocity, double *ElectricField, 
+extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField, 
                                       double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
 	int n;
 	double Ci;
@@ -133,6 +133,12 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
         FluxDiffusive[n+0*Np] = flux_diffusive_x;
         FluxDiffusive[n+1*Np] = flux_diffusive_y;
         FluxDiffusive[n+2*Np] = flux_diffusive_z;
+        FluxAdvective[n+0*Np] = ux*Ci;
+        FluxAdvective[n+1*Np] = uy*Ci;
+        FluxAdvective[n+2*Np] = uz*Ci;
+        FluxElectrical[n+0*Np] = uEPx*Ci;
+        FluxElectrical[n+1*Np] = uEPy*Ci;
+        FluxElectrical[n+2*Np] = uEPz*Ci;
 		
 		// q=0
 		dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
@@ -158,7 +164,7 @@ 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 *FluxDiffusive, double *Velocity, double *ElectricField, 
+extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField, 
                                        double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
 	int n;
 	double Ci;
@@ -197,6 +203,12 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDi
         FluxDiffusive[n+0*Np] = flux_diffusive_x;
         FluxDiffusive[n+1*Np] = flux_diffusive_y;
         FluxDiffusive[n+2*Np] = flux_diffusive_z;
+        FluxAdvective[n+0*Np] = ux*Ci;
+        FluxAdvective[n+1*Np] = uy*Ci;
+        FluxAdvective[n+2*Np] = uz*Ci;
+        FluxElectrical[n+0*Np] = uEPx*Ci;
+        FluxElectrical[n+1*Np] = uEPy*Ci;
+        FluxElectrical[n+2*Np] = uEPz*Ci;
 
 		// q=0
 		dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;

models/IonModel.cpp

@@ -692,7 +692,9 @@ 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);
+	ScaLBL_AllocateDeviceMemory((void **) &FluxDiffusive,  number_ion_species*3*sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &FluxAdvective,  number_ion_species*3*sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &FluxElectrical, 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");
@@ -878,9 +880,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
             
 
             //LB-Ion collison
-            ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
+            ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[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],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
+            ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
                                   rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
             
             if (BoundaryConditionSolid==1){
@@ -934,9 +936,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
             
 
             //LB-Ion collison
-            ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
+            ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[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],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
+            ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
                                   rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
             
             if (BoundaryConditionSolid==1){
@@ -992,6 +994,38 @@ void ScaLBL_IonModel::getIonFluxDiffusive(DoubleArray &IonFlux_x,DoubleArray &Io
 	ScaLBL_Comm->Barrier(); comm.barrier();
 }
 
+void ScaLBL_IonModel::getIonFluxAdvective(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,&FluxAdvective[ic*3*Np+0*Np],IonFlux_x);
+	IonFlux_LB_to_Phys(IonFlux_x,ic);
+	ScaLBL_Comm->Barrier(); comm.barrier();
+
+	ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+1*Np],IonFlux_y);
+	IonFlux_LB_to_Phys(IonFlux_y,ic);
+	ScaLBL_Comm->Barrier(); comm.barrier();
+
+	ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+2*Np],IonFlux_z);
+	IonFlux_LB_to_Phys(IonFlux_z,ic);
+	ScaLBL_Comm->Barrier(); comm.barrier();
+}
+
+void ScaLBL_IonModel::getIonFluxElectrical(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,&FluxElectrical[ic*3*Np+0*Np],IonFlux_x);
+	IonFlux_LB_to_Phys(IonFlux_x,ic);
+	ScaLBL_Comm->Barrier(); comm.barrier();
+
+	ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+1*Np],IonFlux_y);
+	IonFlux_LB_to_Phys(IonFlux_y,ic);
+	ScaLBL_Comm->Barrier(); comm.barrier();
+
+	ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[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){
     //This function write out decomposed data
     DoubleArray PhaseField(Nx,Ny,Nz);
@@ -1048,6 +1082,85 @@ void ScaLBL_IonModel::getIonFluxDiffusive_debug(int timestep){
     }
 }
 
+void ScaLBL_IonModel::getIonFluxAdvective_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,&FluxAdvective[ic*3*Np+0*Np],PhaseField);
+        ScaLBL_Comm->Barrier(); comm.barrier();
+        IonFlux_LB_to_Phys(PhaseField,ic);
+
+        FILE *OUTFILE_X;
+        sprintf(LocalRankFilename,"IonFluxAdvective_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE_X = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE_X);
+        fclose(OUTFILE_X);
+
+        //y-component
+	    ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+1*Np],PhaseField);
+        ScaLBL_Comm->Barrier(); comm.barrier();
+        IonFlux_LB_to_Phys(PhaseField,ic);
+
+        FILE *OUTFILE_Y;
+        sprintf(LocalRankFilename,"IonFluxAdvective_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE_Y = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE_Y);
+        fclose(OUTFILE_Y);
+
+        //z-component
+	    ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+2*Np],PhaseField);
+        ScaLBL_Comm->Barrier(); comm.barrier();
+        IonFlux_LB_to_Phys(PhaseField,ic);
+
+        FILE *OUTFILE_Z;
+        sprintf(LocalRankFilename,"IonFluxAdvective_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE_Z = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE_Z);
+        fclose(OUTFILE_Z);
+    }
+}
+
+void ScaLBL_IonModel::getIonFluxElectrical_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,&FluxElectrical[ic*3*Np+0*Np],PhaseField);
+        ScaLBL_Comm->Barrier(); comm.barrier();
+        IonFlux_LB_to_Phys(PhaseField,ic);
+
+        FILE *OUTFILE_X;
+        sprintf(LocalRankFilename,"IonFluxElectrical_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE_X = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE_X);
+        fclose(OUTFILE_X);
+
+        //y-component
+	    ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+1*Np],PhaseField);
+        ScaLBL_Comm->Barrier(); comm.barrier();
+        IonFlux_LB_to_Phys(PhaseField,ic);
+
+        FILE *OUTFILE_Y;
+        sprintf(LocalRankFilename,"IonFluxElectrical_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE_Y = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE_Y);
+        fclose(OUTFILE_Y);
+
+        //z-component
+	    ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+2*Np],PhaseField);
+        ScaLBL_Comm->Barrier(); comm.barrier();
+        IonFlux_LB_to_Phys(PhaseField,ic);
+
+        FILE *OUTFILE_Z;
+        sprintf(LocalRankFilename,"IonFluxElectrical_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE_Z = 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++){

models/IonModel.h

@@ -86,6 +86,8 @@ public:
     double *FluidVelocityDummy;
     double *ElectricFieldDummy;
     double *FluxDiffusive;
+    double *FluxAdvective;
+    double *FluxElectrical;
 
 private:
 	Utilities::MPI comm;

tests/TestIonModel.cpp

@@ -76,6 +76,8 @@ int main(int argc, char **argv)
         }
         IonModel.getIonConcentration_debug(timestep);
         IonModel.getIonFluxDiffusive_debug(timestep);
+        IonModel.getIonFluxAdvective_debug(timestep);
+        IonModel.getIonFluxElectrical_debug(timestep);
 
         if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
         if (rank==0) printf("*************************************************************\n");