A Smooth Partition of Unity Finite Element Method for Vortex Particle Regularization

TL;DR

This paper introduces a novel smooth finite element method on Cartesian grids for vortex particle regularization, enhancing stability and efficiency in vortex simulations.

Contribution

It develops a new $C^ty$-smooth finite element space with a partition of unity approach and a fictitious domain formulation for general domains.

Findings

Stable and convergent scheme with conservation properties

Numerical results confirm theoretical analysis

Optimal grid-size proportional to square root of particle spacing

Abstract

We present a new class of $C^\infty$-smooth finite element spaces on Cartesian grids, based on a partition of unity approach. We use these spaces to construct smooth approximations of particle fields, i.e., finite sums of weighted Dirac deltas. In order to use the spaces on general domains, we propose a fictitious domain formulation, together with a new high-order accurate stabilization. Stability, convergence, and conservation properties of the scheme are established. Numerical experiments confirm the analysis and show that the Cartesian grid-size $\sigma$ should be taken proportional to the square-root of the particle spacing $h$, resulting in significant speed-ups in vortex methods.

Equations

$$\frac{\partial{\omega}}{\partial{t}}+(\mathbf{u}\cdot\nabla)\omega=0\qquad\mbox{in $\Omega$.}$$

$$\omega (t=0)\approx {\omega }_h(0)=\sum _{i=1}^N{\Gamma }_i\delta (x-x_i(0)),$$

$$\Vert \omega -{\omega }_h{\Vert }_{W^{-(m+1),2}(\Omega )}\le C{h}^{m+1}\Vert \omega {\Vert }_{W^{m+1,2}(\Omega )},$$

$$\mathbf{u}=\mathbf{K}\star \omega ,\mathbf{K}(x)=\frac{1}{2\pi }\frac{(-x_2,x_1{)}^{\top }}{|x{|}^2},$$

$$\Vert \mathbf{K}\star \omega {\Vert }_{W^{1,2}(\Omega )}\le C\Vert \omega {\Vert }_{L^2(\Omega )}.$$

$$\exists M\in \mathbb{N}:\forall x\in \Omega :\mathrm{card}\{i|x\in {\Omega }_i\}\le M.$$

$$\mathrm{supp}{\phi }_i$$

$$\sum _i{\phi }_i(x)$$

$$|{\phi }_i{|}_{W^{k,\infty }({\mathbb{R}}^d)}$$

$$|f|_{W^{k,p}(\mathbb{R}^{d})}:=\begin{cases}\biggl{(}\sum_{|\alpha|=k}\|\partial^{\alpha}f\|_{L^{p}(\mathbb{R}^{d})}^{p}\biggr{)}^{1/p}&p\in[1,\infty),\\ \max_{|\alpha|=k}\|\partial^{\alpha}f\|_{L^{\infty}(\mathbb{R}^{d})}&p=\infty.\end{cases}$$

$$V_{i}^{P}:=\operatorname{span}\bigl{\{}x^{\alpha}\,\bigr{|}\,|\alpha|\leq P\bigr{\}},$$

$$V_{\sigma}^{P}:=\operatorname{span}\bigl{\{}\varphi_{i}v_{i}\,\bigr{|}\,v_{i}\in V_{i}^{P}\bigr{\}}$$

$$|u-v_i{|}_{W^{k,p}(\Omega \cap {\Omega }_i)}\le C{\sigma }^{P+1-k}|u{|}_{W^{P+1,p}(\Omega \cap {\Omega }_i)}\forall k\in {\mathbb{N}}_0,k\le P+1,$$

$$|u-\mathcal{P}u{|}_{W^{k,p}(\Omega )}\le C{\sigma }^{P+1-k}|u{|}_{W^{P+1,p}(\Omega )}\forall k\in {\mathbb{N}}_0,k\le P+1,$$

$$\|\partial^{\alpha}(u-\mathcal{P}u)\|_{L^{p}(\Omega)}^{p}=\bigl{\|}\partial^{\alpha}\sum_{i}(u-v_{i})\varphi_{i}\bigr{\|}_{L^{p}(\Omega)}^{p}=\bigl{\|}\sum_{\beta\leq\alpha}\sum_{i}{\binom{\alpha}{\beta}}\partial^{\beta}\varphi_{i}\partial^{\alpha-\beta}(u-v_{i})\bigr{\|}_{L^{p}(\Omega)}^{p}.$$

$$\biggl{|}\sum_{\beta\leq\alpha}\sum_{i}\begin{pmatrix}\alpha\\ \beta\end{pmatrix}\partial^{\beta}\varphi_{i}\partial^{\alpha-\beta}\bigl{(}u-v_{i}\bigr{)}\biggr{|}^{p}\leq C(\alpha,p)\sum_{\beta\leq\alpha}\biggl{|}\sum_{i}\partial^{\beta}\varphi_{i}\partial^{\alpha-\beta}\bigl{(}u-v_{i}\bigr{)}\biggr{|}^{p}$$

$$\|\partial^{\alpha}(u-\mathcal{P}u)\|_{L^{p}(\Omega)}^{p}\leq C\sum_{\beta\leq\alpha}\ \bigl{\|}\sum_{i}\partial^{\beta}\varphi_{i}\partial^{\alpha-\beta}(u-v_{i})\bigr{\|}_{L^{p}(\Omega)}^{p}.$$

$$\biggl{|}\sum_{i}\partial^{\beta}\varphi_{i}\partial^{\alpha-\beta}(u-v_{i})\biggr{|}^{p}\leq C(M,p)\sum_{i}\bigl{|}\partial^{\beta}\varphi_{i}\partial^{\alpha-\beta}(u-v_{i})\bigr{|}^{p}$$

$$\Vert {\partial }^{\alpha }(u-\mathcal{P}u){\Vert }_{L^p(\Omega )}^p\le C\sum _i\sum _{\beta \le \alpha }\Vert {\partial }^{\beta }{\phi }_i{\partial }^{\alpha -\beta }(u-v_i){\Vert }_{L^p(\Omega \cap {\Omega }_i)}^p.$$

$$\begin{split}\zeta:\mathbb{R}\to[0,K^{-1}],\qquad x&\mapsto\begin{cases}0&\mbox{if $|x|\geq\tfrac{1}{2}$,}\\ K^{-1}\exp{\bigl{(}-\frac{1}{1-4x^{2}}\bigr{)}}&\mbox{else,}\end{cases}\\ K&\approx 0.221\,996\,908\,084\,039\,719,\end{split}$$

$$\zeta :{\mathbb{R}}^d\to [0,K^{-d}],(x_1,\dots ,x_d)\mapsto \prod _{i=1}^d\zeta (x_i),$$

$${\Omega }_i$$

$${\omega }_i$$

$$\varphi_{i}(x):=\bigl{(}\chi_{\omega_{i}}\star\zeta_{\sigma}\bigr{)}(x)=\int_{\omega_{i}}{\zeta_{\sigma}(x-y)}\,{\mathrm{d}}{y}.$$

$${\phi }_i\in C_0^{\infty }({\mathbb{R}}^d)$$

$$\sum _{i\in {\mathbb{Z}}^d}{\phi }_i(x)$$

$$|{\phi }_i{|}_{W^{k,p}({\mathbb{R}}^d)}$$

