/[MITgcm]/MITgcm/pkg/grdchk/grdchk_main.F
ViewVC logotype

Diff of /MITgcm/pkg/grdchk/grdchk_main.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3.6.4 by heimbach, Fri Mar 7 04:01:59 2003 UTC revision 1.6 by heimbach, Mon Sep 16 18:11:58 2002 UTC
# Line 59  c     changed: mlosch@ocean.mit.edu: 09- Line 59  c     changed: mlosch@ocean.mit.edu: 09-
59  c              - added centered difference vs. 1-sided difference option  c              - added centered difference vs. 1-sided difference option
60  c              - improved output format for readability  c              - improved output format for readability
61  c              - added control variable hFacC  c              - added control variable hFacC
 c              heimbach@mit.edu 24-Feb-2003  
 c              - added tangent linear gradient checks  
 c              - fixes for multiproc. gradient checks  
 c              - added more control variables  
62  c      c    
63  c     ==================================================================  c     ==================================================================
64  c     SUBROUTINE grdchk_main  c     SUBROUTINE grdchk_main
# Line 94  c     == local variables == Line 90  c     == local variables ==
90        integer j,  jmin, jmax        integer j,  jmin, jmax
91        integer k        integer k
92    
93          integer jprocs
94          integer proc_grdchk
95        integer icomp        integer icomp
96        integer ichknum        integer ichknum
97        integer icvrec        integer icvrec
# Line 117  c     == local variables == Line 115  c     == local variables ==
115        _RL     localEps        _RL     localEps
116        _RL     grdchk_epsfac        _RL     grdchk_epsfac
117    
       _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)  
       _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)  
   
118  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
119        _RL     g_fc        _RL     g_fc
120        common /g_cost_r/ g_fc        common /g_cost_r/ g_fc
121  #endif  #endif
122    
123          _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
124          _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
125          _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
126    
127  c     == end of interface ==  c     == end of interface ==
128  CEOP  CEOP
129    
# Line 144  c--   initialise variables Line 142  c--   initialise variables
142    
143  c--   Compute the adjoint models' gradients.  c--   Compute the adjoint models' gradients.
144  c--   Compute the unperturbed cost function.  c--   Compute the unperturbed cost function.
145  cph   Gradient via adjoint has already been computed,  c--   Gradient via adjoint has already been computed,
146  cph   and so has unperturbed cost function,  c--   and so has unperturbed cost function,
147  cph   assuming all xx_ fields are initialised to zero.  c--   assuming all xx_ fields are initialised to zero.
148    
149        fcref = fc        fcref = fc
150        ierr_grdchk = 0        ierr_grdchk = 0
# Line 171  cph   assuming all xx_ fields are initia Line 169  cph   assuming all xx_ fields are initia
169           grdchk_epsfac = 1. _d 0           grdchk_epsfac = 1. _d 0
170        end if        end if
171    
       print *, 'ph-grd 3 -------------------------------'  
       print ('(2a)'),  
      &     ' ph-grd 3  proc    #    i    j    k            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'  
 #else  
       print ('(2a)'),  
      &     ' ph-grd 3  proc    #    i    j    k          adj grad',  
      &     '         fd grad      1 - fd/adj'  
 #endif  
   
172  c--   Compute the finite difference approximations.  c--   Compute the finite difference approximations.
173  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
174  c--   gradient checks.  c--   gradient checks.
175          do jprocs = 1,numberOfProcs
176             proc_grdchk = jprocs - 1
177    
178        do icomp = nbeg, nend, nstep           if ( myProcId .eq. proc_grdchk ) then
179    
180           ichknum = (icomp - nbeg)/nstep + 1              do icomp = nbeg, nend, nstep
181    
182           if (ichknum .le. maxgrdchecks ) then                 ichknum = (icomp - nbeg)/nstep + 1
183    
184  c--         Determine the location of icomp on the grid.                 if (ichknum .le. maxgrdchecks ) then
185              if ( myProcId .EQ. grdchkwhichproc ) then  
186                 call grdchk_loc( icomp, ichknum,  c--   Determine the location of icomp on the grid.
187       &              icvrec, itile, jtile, layer,                    call grdchk_loc( icomp, ichknum,
188       &              itilepos, jtilepos, itest, ierr,       &                 icvrec, itile, jtile, layer,
189       &              mythid )       &                 itilepos, jtilepos, itest, ierr,
190              endif       &                 mythid )
             _BARRIER  
