Modelar um fluido não viscoso em movimento usando QUICK-PDE

Nota

As Qiskit Functions são um recurso experimental disponível apenas para usuários do IBM Quantum® Premium Plan, Flex Plan e On-Prem (via IBM Quantum Platform API) Plan. Elas estão em status de lançamento de pré-visualização e sujeitas a alterações.

Estimativa de uso: 50 minutos em um processador Heron r2. (NOTA: Esta é apenas uma estimativa. Seu tempo de execução pode variar.)

Observe que o tempo de execução desta função é geralmente superior a 20 minutos, então você pode querer dividir este tutorial em duas seções: a primeira, na qual você lê o tutorial e inicia os jobs, e a segunda algumas horas depois (dando tempo suficiente para os jobs serem concluídos), para trabalhar com os resultados dos jobs.

Contexto

Este tutorial busca ensinar em nível introdutório como usar a função QUICK-PDE para resolver problemas complexos de múltiplas físicas em QPUs Heron R2 de 156Q usando o H-DES (Hybrid Differential Equation Solver) da ColibriTD. O algoritmo subjacente está descrito no artigo do H-DES. Observe que este solucionador também pode resolver equações não lineares.

Problemas de múltiplas físicas - incluindo dinâmica de fluidos, difusão de calor e deformação de materiais, para citar alguns - podem ser descritos de forma ubíqua por Equações Diferenciais Parciais (PDEs).

Tais problemas são altamente relevantes para várias indústrias e constituem um ramo importante da matemática aplicada. No entanto, resolver PDEs acopladas não lineares multivariáveis com ferramentas clássicas permanece desafiador devido ao requisito de uma quantidade exponencialmente grande de recursos.

Esta função é apropriada para equações com complexidade e variáveis crescentes, e é o primeiro passo para desbloquear possibilidades que antes eram consideradas intratáveis. Para descrever completamente um problema modelado por PDEs, é necessário conhecer as condições iniciais e de contorno. Estas podem fortemente alterar a solução da PDE e o caminho para encontrar sua solução.

Este tutorial ensina você a:

1. Definir os parâmetros da função de condição inicial.
2. Ajustar o número de qubits (usados para codificar a função da equação diferencial), profundidade e número de disparos.
3. Executar QUICK-PDE para resolver a equação diferencial subjacente.

Requisitos

Antes de iniciar este tutorial, certifique-se de ter o seguinte instalado:

• Qiskit SDK v2.0 ou posterior (pip install qiskit)
• Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
• Matplotlib (pip install matplotlib)
• Acesso à função QUICK-PDE. Preencha o formulário para solicitar acesso.

Configuração

Autentique usando sua chave de API e selecione a função da seguinte forma:

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

Passo 1: Definir propriedades do problema a resolver

Este tutorial abrange a experiência do usuário sob duas perspectivas: o problema físico determinado pelas condições iniciais e o componente algorítmico na resolução de um exemplo de dinâmica de fluidos em um computador quântico.

A Dinâmica de Fluidos Computacional (CFD) tem uma ampla gama de aplicações e, portanto, é importante estudar e resolver as PDEs subjacentes. Uma família importante de PDEs são as equações de Navier-Stokes, que são um sistema de equações diferenciais parciais não lineares que descrevem o movimento de fluidos. Elas são altamente relevantes para problemas científicos e aplicações de engenharia.

Sob certas condições, as equações de Navier-Stokes se reduzem à equação de Burgers, uma equação de convecção-difusão que descreve fenômenos que ocorrem em dinâmica de fluidos, dinâmica de gases e acústica não linear, para citar alguns, modelando sistemas dissipativos.

A versão unidimensional da equação depende de duas variáveis: $t \in \mathbb{R}_{\geq 0}$ modelando a dimensão temporal, $x \in \mathbb{R}$ representando a dimensão espacial. A forma geral da equação é chamada equação de Burgers viscosa e se lê:

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

