A fluctuating boundary integral method for Brownian suspensions - - PowerPoint PPT Presentation

A fluctuating boundary integral method for Brownian suspensions

Yuanxun (Bill) Bao Courant Institute, NYU Collaborators: Aleksandar Donev (Courant) Leslie Greengard (Courant) Eric Keaveny (Imperial) Manas Rachh (Yale) SIAM Computational Science and Engineering March 2, 2017

• Y. Bao (CIMS)

Fluctuating BIE

Brownian Dynamics with Hydrodynamic Interactions

⊲ Consider a suspension of Nb rigid bodies with configuration Q = {qβ, θβ}Nb

β=1 consisting of positions and orientations im-

mersed in a Stokes fluid. ⊲ The Ito stochastic equation of Brownian Dynamics (BD) is dQ dt = −N ∂QU + (2kBTN )

1 2 W(t) + (kBT) ∂Q · N ,

where N (Q) is the body mobility matrix, U(Q) is the potential energy, kBT is the temperature, and W(t) is a vector of independent white noise processes. ⊲ Here the stochastic noise amplitude is determined from the fluctuation-dissipation balance: N

1 2

• N

1 2

∗ = N . ⊲ The stochastic drift term ∂Q · N =

j ∂jNij is related to the Ito

interpretation of the noise.

• Y. Bao (CIMS)

Fluctuating BIE

Brownian Dynamics with Hydrodynamic Interactions

⊲ Consider a suspension of Nb rigid bodies with configuration Q = {qβ, θβ}Nb

β=1 consisting of positions and orientations im-

mersed in a Stokes fluid. ⊲ The Ito stochastic equation of Brownian Dynamics (BD) is dQ dt = −N ∂QU + (2kBTN )

1 2 W(t) + (kBT) ∂Q · N ,

where N (Q) is the body mobility matrix, U(Q) is the potential energy, kBT is the temperature, and W(t) is a vector of independent white noise processes. ⊲ Here the stochastic noise amplitude is determined from the fluctuation-dissipation balance: N

1 2

• N

1 2

∗ = N . ⊲ The stochastic drift term ∂Q · N =

j ∂jNij is related to the Ito

interpretation of the noise.

• Y. Bao (CIMS)

Fluctuating BIE

Brownian Dynamics with Hydrodynamic Interactions

⊲ Consider a suspension of Nb rigid bodies with configuration Q = {qβ, θβ}Nb

β=1 consisting of positions and orientations im-

mersed in a Stokes fluid. ⊲ The Ito stochastic equation of Brownian Dynamics (BD) is dQ dt = −N ∂QU + (2kBTN )

1 2 W(t) + (kBT) ∂Q · N ,

where N (Q) is the body mobility matrix, U(Q) is the potential energy, kBT is the temperature, and W(t) is a vector of independent white noise processes. ⊲ Here the stochastic noise amplitude is determined from the fluctuation-dissipation balance: N

1 2

• N

1 2

∗ = N . ⊲ The stochastic drift term ∂Q · N =

j ∂jNij is related to the Ito

interpretation of the noise.

• Y. Bao (CIMS)

Fluctuating BIE

Brownian Dynamics with Hydrodynamic Interactions

⊲ Consider a suspension of Nb rigid bodies with configuration Q = {qβ, θβ}Nb

β=1 consisting of positions and orientations im-

mersed in a Stokes fluid. ⊲ The Ito stochastic equation of Brownian Dynamics (BD) is dQ dt = −N ∂QU + (2kBTN )

1 2 W(t) + (kBT) ∂Q · N ,

where N (Q) is the body mobility matrix, U(Q) is the potential energy, kBT is the temperature, and W(t) is a vector of independent white noise processes. ⊲ Here the stochastic noise amplitude is determined from the fluctuation-dissipation balance: N

1 2

• N

1 2

∗ = N . ⊲ The stochastic drift term ∂Q · N =

j ∂jNij is related to the Ito

interpretation of the noise.

• Y. Bao (CIMS)

Fluctuating BIE

Hydrodynamic Body Mobility Matrix

⊲ The body mobility matrix N (Q) 0 is a symmetric positive semidefinite (SPD) and it includes hydrodynamic interactions and (periodic) boundary conditions. ⊲ For viscous-dominated flows (Re → 0), we can assume steady Stokes flow and solve the Stokes mobility problem, U = N F, where U = {uβ, ωβ}Nb

β=1 collects the linear and angular velocities,

F = {fβ, τ β}Nb

β=1 collects the applied forces and torques.

