Advanced visualization workflows for gyrokinetics simulation data

1 Introduction

The feasibility of tokamaks as reliable, large-scale sources of fusion energy depends on their ability to confine plasma effectively. Tokamaks use strong magnetic fields to establish the required confinement conditions; however, physical processes within the plasma act against confinement. One of the most important of these processes is turbulence. Plasma turbulence drives cross-field transport of heat, particles, and momentum, degrading confinement by reducing confinement time, pressure, and overall performance [1, 2]. Understanding and ultimately mitigating turbulence-driven transport is therefore a central challenge in the design of viable fusion reactors.

Although turbulence can be studied experimentally, acquiring the breadth of measurements needed for a robust understanding remains technically challenging and is often limited by diagnostic access. As a result, much of the field relies on mathematical and numerical investigations based on advanced, high-fidelity simulations. Gyrokinetic modeling is the dominant approach to quantifying turbulent transport in fusion reactors and identifying its principal drivers [2, 3]. Codes such as GKW [4], GENE [5], CGYRO [6], and GS2 [7] have been developed, extensively verified through inter-code benchmarks [8, 9], and are now routinely used to simulate microturbulence under realistic tokamak conditions. Beyond these, surrogate modeling has emerged as an increasingly important complement to full gyrokinetic simulation, particularly in the context of integrated plasma scenario modelling, where reduced or surrogate descriptions can be used inside larger coupled workflows such as JINTRAC-based scenario studies [10]. Paischer et al. [11] recently introduced GyroSwin, a scalable five-dimensional neural surrogate for nonlinear electrostatic flux-tube gyrokinetic distribution functions with adiabatic electrons, reducing the computational cost of fully resolved simulations by three orders of magnitude while remaining physically verifiable against the reference codes within a given input parameter domain.

In a strongly magnetized plasma, the motion of a charged particle can be separated into fast gyromotion around a magnetic field line and slow drift dynamics. In principle, this motion is described by the coupled Vlasov and Maxwell equations. However, the Vlasov equation is six-dimensional, and resolving the cyclotron frequency explicitly is computationally prohibitive. Gyrokinetics removes the fast gyromotion by representing the particle at its gyrocenter and averaging over the gyrophase, reducing the problem to a five-dimensional formulation [2, 4]. The resulting phase-space distribution function describes the gyrocenter density in a five-dimensional phase space containing three spatial and two velocity coordinates. This structure makes direct visualization of the raw simulation output highly challenging [12, 13]. Simulation codes typically post-process the distribution function by computing velocity-space moments, yielding field variables and fluid quantities that can be represented in three-dimensional physical space. These quantities are often reduced further to zero-dimensional (0D) transport outputs, such as fluxes, that can be passed downstream to larger integrated modelling workflows. The spatial quantities considered here include the electrostatic potential ϕ and fluctuations in density δ n and temperature δ T. Although visualizing the full five-dimensional distribution function is an important open problem, here we focus on the three-dimensional post-processed spatial fields rather than on the scalar reduced outputs.

Even within the reduced three-dimensional setting, producing useful visualizations of gyrokinetic data is non-trivial. The simulation domain is defined in magnetic-field-aligned flux-tube coordinates that follow the field line around the torus [14], and mapping these coordinates into real physical space requires a substantial geometric transformation [15]. As a result, gyrokinetics specialists have traditionally worked with only a small subset of the output data, such as one-dimensional time traces or fields examined near the outboard mid-plane. Whole-dataset three-dimensional visualization can provide scientific value, but the required processing has limited its routine use.

Previous work has explored a range of approaches to visualizing gyrokinetic simulation data in three-dimensional space. Crawford et al. [16] introduced interactive volume visualization of plasma flow in toroidal coordinate systems using custom GPU code, demonstrating that existing commercial software packages were inadequate for the non-rectilinear mesh structures produced by gyrokinetic codes. Stantchev et al. [12] developed a confluent volumetric rendering framework that simultaneously encodes two physically meaningful fields, namely, heat flux and electric potential, within a single volumetric representation, providing a direct visual connection between transport and the driving turbulent fluctuations. Jones and Ma [17] addressed the challenge of time-varying multivariate particle data from particle-in-cell gyrokinetic simulations, arguing that traditional x-y plots and contour plots are insufficient for the multi-dimensional quantities of interest and that interactive multidimensional visualization can reveal aspects of the data that were simply not visible before. More recently, Feibush et al. [15] presented three complementary geometric techniques for visualizing GTS simulation data that work directly with the topology of magnetic flux surfaces, producing high-resolution surface models for display in scientific visualization programs (e.g., ParaView).

ParaView is a widely used open-source scientific visualization platform built on VTK [18–20]. It provides a powerful and flexible environment for visualizing large datasets, including support for volume rendering and advanced surface rendering. ParaView’s Python interface (pvpython [18, 19]) allows users to create reproducible visualization pipelines that can be shared and automated. There has also been work demonstrating the value of advanced visualization tools not originally designed for scientific visualization, such as Blender [21–23] and NVIDIA Omniverse [24, 25], in particular when it comes to integrating heterogeneous fusion data [26, 27].

