Interp1PrimFifthOrderWENO_GPU.cu

Go to the documentation of this file.

#include <basic.h>
#include <basic_gpu.h>
#include <arrayfunctions_gpu.h>
#include <mathfunctions.h>
#include <interpolation.h>
#include <mpivars.h>
#include <hypar.h>

#undef  _MINIMUM_GHOSTS_

#define _MINIMUM_GHOSTS_ 3

#ifdef CUDA_VAR_ORDERDING_AOS

__global__
void Interp1PrimFifthOrderWENO_kernel(
    int ngrid_points,
    int ndims,
    int dir,
    int ghosts,
    int nvars,
    int weno_size,
    int offset,
    int stride,
    int upw,
    int uflag,
    const int *dim,
    const double *fC,
    const double *w1,
    const double *w2,
    const double *w3,
    double *fI
)
{
    int p = threadIdx.x + (blockDim.x * blockIdx.x);
    if (p < ngrid_points) {
      int    bounds_inter[GPU_MAX_NDIMS], indexC[GPU_MAX_NDIMS], indexI[GPU_MAX_NDIMS];
      int qm1,qm2,qm3,qp1,qp2;
      const double *ww1, *ww2, *ww3;

      ww1 = w1 + (upw < 0 ? 2*weno_size : 0) + (uflag ? weno_size : 0) + offset;
      ww2 = w2 + (upw < 0 ? 2*weno_size : 0) + (uflag ? weno_size : 0) + offset;
      ww3 = w3 + (upw < 0 ? 2*weno_size : 0) + (uflag ? weno_size : 0) + offset;

      _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
      _ArrayIndexnD_(ndims,p,bounds_inter,indexC,0);
      _ArrayCopy1D_(indexC,indexI,ndims);

      if (upw > 0) {
          indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
          qm3 = qm1 - 2*stride;
          qm2 = qm1 -   stride;
          qp1 = qm1 +   stride;
          qp2 = qm1 + 2*stride;
      } else {
          indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
          qm3 = qm1 + 2*stride;
          qm2 = qm1 +   stride;
          qp1 = qm1 -   stride;
          qp2 = qm1 - 2*stride;
      }

      /* Defining stencil points */
      const double *fm3, *fm2, *fm1, *fp1, *fp2;
      fm3 = (fC+qm3*nvars);
      fm2 = (fC+qm2*nvars);
      fm1 = (fC+qm1*nvars);
      fp1 = (fC+qp1*nvars);
      fp2 = (fC+qp2*nvars);

      /* Candidate stencils and their optimal weights*/
      double f1[GPU_MAX_NVARS], f2[GPU_MAX_NVARS], f3[GPU_MAX_NVARS];

      _ArrayAXBYCZ_(f1,(2*GPU_ONE_SIXTH),fm3,(-7*GPU_ONE_SIXTH) ,fm2,(11*GPU_ONE_SIXTH) ,fm1,nvars);
      _ArrayAXBYCZ_(f2,(-GPU_ONE_SIXTH) ,fm2,(5*GPU_ONE_SIXTH)  ,fm1,(2*GPU_ONE_SIXTH)  ,fp1,nvars);
      _ArrayAXBYCZ_(f3,(2*GPU_ONE_SIXTH),fm1,(5*GPU_ONE_SIXTH)  ,fp1,(-GPU_ONE_SIXTH)   ,fp2,nvars);

      /* calculate WENO weights */
      const double *cur_w1, *cur_w2, *cur_w3;
      cur_w1 = (ww1+p*nvars);
      cur_w2 = (ww2+p*nvars);
      cur_w3 = (ww3+p*nvars);

      _ArrayMultiply3Add1D_((fI+p*nvars),cur_w1,f1,cur_w2,f2,cur_w3,f3,nvars);
    }

    return;
}

extern "C" int gpuInterp1PrimFifthOrderWENO(
    double       *fI,  
    double *fC,  
    double *u,   
    double *x,   
    int    upw,  
    int    dir,  
    void   *s,   
    void   *m,   
    int    uflag 
    )
{
    HyPar           *solver = (HyPar*)          s;
    WENOParameters  *weno   = (WENOParameters*) solver->interp;

    int ghosts = solver->ghosts;
    int ndims  = solver->ndims;
    int nvars  = solver->nvars;
    int *dim   = solver->dim_local;
    int *stride= solver->stride_with_ghosts;

    /* calculate dimension offset */
    int offset = weno->offset[dir];
    int bounds_inter[ndims];
    _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] +=  1;
    int npoints_grid; _ArrayProduct1D_(bounds_inter,ndims,npoints_grid);
    int nblocks = (npoints_grid - 1) / GPU_THREADS_PER_BLOCK + 1;