$$\sum _{i\in {\mathbb{Z}}^d}{\phi }_i(x)=\sum _{i\in {\mathbb{Z}}^d}{\int }_{{\omega }_i}{\zeta }_{\sigma }(x-y)\mathrm{d}y={\int }_{{\mathbb{R}}^d}{\zeta }_{\sigma }(x-y)\mathrm{d}y=1.$$

$$\Vert {\chi }_{{\omega }_i}\star {\partial }^{\alpha }{\zeta }_{\sigma }{\Vert }_{L^p({\mathbb{R}}^d)}\le \Vert {\chi }_{{\omega }_i}{\Vert }_{L^p({\mathbb{R}}^d)}\Vert {\partial }^{\alpha }{\zeta }_{\sigma }{\Vert }_{L^1({\mathbb{R}}^d)}={\sigma }^{d/p-k}\Vert {\partial }^{\alpha }\zeta {\Vert }_{L^1({\mathbb{R}}^d)}.$$

$$\varphi_{i}(x)\equiv\hat{\varphi}\biggl{(}\frac{x-x_{i}}{\sigma}\biggr{)},\quad\hat{\varphi}(x):=\int_{(-\frac{1}{2},\frac{1}{2})^{d}}{\zeta(x-y)}\,{\mathrm{d}}{y}.$$

$$Q_i:=(i_1\sigma ,(i_1+1)\sigma )\times \dots \times (i_d\sigma ,(i_d+1)\sigma ),i\in {\mathbb{Z}}^d.$$

$$\mathcal{B}_{j}^{P}:=\biggl{\{}\underbrace{\biggl{(}\frac{x-x_{j}}{\sigma}\biggr{)}^{\alpha}}_{=:g_{j,\alpha}(x)}\,\biggr{|}\,|\alpha|\leq P\biggr{\}},\qquad j\in\mathbb{Z}^{d}$$

$$g_{j,\alpha}(x)=\bigl{(}\hat{g}_{m,\alpha}\circ\Phi_{i}^{-1}\bigr{)}(x),\qquad x\in Q_{i},\ j\in\mathcal{J}_{i},$$

Full text

{AMS}

65N12, 65N15, 65N30, 65N75, 65N80, 65N85

1 Introduction

Vortex particle methods are numerical schemes for solving the incompressible Navier–Stokes equations. Instead of their formulation in the primitive variables velocity $\mathbf{u}$ and pressure, they make use of the equivalent formulation in terms of the vorticity $\boldsymbol{\omega}=\nabla\times\mathbf{u}$. The core idea is most easily explained in the two-dimensional, inviscid case, where the equation for the scalar vorticity $\omega$ in a bounded domain $\Omega$ then reads:

[TABLE]

Let us for the moment assume the velocity field would be known, in which case we have a linear equation. We then can discretize the initial vorticity field with particles:

[TABLE]

where $\Gamma_{i}$ and $x_{i}$ denote the circulation and position of particle $i$, and $\delta$ is the Dirac delta function. Such particle fields can be seen as quadrature rules for integrating smooth functions $\varphi$ against the vorticity $\omega$ we are aiming to approximate. Let us for example assume we are given a quasi-uniform, shape-regular triangulation of the domain $\Omega$ of mesh-width $h$. If one then applies a quadrature rule of exactness degree $m$ with positive weights to each cell, one obtains a set of quadrature nodes $x_{i}$ with associated weights $w_{i}$. A particle field approximation could then be obtained by setting $\Gamma_{i}:=w_{i}\omega(x_{i})$ in Eq. 2. For such a particle field one can prove error-bounds of the form:

[TABLE]

for $m+1>d/2$, where $d$ is the number of spatial dimensions. Here and throughout this text the symbol $C$ refers to a generic positive constant which is independent of the functions involved.

The reason for choosing this particular discretization is the availability of an analytic solution of Eq. 1 in this case. If one modifies the particles’ positions according to $\frac{{\mathrm{d}}{x_{i}}}{{\mathrm{d}}{t}}=\mathbf{u}(x_{i}(t))$, $i=1,\ldots,N$, the resulting approximation $\omega_{h}$ fulfills the vorticity equation Eq. 1 exactly, i. e., the only error in the approximation comes from the initialization error [9, Appendix A].

In practice, however, the velocity field is of course not known and needs to be retrieved from the vorticity. Let us for simplicity assume that the velocity would vanish at the boundaries. In this case the Helmholtz decomposition theorem tells us that the velocity can be retrieved through the Biot–Savart law without any boundary integral terms:

[TABLE]

where $\star$ denotes convolution. We have the following classical estimate due to Calderón and Zygmund [7]:

[TABLE]

The problem is that the particle approximation $\omega_{h}\not\in L^{2}(\Omega)$ is not smooth enough to apply this estimate; applying the Biot–Savart law to the particle field directly yields a singular velocity field. The question we try to answer in this paper is how to obtain an accurate, smooth approximation $\omega_{\sigma}\in L^{2}(\Omega)$ from the particle field $\omega_{h}$, where $\sigma$ refers to a smoothing length, which will be defined precisely later. This problem is called particle regularization. Once such a smooth approximation has been obtained, one closes the system of equations by setting $\mathbf{u}_{\sigma}:=\mathbf{K}\star\omega_{\sigma}$ and modifying the particle’s positions according to $\frac{{\mathrm{d}}{x_{i}}}{{\mathrm{d}}{t}}=\mathbf{u}_{\sigma}(x_{i}(t),t)$ instead. It is this natural treatment of convection which makes vortex methods so appealing. Given an appropriate choice of $\omega_{\sigma}$ one can show that the resulting method is essentially free of artificial viscosity and conserves mass, circulation, linear momentum, and angular momentum, and the energy of $\mathbf{u}_{\sigma}$ exactly [9, Section 2.6]. When extended to handle physical viscosity, this makes the method particularly attractive for flows at high Reynolds numbers, see for example the recent work by Yokota et al. [25].

The most common approach to the regularization problem is to mollify the particle field with a certain, radially symmetric blob-function $\zeta_{\sigma}$: $\omega_{\sigma}:=\omega_{h}\star\zeta_{\sigma}$, where $\sigma$ denotes the radius of the blob’s core [9, Section 2.3]. Many commonly used blob-functions have infinite support, effectively extending $\omega_{\sigma}$ from $\Omega$ to $\mathbb{R}^{d}$ and blurring the domain’s boundaries. There are approaches to use blob-functions with varying shapes near boundaries [24, 18] or to use image particles outside of the domain [9, Section 4.5.2]. These approaches assume that the boundaries are flat and usually fail in the presence of sharp corners or kinks. Another approach is to interpret the particles’ circulations as weighted function values $\Gamma_{i}=w_{i}\omega(x_{i})$, where the $w_{i}$ are weights from an underlying quadrature rule. While this is strictly speaking only the case during the initialization stage, the approach is then to create a triangulation of the domain using the particles’ positions as grid nodes and to use these values to construct a piece-wise linear approximation [22]. This requires a mesh to be regenerated at every time-step, which is problematic as the particle field gets distorted over time. In Vortex-in-Cell (VIC) schemes one uses interpolation formulas to obtain a grid-based approximation of the vorticity field. In the vicinity of boundaries these formulas need to be specifically adapted to the particular geometry at hand and cannot be used for arbitrary domains [10]. In general the regularization problem causes significant difficulties, and as Cottet and Koumoutsakos point out in the introduction of their book [9]: “To our knowledge there is no completely satisfactory solution for general geometries, in particular because of the need to regularize vortices near the boundary.”

In this work we try to address this problem with the help of a finite element formulation. The non-smooth $W^{-(m+1),2}$-nature of the particle field forces us to use shape functions that are globally $W^{m+1,2}$-smooth, which is not the case for the classical, piecewise linear elements. The partition of unity finite element method (PUFEM) by Melenk and Babuška is a generalization of the classical finite element method (FEM), which can be used to obtain such smooth spaces. Even though already mentioned in their introductory paper [20], there seems to have been little research in this direction. Duarte et al. [12] describe an approach which only works for certain triangulations in two dimensions.

The generation of globally smooth shape functions on general meshes in higher dimensions is a well-known, hard problem. We instead consider simple Cartesian grids, on which the construction of smooth shape functions is easier. We then apply a fictitious domain approach to deal with general geometries. This typically results in instabilities in the cut elements. Under the name ghost penalty Burman [4] presented an effective and accurate stabilization strategy for this problem, which has for example been successfully applied to several other flow problems with cut elements [6, 15, 16, 17]. We use a similar, higher-order approach inspired by Burman and Fernández [5] as well as Cattaneo et al. [8] to achieve accuracy and stability of the resulting discretization.

The rest of this article is structured as follows. In Section 2 we define and construct smooth PUFEM spaces and analyze some of their important properties. In Section 3 we introduce a stabilized variational formulation and prove its stability and convergence. The regularization problem is then treated as a perturbation to this variational formulation. Similar to the approach with blob-functions as mentioned above, the resulting error can be split into regularization and quadrature error parts which need to be carefully balanced. The analysis will show that the smoothing parameter $\sigma$ should be taken proportional to the square-root of the particle-spacing $h$ and that this choice is in a certain sense optimal. In Section 4 we perform numerical experiments, confirming our analysis. As a consequence of the quadratic relation between $h$ and $\sigma$, the computation of the velocity field only has a computational complexity of $\mathcal{O}({h^{-\frac{d}{2}}})$, enabling the use of particle numbers on desktop workstations which were previously only possible on super computers. We finish the article with concluding remarks and acknowledgements.

2 Smooth Partition of Unity Finite Element Spaces

2.1 Basic Theory

We begin this subsection by defining smooth partitions of unity, similar to Melenk’s and Babuška’s theory [20].

Definition 2.1** (Smooth Partition of Unity).**

Let $\Omega\subset\mathbb{R}^{d}$ be an open set and let $\{\Omega_{i}\}$ be an open cover of $\Omega$ satisfying a pointwise overlap condition:

[TABLE]

Let $\{\varphi_{i}\}$ be a Lipschitz partition of unity subordinate to the cover $\{\Omega_{i}\}$ satisfying

[TABLE]

where the $C(k)$ are positive constants and the symbol $|\cdot|_{W^{k,p}(\mathbb{R}^{d})}$ refers to the Sobolev semi-norms:

[TABLE]

*Then, $\{\varphi_{i}\}$ is called a smooth partition of unity subordinate to the cover $\{\Omega_{i}\}$. The sets $\Omega_{i}$ are called patches. *

Using these functions $\{\varphi_{i}\}$, we can define the spaces for the partition of unity finite element method (PUFEM).

Definition 2.2** (PUFEM Spaces).**

Let $\{\varphi_{i}\}$ be a smooth partition of unity subordinate to the open cover $\{\Omega_{i}\}$. For $P,k\in\mathbb{N}_{0}$, $p\in[1,\infty]$ we define polynomial enrichment spaces $V_{i}^{P}\subset W^{k,p}(\Omega\cap\Omega_{i}):$

[TABLE]

and the PUFEM spaces $V_{\sigma}^{P}(\Omega)\subset W^{k,p}(\Omega):$

[TABLE]

*where $\sigma:=\max_{i}\operatorname{diam}{\Omega_{i}}$ refers to the maximum patch diameter. *

{assumption}

We will assume that the shapes of the domain $\Omega$ and the patches $\{\Omega_{i}\}$ are such that we can apply the Bramble–Hilbert lemma [2, Lemma (4.3.8)]. In particular, we will assume that for all $u\in W^{P+1,p}(\Omega\cap\Omega_{i})$, $p\in[1,\infty]$, there exists a $v_{i}\in V_{i}^{P}$ such that:

[TABLE]

where the constant $C$ is independent of $\sigma$ and $u$.

We then have the following estimate, which is a straightforward generalization of the result of Melenk and Babuška [20, Theorem 2.1].

Theorem 2.3**.**

Let $V_{\sigma}^{P}(\Omega)$ be as in Definition 2.2 and let Section 2.1 be fulfilled. Then for any $u\in W^{P+1,p}(\Omega)$, $p\in[1,\infty]$, there exists $\mathcal{P}u\in V_{\sigma}^{P}(\Omega)$ such that:

[TABLE]

*where the constant $C$ is independent of $\sigma$ and $u$. *

Proof 2.4**.**

Here, we will only consider the case $p\in[1,\infty)$; the proof for the case $p=\infty$ is analogous. With $v_{i}\in V_{i}^{P}$ as in Section 2.1 we set $\mathcal{P}u:=\sum_{i}\varphi_{i}v_{i}$. We may then write for any multi-index $\alpha$ with $|\alpha|=k$:

[TABLE]

Considering the absolute value of the expanded derivative on the right, we obtain using Hölder’s inequality:

[TABLE]

and thus:

[TABLE]

Now, using the fact that for every point $x\in\Omega$ there are at most $M$ non-zero terms in the sum over $i$, we obtain by again using Hölder’s inequality:

[TABLE]

After inserting this in the previous relation we may exchange the order of summation. Using the fact that $\varphi_{i}\equiv 0$ outside $\Omega_{i}$ we obtain:

[TABLE]

*Applying Eq. 9, Eq. 13, and collecting all the terms yields the claim. *

2.2 Construction of a Smooth Partition of Unity

In this subsection we are going to construct such a smooth partition of unity using mollification. We will make use of the following two definitions.

Definition 2.5** (Friedrichs’ Mollifier).**

The function:

[TABLE]

is called Friedrichs’ mollifier in one-dimensional space. The constant $K$ was obtained numerically, such that $\|\zeta\|_{L^{1}(\mathbb{R})}=1$. For spatial dimensions greater than $d=1$ we define Friedrichs’ mollifier using the product:

[TABLE]

*where under a slight abuse of notation, we reused the symbol $\zeta$. *

It is well known that $\zeta\in C_{0}^{\infty}(\mathbb{R}^{d})$, and thus also $\zeta\in W^{k,p}(\mathbb{R}^{d}),\ k\in\mathbb{N}_{0},\ p\in[1,\infty]$. Furthermore, we have $\operatorname{supp}{\zeta}=[-\tfrac{1}{2},\tfrac{1}{2}]^{d}$, which leads to the following definition.

Definition 2.6** (Cartesian Grid).**

Given $\sigma>0$, we define Cartesian grid points $x_{i}\in\mathbb{R}^{d}$, $i\in\mathbb{Z}^{d}$, $x_{i}:=(i_{1}\sigma,\ldots,i_{d}\sigma)$. With each grid point we associate a patch $\Omega_{i}$ and a patch-core $\omega_{i}$:

[TABLE]

