Recent Developments in High-Resolution and Adaptive Methods for CFD

1. Recent Developments in High-Resolution and Adaptive Methods for CFD • Joint work with: • Terry Ligocki, Dan Martin, Peter McCorquodale, David Serafini (ANAG / LBNL). • Robert Crockett, Robert Fisher, Richard Klein, Chris McKee (UC Berkeley). • Steve Jardin, Ravi Samtaney (PPPL). • Scott Baden, Greg Balls (UCSD). • Caroline Bono (LBNL / LLNL) • Alex Friedman, Dave Grote, Jean-Luc Vay, Robert Ryne (AFRD / LBNL). Phillip Colella, Applied Numerical Algorithms Group Lawrence Berkeley National Laboratory Funding: US Department of Energy Applied Mathematical Sciences Program, SciDAC Program NASA Computational Technologies Program

2. Topics New methods of potential interest to astrophysical CFD • Godunov methods for MHD • Splitting methods for gravitationally stratified flows • Fast scalable Poisson solvers • AMR for PIC

3. Divergence-Free Constraint for the MHD Equations , The divergence-free constraint is an initial-value constraint: • Computational physicists have been very concerned about the extent to which • This constraint is satisfied numerically. Attempts to deal with this include: • Staggered grids peculiar advection schemes; also a nuisance in AMR, particularly with nonideal effects. • Work in terms of potentials - starting from scratch relative to existing hyperbolic machinery. • Godunov symmetrization (“eight-wave methods”). • Projections.

4. What’s really going on ? (1) Modified equation analysis says that the constraint is not really preserved: In particular, (2) Without the constraint being satisfied, the modified equations are ill-posed: the linearized coefficient matrices can have eigenvector deficiencies.

5. Let’s consider the special case of a isolated propagating plane wave. If we use as initial data where is a right eigenvector of the linearized system, then, to leading order in where is the corresponding eigenvalue. The divergence-free condition for such a plane wave reduces to The equation for the amplitude is

6. If we include the effect of numerical errors in a modified equation, we obtain i.e. the solution has the same regularity in space as the forcing terms. and we lose one derivative. This is what happens if the linearized coefficient matrix of a hyperbolic system has an eigenvector deficiency.

7. Regularization Godunov symmetrization: add a nonlinear multiple of the constraint to the equations in a way that eliminates the deficiencies, while retaining the same smooth solutions if the constraint is satisfied. Parabolic regularization smoothes the deficient field: Divergence filter (Marder): This is not the same as adding a diffusion term to the magnetic field evolution. If the divergence-free condition is satisfied, the regularizing term is identically zero. This issue arises in many other problems with “gauge constraints”: e.g. solid mechanics, Hamilton-Jacobi, electromagnetics, incompressible flow.

8. Evaluating the Possible Remedies We developed an unsplit second-order Godunov method for MHD, and used several of the standard approaches to enforce the constraint. For plane wave solutions, standard co-located finite difference methods without regularization are unstable. Inclusion of a divergence filter stabilizes these methods. For problems with discontinuities, must supplement filter with a more elaborate form of regularization: face-centered magnetic fields used to compute conservative fluxes must be made divergence-free using a projection (true for seven-wave and eight-wave formulations). Often objected to because it turns a local problem into a nonlocal problem. Upwinding + nonconservative terms in eight-wave formulation can lead to loss of accuracy for linear waves nearly parallel to one of the mesh directions - lack of cancellation of errors when wave speed changes sign. Similar failure mechanism in flux tube problems ? How do these considerations effect the staggered methods ?

9. Applications

10. CFD for Gravitationally Stratified Flows Time scales: compressive (acoustic waves), bouyant (gravity waves), advective. Often only interested in the advective scales, which are slower than the other two. Goal: split the equations so that the system is block-diagonalized with respect to the three scales. We will still be solving the full Euler equations.

11. Anelastic Splitting of the Dynamics , Vortical dynamics: velocity, pressure satisfy anelastic equations. , Compressible dynamics: velocity, pressure satisfy wave equation. , Bouyant dynamics: horizontal response to vertical displacements. , ,

12. Anelastic Splitting of the Dynamics Appropriate temporal discretization of each of the components, with explicit predictor-corrector to evaluate the off-diagonal terms. Vortical dynamics: Semi-implicit method, with explicit treatment of advection. , Compressible dynamics: Implicit temporal discretization, or subcycling in time. , Gravity waves: ?

13. Gravity Waves in a Continuously Stratified Medium Asymptotic analysis: aspect ratio << 1; Fr, Ma << 1, Fr ~ Ma , is a self-adjoint second-order differential operator in the vertical direction. Its eigenmodes define an infinite collection of horizontally propagating waves, with wave speeds , Nontrivial coupling to potential velocity:

14. Gravity Waves in a Continuously Stratified Medium • Semi-implicit treatment of gravity waves: • Evolve for one time step using split method described in earlier slide. • Extract fast mode(s) for , intial data, and evolve them for one time step using a stable method (subcycling in time, implicit method). • Correct at the old and new time, and use it to recompute the contribution to the update of the velocities.

