Modelar un fluido no viscoso en movimiento usando QUICK-PDE

Nota

Las funciones de Qiskit son una característica experimental disponible únicamente para usuarios del Plan Premium, Plan Flex y Plan On-Prem (a través de la API de IBM Quantum Platform) de IBM Quantum®. Se encuentran en estado de versión preliminar y están sujetas a cambios.

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

Tenga en cuenta que el tiempo de ejecución de esta función es generalmente superior a 20 minutos, por lo que podría querer dividir este tutorial en dos secciones: la primera, en la que usted lo lee y lanza los trabajos, y la segunda unas horas después (dando tiempo suficiente para que los trabajos se completen), para trabajar con los resultados de los trabajos.

Contexto

Este tutorial busca enseñar a nivel introductorio cómo utilizar la función QUICK-PDE para resolver problemas complejos de multi-física en QPUs Heron R2 de 156Q mediante el uso del H-DES (Hybrid Differential Equation Solver) de ColibriTD. El algoritmo subyacente se describe en el artículo de H-DES. Tenga en cuenta que este solucionador también puede resolver ecuaciones no lineales.

Los problemas de multi-física, incluyendo dinámica de fluidos, difusión de calor y deformación de materiales, por nombrar algunos, pueden describirse de manera ubicua mediante Ecuaciones Diferenciales Parciales (EDPs).

Dichos problemas son altamente relevantes para diversas industrias y constituyen una rama importante de las matemáticas aplicadas. Sin embargo, resolver EDPs acopladas multivariadas no lineales con herramientas clásicas sigue siendo un desafío debido a la necesidad de una cantidad exponencialmente grande de recursos.

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

Este tutorial le enseña cómo:

Definir los parámetros de la función de condición inicial.
Ajustar el número de qubits (utilizado para codificar la función de la ecuación diferencial), la profundidad y el número de disparos.
Ejecutar QUICK-PDE para resolver la ecuación diferencial subyacente.

Requisitos

Antes de comenzar este tutorial, asegúrese de tener instalado lo siguiente:

Qiskit SDK v2.0 o posterior (pip install qiskit)
Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
Matplotlib (pip install matplotlib)
Acceso a la función QUICK-PDE. Complete el formulario para solicitar acceso.

Configuración

Autentíquese usando su clave de API y seleccione la función de la siguiente manera:

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

Paso 1: Establecer las propiedades del problema a resolver

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

La Dinámica de Fluidos Computacional (CFD) tiene un amplio rango de aplicaciones, y por lo tanto es importante estudiar y resolver las EDPs subyacentes. Una familia importante de EDPs son las ecuaciones de Navier-Stokes, que son un sistema de ecuaciones diferenciales parciales no lineales que describen el movimiento de los fluidos. Son altamente relevantes para problemas científicos y aplicaciones de ingeniería.

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

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

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

donde $u(x,t)$ es el campo de velocidad del fluido en una posición $x$ y tiempo $t$ dados, y $\nu$ es la viscosidad del fluido. La viscosidad es una propiedad importante de un fluido que mide su resistencia dependiente de la velocidad al movimiento o la deformación, 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 desarrollar discontinuidades (ondas de choque), debido a la falta de su resistencia interna. En este caso, la ecuación se denomina ecuación de Burgers no viscosa y es un caso especial de ecuación de onda no lineal.

Estrictamente hablando, los flujos no viscosos no ocurren en la naturaleza, pero al modelar flujo aerodinámico, debido al efecto infinitesimal del transporte, utilizar una descripción no viscosa del problema puede ser útil. Sorprendentemente, más del 70% de la teoría aerodinámica trata con flujos no viscosos.

Este tutorial utiliza la ecuación de Burgers no viscosa como ejemplo de CFD para resolver en QPUs de IBM® usando QUICK-PDE, dada por la ecuación:

$\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. Puede ajustar $a$ y $b$ y observar cómo influyen en el proceso de resolución y en la solución.

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

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

De manera predeterminada, el solucionador utiliza parámetros informados por la física, que son parámetros iniciales del circuito para un número de qubits y profundidad dados desde los cuales nuestro solucionador comenzará.

Los disparos también forman parte de los parámetros con un valor predeterminado, ya que su ajuste fino es importante.

Dependiendo de la configuración que esté intentando resolver, los parámetros del algoritmo para lograr soluciones satisfactorias podrían necesitar ser adaptados; por ejemplo, puede requerir más o menos qubits por variable $t$ y $x$ , dependiendo de $a$ y $b$ . Lo siguiente 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 informados por la física podrían dirigir el proceso de optimización en una dirección incorrecta; en ese caso, puede desactivarlos estableciendo la estrategia de initialization en "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  ])}}

Paso 3: Comparar el rendimiento de los algoritmos

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

En el ejemplo de la salida de job_2 (enfoque basado en computación cuántica), solo se optimizan 13 parámetros (12 parámetros del circuito más 1 parámetro de escala) 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 alcanzar una pérdida inferior a 0.0015 después de 28 iteraciones, y optimizando solo unos pocos parámetros clásicos.

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

Graph showing PINN data compared with the HDES-Qiskit function.

Ahora, dado que queremos hacer una comparación justa, debemos utilizar el mismo optimizador en ambos casos. El menor número 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]
   20 5.398521572351456e-02 1.0e+00 9.98e-03  1e-02  1e-02 0:02.3
   40 5.444650724530220e-02 1.0e+00 9.97e-03  1e-02  1e-02 0:05.1
   60 4.447407275438309e-02 1.0e+00 9.95e-03  1e-02  1e-02 0:08.2
   80 2.068969979882240e-02 1.0e+00 9.94e-03  1e-02  1e-02 0:11.7
  120 1.028892211616039e-02 1.0e+00 9.91e-03  1e-02  1e-02 0:20.1
  140 5.140972323715687e-03 1.0e+00 9.90e-03  1e-02  1e-02 0:25.4
  180 3.811701666563749e-03 1.0e+00 9.87e-03  1e-02  1e-02 0:37.4
  200 3.189878538250923e-03 1.0e+00 9.85e-03  1e-02  1e-02 0:44.2
  240 2.547040116041899e-03 1.0e+00 9.83e-03  1e-02  1e-02 0:59.7
  280 2.166548743844032e-03 1.0e+00 9.80e-03  1e-02  1e-02 1:18.0
  300 1.783065614290535e-03 1.0e+00 9.79e-03  1e-02  1e-02 1:28.4
  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

Usted puede hacer lo mismo con sus datos de job_2 y graficar 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 su solución, usted puede ahora elegir qué hacer con ella. A continuación se muestra cómo graficar 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 of the previous code cell

Observe la diferencia en la condición inicial de 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 of the previous code cell

Encuesta sobre el tutorial

Por favor, tómese un minuto para proporcionar comentarios sobre este tutorial. Sus opiniones nos ayudarán a mejorar nuestro contenido y la experiencia del usuario:

Link to survey

Contexto​

Requisitos​

Configuración​

Paso 1: Establecer las propiedades del problema a resolver​

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

Paso 3: Comparar el rendimiento de los algoritmos​

Paso 4: Utilizar el resultado​

Encuesta sobre el tutorial​