Model order reduction for stochastic dynamical systems with continuous symmetries

TL;DR

This paper introduces a new order reduction method for stochastic dynamical systems with continuous symmetries, combining symmetry reduction with standard techniques to efficiently capture coherent patterns.

Contribution

The authors develop a novel symmetry-reduced Dynamically Orthogonal (SDO) scheme that improves efficiency in modeling stochastic systems with symmetries.

Findings

SDO scheme effectively captures patterns with fewer modes

Demonstrated on 1D Korteweg-de Vries and 2D Navier-Stokes equations

Outperforms traditional methods in symmetry-invariant stochastic systems

Abstract

Stochastic dynamical systems with continuous symmetries arise commonly in nature and often give rise to coherent spatio-temporal patterns. However, because of their random locations, these patterns are not well captured by current order reduction techniques and a large number of modes is typically necessary for an accurate solution. In this work, we introduce a new methodology for efficient order reduction of such systems by combining (i) the method of slices, a symmetry reduction tool, with (ii) any standard order reduction technique, resulting in efficient mixed symmetry-dimensionality reduction schemes. In particular, using the Dynamically Orthogonal (DO) equations in the second step, we obtain a novel nonlinear Symmetry-reduced Dynamically Orthogonal (SDO) scheme. We demonstrate the performance of the SDO scheme on stochastic solutions of the 1D Korteweg-de Vries and 2D Navier-Stokes equations.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}=\mathbf{F}(\mathbf{u},t;\omega ),$$

$$⟨{\mathbf{u}}_1,{\mathbf{u}}_2⟩={\int }_D{\mathbf{u}}_1{\mathbf{u}}_2^{\ast }\mathrm{d}\mathbf{x}.$$

$$\mathbf{F}(g\mathbf{u})=g\mathbf{F}(\mathbf{u}).$$

$$\underset{\mathbf{\varphi }}{\min }||g^{-1}(\mathbf{\varphi })\mathbf{u}-{\hat{\mathbf{u}}}^{\prime }||,$$

$$\underset{\mathbf{\varphi }}{\min }||\mathbf{u}-g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime }||,$$

$$\frac{\partial }{\partial {\varphi }_a}||\mathbf{u}-g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime }|{|}^2=\frac{\partial }{\partial {\varphi }_a}⟨\mathbf{u},\mathbf{u}⟩-2\frac{\partial }{\partial {\varphi }_a}⟨\mathbf{u},g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime }⟩+\frac{\partial }{\partial {\varphi }_a}⟨g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime },g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime }⟩=0,$$

$$⟨\mathbf{u},{\mathbf{t}}_a(g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime })⟩=0,$$

$${\mathbf{t}}_a(g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime })=\frac{\partial g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime }}{\partial {\varphi }_a}=\underset{\delta {\varphi }_a\to 0}{\lim }\frac{g(\mathbf{\varphi }+\delta {\varphi }_a){\hat{\mathbf{u}}}^{\prime }-g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime }}{\delta {\varphi }_a}.$$

$${\mathbf{t}}_a(g(\mathbf{\varphi }){\hat{\mathbf{u}}}^{\prime })=g(\mathbf{\varphi }){\mathbf{t}}_a({\hat{\mathbf{u}}}^{\prime }),$$

$${\mathbf{t}}_a({\hat{\mathbf{u}}}^{\prime })=\underset{\delta {\varphi }_a\to 0}{\lim }\frac{g(\delta {\varphi }_a){\hat{\mathbf{u}}}^{\prime }-{\hat{\mathbf{u}}}^{\prime }}{\delta {\varphi }_a}.$$

$$⟨\mathbf{u},g(\mathbf{\varphi }){\mathbf{t}}_a^{\prime }⟩=0\iff ⟨\hat{\mathbf{u}},{\mathbf{t}}_a^{\prime }⟩=0,$$

$${\dot{\varphi }}_a\frac{\partial g(\mathbf{\varphi })\hat{\mathbf{u}}}{\partial {\varphi }_a}+g(\mathbf{\varphi })\frac{\partial \hat{\mathbf{u}}}{\partial t}=\mathbf{F}(g(\mathbf{\varphi })\hat{\mathbf{u}}),$$

$$\frac{\partial g(\mathbf{\varphi })\hat{\mathbf{u}}}{\partial {\varphi }_a}={\mathbf{t}}_a(g(\mathbf{\varphi })\hat{\mathbf{u}})=g(\mathbf{\varphi }){\mathbf{t}}_a(\hat{\mathbf{u}}),$$

$${\mathbf{t}}_a(\hat{\mathbf{u}})=\underset{\delta {\varphi }_a\to 0}{\lim }\frac{g(\delta {\varphi }_a)\hat{\mathbf{u}}-\hat{\mathbf{u}}}{\delta {\varphi }_a}.$$

$$\frac{\partial \hat{\mathbf{u}}}{\partial t}=\mathbf{F}(\hat{\mathbf{u}})-{\dot{\varphi }}_a{\mathbf{t}}_a(\hat{\mathbf{u}}).$$

$$⟨\frac{\partial \hat{\mathbf{u}}}{\partial t},{\mathbf{t}}_a^{\prime }⟩=⟨\mathbf{F}(\hat{\mathbf{u}}),{\mathbf{t}}_a^{\prime }⟩-{\dot{\varphi }}_b⟨{\mathbf{t}}_b(\hat{\mathbf{u}}),{\mathbf{t}}_a^{\prime }⟩=0.$$

$$\dot{\mathbf{\varphi }}={\mathbf{T}}^{-1}\mathbf{f}.$$

$$\hat{\mathbf{u}}(\mathbf{x},t;\omega )=\bar{\mathbf{u}}(\mathbf{x},t)+Y_i(t;\omega ){\hat{\mathbf{u}}}_i(\mathbf{x},t),i=1,...,s,$$

$$⟨\frac{\partial {\hat{\mathbf{u}}}_i}{\partial t},{\hat{\mathbf{u}}}_j⟩=0,i,j=1,...,s,$$

$$\frac{\partial \bar{\mathbf{u}}}{\partial t}=E[\mathbf{F}(\hat{\mathbf{u}})-{\dot{\varphi }}_a{\mathbf{t}}_a(\hat{\mathbf{u}})],$$

$$\frac{\mathrm{d}{Y}_i}{\mathrm{d}t}=⟨\mathbf{F}(\hat{\mathbf{u}})-{\dot{\varphi }}_a{\mathbf{t}}_a(\hat{\mathbf{u}})-E[\mathbf{F}(\hat{\mathbf{u}})-{\dot{\varphi }}_a{\mathbf{t}}_a(\hat{\mathbf{u}})],{\hat{\mathbf{u}}}_i⟩,$$

$$\frac{\partial {\hat{\mathbf{u}}}_i}{\partial t}={\hat{\mathbf{H}}}_i-⟨{\hat{\mathbf{H}}}_i,{\hat{\mathbf{u}}}_j⟩{\hat{\mathbf{u}}}_j,$$

$${\hat{\mathbf{H}}}_i=E[Y_k(\mathbf{F}(\hat{\mathbf{u}})-{\dot{\varphi }}_a{\mathbf{t}}_a(\hat{\mathbf{u}}))]C_{ik}^{-1},$$

$$\dot{\mathbf{\varphi }}={\mathbf{T}}^{-1}\mathbf{f},$$

$$\mathbf{u}(\mathbf{x},t;\omega )=g(\mathbf{\varphi }(t;\omega ))\hat{\mathbf{u}}(\mathbf{x},t;\omega )=g(\mathbf{\varphi }(t;\omega ))\bar{\mathbf{u}}(\mathbf{x},t)+Y_i(t;\omega )g(\mathbf{\varphi }(t;\omega )){\hat{\mathbf{u}}}_i(\mathbf{x},t).$$

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\mu \frac{{\partial }^{3}u}{\partial x^3}=0,$$

$$g(c)u(x)=u(x-c),$$

$$t(u)=\underset{\delta c\to 0}{\lim }\frac{g(\delta c)u-u}{\delta c}=-\frac{\partial u}{\partial x}.$$

$$\frac{\partial \hat{u}}{\partial t}=F(\hat{u})+\dot{c}\frac{\partial \hat{u}}{\partial x},$$

$$F(u)=-u\frac{\partial u}{\partial x}-\mu \frac{{\partial }^{3}u}{\partial x^3},$$

Full text

1 Introduction

Examples of physical systems that can be modeled as stochastic dynamical systems with continuous symmetries abound in the world around us, whether it be pipe flow [15], water waves [36], flame dynamics [27], or nonlinear optics [1], just to name a few. Stochasticity arises from unknown parameters or initial conditions, while continuous symmetry manifests itself as translational and/or rotational invariance due to specific geometry. Such systems are typically described by stochastic partial differential equations (PDEs) with complex nonlinear responses, which makes accurate quantification of their statistical behavior through direct Monte-Carlo simulations a challenge due to the high computational costs involved.