15. Status and Plans • AMR for low Mach number, gravitationally stratified flows • Well-posed initial-boundary value problem • Natural coupling to higher-speed case. • Mixture of 2D and 3D simulations in an AMR setting. • Anisotropic solvers.

16. Analysis-Based Poisson Solvers , Idea: disjoint regions in space are decoupled, modulo analytic functions. Domain decomposition should lead to efficient parallel solvers. • Multigrid: localizes computation, but not communication. • Schwartz domain decomposition: still iterative. • Fast multipole method: localizes computation and communication noniteratively, but cost goes up with the number of dimensions. Real analytic, with rapidly convergent Taylor expansion

17. Method of Local Corrections (Anderson, 1986) MLC is a non-iterative domain decomposition method for computing ( ) (1) Solve local problems on overlapping local patches: (2) Solve a single coarse-grained problem to represent the nonlocal coupling: (3) Compute composite solution as combination of local fields and interpolated corrected global field:

18. Method of Local Corrections Why does this work ? (1) Local problems are independent, and trivially parallel: outside the support of the local charge introduces a small (O(Hp) ) (2) Truncating error. The global solve is a bottleneck, but can be made small relative to the local solutions. (3) The local interactions / local corrections step requires mainly local communications. The correction of the coarse-grid values is done in a way so that the interpolated values are discrete harmonic.

19. Method of Local Corrections • Anderson’s original algorithm was a PPPM method for vortex transport in • 2D - local calculation was a N-body method, coarse calculation a Mehrstellen finite difference method. Baden (1986) used the method as the basis of parallel domain-decomposition algorithm for particle methods. • What is needed to apply this idea to 3D multiresolution gridded calculations? • Efficient grid-based infinite domain solver. • Mehrstellen discretizations play an essential role: • High-order accuracy in regions where the solution is harmonic. • Ref: Balls and Colella, 2002; Baden, Balls, Colella, McCorquodale, 2005.

20. James’ Algorithm for Infinite-Domain Boundary Conditions In 3D, the direct calculation of the surface-surface potential is too expensive (O(N4)). We reduce this to a small fraction of the solver cost using a simplified multipole expansion. The resulting method is 3X faster than domain-doubling method, no self-force issues.

21. Performance of Fast James’ Algorithm In 3D, the direct calculation of the surface-surface potential is too expensive (O(N4)). We reduce this to a small fraction of the solver cost using a simplified multipole expansion. The resulting method is 3X faster than domain-doubling method, no self-force issues.

22. MLC Results 3D MLC calculation of Poisson's equation, with 256^3 mesh broken into 64^3 blocks. Image on right shows solution error through slice at mid-plane.

23. Status and Plans Finish final performance squeeze of two-level MLC. • Eliminate superfluous communications. • Parallelize global solve. • Optimize degree of overlap. Multilevel MLC. • “Non-iterative multigrid” - communication comparable to a single AMR V-cycle. • Fewer terms in multipole expansion - infinite domain solves less expensive. Other problems: heat equation in 3D; cell-centered solvers; nontrivial boundary conditions … Goal: AMR Poisson solvers at cost per grid point = 5X uniform grid FFT.

24. AMR for Electrostatic PIC • Natural adaptive strategy: refine to keep the number of particles per cell fixed, or to resolve large gradients near a boundary. • Progress to date: • Node-centered AMR Poisson solvers, Shortley-Weller treatment of irregular boundaries. • Data structure support for bin-sorted particles on an adaptive grid; PIC interpolation between particles and grid. • Coupling of AMR-PIC to beam dynamics codes. Coupling to QuickPIC in progress.

25. Matching Conditions at Coarse / Fine Boundaries One-way coupling: • Solve on • Solve on , using values interpolated from coarse calculation • Difficulty: Poisson solve delocalizes the error, limiting the accuracy in some cases to that of the coarse grid ( ). Two-way coupling: Well-posed free-boundary problem for Poisson; solution error has correct locality properties.

26. Spurious Self-Forces Due to Refinement Boundaries Possible solutions: • Keep particles away from refinement boundaries. • Use one-way coupling (only coarse Neumann matching at boundary). • Find different particle / grid interaction that does not generate spurious images (cell centering + ?). • Find a different multiresolution Poisson solver that does not generate spurious images.

27. Status and Plans Two-level MLC PIC code implemented; currently testing in adaptive mode. AMR-PIC using multilevel MLC algorithm. New deposition / interpolation algorithms for nodal-point AMR. Open-source AMR cosmology code being developed by Francesco Miniati (compressible gas dynamics and collisionless dark matter)

28. Other Work Open source software: http://davis.lbl.gov/APDEC/software.html http://seesar.lbl.gov/ANAG Volume-of-fluid discretization of moving boundaries - application to low-Mach number flame tracking in supernovae. Low Mach number flames in type 1A SN - LBNL-CCSE / UCSC