Моделювання потоку нев'язкої рідини за допомогою QUICK-PDE

Примітка

Qiskit Functions — це експериментальна функція, доступна лише для користувачів планів IBM Quantum® Premium Plan, Flex Plan та On-Prem (через IBM Quantum Platform API). Вона перебуває у статусі попереднього випуску та може змінюватися.

Орієнтовний час використання: 50 хвилин на процесорі Heron r2. (ПРИМІТКА: Це лише оцінка. Ваш фактичний час виконання може відрізнятися.)

Зверніть увагу, що час виконання цієї функції зазвичай перевищує 20 хвилин, тому Ви можете розділити цей посібник на два розділи: перший, у якому Ви ознайомитеся з матеріалом і запустите завдання, та другий — через кілька годин (надавши достатньо часу для завершення завдань), щоб працювати з результатами завдань.

Передумови

Цей посібник має на меті навчити на вступному рівні, як використовувати функцію QUICK-PDE для розв'язання складних мультифізичних задач на 156-кубітних QPU Heron R2 за допомогою H-DES (гібридного розв'язувача диференціальних рівнянь) від ColibriTD. Базовий алгоритм описано у статті про H-DES. Зверніть увагу, що цей розв'язувач також може розв'язувати нелінійні рівняння.

Мультифізичні задачі — зокрема динаміка рідин, теплова дифузія та деформація матеріалів — можуть бути повсюдно описані за допомогою диференціальних рівнянь у частинних похідних (ДРЧП).

Такі задачі є надзвичайно актуальними для різних галузей промисловості та становлять важливу гілку прикладної математики. Однак розв'язання нелінійних багатовимірних зв'язаних ДРЧП класичними методами залишається складним завданням через необхідність експоненціально великої кількості ресурсів.

Ця функція підходить для рівнянь зі зростаючою складністю та кількістю змінних і є першим кроком до відкриття можливостей, які раніше вважалися нерозв'язними. Для повного опису задачі, змодельованої за допомогою ДРЧП, необхідно знати початкові та граничні умови. Вони можуть суттєво змінити розв'язок ДРЧП та шлях до знаходження цього розв'язку.

Цей посібник навчить Вас:

1. Визначати параметри функції початкових умов.
2. Налаштовувати кількість кубітів (що використовуються для кодування функції диференціального рівняння), глибину та кількість вимірювань.
3. Запускати QUICK-PDE для розв'язання базового диференціального рівняння.

Вимоги

Перед початком роботи з цим посібником переконайтеся, що у Вас встановлено наступне:

• Qiskit SDK v2.0 або новішої версії (pip install qiskit)
• Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
• Matplotlib (pip install matplotlib)
• Доступ до функції QUICK-PDE. Заповніть форму для запиту доступу.

Налаштування

Автентифікуйтесь за допомогою Вашого API-ключа та оберіть функцію наступним чином:

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

Крок 1: Встановлення властивостей задачі для розв'язання

Цей посібник розглядає користувацький досвід з двох перспектив: фізичну задачу, визначену початковими умовами, та алгоритмічну складову розв'язання прикладу з динаміки рідин на квантовому комп'ютері.

Обчислювальна гідродинаміка (CFD) має широкий спектр застосувань, і тому важливо вивчати та розв'язувати відповідні ДРЧП. Важливим сімейством ДРЧП є рівняння Нав'є-Стокса — система нелінійних диференціальних рівнянь у частинних похідних, що описує рух рідин. Вони є надзвичайно важливими для наукових задач та інженерних застосувань.

За певних умов рівняння Нав'є-Стокса зводяться до рівняння Бюргерса — рівняння конвекції-дифузії, що описує явища в динаміці рідин, газодинаміці та нелінійній акустиці, зокрема шляхом моделювання дисипативних систем.

Одновимірна версія рівняння залежить від двох змінних: $t \in \mathbb{R}_{\geq 0}$, що моделює часовий вимір, та $x \in \mathbb{R}$, що представляє просторовий вимір. Загальна форма рівняння називається в'язким рівнянням Бюргерса і має вигляд:

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

