Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ww3 homp #1149

Draft
wants to merge 24 commits into
base: develop
Choose a base branch
from
Draft

Ww3 homp #1149

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
a3e143c
init_bugs2: solve neuman on boundary
aronroland Oct 23, 2023
57ba5f8
Merge branch 'develop' into init_bugs2
aronroland Dec 14, 2023
4cbe073
ug_imp_nmb_acc: Those changes to the number accuracy of the implicit …
aronroland Dec 19, 2023
e0bea58
ug_imp_nmb_acc: add proper settings to limon
aronroland Dec 19, 2023
e60a45e
ug_imp_nmb_acc: add another input file
aronroland Dec 19, 2023
cb7bae8
ww3_homp: add driver of homp explicit
aronroland Dec 20, 2023
cb7dfcd
ww3_homp: update regtest
aronroland Dec 20, 2023
a3b804d
Merge branch 'ug_imp_nmb_acc' into ww3_homp
aronroland Dec 20, 2023
d2009bf
ww3_homp: add explicit and implicit homp
aronroland Dec 20, 2023
c86a567
ww3_homp: consolidiation, extension of omph to implicit solver and so…
aronroland Dec 28, 2023
6f358e2
ww3_homp: some more issues fixed ...
aronroland Dec 28, 2023
bc454a6
Merge branch 'init_bugs2' into ww3_homp
aronroland Dec 28, 2023
f2f6d4d
ww3_homp: clean and merge init_bugs2
aronroland Dec 28, 2023
5482da9
init_bugs2: fix issue #1017 and other issues pointed out by Chris and…
aronroland Jan 10, 2024
9b0083d
Merge branch 'NOAA-EMC:develop' into init_bugs2
thesser1 Jan 11, 2024
e71d161
init_bugs2: fix bug ...
aronroland Jan 16, 2024
57b41c0
Merge branch 'init_bugs2' of https://github.com/erdc/WW3 into init_bugs2
aronroland Jan 16, 2024
b98150e
ww3_homp: revert old regtest
aronroland Jan 16, 2024
712cb55
ww3_homp: add missing end parallel
aronroland Jan 16, 2024
8c6db6c
ww3_homp: regtesting ...
aronroland Jan 16, 2024
6f65346
init_bugs2: remove init emean
aronroland Jan 16, 2024
e517e8d
init_bugs2: thr not defined in db1
aronroland Jan 27, 2024
7f0d537
Merge branch 'develop' into ww3_homp
aronroland Feb 10, 2024
fefdc69
Merge branch 'init_bugs2' into ww3_homp
aronroland Feb 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 180 additions & 26 deletions model/src/w3profsmd_pdlib.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -865,6 +865,120 @@ SUBROUTINE PDLIB_W3XYPUG ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC )
!/ End of W3SPR4 ----------------------------------------------------- /
!/
END SUBROUTINE PDLIB_W3XYPUG

