Accès aux résultats GPU et manipulation des données

Ce script permet d’accéder aux résultats des simulations 2D effectuées sur GPU et de manipuler les données pour une analyse ultérieure. Il est conçu pour être utilisé avec les fichiers de résultats générés par les simulations.

Import des modules nécessaires

import _add_path # for debugging purposes only - must be removed in production

from wolfhece import is_enough
try:
    if not is_enough('2.2.33'):
        raise ImportError("HECE version 2.2.33 or higher is required.")
except ImportError as e:
    print(f"Error: {e}")
    raise

from wolfhece.Results2DGPU import wolfres2DGPU, OneWolfResult, views_2D, WolfArray
from wolfgpu.results_store import ResultsStore, ResultType

from wolfhece.pydownloader import toys_gpu_dataset, GITLAB_EXAMPLE_GPU

import numpy as np

Accès direct depuis le ResultStore

Le ResultStore est la classe qui gère les opérations de lecture et d’écriture des résultats des simulations GPU.

Même si d’autres routines sont disponibles, nous demandons à l’utilisateur de ne pas les utiliser directement, mais de passer par le ResultStore. Il offre en effet une interface fixe pour accéder aux résultats, ce qui permet de garantir la compatibilité avec les futures versions du code, notamment le passage de fichiers “npz” au format “npy”.

En interne, le code GPU utilise un tuilage pour convertir les données de modélisation en un format rectangulaire dense, mieux adapté au stockage et aux calculs sur GPU. En outre, pour le “npz”, seules les valeurs non nulles sont stockées (via le format CSR - Compressed Sparse Row), ce qui permet de réduire la taille des fichiers.

L’utilisateur final voudra habituellement quant à lui disposer des résultats sous forme géoréférencée, ce qui nécessite une conversion des données.

Cette conversion de données n’est pas “gratuite” en termes de temps de calcul. Si on souhaite accéder rapidement aux résultats, le format “npy” est plus adapté, car il stocke les matrices tuilées sans cette étape de conversion CSR.

rs = toys_gpu_dataset('channel_w_archbridge_fully_man004') # Load a GPU dataset from the web and obtain a ResultsStore object

assert isinstance(rs, ResultsStore), "The object should be an instance of ResultsStore"

# An error message will be raised as the bridge_deck is not available in the dataset -- It is normal !

ERROR:root:HTTP error occurred while downloading https://gitlab.uliege.be/HECE/wolfgpu_examples/-/raw/main/channel_w_archbridge_fully_man004/bridge_deck.npy: 404 Client Error: Not Found for url: https://gitlab.uliege.be/HECE/wolfgpu_examples/-/raw/main/channel_w_archbridge_fully_man004/bridge_deck.npy

Obtenir le nombre de résultats disponibles

rs.nb_results # Get the number of results available in the dataset

6

Récupérer les valeurs associées aux sauvegardes

Il est possible de récupérer les séries temporelles des valeurs associées aux résultats. Les valeurs possibles sont dans la classe ResultType et sont les suivantes :

• T : The simulation time (the time inside the simulation)

• LAST_DELTA_T : The last Δt computed at the moment the record is taken.

• NB_ACTIVE_TILES : The number of active tiles in the simulation. An active tile is a tile where there is some water.

• STEP_NUM : The number of the simulation step.

• CLOCK_T : The time spent into the simulation (wall clock time)

• NB_MOSTLY_DRY_MESHES : Number of mostly dry mesh. A mostly dry mesh is a mesh where water depth is not zero but very small (< 1e-5 m).

• NB_WET_MESHES : Number of wet meshes. A wet mesh is a mesh where water depth is above zero.

print(rs.get_named_series(ResultType.T))
print(rs.get_named_series(ResultType.NB_WET_MESHES))

[0.0, 40.60536419227719, 81.21073305234313, 121.81609793752432, 162.4214591383934, 203.02682216465473]
[8466.0, 8466.0, 8466.0, 8466.0, 8466.0, 8466.0]

Récupérer les résultats

• soit le dernier pas disponible

• soit un pas spécifique

res = rs.get_last_result()

# returns [t, dt, n_iter_dry_up_euler, n_iter_dry_up_rk, h, qx, qy]

h = rs.get_last_named_result(ResultType.H)  # Get the last result for the water depth
qx, qy = rs.get_last_named_result([ResultType.QX, ResultType.QY])  # Get the last result for the x-component and y-component of the unit discharge - we can mix ResultType values in a list

