diff --git a/common/Domain.cpp b/common/Domain.cpp index eadc4592..d75ec4a5 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -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 external boundary conditions are applied remove solid - if (BoundaryCondition > 0 && BoundaryCondition !=5 && kproc() == 0){ - if (inlet_layers_z < 4){ - inlet_layers_z=4; - 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); - } - } - for (int k=0; k 0 && BoundaryCondition !=5 && kproc() == nprocz-1){ - if (outlet_layers_z < 4){ - outlet_layers_z=4; - 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); - } - } - for (int k=Nz-outlet_layers_z; k 0 && BoundaryCondition !=5 && kproc() == 0){ +// if (inlet_layers_z < 4){ +// inlet_layers_z=4; +// 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); +// } +// } +// for (int k=0; k 0 && BoundaryCondition !=5 && kproc() == nprocz-1){ +// if (outlet_layers_z < 4){ +// outlet_layers_z=4; +// 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); +// } +// } +// for (int k=Nz-outlet_layers_z; k0)+Psi[ijk]*(id<=0);// get neighbor for phi - 18 //............Compute the Color Gradient................................... - nx = -1.f/18.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)); - nz = -1.f/18.f*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); + //nx = 1.f/6.f*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); + //ny = 1.f/6.f*(m3-m4+0.5*(m7-m8-m9+m10+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[Np+n] = ny; diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index d400df05..89a1c42a 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -6,8 +6,9 @@ #include "common/ReadMicroCT.h" 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), +chargeDen_dummy(0), 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" ); 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 - gamma = 1.0;//time step of LB-Poisson equation - tau = 0.5+k2_inv*gamma; + k2_inv = 4.0;//speed of sound for D3Q7 lattice + tau = 0.5+k2_inv; timestepMax = 100000; tolerance = 1.0e-6;//stopping criterion for obtaining steady-state electricla potential h = 1.0;//resolution; unit: um/lu @@ -36,6 +35,7 @@ void ScaLBL_Poisson::ReadParams(string filename){ analysis_interval = 1000; Vin = 1.0; //Boundary-z (inlet) 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 if (electric_db->keyExists( "timestepMax" )){ @@ -47,12 +47,12 @@ void ScaLBL_Poisson::ReadParams(string filename){ if (electric_db->keyExists( "tolerance" )){ tolerance = electric_db->getScalar( "tolerance" ); } - if (electric_db->keyExists( "gamma" )){ - gamma = electric_db->getScalar( "gamma" ); - } if (electric_db->keyExists( "epsilonR" )){ epsilonR = electric_db->getScalar( "epsilonR" ); } + if (electric_db->keyExists( "DummyChargeDen" )){ + chargeDen_dummy = electric_db->getScalar( "DummyChargeDen" ); + } // Read solid boundary condition specific to Poisson equation BoundaryConditionSolid = 1; @@ -76,7 +76,6 @@ void ScaLBL_Poisson::ReadParams(string filename){ //Re-calcualte model parameters if user updates input epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)] epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity - tau = 0.5+k2_inv*gamma; if (rank==0) printf("***********************************************************************************\n"); 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 if (BoundaryConditionSolid==2){ //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; } label_count[idx] += 1.0; @@ -284,11 +282,10 @@ void ScaLBL_Poisson::Create(){ //........................................................................... ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); 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 **) &Psi, sizeof(double)*Nx*Ny*Nz); ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np); - //ScaLBL_AllocateDeviceMemory((void **) &PoissonSolid, sizeof(double)*Nx*Ny*Nz); //........................................................................... // Update GPU data structures @@ -329,26 +326,13 @@ void ScaLBL_Poisson::Create(){ MPI_Barrier(comm); delete [] neighborList; // copy node ID - ScaLBL_CopyToDevice(dvcID, Mask->id, sizeof(signed char)*Nx*Ny*Nz); - ScaLBL_DeviceBarrier(); + //ScaLBL_CopyToDevice(dvcID, Mask->id, sizeof(signed char)*Nx*Ny*Nz); + //ScaLBL_DeviceBarrier(); //Initialize solid boundary for electric potential ScaLBL_Comm->SetupBounceBackList(Map, Mask->id, Np); 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){ @@ -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, 0, ScaLBL_Comm->LastExterior(), Np); 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){ @@ -418,20 +412,17 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){ //************************************************************************/ // *************ODD TIMESTEP*************// timestep++; - SolveElectricPotentialAAodd(); - //compute electric field - SolveElectricField(); - //perform collision - SolvePoissonAAodd(ChargeDensity); + + SolveElectricPotentialAAodd();//update electric potential + //SolveElectricField(); //deprecated - compute electric field + SolvePoissonAAodd(ChargeDensity);//perform collision ScaLBL_DeviceBarrier(); MPI_Barrier(comm); // *************EVEN TIMESTEP*************// timestep++; - SolveElectricPotentialAAeven(); - //compute electric field - SolveElectricField(); - //perform collision - SolvePoissonAAeven(ChargeDensity); + SolveElectricPotentialAAeven();//update electric potential + //SolveElectricField();//deprecated - compute electric field + SolvePoissonAAeven(ChargeDensity);//perform collision 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; kSendD3Q7AA(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){ for (int k=0;kSendHalo(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(){ // diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index 42c47afb..921963ed 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -28,13 +28,9 @@ public: void Create(); void Initialize(); void Run(double *ChargeDensity); - void SolveElectricPotentialAAodd(); - void SolveElectricPotentialAAeven(); - void SolveElectricField(); - void SolvePoissonAAodd(double *ChargeDensity); - void SolvePoissonAAeven(double *ChargeDensity); void getElectricPotential(int timestep); void getElectricField(int timestep); + void DummyChargeDensity();//for debugging //bool Restart,pBC; int timestep,timestepMax; @@ -43,9 +39,10 @@ public: int BoundaryConditionSolid; double tau; double tolerance; - double k2_inv,gamma; + double k2_inv; double epsilon0,epsilon0_LB,epsilonR,epsilon_LB; double Vin, Vout; + double chargeDen_dummy;//for debugging int Nx,Ny,Nz,N,Np; int rank,nprocx,nprocy,nprocz,nprocs; @@ -66,11 +63,11 @@ public: DoubleArray Psi_host; int *NeighborList; int *dvcMap; - signed char *dvcID; + //signed char *dvcID; double *fq; double *Psi; double *ElectricField; - //double *PoissonSolid; + double *ChargeDensityDummy;// for debugging private: MPI_Comm comm; @@ -85,5 +82,10 @@ private: void AssignSolidBoundary(double *poisson_solid); void Potential_Init(double *psi_init); void ElectricField_LB_to_Phys(DoubleArray &Efield_reg); + void SolveElectricPotentialAAodd(); + void SolveElectricPotentialAAeven(); + //void SolveElectricField(); + void SolvePoissonAAodd(double *ChargeDensity); + void SolvePoissonAAeven(double *ChargeDensity); }; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1a8bfac0..43167b3f 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -47,6 +47,7 @@ ADD_LBPM_TEST( TestTorusEvolve ) ADD_LBPM_TEST( TestTopo3D ) ADD_LBPM_TEST( TestFluxBC ) ADD_LBPM_TEST( TestMap ) +ADD_LBPM_TEST( TestPoissonSolver ) #ADD_LBPM_TEST( TestMRT ) #ADD_LBPM_TEST( TestColorGrad ) #ADD_LBPM_TEST( TestColorGradDFH ) diff --git a/tests/TestPoissonSolver.cpp b/tests/TestPoissonSolver.cpp new file mode 100644 index 00000000..309a03c7 --- /dev/null +++ b/tests/TestPoissonSolver.cpp @@ -0,0 +1,102 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#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