Numerical experiments with a level-set tracking algorithm for - - PowerPoint PPT Presentation

Numerical experiments with a level-set tracking algorithm for generalized diffusion equations

Jacob Cheverie and Leandro Fosque

Mentor: Marianne Korten Department of Mathematics Kansas State University

July 22nd, 2014

Math REU funded by NSF under DMS-1262877.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Introduction

We are analyzing a numerical method introduced by J. Bouillet and J. Etcheverry to compute the solution of the generalized diffusion equation by tracking the evolution of level sets.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Introduction

We are analyzing a numerical method introduced by J. Bouillet and J. Etcheverry to compute the solution of the generalized diffusion equation by tracking the evolution of level sets. The diffusion equation is: ∂u ∂t = ∇2α(u) Where α is non-decreasing, and α(0) = 0.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Introduction

The diffusion equation models the process of diffusion and processes that exhibit similar behaviors, such as information propagation, flow in or through porous media, and temperature-dependent phase changes.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Introduction

Solutions may be only piecewise smooth, and solve the equation in a weak sense:

Rn×(0,T)

uϕt + α(u)∇ϕ = 0 ∀ϕ ∈ C ∞

0 (Rn × (0, T))

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Initial Condition

The initial profile of the concentration, u0, must be given as a function of spatial variables: u0 = u0(x1, x2, · · · , xn) Where x = xˆ e + xˆ e + · · · + xnˆ en is a position vector x ∈ Rn.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Free Boundaries

Assuming u(x, t) discontinuous at the curve S(t)={u = u∗}, we can derive the jump condition. We will later use this jump condition to track the evolution of the free boundaries.

!+ !"

u ! u* u ! u*

S

(x0,t0)

!

Figure 1 : The free boundary and normal vector.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Derivation of the jump condition

We can obtain the Rankine-Hugoniot condition for the jump from the lateral limits of u and ∇α(u) at the sides of the boundary S.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Derivation of the jump condition

We can obtain the Rankine-Hugoniot condition for the jump from the lateral limits of u and ∇α(u) at the sides of the boundary S. Let the test function ϕ(x, t) ∈ C ∞

0 (R × (0, T)) have support

containing a neighborhood (x, t). Let Ω+ ∪ S ∪ Ω− ⊃ suppϕ, and ν(x, t) the unit normal vector to S pointing into Ω+. Then:

Ω+∪Ω−[ϕtu + ∇2 xϕ · α(u)]dxdt = 0

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Derivation of the jump condition Cont.

Since α(u) is smooth where u > u∗ and u < u∗, we can work each side separately and skipping many steps, we get from the divergence theorem:

• S

ϕ[∇xα(u∗) − ∇xα(u∗), −(u∗ − u∗)] · ν(x, t)ds = 0, ∀ϕ ∈ C ∞

• f a neighborhood of (x0, t0).

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Derivation cont.

Therefore, [∇xα(u∗) − ∇xα(u∗), −(u∗ − u∗)] · ν(x, t) =  whence [math stuff] is a tangent vector to S. If the interface S is defined by (ψ(t), t), with ψ : R → Rn Lipschitz, then ν(x, t) =

(ψ′(t),1) (ψ′(t),1)

and, − ψ′(t) = x′ = −∂xα(u∗) − (−∂xα(u∗)) u∗ − u∗ = [−∂xα(u)] [u] Where x′(t) is the velocity of the free boundary.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Riemann Problem

Let α(u) = H(u − c), where H is the Heaviside step function.

c

1

Figure 2 : The step function.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Riemann problem cont.

The data would be uI = a for u < c, and uI = b for u > c. Since α(u) needs to be continuous, we redefine U(x, t) = α(u(x, t)) as,      U = 0 x < s(t) U =

x−s(t) d(t)−s(t)

s(t) ≤ x ≤ d(t) U = 1 x > d(t) and u(x, t) as:      u = a x < s(t) u = c s(t) ≤ x ≤ d(t) u = b x > d(t)

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Riemann problem cont.

Then, (u, U) solves ut = ∇U as ut ≡ 0 = ∇U in s(t) < x < d(t). The jump condition is satisfied along the free boundaries (s(t), t) and (d(t), t).

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

The Method

We introduce a moving grid method that tracks level sets as follows: We split the interval [0, m] in k + 1 equal subintervals [xi−, xi],  ≤ i ≤ k. Then, we discretize the initial data as ui = u(xi), and α(u) as αi = α(ui).

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Discretized form

We move the grid points xi according to the jump condition as

• btained in the Riemann problem.

˙ xi(t) = −  ui+ − ui αi+ − αi xi+ − xi − αi − αi− xi − xi−

• Where αi = α(ui), with i ∈ [, k + ]. Here ui is the exact value of u

at the point xi, and αi the exact value of U at xi.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Discretized Form cont.

Figure 3 : Discretized curve example.

Notice that we are working and discretizing only the first quadrant of the x − y coordinates. Because the diffusion is symmetric, we don’t need to work on the second quadrant.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Runge-Kutta Method

