Using a field that is non-constant in two directions for an EVP

[Joseph Hall]

I am attempting to make some additions to an EVP code in cylindrical geometry. The problem is decomposing waves in a cold ion-electron plasma into eigenmodes to search for a specific surface wave. Currently, we have non-dimensionalized the problem such that there is only one parameter that contains the density and distribution of the plasma in space, call it omega_pnsq. Right now, the problem is symmetric in theta, and the density of the plasma drops off sharply at r=~73, see this plot:

The idea is now to break this cylindrical symmetry by introducing a surface perturbation like this (exaggerated):

I have read on this group that Dedalus does not support multiple Chebyshev bases until v3, and that a coefficient can only be non-constant along the one Chebyshev basis. Hence, I have r coded as a Chebyshev and theta as a Fourier basis. Is there still any way to have this density parameter vary with both r and theta? I followed this and this thread but could not even get it to work when my field was only dependent on r. My code looks like

t_basis = de.Fourier('th', Nt, interval=(0, 2*np.pi))

r_basis = de.Chebyshev('r', Nr, interval=(40, 103))

domain = de.Domain([t_basis, r_basis], grid_dtype=np.complex128)

omega_pnsq = domain.new_field()

th, r = domain.grids(scales=1)

omega_pnsq['g'] = 0.5 * (np.tanh((73.35 - r)/6.2) + 1) # note that this only depends on r

problem.parameters['omega_pnsq'] = omega_pnsq

I then go on to define substitutions using omega_pnsq, which gets fed into what is essentially Maxwell's equations. The error message I get is "[long equation] is non-constant along separable direction 'th'." This problem disappears when I make omega_pnsq a substitution or specify that omega_pnsq.meta['th'] is constant but that's unhelpful since I eventually want it to vary with r and theta as a field. 

summary: Why does the code I posted raise this error, and is there any way to have a parameter/field vary across two bases before the next update?

[Daniel Michael Lecoanet]

Hi,

I think the issue is that you need to tell Dedalus that omega_pnsq is constant in the theta direction. I think the syntax for that is:

omega_pnsq.meta[’th’][‘constant’] = True

Daniel

On May 26, 2021, at 2:09 PM, Joseph Hall <josep...@brown.edu> wrote:

I am attempting to make some additions to an EVP code in cylindrical geometry. The problem is decomposing waves in a cold ion-electron plasma into eigenmodes to search for a specific surface wave. Currently, we have non-dimensionalized the problem such that there is only one parameter that contains the density and distribution of the plasma in space, call it omega_pnsq. Right now, the problem is symmetric in theta, and the density of the plasma drops off sharply at r=~73, see this plot:

<download.png>

The idea is now to break this cylindrical symmetry by introducing a surface perturbation like this (exaggerated):

<download-1.png>

I have read on this group that Dedalus does not support multiple Chebyshev bases until v3, and that a coefficient can only be non-constant along the one Chebyshev basis. Hence, I have r coded as a Chebyshev and theta as a Fourier basis. Is there still any way to have this density parameter vary with both r and theta? I followed this and this thread but could not even get it to work when my field was only dependent on r. My code looks like

t_basis = de.Fourier('th', Nt, interval=(0, 2*np.pi))

r_basis = de.Chebyshev('r', Nr, interval=(40, 103))

domain = de.Domain([t_basis, r_basis], grid_dtype=np.complex128)

omega_pnsq = domain.new_field()

th, r = domain.grids(scales=1)

omega_pnsq['g'] = 0.5 * (np.tanh((73.35 - r)/6.2) + 1) # note that this only depends on r

problem.parameters['omega_pnsq'] = omega_pnsq

I then go on to define substitutions using omega_pnsq, which gets fed into what is essentially Maxwell's equations. The error message I get is "[long equation] is non-constant along separable direction 'th'." This problem disappears when I make omega_pnsq a substitution or specify that omega_pnsq.meta['th'] is constant but that's unhelpful since I eventually want it to vary with r and theta as a field. 

summary: Why does the code I posted raise this error, and is there any way to have a parameter/field vary across two bases before the next update?

--
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/06fe567d-f291-4193-a2a1-192c7b5f2902n%40googlegroups.com.
<download-1.png><download.png>

[Joseph Hall]

That is definitely the reason why it was throwing an error when I used the original density profile. However, won't making omega_pnsq.meta[’th’] a constant negate any asymmetry (which is what I'm trying to investigate)? A quick test shows that the results don't change with meta[’th’] constant, which makes sense. Maybe this can't be done in an EVP but even if it's possible in an IVP it would be useful information.

[Bindesh Tripathi]

