LBPM/hip/FreeLee.hip
2021-11-03 10:55:56 -04:00

3529 lines
187 KiB
Plaintext

#include <math.h>
#include <stdio.h>
#include "hip/hip_runtime.h"
#include "hip/hip_runtime_api.h"
#define STOKES
#define NBLOCKS 1024
#define NTHREADS 256
__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;
double p = 1.0;//NOTE: take initial pressure p=1.0
double chem;
double cg_x,cg_y,cg_z;
//for (n=0; n<Np; n++){
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
if ( n<Np ){
chem = mu_phi[n];//chemical potential
cg_x = ColorGrad[0*Np+n];
cg_y = ColorGrad[1*Np+n];
cg_z = ColorGrad[2*Np+n];
gqbar[0*Np+n] = 0.3333333333333333*p;
gqbar[1*Np+n] = 0.055555555555555555*(p - 0.5*(chem*cg_x+Fx)); //double(100*n)+1.f;
gqbar[2*Np+n] = 0.055555555555555555*(p - 0.5*(-chem*cg_x-Fx)); //double(100*n)+2.f;
gqbar[3*Np+n] = 0.055555555555555555*(p - 0.5*(chem*cg_y+Fy)); //double(100*n)+3.f;
gqbar[4*Np+n] = 0.055555555555555555*(p - 0.5*(-chem*cg_y-Fy)); //double(100*n)+4.f;
gqbar[5*Np+n] = 0.055555555555555555*(p - 0.5*(chem*cg_z+Fz)); //double(100*n)+5.f;
gqbar[6*Np+n] = 0.055555555555555555*(p - 0.5*(-chem*cg_z-Fz)); //double(100*n)+6.f;
gqbar[7*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_x+cg_y)+Fx+Fy)); //double(100*n)+7.f;
gqbar[8*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x-cg_y)-Fx-Fy)); //double(100*n)+8.f;
gqbar[9*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_x-cg_y)+Fx-Fy)); //double(100*n)+9.f;
gqbar[10*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x+cg_y)-Fx+Fy)); //double(100*n)+10.f;
gqbar[11*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_x+cg_z)+Fx+Fz)); //double(100*n)+11.f;
gqbar[12*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x-cg_z)-Fx-Fz)); //double(100*n)+12.f;
gqbar[13*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_x-cg_z)+Fx-Fz)); //double(100*n)+13.f;
gqbar[14*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_x+cg_z)-Fx+Fz)); //double(100*n)+14.f;
gqbar[15*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_y+cg_z)+Fy+Fz)); ; //double(100*n)+15.f;
gqbar[16*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_y-cg_z)-Fy-Fz));; //double(100*n)+16.f;
gqbar[17*Np+n] = 0.0277777777777778*(p-0.5*(chem*(cg_y-cg_z)+Fy-Fz)); ; //double(100*n)+17.f;
gqbar[18*Np+n] = 0.0277777777777778*(p-0.5*(chem*(-cg_y+cg_z)-Fy+Fz));; //double(100*n)+18.f;
}
}
}
__global__ void dvc_ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np)
{
int n;
double p = 1.0;//NOTE: take initial pressure p=1.0
// for (n=0; n<Np; n++){
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
if ( n<Np ){
gqbar[0*Np+n] = 0.3333333333333333*p;
gqbar[1*Np+n] = 0.055555555555555555*(p - 0.5*(Fx)); //double(100*n)+1.f;
gqbar[2*Np+n] = 0.055555555555555555*(p - 0.5*(-Fx)); //double(100*n)+2.f;
gqbar[3*Np+n] = 0.055555555555555555*(p - 0.5*(Fy)); //double(100*n)+3.f;
gqbar[4*Np+n] = 0.055555555555555555*(p - 0.5*(-Fy)); //double(100*n)+4.f;
gqbar[5*Np+n] = 0.055555555555555555*(p - 0.5*(Fz)); //double(100*n)+5.f;
gqbar[6*Np+n] = 0.055555555555555555*(p - 0.5*(-Fz)); //double(100*n)+6.f;
gqbar[7*Np+n] = 0.0277777777777778*(p-0.5*(Fx+Fy)); //double(100*n)+7.f;
gqbar[8*Np+n] = 0.0277777777777778*(p-0.5*(-Fx-Fy)); //double(100*n)+8.f;
gqbar[9*Np+n] = 0.0277777777777778*(p-0.5*(Fx-Fy)); //double(100*n)+9.f;
gqbar[10*Np+n] = 0.0277777777777778*(p-0.5*(-Fx+Fy)); //double(100*n)+10.f;
gqbar[11*Np+n] = 0.0277777777777778*(p-0.5*(Fx+Fz)); //double(100*n)+11.f;
gqbar[12*Np+n] = 0.0277777777777778*(p-0.5*(-Fx-Fz)); //double(100*n)+12.f;
gqbar[13*Np+n] = 0.0277777777777778*(p-0.5*(Fx-Fz)); //double(100*n)+13.f;
gqbar[14*Np+n] = 0.0277777777777778*(p-0.5*(-Fx+Fz)); //double(100*n)+14.f;
gqbar[15*Np+n] = 0.0277777777777778*(p-0.5*(Fy+Fz)); ; //double(100*n)+15.f;
gqbar[16*Np+n] = 0.0277777777777778*(p-0.5*(-Fy-Fz));; //double(100*n)+16.f;
gqbar[17*Np+n] = 0.0277777777777778*(p-0.5*(Fy-Fz)); ; //double(100*n)+17.f;
gqbar[18*Np+n] = 0.0277777777777778*(p-0.5*(-Fy+Fz));; //double(100*n)+18.f;
}
}
}
__global__ void dvc_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){
int idx,n;
double phi;
double nx,ny,nz,cg_mag;
double theta;//anti-diffusion term
double cs2_inv = 4.5;//inverse of speed of sound for D3Q7
double M = 1.0/cs2_inv*(tauM-0.5);//diffusivity (or mobility)
for (idx=start; idx<finish; idx++){
n = Map[idx];
phi = Phi[n];
if (phi > 1.f) phi = 1.0;
if (phi < -1.f) phi = -1.0;
Den[idx] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
//compute unit normal of color gradient
nx = ColorGrad[idx+0*Np];
ny = ColorGrad[idx+1*Np];
nz = ColorGrad[idx+2*Np];
cg_mag = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag_temp = cg_mag;
if (cg_mag==0.0) ColorMag_temp=1.0;
nx = nx/ColorMag_temp;
ny = ny/ColorMag_temp;
nz = nz/ColorMag_temp;
//theta = M*cs2_inv*(1-factor*phi*phi)/W;
theta = 4.5*M*2.0*(1-phi*phi)/W;
//theta = 0; // try more diffusive initial condition
hq[0*Np+idx]=0.3333333333333333*(phi);
hq[1*Np+idx]=0.1111111111111111*(phi+theta*nx);
hq[2*Np+idx]=0.1111111111111111*(phi-theta*nx);
hq[3*Np+idx]=0.1111111111111111*(phi+theta*ny);
hq[4*Np+idx]=0.1111111111111111*(phi-theta*ny);
hq[5*Np+idx]=0.1111111111111111*(phi+theta*nz);
hq[6*Np+idx]=0.1111111111111111*(phi-theta*nz);
}
}
__global__ void dvc_ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi,
double rhoA, double rhoB, int start, int finish, int Np){
int idx,n,nread;
double fq,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 ){
// q=0
fq = hq[n];
phi = fq;
// q=1
nread = neighborList[n];
fq = hq[nread];
phi += fq;
// q=2
nread = neighborList[n+Np];
fq = hq[nread];
phi += fq;
// q=3
nread = neighborList[n+2*Np];
fq = hq[nread];
phi += fq;
// q = 4
nread = neighborList[n+3*Np];
fq = hq[nread];
phi += fq;
// q=5
nread = neighborList[n+4*Np];
fq = hq[nread];
phi += fq;
// q = 6
nread = neighborList[n+5*Np];
fq = hq[nread];
phi += fq;
// save the number densities
Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
// save the phase indicator field
idx = Map[n];
Phi[idx] = 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, double tauM, double W, int start, int finish, int Np){
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, theta;
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++){
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];
theta = 4.5*M*2.0*(1-phi*phi)/W;
/* 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
nr1 = neighborList[n];
nr2 = neighborList[n+Np];
nr3 = neighborList[n+2*Np];
nr4 = neighborList[n+3*Np];
nr5 = neighborList[n+4*Np];
nr6 = neighborList[n+5*Np];
//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 ---------------------------------//
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;
h1 -= (h1 - 0.1111111111111111*nx*theta - phi*(0.1111111111111111 + 0.5*ux))/tauM;
h2 -= (h2 + 0.1111111111111111*nx*theta - phi*(0.1111111111111111 - 0.5*ux))/tauM;
h3 -= (h3 - 0.1111111111111111*ny*theta - phi*(0.1111111111111111 + 0.5*uy))/tauM;
h4 -= (h4 + 0.1111111111111111*ny*theta - phi*(0.1111111111111111 - 0.5*uy))/tauM;
h5 -= (h5 - 0.1111111111111111*nz*theta - phi*(0.1111111111111111 + 0.5*uz))/tauM;
h6 -= (h6 + 0.1111111111111111*nz*theta - phi*(0.1111111111111111 - 0.5*uz))/tauM;
//........................................................................
/*Update the distributions */
// q = 0
hq[n] = h0;
hq[nr2] = h1;
hq[nr1] = h2;
hq[nr4] = h3;
hq[nr3] = h4;
hq[nr6] = h5;
hq[nr5] = h6;
//........................................................................
}
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField( int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
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, theta;
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];
theta = 4.5*M*2.0*(1-phi*phi)/W;
/* 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;
h1 -= (h1 - 0.1111111111111111*nx*theta - phi*(0.1111111111111111 + 0.5*ux))/tauM;
h2 -= (h2 + 0.1111111111111111*nx*theta - phi*(0.1111111111111111 - 0.5*ux))/tauM;
h3 -= (h3 - 0.1111111111111111*ny*theta - phi*(0.1111111111111111 + 0.5*uy))/tauM;
h4 -= (h4 + 0.1111111111111111*ny*theta - phi*(0.1111111111111111 - 0.5*uy))/tauM;
h5 -= (h5 - 0.1111111111111111*nz*theta - phi*(0.1111111111111111 + 0.5*uz))/tauM;
h6 -= (h6 + 0.1111111111111111*nz*theta - phi*(0.1111111111111111 - 0.5*uz))/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
Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
// save the phase indicator field
idx = Map[n];
Phi[idx] = phi;
}
}
}
__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 kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){
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 p;//pressure
double chem;//chemical potential
double phi; //phase field
double rho0;//fluid density
// distribution functions
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m0,m3,m5,m7;
double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
double mm3,mm5,mm7;
double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
double nx,ny,nz;//normal color gradient
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 h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid
//double C,theta;
// double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
// for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
rho0 = Den[n];//load density
// Get the 1D index based on regular data layout
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
//........................................................................
//.................Read Phase Indicator Values............................
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Phi[nn]; // get neighbor for phi - 3
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Phi[nn]; // get neighbor for phi - 4
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Phi[nn]; // get neighbor for phi - 5
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Phi[nn]; // get neighbor for phi - 6
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Phi[nn]; // get neighbor for phi - 7
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Phi[nn]; // get neighbor for phi - 8
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Phi[nn]; // get neighbor for phi - 9
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Phi[nn]; // get neighbor for phi - 10
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Phi[nn]; // get neighbor for phi - 11
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Phi[nn]; // get neighbor for phi - 12
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Phi[nn]; // get neighbor for phi - 13
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Phi[nn]; // get neighbor for phi - 14
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Phi[nn]; // get neighbor for phi - 15
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Phi[nn]; // get neighbor for phi - 16
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Phi[nn]; // get neighbor for phi - 17
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Phi[nn]; // get neighbor for phi - 18
// compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
//........................................................................
nn2x = ijk-2; // neighbor index (get convention)
mm1 = Phi[nn2x]; // get neighbor for phi - 1
mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
//........................................................................
nn2x = ijk+2; // neighbor index (get convention)
mm2 = Phi[nn2x]; // get neighbor for phi - 2
mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
//........................................................................
nn2x = ijk-strideY*2; // neighbor index (get convention)
mm3 = Phi[nn2x]; // get neighbor for phi - 3
mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
//........................................................................
nn2x = ijk+strideY*2; // neighbor index (get convention)
mm4 = Phi[nn2x]; // get neighbor for phi - 4
mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
//........................................................................
nn2x = ijk-strideZ*2; // neighbor index (get convention)
mm5 = Phi[nn2x]; // get neighbor for phi - 5
mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
//........................................................................
nn2x = ijk+strideZ*2; // neighbor index (get convention)
mm6 = Phi[nn2x]; // get neighbor for phi - 6
mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
//........................................................................
nn2x = ijk-strideY*2-2; // neighbor index (get convention)
mm7 = Phi[nn2x]; // get neighbor for phi - 7
mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
//........................................................................
nn2x = ijk+strideY*2+2; // neighbor index (get convention)
mm8 = Phi[nn2x]; // get neighbor for phi - 8
mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
//........................................................................
nn2x = ijk+strideY*2-2; // neighbor index (get convention)
mm9 = Phi[nn2x]; // get neighbor for phi - 9
mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
//........................................................................
nn2x = ijk-strideY*2+2; // neighbor index (get convention)
mm10 = Phi[nn2x]; // get neighbor for phi - 10
mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
//........................................................................
nn2x = ijk-strideZ*2-2; // neighbor index (get convention)
mm11 = Phi[nn2x]; // get neighbor for phi - 11
mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
//........................................................................
nn2x = ijk+strideZ*2+2; // neighbor index (get convention)
mm12 = Phi[nn2x]; // get neighbor for phi - 12
mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
//........................................................................
nn2x = ijk+strideZ*2-2; // neighbor index (get convention)
mm13 = Phi[nn2x]; // get neighbor for phi - 13
mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
//........................................................................
nn2x = ijk-strideZ*2+2; // neighbor index (get convention)
mm14 = Phi[nn2x]; // get neighbor for phi - 14
mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
//........................................................................
nn2x = ijk-strideZ*2-strideY*2; // neighbor index (get convention)
mm15 = Phi[nn2x]; // get neighbor for phi - 15
mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
//........................................................................
nn2x = ijk+strideZ*2+strideY*2; // neighbor index (get convention)
mm16 = Phi[nn2x]; // get neighbor for phi - 16
mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
//........................................................................
nn2x = ijk+strideZ*2-strideY*2; // neighbor index (get convention)
mm17 = Phi[nn2x]; // get neighbor for phi - 17
mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
//........................................................................
nn2x = ijk-strideZ*2+strideY*2; // neighbor index (get convention)
mm18 = Phi[nn2x]; // get neighbor for phi - 18
mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);
//............Compute the Color Gradient...................................
nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//............Compute the Chemical Potential...............................
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;
//............Compute the Mixed Gradient...................................
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));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0
m0 = dist[n];
// q=1
nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
m1 = dist[nr1]; // reading the f1 data into register fq
nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
m2 = dist[nr2]; // reading the f2 data into register fq
// q=3
nr3 = neighborList[n+2*Np]; // neighbor 4
m3 = dist[nr3];
// q = 4
nr4 = neighborList[n+3*Np]; // neighbor 3
m4 = dist[nr4];
// q=5
nr5 = neighborList[n+4*Np];
m5 = dist[nr5];
// q = 6
nr6 = neighborList[n+5*Np];
m6 = dist[nr6];
// q=7
nr7 = neighborList[n+6*Np];
m7 = dist[nr7];
// q = 8
nr8 = neighborList[n+7*Np];
m8 = dist[nr8];
// q=9
nr9 = neighborList[n+8*Np];
m9 = dist[nr9];
// q = 10
nr10 = neighborList[n+9*Np];
m10 = dist[nr10];
// q=11
nr11 = neighborList[n+10*Np];
m11 = dist[nr11];
// q=12
nr12 = neighborList[n+11*Np];
m12 = dist[nr12];
// q=13
nr13 = neighborList[n+12*Np];
m13 = dist[nr13];
// q=14
nr14 = neighborList[n+13*Np];
m14 = dist[nr14];
// q=15
nr15 = neighborList[n+14*Np];
m15 = dist[nr15];
// q=16
nr16 = neighborList[n+15*Np];
m16 = dist[nr16];
// q=17
nr17 = neighborList[n+16*Np];
m17 = dist[nr17];
// q=18
nr18 = neighborList[n+17*Np];
m18 = dist[nr18];
//compute fluid velocity
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);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
//compute pressure
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);
//compute equilibrium distributions
feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) -
0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
(-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) -
0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
(2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) -
0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
(-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) +
rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux +
chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) -
0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
(2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) -
0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
(-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) +
rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy +
chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz)));
feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) -
0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) +
rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz +
chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz))));
feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) -
0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) +
rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz +
chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz))));
feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
(2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)));
feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
(2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) +
2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)));
feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) -
0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) +
2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)));
feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) +
2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)));
feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
(2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)));
feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) -
0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
(2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) +
2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz))));
feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) -
0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) +
2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))));
feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) +
2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)));
feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
(2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)));
feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) -
0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
(2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) +
2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))));
feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) -
0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) +
2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))));
feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) +
2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)));
//------------------------------------------------- BCK collison ------------------------------------------------------------//
// q=0
dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) +
(mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) +
0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz))));
// q = 1
dist[nr2] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) +
(mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) +
0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));
// q=2
dist[nr1] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) +
(mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) +
0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));
// q = 3
dist[nr4] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) +
0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 4
dist[nr3] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) +
0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 5
dist[nr6] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) +
0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 6
dist[nr5] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) +
(mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) +
0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));
// q = 7
dist[nr8] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(0.2222222222222222 + (ux + uy)*(ux + uy) -
0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 8
dist[nr7] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uy)*(ux + uy) +
0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 9
dist[nr10] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) +
0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 10
dist[nr9] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) +
0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 11
dist[nr12] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
(0.2222222222222222 + (ux + uz)*(ux + uz) -
0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 12
dist[nr11] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
+ (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
// q = 13
dist[nr14] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uz)*(ux - uz) +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
(mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
// q= 14
dist[nr13] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
(-0.2222222222222222 - (ux - uz)*(ux - uz) +
0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 15
dist[nr16] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
+ (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));
// q = 16
dist[nr15] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
+ (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));
// q = 17
dist[nr18] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
+ (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));
// q = 18
dist[nr17] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
(-0.2222222222222222 - (uy - uz)*(uy - uz) +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) +
(mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------//
//Update velocity on device
Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy;
Vel[2*Np+n] = uz;
//Update pressure on device
Pressure[n] = p;
//Update chemical potential on device
mu_phi[n] = chem;
//Update color gradient on device
ColorGrad[0*Np+n] = nx;
ColorGrad[1*Np+n] = ny;
ColorGrad[2*Np+n] = nz;
}
}
}
__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 kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){
int n,nn,nn2x,ijk;
double ux,uy,uz;//fluid velocity
double p;//pressure
double chem;//chemical potential
double phi; //phase field
double rho0;//fluid density
// distribution functions
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m0,m3,m5,m7;
double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
double mm3,mm5,mm7;
double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
double nx,ny,nz;//normal color gradient
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 h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid
//double C,theta;
//double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
// for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
rho0 = Den[n];//load density
// Get the 1D index based on regular data layout
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
//........................................................................
//.................Read Phase Indicator Values............................
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Phi[nn]; // get neighbor for phi - 3
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Phi[nn]; // get neighbor for phi - 4
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Phi[nn]; // get neighbor for phi - 5
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Phi[nn]; // get neighbor for phi - 6
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Phi[nn]; // get neighbor for phi - 7
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Phi[nn]; // get neighbor for phi - 8
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Phi[nn]; // get neighbor for phi - 9
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Phi[nn]; // get neighbor for phi - 10
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Phi[nn]; // get neighbor for phi - 11
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Phi[nn]; // get neighbor for phi - 12
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Phi[nn]; // get neighbor for phi - 13
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Phi[nn]; // get neighbor for phi - 14
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Phi[nn]; // get neighbor for phi - 15
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Phi[nn]; // get neighbor for phi - 16
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Phi[nn]; // get neighbor for phi - 17
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Phi[nn]; // get neighbor for phi - 18
// compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
//........................................................................
nn2x = ijk-2; // neighbor index (get convention)
mm1 = Phi[nn2x]; // get neighbor for phi - 1
mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
//........................................................................
nn2x = ijk+2; // neighbor index (get convention)
mm2 = Phi[nn2x]; // get neighbor for phi - 2
mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
//........................................................................
nn2x = ijk-strideY*2; // neighbor index (get convention)
mm3 = Phi[nn2x]; // get neighbor for phi - 3
mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
//........................................................................
nn2x = ijk+strideY*2; // neighbor index (get convention)
mm4 = Phi[nn2x]; // get neighbor for phi - 4
mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
//........................................................................
nn2x = ijk-strideZ*2; // neighbor index (get convention)
mm5 = Phi[nn2x]; // get neighbor for phi - 5
mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
//........................................................................
nn2x = ijk+strideZ*2; // neighbor index (get convention)
mm6 = Phi[nn2x]; // get neighbor for phi - 6
mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
//........................................................................
nn2x = ijk-strideY*2-2; // neighbor index (get convention)
mm7 = Phi[nn2x]; // get neighbor for phi - 7
mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
//........................................................................
nn2x = ijk+strideY*2+2; // neighbor index (get convention)
mm8 = Phi[nn2x]; // get neighbor for phi - 8
mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
//........................................................................
nn2x = ijk+strideY*2-2; // neighbor index (get convention)
mm9 = Phi[nn2x]; // get neighbor for phi - 9
mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
//........................................................................
nn2x = ijk-strideY*2+2; // neighbor index (get convention)
mm10 = Phi[nn2x]; // get neighbor for phi - 10
mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
//........................................................................
nn2x = ijk-strideZ*2-2; // neighbor index (get convention)
mm11 = Phi[nn2x]; // get neighbor for phi - 11
mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
//........................................................................
nn2x = ijk+strideZ*2+2; // neighbor index (get convention)
mm12 = Phi[nn2x]; // get neighbor for phi - 12
mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
//........................................................................
nn2x = ijk+strideZ*2-2; // neighbor index (get convention)
mm13 = Phi[nn2x]; // get neighbor for phi - 13
mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
//........................................................................
nn2x = ijk-strideZ*2+2; // neighbor index (get convention)
mm14 = Phi[nn2x]; // get neighbor for phi - 14
mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
//........................................................................
nn2x = ijk-strideZ*2-strideY*2; // neighbor index (get convention)
mm15 = Phi[nn2x]; // get neighbor for phi - 15
mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
//........................................................................
nn2x = ijk+strideZ*2+strideY*2; // neighbor index (get convention)
mm16 = Phi[nn2x]; // get neighbor for phi - 16
mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
//........................................................................
nn2x = ijk+strideZ*2-strideY*2; // neighbor index (get convention)
mm17 = Phi[nn2x]; // get neighbor for phi - 17
mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
//........................................................................
nn2x = ijk-strideZ*2+strideY*2; // neighbor index (get convention)
mm18 = Phi[nn2x]; // get neighbor for phi - 18
mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);
//............Compute the Color Gradient...................................
nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//............Compute the Chemical Potential...............................
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;
//............Compute the Mixed Gradient...................................
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));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0
m0 = dist[n];
// q=1
m1 = dist[2*Np+n];
// q=2
m2 = dist[1*Np+n];
// q=3
m3 = dist[4*Np+n];
// q = 4
m4 = dist[3*Np+n];
// q=5
m5 = dist[6*Np+n];
// q = 6
m6 = dist[5*Np+n];
// q=7
m7 = dist[8*Np+n];
// q = 8
m8 = dist[7*Np+n];
// q=9
m9 = dist[10*Np+n];
// q = 10
m10 = dist[9*Np+n];
// q=11
m11 = dist[12*Np+n];
// q=12
m12 = dist[11*Np+n];
// q=13
m13 = dist[14*Np+n];
// q=14
m14 = dist[13*Np+n];
// q=15
m15 = dist[16*Np+n];
// q=16
m16 = dist[15*Np+n];
// q=17
m17 = dist[18*Np+n];
// q=18
m18 = dist[17*Np+n];
//compute fluid velocity
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);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
//compute pressure
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);
//compute equilibrium distributions
feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) -
0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
(-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) -
0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
(2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) -
0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
(-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) +
rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux +
chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) -
0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
(2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) -
0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
(-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) +
rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy +
chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz)));
feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) -
0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) +
rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz +
chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz))));
feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) -
0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) +
rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz +
chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz))));
feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
(2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)));
feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
(2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) +
2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)));
feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) -
0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) +
2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)));
feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) +
2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)));
feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
(2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)));
feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) -
0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
(2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) +
2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz))));
feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) -
0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) +
2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))));
feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) +
2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)));
feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
(2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)));
feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) -
0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
(2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) +
2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))));
feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) -
0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) +
2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))));
feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) +
2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)));
//------------------------------------------------- BCK collison ------------------------------------------------------------//
// q=0
dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) +
(mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) +
0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz))));
// q = 1
dist[1*Np+n] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) +
(mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) +
0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));
// q=2
dist[2*Np+n] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) +
(mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) +
0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));
// q = 3
dist[3*Np+n] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) +
0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 4
dist[4*Np+n] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) +
0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 5
dist[5*Np+n] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) +
0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 6
dist[6*Np+n] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) +
(mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) +
0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));
// q = 7
dist[7*Np+n] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(0.2222222222222222 + (ux + uy)*(ux + uy) -
0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 8
dist[8*Np+n] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uy)*(ux + uy) +
0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 9
dist[9*Np+n] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) +
0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 10
dist[10*Np+n] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) +
0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 11
dist[11*Np+n] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
(0.2222222222222222 + (ux + uz)*(ux + uz) -
0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 12
dist[12*Np+n] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
+ (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
// q = 13
dist[13*Np+n] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uz)*(ux - uz) +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
(mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
// q= 14
dist[14*Np+n] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
(-0.2222222222222222 - (ux - uz)*(ux - uz) +
0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 15
dist[15*Np+n] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
+ (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));
// q = 16
dist[16*Np+n] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
+ (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));
// q = 17
dist[17*Np+n] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
+ (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));
// q = 18
dist[18*Np+n] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
(-0.2222222222222222 - (uy - uz)*(uy - uz) +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) +
(mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------//
//Update velocity on device
Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy;
Vel[2*Np+n] = uz;
//Update pressure on device
Pressure[n] = p;
//Update chemical potential on device
mu_phi[n] = chem;
//Update color gradient on device
ColorGrad[0*Np+n] = nx;
ColorGrad[1*Np+n] = ny;
ColorGrad[2*Np+n] = nz;
}
}
}
__global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined(int *neighborList, 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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){
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 p;//pressure
double chem;//chemical potential
double phi; //phase field
double rho0;//fluid density
// distribution functions
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m0,m3,m5,m7;
double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
double mm3,mm5,mm7;
double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
double nx,ny,nz;//normal color gradient
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 h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid
double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
double phi_temp;
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
rho0 = Den[n];//load density
// Get the 1D index based on regular data layout
ijk = Map[n];
phi = Phi[ijk];// load phase field
phi_temp = phi;
if (phi>1.f) phi_temp=1.0;
if (phi<-1.f) phi_temp=-1.0;
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT
//........................................................................
//.................Read Phase Indicator Values............................
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Phi[nn]; // get neighbor for phi - 3
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Phi[nn]; // get neighbor for phi - 4
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Phi[nn]; // get neighbor for phi - 5
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Phi[nn]; // get neighbor for phi - 6
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Phi[nn]; // get neighbor for phi - 7
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Phi[nn]; // get neighbor for phi - 8
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Phi[nn]; // get neighbor for phi - 9
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Phi[nn]; // get neighbor for phi - 10
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Phi[nn]; // get neighbor for phi - 11
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Phi[nn]; // get neighbor for phi - 12
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Phi[nn]; // get neighbor for phi - 13
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Phi[nn]; // get neighbor for phi - 14
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Phi[nn]; // get neighbor for phi - 15
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Phi[nn]; // get neighbor for phi - 16
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Phi[nn]; // get neighbor for phi - 17
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Phi[nn]; // get neighbor for phi - 18
// compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
//........................................................................
nn2x = ijk-2; // neighbor index (get convention)
mm1 = Phi[nn2x]; // get neighbor for phi - 1
mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
//........................................................................
nn2x = ijk+2; // neighbor index (get convention)
mm2 = Phi[nn2x]; // get neighbor for phi - 2
mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
//........................................................................
nn2x = ijk-strideY*2; // neighbor index (get convention)
mm3 = Phi[nn2x]; // get neighbor for phi - 3
mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
//........................................................................
nn2x = ijk+strideY*2; // neighbor index (get convention)
mm4 = Phi[nn2x]; // get neighbor for phi - 4
mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
//........................................................................
nn2x = ijk-strideZ*2; // neighbor index (get convention)
mm5 = Phi[nn2x]; // get neighbor for phi - 5
mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
//........................................................................
nn2x = ijk+strideZ*2; // neighbor index (get convention)
mm6 = Phi[nn2x]; // get neighbor for phi - 6
mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
//........................................................................
nn2x = ijk-strideY*2-2; // neighbor index (get convention)
mm7 = Phi[nn2x]; // get neighbor for phi - 7
mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
//........................................................................
nn2x = ijk+strideY*2+2; // neighbor index (get convention)
mm8 = Phi[nn2x]; // get neighbor for phi - 8
mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
//........................................................................
nn2x = ijk+strideY*2-2; // neighbor index (get convention)
mm9 = Phi[nn2x]; // get neighbor for phi - 9
mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
//........................................................................
nn2x = ijk-strideY*2+2; // neighbor index (get convention)
mm10 = Phi[nn2x]; // get neighbor for phi - 10
mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
//........................................................................
nn2x = ijk-strideZ*2-2; // neighbor index (get convention)
mm11 = Phi[nn2x]; // get neighbor for phi - 11
mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
//........................................................................
nn2x = ijk+strideZ*2+2; // neighbor index (get convention)
mm12 = Phi[nn2x]; // get neighbor for phi - 12
mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
//........................................................................
nn2x = ijk+strideZ*2-2; // neighbor index (get convention)
mm13 = Phi[nn2x]; // get neighbor for phi - 13
mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
//........................................................................
nn2x = ijk-strideZ*2+2; // neighbor index (get convention)
mm14 = Phi[nn2x]; // get neighbor for phi - 14
mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
//........................................................................
nn2x = ijk-strideZ*2-strideY*2; // neighbor index (get convention)
mm15 = Phi[nn2x]; // get neighbor for phi - 15
mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
//........................................................................
nn2x = ijk+strideZ*2+strideY*2; // neighbor index (get convention)
mm16 = Phi[nn2x]; // get neighbor for phi - 16
mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
//........................................................................
nn2x = ijk+strideZ*2-strideY*2; // neighbor index (get convention)
mm17 = Phi[nn2x]; // get neighbor for phi - 17
mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
//........................................................................
nn2x = ijk-strideZ*2+strideY*2; // neighbor index (get convention)
mm18 = Phi[nn2x]; // get neighbor for phi - 18
mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);
//............Compute the Color Gradient...................................
nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//............Compute the Chemical Potential...............................
//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 = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi_temp+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi_temp));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi_temp*(phi_temp+1.0)*(phi_temp-1.0)-kappa*chem;
//............Compute the Mixed Gradient...................................
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));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
//de-noise color gradient and mixed gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
if (C<1.0e-12) nx=ny=nz=0.0;
double mg_mag = sqrt(mgx*mgx+mgy*mgy+mgz*mgz);
if (mg_mag<1.0e-12) mgx=mgy=mgz=0.0;
//maybe you can also de-noise chemical potential ? within the bulk phase chem should be ZERO
if (fabs(chem)<1.0e-12) chem=0.0;
// q=0
m0 = dist[n];
// q=1
nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
m1 = dist[nr1]; // reading the f1 data into register fq
nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
m2 = dist[nr2]; // reading the f2 data into register fq
// q=3
nr3 = neighborList[n+2*Np]; // neighbor 4
m3 = dist[nr3];
// q = 4
nr4 = neighborList[n+3*Np]; // neighbor 3
m4 = dist[nr4];
// q=5
nr5 = neighborList[n+4*Np];
m5 = dist[nr5];
// q = 6
nr6 = neighborList[n+5*Np];
m6 = dist[nr6];
// q=7
nr7 = neighborList[n+6*Np];
m7 = dist[nr7];
// q = 8
nr8 = neighborList[n+7*Np];
m8 = dist[nr8];
// q=9
nr9 = neighborList[n+8*Np];
m9 = dist[nr9];
// q = 10
nr10 = neighborList[n+9*Np];
m10 = dist[nr10];
// q=11
nr11 = neighborList[n+10*Np];
m11 = dist[nr11];
// q=12
nr12 = neighborList[n+11*Np];
m12 = dist[nr12];
// q=13
nr13 = neighborList[n+12*Np];
m13 = dist[nr13];
// q=14
nr14 = neighborList[n+13*Np];
m14 = dist[nr14];
// q=15
nr15 = neighborList[n+14*Np];
m15 = dist[nr15];
// q=16
nr16 = neighborList[n+15*Np];
m16 = dist[nr16];
// q=17
nr17 = neighborList[n+16*Np];
m17 = dist[nr17];
// q=18
nr18 = neighborList[n+17*Np];
m18 = dist[nr18];
//compute fluid velocity
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);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
//compute pressure
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);
//compute equilibrium distributions
feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) -
0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
(-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) -
0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
(2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) -
0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
(-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) +
rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux +
chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) -
0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
(2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) -
0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
(-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) +
rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy +
chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz)));
feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) -
0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) +
rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz +
chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz))));
feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) -
0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) +
rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz +
chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz))));
feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
(2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)));
feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
(2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) +
2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)));
feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) -
0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) +
2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)));
feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) +
2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)));
feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
(2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)));
feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) -
0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
(2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) +
2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz))));
feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) -
0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) +
2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))));
feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) +
2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)));
feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
(2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)));
feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) -
0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
(2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) +
2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))));
feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) -
0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) +
2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))));
feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) +
2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)));
//------------------------------------------------- BCK collison ------------------------------------------------------------//
// q=0
dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) +
(mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) +
0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz))));
// q = 1
dist[nr2] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) +
(mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) +
0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));
// q=2
dist[nr1] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) +
(mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) +
0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));
// q = 3
dist[nr4] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) +
0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 4
dist[nr3] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) +
0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 5
dist[nr6] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) +
0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 6
dist[nr5] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) +
(mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) +
0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));
// q = 7
dist[nr8] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(0.2222222222222222 + (ux + uy)*(ux + uy) -
0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 8
dist[nr7] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uy)*(ux + uy) +
0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 9
dist[nr10] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) +
0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 10
dist[nr9] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) +
0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 11
dist[nr12] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
(0.2222222222222222 + (ux + uz)*(ux + uz) -
0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 12
dist[nr11] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
+ (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
// q = 13
dist[nr14] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uz)*(ux - uz) +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
(mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
// q= 14
dist[nr13] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
(-0.2222222222222222 - (ux - uz)*(ux - uz) +
0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 15
dist[nr16] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
+ (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));
// q = 16
dist[nr15] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
+ (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));
// q = 17
dist[nr18] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
+ (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));
// q = 18
dist[nr17] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
(-0.2222222222222222 - (uy - uz)*(uy - uz) +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) +
(mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(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 = 4.5*M*2.0*(1-phi*phi)/W;
theta = 4.5*M*2.0*(1-phi_temp*phi_temp)/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 - 0.1111111111111111*nx*theta - phi*(0.1111111111111111 + 0.5*ux))/tauM;
// q = 2
hq[nr1] = h2 - (h2 + 0.1111111111111111*nx*theta - phi*(0.1111111111111111 - 0.5*ux))/tauM;
// q = 3
hq[nr4] = h3 - (h3 - 0.1111111111111111*ny*theta - phi*(0.1111111111111111 + 0.5*uy))/tauM;
// q = 4
hq[nr3] = h4 - (h4 + 0.1111111111111111*ny*theta - phi*(0.1111111111111111 - 0.5*uy))/tauM;
// q = 5
hq[nr6] = h5 - (h5 - 0.1111111111111111*nz*theta - phi*(0.1111111111111111 + 0.5*uz))/tauM;
// q = 6
hq[nr5] = h6 - (h6 + 0.1111111111111111*nz*theta - phi*(0.1111111111111111 - 0.5*uz))/tauM;
//........................................................................
//Update velocity on device
Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy;
Vel[2*Np+n] = uz;
//Update pressure on device
Pressure[n] = p;
//Update chemical potential on device
mu_phi[n] = chem;
//Update color gradient on device
ColorGrad[0*Np+n] = nx;
ColorGrad[1*Np+n] = ny;
ColorGrad[2*Np+n] = nz;
}
}
__global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){
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 p;//pressure
double chem;//chemical potential
double phi; //phase field
double rho0;//fluid density
// distribution functions
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m0,m3,m5,m7;
double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
double mm3,mm5,mm7;
double feq0,feq1,feq2,feq3,feq4,feq5,feq6,feq7,feq8,feq9,feq10,feq11,feq12,feq13,feq14,feq15,feq16,feq17,feq18;
double nx,ny,nz;//normal color gradient
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 h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid
double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
double phi_temp;
// for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
rho0 = Den[n];//load density
// Get the 1D index based on regular data layout
ijk = Map[n];
phi = Phi[ijk];// load phase field
phi_temp = phi;
if (phi>1.f) phi_temp=1.0;
if (phi<-1.f) phi_temp=-1.0;
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT
//........................................................................
//.................Read Phase Indicator Values............................
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Phi[nn]; // get neighbor for phi - 3
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Phi[nn]; // get neighbor for phi - 4
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Phi[nn]; // get neighbor for phi - 5
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Phi[nn]; // get neighbor for phi - 6
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Phi[nn]; // get neighbor for phi - 7
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Phi[nn]; // get neighbor for phi - 8
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Phi[nn]; // get neighbor for phi - 9
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Phi[nn]; // get neighbor for phi - 10
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Phi[nn]; // get neighbor for phi - 11
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Phi[nn]; // get neighbor for phi - 12
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Phi[nn]; // get neighbor for phi - 13
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Phi[nn]; // get neighbor for phi - 14
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Phi[nn]; // get neighbor for phi - 15
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Phi[nn]; // get neighbor for phi - 16
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Phi[nn]; // get neighbor for phi - 17
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Phi[nn]; // get neighbor for phi - 18
// compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
//........................................................................
nn2x = ijk-2; // neighbor index (get convention)
mm1 = Phi[nn2x]; // get neighbor for phi - 1
mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
//........................................................................
nn2x = ijk+2; // neighbor index (get convention)
mm2 = Phi[nn2x]; // get neighbor for phi - 2
mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
//........................................................................
nn2x = ijk-strideY*2; // neighbor index (get convention)
mm3 = Phi[nn2x]; // get neighbor for phi - 3
mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
//........................................................................
nn2x = ijk+strideY*2; // neighbor index (get convention)
mm4 = Phi[nn2x]; // get neighbor for phi - 4
mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
//........................................................................
nn2x = ijk-strideZ*2; // neighbor index (get convention)
mm5 = Phi[nn2x]; // get neighbor for phi - 5
mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
//........................................................................
nn2x = ijk+strideZ*2; // neighbor index (get convention)
mm6 = Phi[nn2x]; // get neighbor for phi - 6
mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
//........................................................................
nn2x = ijk-strideY*2-2; // neighbor index (get convention)
mm7 = Phi[nn2x]; // get neighbor for phi - 7
mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
//........................................................................
nn2x = ijk+strideY*2+2; // neighbor index (get convention)
mm8 = Phi[nn2x]; // get neighbor for phi - 8
mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
//........................................................................
nn2x = ijk+strideY*2-2; // neighbor index (get convention)
mm9 = Phi[nn2x]; // get neighbor for phi - 9
mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
//........................................................................
nn2x = ijk-strideY*2+2; // neighbor index (get convention)
mm10 = Phi[nn2x]; // get neighbor for phi - 10
mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
//........................................................................
nn2x = ijk-strideZ*2-2; // neighbor index (get convention)
mm11 = Phi[nn2x]; // get neighbor for phi - 11
mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
//........................................................................
nn2x = ijk+strideZ*2+2; // neighbor index (get convention)
mm12 = Phi[nn2x]; // get neighbor for phi - 12
mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
//........................................................................
nn2x = ijk+strideZ*2-2; // neighbor index (get convention)
mm13 = Phi[nn2x]; // get neighbor for phi - 13
mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
//........................................................................
nn2x = ijk-strideZ*2+2; // neighbor index (get convention)
mm14 = Phi[nn2x]; // get neighbor for phi - 14
mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
//........................................................................
nn2x = ijk-strideZ*2-strideY*2; // neighbor index (get convention)
mm15 = Phi[nn2x]; // get neighbor for phi - 15
mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
//........................................................................
nn2x = ijk+strideZ*2+strideY*2; // neighbor index (get convention)
mm16 = Phi[nn2x]; // get neighbor for phi - 16
mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
//........................................................................
nn2x = ijk+strideZ*2-strideY*2; // neighbor index (get convention)
mm17 = Phi[nn2x]; // get neighbor for phi - 17
mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
//........................................................................
nn2x = ijk-strideZ*2+strideY*2; // neighbor index (get convention)
mm18 = Phi[nn2x]; // get neighbor for phi - 18
mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);
//............Compute the Color Gradient...................................
nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//............Compute the Chemical Potential...............................
//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 = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi_temp+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi_temp));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi_temp*(phi_temp+1.0)*(phi_temp-1.0)-kappa*chem;
//............Compute the Mixed Gradient...................................
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));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
//de-noise color gradient and mixed gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
if (C<1.0e-12) nx=ny=nz=0.0;
double mg_mag = sqrt(mgx*mgx+mgy*mgy+mgz*mgz);
if (mg_mag<1.0e-12) mgx=mgy=mgz=0.0;
//maybe you can also de-noise chemical potential ? within the bulk phase chem should be ZERO
if (fabs(chem)<1.0e-12) chem=0.0;
// q=0
m0 = dist[n];
// q=1
m1 = dist[2*Np+n];
// q=2
m2 = dist[1*Np+n];
// q=3
m3 = dist[4*Np+n];
// q = 4
m4 = dist[3*Np+n];
// q=5
m5 = dist[6*Np+n];
// q = 6
m6 = dist[5*Np+n];
// q=7
m7 = dist[8*Np+n];
// q = 8
m8 = dist[7*Np+n];
// q=9
m9 = dist[10*Np+n];
// q = 10
m10 = dist[9*Np+n];
// q=11
m11 = dist[12*Np+n];
// q=12
m12 = dist[11*Np+n];
// q=13
m13 = dist[14*Np+n];
// q=14
m14 = dist[13*Np+n];
// q=15
m15 = dist[16*Np+n];
// q=16
m16 = dist[15*Np+n];
// q=17
m17 = dist[18*Np+n];
// q=18
m18 = dist[17*Np+n];
//compute fluid velocity
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);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
//compute pressure
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);
//compute equilibrium distributions
feq0 = 0.3333333333333333*p - 0.25*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) -
0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz) - 0.5*(-(nx*ux) - ny*uy - nz*uz)*
(-0.08333333333333333*(rhoA - rhoB)*(ux*ux + uy*uy + uz*uz) + chem*(0.3333333333333333 - 0.5*(ux*ux + uy*uy + uz*uz)));
feq1 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) -
0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx - nx*ux - ny*uy - nz*uz)*
(2*chem*ux*ux - 0.3333333333333333*((-rhoA + rhoB)*ux*ux + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz)));
feq2 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-ux*ux + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) -
0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)) - 0.0625*(nx + nx*ux + ny*uy + nz*uz)*
(-2.*chem*ux*ux + 0.1111111111111111*(-4.*chem + rhoB*(-2.*ux - 1.*ux*ux - 1.*uy*uy - 1.*uz*uz) +
rhoA*(2.*ux + ux*ux + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*ux*ux +
chem*(4.*ux + 2.*ux*ux + 2.*uy*uy + 2.*uz*uz)));
feq3 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) -
0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny - nx*ux - ny*uy - nz*uz)*
(2*chem*uy*uy - 0.3333333333333333*((-rhoA + rhoB)*uy*uy + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz)));
feq4 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uy*uy + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) -
0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.0625*(ny + nx*ux + ny*uy + nz*uz)*
(-2.*chem*uy*uy + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 2.*uy - 1.*uy*uy - 1.*uz*uz) +
rhoA*(ux*ux + 2.*uy + uy*uy + uz*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uy*uy +
chem*(2.*ux*ux + 4.*uy + 2.*uy*uy + 2.*uz*uz)));
feq5 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) -
0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)) - 0.0625*(nx*ux + ny*uy + nz*(-1. + uz))*
(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (2. - 1.*uz)*uz) +
rhoA*(ux*ux + uy*uy + (-2. + uz)*uz)) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz +
chem*(2.*ux*ux + 2.*uy*uy + uz*(-4. + 2.*uz))));
feq6 = 0.05555555555555555*p - 0.08333333333333333*rho0*(-uz*uz + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) -
0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))) - 0.0625*(nz + nx*ux + ny*uy + nz*uz)*
(-2.*chem*uz*uz + 0.1111111111111111*(-4.*chem + rhoB*(-1.*ux*ux - 1.*uy*uy + (-2. - 1.*uz)*uz) +
rhoA*(ux*ux + uy*uy + uz*(2. + uz))) + 0.3333333333333333*((-1.*rhoA + rhoB)*uz*uz +
chem*(2.*ux*ux + 2.*uy*uy + uz*(4. + 2.*uz))));
feq7 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uy)*(ux + uy) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx + ny - nx*ux - ny*uy - nz*uz)*
(2*chem*(ux + uy)*(ux + uy) + 0.3333333333333333*((rhoA - rhoB)*(ux + uy)*(ux + uy) - 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)));
feq8 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(-(nx*(1 + ux)) - ny*(1 + uy) - nz*uz)*
(2*chem*(ux + uy)*(ux + uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uy)*(ux + uy)) +
2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)));
feq9 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) -
0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)) - 0.03125*(nx - nx*ux - ny*(1 + uy) - nz*uz)*
(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) +
2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)));
feq10 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) -
0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uy - 1.*uy*uy +
0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)) - 0.03125*(ny - nx*(1 + ux) - ny*uy - nz*uz)*
(2*chem*(ux - uy)*(ux - uy) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uy)*(ux - uy)) +
2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)));
feq11 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uz)*(ux + uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nx + nz - nx*ux - ny*uy - nz*uz)*
(2*chem*(ux + uz)*(ux + uz) + 0.3333333333333333*((rhoA - rhoB)*(ux + uz)*(ux + uz) - 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)));
feq12 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) -
0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux - 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*(1 + ux)) - ny*uy - nz*(1 + uz))*
(2*chem*(ux + uz)*(ux + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux + uz)*(ux + uz)) +
2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz))));
feq13 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) -
0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))) - 0.03125*(nx - nx*ux - ny*uy - nz*(1 + uz))*
(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) +
2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))));
feq14 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*ux*ux + 2.*ux*uz - 1.*uz*uz +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*(1 + ux) - ny*uy - nz*uz)*
(2*chem*(ux - uz)*(ux - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(ux - uz)*(ux - uz)) +
2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)));
feq15 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(ny + nz - nx*ux - ny*uy - nz*uz)*
(2*chem*(uy + uz)*(uy + uz) + 0.3333333333333333*((rhoA - rhoB)*(uy + uz)*(uy + uz) - 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)));
feq16 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) -
0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy - 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(-(nx*ux) - ny*(1 + uy) - nz*(1 + uz))*
(2*chem*(uy + uz)*(uy + uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy + uz)*(uy + uz)) +
2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))));
feq17 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) -
0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))) - 0.03125*(ny - nx*ux - ny*uy - nz*(1 + uz))*
(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) +
2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))));
feq18 = 0.027777777777777776*p - 0.041666666666666664*rho0*
(-(uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) -
0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*(-0.2222222222222222 - 1.*uy*uy + 2.*uy*uz - 1.*uz*uz +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)) - 0.03125*(nz - nx*ux - ny*(1 + uy) - nz*uz)*
(2*chem*(uy - uz)*(uy - uz) - 0.3333333333333333*(-((rhoA - rhoB)*(uy - uz)*(uy - uz)) +
2*chem*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + 0.1111111111111111*
(4*chem - (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)));
//------------------------------------------------- BCK collison ------------------------------------------------------------//
// q=0
dist[n] = m0 - (m0-feq0)/tau + 0.25*(2*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) +
(mgx*ux + mgy*uy + mgz*uz)*(2*chem*(ux*ux + uy*uy + uz*uz) +
0.3333333333333333*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*uz))));
// q = 1
dist[1*Np+n] = m1 - (m1-feq1)/tau + 0.125*(2*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) +
(mgx*(-1 + ux) + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) +
0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*uz))));
// q=2
dist[2*Np+n] = m2 - (m2-feq2)/tau + 0.125*(2*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) +
(mgx + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(ux*ux) +
0.3333333333333333*((-rhoA + rhoB)*(ux*ux) + 2*chem*(2*ux + ux*ux + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*uz))));
// q = 3
dist[3*Np+n] = m3 - (m3-feq3)/tau + 0.125*(2*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*ux + mgy*(-1 + uy) + mgz*uz)*(-2*chem*(uy*uy) +
0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 4
dist[4*Np+n] = m4 - (m4-feq4)/tau + 0.125*(2*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgy + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uy*uy) +
0.3333333333333333*((-rhoA + rhoB)*(uy*uy) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 5
dist[5*Np+n] = m5 - (m5-feq5)/tau + 0.125*(2*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*ux + mgy*uy + mgz*(-1 + uz))*(-2*chem*(uz*uz) +
0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 6
dist[6*Np+n] = m6 - (m6-feq6)/tau + 0.125*(2*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) +
(mgz + mgx*ux + mgy*uy + mgz*uz)*(-2*chem*(uz*uz) +
0.3333333333333333*((-rhoA + rhoB)*(uz*uz) + 2*chem*(ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + uy*uy + uz*(2 + uz)))));
// q = 7
dist[7*Np+n] = m7 - (m7-feq7)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(0.2222222222222222 + (ux + uy)*(ux + uy) -
0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*(-1 + ux) + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 8
dist[8*Np+n] = m8 - (m8-feq8)/tau + 0.0625*(2*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uy)*(ux + uy) +
0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgx + mgy + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((ux + uy)*(ux + uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uy)*(ux + uy))) + 2*chem*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 9
dist[9*Np+n] = m9 - (m9-feq9)/tau + 0.0625*(2*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) +
0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
(mgy + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))));
// q = 10
dist[10*Np+n] = m10 - (m10-feq10)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) +
0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
(mgx*(1 + ux) + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((ux - uy)*(ux - uy)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uy)*(ux - uy))) + 2*chem*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))));
// q = 11
dist[11*Np+n] = m11 - (m11-feq11)/tau + 0.0625*(-2*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
(0.2222222222222222 + (ux + uz)*(ux + uz) -
0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*(-1 + ux) + mgy*uy + mgz*(-1 + uz))*
(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 12
dist[12*Np+n] = m12 - (m12-feq12)/tau + 0.0625*(2*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))
+ (mgx + mgz + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((ux + uz)*(ux + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux + uz)*(ux + uz))) + 2*chem*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
// q = 13
dist[13*Np+n] = m13 - (m13-feq13)/tau + 0.0625*(2*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uz)*(ux - uz) +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
(mgz + mgx*(-1 + ux) + mgy*uy + mgz*uz)*
(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))));
// q= 14
dist[14*Np+n] = m14 - (m14-feq14)/tau + 0.0625*(2*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
(-0.2222222222222222 - (ux - uz)*(ux - uz) +
0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
(mgx*(1 + ux) + mgy*uy + mgz*(-1 + uz))*
(-2*chem*((ux - uz)*(ux - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((ux - uz)*(ux - uz))) + 2*chem*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))));
// q = 15
dist[15*Np+n] = m15 - (m15-feq15)/tau + 0.0625*(-2*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))
+ (mgx*ux + mgy*(-1 + uy) + mgz*(-1 + uz))*
(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))));
// q = 16
dist[16*Np+n] = m16 - (m16-feq16)/tau + 0.0625*(2*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))
+ (mgy + mgz + mgx*ux + mgy*uy + mgz*uz)*
(-2*chem*((uy + uz)*(uy + uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy + uz)*(uy + uz))) + 2*chem*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))));
// q = 17
dist[17*Np+n] = m17 - (m17-feq17)/tau + 0.0625*(2*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))
+ (mgz + mgx*ux + mgy*(-1 + uy) + mgz*uz)*
(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) +
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))));
// q = 18
dist[18*Np+n] = m18 - (m18-feq18)/tau + 0.0625*(2*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
(-0.2222222222222222 - (uy - uz)*(uy - uz) +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) +
(mgx*ux + mgy*(1 + uy) + mgz*(-1 + uz))*
(-2*chem*((uy - uz)*(uy - uz)) + 0.3333333333333333*
(-((rhoA - rhoB)*((uy - uz)*(uy - uz))) + 2*chem*(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 = 4.5*M*2.0*(1-phi*phi)/W;
theta = 4.5*M*2.0*(1-phi_temp*phi_temp)/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 - 0.1111111111111111*nx*theta - phi*(0.1111111111111111 + 0.5*ux))/tauM;
// q = 2
hq[2*Np+n] = h2 - (h2 + 0.1111111111111111*nx*theta - phi*(0.1111111111111111 - 0.5*ux))/tauM;
// q = 3
hq[3*Np+n] = h3 - (h3 - 0.1111111111111111*ny*theta - phi*(0.1111111111111111 + 0.5*uy))/tauM;
// q = 4
hq[4*Np+n] = h4 - (h4 + 0.1111111111111111*ny*theta - phi*(0.1111111111111111 - 0.5*uy))/tauM;
// q = 5
hq[5*Np+n] = h5 - (h5 - 0.1111111111111111*nz*theta - phi*(0.1111111111111111 + 0.5*uz))/tauM;
// q = 6
hq[6*Np+n] = h6 - (h6 + 0.1111111111111111*nz*theta - phi*(0.1111111111111111 - 0.5*uz))/tauM;
//........................................................................
//Update velocity on device
Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy;
Vel[2*Np+n] = uz;
//Update pressure on device
Pressure[n] = p;
//Update chemical potential on device
mu_phi[n] = chem;
//Update color gradient on device
ColorGrad[0*Np+n] = nx;
ColorGrad[1*Np+n] = ny;
ColorGrad[2*Np+n] = nz;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField_alt(int *neighborList, int *Map, double *hq, double *Den, double *Phi,
double rhoA, double rhoB, int start, int finish, int Np){
int n,idx,nread;
double fq,phi;
// for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
// q=0
fq = hq[n];
phi = fq;
// q=1
nread = neighborList[n];
fq = hq[nread];
phi += fq;
// q=2
nread = neighborList[n+Np];
fq = hq[nread];
phi += fq;
// q=3
nread = neighborList[n+2*Np];
fq = hq[nread];
phi += fq;
// q = 4
nread = neighborList[n+3*Np];
fq = hq[nread];
phi += fq;
// q=5
nread = neighborList[n+4*Np];
fq = hq[nread];
phi += fq;
// q = 6
nread = neighborList[n+5*Np];
fq = hq[nread];
phi += fq;
// save the number densities
Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
// save the phase indicator field
idx = Map[n];
Phi[idx] = phi;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi,
double rhoA, double rhoB, int start, int finish, int Np){
int n,idx;
double fq,phi;
// for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
}
// q=0
fq = hq[n];
phi = fq;
// q=1
fq = hq[2*Np+n];
phi += fq;
// f2 = hq[10*Np+n];
fq = hq[1*Np+n];
phi += fq;
// q=3
fq = hq[4*Np+n];
phi += fq;
// q = 4
fq = hq[3*Np+n];
phi += fq;
// q=5
fq = hq[6*Np+n];
phi += fq;
// q = 6
fq = hq[5*Np+n];
phi += fq;
// save the number densities
Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
// save the phase indicator field
idx = Map[n];
Phi[idx] = phi;
}
}
__global__ void dvc_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){
int n;
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 p;//pressure
// distribution functions
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m0,m3,m5,m7;
// for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
// q=0
m0 = dist[n];
// q=1
nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
m1 = dist[nr1]; // reading the f1 data into register fq
nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
m2 = dist[nr2]; // reading the f2 data into register fq
// q=3
nr3 = neighborList[n+2*Np]; // neighbor 4
m3 = dist[nr3];
// q = 4
nr4 = neighborList[n+3*Np]; // neighbor 3
m4 = dist[nr4];
// q=5
nr5 = neighborList[n+4*Np];
m5 = dist[nr5];
// q = 6
nr6 = neighborList[n+5*Np];
m6 = dist[nr6];
// q=7
nr7 = neighborList[n+6*Np];
m7 = dist[nr7];
// q = 8
nr8 = neighborList[n+7*Np];
m8 = dist[nr8];
// q=9
nr9 = neighborList[n+8*Np];
m9 = dist[nr9];
// q = 10
nr10 = neighborList[n+9*Np];
m10 = dist[nr10];
// q=11
nr11 = neighborList[n+10*Np];
m11 = dist[nr11];
// q=12
nr12 = neighborList[n+11*Np];
m12 = dist[nr12];
// q=13
nr13 = neighborList[n+12*Np];
m13 = dist[nr13];
// q=14
nr14 = neighborList[n+13*Np];
m14 = dist[nr14];
// q=15
nr15 = neighborList[n+14*Np];
m15 = dist[nr15];
// q=16
nr16 = neighborList[n+15*Np];
m16 = dist[nr16];
// q=17
nr17 = neighborList[n+16*Np];
m17 = dist[nr17];
// q=18
nr18 = neighborList[n+17*Np];
m18 = dist[nr18];
//compute fluid velocity
ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(Fx)/3.0);
uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(Fy)/3.0);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(Fz)/3.0);
//compute pressure
p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17);
//------------------------------------------------- BCK collison ------------------------------------------------------------//
// q=0
dist[n] = m0 + 0.5*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) -
(m0 - 0.3333333333333333*p + 0.25*(Fx*ux + Fy*uy + Fz*uz)*
(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz))/
tau;
// q = 1
dist[nr2] = m1 + 0.25*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) -
(m1 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(ux*ux) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) + 0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)))/tau;
// q=2
dist[nr1] = m2 + 0.25*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) -
(m2 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(ux*ux) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) +
0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(ux*ux) +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)))/tau;
// q = 3
dist[nr4] = m3 + 0.25*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) -
(m3 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(uy*uy) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*
(-0.2222222222222222 - 1.*(uy*uy) + 0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
// q = 4
dist[nr3] = m4 + 0.25*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) -
(m4 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(uy*uy) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) +
0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uy*uy) +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
// q = 5
dist[nr6] = m5 + 0.25*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) -
(m5 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) +
0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
// q = 6
dist[nr5] = m6 + 0.25*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) -
(m6 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) +
0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uz*uz) +
0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))))/tau;
// q = 7
dist[nr8] = m7 - 0.125*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(0.2222222222222222 + (ux + uy)*(ux + uy) - 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
- (m7 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux + uy)*(ux + uy)) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) +
0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
// q = 8
dist[nr7] = m8 + 0.125*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))\
- (m8 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux + uy)*(ux + uy)) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) +
0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
// q = 9
dist[nr10] = m9 + 0.125*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))
- (m9 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux - uy)*(ux - uy)) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) +
0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
// q = 10
dist[nr9] = m10 + 0.125*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
- (m10 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux - uy)*(ux - uy)) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) +
0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
// q = 11
dist[nr12] = m11 - 0.125*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
(0.2222222222222222 + (ux + uz)*(ux + uz) - 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
- (m11 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux + uz)*(ux + uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
// q = 12
dist[nr11] = m12 + 0.125*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
- (m12 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux + uz)*(ux + uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;
// q = 13
dist[nr14] = m13 + 0.125*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
- (m13 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux - uz)*(ux - uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;
// q= 14
dist[nr13] = m14 + 0.125*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
(-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
- (m14 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux - uz)*(ux - uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
// q = 15
dist[nr16] = m15 - 0.125*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))\
- (m15 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
// q = 16
dist[nr15] = m16 + 0.125*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))\
- (m16 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) +
0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))))/tau;
// q = 17
dist[nr18] = m17 + 0.125*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))\
- (m17 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) +
0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*
(-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))))/tau;
// q = 18
dist[nr17] = m18 + 0.125*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))\
- (m18 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) +
0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
//----------------------------------------------------------------------------------------------------------------------------------------//
//Update velocity on device
Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy;
Vel[2*Np+n] = uz;
//Update pressure on device
Pressure[n] = p;
}
}
}
__global__ void dvc_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){
int n;
double ux,uy,uz;//fluid velocity
double p;//pressure
// distribution functions
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m0,m3,m5,m7;
// for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
// q=0
m0 = dist[n];
// q=1
m1 = dist[2*Np+n];
// q=2
m2 = dist[1*Np+n];
// q=3
m3 = dist[4*Np+n];
// q = 4
m4 = dist[3*Np+n];
// q=5
m5 = dist[6*Np+n];
// q = 6
m6 = dist[5*Np+n];
// q=7
m7 = dist[8*Np+n];
// q = 8
m8 = dist[7*Np+n];
// q=9
m9 = dist[10*Np+n];
// q = 10
m10 = dist[9*Np+n];
// q=11
m11 = dist[12*Np+n];
// q=12
m12 = dist[11*Np+n];
// q=13
m13 = dist[14*Np+n];
// q=14
m14 = dist[13*Np+n];
// q=15
m15 = dist[16*Np+n];
// q=16
m16 = dist[15*Np+n];
// q=17
m17 = dist[18*Np+n];
// q=18
m18 = dist[17*Np+n];
//compute fluid velocity
ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(Fx)/3.0);
uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(Fy)/3.0);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(Fz)/3.0);
//compute pressure
p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17);
//------------------------------------------------- BCK collison ------------------------------------------------------------//
// q=0
dist[n] = m0 + 0.5*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) -
(m0 - 0.3333333333333333*p + 0.25*(Fx*ux + Fy*uy + Fz*uz)*
(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz))/
tau;
// q = 1
dist[1*Np+n] = m1 + 0.25*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) -
(m1 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(ux*ux) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) +
0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) + 0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)))/tau;
// q=2
dist[2*Np+n] = m2 + 0.25*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux +
0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) -
(m2 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(ux*ux) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) +
0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(ux*ux) +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)))/tau;
// q = 3
dist[3*Np+n] = m3 + 0.25*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) -
(m3 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(uy*uy) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) +
0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)*
(-0.2222222222222222 - 1.*(uy*uy) + 0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
// q = 4
dist[4*Np+n] = m4 + 0.25*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy +
0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) -
(m4 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(uy*uy) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) +
0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uy*uy) +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
// q = 5
dist[5*Np+n] = m5 + 0.25*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) -
(m5 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) +
0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
// q = 6
dist[6*Np+n] = m6 + 0.25*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz +
0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) -
(m6 - 0.05555555555555555*p + 0.08333333333333333*rho0*
(-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) +
0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uz*uz) +
0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))))/tau;
// q = 7
dist[7*Np+n] = m7 - 0.125*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(0.2222222222222222 + (ux + uy)*(ux + uy) - 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
- (m7 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux + uy)*(ux + uy)) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) +
0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
// q = 8
dist[8*Np+n] = m8 + 0.125*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))\
- (m8 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux + uy)*(ux + uy)) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) +
0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
// q = 9
dist[9*Np+n] = m9 + 0.125*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz))
- (m9 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux - uy)*(ux - uy)) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) +
0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) +
0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau;
// q = 10
dist[10*Np+n] = m10 + 0.125*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\
- (m10 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux - uy)*(ux - uy)) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) +
0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) +
0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau;
// q = 11
dist[11*Np+n] = m11 - 0.125*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))*
(0.2222222222222222 + (ux + uz)*(ux + uz) - 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
- (m11 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux + uz)*(ux + uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
// q = 12
dist[12*Np+n] = m12 + 0.125*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
- (m12 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux + uz)*(ux + uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;
// q = 13
dist[13*Np+n] = m13 + 0.125*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))\
- (m13 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux - uz)*(ux - uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) +
0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) +
0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau;
// q= 14
dist[14*Np+n] = m14 + 0.125*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))*
(-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\
- (m14 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((ux - uz)*(ux - uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) +
0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) +
0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau;
// q = 15
dist[15*Np+n] = m15 - 0.125*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))*
(0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))\
- (m15 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) +
0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
// q = 16
dist[16*Np+n] = m16 + 0.125*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))\
- (m16 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) +
0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)*
(-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))))/tau;
// q = 17
dist[17*Np+n] = m17 + 0.125*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)*
(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))\
- (m17 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) +
0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)*
(-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) +
0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))))/tau;
// q = 18
dist[18*Np+n] = m18 + 0.125*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))*
(-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))\
- (m18 - 0.027777777777777776*p + 0.041666666666666664*rho0*
(-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) +
0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))*
(-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) +
0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)))/tau;
//----------------------------------------------------------------------------------------------------------------------------------------//
//Update velocity on device
Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy;
Vel[2*Np+n] = uz;
//Update pressure on device
Pressure[n] = p;
}
}
}
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("hip 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){
dvc_ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init<<<NBLOCKS,NTHREADS >>>( gqbar, Fx, Fy, Fz, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("hip 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,
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);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("hip 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,
double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np)
{
//hipFuncSetCacheConfig((void*)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("hip error in ScaLBL_D3Q7_AAodd_FreeLee_PhaseField: %s \n",hipGetErrorString(err));
}
}
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 tauM, double W, int start, int finish, int Np)
{
//hipFuncSetCacheConfig((void*)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("hip 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((void*)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("hip 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){
//hipFuncSetCacheConfig((void*)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("hip 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,
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){
//hipFuncSetCacheConfig((void*)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("hip error in ScaLBL_D3Q19_AAeven_FreeLeeModel: %s \n",hipGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){
//hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined<<<NBLOCKS,NTHREADS >>>(Map, dist, Den, hq, Phi, mu_phi, Vel, Pressure, ColorGrad,
rhoA, rhoB, tauA, tauB, tauM, kappa, beta, W, Fx, Fy, Fz, strideY, strideZ, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("hip error in ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined: %s \n",hipGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined(int *neighborList, 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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){
//hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined<<<NBLOCKS,NTHREADS >>>(neighborList, Map, dist, hq, Den, Phi, mu_phi, Vel, Pressure, ColorGrad,
rhoA, rhoB, tauA, tauB, tauM, kappa, beta, W, Fx, Fy, Fz, strideY, strideZ, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("hip error in ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined: %s \n",hipGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi,
double rhoA, double rhoB, int start, int finish, int Np){
//hipFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField<<<NBLOCKS,NTHREADS >>>( neighborList, Map, hq, Den, Phi, rhoA, rhoB, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("hip error in ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField: %s \n",hipGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){
//hipFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField<<<NBLOCKS,NTHREADS >>>( Map, hq, Den, Phi, rhoA, rhoB, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("hip error in ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField: %s \n",hipGetErrorString(err));
}
}
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){
//hipFuncSetCacheConfig((void*)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("hip 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,
double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np)
{
//hipFuncSetCacheConfig((void*)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("hip 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){
}