1. Introduction
In recent years, the field of quantum computing (1) has developed into an active and diverse field of research, and significant progress has been made in a number of important areas. For a relatively small number of applications, quantum algorithms have been developed that provide a significant speedup relative to classical methods. Shor’s algorithm for factoring composite integers and Grover’s algorithm for quantum search were key developments in establishing quantum computing. More recently, significant progress has been made in the area of quantum chemistry and quantum physics. Beyond those two fields, only recently have quantum computing applications appeared in other areas of science and engineering, e.g., work in computational electromagnetics (2, 3), mixing in turbulent flow (4), and computational fluid dynamics (5). More general applications have been developed which take advantage of the unique capabilities of quantum computing platforms, e.g., methods for the solution of linear systems of equations (6) and Poisson equation (7).
In recent years significant progress has been made in designing and constructing quantum computers. Currently available quantum computers are relatively smallscale and have become known as noisy intermediatescale quantum (NISQ) computers. These machines have a limited number of qubits (expected to increase to 50–100 in coming years), a limited connectivity between these qubits, a small set of available quantum gates, and typically very little or no quantum error correction.
This chapter describes results of a recent investigation aiming to assess the potential of quantum computing and suitably designed algorithms for future computational fluid dynamics application, particularly for NISQtype quantum hardware. In this work, the quantum circuit model is used for a “universal” or “digital” quantum computer, i.e., work on adiabatic quantum computing is not considered here. In the absence of the required quantum hardware, largescale parallel simulations on parallel classical computers are required in developing such algorithms. In this work the recently developed quantum simulator (5) included in the MΦC multiphysics CFD framework is used (8, 9).
In the near future, the most likely scenario for the introduction of quantum computing hardware is through the quantum coprocessor model, i.e., where a quantum processing unit (QPU) is loosely coupled to a classical computer with one or more CPUs (10). In current designs, the quantum processor requires storage at low temperatures in a cryostat leading to a distinct physical separation between the classical and quantum hardware. Coupling takes place by exchanging classical information. In application of this hybrid quantum/classical approach, the quantum processor acts like a coprocessor with the quantum processor dealing with selected computationally demanding tasks. The quantum processor receives information from the CPU, and this is used to initialize the quantum state in the quantum processor. During the quantum simulation, the quantum state is transformed by application of quantum gates in quantum circuits. Then measurement operations are used to extract classical information from this quantum state, and this is subsequently passed to the CPU. Since in quantum mechanics a measurement leads to the (partial) collapse of the quantum state, in the hybrid classical/quantum approach, typically multiple realizations of the quantum state are needed to obtain classical information with acceptable levels of noise and uncertainty. It is important to recognize that, since initializing a particular quantum state in quantum computer can be a significant challenge, this hybrid approach can only be expected to lead to significant computational speedups in case the quantum simulation is significantly faster for the selected problem than conventional solution methods.
As an example of this hybrid classical/quantum approach, the author introduced a quantum computing application in which the vortexincell method was used to solve the incompressibleflow NavierStokes equations in a regular domain (5). In this algorithm, the Poisson solvers dominating CPU time requirements are based on the quantum computing equivalent of the fast Fourier transform, i.e., the quantum Fourier transform. In this chapter, this algorithm and its application to example flow problems are investigated further. Specifically, the effect of applying an approximate QFT instead of the full QFT is analyzed for different levels of approximation or truncating of rotation gates in the quantum circuit implementation.
The second part of this chapter describes a more recent investigation into the development of quantum algorithms relevant for computational aerodynamics based on modeling at the kinetic level. The key innovation in these developments is the design targeting execution of the algorithm fully on the quantum processor. In particular, at the start of the simulation, multiple quantum states in the quantum processor would be initialized. Then, the quantum algorithms would perform a series of unitary transformations. Only at the end of the simulation would measurements be performed to extract the required classical. A key question this study aims to answer is for which applications this approach is feasible.
This chapter is organized as follows. Section 2 describes key principles of quantum computing relevant to the quantum algorithms for fluid simulations described here. Section 3 describes the hybrid quantum/classical implementation of the vortexincell method along with a number of example applications. The quantum discretevelocity algorithm for kinetic flow modeling is described in Section 4. Conclusions and future research directions are presented in Section 5.
2. Elements of quantum computing relevant in current work
The fundamental unit of quantum computation is the qubit (1). Whereas a classical bit is confined to existing in either the 0 or 1 state, a qubit can be in a state of superposition, i.e., it exists in both states simultaneously. Upon measuring the qubit, the quantum state collapses to either of these two states, and the qubit is no longer in a state of superposition. The state of a qubit is defined through a pair of complex numbers (1). A collection of
$\mathit{nq}$qubits in a coherent state is termed a quantum register of size
$\mathit{nq}$here. Its quantum state is defined by the wave function Ψ> resulting from the tensor product of the quantum states of each qubit in the coherent register. The superposition in this coherent register then creates 2^{nq} different states that can be found upon measurement of the quantum state. In simulating this quantum state on a classical computer, a storage space of 2^{nq} complex numbers is required.
In the present work, the quantum circuit model of quantum computing is used. In this case, the unitary operations on a quantum state allowed by quantum mechanics are represented by a series of quantum (logic) gates acting on the quantum state. A quantum logic gate is an elementary quantum computing device that performs a fixed unitary operation on selected qubits in a fixed period of time. Written in a matrix form, unitary means that the determinant of the transformation equals one.
2.1 Mapping a computational problem onto the quantum state vector
The quantum state Ψ> for a register with
$\mathit{nq}$coherent qubits is represented by a Hilbert space of dimension 2^{nq}. In a quantum computer, different possibilities exist for the physical implementation of qubits, e.g., the “spin” of an electron (with the two possible states being “spinup” and “spindown”) or the plane of polarization of linearly polarized photon has been used. The discrete energy levels in an atom excited by laser pulses present an alternative to electronspin and photonbased qubit implementations.
We will now describe how this quantum state can be used to represent the storage space required for the computational problems of interest here, representing discretized partially differential equations. As a first step, consider a function f discretized on a (regular) mesh with N mesh points. It follows that
${\text{log}}_{2}N$qubits would suffice to create the required number of degrees of freedom in the quantum state vector.
However, it is important to stress that the quantum state vector only represents the likelihood that upon measurement the quantum state collapses into a particular state (1). In other words, with
$\mathit{nq}={\text{log}}_{2}N$, we cannot extract the full information for all N degrees of freedom using a single realization of this quantum state. However, for as long as this classical information is not needed, the quantum state has the required number of degrees of freedom. This superpositionbased principle for storage of discrete data can be extended to multidimensional problems as well. For example, 24 qubits suffice to store a single discretized function on a 256^{3} regular mesh. In application in which multiple variables are to be stored in each mesh point, as in the discretevelocity method discussed in the second part of this chapter, we add further qubits to the quantum register. Specifically, for each qubit added in the coherent register, the number of degrees of freedom is doubled. Once a mapping of the considered computational problem onto the quantum state vector has been designed, calculations are then performed through application of quantum gates as in the quantum circuit model.
2.2 Approximate quantum Fourier transform (AQFT)
The ability to implement the quantum Fourier transform (QFT) efficiently on a quantum computer is of paramount importance for many quantum algorithms. Figure 1 shows the “standard” quantum circuit implementation of the QFT for an example register with 6 qubits. In Figure 1, “H” represents the onequbit Hadamard gate, and the “Rk” gates are controlled rotation gates over an angle defined by index “k,” i.e.,
$2\pi /{2}^{k}$. In this circuit and all subsequent quantum circuits shown in this chapter, the qubit register is represented vertically, with the leftmost qubit in the register at the top. The horizontal direction determines the sequence of gates that are applied to the quantum state. In the QFT example shown, it can be seen that the qubit indices have been reversed in the output state (righthand side) relative to the input, to represent that this standard QFT circuit returns the discrete Fourier transform in bitreversed ordering.
A particular challenge is presented by the controlled rotation gates particularly those involving small angles. The QFT can be implemented approximately by removing all rotation gates with angles smaller than a certain threshold value, resulting in the approximate QFT (AQFT). In particular for faulttolerant implementations, this is desirable as it greatly reduces the gate count. In the following, we define the approximation or “bandlimiting” in the AQFT as follows. The rotation gates are eliminated above a limit value “k,” i.e., for an angle smaller than 2
$\pi /{2}^{k}$, the rotation gate is not included in the quantum circuit.
3. Hybrid quantum/classical vortexincell method
The vortexincell method is a wellstudied hybrid particlemesh method for incompressible flows and is particularly well suited for flows in regular domains such that efficient Poisson solvers can be used. In the present work, the Fourier analysis approach to solving the problem in a fully periodic domain is used, using the QFT for the required discrete Fourier transforms.
The vortexincell (VIC) method solves the incompressibleflow NavierStokes equations, transformed into the Helmholtz equations for vorticity evolution (5):
$\frac{\mathrm{\partial}\overrightarrow{\mathrm{\omega}}}{\mathrm{\partial}t}+\overrightarrow{u}\frac{\mathrm{\partial}\overrightarrow{\mathrm{\omega}}}{\mathrm{\partial}\overrightarrow{x}}=\left(\overrightarrow{\mathrm{\omega}}.\nabla \right)\overrightarrow{u}+\mathrm{\nu}\mathrm{\Delta}\overrightarrow{\mathrm{\omega}};\overrightarrow{\mathrm{\omega}}=\nabla x\overrightarrow{u}$
E1
In simulations using the VIC, the flow evolves through a (large) number of time steps. During each of these time steps, the velocity field is recomputed using solutions of three Poisson problems for stream function
$\overrightarrow{A}$:
$\mathrm{\u2206}\overrightarrow{A}=\overrightarrow{\mathrm{\omega}};\overrightarrow{u}=\nabla x\overrightarrow{A}$
E2
This part of the VIC is of particular interest here, as it represents the part that would be performed by the quantum processor in the quantum coprocessor model.
Figure 2 shows an example of a VIC simulation of two leapfrogging vortex rings, i.e., flow structures of fundamental importance in fluid mechanics. The lower vortex ring is stronger than the ring above it, and it will therefore convect upward faster, leading to the interaction of the vortex rings as shown. The isosurface represents vorticity strength, i.e., a direct indicator of the “strength” of the considered vortex. Results are compared for two different meshes, 128^{3} and 256^{3}, to highlight the dependency of the solution on the chosen mesh size. Also, in the shown simulation, no quantum errors were simulated, and the full QFT was used. If we now replace the QFT with the AQFT, the results shown in Figures 3 and 4 are obtained. In Figure 3, the “k” limit in the QFT is set to five for both meshes, showing that for the finer mesh this leads to unacceptable errors, while the coarsermesh simulation still produces worthwhile results. If the “k” limit for the finer mesh is increased from 5 to 6, i.e., more controlled rotation gates are included in the AQFT circuit, the simulation on the finer mesh can also be made to produce similarly useful results. These example results show what level of approximation in the QFT is tolerable for application of the VIC method. For other QFTbased CFD solvers, a similar sensitivity study would need to be conducted.
4. Quantum algorithm for discretevelocity method
In computational fluid dynamics, the most widely used methods involve solving the NavierStokes equations for a continuum fluid, i.e., where fluid density, velocity components, and energy in each location in the computational domain are to be found from conservation equations. The vortexincell method used in the previous section employs the NavierStokes equations in a transformed form involving vorticity rather than velocity, but importantly it still uses the continuum flow assumption.
An alternative approach to the NavierStokesbased modeling is a description of the flow at a more detailed level, i.e., at the kinetic level (11). Instead of governing equations for mass, momentum, and energy conservation, the flow is now described by the Boltzmann equation governing a particle distribution function in state space (or 3D velocity space for a 3D monatomic gas flow) for each location in the considered domain (10).
For a monatomic gas, the distribution function
$f\left(\overrightarrow{x,}\overrightarrow{c};t\right)$with
$\overrightarrow{x}$defining the coordinates in (physical) space and
$\overrightarrow{c}$the molecular velocity (defined in “velocity space”) is governed by the Boltzmann equation:
$\frac{\mathrm{\partial f}}{\mathrm{\partial t}}+\overrightarrow{\mathrm{c}}\frac{\mathrm{\partial f}}{\mathrm{\partial}\overrightarrow{\mathrm{x}}}=\mathrm{Q}\left(\mathrm{f},\mathrm{f}\right)$
E3
The distribution function defines for each point
$\overrightarrow{x}$the likelihood that molecules have velocity
$\overrightarrow{c}$. On the righthand side of the equation, the collision term
$Q\left(f,f\right)$represents the effect of collisions between particles (12). For each point in space, we can take moments of the distribution function to obtain the local fluid density, velocity components, and energy.
The key advantage of this approach is that noncontinuum flows, i.e., flows for which the density is so low that we cannot assume it to act as a continuum, can also be modeled. However, the main problem is the large computational cost when using a direct discretization approach due to the high dimensionality, i.e., for a 3D flow problem, we have a sixdimensional solution space (or seven when including time).
A further main challenge is the cost of evaluating the collision term. For the freemolecular flows considered here, the collision term can be discarded, and we will use the collisionless Boltzmann equation instead.
In the discretevelocity method used here, the velocity space is discretized using a uniformly spaced Cartesian mesh. A maximum molecular velocity magnitude is defined,
${c}_{\mathit{max}}$, so that the left and right domain boundaries in velocity space are at
${c}_{\mathit{max}}$and
$+{c}_{\mathit{max}}$, respectively, with a uniform spacing Δc separating each point in velocity space. These limits are problem dependent.
The discretevelocity approach used here has a number of characteristics facilitating a quantum implementation:

A uniformly spaced Cartesian mesh is used for the spatial discretization as well as for the discretization of the velocity space.

In case solid objects are present in the computational domain, these are rectangular, and its edges align with the mesh lines in the mesh. Specifically, solid bodies can be defined by “tagging” selected groups of cells in the mesh.

A constant velocityspace discretization is used in each point in space, i.e., the velocityspace boundaries defined earlier as well as the number of discrete velocities are identical in each cell.

The convection part of the Boltzmann equation (i.e., the second term on the lefthand side in the shown equation) along with gassolid interactions determines the time evolution of the distribution function in the absence of interparticle collisions.

The timeintegration method used here is based on the reservoir technique (13), such that during the time integration the convection step always exactly involves the distribution function defined in a cell of computational mesh to move to a cell that is a nearest neighbor. This is commonly referred to as “streaming” of data.
In the examples shown, problems are restricted to twodimensional flows. For this case, the collisionless Boltzmann equation originally defined for 3D can be reduced to two kinetic equations for two reduced distribution functions:
$\frac{\mathrm{\partial f}}{\mathrm{\partial t}}+\overrightarrow{\mathrm{c}}\frac{\mathrm{\partial f}}{\mathrm{\partial}\overrightarrow{\mathrm{x}}}=0;\frac{\mathrm{\partial g}}{\mathrm{\partial t}}+\overrightarrow{\mathrm{c}}\frac{\mathrm{\partial g}}{\mathrm{\partial}\overrightarrow{\mathrm{x}}}=0$
E4
where the velocity and coordinate vectors are now defined in 2D.
Inspired by the previous work on the Dirac equation (14), an efficient quantum circuit implementation of the streaming operations in x and ydirection on a Cartesian 2D mesh can be created as shown in Figure 5. The example shows how a 12qubit register is used for a discretized function with 6 qubits defining the indices of 64 grid points in x and ycoordinate directions. The circuits with the filled circles define streaming to “right” neighbors, i.e., when each control qubits has the 1> state, the target qubit gets negated (“X” symbol used here). Similarly, the left streaming operation employs multiplecontrol NOT gates with target qubits being negated when each control qubit is in state 0>.
For the streaming operations in quantum discretevelocity method, further extensions are needed. First, the quantum register needs to be extended relative to that shown in Figure 5 to account for the storage of two discretized reduced distributions defined in the discretized velocity space. The additional distribution function is accounted for using an additional qubit termed “g” in the following quantum circuit diagrams. The number of additional qubits needed for the discretevelocity mesh clearly depends on the number of discrete velocities used. For example, for a 16 × 16 discretevelocity mesh, we add 8 qubits (four for each direction in state space). The qubits are denoted by the “u0,”…, “v0,”… qubits in the quantum circuits. Finally, to account for solid objects, we use an additional qubit (“BC” in the diagrams) set to 1> to denote a cell within fluid and 0> for a cell within a solid. For a 64 × 64 Cartesian mesh and a 16 × 16 discretevelocity mesh, Figure 6 shows the quantum circuit used to simulate the freemolecular flow around a rectangular body, for which the evolution of the flow field starting from an initial uniform flow is shown in Figure 7. The key feature in the quantum circuits shown in Figure 6 is the extended number of control qubits, i.e., beyond the checks on the qubits corresponding to spatial coordinates, control qubits also involve the “BC” qubit (fluid/solid flag) as well as the qubits related to the discretevelocity indices. This least feature originates from the need to stream only data associated with selected discrete velocities in the used timeintegration method.
4.1 Quantum circuit implementation of specularreflection boundary conditions
A key aspect of the flows simulated by the quantum algorithm described here are gassolid interactions. In this work, specularreflection boundary conditions are assumed. This means that upon hitting a solid surface, a gas particle will bounce off this surface with the wallnormal velocity component effectively getting reversed and the tangentialflow component being preserved. In Figure 7, the time evolution of a Mach 2 freemolecular flow is shown, starting from an initial flow at uniform velocity everywhere.
As can be seen in Figure 7, the specularreflection boundary conditions gradually make the velocity vectors align with the solid wall such that there is no longer a flow into the solid body. This is the physically correct behavior. In the quantum register, the indexing for the velocitymesh data is designed such that a change of the sign of the discrete velocity implies a bit negation of qubits representing the index of the considered discrete velocity. As an example, Figure 8 shows the quantum circuit implementation of the specular reflection for a rectangular body. In this case, only 16 × 16 discrete velocities are used for clarity. It can be seen that control qubits are used to “select” cells in space for which to apply the boundary condition. A further control qubit involves the “BC” qubit representing the solid/fluid flag. The negation operation (“X”) is applied to the qubits representing the velocitymesh index to create the “change of sign” of the considered discretevelocity data. This circuit is shown as an illustration of the quantum algorithm design approach used here.
4.2 Circuit implementations using ancilla qubits
Although the circuits shown in Figure 6 perform the correct convection operations, an important practical constraint needs to be considered. For the quantum computer implementations achieved so far and those foreseen for the near future, the kind of multiqubitcontrolled NOT operations used here cannot be implemented. A small set of native gates will be available which most likely includes NOT, CNOT, and the Toffoli gate.
The solution around this limitation is the introduction of ancilla qubits, which can be regarded as the quantum equivalence of additional workspace in the memory of a classical computer. Then circuits involving multiqubit operations involving a large number of control gates can be transformed into circuits with more qubits and a larger number of gate operations, however now with a smaller number of control qubits. As is typical in this context, we assume that the ancilla qubits are initially in 0>, and since these are to be reused multiple times, the transformed circuits need to reset the ancilla qubits of this state at the end of the operations. For the quantum circuits implementing streaming in positive xdirection considered previously, Figure 9 shows the required circuit transformations for the cases of 1, 2, or 3 ancilla qubits. Increasing the number of qubits will make the circuits more demanding to implement, so in the actual implementation there is an important tradeoff between gate complexity (i.e., within the set of native gates available) and total number of qubits. For the 64 × 64 twodimensional mesh and 16 × 16 discrete velocities considered in Figure 6, the transformed circuits show that with three ancilla qubits, the maximum number of control gates is reduced to four in a fivequbitcontrolled negation gate. However, for most practical hardware implementations, further transformations would be required.
5. Conclusions
This chapter presented two different quantum algorithms with possible applications in computational fluid dynamics. Beyond their very different areas of application, the key differences are the computational model with regard the quantum coprocessor model of quantum computing. The hybrid quantum/classical algorithm for the vortexincell method involves repeated exchanges of information between classical and quantum hardware, i.e., at each time step in the time integration. In contrast, the quantum algorithm implementing a discretevelocity method for kinetic flow modeling can be performed on the quantum processor for the duration of the simulation, with classical information exchange only required at the start and end of the simulation.
This work addressed a number of key challenges that remain to be investigated further. Firstly, the need for further efficient quantum algorithms as well as a further understanding of how to apply the quantum coprocessor model for this type of flow simulations was investigated. Secondly, the measurementbased extraction of classical information fundamentally changes the way quantum algorithms for CFD application will most likely be used. Finally, obtaining detailed information on the full flow field will be a challenge, so applications for which only certain characteristics of the solution are desired would present a good choice for future applications.
Acknowledgments
A part of the simulation results presented were obtained using the EPSRCfunded ARCHIEWeSt High Performance Computer (
This is a syndicated post. Read the original post at Source link .