QUICK-PDE: Uma função Qiskit por ColibriTD

Consulte a referência da API

Visão geral

O solucionador de equações diferenciais parciais (PDE) apresentado aqui faz parte da nossa plataforma Quantum Innovative Computing Kit (QUICK) (QUICK-PDE) e está empacotado como uma função Qiskit. Com a função QUICK-PDE, é possível resolver equações diferenciais parciais específicas de domínios em QPUs do IBM Quantum. Esta função baseia-se no algoritmo descrito no documento de descrição H-DES da ColibriTD's. Este algoritmo pode resolver problemas multifísicos complexos, começando com Dinâmica de Fluidos Computacional (CFD) e Deformação de Materiais (MD), e outros casos de uso em breve.

Para lidar com as equações diferenciais, as soluções experimentais são codificadas como combinações lineares de funções ortogonais (normalmente polinômios de Chebyshev e, mais especificamente, $2^n$ deles, em que $n$ é o número de qubits que codificam sua função), parametrizadas pelos ângulos de um circuito quântico variável (VQC). O ansatz gera um estado que codifica a função, que é avaliada por observáveis cujas combinações permitem a avaliação da função em todos os pontos. Em seguida, é possível avaliar a função de perda na qual as equações diferenciais são codificadas e ajustar os ângulos em um loop híbrido, conforme mostrado a seguir. As soluções de teste se aproximam gradualmente das soluções reais até que você chegue a um resultado satisfatório.

Além desse loop híbrido, você também pode encadear diferentes otimizadores. Isso é útil quando você deseja que um otimizador global encontre um bom conjunto de ângulos e, em seguida, um otimizador mais refinado siga um gradiente até o melhor conjunto de ângulos vizinhos. No caso da dinâmica de fluidos computacional (CFD), a sequência de otimização padrão produz os melhores resultados, mas no caso da deformação de materiais (MD), embora o padrão forneça bons resultados, você pode configurá-lo ainda mais para obter benefícios específicos do problema.

Observe que, para cada variável da função, especificamos o número de qubits (com o qual você pode brincar). Ao empilhar 10 circuitos idênticos e avaliar os 10 observáveis idênticos em qubits diferentes em um circuito grande, é possível mitigar o ruído no processo de otimização do CMA, contando com o método de aprendizagem de ruído, e reduzir significativamente o número de disparos necessários.

Dinâmica de fluidos computacional

A equação de Burgers para fluidos não viscosos modela o fluxo desses fluidos da seguinte forma:

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

$u$ representa o campo de velocidade do fluido. Este caso de uso possui uma condição de contorno temporal: é possível selecionar a condição inicial e, em seguida, permitir que o sistema se estabilize. Atualmente, as únicas condições iniciais aceitas são funções lineares: $ax + b$. A solução analítica é:

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

As equações de Euler sem pressão modelam o fluxo de um fluido compressível e não viscoso com amortecimento da seguinte forma:

$\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$ representa o campo de densidade, $u$ o campo de velocidade e $\mu$ um coeficiente de amortecimento. Em nossa formulação, definimos $\lambda = 1$; portanto, ele não será utilizado como parâmetro a seguir. Este caso de uso possui condições de contorno temporais: $g(0, x) = e^{-x}$ e $u(0, x) = x$. A solução analítica é:

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

Os argumentos para as equações diferenciais da CFD estão em uma grade fixa, como segue:

• $t$ está entre 0 e 0.95, com 41 pontos de amostragem. $x$ está entre 0 e 0.95, com 41 pontos de amostragem.

Deformação do material

Este caso de uso se concentra na deformação hipoelástica por meio do ensaio de tração unidimensional, no qual uma barra fixada no espaço é submetida a tração em sua outra extremidade. Descrevemos o problema da seguinte forma:

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

$\sigma' - b = 0,$

$K$ representa o módulo de volume do material que está sendo distendido, $n$ o expoente de uma lei de potência, $b$ a força por unidade de massa, $\epsilon_0$ o limite de tensão proporcional, $\sigma_0$ o limite de deformação proporcional, $u$ a função de tensão e $\sigma$ a função de deformação. A solução analítica é:

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

onde $\sigma_0 = g(0)$ é a condição de contorno para a deformação em $x=0$.

A barra considerada é de comprimento unitário. Esse caso de uso tem uma condição de limite para a tensão superficial $t$, ou a quantidade de trabalho necessária para esticar a barra.

Os argumentos para as equações diferenciais do MD estão em uma grade fixa, como segue:

• $x$ está entre 0 e 1, com 30 pontos de amostragem.

Referências

A tabela a seguir apresenta estatísticas sobre várias execuções de nossa função.

Introdução

Preencha o formulário para solicitar acesso à função QUICK-PDE. Então, supondo que você já tenha salvo sua conta no seu ambiente local, selecione a função da seguinte maneira:

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

Exemplos

Para começar, experimente um dos exemplos a seguir:

Equação de Burgers para fluidos não viscosos (CFD)

Para a equação de Burgers, quando as condições iniciais são definidas como $u(0,x) = x$, os resultados são os seguintes:

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

Verifique o status da sua carga de trabalho do Qiskit Function ou obtenha os resultados da seguinte maneira:

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

Equação de Euler sem pressão (CFD)

Para a equação de Euler, quando as condições iniciais são definidas como $g(0, x) = e^{-x}$ e $u(0, x) = x$, para um determinado $\mu$ (aqui, $\mu = 0.1$ ) e $\lambda = 1$, os resultados são os seguintes:

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

Deformação do material

O caso de uso de deformação de material requer os parâmetros físicos de seu material e a força aplicada, como 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()

A seguir, apresentamos um exemplo de como obter o valor da função para um conjunto específico de coordenadas:

# 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

Buscar mensagens de erro

Se o status de sua carga de trabalho for ERROR, use job.error_message() para obter a mensagem de erro para ajudar na depuração, como 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)

Obtenha suporte

Para obter suporte, entre em contato com [email protected].

Próximas etapas