A Reduced Order Modeling technique to study bifurcating phenomena: application to the Gross-Pitaevskii equation

TL;DR

This paper introduces a computationally efficient reduced order modeling framework that combines continuation, Newton's method, and hyper-reduction to study bifurcations in nonlinear PDEs, demonstrated on the Gross-Pitaevskii equation.

Contribution

The paper presents a novel ROM-based approach that significantly accelerates bifurcation analysis for nonlinear PDEs, especially for the Gross-Pitaevskii equation.

Findings

60 times faster bifurcation diagram construction compared to full order methods

Effective tracing of steady solution branches with reduced computational cost

Applicable to multi-parameter bifurcation studies

Abstract

We propose a computationally efficient framework to treat nonlinear partial differential equations having bifurcating solutions as one or more physical control parameters are varied. Our focus is on steady bifurcations. Plotting a bifurcation diagram entails computing multiple solutions of a parametrized, nonlinear problem, which can be extremely expensive in terms of computational time. In order to reduce these demanding computational costs, our approach combines a continuation technique and Newton's method with a Reduced Order Modeling (ROM) technique, suitably supplemented with a hyper-reduction method. To demonstrate the effectiveness of our ROM approach, we trace the steady solution branches of a nonlinear Schr\"{o}dinger equation, called Gross-Pitaevskii equation, as one or two physical parameters are varied. In the two parameter study, we show that our approach is 60 times faster in constructing a bifurcation diagram than a standard Full Order Method.

Equations

$$G(X(\mathbf{\mu });\mathbf{\mu })=0,$$

$$⟨G(X_h(\mu );\mu ),Y_h⟩\doteq g(X_h(\mu ),Y_h;\mu )=0,\forall Y_h\in V_h.$$

$$dg[X_h^k(\mu )](\delta X_h,Y_h;\mu )=g(X_h^k(\mu ),Y_h;\mu ),\forall Y_h\in V_h,$$

$$X_h^{k+1}(\mu )=X_h^k(\mu )-\delta X_h$$

$$\mathbb{J}({\vec{X}}_h^k(\mu );\mu )\delta {\vec{X}}_h=G_h({\vec{X}}_h^k(\mu );\mu ),$$

$$\mathbb{J}({\vec{X}}_h^k(\mu );\mu ){)}_{ij}=dg[X_h^k(\mu )](E^j,E^i;\mu ),\text{for all}i,j=1,\dots ,N_h.$$

$$g(X_N(\mu ),Y_N;\mu )=0,\forall Y_N\in V_N,$$

$$dg[X_N^k(\mu )](\delta X_N,Y_N;\mu )=g(X_N^k(\mu ),Y_N;\mu ),\forall Y_N\in V_N.$$

$$X_N^{k+1}(\mu )=X_N^k(\mu )-\delta X_N.$$

$$X_N(\mu )=\sum _{m=1}^N{X}_N^{(m)}(\mu ){\Sigma }^m.$$

$$g\left(\sum _{m=1}^N{X}_N^{(m)}(\mu ){\Sigma }^m,{\Sigma }^n;\mu \right)=0,n=1,\dots ,N.$$

$$(G_N({\vec{X}}_N(\mu );\mu ){)}_n=g\left(\sum _{m=1}^N{X}_N^{(m)}(\mu ){\Sigma }^m,{\Sigma }^n;\mu \right),$$

$$(\mathbb{V}{)}_{jm}={\Sigma }_{(j)}^m$$

$${\mathbb{V}}^T{G}_N(\mathbb{V}{\vec{X}}_N(\mu );\mu )=0.$$

$${\mathbb{J}}_N({\vec{X}}_N^k(\mu );\mu )\delta {\vec{X}}_N=G_N({\vec{X}}_N^k(\mu );\mu ),$$

$${\mathbb{J}}_N({\vec{X}}_N^k(\mu );\mu )={\mathbb{V}}^T\mathbb{J}(\mathbb{V}{\vec{X}}_N^k(\mu );\mu )\mathbb{V}.$$

$${\Phi }_{tot}({\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N)=\prod _{i=1}^N\Phi ({\mathbf{r}}_i),$$

$$N={\int }_D\rho d\mathbf{r},\rho =|\Phi {|}^2,$$

$$i{\partial }_t\Phi =-\frac{1}{2}\Delta \Phi +|\Phi {|}^2\Phi +W(r)\Phi \text{in}D,$$

$$\Phi (\mathbf{r},t)=\varphi (r)\exp (-i\mu t),\varphi (r):\overline{D}\to \mathbb{C},$$

$$G(\varphi ;\mu )\doteq -\frac{1}{2}\Delta \varphi +|\varphi {|}^2\varphi +W(r)\varphi -\mu \varphi =0.$$

$$|m,n⟩:={\varphi }_{m,n}\sim H_m(\sqrt{\Omega }x)H_n(\sqrt{\Omega }y)e^{-\frac{r^2}{2}\Omega },$$

$$a(X,Y)=\frac{1}{2}{\int }_D\nabla X\cdot \nabla Yd\mathbf{x},b(X,Y)=\frac{1}{2}\Omega {\int }_D|\mathbf{r}{|}^{2}X\cdot Yd\mathbf{x},$$

$$d(X,Y;\mu )=\mu {\int }_{D}X\cdot Yd\mathbf{x},n(Y,Z)={\int }_D|Z{|}^{2}Z\cdot Yd\mathbf{x},$$

$$c(X,Y,Z)={\int }_D[2(X\cdot Z)Z+|Z{|}^{2}X]\cdot Yd\mathbf{x},$$

$$a(\delta X_h,Y_h)+b(\delta X_h,Y_h)-d(\delta X_h,Y_h;\mu )+c(\delta X_h,X_h^k,Y_h)=$$

$$a(X_h^k,Y_h)+b(X_h^k,Y_h)-d(X_h^k,Y_h;\mu )+n(Y_h,X_h)\forall Y_h\in V_h.$$

$$({\mathbb{A}}_h{)}_{ij}$$

$$({\mathbb{C}}_h{)}_{ij}$$

$$\mathbb{J}({\vec{X}}_h^k(\mu );\mu )={\mathbb{A}}_h+{\mathbb{B}}_h-{\mathbb{D}}_h+{\mathbb{C}}_h.$$

$$a(\delta X_N,Y_N)+$$

$$a(X_N^k,Y_N)+b(X_N^k,Y_N)-d(X_N^k,Y_N;\mu )+n(Y_N,X_N^k)$$

$${\mathbb{J}}_N({\vec{X}}_N^k(\mu );\mu )={\mathbb{A}}_N+{\mathbb{B}}_N-\mu {\mathbb{D}}_N+{\mathbb{C}}_N$$

Full text

A Reduced Order Modeling technique to study bifurcating phenomena: application to the Gross-Pitaevskii equation

Federico Pichi, Annalisa Quaini, Gianluigi Rozza

Abstract. We propose a computationally efficient framework to treat nonlinear partial differential equations having bifurcating solutions as one or more physical control parameters are varied. Our focus is on steady bifurcations. Plotting a bifurcation diagram entails computing multiple solutions of a parametrized, nonlinear problem, which can be extremely expensive in terms of computational time. In order to reduce these demanding computational costs, our approach combines a continuation technique and Newton’s method with a Reduced Order Modeling (ROM) technique, suitably supplemented with a hyper-reduction method. To demonstrate the effectiveness of our ROM approach, we trace the steady solution branches of a nonlinear Schrödinger equation, called Gross–Pitaevskii equation, as one or two physical parameters are varied. In the two parameter study, we show that our approach is 60 times faster in constructing a bifurcation diagram than a standard Full Order Method.

1 Introduction

We consider the problem of finding a solution $X\in V$ such that

[TABLE]

where $\boldsymbol{\mu}\in\mathcal{D}$ denotes a point in a parameter domain $\mathcal{D}\subset\mathbb{R}^{M}$, $V$ a given functional space, $G(\cdot,\boldsymbol{\mu})$ a given nonlinear functional in $X$. We view problem (1) as the strong form of a nonlinear partial differential equation (PDE) or a system of such equations in which $M$ parameters appear. We are interested in situations in which the solution $X$ of (1) differs in character for parameter vectors $\boldsymbol{\mu}$ in different subregions of the parameter domain $\mathcal{D}$. Such situations occur if $X$ undergoes bifurcations as $\boldsymbol{\mu}$ changes from one subregion to another.

In general, one cannot solve (1) for $X$ so that one instead seeks an approximation of the solution in a $N_{h}$-dimensional subspace $V_{h}\subset V$. Such approximation is usually obtained using a so-called Full Order Method (FOM), like for example the Finite Element Method, which is often expensive, especially if multiple solutions are needed. For this reason, one is interested in finding surrogates methods that are much less costly so that obtaining approximations to the exact solution $X$ of (1) for many choices of the parameter vector $\boldsymbol{\mu}\in\mathcal{D}$ becomes feasible. Such surrogates, called Reduced Order Modeling (ROM) techniques, are constructed using a “few” solutions computed with the FOM.

A practical way to realize a ROM is to organize the computation in two steps:

An offline phase: approximation solutions corresponding to selected representative parameters values/system configurations are computed with a FOM and stored, together with other information about the parametrized problem. This is a computationally expensive step usually performed on high performance computing facilities.

• -

An online phase: the information obtained during the offline phase is used to compute the solution for a newly specified value of the parameters in a short amount of time (ideally in real time), even on a relatively low power device such as a laptop or a smartphone.

These split computational procedures are built in such a way that new parameter dependent quantities are easily and quickly computed online, while representative basis functions for selected parameter values and more demanding quantities are pre-computed offline. Among all the possible ROMs, we choose the Reduced Basis (RB) technique [18, 30, 36].

Recent developments of ROM techniques have focused on the reduction of computational time for a wide range of differential problems [9, 37], while maintaining a prescribed tolerance on error bounds [39, 18, 36]. ROMs in the setting of bifurcating solutions are considered in the early papers [26, 27, 28, 29] for buckling bifurcations in solid mechanics. More recently, in [40] it is shown that a POD approach allows for considerable computational time savings for the analysis of bifurcations in some nonlinear dissipative systems. RB methods have been used to study symmetry breaking bifurcations [15, 35] and Hopf bifurcations [34] for natural convection problems. A RB method for symmetry breaking bifurcations in contraction-expansion channels has been proposed in [33]. In [16], steady bifurcations for both a natural convection problem and contraction-expansion channels are investigated with a localized ROM approach. Yano and Patera [41] introduced a RB method for the stability of flows under perturbations in the forcing term or in the boundary conditions, which is based on a space-time framework that allows for particularly sharp error estimates. Furthermore, in [41] it is shown how a space-time inf-sup constant approaches zero as the computed solutions get close to a bifurcating value. A recent work on ROMs for bifurcating solutions in structural mechanics is [31]. Finally, we would like to mention that machine learning techniques based on sparse optimization have been applied to detect bifurcating branches of solutions in [5, 22] for a two-dimensional laterally heated cavity and Ginzburg-Landau model, respectively.

The methodology we propose to plot a bifurcation diagram makes use of three building blocks: i) a continuation technique to properly follow each solution branch, ii) Newton’s method to deal with the nonlinearity of problem (1), and iii) a RB method to efficiently solve the linearized problem obtained from Newton’s method. The novelty of this approach relies on the combination of well-known and assessed methods to obtain a versatile and global approach to any kind of bifurcation problems modeled by parametric PDEs. We present technical insights to build a ROM approach capable of approximating efficiently and accurately the whole bifurcation diagram with the use of a unique reduced basis. As a concrete setting to illustrate our methodology, we consider the Gross–Pitaevskii equation.

