Machine learning acceleration of simulations of Stokesian suspensions

TL;DR

This paper introduces a machine learning-augmented reduced model that significantly accelerates simulations of Stokesian suspensions, maintaining accuracy across various geometries and boundary conditions.

Contribution

It presents a novel ML-based reduced model that replaces computationally expensive parts of fluid simulations, enabling faster and more flexible analysis of particulate Stokesian flows.

Findings

Model is 10 times faster than traditional schemes.

Accurately reproduces key flow features.

Generalizes to arbitrary geometries without retraining.

Abstract

Particulate Stokesian flows describe the hydrodynamics of rigid or deformable particles in Stokes flows. Due to highly nonlinear fluid-structure interaction dynamics, moving interfaces, and multiple scales, numerical simulations of such flows are challenging and expensive. In this Letter, we propose a generic machine-learning-augmented reduced model for these flows. Our model replaces expensive parts of a numerical scheme with multilayer perceptrons. Given the physical parameters of the particle, our model generalizes to arbitrary geometries and boundary conditions without the need to retrain the regression function. It is 10 times faster than a state-of-the-art numerical scheme having the same number of degrees of freedom and can reproduce several features of the flow quite accurately. We illustrate the performance of our model on integral equation formulation of vesicle suspensions in two dimensions.

Equations

$$\frac{\partial \mathbf{x}}{\partial t}={\mathbf{u}}_{\mathrm{self}}[\gamma ](\mathbf{x})+V[\gamma ]{\mathbf{u}}_{\infty }(\mathbf{x}),$$

$$\frac{\partial \mathbf{x}}{\partial t}=V[{\gamma }^0]{\mathbf{u}}_{\infty }(\mathbf{x});\mathbf{x}(0)={\mathbf{x}}^0\in {\gamma }^0;t\in (0,\Delta t),$$

$$\frac{\partial \mathbf{x}}{\partial t}={\mathbf{u}}_{\mathrm{self}}[{\gamma }^{\ast }](\mathbf{x});\mathbf{x}(0)={\mathbf{x}}^{\ast }\in {\gamma }^{\ast };t\in (0,\Delta t).$$

$$V[\gamma ]{\mathbf{u}}_{\infty }(\mathbf{x})=\sum _{k=1}^{N_f}V[\gamma ]{\varphi }_k(\alpha ){\hat{\mathbf{u}}}_{\infty ,k}=\sum _{k=1}^{N_f}{\Psi }_k[\gamma ](\alpha ){\hat{\mathbf{u}}}_{\infty ,k}.$$

$$\frac{\partial \mathbf{x}}{\partial t}=V[{\gamma }_i^0]{\mathbf{u}}_{\infty }[\gamma ,\Gamma ](\mathbf{x});\mathbf{x}(0)={\mathbf{x}}^0\in {\gamma }_i^0,$$

$$-\mu \Delta \mathbf{u}(\mathbf{x})+\nabla p(\mathbf{x})=0,\text{and}\nabla \cdot \mathbf{u}(\mathbf{x})=0,\mathbf{x}\in \Omega \setminus (\bigcup _i{\gamma }_i),$$

$$\mathbf{u}(\mathbf{x},t)=\frac{\partial \mathbf{x}}{\partial t}(t),\mathbf{x}\in \{{\gamma }_i{\}}_{i=1}^M\mathrm{and}\mathbf{u}(\mathbf{x},t)=\mathbf{U}(\mathbf{x},t),\mathbf{x}\in \{{\Gamma }_k{\}}_{k=0}^{K-1}.$$

$${\mathbf{x}}_s\cdot {\mathbf{u}}_s=0,\mathbf{x}\in \{{\gamma }_i{\}}_{i=1}^M.$$

$$[[\mathbf{T}\mathbf{n}]]=-{\kappa }_b{\mathbf{x}}_{ssss}+(\sigma {\mathbf{x}}_s{)}_s,\mathbf{x}\in \{{\gamma }_i{\}}_{i=1}^M,$$

$$\mathcal{S}[\mathbf{f}](\mathbf{x})={\int }_{\gamma }\mathbf{G}(\mathbf{x}-\mathbf{y})\mathbf{f}d\gamma (\mathbf{y}),$$

$${\mathbf{G}}_{ij}(\mathbf{x}-\mathbf{y})=\frac{1}{4\pi \mu }\left(-{\delta }_{ij}\ln |\mathbf{x}-\mathbf{y}|+\frac{(\mathbf{x}-\mathbf{y}{)}_i(\mathbf{x}-\mathbf{y}{)}_j}{|\mathbf{x}-\mathbf{y}{|}^2}\right).$$

$$\mathcal{B}(\mathbf{x})\mathbf{x}=-{\kappa }_b\mathcal{S}[{\mathbf{x}}_{ssss}]\mathbf{x},$$

$$\mathcal{T}(\mathbf{x})\sigma =\mathcal{S}[\sigma {\mathbf{x}}_s]\mathbf{x},$$

$$\mathcal{P}(\mathbf{x})\mathbf{y}={\mathbf{x}}_s\cdot {\mathbf{y}}_s,$$

$$\mathcal{M}(\mathbf{x})=\mathcal{T}(\mathbf{x}){\left(\mathcal{P}(\mathbf{x})\mathcal{T}(\mathbf{x})\right)}^{-1}\mathcal{P}(\mathbf{x}).$$

$$\frac{\partial \mathbf{x}}{\partial t}=(1-\mathcal{M}(\mathbf{x}))\mathcal{B}(\mathbf{x})\mathbf{x}+(1-\mathcal{M}(\mathbf{x})){\mathbf{u}}_{\infty }(\mathbf{x}).$$

$$\frac{\partial \mathbf{x}}{\partial t}=(1-\mathcal{M}(\mathbf{x}))\mathcal{B}(\mathbf{x})\mathbf{x}+(1-\mathcal{M}(\mathbf{x})){\mathbf{u}}_{\infty }[\gamma ,\Gamma ](\mathbf{x}),$$