Les résultats sont retournés sous la forme de matrices numpy NON géoréférencées (Numpy ne permet pas ce type d’opération).

ATTENTION que les formes sont transposées par rapport à ce que Wolf a l’habitude de faire.

En réalité, la simulation a été définie avec 500 mailles selon X et 21 mailles selon Y.

Dans ces résultats “bruts”, les mailles sont stockées en lignes, ce qui signifie que la première dimension correspond à Y et la seconde à X.

print(type(h), h.shape)
print(type(qx), qx.shape)
print(type(qy), qy.shape)

<class 'numpy.ndarray'> (21, 500)

Graphiques

Il est tout à fait possible de visualiser les résultats obtenus avec des bibliothèques comme Matplotlib ou Seaborn.

import matplotlib.pyplot as plt

plt.imshow(qx, cmap='viridis', origin='lower') # Le stockage sous la forme (21,500) permet d'utiliser imshow directement, tandis que Wolf devrait les transposer (ce qui est la plupart du temps transparent pour l'utilisateur, et sans surcoût).

<matplotlib.image.AxesImage at 0x25eeaf1e2d0>

Récupérer un résultat spécifique

h = rs.get_named_result(ResultType.H, 1)  # Get the first result for the water depth - index is 1-based (like Fortran code)

Travailler avec l’objet “wolfres2DGPU”

“wolfres2DGPU” est un objet qui permet de manipuler les résultats des simulations 2D sur GPU de manière plus complète. Il est basé sur l’objet “Wolfresults_2D” qui permet entre autre de gérer les résultats CPU-Fortran.

En interne, “wolfres2DGPU” utilise le “ResultStore” pour accéder aux résultats, mais il offre des fonctionnalités supplémentaires pour manipuler les données.

url = GITLAB_EXAMPLE_GPU + '/' + 'channel_w_archbridge_fully_man004'
print(url)
res = wolfres2DGPU(url) # Initialisation sur base d'un url - pas possible d'utiliser toys_gpu_dataset ici car il renvoit directement un ResultsStore

ERROR:root:HTTP error occurred while downloading https://gitlab.uliege.be/HECE/wolfgpu_examples/-/raw/main/channel_w_archbridge_fully_man004/bridge_deck.npy: 404 Client Error: Not Found for url: https://gitlab.uliege.be/HECE/wolfgpu_examples/-/raw/main/channel_w_archbridge_fully_man004/bridge_deck.npy

https://gitlab.uliege.be/HECE/wolfgpu_examples/-/raw/main/channel_w_archbridge_fully_man004

Récupérer les temps et les pas de sauvegarde

print(res.get_times_steps())

([0.0, 40.60536419227719, 81.21073305234313, 121.81609793752432, 162.4214591383934, 203.02682216465473], [0, 1000, 2000, 3000, 4000, 5000])

Récupérer un pas spécifique

res.read_oneresult(0) # Here we are 0-based, so we read the first result - more Pythonic
res.read_oneresult(-1) # Here we are 0-based, so we read the last result

Accès aux résultats stockés dans l’objet

Pour diverses raisons difficiles à résumer ici, les différentes matrices de résultats sont stockées dans un dictionnaire Python “myblocks”. En GPU, il n’y a qu’une seule entrée au dictionnaire, mais en CPU, il y peut y en avoir plusieurs (aspects multiblocs obligent).

L’objet interne est “OneWolfResult”. Il contient les matrices de base (hauteur d’eau, débit selon X et débit selon Y) mais également des combinaisons (voir “views_2D”).

Le type de résultat fourni par défaut est la hauteur d’eau.

res.myblocks['block1'].waterdepth  # Access the water depth for the first block (GPU only has one block, but CPU can have several blocks)
res.myblocks['block1'].waterdepth.shape  # Shape of the water depth matrix for the first block

(500, 21)

# res.myblocks['block1'].qx # Access the x-component of the unit discharge for the first block
# res.myblocks['block1'].qy # Access the y-component of the unit discharge for the first block
res.myblocks['block1'].current # Access the current view for the first block

<wolfhece.wolf_array.WolfArray at 0x2ab91d57ed0>

res.myblocks['block1'].current is res.myblocks['block1'].waterdepth # By default, the current view is the water depth, but it can be changed to other views (qx, qy, etc.)

True