Often referred to as a nonlinear Schrödinger equation, the Gross–Pitaevskii equation models certain classes of Bose-Einstein condensates (BECs), a special state of matter formed by identical bosons at ultra-low temperatures. It is well known that the solutions of the Gross–Pitaevskii equation with a parabolic trap in two dimensions exhibit a rich bifurcating behavior [24, 25, 7], which includes symmetry-breaking bifurcations and vortex-bearing states when a (sufficiently strong) rotational angular momentum term is added [14]. The bifurcating behavior becomes even richer for the two-dimensional coupled Gross-Pitaevskii equations [6]. However, for simplicity we stick to the simple Gross-Pitaevskii equation and present a one parameter study (the chemical potential being the only varying parameter) and a two parameter case (varying chemical potential and the normalized trap strength). For both cases, we show that our approach is able to capture the first six solution branches with high accuracy. We stress the need to supplement the RB method with a hyper-reduction technique [4, 8] to enable significant computational time savings with respect to a standard Full Order Method. In particular, we show that if the RB method at the above point iii) is combined with the Discrete Empirical Interpolation Method [8], our approach is up to 60 times faster in constructing a bifurcation diagram than a FOM, making it an ideal tool to study the complex solution behavior of the Gross–Pitaevskii equation and other nonlinear PDE problems.

The work outline is as follows. In Sec. 2, we present the building blocks of our approach to reduce the computational time requested by the construction of the bifurcation diagrams. In Sec. 3, we apply such approach to the Gross–Pitaevskii equation. Numerical results pertaining to the validation of the FOM, reconstruction of bifurcation diagrams with our ROM approach with and without hyper-reduction technique are reported in Sec. 4. Conclusions are drawn in Sec. 5.

2 Numerical approximation of a problem with bifurcations

The nonlinearity in problem (1) can produce a loss of uniqueness for the solution, with multiple solutions branching from a known solution at a bifurcation point. Our aim is to study numerically the associated bifurcation diagrams with contained computational costs. We restrict our attention to steady bifurcations.

The standard assumption for the map $G$ in eq (1) is the continuous differentiability with respect to $X$ and $\boldsymbol{\mu}$. Let $(\bar{X},\bar{\boldsymbol{\mu}})\in V\times\mathcal{D}$ be the known solution, i.e. $G(\bar{X},\bar{\boldsymbol{\mu}})=0$. Let us denote by $D_{X}G(Z,\boldsymbol{\mu}):V\rightarrow V^{\prime}$ and $D_{\boldsymbol{\mu}}G(Z,\boldsymbol{\mu}):\mathcal{D}\rightarrow V^{\prime}$ the partial derivatives of $G$ on a generic point $(Z,\boldsymbol{\mu})\in V\times\mathcal{D}$. A strong assumption usually found in the literature in order to have a local branch of non-singular solutions is that $D_{X}G(\bar{X},\bar{\boldsymbol{\mu}}):V\rightarrow V^{\prime}$ is bijective. Of course, this is not our case: we deal with bifurcation points, which are singularities for the system. Moreover, we do not require $G$ to be affine in $\boldsymbol{\mu}$, because of the nonlinearity of the problem. As we will see later, this forces us to implement a hyper-reduction technique, such as EIM/DEIM [4, 8], to recover efficiency.

2.1 Proposed approach for the branch reconstruction

For simplicity, we present our approach for a scalar parameter $\mu$ (i.e., $M=1$) in parameter domain $\mathcal{D}=[\mu_{i},\mu_{f}]$, although it can easily be extended to a multi-parameter setting (i.e., $M>1$). Our algorithm is based on three building blocks:

A continuation technique: it permits to reconstruct properly the bifurcation diagram by following each branch and providing a suitable initial guess for the nonlinear iterations (next building block). We use a slight modification of the continuation method in [2], which consists in a for loop in $\bar{\mathcal{D}}=[\mu_{1},\dots,\mu_{K}]\subset\mathcal{D}$, a discrete version of the parameter set with cardinality $K$. At each new cycle, we check if a bifurcation occurs, using a threshold $\epsilon_{BIF}$ that controls the norm of the solution. If a bifurcation does occur, we set that bifurcated solution as the initial guess for the next cycle in order to capture the post-bifurcation behavior. 2. 2.

Newton’s method: each point in the bifurcation diagram requires the solution of nonlinear PDE problem for the corresponding value $\mu$. To deal with the nonlinearity, we use the Newton-Kantorovich method [10]. 3. 3.

A Galerkin Finite Element discretization: at each step of the Newton iteration, we have to approximate the solution of a new linear weak formulation. For the space discretization, we use the Galerkin-Finite Element method as discussed in Sec. 2.2. This choice is made also in view of the numerical extension towards the model order reduction.

Algorithm 1 summarizes the three steps for generic nonlinear parametric problem (1).

Let us clarify lines 2-6 in Algorithm 1. In order to decrease the computational cost also during a generic online phase (i.e, also in the case with no hyper-reduction), we avoid applying the continuation method when the bifurcated solutions do not exist yet. To do this, we control the norm of the solution and we keep providing the pre-bifurcation guess, which is different for each of the branches we want to approximate, until the solver converges to a bifurcated solution. At this point, the continuation is enabled and the computed solution becomes the new guess.

We have applied similar strategies for the numerical study of bifurcations arising in different contexts, ranging from fluid dynamics [34, 33, 16, 17] to structural mechanics [31].

2.2 Galerkin Finite Element method for a (generic) nonlinear problem

In this section, we introduce some standard notion for the discretization of generic problem (1) with the Galerkin Finite Element (FE) method.

Let $V_{h}$ be a family of finite dimensional spaces, such that $V_{h}\subset V$. Let $N_{h}=\text{dim}(V_{h})$. We first cast problem (1) in weak form and then for a given parameter $\mu\in\mathcal{D}$ seek $X_{h}(\mu)\in V_{h}$ that satisfies

[TABLE]

To treat the nonlinearity in $G$, we apply the Newton-Kantorovich method [10, 36], which reads as follows. Choose initial guess $X^{0}_{h}(\mu)\in V_{h}$. Then, for every $k=0,1,\dots$:

Step 1: Seek the variation $\delta X_{h}\in V_{h}$ such that

[TABLE]

• -

Step 2: Update the solution

[TABLE]

Step 1 and 2 are repeated until the $L^{2}$-norm of the residual falls below a prescribed tolerance $\epsilon$.

We denote with $\{E^{j}\}_{j=1}^{N_{h}}$ a basis for $V_{h}$. Newton’s method combined with the Galerkin finite element method and applied to problem (2) reads: find $\delta\vec{X}_{h}\in\mathbb{R}^{N_{h}}$ such that

[TABLE]

where the Jacobian matrix in $\mathbb{R}^{N_{h}\times N_{h}}$ is defined as

[TABLE]

2.3 The Reduced Basis method

Once projected onto a suitable Finite Element space, the parametrized discrete problem derived in the previous section leads to a very large nonlinear system which has to be solved for every parameter $\mu\in\mathcal{D}$. The solution computed with a FE method represents the so-called high fidelity approximation, which is computationally expensive. To reduce the computational cost without compromising the accuracy, we choose to use a ROM technique called Reduced Basis (RB) method [18, 30, 36]. Roughly speaking, this method consists in a projection of the high fidelity problem on a subspace of smaller dimension, constructed with some properly chosen basis functions.

RB methods use the offline-online paradigm introduced in Sec. 1. In the offline phase, we explore the parameter space $\mathcal{D}$ in order to construct a basis for the low dimensional manifold, which efficiently approximates the high fidelity space and where the parametrized solutions lie. This entails solving $N_{train}$ times the Galerkin high-fidelity problem associated to $N_{train}$ values of $\mu$ in $\mathcal{D}$. For the numerical results in Sec. 4, we chose $\{\mu^{n}\}_{n=1}^{N_{train}}$ to be an ordered sampling of the interval $\mathcal{D}$. In the online phase the solution is computed through the projection on the low dimensional manifold in an efficient and reliable way for every $\mu\in\mathcal{D}$ we are interested in. The reduced computational cost comes from avoiding to project on the large Finite Element space. To be precise, we want to construct the reduced problem through the projection on a subspace $V_{N}\subset V_{h}$ spanned by a collection of the snapshots, i.e. solutions of the full order problem for selected values of parameter $\mu$, obtained by, e.g., Proper orthogonal decomposition (POD) or Greedy techniques [18, 30, 36].

Let us provide some more details of the online phase for generic problem (1). For a given $\mu\in\mathcal{D}$, we seek $X_{N}(\mu)\in V_{N}$ that satisfies

[TABLE]

where $g(\cdot,\cdot;\mu)$ is defined in (2). Just like for the FOM in Sec. 2.2, we apply the Newton-Kantorovich method. We choose an initial guess $X^{0}_{N}(\mu)\in V_{N}$ and the for every $k=0,1,\dots$:

Step 1: Find the variation $\delta X_{N}\in V_{N}$ such that

[TABLE]

• -

Step 2: Update the solution

[TABLE]

Step 1 and 2 are repeated until the $L^{2}$-norm of the residual falls below a prescribed tolerance $\epsilon$.

Let $\{\Sigma^{m}\}_{m=1}^{N}$ be an orthonormal basis (with respect to the inner product defined on the space $V_{h}$) for $V_{N}$, obtained through POD sampling and the Gram-Schmidt procedure during the offline phase. We remark that this basis will be optimal in the $\ell^{2}$ sense, minimizing over all possible N-dimensional orthonormal bases $W_{N}$ the errors between the snapshots and their projection through $W_{N}$.

Then, $V_{N}=\text{span}\{\Sigma^{1},\dots,\Sigma^{N}\}$ and we can write every $X_{N}(\mu)\in V_{N}$ as

[TABLE]

We denote with $\vec{X}_{N}(\mu)=\{X_{N}^{(m)}(\mu)\}_{m=1}^{N}\in\mathbb{R}^{N}$ the reduced solution vector.

By plugging (8) into (6) and choosing $Y_{N}=\Sigma^{n}\in V_{N}$, for $1\leq n\leq N$, we obtain the following algebraic system

[TABLE]

Let

[TABLE]

be the residual reduced vector. We denote with $\mathbb{V}$ the $N_{h}\times N$ transformation matrix whose elements

[TABLE]

are the nodal evaluation of the m-th basis function at the j-th node. With this new notation, we can rewrite problem (9) as

[TABLE]

