Vlasov-Maxwell solver?

[Zach Williams]

Hello Dedalus Community!

I'm toying with the possibility of using Dedalus (d3 specifically) to solve collisionless kinetic problems in plasma physics, i.e. solving the Vlasov-Maxwell equations. As I understand it, d3 does have the ability to handle integro-differential equations. However, I could see the velocity space dependencies in particular being problematic, as I think Dedalus is only equipped to handle one non-periodic dimension in a given problem (though that may have been a d2 thing, I'm still pretty green to d3 and working on learning). Doing a search for "Vlasov" in the forums reveals only a couple posts from many years ago, so I thought I would check and see if solving these sort of systems is something anyone has thought about recently, or if it's even within the scope of Dedalus to look at these sort of problems. I figure broadly speaking that the answer along the lines of either (1) yes people have thought about this and worked on it, it's can be done or (2) Sounds plausible, but no one's really tried so good luck or (3) Dedalus is not well-suited for these sort of problems.

I know there are plenty of codes out there that solve these sort of problems and I could look elsewhere, but frankly speaking I enjoy working with Dedalus, particular as a tool for onboarding students to computational research. I appreciate any and all thoughts!

Thanks for your time.

Cheers,

Zach

[Zach Williams]

Hi all,

We figured we'd give it a try in the name of exploration, at least in the presumably manageable case of 1 spatial dimension (treated periodically) and 1 velocity dimension (Chebyshev). Just to get our feet wet, we're attempting to model as simple a kinetic case as possible, the two-stream instability, consisting of an initial condition with two counter-streaming Maxwellians (one perturbed to get things going). The behavior of the distribution function for this instability is very well-known, so we thought it would be a nice test case. This is our first crack at doing anything of this sort in Dedalus so I'm experimenting with different input parameters, at this point we're able to get the code to run but it crashes pretty quickly with NaN outputs. I'm sure there's many things we're doing wrong with the implementation, so if anyone could take a look at the attached code we would appreciate any tips or input on how we might do this better, even for this basic scenario. The model is taken from any standard introductory plasma textbook, such as Chen's Introduction to Plasma Physics, see equations 7.1 and 7.3 here: https://farside.ph.utexas.edu/teaching/plasma/Plasmahtml/node90.html for example. We normalized those equations to debye length, plasma frequency, and equilibrium density, if it matters.

Thanks for your time!

Cheers,

Zach

[Daniel Lecoanet]

Hi Zach,

You might want to introduce some diffusion into the equations. The spectral method introduces very little numerical diffusion into the solutions, so often we will introduce some diffusion directly into the equations in order to prevent features from developing at the grid scale. With some diffusion, I think Dedalus would be able to solve 1+1D Vlasov-Poisson equation as you’ve described.

Daniel

-- 
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/735471a5-c936-4051-935a-3ba25f8b325fn%40googlegroups.com.
<kinetic_IVP.py>

[Zach Williams]

Thanks for the quick response, Daniel! That's a good thought. We just tried implementing second- and fourth-order dissipation in the dt(f) equation, as well as fourth-order dissipation in the phi equation (and by this I mean Lap(Lap(quantity))), unfortunately it didn't really have any difference on the crash. That is likely an important feature to keep, but it seems there's more wrong with our implementation yet.

One place I suspect issues may arise is in the velocity boundary conditions. Typically, as I understand it, Dedalus requires boundary conditions for Chebyshev domains, yet for the script included above it is able to run without boundary conditions. When we look at 2D slices of the results, there are some seemingly anomalous structures forming at the edges of the v domain, so it's likely connected to that. We would be happy with the boundary condition of f = 0 at the edges of velocity space, but we can't add such conditions here without raising non-square errors.

Thanks again, any other thoughts or ideas are welcome!

[Piyush Grover]

Zach,

I have solved a Vlasov-like PDE that shows up in flocking problems, specifically Eq. 7 in this paper:https://arxiv.org/abs/1611.02194

I was able to numerically verify the stability results in that paper using this code.

I used similar domain as you, and also played with artificial diffusion (in x) as suggested by Daniel.

P.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/47C1D071-434B-475D-B759-D56D85A9FF15%40northwestern.edu.

[Zach Williams]

Dear Piyush,

Thanks for the reply! That's encouraging to hear that others have had success studying similar systems.

Would you be willing and able to share any of the code used to successfully solved the PDE so I can compare your method of artificial diffusion with ours?

Thanks again,

Zach

[Enrique Rojas Villalba]

Hi Zach, 

I used Dedalus2 to reproduce the electric field energy decay rate of Landau-damped Langmuir waves some time ago. I haven't checked if it's correct, but it looks good enough. I am attaching the jupyter notebook.

Kike 

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/903327f7-6d4f-4213-ac91-77cd7e3e7cf9n%40googlegroups.com.

[Keaton Burns]

Hi Zach,

Along with the above suggestions, you might want to put in an artificial diffusion in v, as well as x, both to help regulate f and to accommodate the f boundary conditions you mentioned at the edge of velocity space. To do this, you would need to add two tau terms to keep the system square. You can also try adding the tau terms and boundary conditions without any diffusion in v, but I think that’s also likely to have stability problems.

Best,

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CACcJjsFMpPJWY%3D4MtkAMkyutkAYKSEhJvQ1vSZF-Sp6%3DuwP%2B2g%40mail.gmail.com.

[Zach Williams]

Hi all,

Thanks for the additional input! As an update, we corrected a couple of mistakes on our end and added numerical dissipation in the form of a laplacian (in x and v) to our distribution function evolution equation. Unfortunately, it still seems to eventually crash, even with pretty large (0.01) amplitude dissipation and velocity space boundary conditions set to zero. My concern is that dissipation this large has a pretty substantial effect on the linear properties, which muddies the purpose we have in mind for our work, so even if dissipation this large fixed it that wouldn't be ideal.

We've tweaked the timestep some, that seems to allow it to run a bit longer but eventually the same problem occurs with quantities (f) exploding to NaNs. I'm none too familiar with a good CFL condition for a kinetic description, so I figured as a baseline that (Lx/Nx)/Lv would give a reasonable timestep lower bound, this might be the problem and we'll explore it more but a couple orders of magnitude change in dt hasn't fixed it yet.

I've attached the updated code in case anyone might take a look and see if there's other suggestions to be made. Note that we're now evaluating only fluctuating quantities and have implicit/explicit separation between linear/nonlinear on the LHS/RHS. The assumed equilibrium consists of two shifted maxwellians, and the initial condition provides a harmonic perturbation to one of them.

Thanks again for all the input so far, it is greatly appreciated!

Cheers,

Zach Williams