Hello,
I’ve tried a few approaches to read ADIOS2 BP format timeseries data from the X-Point Gyrokinetic Code into Paraview for visualization (I’d like to eventually play with VR, and I’d like to represent the timeseries data perhaps through an animation or something). Right now, I’m struggling to conceptualize the best way to associate mesh data (written to an ADIOS2 BP “file” (really a folder) called xgc.mesh.bp). The timeseries data is a variable in files with the naming scheme xgc.3d.XXXXX.bp, where XXXXX refers to the timestep.
For background, XGC uses a simple cylindrical coordinate system (R, Z, phi → toroidal angle), and data is (often) computed on a 3D mesh that is unstructured in R,Z planes and evenly spaced / structured in the phi / toroidal angle direction (think - donut shaped mesh, with unstructured mesh cross sections that are cut out of the donut in even intervals in such a way that the knife passes through the hole in the center in order to get the “poloidal plane” cross section). This data is written out to BPs as a 2D array, with the first dimension containing data associated with all of a poloidal plane’s / cross section’s unstructured mesh points, and the second dimension representing the poloidal plane index (again, these are evenly spaced in phi from 0-2pi).
So, with that, I need to be able to associate this R,Z data with mesh point indices, and compute an evenly spaced phi, then transform to (x,y,z) coordinates and associate these mesh locations to the data in the timeseries BP files for use in paraview. From my reading, I see a few options:
- Write an XML file into the BP folder for each BP (later, if this works, specify it as an attribute written to the BP file) that specifies this data as a VTK “unstructured grid”. It seems like in this case, I have to include the mesh point R,Z,phi or x,y,z locations in a text array of values - which would contribute quite a bit of bloat / duplication (I have ~1500 timestep BP files, and will likely produce more).
- Write some sort of top-level (in the base folder rather than within each BP file) specification file (XML? XDMF?) that I don’t understand well at all, and I’m not sure where I can read up on this.
- Read data and rewrite to a file that VTK inherently understands (probably, with the mesh data saved into each file), which would be impractical since I have very large datasets.
- Write a plugin in python that reads from the xgc.mesh.bp file and computes / loads triangle data as the data is loaded from the timeseries file. I have an example that ChatGPT helped me create below (seems right, but I’m not fully sure), but it’s not clear to me how to leverage the ADIOS2 that Paraview is using under the hood. I’ve attempted to use pip to install adios2 into the correct folder as in here, but it seems that the pip adios2 library is potentially broken (I can’t import adios2 on my miniconda install - I have to use conda to install to get an importable install - nor is it able to be imported from the paraview python interpreter shell). I’ve also tried installing to my miniconda install using conda, then copying the package folder manually to the folder where paraview stores its packages internally, but I can’t import the package properly that way either.
So overall, I’m looking for:
- What’s the best way to get the data into paraview, associating the R,Z poloidal plane mesh with the data contained in xgc.mesh.bp?
- If a plugin is potentially the best method, how would I go about getting a usable adios2 python package working so I can use my below plugin? Does the plugin below seem ok for this use?
Thanks all.
Plugin:
XGC Reader Plugin
from paraview.util.vtkAlgorithm import smproxy, smproperty, smdomain
import numpy as np
import adios2
import vtk
import os
from vtk.util import numpy_support as ns
@smproxy.reader(name="XGCReader", label="XGC BP Reader",
extensions="bp", file_description="Read XGC ADIOS2 timeseries files (xgc.3d..., xgc.f3d...)")
class XGCReader(vtk.vtkAlgorithm):
def __init__(self):
super().__init__()
self.SetNumberOfInputPorts(0)
self.SetNumberOfOutputPorts(1)
self.filename = None
self.rz_mesh = None
self.n_phi = 0
self.phi_vals = None
self.triangles = None
@smproperty.stringvector(name="FileName")
@smdomain.filelist()
def SetFileName(self, fname):
self.filename = fname
self.Modified()
def RequestInformation(self, request, inInfo, outInfo):
return 1
def RequestData(self, request, inInfo, outInfo):
# Extract timestep number to open mesh file
timestep = int(self.filename.split('.')[-2]) # relies on the naming scheme "xgc.{3d,f3d,2d,f2d}.XXXXX.bp"
# Read mesh (once)
if self.rz_mesh is None:
with adios2.open("xgc.mesh.bp", "r") as f:
self.rz_mesh = np.array(f.read("rz")) # shape: [npoints, 2]
with adios2.open(self.filename, "r") as f:
tmp_data = f.read("dpot") # just read an array of data that's going to be present - dpot is NOT present in any
self.n_phi = tmp_data.shape[1] # other type of file, but that usually doesn't have the same layout of data anyway
self.phi_vals = np.linspace(0, 2 * np.pi, self.n_phi, endpoint=False)
self.triangles = np.array(f.read("nd_connect_list")) # shape: [ntris, 3]
n_points_2d = self.rz_mesh.shape[0]
n_phi = self.n_phi
# Build 3D grid
points = vtk.vtkPoints()
for i_phi, phi in enumerate(self.phi_vals):
cos_phi, sin_phi = np.cos(phi), np.sin(phi)
for r, z in self.rz_mesh:
x = r * cos_phi
y = r * sin_phi
points.InsertNextPoint(x, y, z)
# Build triangles
cells = vtk.vtkCellArray()
for i_phi in range(n_phi):
offset = i_phi * n_points_2d
for tri in self.triangles:
ids = [int(idx + offset) for idx in tri]
triangle = vtk.vtkTriangle()
for j in range(3):
triangle.GetPointIDs().SetID(j, ids[j])
cells.InsertNextCell(triangle)
# Read scalar field
with adios2.open(self.filename, "r") as f:
varname = list(f.available_variables.keys())[0]
data = np.array(f.read(varname)) # shape: [npoints_2d, nplanes]
# Transpose to [nplanes, npoints_2d] for easier looping
data = data.T.flatten()
# Build output polydata or unstructured grid
output = vtk.vtkUnstructuredGrid()
output.SetPoints(points)
output.SetCells(vtk.VTK_TRIANGLE, cells)
# Add scalar data
vtk_array = ns.numpy_to_vtk(data, deep=True)
vtk_array.SetName(varname)
output.GetPointData().AddArray(vtk_array)
output.GetPointData().SetActiveScalars(varname)
outInfo.GetInformationObject(0).Set(vtk.vtkDataObject.DATA_OBJECT(), output)
return 1