Modelování proudění neviskózní tekutiny pomocí QUICK-PDE

Poznámka

Qiskit Functions jsou experimentální funkce dostupné pouze uživatelům IBM Quantum® Premium Plan, Flex Plan a On-Prem (prostřednictvím IBM Quantum Platform API) Plan. Jsou ve stavu náhledového vydání a mohou se měnit.

Odhadovaná spotřeba: 50 minut na procesoru Heron r2. (POZNÁMKA: Jedná se pouze o odhad. Skutečný čas běhu se může lišit.)

Vezmi na vědomí, že doba provádění této funkce je obecně delší než 20 minut, takže možná budeš chtít rozdělit tento tutoriál do dvou částí: první, ve které si ho přečteš a spustíš úlohy, a druhé o několik hodin později (s dostatečným časem na dokončení úloh), ve které budeš pracovat s výsledky úloh.

Pozadí

Tento tutoriál si klade za cíl na úvodní úrovni naučit, jak používat funkci QUICK-PDE k řešení složitých multifyzikálních problémů na 156Q Heron R2 QPU pomocí H-DES (Hybrid Differential Equation Solver) od ColibriTD. Základní algoritmus je popsán v článku o H-DES. Vezmi na vědomí, že tento solver dokáže řešit i nelineární rovnice.

Multifyzikální problémy – včetně dynamiky tekutin, tepelné difuze a deformace materiálů, abychom jmenovali alespoň některé – lze všudypřítomně popsat pomocí parciálních diferenciálních rovnic (PDE).

Takové problémy jsou vysoce relevantní pro různá průmyslová odvětví a tvoří důležitou větev aplikované matematiky. Řešení nelineárních multivariátních vázaných PDE klasickými nástroji však zůstává náročné kvůli požadavku exponenciálně velkého množství zdrojů.

Tato funkce je vhodná pro rovnice s rostoucí složitostí a proměnnými a je prvním krokem k odemknutí možností, které byly dříve považovány za neřešitelné. Aby bylo možné plně popsat problém modelovaný PDE, je nutné znát počáteční a okrajové podmínky. Ty mohou výrazně změnit řešení PDE i cestu k nalezení tohoto řešení.

Tento tutoriál tě naučí, jak:

1. Definovat parametry funkce počáteční podmínky.
2. Upravit počet qubitů (použitých pro kódování funkce diferenciální rovnice), hloubku a počet vzorků (shots).
3. Spustit QUICK-PDE k řešení příslušné diferenciální rovnice.

Požadavky

Před zahájením tohoto tutoriálu se ujisti, že máš nainstalované následující:

• Qiskit SDK v2.0 nebo novější (pip install qiskit)
• Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
• Matplotlib (pip install matplotlib)
• Přístup k funkci QUICK-PDE. Vyplň formulář pro žádost o přístup.

Nastavení

Autentizuj se pomocí svého API klíče a vyber funkci takto:

# 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: Nastavení vlastností problému k řešení

Tento tutoriál pokrývá uživatelský pohled ze dvou perspektiv: fyzikálního problému určeného počátečními podmínkami a algoritmické složky při řešení příkladu z dynamiky tekutin na kvantovém počítači.

Výpočetní dynamika tekutin (CFD) má širokou škálu aplikací, a proto je důležité studovat a řešit příslušné PDE. Důležitou rodinou PDE jsou Navier-Stokesovy rovnice, což je soustava nelineárních parciálních diferenciálních rovnic popisujících pohyb tekutin. Jsou vysoce relevantní pro vědecké problémy a inženýrské aplikace.

Za určitých podmínek se Navier-Stokesovy rovnice redukují na Burgersovu rovnici, rovnici konvekce–difuze popisující jevy vyskytující se v dynamice tekutin, dynamice plynů a nelineární akustice, abychom jmenovali alespoň některé, prostřednictvím modelování disipativních systémů.

Jednorozměrná verze rovnice závisí na dvou proměnných: $t \in \mathbb{R}_{\geq 0}$ modelující časovou dimenzi, $x \in \mathbb{R}$ reprezentující prostorovou dimenzi. Obecná forma rovnice se nazývá viskózní Burgersova rovnice a má tvar:

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

