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

add hip version of Lee model

Browse Source
This commit is contained in:
James McClure 2021-03-24 07:57:27 -04:00
parent 32d750a04c
commit f44167d2ef
1 changed files with 794 additions and 689 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

457
hip/FreeLee.cu
View File

@ -1,11 +1,12 @@
#include <math.h> #include <math.h>
#include <stdio.h>
#include "hip/hip_runtime.h" #include "hip/hip_runtime.h"
#define STOKES
#define NBLOCKS 1024 #define NBLOCKS 1024
#define NTHREADS 256 #define NTHREADS 256
#define STOKES
__global__ void dvc_ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np) __global__ void dvc_ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np)
{ {
int n; int n;
@ -186,10 +187,16 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList,
} }
} }
__global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, __global__ void dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
double rhoA, double rhoB, int start, int finish, int Np){ double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
int idx,n;
double fq,phi; int n,idx,nr1,nr2,nr3,nr4,nr5,nr6;
double h0,h1,h2,h3,h4,h5,h6;
double nx,ny,nz,C;
double ux,uy,uz;
double phi;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
double factor = 1.0;
// for (int n=start; n<finish; n++){ // for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1; int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){ for (int s=0; s<S; s++){
@ -197,33 +204,168 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start; n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){ if ( n<finish ){
// q=0
fq = hq[n]; /* load phase indicator field */
phi = fq; idx = Map[n];
phi = Phi[idx];
/* velocity */
ux = Vel[0*Np+n];
uy = Vel[1*Np+n];
uz = Vel[2*Np+n];
/*color gradient */
nx = ColorGrad[0*Np+n];
ny = ColorGrad[1*Np+n];
nz = ColorGrad[2*Np+n];
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C < 1.0e-12) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=1 // q=1
fq = hq[2*Np+n]; nr1 = neighborList[n];
phi += fq; nr2 = neighborList[n+Np];
nr3 = neighborList[n+2*Np];
nr4 = neighborList[n+3*Np];
nr5 = neighborList[n+4*Np];
nr6 = neighborList[n+5*Np];
// f2 = hq[10*Np+n]; //q=0
fq = hq[1*Np+n]; h0 = hq[n];
phi += fq; //q=1
h1 = hq[nr1];
//q=2
h2 = hq[nr2];
//q=3
h3 = hq[nr3];
//q=4
h4 = hq[nr4];
//q=5
h5 = hq[nr5];
//q=6
h6 = hq[nr6];
// q=3 //-------------------------------- BGK collison for phase field ---------------------------------//
fq = hq[4*Np+n]; h0 -= (h0 - 0.3333333333333333*phi)/tauM;
phi += fq; h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//........................................................................
// q = 4 /*Update the distributions */
fq = hq[3*Np+n]; // q = 0
phi += fq; hq[n] = h0;
hq[nr2] = h1;
hq[nr1] = h2;
hq[nr4] = h3;
hq[nr3] = h4;
hq[nr6] = h5;
hq[nr5] = h6;
//........................................................................
// q=5 }
fq = hq[6*Np+n]; }
phi += fq; }
// q = 6
fq = hq[5*Np+n]; __global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField( int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
phi += fq; double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
int idx,n;
double h0,h1,h2,h3,h4,h5,h6;
double nx,ny,nz,C;
double ux,uy,uz;
double phi;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
double factor = 1.0;
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 + start;
if ( n<finish ){
/* load phase indicator field */
idx = Map[n];
phi = Phi[idx];
/* velocity */
ux = Vel[0*Np+n];
uy = Vel[1*Np+n];
uz = Vel[2*Np+n];
/*color gradient */
nx = ColorGrad[0*Np+n];
ny = ColorGrad[1*Np+n];
nz = ColorGrad[2*Np+n];
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C < 1.0e-12) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
h0 = hq[n];
h1 = hq[2*Np+n];
h2 = hq[Np+n];
h3 = hq[4*Np+n];
h4 = hq[3*Np+n];
h5 = hq[6*Np+n];
h6 = hq[5*Np+n];
//-------------------------------- BGK collison for phase field ---------------------------------//
h0 -= (h0 - 0.3333333333333333*phi)/tauM;
h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//........................................................................
/*Update the distributions */
// q = 0
hq[n] = h0;
hq[Np+n] = h1;
hq[2*Np+n] = h2;
hq[3*Np+n] = h3;
hq[4*Np+n] = h4;
hq[5*Np+n] = h5;
hq[6*Np+n] = h6;
//........................................................................
}
}
}
__global__ void dvc_ScaLBL_D3Q7_ComputePhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){
int idx,n;
double h0,h1,h2,h3,h4,h5,h6;
double phi;
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 + start;
if ( n<finish ){
h0 = hq[n];
h1 = hq[1*Np+n];
h2 = hq[2*Np+n];
h3 = hq[3*Np+n];
h4 = hq[4*Np+n];
h5 = hq[5*Np+n];
h6 = hq[6*Np+n];
phi = h0+h1+h2+h3+h4+h5+h6;
// save the number densities // save the number densities
Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
@ -235,8 +377,8 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double
} }
} }
__global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){ int strideY, int strideZ, int start, int finish, int Np){
int n,nn,nn2x,ijk; int n,nn,nn2x,ijk;
@ -256,11 +398,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; //double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid double tau;//position dependent LB relaxation time for fluid
double C,theta; //double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 // double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
// for (int n=start; n<finish; n++){ // for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1; int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){ for (int s=0; s<S; s++){
@ -269,12 +410,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
if ( n<finish ){ if ( n<finish ){
rho0 = Den[n];//load density rho0 = Den[n];//load density
phi = Phi[n];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// Get the 1D index based on regular data layout // Get the 1D index based on regular data layout
ijk = Map[n]; ijk = Map[n];
phi = Phi[ijk];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT // COMPUTE THE COLOR GRADIENT
//........................................................................ //........................................................................
//.................Read Phase Indicator Values............................ //.................Read Phase Indicator Values............................
@ -416,9 +559,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem; chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
//............Compute the Mixed Gradient................................... //............Compute the Mixed Gradient...................................
mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30 mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25; mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25; mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0 // q=0
m0 = dist[n]; m0 = dist[n];
@ -774,62 +917,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)))); 0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------// //----------------------------------------------------------------------------------------------------------------------------------------//
// ----------------------------- compute phase field evolution ----------------------------------------
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//compute surface tension-related parameter
theta = M*4.5*(1-4.0*phi*phi)/W;
//load distributions of phase field
//q=0
h0 = hq[n];
//q=1
h1 = hq[nr1];
//q=2
h2 = hq[nr2];
//q=3
h3 = hq[nr3];
//q=4
h4 = hq[nr4];
//q=5
h5 = hq[nr5];
//q=6
h6 = hq[nr6];
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
// q = 1
hq[nr2] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 2
hq[nr1] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 3
hq[nr4] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 4
hq[nr3] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 5
hq[nr6] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
// q = 6
hq[nr5] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
//........................................................................
//Update velocity on device //Update velocity on device
Vel[0*Np+n] = ux; Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy; Vel[1*Np+n] = uy;
@ -846,12 +933,11 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
} }
} }
__global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){ int strideY, int strideZ, int start, int finish, int Np){
int n,nn,nn2x,ijk; int n,nn,nn2x,ijk;
//int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
double ux,uy,uz;//fluid velocity double ux,uy,uz;//fluid velocity
double p;//pressure double p;//pressure
double chem;//chemical potential double chem;//chemical potential
@ -867,10 +953,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; //double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid double tau;//position dependent LB relaxation time for fluid
double C,theta; //double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 //double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
// for (int n=start; n<finish; n++){ // for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1; int S = Np/NBLOCKS/NTHREADS + 1;
@ -880,12 +966,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
if ( n<finish ){ if ( n<finish ){
rho0 = Den[n];//load density rho0 = Den[n];//load density
phi = Phi[n];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// Get the 1D index based on regular data layout // Get the 1D index based on regular data layout
ijk = Map[n]; ijk = Map[n];
phi = Phi[ijk];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT // COMPUTE THE COLOR GRADIENT
//........................................................................ //........................................................................
//.................Read Phase Indicator Values............................ //.................Read Phase Indicator Values............................
@ -1027,9 +1115,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem; chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
//............Compute the Mixed Gradient................................... //............Compute the Mixed Gradient...................................
mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30 mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25; mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25; mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0 // q=0
m0 = dist[n]; m0 = dist[n];
@ -1091,6 +1179,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0); ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0); uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0); uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
//compute pressure //compute pressure
p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17) p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
+0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz); +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);
@ -1368,62 +1457,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)))); 0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------// //----------------------------------------------------------------------------------------------------------------------------------------//
// ----------------------------- compute phase field evolution ----------------------------------------
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//compute surface tension-related parameter
theta = M*4.5*(1-4.0*phi*phi)/W;
//load distributions of phase field
//q=0
h0 = hq[n];
//q=1
h1 = hq[2*Np+n];
//q=2
h2 = hq[1*Np+n];
//q=3
h3 = hq[4*Np+n];
//q=4
h4 = hq[3*Np+n];
//q=5
h5 = hq[6*Np+n];
//q=6
h6 = hq[5*Np+n];
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
// q = 1
hq[1*Np+n] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 2
hq[2*Np+n] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 3
hq[3*Np+n] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 4
hq[4*Np+n] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 5
hq[5*Np+n] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
// q = 6
hq[6*Np+n] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
//........................................................................
//Update velocity on device //Update velocity on device
Vel[0*Np+n] = ux; Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy; Vel[1*Np+n] = uy;
@ -1436,9 +1469,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
ColorGrad[0*Np+n] = nx; ColorGrad[0*Np+n] = nx;
ColorGrad[1*Np+n] = ny; ColorGrad[1*Np+n] = ny;
ColorGrad[2*Np+n] = nz; ColorGrad[2*Np+n] = nz;
}
} }
}
} }
__global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure, __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,
@ -1970,48 +2003,120 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dis
extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np){ extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np){
dvc_ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init<<<NBLOCKS,NTHREADS >>>( gqbar, mu_phi, ColorGrad, Fx, Fy, Fz, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init: %s \n",hipGetErrorString(err));
}
} }
extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np){ extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np){
dvc_ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init<<<NBLOCKS,NTHREADS >>>( gqbar, Fx, Fy, Fz, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init: %s \n",hipGetErrorString(err));
}
} }
extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad, extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad,
double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){ double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
} dvc_ScaLBL_FreeLeeModel_PhaseField_Init<<<NBLOCKS,NTHREADS >>>(Map, Phi, Den, hq, ColorGrad, rhoA, rhoB, tauM, W, start, finish, Np);
extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, hipError_t err = hipGetLastError();
double rhoA, double rhoB, int start, int finish, int Np){ if (hipSuccess != err){
printf("CUDA error in ScaLBL_FreeLeeModel_PhaseField_Init: %s \n",hipGetErrorString(err));
}
} }
extern "C" void ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np)
double rhoA, double rhoB, int start, int finish, int Np){ {
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField<<<NBLOCKS,NTHREADS >>>(neighborList, Map, hq, Den, Phi, ColorGrad, Vel,
rhoA, rhoB, tauM, W, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in ScaLBL_D3Q7_AAodd_FreeLee_PhaseField: %s \n",hipGetErrorString(err));
}
} }
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, extern "C" void ScaLBL_D3Q7_AAeven_FreeLee_PhaseField( int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField<<<NBLOCKS,NTHREADS >>>( Map, hq, Den, Phi, ColorGrad, Vel, rhoA, rhoB, tauM, W, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in ScaLBL_D3Q7_AAeven_FreeLee_PhaseField: %s \n",hipGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_ComputePhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q7_ComputePhaseField, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q7_ComputePhaseField<<<NBLOCKS,NTHREADS >>>( Map, hq, Den, Phi, rhoA, rhoB, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in ScaLBL_D3Q7_ComputePhaseField: %s \n",hipGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){ int strideY, int strideZ, int start, int finish, int Np){
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel<<<NBLOCKS,NTHREADS >>>(neighborList, Map, dist, Den, Phi, mu_phi, Vel, Pressure, ColorGrad,
rhoA, rhoB, tauA, tauB, kappa, beta, W, Fx, Fy, Fz, strideY, strideZ, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AAodd_FreeLeeModel: %s \n",hipGetErrorString(err));
}
} }
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){ int strideY, int strideZ, int start, int finish, int Np){
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel<<<NBLOCKS,NTHREADS >>>(Map, dist, Den, Phi, mu_phi, Vel, Pressure, ColorGrad,
rhoA, rhoB, tauA, tauB, kappa, beta, W, Fx, Fy, Fz, strideY, strideZ, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AAeven_FreeLeeModel: %s \n",hipGetErrorString(err));
}
} }
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure, extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,
double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){ double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK<<<NBLOCKS,NTHREADS >>>(neighborList, dist, Vel, Pressure,
tau, rho0, Fx, Fy, Fz, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK: %s \n",hipGetErrorString(err));
}
} }
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure,
double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){ double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK<<<NBLOCKS,NTHREADS >>>(dist, Vel, Pressure,
tau, rho0, Fx, Fy, Fz, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK: %s \n",hipGetErrorString(err));
}
} }
extern "C" void ScaLBL_D3Q9_MGTest(int *Map, double *Phi,double *ColorGrad,int strideY, int strideZ, int start, int finish, int Np){
}