QUICK-PDE: una funzione Qiskit di ColibriTD

Consulta la documentazione dell'API

Panoramica

Il risolutore di equazioni differenziali parziali (PDE) qui presentato fa parte della nostra piattaforma Quantum Innovative Computing Kit (QUICK) (QUICK-PDE) ed è confezionato come funzione Qiskit. Con la funzione QUICK-PDE, è possibile risolvere equazioni differenziali parziali specifiche del dominio sulle QPU IBM Quantum. Questa funzione si basa sull'algoritmo descritto nel documento descrittivo H-DES disponibile all'indirizzo ColibriTD's. Questo algoritmo è in grado di risolvere complessi problemi multifisici, a partire dalla fluidodinamica computazionale (CFD) e dalla deformazione dei materiali (MD), con altri casi d'uso in arrivo.

Per affrontare le equazioni differenziali, le soluzioni di prova sono codificate come combinazioni lineari di funzioni ortogonali (tipicamente polinomi di Chebyshev, e più precisamente $2^n$ di essi, dove $n$ è il numero di qubit che codificano la funzione), parametrizzate dagli angoli di un Circuito Quantico Variabile (VQC). L'ansatz genera uno stato che codifica la funzione, valutato da osservabili le cui combinazioni permettono di valutare la funzione in tutti i punti. È quindi possibile valutare la funzione di perdita in cui sono codificate le equazioni differenziali e regolare con precisione gli angoli in un ciclo ibrido, come mostrato di seguito. Le soluzioni di prova si avvicinano gradualmente alle soluzioni reali fino a raggiungere un risultato soddisfacente.

Oltre a questo ciclo ibrido, è possibile concatenare diversi ottimizzatori. Ciò è utile quando si desidera che un ottimizzatore globale trovi un buon insieme di angoli e che un ottimizzatore più preciso segua un gradiente verso il miglior insieme di angoli vicini. Nel caso della fluidodinamica computazionale (CFD), la sequenza di ottimizzazione predefinita produce i risultati migliori; nel caso della deformazione dei materiali (MD), invece, pur fornendo buoni risultati, è possibile configurarla ulteriormente per ottenere vantaggi specifici per il problema.

Si noti che per ogni variabile della funzione viene specificato il numero di qubit (con cui si può giocare). Impilando 10 circuiti identici e valutando le 10 osservabili identiche su qubit diversi in un unico grande circuito, è possibile attenuare il rumore all'interno del processo di ottimizzazione CMA, affidandosi al metodo noise learner, e ridurre significativamente il numero di scatti necessari.

Dinamica dei fluidi computazionale

L'equazione di Burgers per i fluidi non viscosi descrive il flusso di tali fluidi come segue:

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

$u$ rappresenta il campo di velocità del fluido. Questo caso d'uso presenta una condizione al contorno temporale: è possibile selezionare la condizione iniziale e poi lasciare che il sistema si stabilizzi. Attualmente, le uniche condizioni iniziali ammesse sono le funzioni lineari: $ax + b$. La soluzione analitica è:

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

Le equazioni di Eulero senza pressione modellano il flusso di un fluido compressibile e non viscoso con smorzamento come segue:

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

$u\frac{\partial g}{\partial t} + g\frac{\partial u}{\partial t} + u^2\frac{\partial g}{\partial x} + 2gu\frac{\partial u}{\partial x} + \frac{\mu}{(1+t)^{\lambda}} g u = 0,$

$g$ rappresenta il campo di densità, $u$ il campo di velocità e $\mu$ un coefficiente di smorzamento. Nella nostra formulazione, definiamo l' $\lambda = 1$; pertanto, non verrà utilizzata come parametro nel prosieguo. Questo caso d'uso presenta condizioni al contorno temporali: $g(0, x) = e^{-x}$ e $u(0, x) = x$. La soluzione analitica è:

$g(t, x) = \frac{1-\mu}{(1+t)^{1-\mu} - \mu} \exp\!\left(\frac{(\mu-1)\, x}{(1+t)^{1-\mu} - \mu}\right),$

$u(t, x) = \frac{(1-\mu)\, x}{\left((1+t)^{1-\mu} - \mu\right)(1+t)^{\mu}}.$

Le argomentazioni per le equazioni differenziali della CFD sono su una griglia fissa, come segue:

• $t$ è compreso tra 0 e 0.95 con 41 punti di campionamento. $x$ è compreso tra 0 e 0.95 con 41 punti di campionamento.

Deformazione del materiale

Questo caso d'uso si concentra sulla deformazione ipoelastica nell'ambito della prova di trazione unidimensionale, in cui una barra fissata nello spazio viene sottoposta a trazione all'altra estremità. Descriviamo il problema come segue:

$u' - \frac{\sigma}{3K} - \frac{2}{\sqrt{3}}\epsilon_0\left(\frac{\sigma'}{\sigma_0\sqrt{3}}\right)^n = 0,$

$\sigma' - b = 0,$

$K$ rappresenta il modulo di compressibilità del materiale sottoposto a deformazione, $n$ l'esponente di una legge di potenza, $b$ la forza per unità di massa, $\epsilon_0$ il limite di sollecitazione proporzionale, $\sigma_0$ il limite di deformazione proporzionale, $u$ la funzione di sollecitazione e $\sigma$ la funzione di deformazione. La soluzione analitica è:

$\sigma(x) = \sigma_0 - bx,$

