Supervised parallel-in-time algorithm for long-time Lagrangian simulations of stochastic dynamics: Application to hydrodynamics

TL;DR

This paper introduces a supervised parallel-in-time algorithm that accelerates long-time stochastic Lagrangian simulations by leveraging a macroscopic Navier-Stokes model as a predictor, effectively bridging microscopic and macroscopic dynamics.

Contribution

The novel SPASD method uses a low-fidelity macroscopic model to supervise and correct high-dimensional microscopic simulations, enabling efficient long-time stochastic dynamics simulations.

Findings

Accelerates long-time Lagrangian simulations of stochastic dynamics.

Uses a macroscopic Navier-Stokes model as a predictor for microscopic simulations.

Effectively recovers microscopic details and fluctuations through correction.

Abstract

Lagrangian particle methods based on detailed atomic and molecular models are powerful computational tools for studying the dynamics of microscale and nanoscale systems. However, the maximum time step is limited by the smallest oscillation period of the fastest atomic motion, rendering long-time simulations very expensive. To resolve this bottleneck, we propose a supervised parallel-in-time algorithm for stochastic dynamics (SPASD) to accelerate long-time Lagrangian particle simulations. Our method is inspired by bottom-up coarse-graining projections that yield mean-field hydrodynamic behavior in the continuum limit. Here as an example, we use the dissipative particle dynamics (DPD) as the Lagrangian particle simulator that is supervised by its macroscopic counterpart, i.e., the Navier-Stokes simulator. The low-dimensional macroscopic system (here, the Navier-Stokes solver) serves as a predictor to supervise the high-dimensional Lagrangian simulator, in a predictor-corrector type algorithm. The results of the Lagrangian simulation then correct the mean-field prediction and provide the proper microscopic details (e.g., consistent fluctuations, correlations, etc.). The unique feature that sets SPASD apart from other multiscale methods is the use of a low-fidelity macroscopic model as a predictor. The macro-model can be approximate and even inconsistent with the microscale description, but SPASD anticipates the deviation and corrects it internally to recover the true dynamics.

Figures

Equations

$$\frac{\Big{\lVert}\mathcal{G}(\mathscr{F}\{U^{n}_{k+1}\})-\mathcal{G}(\mathscr{F}\{U^{n}_{k}\})\Big{\rVert}_{l_{1}}}{\Big{\lVert}\mathcal{G}(\mathscr{F}\{U^{n}_{k+1}\})\Big{\rVert}_{l_{1}}}\equiv C_{\text{TC}}<C_{\text{Tolerance}},$$

$${\mathrm{T}}^{\text{c}}=N\cdot {\tau }^{\text{c}};$$

$${\mathrm{T}}^{\text{f}}={\tau }^{\text{f}},$$

$$\displaystyle\mathrm{T}^{\text{{SPASD}}}=K\cdot\Big{[}\tau^{\text{f}}+N\cdot(\tau^{\text{c}}+\tau^{\mathscr{F}}+\tau^{\mathscr{R}}+\tau^{\mathscr{P}})\Big{]}.$$

$${\mathrm{T}}_{\text{SD}}^{\text{serial}}=\beta \cdot N\cdot {\tau }^f,$$

$$\beta>\dfrac{K\cdot\Big{[}\tau^{\text{f}}+N\cdot(\tau^{\text{c}}+\tau^{\mathscr{F}}+\tau^{\mathscr{R}}+\tau^{\mathscr{P}})\Big{]}}{N\cdot\tau^{\text{f}}}.$$

$$\displaystyle\beta\gtrapprox K\Big{[}\dfrac{1}{N}+\dfrac{\tau^{\text{c}}}{\tau^{\text{f}}}\Big{]}.$$

$$\frac{d{\mathbf{r}}_i}{dr}=v_i,$$

$$\frac{d{\mathbf{v}}_i}{dt}={\mathbf{F}}_i=\sum _{i\ne j}({\mathbf{F}}_{ij}^C+{\mathbf{F}}_{ij}^D+{\mathbf{F}}_{ij}^R),$$

$${\mathbf{F}}_{ij}^C={\alpha }_{ij}{\omega }_C(r_{ij}){\mathbf{e}}_{ij},$$

$${\mathbf{F}}_{ij}^D=-{\gamma }_{ij}{\omega }_D(r_{ij})({\mathbf{e}}_{ij}\cdot {\mathbf{v}}_{ij}){\mathbf{e}}_{ij},$$

$${\mathbf{F}}_{ij}^R={\sigma }_{ij}{\omega }_R(r_{ij}){\xi }_{ij}\Delta t_f^{-1/2}{\mathbf{e}}_{ij},$$

$${\sigma }_{ij}^2=2{\gamma }_{ij}{k}_{B}T,{\omega }_D(r_{ij})={\omega }_R^2(r_{ij}),$$

$$\frac{\text{D}\mathbf{u}}{\text{D}t}=\nu {\nabla }^2\mathbf{u}+F,$$

$$u(x,t)=\frac{Fd^{2}}{8\nu}\Big{(}1-\Big{(}\frac{2x}{d}^{2}\Big{)}\Big{)}-\sum^{\infty}_{n=0}\frac{4(-1)^{n}Fd^{2}}{\nu\pi^{3}(2n+1)^{3}}\cdot\cos\Big{(}\frac{(2n+1)\pi x}{d}\Big{)}\cdot\exp\Big{(}-\frac{(2n+1)^{2}\pi^{2}\nu t}{d^{2}}\Big{)},$$

$$\displaystyle\epsilon_{l_{2}}\equiv\frac{\Big{\lVert}u^{\text{ref}}-u^{\text{sim}}\Big{\rVert}_{l_{2}}}{\Big{\lVert}u^{\text{ref}}\Big{\rVert}_{l_{2}}},$$

