# Numerical solution of the incompressible Navier-Stokes equations for
chemical mixers via quantum-inspired Tensor Train Finite Element Method

Egor Kornev  Terra Quantum AG, Kornhausstrasse 25, 9000 St. Gallen,
Switzerland  Sergey Dolgov  Terra Quantum AG, Kornhausstrasse 25, 9000 St.
Gallen, Switzerland  Karan Pinto  Terra Quantum AG, Kornhausstrasse 25, 9000
St. Gallen, Switzerland  Markus Pflitsch  Terra Quantum AG, Kornhausstrasse
25, 9000 St. Gallen, Switzerland  Michael Perelshtein<EMAIL_ADDRESS>Terra Quantum AG, Kornhausstrasse 25, 9000 St. Gallen, Switzerland  Artem
Melnikov  Terra Quantum AG, Kornhausstrasse 25, 9000 St. Gallen, Switzerland

###### Abstract

The solution of computational fluid dynamics problems is one of the most
computationally hard tasks, especially in the case of complex geometries and
turbulent flow regimes. We propose to use Tensor Train (TT) methods, which
possess logarithmic complexity in problem size and have great similarities
with quantum algorithms in the structure of data representation. We develop
the Tensor train Finite Element Method – TetraFEM – and the explicit numerical
scheme for the solution of the incompressible Navier-Stokes equation via
Tensor Trains. We test this approach on the simulation of liquids mixing in a
T-shape mixer, which, to our knowledge, was done for the first time using
tensor methods in such non-trivial geometries. As expected, we achieve
exponential compression in memory of all FEM matrices and demonstrate an
exponential speed-up compared to the conventional FEM implementation on dense
meshes. In addition, we discuss the possibility of extending this method to a
quantum computer to solve more complex problems. This paper is based on work
we conducted for Evonik Industries AG.

##  I INTRODUCTION

The modelling of fluid dynamics is a demanding and challenging field [1, 2,
3]. The partial differential equations (PDEs) describing the physical
processes are complex and in an overwhelming number of practical cases, are
not analytically solvable. This means numerical approaches are required [4].
The main limitation of these approaches is the inevitable trade-off between
computational resources and the quality of the solution. Producing
satisfactory results usually requires very dense discretization, which leads
to enormous amounts of memory and processor runtime being consumed by the
algorithm, especially for transient flows [5]. To overcome these challenges,
alternative methods for approximating flow functions with lower computational
complexity are highly sought.

One such promising and novel approach is based on Tensor Networks. This
mathematical tool was first introduced in the context of multi-particle
quantum physics [6]. Independently, it was proposed for the efficient solution
of computational mathematical problems: basic linear algebra [7], including
solution of linear systems of equations [8, 44], optimization [9, 10], and
machine learning [11, 12, 13]. The main advantage of tensor network algorithms
is an exponential reduction in storage and a poly-logarithmic runtime
complexity by the mesh discretization size compared to conventional methods
[14]. For example, this has been demonstrated on a number of PDE problems [15,
16, 45, 17, 47]. In addition, an important feature of tensor networks is that
they can be quite efficiently mapped to a quantum computer [18].

Unfortunately, all of the problems mentioned above were solved in simple
rectangle-type domains in order for the function values on the mesh to be
well-represented as tensors. The purpose of this work is to apply the tensor
network approach to more complex T-mixer-type domains, similarly to Refs. [19,
20].

Here, we want to emphasize why we don’t use quantum PDE solvers despite the
large variety of them [21, 22, 23]. Quantum algorithms usually solve a PDE
problem by reducing it to the solution of a linear system of equations and
applying the famous HHL [24, 25] or more advanced algorithms, e.g. Child’s
[26]. Although these algorithms can theoretically be applied to solving PDEs,
in practice there are many challenges – the main one is the great complexity
and inefficiency of encoding unitary matrices into a quantum computer [27].
Therefore, basically, very small tasks are solved with the help of these
algorithms, and, as Ref. [28] shows, these methods are very poorly scalable.
In this regard, we adhere to a different approach – we first try to solve the
problem using algorithms on tensor networks, which can be realized on a
classical computer, and then try to extend them to a quantum computer. In
addition to the immediate complexity reduction on a classical computer in
comparison to conventional methods, this approach allows us to better
understand which problems can benefit from hybrid quantum computing [29] and
how.

The structure of the paper is as follows. In Section II, we give an overview
of Tensor Trains and the main algorithms. The main ideas behind the finite
element analysis in application to the problem are presented in Section III.
In Section IV, the model equations and the mathematical problem statement are
presented. Section V provides the results of the numerical tests. In Section
VI, we discuss a possible extension of the presented method to a quantum
computer and Section VII concludes the work.

##  II Tensor Trains

Tensor Trains (TTs) are the basic type of Tensor Network on which all basic
linear algebra operations can be efficiently realized [7]. In this section, we
will give a definition of Tensor Trains and describe the main TT-algorithms
used in our implementation, which are Rounding and MatVec operations [7], as
well as solving a linear system of equations in TT format [8] and TT-cross
approximation [30]. In addition, we consider the Quantized Tensor Train format
in its canonical form [31]. Table 1 contains the most important information
about the operations that are used.

[width = 1 ]TT_vector.pdf

