Extracting poloidal-toroidal coefficients of vector fields directly in dedalus

[Andrey Kononov]

Hi all! 

I compute EVP and get a vector field as an eigenfunction. I want to decompose it into poloidal-toroidal amplitudes following this notation: 

Right now I am integrating over a sphere for different radii (code attached) to acquire those coefficients, but I was wandering whether a more direct approach exists. As I understand, in dedalus vector fields on balls/shells are decomposed into spin-weighted spherical harmonics, which are connected to poloidal-toroidal decomposition by simple identites. So perhaps I can get those amplitudes from xi['c'] or something alike? 

Thank you! 

[Semih Tuna]

Hello,

I recently needed this for another application. I haven't worked out the mapping from spin-weighted harmonics to poloidal-toroidal decomposition, because it was convenient to simply extract the these functions from radial components of the vector field and its divergence/curls. 

For example, to get W_lm(r), you can simply define

Nr, Ntheta, Nahi = ...

mesh = ..

comm = ..

dist = d3.Distributor(coords, dtype=dtype, mesh=mesh,

                      comm= comm)

er = dist.VectorField(coords, bases=ball.radial_basis)

er.change_scales(dealias)

er['g'][2] = 1.0

ksi_r = ksi@er

and look at coefficient space data of ksi_r. Since ksi_r is a Dedalus scalar, Dedalus simply gives you coefficients of an expansion in Y_lm, which are related to your W_lm(r). I managed to do this in a parallel-safe manner as below:

radius_arr= ball.radial_basis._radius_grid(dealias)

modes= np.zeros((len(radius_arr),

                  Ntheta,

                  2*Ntheta + 2))

for n, rad_val in enumerate(radius_arr):

    surf= ball.S2_basis(radius= rad_val)

    m_ind, l_ind, n_ind= dist.coeff_layout.local_group_arrays(

        surf.domain(dist),

        scales=1)

    ksi_r_arr= ksi_r(r=rad_val)['c'][:]

    for l in l_arr:

        for m in range(l + 1):

            msk = ((l_ind == l) & (m_ind == m))

            if not np.any(msk):

                continue

            modes[n, l, 2*m]= ksi_r_arr[msk][0]

            modes[n, l, 2*m + 1]= ksi_r_arr[msk][1]

modes array holds W_lm (or W_lm/r).

To get coefficients of Phi and Psi, you can define curl or divergence of ksi and make use of VSH identities. Whatever ends up in the radial component, you can extract its spherical harmonic expansion from coefficient space data ['c'] this way. 

This way, instead of (l x m) volume integrals, you do n= len(rad_arr) Dedalus interpolations, which I found to work faster in my application.

Hope it helps, and let me know if anyone finds a more practical way of doing this.

Semih