Snapshot: Chaotic behavior of the Kelvin-Helmholtz Instability

We investigate the Kelvin-Helmholtz instability setup of Fjordholm et al. [1], where the sharp interface in the initial condition is randomly perturbed. The simulation is done with an entropy stable DGSEM with polynomial degree N=3 and 32 x 32 grid cells. The final time is t = 2.0. Trixi.jl [2] is used with the setup found in examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl.

The following gif is showing 100 simulation results of the density distribution for different random interface perturbations:

Here is the average density distribution of all 100 results:

[1] Ulrik S. Fjordholm, Roger Käppeli, Siddhartha Mishra, Eitan Tadmor. Construction of approximate entropy measure valued solutions for hyperbolic systems of conservation laws, 2014. (https://arxiv.org/abs/1402.0909)
[2] Trixi.jl, a numerical simulation framework for hyperbolic conservation laws written in Julia. (https://github.com/trixi-framework/Trixi.jl)

Thesis Snapshot: Planning of collision-free trajectories for UAVs and optimization using uniform B-splines. (Bachelor thesis by Alexander Bast)

Autonomous vehicles will be an integral part of our environment and will have a significant impact on transportation and other industries. In this context, drones, so-called UAVs will be used for aerial transportation. Due to the high information density, motion planning is a fundamental problem in the autonomous movement of vehicles.
In the future, drones should transport goods or perhaps people in a fast, safe and thus collision-free manner. In doing so, the drone should be able to perform all calculations quickly, accurately and autonomously. Therefore, it makes sense to develop algorithms that meet these requirements.
As part of a bachelor thesis, a method [1] was investigated that uses uniform B-splines to replan a drone’s trajectory in the presence of a potential collision. This method formulates the task of replanning the trajectory as the solution to a non-linear optimization problem.
\begin{align*}
E=\lambda_{ep}(p(t_{ep})-p_{ep})^{2}+\lambda_{ep}(p'(t_{ep})-p’_{ep})^{2} + \lambda_{c} \int_{t_{min}}^{t_{max}}c(p(t))||p'(t)||dt + \sum_{i=2}^{4}\int_{t_{min}}^{t_{max}}\lambda_{q}(p^{(i)}(t))^{2}dt
\end{align*}
With this simplified equation, three cost terms are considered: Endpoint error, collision cost and quadratic derivative cost. With the goal of computing an efficient and collision-free trajectory, an iterative process is used to construct a B-spline curve such that it follows a pre-planned and optimal trajectory. The focus of the analysis was on the selection of the control parameters λep, λc, λq, as well as the number of control points to be optimized and the resulting properties of the optimized trajectory. Using a uniform knot vector, the B-spline is expanded to include control points whose spatial coordinates are determined using nonlinear optimization so that the objective function is minimized. Due to the properties of B-splines, upcoming segments can be locally optimized so that a collision is prevented. The use of a uniform knot vector allows an efficient replanning of the trajectory in real time, which can be performed autonomously by a drone. The following figure shows one result.

[1] Vladyslav Usenko, Lukas von Stumberg, Andrej Pangercic, Daniel Cremers, Real-Time Trajectory Replanning for MAVs using Uniform B-splines and a 3D Circular Buffer, Technical University of Munich, 2017, https://arxiv.org/abs/1703.01416.

Thesis Snapshot: The effect of sea ice seasonality on the two-dimensional energy balance model ‘Klimakoffer’ (Bachelor thesis by Xenia Ibach)

In a previous snapshot, the two-dimensional energy balance model has been already presented:\begin{align*}
C(\hat{r})\frac{\partial T(\hat{r},t)}{\partial t} +A+BT(\hat{r},t)= \nabla\cdot(D(\hat{r})\nabla T(\hat{r},t)) + QS(\hat{r},t)(1-a(\hat{r}))
\end{align*}
The model is based on the assumption of a static surface distribution of the surface types land, ocean, sea ice and snow. However, climate change at an increasing rate is also affecting surface types seasonally as well as permanently over the years. Melting ice and snow masses worldwide are a consequence of this, affecting not only snow and sea ice distribution, but also that of oceans and land surfaces through rising sea ice levels.
Therefore, for more accurate climate modeling with this EBM, it is interesting to modify the model with respect to a time-dependent surface distribution. This modification directly affects the model parameters albedo, solar forcing, and heat capacity, requiring the following reformulation of the model:
\begin{align*}
C(\hat{r},t)\frac{\partial T(\hat{r},t)}{\partial t} +A+BT(\hat{r},t)= \nabla\cdot(D(\hat{r})\nabla T(\hat{r},t)) + QS(\hat{r},t)(1-a(\hat{r},t))
\end{align*}
As part of a bachelor’s thesis, two options were presented for introducing a seasonal sea ice distribution into the EBM by using and modifying simulations from the ‘Klimakoffer‘. Similar approaches are also conceivable for the other surface types.
One possibility is to process documented measurements of sea ice distribution, available for example from institutions such as the NSIDC (National Snow and Sea Ice Data Center) or NASA, so that a monthly change in sea ice distribution at the edge of the sea ice extent at the two poles can be mapped using a specially developed algorithm.

A second option is to directly read in satellite images and categorize the color values of the image’s pixels and map them to the surface types required in the EBM. Here, as an example, a satellite image of the Blue Marble series from 2004 of NASA was read in and categorized on the basis of the color values to the different surface types.

Image from NASA’s Blue Marble Next Generation [3]


[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.
[3] EOS Project Science Office at NASA Goddard Space Flight Center, Blue Marble Next Generation, 2004, https://visibleearth.nasa.gov/collection/1484/blue-marble.

Snapshot: A high order detonation diffraction simulation with reactive Euler equations

We show a detonation diffraction simulation with multiple obstacles in a reflective box using reactive multicomponent Euler equations with Trixi.jl. These kind of detonation waves may result in pressure and density drops close to zero which is numerically difficult not only but especially for high order schemes.

In this simulation we use a reacting model which consists of one reaction and two species. This example has been calculated using the high order DG-FV blending method in Trixi.jl with a HOHQMesh generated mesh whereas the chemical network is solved with KROME.

(Note: Not all of these features are currently available in the main code of Trixi.jl.)

Snapshot: Supersonic turbulence simulations with FLUXO

We show a supersonic turbulence simulation in a periodic box fueled by perpetual kinetic energy injection on the largest three modes. Such setups are subject of active research in astrophysics and are believed to play an important role in the dynamics of galactic clouds and star formation, e.g., [1,2,3]. An initially uniform density distribution (left plot in the video) is gradually driven to a turbulent state till an average velocity of Mach 3.5 is reached (about t=0.19). Locally, speeds can spike above Mach 35 (bottom right plot in the video). The simulation was computed on FLUXO (https://github.com/project-fluxo/fluxo) with a robust invariant domain preserving Discontinuous Galerkin scheme (DG) [4] solving the compressible, quasi-isothermal ($\gamma =1.001$) Euler equations on a uniform grid with $128^3$ degrees-of-freedom.

Smooth blending with a monotone finite volume scheme stabilizes the high-order DG at strongly shocked regions (filaments of high density concentration) displayed by the blending parameter plot (top right in the video). Blue encodes regions of 100% high order DG while red dots indicate focused blending with FV. Hypersonic turbulence regimes are very challenging for any numerical scheme considering the stark gradients in density, bubbles of near-vacuum and the extremely high speeds involved in such simulations.

[1] https://academic.oup.com/mnras/article/436/2/1245/1126116
[2] https://adsabs.harvard.edu/full/1981MNRAS.194..809L]
[3] https://academic.oup.com/mnras/article/480/3/3916/5060766
[4] Rueda-Ramírez, A. M., Pazner, W., & Gassner, G. J. (2022). Subcell limiting strategies for discontinuous Galerkin spectral element methods. arXiv preprint arXiv:2202.00576.

Snapshot: Hypersonic astrophysical jet (Mach~2000) simulation using a hybrid DG/FV solver

Simulation of a hypersonic astrophysical jet with Mach=2156.91 moving through a medium at rest.

The simulation was performed on a 256 × 256 quadrilateral grid with a nodal fourth-order entropy-stable Discontinuous Galerkin (DG) scheme. Local bounds on the density and specific entropy are imposed by blending the DG scheme locally (at the node level) with a second-order Finite Volume (FV) method. Both the DG and the FV methods use the Harten-Lax-van Leer-Einfeldt (HLLE) Riemann solver.

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

The numerical methods used to perform the simulation are described in our preprint:
Rueda-Ramírez, A. M., Pazner, W., Gassner, G.J. (2022). Subcell Limiting Strategies for Discontinuous Galerkin Spectral Element Methods. Submitted. https://arxiv.org/abs/2202.00576.

Snapshot: Immersed Boundary Discontinuous Galerkin Simulation of Flow Past a Cylinder at Re=100

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.

Snapshot: Simulation of a stiff multi-species reacting flow

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

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

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: