Adding permeability simulator

This commit is contained in:
James E McClure 2015-06-15 21:37:07 -04:00
parent 86e3d88d33
commit d28d4c7377
6 changed files with 212 additions and 200 deletions

View File

@ -112,97 +112,6 @@ extern "C" void InitDenColorDistance(char *ID, double *Den, double *Phi, double
}
extern "C" void Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz, int S)
{
int n,N;
// distributions
double f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double vx,vy,vz;
N = Nx*Ny*Nz;
for (n=0; n<N; n++){
if (ID[n] > 0){
//........................................................................
// Registers to store the distributions
//........................................................................
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//........................................................................
f1 = distodd[n];
f3 = distodd[1*N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//.................Compute the velocity...................................
vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
//..................Write the velocity.....................................
vel[n] = vx;
vel[N+n] = vy;
vel[2*N+n] = vz;
//........................................................................
}
}
}
extern "C" void ComputePressureD3Q19(char *ID, double *disteven, double *distodd, double *Pressure,
int Nx, int Ny, int Nz)
{
int n,N;
// distributions
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
N = Nx*Ny*Nz;
for (n=0; n<N; n++){
if (ID[n] > 0){
//........................................................................
// Registers to store the distributions
//........................................................................
f0 = disteven[n];
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//........................................................................
f1 = distodd[n];
f3 = distodd[1*N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//.................Compute the velocity...................................
Pressure[n] = 0.3333333333333333*(f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+
f9+f12+f11+f14+f13+f16+f15+f18+f17);
}
}
}
//*************************************************************************
extern "C" void ColorBC_inlet(double *Phi, double *Den, double *A_even, double *A_odd,

View File

@ -226,3 +226,95 @@ extern "C" void SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, i
}
}
}
extern "C" void ComputeVelocityD319(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz, int S)
{
int n,N;
// distributions
double f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double vx,vy,vz;
N = Nx*Ny*Nz;
for (n=0; n<N; n++){
if (ID[n] > 0){
//........................................................................
// Registers to store the distributions
//........................................................................
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//........................................................................
f1 = distodd[n];
f3 = distodd[1*N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//.................Compute the velocity...................................
vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
//..................Write the velocity.....................................
vel[n] = vx;
vel[N+n] = vy;
vel[2*N+n] = vz;
//........................................................................
}
}
}
extern "C" void ComputePressureD3Q19(char *ID, double *disteven, double *distodd, double *Pressure,
int Nx, int Ny, int Nz)
{
int n,N;
// distributions
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
N = Nx*Ny*Nz;
for (n=0; n<N; n++){
if (ID[n] > 0){
//........................................................................
// Registers to store the distributions
//........................................................................
f0 = disteven[n];
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//........................................................................
f1 = distodd[n];
f3 = distodd[1*N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//.................Compute the velocity...................................
Pressure[n] = 0.3333333333333333*(f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+
f9+f12+f11+f14+f13+f16+f15+f18+f17);
}
}
}

11
cpu/ScaLBL.h Normal file
View File

@ -0,0 +1,11 @@
extern "C" void PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N);
extern "C" void UnpackDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count,double *recvbuf, double *dist, int Nx, int Ny, int Nz);
extern "C" void InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz);
extern "C" void SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz);
extern "C" void MRT(char *ID, double *f_even, double *f_odd, double rlxA, double rlxB, double Fx, double Fy, double Fz,int Nx, int Ny, int Nz);
extern "C" void InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz);

View File