onde $u(x,t)$ é o campo de velocidade do fluido em uma dada posição $x$ e tempo $t$, e $\nu$ é a viscosidade do fluido. A viscosidade é uma propriedade importante de um fluido que mede sua resistência dependente da taxa ao movimento ou deformação e, portanto, desempenha um papel crucial na determinação da dinâmica de um fluido. Quando a viscosidade do fluido é nula ($\nu$ = 0), a equação se torna uma equação de conservação que pode desenvolver descontinuidades (ondas de choque), devido à falta de sua resistência interna. Neste caso, a equação é chamada de equação de Burgers não viscosa e é um caso especial de equação de onda não linear.

Estritamente falando, fluxos não viscosos não ocorrem na natureza, mas ao modelar fluxo aerodinâmico, devido ao efeito infinitesimal do transporte, usar uma descrição não viscosa do problema pode ser útil. Surpreendentemente, mais de 70% da teoria aerodinâmica lida com fluxos não viscosos.

Este tutorial usa a equação de Burgers não viscosa como um exemplo de CFD para resolver em QPUs IBM® usando QUICK-PDE, dada pela equação:

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

A condição inicial para este problema é definida como uma função linear: $u(t=0,x) = ax + b,\text{ com }a,b\in\mathbb{R}$ onde $a$ e $b$ são constantes arbitrárias que influenciam a forma da solução. Você pode ajustar $a$ e $b$ e ver como eles influenciam o processo de resolução e a solução.

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

Passo 2 (se necessário): Otimizar problema para execução em hardware quântico

Por padrão, o solucionador usa parâmetros fisicamente informados, que são parâmetros de circuito iniciais para um dado número de qubits e profundidade a partir dos quais nosso solucionador começará.

Os disparos também fazem parte dos parâmetros com um valor padrão, uma vez que ajustá-los é importante.

Dependendo da configuração que você está tentando resolver, os parâmetros do algoritmo para alcançar soluções satisfatórias podem precisar ser adaptados; por exemplo, pode exigir mais ou menos qubits por variável $t$ e $x$, dependendo de $a$ e $b$. O seguinte ajusta o número de qubits por função por variável, a profundidade por função e o número de disparos.

Você também pode ver como especificar o backend e o modo de execução.

Além disso, parâmetros fisicamente informados podem direcionar o processo de otimização em uma direção errada; nesse caso, você pode desabilitá-lo definindo a estratégia de initialization como "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  ])}}

Etapa 3: Compare o desempenho dos algoritmos

Você pode comparar o processo de convergência da nossa solução (HDES) do job_2 com o desempenho de um algoritmo e solver de redes neurais informadas por física (PINN) (veja o artigo e o repositório GitHub associado).

No exemplo da saída do job_2 (abordagem baseada em computação quântica), apenas 13 parâmetros (12 parâmetros de circuito mais 1 parâmetro de escala) são otimizados com o solver clássico. O processo de convergência é o seguinte:

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

Isso significa que uma perda abaixo de 0,0015 pode ser alcançada após 28 iterações, e com a otimização de apenas alguns parâmetros clássicos.

Agora podemos comparar o mesmo com a solução PINN usando a configuração padrão sugerida pelo artigo, utilizando um otimizador baseado em gradiente. O equivalente do nosso circuito com 13 parâmetros a serem otimizados é a rede neural, que requer pelo menos oito camadas de 20 neurônios, e portanto envolve a otimização de 3021 parâmetros. Então, a perda alvo é alcançada na Etapa 315, perda: 0,0014988397.

Agora, já que queremos fazer uma comparação justa, devemos usar o mesmo otimizador em ambos os casos. O menor número de iterações que encontramos para 12 camadas de 20 neurônios = 4701 parâmetros:

(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

Você pode fazer o mesmo com seus dados do job_2, e plotar uma comparação com a solução PINN.

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

Etapa 4: Use o resultado

Com a sua solução, você pode agora escolher o que fazer com ela. O exemplo a seguir demonstra como plotar o resultado.

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

Observe a diferença da condição inicial para a segunda execução e seu efeito no resultado:

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

Pesquisa do tutorial

Por favor, reserve um minuto para fornecer feedback sobre este tutorial. Suas percepções nos ajudarão a melhorar nossas ofertas de conteúdo e experiência do usuário:

Link to survey