Model order reduction aims at solving this issue by approximating the stochastic solution in terms of a finite sum of deterministic spatial modes multiplied by stochastic scalar coefficients, motivated by the Karhunen-Loève decomposition. In this way, the computation of the stochastic solution is reduced to the evolution of the coefficients and/or the modes, leading to significant computational savings and, in some circumstances, to improved physical understanding of the underlying dynamics. Naturally, there red different possible ways to derive such reduced-order models from the governing equations, and various dimensionality reduction methods have been proposed over the years, such as the Proper Orthogonal Decomposition (POD) combined with Galerkin projection [3, 16], the Polynomial Chaos (PC) expansion [35], or more recently the Dynamically Orthogonal (DO) equations [30].

Stochastic dynamical systems with continuous symmetries often give rise to coherent spatio-temporal patterns [9]. However, due to the fact that individual realizations are invariant along the symmetry directions of the system, the precise location of these spatio-temporal patterns might be subject to large stochastic variability and different realizations might display similar structures at completely different spatial locations. As a result, reduced-order models relying on linear modal decompositions will require a large number of spatial modes to adequately capture the spatio-temporal dynamics of the stochastic solution, which is likely to offset the computational gains. This is because spatial shifts cannot be efficiently represented by a finite linear combination of global spatial modes.

In parallel with dimensionality reduction methods, much work has been done over the last few decades on symmetry reduction, which concerns the removal of continuous symmetries associated with deterministic dynamical systems [8, 11, 28, 4, 33, 13, 20]. After symmetry reduction, the dynamics of the original system along its symmetry directions (i.e. translations and/or rotations) is factored out, so that the resulting symmetry-reduced state is left with non-trivial shape-changing dynamics. Coherent spatio-temporal patterns appearing in the symmetry-reduced solution will therefore remain at the same spatial location while undergoing shape deformations.

In this way, the possibly low-rank structure of the symmetry-reduced state is preserved, which is clearly advantageous for model order reduction purposes. [18, 14] first combined such symmetry reduction techniques with the POD to derive low-dimensional models for deterministic systems governed by PDEs. Their model, however, only described the symmetry-reduced dynamics and no attempt was made to recover the original system state. Such a closure was accomplished shortly thereafter in the same context by [28] with a so-called “reconstruction equation” for the symmetry coordinate, resulting in the first dynamical order reduction framework to take advantage of the continuous symmetries of a system. The same procedure was later adopted by [23] in the context of parametric order reduction using reduced-basis methods. More generally, combinations of nonlinear mappings with reduced-order models have been explored by several authors lately [17, 6, 7, 21, 22].

Similar ideas of reducing the symmetry before performing order reduction have also been applied to stochastic dynamical systems, but these recent studies focused on data reduction [34] or inference [25]. Here, we introduce a new framework for efficient dynamical model order reduction of stochastic dynamical systems with continuous symmetries by combining symmetry reduction with order reduction methods. As we will see, this approach naturally leads to novel nonlinear reduced-order models that efficiently take advantage of the continuous symmetries of the system, leading to much improved accuracy for a given number of modes. This methodology can be applied to any dimensionality reduction method of choice provided the latter preserves the symmetry reduction properties, which is true for techniques like the POD or the DO equations. We will illustrate our approach with the DO equations, which will result in a novel Symmetry-reduced Dynamically Orthogonal (SDO) scheme.

The paper is structured as follows. The general symmetry-dimensionality reduction methodology as well as the derivation of the SDO scheme are presented in Section 2. The performance of the SDO scheme is compared with the standard DO method on stochastic simulations of the 1D Korteweg-de Vries and 2D Navier-Stokes equations in Sections 3 and 4, respectively. Finally, conclusions follow in Section 5.

2 Blending symmetry reduction with dimensionality reduction

Let $(\Omega,\mathcal{B},\mathcal{P})$ be a probability space and $\omega\in\Omega$ indicate an elementary event. Denoting space $\mathbf{x}\in D\subset\mathbb{R}^{n}$ and time $t$, we consider the stochastic partial differential equation

[TABLE]

where $\mathbf{F}$ is a (possibly stochastic and time-dependent) nonlinear differential operator and $\mathbf{u}(\mathbf{x},t;\omega)$ is a random vector field that belongs to the Hilbert space $\mathcal{H}$ of continuous and square-integrable functions with inner product

[TABLE]

In this work, we are interested in dynamical systems that are symmetric under a group $G$ of continuous transformations, that is, the differential operator $\mathbf{F}$ satisfies the following equivariance condition for any group element $g\in G$

[TABLE]

We will restrict ourselves to symmetry groups $G$ that (i) are Lie groups, (ii) preserve inner products and distances. Typically, $G$ will consist of translations along different direction and/or rotations about different axes. When the above relation (3) is satisfied, individual solutions $\mathbf{u}$ of the dynamical system (1) are invariant under $G$, meaning that $g\mathbf{u}$ is also a solution for any $g\in G$.

In this paper, we introduce a new framework for efficient dimensionality reduction of such systems with continuous symmetries. Our methodology comprises the two following steps:

In a first step, we perform symmetry reduction of the dynamical system (1) using the method of slices [28, 13]. In this framework, the original system state $\mathbf{u}$ is decomposed into (i) a stochastic symmetry-reduced state $\hat{\mathbf{u}}$, which is fixed in the physical domain but captures the intrinsic changes in shape of $\mathbf{u}$, and (ii) a finite set of stochastic phase parameters $\boldsymbol{\phi}$ that track the motion of $\mathbf{u}$ along the symmetry directions of the system. The symmetry-reduced state $\hat{\mathbf{u}}$ and phase parameters $\boldsymbol{\phi}$ are defined precisely in Section 2.1, and equations governing their temporal evolution are presented in Section 2.2. 2. 2.

In a second step, we perform dimensionality reduction of the stochastic symmetry-reduced state using any standard model order reduction method of choice. This order reduction step is illustrated in Section 2.3 through the use of the Dynamically Orthogonal (DO) equations [30], which leads to a novel, nonlinear Symmetry-reduced Dynamically Orthogonal (SDO) scheme.

In the following, we adopt a dynamical state-space approach where the system state $\mathbf{u}(t;\omega)$ for given time $t$ and realization $\omega$ is represented by a single point in the infinite dimensional state space $\mathcal{M}$ of all possible solutions.

2.1 Method of slices

Let us first consider a given deterministic state $\mathbf{u}$, for instance a particular realization $\omega_{0}$ at a given time $t_{0}$ of (1). The application of the family of continuous transformations $g\in G$ to $\mathbf{u}$ gives rise to a family of dynamically equivalent states $g\mathbf{u}$ called the group orbit of $\mathbf{u}$, as illustrated in Figure 1(a).

The goal of symmetry reduction is to reduce all these equivalent symmetry copies to a unique representative symmetry-reduced state $\hat{\mathbf{u}}$. For a given state $\mathbf{u}$ and its group orbit, the method of slices defines $\hat{\mathbf{u}}$ by choosing the point on the group orbit of $\mathbf{u}$ that is closest to a fixed template state $\hat{\mathbf{u}}^{\prime}$, so that $\hat{\mathbf{u}}$ overlies as well as possible the template in physical space, as indicated in Figure 1(b). Denoting with $\boldsymbol{\phi}=\phi_{1},...,\phi_{N}$ the $N$ continuous scalar phase parameters of the transformation $G$ (for example the translations amounts and/or rotation angles along different directions) and writing $\mathbf{u}=g(\boldsymbol{\phi})\hat{\mathbf{u}}$ or equivalently $\hat{\mathbf{u}}=g^{-1}(\boldsymbol{\phi})\mathbf{u}$, this condition can be expressed in terms of $\boldsymbol{\phi}$ as

[TABLE]

where we use the $L^{2}$ norm $||\mathbf{u}||^{2}=\langle\mathbf{u},\mathbf{u}\rangle$. As mentioned earlier, we only consider transformations that do not affect the inner product between different states, i.e. $\langle g\mathbf{u}_{1},g\mathbf{u}_{2}\rangle=\langle\mathbf{u}_{1},\mathbf{u}_{2}\rangle$ and $||g\mathbf{u}||=||\mathbf{u}||$. Therefore, the above condition becomes

[TABLE]

from which we can deduce the extremum condition

[TABLE]

where $a=1,...,N$. Strictly speaking, one should also make sure that the second derivative with respect to $\phi_{a}$ be positive in order to have a minimum. Using the distance-preserving property of the transformation, the extremum condition leads to the following slice condition

[TABLE]

where $\mathbf{t}_{a}(g(\boldsymbol{\phi})\hat{\mathbf{u}}^{\prime})$ is the tangent to the group orbit at $g(\boldsymbol{\phi})\hat{\mathbf{u}}^{\prime}$ in direction $\phi_{a}$

[TABLE]

Since $G$ is a Lie group, we have $g(\boldsymbol{\phi}+\delta\phi_{a})\hat{\mathbf{u}}^{\prime}=g(\boldsymbol{\phi})g(\delta\phi_{a})\hat{\mathbf{u}}^{\prime}$, thus we can factor out the group action $g(\boldsymbol{\phi})$ from the above expression to obtain the relation

[TABLE]

where $\mathbf{t}_{a}(\hat{\mathbf{u}}^{\prime})$ is the group orbit tangent at the fixed template $\hat{\mathbf{u}}^{\prime}$ in direction $\phi_{a}$

[TABLE]

As a result, we may now use the distance-preserving property of $g$ to rewrite the slice condition (7) in terms of the symmetry-reduced state $\hat{\mathbf{u}}$