$$\displaystyle\Big{\lVert}u^{\text{ref}}-u^{\text{sim}}\Big{\rVert}_{l_{2}}=\Big{(}\sum_{i}\mathinner{\!\left\lvert u_{i}^{\text{ref}}-u_{i}^{\text{sim}}\right\rvert}^{2}\cdot\Delta x\Big{)}^{1/2},$$

Full text

1 Introduction

Inherent stochastic fluctuations are important in many cellular and subcellular processes in biological systems, as well as in physical systems (soft matter, plasma dynamics, etc.) with timescales across many orders of magnitude. For such systems, a macroscopic deterministic description built on the continuum hypothesis is invalid [1], and therefore Lagrangian methods based on detailed atomic and molecular models are required [2]. However, because the atomistic details are explicitly modeled, the maximum time step is limited by the smallest oscillation period of the fastest atomic motion, which are typically of the order of several femtoseconds [3]. Consequently, the physical time scales accessible to those simulations are too short to address interesting long-time dynamics. This huge disparity in temporal scales, which can also happen at larger values of the smallest timescales, has not been addressed adequately in the computational literature on multiscale physical and biological systems. Consider, for example, the problem of thrombus formation caused by platelet margination and adhesion to the arterial wall, where the characteristic time for platelet activation is a few microseconds. Lagrangian platelet models can only simulate the platelet dynamics up to a few minutes [4, 5], which is far below the time scale for those in vivo processes that usually require at least several hours or even days [6, 7]. With a time step of the order of $10^{-6}~{}s$, a clinically meaningful result for at least one-day time integration may require about $10^{11}$ time steps with explicit platelet models [8], which is computationally prohibitive. Another example is the simulation of protein folding dynamics [9]. In an ab initio molecular dynamics (MD) simulation, the protein molecule and surrounding solvents are treated as classical particles interacting through empirical energy functions. The MD trajectories provide extremely high spatial and temporal resolutions and can identify key intermediates of the folding processes. However, its computational cost hampers the study of large proteins due to their large number of degrees of freedom, greatly increasing the simulation time required to observe a single folding event, which is on the order of tens of microseconds to milliseconds [10]. For timescale on the order of milliseconds, MD simulations are challenging even for the specially designed ASICS machine - the longest continuous-trajectory atomic-scale simulation is a $2.936$ millisecond simulation of NTL9 at $355~{}K$ [11]. Major conformational transitions triggered by enzymatic reactions occur on the millisecond to second time scale or longer [12]. Without algorithmic improvements, long-time integration of protein folding via ab initio MD simulations is a formidable problem.

The complexity due to the huge number of spatio-temporal degrees of freedom needed to resolve atomic or molecular details requires large amounts of computational resources and simulation time to capture the underlying kinetics and dynamics. A typical strategy to leverage parallel computing is dividing the spatial domain into many smaller subdomains so that all subdomains can be processed concurrently by many computing cores. We will refer to this as “conventional spatial decomposition” approach hereafter. Although it is relatively easy-to-implement and very effective, the strong scaling eventually plateaus as the number of unknowns per core becomes too small [13]. The lack of an efficient and stable numerical method remains a serious bottleneck for long-time simulations. Time parallel integration methods can break the bottleneck by decomposing the time domain and solving these subdomains parallel-in-time. Since the inception of time parallel integration methods over 50 years ago [14], many have improved and proposed new parallel-in-time methods. Parareal, proposed by Lions et al. [15], is a popular parallel-in-time scheme and has demonstrated its versatility in various applications. It follows a predictor-corrector paradigm in which the prediction is carried out with an inexpensive solver while the correction is performed with an expensive but parallelizable solver. The prediction is then corrected over iterations, resulting in a refined solution, which then acts as an initial condition in time for the next iteration. The traditional Parareal scheme has been applied to a wide range of problems, including fluid-structure interactions [16], reservoir simulation [17], approximation to Navier-Stokes equations [18], turbulent plasma flow [19], and even to fractional partial differential equations [20]. While most applications employ continuum solvers, a few utilize other methods such as particle solvers [21, 22, 23, 24, 25, 26]. Nevertheless, in most published works the fine and coarse propagators solve the same governing equation in an umbrella model, but with different spatial-temporal discretizations. There are a few exceptions in which the coarse propagator operates on a simplified mathematical model [27, 28, 29, 23]. However, the acceleration of Lagrangian stochastic particle solvers employing parallel-in-time strategy has yet to be developed.

In the present work, we propose a Supervised Parallel-in-time Algorithm for Stochastic Dynamics (SPASD), which aims to significantly accelerate stochastic Lagrangian solvers for long-time simulations. Based on the bottom-up coarse-graining philosophy, stochastic particle models such as dissipative particle dynamics (DPD) converge to continuum macroscopic models in the scale limit [30, 31]. The macroscopic system can then serve as a predictor to supervise the high-dimensional stochastic Lagrangian simulation. Even though the governing equations of the macroscopic model, generally in the form of partial differential equations, are different from those of the microscopic model, the macroscopic model can capture the correct mean-field behavior of the microscopic system in the continuum limit. In particular, an inexpensive continuum solver, also known as the coarse propagator, solves the macroscopic model in serial, and an expensive but parallelizable solver, also known as the fine propagator, resolves the molecular details by performing stochastic microscopic simulations. In the examples, we will show how the coarse propagator can provide a rough prediction of the mean-field hydrodynamics that supervises the fine propagator in the time domain, as well as how the fine propagator corrects the prediction iteratively. For demonstration, we employ DPD to model the stochastic microscopic dynamics and the Navier-Stokes equations to obtain mean-field behavior in the continuum [30, 31]. It is worth noting that SPASD is a scalable and stable parallel-in-time algorithm applicable to many popular stochastic Lagrangian solvers, e.g., MD, DPD, smoothed DPD and Langevin dynamics. Most importantly, the strong scaling of SPASD does not plateau like conventional spatial decomposition methods, and SPASD can be implemented in conjunction with spatial decomposition methods to achieve further speedup.