Il est toutefois possible d’accéder au résultat courant sous le forme d’un WolfArray.

curres_by_index = res[0] # Access the first block result (0-based index)
curres_by_name  = res['block1'] # Access the first block result (1-based index)
type(curres_by_index), type(curres_by_name) # Both should be of type OneWolfResult

print(curres_by_name is curres_by_index) # Both should be the same object
print(np.all(curres_by_index.array == curres_by_name.array)) # But the data should be the same

True
True

Changement de vue

A choisir dans “views_2D”, les vues disponibles sont :

• WATERDEPTH = _(‘Water depth [m]’)

• WATERLEVEL = _(‘Water level [m]’)

• TOPOGRAPHY = _(‘Bottom level [m]’)

• QX = _(‘Discharge X [m2s-1]’)

• QY = _(‘Discharge Y [m2s-1]’)

• QNORM = _(‘Discharge norm [m2s-1]’)

• UX =_(‘Velocity X [ms-1]’)

• UY = _(‘Velocity Y [ms-1]’)

• UNORM = _(‘Velocity norm [ms-1]’)

• HEAD = _(‘Head [m]’)

• FROUDE = _(‘Froude [-]’)

• U_SHEAR = _(‘Shear velocity [ms-1]’)

• SHIELDS_NUMBER = _(‘Shields number - Manning-Strickler’)

• CRITICAL_DIAMETER_SHIELDS = _(‘Critical grain diameter - Shields’)

• CRITICAL_DIAMETER_IZBACH = _(‘Critical grain diameter - Izbach’)

• CRITICAL_DIAMETER_SUSPENSION_50 = _(‘Critical grain diameter - Suspension 50%’)

• CRITICAL_DIAMETER_SUSPENSION_100 = _(‘Critical grain diameter - Suspension 100%’)

res.set_currentview(views_2D.FROUDE)  # Change the current view to the x-component of the unit discharge

res[0].plot_matplotlib()

(<Figure size 640x480 with 1 Axes>, <Axes: >)

Récupérer tous les résultats dans une liste Python

Pour des simulations pas trop lourdes, on peut envisager de récupérer tous les résultats dans une liste Python. Cela permet de manipuler les données plus facilement, mais attention à la mémoire disponible.

res.set_currentview(views_2D.WATERDEPTH)  # Reset the current view to the water depth
all = []
for i in range(res.get_nbresults()):
    res.read_oneresult(i)
    all.append(res.as_WolfArray()) # here using the as_WolfArray method to create a new WolfArray from the current result

print(f"Number of results in the list: {len(all)}")
for i in range(len(all)):
    print('max : ', all[i].array.max(), 'min :', all[i].array.min(), 'median :', np.ma.median(all[i].array), 'mean :', np.ma.mean(all[i].array))

Number of results in the list: 6
max :  8.042151 min : 7.981053 median : 8.007215 mean : 8.017717930545713
max :  8.042157 min : 7.981057 median : 8.007235 mean : 8.017718853354594
max :  8.04216 min : 7.9810605 median : 8.007196 mean : 8.017717930545713
max :  8.042163 min : 7.9810634 median : 8.0072155 mean : 8.017717930545713
max :  8.042169 min : 7.9810696 median : 8.0071945 mean : 8.017717930545713
max :  8.042171 min : 7.981072 median : 8.007233 mean : 8.017717930545713

Création de cartes de danger

Les cartes accessibles sont:

• Danger map - Water depth

• Danger map - Velocity norm

• Danger map - Discharge norm

• Danger map - Water surface elevation

• Danger map - Head

• Danger map - Time of arrival

• Danger map - Time of maximum

dm_h, dm_unorm, dm_qnorm, dm_z, dm_head, dm_toa, dm_tm = res.danger_map() # tuple of WolfArray -  H, U_norm, Q_norm, Z, Head, Time_of_arrival, Time_of_maximum

100%|██████████| 6/6 [00:00<00:00, 103.73it/s]

dm_h.plot_matplotlib()

(<Figure size 640x480 with 1 Axes>, <Axes: >)

Et ensuite ?

La manipulation de résultats et la conversion en WolfArray géoréférencée permet d’effectuer des analyses plus poussées et de visualiser les données de manière plus intuitive.

La combinaison avec des données vectorielles (polygones d’analyse, bâtiments, etc.) permet de croiser les résultats de la simulation avec des données géographiques pour une analyse plus complète.