Modellierung einer strömenden nicht-viskosen Flüssigkeit mit QUICK-PDE

Hinweis

Qiskit Functions sind eine experimentelle Funktion, die nur für Benutzer des IBM Quantum® Premium Plan, Flex Plan und On-Prem (via IBM Quantum Platform API) Plan verfügbar ist. Sie befinden sich im Preview-Release-Status und können sich ändern.

Nutzungsschätzung: 50 Minuten auf einem Heron r2-Prozessor. (HINWEIS: Dies ist nur eine Schätzung. Deine Laufzeit kann variieren.)

Beachte, dass die Ausführungszeit dieser Funktion im Allgemeinen mehr als 20 Minuten beträgt, daher möchtest du dieses Tutorial möglicherweise in zwei Abschnitte unterteilen: den ersten, in dem du es durchliest und die Jobs startest, und den zweiten einige Stunden später (um den Jobs ausreichend Zeit zur Fertigstellung zu geben), um mit den Ergebnissen der Jobs zu arbeiten.

Hintergrund

Dieses Tutorial soll auf einführendem Niveau vermitteln, wie du die QUICK-PDE-Funktion verwenden kannst, um komplexe Multi-Physik-Probleme auf 156Q Heron R2 QPUs unter Verwendung von ColibriTDs H-DES (Hybrid Differential Equation Solver) zu lösen. Der zugrunde liegende Algorithmus wird im H-DES-Paper beschrieben. Beachte, dass dieser Solver auch nicht-lineare Gleichungen lösen kann.

Multi-Physik-Probleme - einschließlich Strömungsdynamik, Wärmediffusion und Material- deformation, um nur einige zu nennen - können allgegenwärtig durch partielle Differentialgleichungen (PDEs) beschrieben werden.

Solche Probleme sind für verschiedene Industrien hochrelevant und stellen einen wichtigen Zweig der angewandten Mathematik dar. Die Lösung nicht-linearer multivariater gekoppelter PDEs mit klassischen Werkzeugen bleibt jedoch aufgrund der Anforderung einer exponentiell großen Menge an Ressourcen herausfordernd.

Diese Funktion ist für Gleichungen mit zunehmender Komplexität und Variablen geeignet und ist der erste Schritt, um Möglichkeiten zu erschließen, die einst als unlösbar galten. Um ein durch PDEs modelliertes Problem vollständig zu beschreiben, ist es notwendig, die Anfangs- und Randbedingungen zu kennen. Diese können die Lösung der PDE und den Weg zur Findung ihrer Lösung stark verändern.

Dieses Tutorial zeigt dir, wie du:

1. Die Parameter der Anfangsbedingungsfunktion definieren.
2. Die Qubit-Anzahl (verwendet zur Codierung der Funktion der Differentialgleichung), Tiefe und Shot-Anzahl anpassen.
3. QUICK-PDE ausführen, um die zugrunde liegende Differentialgleichung zu lösen.

Anforderungen

Stelle vor Beginn dieses Tutorials sicher, dass du Folgendes installiert hast:

• 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ülle das Formular aus, um Zugang anzufordern.

Einrichtung

Authentifiziere dich mit deinem API-Schlüssel und wähle die Funktion wie folgt aus:

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit-ibm-catalog

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

Dieses Tutorial behandelt die Benutzererfahrung aus zwei Perspektiven: das physikalische Problem, das durch die Anfangsbedingungen bestimmt wird, und die algorithmische Komponente zur Lösung eines Strömungsdynamikbeispiels auf einem Quantencomputer.

Computational Fluid Dynamics (CFD) hat ein breites Anwendungsspektrum, und daher ist es wichtig, die zugrunde liegenden PDEs zu studieren und zu lösen. Eine wichtige Familie von PDEs sind die Navier-Stokes-Gleichungen, die ein System nichtlinearer partieller Differentialgleichungen sind, die die Bewegung von Flüssigkeiten beschreiben. Sie sind hochrelevant für wissenschaftliche Probleme und technische Anwendungen.

Unter bestimmten Bedingungen reduzieren sich die Navier-Stokes-Gleichungen auf die Burgers-Gleichung, eine Konvektions-Diffusions-Gleichung, die Phänomene beschreibt, die in Strömungsdynamik, Gasdynamik und nichtlinearer Akustik auftreten, um nur einige zu nennen, indem sie dissipative Systeme modelliert.

Die eindimensionale Version der Gleichung hängt von zwei Variablen ab: $t \in \mathbb{R}_{\geq 0}$ modelliert die zeitliche Dimension, $x \in \mathbb{R}$ repräsentiert die räumliche Dimension. Die allgemeine Form der Gleichung wird die viskose Burgers-Gleichung genannt und lautet:

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

