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

update PoissonSolver and fix numerous bugs

Browse Source
This commit is contained in:
Rex Zhe Li 2020-09-02 11:37:23 -04:00
parent 86a1bb81a1
commit 20c8cc9c3b
9 changed files with 305 additions and 228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
common/Domain.cpp
View File

@ -626,38 +626,38 @@ void Domain::Decomp( const std::string& Filename )
if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6)); if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
//......................................................... //.........................................................
// If external boundary conditions are applied remove solid // If external boundary conditions are applied remove solid
if (BoundaryCondition > 0 && BoundaryCondition !=5 && kproc() == 0){ // if (BoundaryCondition > 0 && BoundaryCondition !=5 && kproc() == 0){
if (inlet_layers_z < 4){ // if (inlet_layers_z < 4){
inlet_layers_z=4; // inlet_layers_z=4;
if(RANK==0){ // if(RANK==0){
printf("NOTE:Non-periodic BC is applied, but the number of Z-inlet layers is not specified (or is smaller than 3 voxels) \n the number of Z-inlet layer is reset to %i voxels, saturated with phase label=%i \n",inlet_layers_z-1,inlet_layers_phase); // printf("NOTE:Non-periodic BC is applied, but the number of Z-inlet layers is not specified (or is smaller than 3 voxels) \n the number of Z-inlet layer is reset to %i voxels, saturated with phase label=%i \n",inlet_layers_z-1,inlet_layers_phase);
} // }
} // }
for (int k=0; k<inlet_layers_z; k++){ // for (int k=0; k<inlet_layers_z; k++){
for (int j=0;j<Ny;j++){ // for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){ // for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i; // int n = k*Nx*Ny+j*Nx+i;
id[n] = inlet_layers_phase; // id[n] = inlet_layers_phase;
} // }
} // }
} // }
} // }
if (BoundaryCondition > 0 && BoundaryCondition !=5 && kproc() == nprocz-1){ // if (BoundaryCondition > 0 && BoundaryCondition !=5 && kproc() == nprocz-1){
if (outlet_layers_z < 4){ // if (outlet_layers_z < 4){
outlet_layers_z=4; // outlet_layers_z=4;
if(RANK==nprocs-1){ // if(RANK==nprocs-1){
printf("NOTE:Non-periodic BC is applied, but the number of Z-outlet layers is not specified (or is smaller than 3 voxels) \n the number of Z-outlet layer is reset to %i voxels, saturated with phase label=%i \n",outlet_layers_z-1,outlet_layers_phase); // printf("NOTE:Non-periodic BC is applied, but the number of Z-outlet layers is not specified (or is smaller than 3 voxels) \n the number of Z-outlet layer is reset to %i voxels, saturated with phase label=%i \n",outlet_layers_z-1,outlet_layers_phase);
} // }
} // }
for (int k=Nz-outlet_layers_z; k<Nz; k++){ // for (int k=Nz-outlet_layers_z; k<Nz; k++){
for (int j=0;j<Ny;j++){ // for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){ // for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i; // int n = k*Nx*Ny+j*Nx+i;
id[n] = outlet_layers_phase; // id[n] = outlet_layers_phase;
} // }
} // }
} // }
} // }
for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-1;k++){ for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-1;k++){
for (int j=1;j<Ny-1;j++){ for (int j=1;j<Ny-1;j++){
for (int i=1;i<Nx-1;i++){ for (int i=1;i<Nx-1;i++){

4
common/ScaLBL.cpp
View File

@ -1443,8 +1443,8 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component){
//...Packing for y face(4,8,9,16,18)................................. //...Packing for y face(4,8,9,16,18).................................
ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,&Aq[Component*7*N],N); ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,&Aq[Component*7*N],N);
MPI_Isend(sendbuf_y, 2*sendCount_y,MPI_DOUBLE,rank_y,sendtag,MPI_COMM_SCALBL,&req1[2]); MPI_Isend(sendbuf_y, sendCount_y,MPI_DOUBLE,rank_y,sendtag,MPI_COMM_SCALBL,&req1[2]);
MPI_Irecv(recvbuf_Y, 2*recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,MPI_COMM_SCALBL,&req2[2]); MPI_Irecv(recvbuf_Y, recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,MPI_COMM_SCALBL,&req2[2]);
//...Packing for Y face(3,7,10,15,17)................................. //...Packing for Y face(3,7,10,15,17).................................
ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,&Aq[Component*7*N],N); ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,&Aq[Component*7*N],N);

4
common/ScaLBL.h
View File

@ -92,10 +92,10 @@ extern "C" void ScaLBL_IonConcentration_Phys(double *Den, double h, int ion_comp
// LBM Poisson solver // LBM Poisson solver
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, 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,
int start, int finish, int Np); int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, 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,
int start, int finish, int Np); 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); extern "C" void ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(int *neighborList,int *Map, double *dist, double *Psi, int start, int finish, int Np);

