/[MITgcm]/MITgcm_contrib/heimbach/OpenAD/code_regress/grdchk_main.F
ViewVC logotype

Annotation of /MITgcm_contrib/heimbach/OpenAD/code_regress/grdchk_main.F

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


Revision 1.3 - (hide annotations) (download)
Wed Oct 8 20:51:52 2008 UTC (16 years, 9 months ago) by utke
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +1 -1 lines
FILE REMOVED
files made part of common changes or obsoleted

1 utke 1.1 C
2 utke 1.3 C $Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_regress/grdchk_main.F,v 1.2 2008/10/03 14:45:54 utke Exp $
3 utke 1.1 C $Name: $
4    
5     #include "AD_CONFIG.h"
6     #include "CPP_OPTIONS.h"
7    
8     CBOI
9     C
10     C !TITLE: GRADIENT CHECK
11     C !AUTHORS: mitgcm developers ( support@mitgcm.org )
12     C !AFFILIATION: Massachussetts Institute of Technology
13     C !DATE:
14     C !INTRODUCTION: gradient check package
15     c \bv
16     c Compare the gradients calculated by the adjoint model with
17     c finite difference approximations.
18     c
19     C !CALLING SEQUENCE:
20     c
21     c the_model_main
22     c |
23     c |-- ctrl_unpack
24     c |-- adthe_main_loop - unperturbed cost function and
25     c |-- ctrl_pack adjoint gradient are computed here
26     c |
27     c |-- grdchk_main
28     c |
29     c |-- grdchk_init
30     c |-- do icomp=... - loop over control vector elements
31     c |
32     c |-- grdchk_loc - determine location of icomp on grid
33     c |
34     c |-- grdchk_getxx - get control vector component from file
35     c | perturb it and write back to file
36     c |-- grdchk_getadxx - get gradient component calculated
37     c | via adjoint
38     c |-- the_main_loop - forward run and cost evaluation
39     c | with perturbed control vector element
40     c |-- calculate ratio of adj. vs. finite difference gradient
41     c |
42     c |-- grdchk_setxx - Reset control vector element
43     c |
44     c |-- grdchk_print - print results
45     c \ev
46     CEOI
47    
48     CBOP
49     C !ROUTINE: grdchk_main
50     C !INTERFACE:
51     subroutine grdchk_main( mythid )
52    
53     C !DESCRIPTION: \bv
54     c ==================================================================
55     c SUBROUTINE grdchk_main
56     c ==================================================================
57     c o Compare the gradients calculated by the adjoint model with
58     c finite difference approximations.
59     c started: Christian Eckert eckert@mit.edu 24-Feb-2000
60     c continued&finished: heimbach@mit.edu: 13-Jun-2001
61     c changed: mlosch@ocean.mit.edu: 09-May-2002
62     c - added centered difference vs. 1-sided difference option
63     c - improved output format for readability
64     c - added control variable hFacC
65     c heimbach@mit.edu 24-Feb-2003
66     c - added tangent linear gradient checks
67     c - fixes for multiproc. gradient checks
68     c - added more control variables
69     c
70     c ==================================================================
71     c SUBROUTINE grdchk_main
72     c ==================================================================
73     C \ev
74    
75     C !USES:
76     implicit none
77    
78     c == global variables ==
79     #include "SIZE.h"
80     #include "EEPARAMS.h"
81     #include "PARAMS.h"
82     #include "grdchk.h"
83     #include "cost.h"
84     #include "ctrl.h"
85     #ifdef ALLOW_TANGENTLINEAR_RUN
86     #include "g_cost.h"
87     #endif
88    
89     C !INPUT/OUTPUT PARAMETERS:
90     c == routine arguments ==
91     integer mythid
92    
93     #ifdef ALLOW_GRDCHK
94     C !LOCAL VARIABLES:
95     c == local variables ==
96     integer myiter
97     _RL mytime
98     integer bi, itlo, ithi
99     integer bj, jtlo, jthi
100     integer i, imin, imax
101     integer j, jmin, jmax
102     integer k
103    
104     integer icomp
105     integer ichknum
106     integer icvrec
107     integer jtile
108     integer itile
109     integer layer
110     integer obcspos
111     integer itilepos
112     integer jtilepos
113     integer icglo
114     integer itest
115     integer ierr
116     integer ierr_grdchk
117     _RL gfd
118     _RL fcref
119     _RL fcpertplus, fcpertminus
120     _RL ratio_ad
121     _RL ratio_ftl
122     _RL xxmemo_ref
123     _RL xxmemo_pert
124     _RL adxxmemo
125     _RL ftlxxmemo
126     _RL localEps
127     _RL grdchk_epsfac
128    
129     _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
130     _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
131     _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
132    
133     CHARACTER*(MAX_LEN_MBUF) msgBuf
134    
135     c == end of interface ==
136     CEOP
137    
138     c-- Set the loop ranges.
139     jtlo = mybylo(mythid)
140     jthi = mybyhi(mythid)
141     itlo = mybxlo(mythid)
142     ithi = mybxhi(mythid)
143     jmin = 1
144     jmax = sny
145     imin = 1
146     imax = snx
147    
148     print *, 'ph-check entering grdchk_main '
149    
150     c-- initialise variables
151     call grdchk_init( mythid )
152    
153     c-- Compute the adjoint models' gradients.
154     c-- Compute the unperturbed cost function.
155     cph Gradient via adjoint has already been computed,
156     cph and so has unperturbed cost function,
157     cph assuming all xx_ fields are initialised to zero.
158    
159     ierr_grdchk = 0
160 utke 1.2 adxxmemo = 0.
161     ftlxxmemo = 0.
162 utke 1.1 #ifdef ALLOW_ADMTLM
163     fcref = objf_state_final(idep,jdep,1,1,1)
164     #else
165     fcref = fc%v
166     #endif
167    
168     print *, 'ph-check fcref = ', fcref
169    
170     do bj = jtlo, jthi
171     do bi = itlo, ithi
172     do k = 1, nr
173     do j = jmin, jmax
174     do i = imin, imax
175     tmpplot1(i,j,k,bi,bj) = 0. _d 0
176     tmpplot2(i,j,k,bi,bj) = 0. _d 0
177     tmpplot3(i,j,k,bi,bj) = 0. _d 0
178     end do
179     end do
180     end do
181     end do
182     end do
183    
184     if ( useCentralDiff ) then
185     grdchk_epsfac = 2. _d 0
186     else
187     grdchk_epsfac = 1. _d 0
188     end if
189    
190 utke 1.2 WRITE(standardmessageunit,'(A)')
191     & 'grad-res -------------------------------'
192     WRITE(standardmessageunit,'(2a)')
193 utke 1.1 & ' grad-res proc # i j k bi bj iobc',
194     & ' fc ref fc + eps fc - eps'
195     #ifdef ALLOW_TANGENTLINEAR_RUN
196 utke 1.2 WRITE(standardmessageunit,'(2a)')
197 utke 1.1 & ' grad-res proc # i j k bi bj iobc',
198     & ' tlm grad fd grad 1 - fd/tlm'
199     #else
200 utke 1.2 WRITE(standardmessageunit,'(2a)')
201 utke 1.1 & ' grad-res proc # i j k bi bj iobc',
202     & ' adj grad fd grad 1 - fd/adj'
203     #endif
204    
205     c-- Compute the finite difference approximations.
206     c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
207     c-- gradient checks.
208    
209     if ( nbeg .EQ. 0 )
210     & call grdchk_get_position( mythid )
211    
212     do icomp = nbeg, nend, nstep
213    
214     ichknum = (icomp - nbeg)/nstep + 1
215    
216     cph(
217     cph-print print *, 'ph-grd _main: nbeg, icomp, ichknum ',
218     cph-print & nbeg, icomp, ichknum
219     cph)
220     if (ichknum .le. maxgrdchecks ) then
221    
222     c-- Determine the location of icomp on the grid.
223     if ( myProcId .EQ. grdchkwhichproc ) then
224     call grdchk_loc( icomp, ichknum,
225     & icvrec, itile, jtile, layer, obcspos,
226     & itilepos, jtilepos, icglo, itest, ierr,
227     & mythid )
228     cph(
229     cph-print print *, 'ph-grd ----- back from loc -----',
230     cph-print & icvrec, itilepos, jtilepos, layer, obcspos
231     cph)
232     endif
233     _BARRIER
234    
235     c******************************************************
236     c-- (A): get gradient component calculated via adjoint
237     c******************************************************
238    
239     c-- get gradient component calculated via adjoint
240     if ( myProcId .EQ. grdchkwhichproc .AND.
241     & ierr .EQ. 0 ) then
242     call grdchk_getadxx( icvrec,
243     & itile, jtile, layer,
244     & itilepos, jtilepos,
245     & adxxmemo, mythid )
246     endif
247     _BARRIER
248    
249     #ifdef ALLOW_TANGENTLINEAR_RUN
250     c******************************************************
251     c-- (B): Get gradient component g_fc from tangent linear run:
252     c******************************************************
253     c--
254     c-- 1. perturb control vector component: xx(i)=1.
255    
256     if ( myProcId .EQ. grdchkwhichproc .AND.
257     & ierr .EQ. 0 ) then
258     localEps = 1. _d 0
259     call grdchk_getxx( icvrec, TANGENT_SIMULATION,
260     & itile, jtile, layer,
261     & itilepos, jtilepos,
262     & xxmemo_ref, xxmemo_pert, localEps,
263     & mythid )
264     endif
265     _BARRIER
266    
267     c--
268     c-- 2. perform tangent linear run
269     mytime = starttime
270     myiter = niter0
271     #ifdef ALLOW_ADMTLM
272     do k=1,4*Nr+1
273     do j=1,sny
274     do i=1,snx
275     g_objf_state_final(i,j,1,1,k) = 0.
276     enddo
277     enddo
278     enddo
279     #else
280     g_fc = 0.
281     #endif
282    
283     call g_the_main_loop( mytime, myiter, mythid )
284     _BARRIER
285     #ifdef ALLOW_ADMTLM
286     ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
287     #else
288     ftlxxmemo = g_fc
289     #endif
290    
291     c--
292     c-- 3. reset control vector
293     if ( myProcId .EQ. grdchkwhichproc .AND.
294     & ierr .EQ. 0 ) then
295     call grdchk_setxx( icvrec, TANGENT_SIMULATION,
296     & itile, jtile, layer,
297     & itilepos, jtilepos,
298     & xxmemo_ref, mythid )
299     endif
300     _BARRIER
301    
302     #endif /* ALLOW_TANGENTLINEAR_RUN */
303    
304    
305     c******************************************************
306     c-- (C): Get gradient via finite difference perturbation
307     c******************************************************
308    
309     c-- get control vector component from file
310     c-- perturb it and write back to file
311     c-- positive perturbation
312     localEps = abs(grdchk_eps)
313     if ( myProcId .EQ. grdchkwhichproc .AND.
314     & ierr .EQ. 0 ) then
315     call grdchk_getxx( icvrec, FORWARD_SIMULATION,
316     & itile, jtile, layer,
317     & itilepos, jtilepos,
318     & xxmemo_ref, xxmemo_pert, localEps,
319     & mythid )
320     endif
321     _BARRIER
322    
323     c-- forward run with perturbed control vector
324     mytime = starttime
325     myiter = niter0
326     call OpenAD_the_main_loop( mytime, myiter, mythid )
327     #ifdef ALLOW_ADMTLM
328     fcpertplus = objf_state_final(idep,jdep,1,1,1)
329     #else
330     fcpertplus = fc%v
331     #endif
332     print *, 'ph-check fcpertplus = ', fcpertplus
333     _BARRIER
334    
335     c-- Reset control vector.
336     if ( myProcId .EQ. grdchkwhichproc .AND.
337     & ierr .EQ. 0 ) then
338     call grdchk_setxx( icvrec, FORWARD_SIMULATION,
339     & itile, jtile, layer,
340     & itilepos, jtilepos,
341     & xxmemo_ref, mythid )
342     endif
343     _BARRIER
344    
345     fcpertminus = fcref
346     print *, 'ph-check fcpertminus = ', fcpertminus
347    
348     if ( useCentralDiff ) then
349    
350     c-- get control vector component from file
351     c-- perturb it and write back to file
352     c-- repeat the proceedure for a negative perturbation
353     if ( myProcId .EQ. grdchkwhichproc .AND.
354     & ierr .EQ. 0 ) then
355     localEps = - abs(grdchk_eps)
356     call grdchk_getxx( icvrec, FORWARD_SIMULATION,
357     & itile, jtile, layer,
358     & itilepos, jtilepos,
359     & xxmemo_ref, xxmemo_pert, localEps,
360     & mythid )
361     endif
362     _BARRIER
363    
364     c-- forward run with perturbed control vector
365     mytime = starttime
366     myiter = niter0
367     call OpenAD_the_main_loop( mytime, myiter, mythid )
368     _BARRIER
369     #ifdef ALLOW_ADMTLM
370     fcpertminus = objf_state_final(idep,jdep,1,1,1)
371     #else
372     fcpertminus = fc%v
373     #endif
374    
375     c-- Reset control vector.
376     if ( myProcId .EQ. grdchkwhichproc .AND.
377     & ierr .EQ. 0 ) then
378     call grdchk_setxx( icvrec, FORWARD_SIMULATION,
379     & itile, jtile, layer,
380     & itilepos, jtilepos,
381     & xxmemo_ref, mythid )
382     endif
383     _BARRIER
384    
385     c-- end of if useCentralDiff ...
386     end if
387    
388     c******************************************************
389     c-- (D): calculate relative differences between gradients
390     c******************************************************
391    
392     if ( myProcId .EQ. grdchkwhichproc .AND.
393     & ierr .EQ. 0 ) then
394    
395     if ( grdchk_eps .eq. 0. ) then
396     gfd = (fcpertplus-fcpertminus)
397     else
398     gfd = (fcpertplus-fcpertminus)
399     & /(grdchk_epsfac*grdchk_eps)
400     endif
401    
402     if ( adxxmemo .eq. 0. ) then
403     ratio_ad = abs( adxxmemo - gfd )
404     else
405     ratio_ad = 1. - gfd/adxxmemo
406     endif
407    
408     if ( ftlxxmemo .eq. 0. ) then
409     ratio_ftl = abs( ftlxxmemo - gfd )
410     else
411     ratio_ftl = 1. - gfd/ftlxxmemo
412     endif
413    
414     tmpplot1(itilepos,jtilepos,layer,itile,jtile)
415     & = gfd
416     tmpplot2(itilepos,jtilepos,layer,itile,jtile)
417     & = ratio_ad
418     tmpplot3(itilepos,jtilepos,layer,itile,jtile)
419     & = ratio_ftl
420    
421     fcrmem ( ichknum ) = fcref
422     fcppmem ( ichknum ) = fcpertplus
423     fcpmmem ( ichknum ) = fcpertminus
424     xxmemref ( ichknum ) = xxmemo_ref
425     xxmempert ( ichknum ) = xxmemo_pert
426     gfdmem ( ichknum ) = gfd
427     adxxmem ( ichknum ) = adxxmemo
428     ftlxxmem ( ichknum ) = ftlxxmemo
429     ratioadmem ( ichknum ) = ratio_ad
430     ratioftlmem ( ichknum ) = ratio_ftl
431    
432     irecmem ( ichknum ) = icvrec
433     bimem ( ichknum ) = itile
434     bjmem ( ichknum ) = jtile
435     ilocmem ( ichknum ) = itilepos
436     jlocmem ( ichknum ) = jtilepos
437     klocmem ( ichknum ) = layer
438     iobcsmem ( ichknum ) = obcspos
439     icompmem ( ichknum ) = icomp
440     ichkmem ( ichknum ) = ichknum
441     itestmem ( ichknum ) = itest
442     ierrmem ( ichknum ) = ierr
443     icglomem ( ichknum ) = icglo
444    
445 utke 1.2 WRITE(standardmessageunit,'(A)')
446     & 'grad-res -------------------------------'
447     WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
448     & ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
449 utke 1.1 & layer,itile,jtile,obcspos,
450     & fcref, fcpertplus, fcpertminus
451     #ifdef ALLOW_TANGENTLINEAR_RUN
452 utke 1.2 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
453     & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
454 utke 1.1 & icompmem(ichknum),itestmem(ichknum),
455     & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
456     & ftlxxmemo, gfd, ratio_ftl
457     WRITE(msgBuf,'(A34,2(1PE24.14,X))')
458     & 'precision_grdchk_result TLM ', fcref, ftlxxmemo
459     CALL PRINT_MESSAGE
460     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
461 utke 1.2 c
462     WRITE(msgBuf,'(A34,1PE24.14)')
463     & ' TLM precision_derivative_cost =', fcref
464     CALL PRINT_MESSAGE
465     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
466     WRITE(msgBuf,'(A34,1PE24.14)')
467     & ' TLM precision_derivative_grad =', ftlxxmemo
468     CALL PRINT_MESSAGE
469     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
470 utke 1.1 #else
471 utke 1.2 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
472     & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
473 utke 1.1 & icompmem(ichknum),itestmem(ichknum),
474     & bimem(ichknum),bjmem(ichknum),obcspos,
475     & adxxmemo, gfd, ratio_ad
476 utke 1.2 c WRITE(msgBuf,'(A34,2(1PE24.14,X))')
477     c & 'precision_grdchk_result ADM ', fcref, adxxmemo
478     c CALL PRINT_MESSAGE
479     c & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
480     WRITE(msgBuf,'(A34,1PE24.14)')
481     & ' ADM precision_derivative_cost =', fcref
482     CALL PRINT_MESSAGE
483     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
484     WRITE(msgBuf,'(A34,1PE24.14)')
485     & ' ADM precision_derivative_grad =', adxxmemo
486 utke 1.1 CALL PRINT_MESSAGE
487     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
488     #endif
489    
490     endif
491    
492     print *, 'ph-grd ierr ---------------------------'
493     print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
494     & ', ichknum = ', ichknum
495    
496     _BARRIER
497    
498     c-- else of if ( ichknum ...
499     else
500     ierr_grdchk = -1
501    
502     c-- end of if ( ichknum ...
503     endif
504    
505     c-- end of do icomp = ...
506     enddo
507    
508     if ( myProcId .EQ. grdchkwhichproc ) then
509     CALL WRITE_REC_XYZ_RL(
510     & 'grd_findiff' , tmpplot1, 1, 0, myThid)
511     CALL WRITE_REC_XYZ_RL(
512     & 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
513     CALL WRITE_REC_XYZ_RL(
514     & 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
515     endif
516    
517     c-- Everyone has to wait for the component to be reset.
518     _BARRIER
519    
520     c-- Print the results of the gradient check.
521     call grdchk_print( ichknum, ierr_grdchk, mythid )
522    
523     #endif /* ALLOW_GRDCHK */
524    
525     end

  ViewVC Help
Powered by ViewVC 1.1.22