[TABLE]

where we have denoted the fixed template tangent $\mathbf{t}_{a}^{\prime}=\mathbf{t}_{a}(\hat{\mathbf{u}}^{\prime})$. Given a state $\mathbf{u}$ and template $\hat{\mathbf{u}}^{\prime}$, the first equality in the above slice condition gives the phase parameters $\boldsymbol{\phi}$ such that the distance between the symmetry-reduced state $\hat{\mathbf{u}}=g^{-1}(\boldsymbol{\phi})\mathbf{u}$ and the template $\hat{\mathbf{u}}^{\prime}$ is minimized. The second equality can be interpreted as an orthogonality condition which states that the symmetry-reduced state $\hat{\mathbf{u}}$ always lies within a hyperplane normal to the $N$ group orbit tangents $\mathbf{t}_{a}^{\prime}$ at the template $\hat{\mathbf{u}}^{\prime}$. This fixed hyperplane thus defines a slice through the full state space containing all symmetry-reduced states and called the symmetry-reduced state space, see Figure 2.

2.2 Dynamics within the symmetry-reduced state space

Consider now that $\mathbf{u}(t;\omega)$ is stochastic and evolves under the governing equation (1), as illustrated in Figure 2. The evolution of $\mathbf{u}(t;\omega)$ in the full state space will correspond a set of stochastic and time-dependent phase parameters $\boldsymbol{\phi}(t;\omega)$ and a stochastic symmetry-reduced state $\hat{\mathbf{u}}(t;\omega)=g^{-1}(\boldsymbol{\phi}(t;\omega))\mathbf{u}(t;\omega)$ evolving in the symmetry-reduced state space in such a way that the slice condition (11) is satisfied at all times and for all realizations. In other words, the symmetry-reduced state is fixed to the template in physical space and captures the shape deformations of $\mathbf{u}(t;\omega)$, while the phase parameters track the motion of $\mathbf{u}(t;\omega)$ along the symmetry directions of the system. In order to find the evolution equations for the symmetry-reduced state $\hat{\mathbf{u}}$ and the phase parameters $\boldsymbol{\phi}$, we first plug the relation $\mathbf{u}=g(\boldsymbol{\phi})\hat{\mathbf{u}}$ into the governing equations

[TABLE]

where repeated indices indicate summation. As was done in (8), we next express the term $\partial g(\boldsymbol{\phi})\hat{\mathbf{u}}/\partial\phi_{a}$ as

[TABLE]

where $\mathbf{t}_{a}(\hat{\mathbf{u}})$ is the group orbit tangent of the time-dependent state $\hat{\mathbf{u}}$ in direction $\phi_{a}$

[TABLE]

Finally, we make use of the equivariance condition (3) to obtain the following evolution equation for the symmetry-reduced state

[TABLE]

It remains to find evolution equations for the phase parameters $\boldsymbol{\phi}$, which is achieved by substituting the above equation into the time derivative of the slice condition (11)

[TABLE]

The above equation can be interpreted as an orthogonality condition that constrains the dynamics of $\hat{\mathbf{u}}$ to remain confined within the symmetry-reduced state space defined by the fixed template $\hat{\mathbf{u}}^{\prime}$. Introducing the stochastic time-dependent matrix $\{\mathbf{T}\}_{ab}=\langle\mathbf{t}_{b}(\hat{\mathbf{u}}),\mathbf{t}_{a}^{\prime}\rangle$ and vector $\{\mathbf{f}\}_{a}=\langle\mathbf{F}(\hat{\mathbf{u}}),\mathbf{t}_{a}^{\prime}\rangle$, one finally gets the following matrix inverse problem for the phase velocity

[TABLE]

Together, equations (15) and (17) govern the evolution of the stochastic symmetry-reduced state $\hat{\mathbf{u}}(t;\omega)$ and the phase parameters $\boldsymbol{\phi}(t;\omega)$ (see also [10] in the deterministic context). In this way, the dynamics of $\mathbf{u}(t;\omega)$ has been separated into shape deformations, reproduced by $\hat{\mathbf{u}}(t;\omega)$ which is fixed in physical space, and motion along the symmetry directions of the system, tracked by $\boldsymbol{\phi}(t;\omega)$. From these two quantities, the evolution of the system in the full state space can be reconstructed exactly as $\mathbf{u}(t;\omega)=g(\boldsymbol{\phi}(t;\omega))\hat{\mathbf{u}}(t;\omega)$, hence no information has been lost so far.

Before moving on to the next step, we note that the choice of template can affect the symmetry reduction procedure in two distinct ways. For a given state $\mathbf{u}$, different templates will naturally lead to different optimal values of $\boldsymbol{\phi}$ for which the symmetry-reduced state $\hat{\mathbf{u}}=g^{-1}(\boldsymbol{\phi})\mathbf{u}$ is closest to the template $\hat{\mathbf{u}}^{\prime}$. Nevertheless, the full state reconstructed from the solution to equations (15) and (17) will be independent of the template. The second and potentially more worrisome effect concerns the conditioning of the matrix $\mathbf{T}$, or worse, the possibility of having singularities in equation (17) whenever the determinant of $\mathbf{T}$ vanishes. This is a consequence of the fact that there exist states $\mathbf{u}$ for which a solution to the minimization problem (4), or equivalently the slice condition (11), ceases to exist (a trivial exemple being the spatially constant state). The number and nature of such problematic states depend on the template function and reflect the finite extent of validity of the symmetry-reduced state space defined by a given template. Therefore, the well-posedness of (17) and the likelihood of encountering singularities are critically tied to the choice of template. We will see in Sections 3.3 and 4.3 that there exist clever choices which almost entirely alleviate this issue.

2.3 Order reduction in the symmetry-reduced state space

We now turn to the second step of our methodology, which consists in applying standard order reduction methods directly to the dynamics of $\hat{\mathbf{u}}$ in the symmetry-reduced space, governed by equations (15) and (17). Because the symmetry-reduced stochastic state is defined in such a way that the distance $||\hat{\mathbf{u}}(\mathbf{x},t;\omega)-\hat{\mathbf{u}}^{\prime}(\mathbf{x})||$ is minimized at all times and for all realizations, the symmetry-reduced realizations $\hat{\mathbf{u}}(\mathbf{x},t;\omega)$ will be grouped together in physical space, as illustrated in Figure 3.

This property is very appealing for low-dimensional models that rely on linear modal decompositions since it implies that the symmetry-reduced realizations $\hat{\mathbf{u}}(\mathbf{x},t;\omega)$ can be approximated much better by such finite-dimensional modal decompositions than their full state space counterparts $\mathbf{u}(\mathbf{x},t;\omega)$. This is inherently connected to the notion of Kolmogorov n-width from approximation theory [19, 24], and we may say that the symmetry-reduced realizations have much lower n-width than their full state space counterparts [23, 6].

Here, as an illustration, we perform dimensionality reduction using the DO framework [30] because of its inherent and desirable ability to deal with strongly transient stochastic responses. Moreover, by considering the DO equations as a computational method for evolving a low-rank matrix representation of the discretized solution, [12] recently showed that DO gives the best possible instantaneous approximation amongst reduced-order models relying on a linear modal decomposition. In general, we note that when choosing an order reduction method, one needs to make sure that the order reduction step preserves the symmetry reduction step, that is, the dynamics of the low-dimensional solution remain confined to the symmetry-reduced state space.

As discussed above, we apply the DO equations directly to the symmetry-reduced state $\hat{\mathbf{u}}$ instead of the original state $\mathbf{u}$, resulting in a new order reduction framework that we call the Symmetry-reduced Dynamically Orthogonal (SDO) equations. Following the standard DO methodology, the stochastic symmetry-reduced state is decomposed into a mean component, and a stochastic part that is projected to a low-dimensional subspace of order $s$ through the following finite-dimensional expansion

[TABLE]

where $\bar{\mathbf{u}}(\mathbf{x},t)$ is the mean of the symmetry-reduced state, $\hat{\mathbf{u}}_{i}(\mathbf{x},t)$ are time-dependent orthonormal modes that capture the principal directions of variance of the symmetry-reduced state, and $Y_{i}(t;\omega)$ are time-dependent stochastic coefficients that characterize the symmetry-reduced stochastic state within this subspace. The redundancy arising from the dependence in time of both the modes and stochastic coefficients is overcome through the DO condition

[TABLE]

which requires that the time variation of the stochastic subspace be orthogonal to itself. Inserting the DO representation (18) into equation (15) governing the evolution of $\hat{\mathbf{u}}(\mathbf{x},t;\omega)$ in the symmetry-reduced state space and using the DO condition (19), one obtains, together with equation (17) for the phase parameters $\boldsymbol{\phi}(t;\omega)$, an explicit set of evolution equations for all unknown quantities. We have the following deterministic PDE for the mean field

[TABLE]

the following stochastic ODEs for the stochastic coefficients

[TABLE]

and the following deterministic PDEs for the modes

[TABLE]

where the deterministic fields $\hat{\mathbf{H}}_{i}$ are defined as

[TABLE]

where $C_{ij}=E[Y_{i}Y_{j}]$ is the covariance matrix. The continuous phase parameters $\boldsymbol{\phi}(t;\omega)$ are evolved according to the stochastic equation (17)

