Interpolation/evaluating expressions

[Daniel Lecoanet]

Hey,

I have some data defined on one basis, and I want to interpolate it onto point in a different basis.  Luckily, we have an Interpolate function which will spectrally interpolate for me.  However, as far as I know, actually using it is complicated.  This is because the evaluator is all mixed in with the problem and timestepping.  But it would be great if I could just do something like

evaluator.vars['u_old_grid'] = u_old_grid

for i in range(len(z)):

  evaluator.vars[zi] = z[i]

  u_new_grid['g'][:,:,i] = evaluator.evaluate("Interpolate(u_old_grid,'z',zi)")

Actually, I'm a bit confused why the evaluator belongs to the problem in the first place -- shouldn't it belong to the domain?  This would make it much easier to use the evaluator when doing data analysis.  You would only need to redefine your basis, opposed to right now when you also need to define a problem.

Daniel

[Keaton Burns]

Hi Daniel,

The evaluator is specifically meant to coordinate the simultaneous evaluation of operator trees when time stepping a problem.  For this sort of case you don’t really need to use it at all, instead I’d just directly build and evaluate Interpolate operator objects:

from dedalus2.data.operators import Interpolate

# Some new grid

new_z = np.random.rand(10)

for new_z_i in new_z:

# Build interpolation operator and evaluate it

op = Interpolate(u, ‘z’, new_z_i).evaluate()

# Get the grid values of the interpolant

op[‘g’]

-Keaton

--
You received this message because you are subscribed to the Google Groups "Dedalus Development" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-dev...@googlegroups.com.
To post to this group, send email to dedal...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DjtcPYoe%2Ba%2BzbrU9UsXJxzAKQWSB79UELW7ocURbRA9gw%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.

[Daniel Lecoanet]

Ah ha!  This is exactly what I was hoping existed.  But how does Interpolate know what 'z' means?  Does it ask u which domain it's defined over, and see if there's a 'z' in that domain?

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/1782B5E0-7C9A-4EE2-A464-1C27BF4D5DEB%40gmail.com.

[Daniel Lecoanet]

Wait, actually, I don't think the domain knows what 'z' means -- only the problem does.  But I guess I can just pass the z_basis instead.

Daniel

[Keaton Burns]

Yeah so sort of both.  In the current code, the axis names get associated with the bases in a domain when you build a problem, so you can use the axis name as a string like I did if you’ve already built a problem with the field’s domain.  Otherwise, you can just directly pass basis objects, as you’ve seen.

A quick note that in the new stuff I’m working on (dedalus2-version2?), you always provide an axis name to each basis on instantiation, so you never have this sort of issue.

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3Dg8ZYiT3yUeaCXTFNHi26DaAAdF_WX-UiGRx-JP-tAoZQ%40mail.gmail.com.

[Daniel Lecoanet]

Hey Keaton,

It seems that I still don't understand how the interpolation works.  Here's my loop where I'm trying to interpolate onto a new grid:

for i in range(shape[2]):

  zi = z[0,0,i]

  u['g'][:,:,i] = Interpolate(u_old,z_basis_full,zi).evaluate()['g'][axslice(2,0,1)].reshape(shape[0],shape[1])

  v['g'][:,:,i] = Interpolate(v_old,z_basis_full,zi).evaluate()['g'][axslice(2,0,1)].reshape(shape[0],shape[1])

where shape = z.shape

This not correctly reproducing the u_old and v_old field.  For example, u is identically zero for x less than about Lx/3, even though this is not the case for u_old.

Any thoughts?

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/F3AE0BA6-EC39-467B-A804-2829E891A512%40gmail.com.

[Keaton Burns]

Hi Daniel,

I think there’s actually a binding problem that occurs when you evaluate an operator and access the data in the same expression.  Basically the field objects get recycled and their data may get cleared when there are no python references to them, which is actually the case when the field gets returned by SomeOperator.evaluate() but isn’t bound to a reference.

Try just separating those expressions into separate lines, something like:

interp = Interpolate(u_old,z_basis_full,zi).evaluate()

v[‘g’][:, :, i] = interp[‘g’][:, :, 0]

I’ve also attached a quick script that sets up a simple example of doing this.