Considering SPASD as a multiscale method, one could argue that SPASD resembles methods such as multigrid, heterogeneous multiscale method, and equation-free approach. However, there are clear distinctions. Categorically, SPASD belongs to the class of concurrent multiscale modeling, for which the macro- and micro-model are solved concurrently. Structurally, SPASD echoes the extended multigrid methods (EMGM) [32], in which the effective problems are solved with different types of models corresponding to different levels of physical representation. To bridge the scales, EMGM employs techniques that are similar to SPASD - interpolation, relaxation and restriction, and relies on iterative refinement for convergence. However, in contrast with SPASD , EMGM is strictly serial in time. There have been some attempts for developing parallel-in-time multigrid methods [33], but the extension to EMGM has yet to be properly explored to the best of our knowledge.

The heterogeneous multiscale method (HMM) views the connection between the micro- and macro-model from a different perspective. In HMM, one has to first guess the form of the macro-model, which can be incomplete [34], e.g. missing the true constitutive law in complex materials. The key is to extract the required information from the micro-model to solve the incomplete macro-model. As expected, the macro-state provides the constraints for setting up the micro-state. In this regard, the difference between HMM and SPASD is easy to spot: the macro-model is complete and serves as a predictor in SPASD . In particular, guessing the wrong values of the macro-model parameters in SPASD is acceptable because the ramification is corrected by the algorithm, whereas guessing the wrong form of the macro-model is likely going to lead to wrong results in HMM [34].

Another well-known multiscale methodology is the equation-free approach (EFA). EFA is intended to be used when a macroscopic description is unavailable in closed form [35]. In EFA, microscopic simulations are performed only in small spatial domains over a short time period, commonly known as patches. The evolution of the macroscopic fields is then macroscopically interpolated across the patches. In contrast, microscopic simulations in SPASD are performed concurrently on the entire temporal domain. Furthermore, we use a macro-model to best describe the complex microscopic dynamics at a coarse level. Due to the reduced representation of the true dynamics, the macro-model has a closed-form, and it is solvable with a fast solver.

In conclusion, treating the macro-model merely as a predictor is unique to SPASD . That is to say, the macro-model can be inaccurate and inconsistent with the microscale description. Because the parallel-in-time algorithm anticipates the deviation from the true (microscopic) dynamics, SPASD can recover the true dynamics, as demonstrated in this work.

The remainder of this paper is organized as follows. In section 2 we describe SPASD and provide an analysis of its theoretical speedup. In Section 3 we demonstrate the implementation of SPASD for a stochastic Lagrangian simulation supervised by its mean-field approximation. Subsequently, we present quantitative results including accuracy, convergence rate and parallel efficiency for an one-dimensional time-dependent problem. Furthermore, we apply SPASD in a two-dimensional cavity flow and analyzed the results. Lastly, the paper ends with a brief summary and discussion in Section 4.

2 Algorithm

Given that SPASD is inspired by the traditional Parareal algorithm [15], we will first briefly summarize the traditional Parareal and then introduce the proposed algorithm, SPASD. The speedup offered by SPASD is determined by the complexity of projection, mapping and filtering operations as described below. In order to retain generality, we present a theoretical speedup in which the walltimes of those operations are represented symbolically. The theoretical speedup is then compared to the speedup of conventional spatial decomposition method in the context of distributed computing.

2.1 Traditional Parareal Algorithm

For the sake of clarity, let us assume that we have reduced the original problem to the following ordinary differential equation (ODE): ${du}/{dt}=f(u),\;t~{}\in~{}[0,T_{\text{final}}]$; see Table 2.1 for the notation. Let $U^{n}_{k}$ be an approximation to the exact solution $u(n\Delta T)$ at time $n\Delta T$ in $k^{\text{th}}$ iteration, where $\Delta T$ denotes the length of a time subdomain. $N=T_{\text{final}}/\Delta T$ is the number of time subdomains. Let $\mathcal{F}$ be the expensive fine propagator that accurately approximates $u(t)$, and $\mathcal{G}$ be a less accurate but inexpensive coarse propagator. $\Delta t_{f}$ and $\Delta t_{c}$ represent the time step of the fine and coarse propagators, respectively. We note that, although $\Delta t_{c}$ is usually equal to the length of a time subdomain $\Delta T$, it can be less than $\Delta T$ as demonstrated in later sections, i.e., $\Delta t_{c}\leq\Delta T$. Summarized in Algorithm 1, the traditional Parareal algorithm [15] places few constraints on the fine and coarse propagators. Typically, the coarse propagator employs the same numerical scheme as the fine solver, but with a coarser spatio-temporal discretization. Alternatively, they can be different numerical schemes with distinct convergence properties [28, 36]. Nevertheless, both the fine and coarse propagators operate on the same underlying governing equations in most applications.

2.2 Supervised Parallel-in-time Algorithm for Stochastic Dynamics (SPASD)

In order to accelerate particle simulations, we use a mean-field approximation as the macroscopic model whose underlying governing equations are different from those of the microscopic model. Summarized in Algorithm 2 and represented graphically in Fig. 1, propagators in SPASD solve their respective sets of governing equations reflecting the model. As a consequence, a solution representing the macroscopic state is different from the one representing the microscopic state. A macroscopic solution from a coarse propagation can be mapped to a microscopic state via a mapping operation (denoted by the mapping operator $\mathscr{R}$). Conversely, a macroscopic solution can be obtained from a microscopic state via a projection operation (denoted by the projection operator $\mathscr{P}$). The projection and mapping processes are the core challenges for virtually all types of multiscale algorithms. There is no universal method bridging the scales, and the realization of those processes should be problem-dependent. During our numerical experiments, we tested some commonly implemented methods in multiscale such as weighted kernel sampling, linearly interpolation, and nearest neighbor interpolation. For our purpose, we found that the choice of method does not impact the refined results, which is partially due to the fact that the information lost in those processes can be regained by the fine propagator.

Noise filtering, denoted by operator $\mathscr{F}$, is crucial to the success of SPASD because in stochastic simulations the noise would accumulate with time if left unattended. The accumulation of noise will cause the refined solution to diverge from the true solution, as demonstrated in Section 3.1.1.

As with all iterative schemes, the refinement process stops when a termination-condition is satisfied. Because the magnitude of the stochastic noise varies from application to application, the idea of a universal condition is imprudent. We thus propose a practical condition that can be tuned based on the application:

[TABLE]

where $C_{\text{Tolerance}}$ is a user-defined tolerance, which should be determined by considering the magnitude of the intrinsic stochastic noise in the system.

2.3 Theoretical Speedup

The total walltime can be expressed as $\mathrm{T}^{\text{{SPASD}}}=K\cdot(\mathrm{T}^{\text{f}}+\mathrm{T}^{\text{c}}+\mathrm{T}^{\mathscr{F}}+\mathrm{T}^{\mathscr{R}}+\mathrm{T}^{\mathscr{P}})$, where $K$ is the number of iterations to convergence. The walltimes taken by the coarse and fine propagators per iteration are denoted by $\mathrm{T}^{c}$ and $\mathrm{T}^{f}$, respectively. The rest of the terms $\mathrm{T}^{\mathscr{F}}$, $\mathrm{T}^{\mathscr{R}}$ and $\mathrm{T}^{\mathscr{P}}$ denote walltime spent on noise-filtering, up-scaling and down-scaling, respectively. In accordance with the notation in Fig. 2, let $N$ be the number of concurrent fine propagations. $\mathrm{T}^{\text{c}}$ and $\mathrm{T}^{\text{f}}$ can therefore be expressed as

[TABLE]

where $\tau^{\text{f}}$ is the walltime taken by the fine solver for one time-subdomain, and $\tau^{\text{c}}$ is the walltime taken by the coarse solver for the same time-subdomain. $\mathrm{T}^{\mathscr{F}}$, $\mathrm{T}^{\mathscr{R}}$ and $\mathrm{T}^{\mathscr{P}}$ can also be expressed similarly as $\mathrm{T}^{\mathscr{F}}=N\cdot\tau^{\mathscr{F}}$, $\mathrm{T}^{\mathscr{R}}=N\cdot\tau^{\mathscr{R}}$ and $\mathrm{T}^{\mathscr{P}}=N\cdot\tau^{\mathscr{P}}$ respectively, where $\tau$ is the walltime taken by each corresponding operation. The total walltime can thus be written as

[TABLE]

In contrast, the same simulation with only the fine solver running in serial would take $\mathrm{T}^{\text{serial}}=N\cdot\tau^{\text{f}}$. When conventional domain decomposition is applied using the same number of cores, the walltime is reduced to

[TABLE]

where $N^{-1}\leq\beta\leq 1$. For methods that scale perfectly linearly, $\beta$ is equal to $N^{-1}$. In the case that additional resources do not result in speedup, $\beta$ is equal to $1$. Comparing the walltime of SPASD in Eq. (4) with conventional domain decomposition in Eq. (5), SPASD offers better speedup when

[TABLE]

Typically, $\tau^{\text{f}}>>\tau^{\text{c}}$, and both $\tau^{\text{f}}$ and $\tau^{\text{c}}$ are much larger than $\tau^{\mathscr{F}}$, $\tau^{\mathscr{R}}$ and $\tau^{\mathscr{P}}$. Therefore, Eq. (6) can be reduced to

[TABLE]

We note that in order to select the number of time subdomains, we must carefully consider the physics of the underlying stochastic dynamics. In particular, each particle model is associated with a set of dynamical properties, which describe the way collective behavior leads to macroscopic observables. For hydrodynamics, this is given by the characteristic time for momentum decorrelation, which can be quantified by the temporal correlation function of momentum, defined as $C(t)=\langle\mathbf{p}(t+\tau)\mathbf{p}(\tau)\rangle$, where $\mathbf{p}(t)=mass\cdot\mathbf{v}(t)$ represents the instantaneous momentum. The computed $C(t)$ of simple isothermal DPD fluid is shown in Fig. 3, normalized by the value of $C(0)$. The decay of $C(t)$ signifies momentum decorrelation over time, and correspondingly we should use a suitable $\Delta T$ in our Parareal algorithm so that momentum is sufficiently decorrelated. From our benchmark tests, $C(t)=10^{-4}$ provides a sufficient decorrelation of momentum and leads to accurate results. Therefore, in the present work, we use $\Delta T=10$, where the value of $C(t)$ is approximately $10^{-4}$ according to Fig. 3.

2.4 Methods

SPASD offers great versatility and robustness in terms of propagator choices. As a framework targeting particle methods, SPASD does not limit the type of particle solvers; that is to say, any particle solver can be accelerated with SPASD. For demonstration, we model mesoscale hydrodynamics with DPD, a stochastic Lagrangian method. DPD can be derived rigorously from coarse-graining of atomistic dynamics [37, 38] and can recover macroscopic hydrodynamics in the continuum limit [39]. For a simple DPD fluid, a suitable low-dimensional model would be the continuum representation of fluid, namely the Navier-Stokes equations. Since there is no constraint on the discretization of the continuum equations, for simplicity, we chose the finite difference (FD) scheme as the coarse propagator.

The concurrent coupling of two independent solvers requires robust on-the-fly data transfer. For the examples shown here, we used a light weight library called Multiscale Universal Interface [40], or MUI for short, to bridge the solvers.

2.4.1 Dissipative Particle Dynamics

In the DPD framework, each particle is represented explicitly by its position and velocity. The time evolution of DPD particles is governed by Newton’s equation of motion [41]:

[TABLE]

