Getting Started

This tutorial will take you step by step through the 3D : Visualization example.

install / prequisites

Make sure you have pygranite setup correctly.

write the script

Create a new python file somewhere (in your virtual environment if you use pipenv) and add the following in it:

#!/usr/bin/env python3

import vedo, numpy, os
import pygranite

from netCDF4 import Dataset

The first line specifies this as a python script (under UNIX systems it allows the script to be run like a program). Further we add the imports for vedo (the visualization library), numpy, os and pygranite. We do as well import the Dataset class from the netcdf library to allow us to fetch our data from *.netcdf files.

writing the data loader

Our first goal is to implement our instance of a DataLoader, which has a barebone structure like this:

class netCDFLoader(pygranite.DataLoader):
    def __init__(self):
        pygranite.DataLoader.__init__(self)  # do not use super()

    def step(self):  # interface method
        pass

    def windfield(self):  # interface method
        pass

    def uplift(self):  # interface method
        pass

    def additionalVolumes(self):  # interface method
        pass

    def constants(self):  # interface method
        pass

The reason this has to be done is that it is usually impossible to store an entire set of windfields on the gpu. The class is used to provide the data to the algorithm such that it can load the data on demand (Note: If your PC has enough main memory (CPU-RAM) to store the entiry dataset, you can even speedup the loading process as well, if you want to compute multiple TrajectorySets on the same dataset).

We want to store a list of the filenames to the dataset, so we adjust our constructor to take an additional argument files, which will be a python list containing the filenames. We also add an additional member _index which will store the index to the next dataset, which is needed next (Note that the underscore indicates (by convention) that the variable is private).

class netCDFLoader(pygranite.DataLoader):
    def __init__(self, files):
        pygranite.DataLoader.__init__(self)  # do not use super()
        self._files = files[:]  # copy
        self._index = 0

    def step(self):  # interface method
        pass

    def windfield(self):  # interface method
        pass

    def uplift(self):  # interface method
        pass

    def additionalVolumes(self):  # interface method
        pass

    def constants(self):  # interface method
        pass

The first method step will be called by the algorithm request an increment to our internal structure (if it returns False the execution will be stopped). In our case this should simply be the case if our _index is less than the available windfields.

    def step(self):  # interface method
        if self._index < len(self._files):
            self._index = (self._index + 1)
            return True
        return False

In windfield we first we have to get our next file:

    def windfield(self):  # interface method
        fn = self._files[self._index]

Now that we have the filename we load the netCDF4.Dataset and extract the required information.

    def windfield(self):  # interface method
        fn = self._files[self._index]
        ds = Dataset(fn, 'r')
        u = ds.variables['u'][:][0] # depends on the files
        v = ds.variables['v'][:][0] # depends on the files
        w = ds.variables['w'][:][0] # depends on the files

And finally return the value of these in form of a list:

    def windfield(self):  # interface method
        fn = self._files[self._index]

        ds = Dataset(fn, 'r')
        u = ds.variables['u'][:][0]
        v = ds.variables['v'][:][0]
        w = ds.variables['w'][:][0]
        
        return [u, v, w]

Note that the windfield data is expected to be addressed like u[z][y][x]. If we are in a 2-dimensional case we simply return a two-dimensional list.

The final loader should look something like this (with an additional reset-method):

class netCDFLoader(pygranite.DataLoader):
    def __init__(self, files):
        pygranite.DataLoader.__init__(self)  # do not use super()
        self._files = files
        self.reset()

    def step(self):  # interface method
        if self._index < len(self._files):
            self._index = (self._index + 1)
            return True
        return False

    def reset(self):
        self._index = 0

    def windfield(self):  # interface method
        fn = self._files[self._index]

        ds = Dataset(fn, 'r')
        u = ds.variables['u'][:][0]
        v = ds.variables['v'][:][0]
        w = ds.variables['w'][:][0]
        
        return [u, v, w]

Note that the loader used in the example file is a little bit different to demonstrate possible other use-cases.

Now we are ready to continue with our main program.

implementing the simulation

With the DataLoader done, the hardest part is done. We start our main script by adding the files location as a variable (to eliminate calls from different directories messing up our relative pathing).

if __name__ == "__main__":
    root = os.path.abspath(os.path.dirname(__file__))  # directory the file is in

Then we initialize our self-written WindfieldLoader with a dataset (example path relative to the file):

    directory = root + "/../data/2020/fields_scaled/"
    loader = netCDFLoader([
        directory + f"{i}.nc" for i in range(len(os.listdir(directory)))
    ])

We then configure our simulation accordingly:

    # simulation properties
    settings = pygranite.IntegratorSettings()
    settings.Space = pygranite.Space.Space3D
    settings.MinimumAliveParticles = 16
    settings.WindfieldMode = pygranite.WindfieldMode.Dynamic
    settings.Integrator = pygranite.Integrator.ClassicRungeKutta
    settings.DeltaT = 1 
    settings.DataTimeDistance = 60
    settings.MaximumSimulationTime = 1000
    settings.GridScale = [25, 25, 25]
    settings.SaveInterval = 1

