OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

and pour the waters of the nile

Browse Source
This commit is contained in:
James E McClure
2018-05-19 19:42:58 -04:00
parent c5994631ee
commit baf501a14a
3 changed files with 271 additions and 217 deletions
Download Patch File Download Diff File Expand all files Collapse all files

181
models/MRTModel.cpp Normal file
View File

@@ -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<Database>( filename );
domain_db = db->getDatabase( "Domain" );
mrt_db = db->getDatabase( "MRT" );
// Color Model parameters
timestepMax = mrt_db->getScalar<int>( "timestepMax" );
tauA = mrt_db->getScalar<double>( "tau" );
Fx = mrt_db->getVector<double>( "F" )[0];
Fy = mrt_db->getVector<double>( "F" )[1];
Fz = mrt_db->getVector<double>( "F" )[2];
Restart = mrt_db->getScalar<bool>( "Restart" );
din = mrt_db->getScalar<double>( "din" );
dout = mrt_db->getScalar<double>( "dout" );
flux = mrt_db->getScalar<double>( "flux" );
// Read domain parameters
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
BoundaryCondition = domain_db->getScalar<int>( "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<Domain>(new Domain(domain_db,comm)); // full domain for analysis
Mask = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // mask domain removes immobile phases
Nx+=2; Ny+=2; Nz += 2;
N = Nx*Ny*Nz;
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
//Averages = std::shared_ptr<TwoPhase> ( 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; i<Nx*Ny*Nz; i++) Dm->id[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<ScaLBL_Communicator>(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);
}

69
models/MRTModel.h Normal file
View File

@@ -0,0 +1,69 @@
/*
* Multi-relaxation time LBM Model
*/
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#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<Database> 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<Domain> Dm; // this domain is for analysis
std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
std::shared_ptr<TwoPhase> Averages;
// input database
std::shared_ptr<Database> db;
std::shared_ptr<Database> domain_db;
std::shared_ptr<Database> color_db;
std::shared_ptr<Database> 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<Database> db0);
};

238
tests/TestPoiseuille.cpp
View File

@@ -8,6 +8,7 @@
#include <fstream>
#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<Database>( FILENAME );
auto domain_db = db->getDatabase( "Domain" );
int Nx = domain_db->getVector<int>( "n" )[0];
int Ny = domain_db->getVector<int>( "n" )[1];
int Nz = domain_db->getVector<int>( "n" )[2];
int nprocx = domain_db->getVector<int>( "nproc" )[0];
int nprocy = domain_db->getVector<int>( "nproc" )[1];
int nprocz = domain_db->getVector<int>( "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<Domain> 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;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
if (i<2) Dm->id[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;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (Dm->id[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_Communicator> 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;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
printf("%i ",Dm->id[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;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
printf("%i ",Dm->id[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");
}
}
// ****************************************************