Parallel multigrid methods for parabolic partial differential - - PowerPoint PPT Presentation

Parallel multigrid methods for parabolic partial differential equations and applications

Feng Wei Yang

Department of Mathematics University of Sussex F.W.Yang@sussex.ac.uk

1 October 2015

Feng Wei Yang Seminar at INI 1 October 2015 1 / 39

Objectives

To solve complex non-linear parabolic systems by applying: Cartesian Grids (FDM) Implicit Schemes (BDFs) Nonlinear Multigrid Method with Full Approximation Scheme (FAS) Adaptive Mesh Refinement (AMR) Adaptive Time-Stepping (ATS) Parallel Technique

Feng Wei Yang Seminar at INI 1 October 2015 2 / 39

An example of the so-called ”stiff” problems

Binary alloy solidification in 3D

• P. Bollada, C.E. Goodyer, P.K. Jimack, A.M. Mullis, F.W. Yang

Journal of Computational Physics, 2015 Feng Wei Yang Seminar at INI 1 October 2015 3 / 39

Outline

Multigrid methods Thin film model from Gaskell et al. Optimal control model with geometric evolution laws for whole cell tracking

Feng Wei Yang Seminar at INI 1 October 2015 4 / 39

Jacobi/Gauss-Seidel iterative methods

Well-known methods Require diagonally-dominant matrices Typically have complexity of O(n2) for general sparse matrices ... Smoothing property

Low frequency of error High frequency of error

S.H. Lui Numerical Analysis of Partial Differential Equations, 2011 Feng Wei Yang Seminar at INI 1 October 2015 5 / 39

Convergence of a typical Jacobi iterative method

source: nkl.cc.u-tokyo.ac.jp Feng Wei Yang Seminar at INI 1 October 2015 6 / 39

Multigrid V-cycle

Finest grid Coarsest grid Grid level 1 Grid level 2 Grid level 3 Grid level 4 x y

Feng Wei Yang Seminar at INI 1 October 2015 7 / 39

Linear multigrid

A linear problem: Au = b, (1) exact error can be obtained as E = u − v, (2) residual can be calculated as: r = b − Av. (3) Error equation: AE = A(u − v) = Au − Av = b − Av = r. (4)

Feng Wei Yang Seminar at INI 1 October 2015 8 / 39

Linear multigrid

Feng Wei Yang Seminar at INI 1 October 2015 9 / 39

Nonlinear multigrid

The Error Equation (4) does not exist in a nonlinear case Full Approximate Scheme (FAS) For problem on coarser grids, a modified RHS is included

Feng Wei Yang Seminar at INI 1 October 2015 10 / 39

Nonlinear multigrid

Feng Wei Yang Seminar at INI 1 October 2015 11 / 39

A nonlinear point-wise smoother

Let’s consider our nonlinear problem: A(v) = f . It can be rewritten as: F(v) = 0. Then the Newton-like nonlinear point-wise smoother at a particular grid point (i, j) ∈ Ω can be the following: vℓ+1,t+1

i,j

= vℓ,t+1

i,j

− F(v) F′(vℓ,t+1

i,j

) .

Feng Wei Yang Seminar at INI 1 October 2015 12 / 39

Domain decomposition and guard cells

boundary

Guard cells as boundary points Guard cells that store values of corresponding grid points on neighboring blocks

Feng Wei Yang Seminar at INI 1 October 2015 13 / 39

Multigrid in parallel

Finest grid Coarsest grid Grid level 1 Grid level 4

Feng Wei Yang Seminar at INI 1 October 2015 14 / 39

Droplet spreading model

∂h ∂t =

∂ ∂x

• h3

3

• ∂p

∂x − Bo ǫ sinα

• + ∂

∂y

• h3

3

• ∂p

∂y

• p =

−△(h) − Π(h) + Boh cos α with boundary conditions: h = h∗ ∂np = 0

• n ∂Ω

Gaskell et al. Int. J. Numer. Meth. Fluids, 45:1161-1186, 2004 Feng Wei Yang Seminar at INI 1 October 2015 15 / 39

Our solver

Cell-centred 2nd order finite difference method PARAMESH library for mesh generation and AMR Fully implicit BDF2 method with adaptive time-stepping MLAT variation of FAS multigrid at each time-step Newton-block 2 × 2 Red-Black (weighted) Gauss-Seidel smoother Full weighting restriction and bilinear interpolation Parallelisation through domain decomposition

Feng Wei Yang Seminar at INI 1 October 2015 16 / 39

Newton-block smoother

Update at a grid point (i, j): hℓ+1,t+1 pℓ+1,t+1

• i,j

= hℓ,t+1 pℓ,t+1

• i,j

−  

∂Fh ∂ht+1

i,j

∂Fh ∂pt+1

i,j

∂Fp ∂ht+1

i,j

∂Fp ∂pt+1

i,j

 

−1 Fh i,j(h, p)

Fp i,j(h, p)

• Feng Wei Yang

Seminar at INI 1 October 2015 17 / 39

Validation

0.2 0.4 0.6 0.8 1 x 10

−5

1 2 3 4 5 t h0(t) 32x32 64x64 128x128 256x256 512x512 1024x1024

Results from Gaskell et al. on the left and our results on the right.

Feng Wei Yang Seminar at INI 1 October 2015 18 / 39

Multigrid linear complexity

10

4

10

5

10

6

10

7

10 10

1

10

2

• No. grid points on the finest grid.

Average CPU time per time step (seconds). CPU time required Line with slope of 1

A log-log plot demonstrating the linear complexity of multigrid.

Feng Wei Yang Seminar at INI 1 October 2015 19 / 39

Multigrid performance

