/[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.33 by jmc, Mon Sep 26 12:56:52 2011 UTC revision 1.40 by heimbach, Wed Aug 22 19:34:51 2012 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "GRDCHK_OPTIONS.h"  #include "GRDCHK_OPTIONS.h"
5    #include "AD_CONFIG.h"
6    
7  CBOI  CBOI
8  C  C
# Line 158  cph   assuming all xx_ fields are initia Line 159  cph   assuming all xx_ fields are initia
159        ierr_grdchk = 0        ierr_grdchk = 0
160        adxxmemo  = 0.        adxxmemo  = 0.
161        ftlxxmemo = 0.        ftlxxmemo = 0.
162  #ifdef ALLOW_ADMTLM  #if (defined  (ALLOW_ADMTLM))
163        fcref = objf_state_final(idep,jdep,1,1,1)        fcref = objf_state_final(idep,jdep,1,1,1)
164    #elif (defined (ALLOW_AUTODIFF_OPENAD))
165          fcref = fc%v
166  #else  #else
167        fcref = fc        fcref = fc
168  #endif  #endif
# Line 229  cph( Line 232  cph(
232  cph-print               print *, 'ph-grd ----- back from loc -----',  cph-print               print *, 'ph-grd ----- back from loc -----',
233  cph-print     &             icvrec, itilepos, jtilepos, layer, obcspos  cph-print     &             icvrec, itilepos, jtilepos, layer, obcspos
234  cph)  cph)
235                else
236                  icvrec   = 0
237                  itile    = 0
238                  jtile    = 0
239                  layer    = 0
240                  obcspos  = 0
241                  itilepos = 0
242                  jtilepos = 0
243                  icglo    = 0
244                  itest    = 0
245              endif              endif
246              _BARRIER  
247    c make sure that all procs have correct file records, so that useSingleCpuIO works
248          CALL GLOBAL_SUM_INT( icvrec , myThid )
249          CALL GLOBAL_SUM_INT( layer , myThid )
250          CALL GLOBAL_SUM_INT( ierr , myThid )
251    
252  c******************************************************  c******************************************************
253  c--   (A): get gradient component calculated via adjoint  c--   (A): get gradient component calculated via adjoint
254  c******************************************************  c******************************************************
255    
256  c--   get gradient component calculated via adjoint  c--   get gradient component calculated via adjoint
             if ( myProcId .EQ. grdchkwhichproc .AND.  
      &           ierr .EQ. 0 ) then  
257                 call grdchk_getadxx( icvrec,                 call grdchk_getadxx( icvrec,
258       &              itile, jtile, layer,       &              itile, jtile, layer,
259       &              itilepos, jtilepos,       &              itilepos, jtilepos,
260       &              adxxmemo, mythid )       &              adxxmemo, ierr, mythid )
             endif  
261  C--   Add a global-sum call so that all proc will get the adjoint gradient  C--   Add a global-sum call so that all proc will get the adjoint gradient
262              _GLOBAL_SUM_RL( adxxmemo, myThid )              _GLOBAL_SUM_RL( adxxmemo, myThid )
 c           _BARRIER  
263    
264  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
265  c******************************************************  c******************************************************
# Line 255  c*************************************** Line 268  c***************************************
268  c--  c--
269  c--   1. perturb control vector component: xx(i)=1.  c--   1. perturb control vector component: xx(i)=1.
270    
             if ( myProcId .EQ. grdchkwhichproc .AND.  
      &           ierr .EQ. 0 ) then  
271                 localEps = 1. _d 0                 localEps = 1. _d 0
272                 call grdchk_getxx( icvrec, TANGENT_SIMULATION,                 call grdchk_getxx( icvrec, TANGENT_SIMULATION,
273       &              itile, jtile, layer,       &              itile, jtile, layer,
274       &              itilepos, jtilepos,       &              itilepos, jtilepos,
275       &              xxmemo_ref, xxmemo_pert, localEps,       &              xxmemo_ref, xxmemo_pert, localEps,
276       &              mythid )       &              ierr, mythid )
             endif  
             _BARRIER  