⊲ At every time step of BD simulation, we need to generate particle velocity in the form of (dropping kBT),

• U = N F + N

1 2 W.

⊲ This talk: How to accurately and efficiently compute the action of N and N

1 2 ?

• Y. Bao (CIMS)

Fluctuating BIE

Hydrodynamic Body Mobility Matrix

⊲ The body mobility matrix N (Q) 0 is a symmetric positive semidefinite (SPD) and it includes hydrodynamic interactions and (periodic) boundary conditions. ⊲ For viscous-dominated flows (Re → 0), we can assume steady Stokes flow and solve the Stokes mobility problem, U = N F, where U = {uβ, ωβ}Nb

β=1 collects the linear and angular velocities,

F = {fβ, τ β}Nb

β=1 collects the applied forces and torques.

⊲ At every time step of BD simulation, we need to generate particle velocity in the form of (dropping kBT),

• U = N F + N

1 2 W.

⊲ This talk: How to accurately and efficiently compute the action of N and N

1 2 ?

• Y. Bao (CIMS)

Fluctuating BIE

Hydrodynamic Body Mobility Matrix

⊲ The body mobility matrix N (Q) 0 is a symmetric positive semidefinite (SPD) and it includes hydrodynamic interactions and (periodic) boundary conditions. ⊲ For viscous-dominated flows (Re → 0), we can assume steady Stokes flow and solve the Stokes mobility problem, U = N F, where U = {uβ, ωβ}Nb

β=1 collects the linear and angular velocities,

F = {fβ, τ β}Nb

β=1 collects the applied forces and torques.

⊲ At every time step of BD simulation, we need to generate particle velocity in the form of (dropping kBT),

• U = N F + N

1 2 W.

⊲ This talk: How to accurately and efficiently compute the action of N and N

1 2 ?

• Y. Bao (CIMS)

Fluctuating BIE

Hydrodynamic Body Mobility Matrix

⊲ The body mobility matrix N (Q) 0 is a symmetric positive semidefinite (SPD) and it includes hydrodynamic interactions and (periodic) boundary conditions. ⊲ For viscous-dominated flows (Re → 0), we can assume steady Stokes flow and solve the Stokes mobility problem, U = N F, where U = {uβ, ωβ}Nb

β=1 collects the linear and angular velocities,

F = {fβ, τ β}Nb

β=1 collects the applied forces and torques.

⊲ At every time step of BD simulation, we need to generate particle velocity in the form of (dropping kBT),

• U = N F + N

1 2 W.

⊲ This talk: How to accurately and efficiently compute the action of N and N

1 2 ?

• Y. Bao (CIMS)

Fluctuating BIE

First-Kind Boundary Integral Formulation

⊲ Let us first ignore Brownian terms and solve a mobility problem to compute N F. ⊲ For simplicity, consider only a single body Ω. The first-kind bound- ary integral equation for the mobility problem, v(q) = u+ω ×q =

• ∂Ω

G(q−q′) µ(q′) dq′ for all q ∈ ∂Ω, (1) along with force and torque balance conditions

• ∂Ω

µ(q) dq = f and

• ∂Ω

q × µ(q) dq = τ, (2) where µ(q ∈ ∂Ω) is the surface traction (single-layer density) and G is the (periodic) Stokeslet. ⊲ Note that one can alternatively use a completed second-kind or a mixed first-second kind formulation for improved conditioning. We only know how to generate Brownian displacements efficiently in the first-kind formulation.

• Y. Bao (CIMS)

Fluctuating BIE

First-Kind Boundary Integral Formulation

⊲ Let us first ignore Brownian terms and solve a mobility problem to compute N F. ⊲ For simplicity, consider only a single body Ω. The first-kind bound- ary integral equation for the mobility problem, v(q) = u+ω ×q =

• ∂Ω

G(q−q′) µ(q′) dq′ for all q ∈ ∂Ω, (1) along with force and torque balance conditions

• ∂Ω

µ(q) dq = f and

• ∂Ω

q × µ(q) dq = τ, (2) where µ(q ∈ ∂Ω) is the surface traction (single-layer density) and G is the (periodic) Stokeslet. ⊲ Note that one can alternatively use a completed second-kind or a mixed first-second kind formulation for improved conditioning. We only know how to generate Brownian displacements efficiently in the first-kind formulation.

• Y. Bao (CIMS)

Fluctuating BIE

First-Kind Boundary Integral Formulation

