SELF_SetInterface_3D.cpp Source File

/*
! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
!
! Maintainers : support@fluidnumerics.com
! Official Repository : https://github.com/FluidNumerics/self/
!
! Copyright © 2024 Fluid Numerics LLC
!
! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in
!    the documentation and/or other materials provided with the distribution.
!
! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from
!    this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
*/

#include "SELF_GPU_Macros.h"

// Element-set interface traces (see SELF_SetInterface_3D_t.f90).
//
// The outbox and ghost buffers are (Np,Np,nvarMax,nFaces) with Np = N+1, column-major,
// so node (i,j), variable v of face f (all 0-based) sits at
//   i + Np*(j + Np*(v + nvarMax*f)).
// The faces of one set are f0 .. f0+nFaces-1; elems[f] and sides[f] are 1-based
// element and side indices in that set's sub-mesh, whose traces have nel elements.

// outbox(i,j,v,f) = trace(i,j,s,e,v): each face's own trace, unflipped
__global__ void GatherSetInterfaceOutbox_kernel(real *outbox, real *trace, int *elems, int *sides,
                                                int f0, int nFaces, int nvar, int nvarMax, int N, int nel){
  uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x;
  uint32_t np = N+1;
  uint32_t total = nFaces*nvar*np*np;

  if(idof < total){
    uint32_t i = idof % np;
    uint32_t j = (idof/np) % np;
    uint32_t v = (idof/(np*np)) % nvar;
    uint32_t f = f0 + idof/(np*np*nvar);
    uint32_t e = elems[f]-1;
    uint32_t s = sides[f]-1;
    outbox[i + np*(j + np*(v + nvarMax*f))] = trace[SCB_3D_INDEX(i,j,s,e,v,N,nel)];
  }
}

// extBoundary(i,j,s,e,v) = ghost(i,j,v,f): the partner trace, already flipped on the host
__global__ void FillSetInterfaceBoundary_kernel(real *extBoundary, real *ghost, int *elems, int *sides,
                                                int f0, int nFaces, int nvar, int nvarMax, int N, int nel){
  uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x;
  uint32_t np = N+1;
  uint32_t total = nFaces*nvar*np*np;

  if(idof < total){
    uint32_t i = idof % np;
    uint32_t j = (idof/np) % np;
    uint32_t v = (idof/(np*np)) % nvar;
    uint32_t f = f0 + idof/(np*np*nvar);
    uint32_t e = elems[f]-1;
    uint32_t s = sides[f]-1;
    extBoundary[SCB_3D_INDEX(i,j,s,e,v,N,nel)] = ghost[i + np*(j + np*(v + nvarMax*f))];
  }
}

extern "C"
{
  void GatherSetInterfaceOutbox_gpu(real *outbox, real *trace, int *elems, int *sides,
                                    int f0, int nFaces, int nvar, int nvarMax, int N, int nel)
  {
    int threads_per_block = 256;
    int total = nFaces*nvar*(N+1)*(N+1);
    int nblocks_x = total/threads_per_block + 1;
    GatherSetInterfaceOutbox_kernel<<<dim3(nblocks_x,1,1), dim3(threads_per_block,1,1), 0, 0>>>(
      outbox, trace, elems, sides, f0, nFaces, nvar, nvarMax, N, nel);
  }

  void FillSetInterfaceBoundary_gpu(real *extBoundary, real *ghost, int *elems, int *sides,
                                    int f0, int nFaces, int nvar, int nvarMax, int N, int nel)
  {
    int threads_per_block = 256;
    int total = nFaces*nvar*(N+1)*(N+1);
    int nblocks_x = total/threads_per_block + 1;
    FillSetInterfaceBoundary_kernel<<<dim3(nblocks_x,1,1), dim3(threads_per_block,1,1), 0, 0>>>(
      extBoundary, ghost, elems, sides, f0, nFaces, nvar, nvarMax, N, nel);
  }
}