Modellare un fluido non viscoso in movimento usando QUICK-PDE

Nota

Le Qiskit Functions sono una funzionalità sperimentale disponibile solo per gli utenti del Piano Premium IBM Quantum®, Piano Flex e Piano On-Prem (tramite IBM Quantum Platform API). Sono in stato di rilascio in anteprima e soggette a modifiche.

Stima di utilizzo: 50 minuti su un processore Heron r2. (NOTA: Questa è solo una stima. Il tempo di esecuzione effettivo potrebbe variare.)

Notate che il tempo di esecuzione di questa funzione è generalmente superiore a 20 minuti, quindi potreste voler dividere questo tutorial in due sezioni: la prima, in cui lo leggete e lanciate i job, e la seconda qualche ora dopo (concedendo ampio tempo per il completamento dei job), per lavorare con i risultati dei job.

Contesto

Questo tutorial cerca di insegnare a livello introduttivo come utilizzare la funzione QUICK-PDE per risolvere problemi multifisici complessi su QPU Heron R2 da 156Q utilizzando l'H-DES (Hybrid Differential Equation Solver) di ColibriTD. L'algoritmo sottostante è descritto nel paper sull'H-DES. Notate che questo risolutore può anche risolvere equazioni non lineari.

I problemi multifisici - tra cui dinamica dei fluidi, diffusione del calore e deformazione dei materiali, solo per citarne alcuni - possono essere descritti ubiquitariamente da Equazioni Differenziali Parziali (PDE).

Tali problemi sono altamente rilevanti per varie industrie e costituiscono un importante ramo della matematica applicata. Tuttavia, risolvere PDE accoppiate multivariabili non lineari con strumenti classici rimane una sfida a causa del requisito di una quantità esponenzialmente grande di risorse.

Questa funzione è appropriata per equazioni con crescente complessità e variabili, ed è il primo passo per sbloccare possibilità che un tempo erano considerate intrattabili. Per descrivere completamente un problema modellato da PDE, è necessario conoscere le condizioni iniziali e al contorno. Queste possono fortemente modificare la soluzione della PDE e il percorso per trovare la loro soluzione.

Questo tutorial vi insegna come:

