Dot product in coefficient space

[Neil Lewis]

Hi There, 

I'm trying to compute the dot product in coefficient space, so that I may output KE spectra. My data is defined on the sphere (a spherical shell). To debug this issue, I'm using Ntheta=16 and Nr=5. 

To do this, I have tried defining the substitutions:

KE = u['c'] @ np.conj(u['c']) (1)

and 

KE = d3.dot(u(r=Ri).evaluate()['c'], np.conj(u(r=Ri).evaluate()['c'])) (2)

based on the answers to this question: https://groups.google.com/g/dedalus-users/c/tFsMV-NEUvo/m/IQv1MEBoCwAJ

However, when I try approach (1) I get the following error: 

Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 15 is different from 5)

and approach (2) yields: 

Both arguments to DotProduct must be Operand

Before I devote a large amount of time to resolving these errors, I wanted to check: can you take the dot product of two vectors in coefficient space? 

I had taken this to be implied by the answer linked to above, but now I'm not sure you can, because, e.g., (1) is so simple that I can't see what I could have done wrong. 

Thanks in advance, 

Neil Lewis 

[Jeffrey S. Oishi]

Hi Neil,

The dot product is a linear operation, so it doesn't matter where you
take it. I'd suggest using something like

KE = u@u

and then asking for KE['c'].

Jeff

> --
> You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/a6ee4721-3769-452f-b9ac-96685723a26bn%40googlegroups.com.

[Neil Lewis]

Hi Jeff, 

I think this will give me what I want, based on my own prior understanding, and additionally the discussion in the previous conversation I linked. 

See, e.g., this message from Ben Brown (copied verbatim from previous discussion):

Ruben,

     re energy spectra: Evan's exactly right (and Daniel wrote to me about the same off-list):

it depends on what you want.  If you want the spectrum of the kinetic energy, then KE['c'] gets you that.  

But if you want the energy spectrum (e.g., as used in turbulence theory), then what you want is dot(u['c'], np.conj(u['c'])), which is not the same.

From Pope's book Turbulent Flows:

The energy spectrum function E(k,t) is:

E(k,t) = int[Φ_ii(k, t)*δ(|k| - k)dk]   (eq 3.166)

with Φ_ii the velocity spectrum tensor, which is given by the fourier transform of the two point correlation.  Since the two point correlation is itself computable from FFTs, then I'm pretty sure

Φ_ii ~ dot(u_i['c'], np.conj(u_i['c'])) (eq 3.160 & 3.163)

up to constant factors, but you should check to make sure everything lines up with what you're comparing against.

Also, in case we forgot to say it elsewhere, one needs to be real careful when computing power spectra where one direction is represented on chebyshevs; coefficients in chebyshevs do not translate in the same way as in spherical harmonics or fourier.  This is why most spherical shell spectra are horizontal spectra only, evaluated at some specified shell depth (like mid-shell).  To do something like this, do:

E = dot(u(r=r_i).evaluate()['c'], np.conj(u(r=r_i).evaluate()['c']))

where r_i = 0.5 or whatever depth you want to evaluate at.

I think this is correct, but haven't tested it.

Also, poloidal/toroidal decompositions are a completely different equation formulation, so the power spectrum of e.g., W and Z from the Magic code's approach are going to be substantially different from your results if you're directly solving for u.  It will take care to convert back and forth between those and the velocity components.  I would probably choose to obtain the velocity components u from your comparison code,  calculate the energy spectrum directly yourself, and compare to that.

--Ben

[Neil Lewis]

Sorry, above I meant to say 'this will not give me what I want'. 

[Jeffrey S. Oishi]

Quite. I didn't read carefully enough. I can't test it right now, but
try using the numpy dot product on your second formulation (rather
than d3.dot).

> To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/2302995e-ee9b-4cf5-b68a-a88caca8d6e8n%40googlegroups.com.

[Neil Lewis]

With the np.dot in the second formulation I get a different error (similar to the error obtained when using @ in the first formulation):

ValueError: shapes (3,16,15,1) and (3,16,15,1) not aligned: 1 (dim 3) != 15 (dim 2)

More generally, might it be the case that one can't take the dot product in coefficient space? 

If this is the case, then I'm wondering if anyone has some advice on how to practically implement Ben's suggestions in the previous thread? While waiting, I'll try to figure this out for myself, and will post back here if I get something that works. 

Thanks, 

Neil 

[Daniel Lecoanet]

Hi Neil,

You don’t want to do a dot product, you want to do a multiplication.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/706c06fa-f509-4071-bd48-817f99c78a4an%40googlegroups.com.

[Neil Lewis]

Hi Daniel, 

By multiplication do you mean with '*'? If Y, then wouldn't that leave me with a vector, whereas KE should be a scalar with dims (horiz_coeff1, horiz_coeff2)? 

Thanks, 

Neil 

[Daniel Lecoanet]

Hi Neil,

After you use [‘c’] everything is in terms of numpy array. When you try to take a dot product, that is a dot product between numpy arrays which is between the last axis of the first array and the second-to-last axis of the second array (https://numpy.org/doc/stable/reference/generated/numpy.dot.html). Here if you take a multiplication you get an array of size (3, 16, 15, 1). If you want to sum over the three components, do np.sum(…. , axis=0).

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/fe39181c-d3cd-409c-8911-d38925fa9632n%40googlegroups.com.