Back to home page

MITgcm
 
 

    

File indexing completed on 2021-11-06 05:17:09 UTC

view on githubraw file Latest commit 1c2b1fa3 on 2021-10-19 18:24:13 UTC
aecc8b0f47 Mart*0001 #include "AUTODIFF_OPTIONS.h"
                0002 #ifdef ALLOW_CTRL
                0003 # include "CTRL_OPTIONS.h"
                0004 #endif
                0005 
                0006 C--   File cg2d_mad.F: Code for hand-written (manual) adjoint of cg2d
                0007 C--    Contents
                0008 C--   CG2D_MAD        - computes the adjoint of cg2d
                0009 C--   CG2D_STORE      - saves input field to or restores from a common block
                0010 
                0011 CBOP
                0012 C     !ROUTINE: CG2D_MAD
                0013 C     !INTERFACE:
                0014       SUBROUTINE CG2D_MAD(
                0015      U                    cg2d_b_ad,
                0016      U                    cg2d_x_ad,
                0017      U                    numIters, nIterMin,
                0018      I                    myThid )
                0019 C     !DESCRIPTION: \bv
                0020 C     *==========================================================*
                0021 C     | SUBROUTINE CG2D_MAD
                0022 C     | o Manual ADjoint version of CG2D, the two-dimensional
                0023 C     |   grid problem conjugate-gradient inverter.
                0024 C     *==========================================================*
                0025 C     | This routine is called from solve_for_pressure_ad that
                0026 C     | is generated by TAF/TAMC. We want a self-adjoint cg2d
                0027 C     | to avoid too many complications, so we provide flow
                0028 C     | directives in cg2d.flow and add the missing AD
                0029 C     | contributions manually (hence the name mad=Manuad ADjoint).
                0030 C     |
                0031 C     | This routine calls cg2d with reversed order of argumgent
                0032 C     | cg2d_x_ad and cg2d_b_ad, and sets dependendcies to zero
                0033 C     | that a non-self-adjoint cg2d would have produced.
                0034 C     | In addition, when runtime parameter cg2dFullAdjoint = T,
                0035 C     | it computes sensitivities to matrix coefficients. This
                0036 C     | makes it the full adjoint of cg2d, only if cg2d converges
                0037 C     | completely. Otherwise it is an excellent approximation.
                0038 C     *==========================================================*
                0039 C     \ev
                0040 
                0041 C     !USES:
                0042       IMPLICIT NONE
                0043 C     === Global data ===
                0044 #include "SIZE.h"
                0045 #include "EEPARAMS.h"
                0046 #include "PARAMS.h"
                0047 #if ( defined NONLIN_FRSURF || defined ALLOW_DEPTH_CONTROL )
                0048 # include "CG2D.h"
                0049 #endif
                0050 #ifdef ALLOW_AUTODIFF_TAMC
                0051 # include "AUTODIFF_PARAMS.h"
                0052 #endif
                0053 
                0054 C     !INPUT/OUTPUT PARAMETERS:
                0055 C     === Routine arguments ===
                0056 C     cg2d_b_ad :: Adjoint of source term (on input = 0)
                0057 C     cg2d_x_ad :: Adjoint of the solution
                0058 C     firstResidual :: the initial residual before any iterations
                0059 C     minResidualSq :: the lowest residual reached (squared)
                0060 C     lastResidual  :: the actual residual reached
                0061 C     numIters  :: Inp: the maximum number of iterations allowed
                0062 C                  Out: the actual number of iterations used
                0063 C     nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
                0064 C                  Out: iteration number corresponding to lowest residual
                0065 C     myThid    :: Thread on which I am working.
                0066       _RL  cg2d_b_ad(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0067       _RL  cg2d_x_ad(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0068 c     _RL  firstResidual, minResidualSq, lastResidual
                0069       INTEGER numIters
                0070       INTEGER nIterMin
                0071       INTEGER myThid
                0072 
                0073 #ifdef ALLOW_AUTODIFF_TAMC
                0074 #if ( defined NONLIN_FRSURF || defined ALLOW_DEPTH_CONTROL )
                0075 C     directly imported from TAF-generated code, make sure that they are
                0076 C     consistent with what is found in S/R update_cg2d_ad
f71d880987 Jean*0077       _RS aW2d_ad(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0078       _RS aS2d_ad(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0079       _RS aC2d_ad(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0080       _RS pW_ad  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0081       _RS pS_ad  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0082       _RS pC_ad  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1c2b1fa3bc Mart*0083       COMMON /ADCG2D_I_RS/
aecc8b0f47 Mart*0084      &     aW2d_ad, aS2d_ad, aC2d_ad, pW_ad, pS_ad, pC_ad
                0085 #endif /* NONLIN_FRSURF or ALLOW_DEPTH_CONTROL */
                0086 
                0087 C     !LOCAL VARIABLES:
                0088 C     === Local variables ===
                0089 C     bi,bj  :: tile indices
                0090 C     i,j    :: Loop counters
                0091 C     cg2d_x :: local copy of cg2d_x to be restored from "tape"
                0092       INTEGER i,j,bi,bj
                0093       _RL  firstResidual, minResidualSq, lastResidual
                0094 #if ( defined NONLIN_FRSURF || defined ALLOW_DEPTH_CONTROL )
                0095       INTEGER numItersFwd, nIterMinFwd
                0096       _RL  recip_cg2dNorm
                0097       _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0098 #endif
                0099       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0100 CEOP
                0101 
                0102 C--   mark output of cg2d as being called from here
                0103       IF ( debugLevel .GE. debLevZero ) THEN
                0104        _BEGIN_MASTER( myThid )
44d3986245 Jean*0105        WRITE(standardMessageUnit,'(A)')
aecc8b0f47 Mart*0106      &      ' Calling cg2d from S/R CG2D_MAD'
                0107        _END_MASTER( myThid )
                0108       ENDIF
                0109 
                0110 #if ( defined NONLIN_FRSURF || defined ALLOW_DEPTH_CONTROL )
                0111 C     make local copies because we need them later
                0112       numItersFwd = numIters
                0113       nIterMinFwd = nIterMin
                0114 #endif
                0115 
                0116 #ifdef ALLOW_SRCG
                0117       IF ( useSRCGSolver ) THEN
                0118 C--   Call the single reduce CG solver
                0119        CALL CG2D_SR(
                0120      U      cg2d_x_ad, cg2d_b_ad,
                0121      O      firstResidual, minResidualSq, lastResidual,
                0122      U      numIters, nIterMin,
                0123      I      myThid )
                0124       ELSE
                0125 #else
                0126       IF (.TRUE.) THEN
                0127 #endif /* ALLOW_SRCG */
                0128 C--   Call the default CG solver
                0129        CALL CG2D(
                0130      U      cg2d_x_ad, cg2d_b_ad,
                0131      O      firstResidual, minResidualSq, lastResidual,
                0132      U      numIters, nIterMin,
                0133      I      myThid )
                0134       ENDIF
                0135 
                0136 CML      CALL ADEXCH_XY_RL( cg2d_b_ad, myThid )
                0137 CML       _EXCH_XY_RL( cg2d_b_ad, myThid )
                0138 
                0139 C     cg2d_x_ad is reset to zero, because we assume that the solution of
                0140 C     the self-adjoint solver does not depend in the intial conditions
                0141 C     so that no sensitivities should be passed back to the caller.
                0142       DO bj=myByLo(myThid),myByHi(myThid)
                0143        DO bi=myBxLo(myThid),myBxHi(myThid)
                0144         DO j=1-OLy,sNy+OLy
                0145          DO i=1-OLx,sNx+OLx
                0146           cg2d_x_ad(i,j,bi,bj) = 0. _d 0
                0147          ENDDO
                0148         ENDDO
                0149        ENDDO
                0150       ENDDO
                0151 
                0152 C--   dump CG2D output
                0153 C     Since there is no clean way in importing myTime into this routine,
                0154 C     we just call it every timestep but only for a higher debugLevel
                0155       IF ( debugLevel .GE. debLevC ) THEN
                0156        _BEGIN_MASTER( myThid )
                0157        WRITE(msgBuf,'(A30,1PE23.14)')
                0158      &       'CG2D_MAD cg2d_init_res =',firstResidual
                0159        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0160        WRITE(msgBuf,'(A37,2I8)')
                0161      &       'CG2D_MAD: cg2d_iters(min,last) =', nIterMin, numIters
                0162        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0163        IF ( minResidualSq.GE.0. ) THEN
                0164         minResidualSq = SQRT(minResidualSq)
                0165         WRITE(msgBuf,'(A30,1PE23.14)')
                0166      &        'CG2D_MAD: cg2d_min_res  =',minResidualSq
                0167         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0168        ENDIF
                0169        WRITE(msgBuf,'(A30,1PE23.14)')
                0170      &       'CG2D_MAD: cg2d_last_res =',lastResidual
                0171        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0172        _END_MASTER( myThid )
                0173       ENDIF
                0174 
                0175 #if ( defined NONLIN_FRSURF || defined ALLOW_DEPTH_CONTROL )
                0176       recip_cg2dNorm = 1.0 _d 0
                0177       IF ( cg2dNorm .NE. 0. _d 0 ) THEN
                0178        recip_cg2dNorm = 1./cg2dNorm
                0179       ENDIF
                0180 
                0181       IF ( cg2dFullAdjoint ) THEN
                0182 C     restore cg2d_x
                0183        CALL CG2D_STORE( cg2d_x, .FALSE., myThid )
                0184 
                0185 C     Compute derivative of coefficient matrix.
                0186        DO bj=myByLo(myThid),myByHi(myThid)
                0187         DO bi=myBxLo(myThid),myBxHi(myThid)
                0188          DO j=1,sNy
                0189           DO i=1,sNx
                0190            aC2d_ad(i,j,bi,bj)   = aC2d_ad(i,j,bi,bj)
                0191      &          - cg2d_b_ad(i,j,bi,bj) * cg2d_x(i,j,bi,bj)
                0192      &          * recip_cg2dNorm
                0193            aW2d_ad(i,  j,bi,bj) = aW2d_ad(i,j,bi,bj)
                0194      &          - cg2d_b_ad(i,j,bi,bj) * cg2d_x(i-1,j,bi,bj)
                0195      &          * recip_cg2dNorm
                0196            aW2d_ad(i+1,j,bi,bj) = aW2d_ad(i+1,j,bi,bj)
                0197      &          - cg2d_b_ad(i,j,bi,bj) * cg2d_x(i+1,j,bi,bj)
                0198      &          * recip_cg2dNorm
                0199            aS2d_ad(i,j,  bi,bj) = aS2d_ad(i,j,bi,bj)
                0200      &          - cg2d_b_ad(i,j,bi,bj) * cg2d_x(i,j-1,bi,bj)
                0201      &          * recip_cg2dNorm
                0202            aS2d_ad(i,j+1,bi,bj) = aS2d_ad(i,j+1,bi,bj)
                0203      &          - cg2d_b_ad(i,j,bi,bj) * cg2d_x(i,j+1,bi,bj)
                0204      &          * recip_cg2dNorm
                0205 C     The sensitivities of preconditioner are something like
                0206 C          pC_ad = pC_ad + cg2d_z_ad * cg2d_r
                0207 C     and similar (cg2d_z = preconditioner * cg2d_r). That is, when the
                0208 C     solver has converged, cg2d_r -> 0 and pC_ad -> 0, unless cg2d_z_ad
                0209 C     is large, which would mean an instability. Therefore we set the
                0210 C     preconditioner adjoint to zero. This also reflects the fact, that
                0211 C     in an ideal world, the preconditioner does not have any effect on
                0212 C     the solution, but only on the rate of convergence.
                0213            pW_ad(i,j,bi,bj) = 0. _d 0
                0214            pS_ad(i,j,bi,bj) = 0. _d 0
                0215            pC_ad(i,j,bi,bj) = 0. _d 0
                0216           ENDDO
                0217          ENDDO
                0218         ENDDO
                0219        ENDDO
                0220       ELSE
                0221 C     Initialise AD fields to zero for an adjoint neglecting the
                0222 C     contributions from the matrix coefficients (default for now).
                0223        DO bj=myByLo(myThid),myByHi(myThid)
                0224         DO bi=myBxLo(myThid),myBxHi(myThid)
                0225          DO j=1-OLy,sNy+OLy
                0226           DO i=1-OLx,sNx+OLx
                0227           aW2d_ad(i,j,bi,bj) = 0. _d 0
                0228           aS2d_ad(i,j,bi,bj) = 0. _d 0
                0229           aC2d_ad(i,j,bi,bj) = 0. _d 0
                0230           pW_ad  (i,j,bi,bj) = 0. _d 0
                0231           pS_ad  (i,j,bi,bj) = 0. _d 0
                0232           pC_ad  (i,j,bi,bj) = 0. _d 0
                0233           ENDDO
                0234          ENDDO
                0235         ENDDO
                0236        ENDDO
                0237       ENDIF
                0238 #endif /* NONLIN_FRSURF or ALLOW_DEPTH_CONTROL */
                0239 #endif /* ALLOW_AUTODIFF_TAMC */
                0240 
                0241       RETURN
                0242       END
                0243 
                0244 C--------1---------2---------3---------4---------5---------6---------7--
                0245 
                0246 CBOP
                0247 C     !ROUTINE: CG2D_STORE
                0248 C     !INTERFACE:
                0249       SUBROUTINE CG2D_STORE(
                0250      U                      cg2d_x,
                0251      I                      doStore,
                0252      I                      myThid )
                0253 C     !DESCRIPTION: \bv
                0254 C     *==========================================================*
                0255 C     | SUBROUTINE CG2D_STORE
                0256 C     | o Hand written tape routine
                0257 C     *==========================================================*
                0258 C     | This routine is called from solve_for_pressure to
                0259 C     | to force TAF to store cg2d_x
                0260 C     *==========================================================*
                0261 C     \ev
                0262 
                0263 C     !USES:
                0264       IMPLICIT NONE
                0265 C     === Global data ===
                0266 #include "SIZE.h"
                0267 #include "EEPARAMS.h"
                0268 #include "PARAMS.h"
                0269 #include "AUTODIFF_PARAMS.h"
                0270 #ifdef ALLOW_AUTODIFF_TAMC
                0271 # include "tamc.h"
                0272 #endif
                0273 
                0274 C     !INPUT/OUTPUT PARAMETERS:
                0275 C     === Routine arguments ===
                0276 C     cg2d_x    :: The solution to be stored and restored
                0277 C     doStore   :: Store or restore based on the logical
                0278 C     myThid    :: Thread on which I am working.
                0279       _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0280       LOGICAL doStore
                0281       INTEGER myThid
                0282 
                0283 #ifdef ALLOW_AUTODIFF_TAMC
                0284 #if ( defined NONLIN_FRSURF || defined ALLOW_DEPTH_CONTROL )
                0285       _RL cg2d_tape(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,nchklev_1)
                0286       COMMON /CG2D_TAPE_RL/ cg2d_tape
                0287 
                0288 C     !LOCAL VARIABLES:
                0289 C     === Local variables ===
                0290 C     bi,bj  :: tile indices
                0291 C     i,j    :: Loop counters
                0292       INTEGER i,j,bi,bj
                0293 CEOP
                0294 
                0295 C--   mark output of cg2d as being called from here
                0296       IF ( debugLevel .GE. debLevZero ) THEN
                0297        _BEGIN_MASTER( myThid )
44d3986245 Jean*0298        WRITE(standardMessageUnit,'(A,L1)')
aecc8b0f47 Mart*0299      &      ' Calling CG2D_STORE with doStore = ', doStore
44d3986245 Jean*0300        WRITE(standardMessageUnit,'(A,I6)')
                0301      &      ' Calling CG2D_STORE with ikey_dynamics=', ikey_dynamics
aecc8b0f47 Mart*0302        _END_MASTER( myThid )
                0303       ENDIF
                0304 
                0305       IF ( doStore ) THEN
                0306 C     Save input cg2d_x to local common block
                0307        DO bj = myByLo(myThid), myByHi(myThid)
                0308         DO bi = myBxLo(myThid), myBxHi(myThid)
                0309          DO j=1-OLy,sNy+OLy
                0310           DO i=1-OLx,sNx+OLx
                0311            cg2d_tape(i,j,bi,bj,ikey_dynamics) = cg2d_x(i,j,bi,bj)
                0312           ENDDO
                0313          ENDDO
                0314         ENDDO
                0315        ENDDO
                0316       ELSE
                0317 C     Restore output cg2d_x from local common block
                0318        DO bj = myByLo(myThid), myByHi(myThid)
                0319         DO bi = myBxLo(myThid), myBxHi(myThid)
                0320          DO j=1-OLy,sNy+OLy
                0321           DO i=1-OLx,sNx+OLx
                0322            cg2d_x(i,j,bi,bj) = cg2d_tape(i,j,bi,bj,ikey_dynamics)
                0323           ENDDO
                0324          ENDDO
                0325         ENDDO
                0326        ENDDO
                0327       ENDIF
                0328 #endif /* NONLIN_FRSURF or ALLOW_DEPTH_CONTROL */
                0329 #endif /* ALLOW_AUTODIFF_TAMC */
                0330 
                0331       RETURN
                0332       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!