191                                
192                      if ( ierr .eq. 0 ) then
193    
194  c******************************************************  c******************************************************
195  c--   (A): get gradient component calculated via adjoint  c--   (A): get gradient component calculated via adjoint
196  c******************************************************  c******************************************************
197                         call grdchk_getadxx( icvrec,
198  c--   get gradient component calculated via adjoint       &                    itile, jtile, layer,
199              if ( myProcId .EQ. grdchkwhichproc .AND.       &                    itilepos, jtilepos,
200       &           ierr .EQ. 0 ) then       &                    adxxmemo, mythid )
201                 call grdchk_getadxx( icvrec,                       _BARRIER
      &              itile, jtile, layer,  
      &              itilepos, jtilepos,  
      &              adxxmemo, mythid )  
             endif  
             _BARRIER  
202    
203  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
204  c******************************************************  c******************************************************
# Line 225  c*************************************** Line 207  c***************************************
207  c--  c--
208  c--   1. perturb control vector component: xx(i)=1.  c--   1. perturb control vector component: xx(i)=1.
209    
210              if ( myProcId .EQ. grdchkwhichproc .AND.                       localEps = 1.
211       &           ierr .EQ. 0 ) then                       call grdchk_getxx( icvrec, TANGENT_SIMULATION,
212                 localEps = 1. _d 0       &                    itile, jtile, layer,
213                 call grdchk_getxx( icvrec, TANGENT_SIMULATION,       &                    itilepos, jtilepos,
214       &              itile, jtile, layer,       &                    xxmemo_ref, xxmemo_pert, localEps,
215       &              itilepos, jtilepos,       &                    mythid )
216       &              xxmemo_ref, xxmemo_pert, localEps,                       _BARRIER
      &              mythid )  
             endif  
             _BARRIER  
217    
218  c--  c--
219  c--   2. perform tangent linear run  c--   2. perform tangent linear run
220              mytime = starttime                       mytime = starttime
221              myiter = niter0                       myiter = niter0
222              g_fc = 0.                       g_fc = 0.
223              call g_the_main_loop( mytime, myiter, mythid )                       call g_the_main_loop( mytime, myiter, mythid )
224              ftlxxmemo = g_fc                       ftlxxmemo = g_fc
             _BARRIER  
225  c--  c--
226  c--   3. reset control vector  c--   3. reset control vector
227              if ( myProcId .EQ. grdchkwhichproc .AND.                       call grdchk_setxx( icvrec, TANGENT_SIMULATION,
228       &           ierr .EQ. 0 ) then       &                    itile, jtile, layer,
229                 call grdchk_setxx( icvrec, TANGENT_SIMULATION,       &                    itilepos, jtilepos,
230       &              itile, jtile, layer,       &                    xxmemo_ref, mythid )
231       &              itilepos, jtilepos,                       _BARRIER
232       &              xxmemo_ref, mythid )  #endif
             endif  
             _BARRIER  
   
 #endif /* ALLOW_TANGENTLINEAR_RUN */  
   
233    
234  c******************************************************  c******************************************************
235  c--   (C): Get gradient via finite difference perturbation  c--   (C): Get gradient via finite difference perturbation
236  c******************************************************  c******************************************************
237    
238  c--   get control vector component from file  c--   get control vector component from file
239  c--   perturb it and write back to file  c--   perturb it and write back to file:
240  c--   positive perturbation  c--   positive perturbation
241              localEps = abs(grdchk_eps)                       localEps = abs(grdchk_eps)
242              if ( myProcId .EQ. grdchkwhichproc .AND.                       call grdchk_getxx( icvrec, FORWARD_SIMULATION,
243       &           ierr .EQ. 0 ) then       &                    itile, jtile, layer,
244                 call grdchk_getxx( icvrec, FORWARD_SIMULATION,       &                    itilepos, jtilepos,
245       &              itile, jtile, layer,       &                    xxmemo_ref, xxmemo_pert, localEps,
246       &              itilepos, jtilepos,       &                    mythid )
247       &              xxmemo_ref, xxmemo_pert, localEps,                       _BARRIER
      &              mythid )  
             endif  
             _BARRIER  
248                                            
249  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
250              mytime = starttime                       mytime = starttime
251              myiter = niter0                       myiter = niter0
252              call the_main_loop( mytime, myiter, mythid )                       call the_main_loop( mytime, myiter, mythid )
253              fcpertplus = fc                       fcpertplus = fc
             _BARRIER  
254                                        
255  c--   Reset control vector.  c--   Reset control vector.
256              if ( myProcId .EQ. grdchkwhichproc .AND.                       call grdchk_setxx( icvrec, FORWARD_SIMULATION,
257       &           ierr .EQ. 0 ) then       &                    itile, jtile, layer,
258                 call grdchk_setxx( icvrec, FORWARD_SIMULATION,       &                    itilepos, jtilepos,
259       &              itile, jtile, layer,       &                    xxmemo_ref, mythid )
260       &              itilepos, jtilepos,                       _BARRIER
      &              xxmemo_ref, mythid )  
             endif  
             _BARRIER  
