64-BitBrainstorm_2026/AMSS-NCKU
8
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add OpenMP parallelization to BSSN RHS hot-path stencil routines

Browse Source 
Enable OpenMP threading for the dominant computational kernels:
- makefile.inc: add -qopenmp to f90appflags
- diff_new.f90: split fderivs/fdderivs into OpenMP interior + serial boundary
- kodiss.f90: split kodis into OpenMP interior + serial boundary
- lopsidediff.f90: add OMP PARALLEL DO COLLAPSE(2) to lopsided
- fmisc.f90: parallelize symmetry_bd bulk array copy
- bssn_rhs.f90: add OMP WORKSHARE to array-syntax operations

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
This commit is contained in:
CGH0S7
2026-02-07 13:58:55 +08:00
parent 09ffdb553d
commit c6e4d4ab71
6 changed files with 162 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

99
AMSS_NCKU_source/diff_new.f90
View File

@@ -997,11 +997,11 @@
fy = ZEO
fz = ZEO
#if 0
do k=1,ex(3)-1
do j=1,ex(2)-1
do i=1,ex(1)-1
#if 0
! x direction
! x direction
if(i+2 <= imax .and. i-2 >= imin)then
!
! f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
@@ -1018,7 +1018,7 @@
! set imax and imin 0
endif
! y direction
! y direction
if(j+2 <= jmax .and. j-2 >= jmin)then
fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
@@ -1029,7 +1029,7 @@
! set jmax and jmin 0
endif
! z direction
! z direction
if(k+2 <= kmax .and. k-2 >= kmin)then
fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
@@ -1040,9 +1040,13 @@
! set kmax and kmin 0
endif
enddo
enddo
enddo
#elif 0
! x direction
if(i+2 <= imax .and. i-2 >= imin)then
do k=1,ex(3)-1
do j=1,ex(2)-1
do i=1,ex(1)-1
!
! f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
! fx(i) = ---------------------------------------------
@@ -1079,8 +1083,32 @@
! set kmax and kmin 0
endif
enddo
enddo
enddo
#else
! for bam comparison
! for bam comparison — split into branch-free interior + serial boundary
! Interior: all stencil points guaranteed in-bounds, no branches needed
!$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
do k=max(3,1),min(ex(3)-1,kmax-2)
do j=max(3,1),min(ex(2)-1,jmax-2)
!DIR$ IVDEP
do i=max(3,1),min(ex(1)-1,imax-2)
fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
enddo
enddo
enddo
!$OMP END PARALLEL DO
! Boundary shell: original branching logic for points near edges
do k=1,ex(3)-1
do j=1,ex(2)-1
do i=1,ex(1)-1
if(i >= 3 .and. i <= imax-2 .and. &
j >= 3 .and. j <= jmax-2 .and. &
k >= 3 .and. k <= kmax-2) cycle
if(i+2 <= imax .and. i-2 >= imin .and. &
j+2 <= jmax .and. j-2 >= jmin .and. &
k+2 <= kmax .and. k-2 >= kmin) then
@@ -1094,10 +1122,10 @@
fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
endif
enddo
enddo
enddo
#endif
enddo
enddo
enddo
return
@@ -1401,10 +1429,10 @@
fxz = ZEO
fyz = ZEO
#if 0
do k=1,ex(3)-1
do j=1,ex(2)-1
do i=1,ex(1)-1
#if 0
!~~~~~~ fxx
if(i+2 <= imax .and. i-2 >= imin)then
!
@@ -1481,9 +1509,48 @@
- (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2)))
elseif(j+1 <= jmax .and. j-1 >= jmin .and. k+1 <= kmax .and. k-1 >= kmin)then
fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
endif
endif
enddo
enddo
enddo
#else
! for bam comparison
! for bam comparison — split into branch-free interior + serial boundary
! Interior: all stencil points guaranteed in-bounds, no branches needed
!$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
do k=max(3,1),min(ex(3)-1,kmax-2)
do j=max(3,1),min(ex(2)-1,jmax-2)
!DIR$ IVDEP
do i=max(3,1),min(ex(1)-1,imax-2)
fxx(i,j,k) = Fdxdx*(-fh(i-2,j,k)+F16*fh(i-1,j,k)-F30*fh(i,j,k) &
-fh(i+2,j,k)+F16*fh(i+1,j,k) )
fyy(i,j,k) = Fdydy*(-fh(i,j-2,k)+F16*fh(i,j-1,k)-F30*fh(i,j,k) &
-fh(i,j+2,k)+F16*fh(i,j+1,k) )
fzz(i,j,k) = Fdzdz*(-fh(i,j,k-2)+F16*fh(i,j,k-1)-F30*fh(i,j,k) &
-fh(i,j,k+2)+F16*fh(i,j,k+1) )
fxy(i,j,k) = Fdxdy*( (fh(i-2,j-2,k)-F8*fh(i-1,j-2,k)+F8*fh(i+1,j-2,k)-fh(i+2,j-2,k)) &
-F8 *(fh(i-2,j-1,k)-F8*fh(i-1,j-1,k)+F8*fh(i+1,j-1,k)-fh(i+2,j-1,k)) &
+F8 *(fh(i-2,j+1,k)-F8*fh(i-1,j+1,k)+F8*fh(i+1,j+1,k)-fh(i+2,j+1,k)) &
- (fh(i-2,j+2,k)-F8*fh(i-1,j+2,k)+F8*fh(i+1,j+2,k)-fh(i+2,j+2,k)))
fxz(i,j,k) = Fdxdz*( (fh(i-2,j,k-2)-F8*fh(i-1,j,k-2)+F8*fh(i+1,j,k-2)-fh(i+2,j,k-2)) &
-F8 *(fh(i-2,j,k-1)-F8*fh(i-1,j,k-1)+F8*fh(i+1,j,k-1)-fh(i+2,j,k-1)) &
+F8 *(fh(i-2,j,k+1)-F8*fh(i-1,j,k+1)+F8*fh(i+1,j,k+1)-fh(i+2,j,k+1)) &
- (fh(i-2,j,k+2)-F8*fh(i-1,j,k+2)+F8*fh(i+1,j,k+2)-fh(i+2,j,k+2)))
fyz(i,j,k) = Fdydz*( (fh(i,j-2,k-2)-F8*fh(i,j-1,k-2)+F8*fh(i,j+1,k-2)-fh(i,j+2,k-2)) &
-F8 *(fh(i,j-2,k-1)-F8*fh(i,j-1,k-1)+F8*fh(i,j+1,k-1)-fh(i,j+2,k-1)) &
+F8 *(fh(i,j-2,k+1)-F8*fh(i,j-1,k+1)+F8*fh(i,j+1,k+1)-fh(i,j+2,k+1)) &
- (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2)))
enddo
enddo
enddo
!$OMP END PARALLEL DO
! Boundary shell: original branching logic for points near edges
do k=1,ex(3)-1
do j=1,ex(2)-1
do i=1,ex(1)-1
if(i >= 3 .and. i <= imax-2 .and. &
j >= 3 .and. j <= jmax-2 .and. &
k >= 3 .and. k <= kmax-2) cycle
if(i+2 <= imax .and. i-2 >= imin .and. &
j+2 <= jmax .and. j-2 >= jmin .and. &
k+2 <= kmax .and. k-2 >= kmin) then
@@ -1518,10 +1585,10 @@
fxz(i,j,k) = Sdxdz*(fh(i-1,j,k-1)-fh(i+1,j,k-1)-fh(i-1,j,k+1)+fh(i+1,j,k+1))
fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
endif
enddo
enddo
enddo
#endif
enddo
enddo
enddo
return