Modele um fluido não viscoso em fluxo usando o QUICK-PDE

Estimativa de uso: 50 minutos em um processador Heron r2. (OBSERVAÇÃO: essa é apenas uma estimativa. Seu tempo de execução pode variar)

Observe que o tempo de execução dessa função geralmente é superior a 20 minutos, portanto, talvez seja melhor dividir este tutorial em duas seções: a primeira, na qual você lê e inicia os trabalhos, e a segunda, algumas horas depois você lê e inicia os trabalhos, e a segunda, algumas horas mais tarde (dando tempo suficiente para que os trabalhos sejam concluídos), para trabalhar com os resultados do (dando tempo suficiente para que os trabalhos sejam concluídos), para trabalhar com os resultados dos trabalhos.

Segundo plano

Este tutorial tem como objetivo ensinar, em um nível introdutório, como usar a função QUICK-PDE para resolver problemas multifísicos complexos em processadores de unidades de processamento de quadrantes (QPUs) Heron R2 d 156Q, utilizando o H-DES (Solucionador Híbrido de Equações Diferenciais) da ColibriTD's. O algoritmo subjacente é descrito no artigo sobre o H-DES. Observe que este solucionador também pode resolver equações não lineares.

Problemas multifísicos - incluindo dinâmica de fluidos, difusão de calor e deformação de materiais deformação de materiais, para citar alguns, podem ser descritos de forma onipresente por Parciais (PDEs).

Esses problemas são altamente relevantes para vários setores e constituem um importante ramo da importante ramo da matemática aplicada. No entanto, a solução de PDEs não lineares pDEs acoplados multivariados não lineares com ferramentas clássicas continua sendo um desafio devido à à necessidade de uma quantidade exponencialmente grande de recursos.

Essa função é apropriada para equações com complexidade e variáveis variáveis, 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. Isso pode alterar fortemente a solução do PDE e o caminho para encontrar sua solução.

Este tutorial ensina você a:

1. Defina os parâmetros da função de condição inicial.
2. Ajuste o número do qubit (usado para codificar a função da equação diferencial), a profundidade e o número de disparos.
3. Execute o QUICK-PDE para resolver a equação diferencial subjacente.

Requisitos

Antes de iniciar este tutorial, verifique se você tem os seguintes itens instalados:

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

Instalação

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

import numpy as np
import matplotlib.pyplot as plt
from qiskit_ibm_catalog import QiskitFunctionsCatalog

catalog = QiskitFunctionsCatalog(
    channel="ibm_quantum_platform",
    instance="INSTANCE_CRN",
    # For `token`, use the 44-character API_KEY you created
    # and saved from the IBM Quantum Platform Home dashboard
    token="YOUR_API_KEY",
)

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

Passo 1: Defina as propriedades do problema a ser resolvido

Este tutorial aborda a experiência do usuário a partir de duas perspectivas: o problema físico problema físico determinado pelas condições iniciais e o componente algorítmico na soluçã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 importante estudar e resolver os 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 parciais não lineares que descrevem o movimento de fluidos. Eles 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, por meio da modelagem de sistemas dissipativos.

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

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

onde$u(x,t)$ é o campo de velocidade do fluido em uma determinada posição$x$ e no momento$t$ , e$\nu$ é a viscosidade do fluido. A viscosidade é uma propriedade importante de um fluido que mede sua resistência ao movimento ou à deformação, dependente da taxa de fluxo, 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 apresentar descontinuidades (ondas de choque), devido à ausência de sua resistência interna. Nesse caso, a equação é chamada de equação de Burgers inviscida e constitui um caso especial da equação de ondas não linear.

A rigor, os fluxos inviscidas não ocorrem na natureza, mas ao modelar o fluxo aerodinâmico, devido ao efeito infinitesimal do transporte, usar uma descrição invisível do problema pode ser útil. Surpreendentemente, mais de 70% da teoria aerodinâmica lida com fluxos inviscidas.

Este tutorial utiliza a equação de Burgers para fluidos não viscosos como exemplo de CFD para resolução em QPUs d IBM®, utilizando o QUICK-PDE, definida da seguinte forma:

$\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{ with }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 observar como eles influenciam o processo de resolução e a solução. Para essa escolha de condição inicial, a equação admite a seguinte solução analítica:

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

o que se mantém válido desde que não tenha se formado nenhum choque, ou seja, por e $t < -1/a$, quando $a < 0$.

job = quick.run(
    use_case="CFD_BURGER",
    physical_parameters={"a": 1.0, "b": 1.0},
)

print(job.result())

Output:

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

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

Por padrão, o solucionador usa parâmetros fisicamente informados, que são parâmetros iniciais do circuito para um determinado número de qubit e profundidade a partir dos quais o solucionador será iniciado.

Os disparos também fazem parte dos parâmetros com um valor padrão, pois é importante fazer o ajuste fino.

Dependendo da configuração que estiver tentando resolver, os parâmetros do algoritmo para para obter soluções satisfatórias podem precisar ser adaptados; por exemplo, ele pode exigir mais ou menos qubits por variável $t$ e $x$, dependendo de $a$ e $b$. A seguir, ajustamos 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, os parâmetros informados fisicamente podem direcionar o processo de otimização para uma direção errada em uma direção errada; nesse caso, você pode desativá-lo definindo a initialization para "RANDOM".

job_2 = quick.run(
    use_case="CFD_BURGER",
    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())

Output:

{'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 os desempenhos dos algoritmos

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

No exemplo da saída job_2 de (abordagem baseada em quantum), apenas 13 parâmetros (12 parâmetros de circuito mais 1 parâmetro de escala) são otimizados com o solucionador 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 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 com a configuração padrão sugerida pelo artigo usando um otimizador baseado em gradiente. O equivalente ao 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, como 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 de job_2 e traçar uma comparação com a solução PINN.

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

Passo 4: Use o resultado

Com sua solução, agora você pode escolher o que fazer com ela. A seguir, demonstramos 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()

Output:

Observe a diferença da condição inicial para a segunda execução e seu efeito sobre o resultado 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()

Output:

Pesquisa tutorial

Reserve um minuto para fazer comentários sobre este tutorial. Suas percepções nos ajudarão a melhorar nossas ofertas de conteúdo e a experiência do usuário:

Link para a pesquisa