Hybrid Octree/Quadtree AMR for Anisotropic Domains Tobin Isaac - - PowerPoint PPT Presentation

Hybrid Octree/Quadtree AMR for Anisotropic Domains

Tobin Isaac

Institute for Computational Engineering & Sciences The University of Texas at Austin

SIAM PP 2014 Portland, Oregon February 19, 2014

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 1 / 23

1

Thin-domain Problems

2

Hybrid Octree/Quadtree AMR

3

Example My collaborators on this work are Carsten Burstedde and Omar Ghattas. These slides can be found at www.ices.utexas.edu/~tisaac/slides.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 2 / 23

Three-dimensional models in thin domains are common

Ocean [Kanarska et al., 2007] Subsurface [Maxwell et al., 2008] Ice [Isaac et al.]

Models with 3D implicit subsystems (e.g. pressure Poisson solves) are becoming more common.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 3 / 23

Many of these problems benefit from local mesh refinement

They benefit from local horizontal mesh refinement [Durand et al., 2009]

Ice-shelf/ice-sheet interface location, with AMR (black) and fixed resolutions (gray).

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 4 / 23

Many of these problems benefit from local mesh refinement

They benefit from local vertical mesh refinement

Boundary singularities, e.g. Dirichlet-Neumann transitions. Spatially varying boundary layer requirements, e.g. in ocean modeling, can be treated with AMR instead of/in addition to hybrid vertical coordinates [Pain et al., 2005].

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 5 / 23

Isotropic mesh refinement has drawbacks

The aspect ratio of the discretization is fixed

It cannot adaptively match the length scales of the solution:

affects accuracy; affects numerics, e.g. smoother convergence rate.

Forcing horizontal refinement when only vertical refinement is desired (and vice versa) increases the mesh size by much more than necessary:

e.g. by a factor of 4 for octrees (8 new octants vs. 2 new layers).

Limits available numerical methods, e.g. geometric multigrid with semi-coarsening.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 6 / 23

Isotropic mesh refinement has drawbacks

There is no natural support for column integrity

Many numerical schemes require degrees of freedom (dofs) to be

• rganized into columns:

for solver convergence, e.g. line-smoothers; for high-throughput computing.

Distributing a mesh based on graph-partitioning or a space-filling curve can split a column of dofs between subproblems:

leads to less-compact subdomains, worse convergence of domain decomposition methods [Farhat et al., 1994].

Ω0 Ω1

vs.

Ω0 Ω1

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 7 / 23

Fully anisotropic AMR has drawbacks

Unstructured anisotropic AMR [Loseille et al., 2010]

Necessary for some problems (e.g. supersonic flows), but metric-tensor based adaptivity is complex and is less scalable than structured or tree-based AMR.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 8 / 23

Fully anisotropic AMR has drawbacks

Hierarchical anisotropic AMR [Williamschen and Groth, 2013]

More complex than isotropic counterpart, cannot adapt to arbitrary target metric tensors without mesh movement.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 9 / 23

The hybrid AMR scheme

partition 0 partition 1 A p4est forest-of-quadtrees to manage columns, with each column stored as a flat, linear binary tree of layers, which guarantees column integrity.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 10 / 23

Some good qualities of isotropic AMR are inherited

Well-shaped partitions quickly generated by a space filling curve,

using weighted partitioning algorithm already present in p4est [Burstedde et al., 2011].

In-place, single-sweep vertical refinement and vertical coarsening algorithms.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 11 / 23

Some good qualities of isotropic AMR are lost

Partitioning has coarser granularity

Load-balancing can be poor, but only if there are columns with Ncolumn Ntotal/P layers.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 12 / 23

Some good qualities of isotropic AMR are lost

Horizontal adaptivity is less flexible

Individual layers cannot be horizontally refined: the mesh size grows quickly with horizontal refinement. needs refinement refinement unnecessary Coarsened columns are combined into least-common-ancestor leaves. level 1 level 2 level 3

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 13 / 23

Mesh quality conditions

Neighboring refinement level is the only free quality condition in isotropic hierarchical AMR, hence the emphasis on 2:1 condition algorithms [Isaac et al., 2012]. In hybrid AMR, the layer aspect ratio is also free.

Maximum and minimum ratios can be enforced at the same time as a 2:1 condition: a single mesh conditioning algorithm.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 14 / 23

What are 2:1 conditions for hybrid AMR?

Some situations that seem fine are wrong C0 finite elements with limited hanging-node stencils

The 2:1 condition is satisfied in terms of adjacency, but there is a hanging node () that depends on another hanging node ().

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 15 / 23

What are 2:1 conditions for hybrid AMR?

