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

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

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


Revision 1.28 - (show annotations) (download)
Fri May 29 06:14:47 2009 UTC (14 years, 11 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint62c, checkpoint62a, checkpoint62, checkpoint62b, checkpoint61q, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.27: +2 -1 lines
Bug fix (spotted by jmc): un-initialized var/

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

  ViewVC Help
Powered by ViewVC 1.1.22