xmsgridtrace

Introduction

xmsgridtrace provides functionality to trace a particle through a vector field on an unstructure grid (xmsgrid). This notebook demonstrates basic functionality of this library.

from shapely.geometry import Polygon
import holoviews as hv
import pandas as pd

import numpy as np
import geoviews as gv
import xarray as xr
import cartopy.crs as ccrs

import hvplot.pandas
import hvplot.xarray

import math
import xmsgridtrace
import xmsgrid

import gridtrace_tools

hv.extension("bokeh")

Create UGrid

First, we will create a simple unstructured grid.

points = [
    (0, 0, 0), (1, 0, 0), (2, 0, 0), 
    (0, 1, 0), (1, 1, 0), (2, 1, 0), 
    (0, 2, 0), (1, 2, 0), (2, 2, 0)
]
cells = [
    xmsgrid.ugrid.UGrid.QUAD, 4, 0, 1, 4, 3, 
    xmsgrid.ugrid.UGrid.QUAD, 4, 1, 2, 5, 4, 
    xmsgrid.ugrid.UGrid.QUAD, 4, 3, 4, 7, 6,
    xmsgrid.ugrid.UGrid.QUAD, 4, 4, 5, 8, 7,
]

ug = xmsgrid.ugrid.UGrid(points, cells)

Display UGrid

%opts Polygons [width=700 height=700]
ug_polygons, ug_lines, ug_points = gridtrace_tools.create_hv_objects_from_ugrid(ug)
ug_polygons = ug_polygons.redim.range(x=(-0.5, 2.5), y=(-0.5, 2.5))
ug_polygons * ug_lines * ug_points

Create Tracer

tracer = xmsgridtrace.gridtrace.GridTrace(
    ugrid=ug,
    vector_multiplier=1,
    max_tracing_time=-1,
    max_tracing_distance=-1,
    min_delta_time=0.01,
    max_change_distance=-1,
    max_change_velocity=-1,
    max_change_direction_in_radians=(math.pi/4)
)

Scalars at Time Steps

This next step will define the vectors on the ugrid. Notice that a "time" is associated with each set of vectors.

scalars = [
    [
        (0, 2, 0), (-0.2, 0, 0), (-2, 0, 0), 
        (0, 0.2, 0), (0, 0, 0), (0, -0.2, 0),
        (2, 0, 0), (0.2, 0, 0), (0, -2, 0),
    ],
    [
        (0, 4, 0), (-0.4, 0, 0), (-4, 0, 0), 
        (0, 0.4, 0), (0, 0, 0), (0, -0.4, 0),
        (4, 0, 0), (0.4, 0, 0), (0, -4, 0),
    ],
]
activity = [True for _ in range(0, len(scalars[0]))]

# Time 0.0
tracer.add_grid_scalars_at_time(
    scalars=scalars[0], 
    scalar_loc="points",
    cell_activity=activity, 
    activity_loc="points",
    time=0.0
)

# Time 20.0
tracer.add_grid_scalars_at_time(
    scalars=scalars[1], 
    scalar_loc="points",
    cell_activity=activity, 
    activity_loc="points",
    time=20.0
)

# draw the vectors
%opts Polygons [width=700 height=700]
vf_scalars = gridtrace_tools.hv_vector_field_from_2_timesteps(
    scalars[0], points
)
ug_polygons * ug_points  * ug_lines * vf_scalars

Trace Point

traces, times = tracer.trace_point((0.5, 0.5, 0), 0.0)
print(traces)
print(times)

[[0.5        0.5        0.        ]
 [0.5        1.25       0.        ]
 [0.54457813 1.33915625 0.        ]
 [0.61632493 1.43549847 0.        ]
 [0.72535406 1.53155337 0.        ]
 [0.88236797 1.61268018 0.        ]
 [0.98873181 1.6331016  0.        ]
 [1.05385039 1.6342606  0.        ]
 [1.1249433  1.56830068 0.        ]
 [1.18950974 1.38634489 0.        ]
 [1.22352421 1.05885901 0.        ]
 [1.22352421 0.90477286 0.        ]
 [1.20053362 0.85080764 0.        ]
 [1.15817907 0.7938777  0.        ]
 [1.08968746 0.74131697 0.        ]
 [0.98966251 0.70663753 0.        ]
 [0.9580615  0.7181798  0.        ]
 [0.92629621 0.77371504 0.        ]
 [0.90239413 0.88917318 0.        ]
 [0.89995173 1.06948757 0.        ]
 [0.91503139 1.09119928 0.        ]
 [0.93816745 1.1127547  0.        ]
 [0.97140029 1.13097896 0.        ]
 [0.99364913 1.13583707 0.        ]
 [1.00715245 1.1364684  0.        ]
 [1.02234471 1.12806558 0.        ]
 [1.03697378 1.0971462  0.        ]
 [1.04673977 1.03713772 0.        ]
 [1.04673977 0.96499504 0.        ]
 [1.03905762 0.95473758 0.        ]
 [1.02764446 0.94488899 0.        ]
 [1.02087912 0.9414954  0.        ]]
(0.0, 0.375, 0.825, 1.3649999999999998, 2.013, 2.7905999999999995, 3.2571599999999994, 3.537095999999999, 3.873019199999999, 4.276127039999999, 4.759856447999998, 5.340331737599998, 6.036902085119998, 6.872786502143997, 7.875847802572797, 9.079521363087355, 9.440623431241724, 9.873945913026965, 10.393932891169255, 11.017917264940003, 11.7666985134649, 12.665236011694777, 13.743481009570628, 14.390428008296139, 14.778596207531445, 15.244398046613812, 15.803360253512654, 16.474114901791264, 17.279020479725595, 18.244907173246794, 19.40397120547223, 20.0)

View Traces

hv_traces = hv.Points(traces).options(color='red', size=4)
ug_polygons * ug_points  * ug_lines * vf_scalars * hv_traces