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

fix a few trivial bugs; add some checkpoint print; still debugging

Browse Source
This commit is contained in:
Rex Zhe Li
2020-08-18 12:40:41 -04:00
parent 771f679f5c
commit 5756d6f138
6 changed files with 65 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
common/ScaLBL.cpp
View File

@@ -1429,37 +1429,37 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component){
ScaLBL_DeviceBarrier();
// Pack the distributions
//...Packing for x face(2,8,10,12,14)................................
ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,&Aq[Component*N],N);
ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,&Aq[Component*7*N],N);
MPI_Isend(sendbuf_x, sendCount_x,MPI_DOUBLE,rank_x,sendtag,MPI_COMM_SCALBL,&req1[0]);
MPI_Irecv(recvbuf_X, recvCount_X,MPI_DOUBLE,rank_X,recvtag,MPI_COMM_SCALBL,&req2[0]);
//...Packing for X face(1,7,9,11,13)................................
ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,&Aq[Component*N],N);
ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,&Aq[Component*7*N],N);
MPI_Isend(sendbuf_X, sendCount_X,MPI_DOUBLE,rank_X,sendtag,MPI_COMM_SCALBL,&req1[1]);
MPI_Irecv(recvbuf_x, recvCount_x,MPI_DOUBLE,rank_x,recvtag,MPI_COMM_SCALBL,&req2[1]);
//...Packing for y face(4,8,9,16,18).................................
ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,&Aq[Component*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_Irecv(recvbuf_Y, 2*recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,MPI_COMM_SCALBL,&req2[2]);
//...Packing for Y face(3,7,10,15,17).................................
ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,&Aq[Component*N],N);
ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,&Aq[Component*7*N],N);
MPI_Isend(sendbuf_Y, sendCount_Y,MPI_DOUBLE,rank_Y,sendtag,MPI_COMM_SCALBL,&req1[3]);
MPI_Irecv(recvbuf_y, recvCount_y,MPI_DOUBLE,rank_y,recvtag,MPI_COMM_SCALBL,&req2[3]);
//...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q19_Pack(6,dvcSendList_z,0,sendCount_z,sendbuf_z,&Aq[Component*N],N);
ScaLBL_D3Q19_Pack(6,dvcSendList_z,0,sendCount_z,sendbuf_z,&Aq[Component*7*N],N);
MPI_Isend(sendbuf_z, sendCount_z,MPI_DOUBLE,rank_z,sendtag,MPI_COMM_SCALBL,&req1[4]);
MPI_Irecv(recvbuf_Z, recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,MPI_COMM_SCALBL,&req2[4]);
//...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q19_Pack(5,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,&Aq[Component*N],N);
ScaLBL_D3Q19_Pack(5,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,&Aq[Component*7*N],N);
//...................................................................................
// Send all the distributions
@@ -1483,33 +1483,33 @@ void ScaLBL_Communicator::RecvD3Q7AA(double *Aq, int Component){
// Unpack the distributions on the device
//...................................................................................
//...Unpacking for x face(2,8,10,12,14)................................
ScaLBL_D3Q7_Unpack(2,dvcRecvDist_x,0,recvCount_x,recvbuf_x,&Aq[Component*N],N);
ScaLBL_D3Q7_Unpack(2,dvcRecvDist_x,0,recvCount_x,recvbuf_x,&Aq[Component*7*N],N);
//...................................................................................
//...Packing for X face(1,7,9,11,13)................................
ScaLBL_D3Q7_Unpack(1,dvcRecvDist_X,0,recvCount_X,recvbuf_X,&Aq[Component*N],N);
ScaLBL_D3Q7_Unpack(1,dvcRecvDist_X,0,recvCount_X,recvbuf_X,&Aq[Component*7*N],N);
//...................................................................................
//...Packing for y face(4,8,9,16,18).................................
ScaLBL_D3Q7_Unpack(4,dvcRecvDist_y,0,recvCount_y,recvbuf_y,&Aq[Component*N],N);
ScaLBL_D3Q7_Unpack(4,dvcRecvDist_y,0,recvCount_y,recvbuf_y,&Aq[Component*7*N],N);
//...................................................................................
//...Packing for Y face(3,7,10,15,17).................................
ScaLBL_D3Q7_Unpack(3,dvcRecvDist_Y,0,recvCount_Y,recvbuf_Y,&Aq[Component*N],N);
ScaLBL_D3Q7_Unpack(3,dvcRecvDist_Y,0,recvCount_Y,recvbuf_Y,&Aq[Component*7*N],N);
//...................................................................................
if (BoundaryCondition > 0){
if (kproc != 0){
//...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,&Aq[Component*N],N);
ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,&Aq[Component*7*N],N);
}
if (kproc != nprocz-1){
//...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,&Aq[Component*N],N);
ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,&Aq[Component*7*N],N);
}
}
else {
//...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,&Aq[Component*N],N);
ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,&Aq[Component*7*N],N);
//...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,&Aq[Component*N],N);
ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,&Aq[Component*7*N],N);
}
//...................................................................................

