From 5756d6f1386291e181a166bc617d1e7e63a2398d Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Tue, 18 Aug 2020 12:40:41 -0400 Subject: [PATCH] fix a few trivial bugs; add some checkpoint print; still debugging --- common/ScaLBL.cpp | 28 ++++++++++----------- cpu/Ion.cpp | 26 +++++++++---------- models/IonModel.cpp | 28 +++++++++++++++------ models/PoissonSolver.cpp | 21 +++++++++------- models/StokesModel.cpp | 2 +- tests/lbpm_electrokinetic_dfh_simulator.cpp | 5 ++++ 6 files changed, 65 insertions(+), 45 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index a77afbca..854c8b49 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -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); } //................................................................................... diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index b04b3c9c..2c7e72a9 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -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; } } diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 14e0a64c..fca2871c 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -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( "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( "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( "BC_Solid" ); + } + // Read domain parameters if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu h = domain_db->getScalar( "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( "BC" ); } - BoundaryConditionSolid = 0; - if (domain_db->keyExists( "BC_Solid" )){ - BoundaryConditionSolid = domain_db->getScalar( "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; icSendD3Q7AA(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; icFirstInterior(), 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; icLastExterior(), 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; icSendD3Q7AA(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; icFirstInterior(), 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; icLastExterior(), Np); diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 27c0cda5..3cf8a6d2 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -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( "epsilonR" ); } + + // Read solid boundary condition specific to Poisson equation + BoundaryConditionSolid = 1; + if (electric_db->keyExists( "BC_Solid" )){ + BoundaryConditionSolid = electric_db->getScalar( "BC_Solid" ); + } + // Read domain parameters if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu h = domain_db->getScalar( "voxel_length" ); @@ -59,10 +66,6 @@ void ScaLBL_Poisson::ReadParams(string filename){ if (domain_db->keyExists( "BC" )){ BoundaryCondition = domain_db->getScalar( "BC" ); } - BoundaryConditionSolid = 1; - if (domain_db->keyExists( "BC_Solid" )){ - BoundaryConditionSolid = domain_db->getScalar( "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; } } diff --git a/models/StokesModel.cpp b/models/StokesModel.cpp index 46f28f45..1b733745 100644 --- a/models/StokesModel.cpp +++ b/models/StokesModel.cpp @@ -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); diff --git a/tests/lbpm_electrokinetic_dfh_simulator.cpp b/tests/lbpm_electrokinetic_dfh_simulator.cpp index bdc15ef4..7cf835da 100644 --- a/tests/lbpm_electrokinetic_dfh_simulator.cpp +++ b/tests/lbpm_electrokinetic_dfh_simulator.cpp @@ -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