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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
JamesEMcclure 2021-06-01 06:43:07 -04:00
commit baedeaa2a3
16 changed files with 6776 additions and 1682 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
analysis/SubPhase.cpp
View File

@ -364,6 +364,7 @@ inline void InterfaceTransportMeasures( double beta, double rA, double rB, doubl
double A1,A2,A3,A4,A5,A6;
double B1,B2,B3,B4,B5,B6;
double nAB,delta;
double phi = (nA-nB)/(nA+nB);
// Instantiate mass transport distributions
// Stationary value - distribution 0
nAB = 1.0/(nA+nB);
@ -402,7 +403,7 @@ inline void InterfaceTransportMeasures( double beta, double rA, double rB, doubl
double uwx = (B1-B2);
double uwy = (B3-B4);
double uwz = (B5-B6);
/*
I.Mn += rA*nA;
I.Mw += rB*nB;
I.Pnx += rA*nA*unx;
@ -413,6 +414,21 @@ inline void InterfaceTransportMeasures( double beta, double rA, double rB, doubl
I.Pwz += rB*nB*uwz;
I.Kn += rA*nA*(unx*unx + uny*uny + unz*unz);
I.Kw += rB*nB*(uwx*uwx + uwy*uwy + uwz*uwz);
*/
if (phi > 0.0){
I.Mn += rA;
I.Pnx += rA*ux;
I.Pny += rA*uy;
I.Pnz += rA*uz;
}
else {
I.Mw += rB;
I.Pwx += rB*ux;
I.Pwy += rB*uy;
I.Pwz += rB*uz;
}
I.Kn += rA*nA*(unx*unx + uny*uny + unz*unz);
I.Kw += rB*nB*(uwx*uwx + uwy*uwy + uwz*uwz);
}
@ -606,8 +622,8 @@ void SubPhase::Full(){
double uy = Vel_y(n);
double uz = Vel_z(n);
if (DelPhi(n) > 1e-3){
// interface region
if (DelPhi(n) > 1e-3 && SDs(n) < 3.0 ){
// film region
double nx = 0.5*(Phi(i+1,j,k)-Phi(i-1,j,k));
double ny = 0.5*(Phi(i,j+1,k)-Phi(i,j-1,k));
double nz = 0.5*(Phi(i,j,k+1)-Phi(i,j,k-1));

167
common/ScaLBL.cpp
View File

@ -989,7 +989,7 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis
}
void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np)
void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC)
{
int idx,i,j,k;
@ -1067,6 +1067,23 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
int *bb_interactions_tmp = new int [local_count];
ScaLBL_AllocateDeviceMemory((void **) &bb_dist, sizeof(int)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &bb_interactions, sizeof(int)*local_count);
int *fluid_boundary_tmp;
double *lattice_weight_tmp;
float *lattice_cx_tmp;
float *lattice_cy_tmp;
float *lattice_cz_tmp;
if(SlippingVelBC==true){
fluid_boundary_tmp = new int [local_count];
lattice_weight_tmp = new double [local_count];
lattice_cx_tmp = new float [local_count];
lattice_cy_tmp = new float [local_count];
lattice_cz_tmp = new float [local_count];
ScaLBL_AllocateDeviceMemory((void **) &fluid_boundary, sizeof(int)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_weight, sizeof(double)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_cx, sizeof(float)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_cy, sizeof(float)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_cz, sizeof(float)*local_count);
}
local_count=0;
for (k=1;k<Nz-1;k++){
@ -1080,36 +1097,78 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
neighbor=Map(i-1,j,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i-1) + (j)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 0.0;
}
bb_dist_tmp[local_count++]=idx + 2*Np;
}
neighbor=Map(i+1,j,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i+1) + (j)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 0.0;
}
bb_dist_tmp[local_count++] = idx + 1*Np;
}
neighbor=Map(i,j-1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j-1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = 0.0;
}
bb_dist_tmp[local_count++]=idx + 4*Np;
}
neighbor=Map(i,j+1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j+1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = 0.0;
}
bb_dist_tmp[local_count++]=idx + 3*Np;
}
neighbor=Map(i,j,k-1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = -1.0;
}
bb_dist_tmp[local_count++]=idx + 6*Np;
}
neighbor=Map(i,j,k+1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 1.0;
}
bb_dist_tmp[local_count++]=idx + 5*Np;
}
}
@ -1127,72 +1186,156 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
neighbor=Map(i-1,j-1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i-1) + (j-1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = 0.0;
}
bb_dist_tmp[local_count++]=idx + 8*Np;
}
neighbor=Map(i+1,j+1,k);
if (neighbor==-1) {
bb_interactions_tmp[local_count] = (i+1) + (j+1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = 0.0;
}
bb_dist_tmp[local_count++]=idx + 7*Np;
}
neighbor=Map(i-1,j+1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i-1) + (j+1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = 0.0;
}
bb_dist_tmp[local_count++]=idx + 10*Np;
}
neighbor=Map(i+1,j-1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i+1) + (j-1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = 0.0;
}
bb_dist_tmp[local_count++]=idx + 9*Np;
}
neighbor=Map(i-1,j,k-1);
if (neighbor==-1) {
bb_interactions_tmp[local_count] = (i-1) + (j)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = -1.0;
}
bb_dist_tmp[local_count++]=idx + 12*Np;
}
neighbor=Map(i+1,j,k+1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i+1) + (j)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 1.0;
}
bb_dist_tmp[local_count++]=idx + 11*Np;
}
neighbor=Map(i-1,j,k+1);
if (neighbor==-1) {
bb_interactions_tmp[local_count] = (i-1) + (j)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 1.0;
}
bb_dist_tmp[local_count++]=idx + 14*Np;
}
neighbor=Map(i+1,j,k-1);
if (neighbor==-1) {
bb_interactions_tmp[local_count] = (i+1) + (j)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = -1.0;
}
bb_dist_tmp[local_count++]=idx + 13*Np;
}
neighbor=Map(i,j-1,k-1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j-1)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = -1.0;
}
bb_dist_tmp[local_count++]=idx + 16*Np;
}
neighbor=Map(i,j+1,k+1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j+1)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = 1.0;
}
bb_dist_tmp[local_count++]=idx + 15*Np;
}
neighbor=Map(i,j-1,k+1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j-1)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = 1.0;
}
bb_dist_tmp[local_count++]=idx + 18*Np;
}
neighbor=Map(i,j+1,k-1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j+1)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = -1.0;
}
bb_dist_tmp[local_count++]=idx + 17*Np;
}
}
@ -1202,10 +1345,24 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
n_bb_d3q19 = local_count; // this gives the d3q19 distributions not part of d3q7 model
ScaLBL_CopyToDevice(bb_dist, bb_dist_tmp, local_count*sizeof(int));
ScaLBL_CopyToDevice(bb_interactions, bb_interactions_tmp, local_count*sizeof(int));
if(SlippingVelBC==true){
ScaLBL_CopyToDevice(fluid_boundary, fluid_boundary_tmp, local_count*sizeof(int));
ScaLBL_CopyToDevice(lattice_weight, lattice_weight_tmp, local_count*sizeof(double));
ScaLBL_CopyToDevice(lattice_cx, lattice_cx_tmp, local_count*sizeof(float));
ScaLBL_CopyToDevice(lattice_cy, lattice_cy_tmp, local_count*sizeof(float));
ScaLBL_CopyToDevice(lattice_cz, lattice_cz_tmp, local_count*sizeof(float));
}
ScaLBL_DeviceBarrier();
delete [] bb_dist_tmp;
delete [] bb_interactions_tmp;
if(SlippingVelBC==true){
delete [] fluid_boundary_tmp;
delete [] lattice_weight_tmp;
delete [] lattice_cx_tmp;
delete [] lattice_cy_tmp;
delete [] lattice_cz_tmp;
}
}
void ScaLBL_Communicator::SolidDirichletD3Q7(double *fq, double *BoundaryValue){
@ -1220,6 +1377,14 @@ void ScaLBL_Communicator::SolidNeumannD3Q7(double *fq, double *BoundaryValue){
ScaLBL_Solid_Neumann_D3Q7(fq,BoundaryValue,bb_dist,bb_interactions,n_bb_d3q7);
}
void ScaLBL_Communicator::SolidSlippingVelocityBCD3Q19(double *fq, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0, double den_scale,double h, double time_conv){
// fq is a D3Q19 distribution
// BoundaryValues is a list of values to assign at bounce-back solid sites
ScaLBL_Solid_SlippingVelocityBC_D3Q19(fq,zeta_potential,ElectricField,SolidGrad,epsilon_LB,tau,rho0,den_scale,h,time_conv,
bb_dist,bb_interactions,fluid_boundary,lattice_weight,lattice_cx,lattice_cy,lattice_cz,n_bb_d3q19,N);
}
void ScaLBL_Communicator::SendD3Q19AA(double *dist){
// NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2

21
common/ScaLBL.h
View File

@ -252,6 +252,14 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(int *Map, double *dist
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined_HigherOrder(int *neighborList, int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined_HigherOrder(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,
double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np);
@ -296,6 +304,12 @@ 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_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
int count, int Np);
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);
@ -375,9 +389,11 @@ public:
void RecvHalo(double *data);
void RecvGrad(double *Phi, double *Gradient);
void RegularLayout(IntArray map, const double *data, DoubleArray &regdata);
void SetupBounceBackList(IntArray &Map, signed char *id, int Np);
void SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC=false);
void SolidDirichletD3Q7(double *fq, double *BoundaryValue);
void SolidNeumannD3Q7(double *fq, double *BoundaryValue);
void SolidSlippingVelocityBCD3Q19(double *fq, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epslion_LB, double tau, double rho0, double den_scale,double h, double time_conv);
// Routines to set boundary conditions
void Color_BC_z(int *Map, double *Phi, double *Den, double vA, double vB);
@ -452,6 +468,9 @@ private:
//......................................................................................
int *bb_dist;
int *bb_interactions;
int *fluid_boundary;
double *lattice_weight;
float *lattice_cx, *lattice_cy, *lattice_cz;
//......................................................................................
};

