Post-Processing and Visualization of Astrophysical Data with PyPLUTO

Speaker 1 [00:00]

So, the last talk of the conference before the closure, congrats, congrats for making it here. It is my pleasure to introduce Giancarlo Mattia, a postdoc at the Max Planck Institute for Astronomy at Heidelberg. He will tell us and entertain us with a Pi-Pluto code for analysis and visualization of numerical simulations.

Speaker 2 [00:23]

Hi, I'm Giancarlo. Thank you very much for staying until the very last talk. Yeah, hooray to you. So the title sounds a little bit like exaggerated, I think, so I changed it in Lessons I Learned While Plotting Data and Building Plotting Routines for a Community Larger Than Me, Myself, and I. This is because PyPluto was a Python package built more than 10 years ago alongside the Pluto code, which is a popular code for astrophysical simulations. Then at some point, the maintainer couldn't do it because he was busy with other stuff related to astronomy. So at some point, I started putting some fixes here and there. And then I asked the maintainer, Varga Vaidya, hey, can I take over? It was 2021. And after five years, the code has been basically completely rewritten. and I had to face the fact that I was not programming anymore just for myself, but also for other people that were relying on the tool. So, I am an astronomer, and most of the work that I do, and most of the time that I invest, has to go to astrophysical questions, like how do jets evolve? Just a disclaimer, not these jets, these jets. Much better. Now, there are two approaches that one can follow, which is through theories, or through equations, and through observations. So, ideally, if I had infinite time and knowledge, I would observe jets propagating through their external medium, or I would just solve the equations of fluids embedded into an electromagnetic field. Unfortunately, life is too short, so these jets move on timescales of millions of years, so good luck trying to see one of those moving. And the models are very hard to interpret, So we have to rely on numerical simulations if we want to discuss the evolution, the temporal evolution of these objects. And I use G-Pluto, which is a code on C, C++, no accelerated on GPU, because historically C++ was a synonym of speed. And Pluto solves the time-dependent differential equation in multidimension, trying to encompass multiphysics of fluids embedded in electromagnetic field. And it's built in a way that it can be massively parallelized and can run almost with ideal scaling on thousands of cores and several GPUs. Because if we want to have some really good results, we need high-resolution grid. This brings us to the issue that we need to deal with large data sets. By large, I mean from the megabyte to the terabyte. I know that some people in the data science community that take data from LinkedIn will just laugh at me and say, this is not large at all. But well, for us, it kind of is, especially because we have bachelor, we have master student, and the storage that we have strongly depends on the university, on the institute that we work. So sometimes we may have to run our post-processing on a laptop, when then gigabyte and terabyte becomes a little bit inaccessible. So the issue of translating a simulation to an output from which we can extract some physical interpretation of what's going on is not just a problem of astrophysics. It's a problem of every simulation-based field. It can go to biology, it can go to economics. If you run a simulation, then you have to be able to interpret the data. And it means that you have to load them from some file or from whatever. You have to manipulate it, and you have to show them in a meaningful way. And that's key. Otherwise, it's really hard to convince other people that you're doing a good job. Now, expectations are always very optimistic. So ideally, we think, okay, our simulation is going to run, no bugs, no issues, no flow. I don't have to fix my code. I'm going to just bulldoze it. And then the plotting time and the analysis time will be at the minimum. And then we have to face reality where the simulations run much more than expected because we face some bugs and issues and fixes and whatever. And then we have to load the data, compute some extravagant quantities. We have to visualise it, but then we decide, okay, we are not good with this and that, so we refine the images, then we have to publish a paper, and then we polish it even more and even more. And in the end, the majority of our time goes into the analysis and the plotting. And that's the first lesson that I learned, that if your analysis pipeline is slow, your science is slow. But that's not on you. It doesn't make you a bad researcher. It just means that the tool that you are using to load the data and to use your data is inefficient. It doesn't say that you're a bad researcher. It just means that that tool can be optimised. How do we solve this problem? Well, for Pluto, for that particular framework, I am trying to build a package that reduces as much as possible the analysis and the plotting time. Just bear in mind that the simulation time in the right plot didn't change because this is something they do on some other portions of my time as Pluto developer, but this has nothing to do with Pluto. If the Pluto simulation is low, then we have to fix it on the Pluto level, of course. And also notice that I didn't reduce the bars to zero, because that would be unrealistic. And we should not aim for unrealistic goals. We should make the process as efficient as possible, but we shouldn't claim that the time will be zero, because that's simply not true. How can we build such a tool? So we tried to follow three keywords. The first is simple. Coding the data, analysing the data, and even doing some small animations like this shouldn't take so many lines of code. Because otherwise the user gets upset and they're like, oh, but this seems so simple, why can't I do it immediately? It has to be flexible, because as astronomers, our ultimate goal sometimes is publishing scientific papers. It shouldn't be, but most of the time it is. So we have to produce in a reasonable amount of time products that are good for a paper. And I don't know if you attended the talk about visualization. It was yesterday morning. And the speaker was saying that scientific articles have slightly different requirements sometimes. And this is kind of true. Sometimes you have to polish your plot basically because you are doing it for a paper, while otherwise you wouldn't have done anything of the sort. And also, since Pluto is a tool that is also used in tutorials for bachelor student, master student, and sometimes people don't have an extensive knowledge of Python and Matloclyp, we built a graphic user interface, because even if you don't really want to do scripting, you can grasp some interpretation of your data, even in a very simple way, but you can just press some buttons and see what you get. Now, let me bring you a little bit in the process of loading, manipulating, and then refining the data. So mostly for legacy reasons, PyPluto has two facade classes, and the first is load. That loads the data. The reason is that when users want to load the data, they don't want to write more than one line of code. They want to load the data, and that's it. So that's what we aimed for. I also tried a small comparison, loading the data in the most direct way possible, just through NumPy, and then using the PyPluto routine that loads the data. And I try to see how slower is PyPluto, because we may have something else, some additional checks and stuff. And notice that when the resolution is slow, do you see the mouse? Yeah. We have a very, very, very bad comparison, which, by the way, this is in units of seconds, so I would challenge anyone to see the difference while you're running your data once. But still, a pretty big difference, but when the resolution increases, and nowadays the resolution can be up to 2,000 or 3,000 cubes, so much more than this. These are all simulations that I ran on a laptop, not on a supercomputer. You see that the speedup has, like, the time needed is almost identical. So why is it so slow at low resolution? Because we do a lot of check. We check that the path that the user put is consistent. We check that the format that the user put is consistent. We load some grid information. We load some additional information of the simulation that comes from the script or file that comes from the code. Because the lesson number two that I learned is faster is good, but breaking users is not. If I build something that is super fast, but then at some point it crashes and the users cannot understand what the hell is going on, then the users aren't just going to fill my email mailbox with, Giancarlo, why doesn't this work? And I'm going to waste so much more time, and they are going to waste so much more time. So the goal was to build something that the users don't really see. It's small, it's compact from the user perspective, but it's also bulletproof. And also one other thing that I really care about is because a lot of the analysis sometimes has to be done on a laptop. And the laptop has a limited RAM, especially now, because RAMs skyrocket their price. So what we're doing is that all the data are loaded with memory mapping, so that if you need just a portion of the data, you don't have to load the entire file and you have to pass only your entire file through the RAM. And this was particularly useful for very high-resolution simulation because it allows you to use your laptop offline with the Deutsche Bahn when your train is late and it's stuck in the middle of nowhere, and then you can still do something. And you can forget that you are very late because Deutsche Bahn is whatever. Now, the next step is, okay, I loaded the file. What do I do with that? So Pluto, but every grid-based code, is based on a grid structure. Right now, we support a different set of coordinates, but we are also going for some more general relativistic framework. And knowing the shape and the size of a cell is very important for average integrated quantities, but also for visualization. Now, this is a jet. That's a very weird jet, right? It's not like the picture that I showed before. Well, because it's in r-theta coordinates. It's exactly how Pluto works, because it works with r-theta if we are in spherical coordinates. This is much more in human understandable framework, because we transformed the frame according to the grid. And now you see that the jet is actually this one, as expected. And also, notice that there's this small vortex, whose size seems comparable to the one of these lobes here in the r-theta, And it's actually not, it's very different because this is a grid feature. And also, we shouldn't forget that all these simulations rely on a temporal structure. We do simulation also because we want to see the temporal evolution of this system. Now this simulation here of this jet ran for like a couple of million time steps. So you have to imagine that this is a grid of 1,000 per whatever hundreds that I had. And this grid has different values at every time step for a million times. Now I don't have the hard disk to fit one million files, that would be unfeasible. So we have to rely on discrete outputs, sometimes the temporal spacing may be not uniform. So we have to aim for loading multiple outputs sometimes at the same time, because this can be used for comparison, to see trends, and to see time variability. This is the same JETS of before, plus other three simulations that we did, and it's very important because I don't want to bore you with the detail, but once the jets go over this shock, they start accelerating. Now, you don't see it because the qualitative evolution of the jet doesn't really change, but quantitatively, if you check the jet velocity, they increase their speed because they are propagating in a less dense environment and blah, blah, blah. This is something you don't see for one snapshot. This is something when you track the temporal evolution of your system. And this brings me to the third lesson that I learned. Plotting simulation data requires understanding the simulation. If you don't know what your simulation is, you don't know what you're plotting. You don't know what you're analysing. It's much harder to infer some quantities, infer some physics, if you don't know what you're doing, right? Now, you loaded the data, we computed whatever based on the grid structure, temporal structure, now we want to see it. I said PyPluto has two classes, and the second is image. Again, for legacy reasons, but also because the users don't want to spend too many, don't want to think about too many classes. With the image, you can create access plot lines, have color mapped plots, field contour lines, simple animation. You can do a lot of things. The image is a really big facade of some other classes that are on the background. And this really looks like, the features look really like a wrapper of Matplotlib. So one could ask, Gercaro, why did you build a Matplotlib wrapper? Well, it's kind of a Matelotip wrapper, but at the same time, I think it's a wrapper that is tailored for the simulations that we are running. Matelotip has to be very general and massive because it has a much broader applicability. We don't use even 10% of all the features that Matelotip has. So, having a wrapper makes us focus on the things that sometimes have a slight ambiguity or some routines that can be simplified because we know what we are looking for, while retaining all the advantages of MatLootLib and making the user feel a little bit safe, especially because figures and axes and subplots are always accessible with the image class. So, if you want to do something that is very extravagant, you are still able to do it. Just if you want to feel safe, then you can. And now I want just to show a little bit of anecdotal experience about how to have a quick visualization and then how to refine it up to journal publication level. So this was 2024, and the referee asked us to plot one quantity for different simulations in polar coordinates. You can see it's the same jet of before. And we wanted to plot the resistivity, which in two of our models is constant, as you can see here, and in two other models was correlating with the jet structure. Now, this was the first thing that we did. Very easy. Few lines of code. You already see that in the constant model is constant, and in the variable models is variable. But I looked at this plot many, many times, so to me it's kind of obvious, but I understand if it's not. Here there's an order of magnitude of difference. Here there's an order of magnitude of difference. This colour map could be summarised in one single value. This as well. But they are the same because there's this scaling. So you can do a direct comparison but it's not super efficient. So the second thing that we did was to polish it a little bit with a couple of few more lines of code. Same quantity, same simulations. We merged into a single view so you don't have to look left, right, top, bottom. You can just look more efficiently. We rescaled the axis because if you notice there were like 25,000, 50,000, 75,000. Now it's more in scale. We used a consistent color scaling along all the simulations. And in this way, you can see clear spatial structure, especially in the variable model on the top. And that allows for a more improved comparison. Is it good? I think so. You can doubt it, but I think it's pretty good. It's not enough for publication level. And this was very similar to what we published. So again, same simulation. We added the units, very important. We used a discrete color map because for that plot, we really wanted to see more the order of magnitude rather than the more refined values. We also added the contour lines of the jet because one thing that you may have asked is, OK, but does the jet exist in the constant models? Because so far you saw just a constant thing, but is the jet there? Yes, it's there. It's there. We plotted it. Then we added black horizontal and vertical lines to just separate a little bit better the simulations. And also, for the sake of this talk, I just put a dark background color, because I think it fits nicer with a dark background slide. And this process is the lesson number four. The simple task, you should aim for effortless. You should aim for trivial. The complex task, you should aim for simple. A complex task will likely not really be trivial, But it can be simple. And since, now I'm sure you're going to say, Giancarlo, this is amazing, and it's script-based. Why do you need a graphics user interface? Well, imagine it's 8 AM, and you are before the time where you take coffee, and you start writing your code. You have no coffee, and then you start putting bugs. Now, unfortunately, you don't really see it very well. Let me just do this then. Let's do a bit impromptu. OK, much better. Now, I highlighted a couple of bugs that now I can challenge you to see. But you can see that there are some user side bugs, like numpy imported as nop, that your finger can slip and just write an extra letter. Or my folder has the quotes that are not opened. So a great user interface helps with fewer user side bugs. Less cognitive load, you don't really have to remember what you should do. You press some buttons and the thing should work. Visual intuition, where's the folder? Well, you can search on a window. Which variable do you show? You will have a window menu. And also, you don't really have to remember the API, because if something is deprecated, we should take care of that when we build a graphic user interface. You shouldn't take care of that. So a graphic user interface works to have an even safer and more controlled work. Let me just show you how it works, how you select a file. This is a video that is also present in the documentation. I'm going to go there in two slides. You select the quantity, the pressure. You select the scale for your color map, and you plot it. And then you decide that you want to change a little bit your color map. you want a divergence, because why not? Ta-da! And you do it. And this is all by clicking button and choosing variables. So, it's very, very intuitive. Also, this, again, anecdotal reference. It was 2021. I started taking over PyPluto, and then I had a really heated discussion with some collaborators, because they wanted to build their own tool that was based on their, like, an obscure language that was not integrated with Python, and I was really fuming because I thought, like, Python is not just a programming language. Python is a community. Because there's a lot of people that use NumPy, there's a lot of people that use Matplotlib, and you should integrate with the people and share the burden and share all the improvements that are happening on a daily, weekly, monthly basis. So we use We use numpy arrays, we use matplotlib, and we use also a lot of other Python packages because we don't want to reinvent the wheel, and we don't want to alienate ourselves from the community. And also, it goes without saying, we try to test our code. And we try to make it Pythonic, we use Pyrefly, for example, we use Ruff, and we use all the type of things that make the code more maintainable because we don't want to rewrite it completely in six months. want the package to stay, and we want the next generation of astronomers and students, when they're going to take over, they're going to have a better package that was the one we took over from. Now, final lesson is, how do people learn to use a tool? Do they read the source code? God, no. Do they read the documentation? Maybe. They look at the example. They mostly look at the example. And that's why we have a documentation. But most important, we have a lot of tests. And these tests cover, I made a small routine that tests, that uses PyTest with this example. It covers more than 70% of the code. So these tests are very extensive. And this is the last lesson, let's say. Examples are the first thing users see and care about. So let me do a little bit of show and don't tell. This is the first test. Very intuitive. You load the data, you load the image, you plot the data, you save the image. End of story. A little bit more complicated, like some subplots, a color bar at the left. You load the data, multiple loading. You create the image, you create the subplots, display, color bar, save the figure. And this has some field lines computed with Matplotlib or with SciPy. So we load the data, we create the image, we compute the field lines, we display the the magnetic field with the colour, string plot. We plot the field lines that we computed with SciPy, and save fig. As you see, we try to make the code compact. Short. You don't have to write too much. This has two insert zooms, one of different quantities. We load the data, display the density on background, zoom, zoom, save fig. End of story. This is an animation, and it's literally like five instructions. You create the image, interactive, it. And that's it. This is a little bit more complicated. But it has a lot of things going on. But I have to be honest, at some point, the code becomes a little bit more cumbersome. But this is a publication figure. So I think it's a little bit more justified. Now, let me tell a little bit about what's the outlook. So what's the next level of simulation? Well, I was discussing with some friends and with the great chair, and we thought, like, what What about something that loads and outputs the data during the simulation? It creates a lot of pros. Easy debug, you avoid large files, you don't have to save things, and for GPU this is really a bottleneck. It becomes also more accessible because you can embed this into a good UI and you're done. The cons is that it's not trivial to build. So me and Francesco, who is the PI, decided to build Karon, which is a Python package for runtime visualization on GPU. And the one on the left, Thibault, is a PhD student at the Max Planck Institute for Nuclear Physics. I don't know what he was thinking at that moment, but I'm sure if he was listening about us talking to Karen, he would have said, wow, it sounds really cool. So let me show a little bit of an example. This is a Kelvin-Helmholtz instability. So this is a fluid instability that happens when you have fluids moving at different velocities. And if you perturb a little bit the interface, your vortices will start to appear. This is a simulation of 512 squared grid cells that I ran on my laptop, which has a GPU. And since we used Jack and the AirPi guy, what happens is that we reached 60 FPS, which means that 60 times per second, we were solving the equation of hydrodynamics and plotting the results at a resolution of 500 squared cells. The resolution is really remarkable because at some point you will see, you really see some secondary vortices that are starting to form. This one. This happens only when the resolution is sufficiently high, which means that we are really reducing the numerical dissipation. The resolution is very, very high and still it is able to run on a laptop. I find it really fascinating. So in conclusion, if you're interested about PyBluetooth processing data from astrophysical simulation, please check the paper. We managed to publish it on the Journal of Open Source Software. By the way, shout out to them. Great journal. Totally recommend. And yeah, about that, thank you very much.

