Boundary integral methods for implicitly defined dynamic interfaces - - PowerPoint PPT Presentation

Boundary integral methods for implicitly defined dynamic interfaces

Richard Tsai

with 


• C. Chen, C. Kublik, Y. Wu

KTH Royal Institute of Technology, Sweden
 and
 The University of Texas at Austin, USA

Celebrating Prof. Yoshikazu Giga’s 60th Birthday

Among the important things

• The existence of a second stomach for sweets
• The importance of “unagi” for survival of summer
• The completion of a meal by soba-yu

The L-solutions of Giga-Sato

Embed the solution graph as a level curve of a higher dimensional function, and look at convergence in Hausdorff distance of the subgraphs:

u(x) x y x phi>0 phi<0

ut +H(t,x,u,ux) = 0 ⇐ ⇒ φt −φyH(t,x,y,−φx φy ) = 0.

Multivalued solutions

−0.4 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −0.4 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 x y

v(y)

Very singular diffusion of Giga

φt + ˜ H(t,x,y,φx,φy) = M|∇φ| ∂ ∂y( φy |φy|).

Numerical solution verifies the entropy condition. (Equal area rule).

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 −0.4 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 −0.4 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

Bunching and anisotropy

−2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5

[Ohtsuka-Giga-T:2013] Computation by T. Ohtsuka

Multivalued solutions for high frequency waves

• tS þ jrxSj2

2 þ V ðxÞ ¼ 0;

[Jin,Liu, Osher-T:2005]

Geometric motion of an interface

Mullins-Sekerka dynamics: Ω(t) Ω(t) vn = ∂u ∂n

• normal velocity of a moving interface

∆u = 0 x / ∈ ∂Ω(t) u = κ, x ∈ ∂Ω(t)

+ far field condition

Solution of linear PDEs by boundary integrals

• Solving for :

γ

• Evaluating the solution:

f(x) = λγ(y) + Z

∂Ω

K(x, y)γ(y)dSy u(x) = Z

∂Ω

˜ K(x, y)γ(y)dSy

∆u = 0, x ∈ Ω u(x) = f(x), x ∈ ∂Ω

Breaking of convexity

Volume preservation

−4 −3 −2 −1 1 2 3 4 −4 −3 −2 −1 1 2 3 4

Merging

Nontrivial surface area-volume relation

Isothermal at infinity.

Closed simple curve Γ : γ(s) = (X(s),Y(s)), parametrized by arclength s. I = Z f ◦γ(s)ds = Z

R2 f(x,y)δ(Γ,x,y)dxdy.

φ<0 φ>0

Explicit curves:

δ(Γ,x,y) =

R δ(x − X(s))δ(y − Y(s))ds. Implicit curves:

δ(Γ,x,y) = δ(φ(x,y))|∇φ(x,y)|.

Surface integral: Z

Rd f(x)δ(φ)|∇φ|dx.

[Peskin, Lai]

Eε,h =

• ∑i,j δε(Γ,xi,j)f(xi,j)h2 −

R

Γ f(γ(s))ds

• Typically Eε,h −

→ 0 if ε = hα,0 ≤ α < 1. (wide support!)

Non-convergence if e.g. δcos

Ch

Relative error E = |Sh − S|/S. Radii: r× > r◦ > r.

(a) Γ

200 400 0.06 0.08 0.1 0.12

(b) δL

h

200 400 0.01 0.012 0.014 0.016

(c) δL

2h

10

10

10

10

10

10

(d)

δcos

2h (dΓ)

10

10

10

10

(e)

δcos

2h (phi)

Theorem

∑i,j∈Ip,q δL

ε(dΓ(xi,j))h2 exact,

if (p,q) = ∇dΓ are relative prime and ε = ε(p,q) =

|p|+|q|

√

p2+q2 h.

ε (p,q) Γ

−1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 Polar plot of ε(θn) θn ε(θn)

[Engquist-Tornberg-T:2005]

Defining a “fat” integral equation

• Proper definition of the equation for computation using

implicit surfaces: 
 
 
 
 Need extension of and the unknown from to

• Parameterization of the interface using nearby level sets
• Averaging over different parameterizations

f

γ

∂Ω {−✏ < < ✏}

f(x) = λρ(y) + Z

∂Ω

Φ(x, y)ρ(y)dSy

fext(x) ⇡ λρext(y) + Z

Rd Φ(x, y)ρext(y) δ✏(φ)|rφ|dy ?

∂Ω = {x : φ(x) = 0}

Advantages of implicit interface boundary integral method

• Inherit certain advantages of boundary integral methods: 


solution on exterior domain problem, easy for both Dirichlet, Neumann, Robin, transmission problems.

• Inherit certain benefits from implicit interface formulation: 

• ne grid for complicated, topologically changing

dynamic interface

• Flexible on the grid geometry used to embed the
• interface. Multi-resolution grid possible.

Inverse problem and shape optimization

Ω

Ω

Ω

Analysis of cancellous bone from micro-CT scan

Compute effective material properties.

Segmented by an MBO type scheme [Esedoglu-T]

LIDARs and 3D scanners

New formulas for integration over implicit surfaces

∂Ω ∂Ωη

y∗ z∗ = y∗ + J ~ ds z = y + ~ ds y

Parametrizing by a nearby level set

y∗ = y ηrd(y) Z

∂Ω

ρ(γ(s))Φ(x, γ(s))ds = Z

∂Ωη

ρ(y∗(sη))Φ(x, y∗(sη))Jηdsη

θ1 θ2 ⌘~ n ds dτ

R2 R1

The Jacobian

A0 = (R1 − η)θ1 · (R2 − η)θ2 Aη = ds dτ = R1θ1 · R2θ2 Aη A0 = R1R2 (R1 − η)(R2 − η) = 1 + (κ1 + κ2)η + O(η2) Jη = 1 + 2Hηη + Gηη2 In fact, we have:

Average the identities

h z∗ = z d(z)rd(z) I0 = Z

∂Ω

ρ(γ(s))Φ(x, γ(s))ds = Z

∂Ωη

ρ(y∗(sη))Φ(x, y∗(sη))Jηdsη = Iη I0 = Z

@Ω

ρ(γ(s))Φ(x, γ(s))ds = Z ✏

−✏

δ✏(η)I⌘dη = Z

Rn ρ(z∗)Φ(x, z∗)δ✏(dΩ(z))J(z; dΩ)dz

The singular values of

Proposition 3. The following identity holds for sufficiently small ✏: ˆ

Σ

g(x)dSx = ˆ

R3 g(PΣ(y)) 12K✏((y))dy,

where K✏ is any symmetric kernel supported in [−✏, ✏] having unit mass.

λ1λ2 = 0 at the boundary of Σ

DPΣ

PΣy =y∗ =y φ(y)rφ(y)

φ(y) = min

x∈Σ |y − x|

[Kublik-T:2015]

The Jacobian and singular values of

Proposition 1. Let 1(x) be the largest singular value of DPΓ(x). The following identity holds for ✏ < ∞, where ∞ is the maximal unsigned curvature of Γ. ˆ

Γ

gds = 1 2⇡ ˆ

R3 g(PΓ(x))1(x)K✏()

• dx.

0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 x K(x)=(1−cos(2 π x))/x

δ✏(x; Γ) λ1 = 0 near the end points of Γ

High order approximation

n Relative Error Order 32 6.2030 ⇥ 10−3

• 64

1.8073 ⇥ 10−4 5.10 128 6.6838 ⇥ 10−6 4.76 256 4.1530 ⇥ 10−7 4.01 512 5.0379 ⇥ 10−8 3.04

Integration over a torus:

Surfaces with boundaries

n Relative Error Order 32 1.1726 ⇥ 10−2

• 64

1.1733 ⇥ 10−3 3.32 128 9.1325 ⇥ 10−4 0.36 256 3.8238 ⇥ 10−4 1.26 512 7.8308 ⇥ 10−5 2.29

3/4 sphere and third-order one-sided finite differencing:

n Relative Error Order 60 5.5078 ⇥ 10−3

• 120

1.1476 ⇥ 10−3 2.63 240 2.3409 ⇥ 10−4 2.29 480 3.7166 ⇥ 10−5 2.66

Summation over point clouds

−1 1 −1 1 −1 −0.5 0.5 1

• 30×30 uniformly distributed point clouds sampling in spherical coordinate the quarter

sphere patch.

• 50 × 50 × 50 uniform Cartesian grid discretizing [−1, 1]3.
• Relative error using = 0.05 = dx: −0.56.
• Relative error using = 0.2 = 4dx: −0.061.

has improved regularity

{dΓN = η} P : x 2 {dΓN = η} 7! x ηrdΓN (x) “interpolates”

ΓN ⊂ Σ

Implicit interface boundary integral equation

Discretized by simple Riemann sum over the grid nodes. Solve the resulting linear system by an iterative solver.

x, y ∈ {−✏ ≤ dΩ ≤ ✏} δ✏(x; dΩ) := δ✏(dΩ(x))J(x; dΩ) x∗ = x dΩ(x)rdΩ(x)

[Kublik-Tanushev-T:2013]

f(x∗) = λρ(x∗) + Z

Rd Φ(x∗, y∗)ρ(y∗)δ✏(y; dΩ)dy

Regularization

• Need regularization of near singularity


• Approximate surfaces by paraboloids
• Approximate weakly in a neighborhood of x on the

paraboloid (the 3D case):

x Ω y

∂Φ ∂n ∂Φ ∂n ∂Φ ∂n (x, y) ' A✏(x, y) for ||x y|| < ε Z

U(x;")

∂Φ ∂ny (x, y)α(y)dS(y) ' Z

˜ U(x;")

A✏(x, y)α(y)dS(y)

Regularization

• sculating paraboloid

Z

U(x;h)

∂Φ ∂ny (x, y)dS(y) ' 1 8πh(κx + κy)|U(x; h)|

O P 0

Condition numbers

h

• Cond. number with the tangent regularization
• Cond. number with the paraboloid regularization

4dx 6.2113 7.8322 2dx 7.3343 6.9390 dx 7.0467 7.4197

dx 2

6.8948 6.7572 h

• Cond. number with the tangent regularization
• Cond. number with the paraboloid regularization

4dx 8.0859 7.8024 2dx 8.2791 7.7919 dx 7.4127 8.1231

dx 2

7.8830 8.1192

N = 503 N = 803

Exterior Neumann Problem for the Helmholtz Equation

dx Re(ue) Re(uKT T ) ErrKT T Re(unew) Errnew

4 256 4 512

dx Im(ue) Im(uKT T ) ErrKT T Im(unew) Errnew

4 256 4 512

k = 2.4048255577 . . . ✏ = 0.145 0 = 0.005 ✏0 = 0.15

     ∆u(x) + k2u(x) = 0 x ∈ ¯ Ωc

∂u ∂n(x) = g(x)

x ∈ ∂Ω lim|x|→∞ |x|

1 2 (

∂ ∂|x| − ik)u(x) = 0

Solutions to Helmholtz equation

Happy Birthday!!