51
cpu/D3Q7BC.cpp
View File

@ -30,6 +30,57 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int
}
}
extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
int count, int Np){
int idx;
int iq,ib,ifluidBC;
double value_b,value_q;
double Ex,Ey,Ez;
double Etx,Ety,Etz;//tangential part of electric field
double E_mag_normal;
double nsx,nsy,nsz;//unit normal solid gradient
double ubx,uby,ubz;//slipping velocity at fluid boundary nodes
float cx,cy,cz;//lattice velocity (D3Q19)
double LB_weight;//lattice weighting coefficient (D3Q19)
double cs2_inv = 3.0;//inverse of cs^2 for D3Q19
double nu_LB = (tau-0.5)/cs2_inv;
for (idx=0; idx<count; idx++){
iq = BounceBackDist_list[idx];
ib = BounceBackSolid_list[idx];
ifluidBC = FluidBoundary_list[idx];
value_b = zeta_potential[ib];//get zeta potential from a solid site
value_q = dist[iq];
//Load electric field and compute its tangential componet
Ex = ElectricField[ifluidBC+0*Np];
Ey = ElectricField[ifluidBC+1*Np];
Ez = ElectricField[ifluidBC+2*Np];
nsx = SolidGrad[ifluidBC+0*Np];
nsy = SolidGrad[ifluidBC+1*Np];
nsz = SolidGrad[ifluidBC+2*Np];
E_mag_normal = Ex*nsx+Ey*nsy+Ez*nsz;//magnitude of electric field in the direction normal to solid nodes
//compute tangential electric field
Etx = Ex - E_mag_normal*nsx;
Ety = Ey - E_mag_normal*nsy;
Etz = Ez - E_mag_normal*nsz;
ubx = -epsilon_LB*value_b*Etx/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
uby = -epsilon_LB*value_b*Ety/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
ubz = -epsilon_LB*value_b*Etz/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
//compute bounce-back distribution
LB_weight = lattice_weight[idx];
cx = lattice_cx[idx];
cy = lattice_cy[idx];
cz = lattice_cz[idx];
dist[iq] = value_q - 2.0*LB_weight*rho0*cs2_inv*(cx*ubx+cy*uby+cz*ubz);
}
}
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];