[TABLE]

where $\{\mathbf{T}\}_{ab}=\langle\mathbf{t}_{b}(\hat{\mathbf{u}}),\mathbf{t}_{a}^{\prime}\rangle$ and $\{\mathbf{f}\}_{a}=\langle\mathbf{F}(\hat{\mathbf{u}}),\mathbf{t}_{a}^{\prime}\rangle$. Altogether, equations (20), (21), (22) and (24) constitute the SDO scheme. They entirely specify the time evolution of the symmetry-reduced state $\hat{\mathbf{u}}(\mathbf{x},t;\omega)$ and the phase parameters $\boldsymbol{\phi}(t;\omega)$. In Appendix A, we show that the dynamics of the reduced-order solution (18) remain confined to the symmetry-reduced state space. Finally, the full state is easily recovered through the group transformation

[TABLE]

The above equation shows that the SDO framework can be interpreted as a nonlinear extension of the usual DO methodology, wherein the mean and modes are allowed further degrees of freedom through the stochastic and time-dependent group transformation $g(\boldsymbol{\phi}(t;\omega))$ along the symmetry directions $G$ of the system. In this sense, it is naturally expected that SDO performs as well or better than DO, which we will see in the following sections where we apply both schemes to concrete one and two-dimensional examples.

3 Application to the KdV equation

We first illustrate our approach with the Korteweg-de Vries (KdV) equation that describes the evolution of weakly nonlinear waves on shallow water surfaces

[TABLE]

where $u(x,t)$ is the surface elevation, $x\in[0,L]$ the space variable and $t$ the time. Stochasticity will be introduced through the initial conditions. The boundary conditions are periodic so that the KdV equation is equivariant under the symmetry group of continuous translations, $G=SO(2)_{x}$. The associated shift operator $g\in G$ writes

[TABLE]

where the continuous phase parameter $c$ represents the shift amount. The tangent to the group orbit at an arbitrary state $u$ is then given by

[TABLE]

3.1 Dynamics within the symmetry-reduced state space

As described in Section 2.2, we now pick a fixed template $\hat{u}^{\prime}(x)$ and we consider the evolution of the symmetry-reduced stochastic state $\hat{u}(x,t;\omega)=g^{-1}(c(t;\omega))u(x,t;\omega)$ such that the distance $||\hat{u}(x,t;\omega)-\hat{u}^{\prime}(x)||$ between the reduced state and the template is minimized at all times $t$ and for all realizations $\omega$. In the case of the KdV equation subject to the symmetry group $G$, the dynamical equations (15) for $\hat{u}$ in the symmetry-reduced state space become

[TABLE]

where $F(u)$ is the differential operator corresponding to the KdV equation in the full state space

[TABLE]

and the stochastic differential equation (17) for the shift amount $c(t;\omega)$ writes

[TABLE]

where $t^{\prime}=t(\hat{u}^{\prime})$ is the group orbit tangent to the fixed template $\hat{u}^{\prime}$.

3.2 Order reduction in the symmetry-reduced state space

In order to formulate the reduced-order SDO equations, let us first approximate the symmetry-reduced stochastic state as a truncated Karhunen-Loeve expansion

[TABLE]

where $\bar{u}$ and $\hat{u}_{i}$ are respectively the symmetry-reduced mean and modes and $Y_{i}$ are the stochastic coefficients of the finite-dimensional expansion. Next, we insert the above representation into the RHS of the evolution equation (29) for the symmetry-reduced state

[TABLE]

where $F_{0}$, $F_{i}$, and $F_{ij}$ are deterministic fields given by

[TABLE]

Now, we substitute the expanded RHS operator (33) into the SDO equations (20), (21) and (22) to find the following evolution equation for the mean

[TABLE]

We then have the following evolution equation for the stochastic coefficients

[TABLE]

while the modes are governed by the following equation

[TABLE]

where the deterministic fields $\hat{H}_{i}$ are defined as

[TABLE]

with $M_{mnk}=E[Y_{m}Y_{n}Y_{k}]$ the matrix of third-order moments. Finally, the time evolution of the shift amount is given by the stochastic equation (31), which can be expanded as

[TABLE]

At any given time, the full stochastic state can then be reconstructed from the reduced-order quantities and the time-integrated stochastic shift amount through the following group transformation

[TABLE]

3.3 Choice of the template

The SDO equations become singular as the denominator in (40) vanishes. A critical issue is therefore to choose a template which ensures that the quantity $\langle t(\hat{u}),t^{\prime}\rangle$ always remains nonzero, for all realizations and at all times. For deterministic scalar systems equivariant under continuous translations, this issue was elegantly addressed by [5] with what they called the first Fourier mode slice, a particular choice of template for which the slice condition (7) is equivalent to setting the phase of the first Fourier mode of the symmetry-reduced state $\hat{u}$ to a fixed value. In this way, $\hat{u}$ is pinned at a specific location in physical space and the method is expected to work as long as the amplitude of the first Fourier mode of $u$ does not vanish.

Consider an arbitrary deterministic state $u$ and its symmetry-reduced counterpart $\hat{u}=g^{-1}(c)u$, where the translation amount $c$ is such that the slice condition $\langle\hat{u},t(\hat{u}^{\prime})\rangle=0$ is satisfied given a fixed template function $\hat{u}^{\prime}$. Introducing the Fourier series decomposition

[TABLE]

where $\tilde{u}(k)$, $k\in\mathbb{Z}$ are the Fourier coefficients of $u$, the symmetry-reduced state is expressed as

[TABLE]

where $|\tilde{u}(k)|$ and $\arg\tilde{u}(k)$ are respectively the modulus and phase angle of $\tilde{u}(k)$. The above relation shows that the Fourier coefficients of $\hat{u}$ are those of $u$ rotated by an amount equal to $2\pi kc/L$. In the first Fourier mode slice method, the $c$ is chosen such that the phase angle of the first Fourier mode of $\hat{u}$ is always a fixed value. Here, we choose this value to be $\pi$ so that the symmetry-reduced state appears centered in the domain. This implies $\arg\tilde{u}(1)+2\pi c/L=\pi$ or equivalently $c=-\arg\tilde{u}(1)L/2\pi+L/2$. As shown in Appendix B, this choice of $c$ corresponds to the following template function

[TABLE]

The first Fourier mode slice can be readily applied to our stochastic SDO equations by choosing (44) as our template function. In this way, all symmetry-reduced realizations $\hat{u}(x,t;\omega)$ have their first Fourier mode phase fixed to $\pi$ and are effectively pinned at the center of the physical domain. Furthermore, equation (40) is well posed as long as the amplitude of the first Fourier mode of $u(x,t;\omega)$ remains finite for all realizations and at all times. This is always true unless one of the realizations has spatial periodicity equal to half that of the domain, which is unlikely to happen in practice.

3.4 Initialization of the SDO quantities

Given an initial ensemble of realizations $u_{0}(x;\omega)$, the initial conditions for the quantities involved in the SDO computation are found by a two-step process, where (i) the symmetry-reduced version $\hat{u}_{0}(x;\omega)$ of each realization is calculated, leading to the initial stochastic translation amount $c_{0}(\omega)$, and (ii) the initial symmetry-reduced mean $\bar{u}_{0}(x)$, modes $\hat{u}_{i0}(x)$ and stochastic coefficients $Y_{i0}(\omega)$ are obtained from a truncated Karhunen-Loeve expansion of $\hat{u}_{0}(x;\omega)$.

Step 1. First, the symmetry-reduced counterparts $\hat{u}_{0}(x;\omega)=g^{-1}(c_{0}(\omega))u_{0}(x;\omega)$ of the initial realizations $u_{0}(x;\omega)$ are calculated from

[TABLE]

where $\tilde{u}_{0}(k;\omega)$ are the Fourier coefficients associated to each realization $u_{0}(x;\omega)$. To bring the symmetry-reduced state $\hat{u}_{0}(x;\omega)$ to the first Fourier mode slice, the initial stochastic translation amount $c_{0}(\omega)$ is chosen such that the phase angle of the first Fourier mode of $\hat{u}_{0}(x;\omega)$ is equal to $\pi$ for all realizations $\omega$, that is $c_{0}(\omega)=-\arg\tilde{u}_{0}(1;\omega)L/2\pi+L/2$.

Step 2. The initial symmetry-reduced realizations $\hat{u}_{0}(x;\omega)$ can then be approximated through a truncated Karhunen-Loeve expansion, giving a set of symmetry-reduced modes and stochastic coefficients from which the SDO modes and coefficients can be initialized. Defining first the initial symmetry-reduced mean $\bar{u}_{0}(x)=E[\hat{u}_{0}(x;\omega)]$, the initial symmetry-reduced orthonormal modes $\hat{u}_{i0}(x)$ are then given by the $s$ most energetic eigenfunctions of the following eigenvalue problem

[TABLE]

where the correlation operator $R(x,y)=E[(\hat{u}_{0}(x;\omega)-\bar{u}_{0}(x))(\hat{u}_{0}(y;\omega)-\bar{u}_{0}(y))]$. Finally, the associated initial stochastic coefficients $Y_{i0}(\omega)$ are obtained by projection of the initial symmetry-reduced realizations to the symmetry-reduced modes

