To't Hauptinholt springen

Een fließende nich-viskose Flüssigkeit mit QUICK-PDE modelleren

Henwies

Qiskit Functions sünd een experimentelle Funktschoon, de bloß för Brukers vun IBM Quantum® Premium Plan, Flex Plan un On-Prem (över IBM Quantum Platform API) Plan verfögbor is. Se sünd in'n Preview-Release-Status un künnt sik ännern.

Bruukschatzung: 50 Minuten op een Heron r2-Prozessor. (HENWIES: Dat is bloß een Schatzung. Dien Looptiet kann variern.)

Merk op, dat de Utföhrtiet vun disse Funktschoon allgemeen mehr as 20 Minuten duert, dorüm wullt du dit Tutorial villicht in twee Afsneden delen: dat eerste, wo du dat dörchleest un de Jobs anfängst, un dat tweete een poor Stünnen later (um de Jobs gnoog Tiet to'n Fardigmaken to geven), um mit de Resultaten vun de Jobs to warken.

Achtergrund

Dit Tutorial is dorop ut, op een inföhren Niveau to wiesen, woans du de QUICK-PDE-Funktschoon bruken deist, um komplexe Multi-Physik-Problemen op 156Q Heron R2 QPUs mit ColibriTDs H-DES (Hybrid Differential Equation Solver) to lösen. De togrund liggende Algorithmus is in dat H-DES-Paper beschreven. Merk op, dat disse Solver ook nich-lineare Gliekungen lösen kann.

Multi-Physik-Problemen – darunner Strömungsdynamik, Wärmtediffuschoon un Material- deformatschoon, um bloß een poor to nömen – künnt allerwegen dörch partielle Differentialgliekungen (PDEs) beschreven warrn.

Söke Problemen sünd för verscheden Industrien hoochrelevant un stellt een wichtigen Tweig vun de anwennte Mathematik dar. Dat Lösen vun nich-lineare multivariiate koppelte PDEs mit klassische Warktüüg blievt swoor wegen de Anforderung vun een exponentiell groot Menge an Ressourcen.

Disse Funktschoon is för Gliekungen mit tonehmende Komplexität un Variabelen passt un is de eerste Schritt, um Möglichkeiten to openen, de fröher as unlösbor gullen. Um een dörch PDEs modelleerd Problem vollständig to beschrieven, is dat nödig, de Anfangs- un Randbedingungen to kennen. De künnt de Lösung vun de PDE un den Weg to't Finnen vun ehr Lösung stark verännnern.

Dit Tutorial wiest di, woans du:

  1. De Parameter vun de Anfangsbedingungsfunktschoon definieerst.
  2. De Qubit-Antall (bruukt för dat Kodieren vun de Funktschoon vun de Differentialgliekung), Deepde un Shot-Antall anpasst.
  3. QUICK-PDE utföhrst, um de togrund liggende Differentialgliekung to lösen.

Anforderungen

Sörg dorför, dat du, ehr du mit dit Tutorial anfängst, dat Folgende installeert hest:

  • Qiskit SDK v2.0 oder höger (pip install qiskit)
  • Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
  • Matplotlib (pip install matplotlib)
  • Togang to de QUICK-PDE-Funktschoon. Vull dat Formular ut, um Togang to anfrodern.

Setup

Authentifizeer di mit dien API-Slötel un wähl de Funktschoon so ut:

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: Egenschapen vun dat to lösende Problem fastleggen

Dit Tutorial behannelt de Brukerörvaring ut twee Perspektiven: dat physikalische Problem, dat dörch de Anfangsbedingungen bestimmt is, un de algorithmische Komponente för dat Lösen vun een Strömungsdynamikbispill op een Quantenrekner.

Computational Fluid Dynamics (CFD) hett een breed Anwendungsspektrum, un dorüm is dat wichtig, de togrund liggende PDEs to studeern un to lösen. Een wichtige Familie vun PDEs sünd de Navier-Stokes-Gliekungen, dat is een System vun nich-lineare partielle Differentialgliekungen, de de Bewegung vun Flüssigkeiten beschrieven. Se sünd hoochrelevant för wetenschaftliche Problemen un technische Anwendungen.

Unner bestimmte Bedingungen reduzeren sik de Navier-Stokes-Gliekungen op de Burgers-Gliekung, een Konvektschoons-Diffuschoons-Gliekung, de Phänomene beschrievt, de in Strömungsdynamik, Gasdynamik un nich-lineare Akustik optreden, um bloß een poor to nömen, indeem se dissipative Systeme modelleert.

De eendimensionale Verschoon vun de Gliekung hängt vun twee Variabelen af: tR0t \in \mathbb{R}_{\geq 0} modelleert de tietliche Dimenschoon, xRx \in \mathbb{R} represendeert de rümliche Dimenschoon. De allgemene Form vun de Gliekung warrt de viskose Burgers-Gliekung nömmt un luudt:

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) dat Fluidgeswindichkeitsfeld an een geven Positschoon xx un Tiet tt is, un ν\nu de Viskosität vun de Flüssigkeit is. Viskosität is een wichtige Egenschap vun een Flüssigkeit, de ehr ratenafhanging Wedderstann gegen Bewegung oder Verformung mett, un dorüm speelt se een entscheidende Rull bi de Bestimmen vun de Dynamik vun een Flüssigkeit. Wenn de Viskosität vun de Flüssigkeit null is (ν\nu = 0), warrt de Gliekung to een Erholtungsgliekung, de Diskontinuitäten (Stotwüllen) ontwikkeln kann, wegen dat Fehlen vun ehr innere Wedderstann. In dat Fall warrt de Gliekung as inviskose Burgers-Gliekung betekent un is een Spezialfall vun een nich-lineare Wüllengliekung.