kde $u(x,t)$ je pole rychlosti tekutiny v dané poloze $x$ a čase $t$ a $\nu$ je viskozita tekutiny. Viskozita je důležitá vlastnost tekutiny, která měří její rychlostně závislý odpor vůči pohybu nebo deformaci, a proto hraje klíčovou roli při určování dynamiky tekutiny. Když je viskozita tekutiny nulová ($\nu$ = 0), rovnice se stává rovnicí zachování, která může vykazovat nespojitosti (rázové vlny) v důsledku absence vnitřního odporu. V tomto případě se rovnice nazývá neviskózní Burgersova rovnice a je zvláštním případem nelineární vlnové rovnice.

Přísně vzato, neviskózní proudění se v přírodě nevyskytuje, ale při modelování aerodynamického proudění může být kvůli nekonečně malému efektu transportu použití neviskózního popisu problému užitečné. Překvapivě více než 70 % aerodynamické teorie se zabývá neviskózním prouděním.

Tento tutoriál používá neviskózní Burgersovu rovnici jako příklad CFD k řešení na IBM® QPU pomocí QUICK-PDE, danou rovnicí:

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

Počáteční podmínka pro tento problém je nastavena na lineární funkci: $u(t=0,x) = ax + b,\text{ with }a,b\in\mathbb{R}$ kde $a$ a $b$ jsou libovolné konstanty ovlivňující tvar řešení. Můžeš upravit $a$ a $b$ a sledovat, jak ovlivňují proces řešení a výsledné řešení.

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 (v případě potřeby): Optimalizace problému pro spuštění na kvantovém hardwaru

Ve výchozím nastavení solver používá fyzikálně informované parametry, což jsou počáteční parametry Circuit pro daný počet qubitů a hloubku, od nichž náš solver začne.

Vzorky (shots) jsou také součástí parametrů s výchozí hodnotou, protože jejich jemné ladění je důležité.

V závislosti na konfiguraci, kterou se snažíš řešit, může být nutné upravit parametry algoritmu k dosažení uspokojivých řešení; například to může vyžadovat více nebo méně qubitů na proměnnou $t$ a $x$ v závislosti na $a$ a $b$. Následující kód upravuje počet qubitů na funkci na proměnnou, hloubku na funkci a počet vzorků.

Můžeš také vidět, jak zadat Backend a režim spuštění.

Navíc fyzikálně informované parametry mohou nasměrovat optimalizační proces špatným směrem; v takovém případě je můžeš zakázat nastavením strategie 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: Porovnání výkonu algoritmů

Můžeš porovnat proces konvergence našeho řešení (HDES) z job_2 s výkonem algoritmu fyzikálně informovaných neuronových sítí (PINN) a příslušného solveru (viz článek a přidružený repozitář na GitHubu).

V příkladu výstupu job_2 (kvantový přístup) je s klasickým solverem optimalizováno pouze 13 parametrů (12 parametrů Circuit plus 1 škálovací parametr). Proces konvergence vypadá takto:

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

To znamená, že ztrátu pod 0,0015 lze dosáhnout po 28 iteracích, a to optimalizací pouze několika klasických parametrů.

Teď můžeme totéž porovnat s řešením PINN s výchozí konfigurací navrženou v článku, která využívá optimizer založený na gradientech. Ekvivalentem našeho Circuit s 13 optimalizovanými parametry je neuronová síť, která vyžaduje alespoň osm vrstev o 20 neuronech, a tedy optimalizaci 3021 parametrů. Cílové ztráty je pak dosaženo v kroku 315, loss: 0.0014988397.

Nyní, protože chceme provést spravedlivé srovnání, měli bychom v obou případech použít stejný optimizer. Nejnižší počet iterací, který jsme nalezli pro 12 vrstev o 20 neuronech = 4701 parametrů:

(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

Totéž můžeš provést se svými daty z job_2 a vykreslit srovnání s řešením PINN.

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

Krok 4: Použití výsledku

Se svým řešením si teď můžeš vybrat, co s ním uděláš. Následující ukázka demonstruje, jak výsledek vykreslit.

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

Všimni si rozdílu počáteční podmínky pro druhý běh a jeho vlivu na výsledek:

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

Anketa k tutoriálu

Věnuj prosím chvilku a poskytni nám zpětnou vazbu k tomuto tutoriálu. Tvé postřehy nám pomohou zlepšit obsah a uživatelský zážitek:

Odkaz na anketu

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.