Saving Eigenvectors

[afra...@wisc.edu]

Hi,

I am familiarizing myself with Dedalus' EVP tools. So far I've gotten the hang of plotting eigenvalues and eigenvectors. However, I would like to save the eigenvectors of an EVP to an h5 file so that I can take a look at them later on without having to re-run the solver. I see that solver.set_state(i) can be used to set the state to the ith eigenvector, and I see that for IVPs you can use something like analysis.add_system(solver.state, layout='g') to save checkpoints. I'm having trouble combining these two in order to save the eigenvectors of an EVP. Is there an easy solution that I'm missing?

Thanks in advance,
Adrian

[Christopher Marcotte]

Hi Adrian,

Wondering if you've figured out a solution to this on your own> If not, I hope this post motivates someone in-the-know to respond, as this would be remarkably useful for me, as well.

[Susan Clark]

Hi Adrian and Christopher,

You can write whatever h5 files you like using something like the following (in which psi has been set to the preferred eigenvector using set_state):

with h5py.File(filename, 'w') as f:

    psi = f.create_dataset("psi", data=psi['g'])

Typically I save my eigenvectors this way, along with any relevant parameters as attributes, i.e.

f.attrs["Pm"] = Pm

This isn't really Dedalus-specific, just an easy way to save your eigenvectors in h5 files. See also the h5py documentation.

Cheers,

Susan

On Mon, Jul 10, 2017 at 11:26 AM, Christopher Marcotte <christophe...@gmail.com> wrote:

Hi Adrian,

Wondering if you've figured out a solution to this on your own> If not, I hope this post motivates someone in-the-know to respond, as this would be remarkably useful for me, as well.

--
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/2d1a8ebe-d7fd-4541-9c1d-b244de88ecd2%40googlegroups.com.

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

[Christopher Marcotte]

Thanks Susan!

I was hoping for something which preserves the `field` type so that I could use the built-in methods to compute inner products (say, between the adjoint and canonical sets of eigenfunctions) accurately. I recognize I could re-compute the weights and then the inner products from the domain grid, but this seems like a duplication of effort and thus prone to errors. I suspect, since you manage with the above method, that no such functionality exists as-is?

[afra...@wisc.edu]

Hi Christopher,

Glad to hear someone else is interested in computing inner products! I have gotten an ugly form of what you're looking for to work, but I haven't yet gotten around to computing inner products. I'm interested in saving eigenvectors, then saving snapshots from nonlinear IVP simulations, then computing inner products between eigenvectors and those snapshots. Please let me know if you wind up working on something in that general direction, as I'm curious to hear how you compute your inner products.

Anyway, here's how I save eigenvectors. It's not pretty, but it works:

        ev = solver.eigenvalues
        vecs = solver.eigenvectors

        analysisc = solver.evaluator.add_file_handler('Eigenmodes_c', iter=0)
        analysisc.add_system(solver.state,layout='c')

        analysisg = solver.evaluator.add_file_handler('Eigenmodes_g', iter=0)
        analysisg.add_system(solver.state,layout='g')

        for i in range(len(ev)):
            if np.isfinite(ev[i]):
                solver.set_state(i)
                analysisc.create_current_file()
                analysisg.create_current_file()
                solver.evaluator.evaluate_handlers([analysisc,analysisg],world_time=kx,wall_time=np.imag(ev[i]),sim_time=np.real(ev[i]),timestep=i,iteration=i)
                #analysis.process(i,i,i,i,i)
            analysisc.set_num += 1
            analysisg.set_num += 1

This saves a separate file/folder for each eigenvector, once for grid space and once for coordinate space. I'd rather set it up to have one file where each iteration is a separate eigenvector, but to be honest I stopped working on making this work nicer as soon as I got it working at all. Please let me know if you have any questions.

Adrian

[afra...@wisc.edu]

By the way, the arguments world_time, wall_time, sim_time, timestep, and iteration within solver.evaluator.evaluate_handlers() need to be specified as something (they were originally written to work with IVPs, not EVPs). I chose to set them to some of the parameters/eigenvalues, but you can set them to anything you want. But you'll get an error if you don't set them to something.

Adrian

[Ben Brown]

