--- MITgcm/pkg/grdchk/grdchk_main.F 2003/10/23 04:41:40 1.10 +++ MITgcm/pkg/grdchk/grdchk_main.F 2011/09/26 12:56:52 1.33 @@ -1,9 +1,7 @@ -C -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/grdchk/grdchk_main.F,v 1.10 2003/10/23 04:41:40 edhill Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/grdchk/grdchk_main.F,v 1.33 2011/09/26 12:56:52 jmc Exp $ C $Name: $ -#include "AD_CONFIG.h" -#include "CPP_OPTIONS.h" +#include "GRDCHK_OPTIONS.h" CBOI C @@ -33,7 +31,7 @@ c | c |-- grdchk_getxx - get control vector component from file c | perturb it and write back to file -c |-- grdchk_getadxx - get gradient component calculated +c |-- grdchk_getadxx - get gradient component calculated c | via adjoint c |-- the_main_loop - forward run and cost evaluation c | with perturbed control vector element @@ -66,7 +64,7 @@ c - added tangent linear gradient checks c - fixes for multiproc. gradient checks c - added more control variables -c +c c ================================================================== c SUBROUTINE grdchk_main c ================================================================== @@ -81,19 +79,20 @@ #include "PARAMS.h" #include "grdchk.h" #include "cost.h" +#include "ctrl.h" #ifdef ALLOW_TANGENTLINEAR_RUN #include "g_cost.h" #endif C !INPUT/OUTPUT PARAMETERS: c == routine arguments == - integer mythid + integer mythid -#ifdef ALLOW_GRADIENT_CHECK +#ifdef ALLOW_GRDCHK C !LOCAL VARIABLES: c == local variables == integer myiter - _RL mytime + _RL mytime integer bi, itlo, ithi integer bj, jtlo, jthi integer i, imin, imax @@ -109,6 +108,7 @@ integer obcspos integer itilepos integer jtilepos + integer icglo integer itest integer ierr integer ierr_grdchk @@ -128,6 +128,8 @@ _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy) _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy) + CHARACTER*(MAX_LEN_MBUF) msgBuf + c == end of interface == CEOP @@ -141,20 +143,26 @@ imin = 1 imax = snx + print *, 'ph-check entering grdchk_main ' + c-- initialise variables call grdchk_init( mythid ) -c-- Compute the adjoint models' gradients. +c-- Compute the adjoint model gradients. c-- Compute the unperturbed cost function. cph Gradient via adjoint has already been computed, cph and so has unperturbed cost function, cph assuming all xx_ fields are initialised to zero. + ierr = 0 ierr_grdchk = 0 -cphadmtlm( + adxxmemo = 0. + ftlxxmemo = 0. +#ifdef ALLOW_ADMTLM + fcref = objf_state_final(idep,jdep,1,1,1) +#else fcref = fc -cphadmtlm fcref = objf_state_final(45,4,1,1) -cphadmtlm) +#endif print *, 'ph-check fcref = ', fcref @@ -178,39 +186,52 @@ grdchk_epsfac = 1. _d 0 end if - print *, 'ph-grd 3 -------------------------------' - print ('(2a)'), - & ' ph-grd 3 proc # i j k fc ref', - & ' fc + eps fc - eps' + WRITE(standardmessageunit,'(A)') + & 'grad-res -------------------------------' + WRITE(standardmessageunit,'(2a)') + & ' grad-res proc # i j k bi bj iobc', + & ' fc ref fc + eps fc - eps' #ifdef ALLOW_TANGENTLINEAR_RUN - print ('(2a)'), - & ' ph-grd 3 proc # i j k tlm grad', - & ' fd grad 1 - fd/tlm' + WRITE(standardmessageunit,'(2a)') + & ' grad-res proc # i j k bi bj iobc', + & ' tlm grad fd grad 1 - fd/tlm' #else - print ('(2a)'), - & ' ph-grd 3 proc # i j k adj grad', - & ' fd grad 1 - fd/adj' + WRITE(standardmessageunit,'(2a)') + & ' grad-res proc # i j k bi bj iobc', + & ' adj grad fd grad 1 - fd/adj' #endif c-- Compute the finite difference approximations. c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep c-- gradient checks. + if ( nbeg .EQ. 0 ) + & call grdchk_get_position( mythid ) + do icomp = nbeg, nend, nstep ichknum = (icomp - nbeg)/nstep + 1 + adxxmemo = 0. +cph( +cph-print print *, 'ph-grd _main: nbeg, icomp, ichknum ', +cph-print & nbeg, icomp, ichknum +cph) if (ichknum .le. maxgrdchecks ) then c-- Determine the location of icomp on the grid. if ( myProcId .EQ. grdchkwhichproc ) then - call grdchk_loc( icomp, ichknum, + call grdchk_loc( icomp, ichknum, & icvrec, itile, jtile, layer, obcspos, - & itilepos, jtilepos, itest, ierr, + & itilepos, jtilepos, icglo, itest, ierr, & mythid ) +cph( +cph-print print *, 'ph-grd ----- back from loc -----', +cph-print & icvrec, itilepos, jtilepos, layer, obcspos +cph) endif _BARRIER - + c****************************************************** c-- (A): get gradient component calculated via adjoint c****************************************************** @@ -223,7 +244,9 @@ & itilepos, jtilepos, & adxxmemo, mythid ) endif - _BARRIER +C-- Add a global-sum call so that all proc will get the adjoint gradient + _GLOBAL_SUM_RL( adxxmemo, myThid ) +c _BARRIER #ifdef ALLOW_TANGENTLINEAR_RUN c****************************************************** @@ -247,20 +270,26 @@ c-- 2. perform tangent linear run mytime = starttime myiter = niter0 -cphadmtlm( +#ifdef ALLOW_ADMTLM + do k=1,4*Nr+1 + do j=1,sny + do i=1,snx + g_objf_state_final(i,j,1,1,k) = 0. + enddo + enddo + enddo +#else g_fc = 0. -cphadmtlm do j=1,sny -cphadmtlm do i=1,snx -cphadmtlm g_objf_state_final(i,j,1,1) = 0. -cphadmtlm enddo -cphadmtlm enddo -cphadmtlm) +#endif + call g_the_main_loop( mytime, myiter, mythid ) -cphadmtlm( - ftlxxmemo = g_fc -cphadmtlm ftlxxmemo = g_objf_state_final(45,4,1,1) -cphadmtlm) _BARRIER +#ifdef ALLOW_ADMTLM + ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1) +#else + ftlxxmemo = g_fc +#endif + c-- c-- 3. reset control vector if ( myProcId .EQ. grdchkwhichproc .AND. @@ -292,17 +321,19 @@ & mythid ) endif _BARRIER - + c-- forward run with perturbed control vector mytime = starttime myiter = niter0 call the_main_loop( mytime, myiter, mythid ) -cphadmtlm( +#ifdef ALLOW_ADMTLM + fcpertplus = objf_state_final(idep,jdep,1,1,1) +#else fcpertplus = fc -cphadmtlm fcpertplus = objf_state_final(45,4,1,1) -cphadmtlm) +#endif + print *, 'ph-check fcpertplus = ', fcpertplus _BARRIER - + c-- Reset control vector. if ( myProcId .EQ. grdchkwhichproc .AND. & ierr .EQ. 0 ) then @@ -314,6 +345,7 @@ _BARRIER fcpertminus = fcref + print *, 'ph-check fcpertminus = ', fcpertminus if ( useCentralDiff ) then @@ -330,14 +362,18 @@ & mythid ) endif _BARRIER - + c-- forward run with perturbed control vector mytime = starttime myiter = niter0 call the_main_loop( mytime, myiter, mythid ) _BARRIER +#ifdef ALLOW_ADMTLM + fcpertminus = objf_state_final(idep,jdep,1,1,1) +#else fcpertminus = fc - +#endif + c-- Reset control vector. if ( myProcId .EQ. grdchkwhichproc .AND. & ierr .EQ. 0 ) then @@ -355,28 +391,28 @@ c-- (D): calculate relative differences between gradients c****************************************************** + if ( grdchk_eps .eq. 0. ) then + gfd = (fcpertplus-fcpertminus) + else + gfd = (fcpertplus-fcpertminus) + & /(grdchk_epsfac*grdchk_eps) + endif + + if ( adxxmemo .eq. 0. ) then + ratio_ad = abs( adxxmemo - gfd ) + else + ratio_ad = 1. - gfd/adxxmemo + endif + + if ( ftlxxmemo .eq. 0. ) then + ratio_ftl = abs( ftlxxmemo - gfd ) + else + ratio_ftl = 1. - gfd/ftlxxmemo + endif + if ( myProcId .EQ. grdchkwhichproc .AND. & ierr .EQ. 0 ) then - if ( grdchk_eps .eq. 0. ) then - gfd = (fcpertplus-fcpertminus) - else - gfd = (fcpertplus-fcpertminus) - & /(grdchk_epsfac*grdchk_eps) - endif - - if ( adxxmemo .eq. 0. ) then - ratio_ad = abs( adxxmemo - gfd ) - else - ratio_ad = 1. - gfd/adxxmemo - endif - - if ( ftlxxmemo .eq. 0. ) then - ratio_ftl = abs( ftlxxmemo - gfd ) - else - ratio_ftl = 1. - gfd/ftlxxmemo - endif - tmpplot1(itilepos,jtilepos,layer,itile,jtile) & = gfd tmpplot2(itilepos,jtilepos,layer,itile,jtile) @@ -406,24 +442,51 @@ ichkmem ( ichknum ) = ichknum itestmem ( ichknum ) = itest ierrmem ( ichknum ) = ierr - - print *, 'ph-grd 3 -------------------------------' - print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ', - & myprocid,ichknum,itilepos,jtilepos,layer, + icglomem ( ichknum ) = icglo + + WRITE(standardmessageunit,'(A)') + & 'grad-res -------------------------------' + WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)') + & ' grad-res ',myprocid,ichknum,itilepos,jtilepos, + & layer,itile,jtile,obcspos, & fcref, fcpertplus, fcpertminus #ifdef ALLOW_TANGENTLINEAR_RUN - print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ', - & myprocid,ichknum,ichkmem(ichknum), + WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)') + & ' grad-res ',myprocid,ichknum,ichkmem(ichknum), & icompmem(ichknum),itestmem(ichknum), + & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum), & ftlxxmemo, gfd, ratio_ftl #else - print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ', - & myprocid,ichknum,ichkmem(ichknum), + WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)') + & ' grad-res ',myprocid,ichknum,ichkmem(ichknum), & icompmem(ichknum),itestmem(ichknum), + & bimem(ichknum),bjmem(ichknum),obcspos, & adxxmemo, gfd, ratio_ad #endif - endif +#ifdef ALLOW_TANGENTLINEAR_RUN + WRITE(msgBuf,'(A34,1PE24.14)') + & ' TLM precision_derivative_cost =', fcref + CALL PRINT_MESSAGE + & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid) + WRITE(msgBuf,'(A34,1PE24.14)') + & ' TLM precision_derivative_grad =', ftlxxmemo + CALL PRINT_MESSAGE + & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid) +#else + WRITE(msgBuf,'(A30,1PE22.14)') + & ' ADM ref_cost_function =', fcref + CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, + & SQUEEZE_RIGHT, myThid ) + WRITE(msgBuf,'(A30,1PE22.14)') + & ' ADM adjoint_gradient =', adxxmemo + CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, + & SQUEEZE_RIGHT, myThid ) + WRITE(msgBuf,'(A30,1PE22.14)') + & ' ADM finite-diff_grad =', gfd + CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, + & SQUEEZE_RIGHT, myThid ) +#endif print *, 'ph-grd ierr ---------------------------' print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp, @@ -435,18 +498,18 @@ else ierr_grdchk = -1 -c-- end of if ( ichknum ... +c-- end of if ( ichknum ... endif -c-- end of do icomp = ... +c-- end of do icomp = ... enddo if ( myProcId .EQ. grdchkwhichproc ) then - CALL WRITE_REC_XYZ_RL( + CALL WRITE_REC_XYZ_RL( & 'grd_findiff' , tmpplot1, 1, 0, myThid) - CALL WRITE_REC_XYZ_RL( + CALL WRITE_REC_XYZ_RL( & 'grd_ratio_ad' , tmpplot2, 1, 0, myThid) - CALL WRITE_REC_XYZ_RL( + CALL WRITE_REC_XYZ_RL( & 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid) endif @@ -456,6 +519,7 @@ c-- Print the results of the gradient check. call grdchk_print( ichknum, ierr_grdchk, mythid ) -#endif /* ALLOW_GRADIENT_CHECK */ +#endif /* ALLOW_GRDCHK */ + return end