Modellieren Sie eine fließende, nicht viskose Flüssigkeit mit QUICK-PDE

Geschätzte Nutzungsdauer: 50 Minuten auf einem Heron r2 Prozessor. (HINWEIS: Dies ist nur eine Schätzung. Ihre Laufzeit kann variieren.)

Beachten Sie, dass die Ausführungszeit dieser Funktion im Allgemeinen mehr als 20 Minuten beträgt, daher sollten Sie dieses Tutorial in zwei Abschnitte aufteilen: den ersten, in dem sie es durchlesen und die Aufträge starten, und den zweiten Abschnitt einige Stunden später (damit Sie genügend Zeit haben, die Aufträge abzuschließen), um mit den Ergebnissen der Aufträge zu arbeiten.

Hintergrund

Dieses Tutorial vermittelt auf einem Einführungsniveau, wie man die QUICK-PDE-Funktion zur Lösung komplexer Mehrphysikprobleme auf 156Q Heron R2 QPUs unter Verwendung des H-DES (Hybrid Differential Equation Solver) ColibriTD's einsetzt. Der zugrunde liegende Algorithmus wird in der H-DES-Veröffentlichung beschrieben. Beachten Sie, dass dieser Solver auch nichtlineare Gleichungen lösen kann.

Multiphysik-Probleme - einschließlich Flüssigkeitsdynamik, Wärmediffusion und Material verformung, um nur einige zu nennen - lassen sich durch partielle Differentialgleichungen (PDEs) beschrieben werden.

Solche Probleme sind für verschiedene Branchen von großer Bedeutung und stellen einen einen wichtigen Zweig der angewandten Mathematik. Das Lösen nichtlinearer gekoppelter multivariater PDEs mit klassischen Werkzeugen bleibt jedoch aufgrund des eine exponentiell große Menge an Ressourcen benötigt.

Diese Funktion eignet sich für Gleichungen mit zunehmender Komplexität und variablen und ist der erste Schritt zur Erschließung von Möglichkeiten, die früher als unlösbar galten. Um ein durch PDEs modelliertes Problem vollständig zu beschreiben, ist es muss man die Anfangs- und Randbedingungen kennen. Diese können stark die Lösung der PDE und den Weg zu ihrer Lösung stark verändern.

In diesem Lernprogramm lernen Sie, wie man:

1. Definieren Sie die Parameter der Anfangsbedingungsfunktion.
2. Stellen Sie die Qubit-Nummer (die zur Codierung der Funktion der Differentialgleichung verwendet wird), die Tiefe und die Schusszahl ein.
3. Führen Sie QUICK-PDE aus, um die zugrunde liegende Differentialgleichung zu lösen.

Anforderungen

Bevor Sie mit diesem Tutorial beginnen, sollten Sie sicherstellen, dass Sie Folgendes installiert haben:

• Qiskit SDK v2.0 oder höher (pip install qiskit)
• Qiskit-Funktionskatalog (pip install qiskit-ibm-catalog)
• Matplotlib (pip install matplotlib)
• Aufruf der Funktion QUICK-PDE. Füllen Sie das Formular aus, um Zugang zu beantragen.

Einrichtung

Authentifizieren Sie sich mit Ihrem API-Schlüssel und wählen Sie die Funktion wie folgt:

import numpy as np
import matplotlib.pyplot as plt
from qiskit_ibm_catalog import QiskitFunctionsCatalog

catalog = QiskitFunctionsCatalog(
    channel="ibm_quantum_platform",
    instance="INSTANCE_CRN",
    # For `token`, use the 44-character API_KEY you created
    # and saved from the IBM Quantum Platform Home dashboard
    token="YOUR_API_KEY",
)

quick = catalog.load("colibritd/quick-pde")

Schritt 1: Eigenschaften des zu lösenden Problems festlegen

Dieses Tutorium behandelt die Benutzererfahrung aus zwei Perspektiven: das physikalische Problem, das durch die Anfangsbedingungen bestimmt wird, und die algorithmische komponente bei der Lösung eines Beispiels aus der Strömungsdynamik auf einem Quantencomputer.

Die numerische Strömungsmechanik (CFD) hat ein breites Anwendungsspektrum, und daher ist es ist es wichtig, die zugrunde liegenden PDEs zu untersuchen und zu lösen. Eine wichtige Familie von PDEs sind die Navier-Stokes-Gleichungen, bei denen es sich um ein System nichtlinearer partieller differentialgleichungen, die die Bewegung von Fluiden beschreiben. Sie sind von großer Bedeutung für wissenschaftliche Probleme und technische Anwendungen.

Unter bestimmten Bedingungen lassen sich die Navier-Stokes-Gleichungen auf die Burgers'sche gleichung, eine Konvektions-Diffusions-Gleichung, die Phänomene beschreibt, die in der strömungsdynamik, Gasdynamik und nichtlineare Akustik, um nur einige zu nennen, durch die Modellierung von dissipative Systeme.

Die eindimensionale Version der Gleichung hängt von zwei Variablen ab: $t \in \mathbb{R}_{\geq 0}$ der Modellierung der zeitlichen Dimension, $x \in \mathbb{R}$ der Darstellung der räumlichen Dimension. Die allgemeine Form der Gleichung wird als die viskose Burgers'sche Gleichung und lautet

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

