/[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.31 - (show annotations) (download)
Tue May 24 22:40:30 2011 UTC (12 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62z, checkpoint62y, checkpoint63, checkpoint63a
Changes since 1.30: +29 -30 lines
include "GRDCHK_OPTIONS.h" instead of "CPP_OPTIONS.h" & "AD_CONFIG.h"

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

  ViewVC Help
Powered by ViewVC 1.1.22