Tensor network reduced order models for wall-bounded flows

We introduce a widely applicable tensor network-based framework for developing reduced order models describing wall-bounded fluid flows. As a paradigmatic example, we consider the incompressible Navier-Stokes equations and the lid-driven cavity in two spatial dimensions. We benchmark our solution against published reference data for low Reynolds numbers and find excellent agreement. In addition, we investigate the short-time dynamics of the flow at high Reynolds numbers for the lid driven and doubly-driven cavities. We represent the velocity components by matrix product states and find that the bond dimension grows logarithmically with simulation time. The tensor network algorithm requires at most a few percent of the number of variables parameterizing the solution obtained by direct numerical simulation, and approximately improves the runtime by an order of magnitude compared to direct numerical simulation on similar hardware. Our approach is readily transferable to other flows, and paves the way towards quantum computational fluid dynamics in complex geometries.


I. INTRODUCTION
Direct numerical simulation (DNS) of the Navier-Stokes equations at large Reynolds numbers would be a highly desirable capability for science and engineering applications.However, it remains an elusive goal due to the extremely large numerical complexity associated with the multiscale nature of turbulence [1,2].The state-ofthe-art method for mitigating this issue in computational fluid dynamics (CFD) is turbulence modelling [3], which continues to be under constant development for improving its accuracy.
A conceptionally different approach to reducing the numerical complexity of DNS is through structure-resolving methodologies [4,5].These methods aim to establish a reduced order model (ROM) of the full system by exploiting correlated structures in the solution.However, identifying suitable modes for building ROMs is difficult and therefore under active investigation [5][6][7][8].
Recently, quantum-inspired tensor network methods have been introduced as a novel paradigm for modelling turbulent flows for diagnostic and predictive purposes [9].The tensor network algorithm in [9] for solving the incompressible Navier-Stokes equation (INSE) approximates the velocity components in matrix product state (MPS) format [10].In the examples studied in [9], the number of variables parameterizing the solution (NVPS) in MPS representation is reduced by over an order of magnitude compared to DNS.The MPS algorithm thus realizes a ROM for the investigated flows.However, the efficient compression reported in [9] only resulted in a computational speedup in a one-dimensional system, but not in two or three spatial dimensions.Furthermore, the examples in [9] are restricted to homogeneous flows with periodic boundary conditions.
Here we show that ROMs based on tensor networks can be extended to wall-bounded flows.We illustrate our approach using the lid-driven cavity in two spatial dimensions, which is a very well-studied problem [11] with tabulated reference solutions [12].We solve the INSE in the streamfunction-vorticity formulation and find that our MPS algorithm reproduces the data in [12] for stationary states at low Reynolds numbers.
As an application of our approach, we assume that the fluid is initially at rest and investigate the short-time dynamics at high Reynolds numbers.We represent the velocity components by MPSs with bond dimension χ and investigate how χ depends on time and grid size.We find that χ grows logarithmically in time and reduces the NVPS compared to direct numerical simulation by about 97%.
As an extension towards more complex flows, we also investigate the doubly-driven cavity where both the top and bottom lids move and find the same qualitative behavior.We compare the runtimes of the MPS and DNS algorithms on similar hardware and at different Reynolds numbers.We find that the MPS algorithm can give rise to significant runtime improvements compared to DNS, peaking at a seventeen-fold speedup in case of the liddriven cavity.
The MPS algorithm in [9] advances the solution to the INSE by solving an optimization problem.More specifically, the continuity equation is combined with the momentum equations via the penalty method [13], and the updated velocity components are obtained by minimizing a single cost function.On the contrary, the MPS algorithm in this work is constructed by emulating the DNS algorithm step-by-step.We achieve this by decomposing the DNS algorithm into four elementary operations (multiplication, addition, matrix-vector operations and solving linear systems of equations) that can be realized in MPS format.It follows that our approach is directly transferable to a broad class of other CFD methodologies and flow geometries.
An important feature of quantum-inspired tensor network algorithms is that they can be ported to a quantum computer [9,14].This transfer can be achieved with quantum circuits of known depth [15] and will provide at least a quadratic speedup over the scaling of the classical tensor network algorithm with the bond dimension [9].
Improved speedups may be achieved by problem-specific quantum circuits [15][16][17] that perform exponentially better than the MPS encoding of flow fields.Our work thus represents a first step towards efficient quantum algorithms for solving CFD problems with boundary conditions.
This paper is organized as follows.The model for the lid-driven cavity in the streamfunction-vorticity formulation is presented in Sec.II.We give a detailed description of the model and the spatial discretization because this forms the foundation for constructing the MPS algorithm.We outline the encoding of flow fields in MPS format and describe how the DNS algorithm can be transformed into MPS format.All technical details are summarized in Appendices.The results are shown in Sec.III and begin with a validation of our tensor network algorithm against previous work.We then consider the short-time dynamics following the quench by the moving lid and analyze the bond dimension as a function of time and grid size.A summary and discussion of our results is provided in Sec.IV.

