/[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.1 - (hide annotations) (download)
Tue Nov 20 15:20:58 2007 UTC (17 years, 8 months ago) by utke
Branch: MAIN
adding regression test code

1 utke 1.1 C
2     C $Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_heat_transport/grdchk_main.F,v 1.3 2007/05/09 21:40:03 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     #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     #ifdef ALLOW_ADMTLM
161     fcref = objf_state_final(idep,jdep,1,1,1)
162     #else
163     fcref = fc%v
164     #endif
165    
166     print *, 'ph-check fcref = ', fcref
167    
168     do bj = jtlo, jthi
169     do bi = itlo, ithi
170     do k = 1, nr
171     do j = jmin, jmax
172     do i = imin, imax
173     tmpplot1(i,j,k,bi,bj) = 0. _d 0
174     tmpplot2(i,j,k,bi,bj) = 0. _d 0
175     tmpplot3(i,j,k,bi,bj) = 0. _d 0
176     end do
177     end do
178     end do
179     end do
180     end do
181    
182     if ( useCentralDiff ) then
183     grdchk_epsfac = 2. _d 0
184     else
185     grdchk_epsfac = 1. _d 0
186     end if
187    
188     print *, 'grad-res -------------------------------'
189     print '(2a)',
190     & ' grad-res proc # i j k bi bj iobc',
191     & ' fc ref fc + eps fc - eps'
192     #ifdef ALLOW_TANGENTLINEAR_RUN
193     print '(2a)',
194     & ' grad-res proc # i j k bi bj iobc',
195     & ' tlm grad fd grad 1 - fd/tlm'
196     #else
197     print '(2a)',
198     & ' grad-res proc # i j k bi bj iobc',
199     & ' adj grad fd grad 1 - fd/adj'
200     #endif
201    
202     c-- Compute the finite difference approximations.
203     c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
204     c-- gradient checks.
205    
206     if ( nbeg .EQ. 0 )
207     & call grdchk_get_position( mythid )
208    
209     do icomp = nbeg, nend, nstep
210    
211     ichknum = (icomp - nbeg)/nstep + 1
212    
213     cph(
214     cph-print print *, 'ph-grd _main: nbeg, icomp, ichknum ',
215     cph-print & nbeg, icomp, ichknum
216     cph)
217     if (ichknum .le. maxgrdchecks ) then
218    
219     c-- Determine the location of icomp on the grid.
220     if ( myProcId .EQ. grdchkwhichproc ) then
221     call grdchk_loc( icomp, ichknum,
222     & icvrec, itile, jtile, layer, obcspos,
223     & itilepos, jtilepos, icglo, itest, ierr,
224     & mythid )
225     cph(
226     cph-print print *, 'ph-grd ----- back from loc -----',
227     cph-print & icvrec, itilepos, jtilepos, layer, obcspos
228     cph)
229     endif
230     _BARRIER
231    
232     c******************************************************
233     c-- (A): get gradient component calculated via adjoint
234     c******************************************************
235    
236     c-- get gradient component calculated via adjoint
237     if ( myProcId .EQ. grdchkwhichproc .AND.
238     & ierr .EQ. 0 ) then
239     call grdchk_getadxx( icvrec,
240     & itile, jtile, layer,
241     & itilepos, jtilepos,
242     & adxxmemo, mythid )
243     endif
244     _BARRIER
245    
246     #ifdef ALLOW_TANGENTLINEAR_RUN
247     c******************************************************
248     c-- (B): Get gradient component g_fc from tangent linear run:
249     c******************************************************
250     c--
251     c-- 1. perturb control vector component: xx(i)=1.
252    
253     if ( myProcId .EQ. grdchkwhichproc .AND.
254     & ierr .EQ. 0 ) then
255     localEps = 1. _d 0
256     call grdchk_getxx( icvrec, TANGENT_SIMULATION,
257     & itile, jtile, layer,
258     & itilepos, jtilepos,
259     & xxmemo_ref, xxmemo_pert, localEps,
260     & mythid )
261     endif
262     _BARRIER
263    
264     c--
265     c-- 2. perform tangent linear run
266     mytime = starttime
267     myiter = niter0
268     #ifdef ALLOW_ADMTLM
269     do k=1,4*Nr+1
270     do j=1,sny
271     do i=1,snx
272     g_objf_state_final(i,j,1,1,k) = 0.
273     enddo
274     enddo
275     enddo
276     #else
277     g_fc = 0.
278     #endif
279    
280     call g_the_main_loop( mytime, myiter, mythid )
281     _BARRIER
282     #ifdef ALLOW_ADMTLM
283     ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
284     #else
285     ftlxxmemo = g_fc
286     #endif
287    
288     c--
289     c-- 3. reset control vector
290     if ( myProcId .EQ. grdchkwhichproc .AND.
291     & ierr .EQ. 0 ) then
292     call grdchk_setxx( icvrec, TANGENT_SIMULATION,
293     & itile, jtile, layer,
294     & itilepos, jtilepos,
295     & xxmemo_ref, mythid )
296     endif
297     _BARRIER
298    
299     #endif /* ALLOW_TANGENTLINEAR_RUN */
300    
301    
302     c******************************************************
303     c-- (C): Get gradient via finite difference perturbation
304     c******************************************************
305    
306     c-- get control vector component from file
307     c-- perturb it and write back to file
308     c-- positive perturbation
309     localEps = abs(grdchk_eps)
310     if ( myProcId .EQ. grdchkwhichproc .AND.
311     & ierr .EQ. 0 ) then
312     call grdchk_getxx( icvrec, FORWARD_SIMULATION,
313     & itile, jtile, layer,
314     & itilepos, jtilepos,
315     & xxmemo_ref, xxmemo_pert, localEps,
316     & mythid )
317     endif
318     _BARRIER
319    
320     c-- forward run with perturbed control vector
321     mytime = starttime
322     myiter = niter0
323     call OpenAD_the_main_loop( mytime, myiter, mythid )
324     #ifdef ALLOW_ADMTLM
325     fcpertplus = objf_state_final(idep,jdep,1,1,1)
326     #else
327     fcpertplus = fc%v
328     #endif
329     print *, 'ph-check fcpertplus = ', fcpertplus
330     _BARRIER
331    
332     c-- Reset control vector.
333     if ( myProcId .EQ. grdchkwhichproc .AND.
334     & ierr .EQ. 0 ) then
335     call grdchk_setxx( icvrec, FORWARD_SIMULATION,
336     & itile, jtile, layer,
337     & itilepos, jtilepos,
338     & xxmemo_ref, mythid )
339     endif
340     _BARRIER
341    
342     fcpertminus = fcref
343     print *, 'ph-check fcpertminus = ', fcpertminus
344    
345     if ( useCentralDiff ) then
346    
347     c-- get control vector component from file
348     c-- perturb it and write back to file
349     c-- repeat the proceedure for a negative perturbation
350     if ( myProcId .EQ. grdchkwhichproc .AND.
351     & ierr .EQ. 0 ) then
352     localEps = - abs(grdchk_eps)
353     call grdchk_getxx( icvrec, FORWARD_SIMULATION,
354     & itile, jtile, layer,
355     & itilepos, jtilepos,
356     & xxmemo_ref, xxmemo_pert, localEps,
357     & mythid )
358     endif
359     _BARRIER
360    
361     c-- forward run with perturbed control vector
362     mytime = starttime
363     myiter = niter0
364     call OpenAD_the_main_loop( mytime, myiter, mythid )
365     _BARRIER
366     #ifdef ALLOW_ADMTLM
367     fcpertminus = objf_state_final(idep,jdep,1,1,1)
368     #else
369     fcpertminus = fc%v
370     #endif
371    
372     c-- Reset control vector.
373     if ( myProcId .EQ. grdchkwhichproc .AND.
374     & ierr .EQ. 0 ) then
375     call grdchk_setxx( icvrec, FORWARD_SIMULATION,
376     & itile, jtile, layer,
377     & itilepos, jtilepos,
378     & xxmemo_ref, mythid )
379     endif
380     _BARRIER
381    
382     c-- end of if useCentralDiff ...
383     end if
384    
385     c******************************************************
386     c-- (D): calculate relative differences between gradients
387     c******************************************************
388    
389     if ( myProcId .EQ. grdchkwhichproc .AND.
390     & ierr .EQ. 0 ) then
391    
392     if ( grdchk_eps .eq. 0. ) then
393     gfd = (fcpertplus-fcpertminus)
394     else
395     gfd = (fcpertplus-fcpertminus)
396     & /(grdchk_epsfac*grdchk_eps)
397     endif
398    
399     if ( adxxmemo .eq. 0. ) then
400     ratio_ad = abs( adxxmemo - gfd )
401     else
402     ratio_ad = 1. - gfd/adxxmemo
403     endif
404    
405     if ( ftlxxmemo .eq. 0. ) then
406     ratio_ftl = abs( ftlxxmemo - gfd )
407     else
408     ratio_ftl = 1. - gfd/ftlxxmemo
409     endif
410    
411     tmpplot1(itilepos,jtilepos,layer,itile,jtile)
412     & = gfd
413     tmpplot2(itilepos,jtilepos,layer,itile,jtile)
414     & = ratio_ad
415     tmpplot3(itilepos,jtilepos,layer,itile,jtile)
416     & = ratio_ftl
417    
418     fcrmem ( ichknum ) = fcref
419     fcppmem ( ichknum ) = fcpertplus
420     fcpmmem ( ichknum ) = fcpertminus
421     xxmemref ( ichknum ) = xxmemo_ref
422     xxmempert ( ichknum ) = xxmemo_pert
423     gfdmem ( ichknum ) = gfd
424     adxxmem ( ichknum ) = adxxmemo
425     ftlxxmem ( ichknum ) = ftlxxmemo
426     ratioadmem ( ichknum ) = ratio_ad
427     ratioftlmem ( ichknum ) = ratio_ftl
428    
429     irecmem ( ichknum ) = icvrec
430     bimem ( ichknum ) = itile
431     bjmem ( ichknum ) = jtile
432     ilocmem ( ichknum ) = itilepos
433     jlocmem ( ichknum ) = jtilepos
434     klocmem ( ichknum ) = layer
435     iobcsmem ( ichknum ) = obcspos
436     icompmem ( ichknum ) = icomp
437     ichkmem ( ichknum ) = ichknum
438     itestmem ( ichknum ) = itest
439     ierrmem ( ichknum ) = ierr
440     icglomem ( ichknum ) = icglo
441    
442     print *, 'grad-res -------------------------------'
443     print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',
444     & myprocid,ichknum,itilepos,jtilepos,
445     & layer,itile,jtile,obcspos,
446     & fcref, fcpertplus, fcpertminus
447     #ifdef ALLOW_TANGENTLINEAR_RUN
448     print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',
449     & myprocid,ichknum,ichkmem(ichknum),
450     & icompmem(ichknum),itestmem(ichknum),
451     & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
452     & ftlxxmemo, gfd, ratio_ftl
453     WRITE(msgBuf,'(A34,2(1PE24.14,X))')
454     & 'precision_grdchk_result TLM ', fcref, ftlxxmemo
455     CALL PRINT_MESSAGE
456     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
457     #else
458     print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',
459     & myprocid,ichknum,ichkmem(ichknum),
460     & icompmem(ichknum),itestmem(ichknum),
461     & bimem(ichknum),bjmem(ichknum),obcspos,
462     & adxxmemo, gfd, ratio_ad
463     WRITE(msgBuf,'(A34,2(1PE24.14,X))')
464     & 'precision_grdchk_result ADM ', fcref, adxxmemo
465     CALL PRINT_MESSAGE
466     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
467     #endif
468    
469     endif
470    
471     print *, 'ph-grd ierr ---------------------------'
472     print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
473     & ', ichknum = ', ichknum
474    
475     _BARRIER
476    
477     c-- else of if ( ichknum ...
478     else
479     ierr_grdchk = -1
480    
481     c-- end of if ( ichknum ...
482     endif
483    
484     c-- end of do icomp = ...
485     enddo
486    
487     if ( myProcId .EQ. grdchkwhichproc ) then
488     CALL WRITE_REC_XYZ_RL(
489     & 'grd_findiff' , tmpplot1, 1, 0, myThid)
490     CALL WRITE_REC_XYZ_RL(
491     & 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
492     CALL WRITE_REC_XYZ_RL(
493     & 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
494     endif
495    
496     c-- Everyone has to wait for the component to be reset.
497     _BARRIER
498    
499     c-- Print the results of the gradient check.
500     call grdchk_print( ichknum, ierr_grdchk, mythid )
501    
502     #endif /* ALLOW_GRDCHK */
503    
504     end

  ViewVC Help
Powered by ViewVC 1.1.22