$u(x) = -\frac{3^{-(3+n)/2}}{2bK(1+n)\,\sigma_0^{n}}\Biggl[3^{(1+n)/2}b^2\sigma_0^n(1+n)x^2 - 2\cdot 3^{(1+n)/2}b\sigma_0^n(1+n)\sigma_0 x - 12\epsilon_0 K\sigma_0^{1+n}$ $- 12b\epsilon_0 K\sigma_0^n x\left(\frac{-bx+\sigma_0}{\sigma_0}\right)^n + 12\epsilon_0 K\sigma_0^{n+1}\left(\frac{-bx+\sigma_0}{\sigma_0}\right)^n - 12\epsilon_0 K \sigma_0^{1+n}\Biggr],$

dove $\sigma_0 = g(0)$ è la condizione al contorno relativa alla deformazione in $x=0$.

La barra considerata è di lunghezza unitaria. Questo caso d'uso prevede una condizione al contorno per la sollecitazione superficiale $t$, ovvero la quantità di lavoro necessaria per allungare la barra.

Le argomentazioni per le equazioni differenziali di MD sono su una griglia fissa, come segue:

• $x$ è compreso tra 0 e 1 con 30 punti di campionamento.

Benchmark

La tabella seguente presenta i dati statistici relativi a diverse esecuzioni della nostra funzione.

Introduzione

Compila il modulo per richiedere l'accesso alla funzione QUICK-PDE. A questo punto, supponendo che tu abbia già salvato il tuo account nell'ambiente locale, seleziona la funzione come segue:

from qiskit_ibm_catalog import QiskitFunctionsCatalog

catalog = QiskitFunctionsCatalog(
    channel="ibm_cloud / ibm_quantum_platform",
    instance="USER_CRN / HGP",
    token="USER_API_KEY / IQP_API_TOKEN",
)

catalog = QiskitFunctionsCatalog(channel="ibm_quantum_platform")

# Verify that you have access to the function
catalog.list()

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

Esempi

Per iniziare, provate uno dei seguenti esempi:

Equazione di Inviscid Burgers (CFD)

Per l'equazione di Burgers, quando le condizioni iniziali sono impostate su $u(0,x) = x$, i risultati sono i seguenti:

# launch the simulation with initial conditions u(0,x) = a*x + b
job = quick.run(
    use_case="CFD_BURGER", physical_parameters={"a": 1.0, "b": 0.0}
)

Per verificare lo stato del carico di lavoro della funzione Qiskit o per ottenere i risultati, procedere come segue:

# Print the ID so you can use it later, if necessary
print(job.job_id)
print(job.status())
solution = job.result()

import numpy as np
import matplotlib.pyplot as plt


def plot_result_3d(result):
    fig = plt.figure()
    ax = fig.add_subplot(projection="3d")

    t, x = np.meshgrid(result["samples"]["t"], result["samples"]["x"])

    ax.plot_surface(
        t,
        x,
        result["functions"]["u"],
        edgecolor="royalblue",
        lw=0.25,
        rstride=26,
        cstride=26,
        alpha=0.3,
    )
    ax.scatter(t, x, result["functions"]["u"], marker=".")
    ax.set(xlabel="t", ylabel="x", zlabel="u(t,x)")

    plt.show()


# Call
plot_result_3d(solution)

Equazione di Eulero senza pressione (CFD)

Per l'equazione di Eulero, quando le condizioni iniziali sono impostate su $g(0, x) = e^{-x}$ e $u(0, x) = x$, per un dato $\mu$ (in questo caso, $\mu = 0.1$ ) e $\lambda = 1$, i risultati sono i seguenti:

# Launches the solving for an arbitrary mu
job = quick.run(use_case="CFD_EULER", physical_parameters={"mu": 0.1})

solution = job.result()


# Colorplot function
def plot_result_2d(result):
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    configs = {
        "g": {"cmap": "viridis", "title": "g(t, x)"},
        "u": {"cmap": "plasma", "title": "u(t, x)"},
    }

    t = result["samples"]["t"]
    x = result["samples"]["x"]

    for ax, (field, cfg) in zip(axes, configs.items()):
        v = result["functions"][field]

        im = ax.contourf(t, x, v, levels=50, cmap=cfg["cmap"])
        fig.colorbar(im, ax=ax, label=cfg["title"])

        ax.set_xlabel("t")
        ax.set_ylabel("x")
        ax.set_title(cfg["title"], fontsize=13, fontweight="bold")

    plt.tight_layout()
    plt.show()


plot_result_2d(solution)

Deformazione del materiale

Il caso d'uso della deformazione del materiale richiede i parametri fisici del materiale e la forza applicata, come segue:

# Select the properties of your material
job = quick.run(
    use_case="MD",
    physical_parameters={
        "t": 12.0,
        "K": 100.0,
        "n": 4.0,
        "b": 10.0,
        "epsilon_0": 0.1,
        "sigma_0": 5.0,
    },
)

# Plot the result
solution = job.result()

_ = plt.figure()
stress_plot = plt.subplot(211)
plt.plot(solution["samples"]["x"], solution["functions"]["u"])
strain_plot = plt.subplot(212)
plt.plot(solution["samples"]["x"], solution["functions"]["sigma"])

plt.show()

Di seguito è riportato un esempio di come ottenere il valore della funzione per un determinato insieme di coordinate:

# u(t=0.2, x=0.7) == 2
assert solution["samples"]["t"][1] == 0.2
assert solution["samples"]["x"][2] == 0.7
assert solution["functions"]["u"][1, 2] == 2

Recupera i messaggi di errore

Se lo stato del carico di lavoro è ERROR, usare job.error_message() per recuperare il messaggio di errore per aiutare il debug, come segue:

job = quick.run(use_case="MD", physical_params={})

print(job.error_message())


# A wrapper can also be used for a more human readable version
def pprint_error(job):
    print("".join(eval(job.error_message())["error"]))


print("___")
pprint_error(job)

Ottenere supporto

Per assistenza, contattare [email protected].

Passi successivi