An illustration of these definitions is given in Fig. 1. It is obvious that the patches $\{\Omega_{i}\}$ form an open cover of $\mathbb{R}^{d}$ with $M$ from Definition 2.1 being equal to $2^{d}$. The patch-cores $\{\omega_{i}\}$ are pairwise disjoint and their closures form a (non-open) cover of $\mathbb{R}^{d}$. Using these definitions, we are now ready to construct smooth partition of unity functions $\{\varphi_{i}\},i\in\mathbb{Z}^{d}$.

Lemma 2.7**.**

For a given $\sigma>0$ and $i\in\mathbb{Z}^{d}$ let $\varphi_{i}$ be the convolution of the characteristic function $\chi_{\omega_{i}}$ of the patch-core $\omega_{i}$ with the scaled Friedrichs’ mollifier $\zeta_{\sigma}(x):=\sigma^{-d}\zeta(x/\sigma)$:

[TABLE]

One then has:

[TABLE]

Proof 2.8**.**

The first property directly follows from the classical properties of mollification [1, sections 2.28 and 2.29]. For the second property we immediately obtain:

[TABLE]

For the last property we obtain with the help of Young’s inequality for convolutions for every multi-index $\alpha$ with $|\alpha|=k$:

[TABLE]

Remark 2.9**.**

There is no closed-form expression for the functions $\{\varphi_{i}\}$ available. However, it is important to notice that we have:

[TABLE]

*Furthermore, $\hat{\varphi}$ inherits the product structure of $\zeta$. In a computer implementation it is thus sufficient to tabulate values for $\hat{\varphi}$ corresponding to the case $d=1$. We can then efficiently approximate $\hat{\varphi}$ using, e. g., cubic Hermite splines. The graph of this function can be seen in Fig. 2. *

2.3 Reference Element and Inverse Estimates

In this subsection we illustrate that the smooth partition of unity constructed in Section 2.2 leads to spaces that can be treated in a manner similar to conventional finite element spaces. In particular, we can subdivide $\mathbb{R}^{d}$ into elements:

[TABLE]

Such an element is for example depicted on the right of Fig. 1. Every $Q_{i}$ may be seen as the image of the reference element $\hat{Q}:=(0,1)^{d}$ under the transformation $\Phi_{i}:\hat{Q}\to Q_{i},\ \hat{x}\mapsto x_{i}+\sigma\hat{x}$. In every element $Q_{i}$ we have a fixed set $\mathcal{J}_{i}$ of $2^{d}$ overlapping patches $\Omega_{j}$. Introducing:

[TABLE]

as bases for the enrichment spaces $V_{j}^{P}$, one quickly sees that within each element $Q_{i}$ the basis functions $g_{j,\alpha}$ can be expressed in terms of mapped reference functions $\hat{g}_{m,\alpha}$:

[TABLE]

where $m$ is the index of the node in the reference element that corresponds to $x_{j}$. Due to Remark 2.9, the same holds true for the partition of unity functions $\varphi_{j}$. This allows us to infer the following classical result, which follows from a scaling argument and the norm-equivalence of finite-dimensional spaces [2, Lemma (4.5.3)].

Lemma 2.10** (Inverse Estimates).**

Let $V_{\sigma}^{P}(\Omega)$, $\sigma>0$, $P\in\mathbb{N}_{0}$ be as in Definition 2.2, with patches as in Definition 2.6. Then, for any element $Q_{i}\subset\Omega$ as defined in Eq. 30 that is completely contained in the domain and every $v_{\sigma}\in V_{\sigma}^{P}(\Omega)$ one has:

[TABLE]

*where the constant $C$ is independent of $\sigma$ and $i$. *

3 Stabilized Variational Formulation

In this section we will introduce a stabilized variational formulation with the aim of mimicking of the $L^{2}(\Omega)$-orthogonal projector onto $V_{\sigma}^{P}(\Omega)$. As the inverse estimates Eq. 33 are not available for elements $Q_{i}$ cut by the boundary $\partial\Omega$, we will employ a fictitious domain approach. In order to ensure coercivity of the resulting bilinear form on the entire fictitious domain, we will add a stabilization term in the cut cells. Once consistency and stability of this formulation have been established, we will model the regularization process as a perturbation to this variational problem.

3.1 Basic Definitions and Conditions

We will restrict ourselves to Hilbert spaces ($p=2$), due to the rich theoretical framework available for this case. We will assume that the shape of the domain $\Omega$ is such that we may apply the Stein extension theorem [1, Section 5.24], i. e., there exists a bounded linear extension operator $\mathcal{E}:W^{k,2}(\Omega)\to W^{k,2}(\mathbb{R}^{d})$ for any natural number $k$. We explicitly wish to include functions that do not vanish on the boundary $\partial\Omega$. For this reason, for any domain $\square\subset\mathbb{R}^{d}$, we will denote by $W^{-k,2}(\square):=W^{k,2}(\square)^{\prime}$ the dual space of $W^{k,2}(\square)$. (Opposed to the convention $W^{-k,2}(\square)=W_{0}^{k,2}(\square)^{\prime}$).

We will need certain geometrical definitions. To this end, let $\sigma>0$ be arbitrary but fixed. We define the fictitious domain $\Omega_{\sigma}$ as the union of all elements that intersect the physical domain $\Omega$. Furthermore we define cut and uncut elements $\Omega_{\sigma}^{\Gamma}$ and $\Omega_{\sigma}^{\circ}$, respectively:

[TABLE]

with $Q_{i}$ as in Eq. 30. Here, we write under a slight abuse of notation $Q_{i}\in\Omega_{\sigma}^{\Gamma}$ if $Q_{i}\subset\Omega_{\sigma}^{\Gamma}$. These domains obviously fulfill $\Omega_{\sigma}^{\circ}\subset\Omega\subset\Omega_{\sigma}$, $\Omega_{\sigma}=\operatorname{int}(\operatorname{clos}\Omega_{\sigma}^{\circ}\cup\operatorname{clos}\Omega_{\sigma}^{\Gamma})$, and $\Omega_{\sigma}^{\circ}\cap\Omega_{\sigma}^{\Gamma}=\emptyset$. Two elements $Q_{i}$ and $Q_{i}^{\prime}$ will be called neighbors if they share at least one node on the Cartesian grid. We will make the following somewhat technical assumption: for every $Q_{i}\in\Omega_{\sigma}^{\Gamma}$ there is a finite sequence of elements $(Q_{i}=Q_{i,1},Q_{i,2},\ldots,Q_{i,K})\subset\Omega_{\sigma}^{\Gamma}$ with the following properties: the number $K$ is bounded independent of $\sigma$, every pair of two subsequent elements are neighbors, and $Q_{i,K}$ has a neighbor in $\Omega_{\sigma}^{\circ}$. This condition means that one can always reach uncut elements from cut elements in a bounded number of steps. For sufficiently fine Cartesian grids this condition is often fulfilled with $K=1$; if necessary it can be enforced by moving additional elements from $\Omega_{\sigma}^{\circ}$ to $\Omega_{\sigma}^{\Gamma}$.

3.2 Introduction of a Higher-order Stabilization Term