SUBROUTINE PDLIB_W3XYPUG_DRIVER ( FACX, FACY, DTG, VGX, VGY, LCALC )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | |
!/ | Aron Roland (BGS IT&E GmbH) |
!/ | Mathieu Dutour-Sikiric (IRB) |
!/ | |
!/ | FORTRAN 90 |
!/ | Last update : 10-Jan-2011 |
!/ +-----------------------------------+
!/
!/ 10-Jan-2008 : Origination. ( version 3.13 )
!/ 10-Jan-2011 : Addition of implicit scheme ( version 3.14.4 )
!/ 06-Feb-2014 : PDLIB parallelization
!/
! 1. Purpose : Explicit advection schemes driver
!
! Propagation in physical space for a given spectral component.
! Gives the choice of scheme on unstructured grid
! Use the geographical parall algorithms for further speed.
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! Local variables.
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! 5. Called by :
!
! W3WAVE Wave model routine.
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
! make the interface between the WAVEWATCH and the WWM code.
!
! 8. Structure :
!
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
!
! 10. Source code :
!/ ------------------------------------------------------------------- /
!/
!
USE CONSTANTS
!
USE W3TIMEMD, only: DSEC21
!
USE W3GDATMD, only: NX, NY, NSPEC, MAPFS, CLATS, &
FLCX, FLCY, NK, NTH, DTH, XFR, &
ECOS, ESIN, SIG, PFMOVE, &
IOBP, IOBPD, &
FSN, FSPSI, FSFCT, FSNIMP, &
GTYPE, UNGTYPE, NBND_MAP, INDEX_MAP
USE YOWNODEPOOL, only: PDLIB_IEN, PDLIB_TRIA
USE W3GDATMD, only: IOBP_LOC, IOBPD_LOC, IOBPA_LOC, IOBDP_LOC
USE YOWNODEPOOL, only: iplg, npa
USE W3WDATMD, only: TIME, VA
USE W3ODATMD, only: TBPI0, TBPIN, FLBPI
USE W3ADATMD, only: CG, CX, CY, ITIME, DW
USE W3IDATMD, only: FLCUR, FLLEV
USE W3GDATMD, only: NSEAL
USE W3ODATMD, only: IAPROC
USE W3DISPMD, only : WAVNU_LOCAL
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
LOGICAL, INTENT(IN) :: LCALC
!/
!/ ------------------------------------------------------------------- /
!/ Local PARAMETERs
INTEGER :: ISP
!/
!/ ------------------------------------------------------------------- /
!
! 1. Preparations --------------------------------------------------- *
! 1.a Set constants
!
#ifdef W3_S
CALL STRACE (IENT, 'W3XYPUG')
#endif
#ifdef W3_DEBUGSOLVER
WRITE(740+IAPROC,*) 'Begin of PDLIB_W3XYPUG_DRIVER'
FLUSH(740+IAPROC)
#endif

#ifdef W3_O MPH
Copy link
Collaborator

@ukmo-ccbunney ukmo-ccbunney Jan 4, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hey @aronroland - just spotted this potential typo?
Should this be W3_OMPH, not W3_O MPH?

Or was it deliberate to disable that OMP directive?

!$OMP PARALLEL DO PRIVATE (ISP)
#endif
DO ISP=1,NSPEC
CALL PDLIB_W3XYPUG ( ISP, FACX, FACX, DTG, VGX, VGY, LCALC )
END DO

!/
!/ End of PDLIB_W3XYPUG ----------------------------------------------------- /
!/
END SUBROUTINE PDLIB_W3XYPUG_DRIVER
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_W3XYPFSN2(ISP, C, LCALC, RD10, RD20, DT, AC)
!/
Expand Down Expand Up @@ -3644,8 +3758,15 @@ SUBROUTINE calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY)

CCOSA = FACX * ECOS(1:NTH)
CSINA = FACX * ESIN(1:NTH)

call print_memcheck(memunit, 'memcheck_____:'//' WW3_JACOBI SECTION 0')

#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (ISP, ITH, IK, CCOS, CSIN, IP, IP_GLOB, CG1, &
!$OMP& CXY, CXYY, FL11, FL12, FL21, FL22, FL31, FL32, &
!$OMP& CRFS, LAMBDA, K, KP, DELTAL, IB1, IB2, IBR, &
!$OMP& I, J, IE, POS, DTK, I1, I2, I3)
#endif
DO ISP = 1, NSPEC

ITH = 1 + MOD(ISP-1,NTH)
Expand Down Expand Up @@ -3747,6 +3868,10 @@ SUBROUTINE calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY)
END DO
END DO ! ISP

#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif

call print_memcheck(memunit, 'memcheck_____:'//' WW3_JACOBI SECTION 1')
#ifdef W3_DEBUGSOLVER
WRITE(740+IAPROC,*) 'sum(VA)=', sum(VA)
Expand Down Expand Up @@ -4470,28 +4595,33 @@ SUBROUTINE calcARRAY_JACOBI_SPECTRAL_1(DTG)
!AR: TODO: check&report if needed ...
LSIG = FLCUR .OR. FLLEV

#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (IP, ITH, IK, ISP, IP_GLOB, ISEA, eSI, CP_SIG, CM_SIG, B_SIG, &
!$OMP& CAS, DMM, CWNB_M2, CWNB_SIG_M2, DWNI_M2, ITH0, eVal, &
!$OMP& CAD, B_THE, CP_THE, CM_THE)
#endif
DO IP = 1, np
IP_glob=iplg(IP)
ISEA=MAPFS(1,IP_glob)
eSI=PDLIB_SI(IP)
IP_glob = iplg(IP)
ISEA = MAPFS(1,IP_glob)
eSI = PDLIB_SI(IP)
IF (FSFREQSHIFT .AND. LSIG) THEN
IF (FreqShiftMethod .eq. 1) THEN
IF (IOBP_LOC(IP).eq.1.and.IOBDP_LOC(IP).eq.1.and.IOBPA_LOC(IP).eq.0) THEN
CALL PROP_FREQ_SHIFT(IP, ISEA, CAS, DMM, DTG)
CP_SIG = MAX(ZERO,CAS)
CM_SIG = MIN(ZERO,CAS)
B_SIG=0
B_SIG = 0
DO ITH=1,NTH
DO IK=1,NK
ISP=ITH + (IK-1)*NTH
ISP = ITH + (IK-1)*NTH
B_SIG(ISP)= CP_SIG(ISP)/DMM(IK-1) - CM_SIG(ISP)/DMM(IK)
END DO
ISP = ITH + (NK-1)*NTH
B_SIG(ISP)= B_SIG(ISP) + CM_SIG(ISP)/DMM(NK) * FACHFA
ISP = ITH + (NK-1)*NTH ! AR: hmm ...
B_SIG(ISP) = B_SIG(ISP) + CM_SIG(ISP)/DMM(NK) * FACHFA
END DO
ASPAR_JAC(:,PDLIB_I_DIAG(IP))=ASPAR_JAC(:,PDLIB_I_DIAG(IP)) + B_SIG(:)*eSI
ASPAR_JAC(:,PDLIB_I_DIAG(IP)) = ASPAR_JAC(:,PDLIB_I_DIAG(IP)) + B_SIG(:) * eSI
ELSE
CAS=0
CAS = 0
END IF
CAS_SIG(:,IP) = CAS
ELSE IF (FreqShiftMethod .eq. 2) THEN
Expand All @@ -4513,7 +4643,7 @@ SUBROUTINE calcARRAY_JACOBI_SPECTRAL_1(DTG)
ELSE
CWNB_M2 = 0
END IF
CWNB_SIG_M2(:,IP)=CWNB_M2
CWNB_SIG_M2(:,IP) = CWNB_M2
END IF
END IF
!
Expand All @@ -4525,7 +4655,7 @@ SUBROUTINE calcARRAY_JACOBI_SPECTRAL_1(DTG)
! CALL PROP_REFRACTION_PR3(ISEA,DTG,CAD, DoLimiterRefraction)
CALL PROP_REFRACTION_PR3(IP,ISEA,DTG,CAD,DoLimiterRefraction)
ELSE
CAD=ZERO
CAD = ZERO
END IF
#ifdef W3_DEBUGREFRACTION
WRITE(740+IAPROC,*) 'refraction IP=', IP, ' ISEA=', ISEA
Expand All @@ -4538,6 +4668,11 @@ SUBROUTINE calcARRAY_JACOBI_SPECTRAL_1(DTG)
ASPAR_JAC(:,PDLIB_I_DIAG(IP))=ASPAR_JAC(:,PDLIB_I_DIAG(IP)) + B_THE(:)*eSI
END IF
END DO

#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif

END SUBROUTINE calcARRAY_JACOBI_SPECTRAL_1
!/ ------------------------------------------------------------------- /
SUBROUTINE calcARRAY_JACOBI_SPECTRAL_2(DTG,ASPAR_DIAG_LOCAL)
Expand Down Expand Up @@ -5546,26 +5681,26 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
INTEGER :: nbIter, ISPnextDir, ISPprevDir
INTEGER :: ISPp1, ISPm1, JP, ICOUNT1, ICOUNT2
! for the exchange
REAL*8 :: CCOS, CSIN, CCURX, CCURY
REAL*8 :: eSum(NSPEC), FRLOCAL
REAL*8 :: eA_THE, eC_THE, eA_SIG, eC_SIG, eSI
REAL*8 :: CAD(NSPEC), CAS(NSPEC), ACLOC(NSPEC)
REAL*8 :: CP_SIG(NSPEC), CM_SIG(NSPEC)
REAL*8 :: eFactM1, eFactP1
REAL*8 :: Sum_Prev, Sum_New, p_is_converged, DiffNew, prop_conv
REAL*8 :: Sum_L2, Sum_L2_GL
REAL :: CCOS, CSIN, CCURX, CCURY
REAL :: eSum(NSPEC), FRLOCAL
REAL :: eA_THE, eC_THE, eA_SIG, eC_SIG, eSI
REAL :: CAD(NSPEC), CAS(NSPEC), ACLOC(NSPEC)
REAL :: CP_SIG(NSPEC), CM_SIG(NSPEC)
REAL :: eFactM1, eFactP1
REAL :: Sum_Prev, Sum_New, p_is_converged, DiffNew, prop_conv
REAL :: Sum_L2, Sum_L2_GL
REAL :: DMM(0:NK2), DAM(NSPEC), DAM2(NSPEC), SPEC(NSPEC)
REAL*8 :: eDiff(NSPEC), eProd(NSPEC), eDiffB(NSPEC)
REAL*8 :: DWNI_M2(NK), CWNB_M2(1-NTH:NSPEC)
REAL :: eDiff(NSPEC), eProd(NSPEC), eDiffB(NSPEC)
REAL :: DWNI_M2(NK), CWNB_M2(1-NTH:NSPEC)
REAL :: VAnew(NSPEC), VFLWN(1-NTH:NSPEC), JAC, JAC2
REAL :: VAAnew(1-NTH:NSPEC+NTH), VAAacloc(1-NTH:NSPEC+NTH)
REAL :: VAinput(NSPEC), VAacloc(NSPEC), ASPAR_DIAG(NSPEC)
REAL :: aspar_diag_local(nspec), aspar_off_diag_local(nspec), b_jac_local(nspec)
REAL*8 :: eDiffSing, eSumPart
REAL :: eDiffSing, eSumPart
REAL :: EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, U10ABS, U10DIR, TAUA, TAUADIR
REAL :: USTAR, USTDIR, TAUWX, TAUWY, CD, Z0, CHARN, FMEANWS, DLWMEAN
REAL*8 :: eVal1, eVal2
REAL*8 :: eVA, eVO, CG2, NEWDAC, NEWAC, OLDAC, MAXDAC
REAL :: eVal1, eVal2
REAL :: eVA, eVO, CG2, NEWDAC, NEWAC, OLDAC, MAXDAC
REAL :: CG1(0:NK+1), WN1(0:NK+1)
LOGICAL :: LCONVERGED(NSEAL), lexist, LLWS(NSPEC)
#ifdef WEIGHTS
Expand Down Expand Up @@ -5641,6 +5776,10 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
CALL ALL_VA_INTEGRAL_PRINT(IMOD, "VA(np) before transform", 0)
CALL ALL_VA_INTEGRAL_PRINT(IMOD, "VA(npa) before transform", 1)
#endif

#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (JSEA, IP, IP_GLOB, ISEA, ISP, ITH, IK, WN1, CG1)
#endif
DO JSEA=1,NSEAL
IP = JSEA
IP_glob = iplg(IP)
Expand All @@ -5656,6 +5795,10 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
VA(ISP,JSEA) = VA(ISP,JSEA) / CG1(IK) * CLATS(ISEA)
END DO
END DO
#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif

VAOLD = VA(1:NSPEC,1:NSEAL)

#ifdef W3_DEBUGSRC
Expand Down Expand Up @@ -5755,6 +5898,13 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL

call print_memcheck(memunit, 'memcheck_____:'//' WW3_PROP SECTION SOLVER LOOP 1')

#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (IP, ISP, IP_GLOB, ISEA, IK, ITH, CG1, WN1, JSEA, eSI, ACLOC, ESUM, &
!$OMP& ASPAR_DIAG, JP, EPROD, ISPprevDir, ISPnextDir, eA_THE, eC_THE, &
!$OMP& ISPm1, ISPp1, eFactM1, eFactP1, eA_SIG, eC_SIG, CAS, CAD, CP_SIG, DMM, &
!$OMP& CWNB_M2, DWNI_M2, p_is_converged, sum_prev, sum_new, diffnew), REDUCTION(+:is_converged)
#endif

DO IP = 1, np

IP_glob = iplg(IP)
Expand Down Expand Up @@ -5806,7 +5956,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
IF (IMEM == 2) THEN
CALL calcARRAY_JACOBI4(IP,DTG,FACX,FACY,VGX,VGY,ASPAR_DIAG_LOCAL,ASPAR_OFF_DIAG_LOCAL,B_JAC_LOCAL)
ASPAR_DIAG(1:NSPEC) = ASPAR_DIAG_LOCAL(1:NSPEC) + ASPAR_DIAG_ALL(1:NSPEC,IP)
esum = B_JAC_LOCAL - ASPAR_OFF_DIAG_LOCAL + B_JAC(1:NSPEC,IP)
esum(1:NSPEC) = B_JAC_LOCAL - ASPAR_OFF_DIAG_LOCAL + B_JAC(1:NSPEC,IP)
ELSEIF (IMEM == 1) THEN
eSum(1:NSPEC) = B_JAC(1:NSPEC,IP)
ASPAR_DIAG(1:NSPEC) = ASPAR_JAC(1:NSPEC,PDLIB_I_DIAG(IP))
Expand Down Expand Up @@ -6013,6 +6163,10 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
#endif
END DO ! IP

#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif

call print_memcheck(memunit, 'memcheck_____:'//' WW3_PROP SECTION SOLVER LOOP 2')

#ifdef W3_DEBUGSOLVERCOH
Expand Down Expand Up @@ -6139,7 +6293,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
#endif
EXIT
END IF
END IF
END IF ! TERMINATE NORM
call print_memcheck(memunit, 'memcheck_____:'//' WW3_PROP SECTION SOLVER LOOP 6')

nbiter = nbiter + 1
Expand Down
6 changes: 2 additions & 4 deletions model/src/w3wavemd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -448,7 +448,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT &
USE W3IOSFMD
#ifdef W3_PDLIB
USE PDLIB_W3PROFSMD, only : APPLY_BOUNDARY_CONDITION_VA
USE PDLIB_W3PROFSMD, only : PDLIB_W3XYPUG, PDLIB_W3XYPUG_BLOCK_IMPLICIT, PDLIB_W3XYPUG_BLOCK_EXPLICIT
USE PDLIB_W3PROFSMD, only : PDLIB_W3XYPUG_DRIVER, PDLIB_W3XYPUG_BLOCK_IMPLICIT, PDLIB_W3XYPUG_BLOCK_EXPLICIT
USE PDLIB_W3PROFSMD, only : ALL_VA_INTEGRAL_PRINT, ALL_VAOLD_INTEGRAL_PRINT, ALL_FIELD_INTEGRAL_PRINT
USE W3PARALL, only : PDLIB_NSEAL, PDLIB_NSEALM
USE yowNodepool, only: npa, iplg, np
Expand Down Expand Up @@ -1816,9 +1816,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT &
#ifdef W3_PDLIB
IF (FLCX .or. FLCY) THEN
IF (.NOT. FSTOTALIMP .AND. .NOT. FSTOTALEXP) THEN
DO ISPEC=1,NSPEC
CALL PDLIB_W3XYPUG ( ISPEC, FACX, FACX, DTG, VGX, VGY, UGDTUPDATE )
END DO
CALL PDLIB_W3XYPUG_DRIVER ( FACX, FACX, DTG, VGX, VGY, UGDTUPDATE )
END IF
END IF
#endif
Expand Down
Loading