Beyond fusion, there is a substantial body of work on multivariate data visualization. Multivariate visualization seeks to combine information from multiple related variables within a single representation, improving discoverability and helping to identify relationships among the variables being represented [28]. High-dimensional scientific datasets pose a general challenge in this regard. Physical display limitations and human visual system constraints prevent direct rendering of structures with more than two or three dimensions, and the most effective strategy depends strongly on the nature and structure of the data [13, 28]. There is therefore no universal solution for multivariate datasets, and the approach must be tailored to the specific physical quantities and use cases of interest.

The aim of this paper is to develop reproducible workflows for visualizing gyrokinetic simulation data in three-dimensional real space and to provide a practical basis for further work on this class of data. We address the lack of tools that span the full pipeline from raw simulation output to publication-ready renders or deployment-ready experiences that natively support multiple gyrokinetic codes. We illustrate the workflow through multivariate visualizations of potential, density, and temperature fluctuations and by applying it across multiple gyrokinetic codes and visualization environments. The remainder of the paper is structured as follows: Section 2 describes the coordinate transformation paths, data processing pipelines, and interactive viewer; Section 3 presents representative visualization results across multiple datasets; and Section 4 discusses the implications, current limitations, and directions for future development.

2 Methods

2.1 Gyrokinetics overview

Gyrokinetics is a reduced kinetic description for strongly magnetized plasmas based on a standard ordering in which fluctuation frequencies remain small compared with the cyclotron frequency, fluctuation amplitudes are small, gyroradii are much smaller than equilibrium scale lengths, and the turbulence is anisotropic, with much shorter perpendicular than parallel scales ( k ⊥ ≫ k ‖ ) [2, 29]. Under this ordering, the fast gyromotion can be averaged away while retaining the low-frequency microturbulence and transport dynamics of interest.

In many local simulations, a common numerical setting is the flux-tube approximation, where the perpendicular domain surrounding a reference magnetic field line is taken to be small compared with the equilibrium variation scale, so equilibrium quantities and gradients are treated as approximately constant across the simulation domain [14, 29]. This motivates the use of field-aligned coordinates: one coordinate follows the magnetic field (often a straight-field-line angle θ), while the perpendicular directions represent radial and binormal displacements (often treated spectrally). In many cases the electrostatic potential ϕ is a primary fluctuating field that mediates the dominant E × B drift and drives cross-field transport, while electromagnetic simulations also evolve magnetic perturbations such as A ‖ and δ B ‖. Fluctuations in density and temperature are obtained either directly as evolved fields in reduced models or, more generally, as velocity-space moments of the perturbed distribution function (e.g., δ n s and δ T s for species s). To visualize ϕ and derived quantities such as density and temperature fluctuations in physical space, the flux-tube domain must be embedded in toroidal geometry and mapped into Cartesian coordinates, as described below.

2.1.1 Coordinate transformation

We employ two related but distinct pipelines to map field-aligned computational data onto real-space Cartesian coordinates for visualization and volumetric output. In the spectral-coordinate pipeline, the perpendicular field is reconstructed in real space using Equation 1, after which the binormal coordinate and toroidal angle are defined by Equations 2, 3, respectively. Parallel closure is imposed through Equation 4, the equilibrium geometry is specified by Equation 5, and the resulting cylindrical coordinates are converted to Cartesian coordinates using Equation 6. In the gridded-format pipeline, the parallel coordinate and radial coordinate are mapped according to Equations 7, 8, respectively, while the geometric toroidal angle is constructed using Equation 9. The analytically shaped flux surface is then defined by Equations 10, 11. Thus, both pipelines involve two principal stages: first, the computational domain is embedded in toroidal geometry; second, the resulting cylindrical coordinates are converted to Cartesian coordinates. The pipelines differ in the native coordinates of the input fields, the source of the geometric embedding, whether equilibrium-derived or analytically shaped, and whether an additional closure is imposed along the parallel coordinate.

In the Pyrokinetics workflow, the coordinate reconstruction summarized below follows the prototype implementation in the Pyrokinetics example script availabe in the Pyrokinetics repository [3, 30]. Fields are reconstructed in perpendicular real space and then embedded using equilibrium or local-geometry flux-surface shapes. A distinctive feature of this path is the explicit enforcement of the parallel closure Equation 4, which provides a phase-consistent periodic extension in θ determined by the safety factor q and the discrete binormal harmonic index n. In contrast, the GKW gridded-format workflow starts from fields already sampled on a real-space computational grid and embeds them using an analytic shaped-torus model. This model supplies a visualization geometry rather than geometry recovered from the simulation, so changes to the shape should be treated as artificial and communication-oriented rather than as scientific statements about the dataset geometry. This route also does not introduce an additional closure in θ beyond the periodicity inherent to the native gridding and any imposed repetition in the binormal direction.

In cases where spectral coordinates are used, flux-tube fields are typically stored in spectral coordinates as ϕ ( k x , k y , θ ) (time index suppressed). After reading input and output data, we reconstruct the perpendicular real-space field via an inverse Fourier transform

ϕ x , y , θ = ∑ k x , k y ϕ k x , k y , θ e i k x x + k y y . (1)

The binormal coordinate is mapped to a 2 π-periodic field-line label α,and the toroidal angle used for embedding is taken (up to sign convention) aswith q evaluated at the reference surface (or locally). A key step specific to this workflow is enforcing periodic closure in θ consistent with the field-line pitch angle (twist): for discrete binormal harmonics n ≡ k y / Δ k y,implemented numerically by appending a boundary slice at θ = π (or equivalently θ + 2 π) with the phase factor in Equation 4. Geometry is obtained from an equilibrium/local-geometry model providing flux-surface shapeswhere x is the local radial coordinate after the perpendicular real-space reconstruction, and ρ ( x ) denotes the corresponding radial coordinate used by the equilibrium or local-geometry model, such as a normalized minor radius or flux-surface label. The 3D embedding then follows by mapping to Cartesian coordinates using the geometric toroidal angle φ = ± ζ:

α = 2 π L y y , α ∈ − π , π , (2)

ζ θ , α = q θ − α , (3)

ϕ x , y , θ + 2 π = ϕ x , y , θ e − 2 π i q n , (4)

R = R ρ , θ , Z = Z ρ , θ , ρ = ρ x , (5)

X = R ⁡ cos ⁡ φ , Y = R ⁡ sin ⁡ φ , Z = Z . (6)

For the GKW gridded-format route, outputs in this format are stored directly on a 3D grid f ( ψ , s , ζ ) (e.g., potential/density/temperature), with optional linear refinement in each direction prior to embedding. The parallel coordinate is mapped to a poloidal-like angleand a minor-radius-like quantity is built from ψ,

θ = 2 π s , (7)

ρ ψ = ε + ρ ⋆ ψ − L x 2 . (8)

The geometric toroidal angle is then constructed explicitly (including magnetic shear and a binormal phase shift) asand the shaped flux surface is given analytically (Miller-like).

φ ψ , s , ζ = q θ 1 + s ̂ ρ − ε − 2 π ρ ⋆ ζ − L y 2 , (9)

R ρ , θ = R 0 1 + α + ρ ⁡ cos θ + sin − 1 δ sin ⁡ θ , (10)

Z ρ , θ = Z 0 + κ ρ ⁡ sin ⁡ θ . (11)

Finally, Cartesian coordinates are obtained via the same cylindrical map Equation 6 (here with φ from Equation 9). Unlike the Pyrokinetics workflow, this GKW gridded-format pipeline does not introduce an additional closure in θ; periodicity is handled through the native gridding and, optionally, by repeating the ζ domain.

2.2 Workflow

Figure 1 summarizes an end-to-end workflow that starts from raw gyrokinetic simulation outputs and ends with publication-ready visualizations. The pipeline was designed for HPC deployment while still supporting workstation rendering, interactive inspection and deployment of custom solutions (e.g., Omniverse-based viewer). The workflow is structured in stages that separate the coordinate transformation and VTK export from the rendering stage, allowing the same processed dataset to be reused across several rendering environments. This also allows users to apply only the stages needed for a given task, such as skipping the rendering stage when only the VTK output is required for further processing or inspection.

The first stage is job setup and generating raw data, which is dependent on the gyrokinetic code being used. Raw data then enter one of two ingestion paths. The GKW gridded format and compatible formats, including GyroSwin predictions, are handled through a direct Python route, whereas outputs from other sources are processed through Pyrokinetics. In both cases, the goal is to map field-aligned quantities into real-space data suitable for visualization. GKW has two output formats relevant to this workflow. One of them is supported by Pyrokinetics, so GKW data can in principle follow the Pyrokinetics route; this possibility is indicated by the dashed connection from GKW to the Pyrokinetics workflow in Figure 1. The GKW and GyroSwin datasets used here use the other output format, which is not compatible with Pyrokinetics, so they are processed through the direct Python route.

This stage can also perform optional analysis-oriented preprocessing. A recurrent feature of gyrokinetic turbulence is the coexistence of zonal flows, which are large-scale axisymmetric E × B structures with k y = 0, and non-zonal turbulent fluctuations. Zonal flows regulate transport by shearing and suppressing turbulent eddies, and are a central element of drift-wave self-organization [1]. Separating the zonal from the turbulent components through binormal Fourier-space filtering is therefore a physically motivated decomposition that we support in the proposed workflows. The data may also be interpolated to improve visual resolution.

After preprocessing, one VTK dataset is written per time step. This step is naturally parallel and can be distributed across many HPC jobs. The resulting VTK series is then passed to a common processing stage, typically in pvpython, where volume representations, isosurfaces, and color mappings are constructed. External geometry, such as tokamak CAD models, can also be loaded at this stage to provide machine context.

