Model a flowing non-viscous fluid using QUICK-PDE

Usage estimate: 50 minutes on ibm_kingston. (NOTE: This is an estimate only. Your runtime may vary.)

Note that the execution time of this function is generally greater than 20 minutes, so you might want to split this tutorial in two sections: the first one, in which you read through it and launch the jobs, and the second one a few hours later (giving ample time for the jobs to complete), to work with the results of the jobs.

Background

This tutorial seeks to teach on an introductory level how to use the QUICK-PDE function to solve complex multi-physics problems on 156Q Heron R2 QPUs by using ColibriTD's H-DES (Hybrid Differential Equation Solver). The underlying algorithm is described in the H-DES paper. Note that this solver can also solve non-linear equations.

Multi-physics problems - including fluids dynamics, heat diffusion, and material deformation, to name a few - can be ubiquitously described by Partial Differential Equations (PDEs).

Such problems are highly relevant for various industries and constitute an important branch of applied mathematics. However, solving non-linear multivariate coupled PDEs with classical tools remains challenging due to the requirement of an exponentially large amount of resources.

This function is appropriate for equations with increasing complexity and variables, and is the first step to unlock possibilities that were once considered intractable. To fully describe a problem modeled by PDEs, it is necessary to know the initial and boundary conditions. These can strongly change the solution of the PDE and the path to find their solution.

In this tutorial, you will learn how to:

1. Define the parameters of the initial condition function.
2. Adjust the qubit number (used to encode the function of the differential equation), depth, and shot number.
3. Run QUICK-PDE to solve the underlying differential equation.

Requirements

Before starting this tutorial, be sure you have the following installed:

• Qiskit SDK v2.0 or later (pip install qiskit)
• Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
• Matplotlib (pip install matplotlib)
• Access to the QUICK-PDE function. Fill out the form to request access.

Setup

Authenticate using your API key and select the function as follows:

import numpy as np
import matplotlib.pyplot as plt
from qiskit_ibm_catalog import QiskitFunctionsCatalog
 
catalog = QiskitFunctionsCatalog(
    channel="ibm_cloud / ibm_quantum / ibm_quantum_platform",
    instance="USER_CRN / HGP",
    token="USER_API_KEY / IQP_API_TOKEN",
)
 
quick = catalog.load("colibritd/quick-pde")

Step 1: Set properties of the problem to solve

This tutorial covers the user experience from two perspectives: the physical problem determined by the initial conditions, and the algorithmic component in solving a fluid dynamics example on a quantum computer.

Computational Fluid Dynamics (CFD) has a broad range of applications, and thus it is important to study and solve the underlying PDEs. An important family of PDEs are the Navier-Stokes equations, which are a system of nonlinear partial differential equations describing the motion of fluids. They are highly relevant for scientific problems and engineering applications.

Under certain conditions, the Navier-Stokes equations reduce to Burgers' equation, a convection–diffusion equation describing phenomena occurring in fluid dynamics, gas dynamics, and nonlinear acoustics, to name a few, by modeling dissipative systems.

The one-dimensional version of the equation depends on two variables: $t \in \mathbb{R}_{\geq 0}$ modeling the temporal dimension, $x \in \mathbb{R}$ representing the spatial dimension. The general form of the equation is called the viscous Burgers' equation and reads:

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

where $u(x,t)$ is the fluid speed field at a given position $x$ and time $t$, and $\nu$ is the viscosity of the fluid. Viscosity is an important property of a fluid that measures its rate-dependent resistance to motion or deformation, and thus it plays a crucial role in the determination of the dynamics of a fluid. When the viscosity of the fluid is null ($\nu$ = 0), the equation becomes a conservation equation that can develop discontinuities (shock waves), due to the lack of its internal resistance. In this case, the equation is called the inviscid Burgers' equation and is a special case of nonlinear wave equation.

Strictly speaking, inviscid flows do not occur in nature, but when modeling aerodynamic flow, due to the infinitesimal effect of transport, using an inviscid description of the problem can be useful. Surprisingly, more than 70% of aerodynamic theory deals with inviscid flows.

In this tutorial we offer the inviscid Burgers' equation as a CFD example to solve on IBM® QPUs using QUICK-PDE, given by the equation: $\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = 0,$

The initial condition for this problem is set to a linear function: $u(t=0,x) = ax + b,\text{ with }a,b\in\mathbb{R}$ where $a$ and $b$ are arbitrary constants influencing the shape of the solution. You can adjust $a$ and $b$ and see how they influence the solving process and the solution.

job = quick.run(
    use_case="cfd",
    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  ])}}

Step 2 (if needed): Optimize problem for quantum hardware execution

By default, the solver uses physically-informed parameters, which are initial circuit parameters for a given qubit number and depth from which our solver will start.

The shots are also part of the parameters with a default value, since fine-tuning them is important.

Depending on the configuration you're trying to solve, the algorithm's parameters to achieve satisfactory solutions might need to be adapted; for example, it can require more or fewer qubits per variable $t$ and $x$, depending on $a$ and $b$. The following adjusts the number of qubits per function per variable, the depth per function, and the number of shots.

You can also see how to specify the backend and the execution mode.

In addition, physically-informed parameters might steer the optimization process in a wrong direction; in that case, you can disable it by setting the initialization strategy to "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())

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

Step 3: Use the result

With your solution, you can now choose what to do with it. The following demonstrates how to plot the result.

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:

Notice the difference of initial condition for the second run and its effect on the result:

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:

Tutorial survey

Please take a minute to provide feedback on this tutorial. Your insights will help us improve our content offerings and user experience:

Link to survey