261    
262              fcpertminus = fcref                       fcpertminus = fcref
263    
264              if ( useCentralDiff ) then                       if ( useCentralDiff ) then
265    
266  c--   get control vector component from file  c--   get control vector component from file
267  c--   perturb it and write back to file  c--   perturb it and write back to file:
268  c--   repeat the proceedure for a negative perturbation  c--   repeat the proceedure for a negative perturbation
269                 if ( myProcId .EQ. grdchkwhichproc .AND.                          localEps = - abs(grdchk_eps)
270       &           ierr .EQ. 0 ) then                          call grdchk_getxx( icvrec, FORWARD_SIMULATION,
271                    localEps = - abs(grdchk_eps)       &                    itile, jtile, layer,
272                    call grdchk_getxx( icvrec, FORWARD_SIMULATION,       &                    itilepos, jtilepos,
273       &                 itile, jtile, layer,       &                    xxmemo_ref, xxmemo_pert, localEps,
274       &                 itilepos, jtilepos,       &                    mythid )
275       &                 xxmemo_ref, xxmemo_pert, localEps,                          _BARRIER
      &                 mythid )  
                endif  
                _BARRIER  
276                                            
277  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
278                 mytime = starttime                          mytime = starttime
279                 myiter = niter0                          myiter = niter0
280                 call the_main_loop( mytime, myiter, mythid )                          call the_main_loop( mytime, myiter, mythid )
281                 _BARRIER                          fcpertminus = fc
                fcpertminus = fc  
282                                        
283  c--   Reset control vector.  c--   Reset control vector.
284                 if ( myProcId .EQ. grdchkwhichproc .AND.                          call grdchk_setxx( icvrec, FORWARD_SIMULATION,
285       &           ierr .EQ. 0 ) then       &                    itile, jtile, layer,
286                    call grdchk_setxx( icvrec, FORWARD_SIMULATION,       &                    itilepos, jtilepos,
287       &                 itile, jtile, layer,       &                    xxmemo_ref, mythid )
288       &                 itilepos, jtilepos,                          _BARRIER
      &                 xxmemo_ref, mythid )  
                endif  
                _BARRIER  
   
 c-- end of if useCentralDiff ...  
             end if  
289    
290                         end if
291    c--
292  c******************************************************  c******************************************************
293  c--   (D): calculate relative differences between gradients  c--   (D): calculate relative differences between gradients
294  c******************************************************  c******************************************************
295    
296              if ( myProcId .EQ. grdchkwhichproc .AND.                       if ( grdchk_eps .eq. 0. ) then
297       &           ierr .EQ. 0 ) then                          gfd = (fcpertplus-fcpertminus)
298                         else
299                            gfd = (fcpertplus-fcpertminus)
300         &                       /(grdchk_epsfac*grdchk_eps)
301                         endif
302    
303                         if ( adxxmemo .eq. 0. ) then
304                            ratio_ad = abs( adxxmemo - gfd )
305                         else
306                            ratio_ad = 1. - gfd/adxxmemo
307                         endif
308                      
309                         if ( ftlxxmemo .eq. 0. ) then
310                            ratio_ftl = abs( ftlxxmemo - gfd )
311                         else
312                            ratio_ftl = 1. - gfd/ftlxxmemo
313                         endif
314                      
315                         tmpplot1(itilepos,jtilepos,layer,itile,jtile) =
316         &                    gfd
317                         tmpplot2(itilepos,jtilepos,layer,itile,jtile) =
318         &                    ratio_ad
319                         tmpplot3(itilepos,jtilepos,layer,itile,jtile) =
320         &                    ratio_ftl
321    
                if ( grdchk_eps .eq. 0. ) then  
                   gfd = (fcpertplus-fcpertminus)  
                else  
                   gfd = (fcpertplus-fcpertminus)  
      &                 /(grdchk_epsfac*grdchk_eps)  
                endif  
