High-order methods for computational astrophysics

Abstract

In computational fluid dynamics, high-order numerical methods have gained quite popularity in the last years due to the need of high fidelity predictions in the simulations. High-order methods are suitable for unsteady flow problems and long-term simulations because they are more efficient when obtaining higher accuracy than low-order methods, and because of their outstanding dissipation and dispersion properties. In the present work, the development and application of three high-order numerical methods, namely, the conservative finite difference (FD) method, the finite volume (FV) method, and the discontinuous Galerkin spectral element method (DGSEM), is presented. These methods are used here for solving three equations systems arising in computational astrophysics on flat spacetimes, specifically, the ideal magnetohydrodynamics (MHD), relativistic hydrodynamics (SRHD) and relativistic magnetohydrodynamics (SRMHD). Our computational framework has been subject to the standard testbench in computational astrophysics. Numerical results of problems having smooth flows, and problems with shock-dominated flows, are also reported.

Finite volume methods are numerical methods based on the weak solution of conservation laws in integral form. Unlike finite volume methods, where cell averages of the solution are evolved in time, in the conservative finite difference schemes only the solution at specific nodal points are considered. This difference offers a high efficiency of finite difference over finite volume methods in two and three dimensional high-order calculations because of the form of the utilized stencils in the reconstruction step. Recently, a lot of effort has been put into the development of efficient high-order accurate reconstruction procedures on structured and unstructured meshes. The most widely used procedure to achieve high-order spatial accuracy in finite volume and conservative finite difference methods is the WENO reconstruction. The basic idea of the WENO schemes is based on an adaptive reconstruction procedure to obtain a higher-order approximation on smooth regions while the scheme remains non-oscillatory near discontinuities. For this reason, the WENO formulation is particularly effective when solving conservation laws containing discontinuities. In this work, the FD and FV methods are extended to very high-order accuracy on regular Cartesian meshes by making use of the arbitrary high-order reconstruction WENO operator. The time discretization is carried out with a strong stability-preserving Runge-Kutta (SSPRK) method. The MHD, SRHD and SRMHD equations are then solved with these two methods for problems having strong shock configurations.

The discontinuous Galerkin (DG) methods combine the ideas of the finite element (FE) and the finite volume methods. From the FE methods, the solution and test functions in the variational formulation of the conservation law are locally represented by polynomials, allowing to be discontinuous at element faces. In order to stabilize the scheme, from the FV methods are borrowed the ideas of using Riemann solvers, which permit to connect a given element with its direct neighboring ones. One special case in the family of DG methods is the DGSEM. In these methods, the domain is decomposed into quadrilateral/hexahedral elements, and the solution and the fluxes are represented by tensor-product basis functions (high-order Lagrangian interpolants). The integrals are approximated by quadrature, and the nodal points, where the solution is computed, are the Gauss-Legendre quadrature points. With these choices, the DG operator has a dimension-by-dimension splitting form, which yields more efficiency due to less operations and less memory consumption. In this work, the DGSEM has been also extended to the equations of computational astrophysics on flat spacetimes, but restricted only to the MHD and SRHD equations.

Because discontinuous solutions form part of the nature of the hyperbolic conservation laws, shock capturing strategies have to be devised, especially for the discontinuous Galerkin method. For the DGSEM, a hybrid DG/FV shock capturing approach is used as the main building block for stabilization of the solution when shocks take place. The hybrid DGSEM/FV is constructed in such a way that, in regions of smooth flows, the DGSEM method is employed, and those parts of the flow having shocks, the DGSEM elements are interpreted as quadrilateral/hexahedral subdomains. In each of these subdomains, the nodal DG solution values are used to build a new local domain composed now of finite volume subcells, which are evolved with a robust finite volume method with third order WENO reconstruction. This new numerical framework for computational astrophysics based on the hybridization of high-order methods brings very promising results.

In der numerischen Strömungsmechanik sind numerische Methoden hoher Ordnung, aufgrund der Notwendigkeit hochgenauer Vorhersagen in den Simulationen, in den letzten Jahren immer populärer geworden. Verfahren hoher Ordnung eignen sich für instationäre Strömungsprobleme und langfristige Simulationen, weil sie beim Erreichen höherer Genauigkeit effizienter sind als Methoden niedriger Ordnung und hervorragende Dissipations- und Dispersionseigenschaften besitzen. In dieser Arbeit wird die Entwicklung und Anwendung dreier numerischer Verfahren hoher Ordnung vorgestellt: der konservativen Finite-Differenzen-Methode (FD), der Finite-Volumen-Methode (FV) und der Discontinuous Galerkin-Spektral-Element-Methode (DGSEM). Diese Methoden werden zur Lösung von drei Gleichungssystemen aus der Computer-Astrophysik auf flachen Raumzeiten verwendet, nämlich der Magnetohydrodynamik (MHD), der relativistischen Hydrodynamik (SRHD) und der relativistischen Magnetohydrodynamik (SRMHD). Unser Berechnungsverfahren wurde den Standard-Testfällen der Computer-Astrophysik unterzogen. Numerische Berechnungen von Problemen mit glatten Strömungsverläufen und Problemen mit starken Stößen werden ebenfalls angegeben.

