Zum Hauptinhalt springe

E fliessendi nöd-viskosi Flüssigkeit mit QUICK-PDE modelliere

Hiiwis

Qiskit Functions sind e experimentelli Funktion, wo nume für Benutzer vom IBM Quantum® Premium Plan, Flex Plan und On-Prem (über IBM Quantum Platform API) Plan verfüegbar isch. Si befinded sich im Preview-Release-Status und chönd sich ändere.

Nutzigsschätzig: 50 Minute uf emene Heron r2-Prozässor. (HIIWIS: Das isch nume ne Schätzig. Dini Laufziit chan variiere.)

Merk, dass d'Usführigsziit vo dere Funktion im Allgemeine meh als 20 Minute betreit, drum wottsch das Tutorial vilicht i zwöi Abschnittiile: dr erscht, wo du 's dureläsisch und d'Jobs startisch, und dr zweit es paar Stund spöter (zum de Jobs gnueg Ziit zur Fertigstellig z'gä), zum mit de Ergäbnis vo de Jobs z'schaffe.

Hingergrund

Das Tutorial will uf ere iführende Stufe zeige, wie du d'QUICK-PDE-Funktion bruuchsch, zum komplexe Multi-Physik-Problem uf 156Q Heron R2 QPUs mit em ColibriTDs H-DES (Hybrid Differential Equation Solver) z'löse. Dr zugrundliegend Algorithmus wird im H-DES-Paper beschribe. Merk, dass dä Solver au nöd-lineari Gliichige löse chan.

Multi-Physik-Problem – drunte Strömigsdynamik, Wärmediffusion und Material- deformation, zum nume es paar z'nenne – chönd allewäg dur partielli Differenzialgliichige (PDEs) beschribe werde.

Sötigi Problem sind für verschideni Industrie hoochwichtig und stelled e wichtige Zwiig vo de agwandte Mathematik dar. S'Löse vo nöd-lineare multivariable gchopplete PDEs mit klassische Wärchzüüg bliibt schwiirig wäge de Aaforderig vo nere exponentiell grosse Mängi a Ressurce.

Dä Funktion isch für Gliichige mit zuenähmender Komplexität und Variable geignet und isch dr erscht Schritt, zum Möglichkeite z'erschliesse, wo emal als unlösbar ggulte händ. Zum es dur PDEs modelliertes Problem vollständig z'beschriibe, isch es nötig, d'Aafangs- und Randbedingige z'chenne. Die chönd d' Lösig vo de PDE und de Wäg zur Findig vo ihrere Lösig staarch verändere.

Das Tutorial zeigt dir, wie du:

  1. D'Parameter vo de Aafangsbedigigsfunktion definiersch.
  2. D'Qubit-Aazaal (bruucht zur Codierig vo de Funktion vo de Differenzialgliichig), Tüüfi und Shot-Aazaal apässisch.
  3. QUICK-PDE usfüehrsch, zum d'zugrundliegendi Differenzialgliichig z'löse.

Aaforderige

Stell sicher, bevor du mit däm Tutorial aafangsch, dass du Folgendes installiert häsch:

  • Qiskit SDK v2.0 oder höcher (pip install qiskit)
  • Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
  • Matplotlib (pip install matplotlib)
  • Zugang zur QUICK-PDE-Funktion. Füll das Formular us, zum Zugang z'afördere.

Setup

Authentifizier dich mit dim API-Schlüssel und 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: Eigeschafte vom z'lösende Problem feschtlege

Das Tutorial behandlet d'Benutzererfahrig us zwöi Perspektive: s' physikalischi Problem, wo dur d'Aafangsbedingige bestimmt wird, und d'algorithmischi Komponänte zum Löse vo nemene Strömigsdynamikbiispil uf emene Quanterechner.

Computational Fluid Dynamics (CFD) het es breits Aawändigsspektrum, und drum isch es wichtig, d'zugrundliegende PDEs z'studiere und z'löse. E wichtigi Familie vo PDEs sind d'Navier-Stokes-Gliichige, das isch es System vo nöd-lineare partielli Differenzialgliichige, wo d'Bewegig vo Flüssigkeite beschribed. Si sind hoochwichtig für wüsseschaftlichi Problem und technischi Aawändige.

Unter bestimmte Bedingige reduzierend sich d'Navier-Stokes-Gliichige uf d'Burgers-Gliichig, e Konvektions-Diffusions-Gliichig, wo Phänomen beschribt, wo i Strömigsdynamik, Gasdynamik und nöd-lineare Akustik uuftrete, zum nume es paar z'nenne, indem si dissipativ System modelliert.

D'eidimensionali Version vo de Gliichig hängt vo zwöi Variable ab: tR0t \in \mathbb{R}_{\geq 0} modelliert d'ziitlichi Dimension, xRx \in \mathbb{R} repräsentiert d'rüümlichi Dimension. D'allgemein Form vo de Gliichig wird d'viskosi Burgers-Gliichig gnännt und luutet:

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