Rendering is chosen according to the scientific or communication objective. ParaView is used for routine scientific analysis, interactive inspection, and reproducible high-fidelity images. Blender [31] is used for communication-oriented outputs, including interactive scene setup, batch rendering of long time series, and multivariate volume visualization with Cycles. Omniverse [32, 33] provides the custom USD-based route used here for interactive local and remote exploration of the processed datasets. Final products are exported as still images, videos, or animations. The rendering stage can be skipped when the goal is to inspect the VTK datasets directly, such as through the Omniverse viewer or by loading them into ParaView for interactive exploration.

The staged structure also makes fidelity and storage trade-offs explicit. Increasing interpolation density and visual detail improves image quality, but also increases VTK file size and total processing cost. The workflow therefore supports low-, mid-, and high-fidelity operation, allowing users to match computational cost to purpose, from rapid exploratory checks to final production rendering.

2.2.1 GKW gridded format

In the direct branch, raw GKW outputs are processed directly with Python scripts on HPC. The scripts ingest post-processed three-dimensional scalar fields on the native field-aligned grid, apply geometric transforms, and perform any required interpolation or filtering. One VTK file is produced per time step, after which the data enter the common downstream stage for VTK processing and rendering in ParaView, Blender, or Omniverse.

One useful feature of the GKW gridded-format workflow is that the visualization-ready output can be shaped by adjusting the Miller-like parameters used in the coordinate transformation. Since this direct route operates on gridded field outputs rather than reading simulation geometry, the geometric presentation is supplied analytically. This kind of shaping should be understood as an artificial visual device, similar to giving a flux tube a finite visible width even though it represents a narrow local domain. It is useful for general-audience communication, where placing the turbulence inside a recognizable tokamak context helps interpretation, but it would not be used for scientific communication about the geometry of the dataset itself.

2.2.2 GyroSwin surrogate model data

GyroSwin predictions can be processed through the same downstream branch as GKW gridded-format data, but they require an additional repadding step before that branch can be used. This follows from the GyroSwin data representation. As described in the GyroSwin paper [11], the electrostatic potential inherited from GKW carries the real-FFT padding used in the original code, so its spatial dimensions do not match those of the associated five-dimensional distribution-function data. For model training and inference, that mismatch is inconvenient because GyroSwin relies on consistent tensor shapes for normalization and for its hierarchical upsampling and downsampling operations. The potential is therefore unpadded during surrogate-data preprocessing.

When these surrogate predictions are returned to the visualization workflow, the padding must be restored so that the predicted potential again matches the spatial layout expected by the GKW-based preprocessing and coordinate-mapping routines. In practice, this is done by repadding the perpendicular planes in Fourier space and then transforming them back to real space before the usual interpolation, geometric mapping, and VTK export are applied. This is the main difference from the GKW gridded-format workflow and allows GyroSwin predictions to be visualized with the same downstream tools and settings as the reference simulations.

2.2.3 Pyrokinetics

The Pyrokinetics branch is the key step that extends the workflow beyond the GKW gridded-format route and opens it to a broader class of gyrokinetic codes. Rather than developing separate preprocessing routines for each code, we use Pyrokinetics as a common interface wherever it can provide the required field data and geometry information. In this work, CGYRO and GENE serve as representative examples, but the same route is intended to apply more generally to any code supported by Pyrokinetics.

In this branch, raw simulation outputs are first processed with Pyrokinetics for spectral-to-real-space reconstruction and geometry-aware embedding before exporting VTK files. As in the direct path, processing is performed per time step and can be parallelized across many HPC jobs. Once the VTK files have been generated, the workflow converges with the direct branch and uses the same downstream tools for processing, rendering, and final media generation.

2.3 Omniverse-based interactive viewer

For interactive inspection of processed datasets, we developed a custom viewer based on NVIDIA Omniverse Kit [32] and Kit-CAE [33]. In this workflow, the viewer sits downstream of the VTK export stage and handles interactive inspection and remote review rather than final offline rendering. The datasets are loaded into a USD stage that stores the data references, volume representations, transfer functions, and rendering settings. It is important to note that production-quality rendering is still possible in Omniverse through the Kit applications, however in this work we leverage Omniverse’s real-time capabilities for interactive exploration and remote review of large datasets through a USD viewer application.

Therefore, the runtime is configured as a USD viewer for CAE data rather than as a general authoring application. Omniverse Kit provides the application framework, OpenUSD provides scene composition, and Kit-CAE provides the VTK delegate and visualization operators used to resolve datasets, fields, and volume representations at run time. Rendering is performed on the server leveraging a RTX GPU and NVIDIA IndeX.