#if defined(GPU_STAT)
    cudaEvent_t start, stop;
    float milliseconds = 0;
    checkCuda( cudaEventCreate(&start) );
    checkCuda( cudaEventCreate(&stop) );

    int weno_memory_accessed = 3*npoints_grid*nvars*sizeof(double);
    int fI_memory_accessed = npoints_grid*nvars*sizeof(double);
    int fC_memory_accessed = 1;
    for (int d=0; d<ndims; d++) {
      if (d == dir) fC_memory_accessed *= (dim[d]+2*ghosts);
      else          fC_memory_accessed *= dim[d];
    }
    fC_memory_accessed *= nvars*sizeof(double);

    checkCuda( cudaEventRecord(start) );
#endif

    Interp1PrimFifthOrderWENO_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
        npoints_grid, ndims, dir, ghosts, nvars, weno->size, offset, stride[dir], upw, uflag,
        solver->gpu_dim_local, fC, weno->w1, weno->w2, weno->w3, fI
    );
    cudaDeviceSynchronize();

#if defined(GPU_STAT)
    checkCuda( cudaEventRecord(stop) );
    checkCuda( cudaEventSynchronize(stop) );
    checkCuda( cudaEventElapsedTime(&milliseconds, start, stop) );

    printf("%-50s GPU time (secs) = %.6f dir = %d bandwidth (GB/s) = %6.2f\n",
            "Interp1PrimFifthOrderWENO", milliseconds*1e-3, dir,
            (1e-6*(weno_memory_accessed+fI_memory_accessed+fC_memory_accessed))/milliseconds);

    checkCuda( cudaEventDestroy(start) );
    checkCuda( cudaEventDestroy(stop) );
#endif

    return 0;
}

#else

__global__
void Interp1PrimFifthOrderWENO_kernel(
    int npoints_grid,
    int npoints_local_wghosts,
    int ndims,
    int dir,
    int ghosts,
    int nvars,
    int weno_size,
    int offset,
    int stride,
    int upw,
    int uflag,
    const int *dim,
    const double *fC,
    const double *w1,
    const double *w2,
    const double *w3,
    double *fI
)
{
    int p = threadIdx.x + (blockDim.x * blockIdx.x);
    if (p < npoints_grid) {
      int    bounds_inter[GPU_MAX_NDIMS], indexC[GPU_MAX_NDIMS], indexI[GPU_MAX_NDIMS];
      int qm1,qm2,qm3,qp1,qp2;
      const double *ww1, *ww2, *ww3;

      ww1 = w1 + (upw < 0 ? 2*weno_size : 0) + (uflag ? weno_size : 0) + offset;
      ww2 = w2 + (upw < 0 ? 2*weno_size : 0) + (uflag ? weno_size : 0) + offset;
      ww3 = w3 + (upw < 0 ? 2*weno_size : 0) + (uflag ? weno_size : 0) + offset;

      _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
      _ArrayIndexnD_(ndims,p,bounds_inter,indexC,0);
      _ArrayCopy1D_(indexC,indexI,ndims);

      if (upw > 0) {
          indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
          qm3 = qm1 - 2*stride;
          qm2 = qm1 -   stride;
          qp1 = qm1 +   stride;
          qp2 = qm1 + 2*stride;
      } else {
          indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
          qm3 = qm1 + 2*stride;
          qm2 = qm1 +   stride;
          qp1 = qm1 -   stride;
          qp2 = qm1 - 2*stride;
      }

      /* Defining stencil points */
      const double *fm3, *fm2, *fm1, *fp1, *fp2;
      /* Candidate stencils and their optimal weights*/
      double f1, f2, f3;

      int l = p;
      for (int j = 0; j < nvars; j++) {
        fm3 = (fC+qm3);
        fm2 = (fC+qm2);
        fm1 = (fC+qm1);
        fp1 = (fC+qp1);
        fp2 = (fC+qp2);

        f1 = (2*GPU_ONE_SIXTH)*fm3[0] + (-7*GPU_ONE_SIXTH)*fm2[0] + (11*GPU_ONE_SIXTH)*fm1[0];
        f2 = ( -GPU_ONE_SIXTH)*fm2[0] + ( 5*GPU_ONE_SIXTH)*fm1[0] + ( 2*GPU_ONE_SIXTH)*fp1[0];
        f3 = (2*GPU_ONE_SIXTH)*fm1[0] + ( 5*GPU_ONE_SIXTH)*fp1[0] + (  -GPU_ONE_SIXTH)*fp2[0];

        qm3 += npoints_local_wghosts;
        qm2 += npoints_local_wghosts;
        qm1 += npoints_local_wghosts;
        qp1 += npoints_local_wghosts;
        qp2 += npoints_local_wghosts;

        /*
        _ArrayAXBYCZ_(f1,(2*GPU_ONE_SIXTH),fm3,(-7*GPU_ONE_SIXTH) ,fm2,(11*GPU_ONE_SIXTH) ,fm1,nvars);
        _ArrayAXBYCZ_(f2,(-GPU_ONE_SIXTH) ,fm2,(5*GPU_ONE_SIXTH)  ,fm1,(2*GPU_ONE_SIXTH)  ,fp1,nvars);
        _ArrayAXBYCZ_(f3,(2*GPU_ONE_SIXTH),fm1,(5*GPU_ONE_SIXTH)  ,fp1,(-GPU_ONE_SIXTH)   ,fp2,nvars);
        */

        /* calculate WENO weights */
        fI[l] = (ww1+l)[0]*f1 + (ww2+l)[0]*f2 + (ww3+l)[0]*f3;
        l += npoints_grid;
        //_ArrayMultiply3Add1D_((fI+p),cur_w1,f1,cur_w2,f2,cur_w3,f3,1);
      }
    }

    return;
}