$${\mathbf{u}}_{\infty }[\gamma ,\Gamma ](\mathbf{x})=\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^M{\int }_{{\gamma }_j}\mathbf{G}(\mathbf{x}-\mathbf{y})\mathbf{f}(\mathbf{y})d{\gamma }_j(\mathbf{y})+\mathcal{W}[\eta ](\mathbf{x}),$$

$$\mathcal{P}(\mathbf{x})\mathcal{T}(\mathbf{x})\sigma (\mathbf{x})=-\mathcal{P}(\mathbf{x})\left(\mathcal{B}(\mathbf{x})\mathbf{x}+{\mathbf{u}}_{\infty }[\gamma ,\Gamma ](\mathbf{x})\right),\mathbf{x}\in {\gamma }_i.$$

$$\mathbf{U}(\mathbf{x})=-\frac{1}{2}\eta (\mathbf{x})+\sum _{j=1}^M{\int }_{{\gamma }_j}\mathbf{G}(\mathbf{x}-\mathbf{y})\mathbf{f}(\mathbf{y})d{\gamma }_j(\mathbf{y})+\mathcal{W}[\eta ](\mathbf{x}),\mathbf{x}\in \Gamma .$$

$$\sigma (\mathbf{x})={\sigma }_{\mathrm{self}}[{\gamma }_i](\mathbf{x})+T[{\gamma }_i]{\mathbf{u}}_{\infty }[\gamma ,\Gamma ],\mathbf{x}\in {\gamma }_i,$$

$$\mathbf{X}=\mathbf{W}\mathbf{\Sigma }{\mathbf{V}}^T$$

$${\mathbf{C}}_{N_p}={\mathbf{W}}_{N_p}^T\mathbf{X}$$

$$\mathcal{J}=\sum _{j=1}^M\Vert {\mathbf{x}}_j^{\Delta t}-{\tilde{\mathbf{x}}}_j^{\Delta t}{\Vert }_2,$$

$$\mathcal{J}=\sum _{j=1}^M\Vert ({\Psi }_k{)}_j-({\tilde{\Psi }}_k{)}_j{\Vert }_2.$$

$$z_2=f(W_1{z}_1+b_1),z_3=f(W_2{z}_2+b_2),\cdots ,z_{K+1}=f(W_K{z}_K+b_K).$$

Full text

*Introduction.—*Particulate Stokesian flows consider the motion of a collection of rigid or deformable particles (e.g., drops, capsules, cells, slender bodies, filaments, active swimmers, possibly elastic or filled by a fluid) that are suspended in a Newtonian fluid and the particle Reynolds number is vanishingly small Lauga and Powers (2009); Shelley and Zhang (2011); Happel and Brenner (2012); Barthés-Biesel (2012). Such flows find many applications in industrial processes, microfluidics, study of complex fluids, bacterial and general active flows. Due to the multiple scales, the strongly nonlinear and nonlocal coupling of the interface deformation to the background flow, and the need for long time horizons, numerical simulations, no matter what the underlying numerical method is, can be extremely expensive. In this letter, we propose a machine-learning-augmented reduced model (MLARM) that consists of three components:

** (i)** We use multilayer perceptrons (MLPs) to approximate several spatial nonlinear operators in our numerical scheme. We choose these operators based on their stiffness properties, computational costs, and ability to generalize to unseen data.

** (ii)** We use high-fidelity (i.e., highly refined in space and time) simulations to train the MLPs. We run these high-fidelity simulations in “burst mode”, single-particle simulations for one time step. These simulations are done only once per particle type (e.g., for various mechanical properties of a particle) and target accuracy (i.e., time step size).

** (iii)** We combine the MLPs with low-fidelity simulations. Although we trained using single-particle, short-horizon, unconfined flows, our method enables us to conduct long-horizon simulations of suspensions with several particles in confined geometries.

The basic idea is to create a regression function that accurately captures the dynamics of the flow using the high-fidelity simulations; and then use this function in the low-fidelity numerical scheme. There exist many schemes for function approximation in high dimensions; we have opted for a multilayer perceptron Hornik et al. (1989). We model different aspects of a flow using separate MLPs. We train a MLP on single particle dynamics in a particle relaxation regime (no imposed flow); and we train another one for the evolution of boundary due to imposed velocity using the mode decomposition of its restriction on the boundary of the particle. Once we have trained the MLPs, the framework can be applied to any flow configuration, that is, we can vary the confined geometry, the number of particles, and the imposed boundary conditions (confined or free-space) without the need to retrain the networks. We demonstrate MLARM’s capability of accurately capturing microscopic and macroscopic flow characteristics to the low-resolution simulations having the same number of degree of freedom for several benchmark problems.

A popular and effective methodology for the mathematical analysis and numerical simulation of Stokesian particulate flows is the boundary integral method (see Rallison and Acrivos (1978); Pozrikidis (1992, 2001); Wang et al. (2013); Freund (2014) for various examples). Our overall method is based on this underlying formulation. To demonstrate our framework, we select a specific particle type, vesicle Seifert (1997); Kantsler and Steinberg (2005); Misbah (2012); Li et al. (2013). (Although the details depend on a particle type, our model is broadly applicable to boundary integral formulations for Stokesian particulate flows). These flows are quite challenging because vesicles are deformable particles that resist bending and stretching, and are filled with a fluid. Vesicles resemble red blood cells and their simulations have been used to help understanding microcirculation and to design microfluidic devices for medical diagnoses and drug delivery systems. Nevertheless, the high computational cost makes simulations of vesicles at realistic concentrations extremely expensive (sometimes they can take several days on large clusters) and limits the ability to perform design optimization and parameter sweeps. Although a numerical scheme based on low resolution discretization in space and time along with a set of correction algorithms enables simulating dense vesicle flows for long time horizons Quaife and Biros (2014); Kabacaoğlu et al. (2018) and solving a design optimization problem Kabacaoğlu and Biros (2018), such simulations still demand even faster algorithms Zhao et al. (2010); Freund (2014); Malhotra et al. (2017). Let us note that although numerous works have used machine learning to tackle computational physics problems there has been little work on Stokesian flows. Most closely related to ours have proposed a reduced order model for unsteady quasi-one dimensional Euler flows Wang et al. (2019). Reduced basis functions are extracted from high-fidelity solutions and deep neural networks are used to map the flow parameters to the coefficients of the reduced basis functions. But, to our knowledge there exist no works on reduced models for Stokesian particulate flows.

