"Terascale Simulation Initiative" PI T. Messacappa, Application Point of Contact: B. Messer (ORNL), TSTT point of Contact: E. D'Azevedo (ORNL)

"Center for Magnetic Reconnection" PI Bhattacharjee, Application Point of Contact: B. Rosner (UofC), TSTT point of Contact: M. Shephard (RPI)

TSTT researchers at ORNL and RPI have active collaborations with two astrophysics SciDAC centers: the Terascale Supernova Initiative (TSI) led by Mezzacappa and the Magnetic Reconnection Center led by Bhattacharjee. The primary focus in both cases is the exploration of alternative discretization strategies for various aspects of the solution process including the neutrino transport problem and hydrodynamics simulations. In addition, we have been exploring a 3D spatial caching strategy to reduce costly evaluations of the scattering kernels for Boltzman transport.

Discontinuous Galerkin Discretizations for Radiative Transfer

TSTT Personnel: Valmor de Almeida (ORNL), Ed D’Azevedo (ORNL)

TSI Personnel: Tony Mezzacappa (ONRL), Bronson Messer (ORNL), Matthias Liebendoerfer (ORNL)

With the TSI team, TSTT researchers at ORNL are exploring new schemes for solving models of radiative transfer in stellar atmospheres representative of supernova type 2 core collapse/explosions. The current solution method, namely, discrete ordinates, in conjunction with finite differences, may require petascale computing power when extended to three dimensions. The TSTT/TSI interaction is aimed at exploring novel schemes that are able to keep computational requirements at the terascale level.

A series of radiative transfer prototype problems have been chosen by the TSI group. These problems are simplified forms of Boltzmann's integro-differential equation for radiative transfer and thus have features significant to supernova collapse/explosion modeling.

The first problem, Milne's axis-symmetric pure scattering of photons in extended stellar atmospheres, has been studied successfully. In this problem the dependent variable is the radiative intensity field (photon distribution function), and the independent variables are the radius of the spherical stellar atmosphere and the polar angle made by the radius and the photon velocity vector. The new method of analysis consists of a discontinuous Galerkin (DG) approach for the approximation of the radiative intensity field in phase space (the Cartesian product of the one-dimensional radius domain and the polar angle domain). The integro-differential equation is treated as a linear hyperbolic PDE at each step of an iterative scheme that corrects the value of the integral scattering kernel. Each iteration is accelerated by evaluating the kernel from the first and zeroth moments of the original problem. The DG method applied to the linear hyperbolic problem at each iteration was implemented explicitly in phase space since the characteristic curves are known a priori. This allows the resulting linear algebraic set of equations to be solved locally and progressively along wave fronts normal to the characteristics. The DG method developed for the Milne's problem was found to be memory efficient (requiring one order of magnitude less storage than the discrete ordinate method) and fast (one order of magnitude faster than DOM). In addition, in view of the adaptive mesh used, the results were significantly more accurate when capturing important features of the solution. For instance, the region of transition from diffusive radiation to streaming was correctly captured, and the outward peaking (a delta function like sharp increase) of the radiation field at the outer boundary of the extended atmosphere was correctly reproduced. The results were obtained for a range of atmosphere sizes including the range of interest to supernova applications. Although the results are exciting, supernova models are significantly more complex, involving the extended set of equations that couple neutrino radiation transfer to relativistic flow, that is radiative hydrodynamics.

Discontinuous Galerkin Discretizations for Hydrodynamics

TSTT Personnel: Mark Shephard (RPI), Joe Flaherty (RPI), Jean-Francios Remacle (RPI)

Astrophysics Personnel: Robert Rosner (UofC), Greg Weirs (UofC), Andrew Seigel (UofC)

The TSTT team at SCOREC (RPI) is working with the astrophysicists at the University of Chicago to develop hp-adaptive discontinuous Galerkin discretizations for flow problems of importance to astrophysics applications. This work is part of a long-standing collaboration with the FLASH center at the UofC and is currently funded by both the ASCI ASAP program and SciDAC. The work outlined in this section was funded predominantly by the ASCI ASAP program, but lays the groundwork for new research that will be funded by SciDAC.

Together with the FLASH team, we determined a set of test problems that would provide a basis of comparison between h-p adaptive DG methods and the currently used Piecewise Parabolic method. With the start of the SciDAC program, we executed these test problems in the SCORE Trellis framework—a flexible, extensible framework for PDE-based application solution. We ran the tests on both uniform meshes and adaptive meshes to demonstrate the increased efficiency over fixed grids.

Based on these and other studies, the h-adaptive techniques in the DG code have been extended in terms of capability and efficiency. Both nonconforming isotropic and conforming refinement strategies are now supported and the conforming refinement includes an initial anisotropic mesh adaptation using an anisotropic error indicator. These tests also showed that the effective application of p-adaptivity requires additional developments in two areas. The first is the appropriate adaptive selection of polynomial order. The second and more difficult task is the construction of a reliable high-order limiting method which works in the general 3D case. This is particularly crucial for hyperbolic equations such as the Euler equations because of shocks and discontinuities. Limiting the solution at the vicinity of discontinuities is always required but usual limiting techniques have the major drawback of reducing the accuracy of the solution in regions where a limiter is not be required. We developed an error estimator for DG based on superconvergence of the solution (in average) at outflow boundaries. Using that important result, we developed a robust shock detector that enables us to turn on/off the limiter procedure when needed. Figure 2 shows an example of the computation of vortex sheets using the detector.

Figure 2

Figure 2 Density contours (left) and adapted mesh (right) for the vortex sheet problem at t=1.

Efforts are now being initiated to integrate DG discretizations, local time stepping and a TSTT mesh interface into the FLASH code. In previous work, a distributed memory, explicit local time stepping technique was developed as the integration scheme used with the DG discretization. This scheme allows us to satisfy the Courant (CFL) stability conditionlocally rather than globally and efficiencies are often 50 to 100 times those of methods using a single time step in a method of lines formulation. This scheme provides accurate resolution near shocks and other discontinuities. It is also useful for resolving layers in viscous problems with small diffusion.

The DG methods extend to viscous problems while conserving the higher order of accuracy, and we think this method is worthwhile for such a class of problems. We are planning to consider an implicit in time integration scheme to enhance the capabilities of the code to solve problems of interest to the FLASH project, such as wave breaking problems or double diffusion. Our experience in using TOPS solvers such as PETSc will help us to provide such capabilities. For DG, a matrix free GMRES type method might be suitable and also very efficient. The compactness of the DG stencil might permit some very simple block diagonal preconditionings (e.g., local Jacobi methods). Using the DG method for implicit time integration is another very interesting possibility that we will explore.

Caching Strategies for Boltzman Transport

TSTT Personnel: Ed D Azevedo (ORNL)

TSI Personnel: Bronson Messer (ORNL)

Another interaction between TSTT and TSI is the use of a 3D spatial cache for avoiding redundant and costly evaluations of scattering kernels for Boltzmann transport. The scattering kernel is a function of density, temperature and species fraction. We have generated trace data over a realistic TSI computation and visualized the trajectories of evaluation in 3D phase space (density, temperature, species fraction). Initial analysis shows 3D interpolation with a small global cache has good hit ratio of over 90%, or about 9 of 10 evaluations of scattering kernels can be computed by 3D interpolation instead of the actual costly computations.