Zum Hauptinhalt springe

Modellierig vun ere fließende nit-viskose Flüssigkeit mit QUICK-PDE

Hinweis

Qiskit Functions sin e experimentelli Funktion, wo nur für Benutzer vum IBM Quantum® Premium Plan, Flex Plan un On-Prem (via IBM Quantum Platform API) Plan verfügbar isch. Se sin im Preview-Release-Status un kenne sich ändere.

Nutzungsschätzung: 50 Minute uf emm Heron r2-Prozessor. (HINWEIS: Des isch nur e Schätzung. Dini Laufzeit kann variiere.)

Beachte, dass d'Ausführungszeit vun dere Funktion im Allgemeine meh als 20 Minute beträgt, drum ka mr des Tutorial vielleicht in zwoi Abschnitt ufteile: d'erschte, wo mr's durchlest un d'Jobs startet, un d'zwote e paar Stunde später (um d'Jobs gnueg Zeit zur Fertigstellung z'gebe), um mit d'Ergebniss vun d'Jobs z'arbeite.

Hintergrund

Des Tutorial isch dazu do, uf einführendem Niveau z'zeige, wie d'QUICK-PDE-Funktion verwendet wird, um komplexe Multi-Physik-Probleme uf 156Q Heron R2 QPUs under Verwendung vum ColibriTDs H-DES (Hybrid Differential Equation Solver) z'löse. D'Grundlage vum Algorithmus isch im H-DES-Paper beschriebe. Beachte, dass derre Solver au nit-lineare Gleichunge löse kann.

Multi-Physik-Probleme – darunter Strömungsdynamik, Wärmediffusion un Materialverformung, um nur e paar z'nenne – kenne überall durch partielle Differentialgleichunge (PDEs) beschriebe werre.

Selle Probleme sin für verschiedeni Industrie hochrelevant un stelle n wichtige Zweig vun d'angewandti Mathematik dar. D'Lösung vun nit-lineare multivariaten gekoppelte PDEs mit klassische Werkzeig bleibt aber aufgrund vum Bedarf an ere exponentiell große Menge an Ressource schwierig.

Disi Funktion isch für Gleichunge mit zunehmender Komplexität un Variable geeignet un isch d'erschte Schritt, um Möglichkeite z'erschließe, wo emol als unlösbar gälte. Um e durch PDEs modelliertes Problem vollständig z'beschreibe, isch es nötig, d'Anfangs- un Randbedingunge z'kenne. Die kenne d'Lösung vun d'PDE un d'Weg zur Findung ihrer Lösung stark verändere.

Des Tutorial zeigt, wie mr:

  1. D'Parameter vun d'Anfangsbedingungsfunktion definiere.
  2. D'Qubit-Anzahl (verwendet zur Codierung vun d'Funktion vun d'Differentialgleichung), Tiefe un Shot-Anzahl anpasse.
  3. QUICK-PDE ausführe, um d'zugrunde liegende Differentialgleichung z'löse.

Anforderunge

Stell vor em Anfang vun diesem Tutorial sicher, dass des Folgende installiert isch:

  • Qiskit SDK v2.0 oder höher (pip install qiskit)
  • Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
  • Matplotlib (pip install matplotlib)
  • Zugang zur QUICK-PDE-Funktion. Füll des Formular us, um Zugang anzufordere.

Setup

Authentifizier dich mit deim API-Schlüssel un wähl d'Funktion so us:

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

Schritt 1: Eigenschafte vum z'lösende Problem festlege

Des Tutorial behandelt d'Benutzererfahrung us zwoi Perspektive: des physikalische Problem, wo durch d'Anfangsbedingunge bestimmt wird, un d'algorithmische Komponente zur Lösung vun emm Strömungsdynamikbeispiel uf emm Quantencomputer.

Computational Fluid Dynamics (CFD) hät e breits Anwendungsspektrum, un drum isch es wichtig, d'zugrunde liegende PDEs z'studiere un z'löse. E wichtigi Familie vun PDEs sin d'Navier-Stokes-Gleichunge, des sin e System nichtlinearer partieller Differentialgleichunge, wo d'Bewegung vun Flüssigkeite beschreibe. Se sin hochrelevant für wissenschaftliche Probleme un technische Anwendunge.

Under bestimmte Bedingunge reduziere sich d'Navier-Stokes-Gleichunge uf d'Burgers-Gleichung, e Konvektions-Diffusions-Gleichung, wo Phänomene beschreibt, wo in Strömungsdynamik, Gasdynamik un nichtlinearer Akustik uffträte, um nur e paar z'nenne, indem se dissipative Systeme modelliert.

D'eindimensionale Version vun d'Gleichung hängt vun zwoi Variable ab: tR0t \in \mathbb{R}_{\geq 0} modelliert d'zeitliche Dimension, xRx \in \mathbb{R} repräsentiert d'räumliche Dimension. D'allgemeine Form vun d'Gleichung wird d'viskose Burgers-Gleichung gnennt un lautet:

ut+uux=ν2u2x,\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial^2 x},

