Anna-Karin Tornberg

Stefan Engblom, Minh Do-Quang, Gustav Amberg et al.

An existing phase-field model of two immiscible fluids with a single soluble surfactant present is discussed in detail. We analyze the well-posedness of the model and provide strong evidence that it is mathematically ill-posed for a large set of physically relevant parameters. As a consequence, critical modifications to the model are suggested that substantially increase the domain of validity. Carefully designed numerical simulations offer informative demonstrations as to the sharpness of our theoretical results and the qualities of the physical model. A fully coupled hydrodynamic test-case demonstrates the potential to capture also non-trivial effects on the overall flow.

Ludvig af Klinteberg, Anna-Karin Tornberg

Boundary integral methods are advantageous when simulating viscous flow around rigid particles, due to the reduction in number of unknowns and straightforward handling of the geometry. In this work we present a fast and accurate framework for simulating spheroids in periodic Stokes flow, which is based on the completed double layer boundary integral formulation. The framework implements a new method known as quadrature by expansion (QBX), which uses surrogate local expansions of the layer potential to evaluate it to very high accuracy both on and off the particle surfaces. This quadrature method is accelerated through a newly developed precomputation scheme. The long range interactions are computed using the spectral Ewald (SE) fast summation method, which after integration with QBX allows the resulting system to be solved in M log M time, where M is the number of particles. This framework is suitable for simulations of large particle systems, and can be used for studying e.g. porous media models.

Chiara Sorgentone, Anna-Karin Tornberg, Petia Vlahovska

We present a highly accurate numerical method based on a boundary integral formulation and the leaky dielectric model to study the dynamics of surfactant-covered drops in the presence of an applied electric field. The method can simulate interacting 3D drops (no axisymmetric simplification) in close proximity, can consider different viscosities, is adaptive in time and able to handle substantial drop deformation. For each drop global representations of the variables based on spherical harmonics expansions are used and the spectral accuracy is achieved by designing specific numerical tools: a specialized quadrature method for the singular and nearly singular integrals that appear in the formulation, a general preconditioner for the implicit treatment of the surfactant diffusion and a reparametrization procedure able to ensure a high-quality representation of the drops also under deformation. Our numerical method is validated against theoretical, numerical and experimental results available in the literature, as well as a new second-order theory developed for a surfactant-laden drop placed in a quadrupole electric field.

Ludvig af Klinteberg, Davoud Saffar Shamshirgar, Anna-Karin Tornberg

We present a spectrally accurate method for the rapid evaluation of free-space Stokes potentials, i.e. sums involving a large number of free space Green's functions. We consider sums involving stokeslets, stresslets and rotlets that appear in boundary integral methods and potential methods for solving Stokes equations. The method combines the framework of the Spectral Ewald method for periodic problems, with a very recent approach to solving the free-space harmonic and biharmonic equations using fast Fourier transforms (FFTs) on a uniform grid. Convolution with a truncated Gaussian function is used to place point sources on a grid. With precomputation of a scalar grid quantity that does not depend on these sources, the amount of oversampling of the grids with Gaussians can be kept at a factor of two, the minimum for aperiodic convolutions by FFTs. The resulting algorithm has a computational complexity of O(N log N) for problems with N sources and targets. Comparison is made with a fast multipole method (FMM) to show that the performance of the new method is competitive.

Ludvig af Klinteberg, Anna-Karin Tornberg

When solving partial differential equations using boundary integral equation methods, accurate evaluation of singular and nearly singular integrals in layer potentials is crucial. A recent scheme for this is quadrature by expansion (QBX), which solves the problem by locally approximating the potential using a local expansion centered at some distance from the source boundary. In this paper we introduce an extension of the QBX scheme in 2D denoted AQBX - adaptive quadrature by expansion - which combines QBX with an algorithm for automated selection of parameters, based on a target error tolerance. A key component in this algorithm is the ability to accurately estimate the numerical errors in the coefficients of the expansion. Combining previous results for flat panels with a procedure for taking the panel shape into account, we derive such error estimates for arbitrarily shaped boundaries in 2D that are discretized using panel-based Gauss-Legendre quadrature. Applying our scheme to numerical solutions of Dirichlet problems for the Laplace and Helmholtz equations, and also for solving these equations, we find that the scheme is able to satisfy a given target tolerance to within an order of magnitude, making it useful for practical applications. This represents a significant simplification over the original QBX algorithm, in which choosing a good set of parameters can be hard.

Michael Siegel, Anna-Karin Tornberg

Accurate evaluation of layer potentials is crucial when boundary integral equation methods are used to solve partial differential equations. Quadrature by expansion (QBX) is a recently introduced method that can offer high accuracy for singular and nearly singular integrals, using truncated expansions to locally represent the potential. The QBX method is typically based on a spherical harmonics expansion which when truncated at order $p$ has $O(p^2)$ terms. This expansion can equivalently be written with $p$ terms, however paying the price that the expansion coefficients will depend on the evaluation/target point. Based on this observation, we develop a target specific QBX method, and apply it to Laplace's equation on multiply connected domains. The method is local in that the QBX expansions only involve information from a neighborhood of the target point. An analysis of the truncation error in the QBX expansions is presented, practical parameter choices are discussed and the method is validated and tested on various problems.

Fredrik Fryklund, Erik Lehto, Anna-Karin Tornberg

We introduce an efficient algorithm, called partition of unity extension or PUX, to construct an extension of desired regularity of a function given on a complex multiply connected domain in $2D$. Function extension plays a fundamental role in extending the applicability of boundary integral methods to inhomogeneous partial differential equations with embedded domain techniques. Overlapping partitions are placed along the boundaries, and a local extension of the function is computed on each patch using smooth radial basis functions; a trivially parallel process. A partition of unity method blends the local extrapolations into a global one, where weight functions impose compact support. The regularity of the extended function can be controlled by the construction of the partition of unity function. We evaluate the performance of the PUX method in the context of solving the Poisson equation on multiply connected domains using a boundary integral method and a spectral solver. With a suitable choice of parameters the error converges as a tenth order method down to $10^{-14}$.

Ludvig af Klinteberg, Anna-Karin Tornberg

In boundary integral methods it is often necessary to evaluate layer potentials on or close to the boundary, where the underlying integral is difficult to evaluate numerically. Quadrature by expansion (QBX) is a new method for dealing with such integrals, and it is based on forming a local expansion of the layer potential close to the boundary. In doing so, one introduces a new quadrature error due to nearly singular integration in the evaluation of expansion coefficients. Using a method based on contour integration and calculus of residues, the quadrature error of nearly singular integrals can be accurately estimated. This makes it possible to derive accurate estimates for the quadrature errors related to QBX, when applied to layer potentials in two and three dimensions. As examples we derive estimates for the Laplace and Helmholtz single layer potentials. These results can be used for parameter selection in practical applications.

Sara Pålsson, Michael Siegel, Anna-Karin Tornberg

Performing highly accurate simulations of droplet systems is a challenging problem. This is primarily due to the interface dynamics which is complicated further by the addition of surfactants. This paper presents a boundary integral method for computing the evolution of surfactant-covered droplets in 2D Stokes flow. The method has spectral accuracy in space and the adaptive time-stepping scheme allows for control of the temporal errors. Previously available semi-analytical solutions (based on conformal-mapping techniques) are extended to include surfactants, and a set of algorithms is introduced to detail their evaluation. These semi-analytical solutions are used to validate and assess the accuracy of the boundary integral method, and it is demonstrated that the presented method maintains its high accuracy even when droplets are in close proximity.

Shriram Srinivasan, Anna-Karin Tornberg

Recently, Gimbutas et al derived an elegant representation for the Green's functions of Stokes flow in a half-space. We present a fast summation method for sums involving these half-space Green's functions (stokeslets, stresslets and rotlets) that consolidates and builds on the work by Klinteberg et al for the corresponding free-space Green's functions. The fast method is based on two main ingredients: The Ewald decomposition and subsequent use of FFTs. The Ewald decomposition recasts the sum into a sum of two exponentially decaying series: one in real-space (short-range interactions) and one in Fourier-space (long-range interactions) with the convergence of each series controlled by a common parameter. The evaluation of short-range interactions is accelerated by restricting computations to neighbours within a specified distance, while the use of FFTs accelerates the computations in Fourier-space thus accelerating the overall sum. We demonstrate that while the method incurs extra costs for the half-space in comparison to the free-space evaluation, greater computational savings is also achieved when compared to their respective direct sums.

Davood Saffar Shamshirgar, Berk Hess, Anna-Karin Tornberg

