Reduced Basis Method for Poisson-Boltzmann Equation Workshop in - - PowerPoint PPT Presentation

Reduced Basis Method for Poisson-Boltzmann Equation

Workshop in Industrial and Applied Mathematics, WIAM16 Cleophas Kweyu, Lihong Feng, Matthias Stein, Peter Benner

September 01, 2016

Partners:

Outline

• 1. Motivation
• 2. Introduction
• 3. Finite Difference Discretization
• 4. Essentials of Reduced Basis Method (RBM)
• 5. Numerical Results
• 6. Conclusions and Outlook

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 2/24

Motivation

Electrostatic Interactions

[Holst ’94]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Complexity of a charged particle in solution surrounded by other charged particles.

Figure: 2-D view of the 3-D Debye-H¨ uckel model.

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 3/24

Introduction

Poisson-Boltzmann Equation (PBE)

[Holst ’94]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

PBE

− ∇.(ǫ(x) ∇u(x)) + ¯ k2(x) sinh(u(x)) = (4πe2

c

kBT )

Nm

• i=1

ziδ(x − xi), in Ω ∈ R3, u(x) = ( e2

c

kBT )

Nm

• i=1

zie−k(d−ai) ǫw(1 + kai)d on ∂Ω, d = |x − xi|, (1) u(∞) = 0. k2 =

8πe2

c I

1000ǫkBT , (I = µ) = 1 2

N

i=1 ciz2 i ,

u(x) = ecψ(x)

kBT .

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 4/24

Introduction

Poisson-Boltzmann Coefficients

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ǫ(x) =

• ǫ1 ≈ 2

if x ∈ Ω1 ǫ2(= ǫ3) ≈ 80 if x ∈ Ω2or Ω3 , ¯ k(x) =

• if x ∈ Ω1or Ω2

√ǫ3k if x ∈ Ω3

Figure: PBE coefficients

Source: Introduction to Molecular Electrostatics with APBS, Robert Konecny

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 5/24

Introduction

Linearized PBE (LPBE)

[Fogolari et al ’99,Holst ’94]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Assumption: ψ(x) ≪ 1.

LPBE

− ∇.(ǫ(x) ∇u(x)) + ¯ k2(x)u(x) = (4πe2

c

kBT )

Nm

• i=1

ziδ(x − xi), (2)

Applications of the PBE and LPBE

potential at the surface of a biomolecule - docking sites, potential outside the molecule - free energy of interaction, free energy of a biomolecule - biomolecular stability.

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 6/24

Finite Difference Discretization

Centered finite differences of LPBE

[Simakov2013]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

− 1 dx2 ǫi+ 1

2 ,j,k(ui+1,j,k − ui,j,k) +

1 dx2 ǫi− 1

2 ,j,k(ui,j,k − ui−1,j,k) −

1 dy2 ǫi,j+ 1

2 ,k(ui,j+1,k − ui,j,k)

+ 1 dy2 ǫi,j− 1

2 ,k(ui,j,k − ui,j−1,k) −

1 dz2 ǫi,j,k+ 1

2 (ui,j,k+1 − ui,j,k) +

1 dz2 ǫi,j,k− 1

2 (ui,j,k − ui,j,k−1)

+ ¯ k2

i,j,kui,j,k = Cqi,j,k.

(3)

(a) Discretization of continuous variables (b) Molecular surfaces and volumes

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 7/24

Essentials of Reduced Basis Method (RBM)

Introduction

[Benner et al ’2015, Eftang ’2011]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Model Reduction: FOM to ROM

Replace FOM AuN (µ) = f N (µ), µ ∈ D, with ROM AuN(µ) = fN(µ), uN(µ) ≈ uN (µ), N ≪ N. RBM is a parametrized model order reduction (PMOR) technique, exploits an offline/online procedure, powerful tools - greedy algorithm and a posteriori error estimation, assumption - typically low dimensional solution manifold, MN = {uN (µ) : µ ∈ D}. (4) RB space V is built upon 4 - generated by greedy algorithm, range(V ) = span{uN (µ1), ..., uN (µN)}, ∀µ1, ..., µN ∈ D. (5)

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 8/24

Essentials of Reduced Basis Method (RBM)

Greedy Algorithm

[Hesthaven et al 2014]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Algorithm 1 Greedy algorithm Input: Training set Ξ ⊂ D including all of µ, i.e., Ξ := {µ1, . . . , µl}. Output: RB basis represented by the projection matrix V .