*MLARM for vesicle flows.—*We assume that fluids in the interior and the exterior of the vesicles are Newtonian with the same viscosity. We consider two-dimensional flows and model them using a quasi-static Stokes approximation scheme Kraus et al. (1996); Danker et al. (2009). Consider a single vesicle in a free-space flow (Fig. 1(a)). Let $\gamma$ denote a vesicle membrane. The points ${\mathbf{x}}(\alpha,t)$ on $\gamma$ are given at $N$ uniformly distributed points $\{\alpha_{k}=2\pi(k-1)/N\}_{k=1}^{N}$ in the parametric domain and time $t$. The velocity of the membrane points is given by

[TABLE]

where $\mathbf{u}_{\mathrm{self}}[\gamma]({\mathbf{x}})$ is the velocity induced by the vesicle itself and $V[\gamma]\mathbf{u}_{\infty}({\mathbf{x}})$ is the matrix-vector multiplication that evaluates the effect of the imposed flow $\mathbf{u}_{\infty}$ on the membrane velocity. We refer to matrices with upper case letters and to vectors with lower case bold letters; $[\gamma]$ denotes dependence (nonlinear and nonlocal) on the interface $\gamma$, and it involves solution of boundary value problems (in our case with integral equations).

Although in our prior work we have found linearly implicit time-stepping schemes to work well Kabacaoğlu et al. (2018), here we use a different method that allows modular MLP training. We use an operator-splitting method Macnamara and Strang (2016) for (1), which brings several advantages when using MLPs (will appear later). The method divides (1) into two problems: advection and relaxation. Given the vesicle membrane $\gamma^{0}$, we first obtain the membrane $\gamma^{*}$ at the intermediate step by evaluating the membrane velocity due to the imposed flow, i.e., solving

[TABLE]

where $\Delta t$ is the time step size. Then, the new membrane position ${\mathbf{x}}^{\Delta t}\in\gamma^{\Delta t}$ is given by solving the relaxation problem

[TABLE]

The computationally expensive steps in this scheme are constructing or applying $V[\gamma]$ and solving the relaxation problem (3), which requires the expensive computation of $\mathbf{u}_{\mathrm{self}}[\gamma]$. Both operators depend nonlinearly on $\gamma$. In MLARM these steps are replaced with function approximations via MLPs as follows. First, for $V[\gamma]$, we need to remove the dependence on $\mathbf{u}_{\infty}$ in order to enable generalization to unseen $\mathbf{u}_{\infty}$ fields. We exploit the fact that $V[\gamma]\mathbf{u}_{\infty}$ depends linearly on $\mathbf{u}_{\infty}$. We decompose $\mathbf{u}_{\infty}$ using an $N_{f}$-term truncated Fourier series, $\mathbf{u}_{\infty}({\mathbf{x}}(\alpha))=\sum_{k=1}^{N_{f}}\phi_{k}(\alpha)\widehat{\mathbf{u}}_{\infty,k}$, where $\phi_{k}$ are Fourier basis vectors and $\widehat{\mathbf{u}}_{\infty,k}$ are the corresponding Fourier coefficients. Then, the term $V[\gamma]\mathbf{u}_{\infty}({\mathbf{x}})$ in (2) becomes

[TABLE]

We approximate $\{\Psi_{k}[\gamma]\}_{k=1}^{N_{f}}$ using $N_{f}$ MLPs, which we term $V$-MLPs. Then, given an unseen $\gamma$ and $\mathbf{u}_{\infty}$, we first compute the $\Psi_{k}$ operators using the $V$-MLPs and then apply the inverse Fourier transform. As a bonus, $\Psi_{k}[\gamma]$ turns out to be linear on certain parameters such as the membrane’s bending rigidity and the time step size, thus, there is no need to retrain the MLPs for these parameters. For the relaxation problem (3) we train another MLP that approximates the nonlinear function $\mathcal{R}(\gamma^{0})=\gamma^{\Delta t}$. We call this one $R$-MLP. The function $\mathcal{R}$ is nonlinear in the problem parameters, therefore, retraining the MLP for different values of bending rigidity and time step size is needed. Choosing a time step size depends on viscous and bending forces on a vesicle. We train several $R$-MLPs for different values of time step size and bending rigidity. Given flow parameters we determine the appropriate time step size and choose the corresponding $R$-MLP. One can also build a parametric reduced model using the trained MLPs Amsallem and Farhat (2008).

For multiple vesicles in confined flows (Fig. 1(b)), only the advection problem (2) changes. For the $i^{th}$ vesicle it becomes

[TABLE]

