Modéliser un fluide non visqueux en écoulement avec QUICK-PDE

Remarque

Les fonctions Qiskit sont une fonctionnalité expérimentale disponible uniquement pour les utilisateurs des plans IBM Quantum® Premium, Flex et On-Prem (via l'API IBM Quantum Platform). Elles sont en version préliminaire et susceptibles d'être modifiées.

Estimation d'utilisation : 50 minutes sur un processeur Heron r2. (REMARQUE : Il ne s'agit que d'une estimation. Votre temps d'exécution peut varier.)

Notez que le temps d'exécution de cette fonction est généralement supérieur à 20 minutes, vous souhaiterez donc peut-être diviser ce tutoriel en deux sections : la première, dans laquelle vous le lisez et lancez les tâches, et la seconde quelques heures plus tard (en laissant suffisamment de temps pour que les tâches se terminent), pour exploiter les résultats des tâches.

Contexte

Ce tutoriel vise à enseigner, à un niveau introductif, comment utiliser la fonction QUICK-PDE pour résoudre des problèmes multi-physiques complexes sur des QPU Heron R2 à 156Q en utilisant le H-DES (Hybrid Differential Equation Solver) de ColibriTD. L'algorithme sous-jacent est décrit dans l'article sur H-DES. Notez que ce solveur peut également résoudre des équations non linéaires.

Les problèmes multi-physiques - incluant la dynamique des fluides, la diffusion thermique et la déformation des matériaux, pour n'en citer que quelques-uns - peuvent être décrits de manière ubiquitaire par des équations aux dérivées partielles (EDP).

Ces problèmes sont très pertinents pour diverses industries et constituent une branche importante des mathématiques appliquées. Cependant, la résolution d'EDP couplées non linéaires multivariées avec des outils classiques reste un défi en raison de la nécessité d'une quantité exponentiellement importante de ressources.

Cette fonction est appropriée pour des équations de complexité et de variables croissantes, et constitue la première étape pour débloquer des possibilités autrefois considérées comme insolubles. Pour décrire complètement un problème modélisé par des EDP, il est nécessaire de connaître les conditions initiales et aux limites. Celles-ci peuvent fortement modifier la solution de l'EDP et le chemin pour trouver leur solution.

Ce tutoriel vous apprend à :

1. Définir les paramètres de la fonction de condition initiale.
2. Ajuster le nombre de qubits (utilisés pour encoder la fonction de l'équation différentielle), la profondeur et le nombre de tirs.
3. Exécuter QUICK-PDE pour résoudre l'équation différentielle sous-jacente.

Prérequis

Avant de commencer ce tutoriel, assurez-vous d'avoir installé les éléments suivants :

• Qiskit SDK v2.0 ou version ultérieure (pip install qiskit)
• Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
• Matplotlib (pip install matplotlib)
• Accès à la fonction QUICK-PDE. Remplissez le formulaire pour demander l'accès.

Configuration

Authentifiez-vous à l'aide de votre clé API et sélectionnez la fonction comme suit :

import numpy as np
import matplotlib.pyplot as plt
from qiskit_ibm_catalog import QiskitFunctionsCatalog

catalog = QiskitFunctionsCatalog(
    channel="ibm_quantum_platform",
    instance="INSTANCE_CRN",
    token="YOUR_API_KEY",  # Use the 44-character API_KEY you created and saved from the IBM Quantum Platform Home dashboard
)

quick = catalog.load("colibritd/quick-pde")

Étape 1 : Définir les propriétés du problème à résoudre

Ce tutoriel couvre l'expérience utilisateur sous deux angles : le problème physique déterminé par les conditions initiales, et la composante algorithmique pour résoudre un exemple de dynamique des fluides sur un ordinateur quantique.

La dynamique des fluides numérique (CFD) possède un large éventail d'applications, et il est donc important d'étudier et de résoudre les EDP sous-jacentes. Une famille importante d'EDP est constituée par les équations de Navier-Stokes, qui forment un système d'équations aux dérivées partielles non linéaires décrivant le mouvement des fluides. Elles sont très pertinentes pour les problèmes scientifiques et les applications d'ingénierie.

Sous certaines conditions, les équations de Navier-Stokes se réduisent à l'équation de Burgers, une équation de convection-diffusion décrivant des phénomènes survenant en dynamique des fluides, dynamique des gaz et acoustique non linéaire, pour n'en citer que quelques-uns, en modélisant des systèmes dissipatifs.

La version unidimensionnelle de l'équation dépend de deux variables : $t \in \mathbb{R}_{\geq 0}$ modélisant la dimension temporelle, $x \in \mathbb{R}$ représentant la dimension spatiale. La forme générale de l'équation est appelée l'équation de Burgers visqueuse et s'écrit :

$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial^2 x},$

