Modelowanie przepływu nielepkiego płynu za pomocą QUICK-PDE

Uwaga

Qiskit Functions to funkcja eksperymentalna dostępna wyłącznie dla użytkowników planów IBM Quantum® Premium Plan, Flex Plan oraz On-Prem (przez IBM Quantum Platform API). Są one w fazie podglądu i mogą ulec zmianie.

Szacowany czas użycia: 50 minut na procesorze Heron r2. (UWAGA: To tylko szacunek. Rzeczywisty czas może się różnić.)

Pamiętaj, że czas wykonania tej funkcji jest zazwyczaj dłuższy niż 20 minut, dlatego warto podzielić ten samouczek na dwie sekcje: pierwszą, w której zapoznasz się z treścią i uruchomisz zadania, oraz drugą — kilka godzin później (dając wystarczająco dużo czasu na zakończenie zadań) — w której przeanalizujesz wyniki.

Tło

Ten samouczek ma na celu wprowadzenie w podstawy korzystania z funkcji QUICK-PDE do rozwiązywania złożonych problemów wielofizycznych na procesorach QPU Heron R2 o 156 qubitach, przy użyciu narzędzia H-DES (Hybrid Differential Equation Solver) firmy ColibriTD. Leżący u podstaw algorytm opisano w artykule o H-DES. Pamiętaj, że solver ten potrafi też rozwiązywać równania nieliniowe.

Problemy wielofizyczne — obejmujące m.in. dynamikę płynów, dyfuzję ciepła i deformację materiałów — można powszechnie opisywać za pomocą Równań Różniczkowych Cząstkowych (PDE).

Takie problemy są niezwykle istotne dla różnych gałęzi przemysłu i stanowią ważną dziedzinę matematyki stosowanej. Jednak rozwiązywanie nieliniowych, wielowymiarowych, sprzężonych PDE za pomocą klasycznych narzędzi pozostaje trudne ze względu na wykładnicze zapotrzebowanie na zasoby obliczeniowe.

Ta funkcja jest odpowiednia dla równań o rosnącej złożoności i liczbie zmiennych, i stanowi pierwszy krok ku możliwościom, które dotychczas uważano za nierozwiązywalne. Aby w pełni opisać problem modelowany przez PDE, konieczna jest znajomość warunków początkowych i brzegowych. Mogą one znacząco wpływać na rozwiązanie PDE oraz na sposób jego znalezienia.

Ten samouczek uczy, jak:

1. Definiować parametry funkcji warunków początkowych.
2. Dostosowywać liczbę qubitów (używanych do kodowania funkcji równania różniczkowego), głębokość i liczbę pomiarów.
3. Uruchamiać QUICK-PDE w celu rozwiązania leżącego u podstaw równania różniczkowego.

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane:

• Qiskit SDK v2.0 lub nowszy (pip install qiskit)
• Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
• Matplotlib (pip install matplotlib)
• Dostęp do funkcji QUICK-PDE. Wypełnij formularz, aby poprosić o dostęp.

Konfiguracja

Uwierzytelnij się za pomocą swojego klucza API i wybierz funkcję w następujący sposób:

# 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")

Krok 1: Ustaw właściwości rozwiązywanego problemu

Ten samouczek obejmuje doświadczenie użytkownika z dwóch perspektyw: problemu fizycznego wyznaczonego przez warunki początkowe oraz komponentu algorytmicznego przy rozwiązywaniu przykładu z dynamiki płynów na komputerze kwantowym.

Obliczeniowa Dynamika Płynów (CFD) ma szeroki zakres zastosowań, dlatego ważne jest badanie i rozwiązywanie leżących u jej podstaw PDE. Ważną rodziną PDE są równania Naviera-Stokesa — układy nieliniowych równań różniczkowych cząstkowych opisujące ruch płynów. Mają one ogromne znaczenie dla problemów naukowych i zastosowań inżynieryjnych.

W pewnych warunkach równania Naviera-Stokesa upraszczają się do równania Burgersa — równania konwekcji-dyfuzji opisującego zjawiska występujące w dynamice płynów, dynamice gazów i nieliniowej akustyce, modelując układy dysypatywne.

Jednowymiarowa wersja tego równania zależy od dwóch zmiennych: $t \in \mathbb{R}_{\geq 0}$ modelującej wymiar czasowy oraz $x \in \mathbb{R}$ reprezentującej wymiar przestrzenny. Ogólna postać równania nosi nazwę lepkiego równania Burgersa i ma postać:

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

gdzie $u(x,t)$ to pole prędkości płynu w danym położeniu $x$ i czasie $t$, a $\nu$ to lepkość płynu. Lepkość jest ważną właściwością płynu, która mierzy jego zależną od tempa odporność na ruch lub odkształcenie, odgrywając kluczową rolę w wyznaczaniu dynamiki płynu. Gdy lepkość płynu wynosi zero ($\nu$ = 0), równanie staje się równaniem zachowania, które może rozwijać nieciągłości (fale uderzeniowe) z powodu braku wewnętrznego oporu. W tym przypadku równanie nosi nazwę nielepkiego równania Burgersa i jest szczególnym przypadkiem nieliniowego równania falowego.

Ściśle rzecz biorąc, przepływy nielepkie nie występują w naturze, ale przy modelowaniu przepływu aerodynamicznego — ze względu na nieskończenie mały efekt transportu — użycie nielepkiego opisu problemu może być przydatne. Co zaskakujące, ponad 70% teorii aerodynamiki dotyczy właśnie przepływów nielepkich.

Ten samouczek używa nielepkiego równania Burgersa jako przykładu CFD do rozwiązania na QPU IBM® za pomocą QUICK-PDE, danego równaniem:

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

Warunek początkowy dla tego problemu jest ustawiony jako funkcja liniowa: $u(t=0,x) = ax + b,\text{ with }a,b\in\mathbb{R}$ gdzie $a$ i $b$ to dowolne stałe wpływające na kształt rozwiązania. Możesz zmieniać $a$ i $b$, obserwując jak wpływają na proces rozwiązywania i na wynik.

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  ])}}

Krok 2 (jeśli potrzeba): Zoptymalizuj problem pod kątem wykonania na sprzęcie kwantowym

Domyślnie solver używa parametrów wyznaczonych fizycznie, czyli początkowych parametrów obwodu dla danej liczby qubitów i głębokości, od których nasz solver rozpocznie działanie.

Pomiary (shots) są również częścią parametrów z domyślną wartością, ponieważ ich precyzyjne dostrojenie jest ważne.

W zależności od konfiguracji, którą próbujesz rozwiązać, parametry algorytmu umożliwiające uzyskanie zadowalających rozwiązań mogą wymagać dostosowania; na przykład może być potrzeba użycia większej lub mniejszej liczby qubitów na zmienną $t$ i $x$, w zależności od $a$ i $b$. Poniżej dostosowywana jest liczba qubitów na funkcję na zmienną, głębokość na funkcję oraz liczba pomiarów.

Możesz też zobaczyć, jak określić Backend i tryb wykonania.

Ponadto, parametry wyznaczone fizycznie mogą kierować proces optymalizacji w złym kierunku; w takim przypadku możesz je wyłączyć, ustawiając strategię initialization na "RANDOM".

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  ])}}

Krok 3: Porównanie wydajności algorytmów

Możesz porównać proces zbieżności naszego rozwiązania (HDES) z zadania job_2 z wydajnością algorytmu i solvera opartego na fizycznie uwarunkowanych sieciach neuronowych (PINN) (zob. artykuł i powiązane repozytorium GitHub).

W przykładzie wyjścia zadania job_2 (podejście kwantowe) tylko 13 parametrów (12 parametrów obwodu oraz 1 parametr skalowania) jest optymalizowanych przez solver klasyczny. Proces zbieżności wygląda następująco:

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

Oznacza to, że stratę poniżej 0,0015 można osiągnąć po 28 iteracjach, optymalizując zaledwie kilka parametrów klasycznych.

Możemy teraz porównać to z rozwiązaniem PINN w domyślnej konfiguracji sugerowanej przez artykuł, używając optymalizatora gradientowego. Odpowiednikiem naszego obwodu z 13 parametrami do optymalizacji jest sieć neuronowa, która wymaga co najmniej ośmiu warstw po 20 neuronów, a zatem wiąże się z optymalizacją 3021 parametrów. Docelowa strata zostaje osiągnięta w kroku 315, loss: 0.0014988397.

Ponieważ zależy nam na rzetelnym porównaniu, powinniśmy użyć tego samego optymalizatora w obu przypadkach. Najniższa liczba iteracji znaleziona dla 12 warstw po 20 neuronów = 4701 parametrów:

(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

Możesz zrobić to samo ze swoimi danymi z zadania job_2 i narysować porównanie z rozwiązaniem PINN.

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

Krok 4: Wykorzystanie wyników

Mając gotowe rozwiązanie, możesz zdecydować, co z nim zrobić. Poniżej pokazano, jak wykreślić wynik.

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()

Zwróć uwagę na różnicę w warunku początkowym dla drugiego uruchomienia i jej wpływ na wynik:

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()

Ankieta dotycząca samouczka

Poświęć chwilę na przekazanie opinii na temat tego samouczka. Twoje spostrzeżenia pomogą nam ulepszyć naszą ofertę treści i doświadczenie użytkownika:

Link do ankiety

Note: This survey is provided by IBM Quantum and relates to the original English content. To give feedback on doQumentation's website, translations, or code execution, please open a GitHub issue.