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

merging electrokinetic & greyscale

Browse Source
This commit is contained in:
James McClure
2021-01-04 14:00:44 -05:00
parent 5c914a5719 d02eff3017
commit 2fc85dead8
22 changed files with 9586 additions and 224 deletions
Download Patch File Download Diff File Expand all files Collapse all files

141
models/GreyscaleModel.cpp
View File

@@ -44,7 +44,7 @@ void ScaLBL_GreyscaleModel::ReadParams(string filename){
din=dout=1.0;
flux=0.0;
dp = 10.0; //unit of 'dp': voxel
CollisionType = 1; //1: IMRT; 2: BGK
CollisionType = 1; //1: IMRT; 2: BGK; 3: MRT
// ---------------------- Greyscale Model parameters -----------------------//
if (greyscale_db->keyExists( "timestepMax" )){
@@ -84,6 +84,9 @@ void ScaLBL_GreyscaleModel::ReadParams(string filename){
if (collision == "BGK"){
CollisionType=2;
}
else if (collision == "MRT"){
CollisionType=3;
}
// ------------------------------------------------------------------------//
//------------------------ Other Domain parameters ------------------------//
@@ -202,9 +205,9 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
for (int idx=0; idx<NLABELS; idx++) label_count[idx]=0;
for (int k=1;k<Nz-1;k++){
for (int j=1;j<Ny-1;j++){
for (int i=1;i<Nx-1;i++){
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
VALUE=id[n];
// Assign the affinity from the paired list
@@ -233,9 +236,9 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
if (NLABELS != PermeabilityList.size()){
ERROR("Error: ComponentLabels and PermeabilityList must be the same length! \n");
}
for (int k=1;k<Nz-1;k++){
for (int j=1;j<Ny-1;j++){
for (int i=1;i<Nx-1;i++){
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
VALUE=id[n];
// Assign the affinity from the paired list
@@ -274,7 +277,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
if (rank==0){
printf("Image resolution: %.5g [um/voxel]\n",Dm->voxel_length);
printf("Component labels: %lu \n",NLABELS);
printf("Number of component labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
POROSITY=PorosityList[idx];
@@ -288,6 +291,56 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
}
}
void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity,double *Permeability,const vector<std::string> &File_poro,const vector<std::string> &File_perm)
{
double *Porosity_host, *Permeability_host;
Porosity_host = new double[N];
Permeability_host = new double[N];
double POROSITY=0.f;
double PERMEABILITY=0.f;
//Initialize a weighted porosity after considering grey voxels
double GreyPorosity_loc=0.0;
GreyPorosity=0.0;
//double label_count_loc = 0.0;
//double label_count_glb = 0.0;
Mask->ReadFromFile(File_poro[0],File_poro[1],Porosity_host);
Mask->ReadFromFile(File_perm[0],File_perm[1],Permeability_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;
POROSITY = Porosity_host[n];
PERMEABILITY = Permeability_host[n];
if (POROSITY<=0.0){
ERROR("Error: Porosity for grey voxels must be 0.0 < Porosity <= 1.0 !\n");
}
else if (PERMEABILITY<=0.0){
ERROR("Error: Permeability for grey voxel must be > 0.0 ! \n");
}
else{
Porosity[idx] = POROSITY;
Permeability[idx] = PERMEABILITY;
GreyPorosity_loc += POROSITY;
//label_count_loc += 1.0;
}
}
}
}
}
GreyPorosity = sumReduce( Dm->Comm, GreyPorosity_loc);
GreyPorosity = GreyPorosity/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (rank==0){
printf("Image resolution: %.5g [um/voxel]\n",Dm->voxel_length);
printf("The weighted porosity, considering both open and grey voxels, is %.3g\n",GreyPorosity);
}
delete [] Porosity_host;
delete [] Permeability_host;
}
void ScaLBL_GreyscaleModel::Create(){
/*
@@ -326,7 +379,6 @@ void ScaLBL_GreyscaleModel::Create(){
neighborSize=18*(Np*sizeof(int));
//...........................................................................
ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Permeability, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Porosity, sizeof(double)*Np);
@@ -334,47 +386,31 @@ void ScaLBL_GreyscaleModel::Create(){
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("Setting up device map and neighbor list \n");
if (rank==0) printf ("Setting up device neighbor list \n");
fflush(stdout);
int *TmpMap;
TmpMap=new int[Np];
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
if (!(idx < 0))
TmpMap[idx] = k*Nx*Ny+j*Nx+i;
}
}
}
// check that TmpMap is valid
for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
int n = TmpMap[idx];
if (n > Nx*Ny*Nz){
printf("Bad value! idx=%i \n");
TmpMap[idx] = Nx*Ny*Nz-1;
}
}
for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
int n = TmpMap[idx];
if (n > Nx*Ny*Nz){
printf("Bad value! idx=%i \n");
TmpMap[idx] = Nx*Ny*Nz-1;
}
}
ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
ScaLBL_DeviceBarrier();
delete [] TmpMap;
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
// initialize phi based on PhaseLabel (include solid component labels)
double *Poros, *Perm;
Poros = new double[Np];
Perm = new double[Np];
AssignComponentLabels(Poros,Perm);
Perm = new double[Np];
if (greyscale_db->keyExists("FileVoxelPorosityMap")){
//NOTE: FileVoxel**Map is a vector, including "file_name, datatype"
auto File_poro = greyscale_db->getVector<std::string>( "FileVoxelPorosityMap" );
auto File_perm = greyscale_db->getVector<std::string>( "FileVoxelPermeabilityMap" );
AssignComponentLabels(Poros,Perm,File_poro,File_perm);
}
else if (greyscale_db->keyExists("PorosityList")){
//initialize voxel porosity and perm from the input list
AssignComponentLabels(Poros,Perm);
}
else {
ERROR("Error: PorosityList or FilenameVoxelPorosityMap cannot be found! \n");
}
ScaLBL_CopyToDevice(Porosity, Poros, Np*sizeof(double));
ScaLBL_CopyToDevice(Permeability, Perm, Np*sizeof(double));
delete [] Poros;
delete [] Perm;
}
@@ -391,6 +427,10 @@ void ScaLBL_GreyscaleModel::Initialize(){
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf("Collision model: BGK.\n");
}
else if (CollisionType==3){
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf("Collision model: MRT.\n");
}
else{
if (rank==0) printf("Unknown collison type! IMRT collision is used.\n");
ScaLBL_D3Q19_GreyIMRT_Init(fq, Np, Den);
@@ -472,6 +512,9 @@ void ScaLBL_GreyscaleModel::Run(){
case 2:
ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
break;
case 3:
ScaLBL_D3Q19_AAodd_Greyscale_MRT(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
default:
ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
@@ -490,6 +533,9 @@ void ScaLBL_GreyscaleModel::Run(){
case 2:
ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
break;
case 3:
ScaLBL_D3Q19_AAodd_Greyscale_MRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
default:
ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
@@ -506,6 +552,9 @@ void ScaLBL_GreyscaleModel::Run(){
case 2:
ScaLBL_D3Q19_AAeven_Greyscale(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
break;
case 3:
ScaLBL_D3Q19_AAeven_Greyscale_MRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
default:
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
@@ -524,6 +573,9 @@ void ScaLBL_GreyscaleModel::Run(){
case 2:
ScaLBL_D3Q19_AAeven_Greyscale(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
break;
case 3:
ScaLBL_D3Q19_AAeven_Greyscale_MRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
default:
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
@@ -597,11 +649,6 @@ void ScaLBL_GreyscaleModel::Run(){
}
}
}
//MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
//MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
//MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
//MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
vax = sumReduce( Mask->Comm, vax_loc);
vay = sumReduce( Mask->Comm, vay_loc);
vaz = sumReduce( Mask->Comm, vaz_loc);