Figure 1: Representation of a $d$-dimensional $n\times\cdots\times n$ tensor
of rank $r$ in a TT format. Shapes with “legs” (TT-cores) represent tensors of
dimension equal to the number of legs. The connection of the tensors along the
leg denotes their convolution by this index.
$\begin{array}[]{clll}\hline\cr\text{ \textnumero }&\text{ Operation }&\text{
Result rank }&\text{ Complexity }\\\ \hline\cr
1&\mathbf{z}=\mathbf{x}\cdot\text{ const
}&\mathrm{r}(\mathbf{z})=\mathrm{r}(\mathbf{x})&O(d\mathrm{r}(\mathbf{x}))\\\
2&\mathbf{z}=\mathbf{x}+\mathbf{y}&\mathrm{r}(\mathbf{z})\leq\mathrm{r}(\mathbf{x})+\mathrm{r}(\mathbf{y})&O\left(nd(\mathrm{r}(\mathbf{x})+\mathrm{r}(\mathbf{y}))^{2}\right)\\\
3&\mathbf{z}=\mathbf{x}\odot\mathbf{y}&\mathrm{r}(\mathbf{z})\leq\mathrm{r}(\mathbf{x})\mathrm{r}(\mathbf{y})&O\left(nd\mathrm{r}^{3}(\mathbf{x})\mathrm{r}^{3}(\mathbf{y})\right)\\\
4&\operatorname{matvec}\mathbf{z}=\mathcal{A}\mathbf{x}&\mathrm{r}(\mathbf{z})\leq\mathrm{r}(\mathcal{A})\mathrm{r}(\mathbf{x})&O\left(nd\mathrm{r}^{3}(\mathcal{A})\mathrm{r}^{3}(\mathbf{x})\right)\\\
5&\operatorname{solve}\mathcal{A}\mathbf{x}=\mathbf{y}&r(\mathbf{x})&O\left(nd\mathrm{r}^{3}(\mathbf{x})\mathrm{r}(\mathcal{A})\right)\\\
6&\mathbf{z}=\operatorname{round}(\mathbf{x},\varepsilon)&\mathrm{r}(\mathbf{z})\leq\mathrm{r}(\mathbf{x})&O\left(nd\mathrm{r}^{3}(\mathbf{x})\right)\\\
\hline\cr\end{array}$ Table 1: $\mathbf{x}$, $\mathbf{y}$ and $\mathbf{z}$ are
tensor train vectors of the same dimensions. Their ranks are
$\mathrm{r}(\mathbf{x})$, $\mathrm{r}(\mathbf{y})$ and
$\mathrm{r}(\mathbf{z})$, respectively. $\mathcal{A}$ is a tensor train matrix
of rank $\mathrm{r}(\mathcal{A})$. (1) - multiplication of a vector by a
constant, (2) - element-wise addition of two vectors, (3) - element-wise
product of two vectors, (4) - matrix-vector multiplication, (5) - solution of
a linear system of equations, (6) - rounding that effectively reduces the rank
of a tensor train with a given precision.

###  II.1 Definition

A Tensor Train can be considered an effective representation of
multidimensional arrays [7]. Such a decomposition is given by

$\displaystyle\mathbf{x}\left(i_{1},\ldots,i_{d}\right)=\sum_{\alpha_{0},\ldots,\alpha_{d-1},\alpha_{d}}$
$\displaystyle
G_{1}\left(\alpha_{0},i_{1},\alpha_{1}\right)G_{2}\left(\alpha_{1},i_{2},\alpha_{2}\right)$
$\displaystyle\ldots G_{d}\left(\alpha_{d-1},i_{d},\alpha_{d}\right),$

where $G_{j}$ is the 3-dimensional tensor called the TT-core. The main
characteristic of such a representation is the rank, $r$, which is equal to
the maximum size among the indices $\alpha_{0},\alpha_{1},\ldots,,\alpha_{d}$
and expresses the correlations and the entanglement in the tensor. In the case
of weak correlations (low rank $r$), this format allows one to store and
perform operations with tensors with logarithmic complexity in their size. The
graphical representation of a tensor train is presented in Fig. 1.

We define a _TT-vector_ as any tensor in the form
$\mathbf{x}\left(i_{1},\ldots,i_{d}\right)$. We implicitly assume that a
vector of length $N=n_{1}\ldots n_{d}$ is treated as a $d$-dimensional tensor
with mode sizes $n_{k}$. This tensor is represented in the TT format as shown
in Fig. 1. Following similar reasoning, we define a _TT-matrix_ as any tensor
of the form $\mathcal{A}\left(i_{1},\ldots,i_{d},j_{1},\ldots,j_{d}\right)$,
where $\left(i_{1},\ldots,i_{d}\right)$ enumerates the rows of $\mathcal{A}$
and $\left(j_{1},\ldots,j_{d}\right)$ enumerates its columns. A graphical
representation of a TT-matrix is presented in Fig. 2 and it reads as
$\mathcal{A}=A_{1}\left(i_{1},j_{1}\right)\ldots
A_{d}\left(i_{d},j_{d}\right)$.

[width = 1 ]TT_matrix.pdf

Figure 2: Graphical representation of a TT-Matrix $\mathcal{A}$. The blue
shapes with legs denote TT-cores. $\left(i_{1},\ldots,i_{d}\right)$ enumerates
the rows of $\mathcal{A}$ and $\left(j_{1},\ldots,j_{d}\right)$ enumerates its
columns.

###  II.2 Algorithms

Basic linear algebra operations, such as multiplication of a TT-vector by a
constant, addition of two TT-vectors, and element-wise multiplication, can be
effectively implemented in TT format (the information is provided in the Table
1) [7]. As one can see, after applying these operations, the rank of the
result can increase significantly, so it needs to be effectively lowered after
each step - the Rounding operation serves this purpose. The Rounding operation
allows one to efficiently reduce the rank of a given TT-vector $\mathbf{x}$
with a specified accuracy. The complexity of this operation is
$O(dnr(\mathbf{x})^{3})$, where $r(\mathbf{x})$ is the initial TT-rank of
$\mathbf{x}$.

The MatVec operation is a multiplication of a TT-vector, $\mathbf{x}$, by a
TT-matrix, $\mathcal{A}$, which gives the TT-vector,
$\mathbf{y}=\mathcal{A}\mathbf{x}$ [7]. The straightforward implementation,
where one collects TT-cores of $\mathbf{y}$ just as a product of $\mathbf{x}$
and $\mathcal{A}$ TT-cores and then performs Rounding to decrease the ranks
has complexity $O(dnr(\mathbf{x})^{3}r(\mathcal{A})^{3})$.

Another important operation we used is the solution of a linear system of
equations $\mathcal{A}\mathbf{x}=\mathbf{b}$ via the Alternating Minimal
Energy (AMEn) algorithm [32]. The AMEn algorithm is a combination of the
single-block Density Matrix Renormalization Group (DMRG) [33] and the steepest
descent [34] algorithms. At each step, the local system of equations is
solved, like in DMRG, but also a corresponding TT-core of the residual
$\mathbf{z}=\mathcal{A}\mathbf{x}-\mathbf{b}$ is efficiently added. This
allows one to vary the ranks of the solution and get rid of the stagnation in
the local minima that are inherent in DMRG. The AMEn algorithm has complexity
$O(dnr(\mathbf{x})^{3}r(\mathcal{A})+dnr(\mathbf{x})^{2}r(\mathcal{A})^{2})$.

Finally, the algorithm that was not included in Table 1 is TT-cross
approximation [30]. It enables one to approximate a wide range of functions in
the TT format with only $O(dnr^{2})$ function evaluations. We use this
algorithm to represent the Jacobian of the domain transformation and initial
conditions in the TT format.