where $t$, $\mathbf{r}_{i}$, $\mathbf{v}_{i}$ and $\mathbf{F}_{i}$ denote time, position, velocity, and force, respectively. The force imposed on particle $i$ is the sum of conservative force $\mathbf{F}^{C}_{ij}$, dissipative force $\mathbf{F}^{D}_{ij}$, and corresponding random force $\mathbf{F}^{R}_{ij}$ from particle $j$ within a radial cutoff $r_{c}$ of $i$. Those pairwise forces are expressed as [42]:

[TABLE]

where $\mathbf{e}_{ij}=\mathbf{r}_{ij}/r_{ij}$ is the unit vector between particles $i$ and $j$, and $\mathbf{v}_{ij}$ is the velocity difference. $\Delta t_{f}$ is the time step of the fine propagator, and $\xi$ is a symmetric Gaussian random variable with zero mean and unit variance [42]. The model parameters $\alpha_{ij}$, $\gamma_{ij}$ and $\sigma_{ij}$ adjust the conservative, dissipative, and random forces, respectively. The corresponding weighting function $\omega_{C}(r_{ij})$, $\omega_{D}(r_{ij})$ and $\omega_{R}(r_{ij})$ regulate the interaction between pairs of particles. The balance between dissipation and thermal fluctuation is maintained by constraints from the fluctuation-dissipation theorem [43]:

[TABLE]

where $k_{B}T$ is the Boltzmann energy unit. A common choice for the weight function is $\omega_{C}(r)=1-r/r_{c}$ and $\omega_{D}(r)=\omega^{2}_{R}(r)=(1-r/r_{c})^{2}$ for $r\leq r_{c}$ and zero for $r>r_{c}$, which is used for the simulations shown in this paper. For simplicity, $k_{B}T$ and the mass of a particle are taken as the energy unit and mass unit, and their values are set to one.

2.4.2 Model Coupling

By applying a Mori-Zwanzig projection operator technique to a DPD system, Español [30] and Marsh et al. [44] formally derived the hydrodynamic equations of a simple DPD fluid using the Green-Kubo formulas, recovering the continuity and Navier-Stokes equations. For more details on the mathematical derivations leading to macroscopic hydrodynamic equation from DPD dynamics, we refer the interested readers to a paper by Español [30] and a book by Zwanzig [31]. Therefore, the continuity and Navier-Stokes equations can be considered as a mean-field representation of the DPD system, and consequently provide macroscopic solution as a rough prediction of the hydrodynamics that can be used to supervise the DPD simulations in the time domain.

For coupling the grid-based Navier-Stokes solver with the Lagrangian DPD solver, we found MUI easy to implement in the sense that codes do not need to be refactored before application. Adding only a few lines of code to existing solvers is sufficient to adhere to MUI’s push-fetch workflow - data pushed by one solver is then fetched by the other. In the context of particle solvers, the relevant data are particle positions and the associated quantities. For continuum solvers, the quantities of interest are associated with nodes of the discretization grid. In a data-transfer, the sender pushes the data into a MUI Interface Layer (IF), which then assists the receiver to fetch and decode these data based on the topological configuration. Because topological configuration in the receiver usually do not exactly overlap with these of the particle solver, a sampler such as Gaussian kernel or nearest neighbor might be used to calculate an appropriate nodal value in the receiver. As for the linking configuration, each of the fine propagators is attached to one IF, and all of those IFs link to the coarse propagator.

3 Numerical Results

In our proposed framework, the macroscopic model is a mean-field approximation to our microscopic model. As a consequence, the physical parameters in the macroscopic model cannot be calculated or deduced precisely from the microscopic model. However, the macroscopic model in SPASD only requires rough estimations of these physical parameters, effects of which can be iteratively “corrected” by the correction term in Algorithm 2. In the example problems below, we also demonstrate empirically that a more accurate parameter estimation in the macroscopic model leads to faster convergence. In addition, we will analyze the accuracy and convergence of SPASD on the example problems.

3.1 One-dimensional time-dependent flow

As a demonstration, we applied SPASD to the one-dimensional time-dependent plane Poiseuille flow with DPD and FD as fine and coarse propagators, respectively. Simulating Poiseuille flow is a well-accepted benchmark test for many Lagrangian methods including DPD because the macroscopic parameter viscosity can be computed empirically from simulations [45]. We refer to the computed viscosity as the true viscosity $\nu_{\text{true}}$, which is incorporated in the true solution that is used in error analysis.

Macroscopically, the pressure-driven Poiseuille flow can be described by the incompressible Navier-Stokes equations:

[TABLE]

where $\nu$ is the viscosity, $F$ denotes a body-force driven by a constant pressure gradient, and $\mathbf{u}$ denotes the velocity field. The exact time-dependent solution to the Navier-Stokes equation is known analytically [46]:

[TABLE]

The true solution $u(\nu_{\text{true}})$ can thus be readily calculated according to Eq. (15). For this problem, we use a periodic box with the size of $30\times 20\times 10$ in DPD units containing $24,000$ particles. The DPD parameters $\alpha$ and $\gamma$ are set to $18.75$ and $4.5$, respectively. All particles are subject to a body-force $F=0.1$ and a radial cutoff $r_{c}=1.58$. Given these parameters, the true viscosity for this system is approximately $0.841$ computed via DPD simulations.

As explained in Section 2.3, the momentum is sufficiently decorrelated beyond $10$ time units according to the temporal correlation function $C(t)$. The communication interval between the propagators $\Delta T$, also known as the time step of the coarse propagator, is thus set to be $10$. With the time step $\Delta t_{f}=0.01$, the fine propagators march forward for $1000$ time steps during the same period. Because the Poiseuille flow is one-dimensional, we use a simple linear interpolation as the projection operation. To find the velocity profile, the simulation domain is divided into slabs of uniform width that are also the size of spatial discretization in the coarse propagator. We then average particle velocities in each slab. The result is a one-dimensional velocity profile that describes the macroscopic state. On the other hand, the macroscopic state is mapped to the particle system with a nearest-neighbor sampler - a particle takes on the velocity of the nearest grid node. The coarse propagation and the filtering of stochastic noise generated by the fine propagator are carried out simultaneously in this example. The viscosity of fluid can naturally create a low-pass filtering effect [47] that damps high-frequency fluctuations as the coarse solver advances. To take advantage of the filtering effect, one time subdomain $\Delta T=10$ is covered by $100$ sequential coarse steps with an effective time step $\Delta t_{c}=0.1$.

3.1.1 Impact of Filtering

As an iterative scheme, the convergence of SPASD greatly relies on the notion that, for a propagation, the stochastic noise in the model should not cause the solution to deviate far from the ensemble solution. To ensure convergence, a filtering operation must be incorporated in the algorithm. If the filtering operation is not enforced, the noise can grow in the temporal space within an iteration and will propagate across iterations. Unchecked noise growth can eventually drive the propagated solution away from the true solution. To demonstrate the impact of filtering, we compared results of the Poiseuille flow simulations without and with partial filtering. For fluids, the macroscopic model has a filtering mechanism build-in: the viscosity of fluid is a natural low-pass filter that damps high-frequency fluctuations. To take advantage of the filtering effect, a time subdomain $\Delta T$ is divided into smaller (effective) time steps $\Delta t_{c}$ which the coarse solver uses.

To demonstrate the impact of filtering, we run a simple example without and with partial-filtering. The same effective time step, subsequently the same Courant-Friedrichs-Lewy condition, are used in both simulations to isolate the effect of filtering. For a time subdomains of $10$ time units, the coarse propagator with an effective time step of $1$ unit filters the solution 10 consecutive times. This approach offers better noise reduction than the case with a subdomain size of $1$ unit and an effective time step of $1$ unit, which offers no filtering. The impact of filtering is apparent in this comparison: the solution in the no-filtering case in Fig. 4LABEL:sub@wofp is substantially nosier than the one in the partial-filtering case in Fig. 4LABEL:sub@wfp. To quantify the deviation from the reference solution, we compute the normalized $l_{2}$ error, expressed as

[TABLE]

where $\Delta x$ denotes the size of spatial discretization, and $i$ is the nodal index. As shown in Fig. 4LABEL:sub@fl2time, the normalized $l_{2}$ error $\epsilon_{l_{2}}$ for the no-filtering case is monotonically increasing with time, which signifies the growth of noise in time within an iteration. The noise also propagates across iterations shown in Fig. 4LABEL:sub@fl2. For the partial-filtering case, the noise accumulates at a much slower rate than the no-filtering case.

3.1.2 Accuracy

To demonstrate the flexibility and robustness of SPASD, we purposely estimated the viscosity in the macroscopic model to be $10\times$ the true viscosity. In other words, the estimated mean-field solution $u(\nu_{\text{est}})=u(\nu_{\text{true}})/10$. To verify the accuracy of SPASD, we compare the true velocity profiles with the ensemble average obtained with multiple SPASD runs. Because the microscopic model is stochastic in nature, it is desirable to obtain an ensemble averaged solution that is less distorted by stochastic noise. In our experiment, the simulated velocity profile is ensemble-averaged from fifty independent simulations. The resulting transient velocity profiles, as well as the evolution of velocity profile at steady state, match those of the true velocity profiles visually as shown in Fig. 5. To quantify the accuracy of SPASD , we computed the normalized $l_{2}$ error $\epsilon_{l_{2}}$ between the reference and simulated solutions. For the one-dimensional Poiseuille plane-flow, the reference solution is the true solution $u(\nu_{\text{true}})$ calculated according to Eq. (15) because the analytical form is known. We found that $\epsilon_{l_{2}}$ decreases over iterations as shown in Fig. 6LABEL:sub@fig:Poise_l2norm. Given that the noise cannot be eliminated by averaging fifty solutions, we accept solutions within $1\%$ error as a threshold denoted as $\epsilon_{\text{Threshold}}$. In this particular example, $\epsilon_{l_{2}}$ falls to $\epsilon_{\text{Threshold}}$ in $20$ iterations. At the threshold, the error has a Gaussian distribution as shown in Fig. 6LABEL:sub@fig:Poise_distri, which matches the noise distribution from the microscopic model [42]. The recovery of the distribution signifies that the SPASD preserved stochastic fluctuations of the microscopic model.

3.2 Convergence

Inspecting $\epsilon_{l_{2}}$ at $T_{\text{final}}=320$ for all iterations in Fig. 7LABEL:sub@fig:rateofconv, SPASD achieves exponential convergence. We repeat the experiment with a more accurate macroscopic viscosity at $2\times$ true viscosity and obtain exponential convergence as well, but only $5$ iterations are needed to reach the same threshold. More importantly, $\epsilon_{l_{2}}$ reaches a minimum level regardless of the estimated viscosity. The minimum error is determined by many factors; in addition to the error intrinsic to Parareal, operations in SPASD such as filtering and projection introduce systematic errors that are difficult to separately quantify.

Lastly, we calculate $C_{\text{TC}}$ from Eq. (1) and plot it in Fig. 7LABEL:sub@fig:convthreshold. Because of the stochastic noise, $C_{\text{TC}}$ varies between runs. We plot the mean and $95\%$ confidence interval for $C_{\text{TC}}$ generated from fifty independent runs. When the parameter estimation is more accurate as in the case of 2x true viscosity, $C_{\text{TC}}$ fluctuates mildly as indicated by the range of the confidence interval and drops sharply to a level that satisfies the termination condition. However, when estimation is rough, the confidence interval is relatively wide without a clear boundary separating iterations above and below $\epsilon_{\text{Threshold}}$. The resulting consequence is discussed in Section 4.

3.3 Efficiency

For particle simulations, spatial domain decomposition is the most utilized acceleration method because the subdomains can be easily assigned to different processors. To compare the performance of SPASD with conventional spatial decomposition (CSD), we apply them separately to the planer Poiseuille flow problem and then calculate the parallel efficiency.