⊲ Let us first ignore Brownian terms and solve a mobility problem to compute N F. ⊲ For simplicity, consider only a single body Ω. The first-kind bound- ary integral equation for the mobility problem, v(q) = u+ω ×q =

• ∂Ω

G(q−q′) µ(q′) dq′ for all q ∈ ∂Ω, (1) along with force and torque balance conditions

• ∂Ω

µ(q) dq = f and

• ∂Ω

q × µ(q) dq = τ, (2) where µ(q ∈ ∂Ω) is the surface traction (single-layer density) and G is the (periodic) Stokeslet. ⊲ Note that one can alternatively use a completed second-kind or a mixed first-second kind formulation for improved conditioning. We only know how to generate Brownian displacements efficiently in the first-kind formulation.

• Y. Bao (CIMS)

Fluctuating BIE

First-Kind Boundary Integral Formulation

⊲ Let us first ignore Brownian terms and solve a mobility problem to compute N F. ⊲ For simplicity, consider only a single body Ω. The first-kind bound- ary integral equation for the mobility problem, v(q) = u+ω ×q =

• ∂Ω

G(q−q′) µ(q′) dq′ for all q ∈ ∂Ω, (1) along with force and torque balance conditions

• ∂Ω

µ(q) dq = f and

• ∂Ω

q × µ(q) dq = τ, (2) where µ(q ∈ ∂Ω) is the surface traction (single-layer density) and G is the (periodic) Stokeslet. ⊲ Note that one can alternatively use a completed second-kind or a mixed first-second kind formulation for improved conditioning. We only know how to generate Brownian displacements efficiently in the first-kind formulation.

• Y. Bao (CIMS)

Fluctuating BIE

First-Kind Boundary Integral Formulation

⊲ Assume that the surface of the body is discretized in some man- ner, and the single-layer operator M is approximated by some quadrature,

• ∂Ω

G(q − q′) µ(q′) dq′ ≡ Mµ → Mλ, where M is a SPD operator with kernel G with r −1 singularity in 3D (log r in 2D), discretized as a SPD mobility matrix M. ⊲ In matrix notation the mobility problem can be written as a saddle- point linear system for the surface forces λ and rigid-body motion U,

• M

−K −K⊤ λ U

• = −
• F
• ,

(3) where K is a simple geometric matrix. ⊲ Using Schur complement to eliminate λ, we get U = N F = (K⊤M−1K)−1F.

• Y. Bao (CIMS)

Fluctuating BIE

First-Kind Boundary Integral Formulation

⊲ Assume that the surface of the body is discretized in some man- ner, and the single-layer operator M is approximated by some quadrature,

• ∂Ω

G(q − q′) µ(q′) dq′ ≡ Mµ → Mλ, where M is a SPD operator with kernel G with r −1 singularity in 3D (log r in 2D), discretized as a SPD mobility matrix M. ⊲ In matrix notation the mobility problem can be written as a saddle- point linear system for the surface forces λ and rigid-body motion U,

• M

−K −K⊤ λ U

• = −
• F
• ,

(3) where K is a simple geometric matrix. ⊲ Using Schur complement to eliminate λ, we get U = N F = (K⊤M−1K)−1F.

• Y. Bao (CIMS)

Fluctuating BIE

First-Kind Boundary Integral Formulation

⊲ Assume that the surface of the body is discretized in some man- ner, and the single-layer operator M is approximated by some quadrature,

• ∂Ω

G(q − q′) µ(q′) dq′ ≡ Mµ → Mλ, where M is a SPD operator with kernel G with r −1 singularity in 3D (log r in 2D), discretized as a SPD mobility matrix M. ⊲ In matrix notation the mobility problem can be written as a saddle- point linear system for the surface forces λ and rigid-body motion U,

• M

−K −K⊤ λ U

• = −
• F
• ,

(3) where K is a simple geometric matrix. ⊲ Using Schur complement to eliminate λ, we get U = N F = (K⊤M−1K)−1F.

• Y. Bao (CIMS)

Fluctuating BIE

Brownian Displacements

⊲ How do we compute the action of N

1 2 ?

More precisely, how to generate a Gaussian random vector with co- variance N ? ⊲ Assume for now we knew the action of M

1 2 on a vector of Gaussian

random variables W, i.e., that we knew how to generate a random “slip” velocity with covariance given by M (periodic Stokeslet). ⊲ Key idea 1: solve the mobility problem with such a random slip

• M

−K −K⊤ λ U

• = −
• M

1 2 W

F

• ,

(4) = ⇒ U = N F + N K⊤M−1M

