[MITgcm-devel] Fwd: TAF issue for "special" parameter lists
Patrick Heimbach
heimbach at MIT.EDU
Fri Feb 24 19:49:02 EST 2006
----- Forwarded message from heimbach at MIT.EDU -----
Date: Fri, 24 Feb 2006 17:57:15 -0500
From: Patrick Heimbach <heimbach at MIT.EDU>
Reply-To: heimbach at MIT.EDU
Subject: TAF issue for "special" parameter lists
To: Ralf Giering <Ralf.Giering at FastOpt.com>
Ralf,
some of my colleagues have recently introduced a construct
into the model which breaks the adjoint.
It's not very clean programming, but works well for them
in pure forward to be able to switch between different
variations of Adams-Bashforth.
The general structure is as following (the underlying code
is calls of S/R gad_calc_rhs by S/R calc_gt, S/R calc_gs
for the case NOT AB_3)
call foo( ..., var, var, ... )
subroutine foo( ..., var1, var2, ... )
u = f1(var1)
v = f2(var2)
var is pure input to S/R foo, not modified.
The reason for doing this is to be able to use the same
argument list for a variation of AB (the AB_3) for which
call gad_calc_rhs( ..., theta, theta, ... ) ! AB_2
is modified into
call gad_calc_rhs( ..., gTnm1, theta, ... ) ! AB_3
Sensitivities for this construct are not
properly accumulated, here's why:
call adfoo( ..., advar, advar, ... )
subroutine adfoo( ..., advar1, advar2, ...)
advar2 = advar2 + adf2(advar2)
advar1 = advar1 + adf1(advar1)
The link of accumulating advar1 and advar2 onto advar
(because they're the same) is missing; strictly:
advar = advar1 + advar2
I have actually hand-coded this variation, i.e.
introduce an auxiliary variable advarh in the call
and then accumulate advarh onto advar, something like
call adfoo( ..., advar, advarh, ... )
advar = advar + advarh
advarh = 0.
It works.
This seems to show a strategy for TAF to properly
handle the construct.
The effect is subtle, not visible at first (and in
gradient checks), but after 13yr integration, the "ECCO"
gradient is messed up.
I am appending the routines where this problem occurs:
calc_gt.F (has the call to "foo" = gad_calc_rhs with
the problematic parameter list)
gad_calc_rhs.F ( = "foo" where theta=Tracer, and theta=TracAB)
adgad_calc_rhs (just to show how adTracer and adTracAB
are treated separately)
adcalc_gt_broken (the code generated by TAF,
it has adtheta twice in param. list)
adcalc_gt_fixed (I introduced on auxiliary variable adthetah,
which after the call gets added onto adtheta).
Hope this is enough info.
I think you'll quickly recognize the general problem.
It's not urgent, we'll re-write with a workaround,
but I think, you can modify TAF in above way to
handle this, provided TAF is able to recognize that
theta is twice in the argument list
(which might not be easy/general though).
Cheers
-Patrick
--
--------------------------------------------------------
Patrick Heimbach Massachusetts Institute of Technology
FON: +1/617/253-5259 EAPS, Room 54-1518
FAX: +1/617/253-4464 77 Massachusetts Avenue
mailto:heimbach at mit.edu Cambridge MA 02139
http://www.mit.edu/~heimbach/ USA
----- End forwarded message -----
--------------------------------------------------------
Patrick Heimbach Massachusetts Institute of Technology
FON: +1/617/253-5259 EAPS, Room 54-1518
FAX: +1/617/253-4464 77 Massachusetts Avenue
mailto:heimbach at mit.edu Cambridge MA 02139
http://www.mit.edu/~heimbach/ USA
-------------- next part --------------
A non-text attachment was scrubbed...
Name: calc_gt.F
Type: text/x-fortran
Size: 7816 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20060224/d98284ef/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: gad_calc_rhs.F
Type: text/x-fortran
Size: 19584 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20060224/d98284ef/attachment-0001.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: adgad_calc_rhs.F
Type: text/x-fortran
Size: 12052 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20060224/d98284ef/attachment-0002.bin>
-------------- next part --------------
subroutine adcalc_gt( bi, bj, imin, imax, jmin, jmax, k, km1, kup,
$ kdown, xa, ya, utrans, vtrans, rtrans, rtranskp1, maskup,
$kappart, myiter, mythid, adutrans, advtrans, adrtrans,
$adrtranskp1, adkappart, adfvert )
...
C----------------------------------------------
C ROUTINE BODY
C----------------------------------------------
calcadvection = tempadvection .and. ( .not. tempmultidimadvec)
iternb = myiter
if (staggertimestep) then
iternb = myiter-1
endif
if (tempforcing .and. ( .not. forcing_in_ab)) then
call adexternal_forcing_t( bi,bj,k )
endif
if (adamsbashforthgt) then
call adadams_bashforth2( bi,bj,k,iternb,adgt,adgtnm1 )
endif
if (tempforcing .and. forcing_in_ab) then
call adexternal_forcing_t( bi,bj,k )
endif
call adgad_calc_rhs( bi,bj,imin,imax,jmin,jmax,k,km1,kup,kdown,xa,
$ya,utrans,vtrans,rtrans,rtranskp1,maskup,diffkht,diffk4t,kappart,
$theta,theta,gad_temperature,tempadvscheme,tempvertadvscheme,
$calcadvection,tempimplvertadv,mythid,adutrans,advtrans,adrtrans,
$adrtranskp1,adkappart,adtheta,adtheta,adfvert,adgt )
end
-------------- next part --------------
subroutine adcalc_gt( bi, bj, imin, imax, jmin, jmax, k, km1, kup,
$ kdown, xa, ya, utrans, vtrans, rtrans, rtranskp1, maskup,
$kappart, myiter, mythid, adutrans, advtrans, adrtrans,
$adrtranskp1, adkappart, adfvert )
...
C----------------------------------------------
C ROUTINE BODY
C----------------------------------------------
cph(
do j = 1-oly, sny+oly
do i = 1-olx, snx+olx
adthetah(i,j,k,bi,bj) = 0.D0
adthetah(i,j,km1,bi,bj) = 0.D0
end do
end do
cph)
calcadvection = tempadvection .and. ( .not. tempmultidimadvec)
iternb = myiter
if (staggertimestep) then
iternb = myiter-1
endif
if (tempforcing .and. ( .not. forcing_in_ab)) then
call adexternal_forcing_t( bi,bj,k )
endif
if (adamsbashforthgt) then
call adadams_bashforth2( bi,bj,k,iternb,adgt,adgtnm1 )
endif
if (tempforcing .and. forcing_in_ab) then
call adexternal_forcing_t( bi,bj,k )
endif
call adgad_calc_rhs( bi,bj,imin,imax,jmin,jmax,k,km1,kup,kdown,xa,
$ya,utrans,vtrans,rtrans,rtranskp1,maskup,diffkht,diffk4t,kappart,
$theta,theta,gad_temperature,tempadvscheme,tempvertadvscheme,
$calcadvection,tempimplvertadv,mythid,adutrans,advtrans,adrtrans,
$adrtranskp1,adkappart,adtheta,adthetah,adfvert,adgt )
cph(
do j = 1-oly, sny+oly
do i = 1-olx, snx+olx
adtheta(i,j,k,bi,bj) = adtheta(i,j,k,bi,bj) +
& adthetah(i,j,k,bi,bj)
adthetah(i,j,k,bi,bj) = 0.D0
adtheta(i,j,km1,bi,bj) = adtheta(i,j,km1,bi,bj) +
& adthetah(i,j,km1,bi,bj)
adthetah(i,j,km1,bi,bj) = 0.D0
end do
end do
cph)
end
More information about the MITgcm-devel
mailing list