The same simulation parameters described in Section 3.1 are used in this benchmark test. To avoid implementing no-slip boundaries, we simulated a reverse Poiseuille flow [48, 49], in which two planar flows are stacked to create periodic boundaries in all directions. At the same time, the $z$-direction is shrunk by one half to preserve the total number of particles. For a general comparison, we simulate a single iteration, and each core was responsible for $10$ time units. When the number of available cores is doubled, the final time $T_{\text{final}}$ is also doubled.

This particular benchmark was run on nodes with AMD Opteron 6274 CPUs. For this example, the performances of CSD and SPASD are similar when fewer cores are available as shown in Fig. 8LABEL:sub@Walltime. As the core-count and $T_{\text{final}}$ increase, their performances start to diverge at approximately $128$ cores. From there, the walltime of CSD grows at a much faster exponential rate than that of SPASD. Two factors contribute to this: First, the size of a spatial subdomain in CSD is limited by the average particle distance. Second, the time spent on communication between nodes grows at a rate that nullifies the benefits of additional cores. SPASD, however, does not suffer from the same issues. Communication between nodes is only invoked at the end of temporal subdomain, and the duration of temporal subdomain can be adjusted at will. This makes SPASD a more scalable algorithm suitable for long-time simulations. For comparison, we used the speed of a single-core simulation, roughly $39$ seconds per time unit, as a basis and computed the speedup of each method. The efficiencies were then calculated as the speedups per core. As expected, SPASD offers better parallel efficiency in the cases of large core counts, as shown in Fig. 8LABEL:sub@Efficiency.

The walltime taken by the coarse propagator and associated operations $\tau^{\text{c}}+\tau^{\mathscr{F}}+\tau^{\mathscr{R}}+\tau^{\mathscr{P}}$, as well as the walltime taken by the fine propagator $\tau^{\text{f}}$, can be readily calculated with a simple linear fit to Eq. (4). For this benchmark case, the walltime of the simulation using SPASD is described by $\mathrm{T}^{\text{{SPASD}}}=197.8+0.08783N$ where $N=T_{\text{final}}/10$. Plugging in the numerical value for $\tau^{\text{f}}$, the walltime of simulation using CSD is reduced to $\mathrm{T}^{\text{serial}}_{\text{SD}}=197.8\beta N$. The scaling factor $\beta$ depends on many factors including the number of parallel cores, distribution topology, and parallel-method efficiency. When few cores are available, $\beta$ is approximately equal to $N^{-1}$ as shown in Fig. 8LABEL:sub@Efficiency, and CSD provides better efficiency and performance. As the number of cores increases, $\beta$ decreases proportionally until around 128 cores at which SPASD and CSD have similar parallel efficiencies. At this core-count, $\beta$ was calculated to be $0.0083$, which agrees with our prediction $\beta_{c}\equiv K(N^{-1}+\tau_{c}/\tau_{f})=0.0084$ given by Eq. (7). From this core-count forward, the flattening of the $\beta$ curve signifies the drop of CSD efficiency.

3.4 Two-dimensional time-dependent flow

To further demonstrate the applicability of SPASD, we simulated a two-dimensional cavity flow. Simulating transient cavity flow is much more complicated than the Poiseuille flow because the convective acceleration and the corner singularity induced by the cavity make it numerically challenging. However, the singularities, including corner singularity [50], crack propagation [51], and moving contact line singularity [52], in partial differential equation (PDE) models (in the limit of $\Delta x\rightarrow 0$) represent rich underlying multiscale physics, and can be captured by particle solvers, e.g., DPD or MD.

In this case, we again use DPD to simulate the lid-driven cavity flow, which is supervised in the temporal domain by the incompressible Navier-Stokes equations. These equations are solved with finite difference on a staggered grid: the pressure is stored at the center of the grid cells, and the velocities are stored at the faces of the grid cells. The second order spatial derivative is approximated with a centered difference scheme. This setup allows pressure and velocity to be solved concurrently in a two-step process. The first step is to compute an intermediate velocity by solving the momentum equation, omitting the effect of pressure. The intermediate velocity is then used to compute new pressure by solving a pressure Poisson equation. The second step is to solve for the new velocity also using the intermediate velocity and pressure. For temporal discretization, we use explicit time stepping.

A fluid region with a dimension of $30\times 30\times 40$ DPD units is surrounded by a solid wall consisted of fixed particles. At a particle density of $8$, we set the DPD parameter $\alpha=9.375$ and $\gamma=13.5$ to counter the effect of compressibility due to the fast moving wall at a velocity of $1.833$. This set of parameters results in a reference Reynolds number $\text{Re}_{\text{ref}}=200$. To create the no-slip cavity lid and wall, we used an effective dissipative coefficient for liquid–solid interactions [53]. For the cavity flow, we used the same projection method as in the Poiseuille flow example, with the consideration that the velocity profile consists of two components. As a result, both $x$ and $y$ components are recorded for each grid node. Because the grid in the coarse propagator spans in two dimensions, the mapping procedure for this example is slightly modified from that of plane-flow. The velocity of a particle is the linear interpolation of two nearest grid values. We used the same filtering techniques as in the prior examples.

Because an analytical solution does not exist for cavity flow, a reference solution is constructed by averaging the velocity fields from fifty independent serial DPD simulations. For the purpose of analysis, instead of particle velocity fields we inspect the equivalent stream function. The stream function can be easily obtained by solving the Poisson equation $-\nabla^{2}\psi=(\nabla\times\bm{u})\cdot\hat{z}$. We then compute the normalized $l_{2}$ error $\epsilon_{l_{2}}$ between simulated and reference stream functions.