2
cpu/D3Q7BC.cpp
View File

@ -11,7 +11,7 @@ extern "C" void ScaLBL_Solid_Dirichlet_D3Q7(double *dist,double *BoundaryValue,i
ib = BounceBackSolid_list[idx]; ib = BounceBackSolid_list[idx];
value_b = BoundaryValue[ib];//get boundary value from a solid site value_b = BoundaryValue[ib];//get boundary value from a solid site
value_q = dist[iq]; value_q = dist[iq];
dist[iq] = -1.0*value_q + value_b*2.0/9.0;//NOTE 2/9 is the speed of sound for D3Q7 lattice dist[iq] = -1.0*value_q + value_b*0.25;//NOTE 0.25 is the speed of sound for D3Q7 lattice
} }
} }

129
cpu/Poisson.cpp
View File

@ -88,10 +88,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
} }
} }
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double tau, double epsilon_LB,double gamma,int start, int finish, int Np){ 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,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
double rho_e;//local charge density double rho_e;//local charge density
double f0,f1,f2,f3,f4,f5,f6; double f0,f1,f2,f3,f4,f5,f6;
int nr1,nr2,nr3,nr4,nr5,nr6; int nr1,nr2,nr3,nr4,nr5,nr6;
@ -102,7 +102,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
//Load data //Load data
rho_e = Den_charge[n]; rho_e = Den_charge[n];
rho_e = gamma*rho_e/epsilon_LB; rho_e = rho_e/epsilon_LB;
idx=Map[n]; idx=Map[n];
psi = Psi[idx]; psi = Psi[idx];
@ -131,51 +131,41 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
nr6 = neighborList[n+5*Np]; nr6 = neighborList[n+5*Np];
f6 = dist[nr6]; f6 = dist[nr6];
//Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
//Ey = (f3-f4)*rlx*4.5; Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice speed of sound
//Ez = (f5-f6)*rlx*4.5; Ez = (f5-f6)*rlx*4.0;
//Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu ElectricField[n+0*Np] = Ex;
//Ey = (f3-f4)*rlx*4.0; ElectricField[n+1*Np] = Ey;
//Ez = (f5-f6)*rlx*4.0; ElectricField[n+2*Np] = Ez;
//ElectricField[n+0*Np] = Ex;
//ElectricField[n+1*Np] = Ey;
//ElectricField[n+2*Np] = Ez;
// q = 0 // q = 0
dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e); dist[n] = f0*(1.0-rlx) + 0.25*(rlx*psi+rho_e);
//dist[n] = f0*(1.0-rlx) + 0.25*(rlx*psi+rho_e);
// q = 1 // q = 1
dist[nr2] = f1*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); dist[nr2] = f1*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//dist[nr2] = f1*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 2 // q = 2
dist[nr1] = f2*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); dist[nr1] = f2*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//dist[nr1] = f2*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 3 // q = 3
dist[nr4] = f3*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); dist[nr4] = f3*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//dist[nr4] = f3*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 4 // q = 4
dist[nr3] = f4*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); dist[nr3] = f4*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//dist[nr3] = f4*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 5 // q = 5
dist[nr6] = f5*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); dist[nr6] = f5*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//dist[nr6] = f5*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 6 // q = 6
dist[nr5] = f6*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e); dist[nr5] = f6*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//dist[nr5] = f6*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//........................................................................ //........................................................................
} }
} }
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, 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,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
double rho_e;//local charge density double rho_e;//local charge density
double f0,f1,f2,f3,f4,f5,f6; double f0,f1,f2,f3,f4,f5,f6;
double rlx=1.0/tau; double rlx=1.0/tau;
@ -185,7 +175,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_c
//Load data //Load data
rho_e = Den_charge[n]; rho_e = Den_charge[n];
rho_e = gamma*rho_e/epsilon_LB; rho_e = rho_e/epsilon_LB;
idx=Map[n]; idx=Map[n];
psi = Psi[idx]; psi = Psi[idx];
@ -198,43 +188,33 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_c
f6 = dist[5*Np+n]; f6 = dist[5*Np+n];
//Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
//Ey = (f3-f4)*rlx*4.5; Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice speed of sound
//Ez = (f5-f6)*rlx*4.5; Ez = (f5-f6)*rlx*4.0;
//Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu ElectricField[n+0*Np] = Ex;
//Ey = (f3-f4)*rlx*4.0; ElectricField[n+1*Np] = Ey;
//Ez = (f5-f6)*rlx*4.0; ElectricField[n+2*Np] = Ez;
//ElectricField[n+0*Np] = Ex;
//ElectricField[n+1*Np] = Ey;
//ElectricField[n+2*Np] = Ez;
// q = 0 // q = 0
dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e); dist[n] = f0*(1.0-rlx) + 0.25*(rlx*psi+rho_e);
//dist[n] = f0*(1.0-rlx) + 0.25*(rlx*psi+rho_e);
// q = 1 // 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.125*(rlx*psi+rho_e);
//dist[1*Np+n] = f1*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 2 // 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.125*(rlx*psi+rho_e);
//dist[2*Np+n] = f2*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 3 // 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.125*(rlx*psi+rho_e);
//dist[3*Np+n] = f3*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 4 // 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.125*(rlx*psi+rho_e);
//dist[4*Np+n] = f4*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 5 // 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.125*(rlx*psi+rho_e);
//dist[5*Np+n] = f5*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 6 // 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.125*(rlx*psi+rho_e);
//dist[6*Np+n] = f6*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//........................................................................ //........................................................................
} }
} }
@ -245,35 +225,13 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, in
int ijk; int ijk;
for (n=start; n<finish; n++){ for (n=start; n<finish; n++){
ijk = Map[n]; ijk = Map[n];
//dist[0*Np+n] = 0.3333333333333333*Psi[n]; dist[0*Np+n] = 0.25*Psi[ijk];
//dist[1*Np+n] = 0.1111111111111111*Psi[n]; dist[1*Np+n] = 0.125*Psi[ijk];
//dist[2*Np+n] = 0.1111111111111111*Psi[n]; dist[2*Np+n] = 0.125*Psi[ijk];
//dist[3*Np+n] = 0.1111111111111111*Psi[n]; dist[3*Np+n] = 0.125*Psi[ijk];
//dist[4*Np+n] = 0.1111111111111111*Psi[n]; dist[4*Np+n] = 0.125*Psi[ijk];
//dist[5*Np+n] = 0.1111111111111111*Psi[n]; dist[5*Np+n] = 0.125*Psi[ijk];
//dist[6*Np+n] = 0.1111111111111111*Psi[n]; dist[6*Np+n] = 0.125*Psi[ijk];
//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.3333333333333333*Psi[ijk];
dist[1*Np+n] = 0.1111111111111111*Psi[ijk];
dist[2*Np+n] = 0.1111111111111111*Psi[ijk];
dist[3*Np+n] = 0.1111111111111111*Psi[ijk];
dist[4*Np+n] = 0.1111111111111111*Psi[ijk];
dist[5*Np+n] = 0.1111111111111111*Psi[ijk];
dist[6*Np+n] = 0.1111111111111111*Psi[ijk];
//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];
} }
} }
@ -368,9 +326,12 @@ extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, s
id = ID[nn]; id = ID[nn];
m18 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 18 m18 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 18
//............Compute the Color Gradient................................... //............Compute the Color Gradient...................................
nx = -1.f/18.f*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); //nx = 1.f/6.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)); //ny = 1.f/6.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)); //nz = 1.f/6.f*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
nx = 1.f/6.f*(m1-m2);//but looks like it needs to multiply another factor of 3
ny = 1.f/6.f*(m3-m4);
nz = 1.f/6.f*(m5-m6);
ElectricField[n] = nx; ElectricField[n] = nx;
ElectricField[Np+n] = ny; ElectricField[Np+n] = ny;

