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

Some updates:(1)add different fq initialization for BGK and IMRT;(2)user can choose collision model

Browse Source
This commit is contained in:
Rex Zhe Li
2020-01-31 15:15:26 -05:00
parent 34c5b223b0
commit 69ee9f79cb
5 changed files with 136 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
common/ScaLBL.h
View File

@@ -46,6 +46,7 @@ extern "C" void ScaLBL_UnpackDenD3Q7(int *list, int count, double *recvbuf, int
extern "C" void ScaLBL_D3Q19_Init(double *Dist, int Np);
extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np);
extern "C" void ScaLBL_D3Q19_Pressure(double *dist, double *press, int Np);
@@ -56,6 +57,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int
extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz);
// GREYSCALE MODEL
extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *Dist, int Np, double Den);
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz,
double *Poros,double *Perm, double *Velocity,double *Pressure);

27
cpu/D3Q19.cpp
View File

@@ -84,6 +84,33 @@ extern "C" void ScaLBL_D3Q19_Init(double *dist, int Np)
}
}
extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den)
{
int n;
for (n=0; n<Np; n++){
dist[n] = Den - 0.6666666666666667;
dist[Np+n] = 0.055555555555555555; //double(100*n)+1.f;
dist[2*Np+n] = 0.055555555555555555; //double(100*n)+2.f;
dist[3*Np+n] = 0.055555555555555555; //double(100*n)+3.f;
dist[4*Np+n] = 0.055555555555555555; //double(100*n)+4.f;
dist[5*Np+n] = 0.055555555555555555; //double(100*n)+5.f;
dist[6*Np+n] = 0.055555555555555555; //double(100*n)+6.f;
dist[7*Np+n] = 0.0277777777777778; //double(100*n)+7.f;
dist[8*Np+n] = 0.0277777777777778; //double(100*n)+8.f;
dist[9*Np+n] = 0.0277777777777778; //double(100*n)+9.f;
dist[10*Np+n] = 0.0277777777777778; //double(100*n)+10.f;
dist[11*Np+n] = 0.0277777777777778; //double(100*n)+11.f;
dist[12*Np+n] = 0.0277777777777778; //double(100*n)+12.f;
dist[13*Np+n] = 0.0277777777777778; //double(100*n)+13.f;
dist[14*Np+n] = 0.0277777777777778; //double(100*n)+14.f;
dist[15*Np+n] = 0.0277777777777778; //double(100*n)+15.f;
dist[16*Np+n] = 0.0277777777777778; //double(100*n)+16.f;
dist[17*Np+n] = 0.0277777777777778; //double(100*n)+17.f;
dist[18*Np+n] = 0.0277777777777778; //double(100*n)+18.f;
}
}
//*************************************************************************
extern "C" void ScaLBL_D3Q19_Swap(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz)
{

40
gpu/D3Q19.cu
View File

@@ -267,6 +267,36 @@ __global__ void dvc_ScaLBL_D3Q19_Init(double *dist, int Np)
}
}
__global__ void dvc_ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den)
{
int n;
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 ){
dist[n] = Den - 0.6666666666666667;
dist[Np+n] = 0.055555555555555555; //double(100*n)+1.f;
dist[2*Np+n] = 0.055555555555555555; //double(100*n)+2.f;
dist[3*Np+n] = 0.055555555555555555; //double(100*n)+3.f;
dist[4*Np+n] = 0.055555555555555555; //double(100*n)+4.f;
dist[5*Np+n] = 0.055555555555555555; //double(100*n)+5.f;
dist[6*Np+n] = 0.055555555555555555; //double(100*n)+6.f;
dist[7*Np+n] = 0.0277777777777778; //double(100*n)+7.f;
dist[8*Np+n] = 0.0277777777777778; //double(100*n)+8.f;
dist[9*Np+n] = 0.0277777777777778; //double(100*n)+9.f;
dist[10*Np+n] = 0.0277777777777778; //double(100*n)+10.f;
dist[11*Np+n] = 0.0277777777777778; //double(100*n)+11.f;
dist[12*Np+n] = 0.0277777777777778; //double(100*n)+12.f;
dist[13*Np+n] = 0.0277777777777778; //double(100*n)+13.f;
dist[14*Np+n] = 0.0277777777777778; //double(100*n)+14.f;
dist[15*Np+n] = 0.0277777777777778; //double(100*n)+15.f;
dist[16*Np+n] = 0.0277777777777778; //double(100*n)+16.f;
dist[17*Np+n] = 0.0277777777777778; //double(100*n)+17.f;
dist[18*Np+n] = 0.0277777777777778; //double(100*n)+18.f;
}
}
}
//*************************************************************************
__global__ void dvc_ScaLBL_D3Q19_Swap_Compact(int *neighborList, double *disteven, double *distodd, int Np, int q){
@@ -2325,7 +2355,15 @@ extern "C" void ScaLBL_D3Q19_Init(double *dist, int Np){
dvc_ScaLBL_D3Q19_Init<<<NBLOCKS,NTHREADS >>>(dist, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AA_Init: %s \n",cudaGetErrorString(err));
printf("CUDA error in ScaLBL_D3Q19_Init: %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den){
dvc_ScaLBL_D3Q19_GreyIMRT_Init<<<NBLOCKS,NTHREADS >>>(dist, Np, Den);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_GreyIMRT_Init: %s \n",cudaGetErrorString(err));
}
}

76
models/GreyscaleModel.cpp
View File

@@ -43,6 +43,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
// ---------------------- Greyscale Model parameters -----------------------//
if (greyscale_db->keyExists( "timestepMax" )){
@@ -77,6 +78,10 @@ void ScaLBL_GreyscaleModel::ReadParams(string filename){
if (greyscale_db->keyExists( "tolerance" )){
tolerance = greyscale_db->getScalar<double>( "tolerance" );
}
auto collision = greyscale_db->getWithDefault<std::string>( "collision", "IMRT" );
if (collision == "BGK"){
CollisionType=2;
}
// ------------------------------------------------------------------------//
//------------------------ Other Domain parameters ------------------------//
@@ -374,7 +379,20 @@ void ScaLBL_GreyscaleModel::Initialize(){
//TODO: for BGK, you need to consider voxel porosity
// for IMRT, the whole set of feq is different
// if in the future you have different collison mode, need to write two set of initialization functions
ScaLBL_D3Q19_Init(fq, Np);
if (CollisionType==1){
ScaLBL_D3Q19_GreyIMRT_Init(fq, Np, Den);
if (rank==0) printf("Collision model: Incompressible MRT.\n");
}
else if (CollisionType==2){
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf("Collision model: BGK.\n");
}
else{
if (rank==0) printf("Unknown collison type! IMRT collision is used.\n");
ScaLBL_D3Q19_GreyIMRT_Init(fq, Np, Den);
CollisionType=1;
greyscale_db->putScalar<std::string>( "collision", "IMRT" );
}
if (Restart == true){
if (rank==0){
@@ -442,8 +460,17 @@ void ScaLBL_GreyscaleModel::Run(){
// *************ODD TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
//ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
switch (CollisionType){
case 1:
ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
case 2:
ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
break;
default:
ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
}
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
// Set BCs
@@ -451,25 +478,52 @@ void ScaLBL_GreyscaleModel::Run(){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
//ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
switch (CollisionType){
case 1:
ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
case 2:
ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
break;
default:
ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// *************EVEN TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
//ScaLBL_D3Q19_AAeven_Greyscale(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
switch (CollisionType){
case 1:
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
case 2:
ScaLBL_D3Q19_AAeven_Greyscale(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
break;
default:
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
}
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
// Set BCs
if (BoundaryCondition == 3){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
//ScaLBL_D3Q19_AAeven_Greyscale(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
switch (CollisionType){
case 1:
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
case 2:
ScaLBL_D3Q19_AAeven_Greyscale(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
break;
default:
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/
if (timestep%analysis_interval==0){

1
models/GreyscaleModel.h
View File

@@ -35,6 +35,7 @@ public:
bool Restart,pBC;
int timestep,timestepMax;
int BoundaryCondition;
int CollisionType;
double tau;
double Den;//constant density
double tolerance;