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

add routine for ion model to read from file

Browse Source
This commit is contained in:
Rex Zhe Li
2020-09-28 17:08:49 -04:00
parent a1f375e86c
commit 4657adbc94
7 changed files with 249 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

147
common/Domain.cpp
View File

@@ -1286,3 +1286,150 @@ void ReadBinaryFile(char *FILENAME, double *Data, size_t N)
File.close();
}
void Domain::ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData)
{
//........................................................................................
// Reading the user-defined input file
// NOTE: so far it only supports BC=0 (periodic) and BC=5 (mixed reflection)
// because if checkerboard or inlet/outlet buffer layers are added, the
// value of the void space is undefined.
// NOTE: if BC=5 is used, where the inlet and outlet layers of the domain are modified,
// user needs to modify the input file accordingly before LBPM simulator read
// the input file.
//........................................................................................
int rank_offset = 0;
int RANK = rank();
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz;
int64_t global_Nx,global_Ny,global_Nz;
int64_t i,j,k,n;
//TODO These offset we may still need them
int64_t xStart,yStart,zStart;
xStart=yStart=zStart=0;
// Read domain parameters
// TODO currently the size of the data is still read from Domain{};
// but user may have a user-specified size
auto size = database->getVector<int>( "n" );
auto SIZE = database->getVector<int>( "N" );
auto nproc = database->getVector<int>( "nproc" );
//TODO currently the funcationality "offset" is disabled as the user-defined input data may have a different size from that of the input domain
if (database->keyExists( "offset" )){
auto offset = database->getVector<int>( "offset" );
xStart = offset[0];
yStart = offset[1];
zStart = offset[2];
}
nx = size[0];
ny = size[1];
nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
global_Nx = SIZE[0];
global_Ny = SIZE[1];
global_Nz = SIZE[2];
nprocs=nprocx*nprocy*nprocz;
double *SegData = NULL;
if (RANK==0){
printf("User-defined input file: %s (data type: %s)\n",Filename.c_str(),Datatype.c_str());
printf("NOTE: currently only BC=0 or 5 supports user-defined input file!\n");
// Rank=0 reads the entire segmented data and distributes to worker processes
printf("Dimensions of the user-defined input file: %ld x %ld x %ld \n",global_Nx,global_Ny,global_Nz);
int64_t SIZE = global_Nx*global_Ny*global_Nz;
if (Datatype == "double"){
printf("Reading input data as double precision floating number\n");
SegData = new double[SIZE];
FILE *SEGDAT = fopen(Filename.c_str(),"rb");
if (SEGDAT==NULL) ERROR("Domain.cpp: Error reading user-defined file!\n");
size_t ReadSeg;
ReadSeg=fread(SegData,8,SIZE,SEGDAT);
if (ReadSeg != size_t(SIZE)) printf("Domain.cpp: Error reading file: %s\n",Filename.c_str());
fclose(SEGDAT);
}
else{
ERROR("Error: User-defined input file only supports double-precision floating number!\n");
}
printf("Read file successfully from %s \n",Filename.c_str());
}
// Get the rank info
int64_t N = (nx+2)*(ny+2)*(nz+2);
// number of sites to use for periodic boundary condition transition zone
//int64_t z_transition_size = (nprocz*nz - (global_Nz - zStart))/2;
//if (z_transition_size < 0) z_transition_size=0;
int64_t z_transition_size = 0;
//char LocalRankFilename[1000];//just for debug
double *loc_id;
loc_id = new double [(nx+2)*(ny+2)*(nz+2)];
// Set up the sub-domains
if (RANK==0){
printf("Decomposing user-defined input file\n");
printf("Distributing subdomains across %i processors \n",nprocs);
printf("Process grid: %i x %i x %i \n",nprocx,nprocy,nprocz);
printf("Subdomain size: %i x %i x %i \n",nx,ny,nz);
printf("Size of transition region: %ld \n", z_transition_size);
for (int kp=0; kp<nprocz; kp++){
for (int jp=0; jp<nprocy; jp++){
for (int ip=0; ip<nprocx; ip++){
// rank of the process that gets this subdomain
int rnk = kp*nprocx*nprocy + jp*nprocx + ip;
// Pack and send the subdomain for rnk
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int64_t x = xStart + ip*nx + i-1;
int64_t y = yStart + jp*ny + j-1;
// int64_t z = zStart + kp*nz + k-1;
int64_t z = zStart + kp*nz + k-1 - z_transition_size;
if (x<xStart) x=xStart;
if (!(x<global_Nx)) x=global_Nx-1;
if (y<yStart) y=yStart;
if (!(y<global_Ny)) y=global_Ny-1;
if (z<zStart) z=zStart;
if (!(z<global_Nz)) z=global_Nz-1;
int64_t nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
int64_t nglobal = z*global_Nx*global_Ny+y*global_Nx+x;
loc_id[nlocal] = SegData[nglobal];
}
}
}
if (rnk==0){
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
UserData[nlocal] = loc_id[nlocal];
}
}
}
}
else{
//printf("Sending data to process %i \n", rnk);
MPI_Send(loc_id,N,MPI_DOUBLE,rnk,15,Comm);
}
// Write the data for this rank data
// NOTE just for debug
//sprintf(LocalRankFilename,"%s.%05i",Filename.c_str(),rnk+rank_offset);
//FILE *ID = fopen(LocalRankFilename,"wb");
//fwrite(loc_id,8,(nx+2)*(ny+2)*(nz+2),ID);
//fclose(ID);
}
}
}
}
else{
// Recieve the subdomain from rank = 0
//printf("Ready to recieve data %i at process %i \n", N,rank);
MPI_Recv(UserData,N,MPI_DOUBLE,0,15,Comm,MPI_STATUS_IGNORE);
}
//Comm.barrier();
MPI_Barrier(Comm);
}

1
common/Domain.h
View File

@@ -178,6 +178,7 @@ public: // Public variables (need to create accessors instead)
void ReadIDs();
void Decomp( const std::string& filename );
void ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData);
void CommunicateMeshHalo(DoubleArray &Mesh);
void CommInit();
int PoreCount();

1
common/ScaLBL.h
View File

@@ -85,6 +85,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Veloci
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np);
extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np);
extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np);

16
cpu/Ion.cpp
View File

@@ -220,6 +220,22 @@ extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit,
}
}
extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np)
{
int n;
double DenInit;
for (n=0; n<Np; n++){
DenInit = Den[n];
dist[0*Np+n] = 0.25*DenInit;
dist[1*Np+n] = 0.125*DenInit;
dist[2*Np+n] = 0.125*DenInit;
dist[3*Np+n] = 0.125*DenInit;
dist[4*Np+n] = 0.125*DenInit;
dist[5*Np+n] = 0.125*DenInit;
dist[6*Np+n] = 0.125*DenInit;
}
}
extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
int n;

33
gpu/Ion.cu
View File

@@ -263,6 +263,27 @@ __global__ void dvc_ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenI
}
}
__global__ void dvc_ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np){
int n;
double DenInit;
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
if (n<Np) {
DenInit = Den[n];
dist[0*Np+n] = 0.25*DenInit;
dist[1*Np+n] = 0.125*DenInit;
dist[2*Np+n] = 0.125*DenInit;
dist[3*Np+n] = 0.125*DenInit;
dist[4*Np+n] = 0.125*DenInit;
dist[5*Np+n] = 0.125*DenInit;
dist[6*Np+n] = 0.125*DenInit;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
int n;
@@ -345,6 +366,18 @@ extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit,
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_Ion_Init_FromFile<<<NBLOCKS,NTHREADS >>>(dist,Den,Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_D3Q7_Ion_Init_FromFile: %s \n",cudaGetErrorString(err));
}
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
//cudaProfilerStart();

55
models/IonModel.cpp
View File

@@ -136,7 +136,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
}
//read initial ion concentration list; INPUT unit [mol/m^3]
//it must be converted to LB unit [mol/lu^3]
if (ion_db->keyExists("IonConcentrationList")){
if (ion_db->keyExists("IonConcentrationList")){
IonConcentration.clear();
IonConcentration = ion_db->getVector<double>( "IonConcentrationList" );
if (IonConcentration.size()!=number_ion_species){
@@ -544,6 +544,35 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
}
}
void ScaLBL_IonModel::AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion)
{
double *Ci_host;
Ci_host = new double[N];
double VALUE=0.f;
Mask->ReadFromFile(File_ion[0],File_ion[1],Ci_host);
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int idx = Map(i,j,k);
if (!(idx < 0)){
int n = k*Nx*Ny+j*Nx+i;
//NOTE: default user-input unit: mol/m^3
VALUE = Ci_host[n]*(h*h*h*1.0e-18);
if (VALUE<0.0){
ERROR("Error: Ion concentration value must be a positive number! \n");
}
else{
Ci[idx] = VALUE;
}
}
}
}
}
delete [] Ci_host;
}
void ScaLBL_IonModel::Create(){
/*
* This function creates the variables needed to run a LBM
@@ -601,8 +630,6 @@ void ScaLBL_IonModel::Create(){
ScaLBL_DeviceBarrier();
delete [] IonSolid_host;
}
}
void ScaLBL_IonModel::Initialize(){
@@ -610,8 +637,26 @@ void ScaLBL_IonModel::Initialize(){
* This function initializes model
*/
if (rank==0) printf ("LB Ion Solver: initializing D3Q7 distributions\n");
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init(&fq[ic*Np*7],&Ci[ic*Np],IonConcentration[ic],Np);
if (ion_db->keyExists("IonConcentrationFile")){
//TODO: Need to figure out how to deal with multi-species concentration initialization
//NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
auto File_ion = ion_db->getVector<std::string>( "IonConcentrationFile" );
double *Ci_host;
Ci_host = new double[number_ion_species*Np];
for (int ic=0; ic<number_ion_species; ic++){
AssignIonConcentration_FromFile(&Ci_host[ic*Np],File_ion);
}
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species*sizeof(double)*Np);
MPI_Barrier(comm);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic*Np*7],&Ci[ic*Np],Np);
}
delete [] Ci_host;
}
else{
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init(&fq[ic*Np*7],&Ci[ic*Np],IonConcentration[ic],Np);
}
}
if (rank==0) printf ("LB Ion Solver: initializing charge density\n");
for (int ic=0; ic<number_ion_species; ic++){

1
models/IonModel.h
View File

@@ -89,5 +89,6 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void AssignSolidBoundary(double *ion_solid);
void AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion);
void IonConcentration_LB_to_Phys(DoubleArray &Den_reg);
};