Finally, we combine Newton’s method and the RB technique. At every iteration $k$ of Newton’s method the problem that has to be solved reads as follows: find $\delta\vec{X}_{N}\in\mathbb{R}^{N}$ such that

[TABLE]

where $\mathbb{J}_{N}$ is the $\mathbb{R}^{N\times N}$ reduced Jacobian matrix

[TABLE]

We remark that eq. (11) involves the degrees of freedom of the high fidelity problem. Because of this, the repeated assembly of the Jacobian compromises the efficiency of the reduced order method during the online phase. As we will see later, this issue can be overcome by adopting an affine recovery technique, which allows a consistent speed-up of the method by interpolating the nonlinear part of the variational form.

3 Application to the Gross–Pitaevskii equation

The Gross–Pitaevskii equation models certain classes of Bose-Einstein condensates. A BEC is a special state of matter formed by an unlimited number of bosons that “condense” into the same energy state at low temperatures. A BEC is formed by cooling a gas of extremely low density, about one-hundred-thousandth the density of normal air, to ultra-low temperatures (close to absolute zero).

A quantum system is the environment to be studied in terms of wave-particle duality (i.e., all particles exhibit a wave nature and viceversa) and it involves the wave-function and its constituents, such as the momentum and wavelength. The Gross–Pitaevskii equation describes the ground state of a quantum system of identical bosons using two simplifications: the Hartree–Fock approximation and the pseudopotential interaction model. In the Hartree–Fock approximation, the total wave-function $\Phi_{tot}$ of a system of $N$ bosons is taken as a product of single-particle functions $\Phi$:

[TABLE]

where $\mathbf{r}_{i}$ is the coordinate of the $i$-th boson. If the single-particle wave-function satisfies the Gross–Pitaevskii equation, the total wave-function minimizes the expectation value (i.e., the probabilistic expected value of the result of an experiment) of the pseudopotential model Hamiltonian under normalization condition:

[TABLE]

where $D\subseteq\mathbb{R}^{2}$ is the domain under consideration and $\rho$ is interpreted as the particle density. The Gross–Pitaevskii equation reads: Find the single-particle wave-function $\Phi(\mathbf{r},t):\overline{D\times\mathbb{R}^{+}}\rightarrow\mathbb{C}$ such that

[TABLE]

where $i$ is the imaginary unit, $r=|\mathbf{r}|=\sqrt{x^{2}+y^{2}}$ is the radial coordinate, and $W(r)=\frac{1}{2}\Omega^{2}r^{2}$ is the external potential, with $\Omega$ being the normalized trap strength, i.e. the ratio of trappings along and transverse to the plane. In this paper, we set $\Omega=0.2$ unless specified otherwise. Notice that we consider a single well potential. Eq. (13) is similar in form to the Ginzburg–Landau equation and is sometimes referred to as a nonlinear Schrödinger equation. Obviously, eq. (13) needs to be supplemented with suitable boundary conditions.

The construction of the steady solution is based on the ansatz:

[TABLE]

where $\mu$ is the chemical potential, which has to satisfy $\mu\geq\Omega$. By plugging (14) into (13), we obtain nonlinear problem

[TABLE]

It is well known that the solutions of the one-dimensional version of problem (15) exhibit a bifurcating behavior [21, 20, 1, 11], which is not particularly rich though. The bifurcations occurring in the two-dimensional problem (15) are far more interesting [24, 25, 7]. Indeed, several secondary bifurcations appear, which include symmetry-breaking bifurcations and vortex-bearing states [14]. The bifurcation diagram plots the number of bosons $N$ in the BEC (12) as a function of the chemical potential $\mu$. When $N\to 0$, the nonlinearity of the problem becomes irrelevant and the states bifurcate from the respective linear limit. Starting from this low-density context, we are interested in exploring the solution modes for greater values of $\mu$, which make the problem strongly nonlinear. Since an arbitrary potential can be approximated as a harmonic potential at the vicinity of a stable equilibrium point, when $N\to 0$ we can decompose the linear eigenfunction $\phi_{m,n}$ in Cartesian form as being proportional to

[TABLE]

where $H_{j}$ is the Hermite polynomial with $j$ being the associated quantum number of the harmonic oscillator. The critical value of the eigenvalue corresponding to linear eigenfunction $\phi_{m,n}$ is $\mu_{crit}=E_{m,n}:=(m+n+1)\Omega$. Notice that the characteristic eigenvalue parameter, i.e. the eigenvalue responsible for the bifurcation, is the chemical potential. Thus, given an initial energy $\mu$ at the linear limit, we increase the chemical potential (and therefore the number of atoms $N$) in order to approach to the strongly nonlinear regime that can lead to the discovery of new states originating from secondary bifurcations.

For the numerical characterization of the stability (and possible classification of the instability) for each state, we refer to [7].

3.1 Finite Element discretization

In this subsection, we apply the method presented in Sec. 2.2 to problem (15). We recall that the solution $\phi$ to eq. (15) is a complex function. Let $\varphi$ and $\psi$ be its real and imaginary part, respectively. Let us introduce:

[TABLE]

The generic $k$-th iteration of the Newton’s method (4) reads: seek $\delta X_{h}=(\delta\varphi_{h},\delta\psi_{h})\in V_{h}$ with $V_{h}\subset(H_{0}^{1})^{2}$, such that

[TABLE]

Note that for simplicity of notation, we have not specified that the solution $X_{h}$ depends on parameter $\mu$.

Let us introduce the finite element discretization matrices:

[TABLE]

The $\mathbb{R}^{N_{h}\times N_{h}}$ Jacobian matrix $\mathbb{J}(\vec{X}_{h}^{k}(\mu);\mu)$ can be written as

