The quadrature of cut elements is crucial for all Finite Element Methods that do not apply boundary-fitted meshes. It should be efficient, accurate, and robust. Various approaches balancing these requirements have been published, with some available as open-source implementations. This work reviews these open-source codes and the methods used. Furthermore, benchmarking examples are developed for 2D and 3D geometries. Implicit and explicit boundary descriptions are available for all models. The different examples test the efficiency, accuracy, versatility, and robustness of the codes. Special focus is set on the influence of the input parameter, which controls the desired quadrature order, on the actual integration error. A detailed comparison of the discussed codes is carried out. The systematic benchmarking allows a conclusive comparison and presents a valuable tool for future code development. All tests are published in an accompanying open-source repository.

In the field of scientific computing, one often finds several alternative software packages (with open or closed source code) for solving a specific problem. These packages sometimes even use alternative methodological approaches, e.g., different numerical discretizations. If one decides to use one of these packages, it is often not clear which one is the best choice. To make an informed decision, it is necessary to measure the performance of the alternative software packages for a suitable set of test problems, i.e. to set up a benchmark. However, setting up benchmarks ad-hoc can become overwhelming as the parameter space expands rapidly. Very often, the design of the benchmark is also not fully set at the start of some project. For instance, adding new libraries, adapting metrics, or introducing new benchmark cases during the project can significantly increase complexity and necessitate laborious re-evaluation of previous results. This paper presents a proven approach that utilizes established Continuous Integration tools and practices to achieve high automation of benchmark execution and reporting. Our use case is the numerical integration (quadrature) on arbitrary domains, which are bounded by implicitly or parametrically defined curves or surfaces in 2D or 3D.

We perform high-order simulations of two-phase flows in capillaries, with and without evaporation. Since a sharp-interface model is used, singularities can arise at the three-phase contact line, where the fluid-fluid interface interacts with the capillary wall. These singularities are especially challenging when a highly accurate, high-order method with very little numerical diffusion is used for the flow solver. In this work, we employ the eXtended Discontinuous Galerkin (XDG) method, which has a very high accuracy but a severe limit regarding e.g., the time-step restriction. To address this challenge and enhance the stability of our numerical method we introduce a novel approach for representing a moving interface in the case of two-phase flows. We propose a global analytical representation of the interface-describing level-set field, defined by a small set of time-dependent parameters. Noteworthy for its simplicity and efficiency, this method effectively addresses the inherent complexity of two-phase flow problems. Furthermore, it significantly improves numerical stability and enables the use of larger time steps, ensuring both reliability and computational efficiency in our simulations. We compare different analytic expressions for level-set representation, including the elliptic function and the fourth-order polynomial, and validate the method against established literature data for capillary rise, both with and without evaporation. These results highlight the effectiveness of our approach in resolving complex interfacial dynamics.

We investigate the nonlinear axisymmetric shape oscillations of inviscid droplets, focusing on the modal coupling between different oscillation modes. The research builds upon the geometrically exact nonlinear theory for droplet shape oscillations by Pl & uuml;macher et al. [Phys. Fluids 32, 067104 (2020)], which allows for large deformations and superimposed initial deformations. The theory rests on the potential flow assumption and the usual momentum jump condition at the droplet surface including the surface tension. Using the unified transform method of A. S. Fokas [Proc. R. Soc. London Ser. A 453, 1411 (1997)], the governing equations are transformed into a system of integrodifferential equations defined on a unit sphere. These equations are then solved numerically with a Galerkin method that uses spherical harmonics and exhibits overall spectral convergence. The frequency shift and the extended time spent in the prolate state within an oscillation period are analyzed for different initial deformation modes and show good agreement with existing literature. The results demonstrate that for large initial deformations, the modal coupling significantly influences the oscillation dynamics, particularly in the interaction between even and odd modes. To study modal coupling in detail, superimposed initial deformations are investigated in which a large-amplitude even mode is superimposed with various individual small-amplitude odd modes. It is shown that the initially excited large-amplitude even mode remains nearly unaffected by the superimposed small-amplitude odd modes. The initially superimposed small-amplitude odd mode, however, is altered by the large-amplitude even mode without being significantly amplified or damped. This nonlinear modal coupling is not present in single even initial deformation modes, where it is shown that an even initial deformation mode cannot excite odd modes.

In this work, a cell agglomeration strategy for the cut cells arising in the eXtended discontinuous Galerkin (XDG) method is presented. Cut cells are a fundamental aspect of unfitted mesh approaches, where complex geometries or interfaces separating subdomains are embedded into structured background grids to facilitate the mesh generation process. In such methods, arbitrary small cells occur due to the intersections of background cells with embedded geometries and lead to discretization difficulties due to their diminutive sizes. Furthermore, temporal evolutions of these geometries may lead to topological changes across different time steps. Both of these issues, that is, small-cut cells and topological changes, can be addressed with a cell agglomeration technique, independent of discretization. However, cell agglomeration encounters significant difficulties in three dimensions due to the complexity of neighborship and issues like cycles and parallel agglomeration chains. The proposed strategy introduces a robust framework that mitigates these problems by incorporating methods for cycle prevention, chain agglomeration, and parallelization. Implemented in the open-source software package BoSSS, this strategy has been successfully tested on multiprocessor systems using dynamic multiphase test cases in both two and three dimensions, enabling simulations that were previously infeasible.

Multigrid methods have been a popular approach for solving linear systems arising from the discretization of partial differential equations (PDEs) for several decades. They are particularly effective for accelerating convergence rates with optimal complexity in terms of both time and space. K-cycle orthonormalization multigrid is a robust variant of the multigrid method that combines the efficiency of multigrid with the robustness of Krylov-type residual minimalizations for problems with strong anisotropies. However, traditional implementations of K-cycle orthonormalization multigrid often rely on bulk-synchronous parallelism, which can limit scalability on modern high-performance computing (HPC) systems. This paper presents a task- parallel variant of the K-cycle orthonormalization multigrid method that leverages asynchronous execution to improve scalability and performance on large-scale parallel systems.

In the sharp-interface modeling of three phase contact lines special care has to be taken with respect to boundary and interface conditions, to avoid singularities. In this work, one such singularity, arising from an incompatibility in a standard model for evaporation/condensation at the contact line, is investigated numerically. By means of employing an extended Discontinuous Galerkin method, established in preliminary works, we study the severity of the singularity in the model. In a first step, the experimental order of convergence is determined in absence of evaporation/condensation and thus without the incompatibility under investigation. The extent of the incompatibility is then measured by computing the convergence orders again, repeating the same simulation with evaporation/condensation present. As one possible model alteration to dispose of the incompatibility causing the singularity and to assess the properties of the numerical method we propose the introduction of a slip condition on the fluid–fluid interface. This part of the study is complemented by an examination of the experimental order of convergence for the test case using the model without and with slip on the interface.

In this work, we present an extended Discontinuous Galerkin method for simulating transient, incompressible two-phase flows, which include heat transfer and thermally driven evaporation at the interface of single-component systems. This expands our previous work to include the consideration of non-material interfaces and a coupling between velocities and temperature gradients at the interface. The phase boundary is represented by the zero-set of a level set function, while effects due to surface tension are treated by the Laplace-Beltrami formulation. This sharp interface model allows for a sub-grid accurate representation of the solution fields. By using compactly supported polynomial solutions, discontinuities at the interface can be sharply represented without employing additional reconstruction schemes. The approach is validated through well-known evaporation test cases. This includes two 1D test cases, known as Stefan and Sucking problem, a 2D film boiling and finally the 3D growth of a vapor bubble, known as Scriven test case.

Nonlinear axisymmetric shape oscillations of a Newtonian drop in a vacuum are investigated using two different theoretical methods, for fundamental interest and for the significance of the oscillations in transport processes across the drop surface. The extended discontinuous Galerkin method is contrasted to the weakly nonlinear theory. While the former allows large drop surface deformation amplitudes to be analyzed with high precision and drop volume errors below 0.11% even at the largest deformations, the latter provides analytical insight into the origin of quasiperiodic time behavior of the oscillations and reveals the oscillation modes coupled in the nonlinear motion. Results from both methods for moderate initial deformation amplitudes at modes of initial drop deformation m=2, 3, and 4 are in excellent agreement, showing the time asymmetry of the oscillation and the decrease of the oscillation frequency with increasing deformation amplitude. The Fourier power spectra for the first oscillation period exhibit decreased dominant frequencies as compared to the linear results as well as the mode coupling as nonlinear effects. The numerical method is used to compute the oscillatory and damping behavior of viscous drops, as well as the interconversion of kinetic and surface energies during the oscillations at strong initial deformations. Published by the American Physical Society 2024

