Tracing field lines or streamlines with yt in an analysator script

This page summarizes how one can use the yt library streamline tracing function with Vlasiator data. This tutorial page is based on the relevant section of the 3D magnetopause identification script.

Installing yt

As in Vorna guide one can pip-install yt, assuming you have the relevant python version/module/environment loaded

export PYTHONUSERBASE=/proj/$USER/.local
export PYTHONPATH=$PYTHONPATH:/proj/$USER/.local
pip3 install yt --cache-dir /proj/$USER/.cache --user

Reading in data

yt can load data on a uniform grid very simply, you need to provide it the data in a dictionary, each component in one entry as a numpy array, as well as a bounding box. Use f.read_variable() for a uniform grid case, or upsample to a fine grid using e.g. f.read_variable_as_fg() to map an AMR grid to the fg grid.

Once you have the components in the dictionary, you load the data into yt with

    yt_dataset = yt.frontends.stream.load_uniform_grid(
        data, # data dictionary
        sizesfs, # number of points per dimension
        bbox=np.array([[boxcoords[0], boxcoords[1]],
                    [boxcoords[2],boxcoords[3]],
                    [boxcoords[4],boxcoords[5]]])) # bounding box (determines coordinate system)

Tracing lines

Prepare a numpy array of seed points in the coordinate system defined by your bbox, when you have that ready you can call the setup function before calling the tracing function itself. The latter takes time and displays so on the command line. You get the streamline_pos (in our example) object that contains the streamline_pos.streamlines object that can be iterated over, each element corresponding to the next seed. streamline_pos.streamlines[s] contains the array of consecutive coordinates along the streamline, which you can use or export. Some cleaning up e.g. of [0,0,0] entries can be done.

Tips

Parallelization

You can pass a nprocs=N argument to yt.load_uniform_grid() if you have called yt.enable_parallelism() before that and your script is called in a parallel instance with N tasks. Streamline tracing scales accordingly, if you have sufficiently many seeds/streamlines to process. Do add a

   from mpi4py import MPI
   MPI.COMM_WORLD.Barrier()

after the call(s) to streamlines_pos.integrate_through_volume() just to make sure you don't break things.

Speeding things up

If you trace thousands of lines for long distances, it can easily take minutes to an hour or more. By default, it traces at "good" resolution/quality. Passing dx=x to yt.visualization.api.Streamlines() where x is a multiple of the grid resolution can speed up things dramatically and in Yann's experience the results are still pretty good (tested looking at four-field junction with a grid resolution of 1e6 m, I got same results going even to dx=6e6 or so.

Of course parallelization might help, though if your script also hopes on threading, good luck, Yann didn't manage to get a script doing both MPI/task parallelization in one section (the above) and threading (pool of threads) elsewhere.