From b05181bfb2d2716d6d5a4de2c83338456326d21f Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 14 Nov 2015 17:03:04 -0500 Subject: [PATCH] Including two domains in lbpm_color_simulator: Full domain Dm and Mask which excludes the solid portion (for LB communications) --- tests/lbpm_color_simulator.cpp | 70 ++++++++++++++++++---------------- 1 file changed, 38 insertions(+), 32 deletions(-) diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index 2948e00c..a3d0c379 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -309,9 +309,15 @@ int main(int argc, char **argv) else Restart=false; NULL_USE(pBC); NULL_USE(velBC); + // Full domain Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + for (i=0; i Averages( new TwoPhase(Dm) ); + // Mask that excludes the solid phase + Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); + InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc, rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z, rank_xy, rank_XY, rank_xY, rank_Xy, rank_xz, rank_XZ, rank_xZ, rank_Xz, @@ -415,8 +421,8 @@ int main(int argc, char **argv) int kstart,kfinish; kstart = 1; kfinish = Nz-1; - if (BoundaryCondition > 0 && Dm.kproc==0) kstart = 4; - if (BoundaryCondition > 0 && Dm.kproc==nprocz-1) kfinish = Nz-4; + if (BoundaryCondition > 0 && Mask.kproc==0) kstart = 4; + if (BoundaryCondition > 0 && Mask.kproc==nprocz-1) kfinish = Nz-4; // Compute the pore volume sum_local = 0.0; @@ -436,7 +442,7 @@ int main(int argc, char **argv) if (rank==0) printf("Media porosity = %f \n",porosity); //......................................................... // If external boundary conditions are applied remove solid - if (BoundaryCondition > 0 && Dm.kproc == 0){ + if (BoundaryCondition > 0 && Mask.kproc == 0){ for (k=0; k<3; k++){ for (j=0;j 0 && Dm.kproc == nprocz-1){ + if (BoundaryCondition > 0 && Mask.kproc == nprocz-1){ for (k=Nz-3; k 0){ @@ -479,11 +485,11 @@ int main(int argc, char **argv) for ( j=0;jSDs(i) -= (1.0); // //....................................................................... // Finalize setup for averaging domain - //Averages->SetupCubes(Dm); + //Averages->SetupCubes(Mask); Averages->UpdateSolid(); //....................................................................... @@ -601,12 +607,12 @@ int main(int argc, char **argv) printf("Setting inlet pressure = %f \n", din); printf("Setting outlet pressure = %f \n", dout); } - if (BoundaryCondition==1 && Dm.kproc == 0) { + if (BoundaryCondition==1 && Mask.kproc == 0) { PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz); ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); } - if (BoundaryCondition==1 && Dm.kproc == nprocz-1){ + if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); } @@ -615,13 +621,13 @@ int main(int argc, char **argv) printf("Setting inlet velocity = %f \n", din); printf("Setting outlet velocity = %f \n", dout); } - if (BoundaryCondition==2 && Dm.kproc == 0) { + if (BoundaryCondition==2 && Mask.kproc == 0) { ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz); //ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); SetPhiSlice_z(Phi,1.0,Nx,Ny,Nz,0); } - if (BoundaryCondition==2 && Dm.kproc == nprocz-1){ + if (BoundaryCondition==2 && Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); //ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); SetPhiSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); @@ -637,11 +643,11 @@ int main(int argc, char **argv) din = 1.0+0.5*dp; dout = 1.0-0.5*dp; // set the initial boundary conditions - if (Dm.kproc == 0) { + if (Mask.kproc == 0) { PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz); ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); } - if (Dm.kproc == nprocz-1){ + if (Mask.kproc == nprocz-1){ PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); } @@ -650,15 +656,15 @@ int main(int argc, char **argv) ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz); ComputeVelocityD3Q19(ID,f_even,f_odd,Velocity,Nx,Ny,Nz); - /* if (BoundaryCondition==1 && Dm.kproc == 0){ + /* if (BoundaryCondition==1 && Mask.kproc == 0){ for (n=Nx*Ny; n<2*Nx*Ny; n++){ - if (Dm.id[n]>0 && 3.0*Pressure[n] != din) printf("Inlet pBC error: %f != %f \n",3.0*Pressure[n],din); + if (Mask.id[n]>0 && 3.0*Pressure[n] != din) printf("Inlet pBC error: %f != %f \n",3.0*Pressure[n],din); } } - if (BoundaryCondition==1 && Dm.kproc == nprocz-1){ + if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ for (n=Nx*Ny*(Nz-2); n0 && 3.0*Pressure[n] != dout) printf("Outlet pBC error: %f != %f \n",3.0*Pressure[n],din); + if (Mask.id[n]>0 && 3.0*Pressure[n] != dout) printf("Outlet pBC error: %f != %f \n",3.0*Pressure[n],din); } } */ @@ -709,10 +715,10 @@ int main(int argc, char **argv) ThreadPool tpool(N_threads); // Create the MeshDataStruct - fillHalo fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1); + fillHalo fillData(Mask.Comm,Mask.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1); std::vector meshData(1); meshData[0].meshName = "domain"; - meshData[0].mesh = std::shared_ptr( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) ); + meshData[0].mesh = std::shared_ptr( new IO::DomainMesh(Mask.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) ); std::shared_ptr PhaseVar( new IO::Variable() ); std::shared_ptr PressVar( new IO::Variable() ); std::shared_ptr SignDistVar( new IO::Variable() ); @@ -822,22 +828,22 @@ int main(int argc, char **argv) DeviceBarrier(); // Pressure boundary conditions - if (BoundaryCondition==1 && Dm.kproc == 0) { + if (BoundaryCondition==1 && Mask.kproc == 0) { PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz); ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); } - if (BoundaryCondition==1 && Dm.kproc == nprocz-1){ + if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); } // Velocity boundary conditions - if (BoundaryCondition==2 && Dm.kproc == 0) { + if (BoundaryCondition==2 && Mask.kproc == 0) { ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz); //ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); SetPhiSlice_z(Phi,1.0,Nx,Ny,Nz,0); } - if (BoundaryCondition==2 && Dm.kproc == nprocz-1){ + if (BoundaryCondition==2 && Mask.kproc == nprocz-1){ ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); //ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); SetPhiSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); @@ -849,11 +855,11 @@ int main(int argc, char **argv) din = 1.0+0.5*dp; dout = 1.0-0.5*dp; // set the initial boundary conditions - if (Dm.kproc == 0) { + if (Mask.kproc == 0) { PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz); ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); } - if (Dm.kproc == nprocz-1){ + if (Mask.kproc == nprocz-1){ PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); } @@ -903,7 +909,7 @@ int main(int argc, char **argv) Averages->PrintComponents(timestep); // ************************************************************************ - int NumberComponents_NWP = ComputeGlobalPhaseComponent(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,Averages->PhaseID,1,Averages->Label_NWP); + int NumberComponents_NWP = ComputeGlobalPhaseComponent(Mask.Nx-2,Mask.Ny-2,Mask.Nz-2,Mask.rank_info,Averages->PhaseID,1,Averages->Label_NWP); printf("Number of non-wetting phase components: %i \n ",NumberComponents_NWP); DeviceBarrier(); CopyToHost(Averages->Phase.get(),Phi,N*sizeof(double));