The basic idea of the ghost penalty method is to control the norm of cut elements by relating them to neighboring uncut elements. In the aforementioned articles [6, 16, 17], for example, this is done by controlling the norms of the gradient-jumps at element boundaries. However, as our PUFEM spaces are globally smooth, they do not contain such jumps. Burman and Fernández [5] and Cattaneo et al. [8] instead use the Brezzi–Pitkäranta stabilization [3]. We will use a higher-order variant of this idea and define the following bilinear form:

[TABLE]

such that $j(u_{\sigma},u_{\sigma})=\sigma^{2(P+1)}|u_{\sigma}|_{W^{P+1,2}(\Omega_{\sigma}^{\Gamma})}^{2}$. We then obtain the following result:

Lemma 3.1**.**

Let $u_{\sigma}\in V_{\sigma}^{P}(\Omega_{\sigma})$. One then has with constants $c$ and $C$ independent of $\sigma$ and the position of $\partial\Omega$ relative to the Cartesian grid:

[TABLE]

Before moving on to the proof of this lemma, let us remark that the stabilization term is necessary. Look for example at the configuration shown in Fig. 3. The partition of unity function corresponding to the node of $Q_{i}$ opposite to $x^{*}$ vanishes on $\Omega_{\sigma}^{\circ}$. Thus its $L^{2}(\Omega_{\sigma})$-norm cannot be controlled by looking at $\Omega_{\sigma}^{\circ}$ only, unless one adds a stabilization term.

Proof 3.2**.**

The second inequality directly follows from the inverse inequalities Eq. 33. For the first inequality, let us first consider the case $K=1$, i. e., a cut element $Q_{i}\in\Omega_{\sigma}^{\Gamma}$ and an associated uncut element $Q_{i}^{\prime}\in\Omega_{\sigma}^{\circ}$ which share a Cartesian grid point $x^{*}$, as for example illustrated in Fig. 3. This configuration can be mapped to one of $P^{2^{d}}_{2}=4^{d}-2^{d}$ reference cases with reference elements $\hat{Q}$ and $\hat{Q}^{\prime}$ using the transformation $\hat{x}=\Phi^{-1}(x):=(x-x^{*})/\sigma$, such that $\hat{x}^{*}=0$. For an arbitrary function $v_{\sigma}\in V_{\sigma}^{P}(Q_{i}\cup Q_{i}^{\prime})$ one obtains with $\hat{v}:=(v_{\sigma}\circ\Phi)\in V_{1}^{P}(\hat{Q}\cup\hat{Q}^{\prime})$:

[TABLE]

We claim that the following expression constitutes a norm on $V_{1}^{P}(\hat{Q}\cup\hat{Q}^{\prime})$:

[TABLE]

It suffices to show that $\|\hat{v}\|_{*}=0\Longrightarrow\hat{v}=0.$ From $\|\hat{v}\|_{L^{2}(\hat{Q}^{\prime})}^{2}=0$ we obtain $\hat{v}\equiv 0$ on $\hat{Q}^{\prime}$ and due to the global smoothness of $\hat{v}$ also $\partial^{\alpha}\hat{v}(\hat{x}^{*})=0$ for all multi-indices $\alpha\in\mathbb{N}_{0}^{d}$. From $|\hat{v}|_{W^{P+1,2}(\hat{Q})}^{2}=0$ we obtain $\partial^{\alpha}\hat{v}\equiv 0,|\alpha|=P+1$ on $\hat{Q}$. Together with $\partial^{\alpha}\hat{v}(\hat{x}^{*})=0$ this implies $\hat{v}\equiv 0$ on $\hat{Q}$ as well. Thus $\|\cdot\|_{*}$ is indeed a norm. After employing the norm-equivalence of finite-dimensional spaces, we can transform back to $Q_{i}\cup Q_{i}^{\prime}$ and obtain:

[TABLE]

*The case $K>1$ with sequences of cells $(Q_{i}=Q_{i,1},Q_{i,2},\ldots,Q_{i,K}=Q_{i}^{\prime})$ follows by induction. Now, summing over all elements and using the finite overlap condition $M=2^{d}$, the claim follows. *

Note that the proof crucially depends on the global smoothness of the spaces $V_{\sigma}^{P}$. In particular, this stabilization does not work with the conventional finite element spaces. As an example consider the case depicted in Fig. 3, and set $v_{\sigma}:=0$ on $Q_{i}^{\prime}$, and $v_{\sigma}:=(x-x^{*})\cdot y$ on $Q_{i}$, where $y\in\mathbb{R}^{d}$ is an arbitrary non-zero vector.

3.3 Stability and Convergence

As described before, we are aiming to mimic the $L^{2}(\Omega)$-orthogonal projector. To this end, we introduce the bilinear forms $a$ and $A$:

[TABLE]

where $\varepsilon>0$ denotes a user-defined stabilization parameter. We define the variational problem as: given any $u$ for which the following makes sense, find $u_{\sigma}\in V_{\sigma}^{P}(\Omega_{\sigma})$ such that:

[TABLE]

We then obtain the following two results.

Theorem 3.3** (Stability).**

The bilinear form $A$ from Eq. 41 fulfills with a constant $C(\varepsilon)$ independent of $\sigma$, $u_{\sigma}$, and the position of $\partial\Omega$ relative to the grid:

[TABLE]

Proof 3.4**.**

For any $u_{\sigma}\in V_{\sigma}^{P}(\Omega_{\sigma})$ one has with the help of Eq. 36:

[TABLE]

Theorem 3.5** (Convergence).**

Let $u\in W^{P+1,2}(\Omega_{\sigma})$. The solution $u_{\sigma}\in V_{\sigma}^{P}(\Omega_{\sigma})$ of the variational problem Eq. 42 then satisfies the following error bound:

[TABLE]

*where the constant $C(\varepsilon)$ is independent of $\sigma$, $u$, and how the boundary $\partial\Omega$ intersects the grid. *

Proof 3.6**.**

According to Theorem 2.3, there exists $\mathcal{P}u\in V_{\sigma}^{P}(\Omega_{\sigma})$ such that:

[TABLE]

We may write:

[TABLE]

For the second term we can apply relation Eq. 46. For the first term we obtain with Theorem 3.3 and the fact that $u_{\sigma}$ solves Eq. 42:

[TABLE]

where we used the Cauchy–Schwarz inequality in the last step. Noting that by the inverse estimates Eq. 33 we have:

[TABLE]

and together with Eq. 46:

[TABLE]

After dividing both sides by $\|\mathcal{P}u-u_{\sigma}\|_{L^{2}(\Omega_{\sigma})}$ we thus obtain:

[TABLE]

*Again applying Eq. 46 to the first term yields the claim. *

3.4 Influence of the Quadrature Error

In vortex methods we are only given a particle field, i. e., a quadrature rule for integrating smooth functions against the underlying vorticity we are aiming to approximate. Furthermore the bilinear form $A$ can usually only be computed approximately, using numerical quadrature. In this subsection we are analyzing the influence of these additional sources of error.

We will assume that the bilinear form $j$ can be computed exactly. This is justified as it is sufficient to perform computations on the reference element $\hat{Q}$, which can be done up to arbitrary precision a priori. As for the bilinear form $a$, we will assume the availability of quadrature rules $I_{m}$ satisfying error bounds of the following form:

[TABLE]

