Practical in situ Visualization on Vector Computers via Vectorized Ray-Casting Method

Abstract

We have been developing VISMO, an in situ visualization library originally developed for plasma simulations. Since this library is designed for scalar computers, it performs poorly on vector computers because such architectures achieve high performance when algorithms provide long loops for vectorization. To address this limitation, we developed a new version of VISMO incorporating a ray-casting-based visualization code optimized for vector computers for scalar fields. Performance evaluation using coupled CFD and Hall MHD simulation codes demonstrates that the new implementation enables practical in situ visualization on vector computers. Furthermore, its vectorization-friendly structure design naturally allows extension of loop length, suggesting that the new implementation can serve as a practical basis for future in situ visualization tools on GPGPU-based systems.

1.  Introduction

Along with theoretical and experimental approaches, numerical simulations play a vital role in elucidating the physics of plasmas. Large-scale plasma simulations now routinely generate massive data with multiple physical variables. Due to the size of modern simulation data, traditional post-processing visualization is becoming increasingly impractical. This is because the approach involves transferring the simulation output to a local personal computer (PC) for visualization, which becomes time-consuming due to the large data size. Furthermore, the CPU performance and memory capacity of the PC are typically insufficient to handle such massive data. Consequently, researchers often must either reduce the data size or restrict visualization to specific regions of interest.

To address this challenge, in situ visualization, in which visualization is performed on the supercomputer simultaneously with the simulation, has been actively studied [1, 2]. By tightly coupling the simulation and visualization codes, this approach eliminates the need for post hoc data transfer and circumvents the limited CPU performance and memory capacity of local PCs. Although in situ visualization suffers from limited interactivity due to batch processing, recent studies aim to alleviate this issue [3, 4].

Currently, many high-performance computing systems rely on specialized accelerators. To exploit these architectures, simulation codes must be optimized for the target architecture. Vector computers are widely used in computational fluid dynamics (CFD) including plasma simulations because these systems enable high-speed execution of CFD codes. Modern vector computers, such as the SX-Aurora TSUBASA vector engine (VE) developed by NEC Corporation, are provided as accelerators. However, visualization codes, typically designed for scalar architectures, often perform poorly on such systems.

Several studies have investigated vectorization of computer graphics and visualization algorithms. Among them, a pioneering study proposed a vectorized ray-tracing method for vector computers [5]. For volume rendering [6], Spalt vectorized the splatting method for volume rendering [7], while Murai et al. developed a visualization library called RVSLIB for the NEC SX-4 [8]. Similarly, Chen et al. vectorized the ray-casting method for unstructured grid visualization on the first Earth Simulator [9].

We have developed VISMO, an in situ visualization library supporting visualization methods such as isosurface, volume rendering, slicing, and streamline visualization [10, 11]. This library was originally developed for particle-in-cell simulations, and it has since been used for other plasma simulations, such as magnetohydrodynamic (MHD) simulations [12], non-plasma fluid simulations, as well as research related to in situ techniques [13]. Since VISMO was designed for scalar computers, visualization methods, including isosurface which is a key visualization method for MHD simulations, do not perform well on vector computers. Therefore, to enable efficient operation on vector architectures, VISMO must be adapted for vector computers.

In this study, we developed a visualization code for scalar fields based on the ray-casting method optimized for vector computers such as SX-Aurora TSUBASA. We integrated it into VISMO and evaluate its performance through coupled simulation codes. The results show that the developed code, VISMO for vector computers, enables practical in situ visualization on vector computers.

This paper is organized as follows. Sections 2 and 3 describe the design of VISMO for vector computers and its rendering characteristics, respectively, and Sec. 4 presents its performance. Discussion, Summary and Conclusion are provided in Secs. 5 and 6, respectively.

2.  Development of VISMO for Vector Computers

We developed VISMO for vector computers to visualize scalar fields defined on a grid in Cartesian coordinates, with a primary focus on the isosurface method. This section details the design and implementation.

2.1  Ray-casting-based visualization of scalar fields

2.1.1  Overview

The VISMO library employs the ray-casting method for scalar-field visualization. The ray-casting method [14] is known to be a technique for implementing volume rendering but it can also be used for isosurface and slicing visualization. Figure 1 illustrates the ray-casting method. The ray-casting method has the property that each ray can be processed independently like the ray-tracing method [5, 8], making it suitable for restructuring into long vector loops. For isosurface rendering, the scalar values at successive sampling points along each ray are compared to determine whether the ray intersects the specified isosurface level—that is, whether it passes through the isosurface. For slicing, intersections with a slicing plane are identified by comparing the positions of successive sampling points. In the case of opaque slices, the intersection point can be directly calculated using the equations of the ray and the slicing plane.

Fig. 1.  Overview of the ray-casting method. As the viewing line passes through the data from the viewpoint through the pixel, the color of the pixel is determined by integrating the colors at the sampling points in the data, taking opacity into account. The color and opacity are obtained by transforming scalar values with the transfer function, which is the relationship between scalar values and color and opacity.

2.1.2  Vectorization strategy

Vector processors achieve high performance when executing long loops. Therefore, our vectorization strategy analyzes the ray-casting algorithm and restructures it to maximize loop length.

Our approach is based on the concept of constructing long vector loops by collecting rays, which was proposed for ray tracing in [5]. In that study, primary rays cast through image pixels and secondary rays spawned by operations such as reflection and refraction are gathered to form a long vector loop of a predetermined length and are then processed together. Specifically, the computations associated with the collected rays, such as intersection tests with scene primitives, are performed using vector operations. The rays that remain to be processed and the newly spawned rays are collected and processed in subsequent vector-processing cycles. This process is repeated until the colors of all pixels are determined. A similar approach has been explored for ray casting in RVSLIB [8], which vectorizes ray casting for volume rendering by grouping sampling points with the same sampling index across multiple rays. However, the details of its implementation are not available in the literature.

We adapt the ray-tracing concept proposed in [5] to the ray-casting method. Unlike ray tracing, the ray-casting method does not generate secondary rays. Each ray is processed as it advances through the data until it exits the data domain or reaches a rendering termination condition. The lifetime of each ray varies depending on the visualization method (e.g., isosurface or volume rendering) and the characteristics of the data. Based on these considerations, we designed the following processing pipeline.

1. 1.   Ray queue construction

   Rays are first cast from the viewpoint through all pixels on the projection plane. Only rays that actually intersect the data domain are queued for processing (see Fig. 2). This eliminates unnecessary computation and increases the efficiency of vector loops.

2. 2.   Vector processing loop

   Two hundred and fifty-six rays, the vector register length of SX-Aurora, are extracted from the head of the ray queue to form the processing set. In other words, multiple rays passing through the data are grouped together to construct loops (see Figs. 2 and 3). Multiple loops for vector processing are created at each sampling point because various computations are performed at each sampling point, such as comparing scalar values and positions at successive sampling points for isosurface and slicing respectively, integrating colors for volume rendering, interpolating normal vectors and advancing ray tips. Other operations are also integrated into vector loops wherever possible.

3. 3.   Dynamic loop reconstruction

   Rays for which computations have been completed (e.g., exit the data domain or hit an opaque isosurface) are removed from the processing set because keeping them in the loops unnecessarily consumes computational resources. Then, remaining rays are shifted forward to fill the gaps, resulting in a continuous processing set, and new rays from the ray queue are added (see Figs. 3 and 4). This dynamic loop reconstruction ensures efficient execution.

Fig. 2.  Construction of ray queue. Rays that require computation are identified in advance and organized to form a queue of rays to be computed, referred to as the ray queue.

Fig. 3.  Loop for vector processing. The loops are constructed from rays extracted from the processing set. As many operations as possible, including advancing the ray tip, are performed through vectorized loops.

Fig. 4.  Loop reconstruction. Rays for which computations have been completed are dynamically removed from the loop, and unprocessed rays are added to maintain loop length to ensure efficient vectorization.

The processes described in steps 2 and 3 are repeated until the ray queue becomes empty. In addition, we created loops for vectorization in other parts of the code as much as possible.

2.2  Use of VHcall

SX-Aurora TSUBASA provides the VHcall function, enabling subroutines to run on the scalar computers called “vector host” (VH). This offloading model is the reverse of the GPGPU offloading model: data are sent from the VE (i.e., the accelerator) to the VH, where subroutines are executed. By using this feature, tasks that are inefficient on the VE can be executed on the VH, improving the overall performance of the program.

In this development, the VE does not handle the PNG image generation task efficiently. Therefore, we used VHcall to execute this subroutine on the VH.

2.3  Integration into VISMO

We integrated the scalar field visualization code developed for vector computers into VISMO with compatibility as a priority. Common components such as the configuration reader were reused. To ensure seamless adoption, both original and vector versions share the same interface names, and conditional compilation selects the appropriate implementation. Users can enable the vector version simply by adjusting compilation options, without modifying their simulation codes.

