Trixi.jl is a numerical simulation framework for adaptive, high-order discretizations of conservation laws. It has a modular architecture that allows users to easily extend its functionality and was designed to be useful to experienced researchers and new users alike. In this talk, we give an overview of Trixi’s current features, present a typical workflow for creating and running a simulation, and show how to add new capabilities for your own research projects.
Abstract: Regarding accuracy and stability of numerical schemes for computational fluid dynamics, the investigation of diffusion/dispersion errors depending on the wave number is of utmost importance. Especially for high order methods, a desired small numerical dissipation competes with robustness and thus has to be carefully analyzed. This wave propagation analysis is often based on pure advection problems. In the literature, various approaches to discretize diffusion terms within a DG scheme have been introduced since the discretization of higher order spatial derivatives within the DG framework is less natural than in case of first order derivatives. In this talk, we will address significant differences in the disspation/dispersion properties for linear advection-diffusion, depending on the specific DG viscous flux discretization which is employed. In addition, results on energy stability of DG viscous flux formulations are dealt with and we show how to formulate the well-known LDG and BR1 fluxes in terms of global upwind SBP operators which complements the derivation and analysis regarding element level SBP properties of the DG scheme.
Abstract: Split form schemes for Euler and Navier-Stokes equations are useful for computation of turbulent flows due to their better robustness. This is because they satisfy additional conservation properties of the governing equations like kinetic energy preservation leading to a reduction in aliasing errors at high orders. Recently, linear stability issues have been pointed out for these schemes for a density wave problem and we investigate this behaviour for some standard split forms. By deriving linearized equations of split form schemes, we show that most existing schemes do not satisfy a perturbation energy equation that holds at the continuous level. A simple modification to the energy flux of some existing schemes is shown to yield a scheme that is consistent with the energy perturbation equation. Numerical tests are given using a discontinuous Galerkin method to demonstrate these results.
Abstract: In this talk, we present a state redistribution method for high order discontinuous Galerkin methods on curvilinear embedded boundary grids. State redistribution relaxes the overly restrictive CFL condition that results from arbitrarily small cut cells when explicit time steppers are used. Thus, the scheme can take time steps that are proportional to the size of cells in the background grid. The discontinuous Galerkin scheme is stabilized by postprocessing the numerical solution after each stage or step of an explicit time stepping method. The advantage of this approach is that it uses only basic mesh information that is already available in many cut cell codes and does not require complex geometric manipulations. We prove that state redistribution is conservative and p-exact. Finally, we solve a number of test problems that demonstrate the encouraging potential of this technique for applications on curvilinear embedded geometries. Numerical experiments reveal that our scheme converges with order $p+1$ in $L_1$ and between $p$ and $p+1$ in $L_\infty$.
Speaker: Dr. Romain Veltz, INRIA, France Date: Thursday, 18th March 2021, 10am (CET) Zoom-Link: Please request via email from Michael Schlottke-Lakemper
Abstract
In this talk, I will first present the basics of bifurcation theory. Then, I will give a panorama of BifurcationKit.jl, a Julia package to perform numerical bifurcation analysis of large dimensional equations (PDE, nonlocal equations, etc) possibly on GPUs using Matrix-Free / Sparse Matrix formulations of the problem. Julia programming language gives access to a rich ecosystem (PDE, GPU, AD, cluster…). Notably, numerical bifurcation analysis can be done entirely on GPU as will be shown in an example.
BifurcationKit incorporates continuation algorithms (PALC, deflated continuation, …) which can be used to perform fully automatic bifurcation diagram computation of stationary states. I will showcase this with the 2d Bratu problem.
Additionally, by leveraging on the above methods, the package can also seek for periodic orbits of Cauchy problems by casting them into an equation of high dimension. It is by now, one of the only software which provides parallel (Standard / Poincaré) shooting methods and finite differences based methods to compute periodic orbits in high dimensions. I will present an application highlighting the ability to fine tune BifurcationKit to get performance.
In a last part, I will describe a mean field model of stochastic spiking neurons described with a 2d measure valued equation. I will present a numerical scheme based on an implicit Finite Volume method. I will then provide some mathematical properties of the mean field concerning well posedness and stationary solutions. Additionally, I will show how BifurcationKit.jl can be used to study numerically the model. Finally, I will conclude on open problems, some of which could hopefully be tackled numerically with Trixi.jl.
Julia has been touted as a programming language especially
well-suited for numerical analysis and scientific computing. However,
while its prevalence is steadily increasing, it has not yet seen
widespread adoption in the computational science or high-performance
computing communities. One of the hurdles is a (perceived) lack of
real-world examples that show how Julia can be used to conduct numerical
simulations and what its advantages and drawbacks are for scientific
applications.
To remediate this, in this talk we discuss the development of a
purely hyperbolic method for self-gravitating gas dynamics within our
Julia-based open source simulation framework Trixi.jl (https://github.com/trixi-framework/Trixi.jl).
In this approach, we reformulate the elliptic gravity problem into a
hyperbolic diffusion problem, which is solved in pseudotime using the
same explicit high-order discontinuous Galerkin method we use for the
flow solution. A key benefit is that in the resulting multi-physics
simulation problem, we can reuse existing hyperbolic solvers while
retaining advanced features such as non-conforming and solution-adaptive
meshes.
Next to presenting numerical results, we will critically examine
our experience with building a multi-physics simulation framework with
Julia. We will discuss its strengths and weaknesses as a programming
language for research software engineering, including an assessment of
how Julia’s claimed benefits hold up against scientific reality, and
give a live demonstration of Julia and Trixi.jl in action.
To make the shown examples reproducible by the audience, the Jupyter notebook used for the live demonstration is available at https://github.com/trixi-framework/talk-2021-julia-adaptive-multi-physics-simulations.
It can be either run from a local Julia/Jupyter installation or in the
cloud via Binder (without having to install Julia locally).
Recently, it was discovered that the entropy-conserving/dissipative high-order split-form discontinuous Galerkin discretizations have robustness issues when trying to solve the simple density wave propagation example for the compressible Euler equations. The issue is related to missing local linear stability, i.e. the stability of the discretization towards perturbations added to a stable base flow. This is strongly related to an anti-diffusion mechanism, that is inherent in entropy-conserving two-point fluxes, which are a key ingredient for the high-order discontinuous Galerkin extension. In this paper, we investigate if pressure equilibrium preservation is a remedy to these recently found local linear stability issues of entropy-conservative/dissipative high-order split-form discontinuous Galerkin methods for the compressible Euler equations. Pressure equilibrium preservation describes the property of a discretization to keep pressure and velocity constant for pure density wave propagation. We present the full theoretical derivation, analysis, and show corresponding numerical results to underline our findings. The source code to reproduce all numerical experiments presented in this article is available online (DOI: 10.5281/zenodo.4054366).
One of the challenges when simulating astrophysical flows with self-gravity is to compute the gravitational forces. In contrast to the hyperbolic compressible Euler equations, the gravity field is described by an elliptic Poisson equation. In [1], we present a purely hyperbolic approach by reformulating the elliptic problem into a hyperbolic diffusion problem, which is solved in pseudotime using the same explicit high-order discontinuous Galerkin method we use for the flow solution.
We start with the Poisson equation for the Newtonian gravitational potential $\phi$,
$$-\vec{\nabla}^2\phi = -4\pi G \rho\label{grav}$$
where G is the universal gravitational constant and $\rho$ is the mass density. Following Nishikawa’s work, we convert the Poisson equation for the gravitational potential into the hyperbolic gravity equations,
where the auxiliary variables $(q_1,q_2)^\intercal\approx\vec{\nabla}\phi$ and $T_r$ is the relaxation time. The steady-state solution of the hyperbolic system is, in fact, the desired solution of the original Poisson problem.
The flow and the gravity solvers operate on a joint hierarchical Cartesian mesh and are two-way coupled in each Runge-Kutta stage via the source terms. A key benefit of our strategy is that it allows the reuse of existing explicit hyperbolic solvers without modifications, while retaining their advanced features such as non-conforming and solution-adaptive grids.
We implemented this purely hyperbolic approach for self-gravitating gas dynamics in Trixi.jl, a tree-based numerical simulation framework for hyperbolic PDEs, which is written in Julia [2]. It was first validated by computing the Jeans gravitational instability, which demonstrates excellent agreement of the numerical results with the analytical solution:
Evolution of the computed kinetic, internal, and potential energies for the Jeans instability using polynomial order N = 3 on a uniform 16 × 16 mesh. The analytical energies are included for reference.
In a second example, we consider a modification of the Sedov blast wave problem that incorporates the effects of gravitational acceleration and which involves strong shocks and complex fluid interactions. We include it to demonstrate the shock capturing and AMR capabilities of Trixi.jl to resolve the cylindrical Sedov blast wave:
Self-gravitating Sedov explosion approximated with polynomial order $N = 3$ on an AMR grid with four levels of possible refinement at $T = 0.5$ (left) and $T = 1.0$ (right). (first and third row) Two dimensional plots of the density and gravitational potential at two times. The white overlay of squares in the density plots shows the AMR grids used by both solvers. (second and last row) One dimensional slices of the solutions along a line from the origin to the edge of the domain in the positive x-direction.
Finally, we show the evolution of the solution together with the adaptive mesh, which is dynamically refined and coarsened to maintain high accuracy while reducing the overall computation time:
Simulation of self-gravitating Sedov blast wave with adaptive mesh refinement, where adaptation performed after every time step of the compressible Euler solver. The total run time compared to a uniformly refined mesh is reduced by a factor of 19, while the solutions are virtually unchanged.
All results in the paper can be easily reproduced by using the resources provided in [3], including step-by-step instructions for how to obtain and install Trixi.jl.
References
[1] Schlottke-Lakemper, Michael and Winters, Andrew R and Ranocha, Hendrik and Gassner, Gregor J, A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics, submitted to J Comp Phys, 2020. arXiv:2008.10593
[2] Schlottke-Lakemper, Michael and Gassner, Gregor J and Ranocha, Hendrik and Winters, Andrew R, Trixi.jl: A tree-based numerical simulation framework for hyperbolic PDEs written in Julia, 2020. doi:10.5281/zenodo.3996439
[3] Schlottke-Lakemper, Michael and Winters, Andrew R and Ranocha, Hendrik and Gassner, Gregor J, Self-gravitating gas dynamics simulations with Trixi.jl, 2020. doi:10.5281/zenodo.3996575
Trixi.jl is a numerical simulation framework for hyperbolic conservation laws written in Julia. A key objective for the framework is to be useful to both scientists and students. Therefore, next to having an extensible design with a fast implementation, Trixi is focused on being easy to use for new or inexperienced users, including the installation and postprocessing procedures. Its features include:
Hierarchical quadtree/octree grid with adaptive mesh refinement
Native support for 2D and 3D simulations
High-order accuracy in space in time
Nodal discontinuous Galerkin spectral element methods
Kinetic energy-preserving and entropy-stable split forms
Entropy-stable shock capturing
Explicit low-storage Runge-Kutta time integration
Square/cubic domains with periodic and Dirichlet boundary conditions
Multiple governing equations:
Compressible Euler equations
Magnetohydrodynamics equations
Hyperbolic diffusion equations for elliptic problems
One of the challenges when simulating astrophysical flows with self-gravity is to compute the gravitational forces. In contrast to the hyperbolic hydrodynamic equations, the gravity field is described by an elliptic Poisson equation. We present a purely hyperbolic approach by reformulating the elliptic problem into a hyperbolic diffusion problem, which is solved in pseudotime using the same explicit high-order discontinuous Galerkin method we use for the flow solution. The flow and the gravity solvers operate on a joint hierarchical Cartesian mesh and are two-way coupled via the source terms. A key benefit of our approach is that it allows the reuse of existing explicit hyperbolic solvers without modifications, while retaining their advanced features such as non-conforming and solution-adaptive grids. By updating the gravitational field in each Runge-Kutta stage of the hydrodynamics solver, high-order convergence is achieved even in coupled multi-physics simulations. After verifying the expected order of convergence for single-physics and multi-physics setups, we validate our approach by a simulation of the Jeans gravitational instability. Furthermore, we demonstrate the full capabilities of our numerical framework by computing a self-gravitating Sedov blast with shock capturing in the flow solver and adaptive mesh refinement for the entire coupled system.