1 2 W = N F + N 1 2 W,

which defines a N

1 2 with the correct covariance:

N

1 2

• N

1 2

∗ = N K⊤M−1M

1 2

• M

1 2

∗ M−1K N = N (K⊤M−1K)N = N (N )−1N = N . (5)

• Y. Bao (CIMS)

Fluctuating BIE

Brownian Displacements

⊲ How do we compute the action of N

1 2 ?

More precisely, how to generate a Gaussian random vector with co- variance N ? ⊲ Assume for now we knew the action of M

1 2 on a vector of Gaussian

random variables W, i.e., that we knew how to generate a random “slip” velocity with covariance given by M (periodic Stokeslet). ⊲ Key idea 1: solve the mobility problem with such a random slip

• M

−K −K⊤ λ U

• = −
• M

1 2 W

F

• ,

(4) = ⇒ U = N F + N K⊤M−1M

1 2 W = N F + N 1 2 W,

which defines a N

1 2 with the correct covariance:

N

1 2

• N

1 2

∗ = N K⊤M−1M

1 2

• M

1 2

∗ M−1K N = N (K⊤M−1K)N = N (N )−1N = N . (5)

• Y. Bao (CIMS)

Fluctuating BIE

Brownian Displacements

⊲ How do we compute the action of N

1 2 ?

More precisely, how to generate a Gaussian random vector with co- variance N ? ⊲ Assume for now we knew the action of M

1 2 on a vector of Gaussian

random variables W, i.e., that we knew how to generate a random “slip” velocity with covariance given by M (periodic Stokeslet). ⊲ Key idea 1: solve the mobility problem with such a random slip

• M

−K −K⊤ λ U

• = −
• M

1 2 W

F

• ,

(4) = ⇒ U = N F + N K⊤M−1M

1 2 W = N F + N 1 2 W,

which defines a N

1 2 with the correct covariance:

N

1 2

• N

1 2

∗ = N K⊤M−1M

1 2

• M

1 2

∗ M−1K N = N (K⊤M−1K)N = N (N )−1N = N . (5)

• Y. Bao (CIMS)

Fluctuating BIE

Brownian Displacements

⊲ How do we compute the action of N

1 2 ?

More precisely, how to generate a Gaussian random vector with co- variance N ? ⊲ Assume for now we knew the action of M

1 2 on a vector of Gaussian

random variables W, i.e., that we knew how to generate a random “slip” velocity with covariance given by M (periodic Stokeslet). ⊲ Key idea 1: solve the mobility problem with such a random slip

• M

−K −K⊤ λ U

• = −
• M

1 2 W

F

• ,

(4) = ⇒ U = N F + N K⊤M−1M

1 2 W = N F + N 1 2 W,

which defines a N

1 2 with the correct covariance:

N

1 2

• N

1 2

∗ = N K⊤M−1M

1 2

• M

1 2

∗ M−1K N = N (K⊤M−1K)N = N (N )−1N = N . (5)

• Y. Bao (CIMS)

Fluctuating BIE

The Single-Layer Mobility Matrix

⊲ We need accurate and efficient routines to compute the action of M and M

1 2 for many bodies.

⊲ Recall that Mµ ≡

• ∂Ω G(q − q′) µ(q′) dq′, where G is a weakly

singular kernel that includes periodic BC effects. ⊲ Key idea 2: Singular quadrature (Alpert in 2D) + Spectral Ewald method to split the Stokeslet into near-field and far-field pieces: G = G(r)

ξ

+ G(w)

ξ

. ⊲ This idea comes from the recent work “Rapid Sampling of Stochastic Displacements in Brownian Dynamics Simulations” by A. M. Fiore,

• F. Balboa Usabiaga, A. Donev and J. W. Swan, to appear in J.
• Chem. Phys., 2017 [?].
• Y. Bao (CIMS)

Fluctuating BIE

The Single-Layer Mobility Matrix

⊲ We need accurate and efficient routines to compute the action of M and M

1 2 for many bodies.

⊲ Recall that Mµ ≡

• ∂Ω G(q − q′) µ(q′) dq′, where G is a weakly

singular kernel that includes periodic BC effects. ⊲ Key idea 2: Singular quadrature (Alpert in 2D) + Spectral Ewald method to split the Stokeslet into near-field and far-field pieces: G = G(r)

ξ

+ G(w)

ξ

. ⊲ This idea comes from the recent work “Rapid Sampling of Stochastic Displacements in Brownian Dynamics Simulations” by A. M. Fiore,