[TABLE]

Note that $s$ not only represents the number of modes and stochastic coefficients of the SDO scheme, but also sets the tolerance on the truncated Karhunen-Loeve decomposition (46) of the initial condition. Therefore, its value should be chosen so that both (i) the stochastic symmetry-reduced initial condition $\hat{u}_{0}(x;\omega)$ and (ii) the subsequent time-evolving stochastic solution $\hat{u}(x,t;\omega)$ are sufficiently well approximated. The accuracy of the initial condition can be monitored by the decay of the eigenvalues $\lambda_{i}$ of the truncated Karhunen-Loeve decomposition (46), which indicate the amount of variance in $\hat{u}_{0}(x;\omega)$ contained along each of the corresponding eigendirections $\hat{u}_{i0}(x)$. Similarly, one can get a rough idea on the accuracy of the stochastic solution by looking at the evolution of the variance in the stochastic coefficients associated with the last few modes. A rigorous approach, however, requires knowledge of the exact solution (either from Monte-Carlo simulations or analytical arguments).

In section 3.6, we will compare numerical results from the SDO framework with those from a DO computation. In the case of DO, the various quantities are initialized following Step 2 directly, i.e. the initial mean $\bar{u}_{0}(x)$, modes $u_{i0}(x)$ and stochastic coefficients $Y_{i0}(\omega)$ are defined from a truncated Karhunen-Loeve expansion of the initial full state-space realizations $u_{0}(x;\omega)$.

3.5 Numerical scheme

The SDO equations for the mean and the modes are implemented using a standard pseudo-spectral method in space with 3/2 antialiasing and a semi-implicit Euler scheme in time, where the third-order derivatives are treated implicitly and the nonlinear terms are treated explicitly. The stochastic coefficients and translation amount are integrated in time with respectively a 4th-order Runge-Kutta and a 2-step Adams-Bashforth scheme, both using a particle method. Note that setting the stochastic translation amount to zero in the SDO equations readily yields the standard DO scheme. For both SDO and DO, we use $L=2\pi$, $\mu=5\cdot 10^{-4}$, 512 Fourier modes, $\Delta t=10^{-4}$ and 1000 Monte-Carlo particles.

3.6 Example with a stochastic soliton

The KdV equation admits a class of solitary wave solutions with shape-dependent propagation speed. In this section, we use the SDO equations to compute the evolution of a stochastic initial condition consisting of such solitons

[TABLE]

where $a(\omega)\sim\mathcal{U}(0.1,0.5)$ is a random variable that affects both the amplitude and the width of the hyperbolic secant profile. Each realization $\omega$ defined by the above initial condition is an exact soliton solution of the KdV equation with propagation velocity $a(\omega)$, which makes this problem challenging for classical order reduction methods as a large number of modes eventually becomes necessary to reproduce faithfully the spatial dispersion of all realisations. The SDO scheme, on the other hand, should not suffer from such issues since the shift amount $c(t;\omega)$ takes care of the spatial translation of the realizations, leaving the symmetry-reduced mean and modes to account for a change in shape that is here nonexistent since the solution consists of solitons. Hence, we use just one mode for the SDO computation, which we initialize according to the procedure described in Section 3.4. For comparison purposes, we also perform a regular DO simulation of the same problem, but using 10 modes instead.

The first two rows in Figure 4 show the symmetry-reduced mean $\bar{u}(x,t)$ and the single symmetry-reduced mode $\hat{u}_{1}(x,t)$ of the SDO solution at final time $t=3$ (left), and their spatio-temporal evolution (right).

As expected, both the symmetry-reduced mean and mode are steady, which reflects the fact that the dynamics of the individual solitons purely consist of translation, which is entirely absorbed in the stochastic shift amount $c(t;\omega)$. The original stochastic state $u(x,t;\omega)$ can be reconstructed from the SDO quantities through the group transformation (41), and we display in the last row of Figure 4 ten full state-space realizations at final time $t=3$ (left) and the spatio-temporal evolution of one of them (right). Each individual soliton realization maintains its initial shape while propagating at a constant amplitude-dependent velocity, which is expected but nonetheless remarkable considering that the SDO computation uses just one single mode.

To put these results into context, Figure 5 shows the mean $\bar{u}(x,t)$, the first two modes $u_{i}(x,t)$ and ten full state-space realizations $u(x,t;\omega)$ of the corresponding DO simulation.

Contrary to the SDO simulation, the DO realizations at final time have not preserved their initial shape although the DO simulation is using ten times as many modes. This is because the DO mean and modes need to account for the translation at random speed of the stochastic soliton, hence they become dispersed over space as time progresses.

This difference in the amount of stochasticity accounted for by the modes is clearly observed in Figure 6, where we plot the time evolution of the energy in the mean $\langle\bar{u},\bar{u}\rangle$ and the variance of the stochastic coefficients $E[Y_{i}^{2}]$ for the SDO (left) and DO (right) solutions.

The stochastic energy of the SDO mode, represented by the variance of the corresponding stochastic coefficient, remains equal to its initial value while that of the DO modes grows drastically over time. Thus, the DO scheme would need a much larger number of modes to achieve a comparable level of accuracy to the SDO simulation, which proves the better efficiency of SDO.

4 Application to the Navier-Stokes equations

In this section, we illustrate our approach with the two-dimensional Navier-Stokes equations

[TABLE]

written in nondimensional variables, with $Re$ the Reynolds number, and defined on the spatial domain $\mathbf{x}\in D=[0,L_{1}]\times[0,L_{2}]$ with periodic boundary conditions. The Navier-Stokes equations are equivariant under translations, rotations and inversion about the origin. Focusing on the group of continuous translations along the two spatial directions $G=\mathrm{SO}(2)_{x_{1}}\times\mathrm{SO}(2)_{x_{2}}$, we define the shift operator $g\in G$ as

[TABLE]

where the continuous phase parameter $\mathbf{c}=(c_{1},c_{2})^{\mathsf{T}}$ represents the shift amounts in the $x_{1}$ and $x_{2}$ directions. The group orbit $g\mathbf{u}$ at an arbitrary state $\mathbf{u}$ thus possesses two tangents, one for each direction $x_{1}$ and $x_{2}$, given by

[TABLE]

4.1 Dynamics within the symmetry-reduced state space

As before, we now pick a fixed template $\hat{\mathbf{u}}^{\prime}(\mathbf{x})$ and we consider the evolution of the symmetry-reduced state $\hat{\mathbf{u}}(\mathbf{x},t;\omega)=g^{-1}(\mathbf{c}(t;\omega))\mathbf{u}(\mathbf{x},t;\omega)$ such that the distance $||\hat{\mathbf{u}}(\mathbf{x},t;\omega)-\hat{\mathbf{u}}^{\prime}(\mathbf{x})||$ between the reduced state and the template is minimized at all times $t$ and for all realizations $\omega$. In the case of the Navier-Stokes equations subject to the symmetry group $G$, the dynamical equations (15) for $\hat{\mathbf{u}}$ in the symmetry-reduced state space become

[TABLE]

where the full state space nonlinear differential operator $\mathbf{F}$ has the form

[TABLE]

with $p$ acting as a Lagrange multiplier to enforce $\nabla\cdot\hat{\mathbf{u}}=0$ since the divergence-free condition on $\mathbf{u}$ implies that $\hat{\mathbf{u}}$ is also divergence-free. Finally, the stochastic differential equation (17) for the shift amount $\mathbf{c}(t;\omega)$ writes

[TABLE]

4.2 Order reduction in the symmetry-reduced state space

In order to formulate the reduced-order SDO equations, let us first approximate the symmetry-reduced stochastic state as a truncated Karhunen-Loeve expansion

[TABLE]

where $\bar{\mathbf{u}}$ and $\hat{\mathbf{u}}_{i}$ are respectively the symmetry-reduced mean and modes and $Y_{i}$ are the stochastic coefficients of the finite-dimensional expansion. Next, we insert the above representation into the RHS of the evolution equation (52) for the symmetry-reduced state to obtain

[TABLE]

where $\mathbf{F}_{0}$, $\mathbf{F}_{i}$ and $\mathbf{F}_{ij}$ are deterministic fields given by

[TABLE]

The stochastic pressure $p$ ensures that $\hat{\mathbf{u}}$ is divergence-free, hence by setting the divergence of the expanded RHS operator (56) to zero (following [31]), one finds that

[TABLE]

where $\Delta p_{0}=\nabla\cdot\mathbf{F}_{0}$, $\Delta p_{i}=\nabla\cdot\mathbf{F}_{i}$ and $\Delta p_{ij}=\nabla\cdot\mathbf{F}_{ij}$ are deterministic pressure fields. Now, we insert the expanded RHS operator (56) governing the evolution of the symmetry-reduced state $\hat{\mathbf{u}}$ into the SDO equations (20), (21) and (22) to find the following evolution equation for the mean

[TABLE]

We then have the following evolution equation for the stochastic coefficients

[TABLE]