In order to find the correct time step, we use the RK 4th order method. k = ˙ x(xj

i)

k = ˙ x(xj

i + ∆t · k/)

k = ˙ x(xj

i + ∆t · k/)

k = ˙ x(xj

i + ∆t · k)

xj+

i

= xj

i + ∆t

 (k + k + k + k) Then we compare it with RK4 using two half-steps. From this comparision we know if it is necessary to shrink the time step.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

R-K Method cont.

Once this system is solved the velocity of each point in the discretized domain is known. Then, given a certain time step h determined by the comparison between the two fourth order Runge-Kutta method, the position values xi will have moved to a new position. When two free boundaries get close enough (xi+ − xi) < δ, where δ is a specified tolerance, the level set xi+ is deleted, as is the corresponding level ui+. The piecewise constant function u(x, h) = ui on [xi, xi+] obtained this way is the approximate solution.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

R-K Method cont.

Once this system is solved the velocity of each point in the discretized domain is known. Then, given a certain time step h determined by the comparison between the two fourth order Runge-Kutta method, the position values xi will have moved to a new position. When two free boundaries get close enough (xi+ − xi) < δ, where δ is a specified tolerance, the level set xi+ is deleted, as is the corresponding level ui+. The piecewise constant function u(x, h) = ui on [xi, xi+] obtained this way is the approximate solution. Movie time.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Radial Solutions in Higher Dimensions

For radial solutions in dimensions n ≥ 2, the expression of U changes: We can obtain the expression for U from the Laplacian for polar coordinates in Rn. 1 rn−1 ∂ ∂r

• rn−1 ∂U

∂r

• = 0

Where U is the redefined α(u) for the discontinuities of u between two consecutive free boundaries.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Radial Solutions in Higher Dimensions cont.

For n = 2 Ur(r, t) =  rln[s(t)/d(t)] For n > 2 Ur(r, t) = ( − n) s(t)n−d(t)n− d(t)n− − s(t)n− r−n The extended U is unique, satisfying U continuous and the jump condition at the free boundaries.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Derivations

Using the divergence form we get the radial jump condition for n dimensions. ˙ r = [− ∂U

∂r ]

[u] For n = 2 we get, ˙ ri = − (ri)(ui+ − ui)

• αi+ − αi

lnri+

ri

− αi − αi− ln ri

ri−

• Jacob Cheverie and Leandro Fosque

Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Derivations cont.

and for, 1 ≤ n = 2 leads to: ˙ ri = −  (ui+ − ui)

• αi+ − αi

r−n

i+ − r−n i

− αi − αi− r−n

i

− r−n

i−

• ( − n)r−n

i

We apply for radial solutions the same method exposed before for the

• ne dimensional case.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Examples

Some preliminaries of the parameters: n is the number of discretization of the x axis. δ is the tolerance of the encounter of xi with xi+. ǫ is the tolerance for the time in the RK process. λ is the distance of spatial grid from origin.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Example: Sin2(x)

• 1

2 3 4 5 6 7 0.2 0.4 0.6 0.8 1.0

Figure 4 : Diffusion of u0 = sin2(x) with α(u) = u.

n = 450, δ = 1.0 × 10−2, ǫ = 1.0 × 10−6, λ = 0.01.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Example: 3 − |x − 3|

• 1

2 3 4 5 6 0.5 1.0 1.5 2.0 2.5 3.0

Figure 5 : Diffusion of u0 = 3 − |x − 3| with α(u) = u.

n = 200, δ = 1.0 × 10−2, ǫ = 1.0 × 10−6, λ = 0.1.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Example: 4 − x2

• 1

2 3 4 5 1 2 3 4

Figure 6 : Diffusion of u0 = 4 − x2 with α(u) = u.

n = 200, δ = 1.0 × 10−9, ǫ = 1.0 × 10−6, λ = 0.1.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Example: e−x2

• 0.5

1.0 1.5 2.0 2.5 3.0 0.2 0.4 0.6 0.8 1.0

Figure 7 : Diffusion of u0 = e−x2 with α(u) = u.

n = 200, δ = 1.0 × 10−9, ǫ = 1.0 × 10−6, λ = 0.1.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

References

Bouillet, J.E. Signed solutions to diffusion-heat conduction equations. Free boundary problems: theory and applications, vol. II (Irsee, 1987), 480-485, Pitman Res. Notes Math. Ser., 186, Longman Sci. Tech., Harlow, 1990. J.E. Bouillet, M.K. Korten, and V. Marquez. Singular Limits and the ”Mesa” Problem. Revista de la Union Matematica Argentina, 41 (1998), 27-39. J.E. Bouillet, J.I. Etcheverry. Numerical Experiments With ut = α(u)xx. Revista de la Union Matematica Argentina, 41 (1998), 15-25.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28

Acknowledgements

The authors are very grateful to their mentor Marianne Korten, to Nathan Albin for invaluable help with the numerical aspect, to the mentors and SUMaR group, and to the KSU math departament.

Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten July 22nd , 2014Math REU funded by / 28