[TABLE]

Next, we will apply the ROM technique presented in Sec. 2.3.

3.2 Application of the RB method

We present the problem that has to be solved in the online phase of the RB method applied to the Gross–Pitaevskii equation. Given $\mu\in\mathcal{D}$, at every iteration $k$ of Newton’s method (7) we seek $\delta X_{N}=(\delta\phi_{N},\delta\psi_{N})\in V_{N}$ such that

[TABLE]

$\forall\,Y_{N}\in V_{N}$.

The reduced Jacobian $\mathbb{J}_{N}(\vec{X}_{N}^{k}(\mu);\mu)\in\mathbb{R}^{N\times N}$ can be written as

[TABLE]

where:

[TABLE]

are the reduced matrices written in terms of finite element matrices (18) and transformation matrix (10).

In general, the time saving promised by the online-offline strategy are enabled by the so called affine decomposition [18], which makes the computations in the online phase independent form the (usually very high) number of degrees of freedom $N_{h}$. Here, however, reduced matrix $\mathbb{C}_{N}$ introduced above depends on $\mu$ through the solution computed at each step of Newton’s method. Thus, we will need an affine-recovery technique called Empirical Interpolation Method (EIM) [4] in order to obtain substantial savings of the computational time during the online phase, as demonstrated in Sec. 4.3.

4 Results

The high fidelity (or full order) approximations used for the results presented in this section were computed with FEniCS [13, 23, 3], while we used RBniCS [38] for the reduced order approximations.

4.1 Validation of the full order method

To validate the Full Order Method for the Gross–Pitaevskii as described in Sec. 3.1, we consider a test proposed in [7]. We approximate the solution to eq. (15) in domain $D=(-12,12)^{2}$, with homogeneous Dirichlet boundary conditions on the entire boundary of $D$. We recall that we set $\Omega=0.2$ and $\mathcal{D}=[0,1.2]$. For the space discretization, we use $\mathbb{P}_{2}$ finite elements and a mesh with 6889 elements.

Fig. 1 displays the FOM bifurcation diagram in the $\mu$-$N$ plane (left) and in the $\mu$-$||\rho||_{\infty}$ plane (right).

The full reconstruction of the bifurcation diagram requires a proper initialization of our algorithm, in particular for Newton’s method. For this purpose, we rely on the linear limit of the system as specified in (16). Indeed, to approximate each branch in Fig. 1, we assign an initial guess proportional to the product of Hermite polynomials $H_{m}$ and $H_{n}$, where $m$ and $n$ chosen according to the value of the critical point $\mu_{crit}=(m+n+1)\Omega$. In the general case (i.e., no analytic information is available), one can recover such an initial guess either from a linearized eigenvalue problem [31] or thorough a deflation method [7, 32].

These diagrams show the first three bifurcation points and the relative non-uniqueness of the solution with respect to the parameter $\mu$. As $\mu$ is increased, the sequence of events is as follows. The ground state $|0,0\rangle$ is the system simplest state. Its linear eigenfunction $\phi_{0,0}$ has corresponding eigenvalue $\mu=\Omega$. The ground state is generically stable, thus no further bifurcations occur from this state [19]. As expected, a unique solution branch departs from $\mu=\Omega$ in Fig. 1. A representative density function for this branch is shown in Fig. 2(a). We see no further bifurcation for $\Omega\leq\mu<2\Omega$. The first interesting events in terms of bifurcation analysis occur for $\mu=2\Omega$, with $n+m=1$: two branches, associated to $|0,1\rangle$ and $|1,0\rangle$, bifurcate from point $(2\Omega,0)$ in the $\mu$-$N$ and $\mu$-$||\rho||_{\infty}$ planes [24, 12]. Indeed, from point $(2\Omega,0)$ in Fig. 1 we observe the two expected branches. Representative density functions for these two branches are reported in Fig. 2(b) and 2(c). The next, more complicated, case of bifurcations emanates from point $(3\Omega,0)$, with $n+m=2$. In Fig. 1, we see that three branches depart from this point, associated to $|1,1\rangle$, $|0,2\rangle$, and $|2,0\rangle$. The corresponding representative densities are shown in Fig. 2(f), 2(e), and 2(d). Finally, all the points without marker in Fig. 1 correspond to the non-physical solution $\phi=0$ that exists since there are no external forces in eq. (15).

Fig. 2 displays the density functions associated to $\mu=1.2$ and all the 6 solution branches in Fig. 1. We observe the richness of density patterns in order of decreasing $N$. In particular, we see the ground state $|0,0\rangle$ in Fig. 2(a), the single charge vortex $|0,1\rangle$ in Fig. 2(b), the 1-dark soliton stripe $|1,0\rangle$ in Fig. 2(c), the dark soliton cross $|1,1\rangle$ in Fig. 2(f), the ring dark soliton $|0,2\rangle$ in Fig. 2(e), the 2-dark soliton stripe $|2,0\rangle$ in Fig. 2(d). Notice that the 6 branches Fig. 1 are related to the first three eigenvalues. For example, the second bifurcation stems from a double eigenvalue and thus we have two branches. This phenomenon is called multiple bifurcations. The stability property of these branches are different for each case, i.e the single charge vortex is always stable while the 1-dark soliton stripe is subject to multiple secondary bifurcations. These properties can be easily studied using standard techniques (see, e.g., [7]). The results in Fig. 1 and 2 are in excellent agreement with the results reported in [7], indicating that the mesh that we use is sufficiently refined for this study.

The overall simulation time required to complete the diagrams in Fig. 1 is roughly 96 minutes with continuation step $\Delta\mu=1.25\cdot 10^{-3}$.

