Skip to content

Commit

Permalink
Include OS7MP in AD-code (MITgcm#814)
Browse files Browse the repository at this point in the history
* show os7md code to autodiff tools, update compatibility check

* use isomip as a test for ad os7mp

* Merge branch 'upstream/master' into os7md_adjoint

* also allow OS7MP in seaice_advection

* remove the SOM-schemes from the list of exclusions

* fix oversight detected by Ou Wang

remove #ifndef ALLOW_AUTODIFF/#endif around ENUM_OS7MP code also in
this routine.

* fix two store directives for som_t/s to reduce memory overhead

- do not store entire fields within bi/bj-loops
- replace a few outdated "byte" parameters by "kind"

* fix comments

* fix and improve error message

* fixing the fix

* un-do changes in input_ad

* Update AD-test isomip.htd

- switch to AdvScheme=7 for both Temp & Salt (so that it get tested)
- reduce T & S diffusivity (by 0.6)
- reduce grdchk_eps to just 1.E-4
- update reference output (from ref. platform villon)

* document including OS7MP in AD-code

---------

Co-authored-by: Jean-Michel Campin <[email protected]>
  • Loading branch information
mjlosch and jm-c authored Mar 29, 2024
1 parent 720a211 commit 598aebf
Show file tree
Hide file tree
Showing 11 changed files with 1,068 additions and 1,372 deletions.
4 changes: 4 additions & 0 deletions doc/tag-index
Original file line number Diff line number Diff line change
@@ -1,6 +1,10 @@
Notes on tags used in MITgcmUV
==============================

o pkg/generic_advdiff:
- adjust 2 SOM (Prather advection) store directives to reduce memory overhead
- include OS7MP (advScheme 7) in AD code ; switch AD test experiment
isomip.htd to use this advection scheme for Temp and Salt.
o pkg/ecco, ctrl:
- put active_read calls within #ifdef ALLOW_AUTODIFF (missing in few places)
and fall back to non-active read (from pkg/rw) if undef;
Expand Down
24 changes: 12 additions & 12 deletions model/src/salt_integrate.F
Original file line number Diff line number Diff line change
Expand Up @@ -217,13 +217,13 @@ SUBROUTINE SALT_INTEGRATE(
ENDDO
#endif /* ALLOW_AUTODIFF */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wFld(:,:,:) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE salt(:,:,:,bi,bj) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE wFld(:,:,:) = comlev1_bibj, key=tkey, kind=isbyte
CADJ STORE salt(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
# ifdef ALLOW_ADAMSBASHFORTH_3
CADJ STORE gsNm(:,:,:,bi,bj,1) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE gsNm(:,:,:,bi,bj,2) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE gsNm(:,:,:,bi,bj,1) = comlev1_bibj, key=tkey, kind=isbyte
CADJ STORE gsNm(:,:,:,bi,bj,2) = comlev1_bibj, key=tkey, kind=isbyte
# else
CADJ STORE gsNm1(:,:,:,bi,bj) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE gsNm1(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
# endif
#endif /* ALLOW_AUTODIFF_TAMC */

Expand All @@ -234,7 +234,7 @@ SUBROUTINE SALT_INTEGRATE(
O kappaRk,
I myThid )
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE kappaRk = comlev1_bibj, key = tkey, byte = isbyte
CADJ STORE kappaRk = comlev1_bibj, key = tkey, kind = isbyte
# endif
#endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */

Expand All @@ -251,8 +251,8 @@ SUBROUTINE SALT_INTEGRATE(
C disable this section of code.
#ifdef GAD_ALLOW_TS_SOM_ADV
# ifdef ALLOW_AUTODIFF_TAMC
C- Note: inefficient storage, storing nSx*nSy times the full som_S array
CADJ STORE som_S = comlev1_bibj, key=tkey, byte=isbyte
C This looks awkward from an F77 point of view, but TAF does the right thing.
CADJ STORE som_S(:,:,:,bi,bj,:) = comlev1_bibj, key=tkey, kind=isbyte
# endif
IF ( saltSOM_Advection ) THEN
# ifdef ALLOW_DEBUG
Expand Down Expand Up @@ -473,16 +473,16 @@ SUBROUTINE SALT_INTEGRATE(
C-- Implicit vertical advection & diffusion

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE gS_loc(:,:,:) = comlev1_bibj, key = tkey, byte = isbyte
CADJ STORE gS_loc(:,:,:) = comlev1_bibj, key = tkey, kind = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef INCLUDE_IMPLVERTADV_CODE
IF ( saltImplVertAdv .OR. implicitDiffusion ) THEN
C to recover older (prior to 2016-10-05) results:
c IF ( saltImplVertAdv ) THEN
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wFld(:,:,:) = comlev1_bibj, key = tkey, byte = isbyte
CADJ STORE salt(:,:,:,bi,bj) = comlev1_bibj, key = tkey, byte = isbyte
CADJ STORE recip_hFac(:,:,:) = comlev1_bibj, key = tkey, byte = isbyte
CADJ STORE wFld(:,:,:) = comlev1_bibj, key = tkey, kind = isbyte
CADJ STORE salt(:,:,:,bi,bj) = comlev1_bibj, key = tkey, kind = isbyte
CADJ STORE recip_hFac(:,:,:) = comlev1_bibj, key = tkey, kind = isbyte
# endif /* ALLOW_AUTODIFF_TAMC */
CALL GAD_IMPLICIT_R(
I saltImplVertAdv, saltVertAdvScheme, GAD_SALINITY,
Expand Down
22 changes: 11 additions & 11 deletions model/src/temp_integrate.F
Original file line number Diff line number Diff line change
Expand Up @@ -225,13 +225,13 @@ SUBROUTINE TEMP_INTEGRATE(
ENDDO
#endif /* ALLOW_AUTODIFF */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wFld(:,:,:) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE wFld(:,:,:) = comlev1_bibj, key=tkey, kind=isbyte
CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
# ifdef ALLOW_ADAMSBASHFORTH_3
CADJ STORE gtNm(:,:,:,bi,bj,1) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE gtNm(:,:,:,bi,bj,2) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE gtNm(:,:,:,bi,bj,1) = comlev1_bibj, key=tkey, kind=isbyte
CADJ STORE gtNm(:,:,:,bi,bj,2) = comlev1_bibj, key=tkey, kind=isbyte
# else
CADJ STORE gtNm1(:,:,:,bi,bj) = comlev1_bibj, key=tkey, byte=isbyte
CADJ STORE gtNm1(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
# endif
#endif /* ALLOW_AUTODIFF_TAMC */

Expand Down Expand Up @@ -259,8 +259,8 @@ SUBROUTINE TEMP_INTEGRATE(
C disable this section of code.
#ifdef GAD_ALLOW_TS_SOM_ADV
# ifdef ALLOW_AUTODIFF_TAMC
C- Note: inefficient storage, storing nSx*nSy times the full som_T array
CADJ STORE som_T = comlev1_bibj, key=tkey, byte=isbyte
C This looks awkward from an F77 point of view, but TAF does the right thing.
CADJ STORE som_T(:,:,:,bi,bj,:) = comlev1_bibj, key=tkey, kind=isbyte
# endif
IF ( tempSOM_Advection ) THEN
# ifdef ALLOW_DEBUG
Expand Down Expand Up @@ -481,16 +481,16 @@ SUBROUTINE TEMP_INTEGRATE(
C-- Implicit vertical advection & diffusion

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE gT_loc (:,:,:) = comlev1_bibj, key = tkey, byte = isbyte
CADJ STORE gT_loc (:,:,:) = comlev1_bibj, key = tkey, kind = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef INCLUDE_IMPLVERTADV_CODE
IF ( tempImplVertAdv .OR. implicitDiffusion ) THEN
C to recover older (prior to 2016-10-05) results:
c IF ( tempImplVertAdv ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE wFld(:,:,:) = comlev1_bibj, key = tkey, byte = isbyte
CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = tkey, byte = isbyte
CADJ STORE recip_hFac(:,:,:) = comlev1_bibj, key = tkey, byte = isbyte
CADJ STORE wFld(:,:,:) = comlev1_bibj, key = tkey, kind = isbyte
CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = tkey, kind = isbyte
CADJ STORE recip_hFac(:,:,:) = comlev1_bibj, key = tkey, kind = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
CALL GAD_IMPLICIT_R(
I tempImplVertAdv, tempVertAdvScheme, GAD_TEMPERATURE,
Expand Down
6 changes: 3 additions & 3 deletions pkg/generic_advdiff/gad_advection.F
Original file line number Diff line number Diff line change
Expand Up @@ -421,11 +421,11 @@ SUBROUTINE GAD_ADVECTION(
CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
I uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
I uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_PPM_NULL_LIMIT .OR.
& advectionScheme.EQ.ENUM_PPM_MONO_LIMIT .OR.
& advectionScheme.EQ.ENUM_PPM_WENO_LIMIT ) THEN
Expand Down Expand Up @@ -642,11 +642,11 @@ SUBROUTINE GAD_ADVECTION(
CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
I vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
I vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_PPM_NULL_LIMIT .OR.
& advectionScheme.EQ.ENUM_PPM_MONO_LIMIT .OR.
& advectionScheme.EQ.ENUM_PPM_WENO_LIMIT ) THEN
Expand Down Expand Up @@ -998,11 +998,11 @@ SUBROUTINE GAD_ADVECTION(
CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
I rTrans, wFld(1-OLx,1-OLy,k), localT3d,
O fVerT(1-OLx,1-OLy,kUp), myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
I rTrans, wFld(1-OLx,1-OLy,k), localT3d,
O fVerT(1-OLx,1-OLy,kUp), myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( usePPMvertAdv .OR. usePQMvertAdv ) THEN
C- copy level from 3d flux data
DO j = 1-OLy,sNy+OLy
Expand Down
8 changes: 0 additions & 8 deletions pkg/generic_advdiff/gad_calc_rhs.F
Original file line number Diff line number Diff line change
Expand Up @@ -285,12 +285,10 @@ SUBROUTINE GAD_CALC_RHS(
CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
I uTrans, uFld, maskLocW, locABT,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
I uTrans, uFld, maskLocW, locABT,
O af, myThid )
#endif
ELSE
STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
ENDIF
Expand Down Expand Up @@ -416,12 +414,10 @@ SUBROUTINE GAD_CALC_RHS(
CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
I vTrans, vFld, maskLocS, locABT,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
I vTrans, vFld, maskLocS, locABT,
O af, myThid )
#endif
ELSE
STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
ENDIF
Expand Down Expand Up @@ -538,12 +534,10 @@ SUBROUTINE GAD_CALC_RHS(
CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
I rTrans, wFld, TracAB,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
I rTrans, wFld, TracAB,
O af, myThid )
#endif
ELSE
STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
ENDIF
Expand Down Expand Up @@ -572,12 +566,10 @@ SUBROUTINE GAD_CALC_RHS(
CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
I rTrans, wFld, TracerN,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
I rTrans, wFld, TracerN,
O af, myThid )
#endif
ELSE
STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
ENDIF
Expand Down
42 changes: 37 additions & 5 deletions pkg/generic_advdiff/gad_check.F
Original file line number Diff line number Diff line change
Expand Up @@ -131,11 +131,43 @@ SUBROUTINE GAD_CHECK( myThid )

C Check compatibility with adjoint
#ifdef ALLOW_AUTODIFF
IF ( tempAdvScheme.EQ.ENUM_OS7MP .OR.
& saltAdvScheme.EQ.ENUM_OS7MP ) THEN
WRITE(msgBuf,'(A,I3,A)') 'GAD_CHECK: advection scheme OS7MP ',
& ENUM_OS7MP,
& ' not yet implemented for adjoint'
IF ( tempAdvScheme.EQ.ENUM_PPM_NULL_LIMIT .OR.
& tempAdvScheme.EQ.ENUM_PPM_MONO_LIMIT .OR.
& tempAdvScheme.EQ.ENUM_PPM_WENO_LIMIT ) THEN
WRITE(msgBuf,'(A,3(I3,A))') 'GAD_CHECK: tempAdvection =',
& ENUM_PPM_NULL_LIMIT, ',', ENUM_PPM_MONO_LIMIT, ', and',
& ENUM_PPM_WENO_LIMIT,
& ' are not yet implemented for adjoint'
CALL PRINT_ERROR( msgBuf, myThid )
errCount = errCount + 1
ENDIF
IF ( tempAdvScheme.EQ.ENUM_PQM_NULL_LIMIT .OR.
& tempAdvScheme.EQ.ENUM_PQM_MONO_LIMIT .OR.
& tempAdvScheme.EQ.ENUM_PQM_WENO_LIMIT ) THEN
WRITE(msgBuf,'(A,3(I3,A))') 'GAD_CHECK: tempAdvection =',
& ENUM_PQM_NULL_LIMIT, ',', ENUM_PQM_MONO_LIMIT, ', and',
& ENUM_PQM_WENO_LIMIT,
& ' are not yet implemented for adjoint'
CALL PRINT_ERROR( msgBuf, myThid )
errCount = errCount + 1
ENDIF
IF ( saltAdvScheme.EQ.ENUM_PPM_NULL_LIMIT .OR.
& saltAdvScheme.EQ.ENUM_PPM_MONO_LIMIT .OR.
& saltAdvScheme.EQ.ENUM_PPM_WENO_LIMIT ) THEN
WRITE(msgBuf,'(A,3(I3,A))') 'GAD_CHECK: saltAdvection =',
& ENUM_PPM_NULL_LIMIT, ',', ENUM_PPM_MONO_LIMIT, ', and',
& ENUM_PPM_WENO_LIMIT,
& ' are not yet implemented for adjoint'
CALL PRINT_ERROR( msgBuf, myThid )
errCount = errCount + 1
ENDIF
IF ( saltAdvScheme.EQ.ENUM_PQM_NULL_LIMIT .OR.
& saltAdvScheme.EQ.ENUM_PQM_MONO_LIMIT .OR.
& saltAdvScheme.EQ.ENUM_PQM_WENO_LIMIT ) THEN
WRITE(msgBuf,'(A,3(I3,A))') 'GAD_CHECK: saltAdvection =',
& ENUM_PQM_NULL_LIMIT, ',', ENUM_PQM_MONO_LIMIT, ', and',
& ENUM_PQM_WENO_LIMIT,
& ' are not yet implemented for adjoint'
CALL PRINT_ERROR( msgBuf, myThid )
errCount = errCount + 1
ENDIF
Expand Down
4 changes: 2 additions & 2 deletions pkg/seaice/seaice_advection.F
Original file line number Diff line number Diff line change
Expand Up @@ -389,11 +389,11 @@ SUBROUTINE SEAICE_ADVECTION(
CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE.,
I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE.,
I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_PPM_NULL_LIMIT .OR.
& advectionScheme.EQ.ENUM_PPM_MONO_LIMIT .OR.
& advectionScheme.EQ.ENUM_PPM_WENO_LIMIT ) THEN
Expand Down Expand Up @@ -603,11 +603,11 @@ SUBROUTINE SEAICE_ADVECTION(
CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE.,
I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE.,
I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
O af, myThid )
#ifndef ALLOW_AUTODIFF
ELSEIF ( advectionScheme.EQ.ENUM_PPM_NULL_LIMIT .OR.
& advectionScheme.EQ.ENUM_PPM_MONO_LIMIT .OR.
& advectionScheme.EQ.ENUM_PPM_WENO_LIMIT ) THEN
Expand Down
4 changes: 2 additions & 2 deletions verification/isomip/code_ad/SIZE.h
Original file line number Diff line number Diff line change
Expand Up @@ -44,8 +44,8 @@ CEOP
PARAMETER (
& sNx = 25,
& sNy = 25,
& OLx = 3,
& OLy = 3,
& OLx = 4,
& OLy = 4,
& nSx = 2,
& nSy = 4,
& nPx = 1,
Expand Down
79 changes: 79 additions & 0 deletions verification/isomip/input_ad.htd/data
Original file line number Diff line number Diff line change
@@ -0,0 +1,79 @@
# ====================
# | Model parameters |
# ====================
#
# Continuous equation parameters
&PARM01
Tref = 30*-1.9,
Sref = 30*34.4,
viscAz=1.E-3,
viscAh=600.0,
no_slip_sides=.FALSE.,
no_slip_bottom=.FALSE.,
diffKhT= 60.0,
diffKhS= 60.0,
#- diffKzT unused when compiled with ALLOW_3D_DIFFKR
#diffKzT=3.E-5,
diffKzS=3.E-5,
bottomDragQuadratic=2.5E-3,
eosType='JMD95Z',
HeatCapacity_Cp = 3974.0,
rhoConst=1030.,
gravity=9.81,
convertFW2Salt = 33.4,
implicitFreeSurface=.TRUE.,
exactConserv=.TRUE.,
hFacMin=0.10,
nonHydrostatic=.FALSE.,
useCDScheme = .TRUE.,
tempAdvScheme = 7,
saltAdvScheme = 7,
#ph(
implicitDiffusion=.TRUE.,
# ivdc_kappa = 7200.,
staggerTimeStep=.TRUE.,
vectorInvariantMomentum=.TRUE.,
nonlinFreeSurf=2,
hFacInf=0.05,
hFacSup=2.0,
#ph)
readBinaryPrec=64,
useSingleCpuIO=.TRUE.,
&

# Elliptic solver parameters
&PARM02
cg2dMaxIters=1000,
cg2dTargetResidual=1.E-13,
&

# Time stepping parameters
&PARM03
nIter0=8640,
nTimeSteps=5,
deltaT=1800.0,
abEps=0.1,
cAdjFreq = 1.,
tauCD = 400000.,
pChkptFreq=0.0,
chkptFreq=0.0,
dumpFreq=0.0,
taveFreq=0.0,
monitorFreq=1.,
monitorSelect=2,
adjMonitorFreq=1800.,
&

# Gridding parameters
&PARM04
usingSphericalPolarGrid=.TRUE.,
ygOrigin = -80.0,
delX=50*0.3,
delY=100*0.1,
delZ=30*30.0,
&

# Input datasets
&PARM05
bathyFile='bathy.box',
&
13 changes: 13 additions & 0 deletions verification/isomip/input_ad.htd/data.grdchk
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
# *******************
# ECCO gradient check
# *******************
&GRDCHK_NML
grdchk_eps = 1.E-4,
iglopos = 20,
jglopos = 24,
kglopos = 20,
# nbeg = 1,
nstep = 1,
nend = 4,
grdchkvarname = "xx_theta",
&
Loading

0 comments on commit 598aebf

Please sign in to comment.