eigenvalue pbm using Chebyshev modes + mapping

[louisalexan...@gmail.com]

Hello Dedalus users,

I'm trying to understand if I can solve an eigenvalue problem using Dedalus.

Basically, I have an ODE of the form:

(1) L_2(\omega,r) d^2/dr^2 f + L_1(\omega,r) d/dr f + L_0(\omega,r) f = 0

that I want to solve for \omega (the eigenvalues) with r \in ]-\infty,\infty[. The L_0,1,2 are real functions of \omega and r (and other parameters).

I'd like to use the mapping:

(2) r -> H s/(1-s^2)

with H a "spreading" parameter. It converts s \in [-1,1] into r \in ]-\infty,\infty[ with a nice clustering near r=0 (point of interest) if you use the classical chebyshev collocation points.

Do you think I can do this with Dedalus?

The only caveat I see is that there is a problem with the mapping (2) at s=-1,+1 where it is singular. So when transforming (1) in s space, I introduce some terms that are singular at s=-1/+1.

If I were to write my own Cheby solver, I think that these terms could be dealt with by using null BCs (Dirichlet) and removing the problematic lines in the Jacobian*Derivative operator matrix of the eigenvalue problem.

But here, Dedalus (nicely) handles everything. Do you see a simple way to do this manually, i.e. go into Dedalus and tell it to remove the problematic lines in the matricial system?

Thanks!
Louis

[Ben Brown]

Louis,

     Boyd has a bunch written on these methods in chapter 17 of "Chebyshev & Fourier Spectral Methods".  I haven't tried these myself yet, but since we don't include +-1 (we're on an interior grid), things might just work.

--Ben

Louis





--
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-users+unsubscribe@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/0e3d7995-1c64-4c8a-94b1-e94ffacb4113%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Keaton Burns]

Hi Louis,

Following up, here’s an example notebook that computes the eigenvalues of the quantum harmonic oscillator over the whole real line using the algebraic mapping described in Boyd, i.e. 

x = H s / sqrt(1 - s^2)

You could use Dedalus substitutions to implement the change of variables, but you may be better off writing out the equations in terms of s and multiplying through by some factors to clear the denominators and give simpler polynomial NCCs.

Best,

-Keaton

To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.

To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAHqBLzybEfuAQw%3D41ArpQOGzZhJJTRth0n%3DHjkzG1_6muPkVEg%40mail.gmail.com.

[louisalexan...@gmail.com]

Hi Keaton,

Your example notebook is great. Exactly what I was looking for.

However, running it in Jupyter notebook I get the error:

"AttributeError: 'EigenvalueSolver' object has no attribute 'solve_sparse'"

I should say that I haven't used Jupyter very much so I may just be missing something trivial. Maybe it's because my Dedalus install is 1.5 years old?

I just checked the class EigenvalueSolver in dedalus/core/solvers.py and it doesn't seem to include the solve_sparse def.

Thanks,
Louis

[Jeffrey S. Oishi]

Hi Louie,

You are using an old version of Dedalus. If you upgrade to the latest version, you'll have that function.

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 post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/101f821f-3358-4c7c-b5ef-71ed537dd8ce%40googlegroups.com.

[louisalexan...@gmail.com]

That's what I thought. Daniel already told me that I should upgrade weeks or months ago...

Can I really just do the:
hg pull
hg update
pip3 install -r requirements.txt
python3 setup.py build_ext --inplace

commands and it ll be updated??

[geoff]

That’s all you need to do.

You might need to use:

pip3 install -r requirements.txt ---upgrade

-Geoff

> --
> 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 post to this group, send email to dedalu...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/d9071585-98fb-4bee-841c-864d12e7c772%40googlegroups.com.

[Jeffrey S. Oishi]

Yes, but you might also need to set $FFTW_PATH and $MPI_PATH to the root of wherever those libraries are installed. You probably don't need to the latter, but you probably do need to do the former. Ask if you get any grief from the process!

J

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/5C824E1E-F171-4235-BC07-CDCB41245615%40gmail.com.

[geoff]

Keaton, Louis, 

I’ve been meaning to get back about this.  

Won’t there be issues with boundary conditions? The equations will be singular at the endpoints. This means there is non need for boundary conditions. It works well in many situations. But dedalus requires you to enter as many bcs as derivatives. 

Sometimes it works with specifying a redundant derivative, sometimes things malfunction. 

To be more specific, 

dx(u) - v = 0

(1+x)*dx(v) - u = F

This system only needs a boundary condition at x = 1., x = -1 is a singular point. But dedalus will complain if you don’t put one in. 

I think the results of Louis’ problem will depend on the details. An alternative would be to put the problem in a large finite domain without singular points at the end points, and put true boundary conditions on the problem. 

Does this make sense?

-Geoff 

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBh4TQfgjmTtBx8-YNZBZdtT-bKfa70Ds9_spqJJfbrfXA%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

<quantum_oscillator.ipynb>

[Keaton Burns]

Boyd has a section about imposing boundary conditions under this remapping (17.8).  The takeaway seems to be that while leaving off boundary conditions at the endpoints often works, adding them in doesn’t seem to introduce any problems, and in some cases it makes the solution a lot better.  Not a strong justification though, and maybe more complicated if your solutions don’t decay, but it seems to work here.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/6381DE0A-FC29-41E4-8D38-8497E01ED57D%40gmail.com.

[geoff]

I just looked at the notebook, and sure enough everything looks great. 

I think the reasoning has to be because of the tau method. 

Just truncating the equations and leaving out the BCs gives a system that is singular at the endpoints in the limit of large resolution. 

Dropping the last rows of the matrix and putting in boundary conditions changes the singular structure of the system before replacing it with BCs. My guess is that dropping the rows from the matrix makes the problem non-singular and hence need BCs.

It’s probably worth tracking this down somewhat. 

The Lane-Emden example problem has always puzzled me a bit for the same reason. 

But we know from experience that we can’t play these same tricks for 2D cylindrical time stepping problems with coordinate singularities. 

Very interesting that the QHO works nicely! 

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBhNTqUwjmw_J9RGDN9foZHUMU3SkdRiyhQO5%2BBEwizcyw%40mail.gmail.com.

[louisalexan...@gmail.com]

Hey,

I have a question about the EVP solver.
What's the pencil parameter?

If the EVP is set up with N modes, I'd think that the EVP should be able to spit out N eigenstates.
Do you set a pencil m < N to speed up the calculations? If yes, how do you know that the m eigenstates you're getting are the ones that matters? Is it based on the largest eigenvalues?

Thanks
Louis

[Jeffrey S. Oishi]

Hi Louis,

The pencil is the mode in the perpendicular direction for fourier bases in other directions. If you have a 1D problem, there's only the 0 pencil. If an EVP has N modes, it will spit out N eigenmodes if you use the solve_dense() method.

If, instead, you use the solve_sparse() method, you can specify the number of modes (and the target eigenvalues) for those modes. This can be < N, and it will be dramatically faster.

Here's an example of solve_sparse() in action:

https://github.com/jsoishi/mri_prefers_3d/blob/master/python/mri.py#L130

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 post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/e46e4ab5-988e-4732-aac3-6940ed84766a%40googlegroups.com.

[louisalexan...@gmail.com]

I see.
I'll think about how to play with the target, which is I guess problem-dependent.
Will the sparse solver spit out the mth eigenstates whose eigenvalues are closest in absolute value to the target?

thanks Jeff!

[Jeffrey S. Oishi]

Yup to both!

--
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 post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/9d44e1bb-82ed-42bd-a967-a6b29ab040c5%40googlegroups.com.

[louisalexan...@gmail.com]

ok.

just out of curiosity, in the EVP, is it always better to multiply the eqns such that the denominators are clean? i.e. is there a problem if one of my equations is:

"problem.add_equation(" u + v/r*C + k*w = 0")"

even when C/r is well defined in 0? (r is the independent variable going from -1 to 1)

should I multiply everything by r so that the denominators are ?

[Jeffrey S. Oishi]

Yes, that's true not just for EVPs, but for all problems. You want to minimize the number of modes you need to resolve non-constant-coefficients because that directly impacts the bandwidth of the matrices you need to solve.

j

--
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 post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/8c5623fb-a5f0-4ea0-baf4-3810b1c2e4d2%40googlegroups.com.

[louisalexan...@gmail.com]

I see.
Thanks Jeff!

[Ben Brown]

Louis,

      To follow up on Jeff's "yes, clean up the denominators": even when using solve_dense(), there seems to be value in cleaning up the denominator and reducing the bandwidth of the matrix.  Upasana Das (copied) found some behaviour with MRI eigenvalues where when she cleaned up the denominator she got good, expected results, but when she didn't she found subtle errors in the eigenvalues.  I don't know if she ever quantified or documented that behavior in their publication, but maybe she'll chime in if they did.

Takeaway: best practices from other problems (e.g., IVP) seem to apply to EVP as well.

--Ben

On Thu, Mar 29, 2018 at 11:35 AM, <louisalexan...@gmail.com> wrote:

I see.
Thanks 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-users+unsubscribe@googlegroups.com.

To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/73537e10-2e9d-4d2e-a7db-fa7cdacb1d86%40googlegroups.com.

[louisalexan...@gmail.com]

Hi Ben,

ok.
I'm trying to reproduce the results in the paper of Mayer and Powell 1992 JFM ("Viscous and inviscid instabilities of a trailing vortex").

I'm getting some errors with the unclean formulation.
I'll let you know if I get it to work with cleaned denominators.

Thanks,
Louis

[Upasana Das]

Hi Ben,

Thanks for tagging me. I do vaguely remember telling you about a difference in result long back, however, I did not quantify it back then so cannot clearly remember what was going wrong. It's still possible that it was some other error on my part.

Since then I have solved various EVPs and have noticed that in my cases (ideal MHD equations with or without rotation) cleaning up the denominator does not matter too much (at least in the growth rates (eigenvalues) of the modes, the quantity I am most interested in). I should point out that in my case the independent variable 'r' never goes to zero and also has fractional negative powers. Nevertheless I do try to multiply the equations by a suitable power-law of 'r' to minimize possible errors (and maybe even reduce run-time) --- although having various combinations of fractional powers in the equations means that this is not super clean.


Regards,

Upasana

[louisalexan...@gmail.com]

Hey,

I'm finally getting back to this.

I'm kind of rusty on the tau method; is it related to how to deal with BCs and what you guys are doing?

From what I gather, for chebyshev EVP, DEDALUS removes the first and last rows of the matrix by default and sub in BCs in the truncated system, correct?
Interestingly, I realized that the endpoints are not even part of the grid, so maybe singularities at the endpoints just don't matter.

In Keaton's QHO there aren't physical BCs but you still enforce left and right null BC for \psi. I guess it doesn't hurt since the solution is supposed to be localized close to 0. Could you put in pretty much anything? (left(psi)=left(psi_s)=0 yields similar results except for one eigenmode that has a mirror symmetry /y=0), again maybe that's just due to the fact that you're interested in eigenmodes well behaved at infinity.

One last question (which I may have already asked...): can you ask DEDALUS to use only even or odd chebyshev modes on (0,1)? There are some parity properties of my variables which I could exploit/should probably enforce to have meaningful eigenstates only.

Thanks
Louis