diff --git a/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice.cu index c869c9f6..ee48fcdf 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice.cu @@ -161,6 +161,87 @@ __device__ void cudaDevice_UpstreamDivAdvFlux(float* scalarField, float* scalarF } //end cudaDevice_UpstreamDivAdvFlux( +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxX(); -------------------------------------------------- +* This is the cuda version of the UpstreamDivAdvFluxX routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_UpstreamDivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float* invD_Jac_d){ + int i,j,k; + int ijk,im1jk,ip1jk; + int iStride,jStride,kStride; + float DscalarDx; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + im1jk = (i-1)*iStride + j*jStride + k*kStride; + ip1jk = (i+1)*iStride + j*jStride + k*kStride; + DscalarDx = ( ( fmaxf(0.0,u_cf[ ip1jk ])*scalarField[ ijk ]+fminf(0.0,u_cf[ ip1jk ])*scalarField[ ip1jk ]) + -( fmaxf(0.0,u_cf[ ijk ])*scalarField[ im1jk ]+fminf(0.0,u_cf[ ijk ])*scalarField[ ijk ]) ); + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDx; + +} //end cudaDevice_UpstreamDivAdvFluxX( + +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxY(); -------------------------------------------------- +* This is the cuda version of the UpstreamDivAdvFluxY routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_UpstreamDivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float* invD_Jac_d){ + int i,j,k; + int ijk,ijm1k,ijp1k; + int iStride,jStride,kStride; + float DscalarDy; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijm1k = i*iStride + (j-1)*jStride + k*kStride; + ijp1k = i*iStride + (j+1)*jStride + k*kStride; + DscalarDy = ( ( fmaxf(0.0,v_cf[ ijp1k ])*scalarField[ ijk ]+fminf(0.0,v_cf[ ijp1k ])*scalarField[ ijp1k ]) + -( fmaxf(0.0,v_cf[ ijk ])*scalarField[ ijm1k ]+fminf(0.0,v_cf[ ijk ])*scalarField[ ijk ]) ); + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDy; + +} //end cudaDevice_UpstreamDivAdvFluxY( + +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxZ(); -------------------------------------------------- +* This is the cuda version of the UpstreamDivAdvFluxZ routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_UpstreamDivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float* invD_Jac_d){ + int i,j,k; + int ijk,ijkm1,ijkp1; + int iStride,jStride,kStride; + float DscalarDz; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijkm1 = i*iStride + j*jStride + (k-1)*kStride; + ijkp1 = i*iStride + j*jStride + (k+1)*kStride; + DscalarDz = ( ( fmaxf(0.0,w_cf[ ijkp1 ])*scalarField[ ijk ]+fminf(0.0,w_cf[ ijkp1 ])*scalarField[ ijkp1 ]) + -( fmaxf(0.0,w_cf[ ijk ])*scalarField[ ijkm1 ]+fminf(0.0,w_cf[ ijk ])*scalarField[ ijk ]) ); + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDz; + +} //end cudaDevice_UpstreamDivAdvFluxZ( + /*----->>>>> __device__ void cudaDevice_SecondDivAdvFlux(); -----------------------------------------------*/ __device__ void cudaDevice_SecondDivAdvFlux(float* scalarField, float* scalarFadv, float* u_cf, float* v_cf, float* w_cf, float* invD_Jac_d){ @@ -315,6 +396,117 @@ __device__ void cudaDevice_HYB34DivAdvFlux(float* scalarField, float* scalarFadv } //end cudaDevice_HYB34DivAdvFlux() +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxX(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB34DivAdvFluxX routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB34DivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,im1jk,ip1jk,im2jk,ip2jk; + int iStride,jStride,kStride; + float DscalarDx; + float one_twelfth; + float flxx_ipf,flxx_imf; + + one_twelfth = 1.0/12.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + im1jk = (i-1)*iStride + j*jStride + k*kStride; + ip1jk = (i+1)*iStride + j*jStride + k*kStride; + im2jk = (i-2)*iStride + j*jStride + k*kStride; + ip2jk = (i+2)*iStride + j*jStride + k*kStride; + + flxx_ipf = one_twelfth * ( 7.0*(scalarField[ ip1jk ]+scalarField[ ijk ])-(scalarField[ ip2jk ]+scalarField[ im1jk ])+ + (1.0-b_hyb_p)*copysign(1.0,u_cf[ ip1jk ])*((scalarField[ ip2jk ]-scalarField[ im1jk ])-3.0*(scalarField[ ip1jk ]-scalarField[ ijk ])) ); + flxx_imf = one_twelfth * ( 7.0*(scalarField[ ijk ]+scalarField[ im1jk ])-(scalarField[ ip1jk ]+scalarField[ im2jk ])+ + (1.0-b_hyb_p)*copysign(1.0,u_cf[ ijk ])*((scalarField[ ip1jk ]-scalarField[ im2jk ])-3.0*(scalarField[ ijk ]-scalarField[ im1jk ])) ); + DscalarDx = u_cf[ ip1jk ]*flxx_ipf - u_cf[ ijk ]*flxx_imf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDx; + +} //end cudaDevice_HYB34DivAdvFluxX() + +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxY(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB34DivAdvFluxY routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB34DivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,ijm1k,ijp1k,ijm2k,ijp2k; + int iStride,jStride,kStride; + float DscalarDy; + float one_twelfth; + float flxy_jpf,flxy_jmf; + + one_twelfth = 1.0/12.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijm1k = i*iStride + (j-1)*jStride + k*kStride; + ijp1k = i*iStride + (j+1)*jStride + k*kStride; + ijm2k = i*iStride + (j-2)*jStride + k*kStride; + ijp2k = i*iStride + (j+2)*jStride + k*kStride; + + flxy_jpf = one_twelfth * ( 7.0*(scalarField[ ijp1k ]+scalarField[ ijk ])-(scalarField[ ijp2k ]+scalarField[ ijm1k ])+ + (1.0-b_hyb_p)*copysign(1.0,v_cf[ ijp1k ])*((scalarField[ ijp2k ]-scalarField[ ijm1k ])-3.0*(scalarField[ ijp1k ]-scalarField[ ijk ])) ); + flxy_jmf = one_twelfth * ( 7.0*(scalarField[ ijk ]+scalarField[ ijm1k ])-(scalarField[ ijp1k ]+scalarField[ ijm2k ])+ + (1.0-b_hyb_p)*copysign(1.0,v_cf[ ijk ])*((scalarField[ ijp1k ]-scalarField[ ijm2k ])-3.0*(scalarField[ ijk ]-scalarField[ ijm1k ])) ); + DscalarDy = v_cf[ ijp1k ]*flxy_jpf - v_cf[ ijk ]*flxy_jmf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDy; + +} //end cudaDevice_HYB34DivAdvFluxY() + +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxZ(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB34DivAdvFluxZ routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB34DivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,ijkm1,ijkp1,ijkm2,ijkp2; + int iStride,jStride,kStride; + float DscalarDz; + float one_twelfth; + float flxz_kpf,flxz_kmf; + + one_twelfth = 1.0/12.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijkm1 = i*iStride + j*jStride + (k-1)*kStride; + ijkp1 = i*iStride + j*jStride + (k+1)*kStride; + ijkm2 = i*iStride + j*jStride + (k-2)*kStride; + ijkp2 = i*iStride + j*jStride + (k+2)*kStride; + + flxz_kpf = one_twelfth * ( 7.0*(scalarField[ ijkp1 ]+scalarField[ ijk ])-(scalarField[ ijkp2 ]+scalarField[ ijkm1 ])+ + (1.0-b_hyb_p)*copysign(1.0,w_cf[ ijkp1 ])*((scalarField[ ijkp2 ]-scalarField[ ijkm1 ])-3.0*(scalarField[ ijkp1 ]-scalarField[ ijk ])) ); + flxz_kmf = one_twelfth * ( 7.0*(scalarField[ ijk ]+scalarField[ ijkm1 ])-(scalarField[ ijkp1 ]+scalarField[ ijkm2 ])+ + (1.0-b_hyb_p)*copysign(1.0,w_cf[ ijk ])*((scalarField[ ijkp1 ]-scalarField[ ijkm2 ])-3.0*(scalarField[ ijk ]-scalarField[ ijkm1 ])) ); + DscalarDz = w_cf[ ijkp1 ]*flxz_kpf - w_cf[ ijk ]*flxz_kmf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDz; + +} //end cudaDevice_HYB34DivAdvFluxZ() + /*----->>>>> __device__ void cudaDevice_HYB56DivAdvFlux(); -------------------------------------------------- * This is the cuda version of the hydro_coreHYB56DivAdvFlx routine from the HYDRO_CORE module */ @@ -382,6 +574,132 @@ __device__ void cudaDevice_HYB56DivAdvFlux(float* scalarField, float* scalarFadv } //end cudaDevice_HYB56DivAdvFlux( +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxX(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB56DivAdvFluxX routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB56DivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,im1jk,ip1jk; + int im2jk,ip2jk; + int im3jk,ip3jk; + int iStride,jStride,kStride; + float DscalarDx; + float one_sixtieth; + float flxx_ipf,flxx_imf; + + one_sixtieth = 1.0/60.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + im1jk = (i-1)*iStride + j*jStride + k*kStride; + ip1jk = (i+1)*iStride + j*jStride + k*kStride; + im2jk = (i-2)*iStride + j*jStride + k*kStride; + ip2jk = (i+2)*iStride + j*jStride + k*kStride; + im3jk = (i-3)*iStride + j*jStride + k*kStride; + ip3jk = (i+3)*iStride + j*jStride + k*kStride; + + flxx_ipf = one_sixtieth * ( 37.0*(scalarField[ ip1jk ]+scalarField[ ijk ])-8.0*(scalarField[ ip2jk ]+scalarField[ im1jk ])+(scalarField[ ip3jk ]+scalarField[ im2jk ])- + (1.0-b_hyb_p)*copysign(1.0,u_cf[ ip1jk ])*((scalarField[ ip3jk ]-scalarField[ im2jk ])-5.0*(scalarField[ ip2jk ]-scalarField[ im1jk ])+10.0*(scalarField[ ip1jk ]-scalarField[ ijk ])) ); + flxx_imf = one_sixtieth * ( 37.0*(scalarField[ ijk ]+scalarField[ im1jk ])-8.0*(scalarField[ ip1jk ]+scalarField[ im2jk ])+(scalarField[ ip2jk ]+scalarField[ im3jk ])- + (1.0-b_hyb_p)*copysign(1.0,u_cf[ ijk ])*((scalarField[ ip2jk ]-scalarField[ im3jk ])-5.0*(scalarField[ ip1jk ]-scalarField[ im2jk ])+10.0*(scalarField[ ijk ]-scalarField[ im1jk ])) ); + + DscalarDx = u_cf[ ip1jk ]*flxx_ipf - u_cf[ ijk ]*flxx_imf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDx; + +} //end cudaDevice_HYB56DivAdvFluxX( + +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxY(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB56DivAdvFluxY routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB56DivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,ijm1k,ijp1k; + int ijm2k,ijp2k; + int ijm3k,ijp3k; + int iStride,jStride,kStride; + float DscalarDy; + float one_sixtieth; + float flxy_jpf,flxy_jmf; + + one_sixtieth = 1.0/60.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijm1k = i*iStride + (j-1)*jStride + k*kStride; + ijp1k = i*iStride + (j+1)*jStride + k*kStride; + ijm2k = i*iStride + (j-2)*jStride + k*kStride; + ijp2k = i*iStride + (j+2)*jStride + k*kStride; + ijm3k = i*iStride + (j-3)*jStride + k*kStride; + ijp3k = i*iStride + (j+3)*jStride + k*kStride; + + flxy_jpf = one_sixtieth * ( 37.0*(scalarField[ ijp1k ]+scalarField[ ijk ])-8.0*(scalarField[ ijp2k ]+scalarField[ ijm1k ])+(scalarField[ ijp3k ]+scalarField[ ijm2k ])- + (1.0-b_hyb_p)*copysign(1.0,v_cf[ ijp1k ])*((scalarField[ ijp3k ]-scalarField[ ijm2k ])-5.0*(scalarField[ ijp2k ]-scalarField[ ijm1k ])+10.0*(scalarField[ ijp1k ]-scalarField[ ijk ])) ); + flxy_jmf = one_sixtieth * ( 37.0*(scalarField[ ijk ]+scalarField[ ijm1k ])-8.0*(scalarField[ ijp1k ]+scalarField[ ijm2k ])+(scalarField[ ijp2k ]+scalarField[ ijm3k ])- + (1.0-b_hyb_p)*copysign(1.0,v_cf[ ijk ])*((scalarField[ ijp2k ]-scalarField[ ijm3k ])-5.0*(scalarField[ ijp1k ]-scalarField[ ijm2k ])+10.0*(scalarField[ ijk ]-scalarField[ ijm1k ])) ); + + DscalarDy = v_cf[ ijp1k ]*flxy_jpf - v_cf[ ijk ]*flxy_jmf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDy; + +} //end cudaDevice_HYB56DivAdvFluxY( + +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxZ(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB56DivAdvFluxZ routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB56DivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,ijkm1,ijkp1; + int ijkm2,ijkp2; + int ijkm3,ijkp3; + int iStride,jStride,kStride; + float DscalarDz; + float one_sixtieth; + float flxz_kpf,flxz_kmf; + + one_sixtieth = 1.0/60.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijkm1 = i*iStride + j*jStride + (k-1)*kStride; + ijkp1 = i*iStride + j*jStride + (k+1)*kStride; + ijkm2 = i*iStride + j*jStride + (k-2)*kStride; + ijkp2 = i*iStride + j*jStride + (k+2)*kStride; + ijkm3 = i*iStride + j*jStride + (k-3)*kStride; + ijkp3 = i*iStride + j*jStride + (k+3)*kStride; + + flxz_kpf = one_sixtieth * ( 37.0*(scalarField[ ijkp1 ]+scalarField[ ijk ])-8.0*(scalarField[ ijkp2 ]+scalarField[ ijkm1 ])+(scalarField[ ijkp3 ]+scalarField[ ijkm2 ])- + (1.0-b_hyb_p)*copysign(1.0,w_cf[ ijkp1 ])*((scalarField[ ijkp3 ]-scalarField[ ijkm2 ])-5.0*(scalarField[ ijkp2 ]-scalarField[ ijkm1 ])+10.0*(scalarField[ ijkp1 ]-scalarField[ ijk ])) ); + flxz_kmf = one_sixtieth * ( 37.0*(scalarField[ ijk ]+scalarField[ ijkm1 ])-8.0*(scalarField[ ijkp1 ]+scalarField[ ijkm2 ])+(scalarField[ ijkp2 ]+scalarField[ ijkm3 ])- + (1.0-b_hyb_p)*copysign(1.0,w_cf[ ijk ])*((scalarField[ ijkp2 ]-scalarField[ ijkm3 ])-5.0*(scalarField[ ijkp1 ]-scalarField[ ijkm2 ])+10.0*(scalarField[ ijk ]-scalarField[ ijkm1 ])) ); + + DscalarDz = w_cf[ ijkp1 ]*flxz_kpf - w_cf[ ijk ]*flxz_kmf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDz; + +} //end cudaDevice_HYB56DivAdvFluxZ( + /*----->>>>> __device__ void cudaDevice_WENO3DivAdvFluxX(); -----------------------------------------------*/ __device__ void cudaDevice_WENO3DivAdvFluxX(float* scalarField, float* scalarFadv,float* u_cf, float* invD_Jac_d){ diff --git a/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice_cu.h b/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice_cu.h index f9234949..359ab3a9 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice_cu.h +++ b/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice_cu.h @@ -50,6 +50,15 @@ __device__ void cudaDevice_calcFaceVelocities(float* hydroFlds_d, float* hydroFa __device__ void cudaDevice_UpstreamDivAdvFlux(float* scalarField, float* scalarFadv, float* u_cf, float* v_cf, float* w_cf, float* invD_Jac_d); +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxX(); -------------------------------------------------- */ +__device__ void cudaDevice_UpstreamDivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxY(); -------------------------------------------------- */ +__device__ void cudaDevice_UpstreamDivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxZ(); -------------------------------------------------- */ +__device__ void cudaDevice_UpstreamDivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float* invD_Jac_d); + /*----->>>>> __device__ void cudaDevice_SecondDivAdvFlux(); -------------------------------------------------- */ __device__ void cudaDevice_SecondDivAdvFlux(float* scalarField, float* scalarFadv, @@ -67,12 +76,30 @@ __device__ void cudaDevice_QUICKDivAdvFlux(float* scalarField, float* scalarFadv __device__ void cudaDevice_HYB34DivAdvFlux(float* scalarField, float* scalarFadv, float* u_cf, float* v_cf, float* w_cf, float b_hyb_p, float* invD_Jac_d); +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxX(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB34DivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float b_hyb_p, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxY(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB34DivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float b_hyb_p, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxZ(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB34DivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float b_hyb_p, float* invD_Jac_d); + /*----->>>>> __device__ void cudaDevice_HYB56DivAdvFlux(); -------------------------------------------------- * This is the cuda version of the HYB56DivAdvFlux routine from the HYDRO_CORE module */ __device__ void cudaDevice_HYB56DivAdvFlux(float* scalarField, float* scalarFadv, float* u_cf, float* v_cf, float* w_cf, float b_hyb_p, float* invD_Jac_d); +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxX(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB56DivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float b_hyb_p, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxY(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB56DivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float b_hyb_p, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxZ(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB56DivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float b_hyb_p, float* invD_Jac_d); + /*----->>>>> __device__ void cudaDevice_WENO3DivAdvFluxX(); -------------------------------------------------- */ __device__ void cudaDevice_WENO3DivAdvFluxX(float* scalarField, float* scalarFadv,float* u_cf, float* invD_Jac_d); diff --git a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu index 5c691ad1..d06a37a1 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu @@ -760,7 +760,7 @@ __global__ void cudaDevice_hydroCoreCommence(int simTime_it, float* hydroFlds_d, fld = &sgstkeScalars_d[fldStride*iFld]; fldBS = &sgstkeScalarsFrhs_d[fldStride*iFld]; // Frhs forcing iwas set to zero, so it can be used here as zero-valued base state if(hydroBCs_d == 1){ //Using LAD BCs - cudaDevice_VerticalAblBCs(iFld, fld, fldBS); + cudaDevice_VerticalAblZeroGradBCs(fld); if(rankXid_d == 0){ cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 0); } @@ -774,7 +774,7 @@ __global__ void cudaDevice_hydroCoreCommence(int simTime_it, float* hydroFlds_d, cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 3); } }else if (hydroBCs_d == 2){ - cudaDevice_VerticalAblBCs(iFld, fld, fldBS); // to apply zero-gradient lower boundary BCs + cudaDevice_VerticalAblZeroGradBCs(fld); if(numProcsX_d==1){ cudaDevice_HorizontalPeriodicXdirBCs(iFld, fld); }//periodic and single rank in X-dir --> implies no MPI exchanges made so perform on-device exchange @@ -902,7 +902,7 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa if((i >= iMin_d)&&(i < iMax_d) && (j >= jMin_d)&&(j < jMax_d) && - (k >= kMin_d)&&(k < kMax_d) ){ + (k >= kMin_d+3)&&(k < kMax_d) ){ // skipping the first 3 vertical levels for(iFld=0; iFld < Nhydro_d; iFld++){ fld = &hydroFlds[fldStride*iFld]; fldFrhs = &hydroFldsFrhs[fldStride*iFld]; @@ -926,24 +926,6 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa } else { // defaults to 1st-order upwinding cudaDevice_UpstreamDivAdvFlux(fld, fldFrhs, u_cf, v_cf, w_cf, invD_Jac_d); } - if(iFld==W_INDX){ - if(dampingLayerSelector_d > 0){ // RAYLEIGH DAMPING ON W ******!!!!!!!! - cudaDevice_topRayleighDampingLayerForcing(fld, fldFrhs, - &rho[0], &rho_BS[0], zPos_d); - } //end if dampingLayerSelector > 0 - if(buoyancySelector_d > 0){ // BUOYANCY SOURCE?SINK OF W ******!!!!!!!! - ijk = i*iStride + j*jStride + k*kStride; - if (moistureSelector_d>0){ - if(moistureNvars_d==1){ - cudaDevice_calcBuoyancyMoistNvar1(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk],&moistScalars[ijk]); - }else if(moistureNvars_d==2){ - cudaDevice_calcBuoyancyMoistNvar2(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk],&moistScalars[ijk],&moistScalars[fldStride+ijk]); - } - }else{ - cudaDevice_calcBuoyancy(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk]); - } - } //end if buoyancySelector > 0 - }//end if iFld==W_INDX }//for iFld if ((turbulenceSelector_d>0) && (TKESelector_d>0)){ // : advection of SGSTKE fields for(iFld=0; iFld < TKESelector_d; iFld++){ @@ -1012,7 +994,93 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa } } } + }//end if in the range of non-halo cells (skipping the first 3 vertical levels) + // lower the order of adection as the surface is approached + if((i >= iMin_d)&&(i < iMax_d) && + (j >= jMin_d)&&(j < jMax_d) && + (k >= kMin_d)&&(k < kMin_d+3) ){ // only first 3 vertical grid levels + for(iFld=0; iFld < Nhydro_d; iFld++){ + fld = &hydroFlds[fldStride*iFld]; + fldFrhs = &hydroFldsFrhs[fldStride*iFld]; + /* Calculate scalar, cell-valued divergence of the advective flux */ + cudaDevice_HYB34DivAdvFluxX(fld, fldFrhs, u_cf, b_hyb_d, invD_Jac_d); + cudaDevice_HYB34DivAdvFluxY(fld, fldFrhs, v_cf, b_hyb_d, invD_Jac_d); + if (k == kMin_d+2) { // hybrid 3rd-4th order + cudaDevice_HYB34DivAdvFluxZ(fld, fldFrhs, w_cf, b_hyb_d, invD_Jac_d); + } else { // 1st-order upwinding + cudaDevice_UpstreamDivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } + }//for iFld + if ((turbulenceSelector_d>0) && (TKESelector_d>0)){ // : advection of SGSTKE fields + for(iFld=0; iFld < TKESelector_d; iFld++){ + fld = &sgstkeScalars[fldStride*iFld]; + fldFrhs = &sgstkeScalarsFrhs[fldStride*iFld]; + TKEAdvSelector_flag = TKEAdvSelector_d; + TKEAdvSelector_b_hyb_flag = TKEAdvSelector_b_hyb_d; + /* Calculate scalar, cell-valued divergence of the advective flux */ + cudaDevice_HYB34DivAdvFluxX(fld, fldFrhs, u_cf, TKEAdvSelector_b_hyb_flag, invD_Jac_d); + cudaDevice_HYB34DivAdvFluxY(fld, fldFrhs, v_cf, TKEAdvSelector_b_hyb_flag, invD_Jac_d); + if (k == kMin_d+2) { // hybrid 3rd-4th order + cudaDevice_HYB34DivAdvFluxZ(fld, fldFrhs, w_cf, TKEAdvSelector_b_hyb_flag, invD_Jac_d); + } else { // 1st-order upwinding + cudaDevice_UpstreamDivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } + } + } + + if ((moistureSelector_d>0) && (moistureNvars_d>0)){ // : advection of moisture fields + for(iFld=0; iFld < moistureNvars_d; iFld++){ + fld = &moistScalars[fldStride*iFld]; + fldFrhs = &moistScalarsFrhs[fldStride*iFld]; + if (iFld==0){ // water vapor + cudaDevice_HYB34DivAdvFluxX(fld, fldFrhs, u_cf, moistureAdvSelectorQv_b_d, invD_Jac_d); + cudaDevice_HYB34DivAdvFluxY(fld, fldFrhs, v_cf, moistureAdvSelectorQv_b_d, invD_Jac_d); + if (k == kMin_d+2) { // hybrid 3rd-4th order + cudaDevice_HYB34DivAdvFluxZ(fld, fldFrhs, w_cf, moistureAdvSelectorQv_b_d, invD_Jac_d); + } else { // 1st-order upwinding + cudaDevice_UpstreamDivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } + } else { // non-qv moisture species (non-oscillatory schemes) + if (moistureAdvSelectorQi_d == 0) { // 1st-order upstream + cudaDevice_UpstreamDivAdvFlux(fld, fldFrhs, u_cf, v_cf, w_cf, invD_Jac_d); + } else { + cudaDevice_WENO3DivAdvFluxX(fld, fldFrhs, u_cf, invD_Jac_d); + cudaDevice_WENO3DivAdvFluxY(fld, fldFrhs, v_cf, invD_Jac_d); + if (k == kMin_d+2) { // 3rd-order WENO + cudaDevice_WENO3DivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } else { // 1st-order upstream + cudaDevice_UpstreamDivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } + } + } + } + } + }//end if in the range of non-halo cells (only first 3 vertical grid levels) + + if((i >= iMin_d)&&(i < iMax_d) && + (j >= jMin_d)&&(j < jMax_d) && + (k >= kMin_d)&&(k < kMax_d) ){ + // W terms + iFld=W_INDX; + fld = &hydroFlds[fldStride*iFld]; + fldFrhs = &hydroFldsFrhs[fldStride*iFld]; + if(dampingLayerSelector_d > 0){ // RAYLEIGH DAMPING ON W ******!!!!!!!! + cudaDevice_topRayleighDampingLayerForcing(fld, fldFrhs, + &rho[0], &rho_BS[0], zPos_d); + } //end if dampingLayerSelector > 0 + if(buoyancySelector_d > 0){ // BUOYANCY SOURCE?SINK OF W ******!!!!!!!! + ijk = i*iStride + j*jStride + k*kStride; + if (moistureSelector_d>0){ + if(moistureNvars_d==1){ + cudaDevice_calcBuoyancyMoistNvar1(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk],&moistScalars[ijk]); + }else if(moistureNvars_d==2){ + cudaDevice_calcBuoyancyMoistNvar2(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk],&moistScalars[ijk],&moistScalars[fldStride+ijk]); + } + }else{ + cudaDevice_calcBuoyancy(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk]); + } + } //end if buoyancySelector > 0 if(coriolisSelector_d > 0){ ijk = i*iStride + j*jStride + k*kStride; cudaDevice_MomentumBS(U_INDX, zPos_d[ijk], &hydroBaseStateFlds[RHO_INDX_BS*fldStride+ijk], &MomBSval[0]); @@ -1031,6 +1099,7 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa &MomBSval[2]); } //end if coriolisSelector_d > 0 }//end if in the range of non-halo cells + if((turbulenceSelector_d > 0) && ((physics_oneRKonly_d==0) || (timeStage==numRKstages))){ cudaDevice_hydroCoreCalcTurbMixing( &hydroFldsFrhs[fldStride*U_INDX],