Reading Time: 13 minutes

La aceleración de GPU en la computación científica a menudo comienza con bibliotecas de alto nivel. Un investigador reemplaza a Numpy con Cupy, usa una FFT acelerada o llama a una rutina de álgebra lineal habilitada para GPU. Esto puede producir una mejora sustancial sin requerir un conocimiento detallado del hardware de GPU.

La programación personalizada del kernel de GPU va un nivel más profundo. En lugar de llamar a una operación predefinida, el desarrollador escribe la función que ejecuta cada subproceso de GPU. Esto proporciona control sobre la indexación de subprocesos, el diseño de datos, las transferencias de memoria, la sincronización, la memoria compartida y el número de operaciones realizadas durante el lanzamiento de cada kernel.

Este control es útil cuando un modelo físico contiene patrones de acceso que no se pueden expresar de manera eficiente a través de operaciones de arreglos ordinarios. Las interacciones de partículas, las plantillas de diferencias finitas, las reglas de colisión, los modelos de fluidos basados en la red, las actualizaciones de materiales y las consultas de geometría son ejemplos comunes.

Un kernel personalizado no es automáticamente más rápido que una biblioteca optimizada. Se vuelve valioso cuando el algoritmo expone suficiente trabajo paralelo y cuando el desarrollador puede reducir el tráfico de memoria, combinar operaciones o adaptar el cálculo a la arquitectura de la GPU.

¿Qué es un kernel de GPU?

Un kernel de GPU es una función ejecutada por muchos subprocesos ligeros. Cada subproceso normalmente maneja una partícula, celda de cuadrícula, elemento, cara, píxel o entrada en una matriz.

El programa host lanza el kernel con una cuadrícula de bloques de subprocesos:

kernel[blocks, threads_per_block](arguments)

Cada subproceso determina qué datos debe procesar desde su bloque e índices de subprocesos. Un índice unidimensional típico es:

index =
    block_index * block_size
    + thread_index

Los hilos se organizan en grupos llamados bloques. Los hilos dentro de un bloque pueden sincronizar e intercambiar datos a través de la memoria compartida. Por lo general, se espera que los bloques separados se ejecuten de forma independiente.

Cuando vale la pena escribir los kernels personalizados

Las bibliotecas de GPU de alto nivel generalmente deberían ser la primera opción. Ya proporcionan multiplicación de matrices optimizada, FFT, reducciones, operaciones escasas, generación de números aleatorios y muchas funciones de elementos.

Un kernel personalizado se vuelve útil cuando:

  • La simulación utiliza un patrón de acceso a la memoria no estándar.
  • Se pueden fusionar varias operaciones pequeñas en una sola pasada de memoria.
  • Los subprocesos vecinos necesitan repetidamente los mismos datos locales.
  • Una regla especializada de partículas, colisión o plantilla domina el tiempo de ejecución.
  • Las matrices intermedias consumen demasiada memoria de dispositivo.
  • La simulación necesita una operación diferenciable personalizada.
  • Una biblioteca existente no puede expresar la lógica de los límites o material de manera eficiente.

Antes de escribir un kernel, perfile el programa. La optimización de una función visualmente compleja que representa solo una pequeña parte del tiempo de ejecución no producirá una aceleración general significativa.

Modelos de programación de GPU

El software científico puede acceder a las GPU a través de varios niveles de abstracción.

Enfoque ejemplos Fuerza principal compensación principal
Descarga basada en directivas OpenMP, OpenACC Aceleración incremental de código C, C++ o Fortran existente Menos control directo sobre los núcleos generados
Programación de kernel orientado a proveedores cuda, cadera Control detallado sobre la ejecución de GPU Mayor coste de implementación y mantenimiento
C++ portátil de rendimiento Sycl, Kokkos, Alpaka Un modelo de fuente para varios backends de hardware La portabilidad no garantiza el mismo rendimiento en todas partes
Arreglos y kernels de GPU de Python Cupy, Numba-cuda, urdimbre de Nvidia Desarrollo rápido con acceso a código GPU compilado Restricciones específicas del marco y comportamiento de compilación

La ENCCS Introducción a los modelos de programación de GPU proporciona una comparación más amplia de enfoques CUDA, HIP, OpenCL, SYCL, Kokkos, OpenMP y relacionados.

cuda

CUDA proporciona acceso directo a la programación de GPU de NVIDIA a través de C++, Fortran y varios enlaces de idioma. Expone cuadrículas, bloques, warps, memoria compartida, flujos, eventos, bibliotecas de dispositivos y herramientas de optimización específicas de hardware.

A menudo se selecciona CUDA cuando el máximo rendimiento específico de NVIDIA y las herramientas maduras son más importantes que la portabilidad del hardware. Los principales costos son la gestión de memoria manual, un modelo de programación de nivel inferior y la dependencia de la plataforma NVIDIA.

Cadera y Rocm

HIP proporciona un modelo de programación C++ similar a CUDA dentro del ecosistema ROCM de AMD. Una gran cantidad de código de estilo CUDA se puede adaptar a HIP, pero la portabilidad de la fuente no significa que un binario se ejecute sin cambios en cada GPU.

Es posible que aún se requiera una compilación, pruebas y ajustes de rendimiento separados para cada arquitectura de destino. Los detalles están disponibles en la Documentación oficial de HIP.

HIP es una opción práctica cuando se debe admitir las GPU AMD o cuando un proyecto quiere reducir la dependencia de un solo proveedor de aceleración mientras se conserva un estilo de programación similar a CUDA.

Granos personalizados de Cupy

Cupy es mejor conocido por sus arreglos compatibles con Numpy y sus funciones numéricas aceleradas. También es compatible con la programación de GPU personalizada.

Los desarrolladores pueden usar:

  • ElementwiseKernel para expresiones personalizadas de elementos
  • ReductionKernel para reducciones personalizadas
  • RawKernel para núcleos escritos en CUDA C++
  • cupyx.jit.rawkernel para las definiciones de kernel JIT de estilo Python

Un kernel crudo simple se puede definir de la siguiente manera:

import cupy as cp

scale_kernel = cp.RawKernel(
    r'''
    extern "C" __global__
    void scale_values(
        const float* input,
        float* output,
        const float factor,
        const int size
    ) {
        int index =
            blockDim.x * blockIdx.x
            + threadIdx.x;

        if (index < size) {
            output[index] =
                factor * input[index];
        }
    }
    ''',
    "scale_values"
)

size = 1_000_000
threads_per_block = 256
blocks = (
    size + threads_per_block - 1
) // threads_per_block

input_values = cp.arange(
    size,
    dtype=cp.float32
)

output_values = cp.empty_like(
    input_values
)

scale_kernel(
    (blocks,),
    (threads_per_block,),
    (
        input_values,
        output_values,
        cp.float32(2.0),
        cp.int32(size)
    )
)

CUPY es útil cuando la mayoría de la aplicación ya usa matrices de GPU y solo las operaciones seleccionadas necesitan kernels personalizados.

numba-cuda

NUMBA-CUDA compila un subconjunto restringido de Python en kernels de GPU. Su modelo de programación sigue de cerca a CUDA C. Los desarrolladores definen una función con un decorador de kernel, calculan un índice de subprocesos e inician la función con una configuración de cuadrícula y bloque.

La documentación activa del proyecto está disponible en numba-cuda. Debido a que su estado de desarrollo puede cambiar, los equipos que crean software de larga duración deben revisar la guía actual de mantenimiento y migración antes de comprometerse con el marco.

un núcleo de partículas numba

El siguiente ejemplo educativo calcula la aceleración gravitacional con las interacciones directas de todos los pares:

import math
import numpy as np
from numba import cuda

@cuda.jit
def gravity_kernel(
    positions,
    masses,
    accelerations,
    gravitational_constant,
    softening_squared
):
    particle = cuda.grid(1)
    particle_count = positions.shape[0]

    if particle >= particle_count:
        return

    px = positions[particle, 0]
    py = positions[particle, 1]
    pz = positions[particle, 2]

    ax = 0.0
    ay = 0.0
    az = 0.0

    for other in range(particle_count):
        if other == particle:
            continue

        dx = positions[other, 0] - px
        dy = positions[other, 1] - py
        dz = positions[other, 2] - pz

        distance_squared = (
            dx * dx
            + dy * dy
            + dz * dz
            + softening_squared
        )

        inverse_distance = (
            1.0
            / math.sqrt(distance_squared)
        )

        inverse_distance_cubed = (
            inverse_distance
            * inverse_distance
            * inverse_distance
        )

        scale = (
            gravitational_constant
            * masses[other]
            * inverse_distance_cubed
        )

        ax += scale * dx
        ay += scale * dy
        az += scale * dz

    accelerations[particle, 0] = ax
    accelerations[particle, 1] = ay
    accelerations[particle, 2] = az


particle_count = 10_000
threads_per_block = 256

blocks = (
    particle_count
    + threads_per_block
    - 1
) // threads_per_block

host_positions = np.random.random(
    (particle_count, 3)
).astype(np.float32)

host_masses = np.ones(
    particle_count,
    dtype=np.float32
)

device_positions = cuda.to_device(
    host_positions
)

device_masses = cuda.to_device(
    host_masses
)

device_accelerations = cuda.device_array(
    (particle_count, 3),
    dtype=np.float32
)

gravity_kernel[
    blocks,
    threads_per_block
](
    device_positions,
    device_masses,
    device_accelerations,
    np.float32(1.0),
    np.float32(1e-4)
)

cuda.synchronize()

accelerations = (
    device_accelerations.copy_to_host()
)

Este kernel tiene complejidad computacional de O(N²). Es adecuado para explicar el mapeo de subprocesos, pero no es un algoritmo de producción eficiente para recuentos de partículas muy grandes.

Los grandes sistemas gravitatorios o de partículas a menudo requieren un método de árbol, un método multipolar rápido, contenedores espaciales o listas de vecinos. Una GPU no puede eliminar el problema de escala de un algoritmo ineficiente.

Deformación de Nvidia

NVIDIA WARP permite a los desarrolladores definir kernels fuertemente tipeados con sintaxis de Python. Warp compila esas funciones en CPU o código CUDA y proporciona primitivas especializadas para geometría, simulación, operaciones escasas, elementos finitos, optimización y diferenciación automática.

La nvidia warp product page y Warp GitHub Repository Incluye ejemplos de partículas, fluidos, mallas, optimización, simulación diferenciable y Programación de GPU basada en azulejos.

El mismo grano de partículas en urdimbre

import numpy as np
import warp as wp

wp.init()

@wp.kernel
def gravity_kernel(
    positions: wp.array(dtype=wp.vec3),
    masses: wp.array(dtype=wp.float32),
    accelerations: wp.array(dtype=wp.vec3),
    gravitational_constant: wp.float32,
    softening_squared: wp.float32
):
    particle = wp.tid()
    particle_count = positions.shape[0]

    position = positions[particle]
    acceleration = wp.vec3(0.0, 0.0, 0.0)

    for other in range(particle_count):
        if other != particle:
            displacement = (
                positions[other]
                - position
            )

            distance_squared = (
                wp.dot(
                    displacement,
                    displacement
                )
                + softening_squared
            )

            inverse_distance = (
                1.0
                / wp.sqrt(distance_squared)
            )

            inverse_distance_cubed = (
                inverse_distance
                * inverse_distance
                * inverse_distance
            )

            acceleration += (
                gravitational_constant
                * masses[other]
                * inverse_distance_cubed
                * displacement
            )

    accelerations[particle] = acceleration


particle_count = 10_000

host_positions = np.random.random(
    (particle_count, 3)
).astype(np.float32)

host_masses = np.ones(
    particle_count,
    dtype=np.float32
)

positions = wp.array(
    host_positions,
    dtype=wp.vec3,
    device="cuda"
)

masses = wp.array(
    host_masses,
    dtype=wp.float32,
    device="cuda"
)

accelerations = wp.zeros(
    particle_count,
    dtype=wp.vec3,
    device="cuda"
)

wp.launch(
    kernel=gravity_kernel,
    dim=particle_count,
    inputs=[
        positions,
        masses,
        accelerations,
        wp.float32(1.0),
        wp.float32(1e-4)
    ],
    device="cuda"
)

wp.synchronize()

Warp determina internamente una configuración de lanzamiento adecuada para un lanzamiento de kernel ordinario. Los usuarios avanzados aún pueden trabajar con dimensiones de bloques, gráficos de comandos, ejecución en mosaico y operaciones de dispositivos especializados cuando sea necesario.

Bucles de zancada

Un kernel no necesita un hilo asignado permanentemente para cada elemento. Un bucle de tracción de cuadrícula permite que cada subproceso procese varias entradas:

from numba import cuda

@cuda.jit
def scale_with_stride(
    input_values,
    output_values,
    factor
):
    index = cuda.grid(1)
    stride = cuda.gridsize(1)

    for position in range(
        index,
        input_values.size,
        stride
    ):
        output_values[position] = (
            factor
            * input_values[position]
        )

Este patrón separa el número de subprocesos lanzados del tamaño total de los datos. Es útil cuando se procesan arreglos muy grandes o se reutiliza una configuración de lanzamiento fijo.

Una introducción práctica a este patrón está disponible en la NUMBA y CUDA de Data Frog Tutorial.

Coalescente de memoria

Los kernels de GPU a menudo están limitados por el ancho de banda de la memoria en lugar del rendimiento aritmético. Los subprocesos dentro de un warp deberían acceder idealmente a las direcciones de memoria cercanas para que el hardware pueda combinar sus solicitudes.

Considere los datos de partículas almacenados como:

particle_0: x, y, z, mass
particle_1: x, y, z, mass
particle_2: x, y, z, mass

Este diseño de matriz de estructuras puede ser conveniente para el código orientado a objetos. Un diseño de estructura de arreglos almacena matrices contiguas separadas:

x_positions[]
y_positions[]
z_positions[]
masses[]

El segundo diseño puede proporcionar un mejor coalescencia cuando cada subproceso lee el mismo campo para una partícula diferente. El mejor diseño aún depende de qué campos se accedan juntos.

Los tipos de vectores no corrigen automáticamente un patrón de acceso deficiente. Los desarrolladores deben inspeccionar las direcciones reales solicitadas por los hilos vecinos.

memoria compartida

La memoria compartida es una pequeña área de memoria de baja latencia accesible por subprocesos en el mismo bloque. Puede reducir las lecturas repetidas de la memoria global.

Una plantilla unidimensional puede cargar un bloque de valores y su halo en la memoria compartida:

import numpy as np
from numba import cuda, float32

BLOCK_SIZE = 256

@cuda.jit
def three_point_stencil(
    input_values,
    output_values
):
    shared = cuda.shared.array(
        shape=BLOCK_SIZE + 2,
        dtype=float32
    )

    local_index = cuda.threadIdx.x
    global_index = cuda.grid(1)
    size = input_values.size

    center = local_index + 1

    if global_index < size:
        shared[center] = (
            input_values[global_index]
        )
    else:
        shared[center] = 0.0

    if local_index == 0:
        left_index = global_index - 1

        shared[0] = (
            input_values[left_index]
            if left_index >= 0
            else 0.0
        )

    if local_index == BLOCK_SIZE - 1:
        right_index = global_index + 1

        shared[BLOCK_SIZE + 1] = (
            input_values[right_index]
            if right_index < size
            else 0.0
        )

    cuda.syncthreads()

    if global_index < size:
        output_values[global_index] = (
            shared[center - 1]
            + shared[center]
            + shared[center + 1]
        ) / 3.0

La llamada a cuda.syncthreads() garantiza que cada subproceso termine de cargar sus valores antes de que comience el cálculo de la plantilla.

La memoria compartida no debe usarse automáticamente. La asignación excesiva de memoria compartida puede reducir la ocupación, agregar el costo de sincronización y hacer que el kernel sea más lento. Se requiere perfilado.

Fusión de grano

Las operaciones de matriz separadas a menudo crean varios lanzamientos de kernel y matrices intermedios:

velocity += dt * acceleration
position += dt * velocity
energy = compute_energy(position, velocity)

Un kernel fusionado puede calcular las tres actualizaciones mientras los valores necesarios permanecen en los registros. Esto reduce la sobrecarga de lanzamiento y el tráfico de memoria global.

La fusión es más beneficiosa cuando las operaciones son simples y están vinculadas a la memoria. Fusionar demasiado trabajo puede aumentar el uso de registros, reducir la ocupación y hacer que el kernel sea difícil de mantener.

Elegir el tamaño del bloque

Los valores como 128 o 256 hilos por bloque son puntos de partida razonables, no como Optima universal.

El mejor tamaño de bloque depende de:

  • Registros utilizados por hilo
  • Memoria compartida utilizada por bloque
  • Divergencia de ramas
  • Mezcla de instrucciones
  • Comportamiento de acceso a la memoria
  • La arquitectura de GPU objetivo

Una regla como lanzar al menos el doble de bloques que los multiprocesadores de transmisión puede ser un experimento inicial útil, pero no garantiza el máximo rendimiento. Las calculadoras de ocupación y las herramientas de perfilado deben guiar la configuración final.

Física diferenciable

La simulación diferenciable calcula cómo cambia una salida con respecto a las entradas como las propiedades del material, las fuerzas, la geometría o las condiciones iniciales.

Warp puede registrar las operaciones del kernel compatible y ejecutar la diferenciación automática en modo inverso. Esto se puede utilizar para problemas inversos, optimización del diseño, estimación de parámetros e integración con flujos de trabajo de aprendizaje automático.

Un patrón simplificado utiliza una cinta:

with wp.Tape() as tape:
    wp.launch(
        kernel=simulation_kernel,
        dim=element_count,
        inputs=[state, parameters],
        outputs=[result],
        device="cuda"
    )

    wp.launch(
        kernel=loss_kernel,
        dim=element_count,
        inputs=[result, target],
        outputs=[loss],
        device="cuda"
    )

tape.backward(loss)

No se garantiza la diferenciación automática para cada kernel. Las sobreescrituras en el lugar, los atómicos no deterministas, el código nativo externo, la lógica discontinua y las operaciones no admitidas pueden requerir reformulación o gradientes personalizados.

Esta capacidad se conecta directamente con aplicaciones como la optimización adjunta y Redes neuronales informadas por la física.

Cuando una GPU puede no ser más rápida

No existe un umbral de tamaño de problema universal en el que una GPU se vuelva más rápida que una CPU. El cruce depende del hardware, la precisión, el movimiento de datos, la estructura del algoritmo, la calidad del compilador y la frecuencia con la que se reutilizan los mismos datos residentes en el dispositivo.

carga de trabajo Adecuación de GPU Consideración principal
Actualización de partículas grandes A menudo fuerte muchas operaciones independientes similares
Plantilla grande regular A menudo fuerte Acceso predecible a la memoria en paralelo
Álgebra lineal densa Fuerte al usar bibliotecas optimizadas Alta intensidad aritmética
pequeña simulación dependiente Lanzamiento y transferencia Los gastos generales pueden dominar
Trasversal de gráficos irregulares Mezclado Divergencia y acceso a la memoria impredecible
Flujo de trabajo enlazado a E/S por lo general limitado La GPU no puede eliminar un cuello de botella de almacenamiento o de red
Algoritmo fuertemente secuencial por lo general débil trabajo independiente insuficiente

El flujo de trabajo correcto es medir la versión de la CPU, la versión de la biblioteca GPU y la versión del kernel personalizado con datos representativos. La Guía de perfiles de rendimiento de computación científica explica cómo identificar el cuello de botella real antes de optimizarlo.

Evitar la transferencia de gastos generales

Las transferencias repetidas entre la memoria del dispositivo y del dispositivo pueden eliminar el beneficio del cálculo de la GPU.

Un bucle ineficiente puede seguir este patrón:

  1. Copiar datos a la GPU.
  2. Ejecute un kernel pequeño.
  3. Copie el resultado a la CPU.
  4. modificarlo en la CPU.
  5. cópielo de nuevo a la GPU.

Un mejor diseño mantiene el estado de simulación en el dispositivo durante muchos pasos y transfiere solo la salida requerida para la visualización, los puntos de control o el análisis.

Los flujos asincrónicos, la memoria del host anclada, la comunicación superpuesta y las funciones de memoria unificada pueden ayudar, pero solo deben introducirse después de que se hayan medido los costos de transferencia ordinarios.

El estudio de caso XLB

El proyecto XLB de Autodesk Research proporciona un ejemplo útil de simulación de GPU nativa de Python. XLB es una biblioteca Boltzmann de código abierto con varios backends computacionales, incluido Nvidia Warp.

En las configuraciones de referencia informadas por Autodesk Research y NVIDIA, el backend warp alcanzó un rendimiento cercano a la implementación comparada de C++/OpenCl FluidX3D para un caso específico de cavidad controlada por TAPA. Una comparación separada informó una aceleración de ocho veces aproximada sobre el backend JAX de XLB en el hardware y las configuraciones seleccionados.

El equipo también demostró un enfoque fuera de núcleo en un grupo GH200 de ocho nodos con un dominio de aproximadamente 50 mil millones de celdas de red. Estos resultados se aplican a las opciones de solución, problema, hardware y la implementación reportada. No deben tratarse como garantías generales de aceleración para los núcleos de Python.

El estudio de caso completo está disponible en Autodesk Research lleva la velocidad de la deformación a la dinámica de fluidos computacional en NVIDIA GH200. El código de proyecto actual está disponible en el repositorio Autodesk XLB.

Cómo comparar correctamente un kernel

Las operaciones de GPU suelen ser asíncronas. Medir solo la duración de la llamada a la función de Python puede informar el tiempo necesario para poner en cola el kernel en lugar del tiempo necesario para ejecutarlo.

Un punto de referencia básico debe:

  1. Ejecute el kernel varias veces para activar la compilación y el calentamiento.
  2. Sincronice antes de iniciar el temporizador.
  3. Ejecute varias iteraciones medidas.
  4. Sincronice antes de detener el temporizador.
  5. Reporte el promedio y la variación entre las repeticiones.
  6. Separe el tiempo de transferencia de datos desde el tiempo de ejecución del kernel.
  7. Verifique que las versiones de CPU y GPU produzcan resultados equivalentes.
import time
from numba import cuda

# Warm-up and JIT compilation
kernel[blocks, threads](*arguments)
cuda.synchronize()

start = time.perf_counter()

for _ in range(100):
    kernel[blocks, threads](*arguments)

cuda.synchronize()

elapsed = time.perf_counter() - start
average = elapsed / 100

print("Average kernel time:", average)

La comparación debe usar un tamaño de simulación realista e incluir el flujo de trabajo completo cuando los costos de transferencia o preprocesamiento importan.