Christopher,

      You can read the saved eigenvectors into a field object and the obtain all of the integration, derivative, interpolation methods etc. that fuels objects normally have. The pseudo code is more or less like so:

domain = dedalus.Domain(stuff)

EV = domain.new_field()

data = h5py.File(stuff to read from psi file)

EV['g'] = data

I'm on my phone, so can't pull a functioning example together, but if you get stuck that could be generated this afternoon. 

Or I may be completely off base in what you're looking for from fields; if so, let us know 

--Ben

On Mon, Jul 10, 2017 at 4:37 AM Christopher Marcotte <christophe...@gmail.com> wrote:

Thanks Susan!

I was hoping for something which preserves the `field` type so that I could use the built-in methods to compute inner products (say, between the adjoint and canonical sets of eigenfunctions) accurately. I recognize I could re-compute the weights and then the inner products from the domain grid, but this seems like a duplication of effort and thus prone to errors. I suspect, since you manage with the above method, that no such functionality exists as-is?

--

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/03bd9ecb-f362-4dbb-bf11-c9d6c741b749%40googlegroups.com.

[Jeffrey S. Oishi]

Hi Christopher,

You're right, as is, there is no way to serialize fields to disk with all of their attributes. It is actually a bit of a tricky problem, though it is in principle possible. Also, as Adrian noted, the output mechanism is definitely geared toward IVPs. Now that we have a significant community of EVP users, we should definitely discuss some additions to make saving eigenvectors better.

I've implemented a somewhat hacky set of functions to do what Ben suggests quasi-automatically, which I've attached here. I am working on polishing it off into something general enough to go in the main Dedalus codebase, but it's not quite there yet, and I'm a bit busy at the moment, so I'm not sure if I'll have time to work on it much. You can use the attached code like so:

import fields_from_file as ff

datafile= "snapshots_s1.h5"

zeta = ff.fields_from_file(datafile,['Fourier','SinCos'],'theta',meta={'y':{'scale':1.0,'parity':-1,'constant':False}})

The meta keyword argument is needed only for SinCos bases. There is currently a limitation that I didn't get around to implementing the calculation of the implementation of the limits (needed for specifying a basis) for Chebyshev bases. This is necessary because HDF5 files from Dedalus store only the grid points, which do not include the end points. Calculating the end points from the grid points for Chebyshev polynomials is trivial. 

I'd just like to point out that we welcome code contributions from anyone. If you have something you think is useful, you can fork the dedalus repository on bitbucket, add your changes, and then issue a pull request. We'll peer review the code and once it passes, add it in. 

J

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAHqBLzzoB%2B5kOo38_16RmMo_mhrX9JwNE-sxMrhUeyoB%3DcGjQw%40mail.gmail.com.

[Christopher Marcotte]

Wow, a lot of discussion.

My interest in saving eigenvectors primarily stems from the siloing of the results from the EVP solver; I simply don't know how else one could reasonably re-use those results after calling the EVP solver from within a function. I am certain some part of this stems from my own python ignorance. 

For my purposes, the most convenient method for saving the eigenfunctions is to use numpy's built-in `savez` and `load` methods. I have a function which takes as input a domain and fields, sets up and solves the eigenproblem using linearization, and then iterates through the N indices of the eigenvalues, calling set_state() appropriately, and copying the coefficients to a pre-allocated numpy matrix, which is returned along with a copy of the similarly ordered eigenvalues. 

I call this function twice -- once for the right eigenfunctions, and once for the left eigenfunctions. The latter use the explicitly defined Jacobian, rather than the parsers linearization.

For computing the inner products, it's most convenient to build a new np.complex128 domain, define a set of fields for the left eigenfunctions of the system (W#) and the right eigenfunctions of the system (V#), and then again iterate through the N modes, copying the amplitudes to the newly defined fields, and then calling (for a two-variable system)

WV[m,n] = de.operators.integrate(W1*V1 + W2*V2, 'x').evaluate()['g'][0]

for each of the N^2 pairings. It's amazingly slow, but perfectly functional, and it avoids some frustrations I ran into with re-using domains.

There are some lingering issues I'm struggling to solve regarding the EVP output, but that's a discussion for another time.