/[MITgcm]/MITgcm/pkg/grdchk/grdchk_main.F
ViewVC logotype

Annotation of /MITgcm/pkg/grdchk/grdchk_main.F

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


Revision 1.8 - (hide annotations) (download)
Tue Jun 24 16:08:45 2003 UTC (20 years, 11 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint51, checkpoint51f_post, checkpoint51d_post, checkpoint51b_pre, branchpoint-genmake2, checkpoint51b_post, checkpoint51c_post, checkpoint51e_post, checkpoint51f_pre, checkpoint51a_post
Branch point for: branch-genmake2
Changes since 1.7: +4 -2 lines
Merging for c51 vs. e34

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

  ViewVC Help
Powered by ViewVC 1.1.22