For remote use, the viewer is coupled to a web application that receives a WebRTC video stream from the Kit session and sends control messages back to the viewer. The browser does not render the data. It provides controls for the live Omniverse session. Scene queries, visibility changes, transfer-function edits, renderer changes, and camera operations are applied directly to the USD stage, with the browser updated from the live response returned by the Kit application. Because the interaction controls an active remote session rather than launching a new offline render, the experience is seamless under suitable network conditions. The main practical limitation is the quality of the network connection, since available bandwidth and latency affect the streamed image quality and responsiveness. The Omniverse viewer application and web application both run on an NVIDIA OVX specification 3.0 server, but the user is able to access the web endpoint from any device in the network, including a laptop or tablet, without needing to install Omniverse locally. This makes it possible to inspect large datasets that are too expensive to load and render on a local workstation, and to share the experience with collaborators in different locations.

The web interface provides the scene hierarchy, transfer-function controls, volume-render settings, and camera controls used for routine inspection. The transfer-function editor can request histograms from the live dataset for the selected volume field, so color and opacity edits are made against the field distribution being viewed. Additional controls allow switching renderer backends, changing resolution-related settings, and restoring predefined or startup camera views. This part of the workflow supports rapid inspection of large three-dimensional fields and remote review of processed gyrokinetic data.

This viewer and web interface are a work in progress, and we plan to extend its feature set in future work. For example, we intend to enable time-series exploration, allow users to upload their own datasets and implement more complex data analysis operations within the viewer. Representative screenshots of the current Omniverse workflow are shown in Figure 8.

3 Results

In this section we present example visualizations produced by the workflows described above. The examples are intended to illustrate the capabilities of the workflows and to demonstrate their application across multiple gyrokinetic codes and visualization environments. They are not intended to provide a comprehensive scientific analysis of the underlying data, which is beyond the scope of this work. The datasets explored here are summarized in Table 1.

TABLE 1. Summary of datasets used in visualization examples. The GKW, GKW MAST-U reshaped, and GYROSWIN datasets are available in the Zenodo record [34].

3.1 GKW data

3.1.1 Single variable visualization

Figures 2–4 show single-field renders of the GKW dataset using both surface and volume representations. The surface views in Figures 2, 3 show the field data available from density and temperature fluctuations to the electrostatic potential and its filtered components (zonal and turbulent). While these views are really good at showing the turbulence structure in a poloidal cross-section they lack the perception of how the data evolves in the toroidal direction.

Figures 3, 4 collect alternative views of the turbulent or full potential field. The contour and slice views in Figure 3 expose interior surfaces and reveal how alternating positive and negative filaments (isosurfaces or contours) extend through the volume. Figure 4 emphasizes the fine grained detail of the turbulence filaments. Together these panels demonstrate that the same processed data can be rendered as opaque surfaces, sliced contours, or semi-transparent volumes depending on what information we want to convey.

Figure 6 focuses on a communication-oriented remapping of the same GKW dataset into a MAST-U context. The GKW field data listed in Table 1 is reshaped analytically with Miller-like parameters chosen to be representative of a MAST-U-like geometry. This additional step is needed because the direct route operates on the gridded GKW output format and does not recover the simulation geometry, whereas the Pyrokinetics route obtains geometry information from the simulation data itself. This example highlights how the GKW workflow can preserve the underlying turbulence content while changing the geometric presentation for outreach and general-audience communication, rather than for scientific communication about the dataset geometry. These figures use some of the same concepts of the scientific visualizations (e.g., segmenting filaments of highest and lowest electrostatic potential) but they present the data with necessary context for audiences less familiar with the science behind the data.

3.1.2 Multivariate and communication-oriented visualization

Figure 5 shows the original GKW case rendered in Blender inside an idealized tokamak CAD context constructed to be consistent with the circular geometry of the original GKW case and adds a higher-resolution detail view of the outboard mid-plane region of the same dataset. These views demonstrate the use of Blender Cycles renderer for multivariate volume rendering. These images show the turbulent electrostatic potential as a semi-transparent volume with a custom shader that maps color to the potential value, opacity to the density fluctuation, and emission to the temperature fluctuation. The result is a single volume that simultaneously encodes three physically meaningful fields, making it possible to visually identify relationships among them.

Figure 6 displays ParaView and Blender renders of the reshaped GKW dataset. The difference between these and the scientific ParaView examples is the use of emissivity in the volume shader to communicate to a general audience that the plasma is a source of energy. This highlights how we can leverage simple features of the rendering environment to produce visual stories that are more effective for communication and outreach, while still being traceable to the same underlying data and scientific content.

Figure 7 presents a full-machine render derived from the same dataset and submitted to the 2025 DiRAC Research Image Competition [36]. This visualization combines elements of storytelling and scientific communication we’ve shown in previous examples, the focus here is on detail and visual impact.

3.1.3 Omniverse-based interactive viewer

Figure 8 documents the browser-based and Omniverse-native interfaces used for interactive inspection. The browser panels show the progression from session entry to operator controls and spectator view, and they make visible the controls available for scene selection, camera placement, colormap editing, and volume transfer-function adjustment. The Kit editor view shows the same content inside the authoring environment, where processed geometry, lighting, and CAE assets can be configured in a shared USD scene.