##  III TENSOR-TRAIN FINITE ELEMENT METHOD (TetraFEM)

In this section, we develop our proposed method of using the TT decomposition
for the solution of PDEs. We use a finite element formulation because it is
easier to handle complex areas and domain meshing.

###  III.1 Finite element discretization of PDEs

The main idea of the finite element method (FEM) is to approximate a target
function, $u(x,y)$, as a sum of $N$ known basis functions
$\\{\varphi_{i}(x,y)\\}_{i=1}^{N}$ each multiplied by coefficients
$\\{\bar{u}_{i}\\}_{i=1}^{N}$,

$u(x,y)\approx u_{N}(x,y):=\sum_{i=1}^{N}\bar{u}_{i}\varphi_{i}(x,y).$ (1)

These basis functions are non-zero only on a small part of the domain, are
differentiable and are of the same type (such as piecewise polynomial).
Usually, these functions are linked with the discretization mesh, for example,

$\varphi_{i}(x,y)=\left\\{\begin{array}[]{rl}1,&\mbox{$(x,y)$ is the $i$-th
mesh node},\\\ 0,&\mbox{$(x,y)$ is the $j$-th mesh node with $j\neq i$,}\\\
&\mbox{polynomial$(x,y)$ between mesh nodes.}\end{array}\right.$

In this paper, we use basis functions that are linear between the mesh nodes,
as shown in Fig. 3.

Partial derivatives of the discrete function follow trivially from Eq. (1):

$\displaystyle\nabla u(x,y)=\sum_{i=1}^{N}\bar{u}_{i}\nabla\varphi_{i}(x,y).$

[width=1]base-functions-overlap-3.pdf

Figure 3: An example of linear basis functions for quadrilateral finite
elements in the 2D case. Each function equals one on its node and zero on
every other node, which allows the coefficients $\bar{u}$ in Eq. (1) to be
interpreted as values of the target function at the nodes. Also, each function
is non-zero only in a small area, which results in the sparseness of the
matrices in Eq. (4) and Eq. (5).

Let us consider a model PDE problem

$\displaystyle-\Delta u$ $\displaystyle=f$ $\displaystyle\mbox{on a domain
}D\subset\mathbb{R}^{2},$ $\displaystyle u$ $\displaystyle=0$
$\displaystyle\mbox{on the boundary }\partial D,$

where $f\in L^{2}(D)$ is a forcing term. The solution $u\in H^{1}_{0}(D)$ can
be sought in the weak formulation

$\iint_{D}(\nabla u(x,y),\nabla\phi(x,y))dxdy=\iint_{D}f(x,y)\phi(x,y)dxdy,$
(2)

for all $\phi\in H^{1}_{0}(D)$. After substituting Eq. (1) into Eq. (2) and
since the basis functions are known, they can be integrated, which results in
a linear system of equations:

$\mathbf{S}\bar{u}=\mathbf{M}\bar{f},$ (3)

where $\bar{f}=\\{\bar{f}_{i}\\}_{i=1}^{N}$ is a vector of coefficients of
$f(x,y)$ in the finite element basis and
$\mathbf{S},\mathbf{M}\in\mathbb{R}^{N\times N}$ are called the stiffness and
mass matrices, respectively. They are populated with elements:

$\displaystyle\mathbf{S}_{ij}=\iint(\nabla\varphi_{i},\nabla\varphi_{j})dxdy,\quad\mathbf{M}_{ij}=\iint\varphi_{i}\varphi_{j}dxdy.$
(4)

Matrices discretizing the first derivatives can be computed similarly:

$\displaystyle(\mathbf{D_{\mu}})_{ij}=\iint(\varphi_{i},\frac{\partial\varphi_{j}}{\partial\mu})dxdy,\quad\mu\in\\{x,y\\}.$
(5)

These matrices are sparse, which follows from the fact that each basis
function is non-zero only in a small area, so the integrals above are also
non-zero only for the neighboring basis functions. This allows us to solve Eq.
(3) efficiently, for example, using iterative methods. Moreover, as we will
see, such matrices are effectively represented in TT format.

###  III.2 Quantized Tensor Train

[width = 1 ]qtt_canonical.jpg

Figure 4: We discretize a two-dimensional function $f(x,y)$ on a uniform grid
and obtain the corresponding tensor $\bar{f}$. We represent this tensor in a
Tensor Train (TT) format and then further compress it into a Quantized Tensor
Train (QTT) format, where each visible index has dimension 2.

The functions we consider in this paper seemingly depend on only two
variables, $x$ and $y$. To compress coefficient vectors, such as $\bar{u}$ or
$\bar{f}$, we use the so-called Quantized Tensor Train (QTT) format [31]. It
is obtained by splitting the index $i=1,\ldots,N$ into artificial indices of
smaller range,

$i=1+(i_{1}-1)+2(i_{2}-1)+\cdots+2^{d-1}(i_{d}-1),$

where $d=\log_{2}N$, such that $i_{k}\in\\{1,2\\}$ and $k=1,\ldots,d$. A
coefficient vector, for example, $\bar{f}$, can now be reshaped into a tensor
with elements $\bar{f}(i_{1},\ldots,i_{d})$. As shown in Fig. 4, this tensor
can be approximated in the TT format as before but the fact that the indices
$i_{1},\ldots,i_{d}$ “quantize” the original index $i$ is referred to by
naming the resulting representation the QTT format of $\bar{f}$. The utility
of the QTT format stems from the fact that the FEM expansion coefficients of
many elementary functions (such as trigonometric, exponential, polynomial and
rational functions) admit rapidly converging, or even exact, QTT
decompositions [35].

[width=1]scheme.pdf

Figure 5: Scheme of the full simulation process via TetraFEM. At first, the
geometry is specified as a set of adjacent subdomains (quadrilaterals with
analytically-defined boundaries). Then, each of the subdomains is transformed
into a rectangle and the mesh is generated. After that, the Jacobian of the
transformation is sampled on the mesh and local finite element matrices (in
each of the subdomains) are assembled in the Tensor Train format. The global
stiffness, mass and partial derivatives matrices are completed as a
concatenation of local ones. Finally, after setting the initial and boundary
conditions, the time-stepping can be performed as a sequence of TT-operations.

In addition, it is experimentally shown that the finite element matrices have
low QTT ranks, and therefore are well-compressible (Fig. 10(a)) [36]. This
means that all matrix-vector operations from Table 1 will have lower
computational complexity than operations with full (even sparse) matrices.

[width=1]transform-2.pdf

Figure 6: The projection of a reference square onto the domain of Eq. (6). The
curved sides of the domain are described with four differentiable analytic
functions.

###  III.3 Domain geometry and mesh

We consider the domain $D$ which consists of $S$ arbitrary curvilinear
quadrilaterals that are connected by their sides as shown in Fig. 7(a). Each
subdomain is specified by a set of four parametrized smooth curves:

$\displaystyle\mathbf{f}_{k}:t\longrightarrow(x,y),\quad t\in[-1,1],\quad
k=1,\dots 4,$
$\displaystyle\mathbf{f}_{0}(-1)=\mathbf{f}_{3}(-1),\quad\mathbf{f}_{1}(-1)=\mathbf{f}_{0}(1),$
$\displaystyle\mathbf{f}_{2}(1)=\mathbf{f}_{1}(1),\quad\mathbf{f}_{3}(1)=\mathbf{f}_{2}(-1).$

The last four equations define the directions of parametrizations and
guarantee that the sides of a quad coincide at the corners. After that, the
differentiable mapping from a reference $[-1,1]\times[-1,1]$ square to a quad
can be constructed:

$\mathbf{N}(\xi,\eta)=-\frac{(1-\xi)(1-\eta)}{4}\mathbf{f}_{0}(-1)-\\\
\frac{(1+\xi)(1-\eta)}{4}\mathbf{f}_{1}(-1)-\frac{(1+\xi)(1+\eta)}{4}\mathbf{f}_{2}(1)-\\\
\frac{(1-\xi)(1+\eta)}{4}\mathbf{f}_{3}(1)+\frac{(1-\eta)}{2}\mathbf{f}_{0}(\xi)+\\\
\frac{(1+\xi)}{2}\mathbf{f}_{1}(\eta)+\frac{(1+\eta)}{2}\mathbf{f}_{2}(\xi)+\frac{(1-\xi)}{2}\mathbf{f}_{3}(\eta).$
(6)

This technique is also known as transfinite interpolation [37]. The Jacobian
of this transformation can be computed as follows:

$\displaystyle J=\begin{bmatrix}J_{11}&J_{12}\\\
J_{21}&J_{22}\end{bmatrix}=\begin{bmatrix}\frac{\partial
x}{\partial\xi}&\frac{\partial x}{\partial\eta}\\\ \frac{\partial
y}{\partial\xi}&\frac{\partial
y}{\partial\eta}\end{bmatrix}=\begin{bmatrix}\frac{\partial\mathbf{N}(\xi,\eta)}{\partial\xi}&\frac{\partial\mathbf{N}(\xi,\eta)}{\partial\eta}\end{bmatrix}.$
(7)

To set the finite element mesh (i.e., the positions of the basis functions), a
uniform Cartesian mesh of size $2^{d/2}\times 2^{d/2}$ is introduced on the
reference square and then projected onto each quadrilateral with Eq. (6) as
shown in Fig. 6. Using this and also the analytic formula for the Jacobian
(Eq. (7)), finite element matrices from Eq. (4) and Eq. (5) are constructed
for each of the $S$ subdomains. So Eq. (4) becomes

$\mathbf{S}_{ij}=\iint\limits_{[-1,1]\times[-1,1]}\Big{(}(J^{T})^{-1}\nabla\Phi_{i},(J^{T})^{-1}\nabla\Phi_{j}\Big{)}|J|dxdy,$
(8)

where $\Phi$ is a $Q1$ basis function on the reference square and $J$ is the
Jacobian. The numerical integration is performed by standard $2\times 2$
Gauss–Legendre quadrature [38]. The matrices $\mathbf{M},\mathbf{D_{x}}$ and
$\mathbf{D_{y}}$ are assembled analogously. To represent the values of the
Jacobian for the arbitrary $\xi$ and $\eta$, the cross-approximation approach
is utilized [30].

This results in the adjacent subdomain meshes having common nodes on the inner
boundaries which can be seen in Fig. 7(b). The coefficient vectors
$\bar{f}_{1},\ldots,\bar{f}_{S}$, corresponding to different subdomains, are
enumerated by an extra index $s=1,\ldots,S$. This allows us to store them all
in a tensor with elements $\bar{f}(i_{1},\ldots,i_{d},s)$ and approximate it
in a $(d+1)$-dimensional QTT format.

###  III.4 Assembly of global FEM matrices

During the assembly, local finite element matrices are generated on each of
the $S$ quadrilateral subdomains. Then, the global matrices are obtained by
so-called subdomain stitching. Let’s denote one of the
$\mathbf{M},\mathbf{S},\mathbf{D_{x}},\mathbf{D_{y}}$ matrices as
$\mathbf{A}$, because the stitching procedure is identical for every type of
matrix.

The idea of stitching is to unite the local basis functions on the borders of
the adjoined subdomains into continuous functions in the global block matrix
by handling the indices during the construction of blocks.

The diagonal block of the global $\mathbf{A}$ is defined as follows:

$\displaystyle\mathbf{A}^{kk}=\mathbf{A}^{k}+c\mathbf{P}^{k},$

where $\mathbf{A}^{k}$ is the local matrix of the $k$-th subdomain,
$\mathbf{P}^{k}$ is the diagonal matrix with ones on the positions
corresponding to the indices of the repeating basis functions and $c$ is the
maximum absolute value across all the main diagonals of $\mathbf{A}^{k}$ for
every $k$.

The construction of non-diagonal blocks is a little more complicated:

$\displaystyle\mathbf{A}^{kl}=\mathbf{P}^{kl}\mathbf{A}^{l}-c\mathbf{P}^{kl}.$

Here, $\mathbf{P}^{kl}_{ij}=1$ if the $i$-th basis function of the $k$-th
subdomain globally coincides with the $j$-th function of $l$-th subdomain and
is zero otherwise. The values are being “pulled” from the neighbor subdomain
matrix.

The final assembly of the global matrix is performed via Kronecker product in
the Tensor Train format:

$\displaystyle\mathbf{A}=\sum_{k,l=1}^{S}\mathbf{A}^{kl}\otimes\mathbf{I}^{kl},$

where $\mathbf{I}^{kl}$ is the $S\times S$ matrix with the only non-zero
element $\mathbf{I}^{kl}_{kl}=1$.

##  IV GOVERNING EQUATIONS AND NUMERICAL SCHEME

###  IV.1 General numerical scheme

To model the flow, the incompressible Navier-Stokes equations are considered:

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}+\boldsymbol{u}\cdot\nabla\boldsymbol{u}=-\frac{1}{\rho}\nabla
p+\nu\nabla^{2}\boldsymbol{u},$ $\displaystyle\nabla\cdot\boldsymbol{u}=0,$