Results from Gaskell et al. on the left and our results on the right.

Feng Wei Yang Seminar at INI 1 October 2015 20 / 39

AMR

AMR with initial condition on the left and final solution on the right.

Feng Wei Yang Seminar at INI 1 October 2015 21 / 39

Animation

Feng Wei Yang Seminar at INI 1 October 2015 22 / 39

Adaptive time-stepping

0.2 0.4 0.6 0.8 1 x 10

−5

1 2 3 4 5 6 7 8 9 10 11 x 10

−7

Time Time step size adaptive time−stepping 1024x1024

Evolution of δt during T = [0, 1 × 10−5].

Feng Wei Yang Seminar at INI 1 October 2015 23 / 39

Adaptive multigrid solver

Cases No. leaf nodes Uniform 10242 1,048,576 AMR 168,480 Cases

• No. time

CPU time steps (seconds) Fixed δt 1000 16721.3 ATS 45 574.4

• F. Yang et al. Advances in Engineering Software, in review, 2015

Feng Wei Yang Seminar at INI 1 October 2015 24 / 39

Optimal control with geometric evolution laws for whole cell tracking

K.N. Blazakis, A. Madzvamuse, C. Reyes-Aldasoro, V. Styles, C. Venkataraman “Whole cell tracking through the optimal control of geometric evolution laws” Journal of Computational Physics, 2015

• F. Yang, C. Venkataraman, V. Styles, A. Madzvamuse

“A robust and efficient adaptive multigrid solver for the optimal control of phase field formulations of geometric evolution laws” in review, 2015 Feng Wei Yang Seminar at INI 1 October 2015 25 / 39

Model objectives

To track the morphology of cells and reconstruct their movements:

• V. Peschetola et al. Cytoskeleton, 2013

Feng Wei Yang Seminar at INI 1 October 2015 26 / 39

What is our signature

Particle tracking: Particle tracking: Particle tracking: The morphology of cells are not considered Manually tracking is slow Automatic tracking algorithms are often flawed

Segmentation is suboptimal for real data Tracking through patten recognition is challenging

Pure geometric math models: Pure geometric math models: Pure geometric math models: Resolution of the data matters Typically no cell-setting is considered It is a complicated procedure to

• btain the results

Computational power and advanced numerical methods have to be included for 3-D real-life cell tracking

Feng Wei Yang Seminar at INI 1 October 2015 27 / 39

Our optimal control model

The volume conserved mean curvature flow:

• V

V V (x x x, t) = (−σH(x x x, t) + η(x x x, t) + λV (t))v v v(x x x, t) on Γ(t), t ∈ (0, T], Γ(0) = Γ0.

The phase-field approximation of the above equation - Allen-Cahn:

     ∂tφ(x x x, t) = △φ(x x x, t) − 1

ǫ2 G ′(φ(x

x x, t)) − 1

ǫ(η(x

x x, t) − λ(t)) in Ω × (0, T], ∇φ · ν ν νΩ = 0 on ∂Ω × (0, T], φ(·, 0) = φ0 in Ω.

Feng Wei Yang Seminar at INI 1 October 2015 28 / 39

Our optimal control model cont.

The objective functional:

J(φ, η) = 1 2

• Ω

(φ(x x x, T) − φobs(x x x))2 dx x x + θ 2 T

• Ω

η(x x x, t)2dx x xdt,

and now we solve the minimisation problem:

minηJ(φ, η), with J given above.

Feng Wei Yang Seminar at INI 1 October 2015 29 / 39

Our optimal control model cont.

The adjoint equation to help computing the derivative of the objective functional:

• ∂tp(x

x x, t) = −△p(x x x, t) + ǫ−2G ′′ (φ (x x x, t))p(x x x, t) in Ω × [0, T), p(x x x, T) = φ(x x x, T) − φobs(x x x) in Ω,

and we update the control as

ηℓ+1 = ηℓ − α

• θηℓ + 1

ǫ pℓ

• .

Feng Wei Yang Seminar at INI 1 October 2015 30 / 39

Numerical challenges

Number of time steps Memory requirement (let’s consider double precision and 100 time steps)

2-D: 5122 requires 0.4 gigabytes 3-D: 5123 requires 215 gigabytes

Feng Wei Yang Seminar at INI 1 October 2015 31 / 39

Two-grid solution strategy

One time step One complete solve for the Allen-Cahn equation from t=(0,T] Intermediate grid(s) Restrict the converged solution

• f ϕ

Fine grid for the Allen-Cahn equation Coarse grid for the adjoint equation One time step One complete solve for the adjoint equation from t=[T,0) Interpolate the computed η Start the next η iteration

Feng Wei Yang Seminar at INI 1 October 2015 32 / 39

Real world example (1)

Feng Wei Yang Seminar at INI 1 October 2015 33 / 39

Real world example (1)

t=0 t=T

Feng Wei Yang Seminar at INI 1 October 2015 34 / 39

Real world example (1)

Feng Wei Yang Seminar at INI 1 October 2015 35 / 39

Real world example (2)

Feng Wei Yang Seminar at INI 1 October 2015 36 / 39

Euler number for topological changes

We compute this Euler number for these time steps with an ”optimized“ control η: X = 1 2π(a − b)

• Ω(a,b)
• −△φ + ∇|∇φ|2 · ∇φ

2|∇φ|2

• dx.
• Q. Du et al. J. Appl. Math., 2005

Feng Wei Yang Seminar at INI 1 October 2015 37 / 39

Real world example (2)

Feng Wei Yang Seminar at INI 1 October 2015 38 / 39

A 3-D example

Feng Wei Yang Seminar at INI 1 October 2015 39 / 39