where $\mathbf{u}_{\infty}$ is the velocity on $\gamma_{i}$ due to its interactions with the other membranes $\gamma$ and the fixed boundaries $\Gamma$. Computing $\mathbf{u}_{\infty}$ requires solving additional integral equations: one for the density on the fixed boundaries and one for the inextensibility constraint. Boundaries do not evolve with time so the related matrices are precomputed and their application can be accelerated with fast multipole methods (FMM) Greengard and Rokhlin (1987). The inextensibility constraint requires solving an equation for Lagrange multiplier (tension). While this equation can be eliminated for a single vesicle, that is not the case for multiple vesicles and confined geometries. The equation resembles (1) and we solve it using $N_{f}+1$ MLPs. Once the tension for all vesicles and the density on the boundaries are obtained, we use FMM to compute $\mathbf{u}_{\infty}[\gamma,\Gamma]$. Then we compute the term $V[\gamma_{i}]\mathbf{u}_{\infty}[\gamma,\Gamma]$ with the approximated $\Psi_{k}[\gamma]\widehat{\mathbf{u}}_{\infty,k}$ using $N_{f}$ MLPs. Lastly, the relaxation problem (3) is solved. Evaluating layer potentials at points close to a boundary (fixed or vesicle) requires special quadratures Pozrikidis (1999); Kropinski (1999); Helsing and Ojala (2008); Zhao et al. (2010); Quaife and Biros (2014). These methods can be quite expensive to be implemented in MLARM. So, we correct near interactions between boundaries by employing a kinematic collision handling that is also used in simulations of emulsions Zinchenko and Davis (2008) and red blood cells Zhao et al. (2010) in 3D. For the minimum arclength spacing $h$ as a threshold, whenever the distance of a membrane point ${\mathbf{x}}$ to another membrane is less than $h$, the point is moved in the $({\mathbf{x}}_{p}-{\mathbf{x}})$ direction, where ${\mathbf{x}}_{p}$ is the projection of ${\mathbf{x}}$ on the other membrane until $\|{\mathbf{x}}-{\mathbf{x}}_{p}\|=h/2$. To have smooth perturbations in vesicle membranes, we also move the neighboring points in the same direction but for smaller amounts. See the Supplemental Material for all the details on the formulation and MLARM.

Overall, the MLPs approximate the computationally expensive terms that depend only on a vesicle shape $\gamma$. Our judicious application of MLPs allows us to train them only for various vesicle shapes, which eases data generation and use the trained MLPs for any flow configuration.

*Training MLPs.—*For the MLPs we use a lower dimensional representation of a vesicle shape based on principal component analysis (PCA). Given a training set of vesicle shapes, we find a reduced modal basis for the shapes using PCA and use it to represent the input and output to the MLPs. Approximately 100,000 distinct shapes make up the library. The outputs are high-fidelity solutions of the related integral equations. For $R$-MLP we solve the relaxation problem (3) for only one time step to obtain $\gamma^{\Delta t}$. The terms approximated by the MLPs linearly depend on translation, rotation and scaling of a vesicle. So, we standardize vesicle shapes so that they have the same center, angular orientation and arclength. The MLPs consist of five to six fully-connected layers with approximately 75,000 to 85,000 parameters. See the Supplemental Material for details on the training data and the MLP’s architectures.

*Numerical experiments.—*We demonstrate MLARM’s capabilities on several benchmark vesicle flows. We measure length in units of the vesicle radius $r_{0}$ (defined as the radius of a circle having the same enclosed area). The flows are characterized by capillary number $\mathrm{Ca}=\mu r_{0}^{3}\dot{\gamma}/\kappa_{b}$ where $\mu$ is the dynamic fluid viscosity and $\dot{\gamma}$ is the imposed shear rate that is varied to adjust Ca. True solutions are obtained using high-fidelity simulations Quaife and Biros (2014). We compare the MLARM simulations and the same-degree-of-freedom (DOF) simulations performed with the numerical scheme that has the same number of DOF as MLARM in terms of accuracy and computation time. We also compare with the same-cost simulations that have the same computation time with the MLARM simulations in terms of accuracy.

First, we consider a vesicle initialized at $y=r_{0}/2$ in a parabolic flow ($\mathbf{u}_{\infty}=r_{0}\dot{\gamma}\left(1-(y/W)^{2}\right){\mathbf{e}}_{x}$) with $W=10r_{0}$. In this flow a vesicle migrates towards the low shear rate region and reaches an equilibrium shape depending on its reduced area and Ca. We reproduce the results for a vesicle of reduced area 0.65 in Kaoui et al. (2009). The equilibrium shape is asymmetric (slipper-like) for all Ca values. Figure 2 shows that the MLARM simulations capture the true equilibrium shapes for all $\mathrm{Ca}<9$ more accurately than the same-DOF simulations do. For $\mathrm{Ca}=9$, the scheme requires finer temporal resolution. For this case, the MLARM simulations are slightly faster than the same-DOF ones. See the Supplemental Material for the vesicles’ equilibrium lateral positions and the comparison with the same-cost simulations. The MLARM simulations are more accurate in capturing the equilibrium lateral positions than the same-cost ones while they have similar accuracy in the equilibrium shapes.

Second, we study cross-streamline migration of a vesicle suspended in a flow with curved flow lines. The setup is the same as in Ghigliotti et al. (2011). The imposed shear rate is $\dot{\gamma}=20/r^{2}$ and the velocity field is $(v_{\theta},v_{r})=(-10/r,0)$ where $r$ is the distance between the vesicle’s center and the origin. The value of Ca depends on $r$ (varying from $\mathrm{Ca}=0.2$ at $r=10r_{0}$ to $\mathrm{Ca}\approx 2.2$ at $r=3r_{0}$). In this flow a vesicle with properties we consider (no viscosity contrast and reduced area of 0.65) migrates towards regions of high shear rates (Fig. 3(a)). We measure migration velocity in the radial direction from $r=10r_{0}$ to $r=3r_{0}$. Figure 3(b) shows that the MLARM simulation captures the migration with the error $|v_{\mathrm{mig}}^{\mathrm{MLA}}-v_{\mathrm{mig}}^{\mathrm{True}}|/v_{\mathrm{mig}}^{\mathrm{True}}=0.06$. We do not present the same-DOF and the same-cost simulations as they are completely wrong.

Third, we consider dilute suspensions in Taylor-Couette flow. Here the inner circle rotates in the counterclockwise direction while the outer one is stationary. The distance between the inner and the outer circles is $10r_{0}$ and Ca = 1.5. In such flows, the force applied by the inner circle balances vesicles’ inward migration. The vesicles eventually organize themselves in a rim with the same interdistance for area fractions of vesicles between approximately 1% and 4% Ghigliotti et al. (2011). We perform simulations with area fractions 2.2% and 4.4% until the inner circle completes 24 rotations. Figure 4 shows that closely and randomly initialized vesicles reach to approximately the same radial positions and are separated uniformly in the azimuthal direction at the equilibrium. The MLARM simulations can capture the vesicles’ spatial order although the individual trajectories are inaccurate. The MLARM simulation becomes less accurate than the same-DOF simulation as area fraction increases since MLARM ignores the near vesicle interactions (see the Supplemental Material for the evolution of the mean radial positions). Note that resolving individual trajectories is not of interest in suspensions because of their sensitivity to (generally unknown) initial conditions.

