/[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.2 by heimbach, Thu Feb 13 23:36:18 2003 UTC revision 1.3.6.3 by heimbach, Tue Mar 4 01:25:20 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 106  c     == local variables == Line 108  c     == local variables ==
108        _RL     gfd        _RL     gfd
109        _RL     fcref        _RL     fcref
110        _RL     fcpertplus, fcpertminus        _RL     fcpertplus, fcpertminus
111        _RL     ratio        _RL     ratio_ad
112          _RL     ratio_ftl
113        _RL     xxmemo_ref        _RL     xxmemo_ref
114        _RL     xxmemo_pert        _RL     xxmemo_pert
115        _RL     adxxmemo        _RL     adxxmemo
116          _RL     ftlxxmemo
117          _RL     localEps
118          _RL     grdchk_epsfac
119    
 cph(  
120        _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _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)        _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
122  cph)        _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
123    
124  Cml(  #ifdef ALLOW_TANGENTLINEAR_RUN
125        _RL     grdchk_epsfac        _RL     g_fc
126  Cml)        common /g_cost_r/ g_fc
127    #endif
128    
129  c     == end of interface ==  c     == end of interface ==
130  CEOP  CEOP
# Line 145  cph   assuming all xx_ fields are initia Line 151  cph   assuming all xx_ fields are initia
151        fcref = fc        fcref = fc
152        ierr_grdchk = 0        ierr_grdchk = 0
153    
 cph(  
154        do bj = jtlo, jthi        do bj = jtlo, jthi
155           do bi = itlo, ithi           do bi = itlo, ithi
156              do k = 1, nr              do k = 1, nr
# Line 153  cph( Line 158  cph(
158                    do i = imin, imax                    do i = imin, imax
159                       tmpplot1(i,j,k,bi,bj) = 0. _d 0                       tmpplot1(i,j,k,bi,bj) = 0. _d 0
160                       tmpplot2(i,j,k,bi,bj) = 0. _d 0                       tmpplot2(i,j,k,bi,bj) = 0. _d 0
161                         tmpplot3(i,j,k,bi,bj) = 0. _d 0
162                    end do                    end do
163                 end do                 end do
164              end do              end do
165           end do           end do
166        end do        end do
 cph)  
167    
168        if ( useCentralDiff ) then        if ( useCentralDiff ) then
169           grdchk_epsfac = 2. _d 0           grdchk_epsfac = 2. _d 0
# Line 170  cph) Line 175  cph)
175        print ('(2a)'),        print ('(2a)'),
176       &     ' ph-grd 3  proc    #    i    j    k            fc ref',       &     ' ph-grd 3  proc    #    i    j    k            fc ref',
177       &     '        fc + eps        fc - eps'       &     '        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)'),        print ('(2a)'),
184       &     ' ph-grd 3  proc    #    i    j    k          adj grad',       &     ' ph-grd 3  proc    #    i    j    k          adj grad',
185       &     '         fd grad      1 - fd/adj'       &     '         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
   
                if (ichknum .le. maxgrdchecks ) then  
197    
198  c--         Determine the location of icomp on the grid.  c--         Determine the location of icomp on the grid.
199                    call grdchk_loc( icomp, ichknum,              if ( myProcId .EQ. grdchkwhichproc ) then
200       &                 icvrec, itile, jtile, layer,                 call grdchk_loc( icomp, ichknum,
201       &                 itilepos, jtilepos, itest, ierr,       &              icvrec, itile, jtile, layer,
202       &                 mythid )       &              itilepos, jtilepos, itest, ierr,
203         &              mythid )
204                endif
205                _BARRIER
206                                
207                    if ( ierr .eq. 0 ) then  c******************************************************
208    c--   (A): get gradient component calculated via adjoint
209  c--   positive perturbation  c******************************************************
                      grdchk_eps = abs(grdchk_eps)  