-Keaton

[Daniel Lecoanet]

Hmm...  It's still now working for me.  I tried doing the following:

for i in range(shape[2]):

  zi = z[0,0,i]

  u_old_interp = Interpolate(u_old,z_basis_full,zi).evaluate()

  v_old_interp = Interpolate(v_old,z_basis_full,zi).evaluate()

  u['g'][:,:,i] = u_old_interp['g'][:,:,0]

  v['g'][:,:,i] = v_old_interp['g'][:,:,0]

u_old_slice = Interpolate(u_old,z_basis_full,Lz_rad/2 + Lz_conv).evaluate()

v_old_slice = Interpolate(v_old,z_basis_full,Lz_rad/2 + Lz_conv).evaluate()

u_slice = Interpolate(u,z_basis,Lz_rad/2 + Lz_conv).evaluate()

v_slice = Interpolate(v,z_basis,Lz_rad/2 + Lz_conv).evaluate()

logger.info(np.max(np.abs(u_old_slice['g'][:,:,0]-u_slice['g'][:,:,0])))

logger.info(np.max(np.abs(v_old_slice['g'][:,:,0]-v_slice['g'][:,:,0])))

and got differences of order unity.

I mean, this should work in general provided that Lz_conv+Lz_rad/2 is in both z_basis and z_basis_full, and u_old, v_old are defined over z_basis_full, and u, v are defined over z_basis, right?

Daniel

--
You received this message because you are subscribed to the Google Groups "Dedalus Development" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-dev...@googlegroups.com.
To post to this group, send email to dedal...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/8F5B77A9-E55E-42E7-9D92-63B6F9EAE7D9%40gmail.com.

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

Interpolate(u_old,z_basis_full,zi).evaluate()

,

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DgXnh7qmNGCKZh0pbPA0YPqDhRx_YWJsy-PnbeRHGJjmA%40mail.gmail.com.

[Daniel Lecoanet]

Ok, so it works in serial, but not in parallel...

Daniel

[Keaton Burns]

Would it be possible to isolate the issue into a standalone script?

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DjnCeABPzj-b_N%3DsNDigbFpz9fu4gCU0jdUa%2Bf7EQFG5A%40mail.gmail.com.

[Daniel Lecoanet]

If you run your script on 4 processors (mesh=[2,2]), half the processors have small errors, and the other half have max errors of 2.  Is that expected?

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/3BF066DD-5289-487A-A9FA-1542B12D2D43%40gmail.com.

[Daniel Lecoanet]

If you run with mesh=[1,4] or mesh=[4,1] then you get the correct answer.

Daniel

[Keaton Burns]

Yeah that won’t work if you’re not looping over global z arrays, so in this case for a 2D mesh.  If you have to do it in parallel, you should loop over the global z array and then just assigning to the local grid for the local points.  When you said you were doing this for analysis I also just assumed you were doing it interactively/serially.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DifToKCmphW%3DdwygrdwQAihFgQyTnxGJ0sMV9qwybhKfQ%40mail.gmail.com.

[Keaton Burns]

Those just get reduced to a 1D mesh=[4]

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3Dge-aLUy0g0jJv1KENZR-_0Vi5k1R7-twnuR8pvMFLF4w%40mail.gmail.com.

[Daniel Lecoanet]

Ok, I see now why you need to go over the full z array for doing this in parallel...  I'll test that out and make sure it works.

On a separate note, did you have an easy example of using GeneralFunction you could point me to?

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/E2700F45-B755-48C7-84AE-CF019B54742E%40gmail.com.

[Keaton Burns]

Here’s the parallelized version of my earlier script.  There are a few subtleties in the conditional (by-process) assignment, since you need to make sure no transposes are attempted by just a subset of processes:

[Daniel Lecoanet]

It's still not working in my code, even when I use the full z basis:

zgrid = z_basis.grid

u_fullz = np.zeros([shape[0],shape[1],zgrid.size])

v_fullz = np.zeros([shape[0],shape[1],zgrid.size])

# Set the ICs