3.2 GyroSwin surrogate model predictions

Figures 2–4 present GyroSwin prediction outputs using the same toroidal mapping and complementary view styles we’ve shown previously for GKW data. Figures 2, 3 collect the surface renders, showing the turbulence detail in a poloidal cross-section and in a toroidal context. Figure 4 finally shows a detailed volume render of the turbulence structures.

3.3 Other gyrokinetic codes

The upcoming subsections show example visualizations of GENE and CGYRO data. These datasets were not chosen because they are particularly interesting in their own right, or because they are directly comparable to the GKW case, but rather because they illustrate the generality of the workflow and its applicability to multiple gyrokinetic codes. The Pyrokinetics branch is designed to be a common interface for any code that can provide the required field data and geometry information, so the GENE and CGYRO examples serve as demonstrations of this generality.

3.3.1 GENE

Figures 3, 4 show that the same visualization pathway can also be applied to GENE outputs. The GENE dataset summarized in Table 1 is drawn from the nonlinear electromagnetic STEP flux-tube case reported by Kennedy et al. [35]. Figure 3 presents the full electrostatic potential field, zonal component, and turbulent component slices side by side on the same mapped geometry. Figure 4 complements this with a cropped turbulent-potential volume rendering in which the more prominent fluctuations were segmented to make their 3D structure more visible.

3.3.2 CGYRO

Figure 4 showcases a detailed crop of a volume rendering of CGYRO outputs. This view was achieved in a similar fashion to the GENE case and the ones visible in the GKW panels of Figure 4.

4 Discussion

4.1 Addressing the visualization gap in gyrokinetics

Advanced three-dimensional visualization of gyrokinetic simulation data has not become routine mainly because the required processing is complex and reproducible tools supporting multiple gyrokinetic codes have been limited. The workflows presented here address this by separating the coordinate transformation from the rendering stage. Once the data have been written as a VTK intermediate product, the same processed dataset can be reused across several rendering environments. This also allows users to apply only the stages needed for a given task.

The staged structure supports a range of fidelity levels. Increasing interpolation density and visual detail improves image quality but increases VTK file size and processing cost; low-fidelity settings are therefore well suited to rapid exploratory work, while full-resolution passes are reserved for final production renders. This makes it practical to incorporate advanced visualization into routine analysis without incurring prohibitive cost at every iteration.

The workflow also allows external tokamak CAD geometry to be loaded at the VTK processing stage, so that fluctuation fields can be rendered within the surrounding machine geometry. This differs from the near-outboard-midplane plots often used in gyrokinetics analysis and makes the data easier to interpret for readers who are not familiar with the coordinate system and physical meaning of the data.

4.2 A common data product for multiple visualization environments

The VTK intermediate format allows the same processed dataset to be used in multiple rendering environments without repeating the upstream processing. In this work, ParaView, Blender, and NVIDIA Omniverse serve different purposes, but all still ingest and operate on the same VTK data.

ParaView is best suited to day-to-day scientific analysis by expert users and to the production of accurate, high-fidelity visualizations for publications and presentations. Its Python scripting interface (pvpython) makes it straightforward to build reproducible rendering pipelines, and its native VTK support means that processed datasets can be loaded with minimal additional configuration. The single-field results for the GKW data were produced primarily in ParaView, and it also provided a natural environment for comparing different filter settings, such as visualizing the full potential field alongside the isolated zonal ( k y = 0 ) and turbulent ( k y ≠ 0 ) components.

Blender is not designed specifically for scientific visualization, but its physically based Cycles renderer is well suited to producing images for dissemination and broader audiences. ParaView also provides established volume-rendering capabilities; the advantage of Blender in this workflow is its flexible material system, which makes it straightforward to build custom shaders that combine emissivity, opacity, and color mapping to encode multiple field variables for communication-oriented multivariate renders.

NVIDIA Omniverse provides real-time interactive visualization within a USD-based environment that can be extended for more complex fusion-visualization tasks. Although it is not a dedicated scientific-visualization tool, its CAE data support and scene management make it useful for local and remote visualization. When paired with specialized hardware (such as the NVIDIA OVX system), it can handle datasets that are too large to inspect frame by frame in a traditional post-processing workflow.

4.3 Multivariate visualization of gyrokinetic fields

Multivariate visualization seeks to represent multiple related variables simultaneously within a single rendering. The introduction noted that there is no universal approach to multivariate datasets and that the optimal strategy depends strongly on the data. In the context of gyrokinetics, the physically meaningful variables are the fluctuating electrostatic potential ϕ, density δ n, and temperature δ T, all of which evolve together and are physically coupled through the underlying gyrokinetic equations. Representing them in a single view creates the possibility of directly examining their spatial correlations, something that is not available from separate scalar plots.

The results also show that Blender’s Cycles renderer can produce multivariate volume visualizations by assigning separate colors, emissivity properties, and opacity mappings to different field variables. This makes it possible to compare the spatial structure of potential, density, and temperature fluctuations in a single view and to identify regions where the fields reinforce or differ from one another.