wobei u(x,t)u(x,t) des Fluidgeschwindigkeitsfeld an ere gegebene Position xx un Zeit tt isch, un ν\nu d'Viskosität vun d'Flüssigkeit isch. Viskosität isch e wichtigi Eigenschaft vun ere Flüssigkeit, wo ihre ratenabhängige Widerstand gege Bewegung oder Verformung misst, un drum spielt se e entscheidende Rolle bi d'Bestimmung vun d'Dynamik vun ere Flüssigkeit. Wenn d'Viskosität vun d'Flüssigkeit null isch (ν\nu = 0), wird d'Gleichung zu ere Erhaltungsgleichung, wo Diskontinuitäte (Stoßwelle) entwickle kann, weil's kein innere Widerstand hät. In dem Fall wird d'Gleichung als inviskose Burgers-Gleichung bezeichnet un isch e Spezialfall vun ere nichtlinearer Wellengleichung.

Streng gnomme träte inviskose Strömunge in d'Natur nit uf, aber bi d'Modellierung aerodynamischer Strömunge kann d'Verwendung vun ere inviskose Beschreibung vum Problem wege dem infinitesimale Transporteffekt nützlich sei. Überraschenderweise befasst sich meh als 70% vun d'aerodynamische Theorie mit inviskose Strömunge.

Des Tutorial verwendet d'inviskose Burgers-Gleichung als CFD-Beispiel zur Lösung uf IBM® QPUs mit QUICK-PDE, gebe durch d'Gleichung:

ut+uux=0\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = 0

D'Anfangsbedingung für des Problem isch uf e lineare Funktion gsetzt: u(t=0,x)=ax+b, mit a,bRu(t=0,x) = ax + b,\text{ mit }a,b\in\mathbb{R} wobei aa un bb beliebige Konstante sin, wo d'Form vun d'Lösung beeinflussen. Mr kenne aa un bb anpasse un luege, wie se d'Lösungsprozess un d'Lösung beeinflussen.

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

Schritt 2 (falls erforderlich): Problem für d'Ausführung uf Quantenhardware optimiere

Standardmäßig verwendet d'Solver physikalisch informierte Parameter, des sin initiale Schaltungsparameter für e gegebeni Qubit-Anzahl un -Tiefe, vun wo unser Solver usgoht.

D'Shots sin au Teil vun d'Parameter mit emm Standardwert, da deren Feinabstimmung wichtig isch.

Je noch d'Konfiguration, wo z'löse versucht wird, müsse d'Parameter vum Algorithmus vielleicht angepasst werre, um zufriedenstellendi Lösunge z'erzielen; zum Beispiel kann es meh oder weniger Qubits pro Variable tt un xx erfordere, je noch aa un bb. Des Folgende passt d'Anzahl vun d'Qubits pro Funktion pro Variable, d'Tiefe pro Funktion un d'Anzahl vun d'Shots an.

Mr kenne au luege, wie mr des Backend un d'Ausführungsmodus spezifiziert.

Darüber hinaus kenne physikalisch informierte Parameter d'Optimierungsprozess in d'falsche Richtung lenke; in dem Fall kenne mr des deaktiviere, indem mr d' initialization-Strategie uf "RANDOM" setzt.

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

Schritt 3: Vergleich vun d'Algorithmusleistunge

Mr kenne d'Konvergenzprozess vun unserer Lösung (HDES) vun job_2 mit d'Leistung vun emm Physics-Informed Neural Networks (PINN)-Algorithmus un -Solver vergleiche (lueg des Paper un des zugehörige GitHub-Repository).

Im Beispiel vun d'Ausgabe vun job_2 (quantenbasierter Ansatz) werre nur 13 Parameter (12 Schaltungsparameter plus 1 Skalierungsparameter) mit dem klassische Solver optimiert. D'Konvergenzprozess verläuft so:

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

Des bedütet, dass e Loss under 0,0015 noch 28 Iteratione erreicht werre kann, un des bi d'Optimierung vun nur wenige klassische Parameter.

Jetzt kenne mir dasselbe mit d'PINN-Lösung mit d'vum Paper vorgeschlagene Standardkonfiguration under Verwendung vun emm gradientenbasierte Optimierer vergleiche. Des Äquivalent vun unserer Schaltung mit 13 z'optimierende Parameter isch des neuronale Netzwerk, wo mindestens acht Schichte mit 20 Neuronen braucht un somit d'Optimierung vun 3021 Parametern beinhaltet. Dann wird d'Ziel-Loss bi Schritt 315, Loss: 0,0014988397, erreicht.

Graph showing PINN data compared with the HDES-Qiskit function.

Jetzt, wo mir en faire Vergleich mache möchte, sollt mir in beide Fäll denselbe Optimierer verwende. D'niedrigste Anzahl vun Iteratione, wo mir für 12 Schichte mit 20 Neuronen = 4701 Parameter gfunde hän:

(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

Mr kenne dasselbe mit de eigene Daten vun job_2 mache un en Vergleich mit d'PINN-Lösung darstelle.

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

Schritt 4: Verwendung vum Ergebnis

Mit der Lösung kenne mr jetzt wähle, was mr damit mache möchte. Des Folgende zeigt, wie mr des Ergebnis darstellt.

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 of the previous code cell

Beachte d'Unterschied vun d'Anfangsbedingung für d'zweite Durchlauf un sini Auswirkung uf des Ergebnis:

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 of the previous code cell

Tutorial-Umfrage

Nimm dir bitte e Minute Zeit, um Feedback zu diesem Tutorial z'gebe. Dini Erkenntnisse helfe uns, unser Inhaltsangebot un unseri Benutzererfahrung z'verbessere:

Link zur Umfrage