/[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.9 - (show annotations) (download)
Thu Oct 2 21:30:22 2003 UTC (20 years, 8 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint51l_post, checkpoint51j_post, checkpoint51l_pre, checkpoint51h_pre, checkpoint51i_post, checkpoint51i_pre, checkpoint51g_post, checkpoint51m_post
Branch point for: tg2-branch
Changes since 1.8: +23 -7 lines
provide ARPACK interface

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

  ViewVC Help
Powered by ViewVC 1.1.22