où $u(x,t)$ est le champ de vitesse du fluide à une position $x$ et un temps $t$ donnés, et $\nu$ est la viscosité du fluide. La viscosité est une propriété importante d'un fluide qui mesure sa résistance au mouvement ou à la déformation dépendante du taux, et joue donc un rôle crucial dans la détermination de la dynamique d'un fluide. Lorsque la viscosité du fluide est nulle ($\nu$ = 0), l'équation devient une équation de conservation qui peut développer des discontinuités (ondes de choc), en raison de l'absence de sa résistance interne. Dans ce cas, l'équation est appelée l'équation de Burgers non visqueuse et constitue un cas particulier d'équation d'onde non linéaire.

À strictement parler, les écoulements non visqueux n'existent pas dans la nature, mais lors de la modélisation d'un écoulement aérodynamique, en raison de l'effet infinitésimal du transport, l'utilisation d'une description non visqueuse du problème peut être utile. De manière surprenante, plus de 70 % de la théorie aérodynamique traite des écoulements non visqueux.

Ce tutoriel utilise l'équation de Burgers non visqueuse comme exemple de CFD à résoudre sur les QPU IBM® en utilisant QUICK-PDE, donnée par l'équation :

$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = 0$

La condition initiale de ce problème est définie comme une fonction linéaire : $u(t=0,x) = ax + b,\text{ with }a,b\in\mathbb{R}$ où $a$ et $b$ sont des constantes arbitraires influençant la forme de la solution. Vous pouvez ajuster $a$ et $b$ et observer comment ils influencent le processus de résolution et la solution.

job = quick.run(
    use_case="cfd",
    physical_parameters={"a": 1.0, "b": 1.0},
)

print(job.result())

{'functions': {'u': array([[1.        , 0.96112378, 0.9230742 , 0.88616096, 0.85058445,
        0.81644741, 0.78376878, 0.75249908, 0.72253689, 0.69374562,
        0.66597013, 0.63905258, 0.61284684, 0.58723093, 0.56211691,
        0.53745752, 0.51324915, 0.48953036, 0.46637547, 0.44388257,
        0.4221554 , 0.40127848, 0.38128488, 0.36211604, 0.34357308,
        0.32525895, 0.30651089, 0.28632252, 0.26325504, 0.23533692],
       [1.2375    , 1.19267729, 1.14850734, 1.10544526, 1.06382155,
        1.02385326, 0.98565757, 0.94926734, 0.91464784, 0.88171402,
        0.85034771, 0.82041411, 0.79177677, 0.76431068, 0.73791248,
        0.71250742, 0.68805224, 0.66453346, 0.64196021, 0.62035121,
        0.59971506, 0.5800232 , 0.56117499, 0.54295419, 0.52497612,
        0.50662498, 0.48698059, 0.4647339 , 0.43809065, 0.40466247],
       [1.475     , 1.4242308 , 1.37394048, 1.32472956, 1.27705866,
        1.23125911, 1.18754636, 1.1460356 , 1.10675879, 1.06968242,
        1.03472529, 1.00177563, 0.9707067 , 0.94139043, 0.91370806,
        0.88755732, 0.86285533, 0.83953655, 0.81754494, 0.79681986,
        0.77727473, 0.75876792, 0.74106511, 0.72379234, 0.70637915,
        0.687991  , 0.66745028, 0.64314527, 0.61292625, 0.57398802],
       [1.7125    , 1.65578431, 1.59937362, 1.54401386, 1.49029576,
        1.43866495, 1.38943515, 1.34280386, 1.29886974, 1.25765082,
        1.21910288, 1.18313715, 1.14963664, 1.11847019, 1.08950364,
        1.06260722, 1.03765842, 1.01453964, 0.99312968, 0.97328851,
        0.95483439, 0.93751264, 0.92095522, 0.90463049, 0.88778219,
        0.86935702, 0.84791997, 0.82155665, 0.78776186, 0.74331358],
       [1.95      , 1.88733782, 1.82480676, 1.76329816, 1.70353287,
        1.6460708 , 1.59132394, 1.53957212, 1.49098069, 1.44561922,
        1.40348046, 1.36449867, 1.32856657, 1.29554994, 1.26529921,
        1.23765712, 1.21246152, 1.18954273, 1.16871442, 1.14975716,
        1.13239406, 1.11625736, 1.10084533, 1.08546864, 1.06918523,
        1.05072304, 1.02838966, 0.99996803, 0.96259746, 0.91263913]])}, 'samples': {'t': array([0.        , 0.03275862, 0.06551724, 0.09827586, 0.13103448,
       0.1637931 , 0.19655172, 0.22931034, 0.26206897, 0.29482759,
       0.32758621, 0.36034483, 0.39310345, 0.42586207, 0.45862069,
       0.49137931, 0.52413793, 0.55689655, 0.58965517, 0.62241379,
       0.65517241, 0.68793103, 0.72068966, 0.75344828, 0.7862069 ,
       0.81896552, 0.85172414, 0.88448276, 0.91724138, 0.95      ]), 'x': array([0.    , 0.2375, 0.475 , 0.7125, 0.95  ])}}

