diff --git a/models/MultiPhysController.cpp b/models/MultiPhysController.cpp index fcfb5403..9b361ad8 100644 --- a/models/MultiPhysController.cpp +++ b/models/MultiPhysController.cpp @@ -2,7 +2,7 @@ ScaLBL_Multiphys_Controller::ScaLBL_Multiphys_Controller(int RANK, int NP, MPI_Comm COMM): rank(RANK),nprocs(NP),Restart(0),timestepMax(0),num_iter_Stokes(0),num_iter_Ion(0), -analysis_interval(0),visualization_interval(0),tolerance(0),comm(COMM) +analysis_interval(0),visualization_interval(0),tolerance(0),time_conv_max(0),comm(COMM) { } @@ -25,6 +25,7 @@ void ScaLBL_Multiphys_Controller::ReadParams(string filename){ analysis_interval = 500; visualization_interval = 10000; tolerance = 1.0e-6; + time_conv_max = 0.0; // load input parameters if (study_db->keyExists( "timestepMax" )){ @@ -135,3 +136,12 @@ vector ScaLBL_Multiphys_Controller::getIonNumIter_PNP_coupling(double Stoke } return num_iter_ion; } + +void ScaLBL_Multiphys_Controller::getTimeConvMax_PNP_coupling(double StokesTimeConv,const vector &IonTimeConv){ + //Return maximum of the time converting factor from Stokes and ion solvers + vector TimeConv; + + TimeConv.assign(IonTimeConv.begin(),IonTimeConv.end()); + TimeConv.insert(TimeConv.begin(),StokesTimeConv); + time_conv_max = *max_element(TimeConv.begin(),TimeConv.end()); +} diff --git a/models/MultiPhysController.h b/models/MultiPhysController.h index f217248f..988f0225 100644 --- a/models/MultiPhysController.h +++ b/models/MultiPhysController.h @@ -27,6 +27,7 @@ public: int getStokesNumIter_PNP_coupling(double StokesTimeConv,const vector &IonTimeConv); vector getIonNumIter_PNP_coupling(double StokesTimeConv,const vector &IonTimeConv); //void getIonNumIter_PNP_coupling(double StokesTimeConv,vector &IonTimeConv,vector &IonTimeMax); + void getTimeConvMax_PNP_coupling(double StokesTimeConv,const vector &IonTimeConv); bool Restart; int timestepMax; @@ -35,6 +36,7 @@ public: int analysis_interval; int visualization_interval; double tolerance; + double time_conv_max; //double SchmidtNum;//Schmidt number = kinematic_viscosity/mass_diffusivity int rank,nprocs; diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index b0dde2c7..96d737bb 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -8,8 +8,11 @@ ScaLBL_Poisson::ScaLBL_Poisson(int RANK, int NP, MPI_Comm COMM): 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),WriteLog(0), -nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),BoundaryConditionSolid(0),Lx(0),Ly(0),Lz(0),comm(COMM) +chargeDen_dummy(0),WriteLog(0),nprocx(0),nprocy(0),nprocz(0), +BoundaryConditionInlet(0),BoundaryConditionOutlet(0),BoundaryConditionSolid(0),Lx(0),Ly(0),Lz(0), +Vin0(0),freqIn(0),t0_In(0),Vin_Type(0),Vout0(0),freqOut(0),t0_Out(0),Vout_Type(0), +TestPeriodic(0),TestPeriodicTime(0),TestPeriodicTimeConv(0),TestPeriodicSaveInterval(0), +comm(COMM) { } @@ -33,10 +36,12 @@ void ScaLBL_Poisson::ReadParams(string filename){ epsilonR = 78.4;//default dielectric constant of water epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity 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] WriteLog = false; + TestPeriodic = false; + TestPeriodicTime = 1.0;//unit: [sec] + TestPeriodicTimeConv = 0.01; //unit [sec/lt] + TestPeriodicSaveInterval = 0.1; //unit [sec] // LB-Poisson Model parameters if (electric_db->keyExists( "timestepMax" )){ @@ -57,6 +62,18 @@ void ScaLBL_Poisson::ReadParams(string filename){ if (electric_db->keyExists( "WriteLog" )){ WriteLog = electric_db->getScalar( "WriteLog" ); } + if (electric_db->keyExists( "TestPeriodic" )){ + TestPeriodic = electric_db->getScalar( "TestPeriodic" ); + } + if (electric_db->keyExists( "TestPeriodicTime" )){ + TestPeriodicTime = electric_db->getScalar( "TestPeriodicTime" ); + } + if (electric_db->keyExists( "TestPeriodicTimeConv" )){ + TestPeriodicTimeConv = electric_db->getScalar( "TestPeriodicTimeConv" ); + } + if (electric_db->keyExists( "TestPeriodicSaveInterval" )){ + TestPeriodicSaveInterval = electric_db->getScalar( "TestPeriodicSaveInterval" ); + } // Read solid boundary condition specific to Poisson equation BoundaryConditionSolid = 1; @@ -65,10 +82,15 @@ void ScaLBL_Poisson::ReadParams(string filename){ } // Read boundary condition for electric potential // BC = 0: normal periodic BC - // BC = 1: fixed inlet and outlet potential - BoundaryCondition = 0; - if (electric_db->keyExists( "BC" )){ - BoundaryCondition = electric_db->getScalar( "BC" ); + // BC = 1: fixed electric potential + // BC = 2: sine/cosine periodic electric potential (need extra input parameters) + BoundaryConditionInlet = 0; + BoundaryConditionOutlet = 0; + if (electric_db->keyExists( "BC_Inlet" )){ + BoundaryConditionInlet = electric_db->getScalar( "BC_Inlet" ); + } + if (electric_db->keyExists( "BC_Outlet" )){ + BoundaryConditionOutlet = electric_db->getScalar( "BC_Outlet" ); } // Read domain parameters @@ -342,15 +364,91 @@ void ScaLBL_Poisson::Create(){ void ScaLBL_Poisson::Potential_Init(double *psi_init){ - if (BoundaryCondition==1){ - if (electric_db->keyExists( "Vin" )){ - Vin = electric_db->getScalar( "Vin" ); - } - if (electric_db->keyExists( "Vout" )){ - Vout = electric_db->getScalar( "Vout" ); - } + //set up default boundary input parameters + Vin0 = Vout0 = 1.0; //unit: [V] + freqIn = freqOut = 50.0; //unit: [Hz] + t0_In = t0_Out = 0.0; //unit: [sec] + Vin_Type = Vout_Type = 1; //1->sin; 2->cos + Vin = 1.0; //Boundary-z (inlet) electric potential + Vout = 1.0; //Boundary-Z (outlet) electric potential + + if (BoundaryConditionInlet>0){ + switch (BoundaryConditionInlet){ + case 1: + if (electric_db->keyExists( "Vin" )){ + Vin = electric_db->getScalar( "Vin" ); + } + if (rank==0) printf("LB-Poisson Solver: inlet boundary; fixed electric potential Vin = %.3g \n",Vin); + break; + case 2: + if (electric_db->keyExists( "Vin0" )){//voltage amplitude; unit: Volt + Vin0 = electric_db->getScalar( "Vin0" ); + } + if (electric_db->keyExists( "freqIn" )){//unit: Hz + freqIn = electric_db->getScalar( "freqIn" ); + } + if (electric_db->keyExists( "t0_In" )){//timestep shift, unit: lt + t0_In = electric_db->getScalar( "t0_In" ); + } + if (electric_db->keyExists( "Vin_Type" )){ + //type=1 -> sine + //tyep=2 -> cosine + Vin_Type = electric_db->getScalar( "Vin_Type" ); + if (Vin_Type>2 || Vin_Type<=0) ERROR("Error: user-input Vin_Type is currently not supported! \n"); + } + if (rank==0){ + if (Vin_Type==1){ + printf("LB-Poisson Solver: inlet boundary; periodic electric potential Vin = %.3g*Sin[2*pi*%.3g*(t+%.3g)] \n",Vin,freqIn,t0_In); + printf(" V0 = %.3g [V], frequency = %.3g [Hz], timestep shift = %.3g [sec] \n",Vin,freqIn,t0_In); + } + else if (Vin_Type==2){ + printf("LB-Poisson Solver: inlet boundary; periodic electric potential Vin = %.3g*Cos[2*pi*%.3g*(t+%.3g)] \n",Vin,freqIn,t0_In); + printf(" V0 = %.3g [V], frequency = %.3g [Hz], timestep shift = %.3g [sec] \n",Vin,freqIn,t0_In); + } + } + break; + } + } + if (BoundaryConditionOutlet>0){ + switch (BoundaryConditionOutlet){ + case 1: + if (electric_db->keyExists( "Vout" )){ + Vout = electric_db->getScalar( "Vout" ); + } + if (rank==0) printf("LB-Poisson Solver: outlet boundary; fixed electric potential Vin = %.3g \n",Vout); + break; + case 2: + if (electric_db->keyExists( "Vout0" )){//voltage amplitude; unit: Volt + Vout0 = electric_db->getScalar( "Vout0" ); + } + if (electric_db->keyExists( "freqOut" )){//unit: Hz + freqOut = electric_db->getScalar( "freqOut" ); + } + if (electric_db->keyExists( "t0_Out" )){//timestep shift, unit: lt + t0_Out = electric_db->getScalar( "t0_Out" ); + } + if (electric_db->keyExists( "Vout_Type" )){ + //type=1 -> sine + //tyep=2 -> cosine + Vout_Type = electric_db->getScalar( "Vout_Type" ); + if (Vout_Type>2 || Vin_Type<=0) ERROR("Error: user-input Vout_Type is currently not supported! \n"); + } + if (rank==0){ + if (Vout_Type==1){ + printf("LB-Poisson Solver: outlet boundary; periodic electric potential Vout = %.3g*Sin[2*pi*%.3g*(t+%.3g)] \n",Vout,freqOut,t0_Out); + printf(" V0 = %.3g [V], frequency = %.3g [Hz], timestep shift = %.3g [sec] \n",Vout,freqOut,t0_Out); + } + else if (Vout_Type==2){ + printf("LB-Poisson Solver: outlet boundary; periodic electric potential Vout = %.3g*Cos[2*pi*%.3g*(t+%.3g)] \n",Vout,freqOut,t0_Out); + printf(" V0 = %.3g [V], frequency = %.3g [Hz], timestep shift = %.3g [sec] \n",Vout,freqOut,t0_Out); + } + } + break; + } } //By default only periodic BC is applied and Vin=Vout=1.0, i.e. there is no potential gradient along Z-axis + if (BoundaryConditionInlet==2) Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,t0_In,Vin_Type,0); + if (BoundaryConditionOutlet==2) Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,t0_Out,Vout_Type,0); double slope = (Vout-Vin)/(Nz-2); double psi_linearized; 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); + if (BoundaryConditionInlet > 0){ + switch (BoundaryConditionInlet){ + case 1: + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + case 2: + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,t0_In,Vin_Type,timestep_from_Study); + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + } + } + if (BoundaryConditionOutlet > 0){ + switch (BoundaryConditionOutlet){ + case 1: + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + case 2: + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,t0_Out,Vout_Type,timestep_from_Study); + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + } } //-------------------------// ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(NeighborList, dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np); } -void ScaLBL_Poisson::SolveElectricPotentialAAeven(){ +void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){ 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); + if (BoundaryConditionInlet > 0){ + switch (BoundaryConditionInlet){ + case 1: + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + case 2: + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,t0_In,Vin_Type,timestep_from_Study); + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + } + } + if (BoundaryConditionOutlet > 0){ + switch (BoundaryConditionOutlet){ + case 1: + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + case 2: + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,t0_Out,Vout_Type,timestep_from_Study); + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + } } //-------------------------// ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np); diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index 74abd775..ebcac179 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -9,6 +9,7 @@ #include #include #include +#include #include "common/ScaLBL.h" #include "common/Communication.h" @@ -16,6 +17,7 @@ #include "analysis/Minkowski.h" #include "ProfilerApp.h" +#define _USE_MATH_DEFINES #ifndef ScaLBL_POISSON_INC #define ScaLBL_POISSON_INC @@ -41,7 +43,8 @@ public: //bool Restart,pBC; int timestep,timestepMax; int analysis_interval; - int BoundaryCondition; + int BoundaryConditionInlet; + int BoundaryConditionOutlet; int BoundaryConditionSolid; double tau; double tolerance; @@ -50,11 +53,18 @@ public: double Vin, Vout; double chargeDen_dummy;//for debugging bool WriteLog; + double Vin0,freqIn,t0_In,Vin_Type; + double Vout0,freqOut,t0_Out,Vout_Type; + bool TestPeriodic; + double TestPeriodicTime;//unit: [sec] + double TestPeriodicTimeConv; //unit [sec/lt] + double TestPeriodicSaveInterval; //unit [sec] int Nx,Ny,Nz,N,Np; int rank,nprocx,nprocy,nprocz,nprocs; double Lx,Ly,Lz; double h;//image resolution + double time_conv;//phys to LB time converting factor; unit=[sec/lt] std::shared_ptr Dm; // this domain is for analysis std::shared_ptr Mask; // this domain is for lbm @@ -97,6 +107,7 @@ private: void SolvePoissonAAodd(double *ChargeDensity); void SolvePoissonAAeven(double *ChargeDensity); void getConvergenceLog(int timestep,double error); + double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int V_type,int time_step); }; #endif diff --git a/tests/TestPoissonSolver.cpp b/tests/TestPoissonSolver.cpp index 32353f65..5683ace1 100644 --- a/tests/TestPoissonSolver.cpp +++ b/tests/TestPoissonSolver.cpp @@ -53,14 +53,36 @@ int main(int argc, char **argv) PoissonSolver.SetDomain(); PoissonSolver.ReadInput(); PoissonSolver.Create(); - PoissonSolver.Initialize(); + if (PoissonSolver.TestPeriodic==true){ + PoissonSolver.Initialize(PoissonSolver.TestPeriodicTimeConv); + } + else { + PoissonSolver.Initialize(0); + } //Initialize dummy charge density for test PoissonSolver.DummyChargeDensity(); - PoissonSolver.Run(PoissonSolver.ChargeDensityDummy); - PoissonSolver.getElectricPotential_debug(1); - PoissonSolver.getElectricField_debug(1); + if (PoissonSolver.TestPeriodic==true){ + if (rank==0) printf("Testing periodic voltage input is enabled. Total test time is %.3g[s], saving data every %.3g[s]; + user-specified time resolution is %.3g[s/lt]\n", + PoissonSolver.TestPeriodicTime,PoissonSolver.TestPeriodicSaveInterval,PoissonSolver.TestPeriodicTimeConv); + int timestep = 0; + while (timestep<(PoissonSolver.TestPeriodicTime/PoissonSolver.TestPeriodicTimeConv)){ + timestep++; + PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,timestep); + if (timestep%(PoissonSolver.TestPeriodicSaveInterval/PoissonSolver.TestPeriodicTimeConv)==0){ + if (rank==0) printf(" Time = %.3g[s]; saving electric potential and field\n",timestep*PoissonSolver.TestPeriodicTimeConv); + PoissonSolver.getElectricPotential_debug(timestep*PoissonSolver.TestPeriodicTimeConv); + PoissonSolver.getElectricField_debug(timestep*PoissonSolver.TestPeriodicTimeConv); + } + } + } + else { + PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,1); + PoissonSolver.getElectricPotential_debug(1); + PoissonSolver.getElectricField_debug(1); + } if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n"); if (rank==0) printf("*************************************************************\n"); diff --git a/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp b/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp index 2b3726a4..93493331 100644 --- a/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp +++ b/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp @@ -80,20 +80,22 @@ int main(int argc, char **argv) IonModel.timestepMax = Study.getIonNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv); IonModel.Initialize(); + // Get maximal time converting factor based on Sotkes and Ion solvers + Study.getTimeConvMax_PNP_coupling(StokesModel.time_conv,IonModel.time_conv); // Initialize LB-Poisson model PoissonSolver.ReadParams(filename); PoissonSolver.SetDomain(); PoissonSolver.ReadInput(); PoissonSolver.Create(); - PoissonSolver.Initialize(); + PoissonSolver.Initialize(Study.time_conv_max); int timestep=0; while (timestep < Study.timestepMax){ timestep++; - PoissonSolver.Run(IonModel.ChargeDensity);//solve Poisson equtaion to get steady-state electrical potental + PoissonSolver.Run(IonModel.ChargeDensity,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 IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential