Programmation du noyau GPU pour la simulation physique personnalisée : écriture de noyaux personnalisés hautes performances en Python

L’accélération des GPU en informatique scientifique commence souvent par des bibliothèques de haut niveau. Un chercheur remplace NumPy par CUPY, utilise une FFT accélérée ou appelle une routine d’algèbre linéaire activée par GPU. Cela peut produire une amélioration substantielle sans nécessiter de connaissances détaillées sur le matériel GPU.

La programmation du noyau GPU personnalisé va d’un niveau plus profondément. Au lieu d’appeler une opération prédéfinie, le développeur écrit la fonction que chaque thread GPU exécute. Cela permet de contrôler l’indexation des threads, la disposition des données, les transferts de mémoire, la synchronisation, la mémoire partagée et le nombre d’opérations effectuées lors de chaque lancement du noyau.

Ce contrôle est utile lorsqu’un modèle physique contient des modèles d’accès qui ne peuvent pas être exprimés efficacement par le biais d’opérations de tableau ordinaires. Les interactions de particules, les pochoirs de différences finies, les règles de collision, les modèles de fluides basés sur le réseau, les mises à jour des matériaux et les requêtes de géométrie sont des exemples courants.

Un noyau personnalisé n’est pas automatiquement plus rapide qu’une bibliothèque optimisée. Cela devient utile lorsque l’algorithme expose suffisamment de travail parallèle et lorsque le développeur peut réduire le trafic mémoire, combiner les opérations ou adapter le calcul à l’architecture GPU.

Qu’est-ce qu’un noyau GPU ?

Un noyau GPU est une fonction exécutée par de nombreux threads légers. Chaque thread gère normalement une particule, une cellule de grille, un élément, une face, un pixel ou une entrée dans un tableau.

Le programme hôte lance le noyau avec une grille de blocs de thread :

kernel[blocks, threads_per_block](arguments)

Chaque thread détermine les données qu’il doit traiter à partir de ses indices de blocs et de threads. Un index unidimensionnel typique est :

index =
    block_index * block_size
    + thread_index

Les threads sont organisés en groupes appelés blocs. Les threads à l’intérieur d’un seul bloc peuvent synchroniser et échanger des données via la mémoire partagée. Les blocs séparés doivent généralement s’exécuter de manière indépendante.

Lorsque les noyaux personnalisés valent la peine d’être écrits

Les bibliothèques de GPU de haut niveau devraient généralement être la première option. Ils fournissent déjà une multiplication matricielle optimisée, des FFT, des réductions, des opérations rares, une génération de nombres aléatoires et de nombreuses fonctions par élément.

Un noyau personnalisé devient utile lorsque :

• La simulation utilise un modèle d’accès à la mémoire non standard.
• Plusieurs petites opérations peuvent être fusionnées en un seul passage sur la mémoire.
• Les threads voisins ont à plusieurs reprises besoin des mêmes données locales.
• Une règle spécialisée de particule, de collision ou de pochoir domine le temps d’exécution.
• Les baies intermédiaires consomment trop de mémoire de périphérique.
• La simulation nécessite une opération différenciable personnalisée.
• Une bibliothèque existante ne peut pas exprimer efficacement la limite ou la logique matérielle requise.

Avant d’écrire un noyau, profilez le programme. L’optimisation d’une fonction visuellement complexe qui ne représente qu’une petite partie du temps d’exécution ne produira pas une accélération globale significative.

Modèles de programmation GPU

Les logiciels scientifiques peuvent accéder aux GPU via plusieurs niveaux d’abstraction.

ENCCS Introduction à la programmation GPU Les modèles fournissent une comparaison plus large des approches CUDA, HIP, OpenCL, SYCL, Kokkos, OpenMP et associées.

Cuda

CUDA fournit un accès direct à la programmation NVIDIA GPU via C++, Fortran et plusieurs liaisons de langage. Il expose les grilles, les blocs, les déformations, la mémoire partagée, les flux, les événements, les bibliothèques de périphériques et les outils d’optimisation spécifiques au matériel.