• F. Balboa Usabiaga, A. Donev and J. W. Swan, to appear in J.
• Chem. Phys., 2017 [?].
• Y. Bao (CIMS)

Fluctuating BIE

The Single-Layer Mobility Matrix

⊲ We need accurate and efficient routines to compute the action of M and M

1 2 for many bodies.

⊲ Recall that Mµ ≡

• ∂Ω G(q − q′) µ(q′) dq′, where G is a weakly

singular kernel that includes periodic BC effects. ⊲ Key idea 2: Singular quadrature (Alpert in 2D) + Spectral Ewald method to split the Stokeslet into near-field and far-field pieces: G = G(r)

ξ

+ G(w)

ξ

. ⊲ This idea comes from the recent work “Rapid Sampling of Stochastic Displacements in Brownian Dynamics Simulations” by A. M. Fiore,

• F. Balboa Usabiaga, A. Donev and J. W. Swan, to appear in J.
• Chem. Phys., 2017 [?].
• Y. Bao (CIMS)

Fluctuating BIE

The Single-Layer Mobility Matrix

⊲ We need accurate and efficient routines to compute the action of M and M

1 2 for many bodies.

⊲ Recall that Mµ ≡

• ∂Ω G(q − q′) µ(q′) dq′, where G is a weakly

singular kernel that includes periodic BC effects. ⊲ Key idea 2: Singular quadrature (Alpert in 2D) + Spectral Ewald method to split the Stokeslet into near-field and far-field pieces: G = G(r)

ξ

+ G(w)

ξ

. ⊲ This idea comes from the recent work “Rapid Sampling of Stochastic Displacements in Brownian Dynamics Simulations” by A. M. Fiore,

• F. Balboa Usabiaga, A. Donev and J. W. Swan, to appear in J.
• Chem. Phys., 2017 [?].
• Y. Bao (CIMS)

Fluctuating BIE

Ewald Splitting of M

⊲ The splitting of G induces a corresponding splitting of M into near- field and far-field pieces: M = M(r) + M(w) = M(r)

Alpert + M(r) trap + M(w),

where M(r)

Alpert is a block-diagonal band-limited matrix whose ele-

ments are (local) Alpert corrections to the trapezoidal rule and

• M(r)

trap

• ij = G(r)

ξ (qi − qj)

and

• M(w)

ij = G(w) ξ

(qi − qj), i = j.

• Y. Bao (CIMS)

Fluctuating BIE

Ewald Splitting of M

⊲ The splitting of G induces a corresponding splitting of M into near- field and far-field pieces: M = M(r) + M(w) = M(r)

Alpert + M(r) trap + M(w),

where M(r)

Alpert is a block-diagonal band-limited matrix whose ele-

ments are (local) Alpert corrections to the trapezoidal rule and

• M(r)

trap

• ij = G(r)

ξ (qi − qj)

and

• M(w)

ij = G(w) ξ

(qi − qj), i = j.

• Y. Bao (CIMS)

Fluctuating BIE

Near-Field Piece of M

Mλ =

• M(r)

Alpert + M(r) trap + M(w)

λ ⊲ M(r)

Alpert: Alpert correction matrix is precomputed for a single body

in some reference configuration, and apply to each body via rotation. ⊲ M(r)

trap: sparse matrix-vector multiplication because the real-space

kernel G(r)

ξ

decays like e−ξ2d2, where d = |qi − qj|. The action of M(r)

trap can be efficiently computed by cell linked-lists

as used in the classical MD method.

• Y. Bao (CIMS)

Fluctuating BIE

Near-Field Piece of M

Mλ =

• M(r)

Alpert + M(r) trap + M(w)

λ ⊲ M(r)

Alpert: Alpert correction matrix is precomputed for a single body

in some reference configuration, and apply to each body via rotation. ⊲ M(r)

trap: sparse matrix-vector multiplication because the real-space

kernel G(r)

ξ

decays like e−ξ2d2, where d = |qi − qj|. The action of M(r)

trap can be efficiently computed by cell linked-lists

as used in the classical MD method.

• Y. Bao (CIMS)

Fluctuating BIE

Far-Field Piece of M

Mλ =

• M(r)

Alpert + M(r) trap + M(w)

λ ⊲ M(w): the wave-space kernel G(w)

ξ

is smooth and regular, G(w)(r) = 1 V

• k=0

eik·(r) H(k, ξ) k2 (I − ˆ kˆ k), (6) where the Hasimoto splitting function H(k, ξ) =

