QUICK-PDE : une fonction Qiskit par l' ColibriTD

Consultez la documentation de l'API

Aperçu

Le solveur d'équations différentielles partielles (PDE) présenté ici fait partie de notre plateforme Quantum Innovative Computing Kit (QUICK) (QUICK-PDE) et est fourni sous forme de fonction Qiskit. La fonction QUICK-PDE vous permet de résoudre des équations aux dérivées partielles spécifiques à un domaine sur les QPU Quantum IBM. Cette fonction est basée sur l'algorithme décrit dans le document de description H-DES d' ColibriTD's. Cet algorithme permet de résoudre des problèmes multiphysiques complexes, à commencer par la dynamique des fluides computationnelle (CFD) et la déformation des matériaux (MD), et d'autres cas d'utilisation seront bientôt disponibles.

Pour résoudre les équations différentielles, les solutions d'essai sont codées sous forme de combinaisons linéaires de fonctions orthogonales (typiquement des polynômes de Chebyshev, et plus précisément $2^n$ où $n$ est le nombre de qubits codant votre fonction), paramétrées par les angles d'un circuit quantique variable (VQC). L'ansatz génère un état codant la fonction, qui est évalué par des observables dont les combinaisons permettent d'évaluer la fonction en tout point. Vous pouvez ensuite évaluer la fonction de perte dans laquelle les équations différentielles sont encodées, et affiner les angles dans une boucle hybride, comme indiqué ci-dessous. Les solutions expérimentales se rapprochent progressivement des solutions réelles jusqu'à ce que vous obteniez un résultat satisfaisant.

En plus de cette boucle hybride, vous pouvez également enchaîner différents optimiseurs. Ceci est utile lorsque vous souhaitez qu'un optimiseur global trouve un bon ensemble d'angles, puis qu'un optimiseur plus fin suive un gradient jusqu'au meilleur ensemble d'angles voisins. Dans le cas de la dynamique des fluides (CFD), la séquence d'optimisation par défaut produit les meilleurs résultats - mais dans le cas de la déformation des matériaux (MD), bien que la séquence par défaut fournisse de bons résultats, vous pouvez la configurer davantage pour obtenir des avantages spécifiques au problème.

Notez que pour chaque variable de la fonction, nous spécifions le nombre de qubits (avec lequel vous pouvez jouer). En empilant 10 circuits identiques et en évaluant les 10 observables identiques sur différents qubits à travers un grand circuit, vous pouvez atténuer le bruit dans le processus d'optimisation CMA, en vous appuyant sur la méthode d'apprentissage du bruit, et réduire de manière significative le nombre de tirs nécessaires.

Dynamique des fluides computationnelle

L'équation de Burgers pour les fluides non visqueux modélise l'écoulement de ces fluides comme suit :

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

$u$ représente le champ de vitesse du fluide. Ce cas d'utilisation comporte une condition aux limites temporelle : vous pouvez choisir la condition initiale, puis laisser le système se stabiliser. Actuellement, les seules conditions initiales admises sont des fonctions linéaires : $ax + b$. La solution analytique est la suivante :

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

Les équations d'Euler sans pression modélisent l'écoulement d'un fluide compressible et non visqueux avec amortissement comme suit :

$\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$ représente le champ de densité, $u$ le champ de vitesse et $\mu$ un coefficient d'amortissement. Dans notre formulation, nous posons $\lambda = 1$; cette variable ne sera donc pas utilisée comme paramètre dans la suite. Ce cas d'utilisation comporte des conditions aux limites temporelles : $g(0, x) = e^{-x}$ et $u(0, x) = x$. La solution analytique est la suivante :

$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}}.$

Les arguments pour les équations différentielles de la CFD sont sur une grille fixe, comme suit :

• $t$ est comprise entre 0 et 0.95 et comporte 41 points d'échantillonnage. $x$ est comprise entre 0 et 0.95 et comporte 41 points d'échantillonnage.

Déformation du matériau

Ce cas d'utilisation porte sur la déformation hypoélastique dans le cadre d'un essai de traction unidimensionnel, au cours duquel une barre fixée dans l'espace est soumise à une traction à son autre extrémité. Nous décrivons le problème comme suit :

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

$\sigma' - b = 0,$

$K$ représente le module de compressibilité du matériau soumis à un étirement, $n$ l'exposant d'une loi de puissance, $b$ la force par unité de masse, $\epsilon_0$ la limite de contrainte proportionnelle, $\sigma_0$ la limite de déformation proportionnelle, $u$ la fonction de contrainte, et $\sigma$ la fonction de déformation. La solution analytique est la suivante :

$\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],$

où $\sigma_0 = g(0)$ correspond à la condition aux limites relative à la déformation au niveau de $x=0$.

La barre considérée est de longueur unitaire. Ce cas d'utilisation a une condition limite pour la contrainte de surface $t$, ou la quantité de travail nécessaire pour étirer la barre.

Les arguments pour les équations différentielles de MD sont sur une grille fixe, comme suit :

• $x$ est comprise entre 0 et 1 et comporte 30 points d'échantillonnage.

Tests de performances

Le tableau suivant présente des statistiques sur les différents parcours de notre fonction.

Premiers pas

Remplissez le formulaire pour demander l'accès à la fonction QUICK-PDE. Ensuite, en partant du principe que vous avez déjà enregistré votre compte dans votre environnement local, sélectionnez la fonction comme suit :

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

Exemples

Pour commencer, essayez l'un des exemples suivants :

Équation d'Inviscid Burgers (CFD)

Pour l'équation de Burgers, lorsque les conditions initiales sont définies comme suit : $u(0,x) = x$, les résultats sont les suivants :

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

Pour vérifier l'état de votre charge de travail Qiskit Function ou obtenir les résultats, procédez comme suit :

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

Équation d'Euler sans pression (CFD)

Pour l'équation d'Euler, lorsque les conditions initiales sont définies comme suit : $g(0, x) = e^{-x}$ et $u(0, x) = x$, pour une valeur donnée de $\mu$ (ici, $\mu = 0.1$ ) et $\lambda = 1$, les résultats sont les suivants :

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

Déformation du matériau

Le cas d'utilisation de la déformation des matériaux nécessite les paramètres physiques de votre matériau et la force appliquée, comme suit :

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

Voici un exemple illustrant comment obtenir la valeur de la fonction pour un ensemble de coordonnées donné :

# 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

Récupérer les messages d'erreur

Si l'état de votre charge de travail est ERROR, utilisez job.error_message() pour récupérer le message d'erreur afin de faciliter le débogage, comme suit :

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)

Obtenir de l'aide

Pour obtenir de l'aide, contactez [email protected].

Etapes suivantes