OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

update Poisson and Stoke solvers to include slipping velocity BC; to be built

Browse Source
This commit is contained in:
Rex Zhe Li
2022-01-18 16:46:09 +11:00
parent 35bea97b2e
commit ee527a42cc
7 changed files with 58 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
common/ScaLBL.h
View File

@@ -296,7 +296,7 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity
* @param finish - lattice node to finish loop * @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure) * @param Np - size of local sub-domain (derived from Domain structure)
*/ */
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, 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, bool UseSlippingVelBC,
int start, int finish, int Np); int start, int finish, int Np);
/** /**
@@ -307,19 +307,22 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *di
* @param Psi - * @param Psi -
* @param ElectricField - electric field * @param ElectricField - electric field
* @param tau - relaxation time * @param tau - relaxation time
* @param epsilon_LB - * @param epsilon_LB - dielectric constant of medium
* @param UseSlippingVelBC - flag indicating the use of Helmholtz-Smoluchowski slipping velocity equation when EDL is too small to resolve
* @param start - lattice node to start loop * @param start - lattice node to start loop
* @param finish - lattice node to finish loop * @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure) * @param Np - size of local sub-domain (derived from Domain structure)
*/ */
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau,
double epsilon_LB, int start, int finish, int Np); double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np);
/** /**
* \brief Poisson-Boltzmann solver / solve electric potential based on AA odd access pattern for D3Q7 * \brief Poisson-Boltzmann solver / solve electric potential based on AA odd access pattern for D3Q7
* @param neighborList - neighbors based on D3Q19 lattice structure * @param neighborList - neighbors based on D3Q19 lattice structure
* @param Map - mapping between sparse and dense representations * @param Map - mapping between sparse and dense representations
* @param dist - D3Q7 distributions * @param dist - D3Q7 distributions
* @param Psi - * @param Psi -
* @param epsilon_LB - dielectric constant of medium
* @param UseSlippingVelBC - flag indicating the use of Helmholtz-Smoluchowski slipping velocity equation when EDL is too small to resolve
* @param start - lattice node to start loop * @param start - lattice node to start loop
* @param finish - lattice node to finish loop * @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure) * @param Np - size of local sub-domain (derived from Domain structure)
@@ -350,10 +353,10 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, in
// LBM Stokes Model (adapted from MRT model) // 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, 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 Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, 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, 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 Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np);
extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np); extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np);

20
cpu/Poisson.cpp
View File

@@ -95,7 +95,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
double *dist, double *Den_charge, double *dist, double *Den_charge,
double *Psi, double *ElectricField, double *Psi, double *ElectricField,
double tau, double epsilon_LB, double tau, double epsilon_LB, bool UseSlippingVelBC,
int start, int finish, int Np) { int start, int finish, int Np) {
int n; int n;
double psi; //electric potential double psi; //electric potential
@@ -109,8 +109,11 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
for (n = start; n < finish; n++) { for (n = start; n < finish; n++) {
//Load data //Load data
rho_e = Den_charge[n]; //rho_e = Den_charge[n];
rho_e = rho_e / epsilon_LB; //rho_e = rho_e / epsilon_LB;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
idx = Map[n]; idx = Map[n];
psi = Psi[idx]; psi = Psi[idx];
@@ -175,8 +178,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist,
double *Den_charge, double *Psi, double *Den_charge, double *Psi,
double *ElectricField, double tau, double *ElectricField, double tau,
double epsilon_LB, int start, double epsilon_LB, bool UseSlippingVelBC,
int finish, int Np) { int start, int finish, int Np) {
int n; int n;
double psi; //electric potential double psi; //electric potential
double Ex, Ey, Ez; //electric field double Ex, Ey, Ez; //electric field
@@ -188,8 +191,11 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist,
for (n = start; n < finish; n++) { for (n = start; n < finish; n++) {
//Load data //Load data
rho_e = Den_charge[n]; //rho_e = Den_charge[n];
rho_e = rho_e / epsilon_LB; //rho_e = rho_e / epsilon_LB;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
idx = Map[n]; idx = Map[n];
psi = Psi[idx]; psi = Psi[idx];

31
cpu/Stokes.cpp
View File

@@ -4,7 +4,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
double *dist, double *Velocity, double *ChargeDensity, double *dist, double *Velocity, double *ChargeDensity,
double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double *ElectricField, double rlx_setA, double rlx_setB, double Gx,
double Gy, double Gz, double rho0, double den_scale, double h, double Gy, double Gz, double rho0, double den_scale, double h,
double time_conv, int start, int finish, int Np) { double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) {
double fq; double fq;
// conserved momemnts // conserved momemnts
double rho, jx, jy, jz; double rho, jx, jy, jz;
@@ -38,13 +38,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
Ey = ElectricField[n + 1 * Np]; Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np]; Ez = ElectricField[n + 2 * Np];
//compute total body force, including input body force (Gx,Gy,Gz) //compute total body force, including input body force (Gx,Gy,Gz)
Fx = Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
Gx + den_scale; //the extra factors at the end necessarily convert unit from phys to LB
rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (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; den_scale;
Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale; den_scale;
// q=0 // q=0
@@ -479,7 +477,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(
int *neighborList, double *dist, double *Velocity, double *ChargeDensity, int *neighborList, double *dist, double *Velocity, double *ChargeDensity,
double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double *ElectricField, double rlx_setA, double rlx_setB, double Gx,
double Gy, double Gz, double rho0, double den_scale, double h, double Gy, double Gz, double rho0, double den_scale, double h,
double time_conv, int start, int finish, int Np) { double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) {
double fq; double fq;
// conserved momemnts // conserved momemnts
double rho, jx, jy, jz; double rho, jx, jy, jz;
@@ -513,12 +511,21 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(
Ex = ElectricField[n + 0 * Np]; Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np]; Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np]; Ez = ElectricField[n + 2 * Np];
//compute total body force, including input body force (Gx,Gy,Gz) //compute total body force, including input body force (Gx,Gy,Gz)
Fx = Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / //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;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and body force induced by external efectric field is reduced to slipping velocity BC.
Fx = (UseSlippingVelBC==1) ? Gx : 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 = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale; den_scale;
Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (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; den_scale;
// q=0 // q=0

18
models/PoissonSolver.cpp
View File

@@ -523,7 +523,7 @@ void ScaLBL_Poisson::Initialize(double time_conv_from_Study){
//} //}
} }
void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study){
//.......create and start timer............ //.......create and start timer............
//double starttime,stoptime,cputime; //double starttime,stoptime,cputime;
@@ -537,13 +537,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
// *************ODD TIMESTEP*************// // *************ODD TIMESTEP*************//
timestep++; timestep++;
SolveElectricPotentialAAodd(timestep_from_Study);//update electric potential SolveElectricPotentialAAodd(timestep_from_Study);//update electric potential
SolvePoissonAAodd(ChargeDensity);//perform collision SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier(); ScaLBL_Comm->Barrier(); comm.barrier();
// *************EVEN TIMESTEP*************// // *************EVEN TIMESTEP*************//
timestep++; timestep++;
SolveElectricPotentialAAeven(timestep_from_Study);//update electric potential SolveElectricPotentialAAeven(timestep_from_Study);//update electric potential
SolvePoissonAAeven(ChargeDensity);//perform collision SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier(); ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/ //************************************************************************/
@@ -724,9 +724,9 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){
ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np);
} }
void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC){
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np);
//TODO: perhaps add another ScaLBL_Comm routine to update Psi values on solid boundary nodes. //TODO: perhaps add another ScaLBL_Comm routine to update Psi values on solid boundary nodes.
//something like: //something like:
//ScaLBL_Comm->SolidDirichletBoundaryUpdates(Psi, Psi_BCLabel, timestep); //ScaLBL_Comm->SolidDirichletBoundaryUpdates(Psi, Psi_BCLabel, timestep);
@@ -739,9 +739,9 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){
//} //}
} }
void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity){ void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC){
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->SolidDirichletAndNeumannD3Q7(fq, Psi, Psi_BCLabel); ScaLBL_Comm->SolidDirichletAndNeumannD3Q7(fq, Psi, Psi_BCLabel);
//if (BoundaryConditionSolid==1){ //if (BoundaryConditionSolid==1){
// ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi); // ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi);

2
models/PoissonSolver.h
View File

@@ -34,7 +34,7 @@ public:
void ReadInput(); void ReadInput();
void Create(); void Create();
void Initialize(double time_conv_from_Study); void Initialize(double time_conv_from_Study);
void Run(double *ChargeDensity, int timestep_from_Study); void Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study);
void getElectricPotential(DoubleArray &ReturnValues); void getElectricPotential(DoubleArray &ReturnValues);
void getElectricPotential_debug(int timestep); void getElectricPotential_debug(int timestep);
void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y, void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y,

14
models/StokesModel.cpp
View File

@@ -593,7 +593,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_StokesMRT( ScaLBL_D3Q19_AAodd_StokesMRT(
NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA,
rlx_setB, Fx, Fy, Fz, rho0, den_scale, h, time_conv, rlx_setB, Fx, Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions // Set boundary conditions
@@ -610,8 +610,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
} }
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity,
ElectricField, rlx_setA, rlx_setB, Fx, Fy, ElectricField, rlx_setA, rlx_setB, Fx, Fy,
Fz, rho0, den_scale, h, time_conv, 0, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
ScaLBL_Comm->LastExterior(), Np); 0, ScaLBL_Comm->LastExterior(), Np);
if (UseSlippingVelBC == true) { if (UseSlippingVelBC == true) {
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19( ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(
@@ -626,8 +626,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_StokesMRT( ScaLBL_D3Q19_AAeven_StokesMRT(
fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx,
Fy, Fz, rho0, den_scale, h, time_conv, ScaLBL_Comm->FirstInterior(), Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions // Set boundary conditions
if (BoundaryCondition == 3) { if (BoundaryCondition == 3) {
@@ -643,8 +643,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
} }
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity,
ElectricField, rlx_setA, rlx_setB, Fx, Fy, ElectricField, rlx_setA, rlx_setB, Fx, Fy,
Fz, rho0, den_scale, h, time_conv, 0, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
ScaLBL_Comm->LastExterior(), Np); 0, ScaLBL_Comm->LastExterior(), Np);
if (UseSlippingVelBC == true) { if (UseSlippingVelBC == true) {
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19( ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(
fq, ZetaPotentialSolid, ElectricField, SolidGrad, epsilon_LB, fq, ZetaPotentialSolid, ElectricField, SolidGrad, epsilon_LB,

2
tests/lbpm_electrokinetic_SingleFluid_simulator.cpp
View File

@@ -94,7 +94,7 @@ int main(int argc, char **argv)
while (timestep < Study.timestepMax){ while (timestep < Study.timestepMax){
timestep++; timestep++;
PoissonSolver.Run(IonModel.ChargeDensity,timestep);//solve Poisson equtaion to get steady-state electrical potental PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential