News

  • Snapshot: Immersed Boundary Discontinuous Galerkin Simulation of Flow Past a Cylinder at Re=100 November 30, 2021

    Simulation of the viscous flow past a cylinder with a discontinuous Galerkin spectral element method on a hierarchical Cartesian mesh with adaptive mesh refinement (AMR). The polynomial degree is N = 3 and the mesh uses five refinement levels. The cylinder boundary is imposed with an immersed boundary method (IBM) based on volume penalization [1], which adds a source term to the compressible Navier-Stokes equations,
    $$
    \partial_t \mathbf{u} + \nabla \cdot \vec{\mathbf{f}} (\mathbf{u},\nabla \mathbf{u}) = \underbrace{\beta (\mathbf{u}_s – \mathbf{u})}_{\text{IBM source}},
    $$
    where $\mathbf{u}$ is the vector of conservative quantities, $\vec{\mathbf{f}}$ is the Navier-Stokes flux, $\beta$ is the penalty parameter, and $\mathbf{u}_s$ is an imposed state. For this simulation, $\mathbf{u}_s$ mimics an isothermal no-slip wall boundary condition.

    The simulation was obtained with the open-source code FLUXO (https://github.com/project-fluxo/fluxo).

    References
    [1] Kou, J., Joshi, S., Hurtado-de-Mendoza, A., Puri, K., Hirsch, C., & Ferrer, E. (2022). Immersed boundary method for high-order flux reconstruction based on volume penalization. Journal of Computational Physics, 448, 110721.

  • New paper published: On the theoretical foundation of overset grid methods for hyperbolic problems: Well-posedness and conservation October 27, 2021

    We use the energy method to study the well-posedness of initial-boundary value problems approximated by overset mesh methods in one and two space dimensions for linear constant-coefficient hyperbolic systems. We show that in one space dimension, for both scalar equations and systems of equations, the problem where one domain partially oversets another is well-posed when characteristic coupling conditions are used. If a system cannot be diagonalized, as is usually the case in multiple space dimensions, then the energy method does not give proper bounds in terms of initial and boundary data. For those problems, we propose a novel penalty approach. We show, by using a global energy that accounts for the energy in the overlap region of the domains, that under well-defined conditions on the coupling matrices the penalized overset domain problems are energy bounded, conservative, well-posed and have solutions equivalent to the original single domain problem.

    Published in Journal of Computational Physics (ScienceDirect)

  • Snapshot: Simulation of a stiff multi-species reacting flow October 19, 2021

    Simulation of an asymmetric two-dimensional stiff multi-species reacting flow for the compressible Euler equations. The computation domain [0,6]x[0,2] of this detonation problem is split into three zones. Zone A [0,0.5]x[0,2] as well as zone B (0.5,6]x[1.2,2] is filled with burnt gas, namely O2,OH and H2O. The last zone C (0.5,6]x[0,1.2) contains unburnt gas, namely H2 and O2.

    The reacting model consists of two reactions H2+O2->2OH, 2OH+H2->2H2O and five species H2,O2,OH,H2O,N2 whereas reaction 1 has a smaller ignition temperature as well as a much faster chemical rate than reaction 2. Species N2 acts as a dilute catalyst here.

    This example has been calculated using Trixi.jl with a P4est mesh using 3072x1024=3.145.728 elements, Dirichlet boundary conditions in x-direction as well as solid wall boundaries in y-direction. The final time is T = 0.1. The convection part is solved with FV whereas the chemical network is solved with KROME. (Note: Not all of these features are currently available in the main code of Trixi.jl.)

    Pressure distribution at final time T=0.1:

    Temperature distribution at final time T=0.1:

    Density distribution at final time T=0.1:

    OH-fraction distribution at final time T=0.1:

  • Snapshot: Numerical simulations of earth’s climate September 30, 2021

    Earth’s climate is changing due to human activity. Therefore, the ability to predict future climate change is of utmost importance if we want to reduce the human interference with the climate system (climate change mitigation) and our vulnerability to the harmful effects of climate change (climate change adaptation).

    It is possible to model the time-dependent local temperatures on the globe using a simple two-dimensional energy balance model [1,2],
    $$ C(\vec{x}) \frac{\partial T}{\partial t} + A + BT = \nabla \cdot \left( D(\vec{x}) \nabla T \right) + Q S(\vec{x},t) (1-a(\vec{x})), $$
    where $T$ is the local temperature; $C(\vec{x})$ is the local heat capacity for a column of air on land, water or ice; the term $A+BT$ models the outgoing long-wave radiation, where $A$ depends on the amount of green-house gases in the atmosphere and $B$ models several feedback effects; $D$ is a latitude-dependent diffusion coefficient; and the last term includes the heat input in the system, where $Q$ is the solar radiation at earth’s position in the solar system, $S$ is the seasonal and latitudinal distribution of insolation, which depends on orbital parameters, such as the eccentricity of earth’s orbit, and the precession and tilt of earth’s rotation axis, and $a(\vec{x})$ is the local albedo, which measures the reflection of solar radiation and depends on the surface cover.

    We developed a julia code that implements the EBM model on the sphere surface and used it to describe the mean temperatures on the globe assuming a constant atmospheric CO$_2$ concentration (of year 1950). The results are presented in the animation.

    [1] North, G.R.; Mengel, J.G.; Short, D.A. (1983). Simple Energy Balance Model Resolving the Seasons and the Continents’ Application to the Astronomical Theory of the Ice Ages. Journal of Geophysical Research.
    [2] Zhuang, K.; North, G.R.; Stevens, Mark J. (2017). A NetCDF version of the two-dimensional energy balance model based on the full multigrid algorithm. SoftwareX.

  • Snapshot: Simulation of a colliding flow setup with DGSEM August 15, 2021
    Simulation of a colliding flow setup with DGSEM (polynomial degree 3,
    compatible sub-cell LGL type FV for shock capturing) of the compressible
    Euler equations, where at x = +/- 64 inflow boundary conditions with
    Mach number Ma = 70 is set (with periodic boundary conditions in y
    direction).

    Density distribution, where the color bar is plotted in logarithmic scale, at
    final simulation time t=25:

    Simulation with AMR (level 3 up to level 9), showing the AMR indicator at
    final time t=25:

    Plot of the DGSEM/FV blending factor for shock capturing at final time
    t=25:

    The combination of the AMR indicator and the shock capturing blending
    factor triggers the final mesh refinement, as can be seen here at the final time t=25:

  • JuliaCon 2021: Adaptive and extendable numerical simulations with Trixi.jl August 9, 2021

    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.

    This talk was given on July 30th, 2021 by Michael Schlottke-Lakemper and Hendrik Ranocha as part of JuliaCon 2021.

    Talk on YouTube: https://www.youtube.com/watch?v=hoViWRAhCBE
    Repository: https://github.com/trixi-framework/talk-2021-juliacon
    Conference agenda entry: https://live.juliacon.org/talk/VAGFD7

  • New paper published: An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. Part II: Subcell finite volume shock capturing August 3, 2021

    The second paper of this series presents two robust entropy stable shock-capturing methods for discontinuous Galerkin spectral element (DGSEM) discretizations of the compressible magneto-hydrodynamics (MHD) equations. Specifically, we use the resistive GLM-MHD equations, which include a divergence cleaning mechanism that is based on a generalized Lagrange multiplier (GLM). For the continuous entropy analysis to hold, and due to the divergence-free constraint on the magnetic field, the GLM-MHD system requires the use of non-conservative terms, which need special treatment.

    Hennemann et al. (2020) [25] recently presented an entropy stable shock-capturing strategy for DGSEM discretizations of the Euler equations that blends the DGSEM scheme with a subcell first-order finite volume (FV) method. Our first contribution is the extension of the method of Hennemann et al. to systems with non-conservative terms, such as the GLM-MHD equations. In our approach, the advective and non-conservative terms of the equations are discretized with a hybrid FV/DGSEM scheme, whereas the visco-resistive terms are discretized only with the high-order DGSEM method. We prove that the extended method is semi-discretely entropy stable on three-dimensional unstructured curvilinear meshes. Our second contribution is the derivation and analysis of a second entropy stable shock-capturing method that provides enhanced resolution by using a subcell reconstruction procedure that is carefully built to ensure entropy stability.

    We provide a numerical verification of the properties of the hybrid FV/DGSEM schemes on curvilinear meshes and show their robustness and accuracy with common benchmark cases, such as the Orszag-Tang vortex and the GEM (Geospace Environmental Modeling) reconnection challenge. Finally, we simulate a space physics application: the interaction of Jupiter’s magnetic field with the plasma torus generated by the moon Io.

    Published in Journal of Computational Physics (ScienceDirect)

  • Snapshot: Interaction of Jupiter’s moon Io with the Jovian magnetosphere July 31, 2021

    Io is the most volcanically active body of the solar system. Furthermore, it is embedded in Jupiter’s magnetic field, the largest and most powerful planetary magnetosphere of the solar system. Due to its strong volcanic activity, Io expels ions and neutrals, which are in turn ionized by ultraviolet and electron impact ionization, forming a plasma torus around Jupiter [1,2]. As Io moves inside the plasma torus, elastic collisions of ions and neutrals inside its atmosphere generate a magnetospheric disturbance that propagates away from Io along the background magnetic field lines at the Alfvén wave speed. This phenomenon creates a pair of Alfvén current tubes that are commonly called Alfvén wings, which have been observed by several flybys [1].

    The figure shows the momentum (non-dimensional) and magnetic field of the plasma that surrounds Io, obtained with the magneto-hydrodynamic (MHD) plasma model of FLUXO*. The Alfvén wings can be observed as disturbances of both the background magnetic and momentum fields. In yellow is the trajectory of the I31 flyby of the Galileo spacecraft, which visited Io in 2001 [1]. The Galileo spacecraft is also depicted (not to scale). A 3D model of the moon’s surface developed by NASA [3] was superposed on the simulation results.

    *FLUXO (www.github.com/project-fluxo/fluxo) is an MPI parallel high-order Discontinuous Galerkin code, which supports unstructured curvilinear hexahedral grids, and is able to perform Adaptive Mesh Refinement (AMR).

    [1] M. Kivelson, K. Khurana, C. Russell, R. Walker, S. Joy, J. Mafi, GALILEO ORBITER AT JUPITER CALIBRATED MAG HIGH RES V1.0, GO-J-MAG-3-RDR-HIGHRES-V1.0, Technical Report, NASA Planetary Data System, 1997.
    [2] J. Saur, F. M. Neubauer, J. E. P. Connerney, Plasma interaction of Io with its plasma torus, 2004.
    [3] https://solarsystem.nasa.gov/resources/2379/io-3d-model/

  • Talk on 2021-07-02: On Wave Propagation Characteristics, Upwind SBP Properties and Energy Stability of DG Viscous Flux Discretizations June 28, 2021

    Speaker: Dr. Sigrun Ortleb, University of Kassel, Germany
    Date & time: Friday, 2nd July 2021, 10 am (CEST)
    Venue: Please request the Zoom meeting link from Michael Schlottke-Lakemper

    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.

  • Snapshot: Evolution of a magnetized torus differentially rotating around a static point mass June 24, 2021

    Evolution of a magnetized torus differentially rotating around a static point mass. [1] The torus gradually forms into an accretion disk with expanding filaments into the surrounding space due to magneto-rotational instabilities similar to the rich structure of the corona of our sun. The simulation was done with a new Discontinuous-Galerkin based MHD solver [2] in FLASH [3]. Fig. 1 and Fig. 2 show the log. scale density cross section in the x-y and x-z plane. Fig 3 highlights the rich structure of the magnetic field (white stream lines) overlayed on top of the log. scale magnetic pressure in the x-z plane.

     

    [1] Machida, Mami, Mitsuru R. Hayashi, and R. Matsumoto.
    “Global simulations of differentially rotating magnetized disks:
    Formation of low-beta filaments and structured coronae.”
    The Astrophysical Journal Letters 532.1 (2000): L67.
    [2] Markert et al. “A Discontinuous Galerkin Method with Sub‑Cell Adaptive Shock Capturing for
    FLASH” (in preparation)
    [3] http://flash.uchicago.edu/

     

  • Snapshot: Gingerbread man simulation with unstructured quadrilateral elements May 25, 2021

    Here we are simulating a gingerbread man, whereby the gingerbread man itself has 861 unstructured quadrilateral elements. For the run parameters, we used polydeg = 10 with the compressible Euler equations and the manufactured (convergence test) solution. The visualization is done on nvisnodes=12 uniform points in each direction, on each element.

    We thank Prof. Dr. Andrew Winters (https://liu.se/en/employee/andwi94) for providing this video/simulation.

  • Talk on 2021-05-21: On a linear stability issue of split form schemes for compressible flows May 19, 2021

    Speaker: Dr. Vikram Singh, Max-Planck-Institute for Meteorology
    Date & time: Friday, 21st May 2021, 10 am (CEST)
    Meeting link: Please request via email from Michael Schlottke-Lakemper

    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.

    Note: This meeting will be recorded.

  • New paper submitted: On the Theoretical Foundation of Overset Grid Methods for Hyperbolic Problems: Well-Posedness and Conservation May 12, 2021

    We use the energy method to study the well-posedness of initial-boundary value problems approximated by overset mesh methods in one and two space dimensions for linear constant-coefficient hyperbolic systems. We show that in one space dimension, for both scalar equations and systems of equations, the problem where one domain partially oversets another is well-posed when characteristic coupling conditions are used. If a system cannot be diagonalized, as is ususally the case in multiple space dimensions, then the energy method does not give proper bounds in terms of initial and boundary data. For those problems, we propose a novel penalty approach. We show, by using a global energy that accounts for the energy in the overlap region of the domains, that under well-defined conditions on the coupling matrices the penalized overset domain problems are energy bounded, conservative, well-posed and have solutions equivalent to the original single domain problem.

    Preprint available at: https://arxiv.org/abs/2105.04664

  • Talk on 2021-05-10: A two-dimensional stabilized discontinuous Galerkin method on curvilinear embedded boundary grids May 8, 2021

    Speaker: Dr. Andrew Giuliani, New York University
    Date & time: Monday, 10th May 2021, 4pm (CEST)/10am (EDT)
    Meeting link: Please request via email from Michael Schlottke-Lakemper

    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$.

  • New paper published: Stability of Discontinuous Galerkin Spectral Element Schemes for Wave Propagation when the Coefficient Matrices have Jumps May 4, 2021

    Abstract:

    We use the behavior of the L2 norm of the solutions of linear hyperbolic equations with discontinuous coefficient matrices as a surrogate to infer stability of discontinuous Galerkin spectral element methods (DGSEM). Although the L2 norm is not bounded by the initial data for homogeneous and dissipative boundary conditions for such systems, the L2 norm is easier to work with than a norm that discounts growth due to the discontinuities. We show that the DGSEM with an upwind numerical flux that satisfies the Rankine-Hugoniot (or conservation) condition has the same energy bound as the partial differential equation does in the L2 norm, plus an added dissipation that depends on how much the approximate solution fails to satisfy the Rankine-Hugoniot jump.

    Accepted in Journal of Scientific Computing.

    arXiv link: https://arxiv.org/abs/2011.11746

  • Snapshot: Astrophysical colliding flow simulation run by a novel Discontinuous Galerkin/Finite Volume (DGFV) blending scheme in FLASH April 29, 2021
    Simulation of an astrophysical colliding flow [1] run by a novel Discontinuous Galerkin/Finite Volume (DGFV) blending scheme [2] which has been implemented in the FLASH code [3].

    The code has following capabilities:

    • fourth-order accurate ideal magneto-hydrodynamics with hyperbolic divergence cleaning [4]
    • octree-based adaptive mesh refinement
    • distributive computing and load balancing
    • multi-species fluid dynamics (N_species > 10)
    • turbulent driving
    • octree-based Poisson solver for self-gravity [5]
    • octree-based radiation physics [6]
    • external gravitional fields
    • sink particles [7]
    • chemical reaction networks [8]

    [1] Weis, Micheal et al. “The Virial Balance of CO-Substructures in Colliding Magnetised Flows” (in preparation)
    [2] Markert, Johannes et al. “A Sub-Element Adaptive Shock Capturing Approach for Discontinuous Galerkin Methods” (submitted)
    [3] Fryxell, Bruce, et al. “FLASH: An adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes.” The Astrophysical Journal Supplement Series 131.1 (2000): 273.
    [4] Markert, Johannes et al. “Flash goes DG” (working title, in preparation)
    [5] Wünsch, Richard, et al. “Tree-based solvers for adaptive mesh refinement code FLASH–I: gravity and optical depths.” Monthly Notices of the Royal Astronomical Society 475.3 (2018): 3393-3418.
    [6] Wünsch, Richard et al. “Tree-based solvers for adaptive mesh refinement code FLASH – II: radiation transport module TreeRay” (submitted)
    [7] Federrath, Christoph, et al. “Modeling collapse and accretion in turbulent gas clouds: implementation and comparison of sink particles in AMR and SPH.” The Astrophysical Journal 713.1 (2010): 269.
    [8] Seifried, D., and S. Walch. “Modelling the chemistry of star-forming filaments–I. H2 and CO chemistry.” Monthly Notices of the Royal Astronomical Society: Letters 459.1 (2016): L11-L15.

     

  • Snapshot: Simulation of a Kelvin-Helmholtz instability using second order Finite Volume schemes and fourth order Discontinuous Galerkin methods April 1, 2021

    We present in-viscid and viscous simulations of a Kelvin-Helmholtz instability using second a order accurate monotoniced-central finite volume (FV) method and a fourth order accurate discontinuous Galerkin (DG) method. The initial condition is given by [1]:

    $$\rho (t=0) = \frac{1}{2}
    + \frac{3}{4} B,
    ~~~~~~~~~
    p (t=0) = 1,~~~~~~~~~~
    $$

    $$
    v_1 (t=0) = \frac{1}{2} \left( B-1 \right),
    ~~~~~~~
    v_2 (t=0) = \frac{1}{10} \sin(2 \pi x),
    $$

    with $$B=\tanh \left( 15 y + 7.5 \right) – \tanh(15y-7.5).$$

    We first present the FV results at end time $t=3.7$, which are computed using a monotoniced-central second order discretization of the Euler equations of gas dynamics on uniform grids.

    The next results use a fourth order DG discretization of the Navier-Stokes equations on uniform grids using $Re=320.000$ at end time $t=3.7$. The highest resolution (4096² DOFs) is a direct numerical simulation (DNS) of the problem, where all scales are resolved.

    It is remarkable that the numerical dissipation of the second order FV scheme causes the in-viscid simulation with 2048² DOFs to look very similar to the viscous DNS solution at $Re=320.000$.

    [1] A.M. Rueda-Ramírez, G.J Gassner (2021). A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations of the Euler Equations. https://arxiv.org/pdf/2102.06017.pdf

  • Talk on 2021-03-18: A tour of BifurcationKit and some results on mean fields of spiking neurons March 12, 2021

    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.

  • New proceedings paper submitted: A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations of the Euler Equations February 25, 2021

    In this paper, we present a positivity-preserving limiter for nodal Discontinuous Galerkin disctretizations of the compressible Euler equations. We use a Legendre-Gauss-Lobatto (LGL) Discontinuous Galerkin Spectral Element Method (DGSEM) and blend it locally with a consistent LGL-subcell Finite Volume (FV) discretization using a hybrid FV/DGSEM scheme that was recently proposed for entropy stable shock capturing. We show that our strategy is able to ensure robust simulations with positive density and pressure when using the standard and the split-form DGSEM. Furthermore, we show the applicability of our FV positivity limiter in extremely under-resolved vortex dominated simulations and in problems with shocks.

    Preprint available at: https://arxiv.org/pdf/2102.06017.pdf

    Evolution of the density for a Sedov blast simulation with periodic boundaries

    Evolution of the density for a Sedov blast simulation with periodic boundaries

  • Snapshot: Accuracy of the LGL-subcell FV scheme with and without inner-element reconstruction February 24, 2021

    We are interested in the accuracy of the finite volume scheme, on a LGL subcell grid induced by a LGL-DGSEM discretization. For comparison, we look at a Kelvin-Helmholtz-Instability test problem, simulated with 4² LGL nodes per element and 32² elements, i.e. 128² spatial degrees of freedom. The high-order DGSEM uses the entropy-conserving split-form powered by the flux of Chandrashekar in the volume, and the HLLC flux at the element interfaces. Positivity is controlled by adding subcell FV where the positivity is not fulfilled with the amount needed. All FV discretizations use the HLLC flux at the element interfaces and at the subcell interfaces as well.  

    The first results show the high-order DGSEM result, which is formally 4th order accurate for smooth problems:

    The next result uses a subcell finite volume approximation on the LGL-subcell grid, directly, i.e. without spatial reconstruction (piecewise constant approximation). The result is very dissipative:

    The last results show the improvement when using a piecewise linear reconstruction with monotonized-central slope limiter on the LGL-subcell grid. The reconstruction is inner-element only, as it does not use element neighbor information. Thus, at the very first and very last LGL subcell, no reconstruction is used. Still, the accuracy is recovered nicely by this simple modification of the subcell FV scheme: 

  • Snapshot: Single-Node Performance Comparison of Four Different Magnetohydrodynamics Codes January 30, 2021

    We compare the runtime performance of four different magnetohydrodynamics codes on a single compute node on the in-house high performance cluster ODIN. A compute node on ODIN consists of 16 cores. We run the ‘3D Alfven wave’ test case up to a fixed simulation time and measure the elapsed wall clock time of each code minus initialization time and input/output operations. For each run the number of cores is successively increased. This allows us to get insights into the scaling behavior (speedup gain wih increasing number of cores) on a single compute node. Furthermore we plot the performance index (PID) over the number of nodes which is a measure of throughput, i.e. how many millions of data points (degrees-of-freedom/DOF) per second are processed by the each code.

    The four contestants are:

    • Flash [1] with Paramesh 4.0 and the Five-wave Bouchut Finite-Volume solver written in Fortran
    • Fluxo [2] an MHD Discontinuous Galerkin Spectral Element code written in Fortran
    • Trixi [3] an MHD Discontinuous Galerkin Spectral Element code written in Julia
    • Nemo an in-house prototype code providing a 2nd order monotonized-central MHD-FV scheme (MCFV) and a hybrid MHD Discontinuous Galerkin Spectral Element / Finite Volume scheme (DGFV) written in Fortran.

    Trixi uses a thread-based parallelization approach while Flash, Fluxo and Nemo
    use the standard MPI method. Furthermore, Nemo also provides OpenMP-based
    parallelization.

    [1] https://flash.uchicago.edu/site/
    [2] https://github.com/project-fluxo/fluxo
    [3] https://github.com/trixi-framework/Trixi.jl

  • Talk on 2021-01-21: Julia for adaptive high-order multi-physics simulations January 25, 2021

    On Wednesday, 27th January 2021, 3:15pm CET, Michael Schlottke-Lakemper will give an online talk on

    Julia for adaptive high-order simulations

    To obtain the Zoom link, please contact the organizers via the official meeting announcement.

    Abstract

    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).

  • New paper submitted: An Entropy Stable Nodal Discontinuous Galerkin Method for the resistive MHD Equations. Part II: Subcell Finite Volume Shock Capturing January 4, 2021

    The second paper of this series presents two robust entropy stable shock-capturing methods for discontinuous Galerkin spectral element (DGSEM) discretizations of the compressible magneto-hydrodynamics (MHD) equations. Specifically, we use the resistive GLM-MHD equations, which include a divergence cleaning mechanism that is based on a generalized Lagrange multiplier (GLM). For the continuous entropy analysis to hold, and due to the divergence-free constraint on the magnetic field, the GLM-MHD system requires the use of non-conservative terms, which need special treatment.

    Hennemann et al. [DOI:10.1016/j.jcp.2020.109935] recently presented an entropy stable shock-capturing strategy for DGSEM discretizations of the Euler equations that blends the DGSEM scheme with a subcell first-order finite volume (FV) method. Our first contribution is the extension of the method of Hennemann et al. to systems with non-conservative terms, such as the GLM-MHD equations. In our approach, the advective and non-conservative terms of the equations are discretized with a hybrid FV/DGSEM scheme, whereas the visco-resistive terms are discretized only with the high-order DGSEM method. We prove that the extended method is entropy stable on three-dimensional unstructured curvilinear meshes. Our second contribution is the derivation and analysis of a second entropy stable shock-capturing method that provides enhanced resolution by using a subcell reconstruction procedure that is carefully built to ensure entropy stability.

    We provide a numerical verification of the properties of the hybrid FV/DGSEM schemes on curvilinear meshes and show their robustness and accuracy with common benchmark cases, such as the Orszag-Tang vortex and the GEM reconnection challenge. Finally, we simulate a space physics application: the interaction of Jupiter’s magnetic field with the plasma torus generated by the moon Io.

    Preprint available at: arXiv:2012.12040

  • Snapshot: Hybrid Discontinuous Galerkin / Finite Volume (DG/FV) simulation of the multicomponent Shock-Bubble Interaction Problem December 24, 2020

    Simulation of the Shock-Bubble Interaction Problem with a hybrid entropy-stable DG/FV method and Adaptive Mesh Refinement (AMR) using Trixi (https://github.com/trixi-framework/Tr…).

    The simulation uses a multicomponent model and computes the spatial operator as a blend of a first-order subcell FV method and a fourth-order DG scheme on a squared domain.

  • New paper published: A Novel Robust Strategy for Discontinuous Galerkin Methods in Computational Fluid Mechanics: Why? When? What? Where? December 15, 2020

    In this paper we will review a recent emerging paradigm shift in the construction and analysis of high order Discontinuous Galerkin (DG) methods applied to approximate solutions of hyperbolic or mixed hyperbolic-parabolic partial differential equations (PDEs) in computational physics. There is a long history using DG methods to approximate the solution of PDEs in computational physics with successful applications in linear wave propagation, like those governed by Maxwell’s equations, incompressible and compressible fluid and plasma dynamics governed by the Navier-Stokes and the Magnetohydrodynamics equations, or as a solver for ordinary differential equations (ODEs), e.g., in structural mechanics. The DG method amalgamates ideas from several existing methods such as the Finite Element Galerkin method (FEM) and the Finite Volume method (FVM) and is specifically applied to problems with advection dominated properties, such as fast moving fluids or wave propagation. In the numerics community, DG methods are infamous for being computationally complex and, due to their high order nature, as having issues with robustness, i.e., these methods are sometimes prone to crashing easily. In this article we will focus on efficient nodal versions of the DG scheme and present recent ideas to restore its robustness, its connections to and influence by other sectors of the numerical community, such as the finite difference community, and further discuss this young, but rapidly developing research topic by highlighting the main contributions and a closing discussion about possible next lines of research.

    Published in: https://www.frontiersin.org/articles/10.3389/fphy.2020.500690/abstract

    Front. Phys. | doi: 10.3389/fphy.2020.500690

  • New paper submitted: A Split-Form, Stable CG/DG-SEM for Wave Propagation Modeled by Linear Hyperbolic Systems December 15, 2020

    We present a hybrid continuous and discontinuous Galerkin spectral element approximation that leverages the advantages of each approach. The continuous Galerkin approximation is used on interior element faces where the equation properties are continuous. A discontinuous Galerkin approximation is used at physical boundaries and if there is a jump in properties at a face. The approximation uses a split form of the equations and two-point fluxes to ensure stability for unstructured quadrilateral/hexahedral meshes with curved elements. The approximation is also conservative and constant state preserving on such meshes. Spectral accuracy is obtained for all examples, which include wave scattering at a discontinuous medium boundary.

    Preprint available at: arXiv:2012.06510

  • Snapshot: Hybrid Discontinuous Galerkin / Finite Volume (DG/FV) simulation of the Orszag-Tang vortex November 27, 2020

    Simulation of the Orszag-Tang vortex test with a hybrid entropy-stable DG/FV method and Adaptive Mesh Refinement (AMR) using FLUXO (https://github.com/project-fluxo/fluxo).

    The simulation uses the GLM-MHD model to control div(B) errors and computes the spatial operator as a blend of a first-order subcell FV method and a fourth-order DG scheme.

    The initial grid has 16×16 elements (64² DOFs), and the maximum resolution is achieved with three refinement levels (equivalent to 512² DOFs).

  • New paper submitted: Stability of Discontinuous Galerkin Spectral Element Schemes for Wave Propagation when the Coefficient Matrices have Jumps November 25, 2020

    We use the behavior of the L2 norm of the solutions of linear hyperbolic equations with discontinuous coefficient matrices as a surrogate to infer stability of discontinuous Galerkin spectral element methods (DGSEM). Although the L2 norm is not bounded by the initial data for homogeneous and dissipative boundary conditions for such systems, the L2 norm is easier to work with than a norm that discounts growth due to the discontinuities. We show that the DGSEM with an upwind numerical flux that satisfies the Rankine-Hugoniot (or conservation) condition has the same energy bound as the partial differential equation does in the L2 norm, plus an added dissipation that depends on how much the approximate solution fails to satisfy the Rankine-Hugoniot jump.

    Preprint available at arXiv:2011.11746

  • New paper submitted: A Sub-Element Adaptive Shock Capturing Approach for Discontinuous Galerkin Methods November 23, 2020

    In this paper, a new strategy for a sub-element based shock capturing for discontinuous Galerkin (DG) approximations is presented. The idea is to interpret a DG element as a collection of data and construct a hierarchy of low to high order discretizations on this set of data, including a first order finite volume scheme up to the full order DG scheme. The different DG discretizations are then blended according to sub-element troubled cell indicators, resulting in a final discretization that adaptively blends from low to high order within a single DG element. The goal is to retain as much high order accuracy as possible, even in simulations with very strong shocks, as e.g. presented in the Sedov test. The framework retains the locality of the standard DG scheme and is hence well suited for a combination with adaptive mesh refinement and parallel computing. The numerical tests demonstrate the sub-element adaptive behavior of the new shock capturing approach and its high accuracy.

    Preprint available at arXiv:2011.03338.

  • Snapshot: 2D Riemann problem computed with a hybrid Discontinuous Galerkin / Finite Volume scheme October 23, 2020

    We present a 2D Riemann problem (configuration 3) computed with a hybrid Discontinuous Galerkin / Finite Volume scheme on a uniform grid with 1024^2 degrees of freedom.

    In this hybrid approach we are blending an 8th order, 4th order and 2nd order DG as well as a 1st order FV scheme.

  • New paper submitted: Preventing pressure oscillations does not fix local linear stability issues of entropy-based split-form high-order schemes September 29, 2020

    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).

    Preprint available at arXiv:2009.13139. Numerical results were obtained with Trixi.jl.

  • Snapshot: Purely hyperbolic gravity simulations with Trixi.jl September 28, 2020

    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,

    $$\frac{\partial}{\partial t}\begin{bmatrix}
    \phi\\[0.1cm]
    q_1\\[0.1cm]
    q_2\\[0.1cm]
    \end{bmatrix}
    +
    \frac{\partial}{\partial x}\begin{bmatrix}
    -q_1\\[0.1cm]
    -\phi/T_r\\[0.1cm]
    0\\[0.1cm]
    \end{bmatrix}
    +
    \frac{\partial}{\partial y}\begin{bmatrix}
    -q_2\\[0.1cm]
    0\\[0.1cm]
    -\phi/T_r\\[0.1cm]
    \end{bmatrix}
    =
    \begin{bmatrix}
    -4\pi G \rho\\[0.1cm]
    -q_1/T_r\\[0.1cm]
    -q_2/T_r\\[0.1cm]
    \end{bmatrix},\label{hyp}$$

    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:

    https://raw.githubusercontent.com/trixi-framework/paper-self-gravitating-gas-dynamics/master/assets/sedov-rho-phi-mesh.gif
    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: A tree-based numerical simulation framework for hyperbolic PDEs written in Julia September 1, 2020

    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
      • Scalar advection
    • Multi-physics simulations
    • Shared-memory parallelization via multithreading
    • Visualization of results with Julia-only tools (2D) or ParaView/VisIt (2D/3D)

    Trixi.jl was initiated by Michael Schlottke-Lakemper and Gregor Gassner (both University of Cologne, Germany). Together with Hendrik Ranocha (KAUST, Saudi Arabia) and Andrew Winters (Linköping University, Sweden), they are the principal developers of Trixi.

    In case of questions, please feel free to create an issue. We are looking forward to feedback and/or potential scientific collaboration.

  • New paper submitted: A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics August 26, 2020

    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.

    Available at arXiv:2008.10593

  • Snapshot: A new approach for approximating conservation laws August 24, 2020

    A new approach for approximating conservation laws has been tested: instead of monitoring the changes of the means of a quantity within a certain volume over time like a finite volume procedure, this method mimics the behavior of simple solutions.

    First the current state is split into waves w, each aligned along one of the eigenvectors of the Fluxjacobian and living on its own grid. Then the corresponding grids are being moved with their eigenvelocity and finally all waves w are being overlayed and normalized yielding the solution after one timestep.

    The plots show the solution of isothermal Euler equations at Time T=0.0000 (initial condition) and T=0.0025 with N=40 gridpoints for each mesh.

  • Snapshot: Higher-order schemes for the MHD equations July 31, 2020

    A robust and easy way of simulating a hyperbolic test case with discontinuities is to use a 1st order finite volume scheme. Using such a method for magnetohydrodynamics (MHD) problems like the Orszag-Tang vortex leads to following results:

    For this and the following examples we used a 4th order time integration scheme and 256 degrees of freedom in each spatial direction.

    A way to generate more accurate results is to increase the order of the scheme, which has to be treated with caution near discontinuities because oscillations may occur. To overcome this issue one could use higher order schemes in smooth regions and lower order schemes in regions with discontinuities. An example for such an approach is our so called DGFV scheme, which blends e.g. a 4th order Discontinuous Galerkin scheme with a 1st order Finite Volume scheme.

    Another way of doing it, is to use a suitable 4th order finite-volume scheme for MHD with a fitting limiter.

  • New paper published: Entropy-Stable p-Nonconforming Discretizations with the Summation-by-Parts Property for the Compressible Navier–Stokes Equations July 16, 2020

    The entropy-conservative/stable, curvilinear, nonconforming, p-refinement algorithm for hyperbolic conservation laws of Del Rey Fernández et al. (2019) is extended from the compressible Euler equations to the compressible Navier–Stokes equations. A simple and flexible coupling procedure with planar interpolation operators between adjoining nonconforming elements is used. Curvilinear volume metric terms are numerically approximated via a minimization procedure and satisfy the discrete geometric conservation law conditions. Distinct curvilinear surface metrics are used on the adjoining interfaces to construct the interface coupling terms, thereby localizing the discrete geometric conservation law constraints to each individual element. The resulting scheme is entropy conservative/stable, element-wise conservative, and freestream preserving. Viscous interface dissipation operators that retain the entropy stability of the base scheme are developed. The accuracy and stability of the resulting numerical scheme are shown to be comparable to those of the original conforming scheme in Carpenter et al. (2014) and Parsani et al. (2016), i.e., this scheme achieves ~p+1/2 convergence on geometrically high-order distorted element grids; this is demonstrated in the context of the viscous shock problem, the Taylor–Green vortex problem at a Reynolds number of Re = 1,600 and a subsonic turbulent flow past a sphere at Re = 2,000.

    sciencedirect
    arxiv

  • New article published in ECCOMAS Newsletter (p.16-20): Split Form Discontinuous Galerkin Methods For Implicit Large Eddy Simulation Of Compressible Turbulence July 16, 2020

    As a teaser we show a numerical demonstration of the capabilities of the split form DG approach. We consider the flow past a plunging SD7003 airfoil at Mach number and cord length based Reynolds number Rec = 40,000. We subdivide the computational domain into 58,490 unstructured curvilinear hexahedral elements and use a polynomial degree N = 7 resulting in a total of about 150 million degrees of freedom.

    The Newsletter is available via: eccomas Newsletter.

  • Snapshot: Variant of the Boss-Bodenheimer test with a self-gravitating high order hydro code June 26, 2020

    Variant of the classical Boss-Bodenheimer Test: We simulate the collapse and fragmentation of a rotating cloud core using a Discontinous Galerkin / Finite Volume hybrid high-order code. The self-gravity of the gas is calculated via the Fast-Multipole-Method.

  • Snapshot: 2D PyOpenGL real-time simulation of the forward-facing step test case May 27, 2020

    WENO methods refers to a class of nonlinear finite volume or finite difference methods which are well suited to approximate solutions of hyperbolic conservation laws with high order accuracy in smooth regions and essentially non-oscillatory transition for solution discontinuities. In the following video we used a 5th order accurate WENO-FD scheme with immersed boundary conditions.

    Here we solve the forward-facing step test case in real-time using Python3/PyOpenGL/GLFW with an Intel Corporation HD Graphics 620 (rev 02). As one can see, it is also possible to add pressure peaks to the simulation by simply clicking with the mouse on the desired field.

  • Snapshot: Smoothed-Particle Hydrodynamics (SPH) simulation of the Numerical Simulation Research Group April 28, 2020

    Smoothed-particle hydrodynamics (SPH) is a computational method used for simulating the mechanics of continuum media, such as solid mechanics and is being increasingly used to model fluid motion as well.
    This method is well-suited for problems dominated by complex boundary dynamics, since SPH is a mesh-free method, as well as for mass conservation problems since the particles themselves represent mass.

    We used a 2D Python/OpenCL SPH code that solves the incompressible Navier-Stokes equations in real time. In this simulation, 52,700 particles were used to smash the people of the Numerical Simulation Research Group against each other.

  • Snapshot: Blast wave simulation computed in FLUXO with an entropy-stable high-order Discontinuous Galerkin/Finite Volume hybrid scheme March 9, 2020

    Blast wave simulation with periodic boundaries computed in FLUXO with an
    entropy-stable high-order Discontinuous Galerkin/Finite Volume scheme.

    The scheme is formulated in Hennemann and Gassner. “A provably entropy stable subcell shock capturing approach for high order split form DG.” Journal of Computational Physics (submitted). The scheme computes the spatial operator as F = (1-a) F_{DG} + a F_{FV}, where a is the blending function

  • Talk: Lea Miko Versbach, Lund University, Sweden, “Implicit DG solvers and multigrid preconditioners” [05.02.2020, 11am] February 4, 2020

    On Wednesday, 5 February 2020 at 11:00 Lea Miko Versbach will talk on the topic “Implicit DG solvers and multigrid preconditioners” 

    Location: Gyrhofstraße 8a (Gebäude 158a), Room 1.105 (1st floor), 50931

    Abstract: In this talk I will give an introduction to implicit DG solvers and the challenges which arise when solving the large nonlinear systems coming from the implicit temporal discretization. We will solve these system using a Jacobian-free Newton-Krylov method. The problem with this method is that the Krylov subspace method converges slowly for discretized PDEs. In order to speed up computations, we need to construct preconditioners. This can be done using iterative methods. Carefully constructed multigrid methods are cheap and effective preconditioners. I will explain more details about the construction of these preconditioners and show some numerical results. 

  • New paper published: Dynamic load balancing for direct-coupled multiphysics simulations January 29, 2020

    High parallel efficiency for large-scale coupled multiphysics simulations requires the computational load to be evenly distributed among all compute cores. For complex applications and massively parallel computations, even minor load imbalances can have a severe impact on the overall performance and resource usage. Exemplarily for a volume-coupled multiphysics simulation, a direct-hybrid method is considered, in which a CFD and a CAA simulation are performed concurrently on the same parallel subdomains. For differing load compositions on each subdomain, accurate computational weights for CFD and CAA cells must be known to determine an efficient domain decomposition. Therefore, a dynamic load balancing scheme is presented, which allows to increase the efficiency of complex coupled simulations with non-trivial domain decompositions. A fully-coupled three-dimensional jet simulation with approximately 300 million degrees of freedom demonstrates the effectiveness of the approach to reduce load imbalances. A detailed performance analysis substantiates the necessity of dynamic load balancing. Furthermore, the results of a strong scaling experiment show the benefit of load balancing to be proportional to the degree of parallelism. In addition, it is shown that the approach allows to attenuate imbalances also for parallel computations on heterogeneous computing hardware. The acoustic field of a chevron nozzle will also be discussed.

    doi:10.1016/j.compfluid.2020.104437

  • Talk: Bernhard Müller, Department of Energy and Process Engineering, Norwegian University of Science and Technology (NTNU), Trondheim, Norway [04.02.2020, 11am] January 28, 2020

    On Tuesday, 4 February 2020 at 11:00 Bernhard Müller will talk on the topic “Immersed Boundary Method for the Compressible Navier-Stokes Equations Using High Order Summation-By-Parts Operators” based on a joint work with M. Ehsan Khalili and Knut Emil Ringstad, NTNU, and Martin Larsson, Sportradar AS, Trondheim.

    Location: Gyrhofstraße 8a (building 158a), Room 1.105 (1st floor), 50931 Cologne

    Abstract:
    An efficient and versatile immersed boundary method (IBM) for simulating compressible viscous flows with complex and moving convex boundaries has been developed. The compressible Navier-Stokes equations are discretized by globally fourth order summation-by-parts (SBP) difference operators with in-built stability properties and the classical fourth order explicit Runge-Kutta method. The proposed Cartesian grid based IBM builds on the ghost point approach in which the solid wall boundary conditions are applied as sharp interface conditions. The interpolation of the flow variables at image points and the solid wall boundary conditions are used to determine the flow variables at three layers of ghost points within the solid body in order to introduce the presence of the body interface in the flow computation and to maintain the overall high order of accuracy of the flow solver. Two different reconstruction procedures, bilinear interpolation and weighted least squares method, are implemented to obtain the values at the ghost points. A robust high order immersed boundary method is achieved by using a hybrid approach in which the first layer of ghost points is treated by using a third order polynomial combined with the weighted least squares method and the second and third layers of ghost points are treated by using bilinear interpolation to find the values at the image points of the corresponding ghost points. After demonstrating the accuracy of the present IBM for low Mach number flow around a circular cylinder, it is applied to simulate flow in the cross-section of the upper airways of a specific obstructive sleep apnea patient. The IBM solver has been further verified and validated for moving boundaries by applying it to a transversely oscillating cylinder in freestream flow and an in-line oscillating cylinder in an initially quiescent fluid. Sound waves generated by the in-line oscillation of the cylinder exhibit both quadrupole and monopole types. The present IBM is also verified for fluid-structure interaction of an elastically mounted circular cylinder in freestream flow at Reynolds number 200, and the rate of energy transferred between the fluid and the structure is investigated.

  • Snapshot: Postdoctoral researcher Dr. Andrés Rueda-Ramírez joins our research group January 24, 2020

    Andrés M. Rueda-Ramírez completed his undergraduate studies in Mechanical Engineering at the National University of Colombia (Universidad Nacional de Colombia, UNAL). After graduating, he worked as a research and teaching assistant at the Research Group on biomechanics at UNAL, where he collaborated with a highly interdisciplinary team in the development of a Finite Element software to simulate bone growing processes.

    AMRR completed his Ph.D. studies in Aerospace Engineering at the Polytechnic University of Madrid (Universidad Politécnica de Madrid). During his Ph.D., AMRR studied and developed p-adaptation algorithms, implicit time-integration schemes, and multigrid solvers for high-order Discontinuous Galerkin Spectral Element Methods (DGSEM).

    AMRR is now a member of the Numerical Simulation Research Group at the University of Cologne, where he joined the development team of the DGSEM code FLUXO. AMRR is currently working on sub-cell shock-capturing schemes and time integrators for the Navier-Stokes and the MHD equations.

  • Talk: Ruediger Pakmor, Max-Planck-Institute for Astrophysics, Garching January 24, 2020

    On Friday, 31 January 2020 at 14:00 Ruediger Pakmor will talk on the topic “The numerical schemes behind the moving mesh code Arepo”

    Location: Weyertal 86-90, Seminar room 1 (ground floor), 50931 Cologne

    Abstract: I will present the schemes behind the moving mesh magnetohydrodynamics code Arepo and discuss in particular the finite volume scheme on a moving mesh and its time integration. I will then show an example of an anisotropic diffusion solver on the unstructured Voronoi mesh used in Arepo and discuss general properties of the Arepo code.

  • Talk: Cedrick Ansorge, Institut für Geophysik und Meteorologie, Universität zu Köln January 21, 2020

    On Tuesday, 28 January 2020 at 11.00 Cedrick Ansorge will talk on the topic “Turbulent Ekman flow as virtual lab in geophysical fluid dynamics”

    Location:  Gyrhofstraße 8a (Gebäude 158a), Room 1.105 (1st floor), 50931 Cologne

    Abstract: The atmospheric boundary layer (ABL) is the lowest part of the atmosphere that is directly linked to the surface through vertical turbulent exchange, typically the lowest 100 to 1000m. There, turbulent mixing is the main vertical transport mechanism for heat, water, momentum and any kind of air constituent. Besides turbulence, the ABL is a multi-physical system comprising also radiative, miro-physical, chemical and other processes on scales from the multi-kilometre range down to the Kolmogorov scale of turbulent motion at the sub-millimetre range in three spatial dimensions. Both the multiphysical complexity and the broad-scale nature are prohibitive for a brute-force approach to numerical modelling of the system. Truncated representations of the ABL are thus inevitable when numerically modeling the ABL.
    I will introduce turbulent Ekman flow–the doubly periodic flow over a flat rotating plate—that physically truncates the ABL to its fluid-mechanical core, the Navier–Stokes equations with appropriate boundary conditions. The governing equations are solved by a highly scalable numerical algorithm that is being used on up to 250,000 compute cores to represent the turbulent flow on grids that routinely use 3 × 230 collocation points in space. The reduced physical complexity allows for high-fidelity turbulence simulations where the entire range of turbulent motion is represented directly–without need for turbulence closure. We can thus study the ABL under conditions for which classical approaches to represent turbulence fail–such as partial or complete laminarization and transitional flows.


     

  • New paper submitted: A Sub-Element Adaptive Shock Capturing Approach for Discontinuous Galerkin Methods December 18, 2019

    In this paper,a new strategy for a sub-element based shock capturing for discontinuous Galerkin (DG) approximations is presented. The idea is to interpret a DG element as a collection of data and construct a hierarchy of low to high order discretisations on this set of data, including a first order finite volume (FV) scheme up to the full order DG scheme. The different DG discretisations are then blended according to sub-element troubled cell indicators, resulting in a final discretisation that adaptively blends from low to high order within a single DG element. The goal is to retain as much high order accuracy as possible, even in simulations with very strong shocks, as e.g. presented in the Sedov test. The framework retains the locality of the standard DG scheme and is hence well suited for a combination with adaptive mesh refinement (AMR) and parallel computing. The numerical tests demonstrate the sub-element adaptive behaviour of the new shock capturing approach and its high accuracy.


     

  • Talk: Almut Gassmann, Leibnitz-Institute of Atmospheric Physics, University Rostock December 12, 2019

    On Friday, 13 December 2019 at 14.00 Almut Gassman will talk on the topic “Entropy production in numerical modeling of a moist atmosphere”

    Location: Weyertal 86 – 90, Seminar room 1 (ground floor), 50931 Cologne

    Abstract: Atmospheric models are based on basic physical laws, like mass conservation, momentum conservation or sometimes energy conservation. In this talk, another law is inspected, namely the second law of thermodynamics which says that the internal entropy production has to be positive. This law has hardly been checked by modellers. But this law determines the direction (but not the strength) of subgrid-scale parameterized fluxes. Furthermore, the choice of numerical schemes for advection is also limited, because the inherent diffusion of those schemes must be downgradient with respect to the gradients of special measures, which are determined by the second law. The formulation of energetically and entropically consistent numerical schemes will be outlined. Examples of failure of contemporary schemes with respect to the second law will be shown. For instance, conventional (non-entropy-consistent) heat flux parameterizations in the mesosphere will amplify waves instead of allowing them to break and dissipate. Or, as another example, entropy-law-inconsistent higher order or TVD advection schemes for temperature might give rise to accelerations of the wind in the wrong direction. When inspecting the different entropy production (or dissipation) rates in the atmosphere, four types of those rates are distinguished, namely (i) dissipation due to friction, (ii) dissipation due to heat fluxes, (iii) dissipation due to mixing of constituents, and (iv) dissipation due to phase changes. The most important among them are the dissipation due to friction and the dissipation due to falling rain, which can be seen as a special case of mixing of moist air and precipitation.


     

  • Snapshot of the month: Supersonic 3D Blob November 15, 2019

    Supersonic 3D Blob Test computed with a ‘High-order Discontinuous Galerkin/Finite-Volume Hybrid Scheme’ on a non-conforming cartesian grid.


     

  • New paper submitted: FLEXI: A high order discontinuous Galerkin framework for hyperbolic-parabolic conservation laws October 25, 2019

    High order (HO) schemes are attractive candidates for the numerical solution of multiscale problems occurring in fluid dynamics and related disciplines. Among the HO discretization variants, discontinuous Galerkin schemes offer a collection of advantageous features which have lead to a strong increase in interest in them and related formulations in the last decade. The methods have matured sufficiently to be of practical use for a range of problems, for example in direct numerical and large eddy simulation of turbulence. However, in order to take full advantage of the potential benefits of these methods, all steps in the simulation chain must be designed and executed with HO in mind. Especially in this area, many commercially available closed-source solutions fall short. In this work, we therefor present the FLEXI framework, a HO consistent, open-source simulation tool chain for solving the compressible Navier-Stokes equations in a high performance computing setting. We describe the numerical algorithms and implementation details and give an overview of the features and capabilities of all parts of the framework. Beyond these technical details, we also discuss the important, but often overlooked issues of code stability, reproducibility and user-friendliness. The benefits gained by developing an open-source framework are discussed, with a particular focus on usability for the open-source community. We close with sample applications that demonstrate the wide range of use cases and the expandability of FLEXI and an overview of current and future developments.

    Link: Paper on arXiv


     

  • Workshop on Efficiency in Computational Science, Cologne, Sep 25th, 2019 September 17, 2019

    On Wednesday, September 25th, 2019, the Workshop on Efficiency in Computational Science will take place at the Department of Mathematics and Computer Science. This workshop brings together researchers in the field of computational science to share and discuss their work where it relates to efficiency. Here, efficiency is understood in a very broad sense, including numerical method development, serial and parallel algorithms, implementation, and hardware aspects. The intention is to provide an informal environment that encourages the exchange of novel ideas and untested approaches. Therefore, the speakers are asked to put an emphasis on work in progress and unsolved issues, and to not restrict themselves to sharing “camera-ready” results only. Ultimately, the goal is to get a fresh perspective on common challenges, to establish new connections across institutional and discipline boundaries, and to identify potential for future scientific collaborations.

    Agenda
    13:00 Welcome
    13:10 Efficiency challenges in adaptive parallel multiphysics simulations
    13:40 Current HPC developments in the TRACE flow solver
    14:10 Towards Large Scale Continual Learning on Modular High Performance Computers
    14:40 Coffee break
    15:20 Vectorization of high-order DG and adaptive linearization
    15:50 Promises and Challenges of Dispersion Relation Preserving Finite Difference Methods
    16:20 Structural modelling for helicopter simulation – or: making small problems even smaller
    16:50 Coffee, discussions & open end

    Full agenda (PDF, 499 KB)


     

  • Viktor Linders from Lund University to visit NumSim group for two weeks September 12, 2019

    From September 16th through September 27th, Dr. Viktor Linders from Lund University is going to visit the NumSim group by invitation of Dr. Michael Schlottke-Lakemper. Together with researchers from RWTH Aachen University, they are working on high-order methods with optimized dispersive properties for aeroacoustics simulations. During his stay in Cologne, Dr. Linders will share some of his previous results on summation-by-parts methods in a talk at the Workshop on Efficiency in Computational Science.