II. MODEL
The setup for the lid-driven cavity in two spatial dimensions is shown in Fig. 1(a).We consider a square box with edge length L, and the upper lid moves with velocity u 0 in x-direction.The x component (y component) of the fluid is denoted by u (v).At t = 0, the fluid is at rest, u = v = 0. We consider a viscous fluid with kinematic viscosity ν and seek solutions to the incompressible Navier-Stokes equations in the streamfunction-vorticity approach [13], The streamfunction ψ and the velocity components u and v are connected via where is the vorticity.Throughout this work we scale time in units of t 0 = L/u 0 , length in terms of L and velocities by u 0 .Solutions to Eq. ( 1) are then characterized by the Reynolds number We discretize the interior of the cavity (excluding boundaries) by a uniform grid with K grid points in each spatial dimension.The computational domain thus comprises K 2 equally spaced points r k with grid spacing h = L/(K + 1) . (5) Each grid point vector r k is uniquely described by a tuple of integers, where k α ∈ {0, . . ., K − 1} is the index of the grid point in the direction êα with α ∈ {x, y}.The one-to-one correspondence in Eq. ( 6) allows us to label discrete function values on the grid by F (r k ) ≡ F k x ,k y .We denote ghost points on the left (bottom) boundary by k x = −1 (k y = −1), and those on the right (top) boundary by The streamfunction ψ must vanish everywhere on the boundary, and all velocity components are zero except for u = u 0 on boundary C t [see Fig. 1(a)].The boundary conditions for ψ, u and v are summarized in Tab.I. We obtain the boundary values for the vorticity w in the standard approach [12] and find (p, q ∈ {0, . . ., K − 1}) The DNS algorithm for solving Eq. ( 1) with the boundary conditions in Tab.I is outlined in Appendix A. For the time integration of Eq. (1a), we use a second-order MacCormack algorithm [13,18,19].Finite-difference operations are realized by sparse matrix-vector multiplications, and we use a preconditioned conjugate gradient algorithm for solving the Poisson equation (1b).The self-consistent solution to the set of Eq. ( 1) is found by iteratively solving Eq. (1b) and Eq.(1a) until convergence is achieved.We begin the description of our MPS algorithm with a discussion of the encoding of discrete functions in MPS format.For this we assume that the number of grid points in each spatial dimension is K = 2 N for an integer N .The binary representation (. ..) 2 of a grid point index k α requires N bits, Ct Cr C b C l u u0 0 0 0 v 0 0 0 0 ψ 0 0 0 0 TABLE I. Dirichlet boundary conditions for velocity fields u, v and the streamfunction ψ on boundaries Cα as indicated in Fig. 1(a).where σ α i ∈ {0, 1}, α ∈ {x, y}, i = 1, . . ., N, and σ α 1 and σ α N are the most and least significant bits, respectively.We approximate a discrete function F by an MPS of bond dimension χ and length 2N , where we introduced The matrices M ωn have dimensions d(n−1)×d(n), where are the internal bonds that are summed over in the product of matrices in Eq. ( 9).These bonds are responsible for describing correlations between different length scales [9,20].The first N matrices in Eq. ( 9) encode the y-components of F , and the remaining N matrices account for the x-components.Note that this encoding employs the scale encoding introduced in [9,20] in each spatial dimension separately.The encoding in Eq. ( 9) thus corresponds to expanding the function f as a sum of product functions, where Y i (X i ) is a function of the y index k y (x index k x ) only.We find that this encoding is more efficient for the cavity geometry than the encoding in [9] where combined scales of all spatial dimensions are considered.Note that the encoding in Eq. ( 9) can be straightforwardly generalized to the case where each spatial dimension is discretized by a different number of grid points.This is of interest for more complex geometries than the square box considered here.
Next we describe how we emulate the DNS algorithm in tensor network format.The DNS algorithm can be broken down into the following elementary operations: (i) addition of flow fields, (ii) multiplication of flow fields, (iii) the algorithm for solving the Poisson equation, and (iv) sparse matrix-vector operations.Sparse matrixvector operations realize finite difference operations on the flow fields, as well as the boundary conditions for the vorticity in Eq. ( 7).
Since all operations (i-iv) can be realized in MPS format, the MPS algorithm for solving Eq. ( 1) can be obtained by replacing each elementary operation in the DNS algorithm by its MPS counterpart.MPSs can be added [10] and multiplied [21], and a Poisson solver in MPS format has been reported in [22].Matrix-vector operations are realized by contracting a matrix product operator (MPO) with an MPS [10], and all MPOs for realizing the finite difference operations and boundary conditions are provided in Appendix C. The numerical complexity of all these operations scales polynomially with the bond dimension χ of the MPS [for details see Appendix B].It follows that the MPS realizations of operations (i-iv) can be numerically more efficient than their standard implementations for sufficiently small χ.
All variables (velocity components u and v, streamfunction ψ and vorticity w) are approximated by an MPS with bond dimension χ.We allow χ to dynamically grow in order to keep the numerical complexity of our algorithm minimal.We achieve this by normalizing the MPSs representing ψ and w to unity, and by inspecting the singular values near the center of these MPSs.We increase the bond dimension if the smallest singular value exceeds a threshold ϵ, which we set to ϵ = 5×10 −8 throughout this work.This choice has been informed by numerical tests, ensuring that all precision targets of the algorithms implementing the elementary operations are met with the smallest possible χ.

