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

add the greyscale effective viscosity back, but by default it is set equal to the normal viscosity

Browse Source
This commit is contained in:
Rex Zhe Li
2020-02-01 14:04:39 -05:00
parent ea8fceda8c
commit 50e4b5a9ba
5 changed files with 77 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
common/ScaLBL.h
View File

@@ -60,16 +60,16 @@ extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int star
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,
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(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 *Pressure);
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz,
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(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 *Pressure);
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz,
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(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);
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz,
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);

48
cpu/Greyscale.cpp
View File

@@ -1,6 +1,6 @@
#include <math.h>
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finish, int Np, double rlx, double Gx, double Gy, double Gz,
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity, double *Pressure){
int n;
// conserved momemnts
@@ -14,7 +14,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finis
double porosity;
double perm;//voxel permeability
double c0, c1; //Guo's model parameters
double mu = (1.0/rlx-0.5)/3.0;//kinematic viscosity
double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
for (int n=start; n<finish; n++){
@@ -42,7 +42,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finis
porosity = Poros[n];
perm = Perm[n];
c0 = 0.5*(1.0+porosity*0.5*mu/perm);
c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
c1 = porosity*0.5*GeoFun/sqrt(perm);
@@ -60,9 +60,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finis
u_mag=sqrt(ux*ux+uy*uy+uz*uz);
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = -porosity*mu/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx;
Fy = -porosity*mu/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy;
Fz = -porosity*mu/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz;
Fx = -porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx;
Fy = -porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy;
Fz = -porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz;
if (porosity==1.0){
Fx=Gx;
Fy=Gy;
@@ -166,7 +166,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finis
}
}
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Gx, double Gy, double Gz,
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity,double *Pressure){
int n;
// conserved momemnts
@@ -181,7 +181,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, in
double porosity;
double perm;//voxel permeability
double c0, c1; //Guo's model parameters
double mu = (1.0/rlx-0.5)/3.0;//kinematic viscosity
double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
int nread;
@@ -264,7 +264,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, in
porosity = Poros[n];
perm = Perm[n];
c0 = 0.5*(1.0+porosity*0.5*mu/perm);
c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
c1 = porosity*0.5*GeoFun/sqrt(perm);
@@ -282,9 +282,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, in
u_mag=sqrt(ux*ux+uy*uy+uz*uz);
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = -porosity*mu/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx;
Fy = -porosity*mu/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy;
Fz = -porosity*mu/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz;
Fx = -porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx;
Fy = -porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy;
Fz = -porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz;
if (porosity==1.0){
Fx=Gx;
Fy=Gy;
@@ -389,7 +389,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, in
}
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int finish, int Np, double rlx, double Gx, double Gy, double Gz,
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity, double Den,double *Pressure){
int n;
double vx,vy,vz,v_mag;
@@ -405,7 +405,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int
double porosity;
double perm;//voxel permeability
double c0, c1; //Guo's model parameters
double mu = (1.0/rlx-0.5)/3.0;//kinematic viscosity
double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
double rlx_setA = rlx;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
@@ -702,7 +702,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int
porosity = Poros[n];
perm = Perm[n];
c0 = 0.5*(1.0+porosity*0.5*mu/perm);
c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
c1 = porosity*0.5*GeoFun/sqrt(perm);
@@ -718,9 +718,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int
u_mag=sqrt(ux*ux+uy*uy+uz*uz);
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = Den*(-porosity*mu/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
Fy = Den*(-porosity*mu/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
Fz = Den*(-porosity*mu/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
Fx = Den*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
Fy = Den*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
Fz = Den*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
if (porosity==1.0){
Fx=Den*Gx;
Fy=Den*Gy;
@@ -850,7 +850,7 @@ 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 Gx, double Gy, double Gz,
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity, double Den,double *Pressure){
int n, nread;
double vx,vy,vz,v_mag;
@@ -866,7 +866,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
double porosity;
double perm;//voxel permeability
double c0, c1; //Guo's model parameters
double mu = (1.0/rlx-0.5)/3.0;//kinematic viscosity
double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
double rlx_setA = rlx;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
@@ -1181,7 +1181,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
porosity = Poros[n];
perm = Perm[n];
c0 = 0.5*(1.0+porosity*0.5*mu/perm);
c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
c1 = porosity*0.5*GeoFun/sqrt(perm);
@@ -1197,9 +1197,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
u_mag=sqrt(ux*ux+uy*uy+uz*uz);
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = Den*(-porosity*mu/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
Fy = Den*(-porosity*mu/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
Fz = Den*(-porosity*mu/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
Fx = Den*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
Fy = Den*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
Fz = Den*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
if (porosity==1.0){
Fx=Den*Gx;
Fy=Den*Gy;

64
gpu/Greyscale.cu
View File

@@ -3,7 +3,7 @@
#define NBLOCKS 1024
#define NTHREADS 256
__global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finish, int Np, double rlx, double Gx, double Gy, double Gz,
__global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity, double *Pressure){
int n;
// conserved momemnts
@@ -17,7 +17,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int f
double porosity;
double perm;//voxel permeability
double c0, c1; //Guo's model parameters
double mu = (1.0/rlx-0.5)/3.0;//kinematic viscosity
double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
int S = Np/NBLOCKS/NTHREADS + 1;
@@ -50,7 +50,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int f
porosity = Poros[n];
perm = Perm[n];
c0 = 0.5*(1.0+porosity*0.5*mu/perm);
c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
c1 = porosity*0.5*GeoFun/sqrt(perm);
@@ -68,9 +68,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int f
u_mag=sqrt(ux*ux+uy*uy+uz*uz);
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = -porosity*mu/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx;
Fy = -porosity*mu/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy;
Fz = -porosity*mu/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz;
Fx = -porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx;
Fy = -porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy;
Fz = -porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz;
if (porosity==1.0){
Fx=Gx;
Fy=Gy;
@@ -176,7 +176,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int f
}
}
__global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Gx, double Gy, double Gz,
__global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity, double *Pressure){
int n;
// conserved momemnts
@@ -191,7 +191,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist
double porosity;
double perm;//voxel permeability
double c0, c1; //Guo's model parameters
double mu = (1.0/rlx-0.5)/3.0;//kinematic viscosity
double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
int S = Np/NBLOCKS/NTHREADS + 1;
@@ -277,7 +277,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist
porosity = Poros[n];
perm = Perm[n];
c0 = 0.5*(1.0+porosity*0.5*mu/perm);
c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
c1 = porosity*0.5*GeoFun/sqrt(perm);
@@ -295,9 +295,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist
u_mag=sqrt(ux*ux+uy*uy+uz*uz);
//Update the body force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = -porosity*mu/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx;
Fy = -porosity*mu/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy;
Fz = -porosity*mu/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz;
Fx = -porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx;
Fy = -porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy;
Fz = -porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz;
if (porosity==1.0){
Fx=Gx;
Fy=Gy;
@@ -402,7 +402,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist
}
}
__global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int finish, int Np, double rlx, double Gx, double Gy, double Gz,
__global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity, double Den, double *Pressure){
int n;
@@ -419,7 +419,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start,
double porosity;
double perm;//voxel permeability
double c0, c1; //Guo's model parameters
double mu = (1.0/rlx-0.5)/3.0;//kinematic viscosity
double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
double rlx_setA = rlx;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
@@ -722,7 +722,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start,
porosity = Poros[n];
perm = Perm[n];
c0 = 0.5*(1.0+porosity*0.5*mu/perm);
c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
c1 = porosity*0.5*GeoFun/sqrt(perm);
@@ -738,9 +738,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start,
u_mag=sqrt(ux*ux+uy*uy+uz*uz);
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = Den*(-porosity*mu/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
Fy = Den*(-porosity*mu/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
Fz = Den*(-porosity*mu/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
Fx = Den*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
Fy = Den*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
Fz = Den*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
if (porosity==1.0){
Fx=Den*Gx;
Fy=Den*Gy;
@@ -873,7 +873,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start,
}
__global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Gx, double Gy, double Gz,
__global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity,double Den, double *Pressure){
int n, nread;
@@ -890,7 +890,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double
double porosity;
double perm;//voxel permeability
double c0, c1; //Guo's model parameters
double mu = (1.0/rlx-0.5)/3.0;//kinematic viscosity
double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
double rlx_setA = rlx;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
@@ -1210,7 +1210,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double
porosity = Poros[n];
perm = Perm[n];
c0 = 0.5*(1.0+porosity*0.5*mu/perm);
c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
c1 = porosity*0.5*GeoFun/sqrt(perm);
@@ -1226,9 +1226,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double
u_mag=sqrt(ux*ux+uy*uy+uz*uz);
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = Den*(-porosity*mu/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
Fy = Den*(-porosity*mu/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
Fz = Den*(-porosity*mu/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
Fx = Den*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
Fy = Den*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
Fz = Den*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
if (porosity==1.0){
Fx=Den*Gx;
Fy=Den*Gy;
@@ -1379,9 +1379,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double
}
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){
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(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 *Pressure){
dvc_ScaLBL_D3Q19_AAeven_Greyscale<<<NBLOCKS,NTHREADS >>>(dist,start,finish,Np,rlx,Fx,Fy,Fz,Poros,Perm,Velocity,Pressure);
dvc_ScaLBL_D3Q19_AAeven_Greyscale<<<NBLOCKS,NTHREADS >>>(dist,start,finish,Np,rlx,rlx_eff,Fx,Fy,Fz,Poros,Perm,Velocity,Pressure);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
@@ -1389,9 +1389,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finis
}
}
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz,double *Poros,double *Perm, double *Velocity,double *Pressure){
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(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 *Pressure){
dvc_ScaLBL_D3Q19_AAodd_Greyscale<<<NBLOCKS,NTHREADS >>>(neighborList,dist,start,finish,Np,rlx,Fx,Fy,Fz,Poros,Perm,Velocity,Pressure);
dvc_ScaLBL_D3Q19_AAodd_Greyscale<<<NBLOCKS,NTHREADS >>>(neighborList,dist,start,finish,Np,rlx,rlx_eff,Fx,Fy,Fz,Poros,Perm,Velocity,Pressure);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
@@ -1399,9 +1399,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, in
}
}
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz,double *Poros,double *Perm, double *Velocity,double Den,double *Pressure){
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(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){
dvc_ScaLBL_D3Q19_AAeven_Greyscale_IMRT<<<NBLOCKS,NTHREADS >>>(dist,start,finish,Np,rlx,Fx,Fy,Fz,Poros,Perm,Velocity,Den,Pressure);
dvc_ScaLBL_D3Q19_AAeven_Greyscale_IMRT<<<NBLOCKS,NTHREADS >>>(dist,start,finish,Np,rlx,rlx_eff,Fx,Fy,Fz,Poros,Perm,Velocity,Den,Pressure);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
@@ -1409,9 +1409,9 @@ 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 Fx, double Fy, double Fz,double *Poros,double *Perm, double *Velocity,double Den,double *Pressure){
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){
dvc_ScaLBL_D3Q19_AAodd_Greyscale_IMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,start,finish,Np,rlx,Fx,Fy,Fz,Poros,Perm,Velocity,Den,Pressure);
dvc_ScaLBL_D3Q19_AAodd_Greyscale_IMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,start,finish,Np,rlx,rlx_eff,Fx,Fy,Fz,Poros,Perm,Velocity,Den,Pressure);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){

29
models/GreyscaleModel.cpp
View File

@@ -14,7 +14,7 @@ void DeleteArray( const TYPE *p )
}
ScaLBL_GreyscaleModel::ScaLBL_GreyscaleModel(int RANK, int NP, MPI_Comm COMM):
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0),Den(0),Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),GreyPorosity(0),
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0),tau_eff(0),Den(0),Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),GreyPorosity(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)
{
SignDist.resize(Nx,Ny,Nz);
@@ -36,6 +36,7 @@ void ScaLBL_GreyscaleModel::ReadParams(string filename){
// set defaults
timestepMax = 100000;
tau = 1.0;
tau_eff = tau;
Den = 1.0;//constant density
tolerance = 0.01;
Fx = Fy = Fz = 0.0;
@@ -52,6 +53,7 @@ void ScaLBL_GreyscaleModel::ReadParams(string filename){
if (greyscale_db->keyExists( "tau" )){
tau = greyscale_db->getScalar<double>( "tau" );
}
tau_eff = greyscale_db->getWithDefault<double>( "tau_eff", tau );
if (greyscale_db->keyExists( "Den" )){
Den = greyscale_db->getScalar<double>( "Den" );
}
@@ -453,6 +455,7 @@ void ScaLBL_GreyscaleModel::Run(){
PROFILE_START("Loop");
auto current_db = db->cloneDatabase();
double rlx = 1.0/tau;
double rlx_eff = 1.0/tau_eff;
double error = 1.0;
double flow_rate_previous = 0.0;
while (timestep < timestepMax && error > tolerance) {
@@ -462,13 +465,13 @@ void ScaLBL_GreyscaleModel::Run(){
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
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);
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;
case 2:
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(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, 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);
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;
}
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
@@ -480,13 +483,13 @@ void ScaLBL_GreyscaleModel::Run(){
}
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);
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;
case 2:
ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, 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);
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;
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
@@ -496,13 +499,13 @@ void ScaLBL_GreyscaleModel::Run(){
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
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);
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;
case 2:
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(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, 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);
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;
}
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
@@ -514,13 +517,13 @@ void ScaLBL_GreyscaleModel::Run(){
}
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);
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, 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);
ScaLBL_D3Q19_AAeven_Greyscale(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, 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);
ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
break;
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);

1
models/GreyscaleModel.h
View File

@@ -37,6 +37,7 @@ public:
int BoundaryCondition;
int CollisionType;
double tau;
double tau_eff;
double Den;//constant density
double tolerance;
double Fx,Fy,Fz,flux;