Some situations that seem fine are wrong C0 finite elements with limited hanging-node stencils

The 2:1 condition is satisfied in terms of adjacency, but there is a hanging node () that depends on another hanging node ().

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 15 / 23

What are 2:1 conditions for hybrid AMR?

Some situations that seem fine are wrong C0 finite elements with limited hanging-node stencils

The 2:1 condition is satisfied in terms of adjacency, but there is a hanging node () that depends on another hanging node ().

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 15 / 23

What are 2:1 conditions for hybrid AMR?

Some situations that seem wrong are fine Cross-refined neighbors

There is no parent/master face, but C0 hanging nodes can be treated as if there were, and face integrals can be calculated as if there were four smaller faces, all using 2:1 interpolation rules.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 16 / 23

What are 2:1 conditions for hybrid AMR?

Some situations that seem wrong are fine Cross-refined neighbors

There is no parent/master face, but C0 hanging nodes can be treated as if there were, and face integrals can be calculated as if there were four smaller faces, all using 2:1 interpolation rules.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 16 / 23

What are 2:1 conditions for hybrid AMR?

Some situations that seem wrong are fine Cross-refined neighbors

There is no parent/master face, but C0 hanging nodes can be treated as if there were, and face integrals can be calculated as if there were four smaller faces, all using 2:1 interpolation rules.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 16 / 23

The full mesh conditioning algorithm

Mesh conditioning (‘‘balance’’)

1 Enforce maximum flatness W/H

by refining columns.

2 Enforce 2:1 between columns.

This does not increase flatness.

3 Enforce minimum flatness by

refining layers.

4 Enforce 2:1 between layers.

This does not decrease flatness, and does not increase maximum flatness. The quadtree/octree ‘‘insulation layer’’ property does not apply. Iterative or all-to-all communication is required in parallel.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 17 / 23

The full mesh conditioning algorithm

Mesh conditioning (‘‘balance’’)

1 Enforce maximum flatness W/H

by refining columns.

2 Enforce 2:1 between columns.

This does not increase flatness.

3 Enforce minimum flatness by

refining layers.

4 Enforce 2:1 between layers.

This does not decrease flatness, and does not increase maximum flatness. The quadtree/octree ‘‘insulation layer’’ property does not apply. Iterative or all-to-all communication is required in parallel.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 17 / 23

The full mesh conditioning algorithm

Mesh conditioning (‘‘balance’’)

1 Enforce maximum flatness W/H

by refining columns.

2 Enforce 2:1 between columns.

This does not increase flatness.

3 Enforce minimum flatness by

refining layers.

4 Enforce 2:1 between layers.

This does not decrease flatness, and does not increase maximum flatness. The quadtree/octree ‘‘insulation layer’’ property does not apply. Iterative or all-to-all communication is required in parallel.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 17 / 23

The full mesh conditioning algorithm

Mesh conditioning (‘‘balance’’)

1 Enforce maximum flatness W/H

by refining columns.

2 Enforce 2:1 between columns.

This does not increase flatness.

3 Enforce minimum flatness by

refining layers.

4 Enforce 2:1 between layers.

This does not decrease flatness, and does not increase maximum flatness. The quadtree/octree ‘‘insulation layer’’ property does not apply. Iterative or all-to-all communication is required in parallel.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 17 / 23

The full mesh conditioning algorithm

Mesh conditioning (‘‘balance’’)

1 Enforce maximum flatness W/H

by refining columns.

2 Enforce 2:1 between columns.

This does not increase flatness.

3 Enforce minimum flatness by

refining layers.

4 Enforce 2:1 between layers.

This does not decrease flatness, and does not increase maximum flatness. The quadtree/octree ‘‘insulation layer’’ property does not apply. Iterative or all-to-all communication is required in parallel.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 17 / 23

The full mesh conditioning algorithm

Mesh conditioning (‘‘balance’’)

1 Enforce maximum flatness W/H

by refining columns.

2 Enforce 2:1 between columns.

This does not increase flatness.

3 Enforce minimum flatness by

refining layers.

4 Enforce 2:1 between layers.

This does not decrease flatness, and does not increase maximum flatness. The quadtree/octree ‘‘insulation layer’’ property does not apply. Iterative or all-to-all communication is required in parallel.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 17 / 23

The full mesh conditioning algorithm

Mesh conditioning (‘‘balance’’)

1 Enforce maximum flatness W/H

by refining columns.

2 Enforce 2:1 between columns.

This does not increase flatness.

3 Enforce minimum flatness by

refining layers.

4 Enforce 2:1 between layers.

This does not decrease flatness, and does not increase maximum flatness. The quadtree/octree ‘‘insulation layer’’ property does not apply. Iterative or all-to-all communication is required in parallel.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 17 / 23