where $\boldsymbol{u}=[u\quad v]^{T}$ stands for the X- and Y-components of
the velocity, $\rho$ is the constant fluid density, $p$ is the pressure and
$\nu$ is the kinematic viscosity. The first equation is called the momentum
equation and the second is the continuity equation.

The possible boundary conditions include:

  1. 1.

_Inlet._

$\boldsymbol{u}=\boldsymbol{u}_{in}$

This boundary condition implies that the inflow is known from the problem
statement and therefore, the velocity on the boundary remains constant for
each time step. The velocity can be constant across the boundary or parabolic,
which is used to eliminate unrealistic gradients at the corners:

$u_{in}(0,y)=u_{0}(1-(1-2y)^{2}),\quad v_{in}(0,y)=0.$

  2. 2.

_Outlet._

$\nu\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{n}}-\boldsymbol{n}p=\boldsymbol{0}$

The outflow means that the fluid is leaving the domain unconstrained. This
approach is called the “do-nothing boundary condition” and is widely used in
computational fluid dynamics.

  3. 3.

_No-slip._

$\boldsymbol{u}=\boldsymbol{0}$

This boundary condition means that the walls are stationary and “sticky.” In
some cases, the wall can be considered “moving” and then the velocity on this
wall has some fixed non-zero value.

  4. 4.

_Natural._

$\frac{\partial p}{\partial\boldsymbol{n}}=0$

The natural boundary condition is used when solving Poisson’s equation for the
pressure. After inspecting the momentum equation, one can notice that a
constant can be added to the pressure without affecting the result.

To impose the boundary conditions, a $\mathbf{Mask}$ matrix is used:

$\mathbf{Mask}_{i,j}=\begin{cases}1&\text{if }i=j\text{ and $\varphi_{i}$ is
interior or outlet}\\\ 0&\text{else}\end{cases}$

The $\mathbf{Mask}$ is diagonal so that after multiplying the solution by the
mask, the essential boundary values become zero. Subsequently, any values can
be inscribed by adding a vector with boundary values (e.g., $u_{in}$ and
$v_{in}$) at the elements corresponding to the essential boundary nodes, and
zeros elsewhere.

For the numerical solution, we chose Chorin’s projection method [39] – the
temporal discretization is explicit Euler:

$\frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\Delta t}=R(t^{n}),$

where $t^{n}=n\Delta t$ is the $n$-th point in time and $R(t)$ is the right
hand side of the momentum equation.

The update of the solution on each time step is performed in three substeps:

  1. 1.

_The predictor step._ The intermediate velocity, $\boldsymbol{u}^{*}$, is
computed without the pressure taken into account:

$\frac{\boldsymbol{u}^{*}-\boldsymbol{u}^{n}}{\Delta
t}=-\boldsymbol{u}^{n}\cdot\nabla\boldsymbol{u}^{n}+\nu\nabla^{2}\boldsymbol{u}^{n}$

  2. 2.

_The pressure step._ Here, Poisson’s equation is solved for the pressure.

$\nabla^{2}p^{n+1}=\frac{\rho}{\Delta t}\nabla\cdot\boldsymbol{u}^{*}$

  3. 3.

_The corrector step._ In this step, the pressure term is finally added to the
equation.

$\frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{*}}{\Delta
t}=-\frac{1}{\rho}\nabla p^{n+1}$

It can be observed that the sum of the predictor and corrector steps results
in the momentum equation. For the convergence of the scheme, it is necessary
to ensure two conditions on the time step. The first one is the
Courant–Friedrichs–Lewy (CFL) condition for the convection term
$\boldsymbol{u}\cdot\nabla\boldsymbol{u}$:

$\Delta t_{1}<\frac{C_{1}\Delta x}{|u|_{\max}},$

where $C_{1}<1$ is the CFL number. For the diffusion term,
$\nu\nabla^{2}\boldsymbol{u}$, the condition is

$\Delta t_{2}<\frac{\Delta x^{2}}{\nu}.$

So the time step $\Delta t$ is chosen such that $\Delta t\leq\min(\Delta
t_{1},\Delta t_{2})$.

[width=1]domain.pdf

(a) The T-mixer geometry, representing an intersection of two straight
microchannels.

[width=1]mesh.pdf

(b) Computational mesh (centers of basis functions).

Figure 7: An example of a computational domain. (a) depicts the geometry
specified as a set of connected quadrilaterals and (b) shows the domain after
assigning the mesh for each subdomain.

###  IV.2 Realization of numerical scheme using Tensor Train operations

After assembling all the matrices, the scheme can be realized as shown in
Algorithm 1.

From this point on, we will be considering the X and Y components of the
velocity $\mathbf{u}=[u,v]$ from Section IV separately as two TT-vectors of
discrete coefficients, $u(i_{1},\ldots,i_{d})=\bar{u}_{i}$ (as defined in Eq.
(1)), and $v(i_{1},\ldots,i_{d})=\bar{v}_{i}$; the pressure $p$ is represented
analogously.

As long as only a limited amount of operations are needed for the basic method
(such as the element-wise sum and product of two vectors, matrix-vector
multiplication and the solving of a linear system), which are all realized in
tensor-train format (Table 1) – after assembling the TT-finite element
matrices, the algorithm could be rewritten as a sequence of operations over
TT-vectors. The notation is as follows: $u$ and $v$ stand for the x- and
y-components of the velocity; $\varepsilon$ is the relative error of the TT-
rounding; TT-matrices are shown in bold; TT-vectors are not in bold.

The solution of linear systems is performed with the AMEn algorithm [32],
which is implemented in a function $amen\\_solve$ in the ttpy package [40].

Algorithm 1 Predictor-corrector time stepping

1:$u=u_{in}.copy()$

2:$v=v_{in}.copy()$

3:for $it=1,nt$ do

4: $\triangleright$ Compute numerical derivatives

5: $Su=solve(\mathbf{M},matvec(\mathbf{S},u)).round(\varepsilon)$

6: $Sv=solve(\mathbf{M},matvec(\mathbf{S},v)).round(\varepsilon)$