Étape 2 (si nécessaire) : Optimiser le problème pour l'exécution sur du matériel quantique

Par défaut, le solveur utilise des paramètres physiquement informés, qui sont des paramètres de circuit initiaux pour un nombre de qubits et une profondeur donnés, à partir desquels notre solveur commencera.

Les tirs font également partie des paramètres avec une valeur par défaut, car leur ajustement fin est important.

Selon la configuration que vous essayez de résoudre, les paramètres de l'algorithme pour obtenir des solutions satisfaisantes peuvent nécessiter une adaptation ; par exemple, cela peut nécessiter plus ou moins de qubits par variable $t$ et $x$, selon $a$ et $b$. Ce qui suit ajuste le nombre de qubits par fonction par variable, la profondeur par fonction et le nombre de tirs.

Vous pouvez également voir comment spécifier le backend et le mode d'exécution.

De plus, les paramètres physiquement informés peuvent orienter le processus d'optimisation dans une mauvaise direction ; dans ce cas, vous pouvez le désactiver en définissant la stratégie d'initialization sur "RANDOM".

job_2 = quick.run(
    use_case="cfd",
    physical_parameters={"a": 0.5, "b": 0.25},
    nb_qubits={"u": {"t": 2, "x": 1}},
    depth={"u": 3},
    shots=[500, 2500, 5000, 10000],
    initialization="RANDOM",
    backend="ibm_kingston",
    mode="session",
)

print(job_2.result())

{'functions': {'u': array([[0.25      , 0.24856543, 0.24687708, 0.2449444 , 0.24277686,
        0.24038389, 0.23777496, 0.23495952, 0.23194702, 0.22874691,
        0.22536866, 0.22182171, 0.21811551, 0.21425952, 0.2102632 ,
        0.20613599, 0.20188736, 0.19752675, 0.19306361, 0.18850741,
        0.18386759, 0.1791536 , 0.17437491, 0.16954096, 0.16466122,
        0.15974512, 0.15480213, 0.1498417 , 0.14487328, 0.13990632],
       [0.36875   , 0.36681313, 0.36457201, 0.36203594, 0.35921422,
        0.35611615, 0.35275103, 0.34912817, 0.34525687, 0.34114643,
        0.33680614, 0.33224532, 0.32747327, 0.32249928, 0.31733266,
        0.31198271, 0.30645873, 0.30077002, 0.29492589, 0.28893564,
        0.28280857, 0.27655397, 0.27018116, 0.26369944, 0.2571181 ,
        0.25044645, 0.24369378, 0.23686941, 0.22998264, 0.22304275],
       [0.4875    , 0.48506084, 0.48226695, 0.47912748, 0.47565158,
        0.47184841, 0.46772711, 0.46329683, 0.45856672, 0.45354594,
        0.44824363, 0.44266894, 0.43683103, 0.43073904, 0.42440212,
        0.41782942, 0.4110301 , 0.4040133 , 0.39678818, 0.38936388,
        0.38174955, 0.37395435, 0.36598742, 0.35785791, 0.34957498,
        0.34114777, 0.33258544, 0.32389713, 0.315092  , 0.30617919],
       [0.60625   , 0.60330854, 0.59996188, 0.59621902, 0.59208895,
        0.58758067, 0.58270318, 0.57746549, 0.57187658, 0.56594545,
        0.55968112, 0.55309256, 0.54618879, 0.53897879, 0.53147158,
        0.52367614, 0.51560147, 0.50725658, 0.49865046, 0.48979211,
        0.48069053, 0.47135472, 0.46179367, 0.45201638, 0.44203186,
        0.4318491 , 0.42147709, 0.41092485, 0.40020136, 0.38931562],
       [0.725     , 0.72155625, 0.71765682, 0.71331056, 0.70852631,
        0.70331293, 0.69767926, 0.69163414, 0.68518643, 0.67834497,
        0.6711186 , 0.66351618, 0.65554655, 0.64721855, 0.63854104,
        0.62952285, 0.62017284, 0.61049986, 0.60051274, 0.59022035,
        0.57963151, 0.56875509, 0.55759992, 0.54617486, 0.53448874,
        0.52255042, 0.51036875, 0.49795257, 0.48531072, 0.47245205]])}, 'samples': {'t': array([0.        , 0.03275862, 0.06551724, 0.09827586, 0.13103448,
       0.1637931 , 0.19655172, 0.22931034, 0.26206897, 0.29482759,
       0.32758621, 0.36034483, 0.39310345, 0.42586207, 0.45862069,
       0.49137931, 0.52413793, 0.55689655, 0.58965517, 0.62241379,
       0.65517241, 0.68793103, 0.72068966, 0.75344828, 0.7862069 ,
       0.81896552, 0.85172414, 0.88448276, 0.91724138, 0.95      ]), 'x': array([0.    , 0.2375, 0.475 , 0.7125, 0.95  ])}}

Étape 3 : Comparer les performances des algorithmes

Vous pouvez comparer le processus de convergence de notre solution (HDES) de job_2 à la performance d'un algorithme et solveur de réseaux de neurones informés par la physique (PINN) (voir l'article et le dépôt GitHub associé).

Dans l'exemple de la sortie de job_2 (approche basée sur le quantique), seuls 13 paramètres (12 paramètres de circuit plus 1 paramètre de mise à l'échelle) sont optimisés avec le solveur classique. Le processus de convergence est le suivant :

optimizers:
   CMA: {'ftarget': np.float64(0.1), 'verb_disp': 10, 'maxiter': 100}
   CMA: {'ftarget': np.float64(0.005), 'verb_disp': 10, 'maxiter': 20}
   CMA: {'ftarget': np.float64(0.0025), 'verb_disp': 10, 'maxiter': 30}
   CMA: {'ftarget': np.float64(0.0005), 'verb_disp': 10, 'maxiter': 10}

500 shots
================== CMA =================
option:  {'ftarget': np.float64(0.1), 'verb_disp': 10, 'maxiter': 100}
0/100, loss: 0.02456641

1000 shots
================== CMA =================
option:  {'ftarget': np.float64(0.005), 'verb_disp': 10, 'maxiter': 20}
0/20, loss: 0.03641833
1/20, loss: 0.02461719
2/20, loss: 0.0283689
3/20, loss: 0.009898383
4/20, loss: 0.04454522
5/20, loss: 0.007019971
6/20, loss: 0.00811147
7/20, loss: 0.01592619
8/20, loss: 0.00764708
9/20, loss: 0.01401516
10/20, loss: 0.01767467
11/20, loss: 0.01220387

5000 shots
================== CMA =================
option:  {'ftarget': np.float64(0.0025), 'verb_disp': 10, 'maxiter': 30}
0/30, loss: 0.01024792
1/30, loss: 0.004343748
2/30, loss: 0.01450951
3/30, loss: 0.008591284
4/30, loss: 0.00266414
5/30, loss: 0.007923613
6/30, loss: 0.02023853
7/30, loss: 0.01031438
8/30, loss: 0.009513116
9/30, loss: 0.008132266
10/30, loss: 0.005787766
11/30, loss: 0.00390582

10000 shots
================== CMA =================
option:  {'ftarget': np.float64(0.0005), 'verb_disp': 10, 'maxiter': 10}
0/10, loss: 0.002386168
1/10, loss: 0.004024823
2/10, loss: 0.001311999
3/10, loss: 0.003433991
4/10, loss: 0.002339664
5/10, loss: 0.002978438
6/10, loss: 0.005458391
7/10, loss: 0.002026701
8/10, loss: 0.00207467
9/10, loss: 0.001947627
final_loss: 0.00151994463476429