Finite-Volumen-Verfahren sind numerische Verfahren, die auf der schwachen Lösung der Erhaltungsgleichungen in integraler Form basieren. Im Gegensatz zu den Finite-Volumen-Verfahren, in denen die Mittelwerte der Lösung in der Zeit entwickelt werden, wird in den konservativen Finite-Differenzen-Verfahren nur die Lösung an bestimmten Knotenpunkten berücksichtigt. Dieser Unterschied führt, aufgrund der Form der im Rekonstruktionsschritt verwendeten Stencils, bei zwei- und dreidimensionalen Berechnungen hoher Ordnung zu einer höheren Effizienz der Finite-Differenzen-Verfahren gegenüber den Finite-Volumen-Verfahren. Das am weitesten verbreitete Verfahren, um räumliche Genauigkeit hoher Ordnung in Finite-Volumen- und konservativen Finite-Differenzen-Methoden zu erreichen ist die WENO-Rekonstruktion. In dieser Arbeit werden die FD- und FV-Methoden durch die Nutzung des WENO-Rekonstruktions-Operators beliebiger Ordnung auf regulären kartesischen Gittern zu einer Genauigkeit sehr hoher Ordnung erweitert. Die Zeitdiskretisierung erfolgt mit einer "Strong Stability-Preserving Runge-Kutta"-Methode (SSPRK). Die MHD-, SRHD- und SRMHD-Gleichungen werden mit diesen beiden Methoden vollständig gelöst und es werden, zum ersten Mal für diese Gleichungen, numerische Berechnungen von Problemen mit Stößen angegeben.

Die Discontinuous Galerkin-Methode (DG) kombiniert die Ideen der Finite-Elemente- (FE) und der Finite-Volumen-Methoden. Wie im FE-Verfahren werden die Lösung und die Testfunktionen in der Variationsformulierung der Erhaltungsgleichung durch Polynome dargestellt, aber mit Zulassung von Unstetigkeiten an Elementgrenzen. Zur Stabilisierung des Schemas, wird von den FV-Methoden die Verwendung von Riemann-Lösern entliehen, welche die Verbindung eines gegebenen Elements mit seinen direkten Nachbarn ermöglichen. Ein Sonderfall in der Familie der DG-Methoden ist die DGSEM. Bei diesen Verfahren wird das Rechengebiet in Vierecks/Hexaeder-Elemente zerlegt, während die Lösung und die Flüsse von Tensorprodukt-Basisfunktionen (Lagrange-Interpolanten hoher Ordnung) dargestellt werden. Die Integrale werden durch Quadratur approximiert und die Knotenpunkte, an denen die Lösung berechnet wird, sind Gauß-Legendre-Quadraturpunkte. Durch diese Wahl hat der DG-Operator eine "Dimension-by-Dimension"-Form, die mehr Effizienz durch weniger Operationen und weniger Speicherbedarf ergibt. In dieser Arbeit wurde die DGSEM auch für die Gleichungen der Computer-Astrophysik auf flachen Raumzeiten erweitert, aber mit Einschränkung auf die MHD- und die SRHD-Gleichungen.

Da diskontinuierliche Lösungen ein Teil der Natur von hyperbolischen Erhaltungsgleichungen sind, müssen Shock-Capturing-Strategien entwickelt werden, insbesondere für das Discontinuous Galerkin-Verfahren. Für die DGSEM wird ein hybrider DG/FV-Ansatz als Hauptbaustein zur Stabilisierung der Lösung verwendet, wenn Stöße auftreten. Das Hybrid DGSEM/FV ist so aufgebaut, dass in Bereichen mit glattem Strömungsverlauf das DGSEM-Verfahren verwendet wird, während in Gebieten mit Stößen die DGSEM-Elemente als Vierecks/Hexahedral-Untergebiete interpretiert werden. In jedem dieser Untergebiete werden die nodalen DG-Lösungswerte verwendet, um ein neues lokales Gebiet aus Finite-Volumen-Unterzellen zu bauen, die mit einem robusten Finite-Volumen-Verfahren mit WENO-Rekonstruktion dritter Ordnung entwickelt werden. Dieses neue numerische Gerüst für die Computer-Astrophysik, das auf der Hybridisierung von Verfahren hoher Ordnung basiert, bringt sehr vielversprechende Ergebnisse.