277    
278  c--  c--
279  c--   2. perform tangent linear run  c--   2. perform tangent linear run
# Line 283  c--   2. perform tangent linear run Line 292  c--   2. perform tangent linear run
292  #endif  #endif
293    
294              call g_the_main_loop( mytime, myiter, mythid )              call g_the_main_loop( mytime, myiter, mythid )
             _BARRIER  
295  #ifdef ALLOW_ADMTLM  #ifdef ALLOW_ADMTLM
296              ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)              ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
297  #else  #else
# Line 292  c--   2. perform tangent linear run Line 300  c--   2. perform tangent linear run
300    
301  c--  c--
302  c--   3. reset control vector  c--   3. reset control vector
             if ( myProcId .EQ. grdchkwhichproc .AND.  
      &           ierr .EQ. 0 ) then  
303                 call grdchk_setxx( icvrec, TANGENT_SIMULATION,                 call grdchk_setxx( icvrec, TANGENT_SIMULATION,
304       &              itile, jtile, layer,       &              itile, jtile, layer,
305       &              itilepos, jtilepos,       &              itilepos, jtilepos,
306       &              xxmemo_ref, mythid )       &              xxmemo_ref, ierr, mythid )
             endif  
             _BARRIER  
307    
308  #endif /* ALLOW_TANGENTLINEAR_RUN */  #endif /* ALLOW_TANGENTLINEAR_RUN */
309    
# Line 312  c--   get control vector component from Line 316  c--   get control vector component from
316  c--   perturb it and write back to file  c--   perturb it and write back to file
317  c--   positive perturbation  c--   positive perturbation
318              localEps = abs(grdchk_eps)              localEps = abs(grdchk_eps)
319              if ( myProcId .EQ. grdchkwhichproc .AND.              call grdchk_getxx( icvrec, FORWARD_SIMULATION,
      &           ierr .EQ. 0 ) then  
                call grdchk_getxx( icvrec, FORWARD_SIMULATION,  
320       &              itile, jtile, layer,       &              itile, jtile, layer,
321       &              itilepos, jtilepos,       &              itilepos, jtilepos,
322       &              xxmemo_ref, xxmemo_pert, localEps,       &              xxmemo_ref, xxmemo_pert, localEps,
323       &              mythid )       &              ierr, mythid )
             endif  
             _BARRIER  
324    
325  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
326              mytime = starttime              mytime = starttime
327              myiter = niter0              myiter = niter0
328    #ifdef ALLOW_AUTODIFF_OPENAD
329                call OpenAD_the_main_loop( mytime, myiter, mythid )
330    #else
331              call the_main_loop( mytime, myiter, mythid )              call the_main_loop( mytime, myiter, mythid )
332  #ifdef ALLOW_ADMTLM  #endif
333    
334    #if (defined (ALLOW_ADMTLM))
335              fcpertplus = objf_state_final(idep,jdep,1,1,1)              fcpertplus = objf_state_final(idep,jdep,1,1,1)
336    #elif (defined (ALLOW_AUTODIFF_OPENAD))
337                fcpertplus = fc%v
338  #else  #else
339              fcpertplus = fc              fcpertplus = fc
340  #endif  #endif
341              print *, 'ph-check fcpertplus = ', fcpertplus              print *, 'ph-check fcpertplus  = ', fcpertplus
             _BARRIER  
342    
343  c--   Reset control vector.  c--   Reset control vector.
344              if ( myProcId .EQ. grdchkwhichproc .AND.              call grdchk_setxx( icvrec, FORWARD_SIMULATION,
      &           ierr .EQ. 0 ) then  
                call grdchk_setxx( icvrec, FORWARD_SIMULATION,  
345       &              itile, jtile, layer,       &              itile, jtile, layer,
346       &              itilepos, jtilepos,       &              itilepos, jtilepos,
347       &              xxmemo_ref, mythid )       &              xxmemo_ref, ierr, mythid )
             endif  
             _BARRIER  