Finally, we consider dense suspensions in a Taylor-Couette flow. The setup is the same as in the previous example. Here vesicles migrate away from the circles and form regions near circles that are free of vesicles (called cell-free layers). We investigate how accurately MLARM can capture these layers. We plot the probability distribution of the distances of the vesicles’ centers to the origin throughout the simulations in Fig. 5 for an area fraction 20%. The MLARM and the same-DOF simulations accurately capture the layers and the same-cost simulation is less accurate. See the Supplemental Material for the results of area fractions 30% and 35% and also the movies. There is no resolution that provides stable and the same computational cost as the MLARM simulation in the 30% and 35% area fraction cases. Here the MLARM simulations are 10$\times$ faster than the same-DOF simulations and 25$\times$ faster than the ground truth simulations. This example clearly shows the effectiveness of the MLARM simulations compared to using just low-resolution simulations.

*Conclusion.—*We propose a reduced model that combines multilayer perceptrons trained with high-fidelity simulations for one time step and a low-fidelity numerical scheme. In various examples of vesicle flows, we compare the MLARM simulations with those given by the numerical scheme with the same degree-of-freedom and also that with the same computational cost. MLARM is faster and also provides simulations that have either similar or better accuracy. One of the limitations of MLARM is that we have not optimized the MLPs’ architectures and training parameters. Such optimization can provide more accurate approximations. Although we presented our analysis in 2D flows (so that we can ensure highly accurate ground truth simulations), the model generalizes as is to 3D vesicle flows and other particles like deformable capsules, drops, filaments and rigid bodies.

I Supplemental Material

I.1 Formulation

In this section, we present the mathematical modeling and the boundary integral equation formulation of vesicle flows in two dimensions.

I.1.1 Governing equations

We assume that there are no external forces and the same Newtonian fluid with viscosity $\mu$ occupies the interior and the exterior of vesicles. $K$ and $M$ are the numbers of fixed boundaries (solid walls) and vesicles, respectively. We denote the membrane (boundary) of the $i^{th}$ vesicle with $\gamma_{i}$. Let $\Gamma_{1},\cdots,\Gamma_{K-1}$ be walls surrounded by another wall $\Gamma_{0}$ (Fig. 1(b)). Since we consider flows in the limit of vanishing Reynolds number, the Stokes equations govern the fluid flow:

[TABLE]

where $\Omega$ is $K$-ly connected domain of interest, ${\mathbf{u}}$ is the fluid velocity field and $p$ is the pressure. We impose the no-slip boundary condition on membranes and the velocity Dirichlet boundary condition on the walls as

[TABLE]

Vesicle membranes are locally inextensible, i.e., each vesicle has a constant arclength. So, (5) is subject to the following constraint

[TABLE]

The subscript $s$ stands for differentiation with respect to the arclength. Finally, the momentum balance on a membrane requires the jump in surface traction to be equal to the membrane force due to bending and tension, i.e.,

[TABLE]

where $\mathbf{T}=-p\mathbf{I}+\mu(\nabla{\mathbf{u}}+\nabla{\mathbf{u}}^{T})$ is the Cauchy stress tensor, ${\mathbf{n}}$ is the outward normal vector on $\gamma_{i}$, $[\![\cdot]\!]$ is the jump across the membrane. The first term on the right hand side is the membrane force due to the bending rigidity $\kappa_{b}$ and the second one is the force due to the tension $\sigma$.

I.1.2 Boundary integral equation formulation

Let us first present the integral equation formulation of (5)-(8) for a single vesicle in a free-space flow, then, extend it to multiple vesicles in confined flows. The membrane force ${\mathbf{f}}=-\kappa_{b}{\mathbf{x}}_{ssss}+(\sigma{\mathbf{x}}_{s})_{s}$ induces velocity. The single layer integral $\mathcal{S}[{\mathbf{f}}]({\mathbf{x}})$ evaluates the velocity at ${\mathbf{x}}$ as

[TABLE]

where the free-space 2D Green’s function for the Stokes equations is

[TABLE]

For convenience, we introduce bending operator $\mathcal{B}$, tension operator $\mathcal{T}$, surface divergence operator $\mathcal{P}$ and stretching operator $\mathcal{M}$ defined for ${\mathbf{x}}\in\gamma$ as follows

[TABLE]

The velocity of the membrane points ${\mathbf{x}}\in\gamma$ in a free-space flow with the velocity field $\mathbf{u}_{\infty}$ is given as

[TABLE]

So, the terms $\mathbf{u}_{\mathrm{self}}[\gamma]({\mathbf{x}})$ and $V[\gamma]$ in the letter (Eq. 1) correspond to $(1-\mathcal{M}({\mathbf{x}}))\mathcal{B}({\mathbf{x}}){\mathbf{x}}$ and $(1-\mathcal{M}({\mathbf{x}}))$, respectively. Here the tension $\sigma$ is eliminated as the stretching operator modifies the membrane velocity to enforce the inextensibility constraint (see Veerapaneni et al. (2009) for its derivation).

For multiple vesicles in confined flows, the velocity of the points ${\mathbf{x}}$ on the $i^{th}$ membrane $\gamma_{i}$ is given by the same equation as (11) with a different composition of the term $\mathbf{u}_{\infty}$, i.e.,

[TABLE]

where $\gamma=\cup_{j=1,j\neq i}\gamma_{j}$ and $\Gamma=\cup_{k=0}\Gamma_{k}$. The term $\mathbf{u}_{\infty}$ is velocity induced by other vesicles in the flow and the walls on the $i^{th}$ membrane instead of the velocity field of the background flow as in (11). The velocity $\mathbf{u}_{\infty}[\gamma,\Gamma]({\mathbf{x}})$ is given as