210    
211  c--   get gradient component calculated via adjoint  c--   get gradient component calculated via adjoint
212                       call grdchk_getadxx( icvrec,              if ( myProcId .EQ. grdchkwhichproc .AND.
213       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
214       &                    itilepos, jtilepos,                 call grdchk_getadxx( icvrec,
215       &                    adxxmemo, mythid )       &              itile, jtile, layer,
216                       _BARRIER       &              itilepos, jtilepos,
217         &              adxxmemo, mythid )
218                endif
219                _BARRIER
220    
221    #ifdef ALLOW_TANGENTLINEAR_RUN
222    c******************************************************
223    c--   (B): Get gradient component g_fc from tangent linear run:
224    c******************************************************
225    c--
226    c--   1. perturb control vector component: xx(i)=1.
227    
228                if ( myProcId .EQ. grdchkwhichproc .AND.
229         &           ierr .EQ. 0 ) then
230                   localEps = 1. _d 0
231                   call grdchk_getxx( icvrec, TANGENT_SIMULATION,
232         &              itile, jtile, layer,
233         &              itilepos, jtilepos,
234         &              xxmemo_ref, xxmemo_pert, localEps,
235         &              mythid )
236                endif
237                _BARRIER
238    
239    c--
240    c--   2. perform tangent linear run
241                mytime = starttime
242                myiter = niter0
243                g_fc = 0.
244                call g_the_main_loop( mytime, myiter, mythid )
245                ftlxxmemo = g_fc
246                _BARRIER
247    c--
248    c--   3. reset control vector
249                if ( myProcId .EQ. grdchkwhichproc .AND.
250         &           ierr .EQ. 0 ) then
251                   call grdchk_setxx( icvrec, TANGENT_SIMULATION,
252         &              itile, jtile, layer,
253         &              itilepos, jtilepos,
254         &              xxmemo_ref, mythid )
255                endif
256                _BARRIER
257    
258    #endif /* ALLOW_TANGENTLINEAR_RUN */
259    
260    
261    c******************************************************
262    c--   (C): Get gradient via finite difference perturbation
263    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                       call grdchk_getxx( icvrec,  c--   positive perturbation
268       &                    itile, jtile, layer,              localEps = abs(grdchk_eps)
269       &                    itilepos, jtilepos,              if ( myProcId .EQ. grdchkwhichproc .AND.
270       &                    xxmemo_ref, xxmemo_pert, mythid )       &           ierr .EQ. 0 ) then
271                       _BARRIER                 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
272         &              itile, jtile, layer,
273         &              itilepos, jtilepos,
274         &              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,              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
   
 c--   repeat the proceedure for a negative perturbation  
                         grdchk_eps = - abs(grdchk_eps)  
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                          call grdchk_getxx( icvrec,  c--   repeat the proceedure for a negative perturbation
303       &                    itile, jtile, layer,                 if ( myProcId .EQ. grdchkwhichproc .AND.
304       &                    itilepos, jtilepos,       &           ierr .EQ. 0 ) then
305       &                    xxmemo_ref, xxmemo_pert, mythid )                    localEps = - abs(grdchk_eps)
306                          _BARRIER                    call grdchk_getxx( icvrec, FORWARD_SIMULATION,
307         &                 itile, jtile, layer,
308         &                 itilepos, jtilepos,
309         &                 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,                 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     reset grdchk_eps to a positive value  c-- end of if useCentralDiff ...
332                          grdchk_eps = abs(grdchk_eps)              end if
333    
334                       end if  c******************************************************
335  c--  c--   (D): calculate relative differences between gradients
336                       if ( grdchk_eps .eq. 0. ) then  c******************************************************
337                          gfd = (fcpertplus-fcpertminus)  
338                       else              if ( myProcId .EQ. grdchkwhichproc .AND.
339                          gfd = (fcpertplus-fcpertminus)       &           ierr .EQ. 0 ) then
      &                       /(grdchk_epsfac*grdchk_eps)  
                      endif  
                       
                      if ( adxxmemo .eq. 0. ) then  
                         ratio = abs( adxxmemo - gfd )  
                      else  
                         ratio = 1. - gfd/adxxmemo  
                      endif  
                     
 cph(  
                      tmpplot1(itilepos,jtilepos,layer,itile,jtile) =  
      &                    gfd  
                      tmpplot2(itilepos,jtilepos,layer,itile,jtile) =  
      &                    ratio  
 cph)  
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  
                      ratiomem  ( ichknum ) = ratio  
   
                      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,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',  
      &                    myprocid,ichknum,itilepos,jtilepos,layer,  
      &                    fcref, fcpertplus, fcpertminus  
                      print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',  
      &                    myprocid,ichknum,ichkmem(ichknum),  
      &                    icompmem(ichknum),itestmem(ichknum),  
      &                    adxxmemo, gfd, ratio  
 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  cph(        if ( myProcId .EQ. grdchkwhichproc ) then
424        CALL WRITE_REC_XYZ_RL( 'grd_findiff' , tmpplot1, 1, 0, myThid)           CALL WRITE_REC_XYZ_RL(
425        CALL WRITE_REC_XYZ_RL( 'grd_ratio' , tmpplot2, 1, 0, myThid)       &        'grd_findiff'   , tmpplot1, 1, 0, myThid)
426  cph)           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.3.6.2  
changed lines
  Added in v.1.3.6.3
  ViewVC Help
Powered by ViewVC 1.1.22