@ -120,104 +120,6 @@ __global__ void dvc_InitDenColorDistance(char *ID, double *Den, double *Phi, do
}
}
}
__global__ void dvc_Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz)
{
int n,N;
// distributions
double f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double vx,vy,vz;
N = Nx*Ny*Nz;
int S = N/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<N && ID[n] > 0){
//........................................................................
// Registers to store the distributions
//........................................................................
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//........................................................................
f1 = distodd[n];
f3 = distodd[1*N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//.................Compute the velocity...................................
vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
//..................Write the velocity.....................................
vel[n] = vx;
vel[N+n] = vy;
vel[2*N+n] = vz;
//........................................................................
}
}
}
__global__ void dvc_ComputePressureD3Q19(char *ID, double *disteven, double *distodd, double *Pressure,
int Nx, int Ny, int Nz)
{
int n,N;
// distributions
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
N = Nx*Ny*Nz;
int S = N/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<N && ID[n] > 0){
//........................................................................
// Registers to store the distributions
//........................................................................
f0 = disteven[n];
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//........................................................................
f1 = distodd[n];
f3 = distodd[1*N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//.................Compute the velocity...................................
Pressure[n] = 0.3333333333333333*(f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+
f9+f12+f11+f14+f13+f16+f15+f18+f17);
}
}
}
//*************************************************************************
__global__ void dvc_ColorBC_inlet(double *Phi, double *Den, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz)
@ -1469,14 +1371,6 @@ extern "C" void InitDenColorDistance(char *ID, double *Den, double *Phi, double
dvc_InitDenColorDistance<<<NBLOCKS,NTHREADS >>>(ID, Den, Phi, Distance, das, dbs, beta, xp, Nx, Ny, Nz);
}
extern "C" void Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz){
dvc_Compute_VELOCITY<<<NBLOCKS,NTHREADS >>>(ID, disteven, distodd, vel, Nx, Ny, Nz);
}
extern "C" void ComputePressureD3Q19(char *ID, double *disteven, double *distodd, double *Pressure,
int Nx, int Ny, int Nz){
dvc_ComputePressureD3Q19<<< NBLOCKS,NTHREADS >>>(ID, disteven, distodd, Pressure, Nx, Ny, Nz);
}
extern "C" void ComputeColorGradient(char *ID, double *phi, double *ColorGrad, int Nx, int Ny, int Nz){
dvc_ComputeColorGradient<<<NBLOCKS,NTHREADS >>>(ID, phi, ColorGrad, Nx, Ny, Nz);
}

View File

@ -234,6 +234,106 @@ __global__ void dvc_SwapD3Q19(char *ID, double *disteven, double *distodd, int
}
}
__global__ void dvc_ComputeVelocityD3Q19(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz)
{
int n,N;
// distributions
double f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double vx,vy,vz;
N = Nx*Ny*Nz;
int S = N/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<N && ID[n] > 0){
//........................................................................
// Registers to store the distributions
//........................................................................
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//........................................................................
f1 = distodd[n];
f3 = distodd[1*N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//.................Compute the velocity...................................
vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
//..................Write the velocity.....................................
vel[n] = vx;
vel[N+n] = vy;
vel[2*N+n] = vz;
//........................................................................
}
}
}
__global__ void dvc_ComputePressureD3Q19(char *ID, double *disteven, double *distodd, double *Pressure,
int Nx, int Ny, int Nz)
{
int n,N;
// distributions
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
N = Nx*Ny*Nz;
int S = N/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<N && ID[n] > 0){
//........................................................................
// Registers to store the distributions
//........................................................................
f0 = disteven[n];
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//........................................................................
f1 = distodd[n];
f3 = distodd[1*N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//.................Compute the velocity...................................
Pressure[n] = 0.3333333333333333*(f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+
f9+f12+f11+f14+f13+f16+f15+f18+f17);
}
}
}
extern "C" void PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N){
int GRID = count / 512 + 1;
dvc_PackDist <<<GRID,512 >>>(q, list, start, count, sendbuf, dist, N);
@ -251,4 +351,12 @@ extern "C" void SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, i
dvc_SwapD3Q19<<<NBLOCKS,NTHREADS >>>(ID, disteven, distodd, Nx, Ny, Nz);
}
extern "C" void ComputeVelocityD3Q19(char *ID, double *disteven, double *distodd, double *vel, int Nx, int \
Ny, int Nz){
dvc_ComputeVelocityD3Q19<<<NBLOCKS,NTHREADS >>>(ID, disteven, distodd, vel, Nx, Ny, Nz);
}
extern "C" void ComputePressureD3Q19(char *ID, double *disteven, double *distodd, double *Pressure,
int Nx, int Ny, int Nz){
dvc_ComputePressureD3Q19<<< NBLOCKS,NTHREADS >>>(ID, disteven, distodd, Pressure, Nx, Ny, Nz);
}

View File

@ -955,8 +955,6 @@ int main(int argc, char **argv)
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
}
int timestep = 0;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
@ -1186,8 +1184,8 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
ComputeVelocityD3Q19(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
MPI_Barrier(MPI_COMM_WORLD);
}
}
//************************************************************************/