[TABLE]

where $\mathcal{W}[\eta]$ is the velocity induced by the walls due to the density $\eta$ on them. Computing (13) requires solving integral equations for the tension $\sigma$ and the density $\eta$. The inextensibility constraint (7) delivers the following equation for the tension

[TABLE]

The Dirichlet velocity boundary condition on the walls provides the following equation for the density $\eta$ on $\Gamma$

[TABLE]

See Rahimian et al. (2010) for the definition of the completed double layer integral $\mathcal{W}[\eta]$. We refer the reader to Rahimian et al. (2010) for the semi-implicit time stepping scheme and Quaife and Biros (2014) for the implicit time stepping scheme to solve (11), (14) and (15). We obtain the true solutions in the present study using the implicit time stepping scheme while our machine-learning-augmented reduced model is built based on the semi-implicit one.

I.2 Machine-learning-augmented reduced model

In this section, we present the details of the machine-learning-augmented reduced model (MLARM). Specifically, we first explain the use of MLPs to obtain tension by solving (14) for multiple vesicles in confined flows. Then, we elaborate on generating a data set, dimensionality reduction using principal component analysis (PCA), multilayer perceptrons’ (MLP) architectures and parameters for training the MLPs. Finally, we present the pseudo-algorithm of MLARM.

I.2.1 Using MLPs to obtain tension for flows of multiple vesicles and solid walls

The way we enforce the inextensibility constraint requires solving (14) to obtain tension for each vesicle at every time step. The solution for the $i^{th}$ membrane can be expressed as

[TABLE]

where $\sigma_{\mathrm{self}}[\gamma_{i}]({\mathbf{x}})=-(\mathcal{P}({\mathbf{x}})\mathcal{T}({\mathbf{x}}))^{-1}\mathcal{P}({\mathbf{x}})\mathcal{B}({\mathbf{x}}){\mathbf{x}}$ and $T[\gamma_{i}]=-(\mathcal{P}({\mathbf{x}})\mathcal{T}({\mathbf{x}}))^{-1}$. So, the first term in (16) is the tension due to a vesicle itself and the second one is the contribution due to the other vesicles and the walls. The form of (16) resembles (1) for the velocity of the membrane points. That is why in MLARM we approximate $\sigma$ with a strategy similar to the one we use for (1). That is, we use a MLP to approximate the mapping between $\gamma_{i}$ and $\sigma_{\mathrm{self}}[\gamma_{i}]({\mathbf{x}})$, and use $N_{f}$ MLPs for the mappings between $\gamma_{i}$ and $\{T[\gamma_{i}]\phi_{k}(s)\}_{k}^{N_{f}}$ where $\phi_{k}(s)$ is the Fourier basis vectors. We call the former $\sigma$-MLP and the latter $T$-MLPs. The input to the MLPs is the PCA coefficients of a vesicle shape $\gamma_{i}$. Since the same PCA basis cannot be used to represent $\sigma_{\mathrm{self}}[\gamma_{i}]({\mathbf{x}})$, the output of $\sigma$-MLP is the Fourier coefficients corresponding to the first 32 modes which accurately reconstruct the term.

I.2.2 Generating a data set

Vesicles show wide variety of shapes in long time horizon simulations of their dense suspensions. In order to have various physical shapes we perform such a simulation of vesicles having the reduced area 0.65 in a confined Taylor-Couette flow (see Kabacaoğlu et al. (2018) for the snapshots of the simulation). We want a library that consists of distinct vesicle shapes by some metric. In order to compare the shapes obtained from the simulation, we, first, standardize the shapes such that their centers, inclination angles111The inclination angle is the angle between the flow direction and the principal axis corresponding to the smallest principal moment of the inertia. The moment of inertia tensor is $J=\frac{1}{4}\int_{\gamma}({\mathbf{r}}\cdot{\mathbf{n}})(|{\mathbf{r}}^{2}|I-{\mathbf{r}}\otimes{\mathbf{r}})\,d\gamma({\mathbf{x}})$ where ${\mathbf{r}}={\mathbf{x}}-{\mathbf{c}}$ and ${\mathbf{c}}$ is the center of the vesicle. and arclengths are the same. Then, we measure the dissimilarity between the standardized shapes based on the Hausdorff metric. Our library consists of $M=100,181$ distinct vesicle shapes. See Fig. 6 for 100 randomly sampled shapes from the library.

Once we have the library, we generate sets of inputs and outputs for the MLPs using a high resolution discretization, i.e., $N=96$ points per vesicle using Algorithm 1. As the approximated terms linearly depend on translation, rotation and scaling of a vesicle, we standardize the input vesicle shapes ${\mathbf{x}}$ so that they have the same center and inclination angle. We also use a standard ordering of the discretization points on vesicle membranes. The first point is the one on the positive $x$-axis, then the other points are equally distributed along the arclength in the counterclockwise direction. The output of $R$-MLP is the solution of (11) with $\mathbf{u}_{\infty}=0$, the bending rigidity $\kappa_{b}=1$ and the time step size $\Delta t=10^{-4}$ only for one time step. This time step size is the largest one with which stable and accurate simulations in a stationary fluid can be performed using the semi-implicit time stepping scheme. For the other MLPs the outputs are the high-fidelity solutions of the related integral equations for vesicles in the library. We have found that $N_{f}=24$ Fourier modes accurately represent the velocity $\mathbf{u}_{\infty}$. Therefore, we use $N_{f}=24$ MLPs to approximate the action of each $V$ and $T$ on the Fourier basis vectors ($V$-MLPs and $T$-MLPs). Each coefficient $k$ has two components, imaginary and real. Additionally, since we represent the vesicles with $N=96$ and each point has two degrees of freedom, $V\phi_{k}$ and $T\phi_{k}$ are of sizes 192 and 96. We reduce this size by subsampling $V\phi_{k}$ to $N=24$ points and $T\phi_{k}$ to $N=48$. The outputs to the $k^{th}$ MLPs for these terms are vectors of size 96 that contain real and imaginary components of the terms. By choosing the subsampling rates as such we can use the same MLPs for both terms.

I.3 Principal component analysis

We perform dimensionality reduction using principal component analysis. Here we describe how to form a reduced basis and represent a vesicle shape using this basis. A column vector ${\mathbf{x}}=[x_{1}\,\,x_{2}\,\,\cdots\,\,x_{N}\,\,y_{1}\,\,y_{2}\,\,\cdots\,\,y_{N}]^{T}$ of size $2N$ contains the $x$ and $y$ coordinates of the membrane points and represents a vesicle shape $\gamma$. Let $\mathbf{X}=[{\mathbf{x}}_{1}\,\,{\mathbf{x}}_{2}\,\,\cdots\,\,{\mathbf{x}}_{M}]$ be a matrix of size $2N$-by-$M$ that stores $M$ vesicle shapes in the library. The singular value decomposition of $\mathbf{X}$ is

[TABLE]

where ${\bm{W}}$ is a $2N$-by-$2N$ matrix, ${\bm{\Sigma}}$ is a $2N$-by-$M$ rectangular diagonal matrix of the singular values of $\mathbf{X}$ and ${\bm{V}}$ is an $M$-by-$M$ matrix. Let ${\mathbf{C}}={\bm{\Sigma}}{\bm{V}}^{T}$ be a $2N$-by-$M$ score matrix whose columns contain the representations of $M$ samples in the orthogonally transformed space. We obtain a truncated ${\mathbf{C}}_{N_{p}}$ matrix of size $N_{p}$-by-$M$ by considering only the first $N_{p}$ largest singular values and the corresponding singular vectors,

[TABLE]

where ${\bm{W}}_{N_{p}}=[w_{1}\,\,w_{2}\,\,\cdots\,\,w_{N_{p}}]$ is an $2N$-by-$N_{p}$ matrix whose columns are the reduced basis vectors. Then the lower dimensional representation of a vesicle ${\mathbf{x}}$ is $c={\bm{W}}^{T}_{N_{p}}{\mathbf{x}}$. We set $N_{p}=16$ as the first 16 PCA modes capture 99% of the total energy, i.e., $(\sum_{i=1}^{16}\sigma_{i}^{2})/(\Sigma_{i=1}^{2N}\sigma_{i}^{2})>99\%$.

I.3.1 MLP architectures

A detailed introduction of multilayer perceptrons is beyond the scope. We only introduce definitions that are helpful to our presentation of the MLPs’ architectures (see Goodfellow et al. (2016) for more information). A MLP has a number of parameters that are found by minimizing a loss function. Let $c$ and $c^{\Delta t}$ be the lower dimensional representations of vesicles ${\mathbf{x}}$ and ${\mathbf{x}}^{\Delta t}$. Then, we find the parameters of $R$-MLP by minimizing the mean squared error

[TABLE]

where $\widetilde{\mathbf{x}}^{\Delta t}_{j}=W_{N_{p}}\widetilde{c}^{\Delta t}_{j}$ correspond to the MLP’s approximation to the $j^{th}$ vesicle in the data set of size $M$. $\|\cdot\|_{2}$ is the Euclidean norm. One can also measure the mean squared error based on the PCA coefficients $c$. We have also used a MLP trained to minimize such an error and obtained similar results. The parameters of the $k^{th}$ $V$-MLP ($T$-MLP as well) minimize the mean squared error

[TABLE]

The parameters of $\sigma$-MLP minimize the mean squared error between the true tension and the one reconstructed using the approximated Fourier coefficients.

All the MLPs we use consist of fully connected layers that are defined as follows. Let $W_{k}$ and $b_{k}$ be the $k^{th}$ fully connected layer’s weight matrix of size $n_{k+1}$-by-$n_{k}$ and bias vector of length $n_{k}$. Then the layer transforms an input vector $z_{k}$ of length $n_{k}$ into a vector $W_{k}z_{k}+b_{k}$ of length $n_{k+1}$. In the MLPs we use batch normalization after every fully connected layer except the last one to shift $W_{k}z_{k}+b_{k}$ to zero mean and scale it to unit variance over the data set. Finally, a nonlinear function $f$ computes an output $z_{k+1}=f(W_{k}z_{k}+b_{k})$. Our choice for $f$ is the so-called leaky rectified linear unit (leaky ReLu), i.e., $f(z)=\beta z$ if $z<0$, $f(z)=z$ if $z\geq 0$, where $\beta$ is a small constant. We set $\beta=0.1$. The function $f$ operates on $W_{k}z_{k}+b_{k}$ elementwise. A MLP consists of a cascade of such transformations. A $K$-layer MLP maps an input $z_{1}$ to an output $z_{K+1}$ as follows

[TABLE]

The length of the vector $z_{k+1}$ is called the width of the $k^{th}$ layer. A MLP’s parameters are the entries of the weight matrices $W_{k}$ and the bias vectors $b_{k}$. We use six fully connected layers in $R$-MLP. The input vector is of length 16 and the widths of the layers are 48, 96, 128, 256, 128 and 16. Therefore, the MLP has 85,920 parameters. $V$-MLPs and $T$-MLPs have the same architecture which consists of five fully connected layers. The widths of the layers are 48, 96, 128, 256 and 96. The MLP has 75,632 parameters. $\sigma$-MLP has five fully connected layers. The widths of the layers are 64, 128, 256, 128 and 64 which make up 83,584 parameters. We normalize the inputs such that they have a unit mean and zero standard deviation.

Finding parameters minimizing the loss functions requires solving a nonconvex optimization problem. We solve this problem using a stochastic gradient descent method (SGDM). SGDM computes the gradients of a loss function with respect to the parameters of a MLP over a randomly sampled subset of the data set which is called a minibatch. In the trainings we use minibatches of size 256. SGDM iterations are terminated either when the error stagnates or SGDM exhausts the entire data set 20 times (i.e., after 20 epochs). SGDM has a parameter called learning rate that scales the size of a step taken in the descent direction. We set the initial learning rate to $5\times 10^{-4}$ and scaled it by a factor of $0.2$ after ten epochs.

I.3.2 MLARM: Pseudo-algorithm

In Algorithm 5 we present the pseudo-algorithm of MLARM only for a vesicle in a free-space flow case due to its simplicity. As mentioned in the letter, we use several correction algorithms in order to improve accuracy of the solution and also stability of MLARM. One of them is the area and length correction algorithm ($\mathtt{correctShape}$). After correcting the errors in a vesicle’s area and arclength, the vesicle might rotate and translate. We remove such artificial rotation and translation using $\mathtt{alignShape}$ algorithm. See Kabacaoğlu et al. (2018) for the details of these algorithms. In the flows with multiple vesicles and solid walls, we compute the pairwise interactions between the boundaries using $N_{LR}$ points to discretize a vesicle membrane. So, the pairwise interactions are computed at a low resolution.

I.4 Additional numerical results

In this section, we explain the setup for the numerical experiments presented in the letter. Then we present results for a single vesicle in a free-space shear flow, the same-cost simulations of a single vesicle in a parabolic flow and discuss how the simulations capture the equilibrium lateral position of a vesicle in this flow. We also present the statistics of vesicles in dilute and dense Taylor-Couette flows and compare the MLARM simulations with the same-cost and the same-DOF ones.

I.4.1 Setup for numerical experiments

The true solutions are obtained using $N=96$ points per vesicle and the time step size $\Delta t=10^{-4}$. In the confined flow examples, we use 192 points to discretize the walls in the MLARM, the same-DOF, and the same-cost simulations and 256 points in the true solution. In the single vesicle examples (parabolic and shear flows, flow with curved flow lines), we use 16 PCA modes in MLARM and $N_{LR}=16$ to solve the equations exactly at every ten time steps. In the confined flow examples, we use 32 PCA modes in order to have more degrees-of-freedom in vesicle shapes and $N_{LR}=32$ to compute the pairwise interactions via FMM. In these examples, we do not solve the equations exactly at all. Accordingly, the number of points per vesicle is $N=16$ in the same-DOF simulations of the single vesicle examples and $N=32$ in those of the confined flow ones. The same-DOF simulations have $\Delta t=10^{-4}$. We determine $N$ and $\Delta t$ for the same-cost simulations by taking several time steps with the numerical scheme for different $N$ and $\Delta t$ and matching the computation time. We set $N=48$ and $\Delta t=3\times 10^{-4}$ in the same-cost simulations as they have the same CPU time as the MLARM simulations for this resolution.

I.4.2 Free-space shear flow

We choose an initial vesicle shape and simulate it in shear flows ($\mathbf{u}_{\infty}=\dot{\gamma}y\mathbf{e}_{x}$) for different Ca values. Since the fluids in the interior and the exterior of a vesicle have the same viscosity, the vesicle tilts to a certain angle and its membrane tank-treads Kraus et al. (1996). Figure 7 shows the equilibrium shapes of the vesicle and compares the MLARM simulations with the same-DOF and the same-cost ones. Both low-resolution simulations and the MLARM simulations are similarly accurate in the equilibrium vesicle shapes. In order to investigate the vesicle’s tank-treading we also present the vesicle’s phase angle as a function of time in Fig. 8, i.e., the angle between its principal axis and the vector from the its center to a specific point on the membrane. The figure shows that although the MLARM simulation captures tank-treading, the vesicle tank-treads faster than the true solution does. Additionally, the same-cost simulation cannot accurately capture the tank-treading velocity either. The error stems from the numerical error introduced by the greater time step size in the same-cost simulation (which is $\Delta t=3\times 10^{-4}$) compared to the same-DOF and the MLARM simulations (which is $\Delta t=10^{-4}$).

I.4.3 Free-space parabolic flow

Figure 2 compares the same-DOF and the MLARM simulations in terms of capturing the equilibrium shapes of a vesicle in parabolic flows with different shear rates. Here we present the equilibrium shapes given by the same-cost simulations as well (Fig. 9). They are more accurate than the same-DOF ones in capturing the equilibrium shapes. However, they are not superior to the MLARM simulations.

We also compare the equilibrium lateral positions of vesicles in these flows (Fig. 9). The MLARM simulations capture the equilibrium lateral positions more accurately than the same-cost ones, especially for high Ca. We omit the same-DOF simulations as they have errors greater than 100%.

I.4.4 Dilute Taylor-Couette flows

Figure 4 shows that vesicles in a dilute Taylor-Couette flow exhibit a spatial order by rotating in a rim with a uniform angular interdistance. We present the evolution of the mean radial position of the vesicles (i.e., the mean of the distances of vesicles’ centers to the inner circle’s center) for area fractions 2.2% and 4.4% in Fig. 10. The figures show that the mean radial position converges as expected. The MLARM and the same-DOF simulations capture the true mean radial position accurately in the 2.2% case. The MLARM simulation has a greater error in the 4.4% case compared to the 2.2% case. The reason is that MLARM ignores near interactions between vesicles and these interactions start dominating the dynamics as area fraction increases. For the denser case, the same-DOF simulation becomes more accurate than the MLARM simulation.

In comparison with the same-cost simulation, the MLARM simulation is more accurate for an area fraction 2.2%. The reason is that the same-cost simulation has a coarser temporal resolution than the MLARM simulation and capturing vesicle migration accurately requires fine temporal resolution. For an area fraction 4.4%, however, the same-cost and the MLARM simulations have similar accuracy.

I.4.5 Dense Taylor-Couette flows

Figure 5 presents the probability distribution of the distance of a vesicles’ centers to the inner circle’s center throughout the simulations at an area fraction 20%. Here we present the results for area fractions 30% and 35% (Fig. 11). The MLARM and the same-DOF simulations accurately capture the cell-free layers whereas there is no resolution that provides stable and the same computational simulations for 30% and 35% area fractions.