where we have defined $\langle\nabla\bar{\mathbf{u}},\hat{\mathbf{u}}_{i}\rangle=(\langle\partial_{x_{1}}\bar{\mathbf{u}},\hat{\mathbf{u}}_{i}\rangle,\langle\partial_{x_{2}}\bar{\mathbf{u}},\hat{\mathbf{u}}_{i}\rangle)^{\mathsf{T}}$, and similarly $\langle\nabla\hat{\mathbf{u}}_{m},\hat{\mathbf{u}}_{i}\rangle=(\langle\partial_{x_{1}}\hat{\mathbf{u}}_{m},\hat{\mathbf{u}}_{i}\rangle,\langle\partial_{x_{2}}\hat{\mathbf{u}}_{m},\hat{\mathbf{u}}_{i}\rangle)^{\mathsf{T}}$. The modes are governed by the following equation

[TABLE]

where the deterministic fields $\hat{\mathbf{H}}_{i}$ are defined as

[TABLE]

with $M_{kmn}=E[Y_{k}Y_{m}Y_{n}]$ the matrix of third-order moments. Finally, the time evolution of the shift amount $\dot{\mathbf{c}}$ is given by the stochastic linear system (54) where the stochastic matrix $\mathbf{T}$ and vector $\mathbf{f}$ can be expanded in functions of $Y_{i}$ and $Y_{i}Y_{j}$ with $i,j=1,...,s$ using that $\mathbf{t}_{a}(\hat{\mathbf{u}})=\mathbf{t}_{a}(\bar{\mathbf{u}})+Y_{i}\mathbf{t}_{a}(\hat{\mathbf{u}}_{i})$, $a=1,2$ and $\mathbf{F}(\hat{\mathbf{u}})=-\nabla p_{0}+\mathbf{F}_{0}+Y_{i}\,(-\nabla p_{i}+\mathbf{F}_{i})+Y_{i}Y_{j}\,(-\nabla p_{ij}+\mathbf{F}_{ij})$.

At any given time, the full stochastic state can be reconstructed from the reduced-order quantities and the time-integrated stochastic shift amount through the group transformation

[TABLE]

4.3 Choice of the template

The first Fourier mode slice is unambiguously defined for scalar fields since there is a unique Fourier coefficient with wavenumber one which defines the symmetry-reduced state. However, ambiguity arises when one considers vector fields with multiples components and therefore multiple Fourier coefficients of wavenumber one. Nevertheless, we show here that for the particular case of fluid flow in two dimensions, one can overcome this issue by defining the first Fourier mode slice based on the scalar vorticity field instead.

Let us introduce an arbitrary deterministic two-dimensional velocity field $\mathbf{u}(\mathbf{x})=(u_{1}(x_{1},x_{2}),u_{2}(x_{2},x_{2}))^{\mathsf{T}}$ and its symmetry-reduced counterpart $\hat{\mathbf{u}}=g^{-1}(\mathbf{c})\mathbf{u}$, where $\mathbf{c}=(c_{1},c_{2})^{\mathsf{T}}$ is such that both slice conditions $\langle\hat{\mathbf{u}},\mathbf{t}_{1}(\hat{\mathbf{u}}^{\prime})\rangle=0$ and $\langle\hat{\mathbf{u}},\mathbf{t}_{2}(\hat{\mathbf{u}}^{\prime})\rangle=0$ are satisfied given a fixed template function $\hat{\mathbf{u}}^{\prime}$. Consider now the scalar vorticity field defined by $\omega=\partial_{x_{1}}u_{2}-\partial_{x_{2}}u_{1}$. Its Fourier series is

[TABLE]

where $\tilde{\omega}=ik_{1}\tilde{u}_{2}-ik_{2}\tilde{u}_{1}$, with $\tilde{u}_{1}$ and $\tilde{u}_{2}$ the Fourier coefficients of $u_{1}$ and $u_{2}$ respectively. The corresponding symmetry-reduced vorticity field is expressed as

[TABLE]

where $|\tilde{\omega}(k_{1},k_{2})|$ and $\arg\tilde{\omega}(k_{1},k_{2})$ are respectively the modulus and phase angle of $\tilde{\omega}$. Similar to the one-dimensional case, the above relation shows that the Fourier coefficients of $\hat{\omega}$ are those of $\omega$ rotated by an amount equal to $2\pi(k_{1}c_{1}/L_{1}+k_{2}c_{2}/L_{2})$. Therefore we can define the first Fourier mode slice by choosing $c_{1}$ and $c_{2}$ such that the phase angles of the first Fourier modes of $\hat{\omega}$ in the $x_{1}$ and $x_{2}$ directions are equal to $\pi$, i.e. $\arg\tilde{\omega}(1,0)+2\pi c_{1}/L_{1}=\pi$ and $\arg\tilde{\omega}(0,1)+2\pi c_{2}/L_{2}=\pi$. In this way, the symmetry-reduced vorticity field will appear centered in the domain. We can now use $\arg\tilde{\omega}(1,0)=\arg\tilde{u}_{2}(1,0)+\pi/2$ and $\arg\tilde{\omega}(0,1)=\arg\tilde{u}_{1}(0,1)-\pi/2$ to express this choice of $c_{1}$ and $c_{2}$ in terms of the Fourier coefficients of the velocity field, resulting in $c_{1}=-\arg\tilde{u}_{2}(1,0)L_{1}/2\pi+L_{1}/4$ and $c_{2}=-\arg\tilde{u}_{1}(0,1)L_{2}/2\pi-L_{2}/4$. Finally, we show in Appendix C that this choice of $\mathbf{c}$ corresponds to the template function

[TABLE]

The first Fourier mode slice defined above can be readily applied to our stochastic SDO equations by choosing (66) as our template function. In this way, all symmetry-reduced realizations $\hat{\mathbf{u}}(\mathbf{x},t;\omega)$ have the phase of their first Fourier modes in the $x$ and $y$ directions fixed to $\pi$ and therefore remain centered in the two-dimensional physical space. Furthermore, it is possible to show that using template (66), the matrix $\mathbf{T}$ appearing in equation (54) is diagonal, with the two diagonal elements given respectively by the amplitudes of the first Fourier modes of $u_{2}$ in the $x$ direction and $u_{1}$ in the $y$ direction. As such, $\mathbf{T}$ is well conditioned so long as the ratio of these two quantities does not get large, and equation (54) is well posed as long as these quantities remain nonzero for all realizations and at all times. This condition is always true unless one of the realizations has spatial periodicity equal to half that of the domain, which is unlikely to happen in practice.

4.4 Initialization of the SDO quantities

Given an initial ensemble of realizations $\mathbf{u}_{0}(\mathbf{x};\omega)$, the initial conditions for the quantities involved in the SDO computation are found by a two-step process similar to Section 3.4, where (i) the symmetry-reduced version $\hat{\mathbf{u}}_{0}(\mathbf{x};\omega)$ of each realization is calculated, leading to the initial stochastic translation amount $\mathbf{c}_{0}(\omega)$, and (ii) the initial symmetry-reduced mean $\bar{\mathbf{u}}_{0}(\mathbf{x})$, modes $\hat{\mathbf{u}}_{i0}(\mathbf{x})$ and stochastic coefficients $Y_{i0}(\omega)$ are obtained from a truncated Karhunen-Loeve expansion of $\hat{\mathbf{u}}_{0}(\mathbf{x};\omega)$.

Step 1. First, the symmetry-reduced counterparts $\hat{\mathbf{u}}_{0}(\mathbf{x};\omega)=g^{-1}(\mathbf{c}_{0}(\omega))u_{0}(\mathbf{x};\omega)$ of the initial realizations $\mathbf{u}_{0}(\mathbf{x};\omega)$ are given by

[TABLE]

where $\tilde{\mathbf{u}}_{0}(k_{1},k_{2};\omega)=(\tilde{u}_{10}(k_{1},k_{2};\omega),\tilde{u}_{20}(k_{1},k_{2};\omega))^{\mathsf{T}}$ are the Fourier coefficients associated to each realization $\mathbf{u}_{0}(\mathbf{x};\omega)=(u_{10}(\mathbf{x};\omega),u_{20}(\mathbf{x};\omega))^{\mathsf{T}}$. To bring the symmetry-reduced state $\hat{\mathbf{u}}_{0}(\mathbf{x};\omega)$ to the first Fourier mode slice, we proceed according to Section 4.3 and choose the initial stochastic translation amount $\mathbf{c}_{0}(\omega)=(c_{10}(\omega),c_{20}(\omega))^{\mathsf{T}}$ as $c_{10}(\omega)=-\arg\tilde{u}_{20}(1,0;\omega)L_{1}/2\pi+L_{1}/4$ and $c_{20}(\omega)=-\arg\tilde{u}_{10}(0,1;\omega)L_{2}/2\pi-L_{2}/4$.

Step 2. The initial symmetry-reduced realizations $\hat{\mathbf{u}}_{0}(\mathbf{x};\omega)$ can then be approximated through a truncated Karhunen-Loeve expansion, giving a set of symmetry-reduced modes and stochastic coefficients from which the SDO modes and coefficients can be initialized. Defining first the initial symmetry-reduced mean $\bar{\mathbf{u}}_{0}(\mathbf{x})=E[\hat{\mathbf{u}}_{0}(\mathbf{x};\omega)]$, the initial symmetry-reduced orthonormal modes $\hat{\mathbf{u}}_{i0}(\mathbf{x})$ are then given by the $s$ most energetic eigenfunctions of the following eigenvalue problem