To test SPASD’s flexibility, we falsely assume that the fluid is much less viscous in the macroscopic model with a Reynolds number $\text{Re}_{\text{est}}=600$ which is $3\times$ the reference Reynolds number $\text{Re}_{\text{ref}}$. It is important to note that the falsely assumed viscosity would lead to an incorrect vortex location and, consequently, vastly different velocity field. As a result, the vortex from our SPASD simulation does not align with the reference vortex in early iterations shown in Fig 9. The mis-alignment is then iteratively corrected. We again calculate $\epsilon_{l_{2}}$ at $T_{\text{final}}$ and plot it against iterations in Fig. 10LABEL:sub@l2error. The rate of convergence is no longer exponential as it was for the Poiseuille plane-flow. This is due to the notion that in two-dimensional space the $l_{2}$ error might not decay in a shortest path. For this particular case, it took 5 iterations to reach the threshold error $\epsilon_{\text{Threshold}}=1\%$ as shown in Fig. 10LABEL:sub@l2error. We also repeated the experiment with an estimated Reynolds number $\text{Re}_{\text{est}}$ that is equal to the reference $\text{Re}_{\text{ref}}$. As expected, it required fewer iterations to converge. The termination criterion was also computed and plotted in Fig. 10LABEL:sub@Conver. It exhibits the same shortcoming as in the Poiseuille flow - no clear separation between iterations above and below $\epsilon_{\text{Threshold}}$. Further discussion on this issue is presented in Section 4.

4 Summary & Discussion

We have developed a supervised parallel-in-time algorithm for stochastic dynamics (SPASD). As an extension to the Parareal algorithm, SPASD uses heterogeneous solvers, e.g., a Navier-Stokes solver as a predictor and the DPD method as a corrector, to accelerate stochastic particle simulations by decomposing the temporal domain. The predictor, also known as the coarse propagator, solves the macroscopic model, which is generally in the form of partial differential equations, in serial. Because the microscopic model converges to continuum macroscopic models in the scale limit, the time evolution of the macroscopic system then serves as an initial condition for all time subdomains. The corrector, also known as the fine propagator, recovers the micro-dynamics with expensive stochastic simulations in parallel. While here we demonstrated SPASD for hydrodynamics, the same framework can be applied to other problems, where long-time integration is hampered by the small characteristic time scale of the micro-dynamics.

To summarize, we first presented the SPASD algorithm and analyzed its theoretical speedup. We then discussed the importance of filtering and demonstrated its effect with an example problem. The accuracy and convergence of SPASD was subsequently tested. We showed that the transient and final solutions match the analytical/reference solutions even when our estimations to the system’s mean-field behavior is far from the true behavior. More importantly, the SPASD algorithm is able to preserve the stochastic fluctuations of the microscopic model. Lastly, we demonstrated that SPASD provides better parallel efficiency and thus better scalability than the conventional domain decomposition method for long-time simulations.

There are three comments and suggestions we would like to make regarding SPASD. First, although the solution converges regardless the value of estimated viscosity, an overly rough estimation requires a large number of iterations. As for the Poiseuille flow example, a rough estimate at $10\times$ true viscosity requires $20$ iterations, whereas only $5$ iterations are needed for a better estimate at $2\times$ the true viscosity. In light of this correlation, we strongly encourage using an informed estimate of the macroscopic parameters in order to maximize the efficiency. Second, if the data from one microscopic state is insufficient to determine a mean-field quantity, the projection operation can be performed with data from multiple states with temporal-proximity. With the additional data, the mean-field quantity can be determined with higher precision. Lastly, we would like to discuss the issue of the termination condition. Given that SPASD is a general parallel-in-time algorithm, the termination condition presented in Section 2.2 should be application-dependent. The stochastic noise is embedded in the refined solution. $C_{\text{Tolerance}}$ thus should be adjusted according to the magnitude of the noise. In the examples, we chose the normalized $l_{2}$ error as the quantity of interest $C_{\text{TC}}$. The termination criterion is satisfied when $C_{\text{TC}}$ falls below a threshold (see Eq. (1)). For the Poiseuille flow, we computed $C_{\text{TC}}$ for all iterations as shown in Fig. 6LABEL:sub@fig:convthreshold. When the parameter estimation is accurate as in the case of $2\times$ reference viscosity, $C_{\text{TC}}$ fluctuates mildly as indicated by the range of confidence interval and drops sharply to a level that satisfies the termination condition. However, when estimation is rough, the confidence interval is relatively large, which might cause the simulation to stop prematurely. In this case, we recommend using noise-filtered solution in the error calculation. Another option is to define a different quantity of interest as the termination criterion with a redefined termination condition. Finally, it is important to mention that the use of filtering helps to damp the high-frequency noises in the mean-field quantities, which only acts as a rough predictor for the fine stochastic solver. Regardless of the choice of filtering algorithms, the final solution obtained with SPASD contains all the statistics of the stochastic fluctuations and their correlations, consistent with unsupervised (serial in time) stochastic dynamics.

In future work, we would like to expand the applications of SPASD in two directions. First, we plan to simulate multi-physics dynamics along with hydrodynamics. For example, we can investigate stochastic non-isothermal systems at the microscale by incorporating the temperature equation in the macro-model. The temperature equation and the Navier-Stokes can be solved independently. For the micro-model, we can employ the energy dissipative particle dynamics, a flavor of DPD that tracks the energy of the system explicitly [49]. Second, we will consider complex fluids that exhibit non-Newtonian behavior while employing a simple fluid model in the macro-scale. For example, we can model blood flow in the microscale, in which the non-Newtonian behavior can be captured with explicit representations of red blood cells and white blood cells [54]. These are currently very expensive simulations and the aim is to test whether SPASD can speed up such simulations by orders of magnitude while maintaining the fidelity of microscale simulations.

Acknowledgements

This work was supported by the DOE PhILMs project (No. DE-SC0019453) and the U.S. Army Research Laboratory under Cooperative Agreement No. W911NF-12-2-0023. An award of computer time was provided by the ASCR Leadership Computing Challenge (ALCC) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. The authors would like to thank Dr. Fangying Song for providing technical support on data analysis and Prof. Yvon Maday for insightful discussions.