Initial device port of sponge BC

.gitignore

@@ -18,3 +18,4 @@ ctags.txt
 reggie_run/
 doc/userguide/userguide.pdf
 src/commit.h
+.vscode/

src/dg/dg.f90

@@ -397,7 +397,13 @@ SUBROUTINE DGTimeDerivative_weakForm(t)
 ! 13. Compute source terms and sponge (in physical space, conversion to reference space inside routines)
 ! TODO: This can be used for latency hiding or not?
 !IF(doCalcSource) CALL CalcSource(Ut,t)
-!IF(doSponge)     CALL Sponge(Ut)
+IF(doSponge) THEN
+!#if USE_GPU
+  CALL Sponge_GPU(d_Ut)
+!#else
+! CALL Sponge(Ut)
+!#endif
+ENDIF
 !IF(doTCSource)   CALL TestcaseSource(Ut)
 
 !! 14. apply Jacobian

src/sponge/sponge.f90

@@ -200,6 +200,7 @@ SUBROUTINE InitSponge
 ! Preparation of the baseflow on each Gauss Point
 SWRITE(UNIT_stdOut,'(A)') '  Initialize Sponge Base Flow...'
 ALLOCATE(SpBaseFlow(PP_nVar,0:PP_N,0:PP_N,0:PP_NZ,nElems))
+!$cuf ALLOCATE( d_SpBaseFlow(PP_nVar,0:PP_N,0:PP_N,0:PP_NZ,nElems) )
 SELECT CASE(SpBaseflowType)
 CASE(SPONGEBASEFLOW_CONSTANT) ! constant baseflow from refstate
   DO iElem=1,nElems
@@ -234,6 +235,11 @@ SUBROUTINE InitSponge
   END IF
 END SELECT
 
+! If using GPU accleration, update base flow array on the device
+!#if USE_GPU
+d_SpBaseFlow = SpBaseFlow
+!#endif
+
 SWRITE(UNIT_stdOut,'(A)')' INIT SPONGE DONE!'
 SWRITE(UNIT_stdOut,'(132("-"))')
 
@@ -364,6 +370,8 @@ SUBROUTINE CalcSpongeRamp()
 nSpongeElems=COUNT(applySponge)
 ALLOCATE(SpongeMat(0:PP_N,0:PP_N,0:PP_NZ,nSpongeElems))
 ALLOCATE(SpongeMap(nSpongeElems))
+!$cuf ALLOCATE( d_SpongeMat(0:PP_N,0:PP_N,0:PP_NZ,nSpongeElems) )
+!$cuf ALLOCATE( d_SpongeMap(nSpongeElems) )
 iSpongeElem=0
 DO iElem=1,nElems
   IF(applySponge(iElem))THEN
@@ -372,6 +380,11 @@ SUBROUTINE CalcSpongeRamp()
   END IF
 END DO
 
+! If using GPU accleration, update SpongeMap on the device
+!#if USE_GPU
+d_spongeMap = spongeMap
+!#endif
+
 ! Calculate the final sponge strength in the sponge region
 SpongeMat=0.
 DO iSpongeElem=1,nSpongeElems
@@ -400,6 +413,11 @@ SUBROUTINE CalcSpongeRamp()
   SpongeMat(:,:,:,iSpongeElem) = damping*sigma(:,:,:)
 END DO !iSpongeElem=1,nSpongeElems
 
+! If using GPU accleration, update SpongeMat on the device
+!#if USE_GPU
+d_SpongeMat = SpongeMat
+!#endif
+
 DEALLOCATE(SpongeShape)
 DEALLOCATE(SpDistance)
 DEALLOCATE(xStart)
@@ -557,7 +575,7 @@ SUBROUTINE Sponge(Ut)
   ELSE
 #endif
     DO k=0,PP_NZ; DO j=0,PP_N; DO i=0,PP_N
-      Ut(:,i,j,k,iElem) = Ut(:,i,j,k,iElem) - SpongeMat(   i,j,k,iSpongeElem) * &
+      Ut(:,i,j,k,iElem) = d_Ut(:,i,j,k,iElem) - SpongeMat(   i,j,k,iSpongeElem) * &
                           (U(:,i,j,k,iElem) - SpBaseFlow(:,i,j,k,iElem))
     END DO; END DO; END DO
 #if FV_ENABLED
@@ -568,6 +586,45 @@ SUBROUTINE Sponge(Ut)
 END SUBROUTINE Sponge
 
 