III. RESULTS
In a first step we validate the MPS algorithm against the tabulated results for the stationary state of the liddriven cavity in [12].We consider a 2 7 × 2 7 grid (N = 7) and Reynolds number Re=1,000.The contours showing the velocity magnitude according to the MPS algorithm and for t/t 0 = 50 are shown in Fig. 1(b).We compare this to the data in [12] in Figs.1(b) and (c).The velocity component u along the vertical line L v [see Fig. 1(a)] according the the MPS algorithm (black solid lines) agrees very well with the data in [12] (red dots).Similarly, we find that our MPS results for v along the horizontal line L h agree very well with [12] as shown in Fig. 1(d).We note that our DNS algorithm is in excellent agreement with the MPS algorithm and with the reference data in [12].
The MPS algorithm dynamically adapts the bond dimension of the MPS representing the flow fields.Initially the fluid it as rest, u = v = 0.This constant velocity field is an MPS with bond dimension χ = 1.However, we find through numerical experiments that the sudden quench induced by the moving lid requires a starting bond dimension of χ = 26.The subsequent evolution of χ with time for the flow in Fig. 1(b) is shown in Fig. 2(a).We find that χ approximately grows logarithmically with time until t/t 0 ≈ 3, and then it stays constant at χ = 38.At t/t 0 ≈ 3, the vortex created by the moving lid has expanded from the top right corner to the whole size of the cavity.While the vortex changes shape until the steady state is reached, the bond dimension stays constant in this regime.
The bond dimension χ is directly related to the NVPS, which is shown in Fig. 2(b) in relation to the total number of grid points K 2 .Initially the NVPS are about 47% of K 2 .For larger times, the NVPS slowly increases to 82% of K 2 .It follows that the MPS format does not result in an efficient compression of the stationary state for Re=1,000.
The situation is completely different in a transient regime at high Reynolds numbers.For this we consider a flow with Re = 24, 000, and Fig. 3 grid at t/t 0 = 3.A magnified view of the vortex forming in the top right corner is shown in Fig. 3(b).It is well known that the lid-driven cavity only exhibits a truly stationary state for Re ≤ 10, 000 [23].For larger Reynolds numbers, the system becomes chaotic and develops random fluctuations that persist for all times.However, we find that at the short times considered here where turbulence has not formed yet, the system is still deterministic.All runs with the same initial conditions give the same result.We expect the onset of turbulence and non-stationary fluctuations at much later times when the vortex has spread to the whole cavity.
Next we investigate the required grid size to correctly represent this transient flow.For this we run the calculation for different grid sizes 2 N × 2 N with N = 8, 9, 10, 11 and 12.The results for N = 9 and N = 10 are shown in Figs.3(c) and (d), respectively.By comparing it to the solution for N = 11 in Fig. 3(b), we find that the flow fields are underresolved on the N = 9, 10 grids.For N = 9 [see Fig. 3(c)], the vortex in the upper right corner is strongly deformed.The amount of deformation is much smaller but still visible for the N = 10 grid [see Fig. 3(d)].On the other hand, increasing the size to N = 12 (not shown) does not result in any significant changes compared with the results for N = 11.We thus conclude that the 2 11 × 2 11 grid is sufficiently large for representing this flow.
The smallest length scale in a fully developed turbulent flow is the Kolmogorov microscale η/L ≈ Re −3/4 [1,2].Although the flow investigated in Fig. 3 is not in the turbulent regime yet, the value of η/L ≈ 5.19 × 10 −4 for Re=24,000 is consistent with the grid point spacing 2 −11 ≈ 4.88 × 10 −4 for the 2 11 × 2 11 grid that resolves this flow.We conclude that the smallest scale according to Kolmogorov theory is excited even in the investigated regime where the flow is still laminar.It follows that η/L gives a reasonable estimate for the required grid size.
The variation of χ with time and for the flow in Fig. 3(a) is shown by the black crosses in Fig. 4(a).At t = 0 we set χ = 40, and after a short initial phase [not shown in Fig. 4(a)] we find that χ grows logarithmically in time.The corresponding NVPS in relation to the total number of grid points is shown by the black crosses in Fig. 4(b).At very short times, the NVPS are only about 1% of K 2 , and thus the MPS format achieves a compression of 99%.For larger times, the NVPS slowly increases to 3.4% of K 2 , corresponding to a compression of 96.6%.
Next we investigate the dependence of the bond dimension on the grid size.We find that for all studied grids (N = 8, 9, 10, 11, 12), χ vs. time has the same qualitative behavior as shown in Fig. 4(a) for N = 11.For each of these curves, we calculate the maximal value χ max at t/t 0 = 3 and the temporally averaged bond dimension χ.The results for χ max and χ are shown by black crosses in Figs.4(a) and (b), respectively.We find that χ max and χ vary with 2N until the grid is fine enough to represent the flow.While χ increases steadily with 2N , χ max first increases then decreases with 2N .We now investigate how the results for the bond dimension in the lid-driven cavity geometry change if we consider a doubly-driven cavity instead, see Fig. 5.The upper lid continues to move at constant velocity u 0 in xdirection.In addition, the bottom lid moves at constant velocity −u 0 in x-direction.The corresponding contours of the velocity magnitude on a 2 11 × 2 11 grid and with Re=24,000 are shown in Fig. 5(b).We find that a second vortex forms in the bottom left corner of the cavity.The corresponding results for the bond dimension as a function of time and grid size are shown by the red circles in Fig. 4. The qualitative behavior of all curves is similar to the lid-driven cavity, but the bond dimension for the doubly-driven cavity is larger than for the lid-driven cavity at each point in time.The NVPS in relation to the total number of grid points grows to about 9% for the doubly-driven cavity, and hence the MPS format still achieves a compression of more than 90%.
The results in Fig. 4 show that the bond dimension only grows logarithmically with simulation time, and that the MPS format achieves an efficient compression of the flow fields.The numerical complexity of the MPS algorithm depends on the bond dimension χ as detailed in Appendix B. While the most costly operation is the multiplication of two MPSs, the algorithm spends the most time on solving the Poisson equation which scales as 2N χ 3 [22].On the other hand, the DNS algorithm can be broken down into sparse matrix-vector multiplications scaling with the total number of grid points K 2 = 2 2N .The exponentially worse scaling of the DNS algorithm with respect to the number of grid points K 2 suggests that the MPS algorithm can give rise to a computational advantage for sufficiently small values of χ.
To address this question we compare the runtimes of the MPS and DNS algorithms.In order to achieve a fair comparison, we implemented the DNS and MPS algorithms in the same programming language (i.e., Matlab [24]), and evaluated all runs on a single node of the ARC facility (Intel Xeon Platinum 8268 CPU @ 2.90GHz) [25].Furthermore, we ensure that the DNS and MPS algorithms solve Eq. ( 1) with the same accuracy (see Appendix B).We find that the MPS algorithm is 5.8 times faster than the DNS algorithm in the case of the lid-driven cavity.The speedup reduces to 3.3 for the doubly-driven cavity since the bond dimension is larger than for the lid-driven cavity at each time step, see Fig. 4(a).
A more comprehensive runtime comparison of the MPS and DNS algorithms at different Reynolds numbers is presented in Fig. 5(a).The grid spacing for each Re is chosen such that it matches the corresponding microscale η/L ≈ Re −3/4 .We show the ratio of the average times T DNS for completing a single iteration of the DNS algorithm and T MPS for completing a single iteration of the MPS algorithm.Since the MPS and DNS algorithms approximately require the same number of iterations, this ratio is also representative of the overall runtime ratio.The MPS algorithm for the lid-driven and doubly-driven cavities runs faster than the DNS algorithm for Re ≥ 9.5 × 10 3 .For a given Reynolds number, the MPS algorithm for the doubly-driven cavity takes more time than in the case of the lid-driven cavity because the former requires larger bond dimensions, see Fig. 5(b).For the largest Reynolds number, the MPS algorithm approximately achieves a seventeen-fold [tenfold] speedup compared with the DNS algorithm for the lid-driven [doubly-driven] cavity.
The speedups shown in Fig. 6 can be qualitatively explained by noting that the DNS algorithm scales like Re 6/4 , whereas the MPS algorithm scales as log Re for fixed bond dimension.However, the bond dimension grows with time and with Reynolds number, and therefore a general scaling of the runtime ratio with Reynolds number is difficult to obtain.At larger simulation times, the runtime advantage of the MPS algorithm may decrease or vanish if the required bond dimension becomes too large.The results in Fig. 6 nevertheless illustrate the tremendous potential of MPS for simulating transient flows.

