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

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

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


Revision 1.1 - (hide annotations) (download)
Tue Aug 1 15:44:29 2006 UTC (20 years, 3 months ago) by utke
Branch: MAIN
gmredi enabled setup

1 utke 1.1 C
2     C $Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_shallow_openad/grdchk_main.F,v 1.2 2006/01/18 23:43:11 utke Exp $
3     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     #ifdef ALLOW_TANGENTLINEAR_RUN
85     #include "g_cost.h"
86     #endif
87    
88     C !INPUT/OUTPUT PARAMETERS:
89     c == routine arguments ==
90     integer mythid
91    
92     #ifdef ALLOW_GRDCHK
93     C !LOCAL VARIABLES:
94     c == local variables ==
95     integer myiter
96     _RL mytime
97     integer bi, itlo, ithi
98     integer bj, jtlo, jthi
99     integer i, imin, imax
100     integer j, jmin, jmax
101     integer k
102    
103     integer icomp
104     integer ichknum
105     integer icvrec
106     integer jtile
107     integer itile
108     integer layer
109     integer obcspos
110     integer itilepos
111     integer jtilepos
112     integer itest
113     integer ierr
114     integer ierr_grdchk
115     _RL gfd
116     _RL fcref
117     _RL fcpertplus, fcpertminus
118     _RL ratio_ad
119     _RL ratio_ftl
120     _RL xxmemo_ref
121     _RL xxmemo_pert
122     _RL adxxmemo
123     _RL ftlxxmemo
124     _RL localEps
125     _RL grdchk_epsfac
126    
127     _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
128     _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
129     _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
130    
131     CHARACTER*(MAX_LEN_MBUF) msgBuf
132    
133     c == end of interface ==
134     CEOP
135    
136     c-- Set the loop ranges.
137     jtlo = mybylo(mythid)
138     jthi = mybyhi(mythid)
139     itlo = mybxlo(mythid)
140     ithi = mybxhi(mythid)
141     jmin = 1
142     jmax = sny
143     imin = 1
144     imax = snx
145    
146     print *, 'ph-check entering grdchk_main '
147    
148     c-- initialise variables
149     call grdchk_init( mythid )
150    
151     c-- Compute the adjoint models' gradients.
152     c-- Compute the unperturbed cost function.
153     cph Gradient via adjoint has already been computed,
154     cph and so has unperturbed cost function,
155     cph assuming all xx_ fields are initialised to zero.
156    
157     ierr_grdchk = 0
158     cphadmtlm(
159     fcref = fc%v
160     cphadmtlm fcref = objf_state_final(45,4,1,1,1)
161     cphadmtlm)
162    
163     print *, 'ph-check fcref = ', fcref
164    
165     do bj = jtlo, jthi
166     do bi = itlo, ithi
167     do k = 1, nr
168     do j = jmin, jmax
169     do i = imin, imax
170     tmpplot1(i,j,k,bi,bj) = 0. _d 0
171     tmpplot2(i,j,k,bi,bj) = 0. _d 0
172     tmpplot3(i,j,k,bi,bj) = 0. _d 0
173     end do
174     end do
175     end do
176     end do
177     end do
178    
179     if ( useCentralDiff ) then
180     grdchk_epsfac = 2. _d 0
181     else
182     grdchk_epsfac = 1. _d 0
183     end if
184    
185     print *, 'grad-res -------------------------------'
186     print ('(2a)'),
187     & ' grad-res proc # i j k fc ref',
188     & ' fc + eps fc - eps'
189     #ifdef ALLOW_TANGENTLINEAR_RUN
190     print ('(2a)'),
191     & ' grad-res proc # i j k tlm grad',
192     & ' fd grad 1 - fd/tlm'
193     #else
194     print ('(2a)'),
195     & ' grad-res proc # i j k adj grad',
196     & ' fd grad 1 - fd/adj'
197     #endif
198    
199     c-- Compute the finite difference approximations.
200     c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
201     c-- gradient checks.
202    
203     if ( nbeg .EQ. 0 ) call grdchk_get_position( mythid )
204    
205     do icomp = nbeg, nend, nstep
206    
207     ichknum = (icomp - nbeg)/nstep + 1
208    
209     if (ichknum .le. maxgrdchecks ) then
210    
211     c-- Determine the location of icomp on the grid.
212     if ( myProcId .EQ. grdchkwhichproc ) then
213     call grdchk_loc( icomp, ichknum,
214     & icvrec, itile, jtile, layer, obcspos,
215     & itilepos, jtilepos, itest, ierr,
216     & mythid )
217     endif
218     _BARRIER
219    
220     c******************************************************
221     c-- (A): get gradient component calculated via adjoint
222     c******************************************************
223    
224     c-- get gradient component calculated via adjoint
225     if ( myProcId .EQ. grdchkwhichproc .AND.
226     & ierr .EQ. 0 ) then
227     call grdchk_getadxx( icvrec,
228     & itile, jtile, layer,
229     & itilepos, jtilepos,
230     & adxxmemo, mythid )
231     endif
232     _BARRIER
233    
234     #ifdef ALLOW_TANGENTLINEAR_RUN
235     c******************************************************
236     c-- (B): Get gradient component g_fc from tangent linear run:
237     c******************************************************
238     c--
239     c-- 1. perturb control vector component: xx(i)=1.
240    
241     if ( myProcId .EQ. grdchkwhichproc .AND.
242     & ierr .EQ. 0 ) then
243     localEps = 1. _d 0
244     call grdchk_getxx( icvrec, TANGENT_SIMULATION,
245     & itile, jtile, layer,
246     & itilepos, jtilepos,
247     & xxmemo_ref, xxmemo_pert, localEps,
248     & mythid )
249     endif
250     _BARRIER
251    
252     c--
253     c-- 2. perform tangent linear run
254     mytime = starttime
255     myiter = niter0
256     cphadmtlm(
257     g_fc = 0.
258     cphadmtlm do j=1,sny
259     cphadmtlm do i=1,snx
260     cphadmtlm g_objf_state_final(i,j,1,1,1) = 0.
261     cphadmtlm g_objf_state_final(i,j,1,1,2) = 0.
262     cphadmtlm enddo
263     cphadmtlm enddo
264     cphadmtlm)
265     call g_the_main_loop( mytime, myiter, mythid )
266     cphadmtlm(
267     ftlxxmemo = g_fc
268     cphadmtlm ftlxxmemo = g_objf_state_final(45,4,1,1,1)
269     cphadmtlm)
270     _BARRIER
271     c--
272     c-- 3. reset control vector
273     if ( myProcId .EQ. grdchkwhichproc .AND.
274     & ierr .EQ. 0 ) then
275     call grdchk_setxx( icvrec, TANGENT_SIMULATION,
276     & itile, jtile, layer,
277     & itilepos, jtilepos,
278     & xxmemo_ref, mythid )
279     endif
280     _BARRIER
281    
282     #endif /* ALLOW_TANGENTLINEAR_RUN */
283    
284    
285     c******************************************************
286     c-- (C): Get gradient via finite difference perturbation
287     c******************************************************
288    
289     c-- get control vector component from file
290     c-- perturb it and write back to file
291     c-- positive perturbation
292     localEps = abs(grdchk_eps)
293     if ( myProcId .EQ. grdchkwhichproc .AND.
294     & ierr .EQ. 0 ) then
295     call grdchk_getxx( icvrec, FORWARD_SIMULATION,
296     & itile, jtile, layer,
297     & itilepos, jtilepos,
298     & xxmemo_ref, xxmemo_pert, localEps,
299     & mythid )
300     endif
301     _BARRIER
302    
303     c-- forward run with perturbed control vector
304     mytime = starttime
305     myiter = niter0
306     call OpenAD_the_main_loop( mytime, myiter, mythid )
307     cphadmtlm(
308     fcpertplus = fc%v
309     cphadmtlm fcpertplus = objf_state_final(45,4,1,1,1)
310     cphadmtlm)
311     print *, 'ph-check fcpertplus = ', fcpertplus
312     _BARRIER
313    
314     c-- Reset control vector.
315     if ( myProcId .EQ. grdchkwhichproc .AND.
316     & ierr .EQ. 0 ) then
317     call grdchk_setxx( icvrec, FORWARD_SIMULATION,
318     & itile, jtile, layer,
319     & itilepos, jtilepos,
320     & xxmemo_ref, mythid )
321     endif
322     _BARRIER
323    
324     fcpertminus = fcref
325     print *, 'ph-check fcpertminus = ', fcpertminus
326    
327     if ( useCentralDiff ) then
328    
329     c-- get control vector component from file
330     c-- perturb it and write back to file
331     c-- repeat the proceedure for a negative perturbation
332     if ( myProcId .EQ. grdchkwhichproc .AND.
333     & ierr .EQ. 0 ) then
334     localEps = - abs(grdchk_eps)
335     call grdchk_getxx( icvrec, FORWARD_SIMULATION,
336     & itile, jtile, layer,
337     & itilepos, jtilepos,
338     & xxmemo_ref, xxmemo_pert, localEps,
339     & mythid )
340     endif
341     _BARRIER
342    
343     c-- forward run with perturbed control vector
344     mytime = starttime
345     myiter = niter0
346     call OpenAD_the_main_loop( mytime, myiter, mythid )
347     _BARRIER
348     fcpertminus = fc%v
349    
350     c-- Reset control vector.
351     if ( myProcId .EQ. grdchkwhichproc .AND.
352     & ierr .EQ. 0 ) then
353     call grdchk_setxx( icvrec, FORWARD_SIMULATION,
354     & itile, jtile, layer,
355     & itilepos, jtilepos,
356     & xxmemo_ref, mythid )
357     endif
358     _BARRIER
359    
360     c-- end of if useCentralDiff ...
361     end if
362    
363     c******************************************************
364     c-- (D): calculate relative differences between gradients
365     c******************************************************
366    
367     if ( myProcId .EQ. grdchkwhichproc .AND.
368     & ierr .EQ. 0 ) then
369    
370     if ( grdchk_eps .eq. 0. ) then
371     gfd = (fcpertplus-fcpertminus)
372     else
373     gfd = (fcpertplus-fcpertminus)
374     & /(grdchk_epsfac*grdchk_eps)
375     endif
376    
377     if ( adxxmemo .eq. 0. ) then
378     ratio_ad = abs( adxxmemo - gfd )
379     else
380     ratio_ad = 1. - gfd/adxxmemo
381     endif
382    
383     if ( ftlxxmemo .eq. 0. ) then
384     ratio_ftl = abs( ftlxxmemo - gfd )
385     else
386     ratio_ftl = 1. - gfd/ftlxxmemo
387     endif
388    
389     tmpplot1(itilepos,jtilepos,layer,itile,jtile)
390     & = gfd
391     tmpplot2(itilepos,jtilepos,layer,itile,jtile)
392     & = ratio_ad
393     tmpplot3(itilepos,jtilepos,layer,itile,jtile)
394     & = ratio_ftl
395    
396     fcrmem ( ichknum ) = fcref
397     fcppmem ( ichknum ) = fcpertplus
398     fcpmmem ( ichknum ) = fcpertminus
399     xxmemref ( ichknum ) = xxmemo_ref
400     xxmempert ( ichknum ) = xxmemo_pert
401     gfdmem ( ichknum ) = gfd
402     adxxmem ( ichknum ) = adxxmemo
403     ftlxxmem ( ichknum ) = ftlxxmemo
404     ratioadmem ( ichknum ) = ratio_ad
405     ratioftlmem ( ichknum ) = ratio_ftl
406    
407     irecmem ( ichknum ) = icvrec
408     bimem ( ichknum ) = itile
409     bjmem ( ichknum ) = jtile
410     ilocmem ( ichknum ) = itilepos
411     jlocmem ( ichknum ) = jtilepos
412     klocmem ( ichknum ) = layer
413     iobcsmem ( ichknum ) = obcspos
414     icompmem ( ichknum ) = icomp
415     ichkmem ( ichknum ) = ichknum
416     itestmem ( ichknum ) = itest
417     ierrmem ( ichknum ) = ierr
418    
419     print *, 'grad-res -------------------------------'
420     print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',
421     & myprocid,ichknum,itilepos,jtilepos,layer,
422     & fcref, fcpertplus, fcpertminus
423     #ifdef ALLOW_TANGENTLINEAR_RUN
424     print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',
425     & myprocid,ichknum,ichkmem(ichknum),
426     & icompmem(ichknum),itestmem(ichknum),
427     & ftlxxmemo, gfd, ratio_ftl
428     WRITE(msgBuf,'(A34,2(1PE24.14,X))')
429     & 'precision_grdchk_result TLM ', fcref, ftlxxmemo
430     CALL PRINT_MESSAGE
431     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
432     #else
433     print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',
434     & myprocid,ichknum,ichkmem(ichknum),
435     & icompmem(ichknum),itestmem(ichknum),
436     & adxxmemo, gfd, ratio_ad
437     WRITE(msgBuf,'(A34,2(1PE24.14,X))')
438     & 'precision_grdchk_result ADM ', fcref, adxxmemo
439     CALL PRINT_MESSAGE
440     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
441     #endif
442    
443     endif
444    
445     print *, 'ph-grd ierr ---------------------------'
446     print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
447     & ', ichknum = ', ichknum
448    
449     _BARRIER
450    
451     c-- else of if ( ichknum ...
452     else
453     ierr_grdchk = -1
454    
455     c-- end of if ( ichknum ...
456     endif
457    
458     c-- end of do icomp = ...
459     enddo
460    
461     if ( myProcId .EQ. grdchkwhichproc ) then
462     CALL WRITE_REC_XYZ_RL(
463     & 'grd_findiff' , tmpplot1, 1, 0, myThid)
464     CALL WRITE_REC_XYZ_RL(
465     & 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
466     CALL WRITE_REC_XYZ_RL(
467     & 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
468     endif
469    
470     c-- Everyone has to wait for the component to be reset.
471     _BARRIER
472    
473     c-- Print the results of the gradient check.
474     call grdchk_print( ichknum, ierr_grdchk, mythid )
475    
476     #endif /* ALLOW_GRDCHK */
477    
478     end

  ViewVC Help
Powered by ViewVC 1.1.22