Speaker 1 [25:13]

so Thank you for that

Speaker 2 [25:16]

that.

Speaker 1 [25:17]

We are all very tired. We have only one question from a very cool guy. Wouldn't building a matplotlib wrapper limit you in case you want later to utilize some visualization feature that you didn't plan for in the beginning? Is it worth the risk?

Speaker 2 [25:33]

I think we try to circumvent this by taking the venture to the fact that we are always, we always have figure and access accessible. So, if you are, this is also more anecdotal prefer. But at the same time, it really works from, it really comes from the feedback. Because the target that we have is broad, but depending on who uses, we can tailor it and make it maybe simpler and a little bit more contained or more flexible. I think we reached a decent compromise by letting the user access the figure and the access to the method. So if you want to do something that is not encompassing the capabilities of PyGluto, you can definitely do it. But at the same time, I think that if you are a pro of matplotlib, you very likely will have your own routine, and you will not rely on a tool that you haven't built yourself. So if I have to choose between experts and beginners, I will always focus on the beginners. And I think this is a very important thing to do, because you will not need help, they Thank you.

Speaker 1 [26:52]

Was I ready? No. But I think we don't have another question. No, we finished all the questions. So that's it. Thanks a lot, Giancarlo. Thank you all for making it this far. And in ten minutes, no, five minutes, there is a closing ceremony downstairs. All right?

Speaker 2 [27:12]

No idea. 16, 12.