[TABLE]

where the correlation matrix is $\mathbf{R}(\mathbf{x},\mathbf{y})=E[(\hat{\mathbf{u}}_{0}(\mathbf{x};\omega)-\bar{\mathbf{u}}_{0}(\mathbf{x}))(\hat{\mathbf{u}}_{0}(\mathbf{y};\omega)-\bar{\mathbf{u}}_{0}(\mathbf{y}))^{\mathsf{T}}]$. Finally, the associated initial stochastic coefficients $Y_{i0}(\omega)$ are obtained by projection of the initial symmetry-reduced realizations to the symmetry-reduced modes

[TABLE]

In section 4.6, we will compare numerical results from the SDO framework with results from the regular DO equations. In the case of DO, the various quantities are initialized following Step 2 directly, i.e. the initial mean $\bar{\mathbf{u}}_{0}(\mathbf{x})$, modes $\mathbf{u}_{i0}(\mathbf{x})$ and stochastic coefficients $Y_{i0}(\omega)$ are defined from a truncated Karhunen-Loeve expansion of the initial full state space realizations $\mathbf{u}_{0}(\mathbf{x};\omega)$.

4.5 Numerical scheme

The SDO equations for the mean and the modes are implemented using a standard pseudo-spectral method in space with 3/2 dealiasing and explicit Euler finite differences in time. The stochastic coefficients and translation amount are integrated in time with respectively a 4th-order Runge-Kutta and a 2-step Adams-Bashforth scheme, both using a particle method. Note that setting the stochastic translation amount to zero in the SDO equations readily yields the standard DO scheme. For both SDO and DO, we use $L_{1}=L_{2}=2\pi$, $64\times 64$ Fourier modes, $\Delta t=0.001$ and 1000 Monte-Carlo particles.

4.6 Example on a stochastically advected vortex

We now use the SDO equations to compute the transient flow of a randomly advected stochastic Lamb-Oseen vortex, described by the following initial velocity field

[TABLE]

where $r=\sqrt{(x_{1}-L_{1}/2)^{2}+(x_{2}-L_{2}/2)^{2}}$ and $\theta=\arctan x_{2}/x_{1}$. In the above initial condition, a Lamb-Oseen vortex centered in the domain with random radial length scale $r_{c}(\omega)\sim\mathcal{N}(0.2,10^{-4})$ is superimposed to a uniform flow of unit magnitude and random direction $\psi(\omega)\sim\mathcal{U}(0,\pi/2)$. We take $\Gamma=10$ and set the Reynolds number to $Re=40$. Each realization $\omega$ defined by the above initial condition therefore undergoes viscous diffusion while being simultaneously advected in a given direction due to the uniform flow. Thus, after a finite time, the realizations will have moved away from their initial position to various locations in the domain.

The SDO simulation is initialized according to the procedure described in Section 4.4. For comparison purposes, we also perform a regular DO simulation of the same problem and we use 6 modes in both computations. Since (70) defines the initial realizations as centered in the domain, they are initially identical to their symmetry-reduced counterparts which are also required to be centered in the domain. It therefore follows that the initial conditions for the SDO and DO mean and modes are the same, and they are displayed in Figure 7 in terms of their vorticity (background color) and velocity (arrows) fields.

Notice that the first two modes mainly represent the variability due to the uniform advection flow of unit magnitude and random direction, whereas the subsequent modes account for the variability due to the random radial length scale of the Lamb-Oseen vortex.

Figure 8 shows the symmetry-reduced mean $\bar{\mathbf{u}}(\mathbf{x},t)$ and symmetry-reduced modes $\hat{\mathbf{u}}_{i}(\mathbf{x},t)$ of the SDO solution at time $t=2.5$, together with one full state space realization reconstructed from the SDO quantities through the group transformation (63).

Compared with the initial condition shown in Figure 7, we observe that the symmetry-reduced mean and modes have remained at their initial location and have merely diffused out due to viscosity. On the other hand, the realization has not only diffused out but has moved away from the center of the domain along the path shown by the white line, thanks to the stochastic shift amount $\mathbf{c}(t;\omega)$ which has tracked the advection due to the stochastic uniform flow.

In the left-hand side of Figure 9, we further display the individual trajectories of 100 realizations, obtained from the time history of the shift amount.

The realizations move away from their initial position at unit velocity and in a random direction uniformly distributed between 0 and $\pi/2$, eventually landing at final time on a circle of radius 2.5. These trajectories are expected given the initial condition (70), but nonetheless remarkable since the SDO equations do not have a priori knowledge of the paths followed by the realizations. The exact radial distance $||\mathbf{c}(t;\omega)||$ of the realizations to their initial position at final time is shown in the right-hand side of Figure 9 and has negligible standard deviation on the order of $10^{-11}$, which demonstrates the numerical accuracy of the stochastic shift amount.

For comparison purposes, Figure 10 shows the mean $\bar{\mathbf{u}}(\mathbf{x},t)$ and modes $\mathbf{u}_{i}(\mathbf{x},t)$ of the DO solution at time $t=2.5$, together with one realization.

Since the standard DO scheme directly performs order reduction in the full state space, the DO mean and modes need to account for the spatial dispersion of the realizations. As a consequence, the modes in Figure 10 do not directly reflect the shape of each individual vortex and a limited number of them is unable to accurately reproduce each realization. On the contrary, the SDO modes are physically more relevant since they directly reveal the underlying non-trivial dynamics of the solution.

Figure 11 shows the statistics of the first four SDO (left) and DO (right) stochastic coefficients at initial and final times, displayed in terms of two-dimensional marginals.

The statistical structure of the SDO solution remains identical over time except for some decline in the variance of $Y_{3}$ and $Y_{4}$ due to viscous diffusion. On the other hand, the statistical structure of the DO solution completely changes over time since the modes have to account for the spatial translation of each realization, thereby obscuring the fact that each individual vortex is merely diffusing besides the advection.

This difference in the amount of stochastic variability accounted for by the modes is clearly apparent in Figure 12, which displays the time evolution of the energy in the mean $\langle\bar{\mathbf{u}},\bar{\mathbf{u}}\rangle$ and the variance of the stochastic coefficients $E[Y_{i}^{2}]$ for the SDO (left) and DO (right) solutions.

The variance of the SDO mean and modes decreases over time since the symmetry-reduced state simply undergoes viscous diffusion. On the other hand, that of their DO counterparts grows over time as the realizations become more and more dispersed. Therefore, the SDO scheme achieves much better accuracy than DO with the same number of modes, as can also be seen from the individual and cumulative distributions of energy in the mean and the modes at final time $t=2.5$ shown in Figure 13.

These simulations are shown as videos in the Supplementary Materials, together with a second case where the direction of the uniform advection flow in the initial condition is uniformly distributed between 0 and $2\pi$, instead of between 0 and $\pi/2$.

5 Conclusions

In this work, we have introduced a novel methodology for efficient order reduction of stochastic dynamical systems with continuous symmetries. This methodology is composed of two steps. In a first step, one performs symmetry reduction of the original dynamical system using the method of slices. In this way, the dynamics of the original system is decoupled into shape deformations, captured by a symmetry-reduced stochastic state fixed in the physical domain, and motion along the symmetry directions of the system, tracked by a set of scalar phase parameters. The second step consists in order reduction of the symmetry-reduced stochastic state, using any standard order reduction method of choice. Since the symmetry-reduced state is fixed in the physical domain, this procedure results in very efficient mixed symmetry-dimensionality reduction schemes. In particular, using the Dynamically Orthogonal (DO) equations to perform the second step, we have obtained a new Symmetry-reduced Dynamically Orthogonal (SDO) scheme that shows much better performance than DO on stochastic solutions of the 1D Korteweg-de Vries and 2D Navier-Stokes equations.

Even though both examples we have presented consider translation invariance and utilize periodic boundary conditions, it is important to recognize that continuous symmetries can arise in bounded domains as well. A prominent example is axisymmetric systems, for instance the flow past an axisymmetric body or in a Taylor-Couette cell. Using axisymmetric coordinates, the rotational symmetry translates to translational invariance in the azimuthal direction, allowing for the direct application of our technique. Furthermore, while we have only considered stochasticity through the initial conditions, we expect our framework to be equally applicable to stochastic operators provided that the stochastic terms (either coefficients or forcing) are statistically homogeneous in the symmetry direction, ensuring that the equivariance condition (3) is satisfied in a statistical sense.

In general, we expect this methodology to work well whenever the stochastic solution is composed of coherent structures or sharp gradients like in shock waves. Possible future work could include the extension of this framework to dynamical systems with self-similar solutions, which can also be symmetry-reduced with the method of slices by rescaling time as well [2, 29, 32]. Finally, one could also envision complementing our method with data assimilation techniques such as Bayesian inference, to which group transformation ideas have already been applied [26].

Acknowledgments

The authors would like to thank Dr. Edouard Boujo for a careful reading of the manuscript, the two anonymous reviewers for their constructive comments, as well as Dr. Sai Ravela, Dr. Mohammad Farazmand and Florian Feppon for helpful discussions related to the subject. This work has been supported through the ONR grant N00014-15-1-2381, the AFOSR grants FA9550-16-1-0231 and W911NF-17-1-0306, and the ARO grant 66710-EG-YIP.

