Modele un fluido no viscoso que fluye utilizando QUICK-PDE

Estimación de uso: 50 minutos en un procesador Heron r2. (NOTA: Esto es sólo una estimación. Su tiempo de ejecución puede variar)

Tenga en cuenta que el tiempo de ejecución de esta función suele ser superior a 20 minutos, por lo que es posible que desee dividir este tutorial en dos secciones: la primera, en la que y lanzar los trabajos, y la segunda unas horas más tarde (dando tiempo a que (dando tiempo suficiente para que los trabajos se completen), para trabajar con los resultados de los trabajos.

En segundo plano

Este tutorial tiene como objetivo enseñar, a nivel introductorio, cómo utilizar la función QUICK-PDE para resolver problemas multifísicos complejos en las unidades de procesamiento de cómputo (QPU) Heron R2 de 156Q, mediante el uso del H-DES (solucionador híbrido de ecuaciones diferenciales) de ColibriTD's. El algoritmo subyacente se describe en el artículo sobre H-DES. Ten en cuenta que este solucionador también puede resolver ecuaciones no lineales.

Los problemas multifísicos, como la dinámica de fluidos, la difusión del calor y la deformación de materiales, por nombrar algunos, pueden describirse mediante ecuaciones Diferenciales parciales (EDP).

Estos problemas son muy relevantes para diversas industrias y constituyen una importante rama de las matemáticas aplicadas rama importante de las matemáticas aplicadas. Sin embargo, la resolución de no lineales acopladas multivariantes con herramientas clásicas sigue siendo un reto debido a la necesidad de una cantidad de recursos exponencialmente grande requiere una cantidad exponencialmente grande de recursos.

Esta función es adecuada para ecuaciones con complejidad y variables, y es el primer paso para desbloquear posibilidades que antes se consideraban consideradas intratables. Para describir completamente un problema modelado por EDP, es necesario conocer las condiciones iniciales y de contorno. Estos pueden cambiar la solución de la EDP y el camino para encontrar su solución.

Este tutorial te enseña cómo:

1. Definir los parámetros de la función de condición inicial.
2. Ajuste el número de qubit (utilizado para codificar la función de la ecuación diferencial), la profundidad y el número de disparo.
3. Ejecute QUICK-PDE para resolver la ecuación diferencial subyacente.

Requisitos

Antes de empezar este tutorial, asegúrate de tener instalado lo siguiente:

• Qiskit SDK v2.0 o posterior (pip install qiskit)
• Catálogo de funciones Qiskit (pip install qiskit-ibm-catalog)
• Matplotlib (pip install matplotlib)
• Acceso a la función QUICK-PDE. Rellena el formulario para solicitar acceso.

Configuración

Autentícate utilizando tu clave API y selecciona la función como se indica a continuación:

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

Paso 1: Establecer las propiedades del problema a resolver

Este tutorial abarca la experiencia del usuario desde dos perspectivas: el problema físico determinado por las condiciones iniciales, y el componente algorítmico algorítmico para resolver un ejemplo de dinámica de fluidos en un ordenador cuántico.

La dinámica de fluidos computacional (CFD) tiene una amplia gama de aplicaciones, por lo que es importante estudiar y resolver las EDP subyacentes importante estudiar y resolver las EDP subyacentes. Una familia importante de EDP son las ecuaciones de Navier-Stokes, un sistema de ecuaciones diferenciales parciales no lineales que describen el movimiento de los fluidos diferenciales parciales no lineales que describen el movimiento de los fluidos. Son muy pertinentes para problemas científicos y aplicaciones de ingeniería.

En determinadas condiciones, las ecuaciones de Navier-Stokes se reducen a la ecuación de Burgers de Burgers, una ecuación de convección-difusión que describe fenómenos la dinámica de fluidos, la dinámica de gases y la acústica no lineal, por citar algunos, modelando sistemas disipativos.

La versión unidimensional de la ecuación depende de dos variables: $t \in \mathbb{R}_{\geq 0}$ modelar la dimensión temporal, $x \in \mathbb{R}$ la representación de la dimensión espacial. La forma general de la ecuación se denomina la ecuación viscosa de Burgers y se lee:

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

donde$u(x,t)$ es el campo de velocidad del fluido en una posición dada$x$ y en un instante dado$t$ , y$\nu$ es la viscosidad del fluido. La viscosidad es una propiedad importante de un fluido que mide su resistencia al movimiento o a la deformación, que depende de la velocidad, y, por lo tanto, desempeña un papel crucial en la determinación de la dinámica de un fluido. Cuando la viscosidad del fluido es nula ( $\nu = 0$ ), la ecuación se convierte en una ecuación de conservación que puede dar lugar a discontinuidades (ondas de choque), debido a la ausencia de su resistencia interna. En este caso, la ecuación se denomina «ecuación de Burgers inviscida» y es un caso especial de la ecuación de ondas no lineal.

En sentido estricto, los flujos invisibles no se dan en la naturaleza, pero al modelizar aerodinámico, debido al efecto infinitesimal del transporte, puede ser útil utilizar una descripción no una descripción no viscosa del problema. Sorprendentemente, más del 70% de la teoría aerodinámica trata de flujos no viscosos.

Este tutorial utiliza la ecuación de Burgers para fluidos no viscosos como ejemplo de CFD para resolverla en las QPU de l IBM®, utilizando QUICK-PDE, y viene dada por:

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

La condición inicial para este problema se establece como una función lineal: $u(t=0,x) = ax + b,\text{ with }a,b\in\mathbb{R},$ donde $a$ y $b$ son constantes arbitrarias que influyen en la forma de la solución. Puedes ajustar $a$ y $b$ y comprobar cómo influyen en el proceso de resolución y en la solución. Para esta elección de condición inicial, la ecuación admite la siguiente solución analítica:

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

lo cual se cumple siempre que no se haya formado ninguna onda de choque; es decir, por ejemp $t < -1/a$, cuando $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  ])}}

Paso 2 (si es necesario): Optimizar el problema para la ejecución en hardware cuántico

Por defecto, el solucionador utiliza parámetros informados físicamente, que son los parámetros iniciales del circuito para un número y profundidad de qubit dados, a partir de los cuales comenzará nuestro solucionador.

Los disparos también forman parte de los parámetros con un valor por defecto, ya que afinarlos es importante.

Dependiendo de la configuración que intente resolver, puede que sea necesario adaptar los parámetros del algoritmo para lograr soluciones satisfactorias parámetros del algoritmo para lograr soluciones satisfactorias; por ejemplo, puede puede requerir más o menos qubits por variable $t$ y $x$, dependiendo de $a$ y $b$. A continuación se ajusta el número de qubits por función por variable, la profundidad por función y el número de disparos.

También puede ver cómo especificar el backend y el modo de ejecución.

Además, los parámetros con información física pueden dirigir el proceso de optimización en la dirección equivocada en una dirección errónea; en ese caso, puede desactivarlo configurando la initialization a "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  ])}}

Paso 3: Comparar el rendimiento de los algoritmos

Puede comparar el proceso de convergencia de nuestra solución (HDES) con job_2 el rendimiento de un algoritmo y solucionador de redes neuronales basadas en la física (PINN) (consulte el artículo y el repositorio GitHub asociado).

En el ejemplo de job_2 la salida de (enfoque basado en la cuántica), solo se optimizan 13 parámetros (12 parámetros de circuito más 1 parámetro de escalado) con el solucionador clásico. El proceso de convergencia es el siguiente:

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

Esto significa que se puede 0.0015 alcanzar una pérdida inferior tras 28 iteraciones y optimizando solo unos pocos parámetros clásicos.

Ahora podemos comparar lo mismo con la solución PINN con la configuración predeterminada sugerida por el artículo utilizando un optimizador basado en gradientes. El equivalente a nuestro circuito con 13 parámetros que optimizar es la red neuronal, que requiere al menos ocho capas de 20 neuronas, lo que implica optimizar 3021 parámetros. A continuación, se alcanza la pérdida objetivo en el paso 315, pérdida: 0.0014988397.

Ahora bien, dado que queremos hacer una comparación justa, deberíamos utilizar el mismo optimizador en ambos casos. El número más bajo de iteraciones que encontramos para 12 capas de 20 neuronas = 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

Puede hacer lo mismo con sus datos de job_2 y trazar una comparación con la solución PINN.

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

Paso 4: Utilizar el resultado

Con tu solución, ahora puedes elegir qué hacer con ella. A continuación se muestra cómo trazar el 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:

Obsérvese la diferencia de condición inicial para la segunda ejecución y su efecto en el 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:

Encuesta tutorial

Tómese un minuto para comentar este tutorial. Su opinión nos ayudará a mejorar nuestra oferta de contenidos y la experiencia de los usuarios:

Enlace a la encuesta