4.2 Tracing bifurcation diagrams with the reduced order approach

4.2.1 One-parameter study

In this section, we present the results obtained with our Reduced Order Method as described in Sec. 2.3 and compare them with the FOM results reported in Sec. 4.1. The only parameter that varies is $\mu$ in interval $[0,1.2]$.

Concerning the construction of the reduced manifold, we employed a training set for the POD with cardinality $N_{train}=160$ for each one of the six branches. Setting the POD tolerance to $10^{-9}$, we obtain a global basis of dimension $N=51$. In the online phase, we reconstruct the reduced bifurcation diagram for all the 961 equally spaced points in $\mathcal{D}=[0,1.2]$ used in the high-fidelity bifurcation diagram shown in Fig. 1. Such points correspond to continuation step $\Delta\mu=1.25\cdot 10^{-3}$.

Fig. 3 shows reduced order errors

[TABLE]

i.e. the difference in absolute value between the branches of the bifurcation diagram computed with FOM and ROM in the $\mu$-$N$ plane (top) and in the $\mu$-$||\rho||_{\infty}$ plane (bottom). In Fig. 3, we see that the largest peaks are associated to the $|1,1\rangle$ branch at $\mu=0.6$. In general, it is expected to have larger errors at the bifurcation points where differentiability with respect to the parameter $\mu$ is lost. We infer that the errors are largest at $\mu=0.6$ due to the more complicated solution structure (compare Fig. 2(f) to the other panels in Fig. 2). From Fig. 3 (bottom), we see the largest error is of the order of $10^{-4}$ also for the infinity norm of the density. However, we observe larger errors over interval $[0.2,1.2]$ for $\mu$, as opposed to localized at $\mu=0.6$ as in Fig. 3 (top).

As further evidence of the accuracy of our ROM approach, we plot in Fig. 4 the difference between $X_{h}$ (FOM solution) and $X_{N}$ (ROM solution) in the $L^{2}$ and $H^{1}_{0}$ norms. Recall that $X_{h}$ and $X_{N}$ consists of both real and imaginary part of the computed solution of the Gross–Pitaevskii equation. As expected, for all the branches the largest errors in both norms occur at the point where each branch departs from the horizontal axis, i.e. $\mu=0.2$ for branch $|0,0\rangle$, $\mu=0.4$ for branches $|1,0\rangle$ and $|0,1\rangle$, and $\mu=0.6$ for branches $|1,1\rangle$, $|0,2\rangle$, and $|2,0\rangle$. In addition, just like in Fig. 3 the largest errors are associated with branch $|1,1\rangle$.

Fig. 5 reports the difference between the density function $\rho$ computed with FOM and ROM in the $L^{2}$ and $H^{1}_{0}$ norms. In the case of the density, the largest errors for each branch occur for $\mu$ larger than the critical value where the branch departs from the horizontal axis. For a better understanding of how the ROM density compares with the FOM density, Fig. 6 displays the difference $\rho_{h}-\rho_{N}$ for $\mu=1.2$. We observe larger errors for the 1- and 2-dark soliton stripe and the dark soliton cross, i.e. for the solutions that do not have central symmetry.

Fig. 3-6 show the ability of our ROM approach to accurately reconstruct bifurcation diagrams as parameter $\mu$ varies. Clearly, it makes sense to set up the ROM machinery if there is a substantial gain in terms of computational time. Because of the nonlinearity in the Gross–Pitaevskii equation, which makes the computations in the online phase dependent on the number of FEM degrees of freedom, the computational speed-up enabled by our ROM approach is only 1.1: it took 86 minutes to generate the data needed for the bifurcation diagrams, using continuation step $\Delta\mu=1.25\cdot 10^{-3}$, accounting only for the online phase computations. Recall that it takes 96 minutes with FOM. These computational savings are not satisfactory, especially if we were to include the cost for the offline phase. Before introducing an affine recovery technique to drastically improve computational efficiency, we present a two-parameter study.

4.2.2 Two-parameter study

In this section, we plot the bifurcation diagram as chemical potential $\mu$ varies in interval $[0,1.2]$ and trap strength $\Omega$ varies in interval $[0.1,0.3]$.

First, we focus on the first bifurcation, i.e. branch $|0,0\rangle$. Fig. 7 show such branch in a two-parameter bifurcation diagram obtained with the Reduced Order Model as $\mu$ and $\Omega$ are varied. As expected from the theory, we see that as $\Omega$ increases the critical value of $\mu$ for the first bifurcation increases linearly. Recall that $\mu_{crit}=\Omega$, for $m=n=0$, and Fig. 7 clearly shows it (see the black dotted line in the highlighted red rectangle). For Fig. 7, we used increment $\Delta\Omega=0.01$ and continuation step $\Delta\mu=1.25\cdot 10^{-3}$. The online phase computations for the graph in Fig. 7 took roughly 164 minutes. The corresponding time required by the FOM is approximately 249 minutes. So with our ROM approach we obtain a speed-up of 1.5, which represents an improvement over the 1.1 speed-up for the one parameter study but it is still not enough to justify the computational costs of the offline phase.

Next, we focus on the first two bifurcations but restrict the attention to $\Omega=0.1,0.2,0.3$. Fig. 8 shows the first three branches, i.e. branches $|0,0\rangle$, $|0,1\rangle$, and $|1,0\rangle$. We see that our ROM approach successfully captures also the critical $\mu$ for the second bifurcation: branches $|0,1\rangle$ and $|1,0\rangle$ depart from $\mu_{crit}=2\Omega$. Just like for Fig. 7, we used continuation step $\Delta\mu=1.25\cdot 10^{-3}$. The time required by the online phase computations to plot the graph in Fig. 8 is roughly 245 minutes.