3.  Rendering Characteristics

VISMO generates the final visualization image by rendering the subdomain assigned to each MPI process and compositing the resulting subimages. Since the image composition requires synchronization across MPI processes, the overall performance is constrained by the process with the highest image rendering workload. Since VISMO employs the ray-casting method for isosurface rendering, the computational cost increases with the number of pixels onto which the subdomain is projected. This implies that, even when the visualization parameters are the same, a larger image size leads to a proportional increase in computational cost. Increasing the number of grid points in the simulation also leads to a drop in rendering performance. Another characteristic is that the computation along a ray terminates once the ray intersects an isosurface. Consequently, when a large isosurface is located close to the projection plane within a subdomain, the rendering computation can be completed relatively quickly. This also implies that the rendering time may vary depending on the viewpoint and viewing direction, even for the same isosurface. Therefore, rendering performance depends not only on the characteristics of the data but also on the viewpoint and viewing direction. As an acceleration technique, VISMO subdivides the subdomain assigned to each process into up to 512 blocks, and determines whether the target isosurfaces exist in each block in advance. As a result, rays can skip blocks that do not contain any part of the target isosurfaces. Accordingly, one of the most computationally demanding cases for VISMO is one in which isosurfaces are distributed throughout the field of view, with gaps between them that allow rays to travel a long distance before either intersecting an isosurface or leaving the subdomain. In such situations, only a limited number of rays can be terminated at an early stage, resulting in a high rendering cost.

The basic rendering characteristics of VISMO were investigated using small synthetic datasets and are described in the following sections. The investigations were conducted on an HPC system equipped with SX-Aurora TSUBASA 10B VEs at University of Hyogo. Because the libpng library was not available on the VEs of this system, the images were output in BMP format without compression, and VHcall was not used.

3.1  Dataset

In this section, the basic rendering characteristics of the newly developed VISMO were investigated using a program that generates small synthetic datasets with two MPI processes. Each MPI process handles a grid of 240³ points, and an identical synthetic dataset is generated in each subdomain. The datasets consist of 2, 16, 54, or 128 spherical isosurfaces. Figure 5 shows the twelve synthetic datasets used in this investigation. The columns correspond to datasets containing 2, 16, 54, and 128 spherical isosurfaces. For each sphere count, the sphere radius is scaled by factors of 1, 1/2, 1/3, and 1/4, respectively. The rows correspond to sphere radii of $R$, $R/2$, and $R/4$, where $R$ denotes the reference radius. This data-generation program was coupled with VISMO, and the generated datasets were visualized to investigate the rendering characteristics of VISMO. Since the visualizations were performed from a symmetric viewpoint with a symmetric viewing direction (see Fig. 5), representative results are reported in this section.

Fig. 5.  Synthetic datasets used for investigation of rendering characteristics. The rows correspond to sphere radii of $R$, $R/2$, and $R/4$, while the columns correspond to datasets containing 2, 16, 54, and 128 spheres, respectively, where $R$ is the reference radius.

3.2  Vectorization efficiency

To investigate the vectorization efficiency, the codes coupled with the newly developed and original VISMO were executed using the single-sphere dataset with the reference radius. The total execution times of the coupled codes were 2.0 and 23.4 s, respectively. The corresponding average vector length and vector operation ratio were 182.4 and 97.4% for the newly developed VISMO, compared with 9.2 and 10.3% for the original one, demonstrating a substantial performance improvement.

For the newly developed VISMO, the preprocessing time required for surveying the ray-skipping blocks was 0.18 s, and the ray-queue construction time, including other preprocessing tasks, was 0.1 s. The ray-casting required 1.28 s, including 0.32 s of dynamic loop reconstruction. The conversion of the rendered image data to BMP format and file output required 0.02 s. These results indicate that the ray-casting stage dominates the overall visualization time. However, the relative contributions of these components depend on the dataset characteristics, viewpoint position, and other visualization parameters.

3.3  Dependence on image size and grid resolution

To evaluate the effect of image resolution on rendering performance, a single-sphere dataset was visualized at image resolutions of 2,048 × 1,024 and 1,024 × 512. The corresponding times required for visualization were 1.61 and 0.56 s, respectively, indicating that the lower-resolution case required approximately 35% of the visualization time of the higher-resolution case. To investigate the effect of grid resolution, the number of grid points was reduced by a factor of eight (240³ to 120³). In a simplified scenario, reducing the grid resolution by a factor of eight decreases the depth of the volume along each axis by a factor of two. In this case, the times required for visualization decreased to 0.96 s, corresponding to approximately 60% of that of the original dataset.