Appendix A DO order reduction preserves symmetry reduction

In this appendix, we show that after DO order reduction of the symmetry-reduced equations, the dynamics of the reduced-order solution remain confined to the symmetry-reduced state space. First, assume that the phase parameters $\boldsymbol{\phi}$ satisfy the stochastic equation (24), which implies

[TABLE]

for all symmetry directions $b=1,...,N$ and all reduced-order realizations (18). We also assume that at the current time instant, the mean and the modes satisfy the slice condition (11), that is $\langle\bar{\mathbf{u}},\mathbf{t}_{a}^{\prime}\rangle=\langle\hat{\mathbf{u}}_{i},\mathbf{t}_{a}^{\prime}\rangle=0$. Then, taking the inner products of the evolution equation (20) for the mean with the template tangents, we find that

[TABLE]

for $b=1,...,N$. Likewise, taking the inner products of the evolution equation (22) for the modes with the template tangents, we have

[TABLE]

where we have used $\langle\hat{\mathbf{u}}_{j},\mathbf{t}_{b}^{\prime}\rangle=0$ and the fact that

[TABLE]

As in (16), equations (72) and (73) show that the mean and modes will satisfy the slice condition (11) throughout time integration, so that the reduced-order solution (18) remains confined to the symmetry-reduced state space. Numerical round-off errors remain negligible in our simulations.

Appendix B First Fourier mode slice for the KdV equation

Here, we justify that the template function (44) corresponds to the first Fourier mode slice defined in Section 3.3, that is the slice condition $\langle\hat{u},t(\hat{u}^{\prime})\rangle=0$ is satisfied when the phase angle of the first Fourier mode of $\hat{u}$ is equal to $\pi$. First, note that the inner product between two generic states $u(x)$ and $v(x)$ can be expressed in terms of their Fourier coefficients $\tilde{u}(k)=a(k)+ib(k)$ and $\tilde{v}(k)=c(k)+id(k)$ as

[TABLE]

where we have used that $\tilde{u}(-k)=\tilde{u}^{*}(k)$ and $\tilde{v}(-k)=\tilde{v}^{*}(k)$ since $u$ and $v$ are real. If we now denote $\tilde{\hat{u}}(k)=a(k)+ib(k)$ and $\tilde{\hat{u}}^{\prime}(k)=c(k)+id(k)$ the Fourier coefficients of the symmetry-reduced state $\hat{u}$ and template $\hat{u}^{\prime}$, respectively, we can use (75) to express the slice condition as

[TABLE]

The first Fourier mode slice is defined by $\arg\tilde{\hat{u}}(1)=\pi$, thus $b(1)=0$, $a(1)<0$ and all other coefficients $a(k)$, $b(k)$ are possibly non-zero. Except for the sign of $a(1)$, these conditions are equivalent to the above slice condition when $c(1)\neq 0$ and all other coefficients $c(k)$, $d(k)$ are zero. Hence a valid choice of template is

[TABLE]

which corresponds to $c(1)=1/2$. With this template, the slice condition (76) reduces to $\langle\hat{u},t(\hat{u}^{\prime})\rangle=-2\pi b(1)=0$, which is equivalent to requiring that $\arg\tilde{\hat{u}}(1)=0$ or $\pi$. This ambiguity between zero and $\pi$ is not a problem in the actual SDO computations since the symmetry-reduced state $\hat{u}$ is initialized with $\arg\tilde{\hat{u}}(1,t_{0};\omega)=\pi$, and continuity of the phase angle during time integration ensures that $\arg\tilde{\hat{u}}(1,t;\omega)$ does not jump to zero at later times. With the phase of its first Fourier mode fixed, the symmetry-reduced stochastic state $\hat{u}(x,t;\omega)$ is effectively pinned at a fixed location in space and this is accomplished with a generic template that is not problem-dependent.

Appendix C First Fourier mode slice for the Navier-Stokes equation

Here, we justify that template (66) corresponds to the first Fourier mode slice defined in Section 4.3 for a two-dimensional velocity field, that is, both slice conditions $\langle\hat{\mathbf{u}},\mathbf{t}_{1}(\hat{\mathbf{u}}^{\prime})\rangle=0$ and $\langle\hat{\mathbf{u}},\mathbf{t}_{2}(\hat{\mathbf{u}}^{\prime})\rangle=0$ are satisfied when the phase angles of the first Fourier modes of $\hat{\omega}$ in the $x$ and $y$ directions are equal to $\pi$. First, note that the inner product between two generic states $\mathbf{u}(\mathbf{x})=(u_{1}(x_{1},x_{2}),u_{2}(x_{1},x_{2}))^{\mathsf{T}}$ and $\mathbf{v}(\mathbf{x})=(v_{1}(x_{1},x_{2}),v_{2}(x_{1},x_{2}))^{\mathsf{T}}$ can be expressed in terms of their Fourier coefficients $\tilde{\mathbf{u}}(\mathbf{k})=(\tilde{u}_{1}(k_{1},k_{2}),\tilde{u}_{2}(k_{1},k_{2}))^{\mathsf{T}}$ and $\tilde{\mathbf{v}}(\mathbf{k})=(\tilde{v}_{1}(k_{1},k_{2}),\tilde{v}_{2}(k_{1},k_{2}))^{\mathsf{T}}$ as

[TABLE]

where $\tilde{u}_{j}(k_{1},k_{2})=a_{j}(k_{1},k_{2})+ib_{j}(k_{1},k_{2})$ and $\tilde{v}_{j}(k_{1},k_{2})=c_{j}(k_{1},k_{2})+id_{j}(k_{1},k_{2})$, $j=1,2$, and we have used the realness of $u_{j}$ and $v_{j}$. If we now denote $\tilde{\hat{\mathbf{u}}}(\mathbf{k})=(\tilde{\hat{u}}_{1}(k_{1},k_{2}),\tilde{\hat{u}}_{2}(k_{1},k_{2}))^{\mathsf{T}}$ and $\tilde{\hat{\mathbf{u}}}^{\prime}(\mathbf{k})=(\tilde{\hat{u}}^{\prime}_{1}(k_{1},k_{2}),\tilde{\hat{u}}^{\prime}_{2}(k_{1},k_{2}))^{\mathsf{T}}$ the Fourier coefficients of the symmetry-reduced state $\hat{\mathbf{u}}$ and template $\hat{\mathbf{u}}^{\prime}$, respectively, we can use (78) to express the slice conditions $\langle\hat{\mathbf{u}},\mathbf{t}_{a}(\hat{\mathbf{u}}^{\prime})\rangle=0$, $a=1,2$ as

[TABLE]

where $\tilde{\hat{u}}_{j}(k_{1},k_{2})=a_{j}(k_{1},k_{2})+ib_{j}(k_{1},k_{2})$ and $\tilde{\hat{u}}^{\prime}_{j}(k_{1},k_{2})=c_{j}(k_{1},k_{2})+id_{j}(k_{1},k_{2})$, $j=1,2$. From Section 4.3, one can show that the first Fourier mode slice corresponds to $\arg\tilde{\hat{u}}_{1}(0,1)=-\pi/2$ and $\arg\tilde{\hat{u}}_{2}(1,0)=\pi/2$, thus $a_{1}(0,1)=a_{2}(1,0)=0$, $b_{1}(0,1)<0$, $b_{2}(1,0)>0$ and all other coefficients $a_{j}(k_{1},k_{2})$, $b_{j}(k_{1},k_{2})$, $j=1,2$ are possibly non-zero. Except for the signs of $b_{1}(0,1)$ and $b_{2}(1,0)$, these conditions are equivalent to the slice conditions (79) when $d_{1}(0,1)\neq 0$, $d_{2}(1,0)\neq 0$ and all other coefficients $c_{j}(k_{1},k_{2})$, $d_{j}(k_{1},k_{2})$, $j=1,2$ are zero. Hence a valid choice of template is

[TABLE]

which corresponds to $d_{1}(0,1)=d_{2}(1,0)=-1/2$. With this template, the slice conditions (79) become $\langle\hat{\mathbf{u}},\mathbf{t}_{1}(\hat{\mathbf{u}}^{\prime})\rangle=-2\pi L_{2}a_{2}(1,0)=0$ and $\langle\hat{\mathbf{u}},\mathbf{t}_{2}(\hat{\mathbf{u}}^{\prime})\rangle=-2\pi L_{1}a_{1}(0,1)=0$, which means requiring that $\arg\tilde{\hat{u}}_{2}(1,0)=\arg\tilde{\hat{u}}_{1}(0,1)=\pm\pi/2$. Same as for the one-dimensional case, this ambiguity between $\pm\pi/2$ is not an issue in the actual SDO computations since the symmetry-reduced state $\hat{\mathbf{u}}$ is initialized according to the first Fourier mode slice with $\arg\tilde{\hat{u}}_{1}(0,1,t_{0};\omega)=-\pi/2$ and $\arg\tilde{\hat{u}}_{2}(1,0,t_{0};\omega)=\pi/2$. Continuity of the phase angles during time integration then ensures that $\arg\tilde{\hat{u}}_{1}(0,1,t;\omega)$ and $\arg\tilde{\hat{u}}_{2}(1,0,t;\omega)$ do not jump to the other root at later times.