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

Merge branch 'VariableDensity'

Browse Source
This commit is contained in:
James E McClure
2017-12-01 11:50:27 -05:00
parent ff0dd164da f3c0cd3aea
commit 1b8ceef7d8
6 changed files with 290 additions and 113 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
common/ScaLBL.h
View File

@@ -77,7 +77,7 @@ extern "C" void ScaLBL_D3Q19_ColorCollide( char *ID, double *disteven, double *d
double alpha, double beta, double Fx, double Fy, double Fz); double alpha, double beta, double Fx, double Fy, double Fz);
extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad, extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad,
double *Velocity, int Nx, int Ny, int Nz, double tau1, double tau2, double *Velocity, int Nx, int Ny, int Nz, double tau1, double tau2, double rho1, double rho2,
double alpha, double beta, double Fx, double Fy, double Fz); double alpha, double beta, double Fx, double Fy, double Fz);
extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A_odd, double *B_even, double *B_odd, extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A_odd, double *B_even, double *B_odd,
@@ -86,10 +86,10 @@ extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A
extern "C" void ScaLBL_ComputePhaseField(char *ID, double *Phi, double *Den, int N); extern "C" void ScaLBL_ComputePhaseField(char *ID, double *Phi, double *Den, int N);
extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, double *A_odd, extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *Vel, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz); double *B_even, double *B_odd, int Nx, int Ny, int Nz);
extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, double *A_odd, extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *Vel, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz); double *B_even, double *B_odd, int Nx, int Ny, int Nz);
extern "C" void ScaLBL_D3Q19_Velocity_BC_z(double *disteven, double *distodd, double uz, extern "C" void ScaLBL_D3Q19_Velocity_BC_z(double *disteven, double *distodd, double uz,

156
cpu/Color.cpp
View File

@@ -103,10 +103,13 @@ extern "C" void ScaLBL_Color_InitDistance(char *ID, double *Den, double *Phi, do
//************************************************************************* //*************************************************************************
extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, double *A_odd, extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *Velocity, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz) double *B_even, double *B_odd, int Nx, int Ny, int Nz)
{ {
int i,j,k,n,N; int i,j,k,n,N;
double na,nb,ux,uy,uz;
double a1,b1,a2,b2;
N = Nx*Ny*Nz; N = Nx*Ny*Nz;
// Fill the inlet with component a // Fill the inlet with component a
for (k=0; k<1; k++){ for (k=0; k<1; k++){
@@ -126,6 +129,7 @@ extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, doub
Den[n] = 1.0; Den[n] = 1.0;
Den[N+n] = 0.0; Den[N+n] = 0.0;
/*
A_even[n] = 0.3333333333333333; A_even[n] = 0.3333333333333333;
A_odd[n] = 0.1111111111111111; A_odd[n] = 0.1111111111111111;
A_even[N+n] = 0.1111111111111111; A_even[N+n] = 0.1111111111111111;
@@ -135,21 +139,81 @@ extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, doub
A_even[3*N+n] = 0.1111111111111111; A_even[3*N+n] = 0.1111111111111111;
B_even[n] = 0.0; B_even[n] = 0.0;
B_odd[n] = 0.0; B_odd[n] = 0.0;
B_even[N+n] = 0.0; B_even[N+n] = 0.0;
B_odd[N+n] = 0.0; B_odd[N+n] = 0.0;
B_even[2*N+n] = 0.0; B_even[2*N+n] = 0.0;
B_odd[2*N+n] = 0.0; B_odd[2*N+n] = 0.0;
B_even[3*N+n] = 0.0; B_even[3*N+n] = 0.0;
*/
na = 1.0;
nb = 0.0;
//....Load the flow velocity...........
ux = Velocity[n];
uy = Velocity[N+n];
uz = Velocity[2*N+n];
A_even[n] = 0.3333333333333333*na;
B_even[n] = 0.3333333333333333*nb;
// Non-Stationary equilibrium distributions
//feq[0] = 0.1111111111111111*(1+4.5*ux);
//feq[1] = 0.1111111111111111*(1-4.5*ux);
//feq[2] = 0.1111111111111111*(1+4.5*uy);
//feq[3] = 0.1111111111111111*(1-4.5*uy);
//feq[4] = 0.1111111111111111*(1+4.5*uz);
//feq[5] = 0.1111111111111111*(1-4.5*uz);
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
a1 = na*(0.1111111111111111*(1+4.5*ux));
b1 = nb*(0.1111111111111111*(1+4.5*ux));
a2 = na*(0.1111111111111111*(1-4.5*ux));
b2 = nb*(0.1111111111111111*(1-4.5*ux));
A_odd[n] = a1;
A_even[N+n] = a2;
B_odd[n] = b1;
B_even[N+n] = b2;
//...............................................
// q = 2
// Cq = {0,1,0}
a1 = na*(0.1111111111111111*(1+4.5*uy));
b1 = nb*(0.1111111111111111*(1+4.5*uy));
a2 = na*(0.1111111111111111*(1-4.5*uy));
b2 = nb*(0.1111111111111111*(1-4.5*uy));
A_odd[N+n] = a1;
A_even[2*N+n] = a2;
B_odd[N+n] = b1;
B_even[2*N+n] = b2;
//...............................................
// q = 4
// Cq = {0,0,1}
a1 = na*(0.1111111111111111*(1+4.5*uz));
b1 = nb*(0.1111111111111111*(1+4.5*uz));
a2 = na*(0.1111111111111111*(1-4.5*uz));
b2 = nb*(0.1111111111111111*(1-4.5*uz));
A_odd[2*N+n] = a1;
A_even[3*N+n] = a2;
B_odd[2*N+n] = b1;
B_even[3*N+n] = b2;
} }
} }
} }
} }
//************************************************************************* //*************************************************************************
extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, double *A_odd, extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *Velocity, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz) double *B_even, double *B_odd, int Nx, int Ny, int Nz)
{ {
int i,j,k,n,N; int i,j,k,n,N;
double na,nb,ux,uy,uz;
double a1,b1,a2,b2;
N = Nx*Ny*Nz; N = Nx*Ny*Nz;
// Fill the outlet with component b // Fill the outlet with component b
for (k=Nz-3; k<Nz-1; k++){ for (k=Nz-3; k<Nz-1; k++){
@@ -161,7 +225,7 @@ extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, doub
Den[n] = 0.0; Den[n] = 0.0;
Den[N+n] = 1.0; Den[N+n] = 1.0;
A_even[n] = 0.0; /* A_even[n] = 0.0;
A_odd[n] = 0.0; A_odd[n] = 0.0;
A_even[N+n] = 0.0; A_even[N+n] = 0.0;
A_odd[N+n] = 0.0; A_odd[N+n] = 0.0;
@@ -176,6 +240,61 @@ extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, doub
B_even[2*N+n] = 0.1111111111111111; B_even[2*N+n] = 0.1111111111111111;
B_odd[2*N+n] = 0.1111111111111111; B_odd[2*N+n] = 0.1111111111111111;
B_even[3*N+n] = 0.1111111111111111; B_even[3*N+n] = 0.1111111111111111;
*/
na = 0.0;
nb = 1.0;
//....Load the flow velocity...........
ux = Velocity[n];
uy = Velocity[N+n];
uz = Velocity[2*N+n];
A_even[n] = 0.3333333333333333*na;
B_even[n] = 0.3333333333333333*nb;
// Non-Stationary equilibrium distributions
//feq[0] = 0.1111111111111111*(1+4.5*ux);
//feq[1] = 0.1111111111111111*(1-4.5*ux);
//feq[2] = 0.1111111111111111*(1+4.5*uy);
//feq[3] = 0.1111111111111111*(1-4.5*uy);
//feq[4] = 0.1111111111111111*(1+4.5*uz);
//feq[5] = 0.1111111111111111*(1-4.5*uz);
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
a1 = na*(0.1111111111111111*(1+4.5*ux));
b1 = nb*(0.1111111111111111*(1+4.5*ux));
a2 = na*(0.1111111111111111*(1-4.5*ux));
b2 = nb*(0.1111111111111111*(1-4.5*ux));
A_odd[n] = a1;
A_even[N+n] = a2;
B_odd[n] = b1;
B_even[N+n] = b2;
//...............................................
// q = 2
// Cq = {0,1,0}
a1 = na*(0.1111111111111111*(1+4.5*uy));
b1 = nb*(0.1111111111111111*(1+4.5*uy));
a2 = na*(0.1111111111111111*(1-4.5*uy));
b2 = nb*(0.1111111111111111*(1-4.5*uy));
A_odd[N+n] = a1;
A_even[2*N+n] = a2;
B_odd[N+n] = b1;
B_even[2*N+n] = b2;
//...............................................
// q = 4
// Cq = {0,0,1}
a1 = na*(0.1111111111111111*(1+4.5*uz));
b1 = nb*(0.1111111111111111*(1+4.5*uz));
a2 = na*(0.1111111111111111*(1-4.5*uz));
b2 = nb*(0.1111111111111111*(1-4.5*uz));
A_odd[2*N+n] = a1;
A_even[3*N+n] = a2;
B_odd[2*N+n] = b1;
B_even[3*N+n] = b2;
} }
} }
@@ -955,8 +1074,8 @@ extern "C" void ScaLBL_D3Q19_ColorCollide( char *ID, double *disteven, double *d
} }
extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad, extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad,
double *Velocity, int Nx, int Ny, int Nz, double tau1, double tau2, double *Velocity, int Nx, int Ny, int Nz, double tau1, double tau2,
double alpha, double beta, double Fx, double Fy, double Fz) double rho1, double rho2, double alpha, double beta, double Fx, double Fy, double Fz)
{ {
int i,j,k,n,nn,N; int i,j,k,n,nn,N;
@@ -969,6 +1088,7 @@ extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, doubl
// additional variables needed for computations // additional variables needed for computations
double rho,jx,jy,jz,C,nx,ny,nz; double rho,jx,jy,jz,C,nx,ny,nz;
double tau,rlx_setA,rlx_setB; double tau,rlx_setA,rlx_setB;
double rho0;
N = Nx*Ny*Nz; N = Nx*Ny*Nz;
char id; char id;
@@ -1097,6 +1217,8 @@ extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, doubl
tau=tau1 + 0.5*(1.0-f1)*(tau2-tau1); tau=tau1 + 0.5*(1.0-f1)*(tau2-tau1);
rlx_setA = 1.f/tau; rlx_setA = 1.f/tau;
rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
// local density
rho0=rho1 + 0.5*(1.0-f1)*(rho2-rho1);
//......No color gradient at z-boundary if pressure BC are set............. //......No color gradient at z-boundary if pressure BC are set.............
// if (pBC && k==0) nx = ny = nz = 0.f; // if (pBC && k==0) nx = ny = nz = 0.f;
@@ -1183,18 +1305,18 @@ extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, doubl
//........................................................................ //........................................................................
//..........Toelke, Fruediger et. al. 2006............... //..........Toelke, Fruediger et. al. 2006...............
if (C == 0.0) nx = ny = nz = 0.0; if (C == 0.0) nx = ny = nz = 0.0;
m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) -alpha*C - m1); m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -alpha*C - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho)- m2); m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2);
m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4);
m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6);
m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8);
m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9);
m10 = m10 + rlx_setA*( - m10); m10 = m10 + rlx_setA*( - m10);
m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(ny*ny-nz*nz)- m11);
m12 = m12 + rlx_setA*( - m12); m12 = m12 + rlx_setA*( - m12);
m13 = m13 + rlx_setA*( (jx*jy/rho) + 0.5*alpha*C*nx*ny - m13); m13 = m13 + rlx_setA*( (jx*jy/rho0) + 0.5*alpha*C*nx*ny - m13);
m14 = m14 + rlx_setA*( (jy*jz/rho) + 0.5*alpha*C*ny*nz - m14); m14 = m14 + rlx_setA*( (jy*jz/rho0) + 0.5*alpha*C*ny*nz - m14);
m15 = m15 + rlx_setA*( (jx*jz/rho) + 0.5*alpha*C*nx*nz - m15); m15 = m15 + rlx_setA*( (jx*jz/rho0) + 0.5*alpha*C*nx*nz - m15);
m16 = m16 + rlx_setB*( - m16); m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17); m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18); m18 = m18 + rlx_setB*( - m18);
@@ -1296,9 +1418,9 @@ extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, doubl
distodd[7*N+n] = f15; distodd[7*N+n] = f15;
distodd[8*N+n] = f17; distodd[8*N+n] = f17;
//...Store the Velocity.......................... //...Store the Velocity..........................
Velocity[n] = jx; Velocity[n] = jx/rho0;
Velocity[N+n] = jy; Velocity[N+n] = jy/rho0;
Velocity[2*N+n] = jz; Velocity[2*N+n] = jz/rho0;
//*************************************************************** //***************************************************************
} // check if n is in the solid } // check if n is in the solid
} // loop over n } // loop over n