wobei $u(x,t)$ das Fluidgeschwindigkeitsfeld an einer gegebenen Position $x$ und Zeit $t$ ist, und $\nu$ die Viskosität der Flüssigkeit ist. Viskosität ist eine wichtige Eigenschaft einer Flüssigkeit, die ihren ratenabhängigen Widerstand gegen Bewegung oder Verformung misst, und daher spielt sie eine entscheidende Rolle bei der Bestimmung der Dynamik einer Flüssigkeit. Wenn die Viskosität der Flüssigkeit null ist ($\nu$ = 0), wird die Gleichung zu einer Erhaltungsgleichung, die Diskontinuitäten (Stoßwellen) entwickeln kann, aufgrund des Fehlens ihres inneren Widerstands. In diesem Fall wird die Gleichung als inviskose Burgers-Gleichung bezeichnet und ist ein Spezialfall einer nichtlinearen Wellengleichung.

Streng genommen treten inviskose Strömungen in der Natur nicht auf, aber bei der Modellierung aerodynamischer Strömungen kann die Verwendung einer inviskosen Beschreibung des Problems aufgrund des infinitesimalen Transporteffekts nützlich sein. Überraschenderweise befasst sich mehr als 70% der aerodynamischen Theorie mit inviskosen Strömungen.

Dieses Tutorial verwendet die inviskose Burgers-Gleichung als CFD-Beispiel zur Lösung auf IBM® QPUs mit QUICK-PDE, gegeben durch die Gleichung:

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

Die Anfangsbedingung für dieses Problem ist auf eine lineare Funktion gesetzt: $u(t=0,x) = ax + b,\text{ mit }a,b\in\mathbb{R}$ wobei $a$ und $b$ beliebige Konstanten sind, die die Form der Lösung beeinflussen. Du kannst $a$ und $b$ anpassen und sehen, wie sie den Lösungsprozess und die 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 die Ausführung auf Quantenhardware optimieren

Standardmäßig verwendet der Solver physikalisch informierte Parameter, die initiale Schaltungsparameter für eine gegebene Qubit-Anzahl und -Tiefe sind, von denen unser Solver ausgehen wird.

Die Shots sind ebenfalls Teil der Parameter mit einem Standardwert, da deren Feinabstimmung wichtig ist.

Abhängig von der Konfiguration, die du zu lösen versuchst, müssen die Parameter des Algorithmus möglicherweise angepasst werden, um zufriedenstellende Lösungen zu erzielen; beispielsweise kann es mehr oder weniger Qubits pro Variable $t$ und $x$ erfordern, abhängig von $a$ und $b$. Das Folgende passt die Anzahl der Qubits pro Funktion pro Variable, die Tiefe pro Funktion und die Anzahl der Shots an.

Du kannst auch sehen, wie du das Backend und den Ausführungsmodus spezifizierst.

Darüber hinaus könnten physikalisch informierte Parameter den Optimierungsprozess in eine falsche Richtung lenken; in diesem Fall kannst du dies deaktivieren, indem du die initialization-Strategie auf "RANDOM" setzen.

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 der Algorithmusleistungen

Du kannst den Konvergenzprozess unserer Lösung (HDES) von job_2 mit der Leistung eines Physics-Informed Neural Networks (PINN)-Algorithmus und -Solvers vergleichen (siehe das Paper und das zugehörige GitHub-Repository).

Im Beispiel der Ausgabe von job_2 (quantenbasierter Ansatz) werden nur 13 Parameter (12 Schaltungsparameter plus 1 Skalierungsparameter) mit dem klassischen Solver optimiert. Der Konvergenzprozess verläuft wie folgt:

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 ein Loss unter 0,0015 nach 28 Iterationen erreicht werden kann, und das bei Optimierung nur weniger klassischer Parameter.

Nun können wir dasselbe mit der PINN-Lösung mit der vom Paper vorgeschlagenen Standardkonfiguration unter Verwendung eines gradientenbasierten Optimierers vergleichen. Das Äquivalent unserer Schaltung mit 13 zu optimierenden Parametern ist das neuronale Netzwerk, das mindestens acht Schichten mit 20 Neuronen benötigt und somit die Optimierung von 3021 Parametern beinhaltet. Dann wird der Ziel-Loss bei Schritt 315, Loss: 0,0014988397, erreicht.

Nun, da wir einen fairen Vergleich durchführen möchten, sollten wir in beiden Fällen denselben Optimierer verwenden. Die niedrigste Anzahl von Iterationen, 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 von job_2 tun und einen Vergleich mit der PINN-Lösung darstellen.

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

Schritt 4: Verwendung des Ergebnisses

Mit deiner Lösung kannst du nun wählen, was du damit tun möchtest. Das Folgende demonstriert, wie du das Ergebnis darstellst.

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 der Anfangsbedingung für den zweiten Durchlauf und seine Auswirkung 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 eine Minute Zeit, um Feedback zu diesem Tutorial zu geben. Deine Erkenntnisse helfen uns, unser Inhaltsangebot und unsere Benutzererfahrung zu verbessern:

Link zur Umfrage