340 likes | 524 Views
Algorithms Overview. Hernan G. Arango Institute of Marine and Coastal Sciences Rutgers University. 2004 ROMS/TOMS European Workshop CNR-ISMAR, Venice, October 18-20. Outline. ROMS/TOMS algorithms status ROMS/TOMS future releases How does one build an adjoint model? Ensemble prediction
E N D
Algorithms Overview Hernan G. Arango Institute of Marine and Coastal Sciences Rutgers University 2004 ROMS/TOMS European Workshop CNR-ISMAR, Venice, October 18-20
Outline • ROMS/TOMS algorithms status • ROMS/TOMS future releases • How does one build an adjoint model? • Ensemble prediction • Variational data assimilation: • Strong constraint 4DVAR • Weak constraint 4DVAR • Final remarks
ROMS/TOMS 2.1 Features • Fasham model revisited (Fennel) • Bio-optical model (up to 84 components), EcoSim (Bissett) • New bottom boundary layer (Blaas); Fixed Styles and Glenn BBL • Sediment model revisited: stratigraphy with Nbed layers (Warner) • Momentum and tracer balances (Crowley) • Time-averaged quadratic terms: <uu>, <uv>, <vv>, <uT>, <vT> • Isobaric Lagrangian trajectories (Warner)
ROMS/TOMS 2.1 Features • Sequential and concurrent coupling with atmospheric models (Moore, Shaffer) • ESMF (initialize, run, finalize) • Atmospheric coupler: Modeling coupling toolkit (MCT, Argonne National Lab) and WRF I/O API • MPI communicator is split between atmosphere and ocean nodes
ROMS/TOMS 2.1 Fixed Bugs • Horizontal viscosity • Parallel periodic boundaries • Tiling in serial applications • Added river mass transport to DU_avg1 and DV_avg1 arrays • MPI parallel bug in restart of floats NetCDF
ROMS/TOMS 2.2 Features • Ice model • Nesting / composed grids • Parallel IO • Improvements to sediment model • Monotonic tracer advection
Serial Versus Parallel NetCDF NCSA IBM P690 Parallel Output Time (0.1 s) 16 Serial 16 Timestep (246 x 240 x 16) Serial Output Time (0.1 s) Parallel 16 16 (Yang, NCSA) Timestep (656 x 640 x 16)
Serial Versus Parallel NetCDF Serial 128 NCAR IBM SP Cluster (WinterHawk II) Output Time (0.1 s) 64 32 16 Timestep (656 x 640 x 16) 128 Output Time (0.1 s) Serial 128 Parallel Timestep (Yang, NCSA) (656 x 640 x 16)
Sediment Model New Features • Suspended-sediment stratification effects in wave boundary layer (neutral currently) • Mechanics for cohesive versus non-cohesive bottom sediments • Gravity-driven transport in bottom boundary layer • Aggregation / dissaggregation • Wetting / drying • Bioturbation in sediment layers • Bedload transport (with wave effects) • Radiation stresses
ROMS/TOMS Adjoint and Data Assimilation Team Hernan G. Arango Boon Chua Bruce D. Cornuelle Emanuele Di Lorenzo Arthur J. Miller Andrew M. Moore Julio Sheinbaum Rutgers University Oregon State University Scripps Institute of Oceanography Georgia Institute of Technology Scripps Institute of Oceanography University of Colorado CICESE
Objectives • To provide the ocean modeling community with analysis and prediction tools that are available in meteorology and Numerical Weather Prediction (NWP), using a community OGCM (ROMS/TOMS). • To build a Generalized Stability Analysis (GSA) platform: eigenmodes, optimal perturbations /singular vectors, forcing singularvectors, stochastic optimals, pseudospectra. • To build an ensemble prediction platform. • To build 4D variational assimilation platforms.
Let’s represent NL ROMS as: • The TL ROMS is derived by considering a small perturbation sto S. A first-order Taylor expansion yields: A is real, non-symmetricPropagator Matrix • The AD ROMS is derived by taking the inner-product with an arbitrary vector , where the inner-product defines an appropriate norm (L2-norm): Overview
How To Build an Adjoint • The ADM can be derived from: • Continuous equations • Discrete equations (A is symmetric; exact) • Hand-coded • Automatic differentiation adjoint compilers (TAMC) • The ADM operator relative to L2-norm can be computed by multiplying each line of the TLM code by the corresponding adjoint variable, and then differentiating with respect the TLM variable. • Use Geiring and Kaminski (1998) transpose TLM operators and recipes. • Non-differentiable algorithms (vertical mixing).
Nonlinear Model DO k=1,N DO i=Istr,Iend+1 FX(i)=0.25_r8*(diff2(i,itrc)+diff2(i-1,itrc))*pmon_u(i)* & (Hz(i,k)+Hz(i-1,k))* & (t(i,k,nrhs,itrc)-t(i-1,k,nrhs,itrc)) END DO DO i=Istr,Iend t(i,k,nnew,itrc)=t(i,k,nnew,itrc)+ & dt*pm(i)*pn(i)*(FX(i+1)-FX(i)) END DO END DO
Tangent Linear Model DO k=1,N DO i=Istr,Iend+1 !! FX(i)=0.25_r8*(diff2(i,itrc)+diff2(i-1,itrc))*pmon_u(i)* !! & (Hz(i,k)+Hz(i-1,k))* !! & (t(i,k,nrhs,itrc)-t(i-1,k,nrhs,itrc)) tl_FX(i)=0.25_r8*(diff2(i,itrc)+diff2(i-1,itrc))*pmon_u(i)* & ((tl_Hz(i,k)+tl_Hz(i-1,k))* & (t(i,k,nrhs,itrc)-t(i-1,k,nrhs,itrc))+ & (Hz(i,k)+Hz(i-1,k))* & (tl_t(i,k,nrhs,itrc) -tl_t(i-1,k,nrhs,itrc))) END DO DO i=Istr,Iend !! t(i,k,nnew,itrc)=t(i,k,nnew,itrc)+ !! & dt*pm(i)*pn(i)*(FX(i+1)-FX(i)) tl_t(i,k,nnew,itrc)=tl_t(i,k,nnew,itrc)+ & dt*pm(i)*pn(i)*(tl_FX(i+1)-tl_FX(i)) END DO END DO
Adjoint Model DO k=1,N DO i=Istr,Iend !! tl_t(i,k,nnew,itrc)=tl_t(i,k,nnew,itrc) + !! & dt*pm(i)*pn(i)*(tl_FX(i+1)-tl_FX(i)) adfac=dt*pm(i)*pn(i)*ad_t(i,k,nnew,itrc) ad_FX(i+1)=ad_FX(i+1)+adfac ad_FX(i )=ad_FX(i ) -adfac END DO DO i=Istr,Iend+1 !! tl_FX(i)=0.25_r8*(diff2(i,itrc)+diff2(i-1,itrc))* pmon_u(i)* !! & ((tl_Hz(i,k)+tl_Hz(i-1,k))* !! & (t(i,k,nrhs,itrc)-t(i-1,k,nrhs,itrc))+ !! & (Hz(i,k)+Hz(i-1,k))* !! & (tl_t(i,k,nrhs,itrc) -tl_t(i-1,k,nrhs,itrc))) adfac =0.25_r8*(diff2(i,itrc)+diff2(i-1,itrc))*pmon_u(i)*ad_FX(i) adfac1=adfac*(t(i,k,nrhs,itrc)-t(i-1,k,nrhs,itrc)) adfac2=adfac*(Hz(i,k)+Hz(i-1,k)) ad_Hz(i ,k)=ad_Hz(i ,k)+adfac1 ad_Hz(i-1,k)=ad_Hz(i-1,k)+adfac1 ad_t(i ,k,nrhs,itrc)=ad_t(i ,k,nrhs,itrc)+adfac2 ad_t(i-1,k,nrhs,itrc)=ad_t(i-1,k,nrhs,itrc) -adfac2 ad_FX(i) =0.0_r8 END DO END DO
Ensemble Prediction • Optimal perturbations / singular vectors and stochastic optimal can also be used to generate ensemble forecasts. • Perturbing the system along the most unstable directions of the state space yields information about the first and second moments of the probability density function (PDF): • ensemble mean • ensemble spread • Adjoint based perturbations excite the full spectrum
Ensemble Prediction For an appropriate forecast skill measure, s
where model, background, observations, inverse background error covariance, inverse observations error covariance • Model solution depends on initial conditions ( ), boundary conditions, and model parameters • Minimize J to produce a best fit between model and observations by adjusting initial conditions, and/or boundary conditions, and/or model parameters. Data Assimilation Overview • Cost Function:
We require the minimum of at which: , , , yielding • AT is the transpose of A, often called the adjoint operator. It can be shown that: The adjoint equation solution provides gradient information Minimization • Perfect model constrained minimization (Lagrange function):
RP: 4D Variational Data Assimilation Platforms (4DVAR) • Strong Constraint (S4DVAR) drivers: • Conventional S4DVAR: outer loop, NL, AD • Incremental S4DVAR: inner and outer loops, NL, TL, AD (Courtier et al., 1994) • Efficient Incremental S4DVAR (Weaver et al., 2003) • Weak Constraint (W4DVAR) - IOM • Indirect Representer Method: inner and outer loops, NL, TL, RP, AD (Egbert et al., 1994; Bennett et al, 1997)
Ipass=1 Outer Loop CALL initial CALL main3d NLM: compute model-observations misfit and cost function CALL ad_initial CALL ad_main3d ADM: compute cost function gradients CALL descent CALL wrt_ini Compute NLM initial conditions using first guess conjugate gradient step size Ipass=2 CALL initial CALL main3d NLM: compute change in cost function Compute NLM initial conditions using refined conjugate gradient step size CALL descent CALL wrt_ini “Conventional” S4DVAR
Incremental S4DVAR CALL initial CALL main3d NLM: compute basic state trajectory and extract model at observations locations Outer Loop Ipass=1 Inner Loop TLM: compute misfit cost function between model (NLM+TLM) and observations CALL tl_initial CALL tl_main3d CALL ad_initial CALL ad_main3d ADM: compute cost function gradients Compute TLM initial conditions using first guess conjugate gradient step size CALL descent CALL tl_wrt_ini Ipass=2 CALL tl_initial CALL tl_main3d TLM: compute change in cost function Compute TLM initial conditions using refined conjugate gradient step size CALL descent CALL tl_wrt_ini Compute NLM new initial conditions (NLM+TLM) CALL ini_adjust CALL wrt_ini
Efficient Incremental S4DVAR CALL initial CALL main3d NLM: compute basic state trajectory and extract model at observations locations Outer Loop CALL ad_initial CALL ad_main3d ADM: compute initial estimate of the gradient CALL ini_descent Initialize conjugate direction as the negative of the gradient (adjoint) solution Inner Loop RPM: compute misfit cost function between model (NLM+TLM) and observations CALL tl_initial CALL tl_main3d CALL ad_initial CALL ad_main3d ADM: compute cost function gradients Compute TLM initial conditions using conjugate gradient step size CALL descent CALL tl_wrt_ini Compute NLM new initial conditions (NLM+TLM) CALL ini_adjust CALL wrt_ini
W4DVAR, IOM nl_roms: compute basic state trajectory nl_roms < nl_roms.in iom_roms: compute first guess and misfitbetween observation and model iom_roms < iom_roms.in Outer Loop Inner Loop ad_roms < ad_roms.in Inner loop, backward (ad_roms) and forward (tl_roms) integrations to compute tl_roms < tl_roms.in IOM components ad_roms: backward integration to compute ad_roms < ad_roms.in iom_roms < iom_roms.in iom_roms: compute
Twin Experiments • Spin-up an idealized, wind-forced double-gyre for 50 years. • Basin dimensions: 1000x2000 km2 • Grid resolution: dx=dy=18.518 km (54x108x4) • Run equilibrium solution for another 5 days and extract observations (true state) daily for each state variable at every spatial grid point. • Initialize 4DVAR algorithms from rest and assimilate observations at day 1. • Force only with the adjoint misfit (model minus observations) terms.
Adjusted Minus Truth Solution Final Adjusted Initial Conditions Vbar Difference RMS = 7.995e-6 Free-surface and Currents Ubar Difference RMS = 1.690e-5 Free-surface Difference RMS = 1.568e-5 S4DVAR
Free-surface Difference Potential Temperature Difference Vbar Difference Free-surface and Currents S4DVAR 3D Double Gyre Adjusted Minus Truth Solution Final Adjusted Initial Conditions Model-Observation Misfit Cost Function Iteration
Adjusted Minus Truth Solution Free-surface and Currents Ubar Difference RMS = 2.960e-2 Vbar Difference RMS = 5.085e-2 Free-surface Difference RMS = 2.136e-3 IOM True Solution Final Adjusted Initial Conditions
Ongoing 4DVAR Applications • Southern California Bight (Cornuelle, Di Lorenzo, Miller) • U.S. East coast (Arango, Moore, Wilkin) • Intra-Americas Sea (Moore, Sheinbaum) • Gulf of Mexico (Moore, Sheinbaum) • East Australia Current (Arango, Wilkin) • Oregon coast (Durski)
Observation Types plus satellite data (SSH, SST) and radar
Timing considerations • SCB – 6 CPU minutes per simulation day per TLM/ADM call on a 833MHz Alpha (78x118x30). • GoM – 17 CPU minutes per simulation day per TLM/ADM call on a 833 MHz Alpha. • IAS – 15 CPU minutes per simulation day per TLM/ADM call on a 833 MHz Alpha. • NENA – 60 CPU minutes per simulation day per TLM/ADM call on a 833 MHz Alpha (384x128x30). • Data assimilation scaling factors: • S4DVAR = 2 • IS4DVAR = 3 inner, 0.5 outer • EIS4DVAR = 2 inner, 0.5 outer • W4DVAR = 2 inner, 2.5 outer.
Final Remarks • Maintenance of TLM, RPM, and ADM models. • Parallelization of TLM, RPM, and ADM models. • Modeling background error covariance. • Training and documentation.
Publications • Moore, A.M., H.G Arango, E. Di Lorenzo, B.D. Cornuelle, A.J. Miller and D. Neilson, 2004: A comprehensive ocean prediction and analysis system based on the tangent linear and adjoint of a regional ocean model, Ocean Modelling, 7, 227-258. http://marine.rutgers.edu/po/Papers/Moore_2004_om.pdf • Arango, H.G., Moore, A.M., E. Di Lorenzo, B.D. Cornuelle, A.J. Miller and D. Neilson, 2003:The ROMS Tangent Linear and Adjoint Models: A comprehensive ocean prediction and analysis system, Rutgers Tech. Report. http://marine.rutgers.edu/po/Papers/roms_adjoint.pdf