Streng nohmen tredt inviskose Strömungen in de Natur nich op, aver bi dat Modelleren vun aerodynamische Strömungen kann dat Bruken vun een inviskose Beschrievung vun dat Problem wegen den infinitesimalen Transporteffekt nützlich ween. Överwunnerlick befatt sik mehr as 70% vun de aerodynamische Theorie mit inviskose Strömungen.

Dit Tutorial bruukt de inviskose Burgers-Gliekung as CFD-Bispill to't Lösen op IBM® QPUs mit QUICK-PDE, geven dörch de Gliekung:

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

De Anfangsbedingung för dit Problem is op een lineare Funktschoon sett: 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 Konstanten sünd, de de Form vun de Lösung beinflusst. Du kannst aa un bb anpassen un sehn, woans se den Lösungsprozess un de Lösung beinflusst.

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 (wenn nödig): Problem för de Utföhren op Quantenhardware optimeren

Standardmäßig bruukt de Solver physikalisch informeerde Parameter, dat sünd initiale Schalltungsparameter för een geven Qubit-Antall un -Deepde, vun de unse Solver utgahn warrt.

De Shots sünd ook Deel vun de Parameter mit een Standardweert, do ehr Fienafstimmen wichtig is.

Afhanging vun de Konfiguratschoon, de du to lösen versökst, möten de Parameter vun den Algorithmus villicht anpasst warrn, um tofredenstellende Lösungen to kriegen; to'n Bispill kann dat mehr oder winniger Qubits pro Variable tt un xx verlangt, afhanging vun aa un bb. Dat Folgende passt de Antall vun de Qubits pro Funktschoon pro Variable, de Deepde pro Funktschoon un de Antall vun de Shots an.

Du kannst ook sehn, woans du dat Backend un den Utföhrmodus spezifizerst.

Daröver rut künnt physikalisch informeerde Parameter den Optimerungsprozess in de verkehrte Richtung lenken; in dat Fall kannst du dat deaktiveren, indeem du de initialization-Strategie op "RANDOM" settst.

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: Algorithmusleistungen verglieken

Du kannst den Konvergenzprozess vun unse Lösung (HDES) vun job_2 mit de Leistung vun een Physics-Informed Neural Networks (PINN)-Algorithmus un -Solver verglieken (kiek dat Paper un dat togehörige GitHub-Repository).

In dat Bispill vun de Utgaav vun job_2 (quantenbaseerde Ansatz) warrt bloß 13 Parameter (12 Schalltungsparameter plus 1 Skalerungsparameter) mit den klassischen Solver optimeert. De Konvergenzprozess verlöppt 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

Dat bedüüdt, dat een Loss unner 0,0015 na 28 Iteratschonen erreicht warrn kann, un dat bi dat Optimeren vun bloß winnig klassische Parameter.

Nu künnt wi datsülvige mit de PINN-Lösung mit de vun dat Paper vörslahn Standardkonfiguratschoon unner Bruuk vun een gradientenbaseerden Optimerer verglieken. Dat Äquivalent vun unse Schalltung mit 13 to optimerende Parameter is dat neuronale Nettwark, dat mindestens acht Schichten mit 20 Neuronen bruukt un so de Optimerung vun 3021 Parameter inslött. Denn warrt de Ziil-Loss bi Schritt 315, Loss: 0,0014988397, erreicht.

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

Nu, wiel wi een fairen Vergliek dörchdragen wüllt, schullt wi in beid Fäll densülvigen Optimerer bruken. De neddrigste Antall vun Iteratschonen, de wi för 12 Schichten mit 20 Neuronen = 4701 Parameter funnen hebbt:

(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 datsülvige mit dien Daten vun job_2 maken un een Vergliek mit de PINN-Lösung darstellen.

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

Schritt 4: Bruuk vun dat Resultat

Mit dien Lösung kannst du nu wählen, wat du dormit maken wullt. Dat Folgende wiest, woans du dat Resultat darstellen deist.

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 op den Unnerscheed vun de Anfangsbedingung för den tweeten Dörchloop un sien Utwarken op dat Resultat:

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-Ümfraag

Nimm di bitte een Minutt Tiet, um Feedback to dit Tutorial to geven. Dien Inzichten helpt uns, unse Inholtsbott un unse Brukerörvaring to verbetern:

Link to de Ümfraag