diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp new file mode 100644 index 00000000..74eae2ff --- /dev/null +++ b/models/MRTModel.cpp @@ -0,0 +1,181 @@ +/* + * Multi-relaxation time LBM Model + */ +#include "models/MRTModel.h" + +ScaLBL_MRTModel::ScaLBL_MRTModel(int RANK, int NP, MPI_Comm 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), +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) +{ + +} +ScaLBL_MRTModel::~ScaLBL_MRTModel(){ + +} + +void ScaLBL_MRTModel::ReadParams(string filename){ + // read the input database + db = std::make_shared( filename ); + domain_db = db->getDatabase( "Domain" ); + mrt_db = db->getDatabase( "MRT" ); + + // Color Model parameters + timestepMax = mrt_db->getScalar( "timestepMax" ); + tauA = mrt_db->getScalar( "tau" ); + Fx = mrt_db->getVector( "F" )[0]; + Fy = mrt_db->getVector( "F" )[1]; + Fz = mrt_db->getVector( "F" )[2]; + Restart = mrt_db->getScalar( "Restart" ); + din = mrt_db->getScalar( "din" ); + dout = mrt_db->getScalar( "dout" ); + flux = mrt_db->getScalar( "flux" ); + + // Read domain parameters + auto L = domain_db->getVector( "L" ); + auto size = domain_db->getVector( "n" ); + auto nproc = domain_db->getVector( "nproc" ); + BoundaryCondition = domain_db->getScalar( "BC" ); + Nx = size[0]; + Ny = size[1]; + Nz = size[2]; + Lx = L[0]; + Ly = L[1]; + Lz = L[2]; + nprocx = nproc[0]; + nprocy = nproc[1]; + nprocz = nproc[2]; + mu=(tau-0.5)/3.0; +} +void ScaLBL_MRTModel::SetDomain(){ + Dm = std::shared_ptr(new Domain(domain_db,comm)); // full domain for analysis + Mask = std::shared_ptr(new Domain(domain_db,comm)); // mask domain removes immobile phases + Nx+=2; Ny+=2; Nz += 2; + N = Nx*Ny*Nz; + for (int i=0; iid[i] = 1; // initialize this way + //Averages = std::shared_ptr ( new TwoPhase(Dm) ); // TwoPhase analysis object + MPI_Barrier(comm); + Dm->CommInit(); + MPI_Barrier(comm); +} + +void ScaLBL_MRTModel::ReadInput(){ + int rank=Dm->rank(); + size_t readID; + //....................................................................... + if (rank == 0) printf("Read input media... \n"); + //....................................................................... + Mask->ReadIDs(); + + sprintf(LocalRankString,"%05d",Dm->rank()); + sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); + sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString); + + // .......... READ THE INPUT FILE ....................................... + //........................................................................... + if (rank == 0) cout << "Reading in signed distance function..." << endl; + //....................................................................... + sprintf(LocalRankString,"%05d",rank); + sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString); + ReadBinaryFile(LocalRankFilename, Averages->SDs.data(), N); + MPI_Barrier(comm); + if (rank == 0) cout << "Domain set." << endl; +} + +void ScaLBL_MRTModel::Create(){ + /* + * This function creates the variables needed to run a LBM + */ + int rank=Mask->rank(); + //......................................................... + // Initialize communication structures in averaging domain + for (int i=0; iid[i] = Mask->id[i]; + Mask->CommInit(); + Np=Mask->PoreCount(); + //........................................................................... + if (rank==0) printf ("Create ScaLBL_Communicator \n"); + // Create a communicator for the device (will use optimized layout) + // ScaLBL_Communicator ScaLBL_Comm(Mask); // original + ScaLBL_Comm = std::shared_ptr(new ScaLBL_Communicator(Mask)); + + int Npad=(Np/16 + 2)*16; + if (rank==0) printf ("Set up memory efficient layout \n"); + Map.resize(Nx,Ny,Nz); Map.fill(-2); + auto neighborList= new int[18*Npad]; + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np); + MPI_Barrier(comm); + //........................................................................... + // MAIN VARIABLES ALLOCATED HERE + //........................................................................... + // LBM variables + if (rank==0) printf ("Allocating distributions \n"); + //......................device distributions................................. + int dist_mem_size = Np*sizeof(double); + int neighborSize=18*(Np*sizeof(int)); + //........................................................................... + ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); + ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); + //........................................................................... + // Update GPU data structures + if (rank==0) printf ("Setting up device map and neighbor list \n"); + // copy the neighbor list + ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); + MPI_Barrier(comm); +} + +void ScaLBL_MRTModel::Initialize(){ + /* + * This function initializes model + */ + if (rank==0) printf ("Initializing distributions \n"); + ScaLBL_D3Q19_Init(fq, Np); +} + +void ScaLBL_MRTModel::Run(){ + double rlx_setA=1.0/tau; + double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); + //.......create and start timer............ + double starttime,stoptime,cputime; + ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + starttime = MPI_Wtime(); + if (rank==0) printf("Beginning AA timesteps...\n"); + if (rank==0) printf("********************************************************\n"); + timestep=0; + while (timestep < timestepMax) { + //************************************************************************/ + timestep++; + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + timestep++; + ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL + ScaLBL_D3Q19_AAeven_MRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_D3Q19_AAeven_MRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + //************************************************************************/ + } + //************************************************************************/ + stoptime = MPI_Wtime(); + // cout << "CPU time: " << (stoptime - starttime) << " seconds" << endl; + cputime = stoptime - starttime; + // cout << "Lattice update rate: "<< double(Nx*Ny*Nz*timestep)/cputime/1000000 << " MLUPS" << endl; + double MLUPS = double(Np*timestep)/cputime/1000000; +// if (rank==0) printf("********************************************************\n"); +// if (rank==0) printf("CPU time = %f \n", cputime); +// if (rank==0) printf("Lattice update rate (per process)= %f MLUPS \n", MLUPS); + MLUPS *= nprocs; +} + +void ScaLBL_MRTModel::VelocityField(double *Vz){ + + int SIZE=Np*sizeof(double); + ScaLBL_D3Q19_Momentum(fq,Velocity, Np); + ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + ScaLBL_CopyToHost(&Vz[0],&Velocity[2*Np],SIZE); + +} diff --git a/models/MRTModel.h b/models/MRTModel.h new file mode 100644 index 00000000..17e0894c --- /dev/null +++ b/models/MRTModel.h @@ -0,0 +1,69 @@ +/* + * Multi-relaxation time LBM Model + */ +#include +#include +#include +#include +#include +#include +#include + +#include "common/Communication.h" +#include "common/MPI_Helpers.h" +#include "ProfilerApp.h" + +class ScaLBL_MRTModel{ +public: + ScaLBL_MRTModel(int RANK, int NP, MPI_Comm COMM); + ~ScaLBL_MRTModel(); + + // functions in they should be run + void ReadParams(string filename); + void ReadParams(std::shared_ptr db0); + void SetDomain(); + void ReadInput(); + void Create(); + void Initialize(); + void Run(); + void VelocityField(double *Vz); + + bool Restart,pBC; + int timestep,timestepMax; + int BoundaryCondition; + double tau,mu; + double Fx,Fy,Fz,flux; + double din,dout; + + int Nx,Ny,Nz,N,Np; + int rank,nprocx,nprocy,nprocz,nprocs; + double Lx,Ly,Lz; + + std::shared_ptr Dm; // this domain is for analysis + std::shared_ptr Mask; // this domain is for lbm + std::shared_ptr ScaLBL_Comm; + std::shared_ptr Averages; + + // input database + std::shared_ptr db; + std::shared_ptr domain_db; + std::shared_ptr color_db; + std::shared_ptr analysis_db; + + IntArray Map; + int *NeighborList; + double *fq; + double *Velocity; + double *Pressure; + +private: + MPI_Comm comm; + + // filenames + char LocalRankString[8]; + char LocalRankFilename[40]; + char LocalRestartFile[40]; + + //int rank,nprocs; + void LoadParams(std::shared_ptr db0); +}; \ No newline at end of file diff --git a/tests/TestPoiseuille.cpp b/tests/TestPoiseuille.cpp index 73909c93..ba23420d 100644 --- a/tests/TestPoiseuille.cpp +++ b/tests/TestPoiseuille.cpp @@ -8,6 +8,7 @@ #include #include "common/ScaLBL.h" #include "common/MPI_Helpers.h" +#include "models/MRTModel.h" //*************************************************************************************** int main(int argc, char **argv) @@ -28,226 +29,35 @@ int main(int argc, char **argv) printf("Running Unit Test: TestPoiseuille \n"); printf("********************************************************\n"); } + + ScaLBL_MRTModel MRT(rank,nprocs,comm); - // BGK Model parameters - double tau,Fx,Fy,Fz; - // Domain variables - int i,j,k,n; - int timestep = 0; - - tau = 1.0; - double mu=(tau-0.5)/3.0; - double rlx_setA=1.0/tau; - double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); - Fx = 0; Fy = 0; - Fz = 1e-3; //1.f; // 1e-3; - // Load inputs - if (rank==0) printf("Loading input database \n"); - auto FILENAME = argv[1]; - auto db = std::make_shared( FILENAME ); - auto domain_db = db->getDatabase( "Domain" ); - int Nx = domain_db->getVector( "n" )[0]; - int Ny = domain_db->getVector( "n" )[1]; - int Nz = domain_db->getVector( "n" )[2]; - int nprocx = domain_db->getVector( "nproc" )[0]; - int nprocy = domain_db->getVector( "nproc" )[1]; - int nprocz = domain_db->getVector( "nproc" )[2]; - if (rank==0){ - printf("********************************************************\n"); - printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz); - printf("********************************************************\n"); - } - - MPI_Barrier(comm); - int kproc = rank/(nprocx*nprocy); - int jproc = (rank-nprocx*nprocy*kproc)/nprocx; - int iproc = rank-nprocx*nprocy*kproc-nprocz*jproc; - - if (rank == 0) { - printf("i,j,k proc=%d %d %d \n",iproc,jproc,kproc); - } - MPI_Barrier(comm); - if (rank == 1){ - printf("i,j,k proc=%d %d %d \n",iproc,jproc,kproc); - printf("\n\n"); - } - - double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz; - - std::shared_ptr Dm( new Domain(domain_db,comm)); - - Nx += 2; - Ny += 2; - Nz += 2; - int N = Nx*Ny*Nz; - - //....................................................................... - // Assign the phase ID field - //....................................................................... - char LocalRankString[8]; - sprintf(LocalRankString,"%05d",rank); - char LocalRankFilename[40]; - sprintf(LocalRankFilename,"ID.%05i",rank); - /* - FILE *IDFILE = fopen(LocalRankFilename,"rb"); - if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx"); - fread(Dm->id,1,N,IDFILE); - fclose(IDFILE); - */ - - // initialize empty domain - for (k=0;kid[n] = 0; - else if (i>Nx-3) Dm->id[n] = 0; - else Dm->id[n]=1; - } - } - } - - Dm->CommInit(); - MPI_Barrier(comm); - - //....................................................................... - // Compute the media porosity - //....................................................................... - double sum; - double sum_local=0.0, porosity; - int Np=0; // number of local pore nodes - for (k=1;kid[n] > 0){ - sum_local+=1.0; - Np++; - } - } - } - } - MPI_Allreduce(&sum_local,&sum,1,MPI_DOUBLE,MPI_SUM,comm); - porosity = sum*iVol_global; - if (rank==0) printf("Media porosity = %f \n",porosity); - - MPI_Barrier(comm); - if (rank == 0) cout << "Domain set." << endl; - if (rank==0) printf ("Create ScaLBL_Communicator \n"); - - // Create a communicator for the device - std::shared_ptr ScaLBL_Comm(new ScaLBL_Communicator(Dm)); - - //...........device phase ID................................................. - if (rank==0) printf ("Copying phase ID to device \n"); - char *ID; - ScaLBL_AllocateDeviceMemory((void **) &ID, N); // Allocate device memory - // Copy to the device - ScaLBL_CopyToDevice(ID, Dm->id, N); - //........................................................................... - - if (rank==0){ - printf("Total domain size = %i \n",N); - printf("Reduced domain size = %i \n",Np); - } - - // LBM variables - if (rank==0) printf ("Allocating distributions \n"); - if (rank==0) printf ("Set up the neighborlist \n"); - int Npad=Np+32; - int neighborSize=18*Npad*sizeof(int); - int *neighborList; - IntArray Map(Nx,Ny,Nz); - neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np); - MPI_Barrier(comm); - - //......................device distributions................................. - int dist_mem_size = Np*sizeof(double); - - int *NeighborList; - // double *f_even,*f_odd; - double * dist; - double * Velocity; - //........................................................................... - ScaLBL_AllocateDeviceMemory((void **) &dist, 19*dist_mem_size); - ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); - ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); - ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); - //........................................................................... - - /* - * AA Algorithm begins here - * - */ - ScaLBL_D3Q19_Init(dist, Np); - - //.......create and start timer............ - double starttime,stoptime,cputime; - - ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - starttime = MPI_Wtime(); - - /************ MAIN ITERATION LOOP (timing communications)***************************************/ - // ScaLBL_Comm->SendD3Q19(dist, &dist[10*Np]); - // ScaLBL_Comm->RecvD3Q19(dist, &dist[10*Np]); - // ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - // - - if (rank==0) printf("Beginning AA timesteps...\n"); - if (rank==0) printf("********************************************************\n"); - - while (timestep < 2000) { - - ScaLBL_Comm->SendD3Q19AA(dist); //READ FROM NORMAL - ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz); - ScaLBL_Comm->RecvD3Q19AA(dist); //WRITE INTO OPPOSITE - ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, 0, ScaLBL_Comm->next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz); - ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - timestep++; - - ScaLBL_Comm->SendD3Q19AA(dist); //READ FORM NORMAL - ScaLBL_D3Q19_AAeven_MRT(dist, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz); - ScaLBL_Comm->RecvD3Q19AA(dist); //WRITE INTO OPPOSITE - ScaLBL_D3Q19_AAeven_MRT(dist, 0, ScaLBL_Comm->next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz); - ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - timestep++; - //************************************************************************/ - - } - //************************************************************************/ - stoptime = MPI_Wtime(); - // cout << "CPU time: " << (stoptime - starttime) << " seconds" << endl; - cputime = stoptime - starttime; - // cout << "Lattice update rate: "<< double(Nx*Ny*Nz*timestep)/cputime/1000000 << " MLUPS" << endl; - double MLUPS = double(Np*timestep)/cputime/1000000; - // if (rank==0) printf("********************************************************\n"); - // if (rank==0) printf("CPU time = %f \n", cputime); - // if (rank==0) printf("Lattice update rate (per process)= %f MLUPS \n", MLUPS); - MLUPS *= nprocs; - - double *Vz; - Vz= new double [Np]; - int SIZE=Np*sizeof(double); - ScaLBL_D3Q19_Momentum(dist,Velocity, Np); - ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - ScaLBL_CopyToHost(&Vz[0],&Velocity[2*Np],SIZE); - - if (rank == 0) printf("Force: %f,%f,%f \n",Fx,Fy,Fz); + MRT.ReadParams(filename); + MRT.SetDomain(); + MRT.Create(); // creating the model will create data structure to match the pore structure and allocate variables + MRT.Initialize(); // initializing the model will set initial conditions for variables + MRT.Run(); + double *Vz; Vz= new double [MRT.Np]; + MRT.VelocityField(Vz); + if (rank == 0) printf("Force: %f,%f,%f \n",MRT.Fx,MRT.Fy,MRT.Fz); + double mu = MRT.mu; + int Nx = MRT.Nx; + int Ny = MRT.Ny; + int Nz = MRT.Nz; + double Fz = MRT.Fz; double vz; double W = 1.f*Nx-4.f; j=Ny/2; k=Nz/2; if (rank == 0) printf("Channel width=%f \n",W); if (rank == 0) printf("ID flag vz analytical\n"); - MPI_Barrier(comm); + if (rank == 0) { for (i=0;iid[n]); - - n = Map(i,j,k); + printf("%i ",MRT.Mask->id[n]); + n = MRT.Map(i,j,k); //printf("%i,%i,%i; %i :",i,j,k,n); if (n<0) {vz =0.f; printf(" b "); } else { vz=Vz[n]; printf(" a "); } @@ -260,14 +70,11 @@ int main(int argc, char **argv) } printf("\n"); } - - if (rank == 1) { for (i=0;iid[n]); - - n = Map(i,j,k); + printf("%i ",MRT.Mask->id[n]); + n = MRT.Map(i,j,k); //printf("%i,%i,%i; %i :",i,j,k,n); if (n<0) {vz =0.f; printf(" b "); } else { vz=Vz[n]; printf(" a "); } @@ -280,9 +87,6 @@ int main(int argc, char **argv) } printf("\n"); } - - - } // ****************************************************