IV. SUMMARY AND DISCUSSION
We have shown that dynamical solutions to the incompressible Navier-Stokes equations for the lid-driven and doubly-driven cavities can be obtained via a tensor network algorithm.Our work extends the results in [9] by showing that the tensor network approach is not restricted to periodic boundary conditions but works equally well for problems with fixed boundary conditions.We achieve this by decomposing a DNS algorithm based on MacCormack's method [13,18,19] into four elementary operations of addition, multiplication, matrix-vector multiplication and solving the Poisson equation.These four operations can be implemented in MPS format and the resulting MPS algorithm automatically builds a ROM characterized by a bond dimension χ.Note that this ROM becomes exact with sufficiently large bond dimen- ) is only taken for t/t0 ≤ 1 (t/t0 ≤ 0.1) due to the large runtimes, and T MPS for the doubly-driven cavity and Re=60.5k is evaluated for t/t0 = 2.1.The grid spacing for each Re is chosen such that it matches the corresponding microscale η/L ≈ Re −3/4 .For data points above (below) the horizontal blue dashed line, the MPS (DNS) algorithm runs faster than its DNS (MPS) counterpart.Solid lines are a guide to the eye.(b) Time-averaged bond dimensions χ(LD) [ χ(DD)] corresponding to the lid-driven [doubly-driven] cavity for different Reynolds numbers.
sion, which distinguishes it from data-driven ROMs [6][7][8] for CFD which lack this guarantee of success.
It is important to note that our approach also applies to other CFD methodologies and flow geometries.For example, the streamfunction-vorticity formulation chosen in this work can be replaced with continuity and momentum equations expressed in terms of velocities and pressure [26].Re-writing this algorithm in terms of tensor network operations follows the same route as presented here.
We run the MPS algorithm on a uniform grid and find that it automatically allocates resources only to those regions in space where they are needed.No a priory knowledge of the flow is required.For example, the NVPS required by MPS to describe the transient regime at large Reynolds number is at most 3% of the total number of gridpoints.This very efficient MPS representation of the flow occurs because the vortex only occupies a small region in space.Very little resources are needed to represent the flow in the large area where the fluid is nearly at rest, see Fig. 3(a).Adding the second vortex in the case of the doubly-driven cavity increases the NVPS to 9%.
A related finding is that the bond dimension of the MPSs representing the flow fields is approximately constant if the grid is fine enough to represent the flow.This feature is related to the known fact that polynomials and Fourier series have efficient MPS representations where the bond dimension is independent of the grid size [27,28].This behavior is also akin to onedimensional quantum systems obeying an area law [29].
We find that the bond dimension of the MPSs repre-senting the transient flows investigated in this work grows logarithmically in time.This slow increase can translate into a runtime advantage of the MPS vs. DNS algorithms if the bond dimension of the initial flow fields is sufficiently small.We find that the MPS algorithm can be significantly faster than the DNS algorithm for simulation times of several units of t 0 = L/u 0 , i.e., the time it takes the lid to traverse the length L of the cavity.In general, our analysis shows that the MPS algorithm will outperform the DNS algorithm at large Reynolds numbers, provided that the required bond dimension is sufficiently small.We anticipate that the maximal bond dimension allowing for a speedup depends on the used hardware and software implementation of the algorithm, which is subject to further study.
Several avenues for further research emerge from here.First, the transient flow example studied in this work may also be efficiently described with adaptive mesh refinement [30,31].In this approach, the grid spacing is dynamically varied in space at the cost of detecting the areas requiring high-resolution grids.It would be interesting to directly compare the performance of these two methods for different flow types, and to establish the differences and similarities between them.
Second, MPS algorithms for CFD may benefit from modern hardware architectures optimized for tensor operations [32].This opens up the exciting prospect of developing tensor network algorithms for technical flows that outperform state-of-the-art CFD algorithms.
Finally, CFD algorithms in tensor network format represent a first step towards solving the Navier-Stokes equations on a quantum computer [15,33,34].Quantum CFD [9,14] promises to enable DNS for analyzing and optimizing technical flows, which would represent a revolutionary improvement of the state-of-the-art [33].Creating and benchmarking quantum CFD algorithms for wall-bounded flows by porting tensor network algorithms to quantum hardware is thus an exciting prospect for future research.and δ fwd x , respectively, [δ fwd y w] p,q = w p,q+1 − w p,q h .(A6b) We update the functions F and G with the predicted solution wt+∆t and obtain F t+∆t and Ḡt+∆t .The solution for the vorticity w t+∆t at t + ∆t is then obtained by approximating the spatial derivatives in Eq. (1a) by first-order accurate backward differences, With the help of the boundary values for w in Eq. ( 7), Eq. (A7) can be evaluated on every point of the inner grid with p, q ∈ {0, . . ., K − 1}.Although the forwardand backward differences in Eqs.(A5)-(A7) are only firstorder accurate in h, the resulting expression for w t+∆t p,q in Eq. (A7) is second-order accurate [13,18,19].
(A8) (iv) The set of equations ( 1) are coupled because the updated streamfunction ψ t+∆t gives rise to new velocity components u t+∆t and v t+∆t via Eq.( 2).We repeat steps (i)-(iii) until a self-consistent solution to Eq. (1) has been found.This results in updated functions ψ t+∆t , w t+∆t , u t+∆t and u t+∆t and completes the time step from t to t + ∆t.We repeat steps (i)-(iv) until the final time is reached.