209
models/PoissonSolver.cpp
View File

@ -6,8 +6,9 @@
#include "common/ReadMicroCT.h" #include "common/ReadMicroCT.h"
ScaLBL_Poisson::ScaLBL_Poisson(int RANK, int NP, MPI_Comm COMM): 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), rank(RANK), nprocs(NP),timestep(0),timestepMax(0),tau(0),k2_inv(0),tolerance(0),h(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), 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),
chargeDen_dummy(0),
nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),BoundaryConditionSolid(0),Lx(0),Ly(0),Lz(0),comm(COMM) nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),BoundaryConditionSolid(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{ {
@ -22,10 +23,8 @@ void ScaLBL_Poisson::ReadParams(string filename){
domain_db = db->getDatabase( "Domain" ); domain_db = db->getDatabase( "Domain" );
electric_db = db->getDatabase( "Poisson" ); electric_db = db->getDatabase( "Poisson" );
k2_inv = 4.5;//speed of sound for D3Q7 lattice k2_inv = 4.0;//speed of sound for D3Q7 lattice
//k2_inv = 4.0;//speed of sound for D3Q7 lattice tau = 0.5+k2_inv;
gamma = 1.0;//time step of LB-Poisson equation
tau = 0.5+k2_inv*gamma;
timestepMax = 100000; timestepMax = 100000;
tolerance = 1.0e-6;//stopping criterion for obtaining steady-state electricla potential tolerance = 1.0e-6;//stopping criterion for obtaining steady-state electricla potential
h = 1.0;//resolution; unit: um/lu h = 1.0;//resolution; unit: um/lu
@ -36,6 +35,7 @@ void ScaLBL_Poisson::ReadParams(string filename){
analysis_interval = 1000; analysis_interval = 1000;
Vin = 1.0; //Boundary-z (inlet) electric potential Vin = 1.0; //Boundary-z (inlet) electric potential
Vout = 1.0; //Boundary-Z (outlet) electric potential Vout = 1.0; //Boundary-Z (outlet) electric potential
chargeDen_dummy = 1.0e-3;//For debugging;unit=[C/m^3]
// LB-Poisson Model parameters // LB-Poisson Model parameters
if (electric_db->keyExists( "timestepMax" )){ if (electric_db->keyExists( "timestepMax" )){
@ -47,12 +47,12 @@ void ScaLBL_Poisson::ReadParams(string filename){
if (electric_db->keyExists( "tolerance" )){ if (electric_db->keyExists( "tolerance" )){
tolerance = electric_db->getScalar<double>( "tolerance" ); tolerance = electric_db->getScalar<double>( "tolerance" );
} }
if (electric_db->keyExists( "gamma" )){
gamma = electric_db->getScalar<double>( "gamma" );
}
if (electric_db->keyExists( "epsilonR" )){ if (electric_db->keyExists( "epsilonR" )){
epsilonR = electric_db->getScalar<double>( "epsilonR" ); epsilonR = electric_db->getScalar<double>( "epsilonR" );
} }
if (electric_db->keyExists( "DummyChargeDen" )){
chargeDen_dummy = electric_db->getScalar<double>( "DummyChargeDen" );
}
// Read solid boundary condition specific to Poisson equation // Read solid boundary condition specific to Poisson equation
BoundaryConditionSolid = 1; BoundaryConditionSolid = 1;
@ -76,7 +76,6 @@ void ScaLBL_Poisson::ReadParams(string filename){
//Re-calcualte model parameters if user updates input //Re-calcualte model parameters if user updates input
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)] epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity
tau = 0.5+k2_inv*gamma;
if (rank==0) printf("***********************************************************************************\n"); if (rank==0) printf("***********************************************************************************\n");
if (rank==0) printf("LB-Poisson Solver: steady-state MaxTimeStep = %i; steady-state tolerance = %.3g \n", timestepMax,tolerance); if (rank==0) printf("LB-Poisson Solver: steady-state MaxTimeStep = %i; steady-state tolerance = %.3g \n", timestepMax,tolerance);
@ -213,7 +212,6 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid)
//NOTE need to convert the user input phys unit to LB unit //NOTE need to convert the user input phys unit to LB unit
if (BoundaryConditionSolid==2){ if (BoundaryConditionSolid==2){
//for BCS=1, i.e. Dirichlet-type, no need for unit conversion //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; AFFINITY = AFFINITY*(h*h*1.0e-12)/epsilon_LB;
} }
label_count[idx] += 1.0; label_count[idx] += 1.0;
@ -284,11 +282,10 @@ void ScaLBL_Poisson::Create(){
//........................................................................... //...........................................................................
ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np); ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
ScaLBL_AllocateDeviceMemory((void **) &dvcID, sizeof(signed char)*Nx*Ny*Nz); //ScaLBL_AllocateDeviceMemory((void **) &dvcID, sizeof(signed char)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &fq, 7*dist_mem_size); ScaLBL_AllocateDeviceMemory((void **) &fq, 7*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz); ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np);
//ScaLBL_AllocateDeviceMemory((void **) &PoissonSolid, sizeof(double)*Nx*Ny*Nz);
//........................................................................... //...........................................................................
// Update GPU data structures // Update GPU data structures
@ -329,26 +326,13 @@ void ScaLBL_Poisson::Create(){
MPI_Barrier(comm); MPI_Barrier(comm);
delete [] neighborList; delete [] neighborList;
// copy node ID // copy node ID
ScaLBL_CopyToDevice(dvcID, Mask->id, sizeof(signed char)*Nx*Ny*Nz); //ScaLBL_CopyToDevice(dvcID, Mask->id, sizeof(signed char)*Nx*Ny*Nz);
ScaLBL_DeviceBarrier(); //ScaLBL_DeviceBarrier();
//Initialize solid boundary for electric potential //Initialize solid boundary for electric potential
ScaLBL_Comm->SetupBounceBackList(Map, Mask->id, Np); ScaLBL_Comm->SetupBounceBackList(Map, Mask->id, Np);
MPI_Barrier(comm); 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;
} }
// 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){ void ScaLBL_Poisson::Potential_Init(double *psi_init){
@ -402,6 +386,16 @@ void ScaLBL_Poisson::Initialize(){
ScaLBL_D3Q7_Poisson_Init(dvcMap, fq, Psi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); 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); ScaLBL_D3Q7_Poisson_Init(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np);
delete [] psi_host; delete [] psi_host;
//extra treatment for halo layer
if (BoundaryCondition==1){
if (Dm->kproc()==0){
ScaLBL_SetSlice_z(Psi,Vin,Nx,Ny,Nz,0);
}
if (Dm->kproc() == nprocz-1){
ScaLBL_SetSlice_z(Psi,Vout,Nx,Ny,Nz,Nz-1);
}
}
} }
void ScaLBL_Poisson::Run(double *ChargeDensity){ void ScaLBL_Poisson::Run(double *ChargeDensity){
@ -418,20 +412,17 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
//************************************************************************/ //************************************************************************/
// *************ODD TIMESTEP*************// // *************ODD TIMESTEP*************//
timestep++; timestep++;
SolveElectricPotentialAAodd();
//compute electric field SolveElectricPotentialAAodd();//update electric potential
SolveElectricField(); //SolveElectricField(); //deprecated - compute electric field
//perform collision SolvePoissonAAodd(ChargeDensity);//perform collision
SolvePoissonAAodd(ChargeDensity);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm); ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// *************EVEN TIMESTEP*************// // *************EVEN TIMESTEP*************//
timestep++; timestep++;
SolveElectricPotentialAAeven(); SolveElectricPotentialAAeven();//update electric potential
//compute electric field //SolveElectricField();//deprecated - compute electric field
SolveElectricField(); SolvePoissonAAeven(ChargeDensity);//perform collision
//perform collision
SolvePoissonAAeven(ChargeDensity);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm); ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/ //************************************************************************/
@ -484,6 +475,75 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
} }
void ScaLBL_Poisson::SolveElectricPotentialAAodd(){
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FROM NORMAL
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
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);
}
void ScaLBL_Poisson::SolveElectricPotentialAAeven(){
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FORM NORMAL
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
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);
}
void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){
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, 0, ScaLBL_Comm->LastExterior(), Np);
if (BoundaryConditionSolid==1){
ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi);
}
else if (BoundaryConditionSolid==2){
ScaLBL_Comm->SolidNeumannD3Q7(fq, Psi);
}
}
void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity){
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, 0, ScaLBL_Comm->LastExterior(), Np);
if (BoundaryConditionSolid==1){
ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi);
}
else if (BoundaryConditionSolid==2){
ScaLBL_Comm->SolidNeumannD3Q7(fq, Psi);
}
}
void ScaLBL_Poisson::DummyChargeDensity(){
double *ChargeDensity_host;
ChargeDensity_host = new double[Np];
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))
ChargeDensity_host[idx] = chargeDen_dummy*(h*h*h*1.0e-18);
}
}
}
ScaLBL_AllocateDeviceMemory((void **) &ChargeDensityDummy, sizeof(double)*Np);
ScaLBL_CopyToDevice(ChargeDensityDummy, ChargeDensity_host, sizeof(double)*Np);
ScaLBL_DeviceBarrier();
delete [] ChargeDensity_host;
}
void ScaLBL_Poisson::getElectricPotential(int timestep){ void ScaLBL_Poisson::getElectricPotential(int timestep){
DoubleArray PhaseField(Nx,Ny,Nz); DoubleArray PhaseField(Nx,Ny,Nz);
@ -528,68 +588,6 @@ void ScaLBL_Poisson::getElectricField(int timestep){
fclose(EZ); fclose(EZ);
} }
void ScaLBL_Poisson::SolveElectricPotentialAAodd(){
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FROM NORMAL
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
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);
}
void ScaLBL_Poisson::SolveElectricField(){
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);
}
void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, 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);
}
}
void ScaLBL_Poisson::SolveElectricPotentialAAeven(){
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FORM NORMAL
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
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);
}
void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity){
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, 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);
}
}
void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){ void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){
for (int k=0;k<Nz;k++){ for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){ for (int j=0;j<Ny;j++){
@ -603,6 +601,19 @@ void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){
} }
} }
//void ScaLBL_Poisson::SolveElectricField(){
// 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);
//
//}
//void ScaLBL_Poisson::getElectricPotential(){ //void ScaLBL_Poisson::getElectricPotential(){
// //

