add output; CPU version compiled; to be tested

common/ScaLBL.h

@@ -49,6 +49,8 @@ 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
@@ -89,6 +91,8 @@ extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit,
 
 extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np);
 
+extern "C" void ScaLBL_IonConcentration_Phys(double *Den, double h, int ion_component, int start, int finish, int Np);
+
 // 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,
@@ -102,10 +106,10 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(double *dist, 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, int start, int finish, int Np);
+                                              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, int start, int finish, int Np);
+                                             double Gx, double Gy, double Gz, double Ex, double Ey, double Ez, 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,

cpu/Ion.cpp

@@ -235,3 +235,13 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity
     }
 }
 
+extern "C" void ScaLBL_IonConcentration_Phys(double *Den, double h, int ion_component, int start, int finish, int Np){
+    //h: resolution [um/lu]
+    int n;
+    double Ci;
+    
+	for (n=start; n<finish; n++){
+            Ci = Den[n+ion_component*Np];
+            Den[n+ion_component*Np] = Ci/(h*h*h*1.0e-18);
+    }
+}

cpu/Stokes.cpp

@@ -1,7 +1,7 @@
 #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, int start, int finish, int Np)
+                                              double Gx, double Gy, double Gz, double Ex_const, double Ey_const, double Ez_const, int start, int finish, int Np)
 {
     double fq;
 	// conserved momemnts
@@ -32,9 +32,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
         
         //Load data
         rhoE = ChargeDensity[n];
-        Ex = ElectricField[n+0*Np];
+        Ex = ElectricField[n+0*Np]+Ex_const;
-        Ey = ElectricField[n+1*Np];
+        Ey = ElectricField[n+1*Np]+Ey_const;
-        Ez = ElectricField[n+2*Np];
+        Ez = ElectricField[n+2*Np]+Ez_const;
         //compute total body force, including input body force (Gx,Gy,Gz)
         Fx = Gx + rhoE*Ex;
         Fy = Gy + rhoE*Ey;
@@ -455,7 +455,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, int start, int finish, int Np)
+                                             double Gx, double Gy, double Gz, double Ex_const, double Ey_const, double Ez_const, int start, int finish, int Np)
 {
     double fq;
 	// conserved momemnts
@@ -487,9 +487,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
 
         //Load data
         rhoE = ChargeDensity[n];
-        Ex = ElectricField[n+0*Np];
+        Ex = ElectricField[n+0*Np]+Ex_const;
-        Ey = ElectricField[n+1*Np];
+        Ey = ElectricField[n+1*Np]+Ey_const;
-        Ez = ElectricField[n+2*Np];
+        Ez = ElectricField[n+2*Np]+Ez_const;
         //compute total body force, including input body force (Gx,Gy,Gz)
         Fx = Gx + rhoE*Ex;
         Fy = Gy + rhoE*Ey;
@@ -955,3 +955,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;
+		//........................................................................
+	}
+}

models/IonModel.cpp

@@ -457,3 +457,22 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
 
 }
 
+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);
+    }
+
+    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);
+    }
+
+}
+

models/IonModel.h

@@ -29,6 +29,7 @@ public:
 	void Create();
 	void Initialize();
 	void Run(double *Velocity, double *ElectricField);
+    void getIonConcentration();
 	
 	//bool Restart,pBC;
 	int timestep,timestepMax;

models/PoissonSolver.cpp

@@ -390,3 +390,14 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
 
 }
 
+void ScaLBL_Poisson::getElectricalPotential(){
+
+        DoubleArray PhaseField(Nx,Ny,Nz);
+	    ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
+        //ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+        FILE *OUTFILE;
+        sprintf(LocalRankFilename,"Electrical_Potential.%05i.raw",rank);
+        OUTFILE = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE);
+        fclose(OUTFILE);
+}

models/PoissonSolver.h

@@ -28,6 +28,7 @@ public:
 	void Create();
 	void Initialize();
 	void Run(double *ChargeDensity);
+    void getElectricalPotential();
 
 	//bool Restart,pBC;
 	int timestep,timestepMax;

models/StokesModel.cpp

@@ -35,6 +35,9 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
 	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;
     //--------------------------------------------------------------------------//
 
 	// Single-fluid Navier-Stokes Model parameters
@@ -55,6 +58,11 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
 		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" );
 	}
@@ -79,6 +87,11 @@ 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*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);
+
 	if (rank==0) printf("*****************************************************\n");
 	if (rank==0) printf("LB Single-Fluid Navier-Stokes Solver: \n");
 	if (rank==0) printf("      Time conversion factor: %.5g [sec/lt]\n", time_conv);
@@ -232,7 +245,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, 
+        ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez, 
                                      ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
         // Set boundary conditions
@@ -248,12 +261,12 @@ 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, 0, ScaLBL_Comm->LastExterior(), Np);
+        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_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, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+        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_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
         // Set boundary conditions
         if (BoundaryCondition == 3){
@@ -268,16 +281,40 @@ 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, 0, ScaLBL_Comm->LastExterior(), Np);
+        ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez, 0, ScaLBL_Comm->LastExterior(), Np);
         ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
         //************************************************************************/
     }
 }
 
-//void ScaLBL_StokesModel::computeVelocity_phys(){
+void ScaLBL_StokesModel::getVelocity(){
-//	ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
+    //get velocity in physical unit [m/sec]
-//	ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+	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;

models/StokesModel.h

@@ -30,12 +30,14 @@ public:
 	void Run();
 	void Run_Lite(double *ChargeDensity, double *ElectricField);
 	void VelocityField();
+    void getVelocity();
 	
 	bool Restart,pBC;
 	int timestep,timestepMax;
 	int BoundaryCondition;
 	double tau,mu;
 	double Fx,Fy,Fz,flux;
+    double Ex,Ey,Ez;
 	double din,dout;
 	double tolerance;
     double nu_phys;

tests/lbpm_electrokinetic_dfh_simulator.cpp

@@ -89,7 +89,9 @@ int main(int argc, char **argv)
             //--------------------------------------------
         }
 
-        //StokesModel.WriteDebug();
+        StokesModel.getVelocity();
+        PoissonSolver.getElectricalPotential();
+        IonModel.getIonConcentration();
 
         PROFILE_STOP("Main");
         PROFILE_SAVE("lbpm_electrokinetic_simulator",1);