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 22 nd , 2014 Math REU funded by NSF under DMS-1262877. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 28

Introduction Solutions may be only piecewise smooth, and solve the equation in a weak sense: � � uϕ t + α ( u ) ∇  ϕ = 0 R n × (0 , T ) 0 ( R n × (0 , T )) ∀ ϕ ∈ C ∞ July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 28

Initial Condition The initial profile of the concentration, u 0 , must be given as a function of spatial variables: u 0 = u 0 ( x 1 , x 2 , · · · , x n ) x ∈ R n . Where � x = x  ˆ e  + x  ˆ e  + · · · + x n ˆ e n is a position vector � July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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 * ! + ! ( x 0 , t 0 ) ! " S u ! u * Figure 1 : The free boundary and normal vector. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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 . July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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: � � Ω + ∪ Ω − [ ϕ t u + ∇ 2 x ϕ · α ( u )] dxdt = 0 July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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: � ϕ [ ∇ x α ( u ∗ ) − ∇ x α ( u ∗ ) , − ( u ∗ − u ∗ )] · ν ( x, t ) ds = 0 , S ∀ ϕ ∈ C ∞ of a neighborhood of ( x 0 , t 0 ). 0 July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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 → R n Lipschitz, then ( ψ ′ ( t ) , 1) ν ( x, t ) = � ( ψ ′ ( t ) , 1) � and, ψ ′ ( t ) = x ′ = − ∂ x α ( u ∗ ) − ( − ∂ x α ( u ∗ )) −  = [ − ∂ x α ( u )] u ∗ − u ∗ [ u ] Where x ′ ( t ) is the velocity of the free boundary. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 28

Riemann Problem Let α ( u ) = H ( u − c ), where H is the Heaviside step function. 1 c Figure 2 : The step function. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 28

Riemann problem cont. The data would be u I = a for u < c , and u I = b for u > c . Since α ( u ) needs to be continuous, we redefine U ( x, t ) = α ( u ( x, t )) as,  U = 0 x < s ( t )   x − s ( t ) U = s ( t ) ≤ x ≤ d ( t ) d ( t ) − s ( 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 )  July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 28

Riemann problem cont. Then, ( u , U ) solves u t = ∇  U as u t ≡ 0 = ∇  U in s ( t ) < x < d ( t ). The jump condition is satisfied along the free boundaries ( s ( t ) , t ) and ( d ( t ) , t ). July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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 [ x i −  , x i ] ,  ≤ i ≤ k . Then, we discretize the initial data as u i = u ( x i ), and α ( u ) as α i = α ( u i ). July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 28

Discretized form We move the grid points x i according to the jump condition as obtained in the Riemann problem.  � α i +  − α i − α i − α i −  � x i ( t ) = − ˙ u i +  − u i x i +  − x i x i − x i −  Where α i = α ( u i ), with i ∈ [ , k +  ]. Here u i is the exact value of u at the point x i , and α i the exact value of U at x i . July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 28

Runge-Kutta Method In order to find the correct time step, we use the RK 4th order method. x ( x j k  = ˙ i ) x ( x j k  = ˙ i + ∆t · k  / ) x ( x j k  = ˙ i + ∆t · k  / ) x ( x j k  = ˙ i + ∆t · k  ) i + ∆t x j +  = x j  ( k  + k  + k  + k  ) i Then we compare it with RK4 using two half-steps. From this comparision we know if it is necessary to shrink the time step. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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 x i will have moved to a new position. When two free boundaries get close enough ( x i +  − x i ) < δ , where δ is a specified tolerance, the level set x i +  is deleted, as is the corresponding level u i +  . The piecewise constant function u ( x, h ) = u i on [ x i , x i +  ] obtained this way is the approximate solution. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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 x i will have moved to a new position. When two free boundaries get close enough ( x i +  − x i ) < δ , where δ is a specified tolerance, the level set x i +  is deleted, as is the corresponding level u i +  . The piecewise constant function u ( x, h ) = u i on [ x i , x i +  ] obtained this way is the approximate solution. Movie time. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 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 R n . 1 ∂ � r n − 1 ∂ U � = 0 r n − 1 ∂ r ∂ r Where U is the redefined α ( u ) for the discontinuities of u between two consecutive free boundaries. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 28

Radial Solutions in Higher Dimensions cont. For n = 2  U r ( r, t ) = rln [ s ( t ) /d ( t )] For n > 2 U r ( 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. July 22 nd , 2014 Math REU funded by Jacob Cheverie and Leandro Fosque Mentor: Marianne Korten / 28