Such error estimates typically arise from the application of quadrature rules of exactness degree $m$ and positive weights to the cells of a quasi-uniform triangulation of the domain $\Omega$ of mesh-width $h$. Note that due to the global smoothness of the PUFEM spaces, these quadrature rules do not need to be aligned with the Cartesian grid. We will write $a_{h}(u_{\sigma},v_{\sigma}):=I_{m}(u_{\sigma}v_{\sigma})$ for the resulting approximate bilinear form. For $u_{\sigma},v_{\sigma}\in V_{\sigma}^{P}(\Omega_{\sigma})$ one then obtains with the help of Hölder’s inequality and the inverse estimates Eq. 33:

[TABLE]

This will require us to couple $h$ and $\sigma$ through a relation like $h=\sigma^{s}$, for some $s>1$. We then obtain coercivity of $A_{h}(u_{\sigma},v_{\sigma}):=a_{h}(u_{\sigma},v_{\sigma})+j(u_{\sigma},v_{\sigma})$:

[TABLE]

where the last constant $C(\varepsilon)$ is independent of $\sigma$, $u$, and the position of $\partial\Omega$ relative to the Cartesian grid, for $\sigma>0$ small enough, $h=\sigma^{s}$, $s>1$.

For the particle field $u_{h}$ we will assume an error bound of the following form:

[TABLE]

which is the typical form arising in vortex methods [9]. Again, the particle field does not in any way need to be aligned to the Cartesian grid. Collecting all of the previous results, we are ready to prove the main result of this article.

Theorem 3.7**.**

Let $h=\sigma^{s}$, $s>1$, and denote $k:=\max\{P,m\}$. Let $u\in W^{k+1,2}(\Omega)$ and let the particle approximation $u_{h}\in W^{-(m+1),2}(\Omega)$ satisfy the error bound Eq. 55. Then for $\sigma>0$ small enough the solution $u_{\sigma}\in V_{\sigma}^{P}(\Omega_{\sigma})$ of the perturbed variational problem:

[TABLE]

satisfies the following error bound:

[TABLE]

*where the constant $C(\varepsilon)$ is independent of $\sigma$, $u$, and the position of $\partial\Omega$ relative to the Cartesian grid. *

Proof 3.8**.**

Let $v_{\sigma}\in V_{\sigma}^{P}(\Omega_{\sigma})$ denote the solution of the unperturbed variational problem Eq. 42, with $u$ extended to $\mathcal{E}u$ by the Stein extension operator. With the help of the coercivity of $A_{h}$ one then obtains:

[TABLE]

Application of the error bounds Eq. 55 and Eq. 53 as well as the inverse estimates Eq. 33 yields:

[TABLE]

and thus:

[TABLE]

Nothing that $\|v_{\sigma}\|_{L^{2}(\Omega_{\sigma})}\leq C(\varepsilon)\|u\|_{W^{P+1,2}(\Omega)}$ we obtain:

[TABLE]

Now, by the triangle inequality:

[TABLE]

*the claim follows by applying Theorem 3.5 and the boundedness of $\mathcal{E}$ to the second term. *

The part $\sigma^{P+1}\|u\|_{W^{k+1,2}(\Omega)}$ is called the smoothing error; $\sigma$ roughly corresponds to the blob-width in conventional vortex particle methods. The second part is called the quadrature error; choosing $s=1+\tfrac{P+1}{m+1}$ balances both terms. In the next section we will illustrate that the choice $P=m$ does not only “feel natural”, but also yields optimal results in a certain sense. In this case we obtain with $s=2$ an overall convergence rate of $\mathcal{O}({\sigma^{P+1}})=\mathcal{O}({h^{\frac{1}{2}(P+1)}})$.

3.5 Optimality of the Smoothed Solution

In this subsection we will assume that we can apply the bilinear form $A$ exactly, i. e., $A_{h}=A$. Furthermore we assume $P=m$. We will show that the smoothed solution $u_{\sigma}$ then satisfies the same asymptotic error bound as $u_{h}$. We will need the following corollary of Theorem 3.5.

Corollary 3.9**.**

The solution operator $\mathcal{S}$ to the problem Eq. 42 is bounded:

[TABLE]

Proof 3.10**.**

*One has $\|\mathcal{S}u\|_{W^{P+1,2}(\Omega_{\sigma})}\leq\|\mathcal{S}u-u\|_{W^{P+1,2}(\Omega_{\sigma})}+\|u\|_{W^{P+1,2}(\Omega_{\sigma})}.$ Using $\mathcal{P}$ from Theorem 2.3, we furthermore obtain: $\|\mathcal{S}u-u\|_{W^{P+1,2}(\Omega_{\sigma})}\leq\|\mathcal{S}u-\mathcal{P}u\|_{W^{P+1,2}(\Omega_{\sigma})}+\|\mathcal{P}u-u\|_{W^{P+1,2}(\Omega_{\sigma})}.$ The second term can be bounded by $C\|u\|_{W^{P+1,2}(\Omega_{\sigma})}$ due to the boundedness of $\mathcal{P}$. The first term can be bounded by first applying the inverse estimates Eq. 33 followed by estimate Eq. 51. *

Theorem 3.11** (Optimality).**

Let the conditions of Theorem 3.7 be fulfilled. Furthermore assume that $A_{h}=A$, $P=m$, and $s=2$. Then the smoothed solution $u_{\sigma}$ fulfills:

[TABLE]

Proof 3.12**.**

Let $\varphi\in W^{P+1,2}(\Omega)$ be arbitrary but fixed. With $\mathcal{P}$ from Theorem 2.3 and the Stein extension operator $\mathcal{E}$ one has:

[TABLE]

For the first term we obtain with the Cauchy–Schwarz inequality, Theorem 2.3, and Theorem 3.7:

[TABLE]

For the second term one has:

[TABLE]

It remains to show that $\|u_{\sigma}\|_{W^{P+1,2}(\Omega_{\sigma})}\leq C(\varepsilon)\|u\|_{W^{P+1,2}(\Omega)}$. To see this, note that we have:

[TABLE]

*Applying inequality Eq. 61 to the second term, collecting all the terms, and noting that $\sigma=\sqrt{h}$ yields the result. *

3.6 Conservation Properties

In the introduction we mentioned the conservation properties of vortex methods as one of their highlights. In this section we make some brief remarks on some of these properties under the assumption that $A_{h}=A$. For brevity, we will focus on the two-dimensional case, but remark that all of the results we present here analogously hold in three-dimensions.

The conserved quantities circulation, linear momentum, and angular momentum are given by $I_{0}=\int_{\Omega}{1\cdot\omega}\,{\mathrm{d}}{x}$, $\mathbf{I}_{1}=\int_{\Omega}{(x_{2},-x_{1})^{\top}\omega}\,{\mathrm{d}}{x}$, and $I_{2}=\int_{\Omega}{|x|^{2}\omega}\,{\mathrm{d}}{x}$, respectively [19, Section 1.7]. Noting that the stabilization term $j$ vanishes if one of its arguments is a polynomial of total degree less than $P$, one obtains for the solution $\omega_{\sigma}$ of Eq. 42 with right-hand side $\omega_{h}$: $(\omega_{\sigma},x^{\alpha})_{L^{2}(\Omega)}=\omega_{h}(x^{\alpha})$ for all $|\alpha|\leq P$. For $P=1$ we consequently conserve $I_{0}$ and $\mathbf{I}_{1}$, for $P=2$ one additionally conserves angular momentum $I_{2}$. This is important, because in vortex methods body forces are often computed using the relation $\mathbf{F}=-\rho\frac{{\mathrm{d}}{\mathbf{I}_{1}}}{{\mathrm{d}}{t}}$, where $\rho$ denotes the fluid’s density.

