Setting an initial condition in spectral space and problems with cfl.compute_dt

[James Penn]

Hi all,

I've just written my first script using Dedalus and am having a couple of problems I can't find documentation to, can you give some advice?

I'm trying to solve the barotropic vorticity equation on a doubly-periodic beta plane.  Taking inspiration from the tutorial notebooks and some of the advice on this list to avoid a singular matrix, I wrote the script attached.

My questions:

1. What is the correct way to set the initial condition in the spectral domain?

   I want the IC to be a random excitation of specific wavenumbers.  I thought this would work:

    zeta = solver.state['zeta']

    k = domain.bases[0].wavenumbers[:, np.newaxis]

    l = domain.bases[1].wavenumbers[np.newaxis, :]

    ksq = k**2 + l**2

    init = np.random.random(ksq.shape) + 1j*np.random.random(ksq.shape)

    init[ksq > 50**2] = 0

    zeta['c'] = init

   but the resulting grid has a single spot in the corners (see attached image of initial cond).

2. cfl.compute_dt() is returning infinity after the first timestep.

   Velocities u and v are diagnostic in my computation, why is this not working?

3. The evolution is slower than an equivalent single-threaded script I wrote that uses numpy fft routines.  Am I doing something silly in defining my problem that is causing it to run slowly?

4. Is there a standard method to add hyperdiffusivity or a high wavenumber filter?

Many thanks in advance,

James

[Ben Brown]

James,

    2. You need to set an initial time step size when you create the CFL object. And a max timstep is often useful as well. See the examples/ivp/2d_rayleigh_benard for a fuller CFL example (I think; not at computer at the moment and can't double-check). 

Right now your velocities don't exist when the first timstep is calculated.  They will for subsequent time steps however.  This is tied into them being algebraic in time (diagnostic) rather than differential.

3. Could you quantify the speed difference you're seeing between your custom single-threaded version and this implementation?  Was your single-threaded version run I the same Python/numpy install as dedalus?  We need a ballpark idea of the difference before we can suggest solutions. 

Thanks!

--Ben

--
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/51553bc9-ca33-449e-98ca-5c110a37f2f9%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[geoff]

James,

I’ve made some changes to your script that highlight some of the features of dedalus. I did this in the hope that other people would also be able to use them.

1. What is the correct way to set the initial condition in the spectral domain?

   I want the IC to be a random excitation of specific wavenumbers.  I thought this would work:

Initial condition are always complicated. In general, you want to satisfy linear constrains of the initial data that are not easy to generate. Suppose you want to set the vorticity or velocity, and want a consistent stream function. Then you have to solve for it. You can use dedalus for this also, but you have to set up the initial problem. 

As you might guess, it can get a lot more complicated depending on the system. Hopefully this example gives you enough information to do this in other cases.

2. cfl.compute_dt() is returning infinity after the first tilmestep.

Like Ben said, you can set a max time step size. But the reason you were having problems was more fundamental than that. 

The initial velocities were zero. You set the “diagnostic” conditions as problem equations. This means you won’t get them without solving. Setting the velocities as substitutions means they will get computed from the stream function whenever you ask for it. They will be non-zero as long as the stream function is non-zero. Which you need to set as an initial condition. 

  Velocities u and v are diagnostic in my computation, why is this not working?

3. The evolution is slower than an equivalent single-threaded script I wrote that uses numpy fft routines.  Am I doing something silly in defining my problem that is causing it to run slowly?

This may have been because you were defining a lot of unnecessary “diagnostic” equations. You don’t need to do this for your problem. Now it only computes what’s needed. Another reason could have been the plotting you were doing each time step. Try it now and see how it goes. We would be very interested uncontrolled speed comparisons. 

4. Is there a standard method to add hyperdiffusivity or a high wavenumber filter?

The way to add hyper diffusion is to define what you want it to mean and include it in your equations. As a philosophy, we want to avoid guessing “standard” parameterisations and including them wiht the code. Reasonable people can disagree what is a good hyper diffusion. We want them all to be able to make it their favourite way.

Spectral filtering is a bit more subtle. We can talk about that more if you find it’s really what you want. In the current version of dedalus it’s a little tricky to make a really good spectral filter. It can be done, but it’s a hassle. We know about this, and are making small changes that will make it much easier. These changes will also allow fractional Laplacians, etc. I don’t want to get into the details of this unless it’s really needed. 

Hope this helps. Let me know if you follow the code in the example. 

[James Penn]

Hi Ben, thanks for your reply.

2. I've changed the setup of the cfl to 

cfl = flow_tools.CFL(solver, initial_dt, max_change=2.0)

and it seems happier now, thanks.

3. For a resolution of N=256 and 1000 iterations my original implementation runs in approx 13 seconds.  The dedalus code as previously attached, but with no plotting performed, runs in 123 seconds (not including the time to create pencil).  I've attached my numpy barotropic vorticity code so others can compare.  Let me know how it runs comparatively on your machine, it may be an issue with how I have installed dedalus (using the install script on the website did not work for me without tweaking, we can follow up on that in another thread if anyone devs are interested on where it failed on my Mac).

Cheers

James

[James Penn]

Geoff,

Many, many thanks for taking the time to do that, it looks worlds better.  I didn't realise you could define substitutions, that makes writing the diagnostic equations much more obvious.  I'll take a full read of the code tomorrow and get back to you with any questions.

I think perhaps in your previous attachment these lines:

psi = solver.state['psi']
psi = init_solver.state['init_psi']

should be something like:

psi = solver.state['psi']
psi['g'] = init_solver.state['init_psi']['g']

otherwise you won't get any evolution?  Then the realtime plotting works for me (I agree that realtime plotting is not the most performant way to run the code!)

James

[geoff]

This took about 19 seconds to run on my laptop. Your other script took about 17 seconds. It’s not clear how the resolutions compare 1-to-1. We are using 3/2 dealiasing, you are using some kind of spectral filter. But it’s clear they are doing about the same. 

See how long this takes for you. That will tell if you have a reasonable install. 

[James Penn]

That runs in about 40s on my computer, 14s for the original numpy script. Far more efficient than my first attempt with dedalus but not sure why I'm still seeing comparatively 3x slower than you are, I'll try a recompile of the libraries.

I adapted the code you sent in your message yesterday, added a file handler and then ran using mpi which all worked as expected and scaled up performance.  It was pleasantly straightforward to do - I'm now working with IT to get dedalus installed on my workstation so that I can run on many more cores.

Really impressive project, now that I've spent a little time in the code I understand why you've arranged things as you have.  I'll post a final version of the barotropic vorticity script when I've implemented the proper initial condition and hyperviscosity.