+!==================================================================================================================================
+!> \brief Apply the sponge to the solution vector on the device (compute contribution to d_Ut).
+!>
+!> GPU accelerated version of "Sponge" subroutine above.
+!> \f$ U_t = U_t - \sigma(x)*(U-U_B) \f$, where \f$ \sigma(x) \f$ is the sponge strength and \f$ U_B \f$ is the base flow.
+!> The operation will be performed in the sponge region only using the sponge mapping. The sponge is already pre-multiplied
+!> by the Jacobian since we are working in the reference space at this point (at the end of DGTimeDerivative_weakForm).
+!==================================================================================================================================
+SUBROUTINE Sponge_GPU(d_Ut)
+  ! MODULES
+  USE MOD_Globals
+  USE MOD_PreProc
+  !$cuf USE MOD_Sponge_Vars, ONLY: d_SpongeMat, d_SpBaseFlow, d_SpongeMap
+  !$cuf USE MOD_DG_Vars     ,ONLY: d_U
+  USE MOD_Mesh_Vars   ,ONLY: nElems
+  IMPLICIT NONE
+  !----------------------------------------------------------------------------------------------------------------------------------
+  ! INPUT/OUTPUT VARIABLES
+  REAL,INTENT(INOUT)  :: d_Ut(PP_nVar,0:PP_N,0:PP_N,0:PP_NZ,nElems) !< DG solution time derivative
+  !$cuf ATTRIBUTES(DEVICE) :: d_Ut (PP_nVar,0:PP_N,0:PP_N,0:PP_NZ,nElems) 
+  !----------------------------------------------------------------------------------------------------------------------------------
+  ! LOCAL VARIABLES
+  INTEGER             :: iElem,iSpongeElem,i,j,k
+  !==================================================================================================================================
+
+  !$cuf kernel do (4) <<< *, * >>>
+  DO iSpongeElem=1,nSpongeElems
+    iElem=d_spongeMap(iSpongeElem)
+
+      DO k=0,PP_NZ; DO j=0,PP_N; DO i=0,PP_N
+        d_Ut(:,i,j,k,iElem) = d_Ut(:,i,j,k,iElem) - d_SpongeMat(   i,j,k,iSpongeElem) * &
+                            (d_U(:,i,j,k,iElem) - d_SpBaseFlow(:,i,j,k,iElem))
+      END DO; END DO; END DO
+
+  END DO
+
+END SUBROUTINE Sponge_GPU
+
+
 !==================================================================================================================================
 !> Deallocate sponge arrays
 !==================================================================================================================================

src/sponge/sponge_vars.f90

@@ -22,13 +22,16 @@ MODULE MOD_Sponge_Vars
 !----------------------------------------------------------------------------------------------------------------------------------
 ! GLOBAL VARIABLES
 !----------------------------------------------------------------------------------------------------------------------------------
-LOGICAL          :: doSponge              !< Turn on to employ sponge regions for reducing reflections at boundaries
-LOGICAL          :: SpongeViz             !< Turn on to write a visualization file of the sponge region and strength
-LOGICAL          :: CalcPruettDamping=.FALSE. !< true if temporally varying, solution adaptive Pruett baseflow is used
-INTEGER          :: nSpongeElems          !< number of elements for which sponge is applied
-INTEGER,ALLOCATABLE :: spongeMap(:)       !< mapping from Elem -> spongElem
-REAL             :: damping               !< Strenght of damping per face
-REAL,ALLOCATABLE :: SpongeMat(:,:,:,:)    !< precomputed sponge functions per DOF and sponge elem
-REAL,ALLOCATABLE,TARGET :: SpBaseFlow(:,:,:,:,:) !< precompute global reference state for whole field
+LOGICAL                             :: doSponge                     !< Turn on to employ sponge regions for reducing reflections at boundaries
+LOGICAL                             :: SpongeViz                    !< Turn on to write a visualization file of the sponge region and strength
+LOGICAL                             :: CalcPruettDamping=.FALSE.    !< true if temporally varying, solution adaptive Pruett baseflow is used
+INTEGER                             :: nSpongeElems                 !< number of elements for which sponge is applied
+INTEGER,ALLOCATABLE                 :: spongeMap(:)                 !< mapping from Elem -> spongElem
+!$cuf INTEGER,ALLOCATABLE,DEVICE    :: d_spongeMap(:)               !< copy of spongeMap on the device
+REAL                                :: damping                      !< Strenght of damping per face
+REAL,ALLOCATABLE                    :: SpongeMat(:,:,:,:)           !< precomputed sponge functions per DOF and sponge elem
+!$cuf REAL,ALLOCATABLE,DEVICE       :: d_SpongeMat(:,:,:,:)         !< copy of SpongeMat on the device
+REAL,ALLOCATABLE,TARGET             :: SpBaseFlow(:,:,:,:,:)        !< precompute global reference state for whole field
+!$cuf REAL,ALLOCATABLE,TARGET,DEVICE:: d_SpBaseFlow(:,:,:,:,:)      !< copy of SpBaseFlow on the device
 !==================================================================================================================================
 END MODULE MOD_Sponge_Vars