322                                            
323                 if ( adxxmemo .eq. 0. ) then                       fcrmem      ( ichknum ) = fcref
324                    ratio_ad = abs( adxxmemo - gfd )                       fcppmem     ( ichknum ) = fcpertplus
325                 else                       fcpmmem     ( ichknum ) = fcpertminus
326                    ratio_ad = 1. - gfd/adxxmemo                       xxmemref    ( ichknum ) = xxmemo_ref
327                 endif                       xxmempert   ( ichknum ) = xxmemo_pert
328                         gfdmem      ( ichknum ) = gfd
329                 if ( ftlxxmemo .eq. 0. ) then                       adxxmem     ( ichknum ) = adxxmemo
330                    ratio_ftl = abs( ftlxxmemo - gfd )                       ftlxxmem    ( ichknum ) = ftlxxmemo
331                         ratioadmem  ( ichknum ) = ratio_ad
332                         ratioftlmem ( ichknum ) = ratio_ftl
333    
334                         irecmem   ( ichknum ) = icvrec
335                         bimem     ( ichknum ) = itile
336                         bjmem     ( ichknum ) = jtile
337                         ilocmem   ( ichknum ) = itilepos
338                         jlocmem   ( ichknum ) = jtilepos
339                         klocmem   ( ichknum ) = layer
340                         icompmem  ( ichknum ) = icomp
341                         ichkmem   ( ichknum ) = ichknum
342                         itestmem  ( ichknum ) = itest
343                         ierrmem   ( ichknum ) = ierr
344                      
345    cph(
346                         print *, 'ph-grd 3 -------------------------------'
347                         print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
348         &                    ichknum,itilepos,jtilepos,layer,
349         &                    fcref, fcpertplus, fcpertminus
350                         print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
351         &                    ichknum,ichkmem(ichknum),
352         &                    icompmem(ichknum),itestmem(ichknum),
353         &                    adxxmemo, gfd, ratio_ad
354                         print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
355         &                    ichknum,ichkmem(ichknum),
356         &                    icompmem(ichknum),itestmem(ichknum),
357         &                    ftlxxmemo, gfd, ratio_ftl
358    cph)
359                      else
360    c
361                         print *, 'ph-grd 3 -------------------------------'
362                         print *, 'ph-grd 3 : ierr = ', ierr,
363         &                                 ', icomp = ', icomp
364                      endif
365                 else                 else
366                    ratio_ftl = 1. - gfd/ftlxxmemo                    ierr_grdchk = -1
367                 endif                 endif
368                                  
369                 tmpplot1(itilepos,jtilepos,layer,itile,jtile)              enddo
      &              = gfd  
                tmpplot2(itilepos,jtilepos,layer,itile,jtile)  
      &              = ratio_ad  
                tmpplot3(itilepos,jtilepos,layer,itile,jtile)  
      &              = ratio_ftl  
   
                fcrmem      ( ichknum ) = fcref  
                fcppmem     ( ichknum ) = fcpertplus  
                fcpmmem     ( ichknum ) = fcpertminus  
                xxmemref    ( ichknum ) = xxmemo_ref  
                xxmempert   ( ichknum ) = xxmemo_pert  
                gfdmem      ( ichknum ) = gfd  
                adxxmem     ( ichknum ) = adxxmemo  
                ftlxxmem    ( ichknum ) = ftlxxmemo  
                ratioadmem  ( ichknum ) = ratio_ad  
                ratioftlmem ( ichknum ) = ratio_ftl  
   
                irecmem   ( ichknum ) = icvrec  
                bimem     ( ichknum ) = itile  
                bjmem     ( ichknum ) = jtile  
                ilocmem   ( ichknum ) = itilepos  
                jlocmem   ( ichknum ) = jtilepos  
                klocmem   ( ichknum ) = layer  
                icompmem  ( ichknum ) = icomp  
                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,  
      &              fcref, fcpertplus, fcpertminus  
 #ifdef ALLOW_TANGENTLINEAR_RUN  
                print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',  
      &              myprocid,ichknum,ichkmem(ichknum),  
      &              icompmem(ichknum),itestmem(ichknum),  
      &              ftlxxmemo, gfd, ratio_ftl  
 #else  
                print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',  
      &              myprocid,ichknum,ichkmem(ichknum),  
      &              icompmem(ichknum),itestmem(ichknum),  
      &              adxxmemo, gfd, ratio_ad  
 #endif  
   
             endif  
   
             print *, 'ph-grd  ierr ---------------------------'  
             print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,  
      &           ', ichknum = ', ichknum  
   
             _BARRIER  
   
 c-- else of if ( ichknum ...  
          else  
             ierr_grdchk = -1  
   
 c-- end of if ( ichknum ...  
370           endif           endif
371    
372  c-- end of do icomp = ...  c--   Everyone has to wait for the component to be reset.
373        enddo           _BARRIER
374    
375        if ( myProcId .EQ. grdchkwhichproc ) then        enddo
          CALL WRITE_REC_XYZ_RL(  
      &        'grd_findiff'   , tmpplot1, 1, 0, myThid)  
          CALL WRITE_REC_XYZ_RL(  
      &        'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)  
          CALL WRITE_REC_XYZ_RL(  
      &        'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)  
       endif  
376    
377  c--   Everyone has to wait for the component to be reset.        CALL WRITE_REC_XYZ_RL( 'grd_findiff'   , tmpplot1, 1, 0, myThid)
378        _BARRIER        CALL WRITE_REC_XYZ_RL( 'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
379          CALL WRITE_REC_XYZ_RL( 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
380    
381  c--   Print the results of the gradient check.  c--   Print the results of the gradient check.
382        call grdchk_print( ichknum, ierr_grdchk, mythid )        call grdchk_print( ichknum, ierr_grdchk, mythid )
Legend:
Removed from v.1.3.6.4  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22