26
cpu/Ion.cpp
View File

@@ -128,22 +128,22 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*Ci;
// q = 1
dist[nr2] = f1*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(ux+uEPx));
dist[nr2] = f1*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(ux+uEPx));
// q=2
dist[nr1] = f2*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(ux+uEPx));
dist[nr1] = f2*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(ux+uEPx));
// q = 3
dist[nr4] = f3*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(uy+uEPy));
dist[nr4] = f3*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(uy+uEPy));
// q = 4
dist[nr3] = f4*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(uy+uEPy));
dist[nr3] = f4*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(uy+uEPy));
// q = 5
dist[nr6] = f5*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(uz+uEPz));
dist[nr6] = f5*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(uz+uEPz));
// q = 6
dist[nr5] = f6*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(uz+uEPz));
dist[nr5] = f6*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(uz+uEPz));
}
}
@@ -183,22 +183,22 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Veloci
dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*Ci;
// q = 1
dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(ux+uEPx));
dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(ux+uEPx));
// q=2
dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(ux+uEPx));
dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(ux+uEPx));
// q = 3
dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(uy+uEPy));
dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(uy+uEPy));
// q = 4
dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(uy+uEPy));
dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(uy+uEPy));
// q = 5
dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(uz+uEPz));
dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(uz+uEPz));
// q = 6
dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(uz+uEPz));
dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(uz+uEPz));
}
@@ -231,7 +231,7 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity
Ci = Den[n+ion_component*Np];
CD = ChargeDensity[n];
CD_tmp = F*IonValence*Ci;
ChargeDensity[n] = CD*(IonValence>0) + CD_tmp;
ChargeDensity[n] = CD*(ion_component>0) + CD_tmp;
}
}

28
models/IonModel.cpp
View File

@@ -36,7 +36,7 @@ void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stoke
//---------------------- Default model parameters --------------------------//
T = 300.0;//temperature; unit [K]
Vt = kb*T/electron_charge;//thermal voltage; unit [Vy]
k2_inv = 4.5;//the inverse of 2nd-rank moment of D3Q7 lattice
k2_inv = 4.5;//speed of sound for D3Q7 lattice
h = 1.0;//resolution; unit: um/lu
tolerance = 1.0e-8;
number_ion_species = 1;
@@ -57,6 +57,8 @@ void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stoke
}
if (ion_db->keyExists( "temperature" )){
T = ion_db->getScalar<int>( "temperature" );
//re-calculate thermal voltage
Vt = kb*T/electron_charge;//thermal voltage; unit [Vy]
}
if (ion_db->keyExists( "number_ion_species" )){
number_ion_species = ion_db->getScalar<int>( "number_ion_species" );
@@ -104,6 +106,12 @@ void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stoke
}
}
//Read solid boundary condition specific to Ion model
BoundaryConditionSolid = 0;
if (ion_db->keyExists( "BC_Solid" )){
BoundaryConditionSolid = ion_db->getScalar<int>( "BC_Solid" );
}
// Read domain parameters
if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
h = domain_db->getScalar<double>( "voxel_length" );
@@ -112,10 +120,6 @@ void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stoke
if (domain_db->keyExists( "BC" )){
BoundaryCondition = domain_db->getScalar<int>( "BC" );
}
BoundaryConditionSolid = 0;
if (domain_db->keyExists( "BC_Solid" )){
BoundaryConditionSolid = domain_db->getScalar<int>( "BC_Solid" );
}
if (rank==0) printf("*****************************************************\n");
if (rank==0) printf("LB Ion Transport Solver: \n");
@@ -128,13 +132,13 @@ void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stoke
switch (BoundaryConditionSolid){
case 0:
if (rank==0) printf("LB Ion Solver: solid boundary: non-flux boundary is assigned");
if (rank==0) printf("LB Ion Solver: solid boundary: non-flux boundary is assigned\n");
break;
case 1:
if (rank==0) printf("LB Ion Solver: solid boundary: Neumann-type surfacen ion concentration is assigned");
if (rank==0) printf("LB Ion Solver: solid boundary: Neumann-type surfacen ion concentration is assigned\n");
break;
default:
if (rank==0) printf("LB Ion Solver: solid boundary: non-flux boundary is assigned");
if (rank==0) printf("LB Ion Solver: solid boundary: non-flux boundary is assigned\n");
break;
}
}
@@ -375,6 +379,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
// *************ODD TIMESTEP*************//
timestep++;
//Update ion concentration and charge density
if (rank==0) printf("timestep=%i; updating ion concentration and charge density\n",timestep);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_IonConcentration(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
@@ -386,7 +391,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
}
//LB-Ion collison
if (rank==0) printf("timestep=%i; execute collision step 1/2\n",timestep);
for (int ic=0; ic<number_ion_species; ic++){
if (rank==0) printf(" ic=%i; execute collsion step 1/2\n",ic);
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
@@ -394,7 +401,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
// Set boundary conditions
/* ... */
if (rank==0) printf("timestep=%i; execute collision step 2/2\n",timestep);
for (int ic=0; ic<number_ion_species; ic++){
if (rank==0) printf(" ic=%i; execute collsion step 2/2\n",ic);
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
}
@@ -409,6 +418,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
// *************EVEN TIMESTEP*************//
timestep++;
//Update ion concentration and charge density
if (rank==0) printf("timestep=%i; updating ion concentration and charge density\n",timestep);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FORM NORMAL
ScaLBL_D3Q7_AAeven_IonConcentration(&fq[ic*Np*7],&Ci[ic*Np],ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
@@ -419,6 +429,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
}
//LB-Ion collison
if (rank==0) printf("timestep=%i; execute collision step 1/2\n",timestep);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
@@ -427,6 +438,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
// Set boundary conditions
/* ... */
if (rank==0) printf("timestep=%i; execute collision step 2/2\n",timestep);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);