wobei$u(x,t)$ das Strömungsgeschwindigkeitsfeld an einer bestimmten Position$x$ und zu einem bestimmten Zeitpunkt$t$ ist und$\nu$ die Viskosität der Flüssigkeit ist. Die Viskosität ist eine wichtige Eigenschaft einer Flüssigkeit, die deren geschwindigkeitsabhängigen Widerstand gegen Bewegung oder Verformung angibt und somit eine entscheidende Rolle bei der Bestimmung der Dynamik einer Flüssigkeit spielt. Wenn die Viskosität der Flüssigkeit null ist ( $\nu = 0$ ), wird die Gleichung zu einer Erhaltungsgleichung, in der sich aufgrund des fehlenden inneren Widerstands Unstetigkeiten (Stoßwellen) bilden können. In diesem Fall wird die Gleichung als viskositätsfreie Burgers- Gleichung bezeichnet und stellt einen Sonderfall der nichtlinearen Wellengleichung dar.

Streng genommen kommen nichtviskose Strömungen in der Natur nicht vor, aber bei der Modellierung der aerodynamischen Strömung kann aufgrund der infinitesimalen Wirkung des Transports eine eine nicht viskose Beschreibung des Problems nützlich sein kann. Überraschenderweise befassen sich mehr als 70 % der aerodynamischen Theorie mit nichtviskosen Strömungen befasst.

Dieses Tutorial verwendet die Burgers-Gleichung für nichtviskose Strömungen als CFD-Beispiel zur Lösung auf QPUs v IBM® s unter Verwendung von QUICK-PDE, die wie folgt lautet:

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

Die Anfangsbedingung für dieses Problem wird durch eine lineare Funktion festgelegt: $u(t=0,x) = ax + b$, $\text{with }a,b\in\mathbb{R}$, wobei $a$ und $b$ beliebige Konstanten sind, die die Form der Lösung beeinflussen. Sie können die Werte für „ $a$ “ und „ $b$ “ anpassen und beobachten, wie sich diese auf den Lösungsprozess und das Ergebnis auswirken. Für diese Wahl der Anfangsbedingung lässt die Gleichung die folgende analytische Lösung zu:

$u(t, x) = \frac{ax + b}{at + 1},$

Dies gilt, solange sich kein Stoß gebildet hat, d. h. für $t < -1/a$, wenn $a < 0$.

job = quick.run(
    use_case="CFD_BURGER",
    physical_parameters={"a": 1.0, "b": 1.0},
)

print(job.result())

Output:

{'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): Optimieren Sie das Problem für die Ausführung auf Quantenhardware

Standardmäßig verwendet der Solver physikalisch informierte Parameter, d. h. anfängliche Schaltkreisparameter für eine bestimmte Qubitanzahl und -tiefe, von denen unser Solver ausgeht.

Die Schüsse gehören ebenfalls zu den Parametern mit einem Standardwert, da die Feinabstimmung wichtig ist.

Je nach der Konfiguration, die Sie zu lösen versuchen, müssen die Parameter des Algorithmus parameter des Algorithmus angepasst werden, um zufriedenstellende Lösungen zu erzielen; zum Beispiel kann er kann mehr oder weniger Qubits pro Variable $t$ und $x$ erfordern, je nach $a$ und $b$. Im Folgenden wird die Anzahl der Qubits pro Funktion pro variable, die Tiefe pro Funktion und die Anzahl der Schüsse.

Sie können auch sehen, wie Sie das Backend und den Ausführungsmodus angeben.

Darüber hinaus können physikalisch informierte Parameter den Optimierungsprozess in eine falsche Richtung lenken in eine falsche Richtung lenken; in diesem Fall können Sie sie deaktivieren, indem Sie die initialization strategie auf "RANDOM" setzen.

job_2 = quick.run(
    use_case="CFD_BURGER",
    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())

Output:

{'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: Vergleichen Sie die Leistungsfähigkeit der Algorithmen

Sie können den Konvergenzprozess unserer Lösung (HDES) mit job_2 der Leistung eines physikalisch informierten neuronalen Netzwerks (PINN) und eines Solvers vergleichen (siehe Artikel und zugehöriges 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 läuft wie folgt ab:

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 Verlust unterhalb von nach 28 Iterationen und mit der Optimierung nur 0.0015 weniger klassischer Parameter erreicht werden kann.

Nun können wir dies mit der PINN-Lösung vergleichen, wobei die in der Veröffentlichung vorgeschlagene Standardkonfiguration mit einem gradientenbasierten Optimierer verwendet wird. Das Äquivalent zu unserer Schaltung mit 13 zu optimierenden Parametern ist das neuronale Netzwerk, das mindestens acht Schichten mit jeweils 20 Neuronen erfordert und somit die Optimierung von 3021 Parametern beinhaltet. Dann wird der Zielverlust in Schritt 315 erreicht, Verlust: 0.0014988397.

Da wir einen fairen Vergleich anstellen wollen, sollten wir in beiden Fällen denselben Optimierer verwenden. Die geringste Anzahl an Iterationen, die wir für 12 Schichten mit jeweils 20 Neuronen gefunden haben = 4701 Parameter:

(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

Sie können dasselbe mit Ihren Daten aus tun job_2 und einen Vergleich mit der PINN-Lösung erstellen.

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

Schritt 4: Verwenden Sie das Ergebnis

Mit Ihrer Lösung können Sie nun entscheiden, was Sie damit machen wollen. Im Folgenden wird gezeigt, wie das Ergebnis dargestellt werden kann.

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:

Beachten Sie den Unterschied der Ausgangsbedingungen für den zweiten Durchlauf und dessen 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()

Output:

Tutorial-Umfrage

Bitte nehmen Sie sich eine Minute Zeit, um uns Ihr Feedback zu diesem Tutorial zu geben. Ihre Erkenntnisse helfen uns, unser Inhaltsangebot und unsere Nutzererfahrung zu verbessern:

Link zur Umfrage