In this paper, we introduce a novel high-order shock tracking method and provide a proof of concept. Our method leverages concepts from implicit shock tracking and extended discontinuous Galerkin methods, primarily designed for solving partial differential equations featuring discontinuities. To address this challenge, we solve a constrained optimization problem aiming at accurately fitting the zero iso-contour of a level set function to the discontinuities. Additionally, we discuss various robustness measures inspired by both numerical experiments and existing literature. Finally, we showcase the capabilities of our method through a series of two-dimensional problems, progressively increasing in complexity.

We present a high-order method that provides numerical integration on volumes, surfaces, and lines defined implicitly by two smooth intersecting level sets. To approximate the integrals, the method maps quadrature rules defined on hypercubes to the curved domains of the integrals. This enables the numerical integration of a wide range of integrands since integration on hypercubes is a well known problem. The mappings are constructed by treating the isocontours of the level sets as graphs of height functions. Numerical experiments with smooth integrands indicate a high-order of convergence for transformed Gauss quadrature rules on domains defined by polynomial, rational, and trigonometric level sets. We show that the approach we have used can be combined readily with adaptive quadrature methods. Moreover, we apply the approach to numerically integrate on difficult geometries without requiring a low-order fallback method.

A numerical simulation of the fluid flow in the gravure printing nip, based on a discontinuous Galerkin algorithm, is used to study the fluid-splitting process and the transition between point and lamella splitting. We study the pressure and shear singularities at the contact point of the printing cylinder and substrate as a function of the variable microscopic residual gap and variations of the printing fluid quantities introduced to the nip. As the hydrodynamic boundary value problem is ill-defined by the nip singularity, we enhance the simulation using renormalization group and algebraic scaling techniques in order to obtain a numerically stable and physically meaningful prediction. Our simulations are compared to analytical results from lubrication theory and to experimental observations on a gravure press.

In this article, we present the foam-dg project, which provides a bridge between OpenFOAM(R) and the high-order DG (discontinuous Galerkin) framework BoSSS. Thanks to the flexibility of the coupling approach, mixed calculations where some parts of the equation system are solved in OpenFOAM(R) and others are solved in BoSSS are easily possible. This is showcased using the convective Cahn-Hilliard equation, where the Cahn-Hilliard part is solved in BoSSS and the Navier-Stokes part is solved in OpenFOAM(R). The obtained results appear reasonable, though the main focus of this paper is to present and document the foam-dg project rather than on quantitative results.

We present a fully coupled solver based on the discontinuous Galerkin method for steady‐state diffusion flames using the low‐Mach approximation of the governing equations with a one‐step kinetic model. The nonlinear equation system is solved with a Newton–Dogleg method and initial estimates for flame calculations are obtained from a flame‐sheet model. Details on the spatial discretization and the nonlinear solver are presented. The method is tested with reactive and nonreactive benchmark cases. Convergence studies are presented, and we show that the expected convergence rates are obtained. The solver for the low‐Mach equations is used for calculating a differentially heated cavity configuration, which is validated against benchmark solutions. Additionally, a two‐dimensional counter diffusion flame is calculated, and the results are compared with the self‐similar one dimensional solution of said configuration.

In this work a solver for two-dimensional, instationary two-phase flows on the basis of the extended discontinuous Galerkin (extended DG/XDG) method is presented. The XDG method adapts the approximation space conformal to the position of the interface. This allows a subcell accurate representation of the incompressible Navier-Stokes equations in their sharp interface formulation. The interface is described as the zero set of a signed-distance level-set function and discretized by a standard DG method. For the interface, resp. level-set, evolution an extension velocity field is used and a two-staged algorithm is presented for its construction on a narrow-band. On the cut-cells a monolithic elliptic extension velocity method is adapted and a fast-marching procedure on the neighboring cells. The spatial discretization is based on a symmetric interior penalty method and for the temporal discretization a moving interface approach is adapted. A cell agglomeration technique is utilized for handling small cut-cells and topology changes during the interface motion. The method is validated against a wide range of typical two-phase surface tension driven flow phenomena in a 2D setting including capillary waves, an oscillating droplet and the rising bubble benchmark.

In this paper, we are going to present a high-order shock fitting approach based on a cut-cell method. We formulate a suitable Constraint Optimization Problem and develop an algorithm aiming to reconstruct the shock front represented by the zero iso-contour of a Level Set function.

In this work, an extended discontinuous Galerkin (extended DG/XDG also called unfitted DG) solver for two‐dimensional flow problems exhibiting moving contact lines is presented. The generalized Navier boundary condition is employed within the XDG discretization for the handling of the moving contact lines. The spatial discretization is based on a symmetric interior penalty method and the numerical treatment of the surface tension force is done via the Laplace–Beltrami formulation. The XDG method adapts the approximation space conformal to the position of the interface and allows a sub‐cell accurate representation within the sharp interface formulation. The interface is described as the zero set of a signed‐distance level‐set function and discretized by a standard DG method. No adaption of the level‐set evolution algorithm is needed for the extension to moving contact line problems. The developed solver is validated against typical two‐dimensional contact line driven flow phenomena including droplet simulations on a wall and the two‐phase Couette flow.

We investigate the molecular origin of shear-thinning in melts of flexible, semiflexible and rigid oligomers with coarse-grained simulations of a sheared melt. Alignment, stretching and tumbling modes or suppression of the latter all contribute to understanding how macroscopic flow properties emerge from the molecular level. By performing simulations of single chains in a shear flow, we identify which of these phenomena are of collective nature and arise through interchain interactions and which are already present in dilute systems. Building upon these microscopic simulations we identify by means of the Irving-Kirkwood formula the corresponding macroscopic stress tensor for a non-Newtonian polymer fluid. Shear-thinning effects in oligomer melts are also demonstrated by macroscopic simulations of a channel flow. The latter have been obtained by the discontinuous Galerkin method approximating macroscopic polymer flows. Our study confirms the influence of microscopic details in the molecular structure of short polymers such as chain flexibility on macroscopic polymer flows.

We present a high-order discontinuous Galerkin (DG) scheme to solve the system of helically symmetric Navier-Stokes equations which are discussed in [28]. In particular, we discretize the helically reduced Navier-Stokes equations emerging from a reduction of the independent variables such that the remaining variables are: t, r, xi with = az+b phi, where r, phi, z are common cylindrical coordinates and t the time. Beside this, all three velocity components are kept non-zero. A new non-singular coordinate eta is introduced which ensures that a mapping of helical solutions into the three-dimensional space is well defined. Using that, periodicity conditions for the helical frame as well as uniqueness conditions at the centerline axis r=0 are derived. In the sector near the axis of the computational domain a change of the polynomial basis is implemented such that all physical quantities are uniquely defined at the centerline. For the temporal integration, we present a semi-explicit scheme of third order where the full spatial operator is splitted into a Stokes operator which is discretized implicitly and an operator for the nonlinear terms which is treated explicitly. Computations are conducted for a cylindrical shell, excluding the centerline axis, and for the full cylindrical domain, where the centerline is included. In all cases we obtain the convergence rates of order O(h(k+1)) that are expected from DG theory. In addition to the first DG discretization of the system of helically invariant NavierStokes equations, the treatment of the central axis, the resulting reduction of the DG space, and the simultaneous use of a semi-explicit time stepper are of particular novelty.

The software package BoSSS serves the discretization of (steady-state or time-dependent) partial differential equations with discontinuous coefficients and/or time-dependent domains by means of an eXtended Discontinuous Galerkin (XDG, resp. DG) method, aka. cut-cell DG, aka. unfitted DG. This work consists of two major parts: First, the XDG method is introduced and a formal notation is developed, which captures important numerical details such as cell-agglomeration and a multigrid framework. In the second part, iterative solvers for extended DG systems are presented and their performance is evaluated.