The rules for coding adjoint and the application of

Document Sample
The rules for coding adjoint and the application of Powered By Docstoc
					                       Lecture Two

The rules for coding adjoint and the application
of 3D variational method to radar data analysis

                         Jidong Gao
                       jdgao@ou.edu
                       November 2006

        Center for Analysis and Prediction of Storms

                  University of Oklahoma
              What is the adjoint?

•   It‟s only a mathematical tool to help you to get
    the gradient of the cost function.

•   In Errico (1997) “What is an adjoint model?”, It
    is said “…the adjoint is used as a tool for
    efficiently determining the optimal solutions.
    Without this tool, the (4D) optimization
    problem (including minimization and
    maximization) could not be solved in a
    reasonable time for application to real-time
    forecasting”. This is a good statement. Then
    what does it mean exactly?
                       The adjoint
•   For a 3D cost function defined as
           1                             1
        J  ( x  xb )T B 1 ( x  xb )  ( H ( x)  y o )T R 1 ( H ( x)  y o )
           2                             2
•   It‟s gradient w.r.t. x is
         x J  ( x  x b )  HT [ H ( x)  y ob ]
•   It involves the HT operator and associated calculation.
    This HT is the „adjoint‟ of H.
•   For a problem of small size, one can explicitly express
    H in a matrix form and perform explicit transpose
    operation and associated caculations. For a large
    problem, this is unrealistic and the sparseness of H
    also makes it unnecessary.
•   Here, the „adjoint‟ coding technique comes to the
    rescue.
•   Mathematical definition of adjoint:

        For any linear forward operator: X=LY, the adjoint code will be: Yad=LTX

•   A simple example for coding adjoint:

forward code:
        do i = 1, N-1
          x(i)=x(i)+a*y(i+1)
        end do
  input: xin, yin; output: xout

the adjoint code:
        do i = 1, N-1
          ad_y(i+1)= a*ad_x(i)
          ad_x(i) = ad_x(i)
        end do
  input: ad_xin,; output: ad_xout, ad_yout.

To verify the correctness of the adjoint code, you can calculate the following two terms:
           xin * ad_xout + yin *ad_yout
           xout *ad_xin
They should be equal.
y    1D variational interpolation scheme



                                        M=6
         1
                                                      x
                                             N=12




A problem: Suppose we have M irregular “observations” of
a sinusoidal function in y, we’d like to determine values at
regular grid points from 1 to N that best represent the
function. First guess: pink curve.
•     We can define a cost function,
                      1                        1
                   J  (T  T b )T (T  T b )  ( H (T )  T o )T ( H (T )  T o )
                      2                        2
    Subroutine FCN(N, T, F) ! subroutine calculating the cost function
    Implicit none
    Use module_const ! Include observational number m=6
    Use module_pass ! Use this to pass forcing term forc(m) to gradient sub.
    !!!! and observation value Tob(6), and its location y(6).
    integer::N, i, j; real::T(N), J, temp ! T is the state vector defined as location vector x
    J = 0.0
    do i=1,n
    J = J+0.5*(T(i)-Tb(i))**2
    enddo
    do j=1, m
      i = max(1, min(n, int (y(j)/dx)+1) )
      temp=T(i)*(x(i+1)-y(j))/dx+T(i+1)*(y(j)-x(i))/dx
      forc(j)= temp-Tob(j)
      J = J+0.5*forc(j) **2
    end do
    return
    end
     Calculation of cost function
              gradident
•    How to calculate the gradient?
    Its not straightforward, because of the reasons stated earlier, and
     because the cost function for obs. is the summation of innovation
     vector in obs. points, not in grid points. But when we do the
     minimization, it is the gradients of cost function with respect to
     the analysis variables at the grid points that are required. We can
     use the adjoint technique to solve this problem. The formulation
     for the gradient is,

            x J  ( x  x b )  HT [ H ( x)  y ob ]
The second term on rhs transfers the gradient of the cost function from the
obs points to the grid points. It can be coded using the adjoint technique.
“Recipes for Adjoint Construction” by Giering and Kaminski 1998) offers a
good tutorial. Linked to at the class web site.
•   Subroutine GRAD(N, t, G)
    Implicit none
•   Use module_const ! Include observational number m=6
•   Use module_pass ! Use this to pass forc(6) from sub FCN
•   !!!! and observation value Tob(6), and its location y(6).
•   integer::N, i, j; real::T(N), J, temp
     G(:) = 0.0
•   do j=1, m
•     i = max(1, min(n, int (y(j)/dx)+1) )
•    G(i) = G(i) + 0.5*forc(j)*(x(j+1)-y(j))/dx
•    G(i+1)=G(i+1)+0.5*forc(j)*(y(j)-x(j))/dx
•   end do

    do i=1,N
      G(i)=G(i)+0.5*(T(i)-Tb(i))
    end do

•   return; end
         storm-scale phenomenon

By storm-scale, we mean
structures like these

We will need to predict
the structure and
evolution of individual
cells

Spatial scale ~ few km
Time scale ~ tens of mins




          For storms, radar is about the only obs platform!
     Radar Data Assimilation

The Radar Problem:

Observe radial velocity and                           Real wind
reflectivity, other dual-
                                          v
                                                         Radial Vr
polarization variables.
                                                  

But NWP model, needs V, T,    y       r               u
                                  
P, qv, qr etc for IC.
                                      x



                              x   y ux  vy
    Vr  u cos   v sin   u  v 
                              r   r    r
NEXRAD Radar provides high spatial and temporal
resolution observations

Assimilation of these radar data into storm-resolving
models is very important for the future of operational
forecasting

Here, I will provide an example using the 3DVAR
method to do single-Doppler velocity retrieval
We are trying to retrieve 3-D wind vectors from single-
Doppler observed reflectivity and radial velocity during a
short time period. A cost-function is defined as follows:

  J=JB + JO + JE+ JD + JS

JB, background field constraint,
Jo, radial wind constraint,
JE, observed tracer constraint (reflectivity),

JD, 3-D anelastic mass continuity equation,
JS, smoothness constraints.
1) Background field constraint

         J   Wu(u  u )2   Wv(v  v )2
           B


            ijk        b     ijk       b

     background can be obtained from other data
     source. A nearby sounding is used here.


2)   Jo, measures the difference between the analyzed radial
          velocity and the observed radial velocity:

          J  1 Wr (Vr V n )2,
              2n
               O
                          rob
  Vr is obtained by interpolating retrieved variables u,, v,and w from
 the grid to observation points, then projecting the winds to the
 radial direction.
3) Observed tracer constraint

    J    = 1  W (E)2,
        E 2 ijkn E

  defines relationship of retrieved variables u, v, w and
  observed tracer h in the simplified equation,


         E  hob  u hob  v hob  w hob
              t       x       y       z
             k 2 hob  kv2 hob  Fm ,
               h H             Z


  Here u, v, and w are velocities to be retrieved.
4)     Anelastic mass continuity constraint

       J   WD(  u   v   w)2
                  x y z
         D


          ijk
      JD imposes a weak anelastic mass conservation on the
      analyzed winds.


5) Smoothness constraint

      J s   wsu(2u)2   wsv (2v)2   wsw(2w)2
           ijk           ijk            ijk
     The spatial smoothness term is used to remove small
     -scale noise and improve the quality of analysis.
          A single observation experiment, no
                   other constraint




horizontal slice at z=9km   Vertical slice at
u-v vector, contour of u    y=42km, w contour
         Same as above, adding JD, JE




horizontal slice at z=9km   Vertical slice at
u-v vector, contour of u    y=42km, w contour
Analysis of May 17, 1981 Arcadia
 Storm, with Cimarron Radar

  Z = 500 m            Z = 9 km




               An analysis with the radial
                velocity obs. only, no bkgd or
                other constraints. We get the
                results as expected

  Z = 4.5 km
Real Radar Data Analysis
Cimarron Radar, cont’d
    with a vertical slice at j = 25
Same as the above, but add JD, JE,
          but without JS



  Z = 500
  m                    Z = 9.0 km




    Z = 4.5 km
Same as above, but for vertical slice j = 25
with bkgd, but without smoothness

Z = 500 m                Z = 4.5 km




                J = 25
16:34 CST, May 17th, 1981, Arcadia, OK storm
        SDVR-Retrieved                  Dual-Doppler




    A                              A




                B                                B
                                 Cima
                                 OK
             Norman,OK                      Norman, OK

Wind vector at Z=500m. Shaded is observed reflectivity field.
      SDVR-Retrieved                Dual-Doppler




Vertical slice through line A-B in the horizontal slice.
                  Comments
•   Single obs. tests illustrates the behavior of
    3DVAR for the analysis of radial wind.

•   When applying 3dvar to radar data,
    smoothness constraints, or background error
    covariances, or filters must be used to
    properly spread the influence of data.

•   With all constraints used, we can obtain
    pretty good analysis of storm structure,
    otherwise we obtain only the radial wind
    component plus a uniform background field.
     Possible future research

–   radar data - essential for storm-scale
    NWP
–   available radar data:
       –   radial velocity
       –   reflectivity
       –   quantities from dual-polarization
           measurements, such as, differential
           reflectivity, differential phase shift, etc
       –   vertical wind profiles
       –   accumulated rainfall
Possible future research (con’d)


   thinning of too dense data consistently
    with the resolution of the analyses


   study of space- and time structure of
    biases between simulated and observed
    data for bias correction
Possible future research (con’d)

–   reflectivity is sensitive to cloud properties
      simulation depends on clouds
       representation in model
      Start from idealized study, to get simulated
       reflectivity from model itself
      study of existing radar simulation models
Possible future research (con’d)

–   code tangent linear and adjoint of
    observation operator
     functional 3DVar minimization
     check speed and quality

–   simulate and assimilate single radar data
     check how 3DVar corrects atmospheric
      fields
Possible future research (con’d)

–   study the impact of reflectivity assimilation to
    forecast
     for different situations (front, strong
      convection...)
     for different resolutions (10 km, 2.5 km)
     check spin up
Thank You for your attention!



     Questions? Send email to:
          jdgao@ou.edu
         or visit Rm 4110