7: $Dxu=solve(\mathbf{M},matvec(\mathbf{D_{x}},u)).round(\varepsilon)$

8: $Dxv=solve(\mathbf{M},matvec(\mathbf{D_{x}},v)).round(\varepsilon)$

9: $Dyu=solve(\mathbf{M},matvec(\mathbf{D_{y}},u)).round(\varepsilon)$

10: $Dyv=solve(\mathbf{M},matvec(\mathbf{D_{y}},v)).round(\varepsilon)$

11: $\triangleright$ Predictor step

12: $u^{*}=u+\Delta t(\nu Su-u\odot Dxu-v\odot Dyu)).round(\varepsilon)$

13: $v^{*}=v+\Delta t(\nu Sv-u\odot Dxv-v\odot Dyv)).round(\varepsilon)$

14: $\triangleright$ Pressure-Poisson equation

15: $Dxu=matvec(\mathbf{D_{x}},u^{*}$)

16: $Dyv=matvec(\mathbf{D_{y}},v^{*}$)

17: $f=\frac{\rho}{\Delta t}(Dxu+Dyv).round(\varepsilon)$

18: $p=solve(\mathbf{S},f)$

19: $\triangleright$ Corrector step

20: $Dxp=solve(\mathbf{M},matvec(\mathbf{D_{x}},p)).round(\varepsilon)$

21: $Dyp=solve(\mathbf{M},matvec(\mathbf{D_{y}},p)).round(\varepsilon)$

22: $u=u^{*}-\frac{\Delta t}{\rho}Dxp$

23: $v=v^{*}-\frac{\Delta t}{\rho}Dyp$

24: $\triangleright$ Apply boundary conditions and update the values

25: $u=matvec(\mathbf{Mask},u)+u_{in}$

26: $v=matvec(\mathbf{Mask},v)+v_{in}$

27:end for

The frequent TT-rounding is used to prevent uncontrolled growth of the ranks
and to keep the operation asymptotics low (see Table 1). The fine tuning of
$\varepsilon$ allows us to control the accuracy of the solution and establish
the optimal trade-off between the error and the computational time.

##  V Numerical experiment – T-Mixer problem

The main problem that we solve here is a mixing of liquids in a rectangular
T-mixer. We took the problem statement and the geometry parameters of the
T-mixer (T-junction) from Ref. [41]:

Parameter | Value
---|---
Inlet width, µm | $265$
Inlet channel length, µm | $900$
Outlet width, µm | $300$
Outlet channel length, µm | $135$

Here, two flows of liquid are entering the junction by the inlets with such
velocities that the flow rates are $4$ ml/h at each channel. We consider only
the case of a mixture of a single liquid – pure water at ambient temperature
with standard parameters. The Reynolds number for this problem is estimated to
be about $Re\approx 10$.

To solve this problem via TetraFEM, we divide the domain into four rectangles:
one for the outlet section, two for the inlet sections, and one for the cross-
section which are aligned as shown in Fig. 7(a). To test the method, we set
the rounding tolerance to $\varepsilon=10^{-3}$ and the discretization
parameter to $d=10$, which means the total number of nodes of the space grid
in one subdomain is $2^{5}\times 2^{5}$. Figs. 8(a) and 8(b) depict the
velocity fields on the outlet in the reached steady state. Fig. 9 shows a
comparison between FEM (numpy realization), TetraFEM and experimental data
[41] for the velocity profiles at the outlet: the curves are almost the same.

[width=1]u.pdf

(a) X component of the velocity.

[width=1]v.pdf

(b) Y component of the velocity.

Figure 8: Values of the established velocities across the domain. The
symmetric laminar flow pattern can be observed, which agrees with the physics
of the problem.

[width=1]1d.pdf

Figure 9: Y-velocity comparison in the fully developed region of the outlet of
the T-mixer. The red and green dashed lines represent conventional (FEM) and
Tensor Train Finite Element Method (TetraFEM) solutions, respectively. The
black line refers to the experimental data from Ref. [41].

Figs. 10(a) and 10(b) depict the number of parameters required to keep the
stiffness matrix and mean time required to perform a step of the simulation
with and without tensor-train format, respectively. It can be seen that the
stiffness matrix in TT format requires significantly less memory compared to
the full format. In addition, the runtime comparison demonstrates an
exponential speed-up of TetraFEM in comparison to the conventional realization
when the discretization exceeds $64$ points per axis in a single subdomain.

[width=1]memory.pdf

(a) The number of parameters required to store the stiffness matrix
$\mathbf{S}$ for the T-mixer geometry. A similar plot is observed for the
other FEM matrices. For all discretization values, the TT format allows
efficient compression of the matrices.

[width=1]time.pdf

(b) The mean duration of a single time step of the T-mixer flow simulation. On
dense meshes, our approach (TetraFEM) demonstrates an exponential speed-up
compared to the conventional FEM implementation.

Figure 10:  Comparison between the full and TT versions of the Finite Element
Method calculation. The X-axis corresponds to the discretization factor
(number of grid points along one axis in a single subdomain). (a) depicts the
number of parameters required to store the stiffness matrix, (b) compares real
times required to perform a single time-step of Algorithm 1. When the number
of discretization points exceeds $64$, the memory of our machine becomes
insufficient to store all the matrices, so the dashed line is obtained by
extrapolation.

We also analyzed the accuracy of the TetraFEM solution depending on the TT-
rounding tolerance $\varepsilon$ (Fig. 12). The relative error of TetraFEM was
calculated with respect to the conventional FEM solution in terms of the
Frobenius norm between the X-component of the velocities. The results mean
that based on the requirements of the method, one can manipulate the TT-
rounding error to get more accurate or faster solutions without altering the
mesh or other parameters of the method.

We performed all the computations on two Intel Xeon 2.2 GHz CPUs and 12 Gb of
RAM.

[width=1]rank.pdf

Figure 11: Dynamics of the mean TT-rank of the pressure vector during
simulation of the T-mixer. As the physical time increases, the flow behavior
becomes more complex, which affects the increase in TT-rank.

[width=1]error.pdf

Figure 12: The dependence of the relative error of our implementation
(TetraFEM) compared to the conventional one (FEM) for various Rounding
tolerances $\varepsilon$. Here, the mesh discretization is equal to $32$.

##  VI Extension to a Quantum Computer

