diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index 07aa3f1d..eace3f3f 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -1175,6 +1175,115 @@ void ScaLBL_Communicator::RecvGrad(double *phi, double *grad){ //................................................................................... } +/** + * + */ +void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component){ + + // NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2 + if (Lock==true){ + ERROR("ScaLBL Error (SendD3Q19): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?"); + } + else{ + Lock=true; + } + // assign tag of 19 to D3Q19 communication + sendtag = recvtag = 7; + ScaLBL_DeviceBarrier(); + // Pack the distributions + //...Packing for x face(2,8,10,12,14)................................ + ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,&Aq[Component*N],N); + + MPI_Isend(sendbuf_x, sendCount_x,MPI_DOUBLE,rank_x,sendtag,MPI_COMM_SCALBL,&req1[0]); + MPI_Irecv(recvbuf_X, recvCount_X,MPI_DOUBLE,rank_X,recvtag,MPI_COMM_SCALBL,&req2[0]); + + //...Packing for X face(1,7,9,11,13)................................ + ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,&Aq[Component*N],N); + + MPI_Isend(sendbuf_X, sendCount_X,MPI_DOUBLE,rank_X,sendtag,MPI_COMM_SCALBL,&req1[1]); + MPI_Irecv(recvbuf_x, recvCount_x,MPI_DOUBLE,rank_x,recvtag,MPI_COMM_SCALBL,&req2[1]); + + //...Packing for y face(4,8,9,16,18)................................. + ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,&Aq[Component*N],N); + + MPI_Isend(sendbuf_y, 2*sendCount_y,MPI_DOUBLE,rank_y,sendtag,MPI_COMM_SCALBL,&req1[2]); + MPI_Irecv(recvbuf_Y, 2*recvCount_Y,MPI_DOUBLE,rank_Y,recvtag,MPI_COMM_SCALBL,&req2[2]); + + //...Packing for Y face(3,7,10,15,17)................................. + ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,&Aq[Component*N],N); + ScaLBL_D3Q19_Pack(3,dvcSendList_Y,sendCount_Y,sendCount_Y,sendbuf_Y,Bq,N); + + MPI_Isend(sendbuf_Y, sendCount_Y,MPI_DOUBLE,rank_Y,sendtag,MPI_COMM_SCALBL,&req1[3]); + MPI_Irecv(recvbuf_y, recvCount_y,MPI_DOUBLE,rank_y,recvtag,MPI_COMM_SCALBL,&req2[3]); + + //...Packing for z face(6,12,13,16,17)................................ + ScaLBL_D3Q19_Pack(6,dvcSendList_z,0,sendCount_z,sendbuf_z,&Aq[Component*N],N); + + MPI_Isend(sendbuf_z, sendCount_z,MPI_DOUBLE,rank_z,sendtag,MPI_COMM_SCALBL,&req1[4]); + MPI_Irecv(recvbuf_Z, recvCount_Z,MPI_DOUBLE,rank_Z,recvtag,MPI_COMM_SCALBL,&req2[4]); + + //...Packing for Z face(5,11,14,15,18)................................ + ScaLBL_D3Q19_Pack(5,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,&Aq[Component*N],N); + + //................................................................................... + // Send all the distributions + MPI_Isend(sendbuf_Z, sendCount_Z,MPI_DOUBLE,rank_Z,sendtag,MPI_COMM_SCALBL,&req1[5]); + MPI_Irecv(recvbuf_z, recvCount_z,MPI_DOUBLE,rank_z,recvtag,MPI_COMM_SCALBL,&req2[5]); + +} + + +void ScaLBL_Communicator::RecvD3Q7AA(double *Aq, int Component){ + + // NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2 + //................................................................................... + // Wait for completion of D3Q19 communication + MPI_Waitall(6,req1,stat1); + MPI_Waitall(6,req2,stat2); + ScaLBL_DeviceBarrier(); + + //................................................................................... + // NOTE: AA Routine writes to opposite + // Unpack the distributions on the device + //................................................................................... + //...Unpacking for x face(2,8,10,12,14)................................ + ScaLBL_D3Q7_Unpack(2,dvcRecvDist_x,0,recvCount_x,recvbuf_x,&Aq[Component*N],N); + //................................................................................... + //...Packing for X face(1,7,9,11,13)................................ + ScaLBL_D3Q7_Unpack(1,dvcRecvDist_X,0,recvCount_X,recvbuf_X,&Aq[Component*N],N); + //................................................................................... + //...Packing for y face(4,8,9,16,18)................................. + ScaLBL_D3Q7_Unpack(4,dvcRecvDist_y,0,recvCount_y,recvbuf_y,&Aq[Component*N],N); + //................................................................................... + //...Packing for Y face(3,7,10,15,17)................................. + ScaLBL_D3Q7_Unpack(3,dvcRecvDist_Y,0,recvCount_Y,recvbuf_Y,&Aq[Component*N],N); + //................................................................................... + + if (BoundaryCondition > 0){ + if (kproc != 0){ + //...Packing for z face(6,12,13,16,17)................................ + ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,&Aq[Component*N],N); + } + if (kproc != nprocz-1){ + //...Packing for Z face(5,11,14,15,18)................................ + ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,&Aq[Component*N],N); + } + } + else { + //...Packing for z face(6,12,13,16,17)................................ + ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,&Aq[Component*N],N); + //...Packing for Z face(5,11,14,15,18)................................ + ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,&Aq[Component*N],N); + } + + //................................................................................... + Lock=false; // unlock the communicator after communications complete + //................................................................................... + +} +/* + */ + void ScaLBL_Communicator::BiSendD3Q7AA(double *Aq, double *Bq){ diff --git a/common/ScaLBL.h b/common/ScaLBL.h index dec8b3d1..23cf6936 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -72,6 +72,17 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz, double *Poros,double *Perm, double *Velocity,double Den,double *Pressure); +// ION TRANSPORT MODEL +extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Velocity, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz); + +extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Velocity, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz); + + +// ION TRANSPORT MODEL +extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz); + +extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz); + // MRT MODEL extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx, @@ -185,6 +196,8 @@ public: int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np); void SendD3Q19AA(double *dist); void RecvD3Q19AA(double *dist); + void SendD3Q7AA(double *fq, int Component); + void RecvD3Q7AA(double *fq, int Component) void BiSendD3Q7AA(double *Aq, double *Bq); void BiRecvD3Q7AA(double *Aq, double *Bq); void TriSendD3Q7AA(double *Aq, double *Bq, double *Cq); diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp new file mode 100644 index 00000000..569c9a0a --- /dev/null +++ b/cpu/Ion.cpp @@ -0,0 +1,123 @@ +extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Velocity, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){ + int n; + // conserved momemnts + double rho,ux,uy,uz,uu; + // non-conserved moments + double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; + + for (int n=start; n 10Np => odd part of dist) + f1 = dist[nr1]; // reading the f1 data into register fq + + nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) + f2 = dist[nr2]; // reading the f2 data into register fq + + // q=3 + nr3 = neighborList[n+2*Np]; // neighbor 4 + f3 = dist[nr3]; + + // q = 4 + nr4 = neighborList[n+3*Np]; // neighbor 3 + f4 = dist[nr4]; + + // q=5 + nr5 = neighborList[n+4*Np]; + f5 = dist[nr5]; + + // q = 6 + nr6 = neighborList[n+5*Np]; + f6 = dist[nr6]; + + rho = f0+f2+f1+f4+f3+f6; + ux = Velocity[n]; + uy = Velocity[n+Np]; + uz = Velocity[n+2*Np]; + uu = 1.5*(ux*ux+uy*uy+uz*uz); + + // q=0 + dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*(1.0-uu); + + // q = 1 + dist[nr2] = f1*(1.0-rlx) + rlx*0.05555555555555555*(rho + 3.0*ux + 4.5*ux*ux - uu) + 0.16666666*Fx; + + // q=2 + dist[nr1] = f2*(1.0-rlx) + rlx*0.05555555555555555*(rho - 3.0*ux + 4.5*ux*ux - uu)- 0.16666666*Fx; + + // q = 3 + dist[nr4] = f3*(1.0-rlx) + + rlx*0.05555555555555555*(rho + 3.0*uy + 4.5*uy*uy - uu) + 0.16666666*Fy; + + // q = 4 + dist[nr3] = f4*(1.0-rlx) + + rlx*0.05555555555555555*(rho - 3.0*uy + 4.5*uy*uy - uu)- 0.16666666*Fy; + + // q = 5 + dist[nr6] = f5*(1.0-rlx) + + rlx*0.05555555555555555*(rho + 3.0*uz + 4.5*uz*uz - uu) + 0.16666666*Fz; + + // q = 6 + dist[nr5] = f6*(1.0-rlx) + + rlx*0.05555555555555555*(rho - 3.0*uz + 4.5*uz*uz - uu) - 0.16666666*Fz; + + + } +} \ No newline at end of file diff --git a/cpu/Poisson.cpp b/cpu/Poisson.cpp new file mode 100644 index 00000000..355e4223 --- /dev/null +++ b/cpu/Poisson.cpp @@ -0,0 +1,123 @@ +extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){ + int n; + // conserved momemnts + double rho,ux,uy,uz,uu; + // non-conserved moments + double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; + + for (int n=start; n 10Np => odd part of dist) + f1 = dist[nr1]; // reading the f1 data into register fq + + nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) + f2 = dist[nr2]; // reading the f2 data into register fq + + // q=3 + nr3 = neighborList[n+2*Np]; // neighbor 4 + f3 = dist[nr3]; + + // q = 4 + nr4 = neighborList[n+3*Np]; // neighbor 3 + f4 = dist[nr4]; + + // q=5 + nr5 = neighborList[n+4*Np]; + f5 = dist[nr5]; + + // q = 6 + nr6 = neighborList[n+5*Np]; + f6 = dist[nr6]; + + rho = f0+f2+f1+f4+f3+f6+f5; + ux = f1-f2; + uy = f3-f4; + uz = f5-f6; + uu = 1.5*(ux*ux+uy*uy+uz*uz); + + + // q=0 + dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*(1.0-uu); + + // q = 1 + dist[nr2] = f1*(1.0-rlx) + rlx*0.05555555555555555*(rho + 3.0*ux + 4.5*ux*ux - uu) + 0.16666666*Fx; + + // q=2 + dist[nr1] = f2*(1.0-rlx) + rlx*0.05555555555555555*(rho - 3.0*ux + 4.5*ux*ux - uu)- 0.16666666*Fx; + + // q = 3 + dist[nr4] = f3*(1.0-rlx) + + rlx*0.05555555555555555*(rho + 3.0*uy + 4.5*uy*uy - uu) + 0.16666666*Fy; + + // q = 4 + dist[nr3] = f4*(1.0-rlx) + + rlx*0.05555555555555555*(rho - 3.0*uy + 4.5*uy*uy - uu)- 0.16666666*Fy; + + // q = 5 + dist[nr6] = f5*(1.0-rlx) + + rlx*0.05555555555555555*(rho + 3.0*uz + 4.5*uz*uz - uu) + 0.16666666*Fz; + + // q = 6 + dist[nr5] = f6*(1.0-rlx) + + rlx*0.05555555555555555*(rho - 3.0*uz + 4.5*uz*uz - uu) - 0.16666666*Fz; + + + } +} \ No newline at end of file diff --git a/models/IonModel.cpp b/models/IonModel.cpp new file mode 100644 index 00000000..38088587 --- /dev/null +++ b/models/IonModel.cpp @@ -0,0 +1,225 @@ +/* + * Multi-relaxation time LBM Model + */ +#include "models/MRT.h" +#include "models/ElectroModel.h" +#include "analysis/distance.h" +#include "common/ReadMicroCT.h" + +ScaLBL_IonModel::ScaLBL_IonModel(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_IonModel::~ScaLBL_IonModel(){ + +} + +void ScaLBL_IonModel::ReadParams(string filename){ + // read the input database + db = std::make_shared( filename ); + domain_db = db->getDatabase( "Domain" ); + ion_db = db->getDatabase( "Ions" ); + + tau = 1.0; + timestepMax = 100000; + tolerance = 1.0e-8; + Fx = Fy = 0.0; + Fz = 1.0e-5; + + // Color Model parameters + if (ion_db->keyExists( "timestepMax" )){ + timestepMax = mrt_db->getScalar( "timestepMax" ); + } + + mu=(tau-0.5)/3.0; +} +void ScaLBL_IonModel::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 + + // domain parameters + Nx = Dm->Nx; + Ny = Dm->Ny; + Nz = Dm->Nz; + Lx = Dm->Lx; + Ly = Dm->Ly; + Lz = Dm->Lz; + + N = Nx*Ny*Nz; + Distance.resize(Nx,Ny,Nz); + Velocity_x.resize(Nx,Ny,Nz); + Velocity_y.resize(Nx,Ny,Nz); + Velocity_z.resize(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); + + rank = Dm->rank(); + nprocx = Dm->nprocx(); + nprocy = Dm->nprocy(); + nprocz = Dm->nprocz(); +} + +void ScaLBL_IonModel::ReadInput(){ + + sprintf(LocalRankString,"%05d",Dm->rank()); + sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); + sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString); + + + if (domain_db->keyExists( "Filename" )){ + auto Filename = domain_db->getScalar( "Filename" ); + Mask->Decomp(Filename); + } + else if (domain_db->keyExists( "GridFile" )){ + // Read the local domain data + auto input_id = readMicroCT( *domain_db, comm ); + // Fill the halo (assuming GCW of 1) + array size0 = { (int) input_id.size(0), (int) input_id.size(1), (int) input_id.size(2) }; + ArraySize size1 = { (size_t) Mask->Nx, (size_t) Mask->Ny, (size_t) Mask->Nz }; + ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 ); + fillHalo fill( comm, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 ); + Array id_view; + id_view.viewRaw( size1, Mask->id ); + fill.copy( input_id, id_view ); + fill.fill( id_view ); + } + else{ + Mask->ReadIDs(); + } + + // Generate the signed distance map + // Initialize the domain and communication + Array id_solid(Nx,Ny,Nz); + // Solve for the position of the solid phase + for (int k=0;kid[n] > 0) id_solid(i,j,k) = 1; + else id_solid(i,j,k) = 0; + } + } + } + // Initialize the signed distance function + for (int k=0;kSDs); + if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); + CalcDist(Distance,id_solid,*Dm); + if (rank == 0) cout << "Domain set." << endl; +} + +void ScaLBL_IonModel::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, number_ion_species*7*dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &Ci, number_ion_species*sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &ChargeDensity, 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_IonModel::Initialize(){ + /* + * This function initializes model + */ + if (rank==0) printf ("Initializing distributions \n"); + ScaLBL_D3Q19_Init(fq, Np); +} + +void ScaLBL_IonModel::Run(double *Velocity){ + + //.......create and start timer............ + double starttime,stoptime,cputime; + ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + starttime = MPI_Wtime(); + if (rank==0) printf("Beginning AA timesteps, timestepMax = %i \n", timestepMax); + if (rank==0) printf("********************************************************\n"); + timestep=0; + double error = 1.0; + double flow_rate_previous = 0.0; + while (timestep < timestepMax && error > tolerance) { + //************************************************************************/ + timestep++; + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + ScaLBL_D3Q7_AAodd_Ion(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + // Set boundary conditions + /* ... */ + ScaLBL_D3Q7_AAodd_Ion(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_D3Q7_AAeven_Ion(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + // Set boundary conditions + /* ... */ + ScaLBL_D3Q7_AAeven_Ion(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + //************************************************************************/ + } + //************************************************************************/ + stoptime = MPI_Wtime(); + if (rank==0) printf("-------------------------------------------------------------------\n"); + // Compute the walltime per timestep + cputime = (stoptime - starttime)/timestep; + // Performance obtained from each node + double MLUPS = double(Np)/cputime/1000000; + + if (rank==0) printf("********************************************************\n"); + if (rank==0) printf("CPU time = %f \n", cputime); + if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); + MLUPS *= nprocs; + if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); + if (rank==0) printf("********************************************************\n"); + +} + diff --git a/models/IonModel.h b/models/IonModel.h new file mode 100644 index 00000000..f6ffad5e --- /dev/null +++ b/models/IonModel.h @@ -0,0 +1,72 @@ +/* + * Multi-relaxation time LBM Model + */ +#include +#include +#include +#include +#include +#include +#include + +#include "common/ScaLBL.h" +#include "common/Communication.h" +#include "common/MPI_Helpers.h" +#include "analysis/Minkowski.h" +#include "ProfilerApp.h" + +class ScaLBL_IonModel{ +public: + ScaLBL_IonModel(int RANK, int NP, MPI_Comm COMM); + ~ScaLBL_IonModel(); + + // 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(double *Velocity); + void VelocityField(); + + bool Restart,pBC; + int timestep,timestepMax; + int BoundaryCondition; + double tau,mu; + double Fx,Fy,Fz,flux; + double din,dout; + double tolerance; + + int number_ion_species; + + 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; + // input database + std::shared_ptr db; + std::shared_ptr domain_db; + std::shared_ptr ion_db; + + IntArray Map; + DoubleArray Distance; + int *NeighborList; + double *fq; + double *Ci; + double *ChargeDensity; + +private: + MPI_Comm comm; + + // filenames + char LocalRankString[8]; + char LocalRankFilename[40]; + char LocalRestartFile[40]; + + //int rank,nprocs; + void LoadParams(std::shared_ptr db0); +}; diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp new file mode 100644 index 00000000..f4f15224 --- /dev/null +++ b/models/PoissonSolver.cpp @@ -0,0 +1,253 @@ +/* + * Multi-relaxation time LBM Model + */ +#include "models/MRT.h" +#include "models/ElectroModel.h" +#include "analysis/distance.h" +#include "common/ReadMicroCT.h" + +ScaLBL_Poisson::ScaLBL_Poisson(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_Poisson::~ScaLBL_Poisson(){ + +} + +void ScaLBL_Poisson::ReadParams(string filename){ + // read the input database + db = std::make_shared( filename ); + domain_db = db->getDatabase( "Domain" ); + mrt_db = db->getDatabase( "MRT" ); + electric_db = db->getDatabase( "Electrochemistry" ); + + tau = 1.0; + timestepMax = 100000; + tolerance = 1.0e-8; + Fx = Fy = 0.0; + Fz = 1.0e-5; + + // Color Model parameters + if (mrt_db->keyExists( "timestepMax" )){ + timestepMax = mrt_db->getScalar( "timestepMax" ); + } + if (mrt_db->keyExists( "tolerance" )){ + tolerance = mrt_db->getScalar( "tolerance" ); + } + if (mrt_db->keyExists( "tau" )){ + tau = mrt_db->getScalar( "tau" ); + } + if (mrt_db->keyExists( "F" )){ + Fx = mrt_db->getVector( "F" )[0]; + Fy = mrt_db->getVector( "F" )[1]; + Fz = mrt_db->getVector( "F" )[2]; + } + if (mrt_db->keyExists( "Restart" )){ + Restart = mrt_db->getScalar( "Restart" ); + } + if (mrt_db->keyExists( "din" )){ + din = mrt_db->getScalar( "din" ); + } + if (mrt_db->keyExists( "dout" )){ + dout = mrt_db->getScalar( "dout" ); + } + if (mrt_db->keyExists( "flux" )){ + flux = mrt_db->getScalar( "flux" ); + } + + // Read domain parameters + if (domain_db->keyExists( "BC" )){ + BoundaryCondition = domain_db->getScalar( "BC" ); + } + + + mu=(tau-0.5)/3.0; +} +void ScaLBL_Poisson::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 + + // domain parameters + Nx = Dm->Nx; + Ny = Dm->Ny; + Nz = Dm->Nz; + Lx = Dm->Lx; + Ly = Dm->Ly; + Lz = Dm->Lz; + + N = Nx*Ny*Nz; + Distance.resize(Nx,Ny,Nz); + Velocity_x.resize(Nx,Ny,Nz); + Velocity_y.resize(Nx,Ny,Nz); + Velocity_z.resize(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); + + rank = Dm->rank(); + nprocx = Dm->nprocx(); + nprocy = Dm->nprocy(); + nprocz = Dm->nprocz(); +} + +void ScaLBL_Poisson::ReadInput(){ + + sprintf(LocalRankString,"%05d",Dm->rank()); + sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); + sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString); + + + if (domain_db->keyExists( "Filename" )){ + auto Filename = domain_db->getScalar( "Filename" ); + Mask->Decomp(Filename); + } + else if (domain_db->keyExists( "GridFile" )){ + // Read the local domain data + auto input_id = readMicroCT( *domain_db, comm ); + // Fill the halo (assuming GCW of 1) + array size0 = { (int) input_id.size(0), (int) input_id.size(1), (int) input_id.size(2) }; + ArraySize size1 = { (size_t) Mask->Nx, (size_t) Mask->Ny, (size_t) Mask->Nz }; + ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 ); + fillHalo fill( comm, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 ); + Array id_view; + id_view.viewRaw( size1, Mask->id ); + fill.copy( input_id, id_view ); + fill.fill( id_view ); + } + else{ + Mask->ReadIDs(); + } + + // Generate the signed distance map + // Initialize the domain and communication + Array id_solid(Nx,Ny,Nz); + // Solve for the position of the solid phase + for (int k=0;kid[n] > 0) id_solid(i,j,k) = 1; + else id_solid(i,j,k) = 0; + } + } + } + // Initialize the signed distance function + for (int k=0;kSDs); + if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); + CalcDist(Distance,id_solid,*Dm); + if (rank == 0) cout << "Domain set." << endl; +} + +void ScaLBL_Poisson::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, 7*dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &Psi, 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_Poisson::Initialize(){ + /* + * This function initializes model + */ + if (rank==0) printf ("Initializing distributions \n"); + ScaLBL_D3Q19_Init(fq, Np); +} + +void ScaLBL_Poisson::Run(double *ChargeDensity){ + + //.......create and start timer............ + double starttime,stoptime,cputime; + ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + starttime = MPI_Wtime(); + if (rank==0) printf("Beginning AA timesteps, timestepMax = %i \n", timestepMax); + if (rank==0) printf("********************************************************\n"); + timestep=0; + double error = 1.0; + double flow_rate_previous = 0.0; + while (timestep < timestepMax && error > tolerance) { + //************************************************************************/ + timestep++; + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + // Set boundary conditions + /* ... */ + ScaLBL_D3Q7_AAodd_Poisson(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_D3Q7_AAeven_Poisson(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + // Set boundary conditions + /* ... */ + ScaLBL_D3Q7_AAeven_Poisson(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz); + ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + //************************************************************************/ + } + //************************************************************************/ + stoptime = MPI_Wtime(); + if (rank==0) printf("-------------------------------------------------------------------\n"); + // Compute the walltime per timestep + cputime = (stoptime - starttime)/timestep; + // Performance obtained from each node + double MLUPS = double(Np)/cputime/1000000; + + if (rank==0) printf("********************************************************\n"); + if (rank==0) printf("CPU time = %f \n", cputime); + if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); + MLUPS *= nprocs; + if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); + if (rank==0) printf("********************************************************\n"); + +} diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h new file mode 100644 index 00000000..63e10df0 --- /dev/null +++ b/models/PoissonSolver.h @@ -0,0 +1,68 @@ +/* + * Multi-relaxation time LBM Model + */ +#include +#include +#include +#include +#include +#include +#include + +#include "common/ScaLBL.h" +#include "common/Communication.h" +#include "common/MPI_Helpers.h" +#include "analysis/Minkowski.h" +#include "ProfilerApp.h" + +class ScaLBL_Poisson{ +public: + ScaLBL_Poisson(int RANK, int NP, MPI_Comm COMM); + ~ScaLBL_Poisson(); + + // 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(double *ChargeDensity); + + bool Restart,pBC; + int timestep,timestepMax; + int BoundaryCondition; + double tau,mu; + double Fx,Fy,Fz,flux; + double din,dout; + double tolerance; + + 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; + // input database + std::shared_ptr db; + std::shared_ptr domain_db; + std::shared_ptr electric_db; + + IntArray Map; + DoubleArray Distance; + int *NeighborList; + double *fq; + double *Psi; + +private: + MPI_Comm comm; + + // filenames + char LocalRankString[8]; + char LocalRankFilename[40]; + char LocalRestartFile[40]; + + //int rank,nprocs; + void LoadParams(std::shared_ptr db0); +}; diff --git a/tests/lbpm_electrokinetic_dfh_simulator.cpp b/tests/lbpm_electrokinetic_dfh_simulator.cpp new file mode 100644 index 00000000..e3765d12 --- /dev/null +++ b/tests/lbpm_electrokinetic_dfh_simulator.cpp @@ -0,0 +1,78 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "models/DFHModel.h" +#include "models/IonModel.h" +#include "models/PoissonSolver.h" + +//#define WRE_SURFACES + +/* + * Simulator for two-phase flow in porous media + * James E. McClure 2013-2014 + */ + +using namespace std; + +//************************************************************************* +// Implementation of Two-Phase Immiscible LBM using CUDA +//************************************************************************* + +int main(int argc, char **argv) +{ + // Initialize MPI + int provided_thread_support = -1; + MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support); + MPI_Comm comm; + MPI_Comm_dup(MPI_COMM_WORLD,&comm); + int rank = comm_rank(comm); + int nprocs = comm_size(comm); + if ( rank==0 && provided_thread_support