wobii u(x,t)u(x,t) s'Fluidgschwündigkeitsfeld a nere ggäbnige Position xx und Ziit tt isch, und ν\nu d'Viskosität vo de Flüssigkeit isch. Viskosität isch e wichtigi Eigeschaft vo nere Flüssigkeit, wo ihre rateabhängig Widerstand gege Bewegig oder Verformig misst, und drum spilt si e entschiidendi Rolle bi de Bestimmig vo de Dynamik vo nere Flüssigkeit. Wenn d' Viskosität vo de Flüssigkeit null isch (ν\nu = 0), wird d'Gliichig zun ere Erhaltigsgliichig, wo Diskontinuitäte (Stosswälle) entwickle chan, wäge däm Fähle vo ihrem innere Widerstand. I däm Fall wird d'Gliichig als inviskosi Burgers-Gliichig bezeichnet und isch e Spezialfall vo nere nöd-lineare Wällegliichig.

Sträng gnoh treted inviskosi Strömige i de Natur nöd uuf, aber bi de Modellierig vo aerodynamische Strömige chan d'Verwendig vo nere inviskose Beschriibig vom Problem wäge däm infinitesimale Transporteffäkt nützlich sii. Überraschenderwis befasst sich meh als 70% vo de aerodynamische Theorie mit inviskose Strömige.

Das Tutorial bruucht d'inviskosi Burgers-Gliichig als CFD-Biispil zur Lösig uf IBM® QPUs mit QUICK-PDE, ggä dur d'Gliichig:

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

D'Aafangsbedingig für das Problem isch uf e lineari Funktion gsetzt: u(t=0,x)=ax+b, mit a,bRu(t=0,x) = ax + b,\text{ mit }a,b\in\mathbb{R} wobii aa und bb beliebigi Konstante sind, wo d'Form vo de Lösig beiflussed. Du chasch aa und bb apasse und luege, wie si de Lösigsprozäss und d'Lösig beiflussed.

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 d'Usfüehrig uf Quantehardware optimiere

Standardmässig bruucht dr Solver physikalisch informierti Parameter, das sind initiali Schaltigsparameter für e ggäbnigi Qubit-Aazaal und -Tüüfi, vo däne use Solver usgoh wird.

D'Shots sind au Tail vo de Parameter mit emene Standardwärt, well ihre Fiiaabstimmig wichtig isch.

Abhängig vo de Konfiguration, wo du z'löse versuesch, müend d'Parameter vom Algorithmus vilicht agpasst werde, zum zufridestellendi Lösige z'erziille; zum Biispil chan's meh oder wäniger Qubits pro Variabli tt und xx verlange, abhängig vo aa und bb. S'Folgende passt d'Aazaal vo de Qubits pro Funktion pro Variabli, d'Tüüfi pro Funktion und d'Aazaal vo de Shots a.

Du chasch au luege, wie du s'Backend und de Usfüehrigsmodus spezifiziersch.

Drüber use chönd physikalisch informierti Parameter de Optimierigsprozäss i d'falschi Richtig leinke; i däm Fall chasch das deaktiviere, indem du d' initialization-Strategie uf "RANDOM" setzisch.

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: Algorithmusleistige vergliiche

Du chasch de Konvergenzprozäss vo eusere Lösig (HDES) vo job_2 mit de Leistig vo nemene Physics-Informed Neural Networks (PINN)-Algorithmus und -Solver vergliiche (lueg s'Paper und s'zuegherig GitHub-Repository).

Im Biispil vo de Usgab vo job_2 (quantebasiert Aasatz) werded nume 13 Parameter (12 Schaltigsparameter plus 1 Skaleringsparameter) mit em klassische Solver optimiert. Dr Konvergenzprozäss verlauft 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

Das bedütet, dass e Loss under 0,0015 nach 28 Iteratione erreicht werde chan, und das bim Optimiere vo nume wänige klassische Parameter.

Jetzt chönd mir s'gliichi mit de PINN-Lösig mit de vom Paper vorgschlagene Standardkonfiguration under Verwendig vo nemene gradientebasierte Optimierer vergliiche. S'Äquivalent vo eusere Schaltig mit 13 z'optimierende Parameter isch s'neuronali Netzwärch, wo mindeschtens acht Schichte mit 20 Neurone bruucht und somit d'Optimierig vo 3021 Parameter iischliesst. Denn wird dr Ziil-Loss bi Schritt 315, Loss: 0,0014988397, erreicht.

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

Jetzt, will mir e faire Vergliich dureführe wänd, söted mir i beide Fäll de gliich Optimierer verwende. D'nidrigscht Aazaal vo Iteratione, wo mir für 12 Schichte mit 20 Neurone = 4701 Parameter gfunde händ:

(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 chasch s'gliichi mit dine Date vo job_2 mache und e Vergliich mit de PINN-Lösig plotte.

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

Schritt 4: Verwendig vom Ergäbnis

Mit dinere Lösig chasch jetzt wähle, was du dermit wottsch mache. S'Folgende zeigt, wie du s'Ergäbnis plottesch.

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

Merk de Underschid vo de Aafangsbedingig für de zweit Durchlauf und siini Uuswirkig uf s'Ergäbnis:

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-Umfrag

Nimm dir bitte e Minutt Ziit, zum Feedback zu däm Tutorial z'gä. Diini Iiblick hälfed eus, eusi Inhaltsbott und eusi Benutzererfahrig z'verbessere:

Link zur Umfrag