Visualization with CUDA-OpenGL Interoperability

Visualization with CUDA-OpenGL Interoperability#

Real-time rendering of simulations using CUDA-OpenGL interoperability. Field data streams directly from GPU memory into the display without intermediate transfers to CPU memory.

Overview#

soliton_solver provides GPU-accelerated visualization using CUDA-OpenGL interoperability:

  • Zero-copy rendering — Field data flows directly from CUDA memory to OpenGL buffers

  • Real-time performance — Interactive exploration at 30–60 FPS on modern GPUs

  • GLFW windowing — Cross-platform window management and input handling

  • Theory-specific visualization — Each theory defines custom colormaps and observables

The visualization backend is implemented in soliton_solver/visualization/gl_backend.py. Each theory provides a theory-specific renderer in render_gl.py.

CUDA-OpenGL Interoperability#

CUDA-OpenGL interop allows CUDA kernels to read/write OpenGL buffer objects directly without CPU-GPU transfers.

CUDA-GL resource mapping#

The workflow is:

  1. Create OpenGL buffer — Pixel buffer object (PBO) allocated via OpenGL

  2. Register with CUDA — PBO is registered as a CUDA graphics resource

  3. Map resource — CUDA gets a device pointer to the PBO memory

  4. Compute on GPU — CUDA kernels write directly to the PBO

  5. Unmap resource — Resource is released for OpenGL to use

  6. Render texture — OpenGL displays the PBO as a texture

Implementation (from gl_backend.py):

class GLBackend:
    def __init__(self, width, height, title="CUDA-OpenGL"):
        # Create PBO
        self.pbo = GL.glGenBuffers(1)
        GL.glBindBuffer(GL.GL_PIXEL_UNPACK_BUFFER, self.pbo)
        GL.glBufferData(
            GL.GL_PIXEL_UNPACK_BUFFER, 
            width * height * 4,  # 4 bytes per RGBA pixel
            None, 
            GL.GL_STREAM_DRAW
        )
        
        # Register PBO with CUDA
        self.resource = _cuda_register_gl_buffer(
            self.pbo, 
            flags=_cuda_gl_write_discard_flag()
        )
    
    def map_pbo(self):
        """Get a CUDA device pointer to the PBO."""
        ptr, nbytes = _cuda_map_resource(self.resource)
        return cuda_array_from_ptr(ptr, (self.height, self.width, 4), dtype=np.uint8)
    
    def unmap_pbo(self):
        """Release the PBO for OpenGL."""
        _cuda_unmap_resource(self.resource)

Zero-copy transfers#

Traditional flow (CPU-GPU copies):

GPU Field Buffer → CPU Memory → OpenGL Texture → Display

CUDA-OpenGL interop (zero-copy):

GPU Field Buffer ← CUDA kernel → PBO ← OpenGL Texture → Display

Benefits:

  • Bandwidth saved — No PCIe transfers of rendered frame data (~50 MB/frame at 1024×1024)

  • Latency reduced — Eliminates GPU→CPU→GPU round-trip

  • Throughput improved — 30–60 FPS typical vs. 5–10 FPS with copies

Rendering pipeline#

The rendering pipeline consists of:

  1. CUDA density computation — Field observables evaluated on GPU

  2. CUDA colormap — Density mapped to RGBA via CUDA kernel

  3. PBO mapping — Pixel buffer mapped from OpenGL to CUDA

  4. CUDA pixel write — Density + colormap written to PBO

  5. PBO unmapping — Resource released to OpenGL

  6. Texture upload — PBO copied to GPU texture (fast, on-device)

  7. OpenGL rendering — Texture displayed as fullscreen quad

Density kernels#

Each visualization mode computes a density field. Examples:

Energy density — Local energy at each grid point:

@cuda.jit
def compute_energy_kernel(en, Field, d1fd1x, p_i, p_f):
    """Compute energy density."""
    x, y = cuda.grid(2)
    if not in_bounds(x, y, p_i):
        return
    en[idx_field(0, x, y, p_i)] = (
        # Theory-specific energy formula
        ...
    )

Skyrmion number — Topological charge density:

@cuda.jit
def compute_skyrmion_number_kernel(density, Field, d1fd1x, p_i, p_f):
    """Compute skyrmion number density."""
    # Integrand of winding number formula
    ...

Field magnitude — Norm of complex/vector fields:

@cuda.jit
def compute_norm_kernel(density, Field, p_i):
    """Compute field magnitude."""
    x, y = cuda.grid(2)
    if not in_bounds(x, y, p_i):
        return
    density[idx_field(0, x, y, p_i)] = math.sqrt(
        Field[idx_field(0, x, y, p_i)]**2 + 
        Field[idx_field(1, x, y, p_i)]**2
    )

Color mapping#

Density values are mapped to RGBA colors. Several colormaps are available:

Jet colormap#

Standard blue-cyan-green-yellow-red colormap. Used for scalar densities.

@cuda.jit
def render_jet_density_to_rgba(pbo_rgba, density, min_val, max_val, p_i):
    """Map density to jet colormap."""
    x, y = cuda.grid(2)
    if not in_bounds(x, y, p_i):
        return
    
    # Normalize density to [0, 1]
    val = (density[idx] - min_val) / (max_val - min_val)
    val = max(0.0, min(1.0, val))  # Clamp to [0, 1]
    
    # Jet colormap interpolation
    if val < 0.125:
        r, g, b = 0.0, 0.0, 0.5 + 0.5 * (val / 0.125)
    elif val < 0.375:
        r, g, b = 0.0, 0.5 * ((val - 0.125) / 0.25), 1.0
    elif val < 0.625:
        r, g, b = 0.5 * ((val - 0.375) / 0.25), 1.0, 1.0 - 0.5 * ((val - 0.375) / 0.25)
    elif val < 0.875:
        r, g, b = 1.0, 1.0 - 0.5 * ((val - 0.625) / 0.25), 0.0
    else:
        r, g, b = 1.0, 0.5 * (1.0 - (val - 0.875) / 0.125), 0.0
    
    pbo_rgba[y, x, 0] = int(r * 255)
    pbo_rgba[y, x, 1] = int(g * 255)
    pbo_rgba[y, x, 2] = int(b * 255)
    pbo_rgba[y, x, 3] = 255  # Alpha

HSV colormap#

Hue-Saturation-Value colormap. Useful for vector fields (magnetization, etc.).

def magnetization_to_hsv(mx, my, mz):
    """Convert magnetization vector to HSV."""
    # Hue from in-plane angle: atan2(mx, my)
    hue = (0.5 + (1.0 / (2 * M_PI)) * atan2(mx, my)) * 360.0
    
    # Saturation from out-of-plane component
    saturation = 0.5 - 0.5 * tanh(3 * (mz - 0.5))
    
    # Value proportional to mz
    value = mz + 1.0
    
    return hue, saturation, value

Advantages:

  • Intuitive — Hue encodes direction, saturation/value encode magnitude

  • Perceptually uniform — Color changes reflect data changes smoothly

  • Vector-friendly — Naturally represents 2D in-plane + 1D perpendicular components

Interactive controls#

The viewer supports real-time interaction via keyboard and mouse.

Keyboard controls#

Typical control scheme (theory-dependent):

Key

Action

n

Toggle arrested Newton flow on/off

k

Toggle velocity arrest on energy increase

F1

Cycle display mode (energy, order parameter, etc.)

v

Enter vortex placement mode

a

Cycle vortex ansatz (Bloch, Neel, etc.)

1–9

Set vortex number (winding number)

o

Save results to disk

Esc

Exit viewer

Mouse interaction#

Vortex placement mode — Click to place vortices at cursor position:

def on_mouse_click(x, y, vortex_type):
    """Place a vortex at (x, y)."""
    # Convert pixel coords to lattice indices
    pxi = int(x * params.xlen / window_width)
    pxj = int(y * params.ylen / window_height)
    
    # Kernel writes vortex configuration
    create_vortex_kernel[grid2d, block2d](
        sim.Field, sim.grid, pxi, pxj, vortex_type, ...
    )

Performance optimization#

Frame rate management#

Balance solver accuracy vs. interactive responsiveness:

run_viewer(sim, params, steps_per_frame=5)

  • steps_per_frame=1 — Rendering dominates; interactive but solver takes few steps

  • steps_per_frame=5 — Good balance; ~10 FPS solver + ~6 fps display

  • steps_per_frame=20 — Solver dominates; slower rendering but faster convergence

Display mode selection#

Different observables have different computational cost:

  • Energy density — Relatively cheap; fast gradient computation

  • Topological charge — Moderately expensive; requires derivatives

  • Magnetic field — Expensive; requires curl computation

Choose the cheapest mode for interactive work; switch to expensive modes for analysis.

Texture format#

RGBA8 (8 bits per channel) vs. RGBA16F (16-bit float):

  • RGBA8 — Default, uses less bandwidth, sufficient for visualization

  • RGBA16F — Higher precision, for detailed feature inspection

Custom visualization#

To add visualization to a new theory, implement render_gl.py:

# soliton_solver/theories/my_theory/render_gl.py

from soliton_solver.visualization.gl_backend import GLBackend
import numpy as np
from numba import cuda

class GLRenderer:
    """Custom renderer for my_theory."""
    
    def __init__(self, width, height):
        self.backend = GLBackend(width, height, title="My Theory")
        self.display_mode = 1
    
    def display_frame(self, sim):
        """Render one frame."""
        # 1. Compute observable (density field)
        grid2d, block2d = launch_2d(sim.p_i_h, threads=(8, 8))
        sim.kernels.compute_my_observable_kernel[grid2d, block2d](
            sim.en, sim.Field, sim.p_i_d, sim.p_f_d
        )
        cuda.synchronize()
        
        # 2. Map PBO for CUDA writing
        pbo_rgba = self.backend.map_pbo()
        
        # 3. Color-map and write to PBO
        my_colormap_kernel[grid2d, block2d](
            pbo_rgba, sim.en, min_val, max_val, sim.p_i_d
        )
        cuda.synchronize()
        
        # 4. Unmap PBO
        self.backend.unmap_pbo()
        
        # 5. Display
        self.backend.upload_and_draw()

def run_viewer(sim, params, *, steps_per_frame=5):
    """Launch interactive viewer."""
    renderer = GLRenderer(params.xlen, params.ylen)
    
    while not renderer.backend.should_close():
        renderer.display_frame(sim)
        
        # Advance solver
        energy, err = sim.step(prev_energy=energy)

Key points:

  • Implement compute_*_kernel for your observables

  • Implement a colormap kernel to map density → RGBA

  • Use GLBackend for the window and PBO management

  • Call map_pbo() / unmap_pbo() around CUDA work

Performance characteristics#

Typical performance on NVIDIA A100:

Grid size

Frame rate

Pixel throughput

512×512

60 FPS

~150 M pixels/sec

1024×1024

30 FPS

~300 M pixels/sec

2048×2048

15 FPS

~600 M pixels/sec

Bottleneck: GPU memory bandwidth and synchronization overhead between CUDA and OpenGL.

Troubleshooting#

Black screen / no display:

  • Check that density kernel writes valid data (not NaN, inf, or all zeros)

  • Verify colormap kernel writes to correct PBO array

Flickering / torn frames:

  • Enable vsync: glfw.swap_interval(1) (default)

  • Increase steps_per_frame if solver is too slow

CUDA-GL errors:

  • Ensure both CUDA and OpenGL operations are synchronized

  • Call cuda.synchronize() before unmapping resources

  • Check that PBO size matches window resolution

Contents