161
gpu/Color.cu
View File

@@ -122,7 +122,7 @@ __global__ void dvc_ScaLBL_Color_InitDistance(char *ID, double *Den, double *Ph
} }
} }
//************************************************************************* //*************************************************************************
__global__ void dvc_ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, double *A_odd, __global__ void dvc_ScaLBL_Color_BC_z(double *Phi, double *Den, double *Velocity, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz) double *B_even, double *B_odd, int Nx, int Ny, int Nz)
{ {
int ij,k,n,N; int ij,k,n,N;
@@ -142,26 +142,48 @@ __global__ void dvc_ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even,
Den[n] = 1.0; Den[n] = 1.0;
Den[N+n] = 0.0; Den[N+n] = 0.0;
A_even[n] = 0.3333333333333333; na = 1.0;
A_odd[n] = 0.1111111111111111; nb = 0.0;
A_even[N+n] = 0.1111111111111111; ux = Velocity[n];
A_odd[N+n] = 0.1111111111111111; uy = Velocity[N+n];
A_even[2*N+n] = 0.1111111111111111; uz = Velocity[2*N+n];
A_odd[2*N+n] = 0.1111111111111111; A_even[n] = 0.3333333333333333*na;
A_even[3*N+n] = 0.1111111111111111; B_even[n] = 0.3333333333333333*nb;
B_even[n] = 0.0; a1 = na*(0.1111111111111111*(1+4.5*ux));
B_odd[n] = 0.0; b1 = nb*(0.1111111111111111*(1+4.5*ux));
B_even[N+n] = 0.0; a2 = na*(0.1111111111111111*(1-4.5*ux));
B_odd[N+n] = 0.0; b2 = nb*(0.1111111111111111*(1-4.5*ux));
B_even[2*N+n] = 0.0;
B_odd[2*N+n] = 0.0; A_odd[n] = a1;
B_even[3*N+n] = 0.0; A_even[N+n] = a2;
B_odd[n] = b1;
B_even[N+n] = b2;
a1 = na*(0.1111111111111111*(1+4.5*uy));
b1 = nb*(0.1111111111111111*(1+4.5*uy));
a2 = na*(0.1111111111111111*(1-4.5*uy));
b2 = nb*(0.1111111111111111*(1-4.5*uy));
A_odd[N+n] = a1;
A_even[2*N+n] = a2;
B_odd[N+n] = b1;
B_even[2*N+n] = b2;
a1 = na*(0.1111111111111111*(1+4.5*uz));
b1 = nb*(0.1111111111111111*(1+4.5*uz));
a2 = na*(0.1111111111111111*(1-4.5*uz));
b2 = nb*(0.1111111111111111*(1-4.5*uz));
A_odd[2*N+n] = a1;
A_even[3*N+n] = a2;
B_odd[2*N+n] = b1;
B_even[3*N+n] = b2;
} }
} }
} }
//************************************************************************* //*************************************************************************
__global__ void dvc_ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, double *A_odd, __global__ void dvc_ScaLBL_Color_BC_Z(double *Phi, double *Den, double *Velocity, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz) double *B_even, double *B_odd, int Nx, int Ny, int Nz)
{ {
int ij,k,n,N; int ij,k,n,N;
@@ -175,21 +197,43 @@ __global__ void dvc_ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even,
Den[n] = 0.0; Den[n] = 0.0;
Den[N+n] = 1.0; Den[N+n] = 1.0;
A_even[n] = 0.0; na = 0.0;
A_odd[n] = 0.0; nb = -1.0;
A_even[N+n] = 0.0; ux = Velocity[n];
A_odd[N+n] = 0.0; uy = Velocity[N+n];
A_even[2*N+n] = 0.0; uz = Velocity[2*N+n];
A_odd[2*N+n] = 0.0; A_even[n] = 0.3333333333333333*na;
A_even[3*N+n] = 0.0; B_even[n] = 0.3333333333333333*nb;
B_even[n] = 0.3333333333333333; a1 = na*(0.1111111111111111*(1+4.5*ux));
B_odd[n] = 0.1111111111111111; b1 = nb*(0.1111111111111111*(1+4.5*ux));
B_even[N+n] = 0.1111111111111111; a2 = na*(0.1111111111111111*(1-4.5*ux));
B_odd[N+n] = 0.1111111111111111; b2 = nb*(0.1111111111111111*(1-4.5*ux));
B_even[2*N+n] = 0.1111111111111111;
B_odd[2*N+n] = 0.1111111111111111; A_odd[n] = a1;
B_even[3*N+n] = 0.1111111111111111; A_even[N+n] = a2;
B_odd[n] = b1;
B_even[N+n] = b2;
a1 = na*(0.1111111111111111*(1+4.5*uy));
b1 = nb*(0.1111111111111111*(1+4.5*uy));
a2 = na*(0.1111111111111111*(1-4.5*uy));
b2 = nb*(0.1111111111111111*(1-4.5*uy));
A_odd[N+n] = a1;
A_even[2*N+n] = a2;
B_odd[N+n] = b1;
B_even[2*N+n] = b2;
a1 = na*(0.1111111111111111*(1+4.5*uz));
b1 = nb*(0.1111111111111111*(1+4.5*uz));
a2 = na*(0.1111111111111111*(1-4.5*uz));
b2 = nb*(0.1111111111111111*(1-4.5*uz));
A_odd[2*N+n] = a1;
A_even[3*N+n] = a2;
B_odd[2*N+n] = b1;
B_even[3*N+n] = b2;
} }
for (k=Nz-1; k<Nz; k++){ for (k=Nz-1; k<Nz; k++){
@@ -953,8 +997,8 @@ __global__ void dvc_ScaLBL_D3Q19_ColorCollide( char *ID, double *disteven, doub
} }
__global__ void dvc_ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad, __global__ void dvc_ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad,
double *Velocity, int Nx, int Ny, int Nz, double tau1, double tau2, double *Velocity, int Nx, int Ny, int Nz, double tau1, double tau2,
double alpha, double beta, double Fx, double Fy, double Fz) double rho1, double rho2, double alpha, double beta, double Fx, double Fy, double Fz)
{ {
char id; char id;
@@ -1098,6 +1142,8 @@ __global__ void dvc_ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven,
rlx_setA = 1.f/tau; rlx_setA = 1.f/tau;
rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
// local density rho0=rho1 + 0.5*(1.0-f1)*(rho2-rho1);
//......No color gradient at z-boundary if pressure BC are set............. //......No color gradient at z-boundary if pressure BC are set.............
// if (pBC && k==0) nx = ny = nz = 0.f; // if (pBC && k==0) nx = ny = nz = 0.f;
// if (pBC && k==Nz-1) nx = ny = nz = 0.f; // if (pBC && k==Nz-1) nx = ny = nz = 0.f;
@@ -1183,22 +1229,23 @@ __global__ void dvc_ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven,
// PERFORM RELAXATION PROCESS // PERFORM RELAXATION PROCESS
//........................................................................ //........................................................................
//..........Toelke, Fruediger et. al. 2006............... //..........Toelke, Fruediger et. al. 2006...............
if (C == 0.0) nx = ny = nz = 0.0; if (C == 0.0) nx = ny = nz = 0.0;
m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) -alpha*C - m1); m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -alpha*C - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho)- m2); m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2);
m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4);
m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6);
m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8);
m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9);
m10 = m10 + rlx_setA*( - m10); m10 = m10 + rlx_setA*( - m10);
m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(ny*ny-nz*nz)- m11);
m12 = m12 + rlx_setA*( - m12); m12 = m12 + rlx_setA*( - m12);
m13 = m13 + rlx_setA*( (jx*jy/rho) + 0.5*alpha*C*nx*ny - m13); m13 = m13 + rlx_setA*( (jx*jy/rho0) + 0.5*alpha*C*nx*ny - m13);
m14 = m14 + rlx_setA*( (jy*jz/rho) + 0.5*alpha*C*ny*nz - m14); m14 = m14 + rlx_setA*( (jy*jz/rho0) + 0.5*alpha*C*ny*nz - m14);
m15 = m15 + rlx_setA*( (jx*jz/rho) + 0.5*alpha*C*nx*nz - m15); m15 = m15 + rlx_setA*( (jx*jz/rho0) + 0.5*alpha*C*nx*nz - m15);
m16 = m16 + rlx_setB*( - m16); m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17); m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18); m18 = m18 + rlx_setB*( - m18);
//.................inverse transformation...................................................... //.................inverse transformation......................................................
f0 = 0.05263157894736842*rho-0.012531328320802*m1+0.04761904761904762*m2; f0 = 0.05263157894736842*rho-0.012531328320802*m1+0.04761904761904762*m2;
f1 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 f1 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2
@@ -1297,9 +1344,9 @@ __global__ void dvc_ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven,
distodd[7*N+n] = f15; distodd[7*N+n] = f15;
distodd[8*N+n] = f17; distodd[8*N+n] = f17;
//...Store the Velocity.......................... //...Store the Velocity..........................
Velocity[n] = jx; Velocity[n] = jx/rho0;
Velocity[N+n] = jy; Velocity[N+n] = jy/rho0;
Velocity[2*N+n] = jz; Velocity[2*N+n] = jz/rho0;
//*************************************************************** //***************************************************************
}// check if n is in the solid }// check if n is in the solid
@@ -1643,8 +1690,8 @@ extern "C" void ColorCollide( char *ID, double *disteven, double *distodd, doubl
} }
extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad, extern "C" void ScaLBL_D3Q19_ColorCollide_gen( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad,
double *Velocity, int Nx, int Ny, int Nz, double tau1, double tau2, double *Velocity, int Nx, int Ny, int Nz, double tau1, double tau2,
double alpha, double beta, double Fx, double Fy, double Fz){ double rho1, double rho2, double alpha, double beta, double Fx, double Fy, double Fz){
dvc_ScaLBL_D3Q19_ColorCollide_gen<<<NBLOCKS,NTHREADS >>>(ID, disteven, distodd, phi, ColorGrad, Velocity, Nx, Ny, Nz, tau1, tau2, dvc_ScaLBL_D3Q19_ColorCollide_gen<<<NBLOCKS,NTHREADS >>>(ID, disteven, distodd, phi, ColorGrad, Velocity, Nx, Ny, Nz, tau1, tau2,
alpha, beta, Fx, Fy, Fz); alpha, beta, Fx, Fy, Fz);
@@ -1672,14 +1719,14 @@ extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A
dvc_ScaLBL_D3Q7_ColorCollideMass<<<NBLOCKS,NTHREADS >>>(ID, A_even, A_odd, B_even, B_odd, Den, Phi, ColorGrad, Velocity, beta, N, pBC); dvc_ScaLBL_D3Q7_ColorCollideMass<<<NBLOCKS,NTHREADS >>>(ID, A_even, A_odd, B_even, B_odd, Den, Phi, ColorGrad, Velocity, beta, N, pBC);
} }
// Pressure Boundary Conditions Functions // Pressure Boundary Conditions Functions
extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, double *A_odd, extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *Vel, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz){ double *B_even, double *B_odd, int Nx, int Ny, int Nz){
int GRID = Nx*Ny / 512 + 1; int GRID = Nx*Ny / 512 + 1;
dvc_ScaLBL_Color_BC_z<<< GRID,512 >>>(Phi, Den, A_even, A_odd, B_even, B_odd, Nx, Ny, Nz); dvc_ScaLBL_Color_BC_z<<< GRID,512 >>>(Phi, Den, Vel, A_even, A_odd, B_even, B_odd, Nx, Ny, Nz);
} }
extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, double *A_odd, extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *Vel, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz){ double *B_even, double *B_odd, int Nx, int Ny, int Nz){
int GRID = Nx*Ny / 512 + 1; int GRID = Nx*Ny / 512 + 1;
dvc_ScaLBL_Color_BC_Z<<< GRID,512 >>>(Phi, Den, A_even, A_odd, B_even, B_odd, Nx, Ny, Nz); dvc_ScaLBL_Color_BC_Z<<< GRID,512 >>>(Phi, Den, Vel, A_even, A_odd, B_even, B_odd, Nx, Ny, Nz);
} }