• 1 + k2

4ξ2

• e−k2/4ξ2.

⊲ We can efficiently compute the action of M(w) in Fourier space by using the Spectral Ewald method of Lindbo/Tornberg [?], M(w) = S†BS, (7) where S is the non-uniform FFT (NUFFT) of Greengard/Lee [?], and B is a SPD block-diagonal matrix (in Fourier space), B(k, ξ) = H(k, ξ) k2 (I − ˆ kˆ k).

• Y. Bao (CIMS)

Fluctuating BIE

Far-Field Piece of M

Mλ =

• M(r)

Alpert + M(r) trap + M(w)

λ ⊲ M(w): the wave-space kernel G(w)

ξ

is smooth and regular, G(w)(r) = 1 V

• k=0

eik·(r) H(k, ξ) k2 (I − ˆ kˆ k), (6) where the Hasimoto splitting function H(k, ξ) =

• 1 + k2

4ξ2

• e−k2/4ξ2.

⊲ We can efficiently compute the action of M(w) in Fourier space by using the Spectral Ewald method of Lindbo/Tornberg [?], M(w) = S†BS, (7) where S is the non-uniform FFT (NUFFT) of Greengard/Lee [?], and B is a SPD block-diagonal matrix (in Fourier space), B(k, ξ) = H(k, ξ) k2 (I − ˆ kˆ k).

• Y. Bao (CIMS)

Fluctuating BIE

Far-Field Random Slip Velocity

⊲ Key Idea 3: random slip velocity with covariance M is generated by M

1 2 W

d.

=

• M(w) 1

2 W(w) +

• M(r) 1

2 W(r),

(8) if both M(w) and M(r) are SPD and W(w)W(r) = 0. Also taken from recent work of Fiore et al. [?] ⊲ The far-field piece M(w) is SPD by construction and we can write M(w) = S†BS =

• S†B

1 2

S†B

1 2

† , (9) so that the wave-space random slip velocity can be generated with a single call to the NUFFT,

• M(w) 1

2 W(w) = S†B 1 2 W(w),

(10) ⊲ This is equivalent to how Brownian displacements are generated in methods like the Fluctuating Immersed Boundary [?] and the fluctu- ating Force Coupling Method [?] by using fluctuating hydrodynamics (putting stochastic forcing on fluid rather than on particles).

• Y. Bao (CIMS)

Fluctuating BIE

Far-Field Random Slip Velocity

⊲ Key Idea 3: random slip velocity with covariance M is generated by M

1 2 W

d.

=

• M(w) 1

2 W(w) +

• M(r) 1

2 W(r),

(8) if both M(w) and M(r) are SPD and W(w)W(r) = 0. Also taken from recent work of Fiore et al. [?] ⊲ The far-field piece M(w) is SPD by construction and we can write M(w) = S†BS =

• S†B

1 2

S†B

1 2

† , (9) so that the wave-space random slip velocity can be generated with a single call to the NUFFT,

• M(w) 1

2 W(w) = S†B 1 2 W(w),

(10) ⊲ This is equivalent to how Brownian displacements are generated in methods like the Fluctuating Immersed Boundary [?] and the fluctu- ating Force Coupling Method [?] by using fluctuating hydrodynamics (putting stochastic forcing on fluid rather than on particles).

• Y. Bao (CIMS)

Fluctuating BIE

Far-Field Random Slip Velocity

⊲ Key Idea 3: random slip velocity with covariance M is generated by M

1 2 W

d.

=

• M(w) 1

2 W(w) +

• M(r) 1

2 W(r),

(8) if both M(w) and M(r) are SPD and W(w)W(r) = 0. Also taken from recent work of Fiore et al. [?] ⊲ The far-field piece M(w) is SPD by construction and we can write M(w) = S†BS =

• S†B

1 2

S†B

1 2

† , (9) so that the wave-space random slip velocity can be generated with a single call to the NUFFT,

• M(w) 1

2 W(w) = S†B 1 2 W(w),

(10) ⊲ This is equivalent to how Brownian displacements are generated in methods like the Fluctuating Immersed Boundary [?] and the fluctu- ating Force Coupling Method [?] by using fluctuating hydrodynamics (putting stochastic forcing on fluid rather than on particles).

• Y. Bao (CIMS)

Fluctuating BIE

Near-Field Random Slip Velocity

⊲ For the near-field random slip velocity, we need to generate the action

• f the square root of a sparse matrix,

M(r) = M(r)

Alpert + M(r) trap