CUDA est souvent sélectionné lorsque des performances maximales spécifiques à NVIDIA et des outils matures sont plus importants que la portabilité matérielle. Les principaux coûts sont la gestion de la mémoire manuelle, un modèle de programmation de niveau inférieur et la dépendance à la plate-forme NVIDIA.

Hanche et ROCM

HIP fournit un modèle de programmation C++ de type CUDA dans l’écosystème ROCM d’AMD. Une grande quantité de code de style CUDA peut être adaptée à la hanche, mais la portabilité de la source ne signifie pas qu’un binaire s’exécute inchangé sur chaque GPU.

Une compilation, des tests et un réglage des performances séparés peuvent toujours être nécessaires pour chaque architecture cible. Les détails sont disponibles dans la documentation officielle de la hanche.

HIP est une option pratique lorsque les GPU AMD doivent être pris en charge ou lorsqu’un projet souhaite réduire la dépendance à un seul fournisseur d’accélérateurs tout en conservant un style de programmation de type CUDA.

Cupy Cernels personnalisés

Cupy est surtout connu pour ses tableaux compatibles avec NumPy et ses fonctions numériques accélérées. Il prend également en charge la programmation GPU personnalisée.

Les développeurs peuvent utiliser :

• ElementwiseKernel pour les expressions personnalisées par élément
• ReductionKernel pour les réductions personnalisées
• RawKernel pour les noyaux écrits en CUDA C++
• cupyx.jit.rawkernel pour les définitions du noyau JIT de style Python

Un simple noyau brut Cupy peut être défini comme suit :

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 est utile lorsque la plupart de l’application utilise déjà des baies de GPU et que seules les opérations sélectionnées ont besoin de noyaux personnalisés.

Numba-cuda

numba-cuda compile un sous-ensemble restreint de python dans les noyaux GPU. Son modèle de programmation suit de près CUDA C. Les développeurs définissent une fonction avec un décorateur de noyau, calculent un index de fil et lancent la fonction avec une configuration de grille et de bloc.

La documentation active du projet est disponible à numba-cuda. Étant donné que son statut de développement peut changer, les équipes créant des logiciels de longue durée doivent revoir les directives de maintenance et de migration actuelles avant de s’engager dans le cadre.

Un noyau de particules Numba

L’exemple éducatif suivant calcule l’accélération gravitationnelle avec des interactions directes toutes paires :

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()
)

Ce noyau a une complexité de calcul de O(N²). Il est adapté pour expliquer la cartographie des threads, mais ce n’est pas un algorithme de production efficace pour les très gros nombres de particules.

Les grands systèmes de gravitation ou de particules nécessitent souvent une méthode d’arbre, une méthode multipôle rapide, des bacs spatiaux ou des listes de voisins. Un GPU ne peut pas supprimer le problème de mise à l’échelle d’un algorithme inefficace.

NVIDIA WARP

NVIDIA WARP permet aux développeurs de définir des noyaux fortement typés avec la syntaxe Python. WARP compile ces fonctions en CPU ou CUDA et fournit des primitives spécialisées pour la géométrie, la simulation, les opérations rares, les éléments finis, l’optimisation et la différenciation automatique.

La page produit NVIDIA WARP et Référentiel GitHub de Warp Inclut des exemples de particules, de fluides, de maillages, d’optimisation, de simulation différenciable et Programmation GPU basée sur des tuiles.

Le même noyau de particules dans WARP

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 détermine en interne une configuration de lancement appropriée pour un lancement de noyau ordinaire. Les utilisateurs avancés peuvent toujours travailler avec des dimensions de blocs, des graphiques de commandes, une exécution en mosaïque et des opérations de périphérique spécialisées si nécessaire.

Boucles de bordure de grille

Un noyau n’a pas besoin d’un thread assigné de façon permanente pour chaque élément. Une boucle Grid-Stride permet à chaque thread de traiter plusieurs entrées :

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]
        )

Ce modèle sépare le nombre de threads lancés de la taille totale des données. Il est utile lors du traitement de très grands tableaux ou de la réutilisation d’une configuration de lancement fixe.

Une introduction pratique à ce modèle est disponible dans le Numba et le CUDA de la grenouille de données Tutoriel.

Coalition de mémoire

Les noyaux GPU sont souvent limités par la bande passante de mémoire plutôt que par le débit arithmétique. Les threads à l’intérieur d’un distorsion doivent idéalement accéder aux adresses mémoire à proximité afin que le matériel puisse combiner leurs demandes.

Considérez les données de particules stockées comme :

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

Cette disposition de tableau de structures peut être pratique pour le code orienté objet. Une mise en page de la structure des tableaux stocke des tableaux contigus distincts :

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

La deuxième mise en page peut fournir une meilleure coalescence lorsque chaque thread lit le même champ pour une particule différente. La meilleure mise en page dépend toujours des champs accessibles ensemble.

Les types de vecteurs ne corrigent pas automatiquement un mauvais modèle d’accès. Les développeurs doivent inspecter les adresses réelles demandées par les threads voisins.

Mémoire partagée

La mémoire partagée est une petite zone de mémoire à faible latence accessible par les threads dans le même bloc. Il peut réduire les lectures répétées de la mémoire globale.

Un pochoir unidimensionnel peut charger un bloc de valeurs et son halo dans la mémoire partagée :

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

L’appel à cuda.syncthreads() garantit que chaque thread finit de charger ses valeurs avant le début du calcul du pochoir.

La mémoire partagée ne doit pas être utilisée automatiquement. Une allocation de mémoire partagée excessive peut réduire l’occupation, augmenter les coûts de synchronisation et ralentir le noyau. Le profilage est requis.

Fusion du noyau

Des opérations de tableaux distinctes créent souvent plusieurs lancements de noyau et des tableaux intermédiaires :

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

Un noyau fusionné peut calculer les trois mises à jour alors que les valeurs nécessaires restent dans les registres. Cela réduit les frais généraux de lancement et le trafic de mémoire globale.

La fusion est plus avantageuse lorsque les opérations sont simples et liées à la mémoire. Fusionner trop de travail peut augmenter l’utilisation des registres, réduire l’occupation et rendre le noyau difficile à entretenir.

Choisir la taille du bloc

Des valeurs telles que 128 ou 256 threads par bloc sont des points de départ raisonnables, et non des optima universels.

La meilleure taille de bloc dépend de :

• Registres utilisés par fil
• Mémoire partagée utilisée par bloc
• Divergence de branche
• mélange d’instructions
• Comportement d’accès à la mémoire
• L’architecture GPU cible

Une règle telle que le lancement d’au moins deux fois plus de blocs que les multiprocesseurs de streaming peut être une expérience initiale utile, mais elle ne garantit pas des performances maximales. Les calculatrices d’occupation et les outils de profilage doivent guider la configuration finale.

Physique différentiable

La simulation différentiable calcule comment une sortie change par rapport aux entrées telles que les propriétés des matériaux, les forces, la géométrie ou les conditions initiales.

WARP peut enregistrer les opérations du noyau pris en charge et exécuter la différenciation automatique en mode inverse. Cela peut être utilisé pour les problèmes inverses, l’optimisation de la conception, l’estimation des paramètres et l’intégration avec les flux de travail d’apprentissage automatique.

Un motif simplifié utilise une bande :

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)

La différenciation automatique n’est pas garantie pour chaque noyau. Les écrasements sur place, les atomes non déterministes, le code natif externe, la logique discontinue et les opérations non pris en charge peuvent nécessiter une reformulation ou des gradients personnalisés.

Cette fonctionnalité est directement connectée à des applications telles que l’optimisation des associations et les réseaux de neurones à base de physique.

Lorsqu’un GPU n’est peut-être pas plus rapide

Il n’y a pas de seuil universel de taille de problème auquel un GPU devient plus rapide qu’un processeur. Le crossover dépend du matériel, de la précision, du mouvement des données, de la structure des algorithmes, de la qualité du compilateur et de la fréquence à laquelle les mêmes données de résidents sont réutilisées.

Le flux de travail correct consiste à mesurer la version du processeur, la version de la bibliothèque GPU et la version du noyau personnalisé avec des données représentatives. Le Guide de profilage de performances de calcul scientifique explique comment identifier le véritable goulot d’étranglement avant de l’optimiser.

Éviter les frais généraux de transfert

Les transferts répétés entre l’hôte et la mémoire de l’appareil peuvent éliminer les avantages du calcul du GPU.

Une boucle inefficace peut suivre ce modèle :

1. Copiez les données dans le GPU.
2. Exécutez un petit noyau.
3. Copiez le résultat sur la CPU.
4. Modifiez-le sur la CPU.
5. Copiez-le dans le GPU.

Une meilleure conception maintient l’état de simulation sur l’appareil pendant de nombreuses étapes de temps et ne transfère que la sortie requise pour la visualisation, le point de contrôle ou l’analyse.

Les flux asynchrones, la mémoire hôte épinglée, la communication qui se chevauche et les fonctionnalités de mémoire unifiée peuvent aider, mais elles ne doivent être introduites qu’après la mesure des coûts de transfert ordinaires.

L’étude de cas XLB

Le projet XLB d’Autodesk Research fournit un exemple utile de simulation GPU native Python. XLB est une bibliothèque de réseau open source Boltzmann avec plusieurs backends informatiques, dont NVIDIA WARP.

Dans les configurations de référence signalées par Autodesk Research et NVIDIA, le backend Warp a atteint des performances proches de la mise en œuvre comparative C++/OpenCL FluidX3D pour un cas de cavité spécifiquement piloté par LID. Une comparaison distincte a rapporté une accélération approximative de huit fois par rapport au backend JAX de XLB sur le matériel et les configurations sélectionnés.

L’équipe a également démontré une approche erronée sur un cluster GH200 à huit nœuds avec un domaine d’environ 50 milliards de cellules de réseau. Ces résultats s’appliquent au solveur, au problème, au matériel et aux choix d’implémentation signalés. Ils ne doivent pas être traités comme des garanties d’accélération générales pour les noyaux Python.

L’étude de cas complète est disponible dans La recherche d’Autodesk apporte une vitesse de distorsion à la dynamique des fluides de calcul sur NVIDIA GH200. Le code de projet actuel est disponible dans le Autodesk XLB Repository.

Comment comparer correctement un noyau

Les opérations GPU sont généralement asynchrones. Mesure uniquement la durée de l’appel de fonctions Python peut indiquer le temps nécessaire pour mettre le noyau en file plutôt que le temps nécessaire pour l’exécuter.

Une référence de base devrait :

1. Exécutez le noyau plusieurs fois pour déclencher la compilation et l’échauffement.
2. Synchronisez avant de démarrer la minuterie.
3. Exécutez plusieurs itérations mesurées.
4. Synchronisez avant d’arrêter la minuterie.
5. Signaler la moyenne et la variation des répétitions.
6. Séparez le temps de transfert de données du temps d’exécution du noyau.
7. Vérifiez que les versions du CPU et du GPU produisent des résultats équivalents.

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 comparaison doit utiliser une taille de simulation réaliste et inclure le flux de travail complet lorsque les coûts de transfert ou de prétraitement sont importants.

La méthodologie de précision du travail peut également comparer le temps d’exécution à l’erreur numérique. La Documentation de référence SCIML fournit des exemples de ce style d’évaluation.

Un flux de travail de développement pratique

1. Implémentez et vérifiez une version claire de la référence CPU.
2. Profilez l’application pour localiser l’opération dominante.
3. Essayez une bibliothèque GPU optimisée avant d’écrire un noyau.
4. Conservez les données fréquemment réutilisées sur l’appareil.
5. Écrivez le noyau correct le plus simple.
6. Validez-le par rapport au résultat du CPU.
7. Mesurez les transferts de mémoire et l’exécution séparément.
8. Inspectez la coalescence, l’occupation, la divergence et l’utilisation des registres.
9. Testez la mémoire partagée ou la fusion uniquement lorsque le profilage les prend en charge.
10. Benchmark plusieurs tailles de problèmes et architectures GPU.

Choisir un cadre

Erreurs courantes de programmation du noyau

• Déplacement des données entre le CPU et le GPU à chaque pas de temps
• Ignorer les threads hors limites
• Utiliser un algorithme inefficace et s’attendre à ce que le matériel corrige sa mise à l’échelle
• Accéder à la mémoire avec une mise en page non fusionnée
• Ajouter de la mémoire partagée sans mesurer si cela aide
• Lancement de nombreux minuscules noyaux au lieu d’envisager la fusion
• Utilisation de registres excessifs ou de mémoire partagée par bloc
• Analyse comparative du code asynchrone sans synchronisation
• Comparaison des sorties sans vérification de la précision numérique
• Utilisation de revendications de performances fixes sur différents GPU et charges de travail
• En supposant que la source portable offre des performances portables
• Application de la différenciation automatique aux opérations sur place non prises en charge

Guides connexes

• scientifique accéléré par GPU Informatique : CUPY, NUMBA et CUDF comparés – Comparez les outils Python de haut niveau pour les charges de travail GPU.
• Accélération GPU pour les simulations FIPY : cupy et Numba Integration Guide — Explorez les composants GPU possibles autour d’un flux de travail Fipy.
• Réseaux de neurones à base de physique pour les simulations scientifiques — Passez en revue la relation entre les modèles physiques et l’apprentissage différentiable.
• Quand utiliser FEM, FVM et FDM Discrétisation avant d’optimiser sa mise en œuvre.
• Profilage des performances : identifier les goulots d’étranglement dans le code scientifique : localiser les limitations du processeur, de la mémoire, du transfert et des E/S.
• Flows Python HPC : depuis un ordinateur portable au superordinateur — Planifiez la progression du développement local vers les grands systèmes.

Conclusion

Les noyaux GPU personnalisés permettent aux développeurs scientifiques de traduire directement des particules, des cellules de la grille, des règles de matériaux et d’autres opérations physiques en fonctions massivement parallèles. Ils permettent de contrôler l’indexation, l’accès à la mémoire, la synchronisation, la fusion du noyau et l’optimisation spécifique à l’appareil.

CUPY propose à la fois des opérations de baies de haut niveau et des interfaces de noyau personnalisées. Numba-CUDA fournit un modèle de type CUDA en Python, tandis que NVIDIA WARP ajoute des primitives de simulation, des opérations de tuiles et une différenciation automatique. CUDA et HIP offrent un contrôle C++ de niveau inférieur, et les cadres portables de performances prennent en charge des objectifs matériels plus larges.

Les gains de performances les plus importants proviennent rarement de la modification de la syntaxe seule. Ils proviennent du choix d’un algorithme parallèle, de la conservation des données sur l’appareil, de la réduction du trafic mémoire, de l’utilisation d’une disposition de données appropriée et de l’élimination des opérations intermédiaires inutiles.

Un noyau personnalisé ne doit être développé qu’après le profilage identifiant un véritable goulot d’étranglement. Il doit être validé par rapport à une solution de référence et étalonné avec la synchronisation, la taille représentative des problèmes et la comptabilisation complète des coûts de transfert de mémoire.

Lorsque ces conditions sont remplies, les outils de noyau basés sur Python peuvent prendre en charge la simulation physique sérieuse sans obliger les chercheurs à déplacer chaque partie de l’application en C++ de bas niveau.