4.4 Visualization of surrogate model predictions

The GyroSwin pathway illustrates an extension of the workflow to surrogate model outputs. Surrogate models for gyrokinetics are an active area of development, and it is useful to be able to inspect their predictions in the same terms as the simulations they are trained on. As described in the methodology, the repadding step restores the spatial layout expected by the coordinate-mapping routines, allowing GyroSwin predictions to pass through the same downstream processing as the reference GKW data.

This makes it straightforward to compare surrogate predictions and reference simulations side by side within the same visualization environment, using the same color scales, camera settings, and geometry. Such visualizations provide a qualitative check on the spatial structure of the predictions and complement the quantitative error metrics usually reported for surrogate models. They are particularly useful when assessing whether physically important features, such as the distinction between zonal and turbulent components, are reproduced.

4.5 Scalability and computational trade-offs

The workflows were developed around the needs of UKAEA and the HPC systems used most frequently at the institution, particularly the Cambridge Service for Data-Driven Discovery (CSD3), the local UKAEA HPC cluster and the local visualization cluster (NVIDIA OVX system). The per-time-step structure of the VTK export stage means that this step is naturally parallel and can be distributed across many independent HPC jobs, which is important given the size of typical gyrokinetics datasets and the number of time steps in a converged run.

The main computational bottleneck is the coordinate transformation combined with the increase in grid resolution required for on-screen visualization. This step is both time-consuming and storage-intensive: the visualization-ready VTK dataset can be orders of magnitude larger than the original simulation output. The current design treats this transformation as a one-time preprocessing cost and prioritizes reproducibility and speed of development. A more scalable solution, such as deferred or progressive loading that avoids materializing the full-resolution dataset at once, would be preferable for routine use and remains future work. The configurable fidelity levels described above partially mitigate this cost in exploratory contexts, but do not eliminate it for high-resolution final outputs.

For simulations approaching billions of elements, rendering every particle, cell, or interpolated sample may provide little additional perceptual benefit once the projected data density exceeds the display pixel density. A native level-of-detail (LOD) strategy is therefore a natural extension of the present workflow. In this setting, given the suitability of volume rendering, it would be particularly useful to have a multiresolution or hierarchical voxelization for adaptive volume visualization dependent on the camera distance and projected data density. This would, in theory, reduce rendering cost and improve remote interactivity while preserving the large-scale turbulent structures needed for scientific exploration. Such an approach is not yet implemented or benchmarked in the current workflows.

4.6 Performance analysis

To quantify the cost of the staged workflow, we instrumented a representative GKW run using the same stage boundaries shown in Figure 1. The benchmark used the sampled time steps 3, 201, 399, 597, and 798, and reports three processing stages: Stage 1, raw GKW output to VTKHDF; Stage 2, VTKHDF to OpenVDB; and Stage 3, rendering to PNG. Stage 3 was measured with both the Blender/Cycles branch, which renders from VDB volumes, and the ParaView branch, which renders directly from VTKHDF. Each benchmark job recorded wall time, CPU usage, peak resident memory, artifact sizes, temporary storage growth, SLURM allocation data, node specifications, and success state. GPU rendering was additionally sampled with nvidia-smi to record GPU utilization and VRAM use. The hardware used is summarized in Table 2.

TABLE 2. Hardware and network contexts used for the performance measurements. The server-side workflow stages were run on a local UKAEA HPC cluster; the Omniverse workflows were hosted on a local UKAEA OVX system and remote-viewing values are from a separate browser-based point test.

In practical terms, each benchmark case was run as a small batch workflow. For a selected simulation time step, the raw data were first converted to VTKHDF, the resulting VTKHDF file was then converted to VDB for the Blender path, and finally the data were rendered to PNG using both Blender and ParaView. This sequence was repeated for the five sampled time steps and for the low, mid, and high settings listed in Table 3. The bars in Figure 9 therefore show the average cost of one frame-level task, while the whiskers show the spread across the five sampled time steps.

TABLE 3. Processing stages and fidelity settings used in the performance sweep that map to the workflow in Figure 1.

The benchmark measurements are per-timestep measurements: each data point corresponds to loading and processing one sampled time step, not to processing a full simulation time series in serial. This is important for interpreting the wall times. Since the transformation and rendering for each time step are independent, the same workflow can be scaled across a larger HPC allocation by distributing many time steps across independent jobs. The per-stage values therefore characterize the cost and resource requirement of each frame-level task, while total campaign wall time depends on the available compute resources.

Figure 9 and Table 4 summarize how the main server-side cost changes with the low, mid, and high setting for each processing stage defined in Table 3. The figure gives the visual comparison, while the table gives the corresponding elapsed time, resident memory, and disk-usage values. Within each figure panel, grouped bars compare Stage 1, Stage 2, Stage 3 rendering with Blender, and Stage 3 rendering with ParaView. Each bar gives the mean over the five sampled time steps, with the whisker spanning the minimum and maximum values.