⊲ For sparse matrices, the principal square root can be efficiently com- puted by a Krylov Lanczos method of Chow/Saad [?]. ⊲ In general, M(r)

Alpert is not symmetric, so M(r) is not SPD strictly

• speaking. Nevertheless, we find that symmetrizing M(r)

Alpert preserves

the order of accuracy of Alpert quadrature, and the Krylov Lanczos iteration is rather insensitive to any small negative eigenvalues.

• Y. Bao (CIMS)

Fluctuating BIE

Near-Field Random Slip Velocity

⊲ For the near-field random slip velocity, we need to generate the action

• f the square root of a sparse matrix,

M(r) = M(r)

Alpert + M(r) trap

⊲ For sparse matrices, the principal square root can be efficiently com- puted by a Krylov Lanczos method of Chow/Saad [?]. ⊲ In general, M(r)

Alpert is not symmetric, so M(r) is not SPD strictly

• speaking. Nevertheless, we find that symmetrizing M(r)

Alpert preserves

the order of accuracy of Alpert quadrature, and the Krylov Lanczos iteration is rather insensitive to any small negative eigenvalues.

• Y. Bao (CIMS)

Fluctuating BIE

Near-Field Random Slip Velocity

⊲ For the near-field random slip velocity, we need to generate the action

• f the square root of a sparse matrix,

M(r) = M(r)

Alpert + M(r) trap

⊲ For sparse matrices, the principal square root can be efficiently com- puted by a Krylov Lanczos method of Chow/Saad [?]. ⊲ In general, M(r)

Alpert is not symmetric, so M(r) is not SPD strictly

• speaking. Nevertheless, we find that symmetrizing M(r)

Alpert preserves

the order of accuracy of Alpert quadrature, and the Krylov Lanczos iteration is rather insensitive to any small negative eigenvalues.

• Y. Bao (CIMS)

Fluctuating BIE

Block-Diagonal Preconditioners

⊲ To mitigate the inherent ill-conditioning of M due to the use of a first-kind boundary integral formulation, we apply a block-diagonal preconditioner, i.e., we simply neglect all hydrodynamic interactions between distinct bodies in the preconditioner, both when solving the saddle-point mobility problem using GMRES, and in the Lanczos iteration for generating

• M(r) 1

2 W(r).

⊲ Both preconditioners can be precomputed using LAPACK for a sin- gle body, and then applied to many bodies via two fast vector rota- tions per body. ⊲ GMRES and Lanczos converge in a constant number of iterations, growing only weakly with packing density.

• Y. Bao (CIMS)

Fluctuating BIE

Block-Diagonal Preconditioners

⊲ To mitigate the inherent ill-conditioning of M due to the use of a first-kind boundary integral formulation, we apply a block-diagonal preconditioner, i.e., we simply neglect all hydrodynamic interactions between distinct bodies in the preconditioner, both when solving the saddle-point mobility problem using GMRES, and in the Lanczos iteration for generating

• M(r) 1

2 W(r).

⊲ Both preconditioners can be precomputed using LAPACK for a sin- gle body, and then applied to many bodies via two fast vector rota- tions per body. ⊲ GMRES and Lanczos converge in a constant number of iterations, growing only weakly with packing density.

• Y. Bao (CIMS)

Fluctuating BIE

Block-Diagonal Preconditioners

⊲ To mitigate the inherent ill-conditioning of M due to the use of a first-kind boundary integral formulation, we apply a block-diagonal preconditioner, i.e., we simply neglect all hydrodynamic interactions between distinct bodies in the preconditioner, both when solving the saddle-point mobility problem using GMRES, and in the Lanczos iteration for generating

• M(r) 1

2 W(r).

⊲ Both preconditioners can be precomputed using LAPACK for a sin- gle body, and then applied to many bodies via two fast vector rota- tions per body. ⊲ GMRES and Lanczos converge in a constant number of iterations, growing only weakly with packing density.

• Y. Bao (CIMS)

Fluctuating BIE

Numerical Results

⊲ This proof-of-concept algorithm/implementation is in 2D only, but the main ideas can be carried over to 3D in principle (but with some technical difficulties that need to be overcome!).

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 φ = 0.25, Nbody = 100

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

φ = 0.5, Nbody = 100

Figure: Random configurations of 100 disks with packing ratio φ = 0.25 (low density) and φ = 0.5 (moderately high density)

• Y. Bao (CIMS)

Fluctuating BIE

Accuracy of N F

32 64 128

• Num. pts. per body

