From 5109c08bc378c9a965b445f062b2beeb38a0ac65 Mon Sep 17 00:00:00 2001 From: Mark Berrill Date: Wed, 7 Oct 2020 12:31:49 -0400 Subject: [PATCH] Initial hip addition --- CMakeLists.txt | 26 +- {gpu => cuda}/BGK.cu | 0 {gpu => cuda}/Color.cu | 0 {gpu => cuda}/CudaExtras.cu | 0 {gpu => cuda}/D3Q19.cu | 0 {gpu => cuda}/D3Q7.cu | 0 {gpu => cuda}/Extras.cu | 0 {gpu => cuda}/MRT.cu | 0 {gpu => cuda}/dfh.cu | 0 {gpu => cuda}/exe/CMakeLists.txt | 0 {gpu => cuda}/exe/lb1_MRT-swap.cu | 0 {gpu => cuda}/exe/lb1_MRT.cu | 0 {gpu => cuda}/exe/lb1_MRT_mpi.cpp | 0 {gpu => cuda}/exe/lb1_MRT_mpi.cu | 0 {gpu => cuda}/exe/lb2_Color.cpp | 0 {gpu => cuda}/exe/lb2_Color.cu | 0 {gpu => cuda}/exe/lb2_Color_mpi.cpp | 0 {gpu => cuda}/exe/lb2_Color_pBC_wia_mpi.cpp | 0 hip/BGK.hip | 311 ++ hip/CMakeLists.txt | 9 + hip/Color.hip | 4131 +++++++++++++++++++ hip/CudaExtras.hip | 34 + hip/D3Q19.hip | 2645 ++++++++++++ hip/D3Q7.hip | 246 ++ hip/Extras.hip | 62 + hip/MRT.hip | 310 ++ hip/dfh.hip | 1508 +++++++ 27 files changed, 9278 insertions(+), 4 deletions(-) rename {gpu => cuda}/BGK.cu (100%) rename {gpu => cuda}/Color.cu (100%) rename {gpu => cuda}/CudaExtras.cu (100%) rename {gpu => cuda}/D3Q19.cu (100%) rename {gpu => cuda}/D3Q7.cu (100%) rename {gpu => cuda}/Extras.cu (100%) rename {gpu => cuda}/MRT.cu (100%) rename {gpu => cuda}/dfh.cu (100%) rename {gpu => cuda}/exe/CMakeLists.txt (100%) rename {gpu => cuda}/exe/lb1_MRT-swap.cu (100%) rename {gpu => cuda}/exe/lb1_MRT.cu (100%) rename {gpu => cuda}/exe/lb1_MRT_mpi.cpp (100%) rename {gpu => cuda}/exe/lb1_MRT_mpi.cu (100%) rename {gpu => cuda}/exe/lb2_Color.cpp (100%) rename {gpu => cuda}/exe/lb2_Color.cu (100%) rename {gpu => cuda}/exe/lb2_Color_mpi.cpp (100%) rename {gpu => cuda}/exe/lb2_Color_pBC_wia_mpi.cpp (100%) create mode 100644 hip/BGK.hip create mode 100644 hip/CMakeLists.txt create mode 100644 hip/Color.hip create mode 100644 hip/CudaExtras.hip create mode 100644 hip/D3Q19.hip create mode 100644 hip/D3Q7.hip create mode 100644 hip/Extras.hip create mode 100644 hip/MRT.hip create mode 100644 hip/dfh.hip diff --git a/CMakeLists.txt b/CMakeLists.txt index 1e7eeaea..33528b62 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -104,7 +104,7 @@ IF ( USE_DOXYGEN ) ADD_DEPENDENCIES( doc latex_docs doxygen ) ELSE() SET( USE_DOXYGEN 0 ) - ENDIF() + ENDIF()lbpm-wia ENDIF() @@ -123,11 +123,26 @@ IF ( USE_CUDA ) ADD_DEFINITIONS( -DUSE_CUDA ) ENABLE_LANGUAGE( CUDA ) ELSEIF ( USE_HIP ) - FIND_PACKAGE( HIP ) - MESSAGE( FATAL_ERROR "STOP" ) + IF ( NOT DEFINED HIP_PATH ) + IF ( NOT DEFINED ENV{HIP_PATH} ) + SET( HIP_PATH "/opt/rocm/hip" CACHE PATH "Path to which HIP has been installed" ) + ELSE() + SET( HIP_PATH $ENV{HIP_PATH} CACHE PATH "Path to which HIP has been installed" ) + ENDIF() + ENDIF() + SET( CMAKE_MODULE_PATH "${HIP_PATH}/cmake" ${CMAKE_MODULE_PATH} ) + FIND_PACKAGE( HIP REQUIRED ) + FIND_PACKAGE( CUDA QUIET ) + MESSAGE( "HIP Found") + MESSAGE( " HIP version: ${HIP_VERSION_STRING}") + MESSAGE( " HIP platform: ${HIP_PLATFORM}") + MESSAGE( " HIP Include Path: ${HIP_INCLUDE_DIRS}") + MESSAGE( " HIP Libraries: ${HIP_LIBRARIES}") + ADD_DEFINITIONS( -DUSE_HIP ) ENDIF() + # Configure external packages IF ( NOT ONLY_BUILD_DOCS ) CONFIGURE_MPI() # MPI must be before other libraries @@ -161,7 +176,10 @@ IF ( NOT ONLY_BUILD_DOCS ) ADD_PACKAGE_SUBDIRECTORY( StackTrace ) ADD_PACKAGE_SUBDIRECTORY( models ) IF ( USE_CUDA ) - ADD_PACKAGE_SUBDIRECTORY( gpu ) + ADD_PACKAGE_SUBDIRECTORY( cuda ) + ELSEIF ( USE_HIP ) + ADD_SUBDIRECTORY( gpu ) + SET( LBPM_LIBRARIES lbpm-hip lbpm-wia ) ELSE() ADD_PACKAGE_SUBDIRECTORY( cpu ) ENDIF() diff --git a/gpu/BGK.cu b/cuda/BGK.cu similarity index 100% rename from gpu/BGK.cu rename to cuda/BGK.cu diff --git a/gpu/Color.cu b/cuda/Color.cu similarity index 100% rename from gpu/Color.cu rename to cuda/Color.cu diff --git a/gpu/CudaExtras.cu b/cuda/CudaExtras.cu similarity index 100% rename from gpu/CudaExtras.cu rename to cuda/CudaExtras.cu diff --git a/gpu/D3Q19.cu b/cuda/D3Q19.cu similarity index 100% rename from gpu/D3Q19.cu rename to cuda/D3Q19.cu diff --git a/gpu/D3Q7.cu b/cuda/D3Q7.cu similarity index 100% rename from gpu/D3Q7.cu rename to cuda/D3Q7.cu diff --git a/gpu/Extras.cu b/cuda/Extras.cu similarity index 100% rename from gpu/Extras.cu rename to cuda/Extras.cu diff --git a/gpu/MRT.cu b/cuda/MRT.cu similarity index 100% rename from gpu/MRT.cu rename to cuda/MRT.cu diff --git a/gpu/dfh.cu b/cuda/dfh.cu similarity index 100% rename from gpu/dfh.cu rename to cuda/dfh.cu diff --git a/gpu/exe/CMakeLists.txt b/cuda/exe/CMakeLists.txt similarity index 100% rename from gpu/exe/CMakeLists.txt rename to cuda/exe/CMakeLists.txt diff --git a/gpu/exe/lb1_MRT-swap.cu b/cuda/exe/lb1_MRT-swap.cu similarity index 100% rename from gpu/exe/lb1_MRT-swap.cu rename to cuda/exe/lb1_MRT-swap.cu diff --git a/gpu/exe/lb1_MRT.cu b/cuda/exe/lb1_MRT.cu similarity index 100% rename from gpu/exe/lb1_MRT.cu rename to cuda/exe/lb1_MRT.cu diff --git a/gpu/exe/lb1_MRT_mpi.cpp b/cuda/exe/lb1_MRT_mpi.cpp similarity index 100% rename from gpu/exe/lb1_MRT_mpi.cpp rename to cuda/exe/lb1_MRT_mpi.cpp diff --git a/gpu/exe/lb1_MRT_mpi.cu b/cuda/exe/lb1_MRT_mpi.cu similarity index 100% rename from gpu/exe/lb1_MRT_mpi.cu rename to cuda/exe/lb1_MRT_mpi.cu diff --git a/gpu/exe/lb2_Color.cpp b/cuda/exe/lb2_Color.cpp similarity index 100% rename from gpu/exe/lb2_Color.cpp rename to cuda/exe/lb2_Color.cpp diff --git a/gpu/exe/lb2_Color.cu b/cuda/exe/lb2_Color.cu similarity index 100% rename from gpu/exe/lb2_Color.cu rename to cuda/exe/lb2_Color.cu diff --git a/gpu/exe/lb2_Color_mpi.cpp b/cuda/exe/lb2_Color_mpi.cpp similarity index 100% rename from gpu/exe/lb2_Color_mpi.cpp rename to cuda/exe/lb2_Color_mpi.cpp diff --git a/gpu/exe/lb2_Color_pBC_wia_mpi.cpp b/cuda/exe/lb2_Color_pBC_wia_mpi.cpp similarity index 100% rename from gpu/exe/lb2_Color_pBC_wia_mpi.cpp rename to cuda/exe/lb2_Color_pBC_wia_mpi.cpp diff --git a/hip/BGK.hip b/hip/BGK.hip new file mode 100644 index 00000000..f3e746af --- /dev/null +++ b/hip/BGK.hip @@ -0,0 +1,311 @@ +#include + +#define NBLOCKS 1024 +#define NTHREADS 256 + +__global__ void dvc_ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){ + int n; + // conserved momemnts + double rho,ux,uy,uz,uu; + // non-conserved moments + double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s 10Np => odd part of dist) + f1 = dist[nr1]; // reading the f1 data into register fq + + nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) + f2 = dist[nr2]; // reading the f2 data into register fq + + // q=3 + nr3 = neighborList[n+2*Np]; // neighbor 4 + f3 = dist[nr3]; + + // q = 4 + nr4 = neighborList[n+3*Np]; // neighbor 3 + f4 = dist[nr4]; + + // q=5 + nr5 = neighborList[n+4*Np]; + f5 = dist[nr5]; + + // q = 6 + nr6 = neighborList[n+5*Np]; + f6 = dist[nr6]; + + // q=7 + nr7 = neighborList[n+6*Np]; + f7 = dist[nr7]; + + // q = 8 + nr8 = neighborList[n+7*Np]; + f8 = dist[nr8]; + + // q=9 + nr9 = neighborList[n+8*Np]; + f9 = dist[nr9]; + + // q = 10 + nr10 = neighborList[n+9*Np]; + f10 = dist[nr10]; + + // q=11 + nr11 = neighborList[n+10*Np]; + f11 = dist[nr11]; + + // q=12 + nr12 = neighborList[n+11*Np]; + f12 = dist[nr12]; + + // q=13 + nr13 = neighborList[n+12*Np]; + f13 = dist[nr13]; + + // q=14 + nr14 = neighborList[n+13*Np]; + f14 = dist[nr14]; + + // q=15 + nr15 = neighborList[n+14*Np]; + f15 = dist[nr15]; + + // q=16 + nr16 = neighborList[n+15*Np]; + f16 = dist[nr16]; + + // q=17 + //fq = dist[18*Np+n]; + nr17 = neighborList[n+16*Np]; + f17 = dist[nr17]; + + // q=18 + nr18 = neighborList[n+17*Np]; + f18 = dist[nr18]; + + rho = f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17; + ux = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14; + uy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18; + uz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18; + uu = 1.5*(ux*ux+uy*uy+uz*uz); + + // q=0 + dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*(1.0-uu); + + // q = 1 + dist[nr2] = f1*(1.0-rlx) + rlx*0.05555555555555555*(rho + 3.0*ux + 4.5*ux*ux - uu) + 0.16666666*Fx; + + // q=2 + dist[nr1] = f2*(1.0-rlx) + rlx*0.05555555555555555*(rho - 3.0*ux + 4.5*ux*ux - uu)- 0.16666666*Fx; + + // q = 3 + dist[nr4] = f3*(1.0-rlx) + + rlx*0.05555555555555555*(rho + 3.0*uy + 4.5*uy*uy - uu) + 0.16666666*Fy; + + // q = 4 + dist[nr3] = f4*(1.0-rlx) + + rlx*0.05555555555555555*(rho - 3.0*uy + 4.5*uy*uy - uu)- 0.16666666*Fy; + + // q = 5 + dist[nr6] = f5*(1.0-rlx) + + rlx*0.05555555555555555*(rho + 3.0*uz + 4.5*uz*uz - uu) + 0.16666666*Fz; + + // q = 6 + dist[nr5] = f6*(1.0-rlx) + + rlx*0.05555555555555555*(rho - 3.0*uz + 4.5*uz*uz - uu) - 0.16666666*Fz; + + // q = 7 + dist[nr8] = f7*(1.0-rlx) + + rlx*0.02777777777777778*(rho + 3.0*(ux+uy) + 4.5*(ux+uy)*(ux+uy) - uu) + 0.08333333333*(Fx+Fy); + + // q = 8 + dist[nr7] = f8*(1.0-rlx) + + rlx*0.02777777777777778*(rho - 3.0*(ux+uy) + 4.5*(ux+uy)*(ux+uy) - uu) - 0.08333333333*(Fx+Fy); + + // q = 9 + dist[nr10] = f9*(1.0-rlx) + + rlx*0.02777777777777778*(rho + 3.0*(ux-uy) + 4.5*(ux-uy)*(ux-uy) - uu) + 0.08333333333*(Fx-Fy); + + // q = 10 + dist[nr9] = f10*(1.0-rlx) + + rlx*0.02777777777777778*(rho - 3.0*(ux-uy) + 4.5*(ux-uy)*(ux-uy) - uu) - 0.08333333333*(Fx-Fy); + + // q = 11 + dist[nr12] = f11*(1.0-rlx) + + rlx*0.02777777777777778*(rho + 3.0*(ux+uz) + 4.5*(ux+uz)*(ux+uz) - uu) + 0.08333333333*(Fx+Fz); + + // q = 12 + dist[nr11] = f12*(1.0-rlx) + + rlx*0.02777777777777778*(rho - 3.0*(ux+uz) + 4.5*(ux+uz)*(ux+uz) - uu) - 0.08333333333*(Fx+Fz); + + // q = 13 + dist[nr14] = f13*(1.0-rlx) + + rlx*0.02777777777777778*(rho + 3.0*(ux-uz) + 4.5*(ux-uz)*(ux-uz) - uu) + 0.08333333333*(Fx-Fz); + + // q= 14 + dist[nr13] = f14*(1.0-rlx) + + rlx*0.02777777777777778*(rho - 3.0*(ux-uz) + 4.5*(ux-uz)*(ux-uz) - uu)- 0.08333333333*(Fx-Fz); + + // q = 15 + dist[nr16] = f15*(1.0-rlx) + + rlx*0.02777777777777778*(rho + 3.0*(uy+uz) + 4.5*(uy+uz)*(uy+uz) - uu) + 0.08333333333*(Fy+Fz); + + // q = 16 + dist[nr15] = f16*(1.0-rlx) + + rlx*0.02777777777777778*(rho - 3.0*(uy+uz) + 4.5*(uy+uz)*(uy+uz) - uu) - 0.08333333333*(Fy+Fz); + + // q = 17 + dist[nr18] = f17*(1.0-rlx) + + rlx*0.02777777777777778*(rho + 3.0*(uy-uz) + 4.5*(uy-uz)*(uy-uz) - uu) + 0.08333333333*(Fy-Fz); + + // q = 18 + dist[nr17] = f18*(1.0-rlx) + + rlx*0.02777777777777778*(rho - 3.0*(uy-uz) + 4.5*(uy-uz)*(uy-uz) - uu) - 0.08333333333*(Fy-Fz); + } + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){ + + dvc_ScaLBL_D3Q19_AAeven_BGK<<>>(dist,start,finish,Np,rlx,Fx,Fy,Fz); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_BGK: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){ + dvc_ScaLBL_D3Q19_AAodd_BGK<<>>(neighborList,dist,start,finish,Np,rlx,Fx,Fy,Fz); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_BGK: %s \n",hipGetErrorString(err)); + } +} diff --git a/hip/CMakeLists.txt b/hip/CMakeLists.txt new file mode 100644 index 00000000..38ef7c27 --- /dev/null +++ b/hip/CMakeLists.txt @@ -0,0 +1,9 @@ +SET( HIP_SEPERABLE_COMPILATION ON ) +SET_SOURCE_FILES_PROPERTIES( BGK.hip Color.hip CudaExtras.hip D3Q19.hip D3Q7.hip dfh.hip Extras.hip MRT.hip PROPERTIES HIP_SOURCE_PROPERTY_FORMAT 1 ) +HIP_ADD_LIBRARY( lbpm-hip BGK.hip Color.hip CudaExtras.hip D3Q19.hip D3Q7.hip dfh.hip Extras.hip MRT.hip SHARED HIPCC_OPTIONS ${HIP_HIPCC_OPTIONS} HCC_OPTIONS ${HIP_HCC_OPTIONS} NVCC_OPTIONS ${HIP_NVCC_OPTIONS} ${HIP_NVCC_FLAGS} ) +TARGET_LINK_LIBRARIES( lbpm-hip /opt/rocm-3.3.0/lib/libhip_hcc.so ) +TARGET_LINK_LIBRARIES( lbpm-wia lbpm-hip ) +ADD_DEPENDENCIES( lbpm-hip copy-include ) + + + diff --git a/hip/Color.hip b/hip/Color.hip new file mode 100644 index 00000000..b802ab1f --- /dev/null +++ b/hip/Color.hip @@ -0,0 +1,4131 @@ +#include +#include +#include "hip/hip_runtime.h" + +#define NBLOCKS 1024 +#define NTHREADS 256 + +__global__ void dvc_ScaLBL_Color_Init(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz) +{ + //int i,j,k; + int n,N; + char id; + N = Nx*Ny*Nz; + + int S = N/NBLOCKS/NTHREADS + 1; + for (int s=0; s 0){ + + // Retrieve the color gradient + nx = ColorGrad[n]; + ny = ColorGrad[N+n]; + nz = ColorGrad[2*N+n]; + //...........Normalize the Color Gradient................................. + C = sqrt(nx*nx+ny*ny+nz*nz); + nx = nx/C; + ny = ny/C; + nz = nz/C; + //......No color gradient at z-boundary if pressure BC are set............. + // if (pBC && k==0) nx = ny = nz = 0.f; + // if (pBC && k==Nz-1) nx = ny = nz = 0.f; + //........................................................................ + // READ THE DISTRIBUTIONS + // (read from opposite array due to previous swap operation) + //........................................................................ + f2 = distodd[n]; + f4 = distodd[N+n]; + f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; + //........................................................................ + f0 = disteven[n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + f17 = disteven[9*N+n]; + //........................................................................ + // PERFORM RELAXATION PROCESS + //........................................................................ + //....................compute the moments............................................... + rho = f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17; + m1 = -30*f0-11*(f2+f1+f4+f3+f6+f5)+8*(f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18 +f17); + m2 = 12*f0-4*(f2+f1 +f4+f3+f6 +f5)+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17; + jx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14; + m4 = 4*(-f1+f2)+f7-f8+f9-f10+f11-f12+f13-f14; + jy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18; + m6 = -4*(f3-f4)+f7-f8-f9+f10+f15-f16+f17-f18; + jz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18; + m8 = -4*(f5-f6)+f11-f12-f13+f14+f15-f16-f17+f18; + m9 = 2*(f1+f2)-f3-f4-f5-f6+f7+f8+f9+f10+f11+f12+f13+f14-2*(f15+f16+f17+f18); + m10 = -4*(f1+f2)+2*(f4+f3+f6+f5)+f8+f7+f10+f9+f12+f11+f14+f13-2*(f16+f15+f18+f17); + m11 = f4+f3-f6-f5+f8+f7+f10+f9-f12-f11-f14-f13; + m12 = -2*(f4+f3-f6-f5)+f8+f7+f10+f9-f12-f11-f14-f13; + m13 = f8+f7-f10-f9; + m14 = f16+f15-f18-f17; + m15 = f12+f11-f14-f13; + m16 = f7-f8+f9-f10-f11+f12-f13+f14; + m17 = -f7+f8+f9-f10+f15-f16+f17-f18; + m18 = f11-f12-f13+f14-f15+f16+f17-f18; + //..........Toelke, Fruediger et. al. 2006............... + if (C == 0.0) nx = ny = nz = 1.0; + m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) -alpha*C - m1); + m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho)- m2); + m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); + m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); + m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); + m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); + m10 = m10 + rlx_setA*(-0.5*((2*jx*jx-jy*jy-jz*jz)/rho) - m10); + m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); + m12 = m12 + rlx_setA*( -0.5*((jy*jy-jz*jz)/rho) - m12); + m13 = m13 + rlx_setA*( (jx*jy/rho) + 0.5*alpha*C*nx*ny - m13); + m14 = m14 + rlx_setA*( (jy*jz/rho) + 0.5*alpha*C*ny*nz - m14); + m15 = m15 + rlx_setA*( (jx*jz/rho) + 0.5*alpha*C*nx*nz - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); + //.................inverse transformation...................................................... + f0 = 0.05263157894736842*rho-0.012531328320802*m1+0.04761904761904762*m2; + f1 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(jx-m4)+0.0555555555555555555555555*(m9-m10); + f2 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(m4-jx)+0.0555555555555555555555555*(m9-m10); + f3 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(jy-m6)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m11-m12); + f4 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(m6-jy)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m11-m12); + f5 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(jz-m8)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m12-m11); + f6 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(m8-jz)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m12-m11); + f7 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2+0.1*(jx+jy)+0.025*(m4+m6) + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12+0.25*m13+0.125*(m16-m17); + f8 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2-0.1*(jx+jy)-0.025*(m4+m6) + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12+0.25*m13+0.125*(m17-m16); + f9 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2+0.1*(jx-jy)+0.025*(m4-m6) + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12-0.25*m13+0.125*(m16+m17); + f10 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2+0.1*(jy-jx)+0.025*(m6-m4) + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12-0.25*m13-0.125*(m16+m17); + f11 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jx+jz)+0.025*(m4+m8) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12+0.25*m15+0.125*(m18-m16); + f12 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2-0.1*(jx+jz)-0.025*(m4+m8) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12+0.25*m15+0.125*(m16-m18); + f13 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jx-jz)+0.025*(m4-m8) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12-0.25*m15-0.125*(m16+m18); + f14 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jz-jx)+0.025*(m8-m4) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12-0.25*m15+0.125*(m16+m18); + f15 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jy+jz)+0.025*(m6+m8) + -0.0555555555555555555555555*m9-0.02777777777777778*m10+0.25*m14+0.125*(m17-m18); + f16 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2-0.1*(jy+jz)-0.025*(m6+m8) + -0.0555555555555555555555555*m9-0.02777777777777778*m10+0.25*m14+0.125*(m18-m17); + f17 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jy-jz)+0.025*(m6-m8) + -0.0555555555555555555555555*m9-0.02777777777777778*m10-0.25*m14+0.125*(m17+m18); + f18 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jz-jy)+0.025*(m8-m6) + -0.0555555555555555555555555*m9-0.02777777777777778*m10-0.25*m14-0.125*(m17+m18); + //....................................................................................................... + // incorporate external force + f1 += 0.16666666*Fx; + f2 -= 0.16666666*Fx; + f3 += 0.16666666*Fy; + f4 -= 0.16666666*Fy; + f5 += 0.16666666*Fz; + f6 -= 0.16666666*Fz; + f7 += 0.08333333333*(Fx+Fy); + f8 -= 0.08333333333*(Fx+Fy); + f9 += 0.08333333333*(Fx-Fy); + f10 -= 0.08333333333*(Fx-Fy); + f11 += 0.08333333333*(Fx+Fz); + f12 -= 0.08333333333*(Fx+Fz); + f13 += 0.08333333333*(Fx-Fz); + f14 -= 0.08333333333*(Fx-Fz); + f15 += 0.08333333333*(Fy+Fz); + f16 -= 0.08333333333*(Fy+Fz); + f17 += 0.08333333333*(Fy-Fz); + f18 -= 0.08333333333*(Fy-Fz); + //*********** WRITE UPDATED VALUES TO MEMORY ****************** + // Write the updated distributions + //....EVEN..................................... + disteven[n] = f0; + disteven[N+n] = f2; + disteven[2*N+n] = f4; + disteven[3*N+n] = f6; + disteven[4*N+n] = f8; + disteven[5*N+n] = f10; + disteven[6*N+n] = f12; + disteven[7*N+n] = f14; + disteven[8*N+n] = f16; + disteven[9*N+n] = f18; + //....ODD...................................... + distodd[n] = f1; + distodd[N+n] = f3; + distodd[2*N+n] = f5; + distodd[3*N+n] = f7; + distodd[4*N+n] = f9; + distodd[5*N+n] = f11; + distodd[6*N+n] = f13; + distodd[7*N+n] = f15; + distodd[8*N+n] = f17; + + //...Store the Velocity.......................... + Velocity[n] = jx; + Velocity[N+n] = jy; + Velocity[2*N+n] = jz; + /* Velocity[3*n] = jx; + Velocity[3*n+1] = jy; + Velocity[3*n+2] = jz; + */ //...Store the Color Gradient.................... + // ColorGrad[3*n] = nx*C; + // ColorGrad[3*n+1] = ny*C; + // ColorGrad[3*n+2] = nz*C; + //............................................... + //*************************************************************** + } // check if n is in the solid + } // loop over n +} + +__global__ void +__launch_bounds__(512,2) +dvc_ScaLBL_D3Q19_ColorCollide( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad, + double *Velocity, int Nx, int Ny, int Nz, double rlx_setA, double rlx_setB, + double alpha, double beta, double Fx, double Fy, double Fz) +{ + + int i,j,k,n,nn,N; + // distributions + double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9; + double f10,f11,f12,f13,f14,f15,f16,f17,f18; + + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + // additional variables needed for computations + double rho,jx,jy,jz,C,nx,ny,nz; + char id; + + N = Nx*Ny*Nz; + + int S = N/NBLOCKS/NTHREADS + 1; + for (int s=0; s 0){ + + //.......Back out the 3-D indices for node n.............. + k = n/(Nx*Ny); + j = (n-Nx*Ny*k)/Nx; + i = n-Nx*Ny*k-Nx*j; + //........................................................................ + //........Get 1-D index for this thread.................... + // n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x; + //........................................................................ + // COMPUTE THE COLOR GRADIENT + //........................................................................ + //.................Read Phase Indicator Values............................ + //........................................................................ + nn = n-1; // neighbor index (get convention) + if (i-1<0) nn += Nx; // periodic BC along the x-boundary + f1 = phi[nn]; // get neighbor for phi - 1 + //........................................................................ + nn = n+1; // neighbor index (get convention) + if (!(i+10)) delta=0; + a1 = na*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nb*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = na*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nb*(0.1111111111111111*(1-4.5*ux))+delta; + + A_odd[n] = a1; + A_even[N+n] = a2; + B_odd[n] = b1; + B_even[N+n] = b2; + //............................................... + // q = 2 + // Cq = {0,1,0} + delta = beta*na*nb*nab*0.1111111111111111*ny; + if (!(na*nb*nab>0)) delta=0; + a1 = na*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nb*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = na*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nb*(0.1111111111111111*(1-4.5*uy))+delta; + + A_odd[N+n] = a1; + A_even[2*N+n] = a2; + B_odd[N+n] = b1; + B_even[2*N+n] = b2; + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*na*nb*nab*0.1111111111111111*nz; + if (!(na*nb*nab>0)) delta=0; + a1 = na*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nb*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = na*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nb*(0.1111111111111111*(1-4.5*uz))+delta; + + A_odd[2*N+n] = a1; + A_even[3*N+n] = a2; + B_odd[2*N+n] = b1; + B_even[3*N+n] = b2; + + } + } +} + +//************************************************************************* +__global__ void dvc_DensityStreamD3Q7(char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity, + double beta, int Nx, int Ny, int Nz, bool pBC) +{ + char id; + + int idx; + int in,jn,kn,n,nn,N; + int q,Cqx,Cqy,Cqz; + // int sendLoc; + + double na,nb; // density values + double ux,uy,uz; // flow velocity + double nx,ny,nz,C; // color gradient components + double a1,a2,b1,b2; + double sp,delta; + double feq[6]; // equilibrium distributions + // Set of Discrete velocities for the D3Q19 Model + int D3Q7[3][3]={{1,0,0},{0,1,0},{0,0,1}}; + N = Nx*Ny*Nz; + + int S = N/NBLOCKS/NTHREADS + 1; + for (int s=0; s 0 && na+nb > 0.0){ + //.......Back out the 3-D indices for node n.............. + int k = n/(Nx*Ny); + int j = (n-Nx*Ny*k)/Nx; + int i = n-Nx*Ny*k-Nx*j; + //.....Load the Color gradient......... + nx = ColorGrad[n]; + ny = ColorGrad[N+n]; + nz = ColorGrad[2*N+n]; + C = sqrt(nx*nx+ny*ny+nz*nz); + if (C == 0.0) C=1.0; + nx = nx/C; + ny = ny/C; + nz = nz/C; + //....Load the flow velocity........... + ux = Velocity[n]; + uy = Velocity[N+n]; + uz = Velocity[2*N+n]; + //....Instantiate the density distributions + // Generate Equilibrium Distributions and stream + // Stationary value - distribution 0 + // Den[2*n] += 0.3333333333333333*na; + // Den[2*n+1] += 0.3333333333333333*nb; + Den[2*n] += 0.3333333333333333*na; + Den[2*n+1] += 0.3333333333333333*nb; + // Non-Stationary equilibrium distributions + feq[0] = 0.1111111111111111*(1+3*ux); + feq[1] = 0.1111111111111111*(1-3*ux); + feq[2] = 0.1111111111111111*(1+3*uy); + feq[3] = 0.1111111111111111*(1-3*uy); + feq[4] = 0.1111111111111111*(1+3*uz); + feq[5] = 0.1111111111111111*(1-3*uz); + // Construction and streaming for the components + for (idx=0; idx<3; idx++){ + // Distribution index + q = 2*idx; + // Associated discrete velocity + Cqx = D3Q7[idx][0]; + Cqy = D3Q7[idx][1]; + Cqz = D3Q7[idx][2]; + // Generate the Equilibrium Distribution + a1 = na*feq[q]; + b1 = nb*feq[q]; + a2 = na*feq[q+1]; + b2 = nb*feq[q+1]; + // Recolor the distributions + if (C > 0.0){ + sp = nx*double(Cqx)+ny*double(Cqy)+nz*double(Cqz); + //if (idx > 2) sp = 0.7071067811865475*sp; + //delta = sp*min( min(a1,a2), min(b1,b2) ); + delta = na*nb/(na+nb)*0.1111111111111111*sp; + //if (a1>0 && b1>0){ + a1 += beta*delta; + a2 -= beta*delta; + b1 -= beta*delta; + b2 += beta*delta; + } + + // .......Get the neighbor node.............. + //nn = n + Stride[idx]; + in = i+Cqx; + jn = j+Cqy; + kn = k+Cqz; + + // Adjust for periodic BC, if necessary + // if (in<0) in+= Nx; + // if (jn<0) jn+= Ny; + // if (kn<0) kn+= Nz; + // if (!(in 0){ + // Get the density value (Streaming already performed) + Na = Den[n]; + Nb = Den[N+n]; + Phi[n] = (Na-Nb)/(Na+Nb); + } + } + } + //................................................................... +} + +__global__ void dvc_ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, int Nz, int Slice) +{ + int n = Slice*Nx*Ny + blockIdx.x*blockDim.x + threadIdx.x; + if (n < (Slice+1)*Nx*Ny){ + Phi[n] = value; + } +} + + + +__global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ + int ijk,nn,n; + double fq; + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + double m3,m5,m7; + double nA,nB; // number density + double a1,b1,a2,b2,nAB,delta; + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + double phi,tau,rho0,rlx_setA,rlx_setB; + + const double mrt_V1=0.05263157894736842; + const double mrt_V2=0.012531328320802; + const double mrt_V3=0.04761904761904762; + const double mrt_V4=0.004594820384294068; + const double mrt_V5=0.01587301587301587; + const double mrt_V6=0.0555555555555555555555555; + const double mrt_V7=0.02777777777777778; + const double mrt_V8=0.08333333333333333; + const double mrt_V9=0.003341687552213868; + const double mrt_V10=0.003968253968253968; + const double mrt_V11=0.01388888888888889; + const double mrt_V12=0.04166666666666666; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + Aq[1*Np+n] = a1; + Bq[1*Np+n] = b1; + Aq[2*Np+n] = a2; + Bq[2*Np+n] = b2; + + //............................................... + // q = 2 + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + Aq[3*Np+n] = a1; + Bq[3*Np+n] = b1; + Aq[4*Np+n] = a2; + Bq[4*Np+n] = b2; + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + Aq[5*Np+n] = a1; + Bq[5*Np+n] = b1; + Aq[6*Np+n] = a2; + Bq[6*Np+n] = b2; + //............................................... + + } + } +} + + +__global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, + double *Phi, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ + + int n,nn,ijk,nread; + int nr1,nr2,nr3,nr4,nr5,nr6; + int nr7,nr8,nr9,nr10; + int nr11,nr12,nr13,nr14; + //int nr15,nr16,nr17,nr18; + double fq; + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + double m3,m5,m7; + double nA,nB; // number density + double a1,b1,a2,b2,nAB,delta; + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + double phi,tau,rho0,rlx_setA,rlx_setB; + + const double mrt_V1=0.05263157894736842; + const double mrt_V2=0.012531328320802; + const double mrt_V3=0.04761904761904762; + const double mrt_V4=0.004594820384294068; + const double mrt_V5=0.01587301587301587; + const double mrt_V6=0.0555555555555555555555555; + const double mrt_V7=0.02777777777777778; + const double mrt_V8=0.08333333333333333; + const double mrt_V9=0.003341687552213868; + const double mrt_V10=0.003968253968253968; + const double mrt_V11=0.01388888888888889; + const double mrt_V12=0.04166666666666666; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s even part of dist) + //fq = dist[nread]; // reading the f2 data into register fq + nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) + fq = dist[nr2]; // reading the f2 data into register fq + rho += fq; + m1 -= 11.0*(fq); + m2 -= 4.0*(fq); + jx -= fq; + m4 += 4.0*(fq); + m9 += 2.0*(fq); + m10 -= 4.0*(fq); + + // q=3 + //nread = neighborList[n+2*Np]; // neighbor 4 + //fq = dist[nread]; + nr3 = neighborList[n+2*Np]; // neighbor 4 + fq = dist[nr3]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy = fq; + m6 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 = fq; + m12 = -2.0*fq; + + // q = 4 + //nread = neighborList[n+3*Np]; // neighbor 3 + //fq = dist[nread]; + nr4 = neighborList[n+3*Np]; // neighbor 3 + fq = dist[nr4]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy -= fq; + m6 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 += fq; + m12 -= 2.0*fq; + + // q=5 + //nread = neighborList[n+4*Np]; + //fq = dist[nread]; + nr5 = neighborList[n+4*Np]; + fq = dist[nr5]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz = fq; + m8 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + + // q = 6 + //nread = neighborList[n+5*Np]; + //fq = dist[nread]; + nr6 = neighborList[n+5*Np]; + fq = dist[nr6]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz -= fq; + m8 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + // q=7 + //nread = neighborList[n+6*Np]; + //fq = dist[nread]; + nr7 = neighborList[n+6*Np]; + fq = dist[nr7]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 = fq; + m16 = fq; + m17 = -fq; + + // q = 8 + //nread = neighborList[n+7*Np]; + //fq = dist[nread]; + nr8 = neighborList[n+7*Np]; + fq = dist[nr8]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 += fq; + m16 -= fq; + m17 += fq; + + // q=9 + //nread = neighborList[n+8*Np]; + //fq = dist[nread]; + nr9 = neighborList[n+8*Np]; + fq = dist[nr9]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 += fq; + m17 += fq; + + // q = 10 + //nread = neighborList[n+9*Np]; + //fq = dist[nread]; + nr10 = neighborList[n+9*Np]; + fq = dist[nr10]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 -= fq; + m17 -= fq; + + // q=11 + //nread = neighborList[n+10*Np]; + //fq = dist[nread]; + nr11 = neighborList[n+10*Np]; + fq = dist[nr11]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 = fq; + m16 -= fq; + m18 = fq; + + // q=12 + //nread = neighborList[n+11*Np]; + //fq = dist[nread]; + nr12 = neighborList[n+11*Np]; + fq = dist[nr12]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 += fq; + m16 += fq; + m18 -= fq; + + // q=13 + //nread = neighborList[n+12*Np]; + //fq = dist[nread]; + nr13 = neighborList[n+12*Np]; + fq = dist[nr13]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 -= fq; + m18 -= fq; + + // q=14 + //nread = neighborList[n+13*Np]; + //fq = dist[nread]; + nr14 = neighborList[n+13*Np]; + fq = dist[nr14]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 += fq; + m18 += fq; + + // q=15 + nread = neighborList[n+14*Np]; + fq = dist[nread]; + //fq = dist[17*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 = fq; + m17 += fq; + m18 -= fq; + + // q=16 + nread = neighborList[n+15*Np]; + fq = dist[nread]; + //fq = dist[8*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 += fq; + m17 -= fq; + m18 += fq; + + // q=17 + //fq = dist[18*Np+n]; + nread = neighborList[n+16*Np]; + fq = dist[nread]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 += fq; + m18 += fq; + + // q=18 + nread = neighborList[n+17*Np]; + fq = dist[nread]; + //fq = dist[9*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 -= fq; + m18 -= fq; + + //........................................................................ + //..............carry out relaxation process.............................. + //..........Toelke, Fruediger et. al. 2006................................ + if (C == 0.0) nx = ny = nz = 0.0; + m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -19*alpha*C - m1); + m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); + m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); + m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); + m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); + m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); + m10 = m10 + rlx_setA*( - m10); + m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); + m12 = m12 + rlx_setA*( - m12); + m13 = m13 + rlx_setA*( (jx*jy/rho0) + 0.5*alpha*C*nx*ny - m13); + m14 = m14 + rlx_setA*( (jy*jz/rho0) + 0.5*alpha*C*ny*nz - m14); + m15 = m15 + rlx_setA*( (jx*jz/rho0) + 0.5*alpha*C*nx*nz - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); + //.................inverse transformation...................................................... + + // q=0 + fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; + dist[n] = fq; + + // q = 1 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10)+0.16666666*Fx; + //nread = neighborList[n+Np]; + dist[nr2] = fq; + + // q=2 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10) - 0.16666666*Fx; + //nread = neighborList[n]; + dist[nr1] = fq; + + // q = 3 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) + 0.16666666*Fy; + //nread = neighborList[n+3*Np]; + dist[nr4] = fq; + + // q = 4 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) - 0.16666666*Fy; + //nread = neighborList[n+2*Np]; + dist[nr3] = fq; + + // q = 5 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) + 0.16666666*Fz; + //nread = neighborList[n+5*Np]; + dist[nr6] = fq; + + // q = 6 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) - 0.16666666*Fz; + //nread = neighborList[n+4*Np]; + dist[nr5] = fq; + + // q = 7 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17) + 0.08333333333*(Fx+Fy); + //nread = neighborList[n+7*Np]; + dist[nr8] = fq; + + // q = 8 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11 + +mrt_V12*m12+0.25*m13+0.125*(m17-m16) - 0.08333333333*(Fx+Fy); + //nread = neighborList[n+6*Np]; + dist[nr7] = fq; + + // q = 9 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17) + 0.08333333333*(Fx-Fy); + //nread = neighborList[n+9*Np]; + dist[nr10] = fq; + + // q = 10 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17)- 0.08333333333*(Fx-Fy); + //nread = neighborList[n+8*Np]; + dist[nr9] = fq; + + // q = 11 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12+0.25*m15+0.125*(m18-m16) + 0.08333333333*(Fx+Fz); + //nread = neighborList[n+11*Np]; + dist[nr12] = fq; + + // q = 12 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+ + mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18) - 0.08333333333*(Fx+Fz); + //nread = neighborList[n+10*Np]; + dist[nr11]= fq; + + // q = 13 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15-0.125*(m16+m18) + 0.08333333333*(Fx-Fz); + //nread = neighborList[n+13*Np]; + dist[nr14] = fq; + + // q= 14 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15+0.125*(m16+m18) - 0.08333333333*(Fx-Fz); + //nread = neighborList[n+12*Np]; + dist[nr13] = fq; + + + // q = 15 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18) + 0.08333333333*(Fy+Fz); + nread = neighborList[n+15*Np]; + dist[nread] = fq; + + // q = 16 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17)- 0.08333333333*(Fy+Fz); + nread = neighborList[n+14*Np]; + dist[nread] = fq; + + + // q = 17 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) + -mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18) + 0.08333333333*(Fy-Fz); + nread = neighborList[n+17*Np]; + dist[nread] = fq; + + // q = 18 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6) + -mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18) - 0.08333333333*(Fy-Fz); + nread = neighborList[n+16*Np]; + dist[nread] = fq; + + // write the velocity + ux = jx / rho0; + uy = jy / rho0; + uz = jz / rho0; + Velocity[n] = ux; + Velocity[Np+n] = uy; + Velocity[2*Np+n] = uz; + + // Instantiate mass transport distributions + // Stationary value - distribution 0 + nAB = 1.0/(nA+nB); + Aq[n] = 0.3333333333333333*nA; + Bq[n] = 0.3333333333333333*nB; + + //............................................... + // q = 0,2,4 + // Cq = {1,0,0}, {0,1,0}, {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nx; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + // q = 1 + //nread = neighborList[n+Np]; + Aq[nr2] = a1; + Bq[nr2] = b1; + // q=2 + //nread = neighborList[n]; + Aq[nr1] = a2; + Bq[nr1] = b2; + + //............................................... + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + // q = 3 + //nread = neighborList[n+3*Np]; + Aq[nr4] = a1; + Bq[nr4] = b1; + // q = 4 + //nread = neighborList[n+2*Np]; + Aq[nr3] = a2; + Bq[nr3] = b2; + + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + // q = 5 + //nread = neighborList[n+5*Np]; + Aq[nr6] = a1; + Bq[nr6] = b1; + // q = 6 + //nread = neighborList[n+4*Np]; + Aq[nr5] = a2; + Bq[nr5] = b2; + //............................................... + } + } +} + +__global__ void dvc_ScaLBL_D3Q19_AAodd_ColorMomentum(int *neighborList, double *dist, double *Den, + double *Velocity, double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int start, int finish, int Np){ + + int n,nread; + double fq; + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + double nA,nB; // number density + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + double phi,tau,rho0,rlx_setA,rlx_setB; + + const double mrt_V1=0.05263157894736842; + const double mrt_V2=0.012531328320802; + const double mrt_V3=0.04761904761904762; + const double mrt_V4=0.004594820384294068; + const double mrt_V5=0.01587301587301587; + const double mrt_V6=0.0555555555555555555555555; + const double mrt_V7=0.02777777777777778; + const double mrt_V8=0.08333333333333333; + const double mrt_V9=0.003341687552213868; + const double mrt_V10=0.003968253968253968; + const double mrt_V11=0.01388888888888889; + const double mrt_V12=0.04166666666666666; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s 10Np => odd part of dist) + fq = dist[nread]; // reading the f1 data into register fq + //fp = dist[10*Np+n]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jx = fq; + m4 = -4.0*fq; + m9 = 2.0*fq; + m10 = -4.0*fq; + + // f2 = dist[10*Np+n]; + nread = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) + fq = dist[nread]; // reading the f2 data into register fq + //fq = dist[Np+n]; + rho += fq; + m1 -= 11.0*(fq); + m2 -= 4.0*(fq); + jx -= fq; + m4 += 4.0*(fq); + m9 += 2.0*(fq); + m10 -= 4.0*(fq); + + // q=3 + nread = neighborList[n+2*Np]; // neighbor 4 + fq = dist[nread]; + //fq = dist[11*Np+n]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy = fq; + m6 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 = fq; + m12 = -2.0*fq; + + // q = 4 + nread = neighborList[n+3*Np]; // neighbor 3 + fq = dist[nread]; + //fq = dist[2*Np+n]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy -= fq; + m6 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 += fq; + m12 -= 2.0*fq; + + // q=5 + nread = neighborList[n+4*Np]; + fq = dist[nread]; + //fq = dist[12*Np+n]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz = fq; + m8 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + + // q = 6 + nread = neighborList[n+5*Np]; + fq = dist[nread]; + //fq = dist[3*Np+n]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz -= fq; + m8 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + // q=7 + nread = neighborList[n+6*Np]; + fq = dist[nread]; + //fq = dist[13*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 = fq; + m16 = fq; + m17 = -fq; + + // q = 8 + nread = neighborList[n+7*Np]; + fq = dist[nread]; + //fq = dist[4*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 += fq; + m16 -= fq; + m17 += fq; + + // q=9 + nread = neighborList[n+8*Np]; + fq = dist[nread]; + //fq = dist[14*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 += fq; + m17 += fq; + + // q = 10 + nread = neighborList[n+9*Np]; + fq = dist[nread]; + //fq = dist[5*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 -= fq; + m17 -= fq; + + // q=11 + nread = neighborList[n+10*Np]; + fq = dist[nread]; + //fq = dist[15*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 = fq; + m16 -= fq; + m18 = fq; + + // q=12 + nread = neighborList[n+11*Np]; + fq = dist[nread]; + //fq = dist[6*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 += fq; + m16 += fq; + m18 -= fq; + + // q=13 + nread = neighborList[n+12*Np]; + fq = dist[nread]; + //fq = dist[16*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 -= fq; + m18 -= fq; + + // q=14 + nread = neighborList[n+13*Np]; + fq = dist[nread]; + //fq = dist[7*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 += fq; + m18 += fq; + + // q=15 + nread = neighborList[n+14*Np]; + fq = dist[nread]; + //fq = dist[17*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 = fq; + m17 += fq; + m18 -= fq; + + // q=16 + nread = neighborList[n+15*Np]; + fq = dist[nread]; + //fq = dist[8*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 += fq; + m17 -= fq; + m18 += fq; + + // q=17 + //fq = dist[18*Np+n]; + nread = neighborList[n+16*Np]; + fq = dist[nread]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 += fq; + m18 += fq; + + // q=18 + nread = neighborList[n+17*Np]; + fq = dist[nread]; + //fq = dist[9*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 -= fq; + m18 -= fq; + + //........................................................................ + //..............carry out relaxation process.............................. + //..........Toelke, Fruediger et. al. 2006................................ + if (C == 0.0) nx = ny = nz = 0.0; + m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -alpha*C - m1); + m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); + m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); + m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); + m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); + m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); + m10 = m10 + rlx_setA*( - m10); + m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); + m12 = m12 + rlx_setA*( - m12); + m13 = m13 + rlx_setA*( (jx*jy/rho0) + 0.5*alpha*C*nx*ny - m13); + m14 = m14 + rlx_setA*( (jy*jz/rho0) + 0.5*alpha*C*ny*nz - m14); + m15 = m15 + rlx_setA*( (jx*jz/rho0) + 0.5*alpha*C*nx*nz - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); + //.................inverse transformation...................................................... + + // q=0 + fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; + dist[n] = fq; + + // q = 1 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10)+0.16666666*Fx; + nread = neighborList[n+Np]; + dist[nread] = fq; + + // q=2 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10) - 0.16666666*Fx; + nread = neighborList[n]; + dist[nread] = fq; + + // q = 3 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) + 0.16666666*Fy; + nread = neighborList[n+3*Np]; + dist[nread] = fq; + + // q = 4 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) - 0.16666666*Fy; + nread = neighborList[n+2*Np]; + dist[nread] = fq; + + // q = 5 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) + 0.16666666*Fz; + nread = neighborList[n+5*Np]; + dist[nread] = fq; + + // q = 6 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) - 0.16666666*Fz; + nread = neighborList[n+4*Np]; + dist[nread] = fq; + + // q = 7 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17) + 0.08333333333*(Fx+Fy); + nread = neighborList[n+7*Np]; + dist[nread] = fq; + + // q = 8 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11 + +mrt_V12*m12+0.25*m13+0.125*(m17-m16) - 0.08333333333*(Fx+Fy); + nread = neighborList[n+6*Np]; + dist[nread] = fq; + + // q = 9 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17) + 0.08333333333*(Fx-Fy); + nread = neighborList[n+9*Np]; + dist[nread] = fq; + + // q = 10 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17)- 0.08333333333*(Fx-Fy); + nread = neighborList[n+8*Np]; + dist[nread] = fq; + + // q = 11 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12+0.25*m15+0.125*(m18-m16) + 0.08333333333*(Fx+Fz); + nread = neighborList[n+11*Np]; + dist[nread] = fq; + + // q = 12 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+ + mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18) - 0.08333333333*(Fx+Fz); + nread = neighborList[n+10*Np]; + dist[nread]= fq; + + // q = 13 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15-0.125*(m16+m18) + 0.08333333333*(Fx-Fz); + nread = neighborList[n+13*Np]; + dist[nread] = fq; + + // q= 14 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15+0.125*(m16+m18) - 0.08333333333*(Fx-Fz); + nread = neighborList[n+12*Np]; + dist[nread] = fq; + + + // q = 15 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18) + 0.08333333333*(Fy+Fz); + nread = neighborList[n+15*Np]; + dist[nread] = fq; + + // q = 16 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17)- 0.08333333333*(Fy+Fz); + nread = neighborList[n+14*Np]; + dist[nread] = fq; + + + // q = 17 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) + -mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18) + 0.08333333333*(Fy-Fz); + nread = neighborList[n+17*Np]; + dist[nread] = fq; + + // q = 18 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6) + -mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18) - 0.08333333333*(Fy-Fz); + nread = neighborList[n+16*Np]; + dist[nread] = fq; + + // write the velocity + ux = jx / rho0; + uy = jy / rho0; + uz = jz / rho0; + Velocity[n] = ux; + Velocity[Np+n] = uy; + Velocity[2*Np+n] = uz; + } + } +} + +__global__ void dvc_ScaLBL_D3Q19_AAeven_ColorMomentum(double *dist, double *Den, double *Velocity, + double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int start, int finish, int Np){ + int n; + double fq; + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + double nA,nB; // number density + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + double phi,tau,rho0,rlx_setA,rlx_setB; + + const double mrt_V1=0.05263157894736842; + const double mrt_V2=0.012531328320802; + const double mrt_V3=0.04761904761904762; + const double mrt_V4=0.004594820384294068; + const double mrt_V5=0.01587301587301587; + const double mrt_V6=0.0555555555555555555555555; + const double mrt_V7=0.02777777777777778; + const double mrt_V8=0.08333333333333333; + const double mrt_V9=0.003341687552213868; + const double mrt_V10=0.003968253968253968; + const double mrt_V11=0.01388888888888889; + const double mrt_V12=0.04166666666666666; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + Aq[1*Np+n] = a1; + Bq[1*Np+n] = b1; + Aq[2*Np+n] = a2; + Bq[2*Np+n] = b2; + + //............................................... + // q = 2 + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + Aq[3*Np+n] = a1; + Bq[3*Np+n] = b1; + Aq[4*Np+n] = a2; + Bq[4*Np+n] = b2; + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + Aq[5*Np+n] = a1; + Bq[5*Np+n] = b1; + Aq[6*Np+n] = a2; + Bq[6*Np+n] = b2; + //............................................... + + } + } +} + +__global__ void dvc_ScaLBL_D3Q19_AAodd_ColorMass(int *neighborList, double *Aq, double *Bq, double *Den, + double *Velocity, double *ColorGrad, double beta, int start, int finish, int Np){ + + int n,nread; + // non-conserved moments + double nA,nB; // number density + double a1,b1,a2,b2,nAB,delta; + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + // q = 1 + nread = neighborList[n+Np]; + Aq[nread] = a1; + Bq[nread] = b1; + // q=2 + nread = neighborList[n]; + Aq[nread] = a2; + Bq[nread] = b2; + + //............................................... + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + // q = 3 + nread = neighborList[n+3*Np]; + Aq[nread] = a1; + Bq[nread] = b1; + // q = 4 + nread = neighborList[n+2*Np]; + Aq[nread] = a2; + Bq[nread] = b2; + + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + // q = 5 + nread = neighborList[n+5*Np]; + Aq[nread] = a1; + Bq[nread] = b1; + // q = 6 + nread = neighborList[n+4*Np]; + Aq[nread] = a2; + Bq[nread] = b2; + //............................................... + } + } +} + +__global__ void dvc_ScaLBL_D3Q7_AAodd_PhaseField(int *neighborList, int *Map, double *Aq, double *Bq, + double *Den, double *Phi, int start, int finish, int Np){ + int idx,n,nread; + double fq,nA,nB; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s 1.f){ + nA = 1.0; nB = 0.f; + } + else if (phi < -1.f){ + nB = 1.0; nA = 0.f; + } + else{ + nA=0.5*(phi+1.f); + nB=0.5*(1.f-phi); + } + Den[idx] = nA; + Den[Np+idx] = nB; + + Aq[idx]=0.3333333333333333*nA; + Aq[Np+idx]=0.1111111111111111*nA; + Aq[2*Np+idx]=0.1111111111111111*nA; + Aq[3*Np+idx]=0.1111111111111111*nA; + Aq[4*Np+idx]=0.1111111111111111*nA; + Aq[5*Np+idx]=0.1111111111111111*nA; + Aq[6*Np+idx]=0.1111111111111111*nA; + + Bq[idx]=0.3333333333333333*nB; + Bq[Np+idx]=0.1111111111111111*nB; + Bq[2*Np+idx]=0.1111111111111111*nB; + Bq[3*Np+idx]=0.1111111111111111*nB; + Bq[4*Np+idx]=0.1111111111111111*nB; + Bq[5*Np+idx]=0.1111111111111111*nB; + Bq[6*Np+idx]=0.1111111111111111*nB; + } + } +} + +extern "C" void ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, int Nz, int Slice){ + int GRID = Nx*Ny / 512 + 1; + dvc_ScaLBL_SetSlice_z<<>>(Phi,value,Nx,Ny,Nz,Slice); +} + +extern "C" void ScaLBL_Color_Init(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz){ + dvc_ScaLBL_Color_Init<<>>(ID, Den, Phi, das, dbs, Nx, Ny, Nz); +} + +extern "C" void ScaLBL_Color_InitDistance(char *ID, double *Den, double *Phi, double *Distance, + double das, double dbs, double beta, double xp, int Nx, int Ny, int Nz){ + + dvc_ScaLBL_Color_InitDistance<<>>(ID, Den, Phi, Distance, das, dbs, beta, xp, Nx, Ny, Nz); +} + +extern "C" void ScaLBL_D3Q19_ColorGradient(char *ID, double *phi, double *ColorGrad, int Nx, int Ny, int Nz){ + dvc_ScaLBL_D3Q19_ColorGradient<<>>(ID, phi, ColorGrad, Nx, Ny, Nz); +} + +extern "C" void ColorCollide( char *ID, double *disteven, double *distodd, double *ColorGrad, + double *Velocity, int Nx, int Ny, int Nz,double rlx_setA, double rlx_setB, + double alpha, double beta, double Fx, double Fy, double Fz, bool pBC){ + dvc_ColorCollide<<>>( ID, disteven, distodd, ColorGrad,Velocity, Nx, Ny, Nz,rlx_setA, rlx_setB, + alpha, beta, Fx, Fy, Fz, pBC); + +} + +extern "C" void ScaLBL_D3Q19_ColorCollide( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad, + double *Velocity, int Nx, int Ny, int Nz,double rlx_setA, double rlx_setB, + double alpha, double beta, double Fx, double Fy, double Fz){ + dvc_ScaLBL_D3Q19_ColorCollide<<>>(ID, disteven, distodd, phi, ColorGrad, Velocity, Nx, Ny, Nz, rlx_setA, rlx_setB, + alpha, beta, Fx, Fy, Fz); + +} + +extern "C" void DensityStreamD3Q7(char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity, + double beta, int Nx, int Ny, int Nz, bool pBC){ + + dvc_DensityStreamD3Q7<<>>(ID, Den, Copy, Phi, ColorGrad, Velocity, beta, Nx, Ny, Nz, pBC); +} + +extern "C" void ScaLBL_ComputePhaseField(char *ID, double *Phi, double *Den, int N){ + dvc_ScaLBL_ComputePhaseField<<>>(ID, Phi, Den, N); +} +extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A_odd, double *B_even, double *B_odd, + double *Den, double *Phi, double *ColorGrad, double *Velocity, double beta, int N, bool pBC){ + dvc_ScaLBL_D3Q7_ColorCollideMass<<>>(ID, A_even, A_odd, B_even, B_odd, Den, Phi, ColorGrad, Velocity, beta, N, pBC); +} +// Pressure Boundary Conditions Functions + +extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_Color, hipFuncCachePreferL1); + + dvc_ScaLBL_D3Q19_AAeven_Color<<>>(Map, dist, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_Color: %s \n",hipGetErrorString(err)); + } + +} + +extern "C" void ScaLBL_D3Q19_AAodd_Color(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, + double *Phi, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_Color, hipFuncCachePreferL1); + + dvc_ScaLBL_D3Q19_AAodd_Color<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, Vel, + rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_Color: %s \n",hipGetErrorString(err)); + } + hipProfilerStop(); +} + +extern "C" void ScaLBL_D3Q7_AAodd_PhaseField(int *NeighborList, int *Map, double *Aq, double *Bq, + double *Den, double *Phi, int start, int finish, int Np){ + + dvc_ScaLBL_D3Q7_AAodd_PhaseField<<>>(NeighborList, Map, Aq, Bq, Den, Phi, start, finish, Np); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_AAodd_PhaseField: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq, double *Den, double *Phi, + int start, int finish, int Np){ + + dvc_ScaLBL_D3Q7_AAeven_PhaseField<<>>(Map, Aq, Bq, Den, Phi, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_AAeven_PhaseField: %s \n",hipGetErrorString(err)); + } + +} + +extern "C" void ScaLBL_D3Q19_Gradient(int *Map, double *Phi, double *ColorGrad, int start, int finish, int Np, + int Nx, int Ny, int Nz){ + + int strideY=Nx; + int strideZ=Nx*Ny; + dvc_ScaLBL_D3Q19_Gradient<<>>(Map, Phi, ColorGrad, start, finish, Np, strideY, strideZ); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_ColorGrad: %s \n",hipGetErrorString(err)); + } + +} + +extern "C" void ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np){ + dvc_ScaLBL_PhaseField_Init<<>>(Map, Phi, Den, Aq, Bq, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_PhaseField_Init: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_ColorMomentum(double *dist, double *Den, double *Vel, + double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_ColorMomentum, hipFuncCachePreferL1); + + dvc_ScaLBL_D3Q19_AAeven_ColorMomentum<<>>(dist, Den, Vel, ColorGrad, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_ColorMomentum: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAodd_ColorMomentum(int *d_neighborList, double *dist, double *Den, double *Vel, + double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_ColorMomentum, hipFuncCachePreferL1); + + dvc_ScaLBL_D3Q19_AAodd_ColorMomentum<<>>(d_neighborList, dist, Den, Vel, ColorGrad, + rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_ColorMomentum: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_ColorMass(double *Aq, double *Bq, double *Den, double *Vel, + double *ColorGrad, double beta, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_ColorMass, hipFuncCachePreferL1); + + dvc_ScaLBL_D3Q19_AAeven_ColorMass<<>>(Aq, Bq, Den, Vel, ColorGrad, beta, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_Color: %s \n",hipGetErrorString(err)); + } + +} + +extern "C" void ScaLBL_D3Q19_AAodd_ColorMass(int *d_neighborList, double *Aq, double *Bq, double *Den, double *Vel, + double *ColorGrad, double beta, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_ColorMass, hipFuncCachePreferL1); + + dvc_ScaLBL_D3Q19_AAodd_ColorMass<<>>(d_neighborList, Aq, Bq, Den, Vel, ColorGrad, beta, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_Color: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_Color_BC_z(int *list, int *Map, double *Phi, double *Den, double vA, double vB, int count, int Np){ + int GRID = count / 512 + 1; + dvc_ScaLBL_Color_BC_z<<>>(list, Map, Phi, Den, vA, vB, count, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_Color_BC_z: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_Color_BC_Z(int *list, int *Map, double *Phi, double *Den, double vA, double vB, int count, int Np){ + int GRID = count / 512 + 1; + dvc_ScaLBL_Color_BC_Z<<>>(list, Map, Phi, Den, vA, vB, count, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_Color_BC_Z: %s \n",hipGetErrorString(err)); + } +} + + + diff --git a/hip/CudaExtras.hip b/hip/CudaExtras.hip new file mode 100644 index 00000000..5be72e5a --- /dev/null +++ b/hip/CudaExtras.hip @@ -0,0 +1,34 @@ +// Basic hip functions callable from C/C++ code +#include "hip/hip_runtime.h" + +extern "C" void dvc_AllocateDeviceMemory(void** address, size_t size){ + hipMalloc(address,size); + hipMemset(*address,0,size); +} + +extern "C" void dvc_CopyToDevice(void* dest, void* source, size_t size){ + hipMemcpy(dest,source,size,hipMemcpyHostToDevice); +} + + +extern "C" void dvc_CopyToHost(void* dest, void* source, size_t size){ + hipMemcpy(dest,source,size,hipMemcpyDeviceToHost); +} + +extern "C" void dvc_Barrier(){ + hipDeviceSynchronize(); +} +/* +#if __CUDA_ARCH__ < 600 +__device__ double atomicAdd(double* address, double val) { +unsigned long long int* address_as_ull = (unsigned long long int*)address; unsigned long long int old = *address_as_ull, assumed; + do { + assumed = old; + old = atomicCAS(address_as_ull, assumed, __double_as_longlong(val + __longlong_as_double(assumed))); + // Note: uses integer comparison to avoid hang in case of NaN (since NaN != NaN) + } + while (assumed != old); return __longlong_as_double(old); +} + +#endif +*/ diff --git a/hip/D3Q19.hip b/hip/D3Q19.hip new file mode 100644 index 00000000..1bcfa04d --- /dev/null +++ b/hip/D3Q19.hip @@ -0,0 +1,2645 @@ +#include +#include "hip/hip_runtime.h" +#include + +#define NBLOCKS 1024 +#define NTHREADS 256 + +/* +1. constants that are known at compile time should be defined using preprocessor macros (e.g. #define) or via C/C++ const variables at global/file scope. +2. Usage of __constant__ memory may be beneficial for programs who use certain values that don't change for the duration of the kernel and for which certain access patterns are present (e.g. all threads access the same value at the same time). This is not better or faster than constants that satisfy the requirements of item 1 above. +3. If the number of choices to be made by a program are relatively small in number, and these choices affect kernel execution, one possible approach for additional compile-time optimization would be to use templated code/kernels + */ + +__constant__ __device__ double mrt_V1=0.05263157894736842; +__constant__ __device__ double mrt_V2=0.012531328320802; +__constant__ __device__ double mrt_V3=0.04761904761904762; +__constant__ __device__ double mrt_V4=0.004594820384294068; +__constant__ __device__ double mrt_V5=0.01587301587301587; +__constant__ __device__ double mrt_V6=0.0555555555555555555555555; +__constant__ __device__ double mrt_V7=0.02777777777777778; +__constant__ __device__ double mrt_V8=0.08333333333333333; +__constant__ __device__ double mrt_V9=0.003341687552213868; +__constant__ __device__ double mrt_V10=0.003968253968253968; +__constant__ __device__ double mrt_V11=0.01388888888888889; +__constant__ __device__ double mrt_V12=0.04166666666666666; + + +// functionality for parallel reduction in Flux BC routines -- probably should be re-factored to another location +// functions copied from https://devblogs.nvidia.com/parallelforall/faster-parallel-reductions-kepler/ + +//__shared__ double Transform[722]= +// {}; + +#if !defined(__CUDA_ARCH__) || __CUDA_ARCH__ >= 600 +#else +__device__ double atomicAdd(double* address, double val) { + unsigned long long int* address_as_ull = (unsigned long long int*)address; + unsigned long long int old = *address_as_ull, assumed; + + do { + assumed = old; + old = atomicCAS(address_as_ull, assumed, __double_as_longlong(val+__longlong_as_double(assumed))); + } while (assumed != old); + return __longlong_as_double(old); +} +#endif + +using namespace cooperative_groups; +__device__ double reduce_sum(thread_group g, double *temp, double val) +{ + int lane = g.thread_rank(); + + // Each iteration halves the number of active threads + // Each thread adds its partial sum[i] to sum[lane+i] + for (int i = g.size() / 2; i > 0; i /= 2) + { + temp[lane] = val; + g.sync(); // wait for all threads to store + if(lane 0; offset /= 2) + val += __shfl_down_sync(0xFFFFFFFF, val, offset, 32); + return val; +} + +__inline__ __device__ +double blockReduceSum(double val) { + + static __shared__ double shared[32]; // Shared mem for 32 partial sums + int lane = threadIdx.x % warpSize; + int wid = threadIdx.x / warpSize; + + val = warpReduceSum(val); // Each warp performs partial reduction + + if (lane==0) shared[wid]=val; // Write reduced value to shared memory + + __syncthreads(); // Wait for all partial reductions + + //read from shared memory only if that warp existed + val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0; + + if (wid==0) val = warpReduceSum(val); //Final reduce within first warp + + return val; +} + +__global__ void deviceReduceKernel(double *in, double* out, int N) { + double sum = 0; + //reduce multiple elements per thread + for (int i = blockIdx.x * blockDim.x + threadIdx.x; + i < N; + i += blockDim.x * gridDim.x) { + sum += in[i]; + } + sum = blockReduceSum(sum); + if (threadIdx.x==0) + out[blockIdx.x]=sum; +} + +__global__ void dvc_ScaLBL_D3Q19_Pack(int q, int *list, int start, int count, double *sendbuf, double *dist, int N){ + //.................................................................................... + // Pack distribution q into the send buffer for the listed lattice sites + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int idx,n; + idx = blockIdx.x*blockDim.x + threadIdx.x; + if (idx 0 ){ + f_even[n] = 0.3333333333333333; + f_odd[n] = 0.055555555555555555; //double(100*n)+1.f; + f_even[N+n] = 0.055555555555555555; //double(100*n)+2.f; + f_odd[N+n] = 0.055555555555555555; //double(100*n)+3.f; + f_even[2*N+n] = 0.055555555555555555; //double(100*n)+4.f; + f_odd[2*N+n] = 0.055555555555555555; //double(100*n)+5.f; + f_even[3*N+n] = 0.055555555555555555; //double(100*n)+6.f; + f_odd[3*N+n] = 0.0277777777777778; //double(100*n)+7.f; + f_even[4*N+n] = 0.0277777777777778; //double(100*n)+8.f; + f_odd[4*N+n] = 0.0277777777777778; //double(100*n)+9.f; + f_even[5*N+n] = 0.0277777777777778; //double(100*n)+10.f; + f_odd[5*N+n] = 0.0277777777777778; //double(100*n)+11.f; + f_even[6*N+n] = 0.0277777777777778; //double(100*n)+12.f; + f_odd[6*N+n] = 0.0277777777777778; //double(100*n)+13.f; + f_even[7*N+n] = 0.0277777777777778; //double(100*n)+14.f; + f_odd[7*N+n] = 0.0277777777777778; //double(100*n)+15.f; + f_even[8*N+n] = 0.0277777777777778; //double(100*n)+16.f; + f_odd[8*N+n] = 0.0277777777777778; //double(100*n)+17.f; + f_even[9*N+n] = 0.0277777777777778; //double(100*n)+18.f; + } + else{ + for(int q=0; q<9; q++){ + f_even[q*N+n] = -1.0; + f_odd[q*N+n] = -1.0; + } + f_even[9*N+n] = -1.0; + } + } + } +} + +__global__ void dvc_ScaLBL_D3Q19_AA_Init(double *f_even, double *f_odd, int Np) +{ + int n; + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s 10Np => odd part of dist) + fq = dist[nread]; // reading the f1 data into register fq + //fp = dist[10*Np+n]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jx = fq; + m4 = -4.0*fq; + m9 = 2.0*fq; + m10 = -4.0*fq; + + // f2 = dist[10*Np+n]; + nread = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) + fq = dist[nread]; // reading the f2 data into register fq + //fq = dist[Np+n]; + rho += fq; + m1 -= 11.0*(fq); + m2 -= 4.0*(fq); + jx -= fq; + m4 += 4.0*(fq); + m9 += 2.0*(fq); + m10 -= 4.0*(fq); + + // q=3 + nread = neighborList[n+2*Np]; // neighbor 4 + fq = dist[nread]; + //fq = dist[11*Np+n]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy = fq; + m6 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 = fq; + m12 = -2.0*fq; + + // q = 4 + nread = neighborList[n+3*Np]; // neighbor 3 + fq = dist[nread]; + //fq = dist[2*Np+n]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy -= fq; + m6 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 += fq; + m12 -= 2.0*fq; + + // q=5 + nread = neighborList[n+4*Np]; + fq = dist[nread]; + //fq = dist[12*Np+n]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz = fq; + m8 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + + // q = 6 + nread = neighborList[n+5*Np]; + fq = dist[nread]; + //fq = dist[3*Np+n]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz -= fq; + m8 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + // q=7 + nread = neighborList[n+6*Np]; + fq = dist[nread]; + //fq = dist[13*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 = fq; + m16 = fq; + m17 = -fq; + + // q = 8 + nread = neighborList[n+7*Np]; + fq = dist[nread]; + //fq = dist[4*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 += fq; + m16 -= fq; + m17 += fq; + + // q=9 + nread = neighborList[n+8*Np]; + fq = dist[nread]; + //fq = dist[14*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 += fq; + m17 += fq; + + // q = 10 + nread = neighborList[n+9*Np]; + fq = dist[nread]; + //fq = dist[5*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 -= fq; + m17 -= fq; + + // q=11 + nread = neighborList[n+10*Np]; + fq = dist[nread]; + //fq = dist[15*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 = fq; + m16 -= fq; + m18 = fq; + + // q=12 + nread = neighborList[n+11*Np]; + fq = dist[nread]; + //fq = dist[6*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 += fq; + m16 += fq; + m18 -= fq; + + // q=13 + nread = neighborList[n+12*Np]; + fq = dist[nread]; + //fq = dist[16*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 -= fq; + m18 -= fq; + + // q=14 + nread = neighborList[n+13*Np]; + fq = dist[nread]; + //fq = dist[7*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 += fq; + m18 += fq; + + // q=15 + nread = neighborList[n+14*Np]; + fq = dist[nread]; + //fq = dist[17*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 = fq; + m17 += fq; + m18 -= fq; + + // q=16 + nread = neighborList[n+15*Np]; + fq = dist[nread]; + //fq = dist[8*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 += fq; + m17 -= fq; + m18 += fq; + + // q=17 + //fq = dist[18*Np+n]; + nread = neighborList[n+16*Np]; + fq = dist[nread]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 += fq; + m18 += fq; + + // q=18 + nread = neighborList[n+17*Np]; + fq = dist[nread]; + //fq = dist[9*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 -= fq; + m18 -= fq; + + //..............incorporate external force................................................ + //..............carry out relaxation process............................................... + m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) - m1); + m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho) - m2); + m4 = m4 + rlx_setB*((-0.6666666666666666*jx) - m4); + m6 = m6 + rlx_setB*((-0.6666666666666666*jy) - m6); + m8 = m8 + rlx_setB*((-0.6666666666666666*jz) - m8); + m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) - m9); + m10 = m10 + rlx_setA*(-0.5*((2*jx*jx-jy*jy-jz*jz)/rho) - m10); + m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) - m11); + m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho) - m12); + m13 = m13 + rlx_setA*((jx*jy/rho) - m13); + m14 = m14 + rlx_setA*((jy*jz/rho) - m14); + m15 = m15 + rlx_setA*((jx*jz/rho) - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); + //....................................................................................................... + //.................inverse transformation...................................................... + + // q=0 + fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; + dist[n] = fq; + + // q = 1 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10)+0.16666666*Fx; + nread = neighborList[n+Np]; + dist[nread] = fq; + + // q=2 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10) - 0.16666666*Fx; + nread = neighborList[n]; + dist[nread] = fq; + + // q = 3 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) + 0.16666666*Fy; + nread = neighborList[n+3*Np]; + dist[nread] = fq; + + // q = 4 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) - 0.16666666*Fy; + nread = neighborList[n+2*Np]; + dist[nread] = fq; + + // q = 5 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) + 0.16666666*Fz; + nread = neighborList[n+5*Np]; + dist[nread] = fq; + + // q = 6 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) - 0.16666666*Fz; + nread = neighborList[n+4*Np]; + dist[nread] = fq; + + // q = 7 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+mrt_V7*m9+mrt_V11*m10+ + mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17) + 0.08333333333*(Fx+Fy); + + nread = neighborList[n+7*Np]; + dist[nread] = fq; + + // q = 8 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11 + +mrt_V12*m12+0.25*m13+0.125*(m17-m16) - 0.08333333333*(Fx+Fy); + nread = neighborList[n+6*Np]; + dist[nread] = fq; + + // q = 9 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+mrt_V7*m9+mrt_V11*m10+ + mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17) + 0.08333333333*(Fx-Fy); + nread = neighborList[n+9*Np]; + dist[nread] = fq; + + // q = 10 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+mrt_V7*m9+mrt_V11*m10+ + mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17)- 0.08333333333*(Fx-Fy); + nread = neighborList[n+8*Np]; + dist[nread] = fq; + + // q = 11 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12+0.25*m15+0.125*(m18-m16) + 0.08333333333*(Fx+Fz); + nread = neighborList[n+11*Np]; + dist[nread] = fq; + + // q = 12 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+ + mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18) - 0.08333333333*(Fx+Fz); + nread = neighborList[n+10*Np]; + dist[nread]= fq; + + // q = 13 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15-0.125*(m16+m18) + 0.08333333333*(Fx-Fz); + nread = neighborList[n+13*Np]; + dist[nread] = fq; + + // q= 14 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15+0.125*(m16+m18) - 0.08333333333*(Fx-Fz); + nread = neighborList[n+12*Np]; + dist[nread] = fq; + + + // q = 15 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18) + 0.08333333333*(Fy+Fz); + nread = neighborList[n+15*Np]; + dist[nread] = fq; + + // q = 16 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17)- 0.08333333333*(Fy+Fz); + nread = neighborList[n+14*Np]; + dist[nread] = fq; + + + // q = 17 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) + -mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18) + 0.08333333333*(Fy-Fz); + nread = neighborList[n+17*Np]; + dist[nread] = fq; + + // q = 18 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6) + -mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18) - 0.08333333333*(Fy-Fz); + nread = neighborList[n+16*Np]; + dist[nread] = fq; + + } + } +} + + +//__launch_bounds__(512,1) +__global__ void +dvc_ScaLBL_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx, double Fy, double Fz) { + + int n; + double fq; + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s 0){ + //.......Back out the 3-D indices for node n.............. + k = n/(Nx*Ny); + j = (n-Nx*Ny*k)/Nx; + i = n-Nx*Ny*k-Nx*j; + //........................................................................ + // Retrieve even distributions from the local node (swap convention) + // f0 = disteven[n]; // Does not particupate in streaming + f1 = distodd[n]; + f3 = distodd[N+n]; + f5 = distodd[2*N+n]; + f7 = distodd[3*N+n]; + f9 = distodd[4*N+n]; + f11 = distodd[5*N+n]; + f13 = distodd[6*N+n]; + f15 = distodd[7*N+n]; + f17 = distodd[8*N+n]; + //........................................................................ + + //........................................................................ + // Retrieve odd distributions from neighboring nodes (swap convention) + //........................................................................ + nn = n+1; // neighbor index (pull convention) + if (!(i+1 0.0){ + distodd[n] = f2; + disteven[N+nn] = f1; + } + //} + //........................................................................ + nn = n+Nx; // neighbor index (pull convention) + if (!(j+1 0.0){ + distodd[N+n] = f4; + disteven[2*N+nn] = f3; + // } + } + //........................................................................ + nn = n+Nx*Ny; // neighbor index (pull convention) + if (!(k+1 0.0){ + distodd[2*N+n] = f6; + disteven[3*N+nn] = f5; + // } + } + //........................................................................ + nn = n+Nx+1; // neighbor index (pull convention) + if (!(i+1 0.0){ + distodd[3*N+n] = f8; + disteven[4*N+nn] = f7; + // } + } + //........................................................................ + nn = n-Nx+1; // neighbor index (pull convention) + if (!(i+1 0.0){ + distodd[4*N+n] = f10; + disteven[5*N+nn] = f9; + // } + } + //........................................................................ + nn = n+Nx*Ny+1; // neighbor index (pull convention) + if (!(i+1 0.0){ + distodd[5*N+n] = f12; + disteven[6*N+nn] = f11; + // } + } + //........................................................................ + nn = n-Nx*Ny+1; // neighbor index (pull convention) + if (!(i+1 0.0){ + distodd[6*N+n] = f14; + disteven[7*N+nn] = f13; + // } + } + //........................................................................ + nn = n+Nx*Ny+Nx; // neighbor index (pull convention) + if (!(j+1 0.0){ + distodd[7*N+n] = f16; + disteven[8*N+nn] = f15; + // } + } + //........................................................................ + nn = n-Nx*Ny+Nx; // neighbor index (pull convention) + if (!(j+1 0.0){ + distodd[8*N+n] = f18; + disteven[9*N+nn] = f17; + // } + } + //........................................................................ + + } + } + } +} + + +__global__ void dvc_ScaLBL_D3Q19_Momentum(double *dist, double *vel, int N) +{ + int n; + // distributions + double f1,f2,f3,f4,f5,f6,f7,f8,f9; + double f10,f11,f12,f13,f14,f15,f16,f17,f18; + double vx,vy,vz; + + int S = N/NBLOCKS/NTHREADS + 1; + for (int s=0; s 0 ){ + f_even[n] = 0 + 0.01*0; + f_odd[n] = 0+ 0.01*1; //double(100*n)+1.f; + f_even[N+n] = 1+ 0.01*2; //double(100*n)+2.f; + f_odd[N+n] = 1+ 0.01*3; //double(100*n)+3.f; + f_even[2*N+n] = 2+ 0.01*4; //double(100*n)+4.f; + f_odd[2*N+n] = 2+ 0.01*5; //double(100*n)+5.f; + f_even[3*N+n] = 3+ 0.01*6; //double(100*n)+6.f; + f_odd[3*N+n] = 3+ 0.01*7; //double(100*n)+7.f; + f_even[4*N+n] = 4+ 0.01*8; //double(100*n)+8.f; + f_odd[4*N+n] = 4+ 0.01*9; //double(100*n)+9.f; + f_even[5*N+n] = 5+ 0.01*10; //double(100*n)+10.f; + f_odd[5*N+n] = 5+ 0.01*11; //double(100*n)+11.f; + f_even[6*N+n] = 6+ 0.01*12; //double(100*n)+12.f; + f_odd[6*N+n] = 6+ 0.01*13; //double(100*n)+13.f; + f_even[7*N+n] = 7+ 0.01*14; //double(100*n)+14.f; + f_odd[7*N+n] = 7+ 0.01*15; //double(100*n)+15.f; + f_even[8*N+n] = 8+ 0.01*16; //double(100*n)+16.f; + f_odd[8*N+n] = 8+ 0.01*17; //double(100*n)+17.f; + f_even[9*N+n] = 9+ 0.01*18; //double(100*n)+18.f; + } + else{ + for(int q=0; q<9; q++){ + f_even[q*N+n] = -1.0; + f_odd[q*N+n] = -1.0; + } + f_even[9*N+n] = -1.0; + } + } + } +} + + +//************************************************************************* + +//extern "C" void ScaLBL_D3Q19_MapRecv(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count, +// int *d3q19_recvlist, int Nx, int Ny, int Nz){ +// int GRID = count / 512 + 1; +// dvc_ScaLBL_D3Q19_Unpack <<>>(q, Cqx, Cqy, Cqz, list, start, count, d3q19_recvlist, Nx, Ny, Nz); +//} + +extern "C" void ScaLBL_D3Q19_Pack(int q, int *list, int start, int count, double *sendbuf, double *dist, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_Pack <<>>(q, list, start, count, sendbuf, dist, N); +} + +extern "C" void ScaLBL_D3Q19_Unpack(int q, int *list, int start, int count, double *recvbuf, double *dist, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_Unpack <<>>(q, list, start, count, recvbuf, dist, N); +} +//************************************************************************* + +extern "C" void ScaLBL_D3Q19_AA_Init(double *f_even, double *f_odd, int Np){ + dvc_ScaLBL_D3Q19_AA_Init<<>>(f_even, f_odd, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AA_Init: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_Init(double *dist, int Np){ + dvc_ScaLBL_D3Q19_Init<<>>(dist, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AA_Init: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_Swap(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz){ + dvc_ScaLBL_D3Q19_Swap<<>>(ID, disteven, distodd, Nx, Ny, Nz); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Swap: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_Swap_Compact(int *neighborList, double *disteven, double *distodd, int Np) +{ + + const int Q = 9; + // hipStream_t streams[Q]; + // Launch the swap operation as different kernels + for (int q=0; q>>(neighborList, disteven, distodd, Np, q); + } + // cpu should wait for all kernels to finish (to avoid launch of dependent kernels) + //hipDeviceSynchronize(); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Swap: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_Compact(char * ID, double *d_dist, int Np) { + hipFuncSetCacheConfig(dvc_ScaLBL_AAeven_Compact, hipFuncCachePreferL1); + dvc_ScaLBL_AAeven_Compact<<>>(ID, d_dist, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Init: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAodd_Compact(char * ID, int *d_neighborList, double *d_dist, int Np) { + hipFuncSetCacheConfig(dvc_ScaLBL_AAodd_Compact, hipFuncCachePreferL1); + dvc_ScaLBL_AAodd_Compact<<>>(ID,d_neighborList, d_dist,Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Init: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np){ + + dvc_ScaLBL_D3Q19_Momentum<<>>(dist, vel, Np); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Velocity: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_Pressure(double *fq, double *Pressure, int Np){ + dvc_ScaLBL_D3Q19_Pressure<<< NBLOCKS,NTHREADS >>>(fq, Pressure, Np); +} + +extern "C" void ScaLBL_D3Q19_Velocity_BC_z(double *disteven, double *distodd, double uz,int Nx, int Ny, int Nz){ + int GRID = Nx*Ny / 512 + 1; + dvc_D3Q19_Velocity_BC_z<<>>(disteven,distodd, uz, Nx, Ny, Nz); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Velocity_BC_z: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_Velocity_BC_Z(double *disteven, double *distodd, double uz, int Nx, int Ny, int Nz, int outlet){ + int GRID = Nx*Ny / 512 + 1; + dvc_D3Q19_Velocity_BC_Z<<>>(disteven, distodd, uz, Nx, Ny, Nz, outlet); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Velocity_BC_Z: %s \n",hipGetErrorString(err)); + } +} + +extern "C" double ScaLBL_D3Q19_Flux_BC_z(double *disteven, double *distodd, double flux,int Nx, int Ny, int Nz){ + + int GRID = Nx*Ny / 512 + 1; + + // IMPORTANT -- this routine may fail if Nx*Ny > 512*512 + if (Nx*Ny > 512*512){ + printf("WARNING (ScaLBL_D3Q19_Flux_BC_z): CUDA reduction operation may fail if Nx*Ny > 512*512"); + } + + // Allocate memory to store the sums + double din; + double sum[1]; + double *dvcsum; + int sharedBytes = NTHREADS*sizeof(double); + hipMalloc((void **)&dvcsum,sizeof(double)*Nx*Ny); + hipMemset(dvcsum,0,sizeof(double)*Nx*Ny); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Flux_BC_z (memory allocation): %s \n",hipGetErrorString(err)); + } + + // compute the local flux and store the result + dvc_D3Q19_Flux_BC_z<<>>(disteven, distodd, flux, dvcsum, Nx, Ny, Nz); + + err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Flux_BC_z (flux calculation, step 1): %s \n",hipGetErrorString(err)); + } + + // Now read the total flux + hipMemcpy(&sum[0],dvcsum,sizeof(double),hipMemcpyDeviceToHost); + din=sum[0]; + + err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Flux_BC_z (flux calculation, step 2): %s \n",hipGetErrorString(err)); + } + + // free the memory needed for reduction + hipFree(dvcsum); + + return din; +} + + +extern "C" void ScaLBL_D3Q19_AAeven_Pressure_BC_z(int *list, double *dist, double din, int count, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAeven_Pressure_BC_z<<>>(list, dist, din, count, N); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_Pressure_BC_z (kernel): %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_Pressure_BC_Z(int *list, double *dist, double dout, int count, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAeven_Pressure_BC_Z<<>>(list, dist, dout, count, N); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_Pressure_BC_Z (kernel): %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAodd_Pressure_BC_z(int *neighborList, int *list, double *dist, double din, int count, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAodd_Pressure_BC_z<<>>(neighborList, list, dist, din, count, N); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_Pressure_BC_z (kernel): %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAodd_Pressure_BC_Z(int *neighborList, int *list, double *dist, double dout, int count, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAodd_Pressure_BC_Z<<>>(neighborList, list, dist, dout, count, N); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_Pressure_BC_Z (kernel): %s \n",hipGetErrorString(err)); + } +} + + +extern "C" double ScaLBL_D3Q19_AAeven_Flux_BC_z(int *list, double *dist, double flux, double area, + int count, int N){ + + int GRID = count / 512 + 1; + + // IMPORTANT -- this routine may fail if Nx*Ny > 512*512 + if (count > 512*512){ + printf("WARNING (ScaLBL_D3Q19_Flux_BC_Z): CUDA reduction operation may fail if count > 512*512"); + } + + // Allocate memory to store the sums + double din; + double sum[1]; + double *dvcsum; + hipMalloc((void **)&dvcsum,sizeof(double)*count); + hipMemset(dvcsum,0,sizeof(double)*count); + int sharedBytes = 512*sizeof(double); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_Flux_BC_z (memory allocation): %s \n",hipGetErrorString(err)); + } + + // compute the local flux and store the result + dvc_ScaLBL_D3Q19_AAeven_Flux_BC_z<<>>(list, dist, flux, area, dvcsum, count, N); + err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_Flux_BC_z (kernel): %s \n",hipGetErrorString(err)); + } + + // Now read the total flux + hipMemcpy(&sum[0],dvcsum,sizeof(double),hipMemcpyDeviceToHost); + din=sum[0]; + err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_Flux_BC_z (reduction): %s \n",hipGetErrorString(err)); + } + + // free the memory needed for reduction + hipFree(dvcsum); + + return din; +} + +extern "C" double ScaLBL_D3Q19_AAodd_Flux_BC_z(int *neighborList, int *list, double *dist, double flux, + double area, int count, int N){ + + int GRID = count / 512 + 1; + + // IMPORTANT -- this routine may fail if Nx*Ny > 512*512 + if (count > 512*512){ + printf("WARNING (ScaLBL_D3Q19_AAodd_Flux_BC_z): CUDA reduction operation may fail if count > 512*512"); + } + + // Allocate memory to store the sums + double din; + double sum[1]; + double *dvcsum; + hipMalloc((void **)&dvcsum,sizeof(double)*count); + hipMemset(dvcsum,0,sizeof(double)*count); + int sharedBytes = 512*sizeof(double); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_Flux_BC_z (memory allocation): %s \n",hipGetErrorString(err)); + } + + // compute the local flux and store the result + dvc_ScaLBL_D3Q19_AAodd_Flux_BC_z<<>>(neighborList, list, dist, flux, area, dvcsum, count, N); + err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_Flux_BC_z (kernel): %s \n",hipGetErrorString(err)); + } + // Now read the total flux + hipMemcpy(&sum[0],dvcsum,sizeof(double),hipMemcpyDeviceToHost); + din=sum[0]; + err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_Flux_BC_z (reduction): %s \n",hipGetErrorString(err)); + } + + // free the memory needed for reduction + hipFree(dvcsum); + + return din; +} + +extern "C" double ScaLBL_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double flux, int Nx, int Ny, int Nz, int outlet){ + + int GRID = Nx*Ny / 512 + 1; + + // IMPORTANT -- this routine may fail if Nx*Ny > 512*512 + if (Nx*Ny > 512*512){ + printf("WARNING (ScaLBL_D3Q19_Flux_BC_Z): CUDA reduction operation may fail if Nx*Ny > 512*512"); + } + + // Allocate memory to store the sums + double dout; + double sum[1]; + double *dvcsum; + hipMalloc((void **)&dvcsum,sizeof(double)*Nx*Ny); + hipMemset(dvcsum,0,sizeof(double)*Nx*Ny); + + // compute the local flux and store the result + dvc_D3Q19_Flux_BC_Z<<>>(disteven, distodd, flux, dvcsum, Nx, Ny, Nz, outlet); + + // Now read the total flux + hipMemcpy(&sum[0],dvcsum,sizeof(double),hipMemcpyDeviceToHost); + + // free the memory needed for reduction + + dout = sum[0]; + + hipFree(dvcsum); + + return dout; + +} + +extern "C" double deviceReduce(double *in, double* out, int N) { + int threads = 512; + int blocks = min((N + threads - 1) / threads, 1024); + + double sum = 0.f; + deviceReduceKernel<<>>(in, out, N); + deviceReduceKernel<<<1, 1024>>>(out, out, blocks); + return sum; +} + +// +//extern "C" void ScaLBL_D3Q19_Pressure_BC_Z(int *list, double *dist, double dout, int count, int Np){ +// int GRID = count / 512 + 1; +// dvc_ScaLBL_D3Q19_Pressure_BC_Z<<>>(disteven, distodd, dout, Nx, Ny, Nz, outlet); +//} + +extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx, + double Fy, double Fz){ + + dvc_ScaLBL_AAeven_MRT<<>>(dist,start,finish,Np,rlx_setA,rlx_setB,Fx,Fy,Fz); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_MRT: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAodd_MRT(int *neighborlist, double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx, + double Fy, double Fz){ + + dvc_ScaLBL_AAodd_MRT<<>>(neighborlist,dist,start,finish,Np,rlx_setA,rlx_setB,Fx,Fy,Fz); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_MRT: %s \n",hipGetErrorString(err)); + } +} + diff --git a/hip/D3Q7.hip b/hip/D3Q7.hip new file mode 100644 index 00000000..16863fec --- /dev/null +++ b/hip/D3Q7.hip @@ -0,0 +1,246 @@ +// GPU Functions for D3Q7 Lattice Boltzmann Methods + +#define NBLOCKS 560 +#define NTHREADS 128 + +__global__ void dvc_ScaLBL_Scalar_Pack(int *list, int count, double *sendbuf, double *Data, int N){ + //.................................................................................... + // Pack distribution q into the send buffer for the listed lattice sites + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int idx,n; + idx = blockIdx.x*blockDim.x + threadIdx.x; + if (idx 0){ + value = Den[n]; + f_even[n] = 0.3333333333333333*value; + f_odd[n] = 0.1111111111111111*value; //double(100*n)+1.f; + f_even[N+n] = 0.1111111111111111*value; //double(100*n)+2.f; + f_odd[N+n] = 0.1111111111111111*value; //double(100*n)+3.f; + f_even[2*N+n] = 0.1111111111111111*value; //double(100*n)+4.f; + f_odd[2*N+n] = 0.1111111111111111*value; //double(100*n)+5.f; + f_even[3*N+n] = 0.1111111111111111*value; //double(100*n)+6.f; + } + else{ + for(int q=0; q<3; q++){ + f_even[q*N+n] = -1.0; + f_odd[q*N+n] = -1.0; + } + f_even[3*N+n] = -1.0; + } + } + } +} + +//************************************************************************* +__global__ void dvc_ScaLBL_D3Q7_Swap(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz) +{ + int i,j,k,n,nn,N; + // distributions + double f1,f2,f3,f4,f5,f6; + char id; + N = Nx*Ny*Nz; + + int S = N/NBLOCKS/NTHREADS + 1; + for (int s=0; s 0){ + //.......Back out the 3-D indices for node n.............. + k = n/(Nx*Ny); + j = (n-Nx*Ny*k)/Nx; + i = n-Nx*Ny*k-Nx*j; + //........................................................................ + // Retrieve even distributions from the local node (swap convention) + // f0 = disteven[n]; // Does not particupate in streaming + f1 = distodd[n]; + f3 = distodd[N+n]; + f5 = distodd[2*N+n]; + //........................................................................ + + //........................................................................ + // Retrieve odd distributions from neighboring nodes (swap convention) + //........................................................................ + nn = n+1; // neighbor index (pull convention) + if (!(i+1 0 ){ + // Read the distributions + f0 = disteven[n]; + f2 = disteven[N+n]; + f4 = disteven[2*N+n]; + f6 = disteven[3*N+n]; + f1 = distodd[n]; + f3 = distodd[N+n]; + f5 = distodd[2*N+n]; + // Compute the density + Den[n] = f0+f1+f2+f3+f4+f5+f6; + } + } + } +} + +extern "C" void ScaLBL_D3Q7_Unpack(int q, int *list, int start, int count, double *recvbuf, double *dist, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q7_Unpack <<>>(q, list, start, count, recvbuf, dist, N); +} + +extern "C" void ScaLBL_Scalar_Pack(int *list, int count, double *sendbuf, double *Data, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_Scalar_Pack <<>>(list, count, sendbuf, Data, N); +} + +extern "C" void ScaLBL_Scalar_Unpack(int *list, int count, double *recvbuf, double *Data, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_Scalar_Unpack <<>>(list, count, recvbuf, Data, N); +} +extern "C" void ScaLBL_PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_PackDenD3Q7 <<>>(list, count, sendbuf, number, Data, N); +} + +extern "C" void ScaLBL_UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N){ + int GRID = count / 512 + 1; + dvc_ScaLBL_UnpackDenD3Q7 <<>>(list, count, recvbuf, number, Data, N); +} + +extern "C" void ScaLBL_D3Q7_Init(char *ID, double *f_even, double *f_odd, double *Den, int Nx, int Ny, int Nz){ + dvc_ScaLBL_D3Q7_Init <<>>(ID, f_even, f_odd, Den, Nx, Ny, Nz); +} + +extern "C" void ScaLBL_D3Q7_Swap(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz){ + dvc_ScaLBL_D3Q7_Swap <<>>(ID, disteven, distodd, Nx, Ny, Nz); +} + +extern "C" void ScaLBL_D3Q7_Density(char *ID, double *disteven, double *distodd, double *Den, + int Nx, int Ny, int Nz){ + dvc_ScaLBL_D3Q7_Density <<>>(ID, disteven, distodd, Den, Nx, Ny, Nz); +} + diff --git a/hip/Extras.hip b/hip/Extras.hip new file mode 100644 index 00000000..6d2a65e3 --- /dev/null +++ b/hip/Extras.hip @@ -0,0 +1,62 @@ +// Basic hip functions callable from C/C++ code +#include "hip/hip_runtime.h" +#include + +extern "C" int ScaLBL_SetDevice(int rank){ + int n_devices; + //int local_rank = atoi(getenv("MV2_COMM_WORLD_LOCAL_RANK")); + hipGetDeviceCount(&n_devices); + //int device = local_rank % n_devices; + int device = rank % n_devices; + hipSetDevice(device); + if (rank < n_devices) printf("MPI rank=%i will use GPU ID %i / %i \n",rank,device,n_devices); + return device; +} + +extern "C" void ScaLBL_AllocateDeviceMemory(void** address, size_t size){ + hipMalloc(address,size); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("Error in hipMalloc: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_FreeDeviceMemory(void* pointer){ + hipFree(pointer); +} + +extern "C" void ScaLBL_CopyToDevice(void* dest, const void* source, size_t size){ + hipMemcpy(dest,source,size,hipMemcpyHostToDevice); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("Error in hipMemcpy (host->device): %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_AllocateZeroCopy(void** address, size_t size){ + //hipMallocHost(address,size); + hipMalloc(address,size); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("Error in hipMallocHost: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_CopyToZeroCopy(void* dest, const void* source, size_t size){ + hipMemcpy(dest,source,size,hipMemcpyHostToDevice); + hipError_t err = hipGetLastError(); + //memcpy(dest, source, size); + +} + +extern "C" void ScaLBL_CopyToHost(void* dest, const void* source, size_t size){ + hipMemcpy(dest,source,size,hipMemcpyDeviceToHost); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("Error in hipMemcpy (device->host): %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_DeviceBarrier(){ + hipDeviceSynchronize(); +} diff --git a/hip/MRT.hip b/hip/MRT.hip new file mode 100644 index 00000000..671e2801 --- /dev/null +++ b/hip/MRT.hip @@ -0,0 +1,310 @@ +//************************************************************************* +// CUDA kernels for single-phase ScaLBL_D3Q19_MRT code +// James McClure +//************************************************************************* +#include "hip/hip_runtime.h" + +#define NBLOCKS 560 +#define NTHREADS 128 + +__global__ void INITIALIZE(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz) +{ + int n,N; + N = Nx*Ny*Nz; + int S = N/NBLOCKS/NTHREADS + 1; + for (int s=0; s 0){ + f_even[n] = 0.3333333333333333; + f_odd[n] = 0.055555555555555555; //double(100*n)+1.f; + f_even[N+n] = 0.055555555555555555; //double(100*n)+2.f; + f_odd[N+n] = 0.055555555555555555; //double(100*n)+3.f; + f_even[2*N+n] = 0.055555555555555555; //double(100*n)+4.f; + f_odd[2*N+n] = 0.055555555555555555; //double(100*n)+5.f; + f_even[3*N+n] = 0.055555555555555555; //double(100*n)+6.f; + f_odd[3*N+n] = 0.0277777777777778; //double(100*n)+7.f; + f_even[4*N+n] = 0.0277777777777778; //double(100*n)+8.f; + f_odd[4*N+n] = 0.0277777777777778; //double(100*n)+9.f; + f_even[5*N+n] = 0.0277777777777778; //double(100*n)+10.f; + f_odd[5*N+n] = 0.0277777777777778; //double(100*n)+11.f; + f_even[6*N+n] = 0.0277777777777778; //double(100*n)+12.f; + f_odd[6*N+n] = 0.0277777777777778; //double(100*n)+13.f; + f_even[7*N+n] = 0.0277777777777778; //double(100*n)+14.f; + f_odd[7*N+n] = 0.0277777777777778; //double(100*n)+15.f; + f_even[8*N+n] = 0.0277777777777778; //double(100*n)+16.f; + f_odd[8*N+n] = 0.0277777777777778; //double(100*n)+17.f; + f_even[9*N+n] = 0.0277777777777778; //double(100*n)+18.f; + } + else{ + for(int q=0; q<9; q++){ + f_even[q*N+n] = -1.0; + f_odd[q*N+n] = -1.0; + } + f_even[9*N+n] = -1.0; + } + } + } +} + +__global__ void Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz) +{ + int n,N; + // distributions + double f1,f2,f3,f4,f5,f6,f7,f8,f9; + double f10,f11,f12,f13,f14,f15,f16,f17,f18; + double vx,vy,vz; + + N = Nx*Ny*Nz; + int S = N/NBLOCKS/NTHREADS + 1; + // S - number of threadblocks per grid block + for (int s=0; s 0){ + //........................................................................ + // Registers to store the distributions + //........................................................................ + f2 = disteven[N+n]; + f4 = disteven[2*N+n]; + f6 = disteven[3*N+n]; + f8 = disteven[4*N+n]; + f10 = disteven[5*N+n]; + f12 = disteven[6*N+n]; + f14 = disteven[7*N+n]; + f16 = disteven[8*N+n]; + f18 = disteven[9*N+n]; + //........................................................................ + f1 = distodd[n]; + f3 = distodd[1*N+n]; + f5 = distodd[2*N+n]; + f7 = distodd[3*N+n]; + f9 = distodd[4*N+n]; + f11 = distodd[5*N+n]; + f13 = distodd[6*N+n]; + f15 = distodd[7*N+n]; + f17 = distodd[8*N+n]; + //.................Compute the velocity................................... + vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14; + vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18; + vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18; + //..................Write the velocity..................................... + vel[n] = vx; + vel[N+n] = vy; + vel[2*N+n] = vz; + //........................................................................ + + } + } + } +} + +//************************************************************************* +__global__ void +__launch_bounds__(512,2) +D3Q19_MRT(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz, + double rlx_setA, double rlx_setB, double Fx, double Fy, double Fz) +{ + + int n,N; + // distributions + double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9; + double f10,f11,f12,f13,f14,f15,f16,f17,f18; + + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + + N = Nx*Ny*Nz; + + char id; + int S = N/NBLOCKS/NTHREADS + 1; + // S - number of threadblocks per grid block + for (int s=0; s 0){ + //........................................................................ + // Registers to store the distributions - read based on swap convention + //........................................................................ + f2 = distodd[n]; + f4 = distodd[N+n]; + f6 = distodd[2*N+n]; + f8 = distodd[3*N+n]; + f10 = distodd[4*N+n]; + f12 = distodd[5*N+n]; + f14 = distodd[6*N+n]; + f16 = distodd[7*N+n]; + f18 = distodd[8*N+n]; + //........................................................................ + f0 = disteven[n]; + f1 = disteven[N+n]; + f3 = disteven[2*N+n]; + f5 = disteven[3*N+n]; + f7 = disteven[4*N+n]; + f9 = disteven[5*N+n]; + f11 = disteven[6*N+n]; + f13 = disteven[7*N+n]; + f15 = disteven[8*N+n]; + f17 = disteven[9*N+n]; + //........................................................................ + //....................compute the moments............................................... + rho = f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17; + m1 = -30*f0-11*(f2+f1+f4+f3+f6+f5)+8*(f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18 +f17); + m2 = 12*f0-4*(f2+f1 +f4+f3+f6 +f5)+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17; + jx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14; + m4 = 4*(-f1+f2)+f7-f8+f9-f10+f11-f12+f13-f14; + jy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18; + m6 = -4*(f3-f4)+f7-f8-f9+f10+f15-f16+f17-f18; + jz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18; + m8 = -4*(f5-f6)+f11-f12-f13+f14+f15-f16-f17+f18; + m9 = 2*(f1+f2)-f3-f4-f5-f6+f7+f8+f9+f10+f11+f12+f13+f14-2*(f15+f16+f17+f18); + m10 = -4*(f1+f2)+2*(f4+f3+f6+f5)+f8+f7+f10+f9+f12+f11+f14+f13-2*(f16+f15+f18+f17); + m11 = f4+f3-f6-f5+f8+f7+f10+f9-f12-f11-f14-f13; + m12 = -2*(f4+f3-f6-f5)+f8+f7+f10+f9-f12-f11-f14-f13; + m13 = f8+f7-f10-f9; + m14 = f16+f15-f18-f17; + m15 = f12+f11-f14-f13; + m16 = f7-f8+f9-f10-f11+f12-f13+f14; + m17 = -f7+f8+f9-f10+f15-f16+f17-f18; + m18 = f11-f12-f13+f14-f15+f16+f17-f18; + //..............incorporate external force................................................ + //jx += 0.5*Fx; + //jy += 0.5*Fy; + //jz += 0.5*Fz; + //..............carry out relaxation process............................................... + m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) - m1); + m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho) - m2); + m4 = m4 + rlx_setB*((-0.6666666666666666*jx) - m4); + m6 = m6 + rlx_setB*((-0.6666666666666666*jy) - m6); + m8 = m8 + rlx_setB*((-0.6666666666666666*jz) - m8); + m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) - m9); + m10 = m10 + rlx_setA*(-0.5*((2*jx*jx-jy*jy-jz*jz)/rho) - m10); + m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) - m11); + m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho) - m12); + m13 = m13 + rlx_setA*((jx*jy/rho) - m13); + m14 = m14 + rlx_setA*((jy*jz/rho) - m14); + m15 = m15 + rlx_setA*((jx*jz/rho) - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); + //.................inverse transformation...................................................... + f0 = 0.05263157894736842*rho-0.012531328320802*m1+0.04761904761904762*m2; + f1 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(jx-m4)+0.05555555555555555*(m9-m10); + f2 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(m4-jx)+0.05555555555555555*(m9-m10); + f3 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(jy-m6)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m11-m12); + f4 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(m6-jy)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m11-m12); + f5 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(jz-m8)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m12-m11); + f6 = 0.05263157894736842*rho-0.004594820384294068*m1-0.01587301587301587*m2 + +0.1*(m8-jz)+0.02777777777777778*(m10-m9)+0.08333333333333333*(m12-m11); + f7 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2+0.1*(jx+jy)+0.025*(m4+m6) + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12+0.25*m13+0.125*(m16-m17); + f8 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2-0.1*(jx+jy)-0.025*(m4+m6) + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12+0.25*m13+0.125*(m17-m16); + f9 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2+0.1*(jx-jy)+0.025*(m4-m6) + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12-0.25*m13+0.125*(m16+m17); + f10 = 0.05263157894736842*rho+0.003341687552213868*m1+0.003968253968253968*m2+0.1*(jy-jx)+0.025*(m6-m4) + +0.02777777777777778*m9+0.01388888888888889*m10+0.08333333333333333*m11 + +0.04166666666666666*m12-0.25*m13-0.125*(m16+m17); + f11 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jx+jz)+0.025*(m4+m8) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12+0.25*m15+0.125*(m18-m16); + f12 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2-0.1*(jx+jz)-0.025*(m4+m8) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12+0.25*m15+0.125*(m16-m18); + f13 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jx-jz)+0.025*(m4-m8) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12-0.25*m15-0.125*(m16+m18); + f14 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jz-jx)+0.025*(m8-m4) + +0.02777777777777778*m9+0.01388888888888889*m10-0.08333333333333333*m11 + -0.04166666666666666*m12-0.25*m15+0.125*(m16+m18); + f15 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jy+jz)+0.025*(m6+m8) + -0.05555555555555555*m9-0.02777777777777778*m10+0.25*m14+0.125*(m17-m18); + f16 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2-0.1*(jy+jz)-0.025*(m6+m8) + -0.05555555555555555*m9-0.02777777777777778*m10+0.25*m14+0.125*(m18-m17); + f17 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jy-jz)+0.025*(m6-m8) + -0.05555555555555555*m9-0.02777777777777778*m10-0.25*m14+0.125*(m17+m18); + f18 = 0.05263157894736842*rho+0.003341687552213868*m1 + +0.003968253968253968*m2+0.1*(jz-jy)+0.025*(m8-m6) + -0.05555555555555555*m9-0.02777777777777778*m10-0.25*m14-0.125*(m17+m18); + //....................................................................................................... + // incorporate external force + f1 += 0.16666666*Fx; + f2 -= 0.16666666*Fx; + f3 += 0.16666666*Fy; + f4 -= 0.16666666*Fy; + f5 += 0.16666666*Fz; + f6 -= 0.16666666*Fz; + f7 += 0.08333333333*(Fx+Fy); + f8 -= 0.08333333333*(Fx+Fy); + f9 += 0.08333333333*(Fx-Fy); + f10 -= 0.08333333333*(Fx-Fy); + f11 += 0.08333333333*(Fx+Fz); + f12 -= 0.08333333333*(Fx+Fz); + f13 += 0.08333333333*(Fx-Fz); + f14 -= 0.08333333333*(Fx-Fz); + f15 += 0.08333333333*(Fy+Fz); + f16 -= 0.08333333333*(Fy+Fz); + f17 += 0.08333333333*(Fy-Fz); + f18 -= 0.08333333333*(Fy-Fz); + //....................................................................................................... + // Write data based on un-swapped convention + disteven[n] = f0; + disteven[N+n] = f2; + disteven[2*N+n] = f4; + disteven[3*N+n] = f6; + disteven[4*N+n] = f8; + disteven[5*N+n] = f10; + disteven[6*N+n] = f12; + disteven[7*N+n] = f14; + disteven[8*N+n] = f16; + disteven[9*N+n] = f18; + + distodd[n] = f1; + distodd[N+n] = f3; + distodd[2*N+n] = f5; + distodd[3*N+n] = f7; + distodd[4*N+n] = f9; + distodd[5*N+n] = f11; + distodd[6*N+n] = f13; + distodd[7*N+n] = f15; + distodd[8*N+n] = f17; + //....................................................................................................... + } + } + } +} + +extern "C" void ScaLBL_D3Q19_MRT(char *ID, double *f_even, double *f_odd, double rlxA, double rlxB, double Fx, double Fy, double Fz,int Nx, int Ny, int Nz) +{ + D3Q19_MRT <<< NBLOCKS,NTHREADS>>> (ID, f_even, f_odd, Nx, Ny, Nz, rlxA, rlxB, Fx, Fy, Fz); +} + + diff --git a/hip/dfh.hip b/hip/dfh.hip new file mode 100644 index 00000000..37f91498 --- /dev/null +++ b/hip/dfh.hip @@ -0,0 +1,1508 @@ +#include +#include +#include "hip/hip_runtime.h" + +#define NBLOCKS 1024 +#define NTHREADS 256 + +#if !defined(__CUDA_ARCH__) || __CUDA_ARCH__ >= 600 +#else +__device__ double atomicAdd(double* address, double val) { + unsigned long long int* address_as_ull = (unsigned long long int*)address; + unsigned long long int old = *address_as_ull, assumed; + + do { + assumed = old; + old = atomicCAS(address_as_ull, assumed, __double_as_longlong(val+__longlong_as_double(assumed))); + } while (assumed != old); + return __longlong_as_double(old); +} +#endif + +__global__ void dvc_ScaLBL_Gradient_Unpack(double weight, double Cqx, double Cqy, double Cqz, + int *list, int start, int count, double *recvbuf, double *phi, double *grad, int N){ + //.................................................................................... + // Unpack distribution from the recv buffer + // Distribution q matche Cqx, Cqy, Cqz + // swap rule means that the distributions in recvbuf are OPPOSITE of q + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int n,idx; + double value, tmp; + idx = blockIdx.x*blockDim.x + threadIdx.x; + if (idx 0.f){ + nA = 1.0; nB = 0.f; + } + else{ + nB = 1.0; nA = 0.f; + } + Den[idx] = nA; + Den[Np+idx] = nB; + + Aq[idx]=0.3333333333333333*nA; + Aq[Np+idx]=0.1111111111111111*nA; + Aq[2*Np+idx]=0.1111111111111111*nA; + Aq[3*Np+idx]=0.1111111111111111*nA; + Aq[4*Np+idx]=0.1111111111111111*nA; + Aq[5*Np+idx]=0.1111111111111111*nA; + Aq[6*Np+idx]=0.1111111111111111*nA; + + Bq[idx]=0.3333333333333333*nB; + Bq[Np+idx]=0.1111111111111111*nB; + Bq[2*Np+idx]=0.1111111111111111*nB; + Bq[3*Np+idx]=0.1111111111111111*nB; + Bq[4*Np+idx]=0.1111111111111111*nB; + Bq[5*Np+idx]=0.1111111111111111*nB; + Bq[6*Np+idx]=0.1111111111111111*nB; + } + } +} + + +// LBM based on density functional hydrodynamics +__global__ void dvc_ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + double *Gradient, double *SolidForce, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int start, int finish, int Np){ + int nn,n; + double fq; + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + double m3,m5,m7; + double nA,nB; // number density + double a1,b1,a2,b2,nAB,delta; + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + double phi,tau,rho0,rlx_setA,rlx_setB; + double force_x,force_y,force_z; + + const double mrt_V1=0.05263157894736842; + const double mrt_V2=0.012531328320802; + const double mrt_V3=0.04761904761904762; + const double mrt_V4=0.004594820384294068; + const double mrt_V5=0.01587301587301587; + const double mrt_V6=0.0555555555555555555555555; + const double mrt_V7=0.02777777777777778; + const double mrt_V8=0.08333333333333333; + const double mrt_V9=0.003341687552213868; + const double mrt_V10=0.003968253968253968; + const double mrt_V11=0.01388888888888889; + const double mrt_V12=0.04166666666666666; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + Aq[1*Np+n] = a1; + Bq[1*Np+n] = b1; + Aq[2*Np+n] = a2; + Bq[2*Np+n] = b2; + + //............................................... + // q = 2 + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + Aq[3*Np+n] = a1; + Bq[3*Np+n] = b1; + Aq[4*Np+n] = a2; + Bq[4*Np+n] = b2; + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + Aq[5*Np+n] = a1; + Bq[5*Np+n] = b1; + Aq[6*Np+n] = a2; + Bq[6*Np+n] = b2; + //............................................... + } + } +} + + +__global__ void dvc_ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, + double *Phi, double *Gradient, double *SolidForce, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int start, int finish, int Np){ + + int n,nn,nread; + int nr1,nr2,nr3,nr4,nr5,nr6; + int nr7,nr8,nr9,nr10; + int nr11,nr12,nr13,nr14; + //int nr15,nr16,nr17,nr18; + double fq; + // conserved momemnts + double rho,jx,jy,jz; + // non-conserved moments + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + double m3,m5,m7; + double nA,nB; // number density + double a1,b1,a2,b2,nAB,delta; + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + double phi,tau,rho0,rlx_setA,rlx_setB; + double force_x,force_y,force_z; + + const double mrt_V1=0.05263157894736842; + const double mrt_V2=0.012531328320802; + const double mrt_V3=0.04761904761904762; + const double mrt_V4=0.004594820384294068; + const double mrt_V5=0.01587301587301587; + const double mrt_V6=0.0555555555555555555555555; + const double mrt_V7=0.02777777777777778; + const double mrt_V8=0.08333333333333333; + const double mrt_V9=0.003341687552213868; + const double mrt_V10=0.003968253968253968; + const double mrt_V11=0.01388888888888889; + const double mrt_V12=0.04166666666666666; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s even part of dist) + //fq = dist[nread]; // reading the f2 data into register fq + nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) + fq = dist[nr2]; // reading the f2 data into register fq + rho += fq; + m1 -= 11.0*(fq); + m2 -= 4.0*(fq); + jx -= fq; + m4 += 4.0*(fq); + m9 += 2.0*(fq); + m10 -= 4.0*(fq); + + // q=3 + //nread = neighborList[n+2*Np]; // neighbor 4 + //fq = dist[nread]; + nr3 = neighborList[n+2*Np]; // neighbor 4 + fq = dist[nr3]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy = fq; + m6 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 = fq; + m12 = -2.0*fq; + + // q = 4 + //nread = neighborList[n+3*Np]; // neighbor 3 + //fq = dist[nread]; + nr4 = neighborList[n+3*Np]; // neighbor 3 + fq = dist[nr4]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jy -= fq; + m6 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 += fq; + m12 -= 2.0*fq; + + // q=5 + //nread = neighborList[n+4*Np]; + //fq = dist[nread]; + nr5 = neighborList[n+4*Np]; + fq = dist[nr5]; + rho += fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz = fq; + m8 = -4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + + // q = 6 + //nread = neighborList[n+5*Np]; + //fq = dist[nread]; + nr6 = neighborList[n+5*Np]; + fq = dist[nr6]; + rho+= fq; + m1 -= 11.0*fq; + m2 -= 4.0*fq; + jz -= fq; + m8 += 4.0*fq; + m9 -= fq; + m10 += 2.0*fq; + m11 -= fq; + m12 += 2.0*fq; + + // q=7 + //nread = neighborList[n+6*Np]; + //fq = dist[nread]; + nr7 = neighborList[n+6*Np]; + fq = dist[nr7]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 = fq; + m16 = fq; + m17 = -fq; + + // q = 8 + //nread = neighborList[n+7*Np]; + //fq = dist[nread]; + nr8 = neighborList[n+7*Np]; + fq = dist[nr8]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 += fq; + m16 -= fq; + m17 += fq; + + // q=9 + //nread = neighborList[n+8*Np]; + //fq = dist[nread]; + nr9 = neighborList[n+8*Np]; + fq = dist[nr9]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jy -= fq; + m6 -= fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 += fq; + m17 += fq; + + // q = 10 + //nread = neighborList[n+9*Np]; + //fq = dist[nread]; + nr10 = neighborList[n+9*Np]; + fq = dist[nr10]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jy += fq; + m6 += fq; + m9 += fq; + m10 += fq; + m11 += fq; + m12 += fq; + m13 -= fq; + m16 -= fq; + m17 -= fq; + + // q=11 + //nread = neighborList[n+10*Np]; + //fq = dist[nread]; + nr11 = neighborList[n+10*Np]; + fq = dist[nr11]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 = fq; + m16 -= fq; + m18 = fq; + + // q=12 + //nread = neighborList[n+11*Np]; + //fq = dist[nread]; + nr12 = neighborList[n+11*Np]; + fq = dist[nr12]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 += fq; + m16 += fq; + m18 -= fq; + + // q=13 + //nread = neighborList[n+12*Np]; + //fq = dist[nread]; + nr13 = neighborList[n+12*Np]; + fq = dist[nr13]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx += fq; + m4 += fq; + jz -= fq; + m8 -= fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 -= fq; + m18 -= fq; + + // q=14 + //nread = neighborList[n+13*Np]; + //fq = dist[nread]; + nr14 = neighborList[n+13*Np]; + fq = dist[nr14]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jx -= fq; + m4 -= fq; + jz += fq; + m8 += fq; + m9 += fq; + m10 += fq; + m11 -= fq; + m12 -= fq; + m15 -= fq; + m16 += fq; + m18 += fq; + + // q=15 + nread = neighborList[n+14*Np]; + fq = dist[nread]; + //fq = dist[17*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 = fq; + m17 += fq; + m18 -= fq; + + // q=16 + nread = neighborList[n+15*Np]; + fq = dist[nread]; + //fq = dist[8*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 += fq; + m17 -= fq; + m18 += fq; + + // q=17 + //fq = dist[18*Np+n]; + nread = neighborList[n+16*Np]; + fq = dist[nread]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy += fq; + m6 += fq; + jz -= fq; + m8 -= fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 += fq; + m18 += fq; + + // q=18 + nread = neighborList[n+17*Np]; + fq = dist[nread]; + //fq = dist[9*Np+n]; + rho += fq; + m1 += 8.0*fq; + m2 += fq; + jy -= fq; + m6 -= fq; + jz += fq; + m8 += fq; + m9 -= 2.0*fq; + m10 -= 2.0*fq; + m14 -= fq; + m17 -= fq; + m18 -= fq; + + //........................................................................ + //..............carry out relaxation process.............................. + //..........Toelke, Fruediger et. al. 2006................................ + if (C == 0.0) nx = ny = nz = 0.0; + m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -alpha*C - m1); + m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); + m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); + m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6); + m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8); + m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9); + m10 = m10 + rlx_setA*( - m10); + m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) + 0.5*alpha*C*(ny*ny-nz*nz)- m11); + m12 = m12 + rlx_setA*( - m12); + m13 = m13 + rlx_setA*( (jx*jy/rho0) + 0.5*alpha*C*nx*ny - m13); + m14 = m14 + rlx_setA*( (jy*jz/rho0) + 0.5*alpha*C*ny*nz - m14); + m15 = m15 + rlx_setA*( (jx*jz/rho0) + 0.5*alpha*C*nx*nz - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); + + // assign force with wetting BC + force_x = alpha*(nA-nB)*SolidForce[n] + Fx; + force_y = alpha*(nA-nB)*SolidForce[n+Np] + Fy; + force_z = alpha*(nA-nB)*SolidForce[n+2*Np] + Fz; + + //.................inverse transformation...................................................... + // q=0 + fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; + dist[n] = fq; + + // q = 1 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10)+0.16666666*force_x; + //nread = neighborList[n+Np]; + dist[nr2] = fq; + + // q=2 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10) - 0.16666666*force_x; + //nread = neighborList[n]; + dist[nr1] = fq; + + // q = 3 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) + 0.16666666*force_y; + //nread = neighborList[n+3*Np]; + dist[nr4] = fq; + + // q = 4 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) - 0.16666666*force_y; + //nread = neighborList[n+2*Np]; + dist[nr3] = fq; + + // q = 5 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) + 0.16666666*force_z; + //nread = neighborList[n+5*Np]; + dist[nr6] = fq; + + // q = 6 + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) - 0.16666666*force_z; + //nread = neighborList[n+4*Np]; + dist[nr5] = fq; + + // q = 7 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17) + 0.08333333333*(force_x+force_y); + //nread = neighborList[n+7*Np]; + dist[nr8] = fq; + + // q = 8 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11 + +mrt_V12*m12+0.25*m13+0.125*(m17-m16) - 0.08333333333*(force_x+force_y); + //nread = neighborList[n+6*Np]; + dist[nr7] = fq; + + // q = 9 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17) + 0.08333333333*(force_x-force_y); + //nread = neighborList[n+9*Np]; + dist[nr10] = fq; + + // q = 10 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+ + mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17)- 0.08333333333*(force_x-force_y); + //nread = neighborList[n+8*Np]; + dist[nr9] = fq; + + // q = 11 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12+0.25*m15+0.125*(m18-m16) + 0.08333333333*(force_x+force_z); + //nread = neighborList[n+11*Np]; + dist[nr12] = fq; + + // q = 12 + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+ + mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18) - 0.08333333333*(force_x+force_z); + //nread = neighborList[n+10*Np]; + dist[nr11]= fq; + + // q = 13 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15-0.125*(m16+m18) + 0.08333333333*(force_x-force_z); + //nread = neighborList[n+13*Np]; + dist[nr14] = fq; + + // q= 14 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4) + +mrt_V7*m9+mrt_V11*m10-mrt_V8*m11 + -mrt_V12*m12-0.25*m15+0.125*(m16+m18) - 0.08333333333*(force_x-force_z); + //nread = neighborList[n+12*Np]; + dist[nr13] = fq; + + + // q = 15 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18) + 0.08333333333*(force_y+force_z); + nread = neighborList[n+15*Np]; + dist[nread] = fq; + + // q = 16 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8) + -mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17)- 0.08333333333*(force_y+force_z); + nread = neighborList[n+14*Np]; + dist[nread] = fq; + + + // q = 17 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8) + -mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18) + 0.08333333333*(force_y-force_z); + nread = neighborList[n+17*Np]; + dist[nread] = fq; + + // q = 18 + fq = mrt_V1*rho+mrt_V9*m1 + +mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6) + -mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18) - 0.08333333333*(force_y-force_z); + nread = neighborList[n+16*Np]; + dist[nread] = fq; + + // write the velocity + ux = (jx + force_x) / rho0; + uy = (jy + force_y) / rho0; + uz = (jz + force_z) / rho0; + //Velocity[n] = ux; + //Velocity[Np+n] = uy; + //Velocity[2*Np+n] = uz; + + // Instantiate mass transport distributions + // Stationary value - distribution 0 + nAB = 1.0/(nA+nB); + Aq[n] = 0.3333333333333333*nA; + Bq[n] = 0.3333333333333333*nB; + + //............................................... + // q = 0,2,4 + // Cq = {1,0,0}, {0,1,0}, {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nx; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + // q = 1 + //nread = neighborList[n+Np]; + Aq[nr2] = a1; + Bq[nr2] = b1; + // q=2 + //nread = neighborList[n]; + Aq[nr1] = a2; + Bq[nr1] = b2; + + //............................................... + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + // q = 3 + //nread = neighborList[n+3*Np]; + Aq[nr4] = a1; + Bq[nr4] = b1; + // q = 4 + //nread = neighborList[n+2*Np]; + Aq[nr3] = a2; + Bq[nr3] = b2; + + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + // q = 5 + //nread = neighborList[n+5*Np]; + Aq[nr6] = a1; + Bq[nr6] = b1; + // q = 6 + //nread = neighborList[n+4*Np]; + Aq[nr5] = a2; + Bq[nr5] = b2; + //............................................... + } + } +} + +__global__ void dvc_ScaLBL_D3Q7_AAodd_DFH(int *neighborList, double *Aq, double *Bq, + double *Den, double *Phi, int start, int finish, int Np){ + int n,nread; + double fq,nA,nB; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s>>(weight, Cqx, Cqy, Cqz, list, start, count, recvbuf, phi, grad, N); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_Gradient_Unpack: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np){ + dvc_ScaLBL_DFH_Init<<>>(Phi, Den, Aq, Bq, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_DFH_Init: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + double *Gradient, double *SolidForce, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_DFH, hipFuncCachePreferL1); + + dvc_ScaLBL_D3Q19_AAeven_DFH<<>>(neighborList, dist, Aq, Bq, Den, Phi, Gradient, SolidForce, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_DFH: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, + double *Phi, double *Gradient, double *SolidForce, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double Fx, double Fy, double Fz, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_DFH, hipFuncCachePreferL1); + + dvc_ScaLBL_D3Q19_AAodd_DFH<<>>(neighborList,dist, Aq, Bq, Den, Phi, Gradient, + SolidForce, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, start, finish, Np); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_DFH: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q7_AAodd_DFH(int *NeighborList, double *Aq, double *Bq, + double *Den, double *Phi, int start, int finish, int Np){ + + dvc_ScaLBL_D3Q7_AAodd_DFH<<>>(NeighborList, Aq, Bq, Den, Phi, start, finish, Np); + + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_AAodd_DFH: %s \n",hipGetErrorString(err)); + } +} + +extern "C" void ScaLBL_D3Q7_AAeven_DFH(double *Aq, double *Bq, double *Den, double *Phi, + int start, int finish, int Np){ + + dvc_ScaLBL_D3Q7_AAeven_DFH<<>>(Aq, Bq, Den, Phi, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_AAeven_DFH: %s \n",hipGetErrorString(err)); + } + +} + +extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *neighborList, double *Phi, double *ColorGrad, int start, int finish, int Np){ + + dvc_ScaLBL_D3Q19_Gradient_DFH<<>>(neighborList, Phi, ColorGrad, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Gradient_DFH: %s \n",hipGetErrorString(err)); + } + +}