348    
349              fcpertminus = fcref              fcpertminus = fcref
350              print *, 'ph-check fcpertminus = ', fcpertminus              print *, 'ph-check fcpertminus = ', fcpertminus
# Line 352  c--   Reset control vector. Line 354  c--   Reset control vector.
354  c--   get control vector component from file  c--   get control vector component from file
355  c--   perturb it and write back to file  c--   perturb it and write back to file
356  c--   repeat the proceedure for a negative perturbation  c--   repeat the proceedure for a negative perturbation
357                 if ( myProcId .EQ. grdchkwhichproc .AND.                 localEps = - abs(grdchk_eps)
358       &           ierr .EQ. 0 ) then                 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
                   localEps = - abs(grdchk_eps)  
                   call grdchk_getxx( icvrec, FORWARD_SIMULATION,  
359       &                 itile, jtile, layer,       &                 itile, jtile, layer,
360       &                 itilepos, jtilepos,       &                 itilepos, jtilepos,
361       &                 xxmemo_ref, xxmemo_pert, localEps,       &                 xxmemo_ref, xxmemo_pert, localEps,
362       &                 mythid )       &                 ierr, mythid )
                endif  
                _BARRIER  
363    
364  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
365                 mytime = starttime                 mytime = starttime
366                 myiter = niter0                 myiter = niter0
367    #ifdef ALLOW_AUTODIFF_OPENAD
368                   call OpenAD_the_main_loop( mytime, myiter, mythid )
369    #else
370                 call the_main_loop( mytime, myiter, mythid )                 call the_main_loop( mytime, myiter, mythid )
371                 _BARRIER  #endif
372  #ifdef ALLOW_ADMTLM  
373    #if (defined (ALLOW_ADMTLM))
374                 fcpertminus = objf_state_final(idep,jdep,1,1,1)                 fcpertminus = objf_state_final(idep,jdep,1,1,1)
375    #elif (defined (ALLOW_AUTODIFF_OPENAD))
376                   fcpertminus = fc%v
377  #else  #else
378                 fcpertminus = fc                 fcpertminus = fc
379  #endif  #endif
380    
381  c--   Reset control vector.  c--   Reset control vector.
382                 if ( myProcId .EQ. grdchkwhichproc .AND.                 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
      &           ierr .EQ. 0 ) then  
                   call grdchk_setxx( icvrec, FORWARD_SIMULATION,  
383       &                 itile, jtile, layer,       &                 itile, jtile, layer,
384       &                 itilepos, jtilepos,       &                 itilepos, jtilepos,
385       &                 xxmemo_ref, mythid )       &                 xxmemo_ref, ierr, mythid )
                endif  
                _BARRIER  
386    
387  c-- end of if useCentralDiff ...  c-- end of if useCentralDiff ...
388              end if              end if
# Line 412  c*************************************** Line 412  c***************************************
412    
413              if ( myProcId .EQ. grdchkwhichproc .AND.              if ( myProcId .EQ. grdchkwhichproc .AND.
414       &           ierr .EQ. 0 ) then       &           ierr .EQ. 0 ) then
   