10 -14 10 -12 10 -10 10 -8 10 -6 10 -4 10 -2 Relative error in U = N F Area fraction φ = 0.25 1st-kind, 4th-order 1st-kind, 8th-order 2nd-kind 32 64 128

• Num. pts. per body

10 -11 10 -9 10 -7 10 -5 10 -3 Relative error in U = N F Area fraction φ = 0.5 1st-kind, 4th-order 1st-kind, 8th-order 2nd-kind

Figure: Accuracy of 1st- and 2nd-kind (spectral in 2D!) mobility solvers for dilute and dense hard-disk suspensions. While the 2nd kind gives spectral accuracy and converges faster with number of DOFs, the first kind is more accurate for low resolutions especially at higher densities (but what about 3D?).

• Y. Bao (CIMS)

Fluctuating BIE

Convergence and robustness (2D specific!)

10 20 30 40 50 60 Number of GMRES iter. 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 Relative residual 1st- vs 2nd-kind Mobility problem, GMRES (φ = 0.5) ξ ≈ 20, 1st-kind ξ ≈ 20, 2nd-kind ξ ≈ 100, 1st-kind ξ ≈ 100, 2nd-kind 10 20 30 40 50 Number of Lanczos iter. 10 -8 10 -6 10 -4 10 -2 10 0 Relative residual (M(r))

1 2 W, Lanczos iter., (φ = 0.5)

ξ ≈ 20 ξ ≈ 30 ξ ≈ 40 ξ ≈ 50 ξ ≈ 60 ξ ≈ 70 ξ ≈ 80 ξ ≈ 90 ξ ≈ 100

Figure: We expect much better scaling in 3D due to faster decay of Stokeslet.

• Y. Bao (CIMS)

Fluctuating BIE

Efficiency and Scaling

20 40 60 80 100 ξ 10 -1 10 0 10 1 10 2 10 3 10 4 Timing (in seconds) φ = 0.5, Nbody = 100, Alpert = 4 export StokesletRsum export StressletRsum

• Mob. GMRES, 1st-kind
• Mob. GMRES, 2nd-kind

generate (M(r))

1 2 W(r)

generate (M(w))

1 2 W(w)

Total 1st-kind Total 2nd-kind

Nbody

200 400 600 800

Timing (in seconds)

100 200 300 400

Timing of FBEM (φ = 0.5, ξ ≈ 80, norder = 4)

FBEM, 1st-kind FBEM, 2nd-kind

Figure: Left: Optimal Ewald splitting parameter, specific to our 2D Matlab

• implementation. In 3D one expects to see an optimal splitting parameter (see

Fiore et al. [?]). Right: Linear scaling of the algorithm with the number of bodies.

• Y. Bao (CIMS)

Fluctuating BIE

Conclusion

⊲ Ewald (Hasimoto) splitting can be used to accelerate both deter- ministic and stochastic simulations in periodic domains. ⊲ Key is to ensure both far-field and near-field are (essentially) SPD so one piece is generated using FFTs and the other using iterative methods. ⊲ Using these principles we have constructed a linear-scaling fluctu- ating boundary element method for Brownian suspensions. ⊲ Can a similar idea be used with boundary integral methods based on grid-free fast multipole methods (e.g., unbounded domains)?

• Y. Bao (CIMS)

Fluctuating BIE

References

• A. M. Fiore, F. Balboa Usabiaga, A. Donev, and J. W. Swan.

Rapid Sampling of Stochastic Displacements in Brownian Dynamics Simulations. nov 2016. Dag Lindbo and Anna-Karin Tornberg. Spectrally accurate fast summation for periodic stokes potentials. Journal of Computational Physics, 229(23):8994–9010, 2010.

• L. Greengard and J. Lee.

Accelerating the nonuniform fast fourier transform. SIAM Review, 46(3):443–454, 2004.

• S. Delong, F. Balboa Usabiaga, R. Delgado-Buscalioni, B. E. Griffith, and A. Donev.

Brownian Dynamics without Green’s Functions.

• J. Chem. Phys., 140(13):134110, 2014.

Software available at https://github.com/stochasticHydroTools/FIB. Eric E. Keaveny. Fluctuating force-coupling method for simulations of colloidal suspensions.

• J. Comp. Phys., 269(0):61 – 79, 2014.

Edmond Chow and Yousef Saad. Preconditioned krylov subspace methods for sampling multivariate gaussian distributions. SIAM Journal on Scientific Computing, 36(2):A588–A608, 2014.

• Y. Bao (CIMS)

Fluctuating BIE