4 Numerical Experiments

We have now established the necessary results to return to our original motivation. Given a particle approximation $\omega_{h}\in W^{-(m+1),2}(\Omega)$ of the vorticity $\omega$ that satisfies an error-bound of the form $\|\omega_{h}-\omega\|_{W^{-(m+1),2}(\Omega)}\leq Ch^{m+1}\|\omega\|_{W^{m+1,2}(\Omega)}$, we want to obtain a smooth approximation $\omega_{\sigma}$, such that we can compute the corresponding induced velocity field using the Biot–Savart law $\mathbf{u}_{\sigma}=\mathbf{K}\star\omega_{\sigma}$. One can then use this approximate velocity field to advance $\omega_{h}$ in time by convecting the particles according to $\frac{{\mathrm{d}}{x_{i}}}{{\mathrm{d}}{t}}(t)=\mathbf{u}_{\sigma}(x_{i}(t),t)$.

In Section 2 we introduced the spaces $V_{\sigma}^{P}(\Omega)$ that can be used as test-spaces for the particle field $\omega_{h}$. In Section 3 we modeled the regularization problem as a perturbation to a stabilized $L^{2}$-projection onto the spaces $V_{\sigma}^{P}(\Omega)$. The analysis indicated that one should choose $P=m$ and $\sigma=\sqrt{h}$, resulting in an a-priori error estimate of $\|\omega_{\sigma}-\omega\|_{L^{2}(\Omega)}\leq C\sigma^{P+1}\|\omega\|_{W^{P+1,2}(\Omega)}$. The Calderón–Zygmund inequality Eq. 5 then tells us that one may expect $\|\mathbf{u}_{\sigma}-\mathbf{u}\|_{W^{1,2}(\Omega)}\leq C\sigma^{P+1}\|\omega\|_{W^{P+1,2}(\Omega)}$ for the resulting velocity field $\mathbf{u}_{\sigma}:=\mathbf{K}\star\omega_{\sigma}$. This analogously holds in the three-dimensional case.

In this section we perform several numerical experiments. We will first describe the experimental setup. We then perform experiments on a scalar particle field and confirm the results of our analysis. In particular, the experiments will show that the common practice of choosing $\sigma$ proportional to $h$ instead of $\sqrt{h}$ does not lead to convergent schemes. We will then illustrate the practicality of our scheme, by approximating a vector-valued vorticity field, computing its induced velocity field, and measuring the error. We finish this section with experiments on the condition number of the resulting systems and its dependence of the stabilization parameter $\varepsilon$.

4.1 Setup

We define our computational test domain as $\Omega=(-\tfrac{1}{2},\tfrac{1}{2})^{3}$. While this is one of the simplest cases for mesh-based methods, due to its sharp corners and edges it is one of the hardest for conventional vortex blob methods. In order to obtain quadrature rules which are not aligned to the Cartesian grid, the mesh generator Gmsh [13] was used to obtain a tetrahedral mesh of the domain, consisting of 24 tetrahedra with maximum edge-length $h=1$. The quadrature rules are obtained by applying the mid-point rule to this mesh and its subsequent uniform refinements from level $l=0$ down to level $l=8$, corresponding to $h=2^{-8}\approx 0.004$ and $N=402\,653\,184$ quadrature nodes.

Preliminary experiments showed good results for a stabilization parameter of $\varepsilon=0.001$. Unless explicitly stated otherwise, we will use this value for all of our computations. We will use degree $P=1$ for the PUFEM spaces, set $\sigma:=Ch^{1/s}$, and experiment on various choices of $C$ and $s=1,2$. For the integration of the bilinear form $A_{h}$ we use the following approach: if in a pair of basis functions one of them has cut support, we use the same quadrature rule as for the particle field. Otherwise precomputed values from the reference element $\hat{Q}$ are used. The resulting systems of equations are solved using the conjugate gradient method, where we apply a simple diagonal scaling as preconditioner. The iteration was stopped when a relative residual of $10^{-12}$ was reached. This was usually the case after less than 100 iterations, with some exceptions for coarse refinement levels $l$ and the case $C=0.5,s=1$.

4.2 Scalar Particle Field

The common practice to choose the smoothing length $\sigma$ proportional to $h$ may in special cases be justified with the analysis of Cottet and Koumoutsakos [9, Section 2.6]. They assume that the quadrature rules used are of infinite order, essentially corresponding to the case $m=\infty$. Such rules, however, typically only exist in very special cases, such as a cube with periodic or zero boundary conditions. To show that this approach does not work in a more general setting, we aim to approximate the following function:

[TABLE]

This function does not vanish at the boundary. The application of conventional blob-methods would thus blur the boundaries and lead to only slowly converging schemes. We define the particle field as $u_{h}:=\sum_{i=1}^{N}w_{i}u(x_{i})\delta(x-x_{i})$, with $\delta$ denoting the Dirac Delta, and $x_{i}$ and $w_{i}$ being the positions and weights of the mid-point quadrature rule applied to the tetrahedra of the mesh at various refinement levels $l$.

Figure 5 shows the error $\|u-u_{\sigma}\|_{L^{2}(\Omega)}$ for $\sigma=Ch$ for various choices of $C$ at different refinement levels. Choosing $C=0.5$ results in approximations with large errors, which do not decrease significantly under mesh refinement. The case $l=8$ was not computed due to the large memory requirements. The other curves exhibit similar behavior: in the beginning and intermediate stages the error decreases, however, only at an approximately linear, not quadratic rate. This rate further decreases and approaches zero under mesh refinement, confirming the predicted bound of the quadrature error $\mathcal{O}({h^{m+1}\sigma^{-(m+1)}})=\mathcal{O}({1})$. Choosing larger values $C$ somewhat delays but does not prevent this effect, at the cost of larger errors on coarse refinement levels.

Figure 5 shows the corresponding error for the case $\sigma=C\sqrt{h}$. All choices of $C$ lead to convergent schemes which approach the predicted convergence rate of $\mathcal{O}({h})$. In our experiments, smaller choices of $C$ lead to smaller errors; however choosing $C$ too small causes larger errors in the coarser cases. In our test case a choice somewhere between $C=0.25$ and $C=0.5$ seems to be optimal.

4.3 Vector-valued Particle Field and Velocity Evaluation

In this section we show that our scheme can drastically reduce the cost of the computationally most expensive part of vortex methods, the velocity evaluation. To this end, we prescribe:

[TABLE]

This velocity field is smooth and fulfills $\nabla\cdot\mathbf{u}\equiv 0$. It was chosen such that it vanishes at the boundaries, so that it can be retrieved from the vorticity field $\boldsymbol{\omega}:=\nabla\times\mathbf{u}$ through the Biot–Savart law without any boundary integral terms:

[TABLE]

Analogous to the previous section, we define the particle approximation:

[TABLE]