21
models/PoissonSolver.cpp
View File

@@ -22,7 +22,7 @@ void ScaLBL_Poisson::ReadParams(string filename){
domain_db = db->getDatabase( "Domain" );
electric_db = db->getDatabase( "Poisson" );
k2_inv = 4.5;//the inverse of 2nd-rank moment of D3Q7 lattice
k2_inv = 4.5;//speed of sound for D3Q7 lattice
gamma = 0.3;//time step of LB-Poisson equation
tau = 0.5+k2_inv*gamma;
timestepMax = 100000;
@@ -30,7 +30,7 @@ void ScaLBL_Poisson::ReadParams(string filename){
h = 1.0;//resolution; unit: um/lu
epsilon0 = 8.85e-12;//electrical permittivity of vaccum; unit:[C/(V*m)]
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilonR = 78.4;//default dielectric constant for water
epsilonR = 78.4;//default dielectric constant of water
epsilon_LB = epsilon0_LB*epsilonR;//electrical permittivity
analysis_interval = 1000;
@@ -50,6 +50,13 @@ void ScaLBL_Poisson::ReadParams(string filename){
if (electric_db->keyExists( "epsilonR" )){
epsilonR = electric_db->getScalar<double>( "epsilonR" );
}
// Read solid boundary condition specific to Poisson equation
BoundaryConditionSolid = 1;
if (electric_db->keyExists( "BC_Solid" )){
BoundaryConditionSolid = electric_db->getScalar<int>( "BC_Solid" );
}
// Read domain parameters
if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
h = domain_db->getScalar<double>( "voxel_length" );
@@ -59,10 +66,6 @@ void ScaLBL_Poisson::ReadParams(string filename){
if (domain_db->keyExists( "BC" )){
BoundaryCondition = domain_db->getScalar<int>( "BC" );
}
BoundaryConditionSolid = 1;
if (domain_db->keyExists( "BC_Solid" )){
BoundaryConditionSolid = domain_db->getScalar<int>( "BC_Solid" );
}
//Re-calcualte model parameters if user updates input
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
@@ -76,13 +79,13 @@ void ScaLBL_Poisson::ReadParams(string filename){
switch (BoundaryConditionSolid){
case 1:
if (rank==0) printf("LB-Poisson Solver: solid boundary: Dirichlet-type surfacen potential is assigned");
if (rank==0) printf("LB-Poisson Solver: solid boundary: Dirichlet-type surfacen potential is assigned\n");
break;
case 2:
if (rank==0) printf("LB-Poisson Solver: solid boundary: Neumann-type surfacen charge density is assigned");
if (rank==0) printf("LB-Poisson Solver: solid boundary: Neumann-type surfacen charge density is assigned\n");
break;
default:
if (rank==0) printf("LB-Poisson Solver: solid boundary: Dirichlet-type surfacen potential is assigned");
if (rank==0) printf("LB-Poisson Solver: solid boundary: Dirichlet-type surfacen potential is assigned\n");
break;
}
}

2
models/StokesModel.cpp
View File

@@ -86,7 +86,7 @@ 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]
time_conv = (h*h*1.0e-12)*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);

5
tests/lbpm_electrokinetic_dfh_simulator.cpp
View File

@@ -78,8 +78,13 @@ int main(int argc, char **argv)
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
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
if (rank==0) printf("timestep=%i; running Ion model\n",timestep);
IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential