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

fix miscellaneous bugs and update the data structure of electric potential

Browse Source
This commit is contained in:
Rex Zhe Li
2020-08-28 11:15:55 -04:00
parent 59ffd7bfd6
commit aa26fcafda
12 changed files with 1027 additions and 206 deletions
Download Patch File Download Diff File Expand all files Collapse all files

33
common/ScaLBL.cpp
View File

@@ -2053,3 +2053,36 @@ void ScaLBL_Communicator::PrintD3Q19(){
delete [] TempBuffer;
}
void ScaLBL_Communicator::D3Q7_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time){
if (kproc == 0) {
if (time%2==0){
ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(dvcSendList_z, fq, Vin, sendCount_z, N);
}
else{
ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_z(neighborList, dvcSendList_z, fq, Vin, sendCount_z, N);
}
}
}
void ScaLBL_Communicator::D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time){
if (kproc == nprocz-1){
if (time%2==0){
ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_Z(dvcSendList_Z, fq, Vout, sendCount_Z, N);
}
else{
ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(neighborList, dvcSendList_Z, fq, Vout, sendCount_Z, N);
}
}
}
void ScaLBL_Communicator::Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin){
if (kproc == 0) {
ScaLBL_Poisson_D3Q7_BC_z(dvcSendList_z, Map, Psi, Vin, sendCount_z);
}
}
void ScaLBL_Communicator::Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout){
if (kproc == nprocz-1){
ScaLBL_Poisson_D3Q7_BC_Z(dvcSendList_Z, Map, Psi, Vout, sendCount_Z);
}
}

40
common/ScaLBL.h
View File

@@ -46,11 +46,8 @@ extern "C" void ScaLBL_UnpackDenD3Q7(int *list, int count, double *recvbuf, int
extern "C" void ScaLBL_D3Q19_Init(double *Dist, int Np);
extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np);
extern "C" void ScaLBL_D3Q19_Momentum_Phys(double *dist, double *vel, double h, double time_conv, int Np);
extern "C" void ScaLBL_D3Q19_Pressure(double *dist, double *press, int Np);
// BGK MODEL
@@ -95,21 +92,30 @@ extern "C" void ScaLBL_IonConcentration_Phys(double *Den, double h, int ion_comp
// LBM Poisson solver
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, double *Den_charge, double *Psi, double *ElectricField, 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,double gamma,
int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, double *Psi, double *ElectricField, 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,double gamma,
int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Poisson_Init(double *dist, 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_AAeven_Poisson_ElectricPotential(int *Map, double *dist, double *Psi, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Poisson_getElectricField(double *dist, double *ElectricField, double tau, int Np);
extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
int strideY, int strideZ,int start, int finish, int Np);
// LBM Stokes Model (adapted from MRT model)
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz, double Ex, double Ey, double Ez, int start, int finish, int Np);
double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz, double Ex, double Ey, double Ez, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np);
// MRT MODEL
extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx,
@@ -190,6 +196,18 @@ extern "C" void ScaLBL_Solid_Dirichlet_D3Q7(double *dist,double *BoundaryValue,i
extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int *BounceBackDist_list,int *BounceBackSolid_list,int N);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np);
extern "C" void ScaLBL_Poisson_D3Q7_BC_z(int *list, int *Map, double *Psi, double Vin, int count);
extern "C" void ScaLBL_Poisson_D3Q7_BC_Z(int *list, int *Map, double *Psi, double Vout, int count);
class ScaLBL_Communicator{
public:
//......................................................................................
@@ -249,6 +267,10 @@ public:
void D3Q19_Reflection_BC_z(double *fq);
void D3Q19_Reflection_BC_Z(double *fq);
double D3Q19_Flux_BC_z(int *neighborList, double *fq, double flux, int time);
void D3Q7_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time);
void D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time);
void Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin);
void Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout);
// Debugging and unit testing functions
void PrintD3Q19();

107
cpu/D3Q7BC.cpp
View File

@@ -30,3 +30,110 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){
for (int idx=0; idx<count; idx++){
int n = list[idx];
double f0 = dist[n];
double f1 = dist[2*Np+n];
double f2 = dist[1*Np+n];
double f3 = dist[4*Np+n];
double f4 = dist[3*Np+n];
double f6 = dist[5*Np+n];
//...................................................
double f5 = Vin - (f0+f1+f2+f3+f4+f6);
dist[6*Np+n] = f5;
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np){
for (int idx=0; idx<count; idx++){
int n = list[idx];
double f0 = dist[n];
double f1 = dist[2*Np+n];
double f2 = dist[1*Np+n];
double f3 = dist[4*Np+n];
double f4 = dist[3*Np+n];
double f5 = dist[6*Np+n];
//...................................................
double f6 = Vout - (f0+f1+f2+f3+f4+f5);
dist[5*Np+n] = f6;
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np){
int nread,nr5;
for (int idx=0; idx<count; idx++){
int n = list[idx];
double f0 = dist[n];
nread = d_neighborList[n];
double f1 = dist[nread];
nread = d_neighborList[n+2*Np];
double f3 = dist[nread];
nread = d_neighborList[n+Np];
double f2 = dist[nread];
nread = d_neighborList[n+3*Np];
double f4 = dist[nread];
nread = d_neighborList[n+5*Np];
double f6 = dist[nread];
// Unknown distributions
nr5 = d_neighborList[n+4*Np];
double f5 = Vin - (f0+f1+f2+f3+f4+f6);
dist[nr5] = f5;
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np){
int nread,nr6;
for (int idx=0; idx<count; idx++){
int n = list[idx];
double f0 = dist[n];
nread = d_neighborList[n];
double f1 = dist[nread];
nread = d_neighborList[n+2*Np];
double f3 = dist[nread];
nread = d_neighborList[n+4*Np];
double f5 = dist[nread];
nread = d_neighborList[n+Np];
double f2 = dist[nread];
nread = d_neighborList[n+3*Np];
double f4 = dist[nread];
// unknown distributions
nr6 = d_neighborList[n+5*Np];
double f6 = Vout - (f0+f1+f2+f3+f4+f5);
dist[nr6] = f6;
}
}
extern "C" void ScaLBL_Poisson_D3Q7_BC_z(int *list, int *Map, double *Psi, double Vin, int count)
{
int idx,n,nm;
for (idx=0; idx<count; idx++){
n = list[idx];
nm = Map[n];
Psi[nm] = Vin;
}
}
extern "C" void ScaLBL_Poisson_D3Q7_BC_Z(int *list, int *Map, double *Psi, double Vout, int count)
{
int idx,n,nm;
for (idx=0; idx<count; idx++){
n = list[idx];
nm = Map[n];
Psi[nm] = Vout;
}
}

350
cpu/Poisson.cpp
View File

@@ -1,19 +1,110 @@
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,
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){
int n;
double psi;//electric potential
double Ex,Ey,Ez;//electrical field
double fq;
int nread;
int idx;
for (n=start; n<finish; n++){
// q=0
fq = dist[n];
psi = fq;
// q=1
nread = neighborList[n];
fq = dist[nread];
psi += fq;
// q=2
nread = neighborList[n+Np];
fq = dist[nread];
psi += fq;
// q=3
nread = neighborList[n+2*Np];
fq = dist[nread];
psi += fq;
// q = 4
nread = neighborList[n+3*Np];
fq = dist[nread];
psi += fq;
// q=5
nread = neighborList[n+4*Np];
fq = dist[nread];
psi += fq;
// q = 6
nread = neighborList[n+5*Np];
fq = dist[nread];
psi += fq;
idx=Map[n];
Psi[idx] = psi;
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *dist, double *Psi, int start, int finish, int Np){
int n;
double psi;//electric potential
double fq;
int idx;
for (n=start; n<finish; n++){
// q=0
fq = dist[n];
psi = fq;
// q=1
fq = dist[2*Np+n];
psi += fq;
// q=2
fq = dist[1*Np+n];
psi += fq;
// q=3
fq = dist[4*Np+n];
psi += fq;
// q=4
fq = dist[3*Np+n];
psi += fq;
// q=5
fq = dist[6*Np+n];
psi += fq;
// q=6
fq = dist[5*Np+n];
psi += fq;
idx=Map[n];
Psi[idx] = psi;
}
}
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,double gamma,int start, int finish, int Np){
int n;
double psi;//electric potential
double Ex,Ey,Ez;//electric field
double rho_e;//local charge density
double f0,f1,f2,f3,f4,f5,f6;
int nr1,nr2,nr3,nr4,nr5,nr6;
double rlx=1.0/tau;
int idx;
for (n=start; n<finish; n++){
//Load data
rho_e = Den_charge[n];
rho_e = gamma*rho_e/epsilon_LB;
idx=Map[n];
psi = Psi[idx];
// q=0
f0 = dist[n];
@@ -40,53 +131,63 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, doubl
nr6 = neighborList[n+5*Np];
f6 = dist[nr6];
psi = f0+f2+f1+f4+f3+f6+f5;
Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
Ey = (f3-f4)*rlx*4.5;
Ez = (f5-f6)*rlx*4.5;
ElectricField[n+0*Np] = Ex;
ElectricField[n+1*Np] = Ey;
ElectricField[n+2*Np] = Ez;
Psi[n] = psi;
//Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
//Ey = (f3-f4)*rlx*4.5;
//Ez = (f5-f6)*rlx*4.5;
//Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
//Ey = (f3-f4)*rlx*4.0;
//Ez = (f5-f6)*rlx*4.0;
//ElectricField[n+0*Np] = Ex;
//ElectricField[n+1*Np] = Ey;
//ElectricField[n+2*Np] = Ez;
// q = 0
dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e);
//dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e);
dist[n] = f0*(1.0-rlx) + 0.25*(rlx*psi+rho_e);
// q = 1
dist[nr2] = f1*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
//dist[nr2] = f1*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
dist[nr2] = f1*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 2
dist[nr1] = f2*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
//dist[nr1] = f2*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
dist[nr1] = f2*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 3
dist[nr4] = f3*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
//dist[nr4] = f3*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
dist[nr4] = f3*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 4
dist[nr3] = f4*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
//dist[nr3] = f4*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
dist[nr3] = f4*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 5
dist[nr6] = f5*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
//dist[nr6] = f5*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
dist[nr6] = f5*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// q = 6
dist[nr5] = f6*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
//dist[nr5] = f6*(1.0-rlx) + 0.1111111111111111*(rlx*psi+rho_e);
dist[nr5] = f6*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//........................................................................
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, double *Psi, double *ElectricField, 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,double gamma,int start, int finish, int Np){
int n;
double psi;//electric potential
double Ex,Ey,Ez;//electrical field
double Ex,Ey,Ez;//electric field
double rho_e;//local charge density
double f0,f1,f2,f3,f4,f5,f6;
double rlx=1.0/tau;
int idx;
for (n=start; n<finish; n++){
//Load data
rho_e = Den_charge[n];
rho_e = gamma*rho_e/epsilon_LB;
idx=Map[n];
psi = Psi[idx];
f0 = dist[n];
f1 = dist[2*Np+n];
@@ -97,49 +198,208 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, dou
f6 = dist[5*Np+n];
psi = f0+f2+f1+f4+f3+f6+f5;
Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
Ey = (f3-f4)*rlx*4.5;
Ez = (f5-f6)*rlx*4.5;
ElectricField[n+0*Np] = Ex;
ElectricField[n+1*Np] = Ey;
ElectricField[n+2*Np] = Ez;
Psi[n] = psi;
//Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
//Ey = (f3-f4)*rlx*4.5;
//Ez = (f5-f6)*rlx*4.5;
//Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
//Ey = (f3-f4)*rlx*4.0;
//Ez = (f5-f6)*rlx*4.0;
//ElectricField[n+0*Np] = Ex;
//ElectricField[n+1*Np] = Ey;
//ElectricField[n+2*Np] = Ez;
// q = 0
dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e);
//dist[n] = f0*(1.0-rlx) + 0.3333333333333333*(rlx*psi+rho_e);
dist[n] = f0*(1.0-rlx) + 0.25*(rlx*psi+rho_e);
// 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.1111111111111111*(rlx*psi+rho_e);
dist[1*Np+n] = f1*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// 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.1111111111111111*(rlx*psi+rho_e);
dist[2*Np+n] = f2*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// 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.1111111111111111*(rlx*psi+rho_e);
dist[3*Np+n] = f3*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// 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.1111111111111111*(rlx*psi+rho_e);
dist[4*Np+n] = f4*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// 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.1111111111111111*(rlx*psi+rho_e);
dist[5*Np+n] = f5*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
// 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.1111111111111111*(rlx*psi+rho_e);
dist[6*Np+n] = f6*(1.0-rlx) + 0.125*(rlx*psi+rho_e);
//........................................................................
}
}
extern "C" void ScaLBL_D3Q7_Poisson_Init(double *dist, int Np)
extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np)
{
int n;
for (n=0; n<Np; n++){
dist[0*Np+n] = 0.3333333333333333;
dist[1*Np+n] = 0.1111111111111111;
dist[2*Np+n] = 0.1111111111111111;
dist[3*Np+n] = 0.1111111111111111;
dist[4*Np+n] = 0.1111111111111111;
dist[5*Np+n] = 0.1111111111111111;
dist[6*Np+n] = 0.1111111111111111;
int ijk;
for (n=start; n<finish; n++){
ijk = Map[n];
//dist[0*Np+n] = 0.3333333333333333*Psi[n];
//dist[1*Np+n] = 0.1111111111111111*Psi[n];
//dist[2*Np+n] = 0.1111111111111111*Psi[n];
//dist[3*Np+n] = 0.1111111111111111*Psi[n];
//dist[4*Np+n] = 0.1111111111111111*Psi[n];
//dist[5*Np+n] = 0.1111111111111111*Psi[n];
//dist[6*Np+n] = 0.1111111111111111*Psi[n];
//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.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];
}
}
extern "C" void ScaLBL_D3Q7_Poisson_getElectricField(double *dist, double *ElectricField, double tau, int Np){
int n;
// distributions
double f1,f2,f3,f4,f5,f6;
double Ex,Ey,Ez;
double rlx=1.0/tau;
for (n=0; n<Np; n++){
//........................................................................
// Registers to store the distributions
//........................................................................
f1 = dist[Np+n];
f2 = dist[2*Np+n];
f3 = dist[3*Np+n];
f4 = dist[4*Np+n];
f5 = dist[5*Np+n];
f6 = dist[6*Np+n];
//.................Compute the Electric Field...................................
//Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
//Ey = (f3-f4)*rlx*4.5;
//Ez = (f5-f6)*rlx*4.5;
Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
Ey = (f3-f4)*rlx*4.0;
Ez = (f5-f6)*rlx*4.0;
//..................Write the Electric Field.....................................
ElectricField[0*Np+n] = Ex;
ElectricField[1*Np+n] = Ey;
ElectricField[2*Np+n] = Ez;
//........................................................................
}
}
extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
int strideY, int strideZ,int start, int finish, int Np){
int n,nn;
int ijk;
int id;
// distributions
double m1,m2,m3,m4,m5,m6,m7,m8,m9;
double m10,m11,m12,m13,m14,m15,m16,m17,m18;
double nx,ny,nz;
for (n=start; n<finish; n++){
// Get the 1D index based on regular data layout
ijk = Map[n];
// COMPUTE THE COLOR GRADIENT
//........................................................................
//.................Read Phase Indicator Values............................
//........................................................................
nn = ijk-1; // neighbor index (get convention)
id = ID[nn];
m1 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 1
//........................................................................
nn = ijk+1; // neighbor index (get convention)
id = ID[nn];
m2 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 2
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
id = ID[nn];
m3 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 3
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
id = ID[nn];
m4 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 4
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
id = ID[nn];
m5 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 5
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
id = ID[nn];
m6 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 6
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
id = ID[nn];
m7 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 7
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
id = ID[nn];
m8 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 8
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
id = ID[nn];
m9 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 9
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
id = ID[nn];
m10 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 10
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
id = ID[nn];
m11 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 11
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
id = ID[nn];
m12 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 12
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
id = ID[nn];
m13 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 13
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
id = ID[nn];
m14 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 14
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
id = ID[nn];
m15 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 15
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
id = ID[nn];
m16 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 16
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
id = ID[nn];
m17 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+Psi[ijk]*(id<=0);// get neighbor for phi - 17
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
id = ID[nn];
m18 = SolidBC==1 ? Psi[nn] : Psi[nn]*(id>0)+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));
ElectricField[n] = nx;
ElectricField[Np+n] = ny;
ElectricField[2*Np+n] = nz;
}
}

162
cpu/Stokes.cpp
View File

@@ -1,7 +1,6 @@
#include <stdio.h>
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz, double Ex_const, double Ey_const, double Ez_const, int start, int finish, int Np)
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np)
{
double fq;
// conserved momemnts
@@ -32,13 +31,13 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
//Load data
rhoE = ChargeDensity[n];
Ex = ElectricField[n+0*Np]+Ex_const;
Ey = ElectricField[n+1*Np]+Ey_const;
Ez = ElectricField[n+2*Np]+Ez_const;
Ex = ElectricField[n+0*Np];
Ey = ElectricField[n+1*Np];
Ez = ElectricField[n+2*Np];
//compute total body force, including input body force (Gx,Gy,Gz)
Fx = Gx + rhoE*Ex;
Fy = Gy + rhoE*Ey;
Fz = Gz + rhoE*Ez;
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
// q=0
fq = dist[n];
@@ -311,9 +310,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
m18 -= fq;
// write the velocity
ux = jx / rho;
uy = jy / rho;
uz = jz / rho;
ux = jx / rho0;
uy = jy / rho0;
uz = jz / rho0;
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
@@ -326,18 +325,18 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
//..............incorporate external force................................................
//..............carry out relaxation process...............................................
m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho) - m2);
m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0) - m2);
m4 = m4 + rlx_setB*((-0.6666666666666666*jx) - m4);
m6 = m6 + rlx_setB*((-0.6666666666666666*jy) - m6);
m8 = m8 + rlx_setB*((-0.6666666666666666*jz) - m8);
m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) - m9);
m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) - m9);
m10 = m10 + rlx_setA*(-0.5*((2*jx*jx-jy*jy-jz*jz)/rho) - m10);
m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) - m11);
m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho) - m12);
m13 = m13 + rlx_setA*((jx*jy/rho) - m13);
m14 = m14 + rlx_setA*((jy*jz/rho) - m14);
m15 = m15 + rlx_setA*((jx*jz/rho) - m15);
m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) - m11);
m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho0) - m12);
m13 = m13 + rlx_setA*((jx*jy/rho0) - m13);
m14 = m14 + rlx_setA*((jy*jz/rho0) - m14);
m15 = m15 + rlx_setA*((jx*jz/rho0) - m15);
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
@@ -454,8 +453,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, do
}
}
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz, double Ex_const, double Ey_const, double Ez_const, int start, int finish, int Np)
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np)
{
double fq;
// conserved momemnts
@@ -487,13 +485,13 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
//Load data
rhoE = ChargeDensity[n];
Ex = ElectricField[n+0*Np]+Ex_const;
Ey = ElectricField[n+1*Np]+Ey_const;
Ez = ElectricField[n+2*Np]+Ez_const;
Ex = ElectricField[n+0*Np];
Ey = ElectricField[n+1*Np];
Ez = ElectricField[n+2*Np];
//compute total body force, including input body force (Gx,Gy,Gz)
Fx = Gx + rhoE*Ex;
Fy = Gy + rhoE*Ey;
Fz = Gz + rhoE*Ez;
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
// q=0
fq = dist[n];
@@ -803,27 +801,27 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
m18 -= fq;
// write the velocity
ux = jx / rho;
uy = jy / rho;
uz = jz / rho;
ux = jx / rho0;
uy = jy / rho0;
uz = jz / rho0;
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
//..............incorporate external force................................................
//..............carry out relaxation process...............................................
m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho) - m2);
m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0) - m2);
m4 = m4 + rlx_setB*((-0.6666666666666666*jx) - m4);
m6 = m6 + rlx_setB*((-0.6666666666666666*jy) - m6);
m8 = m8 + rlx_setB*((-0.6666666666666666*jz) - m8);
m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) - m9);
m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) - m9);
m10 = m10 + rlx_setA*(-0.5*((2*jx*jx-jy*jy-jz*jz)/rho) - m10);
m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) - m11);
m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho) - m12);
m13 = m13 + rlx_setA*((jx*jy/rho) - m13);
m14 = m14 + rlx_setA*((jy*jz/rho) - m14);
m15 = m15 + rlx_setA*((jx*jz/rho) - m15);
m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) - m11);
m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho0) - m12);
m13 = m13 + rlx_setA*((jx*jy/rho0) - m13);
m14 = m14 + rlx_setA*((jy*jz/rho0) - m14);
m15 = m15 + rlx_setA*((jx*jz/rho0) - m15);
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
@@ -955,47 +953,47 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
}
}
extern "C" void ScaLBL_D3Q19_Momentum_Phys(double *dist, double *vel, double h, double time_conv, int Np)
{
//h: resolution [um/lu]
//time_conv: time conversion factor [sec/lt]
int n;
// distributions
double f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double vx,vy,vz;
for (n=0; n<Np; n++){
//........................................................................
// Registers to store the distributions
//........................................................................
f2 = dist[2*Np+n];
f4 = dist[4*Np+n];
f6 = dist[6*Np+n];
f8 = dist[8*Np+n];
f10 = dist[10*Np+n];
f12 = dist[12*Np+n];
f14 = dist[14*Np+n];
f16 = dist[16*Np+n];
f18 = dist[18*Np+n];
//........................................................................
f1 = dist[Np+n];
f3 = dist[3*Np+n];
f5 = dist[5*Np+n];
f7 = dist[7*Np+n];
f9 = dist[9*Np+n];
f11 = dist[11*Np+n];
f13 = dist[13*Np+n];
f15 = dist[15*Np+n];
f17 = dist[17*Np+n];
//.................Compute the velocity...................................
vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
//..................Write the velocity.....................................
vel[0*Np+n] = vx*(h*1.0e-6)/time_conv;
vel[1*Np+n] = vy*(h*1.0e-6)/time_conv;
vel[2*Np+n] = vz*(h*1.0e-6)/time_conv;
//........................................................................
}
}
//extern "C" void ScaLBL_D3Q19_Momentum_Phys(double *dist, double *vel, double h, double time_conv, int Np)
//{
// //h: resolution [um/lu]
// //time_conv: time conversion factor [sec/lt]
// int n;
// // distributions
// double f1,f2,f3,f4,f5,f6,f7,f8,f9;
// double f10,f11,f12,f13,f14,f15,f16,f17,f18;
// double vx,vy,vz;
//
// for (n=0; n<Np; n++){
// //........................................................................
// // Registers to store the distributions
// //........................................................................
// f2 = dist[2*Np+n];
// f4 = dist[4*Np+n];
// f6 = dist[6*Np+n];
// f8 = dist[8*Np+n];
// f10 = dist[10*Np+n];
// f12 = dist[12*Np+n];
// f14 = dist[14*Np+n];
// f16 = dist[16*Np+n];
// f18 = dist[18*Np+n];
// //........................................................................
// f1 = dist[Np+n];
// f3 = dist[3*Np+n];
// f5 = dist[5*Np+n];
// f7 = dist[7*Np+n];
// f9 = dist[9*Np+n];
// f11 = dist[11*Np+n];
// f13 = dist[13*Np+n];
// f15 = dist[15*Np+n];
// f17 = dist[17*Np+n];
// //.................Compute the velocity...................................
// vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
// vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
// vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
// //..................Write the velocity.....................................
// vel[0*Np+n] = vx*(h*1.0e-6)/time_conv;
// vel[1*Np+n] = vy*(h*1.0e-6)/time_conv;
// vel[2*Np+n] = vz*(h*1.0e-6)/time_conv;
// //........................................................................
// }
//}

31
models/IonModel.cpp
View File

@@ -364,11 +364,16 @@ void ScaLBL_IonModel::Initialize(){
void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//Input parameter:
//1. Velocity is from StokesModel
//2. ElectricField is from Poisson model
//LB-related parameter
vector<double> rlx(tau.begin(),tau.end());
for (double item : rlx){
item = 1.0/item;
}
//.......create and start timer............
//double starttime,stoptime,cputime;
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
@@ -462,7 +467,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
}
void ScaLBL_IonModel::getIonConcentration(){
//TODO this ruin the ion concentration on device
//need to do something similar to electric field
void ScaLBL_IonModel::getIonConcentration(int timestep){
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_IonConcentration_Phys(Ci, h, ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_IonConcentration_Phys(Ci, h, ic, 0, ScaLBL_Comm->LastExterior(), Np);
@@ -474,7 +481,7 @@ void ScaLBL_IonModel::getIonConcentration(){
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
FILE *OUTFILE;
sprintf(LocalRankFilename,"Ion%02i.%05i.raw",ic+1,rank);
sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE);
@@ -482,3 +489,23 @@ void ScaLBL_IonModel::getIonConcentration(){
}
//void ScaLBL_IonModel::getIonConcentration(){
// for (int ic=0; ic<number_ion_species; ic++){
// ScaLBL_IonConcentration_Phys(Ci, h, ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
// ScaLBL_IonConcentration_Phys(Ci, h, ic, 0, ScaLBL_Comm->LastExterior(), Np);
// }
//
// DoubleArray PhaseField(Nx,Ny,Nz);
// for (int ic=0; ic<number_ion_species; ic++){
// ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
// ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//
// FILE *OUTFILE;
// sprintf(LocalRankFilename,"Ion%02i.%05i.raw",ic+1,rank);
// OUTFILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,OUTFILE);
// fclose(OUTFILE);
// }
//
//}

3
models/IonModel.h
View File

@@ -29,7 +29,7 @@ public:
void Create();
void Initialize();
void Run(double *Velocity, double *ElectricField);
void getIonConcentration();
void getIonConcentration(int timestep);
//bool Restart,pBC;
int timestep,timestepMax;
@@ -40,6 +40,7 @@ public:
double kb,electron_charge,T,Vt;
double k2_inv;
double tolerance;
double Ex,Ey,Ez;
int number_ion_species;
vector<double> IonDiffusivity;//User input unit [m^2/sec]

383
models/PoissonSolver.cpp
View File

@@ -7,7 +7,7 @@
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),
epsilon0(0),epsilon0_LB(0),epsilonR(0),epsilon_LB(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),
nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),BoundaryConditionSolid(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{
@@ -22,17 +22,20 @@ 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
gamma = 0.3;//time step of LB-Poisson equation
//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;
timestepMax = 100000;
tolerance = 1.0e-6;//stopping criterion for obtaining steady-state electricla potential
h = 1.0;//resolution; unit: um/lu
epsilon0 = 8.85e-12;//electrical permittivity of vaccum; unit:[C/(V*m)]
epsilon0 = 8.85e-12;//electric permittivity of vaccum; unit:[C/(V*m)]
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilonR = 78.4;//default dielectric constant of water
epsilon_LB = epsilon0_LB*epsilonR;//electrical permittivity
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
// LB-Poisson Model parameters
if (electric_db->keyExists( "timestepMax" )){
@@ -56,20 +59,23 @@ void ScaLBL_Poisson::ReadParams(string filename){
if (electric_db->keyExists( "BC_Solid" )){
BoundaryConditionSolid = electric_db->getScalar<int>( "BC_Solid" );
}
// Read boundary condition for electric potentiona
// BC = 0: normal periodic BC
// BC = 1: fixed inlet and outlet potential
BoundaryCondition = 0;
if (electric_db->keyExists( "BC" )){
BoundaryCondition = electric_db->getScalar<int>( "BC" );
}
// Read domain parameters
if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
h = domain_db->getScalar<double>( "voxel_length" );
}
BoundaryCondition = 0;
if (domain_db->keyExists( "BC" )){
BoundaryCondition = domain_db->getScalar<int>( "BC" );
}
//Re-calcualte model parameters if user updates input
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilon_LB = epsilon0_LB*epsilonR;//electrical permittivity
epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity
tau = 0.5+k2_inv*gamma;
if (rank==0) printf("***********************************************************************************\n");
@@ -202,13 +208,13 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid)
AFFINITY=0.f;
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
//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
AFFINITY = AFFINITY*(h*h*1.0e-12);
//TODO maybe there is a factor of gamm missing here ?
AFFINITY = AFFINITY*(h*h*1.0e-12)/epsilon_LB;
}
label_count[idx] += 1.0;
idx = NLABELS;
@@ -244,7 +250,6 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid)
}
}
void ScaLBL_Poisson::Create(){
/*
* This function creates the variables needed to run a LBM
@@ -260,6 +265,7 @@ void ScaLBL_Poisson::Create(){
// Create a communicator for the device (will use optimized layout)
// ScaLBL_Communicator ScaLBL_Comm(Mask); // original
ScaLBL_Comm = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
ScaLBL_Comm_Regular = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
int Npad=(Np/16 + 2)*16;
if (rank==0) printf ("LB-Poisson Solver: Set up memory efficient layout \n");
@@ -277,37 +283,125 @@ void ScaLBL_Poisson::Create(){
int neighborSize=18*(Np*sizeof(int));
//...........................................................................
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 **) &fq, 7*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Np);
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);
//ScaLBL_AllocateDeviceMemory((void **) &PoissonSolid, sizeof(double)*Nx*Ny*Nz);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("LB-Poisson Solver: Setting up device map and neighbor list \n");
fflush(stdout);
int *TmpMap;
TmpMap=new int[Np];
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
if (!(idx < 0))
TmpMap[idx] = k*Nx*Ny+j*Nx+i;
}
}
}
// check that TmpMap is valid
for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
auto n = TmpMap[idx];
if (n > Nx*Ny*Nz){
printf("Bad value! idx=%i \n", n);
TmpMap[idx] = Nx*Ny*Nz-1;
}
}
for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
auto n = TmpMap[idx];
if ( n > Nx*Ny*Nz ){
printf("Bad value! idx=%i \n",n);
TmpMap[idx] = Nx*Ny*Nz-1;
}
}
ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
ScaLBL_DeviceBarrier();
delete [] TmpMap;
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
ScaLBL_DeviceBarrier();
MPI_Barrier(comm);
delete [] neighborList;
// copy node ID
ScaLBL_CopyToDevice(dvcID, Mask->id, sizeof(signed char)*Nx*Ny*Nz);
ScaLBL_DeviceBarrier();
//Initialize solid boundary for electrical potential
//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;
//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){
if (BoundaryCondition==1){
if (electric_db->keyExists( "Vin" )){
Vin = electric_db->getScalar<double>( "Vin" );
}
if (electric_db->keyExists( "Vout" )){
Vout = electric_db->getScalar<double>( "Vout" );
}
}
//By default only periodic BC is applied and Vin=Vout=1.0, i.e. there is no potential gradient along Z-axis
double slope = (Vout-Vin)/(Nz-2);
double psi_linearized;
for (int k=0;k<Nz;k++){
if (k==0 || k==1){
psi_linearized = Vin;
}
else if (k==Nz-1 || k==Nz-2){
psi_linearized = Vout;
}
else{
psi_linearized = slope*(k-1)+Vin;
}
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
if (Mask->id[n]>0){
psi_init[n] = psi_linearized;
}
}
}
}
}
void ScaLBL_Poisson::Initialize(){
/*
* This function initializes model
*/
if (rank==0) printf ("LB-Poisson Solver: initializing D3Q7 distributions\n");
ScaLBL_D3Q7_Poisson_Init(fq, Np);
//NOTE the initialization involves two steps:
//1. assign solid boundary value (surface potential or surface change density)
//2. Initialize electric potential for pore nodes
double *psi_host;
psi_host = new double [Nx*Ny*Nz];
AssignSolidBoundary(psi_host);//step1
Potential_Init(psi_host);//step2
ScaLBL_CopyToDevice(Psi, psi_host, Nx*Ny*Nz*sizeof(double));
ScaLBL_DeviceBarrier();
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;
}
void ScaLBL_Poisson::Run(double *ChargeDensity){
@@ -325,30 +419,83 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
// *************ODD TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
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
/* ... */
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->SolidDirichletD3Q7(fq, PoissonSolid);
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);
//compute electric field
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);
//perform collision
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, 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);
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// *************EVEN TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FORM NORMAL
ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
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
/* ... */
ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->SolidDirichletD3Q7(fq, PoissonSolid);
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);
//compute electric field
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);
//perform collision
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, 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);
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/
// Check convergence of steady-state solution
if (timestep%analysis_interval==0){
ScaLBL_Comm->RegularLayout(Map,Psi,Psi_host);
//ScaLBL_Comm->RegularLayout(Map,Psi,Psi_host);
ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
double count_loc=0;
double count;
double psi_avg;
@@ -393,14 +540,180 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
}
void ScaLBL_Poisson::getElectricalPotential(){
void ScaLBL_Poisson::getElectricPotential(int timestep){
DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
//ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
ScaLBL_CopyToHost(PhaseField.data(),Psi,sizeof(double)*Nx*Ny*Nz);
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
FILE *OUTFILE;
sprintf(LocalRankFilename,"Electrical_Potential.%05i.raw",rank);
sprintf(LocalRankFilename,"Electric_Potential_Time_%i.%05i.raw",timestep,rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE);
}
void ScaLBL_Poisson::getElectricField(int timestep){
//ScaLBL_D3Q7_Poisson_getElectricField(fq,ElectricField,tau,Np);
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map,&ElectricField[0*Np],PhaseField);
ElectricField_LB_to_Phys(PhaseField);
FILE *EX;
sprintf(LocalRankFilename,"ElectricField_X_Time_%i.%05i.raw",timestep,rank);
EX = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,EX);
fclose(EX);
ScaLBL_Comm->RegularLayout(Map,&ElectricField[1*Np],PhaseField);
ElectricField_LB_to_Phys(PhaseField);
FILE *EY;
sprintf(LocalRankFilename,"ElectricField_Y_Time_%i.%05i.raw",timestep,rank);
EY = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,EY);
fclose(EY);
ScaLBL_Comm->RegularLayout(Map,&ElectricField[2*Np],PhaseField);
ElectricField_LB_to_Phys(PhaseField);
FILE *EZ;
sprintf(LocalRankFilename,"ElectricField_Z_Time_%i.%05i.raw",timestep,rank);
EZ = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,EZ);
fclose(EZ);
}
void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){
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)){
Efield_reg(i,j,k) = Efield_reg(i,j,k)/(h*1.0e-6);
}
}
}
}
}
//void ScaLBL_Poisson::getElectricPotential(){
//
// DoubleArray PhaseField(Nx,Ny,Nz);
// ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
// //ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// FILE *OUTFILE;
// sprintf(LocalRankFilename,"Electric_Potential.%05i.raw",rank);
// OUTFILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,OUTFILE);
// fclose(OUTFILE);
//}
//old version where Psi is of size Np
//void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid)
//{
// size_t NLABELS=0;
// signed char VALUE=0;
// double AFFINITY=0.f;
//
// auto LabelList = electric_db->getVector<int>( "SolidLabels" );
// auto AffinityList = electric_db->getVector<double>( "SolidValues" );
//
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: LB-Poisson Solver: SolidLabels and SolidValues must be the same length! \n");
// }
//
// double label_count[NLABELS];
// double label_count_global[NLABELS];
// // Assign the labels
//
// for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=Mask->id[n];
// AFFINITY=0.f;
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// //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;
// idx = NLABELS;
// //Mask->id[n] = 0; // set mask to zero since this is an immobile component
// }
// }
// poisson_solid[n] = AFFINITY;
// }
// }
// }
//
// for (size_t idx=0; idx<NLABELS; idx++)
// label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
//
// if (rank==0){
// printf("LB-Poisson Solver: number of Poisson solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// VALUE=LabelList[idx];
// AFFINITY=AffinityList[idx];
// double volume_fraction = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
// switch (BoundaryConditionSolid){
// case 1:
// printf(" label=%d, surface potential=%.3g [V], volume fraction=%.2g\n",VALUE,AFFINITY,volume_fraction);
// break;
// case 2:
// printf(" label=%d, surface charge density=%.3g [C/m^2], volume fraction=%.2g\n",VALUE,AFFINITY,volume_fraction);
// break;
// default:
// printf(" label=%d, surface potential=%.3g [V], volume fraction=%.2g\n",VALUE,AFFINITY,volume_fraction);
// break;
// }
// }
// }
//}
// old version where Psi is of size Np
//void ScaLBL_Poisson::Potential_Init(double *psi_init){
//
// if (BoundaryCondition==1){
// if (electric_db->keyExists( "Vin" )){
// Vin = electric_db->getScalar<double>( "Vin" );
// }
// if (electric_db->keyExists( "Vout" )){
// Vout = electric_db->getScalar<double>( "Vout" );
// }
// }
// //By default only periodic BC is applied and Vin=Vout=1.0, i.e. there is no potential gradient along Z-axis
// double slope = (Vout-Vin)/(Nz-2);
// double psi_linearized;
// for (int k=0;k<Nz;k++){
// if (k==0 || k==1){
// psi_linearized = Vin;
// }
// else if (k==Nz-1 || k==Nz-2){
// psi_linearized = Vout;
// }
// else{
// psi_linearized = slope*(k-1)+Vin;
// }
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int idx = Map(i,j,k);
// if (!(idx < 0)){
// psi_init[idx] = psi_linearized;
// }
// }
// }
// }
//}

