Modéliser un fluide non visqueux en écoulement à l'aide de QUICK-PDE

Estimation d'utilisation : 50 minutes sur un processeur Heron r2. (NOTE : Il s'agit uniquement d'une estimation. Votre durée d'exécution peut varier.)

Notez que le temps d'exécution de cette fonction est généralement supérieur à 20 minutes, il est donc conseillé de diviser ce tutoriel en deux parties : la première, au cours de laquelle vous lisez et lancez les travaux, et la seconde, quelques heures plus tard, au cours de laquelle vous lancez les travaux la première, au cours de 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 travailler avec les résultats des tâches.

Arrière-plan

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

Les problèmes multi-physiques - y compris la dynamique des fluides, la diffusion de la chaleur et la déformation des matériaux, pour n'en citer que quelques-uns - peuvent être décrits par des équations aux Par des équations différentielles partielles (EDP).

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

Cette fonction est adaptée aux équations de plus en plus complexes et aux variables variables, et constitue la première étape pour débloquer des possibilités qui étaient auparavant considérées comme irréalisables. Pour décrire complètement un problème modélisé par des EDP, il est nécessaire de connaître les conditions initiales et les conditions aux limites de connaître les conditions initiales et les conditions aux limites. Ceux-ci peuvent fortement modifier la solution de l'EDP et le chemin pour trouver leur solution modifier fortement 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. Ajustez le nombre de qubits (utilisé pour coder la fonction de l'équation différentielle), la profondeur et le nombre de tirs.
3. Exécutez QUICK-PDE pour résoudre l'équation différentielle sous-jacente.

Exigences

Avant de commencer ce tutoriel, assurez-vous que les éléments suivants sont installés :

• Qiskit SDK v2.0 ou plus récent (pip install qiskit)
• Catalogue de fonctions Qiskit (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",
    # For `token`, use the 44-character API_KEY you created
    # and saved from the IBM Quantum Platform Home dashboard
    token="YOUR_API_KEY",
)

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 de l'utilisateur de deux points de vue : le problème physique déterminé par les conditions initiales et l'algorithme problème physique déterminé par les conditions initiales, et la composante algorithmique algorithmique dans la résolution d'un exemple de dynamique des fluides sur un ordinateur quantique.

La dynamique des fluides numérique (CFD) a un large éventail d'applications, et il est donc important d'étudier et de résoudre les EDP sous-jacentes il est donc important d'étudier et de résoudre les EDP sous-jacentes. Une famille importante d'EDP sont les équations de Navier-Stokes, qui sont un système d'équations aux non linéaires décrivant le mouvement des fluides. Ils sont très pertinents pour les problèmes scientifiques et les applications d'ingénierie 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 de Burgers, une équation de convection-diffusion décrivant des phénomènes se produisant dans la dynamique des fluides dynamique des fluides, la dynamique des gaz et l'acoustique non linéaire, pour n'en citer que quelques-uns, en modélisant des systèmes dissipatifs systèmes dissipatifs.

La version unidimensionnelle de l'équation dépend de deux variables : $t \in \mathbb{R}_{\geq 0}$ la modélisation de 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 se lit comme suit :

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

où$u(x,t)$ représente le champ de vitesse du fluide à une position donnée$x$ et à un instant donné$t$ , 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, fonction de la vitesse, et qui, de ce fait, joue 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 pouvant donner lieu à des discontinuités (ondes de choc), en raison de l'absence de résistance interne. Dans ce cas, l'équation est appelée « équation de Burgers sans viscosité » et constitue un cas particulier de l'équation d'onde non linéaire.

À proprement parler, les écoulements inviscides n'existent pas dans la nature l'écoulement aérodynamique, en raison de l'effet infinitésimal du transport, l'utilisation d'une description description inviscide du problème peut s'avérer utile. Il est surprenant de constater que plus de 70% de la théorie aérodynamique traite des écoulements inviscides.

Ce tutoriel utilise l'équation de Burgers sans viscosité comme exemple de simulation CFD à résoudre sur les QPU d' IBM® s à l'aide de QUICK-PDE, telle qu'elle est définie ci-dessous :

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

La condition initiale pour ce problème est définie par 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 qui influencent la forme de la solution. Vous pouvez modifier les paramètres $a$ et $b$ et observer leur influence sur le processus de résolution et la solution obtenue. Pour ce choix de condition initiale, l'équation admet la solution analytique suivante :

$u(t, x) = \frac{ax + b}{at + 1},$

ce qui est vrai tant qu'aucun choc ne s'est formé, c'est-à-dire, par ex $t < -1/a$, lorsque $a < 0$.

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

print(job.result())

Output:

{'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 initiaux du circuit pour un nombre de qubits et une profondeur donnés, à partir desquels notre solveur démarre.

Les tirs font également partie des paramètres avec une valeur par défaut, car il est important de les affiner.

En fonction de la configuration que vous essayez de résoudre, il peut être nécessaire d'adapter les paramètres de l'algorithme pour obtenir des solutions satisfaisantes de l'algorithme pour obtenir des solutions satisfaisantes peuvent devoir être adaptés; par exemple, ils peuvent nécessiter plus ou moins de qubits par variable en fonction de la configuration peut nécessiter plus ou moins de qubits par variable $t$ et $x$, en fonction de $a$ et de $b$. La procédure suivante permet d'ajuster le nombre de qubits par fonction et par variable, la profondeur par fonction et le nombre de qubits par variable 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.

En outre, les paramètres à information physique peuvent orienter le processus d'optimisation dans une mauvaise direction dans une mauvaise direction; dans ce cas, vous pouvez les désactiver en réglant la initialization à "RANDOM".

job_2 = quick.run(
    use_case="CFD_BURGER",
    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())

Output:

{'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) à job_2 la performance d'un algorithme et d'un solveur de réseaux neuronaux informés par la physique (PINN) (voir l 'article et le référentiel GitHub associé).

Dans l'exemple de sortie job_2 de (approche quantique), seuls 13 paramètres (12 paramètres de circuit plus 1 paramètre d'échelle) sont optimisés avec le solveur classique. Le processus de convergence se déroule comme suit :

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 optimisant seulement quelques paramètres classiques.

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

Maintenant, puisque nous voulons faire une comparaison équitable, nous devrions utiliser le même optimiseur dans les deux cas. Le nombre d'itérations le plus faible 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 provenant 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 : Utilisez le résultat

Avec votre solution, vous pouvez maintenant choisir ce qu'il faut en faire. L'exemple suivant 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()

Output:

Remarquez la différence de conditions initiales pour la deuxième série et son effet sur le résultat 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()

Output:

Enquête tutorielle

Veuillez prendre une minute pour nous faire part de vos commentaires sur ce tutoriel. Vos commentaires nous aideront à améliorer nos offres de contenu et l'expérience des utilisateurs :

Lien vers l'enquête