The smooth particle mesh Ewald (SPME) method is an FFT based method for the fast evaluation of electrostatic interactions under periodic boundary conditions. A highly optimized implementation of this method is available in GROMACS, a widely used software for molecular dynamics simulations. In this article, we compare a more recent method from the same family of methods, the spectral Ewald (SE) method, to the SPME method in terms of performance and efficiency. We consider serial and parallel implementations of both methods for single and multiple core computations on a desktop machine as well as the Beskow supercomputer at KTH Royal Institute of Technology. The implementation of the SE method has been well optimized, however not yet comparable to the level of the SPME implementation that has been improved upon for many years. We show that the SE method is very efficient whenever used to achieve high accuracy and that it already at this level of optimization can be competitive for low accuracy demands.

Erik Boström, Anna-Karin Tornberg, Ludvig af Klinteberg

Fast Ewald summation efficiently evaluates Coulomb interactions and is widely used in molecular dynamics simulations. It is based on a split into a short-range and a long-range part, where evaluation of the latter is accelerated using the fast Fourier transform (FFT). The accuracy and computational cost depend critically on the mollifier in the kernel split and the window function used in the spreading and interpolation steps that enable the use of the FFT. The first prolate spheroidal wavefunction (PSWF) has optimal concentration in real and Fourier space simultaneously, and is used when defining both a mollifier and a window function. We provide a complete description of the method and derive rigorous error estimates. In addition, we obtain closed-form approximations of the Fourier truncation and aliasing errors, yielding explicit parameter choices for the achieved error to closely match the prescribed tolerance. Numerical experiments confirm the analysis: PSWF-based Ewald summation achieves a given accuracy with significantly fewer Fourier modes and smaller window supports than Gaussian- and B-spline-based approaches, providing a superior alternative to existing Ewald methods for particle simulations.

Dag Lindbo, Anna-Karin Tornberg

A new method for interface tracking is presented. The interface representation, based on domain decomposition, provides the interface location explicitly, yet is Eulerian. This allows for well established finite difference methods on uniform grids to be used for the numerics of advecting the interface and other computations. CFL-stable and second order accurate explicit and implicit time-stepping methods are derived. Numerical results are given to substantiate stated convergence properties, as well as convergence in interface curvature and mass conservation. Our method is applied to a boundary integral formulation for Stokes flow, and the resulting integrals are analyzed and treated numerically to second order accuracy. Finally, we embed our method in the familiar immersed boundary- and immersed interface methods for two-phase Navier-Stokes flow.

Chiara Sorgentone, Anna-Karin Tornberg

The presence of surfactants alters the dynamics of viscous drops immersed in an ambient viscous fluid. This is specifically true at small scales, such as in applications of droplet based microfluidics, where the interface dynamics become of increased importance. At such small scales, viscous forces dominate and inertial effects are often negligible. Considering Stokes flow, a numerical method based on a boundary integral formulation is presented for simulating 3D drops covered by an insoluble surfactant. The method is able to simulate drops with different viscosities and close interactions, automatically controlling the time step size and maintaining high accuracy also when substantial drop deformation appears. To achieve this, the drop surfaces as well as the surfactant concentration on each surface are represented by spherical harmonics expansions. A novel reparameterization method is introduced to ensure a high-quality representation of the drops also under deformation, specialized quadrature methods for singular and nearly singular integrals that appear in the formulation are evoked and the adaptive time stepping scheme for the coupled drop and surfactant evolution is designed with a preconditioned implicit treatment of the surfactant diffusion.

Davoud Saffar Shamshirgar, Anna-Karin Tornberg

We present a fast and spectrally accurate method for efficient computation of the three dimensional Coulomb potential with periodicity in one direction. The algorithm is FFT-based and uses the so-called Ewald decomposition, which is naturally most efficient for the triply periodic case. In this paper, we show how to extend the triply periodic Spectral Ewald method to the singly periodic case, such that the cost of computing the singly periodic potential is only marginally larger than the cost of computing the potential for the corresponding triply periodic system. In the Fourier space contribution of the Ewald decomposition, a Fourier series is obtained in the periodic direction with a Fourier integral over the non periodic directions for each discrete wave number. We show that upsampling to resolve the integral is only needed for modes with small wave numbers. For the zero wave number, this Fourier integral has a singularity. For this mode, we effectively need to solve a free-space Poisson equation in two dimensions. A very recent idea by Vico et al. makes it possible to use FFTs to solve this problem, allowing us to unify the treatment of all modes. An adaptive 3D FFT can be established to apply different upsampling rates locally. The computational cost for other parts of the algorithm is essentially unchanged as compared to the triply periodic case, in total yielding only a small increase in both computational cost and memory usage for this singly periodic case.