415                 tmpplot1(itilepos,jtilepos,layer,itile,jtile)                 tmpplot1(itilepos,jtilepos,layer,itile,jtile)
416       &              = gfd       &              = gfd
417                 tmpplot2(itilepos,jtilepos,layer,itile,jtile)                 tmpplot2(itilepos,jtilepos,layer,itile,jtile)
418       &              = ratio_ad       &              = ratio_ad
419                 tmpplot3(itilepos,jtilepos,layer,itile,jtile)                 tmpplot3(itilepos,jtilepos,layer,itile,jtile)
420       &              = ratio_ftl       &              = ratio_ftl
421                endif
422    
423                if ( ierr .EQ. 0 ) then
424                 fcrmem      ( ichknum ) = fcref                 fcrmem      ( ichknum ) = fcref
425                 fcppmem     ( ichknum ) = fcpertplus                 fcppmem     ( ichknum ) = fcpertplus
426                 fcpmmem     ( ichknum ) = fcpertminus                 fcpmmem     ( ichknum ) = fcpertminus
# Line 443  c*************************************** Line 444  c***************************************
444                 itestmem  ( ichknum ) = itest                 itestmem  ( ichknum ) = itest
445                 ierrmem   ( ichknum ) = ierr                 ierrmem   ( ichknum ) = ierr
446                 icglomem  ( ichknum ) = icglo                 icglomem  ( ichknum ) = icglo
447                endif
448    
449                if ( myProcId .EQ. grdchkwhichproc .AND.
450         &           ierr .EQ. 0 ) then
451    
452                 WRITE(standardmessageunit,'(A)')                 WRITE(standardmessageunit,'(A)')
453       &             'grad-res -------------------------------'       &             'grad-res -------------------------------'
# Line 465  c*************************************** Line 470  c***************************************
470  #endif  #endif
471              endif              endif
472  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
473              WRITE(msgBuf,'(A34,1PE24.14)')              WRITE(msgBuf,'(A30,1PE22.14)')
474       &              ' TLM  precision_derivative_cost =', fcref       &              ' TLM  ref_cost_function      =', fcref
475              CALL PRINT_MESSAGE              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
476       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)       &                          SQUEEZE_RIGHT, myThid )
477              WRITE(msgBuf,'(A34,1PE24.14)')              WRITE(msgBuf,'(A30,1PE22.14)')
478       &              ' TLM  precision_derivative_grad =', ftlxxmemo       &              ' TLM  tangent-lin_grad       =', ftlxxmemo
479              CALL PRINT_MESSAGE              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
480       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)       &                          SQUEEZE_RIGHT, myThid )
481                WRITE(msgBuf,'(A30,1PE22.14)')
482         &              ' TLM  finite-diff_grad       =', gfd
483                CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
484         &                          SQUEEZE_RIGHT, myThid )
485  #else  #else
486              WRITE(msgBuf,'(A30,1PE22.14)')              WRITE(msgBuf,'(A30,1PE22.14)')
487       &              ' ADM  ref_cost_function      =', fcref       &              ' ADM  ref_cost_function      =', fcref
# Line 492  c*************************************** Line 501  c***************************************
501              print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,              print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,
502       &           ', ichknum = ', ichknum       &           ', ichknum = ', ichknum
503    
             _BARRIER  
   
504  c-- else of if ( ichknum ...  c-- else of if ( ichknum ...
505           else           else
506              ierr_grdchk = -1              ierr_grdchk = -1
# Line 504  c-- end of if ( ichknum ... Line 511  c-- end of if ( ichknum ...
511  c-- end of do icomp = ...  c-- end of do icomp = ...
512        enddo        enddo
513    
514        if ( myProcId .EQ. grdchkwhichproc ) then        if (myProcId .EQ. grdchkwhichproc .AND. .NOT.useSingleCpuIO) then
515           CALL WRITE_REC_XYZ_RL(           CALL WRITE_REC_XYZ_RL(
516       &        'grd_findiff'   , tmpplot1, 1, 0, myThid)       &        'grd_findiff'   , tmpplot1, 1, 0, myThid)
517           CALL WRITE_REC_XYZ_RL(           CALL WRITE_REC_XYZ_RL(
# Line 514  c-- end of do icomp = ... Line 521  c-- end of do icomp = ...
521        endif        endif
522    
523  c--   Everyone has to wait for the component to be reset.  c--   Everyone has to wait for the component to be reset.
524        _BARRIER  c     _BARRIER
525    
526  c--   Print the results of the gradient check.  c--   Print the results of the gradient check.
527        call grdchk_print( ichknum, ierr_grdchk, mythid )        call grdchk_print( ichknum, ierr_grdchk, mythid )
Legend:
Removed from v.1.33  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.22