//************************************************************************* // Lattice Boltzmann Simulator for Single Phase Flow in Porous Media // James E. McCLure //************************************************************************* #include #include #include #include "Domain.h" #include "D3Q19.h" #include "D3Q7.h" #include "Array.h" #include "Extras.h" #include using namespace std; extern void GlobalFlipInitD3Q19(double *dist_even, double *dist_odd, int Nx, int Ny, int Nz, int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz) { // Set of Discrete velocities for the D3Q19 Model static int D3Q19[18][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}, {1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1}, {0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}}; int q,i,j,k,n,N; int Cqx,Cqy,Cqz; // Discrete velocity int x,y,z; // Global indices int xn,yn,zn; // Global indices of neighbor int X,Y,Z; // Global size X = Nx*nprocx; Y = Ny*nprocy; Z = Nz*nprocz; N = (Nx+2)*(Ny+2)*(Nz+2); // size of the array including halo for (k=0; k> nprocx; domain >> nprocy; domain >> nprocz; domain >> Nx; domain >> Ny; domain >> Nz; domain >> nspheres; domain >> Lx; domain >> Ly; domain >> Lz; //....................................................................... // ************************************************************** } // ************************************************************** // Broadcast simulation parameters from rank 0 to all other procs MPI_Barrier(MPI_COMM_WORLD); //................................................. MPI_Bcast(&Nx,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&Ny,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&Nz,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&nBlocks,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&nthreads,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(×tepMax,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&Nx,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&Ny,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&Nz,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&nprocx,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&nprocy,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&nprocz,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&nspheres,1,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&Lx,1,MPI_DOUBLE,0,MPI_COMM_WORLD); MPI_Bcast(&Ly,1,MPI_DOUBLE,0,MPI_COMM_WORLD); MPI_Bcast(&Lz,1,MPI_DOUBLE,0,MPI_COMM_WORLD); //................................................. MPI_Barrier(MPI_COMM_WORLD); // ************************************************************** // ************************************************************** if (nprocs != nprocx*nprocy*nprocz){ printf("Fatal error in processor number! \n"); printf("nprocx = %i \n",nprocx); printf("nprocy = %i \n",nprocy); printf("nprocz = %i \n",nprocz); abort(); } if (rank==0){ printf("********************************************************\n"); printf("Sub-domain size = %i x %i x %i\n",Nz,Nz,Nz); printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz); printf("********************************************************\n"); } MPI_Barrier(MPI_COMM_WORLD); kproc = rank/(nprocx*nprocy); jproc = (rank-nprocx*nprocy*kproc)/nprocx; iproc = rank-nprocx*nprocy*kproc-nprocz*jproc; //.......................................... // set up the neighbor ranks //.......................................... i=iproc; j=jproc; k =kproc; i+=1; j+=0; k+=0; if (i<0) i+=nprocx; if (j<0) j+=nprocy; if (k<0) k+=nprocz; if (!(i 0.0){ id[n] = 1; sum++; } } } } sum_local = 1.0*sum; MPI_Allreduce(&sum_local,&porosity,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); porosity = porosity*iVol_global; if (rank==0) printf("Media porosity = %f \n",porosity); */ //........................................................................... MPI_Barrier(MPI_COMM_WORLD); if (rank == 0) cout << "Domain set." << endl; //........................................................................... //........................................................................... // Set up MPI communication structurese if (rank==0) printf ("Setting up communication control structures \n"); //...................................................................................... // Get the actual D3Q19 communication counts (based on location of solid phase) // Discrete velocity set symmetry implies the sendcount = recvcount int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z; int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ; int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ; sendCount_x = sendCount_y = sendCount_z = sendCount_X = sendCount_Y = sendCount_Z = 0; sendCount_xy = sendCount_yz = sendCount_xz = sendCount_Xy = sendCount_Yz = sendCount_xZ = 0; sendCount_xY = sendCount_yZ = sendCount_Xz = sendCount_XY = sendCount_YZ = sendCount_XZ = 0; //...................................................................................... for (k=0; k