Modellierung von aner strömenden nicht-viskosen Flüssigkeit mit QUICK-PDE
Qiskit Functions san a experimentelles Feature, das nur für Nutzerinnen und Nutzer von IBM Quantum® Premium Plan, Flex Plan und On-Prem (über IBM Quantum Platform API) Plan zugänglich is. Es befind sich im Preview-Release-Status und kann sich noch ändern.
Geschätzter Aufwand: 50 Minuten auf'm Heron-r2-Prozessor. (HINWEIS: Das is nur a Schätzung. Deine tatsächliche Laufzeit kann abweichen.)
Beachte, dass die Ausführungszeit dieser Funktion im Allgemeinen über 20 Minuten liegt, also willst du dieses Tutorial vielleicht in zwei Abschnitte aufteilen: zuerst liest du's durch und startest die Jobs, und a paar Stunden später (um den Jobs genug Zeit zum Fertigwerden zu lassen) arbeitest du dann mit den Ergebnissen.
Hintergrund
In diesem Tutorial lernst du auf Einsteigerniveau, wie du die QUICK-PDE-Funktion verwendest, um komplexe Multiphysik-Probleme auf 156Q-Heron-R2-QPUs zu lösen — mithilfe von ColibriTDs H-DES (Hybrid Differential Equation Solver). Der zugrundeliegende Algorithmus is im H-DES-Paper beschrieben. Beachte, dass dieser Solver auch nichtlineare Gleichungen lösen kann.
Multiphysik-Probleme — darunter Fluiddynamik, Wärmediffusion und Materialverformung, um nur ein paar zu nennen — lassen sich allgemein durch partielle Differentialgleichungen (PDEs) beschreiben.
Solche Probleme san in vielen Branchen hoch relevant und bilden a wichtigen Zweig der angewandten Mathematik. Das Lösen von nichtlinearen multivariaten gekoppelten PDEs mit klassischen Werkzeugen is allerdings nach wie vor a Herausforderung, weil dafür exponentiell viele Ressourcen benötigt werden.
Diese Funktion eignet sich für Gleichungen mit zunehmender Komplexität und mehr Variablen und is a erster Schritt, um Möglichkeiten zu erschließen, die einmal als unlösbar galten. Um a Problem, das durch PDEs modelliert wird, vollständig zu beschreiben, muss man die Anfangs- und Randbedingungen kennen. Die können die Lösung der PDE und den Weg dorthin stark beeinflussen.
In diesem Tutorial lernst du:
- Wie du die Parameter der Anfangsbedingungsfunktion definierst.
- Wie du die Qubit-Anzahl (zum Kodieren der Differentialgleichungsfunktion), die Tiefe und die Shot-Anzahl anpasst.
- Wie du QUICK-PDE ausführst, um die zugrundeliegende Differentialgleichung zu lösen.
Voraussetzungen
Bevor du mit dem Tutorial anfängst, stell sicher, dass du Folgendes installiert hast:
- Qiskit SDK v2.0 oder neuer (
pip install qiskit) - Qiskit Functions Catalog (
pip install qiskit-ibm-catalog) - Matplotlib (
pip install matplotlib) - Zugang zur QUICK-PDE-Funktion. Füll das Formular für den Zugangsantrag aus.
Einrichtung
Authentifizier dich mit deinem API-Schlüssel und wähl die Funktion wie folgt aus:
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: Eigenschaften des zu lösenden Problems festlegen
In diesem Tutorial wird die Nutzererfahrung aus zwei Perspektiven beleuchtet: dem physikalischen Problem, das durch die Anfangsbedingungen bestimmt wird, und der algorithmischen Komponente beim Lösen eines Fluiddynamik-Beispiels auf'm Quantencomputer.
Computational Fluid Dynamics (CFD) hat a breites Anwendungsspektrum, daher is es wichtig, die zugrundeliegenden PDEs zu untersuchen und zu lösen. Eine bedeutende Klasse von PDEs san die Navier-Stokes-Gleichungen — a System nichtlinearer partieller Differentialgleichungen, die die Bewegung von Flüssigkeiten beschreiben. Sie san für wissenschaftliche Probleme und technische Anwendungen von großer Bedeutung.
Unter bestimmten Bedingungen reduzieren sich die Navier-Stokes-Gleichungen auf die Burgers-Gleichung, a Konvektions-Diffusions-Gleichung, die unter anderem Phänomene in der Fluiddynamik, der Gasdynamik und der nichtlinearen Akustik beschreibt, indem sie dissipative Systeme modelliert.
Die eindimensionale Version der Gleichung hängt von zwei Variablen ab: modelliert die Zeitdimension, steht für die Raumdimension. Die allgemeine Form der Gleichung heißt viskose Burgers-Gleichung und lautet:
wobei das Strömungsgeschwindigkeitsfeld bei einer bestimmten Position und Zeit is, und die Viskosität der Flüssigkeit darstellt. Viskosität is a wichtige Eigenschaft einer Flüssigkeit, die ihren geschwindigkeitsabhängigen Widerstand gegen Bewegung oder Verformung misst, und spielt daher a entscheidende Rolle bei der Bestimmung der Dynamik einer Flüssigkeit. Wenn die Viskosität null is ( = 0), wird die Gleichung zu a Erhaltungsgleichung, die Unstetigkeiten (Stoßwellen) entwickeln kann, weil der innere Widerstand fehlt. In diesem Fall nennt man die Gleichung die inviskide Burgers-Gleichung, und sie is a Sonderfall der nichtlinearen Wellengleichung.
Streng genommen kommen inviskide Strömungen in der Natur nicht vor, aber bei der Modellierung von aerodynamischen Strömungen kann eine inviskide Beschreibung des Problems nützlich sein — wegen des infinitesimalen Transporteffekts. Überraschenderweise befasst sich mehr als 70 % der Aerodynamiktheorie mit inviskiden Strömungen.
Dieses Tutorial verwendet die inviskide Burgers-Gleichung als CFD-Beispiel, das auf IBM®-QPUs mit QUICK-PDE gelöst wird, gegeben durch die Gleichung:
Die Anfangsbedingung für dieses Problem is auf eine lineare Funktion gesetzt: wobei und beliebige Konstanten san, die die Form der Lösung beeinflussen. Du kannst und anpassen und schauen, wie sie den Lösungsprozess und das Ergebnis 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 nötig): Problem für die Ausführung auf Quantenhardware optimieren
Standardmäßig verwendet der Solver physikalisch informierte Parameter — das san anfängliche Schaltkreisparameter für eine bestimmte Qubit-Anzahl und Tiefe, von der unser Solver startet.
Die Shots san ebenfalls Teil der Parameter mit a Standardwert, da ihre Feinabstimmung wichtig is.
Abhängig von der Konfiguration, die du lösen willst, müssen die Algorithmusparameter für zufriedenstellende Lösungen möglicherweise angepasst werden; zum Beispiel kann's je nach und mehr oder weniger Qubits pro Variable und brauchen. Das Folgende passt die Anzahl der Qubits pro Funktion pro Variable, die Tiefe pro Funktion und die Shot-Anzahl an.
Du kannst auch schauen, wie du das Backend und den Ausführungsmodus festlegst.
Außerdem können physikalisch informierte Parameter den Optimierungsprozess in die falsche Richtung lenken; in dem Fall kannst du das deaktivieren, indem du die initialization-Strategie auf "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: Algorithmusleistung vergleichen
Du kannst den Konvergenzprozess unserer Lösung (HDES) von job_2 mit der Leistung eines physikalisch informierten neuronalen Netzes (PINN) und dessen Solver vergleichen (sieh das Paper und das dazugehörige GitHub-Repository).
Im Beispiel der Ausgabe von job_2 (quantenbasierter Ansatz) werden nur 13 Parameter (12 Schaltkreisparameter plus 1 Skalierungsparameter) mit dem klassischen Solver optimiert. Der Konvergenzprozess is der Folgende:
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
Das bedeutet, dass nach 28 Iterationen und mit der Optimierung von nur a paar klassischen Parametern a Loss unter 0,0015 erreicht werden kann.
Jetzt können wir das mit der PINN-Lösung mit der vom Paper vorgeschlagenen Standardkonfiguration und einem gradientenbasierten Optimizer vergleichen. Das Äquivalent zu unserem Schaltkreis mit 13 zu optimierenden Parametern is das neuronale Netz, das mindestens acht Schichten mit je 20 Neuronen braucht und somit die Optimierung von 3021 Parametern umfasst. Dabei wird der Ziel-Loss bei Schritt 315, Loss: 0,0014988397 erreicht.
Da wir a fairen Vergleich anstellen wollen, solln wir in beiden Fällen denselben Optimizer verwenden. Die niedrigste Iterationszahl, die wir für 12 Schichten mit 20 Neuronen = 4701 Parameter gefunden haben:
(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
Du kannst dasselbe mit deinen Daten aus job_2 machen und a Vergleich mit der PINN-Lösung plotten.
# check the loss function and compare between the two approaches
print(job_2.logs())
Schritt 4: Das Ergebnis verwenden
Mit deiner Lösung kannst du jetzt entscheiden, was du damit machen willst. Das Folgende zeigt, wie du das Ergebnis plottest.
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()
Beachte den Unterschied in der Anfangsbedingung für den zweiten Durchlauf und seinen Einfluss auf das 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()
Tutorial-Umfrage
Bitte nimm dir a Minüterl Zeit und gib uns Feedback zu diesem Tutorial. Deine Rückmeldung hilft uns, unsere Inhalte und die Nutzererfahrung zu verbessern: