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

re-factoring flux BC

Browse Source
This commit is contained in:
James E McClure 2017-12-19 13:00:27 -05:00
parent bee3713349
commit 25c3b89abe
3 changed files with 25 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
cpu/D3Q19.cpp
View File

@ -282,7 +282,7 @@ extern "C" void ScaLBL_D3Q19_Swap_Compact(int *neighborList, double *disteven, d
}
extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID, double *disteven, double *distodd, double flux,
extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID, double *disteven, double *distodd, double Q,
int Nx, int Ny, int Nz){
// Note that this routine assumes the distributions are stored "opposite"
// odd distributions in disteven and even distributions in distodd.
@ -293,7 +293,6 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID, double *disteven, double *di
double din = 0.f;
N = Nx*Ny*Nz;
double A = 1.f*double((Nx-2)*(Ny-2));
double sum = 0.f;
char id;
for (n=Nx*Ny; n<2*Nx*Ny; n++){
@ -331,11 +330,11 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID, double *disteven, double *di
sum += (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f6+f12+f13+f16+f17));
}
}
din = sum/(A*(1.0-flux));
din = sum/(1.0-Q);
return din;
}
extern "C" double ScaLBL_D3Q19_Flux_BC_Z(char *ID, double *disteven, double *distodd, double flux,
extern "C" double ScaLBL_D3Q19_Flux_BC_Z(char *ID, double *disteven, double *distodd, double Q,
int Nx, int Ny, int Nz, int outlet){
// Note that this routine assumes the distributions are stored "opposite"
// odd distributions in disteven and even distributions in distodd.
@ -347,7 +346,6 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_Z(char *ID, double *disteven, double *dis
N = Nx*Ny*Nz;
// Loop over the boundary - threadblocks delineated by start...finish
double A = 1.f*double((Nx-2)*(Ny-2));
double sum = 0.f;
char id;
for (n=outlet; n<N-Nx*Ny; n++){
@ -380,7 +378,7 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_Z(char *ID, double *disteven, double *dis
sum += (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f5+f11+f14+f15+f18));
}
}
dout = sum/(A*(1.0+flux));
dout = sum/(1.0+Q);
return dout;
}

15
gpu/D3Q19.cu
View File

@ -547,7 +547,7 @@ __global__ void dvc_D3Q19_Velocity_BC_Z(double *disteven, double *distodd, doubl
}
}
__global__ void dvc_D3Q19_Flux_BC_z(double *disteven, double *distodd, double flux, double *dvcsum,
__global__ void dvc_D3Q19_Flux_BC_z(char *ID, double *disteven, double *distodd, double flux, double *dvcsum,
int Nx, int Ny, int Nz){
// Note that this routine assumes the distributions are stored "opposite"
// odd distributions in disteven and even distributions in distodd.
@ -557,7 +557,7 @@ __global__ void dvc_D3Q19_Flux_BC_z(double *disteven, double *distodd, double fl
double f10,f12,f13,f16,f17;
//double A = 1.f*double(Nx*Ny);
double factor = 1.f/(double((Nx-2)*(Ny-2))*(1.0-flux));
double factor = 1.f/((1.0-flux));
double sum = 0.f;
@ -565,7 +565,8 @@ __global__ void dvc_D3Q19_Flux_BC_z(double *disteven, double *distodd, double fl
n = Nx*Ny + blockIdx.x*blockDim.x + threadIdx.x;
if (n < 2*Nx*Ny){
char id=ID[n];
if (id > 0){
//........................................................................
f2 = distodd[n];
f4 = distodd[N+n];
@ -592,6 +593,7 @@ __global__ void dvc_D3Q19_Flux_BC_z(double *disteven, double *distodd, double fl
//sum = (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f6+f12+f13+f16+f17))/(A*(1.0-flux));
sum = factor*(f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f6+f12+f13+f16+f17));
//localsum[n]=sum;
}
}
//atomicAdd(dvcsum, sum);
@ -607,7 +609,7 @@ __global__ void dvc_D3Q19_Flux_BC_z(double *disteven, double *distodd, double fl
atomicAdd(dvcsum, sum);
}
__global__ void dvc_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double flux, double *dvcsum,
__global__ void dvc_D3Q19_Flux_BC_Z(char *ID, double *disteven, double *distodd, double flux, double *dvcsum,
int Nx, int Ny, int Nz, int outlet){
int n,N;
// distributions
@ -617,11 +619,13 @@ __global__ void dvc_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double fl
N = Nx*Ny*Nz;
n = outlet + blockIdx.x*blockDim.x + threadIdx.x;
double factor = 1.f/(double((Nx-2)*(Ny-2))*(1.0+flux));
double factor = 1.f/(1.0+flux);
double sum = 0.f;
// Loop over the boundary - threadblocks delineated by start...finish
if ( n<N-Nx*Ny ){
char id=ID[n];
if (id>0){
//........................................................................
// Read distributions from "opposite" memory convention
//........................................................................
@ -650,6 +654,7 @@ __global__ void dvc_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double fl
//sum = (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f5+f11+f14+f15+f18))/(A*(1.0+flux));
sum = factor*(f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f5+f11+f14+f15+f18));
//localsum[n]=sum;
}
}
//sum = warpReduceSum(sum);

25
tests/lbpm_color_macro_simulator.cpp
View File

@ -501,17 +501,18 @@ int main(int argc, char **argv)
for ( j=0;j<Ny;j++){
for ( i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
if (Dm.kproc==0 && k==0) id[n]=1;
if (Dm.kproc==0 && k==1) id[n]=1;
if (Dm.kproc==nprocz-1 && k==Nz-2) id[n]=2;
if (Dm.kproc==nprocz-1 && k==Nz-1) id[n]=2;
if (id[n] > 0){
if (Dm.kproc==0 && k==0) id[n]=1;
if (Dm.kproc==0 && k==1) id[n]=1;
if (Dm.kproc==nprocz-1 && k==Nz-2) id[n]=2;
if (Dm.kproc==nprocz-1 && k==Nz-1) id[n]=2;
}
Mask.id[n] = id[n];
}
}
}
}
// Initialize communication structures in averaging domain
for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
Mask.CommInit(comm);
@ -711,10 +712,7 @@ int main(int argc, char **argv)
double Area=double(Dm.nprocx*Dm.nprocy*(Dm.Nx-2)*(Dm.Ny-2));
if (rank==0){
printf("Using flux boundary condition \n");
printf(" flux = %f \n",flux);
printf(" area = %f \n",Area);
printf(" Q = %f \n",flux*Area);
printf(" dsw /dt = %f (expected)\n",flux/(porosity*double(Dm.nprocz*(Dm.Nz-2))));
printf(" Volumetric flux: Q = %f \n",flux);
printf(" outlet pressure: %f \n",dout);
}
@ -723,14 +721,14 @@ int main(int argc, char **argv)
}
double tmpdin=din;
MPI_Allreduce(&tmpdin,&din,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
din=din/(1.0*Dm.nprocx*Dm.nprocy);
if (pBC && Dm.kproc == 0){
if (rank==0) printf("Flux = %.3e, Computed inlet pressure: %f \n",flux,din);
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
}
if (pBC && Dm.kproc == nprocz-1){
@ -738,6 +736,8 @@ int main(int argc, char **argv)
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
@ -955,15 +955,12 @@ int main(int argc, char **argv)
// Set flux boundary condition
if (BoundaryCondition==4){
din=0.f;
double Area=double(Dm.nprocx*Dm.nprocy*(Dm.Nx-2)*(Dm.Ny-2));
if (pBC && Dm.kproc == 0){
din = ScaLBL_D3Q19_Flux_BC_z(ID,f_even,f_odd,flux,Nx,Ny,Nz);
}
double tmpdin=din;
MPI_Allreduce(&tmpdin,&din,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
din=din/(1.0*Dm.nprocx*Dm.nprocy);
if (pBC && Dm.kproc == 0){
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);