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