Appendix B: MPS Algorithms
Table II outlines the MPS algorithms for realising the required elementary operations as well as their scaling with the bond dimension χ.All these algorithms have in common that they are variational in nature.The desired MPS for representing the target, i.e., the sum or product of MPSs or the solution to the Poisson equation, is found by minimising a cost function.These cost functions are quadratic in the variables and hence efficient and reliable methods for finding optimal solutions exist.We employ single-site DMRG-like [10] sweeps where each tensor in the MPS is sequentially optimised until overall convergence has been achieved.
In order to make the results of the MPS algorithm comparable to the DNS results, we impose the same accuracy goal for solving the Poisson equation and the same convergence criterion for solving Eq. ( 1) in both algorithms.Here we show how the required finite difference operations can be created in the MPO-MPS formalism.We denote an MPO by Q and its contraction with an MPS f as Qf .A generic MPO with bond dimension D can be written as [35] where A is a 1×D row vector, C is a D ×1 column vector, and B[k] with k ∈ {1, . . ., 2N } are D × D matrices whose matrix elements are 2 × 2 matrices.Any 2 × 2 matrix can be expanded in terms of the following four operators, For convenience, we also introduce the identity matrix When multiplying the matrices B[k] in Eq. (C1), we take the outer product of the matrix-valued matrix elements.
In order to illustrate this notation, we consider the following example for N = 1, A = (1, 0) , (C4) where Q fwd y is given in Eq. (C11) and the MPS f u0 is defined in Eq. (C23).The MPO Q Ct is defined in Eq. (C22).In Eq. (C28), f fwd wy is the MPS representing δ fwd y w. • Backward-forward-difference of w in y-direction: where Q C b is defined in Eq. (C20) and Q b extracts the values of a function on the line k y = 0, with (C30e)

Central differences
Here we provide the MPO representations for the central differences in Eq. (A1).
• Second-order accurate approximation of the first derivative in x-direction: with (C32e) • Second-order accurate approximation of the first derivative in y-direction: with A = (1/2, 0, 0) /h , (C33b) C = (0, 1, −1) t . (C33e) FIG. 1. (Color online) (a) Setup of the square lid-driven cavity with edge length L. The upper lid moves at constant velocity u0 in x-direction.Ct, Cr, C b , C l are the top, right, bottom and left boundaries, respectively.Lv (L h ) denotes a vertical (horizontal) line through the center of the cavity.(b) Contour plot of the velocity magnitude s = √ u 2 + v 2 at t = 50 for Re=1000 and evaluated with the tensor network algorithm.(c) Comparison of the tensor network solution for the x-component u of the velocity along Lv (black solid line) with the reference values in [12] (red dots).(d) Comparison of the tensor network solution for the y-component v of the velocity along L h (black solid line) with the reference values in [12] (red dots).

FIG. 2 .
FIG. 2. (a) Bond dimension χ versus time on a logarithmic scale and for the flow in Fig. 1(b).(b) The ratio between the NVPS and the total number of grid points K 2 in percent and as a function of time.Solid lines are a guide to the eye.

11 s
FIG. 3. (Color online) (a) Contour plot of the velocity magnitude s = √ u 2 + v 2 at t/t0 = 3 for the flow configuration shown in Fig. 1(a).The grid size is K 2 = 2 11 × 2 11 , Re=24,000 and results are obtained with the MPS algorithm.(b) Same as in (a) but focussing on the region where the initial vortex forms.(c) Same as in (b) but for K 2 = 2 9 × 2 9 .(d) Same as in (b) but for K 2 = 2 10 × 2 10 .