Although the ray-casting stage is the most computationally demanding part of the visualization process, the reductions in visualization time were not directly proportional to the decreases in image size or grid resolution. Nevertheless, both experiments demonstrated that reducing the image resolution or grid resolution shortened the visualization time.

3.4  Dependence on isosurface distribution

The relationship between isosurface distribution and the time required for visualization was investigated, and the results are summarized in Table 1. For the datasets with the reference radius, the longest visualization time was observed in the case containing 16 spheres. This is likely because the arrangement produced relatively wide gaps along the viewing direction, allowing rays to penetrate deeply into the subdomain before intersecting an isosurface or exiting it. In contrast, the 54 and 128 sphere cases exhibited shorter visualization times, despite containing a larger number of isosurface spheres. Although gaps were still present, a larger fraction of rays intersected an isosurface before reaching the opposite side, allowing them to be terminated after traveling a shorter distance and reducing the overall rendering cost.

Table 1. Visualization times (s) for sphere datasets.

Radius	2 spheres	16 spheres	54 spheres	128 spheres
$R$	1.61	2.06	1.65	1.39
$R/2$	0.88	1.14	2.32	3.66
$R/4$	0.66	1.46	1.71	4.52

For the datasets with reference radii of one-half and one-quarter of the baseline value, the visualization time increased monotonically with the number of spheres. In these cases, a larger number of spheres resulted in more widespread isosurface distributions throughout the domain with relatively wide gaps between the isosurfaces. Consequently, the effectiveness of block-skipping was reduced, while many rays still traversed regions that did not contain isosurfaces, leading to an increase in rendering time. It is also noteworthy that, for a given number of spheres, a smaller reference radius did not necessarily result in a shorter time.

These results indicate that VISMO requires longer visualization times when isosurfaces are distributed throughout the domain, with gaps between them that allow rays to penetrate deep into the subdomain before intersecting an isosurface or exiting it. In contrast, high visualization performance can be achieved when rays terminate at relatively short distances due to localized isosurfaces or effective block skipping.

4.  Performance Evaluation

We evaluated the isosurface-rendering performance of VISMO for vector computers, the most important visualization method for our application, to determine whether it is sufficiently fast for practical in-situ visualization. Two types of simulations were prepared for the performance evaluation: a small-scale CFD simulation and a large-scale Hall MHD simulation.

4.1  Sample visualization image

To demonstrate the visualization images that can be generated with VISMO for vector computers, two sample images are shown in Figs. 6 and 7. Figure 6 was generated on the HPC system with SX-Aurora TSUBASA 10B VEs at University of Hyogo, whereas Fig. 7 was generated on the Plasma Simulator built with SX-Aurora TSUBASA A412-8 at the National Institute of Fusion Science (NIFS). An isosurface of the “Smoke Ring” CFD code [15] is shown in Fig. 6. Two isosurfaces of Hall MHD turbulence [16] are rendered with their levels set to the mean plus four times the standard deviation ($\mu +4\sigma$) for enstrophy and current density.

Fig. 6.  Visualization image of CFD simulation code “Smoke Ring”. An isosurface is visualized, corresponding to the enstrophy.

Fig. 7.  Visualization image of Hall MHD turbulence simulation. Two isosurfaces are visualized, corresponding to the enstrophy (green) and the current density (gray), respectively.

4.2  CFD simulation

To evaluate the performance of the newly developed VISMO in a small-scale but practical CFD application, the publicly available CFD simulation code “Smoke Ring” [15] was modified and coupled with it. In this code, the Navier-Stokes equations are solved in Cartesian coordinates. A second-order central difference scheme is employed for spatial discretization, and a fourth-order Runge-Kutta method is used for time integration. The simulation was conducted using two MPI processes, each handling a grid of 240³ points, resulting in a total computational domain of 480 × 240 × 240. The simulation was executed for 3,000 time steps, and an isosurface image of enstrophy with a resolution of 2,048 × 1,024 was generated every 1,000 steps. This performance evaluation was conducted on the HPC system at University of Hyogo using two VEs. The total simulation time was 881.1 s, while the visualization time required for three images was approximately 2.5 s, partly owing to the relatively simple spatial distribution of the isosurface. The original VISMO required 38.4 s, indicating a speedup of more than 15 times through vectorization.