8
tests/TestBubble.cpp
View File

@@ -601,12 +601,12 @@ int main(int argc, char **argv)
} }
if (pBC && kproc == 0) { if (pBC && kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
} }
if (pBC && kproc == nprocz-1){ if (pBC && kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
} }
//........................................................................... //...........................................................................
@@ -710,12 +710,12 @@ int main(int argc, char **argv)
// Pressure boundary conditions // Pressure boundary conditions
if (pBC && kproc == 0) { if (pBC && kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
} }
if (pBC && kproc == nprocz-1){ if (pBC && kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
} }
//................................................................................... //...................................................................................

48
tests/lbpm_color_macro_simulator.cpp
View File

@@ -136,7 +136,8 @@ int main(int argc, char **argv)
double D = 1.0; // reference length for non-dimensionalization double D = 1.0; // reference length for non-dimensionalization
// Color Model parameters // Color Model parameters
int timestepMax; int timestepMax;
double tau1, tau2, Fx,Fy,Fz,tol,err; double tau1, tau2, rho1,rho2;
double Fx,Fy,Fz,tol,err;
double alpha, beta; double alpha, beta;
double das, dbs, phi_s; double das, dbs, phi_s;
double din,dout; double din,dout;
@@ -163,8 +164,10 @@ int main(int argc, char **argv)
ifstream input("Color.in"); ifstream input("Color.in");
if (input.is_open()){ if (input.is_open()){
// Line 1: model parameters (tau, alpha, beta, das, dbs) // Line 1: model parameters (tau, alpha, beta, das, dbs)
input >> tau1; // Viscosity parameter input >> tau1; // Viscosity non-wetting
input >> tau2; // Viscosity parameter input >> tau2; // Viscosity wetting
input >> rho1; // density non-wetting
input >> rho2; // density wetting
input >> alpha; // Surface Tension parameter input >> alpha; // Surface Tension parameter
input >> beta; // Width of the interface input >> beta; // Width of the interface
input >> phi_s; // value of phi at the solid surface input >> phi_s; // value of phi at the solid surface
@@ -192,6 +195,7 @@ int main(int argc, char **argv)
// Print warning // Print warning
printf("WARNING: No input file provided (Color.in is missing)! Default parameters will be used. \n"); printf("WARNING: No input file provided (Color.in is missing)! Default parameters will be used. \n");
tau1 = tau2 = 1.0; tau1 = tau2 = 1.0;
rho1 = rho2 = 1.0;
alpha=0.005; alpha=0.005;
beta= 0.9; beta= 0.9;
Fx = Fy = Fz = 0.0; Fx = Fy = Fz = 0.0;
@@ -235,6 +239,8 @@ int main(int argc, char **argv)
//................................................. //.................................................
MPI_Bcast(&tau1,1,MPI_DOUBLE,0,comm); MPI_Bcast(&tau1,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&tau2,1,MPI_DOUBLE,0,comm); MPI_Bcast(&tau2,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&rho1,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&rho2,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm); MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm); MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&das,1,MPI_DOUBLE,0,comm); MPI_Bcast(&das,1,MPI_DOUBLE,0,comm);
@@ -292,6 +298,8 @@ int main(int argc, char **argv)
printf("********************************************************\n"); printf("********************************************************\n");
printf("tau (non-wetting) = %f \n", tau1); printf("tau (non-wetting) = %f \n", tau1);
printf("tau (wetting) = %f \n", tau2); printf("tau (wetting) = %f \n", tau2);
printf("density (non-wetting) = %f \n", rho1);
printf("density (wetting) = %f \n", rho2);
printf("alpha = %f \n", alpha); printf("alpha = %f \n", alpha);
printf("beta = %f \n", beta); printf("beta = %f \n", beta);
printf("das = %f \n", das); printf("das = %f \n", das);
@@ -647,13 +655,13 @@ int main(int argc, char **argv)
} }
if (BoundaryCondition==1 && Mask.kproc == 0) { if (BoundaryCondition==1 && Mask.kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
@@ -663,13 +671,13 @@ int main(int argc, char **argv)
} }
if (BoundaryCondition==2 && Mask.kproc == 0) { if (BoundaryCondition==2 && Mask.kproc == 0) {
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (BoundaryCondition==2 && Mask.kproc == nprocz-1){ if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
@@ -686,13 +694,13 @@ int main(int argc, char **argv)
// set the initial boundary conditions // set the initial boundary conditions
if (Mask.kproc == 0) { if (Mask.kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (Mask.kproc == nprocz-1){ if (Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
} }
@@ -720,7 +728,7 @@ int main(int argc, char **argv)
if (pBC && Dm.kproc == 0){ if (pBC && Dm.kproc == 0){
if (rank==0) printf("Flux = %.3e, Computed inlet pressure: %f \n",flux,din); if (rank==0) printf("Flux = %.3e, Computed inlet pressure: %f \n",flux,din);
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
@@ -728,7 +736,7 @@ int main(int argc, char **argv)
if (pBC && Dm.kproc == nprocz-1){ if (pBC && Dm.kproc == nprocz-1){
// if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout); // if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout);
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
} }
@@ -828,7 +836,7 @@ int main(int argc, char **argv)
// Fused Color Gradient and Collision // Fused Color Gradient and Collision
//************************************************************************* //*************************************************************************
ScaLBL_D3Q19_ColorCollide_gen( ID,f_even,f_odd,Phi,ColorGrad, ScaLBL_D3Q19_ColorCollide_gen( ID,f_even,f_odd,Phi,ColorGrad,
Velocity,Nx,Ny,Nz,tau1,tau2,alpha,beta,Fx,Fy,Fz); Velocity,Nx,Ny,Nz,tau1,tau2,rho1,rho2,alpha,beta,Fx,Fy,Fz);
//************************************************************************* //*************************************************************************
ScaLBL_DeviceBarrier(); ScaLBL_DeviceBarrier();
@@ -899,24 +907,24 @@ int main(int argc, char **argv)
// Pressure boundary conditions // Pressure boundary conditions
if (BoundaryCondition==1 && Mask.kproc == 0) { if (BoundaryCondition==1 && Mask.kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
// Velocity boundary conditions // Velocity boundary conditions
if (BoundaryCondition==2 && Mask.kproc == 0) { if (BoundaryCondition==2 && Mask.kproc == 0) {
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (BoundaryCondition==2 && Mask.kproc == nprocz-1){ if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
@@ -928,12 +936,12 @@ int main(int argc, char **argv)
// set the initial boundary conditions // set the initial boundary conditions
if (Mask.kproc == 0) { if (Mask.kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (Mask.kproc == nprocz-1){ if (Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
} }
@@ -953,14 +961,14 @@ int main(int argc, char **argv)
if (pBC && Dm.kproc == 0){ if (pBC && Dm.kproc == 0){
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (pBC && Dm.kproc == nprocz-1){ if (pBC && Dm.kproc == nprocz-1){
// if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout); // if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout);
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
} }

24
tests/lbpm_color_simulator.cpp
View File

@@ -637,13 +637,13 @@ int main(int argc, char **argv)
} }
if (BoundaryCondition==1 && Mask.kproc == 0) { if (BoundaryCondition==1 && Mask.kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
@@ -653,13 +653,13 @@ int main(int argc, char **argv)
} }
if (BoundaryCondition==2 && Mask.kproc == 0) { if (BoundaryCondition==2 && Mask.kproc == 0) {
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (BoundaryCondition==2 && Mask.kproc == nprocz-1){ if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
// Set dynamic pressure boundary conditions // Set dynamic pressure boundary conditions
@@ -675,12 +675,12 @@ int main(int argc, char **argv)
// set the initial boundary conditions // set the initial boundary conditions
if (Mask.kproc == 0) { if (Mask.kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (Mask.kproc == nprocz-1){ if (Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
} }
@@ -849,24 +849,24 @@ int main(int argc, char **argv)
// Pressure boundary conditions // Pressure boundary conditions
if (BoundaryCondition==1 && Mask.kproc == 0) { if (BoundaryCondition==1 && Mask.kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (BoundaryCondition==1 && Mask.kproc == nprocz-1){ if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
// Velocity boundary conditions // Velocity boundary conditions
if (BoundaryCondition==2 && Mask.kproc == 0) { if (BoundaryCondition==2 && Mask.kproc == 0) {
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,0);
} }
if (BoundaryCondition==2 && Mask.kproc == nprocz-1){ if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1);
} }
@@ -878,12 +878,12 @@ int main(int argc, char **argv)
// set the initial boundary conditions // set the initial boundary conditions
if (Mask.kproc == 0) { if (Mask.kproc == 0) {
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz); ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
} }
if (Mask.kproc == nprocz-1){ if (Mask.kproc == nprocz-1){
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2)); ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz); ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
} }
} }