From becd898547fa0ef3ce3fd6e38c1283e7cd0bda4c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 20 Nov 2017 22:19:45 -0500 Subject: [PATCH 1/7] Fixed BC bug --- cpu/D3Q19.cpp | 130 +++++++++++++++++++++++++++----------------------- 1 file changed, 69 insertions(+), 61 deletions(-) diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp index 52473181..a5e6b5fe 100644 --- a/cpu/D3Q19.cpp +++ b/cpu/D3Q19.cpp @@ -293,7 +293,7 @@ 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*Ny); + double A = 1.f*double((Nx-2)*(Ny-2)); double sum = 0.f; char id; for (n=Nx*Ny; n<2*Nx*Ny; n++){ @@ -303,26 +303,26 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID, double *disteven, double *di // Read distributions from "opposite" memory convention //........................................................................ //........................................................................ - f1 = distodd[n]; - f3 = distodd[N+n]; - f5 = distodd[2*N+n]; - f7 = distodd[3*N+n]; - f9 = distodd[4*N+n]; - f11 = distodd[5*N+n]; - f13 = distodd[6*N+n]; - f15 = distodd[7*N+n]; - f17 = distodd[8*N+n]; + f2 = distodd[n]; + f4 = distodd[N+n]; + f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; //........................................................................ f0 = disteven[n]; - f2 = disteven[N+n]; - f4 = disteven[2*N+n]; - f6 = disteven[3*N+n]; - f8 = disteven[4*N+n]; - f10 = disteven[5*N+n]; - f12 = disteven[6*N+n]; - f14 = disteven[7*N+n]; - f16 = disteven[8*N+n]; - f18 = disteven[9*N+n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + f17 = disteven[9*N+n]; //................................................... // Determine the outlet flow velocity @@ -402,26 +402,26 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distodd, do // Read distributions from "opposite" memory convention //........................................................................ //........................................................................ - f1 = distodd[n]; - f3 = distodd[N+n]; - f5 = distodd[2*N+n]; - f7 = distodd[3*N+n]; - f9 = distodd[4*N+n]; - f11 = distodd[5*N+n]; - f13 = distodd[6*N+n]; - f15 = distodd[7*N+n]; - f17 = distodd[8*N+n]; + f2 = distodd[n]; + f4 = distodd[N+n]; + f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; //........................................................................ f0 = disteven[n]; - f2 = disteven[N+n]; - f4 = disteven[2*N+n]; - f6 = disteven[3*N+n]; - f8 = disteven[4*N+n]; - f10 = disteven[5*N+n]; - f12 = disteven[6*N+n]; - f14 = disteven[7*N+n]; - f16 = disteven[8*N+n]; - f18 = disteven[9*N+n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + f17 = disteven[9*N+n]; //................................................... //........Determine the inlet flow velocity......... // uz = -1 + (f0+f3+f4+f1+f2+f7+f8+f10+f9 @@ -447,11 +447,19 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distodd, do f15 = f16 + 0.16666666666666678*(uy+uz)-Cyz; f18 = f17 + 0.16666666666666678*(uz-uy)+Cyz; //........Store in "opposite" memory location.......... - distodd[2*N+n] = f5; + /* distodd[2*N+n] = f5; distodd[5*N+n] = f11; disteven[7*N+n] = f14; distodd[7*N+n] = f15; disteven[9*N+n] = f18; + + */ + + disteven[3*N+n] = f5; + disteven[6*N+n] = f11; + distodd[6*N+n] = f14; + disteven[8*N+n] = f15; + distodd[8*N+n] = f18; /* printf("Site=%i\n",n); printf("ux=%f, uy=%f, uz=%f\n",ux,uy,uz); @@ -480,26 +488,26 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distodd, do //........................................................................ // Read distributions from "opposite" memory convention //........................................................................ - f1 = distodd[n]; - f3 = distodd[N+n]; - f5 = distodd[2*N+n]; - f7 = distodd[3*N+n]; - f9 = distodd[4*N+n]; - f11 = distodd[5*N+n]; - f13 = distodd[6*N+n]; - f15 = distodd[7*N+n]; - f17 = distodd[8*N+n]; + f2 = distodd[n]; + f4 = distodd[N+n]; + f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; //........................................................................ f0 = disteven[n]; - f2 = disteven[N+n]; - f4 = disteven[2*N+n]; - f6 = disteven[3*N+n]; - f8 = disteven[4*N+n]; - f10 = disteven[5*N+n]; - f12 = disteven[6*N+n]; - f14 = disteven[7*N+n]; - f16 = disteven[8*N+n]; - f18 = disteven[9*N+n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + f17 = disteven[9*N+n]; //........Determine the outlet flow velocity......... // uz = 1 - (f0+f3+f4+f1+f2+f7+f8+f10+f9+ // 2*(f6+f16+f17+f12+f13))/dout; @@ -531,11 +539,11 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distodd, do f17 = f18 - 0.16666666666666678*(uz-uy)-Cyz; //........Store in "opposite" memory location.......... - disteven[3*N+n] = f6; - disteven[6*N+n] = f12; - distodd[6*N+n] = f13; - disteven[8*N+n] = f16; - distodd[8*N+n] = f17; + distodd[2*N+n] = f6; + distodd[5*N+n] = f12; + disteven[7*N+n] = f13; + distodd[7*N+n] = f16; + disteven[9*N+n] = f17; //................................................... } From 4261373574c5ee2115e85cc1ed44bc23bcd88f8a Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Tue, 21 Nov 2017 14:53:45 +1100 Subject: [PATCH 2/7] fix a little bug in SetSlice for the 'Phi' values to be consistent with the corresponding bounary conditions --- tests/lbpm_color_macro_simulator.cpp | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/tests/lbpm_color_macro_simulator.cpp b/tests/lbpm_color_macro_simulator.cpp index 03ba424b..5a1ef5ae 100644 --- a/tests/lbpm_color_macro_simulator.cpp +++ b/tests/lbpm_color_macro_simulator.cpp @@ -648,13 +648,13 @@ int main(int argc, char **argv) if (BoundaryCondition==1 && Mask.kproc == 0) { ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,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,0); } if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,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-1); } if (rank==0 && BoundaryCondition==2){ @@ -664,13 +664,13 @@ int main(int argc, char **argv) if (BoundaryCondition==2 && Mask.kproc == 0) { ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,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,0); } if (BoundaryCondition==2 && Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,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-1); } // Set dynamic pressure boundary conditions @@ -687,13 +687,13 @@ int main(int argc, char **argv) if (Mask.kproc == 0) { ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,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,0); } if (Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,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-1); } } @@ -721,7 +721,7 @@ int main(int argc, char **argv) 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,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,0); } @@ -729,7 +729,7 @@ int main(int argc, char **argv) // if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,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-1); } } @@ -900,22 +900,24 @@ int main(int argc, char **argv) if (BoundaryCondition==1 && Mask.kproc == 0) { ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); } if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); } // Velocity boundary conditions if (BoundaryCondition==2 && Mask.kproc == 0) { ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,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,0); } if (BoundaryCondition==2 && Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,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-1); } if (BoundaryCondition==3){ @@ -927,10 +929,12 @@ int main(int argc, char **argv) if (Mask.kproc == 0) { ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); } if (Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); } } @@ -950,14 +954,14 @@ int main(int argc, char **argv) if (pBC && Dm.kproc == 0){ ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,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,0); } if (pBC && Dm.kproc == nprocz-1){ // if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,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-1); } } From 1d645183c0e206dcbda39b0abc1e9a38b298066e Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Tue, 21 Nov 2017 14:58:51 +1100 Subject: [PATCH 3/7] lbpm_color_simulator: fix a little bug in SetSlice for the 'Phi' values to be consistent with the corresponding bounary conditions --- tests/lbpm_color_simulator.cpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index c808e668..73181948 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -638,11 +638,13 @@ int main(int argc, char **argv) if (BoundaryCondition==1 && Mask.kproc == 0) { ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); } if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); } if (rank==0 && BoundaryCondition==2){ @@ -651,16 +653,14 @@ int main(int argc, char **argv) } if (BoundaryCondition==2 && Mask.kproc == 0) { ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz); - //ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,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,0); } if (BoundaryCondition==2 && Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); - //ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,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-1); } // Set dynamic pressure boundary conditions double dp, slope; @@ -676,10 +676,12 @@ int main(int argc, char **argv) if (Mask.kproc == 0) { ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); } if (Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); } } @@ -848,10 +850,12 @@ int main(int argc, char **argv) if (BoundaryCondition==1 && Mask.kproc == 0) { ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); } if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); } // Velocity boundary conditions @@ -875,10 +879,12 @@ int main(int argc, char **argv) if (Mask.kproc == 0) { ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); } if (Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); + ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); } } From 250551b2c7231e27dfecfe191bcaf59939434b2f Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Tue, 21 Nov 2017 15:18:27 +1100 Subject: [PATCH 4/7] D3Q19.cu: fix the bug for read/write associated with the opposite memory location of distribution functions --- gpu/D3Q19.cu | 172 +++++++++++++++++++++++++-------------------------- 1 file changed, 86 insertions(+), 86 deletions(-) diff --git a/gpu/D3Q19.cu b/gpu/D3Q19.cu index b125a2ca..b137eeef 100644 --- a/gpu/D3Q19.cu +++ b/gpu/D3Q19.cu @@ -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*Ny)*(1.0-flux)); + double factor = 1.f/(double((Nx-2)*(Ny-2))*(1.0-flux)); double sum = 0.f; @@ -567,26 +567,26 @@ __global__ void dvc_D3Q19_Flux_BC_z(double *disteven, double *distodd, double fl if (n < 2*Nx*Ny){ //........................................................................ - f1 = distodd[n]; - f3 = distodd[N+n]; -// f5 = distodd[2*N+n]; - f7 = distodd[3*N+n]; - f9 = distodd[4*N+n]; -// f11 = distodd[5*N+n]; - f13 = distodd[6*N+n]; -// f15 = distodd[7*N+n]; - f17 = distodd[8*N+n]; + f2 = distodd[n]; + f4 = distodd[N+n]; + f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + f12 = distodd[5*N+n]; + //f14 = distodd[6*N+n]; + f16 = distodd[7*N+n]; + //f18 = distodd[8*N+n]; //........................................................................ f0 = disteven[n]; - f2 = disteven[N+n]; - f4 = disteven[2*N+n]; - f6 = disteven[3*N+n]; - f8 = disteven[4*N+n]; - f10 = disteven[5*N+n]; - f12 = disteven[6*N+n]; -// f14 = disteven[7*N+n]; - f16 = disteven[8*N+n]; -// f18 = disteven[9*N+n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + //f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + //f11 = disteven[6*N+n]; + f13 = disteven[7*N+n]; + //f15 = disteven[8*N+n]; + f17 = disteven[9*N+n]; //................................................... // compute local sum to determine the density value to set pressure //sum = (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f6+f12+f13+f16+f17))/(A*(1.0-flux)); @@ -617,7 +617,7 @@ __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*Ny)*(1.0+flux)); + double factor = 1.f/(double((Nx-2)*(Ny-2))*(1.0+flux)); double sum = 0.f; // Loop over the boundary - threadblocks delineated by start...finish @@ -625,26 +625,26 @@ __global__ void dvc_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double fl //........................................................................ // Read distributions from "opposite" memory convention //........................................................................ - f1 = distodd[n]; - f3 = distodd[N+n]; - f5 = distodd[2*N+n]; - f7 = distodd[3*N+n]; - f9 = distodd[4*N+n]; - f11 = distodd[5*N+n]; -// f13 = distodd[6*N+n]; - f15 = distodd[7*N+n]; -// f17 = distodd[8*N+n]; + f2 = distodd[n]; + f4 = distodd[N+n]; + //f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + //f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + //f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; //........................................................................ f0 = disteven[n]; - f2 = disteven[N+n]; - f4 = disteven[2*N+n]; -// f6 = disteven[3*N+n]; - f8 = disteven[4*N+n]; - f10 = disteven[5*N+n]; -// f12 = disteven[6*N+n]; - f14 = disteven[7*N+n]; -// f16 = disteven[8*N+n]; - f18 = disteven[9*N+n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + //f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + //f17 = disteven[9*N+n]; // Local sum (based on the consistency condition) //sum = (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f5+f11+f14+f15+f18))/(A*(1.0+flux)); @@ -682,26 +682,26 @@ __global__ void dvc_ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distod // Read distributions from "opposite" memory convention //........................................................................ //........................................................................ - f1 = distodd[n]; - f3 = distodd[N+n]; - f5 = distodd[2*N+n]; - f7 = distodd[3*N+n]; - f9 = distodd[4*N+n]; - f11 = distodd[5*N+n]; - f13 = distodd[6*N+n]; - f15 = distodd[7*N+n]; - f17 = distodd[8*N+n]; - //........................................................................ - f0 = disteven[n]; - f2 = disteven[N+n]; - f4 = disteven[2*N+n]; - f6 = disteven[3*N+n]; - f8 = disteven[4*N+n]; - f10 = disteven[5*N+n]; - f12 = disteven[6*N+n]; - f14 = disteven[7*N+n]; - f16 = disteven[8*N+n]; - f18 = disteven[9*N+n]; + f2 = distodd[n]; + f4 = distodd[N+n]; + f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; + //........................................................................ + f0 = disteven[n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + f17 = disteven[9*N+n]; //................................................... //........Determine the inlet flow velocity......... // uz = -1 + (f0+f3+f4+f1+f2+f7+f8+f10+f9 @@ -727,11 +727,11 @@ __global__ void dvc_ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distod f15 = f16 + 0.16666666666666678*(uy+uz)-Cyz; f18 = f17 + 0.16666666666666678*(uz-uy)+Cyz; //........Store in "opposite" memory location.......... - distodd[2*N+n] = f5; - distodd[5*N+n] = f11; - disteven[7*N+n] = f14; - distodd[7*N+n] = f15; - disteven[9*N+n] = f18; + disteven[3*N+n] = f5; + disteven[6*N+n] = f11; + distodd[6*N+n] = f14; + disteven[8*N+n] = f15; + distodd[8*N+n] = f18; } } @@ -753,26 +753,26 @@ __global__ void dvc_ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distod //........................................................................ // Read distributions from "opposite" memory convention //........................................................................ - f1 = distodd[n]; - f3 = distodd[N+n]; - f5 = distodd[2*N+n]; - f7 = distodd[3*N+n]; - f9 = distodd[4*N+n]; - f11 = distodd[5*N+n]; - f13 = distodd[6*N+n]; - f15 = distodd[7*N+n]; - f17 = distodd[8*N+n]; + f2 = distodd[n]; + f4 = distodd[N+n]; + f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; //........................................................................ f0 = disteven[n]; - f2 = disteven[N+n]; - f4 = disteven[2*N+n]; - f6 = disteven[3*N+n]; - f8 = disteven[4*N+n]; - f10 = disteven[5*N+n]; - f12 = disteven[6*N+n]; - f14 = disteven[7*N+n]; - f16 = disteven[8*N+n]; - f18 = disteven[9*N+n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + f17 = disteven[9*N+n]; //........Determine the outlet flow velocity......... // uz = 1 - (f0+f3+f4+f1+f2+f7+f8+f10+f9+ // 2*(f6+f16+f17+f12+f13))/dout; @@ -804,11 +804,11 @@ __global__ void dvc_ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distod f17 = f18 - 0.16666666666666678*(uz-uy)-Cyz; //........Store in "opposite" memory location.......... - disteven[3*N+n] = f6; - disteven[6*N+n] = f12; - distodd[6*N+n] = f13; - disteven[8*N+n] = f16; - distodd[8*N+n] = f17; + distodd[2*N+n] = f6; + distodd[5*N+n] = f12; + disteven[7*N+n] = f13; + distodd[7*N+n] = f16; + disteven[9*N+n] = f17; //................................................... } } From 9c7605403ff02ba0b95148c00c167089347bbda0 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Tue, 21 Nov 2017 16:11:52 +1100 Subject: [PATCH 5/7] D3Q19.cpp: for flux BC comment out the unncessary read of the distribution functions --- common/ScaLBL.h | 2 +- cpu/D3Q19.cpp | 69 ++++++++++++++++++++++---------------------- tests/TestFluxBC.cpp | 2 +- 3 files changed, 37 insertions(+), 36 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 06afeb4d..a8dd8dd8 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -62,7 +62,7 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distodd, do extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID, double *disteven, double *distodd, double flux, int Nx, int Ny, int Nz); -extern "C" double ScaLBL_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double flux, +extern "C" double ScaLBL_D3Q19_Flux_BC_Z(char *ID, double *disteven, double *distodd, double flux, int Nx, int Ny, int Nz, int outlet); extern "C" void ScaLBL_Color_Init(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz); diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp index a5e6b5fe..5d08cc8f 100644 --- a/cpu/D3Q19.cpp +++ b/cpu/D3Q19.cpp @@ -309,19 +309,19 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID, double *disteven, double *di f8 = distodd[3*N+n]; f10 = distodd[4*N+n]; f12 = distodd[5*N+n]; - f14 = distodd[6*N+n]; + //f14 = distodd[6*N+n]; f16 = distodd[7*N+n]; - f18 = distodd[8*N+n]; + //f18 = distodd[8*N+n]; //........................................................................ f0 = disteven[n]; f1 = disteven[N+n]; f3 = disteven[2*N+n]; - f5 = disteven[3*N+n]; + //f5 = disteven[3*N+n]; f7 = disteven[4*N+n]; f9 = disteven[5*N+n]; - f11 = disteven[6*N+n]; + //f11 = disteven[6*N+n]; f13 = disteven[7*N+n]; - f15 = disteven[8*N+n]; + //f15 = disteven[8*N+n]; f17 = disteven[9*N+n]; //................................................... @@ -335,7 +335,7 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID, double *disteven, double *di return din; } -extern "C" double ScaLBL_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double flux, +extern "C" double ScaLBL_D3Q19_Flux_BC_Z(char *ID, double *disteven, double *distodd, double flux, 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. @@ -344,40 +344,41 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_Z(double *disteven, double *distodd, doub double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9; double f10,f11,f12,f13,f14,f15,f16,f17,f18; double dout = 0.f; - N = Nx*Ny*Nz; // Loop over the boundary - threadblocks delineated by start...finish - double A = 1.f*double(Nx*Ny); + double A = 1.f*double((Nx-2)*(Ny-2)); double sum = 0.f; + char id; for (n=outlet; n0){ + //........................................................................ + // Read distributions from "opposite" memory convention + //........................................................................ + f2 = distodd[n]; + f4 = distodd[N+n]; + //f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + //f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + //f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; + //........................................................................ + f0 = disteven[n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + //f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + //f17 = disteven[9*N+n]; - //........................................................................ - // Read distributions from "opposite" memory convention - //........................................................................ - f1 = distodd[n]; - f3 = distodd[N+n]; - f5 = distodd[2*N+n]; - f7 = distodd[3*N+n]; - f9 = distodd[4*N+n]; - f11 = distodd[5*N+n]; - f13 = distodd[6*N+n]; - f15 = distodd[7*N+n]; - f17 = distodd[8*N+n]; - //........................................................................ - f0 = disteven[n]; - f2 = disteven[N+n]; - f4 = disteven[2*N+n]; - f6 = disteven[3*N+n]; - f8 = disteven[4*N+n]; - f10 = disteven[5*N+n]; - f12 = disteven[6*N+n]; - f14 = disteven[7*N+n]; - f16 = disteven[8*N+n]; - f18 = disteven[9*N+n]; - - sum += (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f5+f11+f14+f15+f18)); - + sum += (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f5+f11+f14+f15+f18)); + } } dout = sum/(A*(1.0+flux)); return dout; diff --git a/tests/TestFluxBC.cpp b/tests/TestFluxBC.cpp index ac5e49d7..af5b4360 100644 --- a/tests/TestFluxBC.cpp +++ b/tests/TestFluxBC.cpp @@ -105,7 +105,7 @@ int main (int argc, char **argv) } if (pBC && Dm.kproc == nprocz-1){ - dout = ScaLBL_D3Q19_Flux_BC_Z(f_even,f_odd,flux,Nx,Ny,Nz,Nx*Ny*(Nz-2)); + dout = ScaLBL_D3Q19_Flux_BC_Z(ID,f_even,f_odd,flux,Nx,Ny,Nz,Nx*Ny*(Nz-2)); printf("Computed outlet pressure: %f \n", dout); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); } From 0efeda94cc19272d39f48e4db2f5e38f4192d564 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Thu, 23 Nov 2017 16:07:24 +1100 Subject: [PATCH 6/7] fill zero velocity values for the solid nodes --- cpu/D3Q19.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp index 5d08cc8f..9428913f 100644 --- a/cpu/D3Q19.cpp +++ b/cpu/D3Q19.cpp @@ -714,6 +714,11 @@ extern "C" void ScaLBL_D3Q19_Velocity(char *ID, double *disteven, double *distod distodd[q*N+n] = -1.0; } disteven[9*N+n] = -1.0; + + //For ID[n]<0 - solid nodes + vel[n] = 0.0; + vel[N+n] = 0.0; + vel[2*N+n] = 0.0; } } } From c14fa7cbbfb572a9b0e403944f425e5c62bb6d98 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Thu, 23 Nov 2017 17:12:38 +1100 Subject: [PATCH 7/7] TestFluxBC.cpp: update the test and looks like our flux BC works ! --- tests/TestFluxBC.cpp | 31 ++++++++++++++++++++++++------- 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/tests/TestFluxBC.cpp b/tests/TestFluxBC.cpp index af5b4360..73a75804 100644 --- a/tests/TestFluxBC.cpp +++ b/tests/TestFluxBC.cpp @@ -70,6 +70,15 @@ int main (int argc, char **argv) if (rank==0) printf ("Copying phase ID to device \n"); char *ID; ScaLBL_AllocateDeviceMemory((void **) &ID, N); // Allocate device memory + // Don't compute in the halo + for (k=0;k 1e-14){ @@ -145,10 +158,14 @@ int main (int argc, char **argv) for (j=0; j