TABLE 4. Values are computed over sampled time steps 3, 201, 399, 597, and 798. Elapsed time is reported as mean with minimum-maximum range; maximum RSS is the largest resident-memory value observed across the five samples; mean disk usage is the mean output artifact size produced by that stage for one time step. Same data can be seen in Figure 9.

TABLE 5. Representative Omniverse remote rendering measurements. FPS and GPU use are point observations from the captured session; the ping result is a network round-trip time measurement, not a full application-level interaction-latency measurement.

The resource allocation differs between the processing and rendering stages. Stage 1 and Stage 2 were CPU-only jobs and requested 64 CPUs; at the highest setting, their mean allocated CPU use was 1.90% and 44.60%, respectively. The rendering jobs ran on visualization nodes with NVIDIA A40 GPUs. Stage 3 with Blender/Cycles requested 32 CPUs and five A40 GPUs, while Stage 3 with ParaView requested 32 CPUs and one A40 GPU. At the highest setting, the sampled GPU utilization reached 100% for both renderers, with peak VRAM growth of 30,914 MiB for Blender/Cycles and 5930 MiB for ParaView.

Table 4 reports the corresponding numerical values, including the mean disk usage for the output artifact produced by each stage. The fidelity sweep confirms that the controls change both runtime and memory. For Stage 1, increasing the interpolation parameter from 5 to 50 raised the VTKHDF output from 1.58 GiB to 13.23 GiB per sampled timestep and increased conversion memory usage from 4.24 GiB to 34.72 GiB. For Stage 2, increasing the VDB sampling from 250 to 1500 raised the mean VDB output from 62 MiB to 11.41 GiB. For Stage 3, increasing the output image size from 1280 × 720 to 3840 × 2160 raised the Blender/Cycles render time from 932 s to 3394 s per timestep and the ParaView render time from 35 s to 145 s per timestep. These results show that the lower-fidelity modes are useful for rapid inspection, while the highest settings are production-rendering settings.

The Omniverse branch was characterized separately because it is a live remote visualization mode rather than an offline image-production stage. In the representative browser-based WebRTC session, the application streamed a 1080p viewport using the RTX - Real-Time 2.0 renderer, with the volume resolution parameter set to 1500. The client accessed the streaming application on Google Chrome web-browser over VPN from a different network, so the network results reflect a realistic remote-access configuration rather than a same-machine or local-desktop test. The current Omniverse viewer does not yet support time-series exploration, so these measurements describe loading and interacting with the data for a single timestep. Client-side performance analytics were measured using Google Chrome’s DevTools, while server-side analytics were measured with the ping and nvidia-smi commands.

These Omniverse measurements (see Table 5) reveal that the streaming application is able to serve the client with 1080p video and keep interactivity real-time. The observed 60.0 FPS point measurement at 1080p is above the 5–10 FPS range often sufficient for qualitative scientific inspection, while the VPN remote-access path had 0.0% packet loss and a mean round-trip time of 27.6 m. The result should nevertheless be interpreted as a representative point characterization rather than a full benchmark distribution: longer instrumented sessions would be required to report robust mean, percentile, and minimum FPS values during extended camera motion. Moreover, the reported values are expected to degrade on worse network conditions on the client-side.

Several limitations remain. The present workflows operate on three-dimensional field data and velocity-space moments, not on the full five-dimensional gyrokinetic distribution function. They are also primarily post-processing workflows rather than fully in-situ visualization systems. The Omniverse viewer currently supports remote inspection of an individual loaded timestep but not time-series playback. Finally, users still need familiarity with gyrokinetic parameters to obtain meaningful outputs; future ParaView and Blender plugins could expose the most important controls through simpler interfaces.

4.7 Conclusion

This paper presents reproducible, scalable, modular and general workflows for advanced visualization of gyrokinetic simulation data and addresses a practical gap between simulation output and interpretable three-dimensional views in real space. The proposed preprocessing pipeline supports multiple gyrokinetic codes through either a GKW-compatible route or a Pyrokinetics-based route, enabling magnetic-field-aligned data to be transformed into a common visualization-ready representation - VTK dataset.

By writing the transformed data to a common VTK representation, the workflows decouple data preparation from rendering and make the same processed datasets reusable across ParaView, Blender, and NVIDIA Omniverse. The examples shown here demonstrate that this common pathway can support single-field views, multivariate renderings of coupled gyrokinetic quantities, and qualitative inspection of surrogate-model predictions alongside reference simulations. Taken together, these workflows provide a practical and extensible foundation for routine visualization of gyrokinetic field data, while also establishing a base for future work on more interactive, scalable, and higher-dimensional visualization capabilities.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2026.1924612/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Didier Mazon, CEA Cadarache, France

Reviewed by: Cristobal Navarro, Austral University of Chile, Chile

Nuno Cruz, University of Lisbon, Portugal