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 ( n 2 ) for general sparse matrices ... Smoothing property High frequency of error Low 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 Grid level 4 Grid level 3 Grid level 2 Coarsest grid Grid level 1 y x 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 (4) = b − Av = r . 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: F ( v ) v ℓ +1 , t +1 = v ℓ, t +1 − . i , j i , j F ′ ( v ℓ, t +1 ) i , j Feng Wei Yang Seminar at INI 1 October 2015 12 / 39

Domain decomposition and guard cells Guard cells that store values of corresponding grid points on neighboring blocks Guard cells as boundary points boundary Feng Wei Yang Seminar at INI 1 October 2015 13 / 39

Multigrid in parallel Finest grid Grid level 4 Coarsest grid Grid level 1 Feng Wei Yang Seminar at INI 1 October 2015 14 / 39

Droplet spreading model ∂ h � � �� � � �� h 3 h 3 ∂ p ∂ p ∂ ∂ x − B o + ∂ ∂ t = ǫ sin α ∂ x 3 ∂ y 3 ∂ y p = −△ ( h ) − Π( h ) + B o h cos α with boundary conditions: h = h ∗ ∂ n p = 0 on ∂ Ω 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 2 nd 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 ): − 1 � F h i , j ( h , p ) � h ℓ +1 , t +1 � h ℓ, t +1 ∂ F h ∂ F h   � � � ∂ h t +1 ∂ p t +1 i , j i , j = − p ℓ +1 , t +1 p ℓ, t +1  ∂ F p ∂ F p  F p i , j ( h , p ) ∂ h t +1 ∂ p t +1 i , j i , j i , j i , j Feng Wei Yang Seminar at INI 1 October 2015 17 / 39

Validation 32x32 64x64 5 128x128 256x256 512x512 1024x1024 4 3 h 0 (t) 2 1 0 0 0.2 0.4 0.6 0.8 1 t − 5 x 10 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 2 10 CPU time required Line with slope of 1 Average CPU time per time step (seconds). 1 10 0 10 4 5 6 7 10 10 10 10 No. grid points on the finest grid. 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 − 7 x 10 11 adaptive time − stepping 1024x1024 10 9 8 7 Time step size 6 5 4 3 2 1 0 0.2 0.4 0.6 0.8 1 Time − 5 x 10 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 1024 2 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 Pure geometric math models: Pure geometric math models: Pure geometric math models: Particle tracking: Particle tracking: Particle tracking: Resolution of the data matters The morphology of cells are not Typically no cell-setting is considered considered Manually tracking is slow It is a complicated procedure to Automatic tracking algorithms obtain the results are often flawed Segmentation is suboptimal Computational power and for real data advanced numerical methods Tracking through patten have to be included for 3-D recognition is challenging 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 ( x x , t ) = ( − σ H ( x x , t ) + η ( x x , t ) + λ V ( t )) v v ( x x , t ) on Γ( t ) , t ∈ (0 , T ] , V V x x x v x Γ(0) = Γ 0 . The phase-field approximation of the above equation - Allen-Cahn:  x , t ) − 1 x , t )) − 1 ∂ t φ ( x x x , t ) = △ φ ( x x ǫ 2 G ′ ( φ ( x x ǫ ( η ( x x x , t ) − λ ( t )) in Ω × (0 , T ] ,   ∇ φ · ν ν ν Ω = 0 on ∂ Ω × (0 , T ] , = φ 0 in Ω .  φ ( · , 0)  Feng Wei Yang Seminar at INI 1 October 2015 28 / 39

Our optimal control model cont. The objective functional: � T J ( φ, η ) = 1 � x + θ � x )) 2 dx x , t ) 2 dx ( φ ( x x , T ) − φ obs ( x η ( x x x x x xdt , x 2 2 Ω 0 Ω 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: x , t ) + ǫ − 2 G ′′ ( φ ( x � ∂ t p ( x x x , t ) = −△ p ( x x x x , t )) p ( x x x , t ) in Ω × [0 , T ) , p ( x x , T ) = φ ( x 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: 512 2 requires 0.4 gigabytes 3-D: 512 3 requires 215 gigabytes Feng Wei Yang Seminar at INI 1 October 2015 31 / 39

Two-grid solution strategy One complete solve for the Allen-Cahn equation from t=(0,T] One time step Start the next η iteration Restrict the converged solution of ϕ Fine grid for the Allen-Cahn equation One complete solve for the adjoint equation from t=[T,0) Intermediate grid(s) Interpolate the computed η One time step Coarse grid for the adjoint equation 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