py_helmholtz

Compute streamfunction, velocity potential and helmholtz decomposition from non-global wind data.

Implemented following:

Li, Zhijin and Chao, Yi and McWilliams, James C., 2006: Computation of the Streamfunction and Velocity Potential for Limited and Irregular Domains. Monthly Weather Review, 3384-3394.

1. To solve for streamfunction (psi) and velocity potential (chi) from u- and v- winds, on a uniform grid (uniform dx and dy everywhere):

u and v are given in even grid (n x m).

streamfunction (psi) and velocity potential (chi) are defined on a dual grid ((n+1) x (m+1)), where psi and chi are defined on the 4 corners of u and v.

Define:

u = u_chi + u_psi
v = v_chi + v_psi

u_psi = -dpsi/dy
v_psi = dpsi/dx
u_chi = dchi/dx
v_chi = dchi/dy

Define 2 2x2 kernels:

k_x = |-0.5 0.5|
      |-0.5 0.5| / dx

k_y = |-0.5 -0.5|
      |0.5   0.5| / dy

Then u_chi = chi \bigotimes k_x where \bigotimes is cross-correlation,

Similarly:

v_chi = chi \bigotimes k_y
u_psi = psi \bigotimes -k_y
v_psi = psi \bigotimes k_x

Define cost function J = (uhat - u)**2 + (vhat - v)**2

Gradients of chi and psi:

dJ/dchi = (uhat - u) du_chi/dchi + (vhat - v) dv_chi/dchi
dJ/dpsi = (uhat - u) du_psi/dpsi + (vhat - v) dv_psi/dpsi

(uhat - u) du_chi/dchi = (uhat - u) \bigotimes Rot180(k_x) = (uhat - u) \bigotimes -k_x
(vhat - v) dv_chi/dchi = (vhat - v) \bigotimes Rot180(k_y) = (vhat - v) \bigotimes -k_y
(uhat - u) du_psi/dpsi = (uhat - u) \bigotimes k_y
(vhat - v) dv_psi/dpsi = (vhat - v) \bigotimes Rot180(k_x) = (vhat - v) \bigotimes -k_x

Rot180() is a 180 degree rotation, and for k_x and k_y, it's the same as the inverting the sign. This Rot180() process is similar as the error back-propagation process in the convolutional neural network training process.

Add the regularization term:

J = (uhat - u)**2 + (vhat - v)**2 + lambda(chi**2 + psi**2)

2. To solve for streamfunction and velocity potential from u- and v- winds on irregular grid (e.g. mercator):

Use similar definition of cost function and gradients, except that the computation of component winds and derivatives are performed on steps for NE, NW, SE, SW qudarants:

u_chi =  0.5*((vp[:-1,1:]-vp[:-1,:-1])/dx_n + (vp[1:,1:]-vp[1:,:-1])/dx_s)
v_chi =  0.5*((vp[1:,:-1]-vp[:-1,:-1])/dy_w + (vp[1:,1:]-vp[:-1,1:])/dy_e)
u_psi = -0.5*((sf[1:,:-1]-sf[:-1,:-1])/dy_w + (sf[1:,1:]-sf[:-1,1:])/dy_e)
v_psi =  0.5*((sf[:-1,1:]-sf[:-1,:-1])/dx_n + (sf[1:,1:]-sf[1:,:-1])/dx_s)

du_chi/dchi, dv_chi/dchi, du_psi/dpsi and dv_psi/dpsi are also compose
by 4 quadrants, see code for details.

3. The Wind2D class is designed for computation of netcdf data via the CDAT interface.

CDAT: https://cdat.llnl.gov/index.html

4. Example

# read in some wind data as `u` and `v` via cdms

u=u(latitude=(5,50),longitude=(100,180))
v=v(latitude=(5,50),longitude=(100,180))

# create an wind obj, optimization will use gradient descent
w1=Wind2D(u,v,'GD')

# compute streamfunction and velocity potential
sf1,vp1=w1.getSFVP()

# get irrotational and non-divergent components
uchi1,vchi1,upsi1,vpsi1=w1.helmholtz()
uhat1=uchi1+upsi1
vhat1=vchi1+vpsi1

# create an wind obj, optimization will use scipy.optimize (recommended)
w2=Wind2D(u,v,'optimize')
sf2,vp2=w2.getSFVP()
uchi2,vchi2,upsi2,vpsi2=w2.helmholtz()
uhat2=uchi2+upsi2
vhat2=vchi2+vpsi2

# recompute and interpolate to the same grid of u and v
sf3,vp3=w2.getSFVP(interp=True)