extern "C" int gpuInterp1PrimFifthOrderWENO(
    double       *fI,  
    double *fC,  
    double *u,   
    double *x,   
    int    upw,  
    int    dir,  
    void   *s,   
    void   *m,   
    int    uflag 
    )
{
    HyPar           *solver = (HyPar*)          s;
    WENOParameters  *weno   = (WENOParameters*) solver->interp;

    int ghosts = solver->ghosts;
    int ndims  = solver->ndims;
    int nvars  = solver->nvars;
    int *dim   = solver->dim_local;
    int *stride= solver->stride_with_ghosts;

    /* calculate dimension offset */
    int offset = weno->offset[dir];
    int bounds_inter[ndims];
    _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] +=  1;
    int npoints_grid; _ArrayProduct1D_(bounds_inter,ndims,npoints_grid);
    int nblocks = (npoints_grid - 1) / GPU_THREADS_PER_BLOCK + 1;

#if defined(GPU_STAT)
    cudaEvent_t start, stop;
    float milliseconds = 0;
    checkCuda( cudaEventCreate(&start) );
    checkCuda( cudaEventCreate(&stop) );


    int weno_memory_accessed = 3*npoints_grid*nvars*sizeof(double);
    int fI_memory_accessed = npoints_grid*nvars*sizeof(double);
    int fC_memory_accessed = 1;
    for (int d=0; d<ndims; d++) {
      if (d == dir) fC_memory_accessed *= (dim[d]+2*ghosts);
      else          fC_memory_accessed *= dim[d];
    }
    fC_memory_accessed *= nvars*sizeof(double);

    checkCuda( cudaEventRecord(start) );
#endif

    Interp1PrimFifthOrderWENO_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
        npoints_grid, solver->npoints_local_wghosts, ndims, dir, ghosts, nvars, weno->size, offset, stride[dir], upw, uflag,
        solver->gpu_dim_local, fC, weno->w1, weno->w2, weno->w3, fI
    );
#if defined(GPU_STAT)
    checkCuda( cudaEventRecord(stop) );
    checkCuda( cudaEventSynchronize(stop) );
#endif

    cudaDeviceSynchronize();

#if defined(GPU_STAT)
    checkCuda( cudaEventElapsedTime(&milliseconds, start, stop) );

    printf("%-50s GPU time (secs) = %.6f dir = %d bandwidth (GB/s) = %6.2f\n",
            "Interp1PrimFifthOrderWENO2", milliseconds*1e-3, dir,
            (1e-6*(weno_memory_accessed+fI_memory_accessed+fC_memory_accessed))/milliseconds);

    checkCuda( cudaEventDestroy(start) );
    checkCuda( cudaEventDestroy(stop) );
#endif

    return 0;
}

#endif