11
models/PoissonSolver.h
View File

@@ -28,7 +28,8 @@ public:
void Create();
void Initialize();
void Run(double *ChargeDensity);
void getElectricalPotential();
void getElectricPotential(int timestep);
void getElectricField(int timestep);
//bool Restart,pBC;
int timestep,timestepMax;
@@ -39,6 +40,7 @@ public:
double tolerance;
double k2_inv,gamma;
double epsilon0,epsilon0_LB,epsilonR,epsilon_LB;
double Vin, Vout;
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
@@ -48,6 +50,7 @@ public:
std::shared_ptr<Domain> Dm; // this domain is for analysis
std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
// input database
std::shared_ptr<Database> db;
std::shared_ptr<Database> domain_db;
@@ -57,10 +60,12 @@ public:
DoubleArray Distance;
DoubleArray Psi_host;
int *NeighborList;
int *dvcMap;
signed char *dvcID;
double *fq;
double *Psi;
double *ElectricField;
double *PoissonSolid;
//double *PoissonSolid;
private:
MPI_Comm comm;
@@ -73,4 +78,6 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void AssignSolidBoundary(double *poisson_solid);
void Potential_Init(double *psi_init);
void ElectricField_LB_to_Phys(DoubleArray &Efield_reg);
};

90
models/StokesModel.cpp
View File