The full mesh conditioning algorithm

Mesh conditioning (‘‘balance’’)

1 Enforce maximum flatness W/H

by refining columns.

2 Enforce 2:1 between columns.

This does not increase flatness.

3 Enforce minimum flatness by

refining layers.

4 Enforce 2:1 between layers.

This does not decrease flatness, and does not increase maximum flatness. The quadtree/octree ‘‘insulation layer’’ property does not apply. Iterative or all-to-all communication is required in parallel.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 17 / 23

The full mesh conditioning algorithm

Mesh conditioning (‘‘balance’’)

1 Enforce maximum flatness W/H

by refining columns.

2 Enforce 2:1 between columns.

This does not increase flatness.

3 Enforce minimum flatness by

refining layers.

4 Enforce 2:1 between layers.

This does not decrease flatness, and does not increase maximum flatness. The quadtree/octree ‘‘insulation layer’’ property does not apply. Iterative or all-to-all communication is required in parallel.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 17 / 23

Implementation and availability

Almost the full forest-of-quadtree/octree API is reproduced. method additional communication refine +1 allgather coarsen +1 allgather partition +1 allgather, +1 isend/irecv balance +? allgather, +? isend/irecv lnodes (FE mesh) +1 allgather, +1 isend/irecv ghost +1 isend/irecv load/save +1 parallel read/write vtk Available as example/p6est (halfway between quadtree and octree) in branch p6est on github.com/tisaac/p4est.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 18 / 23

Example Problem

Stokes flow in a thin, periodic domain

Stress-free top surface and normal Dirichlet / tangential Robin-type boundary conditions on the bottom surface, partitioned across P = 7 processes. Octree partitioning Hybrid partitioning

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 19 / 23

Stokes flow in a thin domain

Krylov method convergence

The Stokes preconditioner has some block-Jacobi (non-overlapping DD)

• components. As the domain gets thinner, dofs become more tightly coupled to

vertical neighbors, hence the importance of column integrity.

20 40 60 80 100 10−16 10−8 100 Krylov iteration

Octree partitioning: rk/r0

20 40 60 80 100 10−16 10−8 100 Krylov iteration

Hybrid partitioning: rk/r0

H :W = 1:1 H :W = 1:10 H :W = 1:100

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 20 / 23

References I

Carsten Burstedde, Lucas C. Wilcox, and Omar Ghattas. p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees. SIAM Journal on Scientific Computing, 33(3):1103--1133, 2011. doi: 10.1137/100791634. Gael Durand, Olivier Gagliardini, Thomas Zwinger, Emmanuel Le Meur, and Richard CA Hindmarsh. Full stokes modeling of marine ice sheets: influence of the grid size. Annals of glaciology, 50(52):109--114, 2009. Charbel Farhat, Jan Mandel, and Francois Xavier Roux. Optimal convergence properties of the feti domain decomposition method. Computer methods in applied mechanics and engineering, 115(3):365--385, 1994. Tobin Isaac, Georg Stadler, and Omar Ghattas. Solvers for higher-order finite elements on nonconforming and anisotropic meshes, applied to ice sheet

• dynamics. In preparation.
• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 21 / 23

References II

Tobin Isaac, Carsten Burstedde, and Omar Ghattas. Low-cost parallel algorithms for 2:1 octree balance. In Proceedings of the 26th IEEE International Parallel & Distributed Processing Symposium. IEEE, 2012. Y Kanarska, A Shchepetkin, and JC McWilliams. Algorithm for non-hydrostatic dynamics in the regional oceanic modeling system. Ocean Modelling, 18(3): 143--174, 2007. Adrien Loseille, Alain Dervieux, and Frédéric Alauzet. Fully anisotropic goal-oriented mesh adaptation for 3d steady euler equations. Journal of Computational Physics, 229(8):2866--2897, 2010. Reed M Maxwell, Steven F Carle, and Andrew FB Tompson. Contamination, risk, and heterogeneity: on the effectiveness of aquifer remediation. Environmental Geology, 54(8):1771--1786, 2008. CC Pain, MD Piggott, AJH Goddard, F Fang, GJ Gorman, DP Marshall, MD Eaton, PW Power, and CRE De Oliveira. Three-dimensional unstructured mesh ocean

• modelling. Ocean Modelling, 10(1):5--33, 2005.
• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 22 / 23

References III

MJ Williamschen and CPT Groth. Parallel Anisotropic Block-Based Adaptive Mesh Refinement Algorithm For Three-Dimensional Flows. PhD thesis, University of Toronto, 2013.

• T. Isaac (ICES, UT Austin)

p6est SIAM PP 2014 23 / 23