2925
cpu/FreeLee.cpp
View File

File diff suppressed because it is too large Load Diff

68
cuda/D3Q7BC.cu
View File

@ -38,6 +38,57 @@ __global__ void dvc_ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValu
}
}
__global__ void dvc_ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
int count, int Np)
{
int idx;
int iq,ib,ifluidBC;
double value_b,value_q;
double Ex,Ey,Ez;
double Etx,Ety,Etz;//tangential part of electric field
double E_mag_normal;
double nsx,nsy,nsz;//unit normal solid gradient
double ubx,uby,ubz;//slipping velocity at fluid boundary nodes
float cx,cy,cz;//lattice velocity (D3Q19)
double LB_weight;//lattice weighting coefficient (D3Q19)
double cs2_inv = 3.0;//inverse of cs^2 for D3Q19
double nu_LB = (tau-0.5)/cs2_inv;
idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
iq = BounceBackDist_list[idx];
ib = BounceBackSolid_list[idx];
ifluidBC = FluidBoundary_list[idx];
value_b = zeta_potential[ib];//get zeta potential from a solid site
value_q = dist[iq];
//Load electric field and compute its tangential componet
Ex = ElectricField[ifluidBC+0*Np];
Ey = ElectricField[ifluidBC+1*Np];
Ez = ElectricField[ifluidBC+2*Np];
nsx = SolidGrad[ifluidBC+0*Np];
nsy = SolidGrad[ifluidBC+1*Np];
nsz = SolidGrad[ifluidBC+2*Np];
E_mag_normal = Ex*nsx+Ey*nsy+Ez*nsz;//magnitude of electric field in the direction normal to solid nodes
//compute tangential electric field
Etx = Ex - E_mag_normal*nsx;
Ety = Ey - E_mag_normal*nsy;
Etz = Ez - E_mag_normal*nsz;
ubx = -epsilon_LB*value_b*Etx/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
uby = -epsilon_LB*value_b*Ety/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
ubz = -epsilon_LB*value_b*Etz/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
//compute bounce-back distribution
LB_weight = lattice_weight[idx];
cx = lattice_cx[idx];
cy = lattice_cy[idx];
cz = lattice_cz[idx];
dist[iq] = value_q - 2.0*LB_weight*rho0*cs2_inv*(cx*ubx+cy*uby+cz*ubz);
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np)
{
int idx,n;
@ -410,6 +461,23 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValue, i
}
}
extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
int count, int Np){
int GRID = count / 512 + 1;
dvc_ScaLBL_Solid_SlippingVelocityBC_D3Q19<<<GRID,512>>>(dist, zeta_potential, ElectricField, SolidGrad,
epsilon_LB, tau, rho0, den_scale, h, time_conv,
BounceBackDist_list, BounceBackSolid_list, FluidBoundary_list,
lattice_weight, lattice_cx, lattice_cy, lattice_cz,
count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_Solid_SlippingVelocityBC_D3Q19 (kernel): %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z<<<GRID,512>>>(list, dist, Vin, count, Np);

4868
cuda/FreeLee.cu
View File

File diff suppressed because it is too large Load Diff

72
models/ColorModel.cpp
View File

@ -2102,9 +2102,79 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
return(TOTAL_MASS_CHANGE);
}
void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){
int Np = M.Np;
double dA, dB, phi;
double mass_a, mass_b, mass_a_global, mass_b_global;
double *Aq_tmp, *Bq_tmp;
double *Vel_x, *Vel_y, *Vel_z, *Phase;
Aq_tmp = new double [7*Np];
Bq_tmp = new double [7*Np];
ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
if (dA > 1.0){
double mass_change = dA - 1.0;
Aq_tmp[n] -= 0.333333333333333*mass_change;
Aq_tmp[n+Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
}
if (dB > 1.0){
double mass_change = dB - 1.0;
Bq_tmp[n] -= 0.333333333333333*mass_change;
Bq_tmp[n+Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
}
}
for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){
dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
if (dA > 1.0){
double mass_change = dA - 1.0;
Aq_tmp[n] -= 0.333333333333333*mass_change;
Aq_tmp[n+Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
}
if (dB > 1.0){
double mass_change = dB - 1.0;
Bq_tmp[n] -= 0.333333333333333*mass_change;
Bq_tmp[n+Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
}
}
ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
}
double FlowAdaptor::MoveInterface(ScaLBL_ColorModel &M){
double INTERFACE_CUTOFF = M.color_db->getWithDefault<double>( "move_interface_cutoff", 0.975 );
double INTERFACE_CUTOFF = M.color_db->getWithDefault<double>( "move_interface_cutoff", 0.1 );
double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault<double>( "move_interface_factor", 10.0 );
ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) );

1
models/ColorModel.h
View File

@ -114,6 +114,7 @@ public:
~FlowAdaptor();
double MoveInterface(ScaLBL_ColorModel &M);
double UpdateFractionalFlow(ScaLBL_ColorModel &M);
void Flatten(ScaLBL_ColorModel &M);
DoubleArray phi;
DoubleArray phi_t;
private:

33
models/FreeLeeModel.cpp
View File

@ -59,6 +59,30 @@ void ScaLBL_FreeLeeModel::getVelocity(DoubleArray &Vel_x, DoubleArray &Vel_y, Do
ScaLBL_Comm->Barrier(); comm.barrier();
}
void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray &regdata){
// Gets data (in optimized layout) from the HOST and stores in regular layout
// Primarly for debugging
int i,j,k,idx;
int n;
// initialize the array
regdata.fill(0.f);
double value;
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
n=k*Nx*Ny+j*Nx+i;
idx=Map(i,j,k);
if (!(idx<0)){
value=data[idx];
regdata(i,j,k)=value;
}
}
}
}
}
void ScaLBL_FreeLeeModel::ReadParams(string filename){
// read the input database
db = std::make_shared<Database>( filename );
@ -574,27 +598,28 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels_ChemPotential_ColorGrad()
DoubleArray PhaseField(Nx,Ny,Nz);
FILE *OUTFILE;
ScaLBL_Comm->RegularLayout(Map,mu_phi_host,PhaseField);
getData_RegularLayout(mu_phi_host,PhaseField);
sprintf(LocalRankFilename,"Chem_Init.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE);
ScaLBL_Comm->RegularLayout(Map,&ColorGrad_host[0],PhaseField);
getData_RegularLayout(&ColorGrad_host[0],PhaseField);
FILE *CGX_FILE;
sprintf(LocalRankFilename,"Gradient_X_Init.%05i.raw",rank);
CGX_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CGX_FILE);
fclose(CGX_FILE);
ScaLBL_Comm->RegularLayout(Map,&ColorGrad_host[Np],PhaseField);
getData_RegularLayout(&ColorGrad_host[Np],PhaseField);
FILE *CGY_FILE;
sprintf(LocalRankFilename,"Gradient_Y_Init.%05i.raw",rank);
CGY_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CGY_FILE);
fclose(CGY_FILE);
ScaLBL_Comm->RegularLayout(Map,&ColorGrad_host[2*Np],PhaseField);
getData_RegularLayout(&ColorGrad_host[2*Np],PhaseField);
FILE *CGZ_FILE;
sprintf(LocalRankFilename,"Gradient_Z_Init.%05i.raw",rank);
CGZ_FILE = fopen(LocalRankFilename,"wb");