@@ -7,7 +7,7 @@
ScaLBL_StokesModel::ScaLBL_StokesModel(int RANK, int NP, MPI_Comm COMM):
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),mu(0),h(0),nu_phys(0),time_conv(0),tolerance(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),mu(0),h(0),nu_phys(0),rho_phys(0),rho0(0),den_scale(0),time_conv(0),tolerance(0),
Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{
@@ -27,17 +27,17 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
//-------------------------------------------------------------------//
//---------------------- Default model parameters --------------------------//
rho_phys = 1000.0; //by default use water density; unit [kg/m^3]
nu_phys = 1.004e-6;//by default use water kinematic viscosity at 20C; unit [m^2/sec]
h = 1.0;//image resolution;[um]
tau = 1.0;
mu = (tau-0.5)/3.0;//LB kinematic viscosity;unit [lu^2/lt]
time_conv = h*h*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
rho0 = 1.0;//LB density
den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
tolerance = 1.0e-8;
Fx = Fy = 0.0;
Fz = 1.0e-5;
//Body electric field [V/lu]
Ex = Ey = 0.0;
Ez = 1.0e-3;
//--------------------------------------------------------------------------//
// Read domain parameters
@@ -59,19 +59,20 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
if (stokes_db->keyExists( "tau" )){
tau = stokes_db->getScalar<double>( "tau" );
}
if (stokes_db->keyExists( "rho0" )){
rho0 = stokes_db->getScalar<double>( "rho0" );
}
if (stokes_db->keyExists( "nu_phys" )){
nu_phys = stokes_db->getScalar<double>( "nu_phys" );
}
if (stokes_db->keyExists( "rho_phys" )){
rho_phys = stokes_db->getScalar<double>( "rho_phys" );
}
if (stokes_db->keyExists( "F" )){
Fx = stokes_db->getVector<double>( "F" )[0];
Fy = stokes_db->getVector<double>( "F" )[1];
Fz = stokes_db->getVector<double>( "F" )[2];
}
if (stokes_db->keyExists( "ElectricField" )){//NOTE user-input has physical unit [V/m]
Ex = stokes_db->getVector<double>( "ElectricField" )[0];
Ey = stokes_db->getVector<double>( "ElectricField" )[1];
Ez = stokes_db->getVector<double>( "ElectricField" )[2];
}
if (stokes_db->keyExists( "Restart" )){
Restart = stokes_db->getScalar<bool>( "Restart" );
}
@@ -88,10 +89,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*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);
Ez = Ez*(h*1.0e-6);
den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
if (rank==0) printf("*****************************************************\n");
if (rank==0) printf("LB Single-Fluid Navier-Stokes Solver: \n");
@@ -246,7 +244,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
while (timestep < timestepMax) {
//************************************************************************/
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez,
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
@@ -262,12 +260,14 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
if (BoundaryCondition == 3){
@@ -282,41 +282,87 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, Ex, Ey, Ez, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/
}
}
void ScaLBL_StokesModel::getVelocity(){
void ScaLBL_StokesModel::getVelocity(int timestep){
//get velocity in physical unit [m/sec]
ScaLBL_D3Q19_Momentum_Phys(fq, Velocity, h, time_conv, Np);
ScaLBL_D3Q19_Momentum(fq, Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
Velocity_LB_to_Phys(PhaseField);
FILE *VELX_FILE;
sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);
sprintf(LocalRankFilename,"Velocity_X_Time_%i.%05i.raw",timestep,rank);
VELX_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,VELX_FILE);
fclose(VELX_FILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
Velocity_LB_to_Phys(PhaseField);
FILE *VELY_FILE;
sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank);
sprintf(LocalRankFilename,"Velocity_Y_Time_%i.%05i.raw",timestep,rank);
VELY_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,VELY_FILE);
fclose(VELY_FILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
Velocity_LB_to_Phys(PhaseField);
FILE *VELZ_FILE;
sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank);
sprintf(LocalRankFilename,"Velocity_Z_Time_%i.%05i.raw",timestep,rank);
VELZ_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,VELZ_FILE);
fclose(VELZ_FILE);
}
void ScaLBL_StokesModel::Velocity_LB_to_Phys(DoubleArray &Vel_reg){
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)){
Vel_reg(i,j,k) = Vel_reg(i,j,k)*(h*1.0e-6)/time_conv;
}
}
}
}
}
//void ScaLBL_StokesModel::getVelocity(){
// //get velocity in physical unit [m/sec]
// ScaLBL_D3Q19_Momentum_Phys(fq, Velocity, h, time_conv, Np);
// ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//
// DoubleArray PhaseField(Nx,Ny,Nz);
// ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
// FILE *VELX_FILE;
// sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);
// VELX_FILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,VELX_FILE);
// fclose(VELX_FILE);
//
// ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
// FILE *VELY_FILE;
// sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank);
// VELY_FILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,VELY_FILE);
// fclose(VELY_FILE);
//
// ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
// FILE *VELZ_FILE;
// sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank);
// VELZ_FILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,VELZ_FILE);
// fclose(VELZ_FILE);
//
//}
void ScaLBL_StokesModel::Run(){
double rlx_setA=1.0/tau;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);

7
models/StokesModel.h
View File

@@ -30,19 +30,21 @@ public:
void Run();
void Run_Lite(double *ChargeDensity, double *ElectricField);
void VelocityField();
void getVelocity();
void getVelocity(int timestep);
bool Restart,pBC;
int timestep,timestepMax;
int BoundaryCondition;
double tau,mu;
double rho0;
double Fx,Fy,Fz,flux;
double Ex,Ey,Ez;
double din,dout;
double tolerance;
double nu_phys;
double rho_phys;
double time_conv;
double h;//image resolution
double den_scale;//scale factor for density
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
@@ -78,4 +80,5 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void Velocity_LB_to_Phys(DoubleArray &Vel_reg);
};

16
tests/lbpm_electrokinetic_dfh_simulator.cpp
View File

@@ -78,14 +78,17 @@ int main(int argc, char **argv)
while (timestep < Study.timestepMax){
timestep++;
if (rank==0) printf("timestep=%i; running Poisson solver\n",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);
//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);
//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
@@ -94,9 +97,10 @@ int main(int argc, char **argv)
//--------------------------------------------
}
StokesModel.getVelocity();
PoissonSolver.getElectricalPotential();
IonModel.getIonConcentration();
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");