Application of MatrixFree Methods in Ocean and Climate Modeling

1
Application of Matrix-Free Methods in Ocean and
Climate Modeling

• Samar Khatiwala
• Columbia University

Collaborators David Keyes (CU/APAM), Tim Kelley
(NCSU/Math), Tim Merlis
(CU undergrad), Mark Cane (CU/LDEO)
2
Solving A x b

• 2 classes of methods to do this
• Direct methods
• - storage and direct manipulation of the
  elements of A
• - exact solution
• - good for small, dense problems
• Iterative methods
• - approximate solution
• - good for large, sparse problems
• - KRYLOV SUBSPACE METHODS
• - DONT NEED TO KNOW A
• - ONLY REQUIRE ABILITY TO COMPUTE
• MATRIX-VECTOR PRODUCTS, i.e., Ax

Ax
Matrix-Free
x
yAx
3
Some applications

• Finding fixed points and limit cycles of large
  dynamical systems (e.g., ocean GCM)
• Linear stability and bifurcation analysis
• Parameter sensitivity
• Obtaining compact/reduced dimension
  representations of complex models
• Probing and analyzing complex systems defined
  through a (legacy) time-stepper
• Empirical parameterization of unresolved physics
  (e.g., spatial variations in eddy-diffusivity)

4
Efficient computation of equilibrium solutions of
seasonally forced ocean GCMs

• Ocean and climate models driven by
  time-independent or periodic (seasonal) forcing
• Very desirable to obtain equilibrium solutions
  (e.g., for initializing climate change
  simulations, parameter sensitivity studies, etc)

5
Current practice

• Direct forward integration until transients have
  died out
• If the transients disappear after a few periods,
  this is the simplest and most efficient approach
• HOWEVER
• Time scale for dynamical adjustment of the ocean
  is very long (several 1000s of years)
• Serious computational challenge in climate
  modeling

6
Matrix-Free Newton-Krylov Approach

• GCM time stepper
• u(t) F(u(0),t)

F RN RN
state transition function
model state at time t
N O(106)
Steady-state formulation F(u) u - F(u,T) 0
- Nonlinear algebraic system - T is the (known)
forcing period - System is nonautonomous so
solution is unique
7
How to solve F(u)0?

• Newtons method
• Starting with initial guess u0
• repeat for k0,1,
• (1) solve J(uk)duk -F(uk), J?F/?u
• (2) uk1 uk duk

Linear system of equations!

• - J is a DENSE matrix, impossible to store or
  compute
• Use a Krylov method to solve the linear system
• Only need to be able to evaluate F, i.e.,
• integrate the GCM

8
Initial results from an ocean GCM

• MIT OGCM, sector configuration (2 resolution)
• Seasonally forced
• N 40000

Recomputation of approximate preconditioner
9
Key challenge Matrix-free preconditioning

• PC via time stepper map which has better
  preconditioning than original problem (EV
  clustering)
• PC using sparse representation of model
  advection-diffusion operator and fast (Krylov)
  methods for matrix exponential calculation
• AD implicit time stepping of coarse grained
  problem
• Graph coloring methods for efficient Jacobian
  evaluation
• Multigrid, coarse-graining
• Problem scaling, globalization of Newton

10
Linear stability analysis of a global ocean GCM
MIT GCM Wrapper ARPACK
11
Temperature eigenvectors associated with
smallest Eigenvalue (gravest eigenmode)