де $u(x,t)$ — поле швидкості рідини у заданій позиції $x$ та момент часу $t$, а $\nu$ — в'язкість рідини. В'язкість — це важлива властивість рідини, що вимірює її залежний від швидкості опір руху або деформації, і тому вона відіграє вирішальну роль у визначенні динаміки рідини. Коли в'язкість рідини дорівнює нулю ($\nu$ = 0), рівняння стає рівнянням збереження, яке може розвивати розриви (ударні хвилі) через відсутність внутрішнього опору. У цьому випадку рівняння називається нев'язким рівнянням Бюргерса і є окремим випадком нелінійного хвильового рівняння.

Строго кажучи, нев'язкі потоки не існують у природі, але при моделюванні аеродинамічного потоку, через нескінченно малий вплив переносу, використання нев'язкого опису задачі може бути корисним. Дивовижно, але понад 70% аеродинамічної теорії має справу з нев'язкими потоками.

У цьому посібнику використовується нев'язке рівняння Бюргерса як приклад обчислювальної гідродинаміки для розв'язання на QPU IBM® за допомогою QUICK-PDE, задане рівнянням:

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

Початкова умова для цієї задачі задається лінійною функцією: $u(t=0,x) = ax + b,\text{ with }a,b\in\mathbb{R}$ де $a$ та $b$ — довільні константи, що впливають на форму розв'язку. Ви можете змінювати $a$ та $b$ і спостерігати, як вони впливають на процес розв'язання та сам розв'язок.

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

Крок 2 (за потреби): Оптимізація задачі для виконання на квантовому обладнанні

За замовчуванням розв'язувач використовує фізично обґрунтовані параметри — початкові параметри схеми для заданої кількості кубітів та глибини, від яких наш розв'язувач починає роботу.

Кількість вимірювань також є частиною параметрів зі значенням за замовчуванням, оскільки їх точне налаштування є важливим.

Залежно від конфігурації, яку Ви намагаєтеся розв'язати, параметри алгоритму для досягнення задовільних розв'язків можуть потребувати адаптації; наприклад, це може вимагати більше або менше кубітів на змінну $t$ та $x$, залежно від $a$ та $b$. Нижче наведено налаштування кількості кубітів на функцію на змінну, глибини на функцію та кількості вимірювань.

Ви також можете побачити, як вказати бекенд та режим виконання.

Крім того, фізично обґрунтовані параметри можуть направити процес оптимізації в неправильному напрямку; у такому випадку Ви можете вимкнути їх, встановивши стратегію initialization на "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  ])}}

Крок 3: Порівняння продуктивності алгоритмів

Ви можете порівняти процес збіжності нашого розв'язку (HDES) з job_2 із продуктивністю алгоритму фізично-інформованих нейронних мереж (PINN) та розв'язувача (див. статтю та пов'язаний репозиторій GitHub).

У прикладі виводу job_2 (квантовий підхід) лише 13 параметрів (12 параметрів схеми плюс 1 параметр масштабування) оптимізуються за допомогою класичного розв'язувача. Процес збіжності виглядає наступним чином:

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

Це означає, що значення функції втрат нижче 0.0015 може бути досягнуто після 28 ітерацій і при оптимізації лише кількох класичних параметрів.

Тепер ми можемо порівняти те саме з розв'язком PINN із конфігурацією за замовчуванням, запропонованою у статті, з використанням градієнтного оптимізатора. Еквівалентом нашої схеми з 13 параметрами для оптимізації є нейронна мережа, яка потребує щонайменше вісім шарів по 20 нейронів і, таким чином, передбачає оптимізацію 3021 параметра. Тоді цільове значення функції втрат досягається на кроці 315, loss: 0.0014988397.

Оскільки ми хочемо провести справедливе порівняння, нам слід використовувати однаковий оптимізатор в обох випадках. Найменша кількість ітерацій, яку ми знайшли для 12 шарів по 20 нейронів = 4701 параметр:

(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

Ви можете зробити те саме зі своїми даними з job_2 та побудувати графік порівняння з розв'язком PINN.

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

Крок 4: Використання результату

Маючи Ваш розв'язок, Ви тепер можете обрати, що з ним робити. Нижче продемонстровано, як побудувати графік результату.

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

Зверніть увагу на різницю початкової умови для другого запуску та її вплив на результат:

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

Опитування щодо посібника

Будь ласка, приділіть хвилину, щоб залишити відгук про цей посібник. Ваші думки допоможуть нам покращити наші навчальні матеріали та досвід користувачів:

Link to survey