1
models/FreeLeeModel.h
View File

@ -84,6 +84,7 @@ public:
void getPhase(DoubleArray &PhaseValues);
void getPotential(DoubleArray &PressureValues, DoubleArray &MuValues);
void getVelocity(DoubleArray &Vx, DoubleArray &Vy, DoubleArray &Vz);
void getData_RegularLayout(const double *data, DoubleArray &regdata);
DoubleArray SignDist;

4
models/IonModel.cpp
View File

@ -831,8 +831,8 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
if (BoundaryConditionSolid==1){
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic*Np*7], IonSolid);
ScaLBL_Comm->Barrier(); comm.barrier();
}
ScaLBL_Comm->Barrier(); comm.barrier();
// *************EVEN TIMESTEP*************//
timestep++;
@ -875,8 +875,8 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
if (BoundaryConditionSolid==1){
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic*Np*7], IonSolid);
ScaLBL_Comm->Barrier(); comm.barrier();
}
ScaLBL_Comm->Barrier(); comm.barrier();
}
}

1
models/PoissonSolver.cpp
View File

@ -114,7 +114,6 @@ void ScaLBL_Poisson::ReadParams(string filename){
h = domain_db->getScalar<double>( "voxel_length" );
}
//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

216
models/StokesModel.cpp
View File

@ -8,6 +8,7 @@
ScaLBL_StokesModel::ScaLBL_StokesModel(int RANK, int NP, const Utilities::MPI& 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),rho_phys(0),rho0(0),den_scale(0),time_conv(0),tolerance(0),
epsilon0(0),epsilon0_LB(0),epsilonR(0),epsilon_LB(0),UseSlippingVelBC(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)
{
@ -38,6 +39,12 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
tolerance = 1.0e-8;
Fx = Fy = 0.0;
Fz = 1.0e-5;
//Stokes solver also needs the following parameters for slipping velocity BC
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;//electric permittivity
UseSlippingVelBC = false;
//--------------------------------------------------------------------------//
// Read domain parameters
@ -85,12 +92,19 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
if (stokes_db->keyExists( "flux" )){
flux = stokes_db->getScalar<double>( "flux" );
}
if (stokes_db->keyExists( "UseElectroosmoticVelocityBC" )){
UseSlippingVelBC = stokes_db->getScalar<bool>( "UseElectroosmoticVelocityBC" );
}
if (stokes_db->keyExists( "epsilonR" )){
epsilonR = stokes_db->getScalar<double>( "epsilonR" );
}
// 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]
den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity
}
void ScaLBL_StokesModel::ReadParams(string filename){
@ -100,7 +114,6 @@ void ScaLBL_StokesModel::ReadParams(string filename){
db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" );
stokes_db = db->getDatabase( "Stokes" );
//---------------------- Default model parameters --------------------------//
rho_phys = 1000.0; //by default use water density; unit [kg/m^3]
@ -114,6 +127,12 @@ void ScaLBL_StokesModel::ReadParams(string filename){
tolerance = 1.0e-8;
Fx = Fy = 0.0;
Fz = 1.0e-5;
//Stokes solver also needs the following parameters for slipping velocity BC
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;//electric permittivity
UseSlippingVelBC = false;
//--------------------------------------------------------------------------//
// Read domain parameters
@ -161,12 +180,19 @@ void ScaLBL_StokesModel::ReadParams(string filename){
if (stokes_db->keyExists( "flux" )){
flux = stokes_db->getScalar<double>( "flux" );
}
if (stokes_db->keyExists( "UseElectroosmoticVelocityBC" )){
UseSlippingVelBC = stokes_db->getScalar<bool>( "UseElectroosmoticVelocityBC" );
}
if (stokes_db->keyExists( "epsilonR" )){
epsilonR = stokes_db->getScalar<double>( "epsilonR" );
}
// 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]
den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity
}
void ScaLBL_StokesModel::SetDomain(){
@ -258,6 +284,159 @@ void ScaLBL_StokesModel::ReadInput(){
if (rank == 0) cout << " Domain set." << endl;
}
void ScaLBL_StokesModel::AssignZetaPotentialSolid(double *zeta_potential_solid)
{
size_t NLABELS=0;
signed char VALUE=0;
double AFFINITY=0.f;
auto LabelList = stokes_db->getVector<int>( "SolidLabels" );
auto AffinityList = stokes_db->getVector<double>( "ZetaPotentialSolidList" );
NLABELS=LabelList.size();
if (NLABELS != AffinityList.size()){
ERROR("Error: LB Stokes Solver: SolidLabels and ZetaPotentialSolidList must be the same length! \n");
}
double label_count[NLABELS];
double label_count_global[NLABELS];
for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
// Assign the labels
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++){
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];//no need to convert unit for zeta potential (i.e. volt)
label_count[idx] += 1.0;
idx = NLABELS;
}
}
zeta_potential_solid[n] = AFFINITY;
}
}
}
for (size_t idx=0; idx<NLABELS; idx++)
label_count_global[idx]=Dm->Comm.sumReduce( label_count[idx]);
if (rank==0){
printf("LB Stokes Solver: number of 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);
printf(" label=%d, zeta potential=%.3g [V], volume fraction=%.2g\n",VALUE,AFFINITY,volume_fraction);
}
}
}
void ScaLBL_StokesModel::AssignSolidGrad(double *solid_grad)
{
double *Dst;
Dst = new double [3*3*3];
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
}
}
}
//implement a D3Q19 lattice
double w_face = 1.0/18.0;
double w_edge = 0.5*w_face;
double w_corner = 0.0;
//local
Dst[13] = 0.f;
//faces
Dst[4] = w_face;
Dst[10] = w_face;
Dst[12] = w_face;
Dst[14] = w_face;
Dst[16] = w_face;
Dst[22] = w_face;
// corners
Dst[0] = w_corner;
Dst[2] = w_corner;
Dst[6] = w_corner;
Dst[8] = w_corner;
Dst[18] = w_corner;
Dst[20] = w_corner;
Dst[24] = w_corner;
Dst[26] = w_corner;
// edges
Dst[1] = w_edge;
Dst[3] = w_edge;
Dst[5] = w_edge;
Dst[7] = w_edge;
Dst[9] = w_edge;
Dst[11] = w_edge;
Dst[15] = w_edge;
Dst[17] = w_edge;
Dst[19] = w_edge;
Dst[21] = w_edge;
Dst[23] = w_edge;
Dst[25] = w_edge;
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)){
double phi_x = 0.f;
double phi_y = 0.f;
double phi_z = 0.f;
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
double weight= Dst[index];
int idi=i+ii-1;
int idj=j+jj-1;
int idk=k+kk-1;
if (idi < 0) idi=0;
if (idj < 0) idj=0;
if (idk < 0) idk=0;
if (!(idi < Nx)) idi=Nx-1;
if (!(idj < Ny)) idj=Ny-1;
if (!(idk < Nz)) idk=Nz-1;
int nn = idk*Nx*Ny + idj*Nx + idi;
double vec_x = double(ii-1);
double vec_y = double(jj-1);
double vec_z = double(kk-1);
double GWNS = double(Mask->id[nn]);
//Since the solid unit normal vector is wanted, treat
//wet node as 0.0 and solid node as 1.0
GWNS = (GWNS>0.0) ? 0.0:1.0;
phi_x += GWNS*weight*vec_x;
phi_y += GWNS*weight*vec_y;
phi_z += GWNS*weight*vec_z;
}
}
}
//solid_grad normalization
double phi_mag=sqrt(phi_x*phi_x+phi_y*phi_y+phi_z*phi_z);
if (phi_mag==0.0) phi_mag=1.0;
solid_grad[idx+0*Np] = phi_x/phi_mag;
solid_grad[idx+1*Np] = phi_y/phi_mag;
solid_grad[idx+2*Np] = phi_z/phi_mag;
}
}
}
}
}
void ScaLBL_StokesModel::Create(){
/*
* This function creates the variables needed to run a LBM
@ -301,6 +480,26 @@ void ScaLBL_StokesModel::Create(){
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
comm.barrier();
if (UseSlippingVelBC==true){
ScaLBL_Comm->SetupBounceBackList(Map, Mask->id.data(), Np,1);
comm.barrier();
//For slipping velocity BC, need zeta potential and solid unit normal vector
ScaLBL_AllocateDeviceMemory((void **) &ZetaPotentialSolid, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &SolidGrad, sizeof(double)*3*Np); //unit normal vector of solid nodes
double *ZetaPotentialSolid_host;
ZetaPotentialSolid_host = new double[Nx*Ny*Nz];
AssignZetaPotentialSolid(ZetaPotentialSolid_host);
double *SolidGrad_host;
SolidGrad_host = new double[3*Np];
AssignSolidGrad(SolidGrad_host);
ScaLBL_CopyToDevice(ZetaPotentialSolid, ZetaPotentialSolid_host, Nx*Ny*Nz*sizeof(double));
ScaLBL_CopyToDevice(SolidGrad, SolidGrad_host, 3*Np*sizeof(double));
ScaLBL_Comm->Barrier();
delete [] ZetaPotentialSolid_host;
delete [] SolidGrad_host;
}
}
void ScaLBL_StokesModel::Initialize(){
@ -324,6 +523,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
timestep = 0;
while (timestep < timestepMax) {
//************************************************************************/
//**************ODD TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
@ -344,8 +544,14 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
}
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);
if (UseSlippingVelBC==true){
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(fq, ZetaPotentialSolid, ElectricField, SolidGrad,
epsilon_LB, 1.0/rlx_setA, rho0, den_scale, h, time_conv);
}
ScaLBL_Comm->Barrier(); comm.barrier();
//**************EVEN TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
@ -366,6 +572,10 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
}
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);
if (UseSlippingVelBC==true){
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(fq, ZetaPotentialSolid, ElectricField, SolidGrad,
epsilon_LB, 1.0/rlx_setA, rho0, den_scale, h, time_conv);
}
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
}

6
models/StokesModel.h
View File

@ -51,6 +51,8 @@ public:
double time_conv;
double h;//image resolution
double den_scale;//scale factor for density
double epsilon0,epsilon0_LB,epsilonR,epsilon_LB;//Stokes solver also needs this for slipping velocity BC
bool UseSlippingVelBC;
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
@ -70,6 +72,8 @@ public:
double *fq;
double *Velocity;
double *Pressure;
double *ZetaPotentialSolid;
double *SolidGrad;
//Minkowski Morphology;
DoubleArray Velocity_x;
@ -88,5 +92,7 @@ private:
void LoadParams(std::shared_ptr<Database> db0);
void Velocity_LB_to_Phys(DoubleArray &Vel_reg);
vector<double> computeElectricForceAvg(double *ChargeDensity, double *ElectricField);
void AssignSolidGrad(double *solid_grad);
void AssignZetaPotentialSolid(double *zeta_potential_solid);
};
#endif

2
tests/lbpm_color_simulator.cpp
View File

@ -100,6 +100,8 @@ int main( int argc, char **argv )
Adapt.MoveInterface(ColorModel);
skip_time += ANALYSIS_INTERVAL;
}
//Adapt.Flatten(ColorModel);
}
ColorModel.WriteDebug();
}