for i in range(zgrid.size):

  zi = zgrid[i]

  u_old_interp = Interpolate(u_old,z_basis_full,zi).evaluate()

  v_old_interp = Interpolate(v_old,z_basis_full,zi).evaluate()

  u_fullz[:,:,i] = u_old_interp['g'][:,:,0]

  v_fullz[:,:,i] = v_old_interp['g'][:,:,0]

u['g'] = u_fullz[:,:,slices[2]]

v['g'] = v_fullz[:,:,slices[2]]

logger.info(Lz_rad/2 + Lz_conv)

u_old_slice = Interpolate(u_old,z_basis_full,Lz_rad/2 + Lz_conv).evaluate()

v_old_slice = Interpolate(v_old,z_basis_full,Lz_rad/2 + Lz_conv).evaluate()

u_slice = Interpolate(u,z_basis,Lz_rad/2 + Lz_conv).evaluate()

v_slice = Interpolate(v,z_basis,Lz_rad/2 + Lz_conv).evaluate()

logger.info(np.max(np.abs(u_old_slice['g'][:,:,0]-u_slice['g'][:,:,0])))

logger.info(np.max(np.abs(v_old_slice['g'][:,:,0]-v_slice['g'][:,:,0])))

I'll see if I can break your test problem in the same way.

Daniel

[Daniel Lecoanet]

Hmm... when I do this with your test problem, it seems to work fine...

Daniel

[Daniel Lecoanet]

Your method doesn't work on my problem either.  I'll try to find a minimal version of my problem which breaks.

Daniel

On Wed, Nov 5, 2014 at 7:49 PM, Keaton Burns <keaton...@gmail.com> wrote:

Here’s the parallelized version of my earlier script.  There are a few subtleties in the conditional (by-process) assignment, since you need to make sure no transposes are attempted by just a subset of processes:

--
You received this message because you are subscribed to the Google Groups "Dedalus Development" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-dev...@googlegroups.com.
To post to this group, send email to dedal...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/A2A9AADE-7D4A-44C2-8BAF-D3722C97ABF8%40gmail.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3Dg%2BnGusCnpZoBBoS%2BB5VZSu%3DR5aw_w_5EoJoRvbLHS4kA%40mail.gmail.com.

[Daniel Lecoanet]

Ok, if you replace

zbc = dev.basis.Chebyshev(nz, interval=[0,1])

with

zbc1 = dev.basis.Chebyshev(nz, interval=[0,0.5])

zbc2 = dev.basis.Chebyshev(nz, interval=[0.5,1])

zbc = dev.basis.Compound((zbc1,zbc2))

your script no longer works on 4 processors (mesh=[2,2]).

Daniel

[Keaton Burns]

Ah yes this is a bug I fixed with the new parser.  The functionals don’t work completely right on compound bases — they only set the first subbasis's data to the computed constant.  I’ll try to make a quick patch for this.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3Djx0V%2BVn_TO80FUmgpDV0WJ4N9j3yY6iRhD9Yr2NJM1gw%40mail.gmail.com.

[Daniel Lecoanet]

Ah ha!  Yeah, it would be great if you could fix that up...

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/754A090D-E8DF-4E89-8598-594DA74F5875%40gmail.com.

[Keaton Burns]

Just pushed a patch to the main repo.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3Dhpjyjo6_ZzmJtfOodyR7nZX2Ufn2_PGcpcD%2BZ35272vA%40mail.gmail.com.

[Daniel Lecoanet]

Still a problem... now with dealiasing.  If you run your test problem on 16 processors (mesh=[4,4]), set nz=32*3/2, and dealias all the z-bases by 2/3, the interpolation doesn't work properly.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/895C5700-A859-4000-9B9F-82FF890AD253%40gmail.com.

[Keaton Burns]

Patched the patch!

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/CAJoYf%3DjHB9e3wEsq8CpuwCvR2_ZKxD-s4hW-vGtHYcz7EvfLdQ%40mail.gmail.com.

[geoff]

not to say that we all aren’t some of those things. I’m just definitely not “move fast”.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/2EC2BB62-1E7F-49D1-A217-F447FD45CDFA%40gmail.com.

[Daniel Lecoanet]

Looks like it's working!

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-dev/A815759D-4864-4BB5-9E9F-B044184EFB97%40gmail.com.