This section outlines how utilizing a quantum computer can provide further
benefits to this method. As we see, we can realize the whole numerical scheme
as operations with Tensor Trains. However, if we observe the dynamics of TT-
ranks during the iterations, we will see that they increase over time (Fig.
11). This happens because the flow becomes more complex over time – there are
more correlations in the solution vector during the iterations. Thus, if we
want to solve problems with more complex flows or geometric shapes, or solve
the problem more precisely, then we will inevitably face an increase in TT-
ranks and, accordingly, the complexity of the entire algorithm.

A quantum computer can help us resolve this problem. Indeed, a quantum
computer can do the following: multiply an initial state-vector, $X_{0}$, by
matrices:

$X=V_{n}V_{n-1}...V_{2}V_{1}X_{0},$

where $V_{i}$ are unitary matrices. Our assumption is that if we have a
numerical scheme written as operations with Tensor Trains (as in our case),
then we will be able to efficiently create a $V_{i}$ matrix for each iteration
step. By efficiency, we mean that a quantum circuit, which corresponds to each
$V_{i}$ decomposed into single and two-qubit gates, is shallow and does not
require full connectivity between the qubits. For example, any Tensor Train
can be efficiently encoded into a quantum circuit [18, 42, 43, 46]. We
conjecture that other operations can also be extended, for example,
multiplication by a TT-matrix seems to be the most native operation we are
currently working on. The most important thing is that if we manage to
implement all operations with tensor networks on a quantum computer, then we
will be able to solve problems with any TT-ranks since the complexity of
quantum computer algorithms will not depend on them at all.

##  VII CONCLUSION

To summarize, we propose to use the Tensor Train Finite Element Method for the
solution of the Navier-Stokes equation in complex geometries. For this, we
divide a complex domain into curvilinear quadrilaterals, transform them into
quads, and generate local stiffness, mass and first-derivative matrices. Then,
we assemble the global FEM matrices taking into account the connection of
subdomains on the borders. We represent these matrices and initial vectors in
the Quantized Tensor Train format in order to implement an explicit iterative
scheme for the numerical solution of Navier-Stokes in the TT format. This
allows us to exponentially reduce the memory consumption and gain an
exponential speed up in comparison to the conventional Finite Element Method.
In addition, the results show that the equations are solved correctly since
the solutions coincide with the experimental data. Moreover, we can control
the solution accuracy and lower it, if desired, in return for a higher speed.

We plan to develop an extended method for curvilinear domains with analytic
boundaries in 3D domains, wrapping it in a multi-purpose Finite Element
toolbox. It will also possess a more accurate connection of sub-domains and
faster reassembling of the convection term matrix $-(\vec{v}\cdot\nabla)$ on
each iteration step. This efficient and flexible toolbox will allow us to
solve PDEs in more complex geometries much more efficiently, which can be
applied to tasks in various fields requiring huge discretization that modern
computers cannot handle.

In addition, we discuss how to extend this approach to a quantum computer. To
do this, one needs to develop the methods that implement all operations with
Tensor Trains on a quantum computer. This will allow us to tackle more complex
flows and achieve better accuracy since the complexity will no longer depend
on the ranks of the solution.