La metodología de precisión de trabajo también puede comparar el tiempo de ejecución con el error numérico. La Documentación de referencia de SCML proporciona ejemplos de este estilo de evaluación.

un flujo de trabajo de desarrollo práctico

  1. Implemente y verifique una versión de referencia de CPU clara.
  2. Perfil de la aplicación para localizar la operación dominante.
  3. Pruebe una biblioteca de GPU optimizada antes de escribir un kernel.
  4. Mantenga los datos reutilizados con frecuencia en el dispositivo.
  5. Escriba el kernel correcto más simple.
  6. Validarlo contra el resultado de la CPU.
  7. Mida las transferencias de memoria y la ejecución por separado.
  8. Inspeccione el uso de coalescencia, ocupación, divergencia y registro.
  9. Pruebe la memoria compartida o la fusión solo cuando la creación de perfiles los admita.
  10. Compare varios tamaños de problemas y arquitecturas de GPU.

Elegir un marco

Requisito Posible punto de partida
Flujo de trabajo de GPU de estilo numpy existente CUPY con operaciones integradas o kernels personalizados
Pequeña cantidad de kernels de Python de estilo CUDA Numba-CUDA, después de revisar su estado de soporte actual
Simulación, geometría y kernels diferenciables Deformación de Nvidia
Control máximo específico de NVIDIA cudá C++
Objetivo de GPU AMD con fuente CUDA Cadera y Rocm
C++ portátil en varios backends SYCL, Kokkos u otra capa de rendimiento-portabilidad

Errores comunes de programación del kernel

  • Mover datos entre la CPU y la GPU dentro de cada paso del tiempo
  • Ignorar hilos fuera de los límites
  • Usando un algoritmo ineficiente y esperando hardware para arreglar su escala
  • Acceder a la memoria con un diseño no coalescido
  • Agregar memoria compartida sin medir si ayuda
  • lanzando muchos granos pequeños en lugar de considerar la fusión
  • Uso de registros excesivos o memoria compartida por bloque
  • Benchmarking código asíncrono sin sincronización
  • Comparación de salidas sin comprobar la precisión numérica
  • Uso de reclamos de rendimiento fijo en diferentes GPU y cargas de trabajo
  • Suponiendo que la fuente portátil proporcione un rendimiento portátil
  • Aplicación de diferenciación automática a operaciones in situ no soportadas

Guías relacionadas

Conclusión

Los kernels de GPU personalizados permiten a los desarrolladores científicos traducir partículas, celdas de cuadrícula, reglas de materiales y otras operaciones físicas directamente en funciones masivamente paralelas. Proporcionan control sobre la indexación, el acceso a la memoria, la sincronización, la fusión del kernel y la optimización específica del dispositivo.

CUPY ofrece operaciones de matriz de alto nivel e interfaces de kernel personalizadas. Numba-Cuda proporciona un modelo similar a CUDA en Python, mientras que Nvidia Warp agrega primitivas de simulación, operaciones de mosaico y diferenciación automática. CUDA y HIP proporcionan control C++ de nivel inferior, y los marcos portátiles de rendimiento admiten objetivos de hardware más amplios.

Las mayores ganancias de rendimiento rara vez provienen de cambiar la sintaxis por sí sola. Provienen de elegir un algoritmo paralelo, mantener datos en el dispositivo, reducir el tráfico de memoria, utilizar un diseño de datos adecuado y eliminar operaciones intermedias innecesarias.

Un kernel personalizado debe desarrollarse solo después de que la creación de perfiles identifique un cuello de botella real. Debe validarse contra una solución de referencia y compararse con la sincronización, los tamaños de problemas representativos y la contabilidad completa de los costos de transferencia de memoria.

Cuando se cumplen esas condiciones, las herramientas de kernel basadas en Python pueden respaldar una simulación física seria sin obligar a los investigadores a mover cada parte de la aplicación a C++ de bajo nivel.