(Thanks, Daniel, for copying me on this.)

Hi Joseph,

I recently leveraged Dedalus to perform a 2D EVP calculations -- although it required to write down many equations (on the order of 50-100), coupled along the Fourier direction. I believe, the current implementation of Dedalus does not do this job for you automatically. So, we do have to combine a bit of analytical work with numerics. The idea basically is to create a 1D EVP but promote the state vector to a higher dimensional one where the state vector is X = [ f_m(r) for m in range(1, N) ] and f_m(r) is the eigenfunction at Fourier mode number m and N is the maximum number of Fourier modes that you wish to couple (this coupling is necessary when you Fourier transform analytically the product of background term that depends both on theta and r and the other term which is your desired state vector.)

For example, if you do not have a theta dependence on the background field, then the perturbations written as f_m(r)*exp(i*m*theta)*exp(i*lambda*t)* immediately yields a 1D-EVP:

dt(M(f_m)) = L(f_m) where f_m is the eigenfunction of the linear operator L (that can have terms like d/dr and radial dependent terms), M is another linear operator. This is straightforward.

__________________________________________________

Now, bring in the background field that has both theta and r dependence. Linearizing your equations to keep first order terms in state vectors (i.e., f_m's), you will this time have coupling between different f_m's (m is the Fourier mode number) as the background field (let's call it B, for brevity and background, instead of omega_pnsq) also has m dependance. You will probably have something like:

dt(M(f_m)) = L(f_m)  + \sum_j {A_j * f_{m-j}} where A_j is another linear operator (some algebraic function of the background B; here j in A_j represents the jth Fourier component of A (which is eventually of B)). Please note that this is still a linear equation, but the only complexity is that it couples many f_m's together (thus you would get the 2D eigenmodes later on; think of them as different modes of vibration of a 2D membrane like that of a drum).

Set you background fields in different variables named as B_ij where ij represents the Fourier mode number:

for jj in range(0, N+1):

    nccvar = domain.new_field()

    nccvar.set_scales('1', keep_data=False)

    nccvar['c'] = feed_your_background_data_for_m_equals_to_jj_but_only_the_radial_dependence

    problem.parameters['B_jj'] = ncc_var

Now, the task is to write equations like below (which can be done via a python script to write strings that represent all of your coupled but still linear equations and then pass them onto one of the greatest friends -- Dedalus :) ).

for m in range(1,N):

    problem.add_equation("dt(M(f_m)) - L(f_m) - \sum_j {A_j* f_{m-j}} = 0")       ...(1)

    problem.add_bc("write your BCs for each f_m")

Now, build an eigenvalue solver and proceed with the usual 1D EVP. The only difference between the usual 1D EVP (where the background does not depend on theta) and 2D (as done here) is that the earlier one has problem vars = [f_{one_particular_m_only}] and the latter one has problem_vars = [f_1, f_2, ..., f_N]. We still employ the same 1D Chebyshev basis, in this method. Maybe Dedalus v3 has interesting features to deal with such problems.

**Note: Since the coupling between the background and the state vectors brings in the full summation in j (in equation 1 above), you will end up with one particular annoying term: A_m* f_0, which unless f_0=0, you would have to get rid of this term by solving an associated linear BVP as discussed briefly in https://groups.google.com/g/dedalus-users/c/PeQLSCUArRU**

BTW, I noticed that you have created 2D domain and I was wondering why you have np.complex128 for the eigenfunctions (I am assuming you are working with real observables). Since the method employed here is 1D EVP but for the appended state vectors, i.e., the new state vector =  X = [f_1(r), f_2(r), ...f_N(r)], there should be 1D domain and thus it is now required to have np.complex128 data type for the eigenfunctions as you have imposed as they are Fourier-transformed materials.

Hope this helps. Please let us know if you happen to have some more queries.

Thanks a lot,

Bindesh

UW-Madison

 

[Daniel Michael Lecoanet]

Hi,

Yes, Bindesh is right that there is no native support for 2D EVPs in Dedalus. Writing out multiple equations is one way to deal with this. However, if you’re just interested in the fastest growing eigenmode (or the most slowly decaying eigenmode) then you could solve the linearized 2D problem as an IVP. Provided there is a unique fastest growing mode, you should be able to measure the growth rate and eigenfunction from the IVP. I like to call this the power method.

Daniel

On May 26, 2021, at 7:37 PM, 'Bindesh Tripathi' via Dedalus Users <dedalu...@googlegroups.com> wrote:

(Thanks, Daniel, for copying me on this.)

Hi Joseph,

I recently leveraged Dedalus to perform a 2D EVP calculations -- although it required to write down many equations (on the order of 50-100), coupled along the Fourier direction. I believe, the current implementation of Dedalus does not do this job for you automatically. So, we do have to combine a bit of analytical work with numerics. The idea basically is to create a 1D EVP but promote the state vector to a higher dimensional one where the state vector is X = [ f_m(r) for m in range(1, N) ] and f_m(r) is the eigenfunction at Fourier mode number m and N is the maximum number of Fourier modes that you wish to couple (this coupling is necessary when you Fourier transform analytically the product of background term that depends both on theta and r and the other term which is your desired state vector.)

For example, if you do not have a theta dependence on the background field, then the perturbations written as f_m(r)*exp(i*m*theta)*exp(i*lambda*t)* immediately yields a 1D-EVP:

dt(M(f_m)) = L(f_m) where f_m is the eigenfunction of the linear operator L (that can have terms like d/dr and radial dependent terms), M is another linear operator. This is straightforward.

__________________________________________________

Now, bring in the background field that has both theta and r dependence. Linearizing your equations to keep first order terms in state vectors (i.e., f_m's), you will this time have coupling between different f_m's (m is the Fourier mode number) as the background field (let's call it B, for brevity and background, instead of omega_pnsq) also has m dependance. You will probably have something like:

dt(M(f_m)) = L(f_m)  + \sum_j {A_j * f_{m-j}} where A_j is another linear operator (some algebraic function of the background B; here j in A_j represents the jth Fourier component of A (which is eventually of B)). Please note that this is still a linear equation, but the only complexity is that it couples many f_m's together (thus you would get the 2D eigenmodes later on; think of them as different modes of vibration of a 2D membrane like that of a drum).

Set you background fields in different variables named as B_ij where ij represents the Fourier mode number:

for jj in range(0, N+1):

    nccvar = domain.new_field()

    nccvar.set_scales('1', keep_data=False)

    nccvar['c'] = feed_your_background_data_for_m_equals_to_jj_but_only_the_radial_dependence

    problem.parameters['B_jj'] = ncc_var

Now, the task is to write equations like below (which can be done via a python script to write strings that represent all of your coupled but still linear equations and then pass them onto one of the greatest friends -- Dedalus :) ).

for m in range(1,N):

    problem.add_equation("dt(M(f_m)) - L(f_m) - \sum_j {A_j* f_{m-j}} = 0")       ...(1)

    problem.add_bc("write your BCs for each f_m")

Now, build an eigenvalue solver and proceed with the usual 1D EVP. The only difference between the usual 1D EVP (where the background does not depend on theta) and 2D (as done here) is that the earlier one has problem vars = [f_{one_particular_m_only}] and the latter one has problem_vars = [f_1, f_2, ..., f_N].We still employ the same 1D Chebyshev basis, in this method. Maybe Dedalus v3 has interesting features to deal with such problems.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/13e203f7-38af-4818-8f13-7a05c8fa134dn%40googlegroups.com.

[Ben Brown]

Following Daniel's note:

      for IVPs, put the term that's non-constant in 2-D (or more) on the RHS.

Everything is fine and fair there, no concerns about coupled directions etc.

You can do fancier things by splitting into a chebyshev-direction variation (LHS) and an all-the-rest variation (RHS) but whether or not that does anything for you depends on your problem.

See Eric Hester's (& friends) work on immersed boundary conditions for a logical limit (in my mind) of these techniques.

--Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/A4125457-D928-4F52-B296-2180AACE5984%40northwestern.edu.

[Joseph Hall]

Thank you everyone for all the information and support! It seems like the best thing to do will be keep working on the IVP and return to the EVP later, but I will update if anything comes up.

- Joe

[Joseph Hall]

This is a pretty small question so I didn't want to start a new thread, but what does the method .T do? For example in the TC flow notebook where they plot the output:

pcolormesh((r[0]*ones_like(z)).T, (z*ones_like(r)).T, v['g'].T, cmap='PuOr') 

quiver((r[0]*ones_like(z)).T, (z*ones_like(r)).T,u['g'].T, w['g'].T, width=0.005)

When I call my grid fields with this argument, the output is qualitatively different and yet I can't seem to find any forum post or documentation on what this method actually is. Out of curiosity, what does it do? The only thing that comes to mind is transpose, but I can't quite tell if that's true from my own code.

[Bindesh Tripathi]

It should transpose, as you correctly mentioned.

[Ben Brown]

Yep.  It’s the numpy array intrinsic. 

—Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/9291415a-0c05-4d5b-b888-26f07bce8cd4n%40googlegroups.com.