18
models/PoissonSolver.h
View File

@ -28,13 +28,9 @@ public:
void Create(); void Create();
void Initialize(); void Initialize();
void Run(double *ChargeDensity); void Run(double *ChargeDensity);
void SolveElectricPotentialAAodd();
void SolveElectricPotentialAAeven();
void SolveElectricField();
void SolvePoissonAAodd(double *ChargeDensity);
void SolvePoissonAAeven(double *ChargeDensity);
void getElectricPotential(int timestep); void getElectricPotential(int timestep);
void getElectricField(int timestep); void getElectricField(int timestep);
void DummyChargeDensity();//for debugging
//bool Restart,pBC; //bool Restart,pBC;
int timestep,timestepMax; int timestep,timestepMax;
@ -43,9 +39,10 @@ public:
int BoundaryConditionSolid; int BoundaryConditionSolid;
double tau; double tau;
double tolerance; double tolerance;
double k2_inv,gamma; double k2_inv;
double epsilon0,epsilon0_LB,epsilonR,epsilon_LB; double epsilon0,epsilon0_LB,epsilonR,epsilon_LB;
double Vin, Vout; double Vin, Vout;
double chargeDen_dummy;//for debugging
int Nx,Ny,Nz,N,Np; int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs; int rank,nprocx,nprocy,nprocz,nprocs;
@ -66,11 +63,11 @@ public:
DoubleArray Psi_host; DoubleArray Psi_host;
int *NeighborList; int *NeighborList;
int *dvcMap; int *dvcMap;
signed char *dvcID; //signed char *dvcID;
double *fq; double *fq;
double *Psi; double *Psi;
double *ElectricField; double *ElectricField;
//double *PoissonSolid; double *ChargeDensityDummy;// for debugging
private: private:
MPI_Comm comm; MPI_Comm comm;
@ -85,5 +82,10 @@ private:
void AssignSolidBoundary(double *poisson_solid); void AssignSolidBoundary(double *poisson_solid);
void Potential_Init(double *psi_init); void Potential_Init(double *psi_init);
void ElectricField_LB_to_Phys(DoubleArray &Efield_reg); void ElectricField_LB_to_Phys(DoubleArray &Efield_reg);
void SolveElectricPotentialAAodd();
void SolveElectricPotentialAAeven();
//void SolveElectricField();
void SolvePoissonAAodd(double *ChargeDensity);
void SolvePoissonAAeven(double *ChargeDensity);
}; };

1
tests/CMakeLists.txt
View File

@ -47,6 +47,7 @@ ADD_LBPM_TEST( TestTorusEvolve )
ADD_LBPM_TEST( TestTopo3D ) ADD_LBPM_TEST( TestTopo3D )
ADD_LBPM_TEST( TestFluxBC ) ADD_LBPM_TEST( TestFluxBC )
ADD_LBPM_TEST( TestMap ) ADD_LBPM_TEST( TestMap )
ADD_LBPM_TEST( TestPoissonSolver )
#ADD_LBPM_TEST( TestMRT ) #ADD_LBPM_TEST( TestMRT )
#ADD_LBPM_TEST( TestColorGrad ) #ADD_LBPM_TEST( TestColorGrad )
#ADD_LBPM_TEST( TestColorGradDFH ) #ADD_LBPM_TEST( TestColorGradDFH )

102
tests/TestPoissonSolver.cpp Normal file
View File

@ -0,0 +1,102 @@
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include <math.h>
#include "models/PoissonSolver.h"
using namespace std;
//********************************************************
// Test lattice-Boltzmann solver of Poisson equation
//********************************************************
int main(int argc, char **argv)
{
// Initialize MPI
int provided_thread_support = -1;
MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
MPI_Comm comm;
MPI_Comm_dup(MPI_COMM_WORLD,&comm);
int rank = comm_rank(comm);
int nprocs = comm_size(comm);
if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
}
// Limit scope so variables that contain communicators will free before MPI_Finialize
{
if (rank == 0){
printf("********************************************************\n");
printf("Running Test for LB-Poisson Solver \n");
printf("********************************************************\n");
}
//PROFILE_ENABLE_TRACE();
//PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START("Main");
Utilities::setErrorHandlers();
auto filename = argv[1];
ScaLBL_Poisson PoissonSolver(rank,nprocs,comm);
// Initialize LB-Poisson model
PoissonSolver.ReadParams(filename);
PoissonSolver.SetDomain();
PoissonSolver.ReadInput();
PoissonSolver.Create();
PoissonSolver.Initialize();
PoissonSolver.getElectricPotential(0);
//Initialize dummy charge density for test
PoissonSolver.DummyChargeDensity();
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy);
PoissonSolver.getElectricPotential(1);
PoissonSolver.getElectricField(1);
//int timestep=0;
//while (timestep < Study.timestepMax){
//
// 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);
// 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);
// IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
// //IonModel.getIonConcentration(timestep);
//
//
// timestep++;//AA operations
// //--------------------------------------------
// //potentially leave analysis module for future
// //--------------------------------------------
//}
//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");
PROFILE_STOP("Main");
PROFILE_SAVE("TestPoissonSolver",1);
// ****************************************************
MPI_Barrier(comm);
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Comm_free(&comm);
MPI_Finalize();
}