1: Choose µ∗ ∈ Ξ arbitrarily 2: Solve FOM for uN (µ∗) 3: S1 = {µ∗}, V1 = [uN (µ∗)], N = 1 4: while max

µ∈Ξ ∆N(µ) ≥ ǫ do

5:

µ∗ = arg max

µ∈Ξ ∆N(µ)

6:

Solve FOM for uN (µ∗)

7:

SN+1 = SN ∪ µ∗, VN+1 = [VN uN (µ∗)]

8:

Orthonormalize the columns of VN+1

9:

N = N + 1

10: end while

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 9/24

Essentials of Reduced Basis Method (RBM)

Computational complexity of the Reduced order Model (ROM)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Nonaffine parameter dependence

(A1 + µA2)uN (µ) = ρ + b(µ), µ ∈ D. (6)

Consider the reduced order model (ROM);

( ˆ A1

• N×N

+µ ˆ A2

• N×N

) uN(µ)

N×1

= ˆ ρ

• N×1

+ V T

• N×N

b(µ)

• N×1

, (7) where ˆ A1 = V TA1V , ˆ A2 = V TA2V , ˆ ρ = V Tρ, and N ≪ N. matrix-vector products require 2NN flops, full evaluation of b(µ).

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 10/24

Essentials of Reduced Basis Method (RBM)

Discrete Empirical Interpolation Method (DEIM)

[Chaturantabut 2010]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 10 15 20 10−15 10−5 105 Number of singular values Singular values

Figure: Decay of singular values

Compute snapshot matrix F = [b(µ1), . . . , b(µl)] ∈ RN×l, apply SVD to F : F = UFΣW T, UF ∈ RN×l, Σ ∈ Rl×l, and W ∈ Rl×l, Σ = diag(σ1, . . . , σl) s.t, σ1 > . . . > σl ≥ 0,

l

• i=r+1

σi

l

• 1=1

σi

< ǫsvd, ǫsvd = 10−13.

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 11/24

Essentials of Reduced Basis Method (RBM)

Discrete Empirical Interpolation Method (DEIM)

[Volkwein 2010]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Select basis set {uF

i }r i=1 of rank r from UF which solves,

arg min

{uF

i }r i=1

l

j=1 xj − r i=1xj, uF i uF i 2 2,

s.t.ui, uj = δij, DEIM determines UFc(µ) s.t, b(µ) ≈ UFc(µ), c(µ) ∈ Rr, determine c(µ) by selecting r rows from b(µ) = UFc(µ), suppose PTU is nonsingular, for P = [e℘1, . . . , e℘r ] ∈ RN×r, then, PTb(µ) = PTUFc(µ) = ⇒ c(µ) = (PTUF)−1PTb(µ), (8) ∴ b(µ) ≈ UF(PTUF)−1PTb(µ). (9) ROM with DEIM approximation becomes, ( ˆ A1

• N×N

+µ ˆ A2

• N×N

) uN(µ)

N×1

= ˆ ρ

• N×1

+ V TUF(PTUF)−1

• N×r

PTb(µ)

r×1

. (10)

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 12/24

Essentials of Reduced Basis Method (RBM)

Discrete Empirical Interpolation Method (DEIM)

[Chaturantabut 2010]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Algorithm 2 DEIM algorithm Input: Basis {uF

i }r i=1 for F.

Output: DEIM basis UF and indices ℘ = [℘1, . . . , ℘r]T ∈ Rr.

1: ℘1 = arg max{|uF

1 |},

2: UF = [uF

1 ], P = [e℘1],

℘ = [℘1].

3: for i = 2 to r do 4:

Solve (PTUF)α = PTuF

i for α, where α = (α1, . . . , αi−1)T,

5:

r = uF

i − UFα,

6:

℘i = arg max{|r|},

7:

UF ← [UF uF

i ], P ← [P e℘i],

℘ ← ℘ ℘i

• .

8: end for

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 13/24

Essentials of Reduced Basis Method (RBM)

DEIM Approximation Error

[Feng et al 2016, Wirtz et al 2014]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

DEIM error is given by, eDEIM = b(µ) − ˜ b(µ) = Π2(I − Π)b(µ), (11) Π and Π2 are oblique projectors defined as, Π = UF(PTUF)−1PT, (12) Π2 = (I − Π)˜ UF(˜ PT(I − Π)˜ UF)−1 ˜ PT, (13) ˜ UF = U∗

F(:, r + 1 : r∗) and ˜

P = P∗(:, r + 1 : r∗) such that U∗

F = [UF, ˜

UF] and P∗ = [P, ˜ P], b(µ) = U∗

F((P∗)TU∗ F)−1(P∗)Tb(µ).

(14)

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 14/24

Essentials of Reduced Basis Method (RBM)

A Posteriori Error Estimation

[Quarteroni et al 2016]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Compute residual due to DEIM, rDEIM

N

(uN; µ) = (ρ + ˜ b(µ)) − AN (µ)uN(µ), (15) general residual becomes, rN(uN; µ) = (ρ + b(µ)) − AN (µ)uN(µ) = (ρ + ˜ b(µ)) − AN (µ)uN(µ) + b(µ) − ˜ b(µ) = rDEIM

N

(uN; µ) + b(µ) − ˜ b(µ)

• :=eDEIM

. (16) a posteriori error can be derived from 16 by, rN(uN; µ) = AN (µ)uN (µ) − AN (µ)uN(µ) = AN (µ)e(µ) (17)

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 15/24

Essentials of Reduced Basis Method (RBM)

A Posteriori Error Estimation

[Quarteroni ’2015]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

The error e(µ) := uN (µ) − uN(µ) is given by e(µ) = (AN (µ))−1rN(uN; µ), (18)

• btain an upper bound for the 2-norm of the error,

e(µ)2 ≤ (AN )−1(µ)2rN(uN; µ)2 = rN(uN; µ)2 σmin(AN (µ)) =: ∆N(µ), (19) where σmin(AN (µ)) is the smallest singular value of AN (µ), in our case the a posteriori error is, e(µ)2 ≈ rN(uN; µ)2 = ∆N(µ). (20)

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 16/24

Numerical Results

Finite Difference Results

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 4 6 8 10 10−12 10−8 10−4 100 Iteration number

Figure: Relative residual for PCG

Consider LPBE (1), parameter domain µ ∈ D = [0.05, 0.15], physical domain Ω = 60˚ A × 60˚ A × 60˚ A, dimension N = 2, 146, 689, PQR file, Cubic B-spline interpolation (basis spline), PCG with algebraic multigrid v-cycle preconditioner.

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 17/24

Numerical Results

Finite Difference Results

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Computational time to solve uN (µ) is ≈ 50 seconds on average.

Figure: uN (µ) at µ = 0 Figure: uN (µ) at µ = 0.05

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 18/24

Numerical Results

Finite Difference Results

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure: uN (µ) at µ = 0.15 Figure: uN (µ) at µ = 0.5

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 19/24

Numerical Results

Reduced Basis Results

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Computational time to solve uN(µ) is ≈ 0.065 seconds on average. True error = uN(µ) − uN(µ)2, ∆max

N

(µ) = max

µ∈Ξ rN(uN; µ)2, Relative ∆max N

(µ) =

∆max

N

(µ) uN(µ∗)2 . µ∗ = arg max µ∈Ξ rN(uN; µ)2.

True error Maximal error 1 2 3 4 5 6 10−5 10−3 10−1 101 103 Reduced Dimension N

(a) Maximal versus true error

1 2 3 4 5 6 10−9 10−7 10−5 10−3 10−1 Reduced Dimension N

(b) Relative ∆max

N

(µ) vs true error Figure: Comparison between true error and maximal error Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 20/24

Numerical Results

Reduced Basis Results

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

True error Error estimate 5 10 15 20 10−8 10−6 10−4 Parameter (µ) sample size

Figure: Error estimate versus true error

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 21/24

Numerical Results

Error analysis between FDM and RBM

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure: Absolute error at µ = 0.05101

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 22/24

Conclusions and Outlook

Conclusions and Outlook

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Conclusions

Applied RBM to LPBE with ionic strength as meaningful parameter. RBM reduces the high dimensional FOM by a factor of ≈ 360, 000 and computational time by a factor of approximately over 800, DEIM error costly in online stage, error estimator provided fast convergence to the RB approximation.

Outlook

Develop an efficient error estimator, reduce DEIM error costs.

Thank you for your attention!

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 23/24

RBM Summary

Conclusions and Outlook

[Quarteroni et al 2016]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure: RB workflow

Cleophas Kweyu, kweyu@mpi-magdeburg.mpg.de WIAM16, August 31- September 2, 2016 24/24