The visualization time required for a visualization image corresponds to approximately 3 simulation time steps. Therefore, if a visualization overhead comparable to the CFD computation cost is acceptable, it may be possible to generate approximately 1,000 images over 3,000 steps of the simulation. Note that the images were output in BMP format without image compression. If image compression had been used, the time required for visualization would likely have been longer.

4.3  Hall MHD simulation

4.3.1  Performance improvements

Performance was measured using the newly developed and original VISMO, with and without VHcall on the Plasma Simulator at NIFS using 16 VEs (64 MPI processes, 2 OpenMP threads each). Both codes were coupled with the Hall MHD code [12, 17], and the data grid size and image resolution were 1024³ and 1024², respectively. Both codes generated three images: two with different isosurfaces, and one with both isosurfaces rendered simultaneously. We measured the times required for visualization, along with the vectorization ratio and the average vector length of the ray-casting subroutines, using VHcall. Table 2 shows the time required for the visualization. VISMO for vector computers achieved more than a 12-fold speedup, with a vectorization ratio of 85.4% and an average vector length of 196.8. VHcall reduced approximately 2 to 3 seconds by offloading PNG generation.

Table 2. Times required for visualization (s).

VHcall	developed code	original code
true	8.79	130.95
false	10.86	134.32

Both metrics showed significant improvement. It was confirmed that the long vector loops formed by 256 rays contribute to the performance gain, reducing the visualization time to one-twelfth or less. It should be noted that these values may vary depending on the characteristics of the data and the visualization parameters, although they are very difficult to quantify (see Sec. 3).

4.3.2  Rendering cost relative to simulation

In this section, to evaluate relative rendering cost compared to simulation, simulation and visualization times were measured for four different cases, varying the data size, image resolution, number of isosurfaces rendered, number of vector hosts, number of MPI processes and OpenMP threads, as well as whether VHcall was used. The simulation code [12, 17] coupled with VISMO for vector computers is a three-dimensional incompressible Hall MHD turbulence code using the pseudo-spectral method and the Runge-Kutta-Gill scheme. These measurements were also performed on the Plasma Simulator at NIFS. The results are summarized in Table 3.

Table 3. Simulation and visualization performance.

	Case 1	Case 2	Case 3	Case 4
Grid size of data	1,0243	2,0483	2,0483	4,0963
Image resolution	1,0242	1,0242	2,0482	2,0482
Number of vector hosts	1	32	32	256
Number of MPI processes	64	2,048	1,024	8,192
Number of threads per process	1	1	2	2
VHcall	True	True	False	False
Number of images (/step)	3	3	18	33
Simulation time (s/step)	17.51	16.13	24.22	32.47
Visualization time (s/step)	69.35	17.59	167.49	158.33
Time per image (s/image)	22.74	8.67	9.30	4.79

In cases 1 and 2, three isosurface visualization images were generated per step. In case 3, a total of 18 isosurface images were generated per step (each of the three images was rendered from 6 different viewpoints, resulting in a total of 18 images), and in case 4, 33 isosurface images were generated per step (each of the three images was rendered from 11 different viewpoints, resulting in a total of 33 images). In cases 1 and 2, one isosurface is rendered in each image. In cases 3 and 4, two isosurfaces are rendered in each image like the sample image in the previous section (Fig. 7). It is worth noting that in situ visualization commonly involves generating images from various viewpoints, which naturally increases the rendering cost and must be taken into account when estimating overall runtime.

5.  Discussion

The performance evaluation results for the Hall MHD simulation show that the computational cost of visualization is higher than that of a single simulation time step, whereas it is comparable to or lower than the simulation time when limited to a single visualization image. In practice, visualization is not performed at every time step, as changes between successive steps are often negligible. Instead, visualization is typically carried out once every 10 to 100 simulation steps. Therefore, the performance achieved here is sufficiently practical, provided that an excessive number of visualization images are not required. In particular, for large data such as 4,096³, it is not practical to store a large volume of raw simulation data. In such cases, saving a limited number of visualization images with various settings (e.g., isosurface levels, viewpoints) is more convenient for data analysis, even if it consumes computation time to some extent. In contrast, when using the original code, the processing time is more than ten times longer, resulting in a rendering capacity of less than one-tenth that of the developed code, which indicates that it is insufficient for practical in-situ visualization.

The performance evaluation using a small-scale CFD simulation demonstrated that the new implementation can achieve practical visualization performance not only for Hall MHD simulations but also for other simulations. In addition, the investigation of rendering characteristics revealed that the visualization performance of VISMO is influenced by the spatial distribution of the target isosurfaces. The investigation also showed that the time used for ray-queue construction and dynamic loop reconstruction was limited compared with that for the ray-casting computation, indicating that this implementation introduces only a modest additional cost. These results indicate that, although the visualization performance depends on the distribution of the target isosurfaces, VISMO is applicable not only to Hall MHD simulations but also to various CFD and plasma simulations.

The present implementation extracts 256 rays to create a vectorization loop but extending the loop length is straightforward if memory permits. This design is expected to serve as a promising foundation for efficient visualization on GPGPU-based computers, such as the new Plasma Simulator built with AMD Instinct MI300A.

6.  Summary and Conclusion

We developed VISMO for vector computers, incorporating scalar-field visualization methods optimized for such systems, with a primary focus on isosurface visualization. The experiments with synthetic datasets clarified that rendering performance is influenced by the spatial distribution of isosurfaces and the distance traversed by rays within the subdomain. The performance evaluations using both the synthetic dataset program and the CFD simulation code coupled with VISMO for vector computers demonstrated speedups of more than one order of magnitude. In the case of the large-scale Hall MHD simulation, the new implementation achieved a vectorization ratio of 85.4%, an average vector length of 196.8, and a practical rendering speed for SX-Aurora TSUBASA. A single visualization image requires computation time comparable to one simulation step, and performing visualization intermittently imposes only a modest overhead. Overall, our results demonstrate that VISMO for vector computers eliminates a major performance constraint, thereby enabling practical in-situ visualization on SX-Aurora TSUBASA-based vector computers.

The improved performance enables a variety of in situ visualizations, including visualization of full-data rendering and the tracking of local structure formation, contributing to deeper understanding of plasma phenomena. Future work includes optimizing vector field visualization methods for vector computers and adapting VISMO to GPGPU-based platforms.

Acknowledgements

The authors would like to thank Prof. Akira Kageyama at Kobe University for his valuable advice on modifying the Smoke Ring code and coupling it with VISMO. This work was supported by JSPS KAKENHI Grant Number JP24K06893. The numerical simulations were performed on an NEC SX-Aurora TSUBASA A412-8 Plasma Simulator at the National Institute of Fusion Science (NIFS), Japan, with the support and under the auspices of the NIFS Collaboration Research program (F23KISS030). This research was partially supported by the “Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures” (JHPCN jh230004, jh240004, jh250013) and by the NINS program of Promoting Research by Networking among Institutions (Grant Number 01422301).

References

• [1]   K.-L.  Ma, IEEE Comput. Graphics Appl. 29, 14 (2009).
• [2]   H.  Yu et al., IEEE Comput. Graphics Appl. 30, 45 (2010).
• [3]   J.  Ahrens et al., In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, (New Orleans, LA, USA, 2014) 424.
• [4]   A.  Kageyama and  N.  Sakamoto, PeerJ Comput. Sci. 6, e305 (2020).
• [5]   D.J.  Plunkett and  M.J.  Bailey, IEEE Comput. Graphics Appl. 5, 52 (1985).
• [6]   R.A.  Drebin et al., In Proceedings of the 15th annual conference on Computer graphics and interactive techniques, (Vancouver, BC, Canada, 1988) 65.
• [7]   A.  Spalt, Proceedings of EDU and COMPUGRAPHICS (ACM, Portugal, 1993) 154.
• [8]   T.  Takei et al., J. Vis. Soc. Jpn. 20(Suppl 1), 115 (2000) (in Japanese).
• [9]   L.  Chen et al., Parallel Comput., 29, 355 (2003).
• [10]   N.  Ohno and  H.  Ohtani, Plasma Fusion Res. 9, 3401071 (2014).
• [11]   N.  Ohno and  A.  Kageyama, Earth Planets Space, 73, 158 (2021).
• [12]   H.  Miura, Plasma 7, 793 (2024).
• [13]   A.  Kageyama et al., Plasma Fusion Res. 15, 1401065 (2020).
• [14]   M.  Levoy, IEEE Comput. Graphics Appl. 8, 29 (1988).
• [15]  Smoke Ring. https://github.com/akageyama/smoke-ring, (accessed 2026-06-29).
• [16]   H.  Miura and  K.  Araki, Phys. Plasmas, 21, 072313 (2014).
• [17]   H.  Miura and  F.  Hamba, J. Comput. Phys. 448, 110692 (2022).