## References

  * [1] Huiuk Yi, Minsik Kim, Dongkil Lee, and Jongmyung Park.  Applications of computational fluid dynamics for mine ventilation in mineral development.  Energies, 15(22), 2022.
  * [2] Alessandro Candreva, Giuseppe De Nisco1, Maurizio Lodi Rizzini, Fabrizio D’Ascenzo, Gaetano Maria De Ferrari, Diego Gallo, Umberto Morbiducci, and Claudio Chiastra.  Current and future applications of computational fluid dynamics in coronary artery disease.  RCM, 23(11):377–null, 2022.
  * [3] Bin Han, Xiaoliang Zhu, Bao-Wen Yang, Shenghui Liu, and Aiguo Liu.  Verification and validation of CFD and its application in PWR fuel assembly.  Progress in Nuclear Energy, 154:104485, 2022.
  * [4] Sergey V Ershkov, Evgeniy Yu Prosviryakov, Natalya V Burmasheva, and Victor Christianto.  Towards understanding the algorithms for solving the Navier–Stokes equations.  Fluid Dynamics Research, 53(4):044501, 2021.
  * [5] Terry L Holst.  Numerical solution of the Navier-Stokes equations about three-dimensional configurations-a survey.  In NASA Conference for Supercomputing in Aerospace, pages 281–298, 1987.
  * [6] Román Orús.  Tensor networks for complex quantum systems.  Nature Reviews Physics, 1(9):538–550, 2019.
  * [7] Ivan Oseledets.  Tensor-train decomposition.  SIAM J. Scientific Computing, 33:2295–2317, 01 2011.
  * [8] Sergey Dolgov and Dmitry Savostyanov.  Alternating minimal energy methods for linear systems in higher dimensions.  SIAM Journal on Scientific Computing, 36(5):1–24, September 2014\.
  * [9] Konstantin Sozykin, Andrei Chertkov, Roman Schutski, Anh-Huy Phan, Andrzej Cichocki, and Ivan Oseledets.  TTOpt: A maximum volume quantized tensor train-based optimization and its application to reinforcement learning.  arXiv preprint arXiv:2205.00293, 2022.
  * [10] Dmitry Morozov, Artem Melnikov, Vishal Shete, and Michael Perelshtein.  Protein-protein docking using a tensor train black-box optimization method.  arXiv preprint arXiv:2302.03410, 2023.
  * [11] Alexander Novikov, Dmitrii Podoprikhin, Anton Osokin, and Dmitry P Vetrov.  Tensorizing neural networks.  Advances in neural information processing systems, 28, 2015.
  * [12] Asel Sagingalieva, Andrii Kurkin, Artem Melnikov, Daniil Kuhmistrov, Michael Perelshtein, Alexey Melnikov, Andrea Skolik, and David Von Dollen.  Hyperparameter optimization of hybrid quantum neural networks for car classification.  arXiv preprint arXiv:2205.04878, 2022.
  * [13] A Naumov, Ar Melnikov, V Abronin, F Oxanichenko, K Izmailov, M Pflitsch, A Melnikov, and M Perelshtein.  Tetra-AML: Automatic machine learning via tensor networks.  arXiv preprint arXiv:2303.16214, 2023.
  * [14] Jiri Blazek.  Computational fluid dynamics: principles and applications.  Butterworth-Heinemann, 2015.
  * [15] Sergey V Dolgov, Boris N Khoromskij, and Ivan V Oseledets.  Fast solution of parabolic problems in the tensor train/quantized tensor train format with initial application to the Fokker–Planck equation.  SIAM Journal on Scientific Computing, 34(6):A3016–A3038, 2012.
  * [16] S. Dolgov and J. W. Pearson.  Preconditioners and tensor product solvers for optimal control problems from chemotaxis.  SIAM Journal on Scientific Computing, 41(6):B1228–B1253, 2019.
  * [17] Nikita Gourianov, Michael Lubasch, Sergey Dolgov, Quincy Y van den Berg, Hessam Babaee, Peyman Givi, Martin Kiffner, and Dieter Jaksch.  A quantum-inspired approach to exploit turbulence structures.  Nature Computational Science, 2(1):30–37, 2022.
  * [18] C. Schön, E. Solano, F. Verstraete, J. I. Cirac, and M. M. Wolf.  Sequential generation of entangled multiqubit states.  Phys. Rev. Lett., 95:110503, Sep 2005.
  * [19] Alexandra Bünger, Sergey Dolgov, and Martin Stoll.  A low-rank tensor method for pde-constrained optimization with isogeometric analysis.  SIAM Journal on Scientific Computing, 42(1):A140–A161, 2020.
  * [20] Larisa Markeeva, I Tsybulin, and I Oseledets.  QTT-isogeometric solver in two dimensions.  Journal of Computational Physics, 424:109835, 2021.
  * [21] Frank Gaitan.  Finding flows of a Navier–Stokes fluid through quantum computing.  npj Quantum Information, 6(1):1–6, 2020.
  * [22] Shengbin Wang, Zhimin Wang, Wendong Li, Lixin Fan, Guolong Cui, Zhiqiang Wei, and Yongjian Gu.  A quantum Poisson solver implementable on NISQ devices.  arXiv preprint arXiv:2005.00256, 2020.
  * [23] Andrew M Childs, Jin-Peng Liu, and Aaron Ostrander.  High-precision quantum algorithms for partial differential equations.  Quantum, 5:574, 2021.
  * [24] Aram W Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Physical review letters, 103(15):150502, 2009.
  * [25] MR Perelshtein, AI Pakhomchik, AA Melnikov, AA Novikov, A Glatz, GS Paraoanu, VM Vinokur, and GB Lesovik.  Solving large-scale linear systems of equations by a quantum hybrid algorithm.  Annalen der Physik, 534(7):2200082, 2022.
  * [26] Andrew M Childs, Robin Kothari, and Rolando D Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing, 46(6):1920–1950, 2017.
  * [27] Anna M Krol, Aritra Sarkar, Imran Ashraf, Zaid Al-Ars, and Koen Bertels.  Efficient decomposition of unitary matrices in quantum circuit compilers.  Applied Sciences, 12(2):759, 2022.
  * [28] Leigh Lapworth.  A hybrid quantum-classical CFD methodology with benchmark HHL solutions, 2022.
  * [29] Michael Perelshtein, Asel Sagingalieva, Karan Pinto, Vishal Shete, Alexey Pakhomchik, Artem Melnikov, Florian Neukart, Georg Gesek, Alexey Melnikov, and Valerii Vinokur.  Practical application-specific advantage through hybrid quantum computing.  2022\.
  * [30] Ivan Oseledets and Eugene Tyrtyshnikov.  TT-cross approximation for multidimensional arrays.  Linear Algebra and its Applications, 432(1):70–88, 2010.
  * [31] B. Khoromskij.  $O(d\log N)$-quantics approximation of $N-d$ tensors in high-dimensional numerical modeling.  Constructive Approximation - CONSTR APPROX, 34, 01 2009.
  * [32] Sergey V. Dolgov and Dmitry V. Savostyanov.  Alternating minimal energy methods for linear systems in higher dimensions.  SIAM Journal on Scientific Computing, 36(5):A2248–A2271, jan 2014\.
  * [33] Steven R. White.  Density matrix formulation for quantum renormalization groups.  Phys. Rev. Lett., 69:2863–2866, Nov 1992.
  * [34] Yousef Saad.  Iterative Methods for Sparse Linear Systems.  Other Titles in Applied Mathematics. SIAM, second edition, 2003.
  * [35] Ivan Oseledets.  Constructive representation of functions in low-rank tensor formats.  Constructive Approximation, 37, 09 2010.
  * [36] L. Markeeva, I. Tsybulin, and I. Oseledets.  QTT-isogeometric solver in two dimensions.  2018\.
  * [37] William J. Gordon and Linda C. Thiel.  Transfinite mappings and their application to grid generation.  Applied Mathematics and Computation, 10-11:171–233, 1982.
  * [38] Howard Elman, David Silvester, and Andrew Wathen.  Finite elements and fast iterative solvers: With applications in incompressible fluid dynamics.  01 2006.
  * [39] Chorin A. J.  The numerical solution of the Navier-Stokes equations for an incompressible fluid.  Bulletin of the American Mathematical Society, 73(6):928–931, 1967\.
  * [40] Ivan Oseledets.  ttpy v.1.2.1, 2022.
  * [41] Pamela Vocale, Giacomo Puccetti, Beatrice Pulvirenti, and Gian Luca Morini.  Experimental and numerical analysis of single phase flow in a micro T-junction.  4th Micro and Nano Flows Conference (MNF2014), 2014.
  * [42] Shi-Ju Ran.  Encoding of matrix product states into quantum circuits of one- and two-qubit gates.  Phys. Rev. A, 101:032310, Mar 2020.
  * [43] Artem A. Melnikov, Alena A. Termanova, Sergey V. Dolgov, Florian Neukart, and Michael R. Perelshtein.  QPrep: Library for quantum state preparation using tensor networks.
  * [44] Michael Lubasch, Pierre Moinier, Dieter Jaksch.  Multigrid renormalization.  Journal of Computational Physics, 372, 587-602, 1 Nov 2018
  * [45] Michael Lubasch, Jaewoo Joo, Pierre Moinier, Martin Kiffner, Dieter Jaksch.  Variational quantum algorithms for nonlinear problems.  Phys. Rev. A, 101:10301, 6 Jan 2020
  * [46] Matan Ben Dov, David Shnaiderov, Adi Makmal, Emanuele G. Dalla Torre.  Approximate encoding of quantum states using shallow circuits.  arXiv preprint arXiv:2207.00028, 2022.
  * [47] Gourianov N.  Exploiting the structure of turbulence with tensor networks.  University of Oxford, 2022.