4.3 Hyper-reduction techniques

In this section, we use two affine-recovery techniques to reconstruct branch $|0,1\rangle$ for the one parameter study. These techniques are called Empirical Interpolation Method (EIM) [4] and Discrete Empirical Interpolation Method (DEIM) [8] and they makes the online phase computations independent from the number of degrees of freedom of the finite element method.

As we have remarked earlier, if the affine dependence assumption is not fulfilled, the speedup of the online reconstruction can be limited, thereby compromising the whole methodology. With the hyper-reduction techniques one aims at approximating a general parametrized function $g:D\times\mathcal{D}\to\mathbb{R}$ by a sum of affine terms:

[TABLE]

where the basis functions $h_{q}$ are obtained by means of a linear combination of $Q$ snapshots $\{g_{\mu_{q}}\}_{q=1}^{Q}$ and the sample points are chosen through a Greedy approach. In (22), the coefficients of such expansion $c_{q}(\mu)$ are found by solving

[TABLE]

in some particular points $\{\mathbf{x}_{j}\}_{j=1}^{Q}$ of the domain $D$ called magic points. Hence, Empirical Interpolation strategies provide a discrete version of $\texttt{I}[g_{\mu}](x)$ as

[TABLE]

The main difference between EIM and DEIM is in the construction of $\mathbb{H}$. In fact, EIM embeds the construction of the basis inside the Greedy procedure, while DEIM exploits a POD on a set of snapshots. Moreover, the DEIM strategy starts with discretizing the nonlinearity, while EIM constructs the set of the magic points and the basis functions before the discretization step.

Let us express residual (4) as

[TABLE]

where parameter affine dependency is guaranteed by the forms $\theta_{q}^{g}(\mu)$. The reduced residual (11) can then be expressed as:

[TABLE]

The hyper-reduction techniques described above provide the following affine approximation of the reduced residual

[TABLE]

where $\{h^{q}\}_{q=1}^{Q_{g}}$ represent a suitable basis and $c_{q}$ the interpolation coefficients. The Jacobian matrix $\mathbb{J}_{N}(\vec{X}_{N};\mu)$ is assembled in a similar fashion.

The purpose of this section is two-fold: show that the branch is reconstructed accurately and report on the substantial computational time savings enable by the hyper-reduction strategies. Fig. 9 shows reduced order errors $E_{N}$ (right) and $E_{\rho}$ (right) for branch $|0,1\rangle$. We see that the values of $E_{N}$ for EIM and DEIM are comparable over the entire $[0,1.2]$ interval, with a peak at $\mu=1.2$ of the order of $10^{-5}$. Error $E_{\rho}$ is slightly larger for EIM, with a peak of the order of $10^{-4}$.

Finally, Fig. 10 shows the difference between FOM solution and ROM solution computed with EIM (left) and DEIM (right) in the $L^{2}$ and $H^{1}_{0}$ norms, again for branch $|0,1\rangle$. We observe slightly larger error peaks than in the case of ROM with no affine-recovery technique: compare Fig. 10 with Fig. 4(b). These increased errors are the price to pay for a considerable computational speed-up. With EIM, it takes only 55 s to construct branch $|0,1\rangle$ while it takes 246 s with FOM. So, our ROM approach coupled with EIM is almost five times faster than the FOM. As for DEIM, the computational time savings are even better: it takes only 7 s for the construction of branch $|0,1\rangle$, with corresponds to a speed-up factor of almost 32. This drastic reduction of the computational time allowed by the DEIM is expected [8].

The hyper-reduction techniques is especially effective in reducing the computational time needed for the two-parameter study presented in Sec. 4.2.2. Our ROM approach with DEIM takes only 6 minutes to reconstruct the graph in Fig. 8. This is a speed-up factor of 40 with respect to our ROM approach without DEIM and a factor 60 with respect to the FOM.

5 Conclusions and perspectives

We have exploited reduced order methods to drastically reduce the computational time required to trace a bifurcation diagram. We proposed a combination of different techniques to overcome the curse of dimensionality, mostly when more than one physical parameter varies. In particular we used the standard Finite Element method as an high fidelity solver for the offline phase and a Proper orthogonal decomposition as a reduction technique to construct a basis for the approximation manifold. Moreover, we implemented a simple continuation method to prevent the divergence of the Newton-Kantorovich method.

To demonstrate the effectiveness of our approach, we studied bifurcating phenomena in quantum mechanics described by the Gross-Pitaevskii equation. We were able to trace the full bifurcation diagram with high accuracy even when dealing with a multi-parameter context. We also highlighted the need for an affine recovery technique (EIM/DEIM) in order to obtain important computational time savings.

There are several ways in which this work could be expanded. One can consider different reduction strategies, such as POD-Greedy for the 2-parameters test case, in order to understand how the bifurcating parameter varies and then discover new branches at a reduced computational cost. In addition, a deflation method could be implemented and paired with a more involved and smart continuation method. Finally, an extension to the 3D version of the same model or to multi-component systems could be investigated.

Acknowledgements

This work was supported by European Union Funding for Research and Innovation through the European Research Council (project H2020 ERC CoG 2015 AROMA-CFD project 681447, P.I. Prof. G. Rozza). This work was also partially supported by NSF through grant DMS-1620384 (Prof. A. Quaini). We acknowledge fruitful conversations with Dr. S. Garner. We thank Dr. F. Ballarin (SISSA) for his great help with the RBniCS software library and precious discussion.