FIG. 4 .
FIG. 4.Analysis of the bond dimension as a function of time and grid size.Black crosses [red circles] correspond to the flow in Fig.3(a) [Fig.5], and solid lines are a guide to the eye.(a) Bond dimension χ versus time on a logarithmic scale.(b) The ratio between the NVPS and the total number of grid points K 2 in percent and as a function of time.(c) Bond dimension χmax at t/t0 = 3 as a function of grid size.(d) Temporally averaged bond dimension χ as a function of grid size.

FIG. 5 .
FIG. 5. (Color online) (a) Boundary conditions corresponding to the doubly-driven cavity where the upper [bottom] lid moves at constant velocity u0 [−u0] in x-direction.(b) Contour plot of the velocity magnitude s = √ u 2 + v 2 at t/t0 = 3 for the doubly-driven cavity on a K 2 = 2 11 × 2 11 grid with Re=24,000 and evaluated with the MPS algorithm.

7 FIG. 6 .
FIG. 6. (Color online) (a) Ratio of the average times T MPS for completing a single iteration of the MPS algorithm and T DNS for completing a single iteration of the DNS algorithm as a function of Reynolds number Re. Black crosses [red circles] correspond to the lid-driven [doubly-driven] cavity.Averages are taken up to t/t0 = 3. T DNS for Re=24k (Re=60.5k) is only taken for t/t0 ≤ 1 (t/t0 ≤ 0.1) due to the large runtimes, and T MPS for the doubly-driven cavity and Re=60.5k is evaluated for t/t0 = 2.1.The grid spacing for each Re is chosen such that it matches the corresponding microscale η/L ≈ Re −3/4 .For data points above (below) the horizontal blue dashed line, the MPS (DNS) algorithm runs faster than its DNS (MPS) counterpart.Solid lines are a guide to the eye.(b) Time-averaged bond dimensions χ(LD) [ χ(DD)] corresponding to the lid-driven [doubly-driven] cavity for different Reynolds numbers.

TABLE II .
Overview of the algorithms for realising the building blocks of the DNS algorithm in MPS format.The last column indicates the scaling of the operation with the bond dimension of the MPSs and MPOs.