Our simulation space is 3D (as its a 3D windfield loader), we want to check each time if at least 16 particles are alive. We have cartesian coordinates, we want to interpolate between the windfields over time, use the default integrator and set the correct timings. The GridSize property specifies the distance between the datapoints in the windfields (in this case 25 (meters)).

Now we have to generate the starting positions of the particles (here done inside a specified spawn box):

    spawn_min = numpy.array([2400, 2530, 150]) # the spawn domain
    spawn_max = numpy.array([2520, 2785, 170])

    inp = numpy.random.rand(512, 3)
    for i in range(inp.shape[0]):
        inp[i] = ((spawn_max - spawn_min) * inp[i] + spawn_min)
    start_set = pygranite.TrajectorySet(inp)

Finally we pass everything to an integrator instance and compute the trajectories:

    # computation
    integrator = pygranite.TrajectoryIntegrator(settings, loader, start_set) 
    result_set = integrator.compute()

visiualizing trajectories

The result_set is a TrajectorySet instance providing some basic operations (and resuse of final positions - to pause computation for example). It provides a method cloud() returning all trajectory points in order, which makes it possible to color them accordingly:

    points = result_set.cloud()

    colors = [
        i / result_set.numberTrajectories() \
            for i in range(result_set.numberTrajectories()) \
                for _ in range(result_set.lengthTrajectories())
    ]

For the visualization we create a domain box (depending on the provided data) and pass the points to the vedo library:

    # visualization using vedo
    domain = [settings.GridScale[0] * 320, settings.GridScale[1] * 192, settings.GridScale[1] * 120]
    world = vedo.Box(list(map(lambda x: x / 2, domain)),  # box origin = box center
                     size=tuple(domain)).wireframe().color([1, 0, 0]).lineWidth(4).flat()
    pts = vedo.Points(points).legend("trajectory points") \
        .cmap("viridis", input_array=numpy.array(colors))
    vedo.show(pts, world)

The final script should look something like this:

#!/usr/bin/env python3

import vedo, numpy, os
import pygranite

from netCDF4 import Dataset


class netCDFLoader(pygranite.DataLoader):
    def __init__(self, files):
        pygranite.DataLoader.__init__(self)  # do not use super()
        self._files = files
        self.reset()

    def step(self):  # interface method
        if self._index < len(self._files):
            self._index = (self._index + 1)
            return True
        return False

    def reset(self):
        self._index = 0

    def windfield(self):  # interface method
        fn = self._files[self._index]

        ds = Dataset(fn, 'r')
        u = ds.variables['u'][:][0]
        v = ds.variables['v'][:][0]
        w = ds.variables['w'][:][0]
        
        return [u, v, w]


if __name__ == "__main__":
    root = os.path.abspath(os.path.dirname(__file__))

    # initialize windfield loader
    directory = root + "/../data/2020/fields_scaled/"
    loader = netCDFLoader([
        directory + f"{i}.nc" for i in range(len(os.listdir(directory)))
    ])

    # simulation properties
    settings = pygranite.IntegratorSettings()
    settings.Space = pygranite.Space.Space3D
    settings.MinimumAliveParticles = 0
    settings.WindfieldMode = pygranite.WindfieldMode.Dynamic
    settings.Integrator = pygranite.Integrator.ClassicRungeKutta
    settings.DeltaT = 1 
    settings.DataTimeDistance = 60
    settings.MaximumSimulationTime = 1000
    settings.GridScale = [25, 25, 25]
    settings.SaveInterval = 1

    # generate start positions
    spawn_min = numpy.array([2400, 2530, 150])
    spawn_max = numpy.array([2520, 2785, 170])

    inp = numpy.random.rand(512, 3)
    for i in range(inp.shape[0]):
        inp[i] = ((spawn_max - spawn_min) * inp[i] + spawn_min)
    start_set = pygranite.TrajectorySet(inp)

    # computation
    integrator = pygranite.TrajectoryIntegrator(settings, loader, start_set) 
    result_set = integrator.compute()

    # visualization using vedo
    points = result_set.cloud() # .trajectory(0) # .cloud() 

    colors = [
        i / result_set.numberTrajectories() \
            for i in range(result_set.numberTrajectories()) \
                for _ in range(result_set.lengthTrajectories())
    ]


    domain = [settings.GridScale[0] * 320, settings.GridScale[1] * 192, settings.GridScale[1] * 120]
    world = vedo.Box(list(map(lambda x: x / 2, domain)),  # box origin = box center
                     size=tuple(domain)).wireframe().color([1, 0, 0]).lineWidth(4).flat()
    pts = vedo.Points(points).legend("trajectory points") \
        .cmap("viridis", input_array=numpy.array(colors))
    vedo.show(pts, world)