Cela signifie qu'une perte inférieure à 0,0015 peut être atteinte après 28 itérations, et en n'optimisant que quelques paramètres classiques.

Nous pouvons maintenant comparer cela à la solution PINN avec la configuration par défaut suggérée par l'article en utilisant un optimiseur basé sur le gradient. L'équivalent de notre circuit avec 13 paramètres à optimiser est le réseau de neurones, qui nécessite au moins huit couches de 20 neurones, et implique donc l'optimisation de 3021 paramètres. Ensuite, la perte cible est atteinte à l'étape 315, perte : 0,0014988397.

Maintenant, puisque nous souhaitons effectuer une comparaison équitable, nous devrions utiliser le même optimiseur dans les deux cas. Le plus faible nombre d'itérations que nous avons trouvé pour 12 couches de 20 neurones = 4701 paramètres :

(10_w,20)-aCMA-ES (mu_w=5.9,w_1=27%) in dimension 4701 (seed=351961)
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1     20 5.398521572351456e-02 1.0e+00 9.98e-03  1e-02  1e-02 0:02.3
    2     40 5.444650724530220e-02 1.0e+00 9.97e-03  1e-02  1e-02 0:05.1
    3     60 4.447407275438309e-02 1.0e+00 9.95e-03  1e-02  1e-02 0:08.2
    4     80 2.068969979882240e-02 1.0e+00 9.94e-03  1e-02  1e-02 0:11.7
    6    120 1.028892211616039e-02 1.0e+00 9.91e-03  1e-02  1e-02 0:20.1
    7    140 5.140972323715687e-03 1.0e+00 9.90e-03  1e-02  1e-02 0:25.4
    9    180 3.811701666563749e-03 1.0e+00 9.87e-03  1e-02  1e-02 0:37.4
   10    200 3.189878538250923e-03 1.0e+00 9.85e-03  1e-02  1e-02 0:44.2
   12    240 2.547040116041899e-03 1.0e+00 9.83e-03  1e-02  1e-02 0:59.7
   14    280 2.166548743844032e-03 1.0e+00 9.80e-03  1e-02  1e-02 1:18.0
   15    300 1.783065614290535e-03 1.0e+00 9.79e-03  1e-02  1e-02 1:28.4
   16    320 2.045844215899706e-03 1.0e+00 9.78e-03  1e-02  1e-02 1:39.8
Stopping early: loss 0.001405 <= target 0.0015
CMA-ES finished. Best loss: 0.001404788694344461

Vous pouvez faire de même avec vos données de job_2, et tracer une comparaison avec la solution PINN.

# check the loss function and compare between the two approaches
print(job_2.logs())

Étape 4 : Utiliser le résultat

Avec votre solution, vous pouvez maintenant choisir ce que vous souhaitez en faire. Ce qui suit montre comment tracer le résultat.

solution = job.result()

# Plot the solution of the second simulation job_2
_ = plt.figure()
ax = plt.axes(projection="3d")

# plot the solution using the 3d plotting capabilities of pyplot
t, x = np.meshgrid(solution["samples"]["t"], solution["samples"]["x"])
ax.plot_surface(
    t,
    x,
    solution["functions"]["u"],
    edgecolor="royalblue",
    lw=0.25,
    rstride=26,
    cstride=26,
    alpha=0.3,
)
ax.scatter(t, x, solution, marker=".")
ax.set(xlabel="t", ylabel="x", zlabel="u(t,x)")

plt.show()

Remarquez la différence de condition initiale pour la seconde exécution et son effet sur le résultat :

solution_2 = job_2.result()

# Plot the solution of the second simulation job_2
_ = plt.figure()
ax = plt.axes(projection="3d")

# plot the solution using the 3d plotting capabilities of pyplot
t, x = np.meshgrid(solution_2["samples"]["t"], solution_2["samples"]["x"])
ax.plot_surface(
    t,
    x,
    solution_2["functions"]["u"],
    edgecolor="royalblue",
    lw=0.25,
    rstride=26,
    cstride=26,
    alpha=0.3,
)
ax.scatter(t, x, solution_2, marker=".")
ax.set(xlabel="t", ylabel="x", zlabel="u(t,x)")

plt.show()

Enquête sur le tutoriel

Veuillez prendre une minute pour donner votre avis sur ce tutoriel. Vos retours nous aideront à améliorer nos contenus et l'expérience utilisateur :

Lien vers l'enquête