Definire i parametri della funzione di condizione iniziale.
Regolare il numero di qubit (utilizzati per codificare la funzione dell'equazione differenziale), la profondità e il numero di shot.
Eseguire QUICK-PDE per risolvere l'equazione differenziale sottostante.

Requisiti

Prima di iniziare questo tutorial, assicuratevi di avere installato quanto segue:

Qiskit SDK v2.0 o successivo (pip install qiskit)
Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
Matplotlib (pip install matplotlib)
Accesso alla funzione QUICK-PDE. Compilate il modulo per richiedere l'accesso.

Configurazione

Autenticatevi utilizzando la vostra chiave API e selezionate la funzione come segue:

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

Step 1: Impostare le proprietà del problema da risolvere

Questo tutorial copre l'esperienza utente da due prospettive: il problema fisico determinato dalle condizioni iniziali, e la componente algoritmica nel risolvere un esempio di dinamica dei fluidi su un computer quantistico.

La Computational Fluid Dynamics (CFD) ha una vasta gamma di applicazioni, ed è quindi importante studiare e risolvere le PDE sottostanti. Una famiglia importante di PDE sono le equazioni di Navier-Stokes, che sono un sistema di equazioni differenziali parziali non lineari che descrivono il moto dei fluidi. Sono altamente rilevanti per problemi scientifici e applicazioni ingegneristiche.

In determinate condizioni, le equazioni di Navier-Stokes si riducono all'equazione di Burgers, un'equazione di convezione-diffusione che descrive fenomeni che si verificano nella dinamica dei fluidi, nella dinamica dei gas e nell'acustica non lineare, solo per citarne alcuni, modellando sistemi dissipativi.

La versione unidimensionale dell'equazione dipende da due variabili: $t \in \mathbb{R}_{\geq 0}$ che modella la dimensione temporale, $x \in \mathbb{R}$ che rappresenta la dimensione spaziale. La forma generale dell'equazione è chiamata equazione di Burgers viscosa e si legge:

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

dove $u(x,t)$ è il campo di velocità del fluido in una data posizione $x$ e tempo $t$ , e $\nu$ è la viscosità del fluido. La viscosità è una proprietà importante di un fluido che misura la sua resistenza dipendente dalla velocità al movimento o alla deformazione, e quindi gioca un ruolo cruciale nella determinazione della dinamica di un fluido. Quando la viscosità del fluido è nulla ( $\nu$ = 0), l'equazione diventa un'equazione di conservazione che può sviluppare discontinuità (onde d'urto), a causa della mancanza di resistenza interna. In questo caso, l'equazione è chiamata equazione di Burgers non viscosa ed è un caso speciale di equazione d'onda non lineare.

Rigorosamente parlando, i flussi non viscosi non si verificano in natura, ma quando si modella il flusso aerodinamico, a causa dell'effetto infinitesimale del trasporto, l'utilizzo di una descrizione non viscosa del problema può essere utile. Sorprendentemente, più del 70% della teoria aerodinamica si occupa di flussi non viscosi.

Questo tutorial utilizza l'equazione di Burgers non viscosa come esempio di CFD da risolvere sulle QPU IBM® utilizzando QUICK-PDE, data dall'equazione:

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

La condizione iniziale per questo problema è impostata su una funzione lineare: $u(t=0,x) = ax + b,\text{ con }a,b\in\mathbb{R}$ dove $a$ e $b$ sono costanti arbitrarie che influenzano la forma della soluzione. Potete regolare $a$ e $b$ e vedere come influenzano il processo di risoluzione e la soluzione.

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

Step 2 (se necessario): Ottimizzare il problema per l'esecuzione su hardware quantistico

Per impostazione predefinita, il risolutore utilizza parametri fisicamente informati, che sono parametri iniziali del circuito per un dato numero di qubit e profondità da cui il nostro risolutore partirà.

Anche gli shot fanno parte dei parametri con un valore predefinito, poiché la loro regolazione fine è importante.

A seconda della configurazione che state cercando di risolvere, i parametri dell'algoritmo per ottenere soluzioni soddisfacenti potrebbero dover essere adattati; ad esempio, può richiedere più o meno qubit per variabile $t$ e $x$ , a seconda di $a$ e $b$ . Quanto segue regola il numero di qubit per funzione per variabile, la profondità per funzione e il numero di shot.

Potete anche vedere come specificare il backend e la modalità di esecuzione.

Inoltre, i parametri fisicamente informati potrebbero indirizzare il processo di ottimizzazione in una direzione sbagliata; in tal caso, potete disabilitarli impostando la strategia di initialization su "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  ])}}

Step 3: Confrontare le prestazioni degli algoritmi

Potete confrontare il processo di convergenza della nostra soluzione (HDES) di job_2 con le prestazioni di un algoritmo e risolutore basato su reti neurali informate dalla fisica (PINN) (consultate il paper e il repository GitHub associato).

Nell'esempio dell'output di job_2 (approccio basato su quantum), vengono ottimizzati solo 13 parametri (12 parametri del circuito più 1 parametro di scala) con il risolutore classico. Il processo di convergenza è il seguente:

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

Ciò significa che una loss inferiore a 0.0015 può essere raggiunta dopo 28 iterazioni, e ottimizzando solo pochi parametri classici.

Ora possiamo confrontare lo stesso con la soluzione PINN con la configurazione predefinita suggerita dal paper utilizzando un ottimizzatore basato sul gradiente. L'equivalente del nostro circuito con 13 parametri da ottimizzare è la rete neurale, che richiede almeno otto strati di 20 neuroni, e quindi comporta l'ottimizzazione di 3021 parametri. Poi, la loss target viene raggiunta allo Step 315, loss: 0.0014988397.

Graph showing PINN data compared with the HDES-Qiskit function.

Ora, poiché vogliamo fare un confronto equo, dovremmo usare lo stesso ottimizzatore in entrambi i casi. Il numero più basso di iterazioni che abbiamo trovato per 12 strati di 20 neuroni = 4701 parametri:

(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]
   20 5.398521572351456e-02 1.0e+00 9.98e-03  1e-02  1e-02 0:02.3
   40 5.444650724530220e-02 1.0e+00 9.97e-03  1e-02  1e-02 0:05.1
   60 4.447407275438309e-02 1.0e+00 9.95e-03  1e-02  1e-02 0:08.2
   80 2.068969979882240e-02 1.0e+00 9.94e-03  1e-02  1e-02 0:11.7
  120 1.028892211616039e-02 1.0e+00 9.91e-03  1e-02  1e-02 0:20.1
  140 5.140972323715687e-03 1.0e+00 9.90e-03  1e-02  1e-02 0:25.4
  180 3.811701666563749e-03 1.0e+00 9.87e-03  1e-02  1e-02 0:37.4
  200 3.189878538250923e-03 1.0e+00 9.85e-03  1e-02  1e-02 0:44.2
  240 2.547040116041899e-03 1.0e+00 9.83e-03  1e-02  1e-02 0:59.7
  280 2.166548743844032e-03 1.0e+00 9.80e-03  1e-02  1e-02 1:18.0
  300 1.783065614290535e-03 1.0e+00 9.79e-03  1e-02  1e-02 1:28.4
  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

Potete fare lo stesso con i vostri dati da job_2 e tracciare un confronto con la soluzione PINN.

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

Step 4: Utilizzare il risultato

Con la vostra soluzione, potete ora scegliere cosa farne. Il seguente esempio dimostra come tracciare il risultato.

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 of the previous code cell

Notate la differenza della condizione iniziale per la seconda esecuzione e il suo effetto sul risultato:

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 of the previous code cell

Sondaggio sul tutorial

Per favore, dedicate un minuto a fornire un feedback su questo tutorial. Le vostre osservazioni ci aiuteranno a migliorare i nostri contenuti e l'esperienza utente:

Link to survey

Contesto​

Requisiti​

Configurazione​

Step 1: Impostare le proprietà del problema da risolvere​

Step 2 (se necessario): Ottimizzare il problema per l'esecuzione su hardware quantistico​

Step 3: Confrontare le prestazioni degli algoritmi​

Step 4: Utilizzare il risultato​

Sondaggio sul tutorial​