Experiments in the previous section suggested a choice of $\sigma=C\sqrt{h}$, with $C$ between $0.25$ and $0.5$. We consequently choose $C=0.375$ and obtain after applying the method to each component a smoothed approximation $\boldsymbol{\omega}_{\sigma}$ with an anticipated convergence rate of $\mathcal{O}({\sigma^{2}})$ in the $L^{2}$-norm. In order to evaluate the Biot–Savart law for this vorticity field, we chose the coarsest level $l$ such that the corresponding mesh width $2^{-l}$ is smaller than $\sigma$. We then compute the orthogonal projection of $\boldsymbol{\omega}_{\sigma}$ onto the standard finite element space of piecewise linear functions on that level. The Biot–Savart integral can then be computed by summing over the tetrahedra, for which Suh published analytic formulas [23]. We couple these formulas with a fast multipole method [14, 11] for the far-field evaluation. The resulting velocity field is approximated by taking the nodal interpolation onto the standard finite element space of piecewise quadratics to obtain an approximate velocity field $\mathbf{u}_{\sigma}$.

Most conventional schemes apply the fast multipole method directly to the particle field, leading to a complexity of $\mathcal{O}({N})=\mathcal{O}({h^{-d}})$, with a large hidden constant. Note that in our case the method is applied to the coarser smoothed approximation, leading to a complexity of only $\mathcal{O}({h^{-\frac{d}{2}}})$.

Fig. 7 shows the $L^{2}$-errors in the approximate smoothed vorticity field $\boldsymbol{\omega}_{\sigma}$ and the velocity field $\mathbf{u}_{\sigma}$. The smoothed vorticity field converges at a rate of $\mathcal{O}({\sigma^{2}})=\mathcal{O}({h})$ as expected. With the Calderón–Zygmund inequality Eq. 5 we obtain that the same error bound holds for the velocity in the $W^{1,2}$-norm. As the results indicate, in the $L^{2}$-norm the error seems to reduce by one power in $\sigma$ faster, resulting in a rate of $\mathcal{O}({h^{1.5}})$.

4.4 System Condition Number

In this section we investigate the effect of the stabilization parameter $\varepsilon$ on the condition number of the system matrix. In Section 4.2 we observed instabilities on the coarse levels $l$ in the case $C=0.25,s=2$. We therefore chose this particular configuration for our experiments. We used the following set of functions as our basis for the spaces $V_{\sigma}^{P}(\Omega_{\sigma})$:

[TABLE]

We may assign a numbering $\mathcal{I}=\{1,\ldots,n\}$ to this set, and subsequently refer to its members as $\mathcal{B}_{\sigma}^{P}\ni\psi_{k}$, $k\in\mathcal{I}$. We can then define the system matrix $\mathsf{A_{h}}\in\mathbb{R}^{n\times n}$ via the relation:

[TABLE]

where $\mathsf{e_{k}}\in\mathbb{R}^{n}$ refers to the $k$-th Cartesian basis vector, and the approximate bilinear form $a_{h}$ is defined as described in the numerical setup (Section 4.1). We are then interested in the condition number of the diagonally scaled matrix $\mathsf{D^{-1}A_{h}}$, where $\mathsf{D}:=\operatorname{diag}\mathsf{A_{h}}$.

Figure 7 shows the condition number of $\mathsf{D^{-1}A_{h}}$ for various refinements levels $l$ as a function of $\varepsilon$. In the case $l=2$ the quadrature error is so large that the resulting matrix $\mathsf{A_{h}}$ ceased being positive definite for $\varepsilon=10^{-3}$, and is even singular for $\varepsilon=0$. This explains the large error observed in Section 4.2 for this case. But even then a sufficiently large choice of $\varepsilon$ results in a well conditioned system. For the finer refinement levels a choice of $\varepsilon$ between $10^{-3}$ and $10^{-1}$ seems to be optimal and reduces the matrix’ condition number below 100. The effect becomes slightly less pronounced with increasing $l$. We can thus conclude that for such a choice of $\varepsilon$ the stabilization removes the ill-conditioning of the system, especially in the presence of moderate quadrature errors.

5 Conclusions and Outlook

We have presented a new method to tackle the particle regularization problem, based on a stabilized fictitious domain formulation with smooth shape-functions. Our approach enjoys all the benefits of the conventional blob-methods: the resulting smoothed approximations are $C^{\infty}$ functions and conserve all moments up to order $P$. On top of that, our approach can accurately handle general geometries. The evaluation of the smoothed approximations is cheap and straightforward and does not require a summation over all particles as in the case of blob functions. The fact that we can only achieve a convergence rate of $\mathcal{O}(h^{\frac{m+1}{2}})$ as opposed to $\mathcal{O}(h^{m+1})$ might seem disappointing, but is intrinsic to the smoothing problem at hand. This can be illustrated in a simple one-dimensional example: given an interval of length $h$ on the real line, the $m$-node Gaussian quadrature rule will have an error bound of $\mathcal{O}(h^{2m})$. With $m$ function values, however, we can only construct an interpolation polynomial of degree $m-1$, having the halved error bound $\mathcal{O}(h^{m})$. Theorem 3.11 shows that the smoothed approximation is essentially just as accurate as the particle field. This also means that it can be used to reinitialize overly distorted particle fields. Furthermore, this means that the smoothed vorticity field has much greater length-scales than the particle spacing.

As a consequence the velocity evaluation – usually the most expensive part of vortex methods – can be drastically sped up. In our numerical experiments we gave an example of a simple mesh-based scheme for this, which was chosen because of its simplicity. A disadvantage of this approach is that the resulting velocity approximation ceases being divergence-free. For future research it would be interesting to make use of the fact that in three-dimensional space the Biot–Savart law is of the form $\mathbf{u}=\nabla\times(G\star\boldsymbol{\omega})$, where $G(x)=(4\pi|x|)^{-1}$ denotes the fundamental solution of the Laplacian. The “curl spaces” \nabla\times\bigl{(}V_{\sigma}^{P}(\Omega)\bigr{)}^{3} would thus be a more natural choice for approximating the velocity, while also being divergence-free in the strong, pointwise sense.

It is not clear whether the exact variational formulation Eq. 42 is actually unstable without stabilization. A result by Reusken [21, Theorem 5] indicates that a rescaling might be sufficient to achieve a stable formulation. On the other hand, this result assumes that the bilinear form $a$ can be computed exactly. The experiments of Section 4.4 suggest that stabilization is especially beneficial in the presence of quadrature errors.

Our current approach uses a uniform grid size $\sigma$ and a fixed polynomial degree $P$. The partition of unity approach, however, is general enough to be extended to adaptive grids and varying polynomial degrees. For the future, experiments with $\sigma$-, $P$-, or $\sigma P$-adaptive schemes are another interesting field for further research.

We believe the stabilized fictitious domain approach with smooth shape functions is not only useful for vortex particle regularization but also for other problems which require higher degrees of smoothness, such as certain problems from linear elasticity like the Kirchhoff–Love thin-plate theory.

The source code of the software used to obtain the results of this article can be obtained from the authors upon request.

Acknowledgements

We would like to thank the reviewers for their fruitful comments and suggestions. The first author would also like to express his gratitude to Sven Groß, Arnold Reusken, and all the members of the DROPS team at the Lehrstuhl für Numerische Mathematik (LNM) at RWTH Aachen University, Germany, with whom he previously worked on two-phase flows and the ghost penalty stabilization. Without their support and the knowledge received during that time, this work would not have been possible.

Last but not least, the first author receives the MEXT scholarship of the Japanese Ministry of Education and was supported by the Keio Leading Edge Laboratory of Science and Technology (no grant numbers allotted). Without their financial support this research would have been impossible to conduct.