text:
\section{Introduction}\label{Introduction}

Free surface simulations can provide highly valuable quantitative and qualitative insights in a large variety of engineering applications. The first numerical models were developed in the 1970's and the difficulties of this subject are reflected in the continued high level of activity within the research field. Several different approaches are used to model the free surface \cite{Worner.2012}: moving meshes, marker points, level set, phase fields and volume of fluid. The first two methods are not optimal for large interface deformation and are difficult to implement. The volume of fluid and level set approaches do not suffer from these problems and are the most commonly used interface methods for practical engineering applications. A prominent weakness of the level set is its lack of strict mass conservation which is caused by the fundamental problem that the signed distance function is not a passive tracer field. However, in practice it allows for good predictions if the discretization is fine enough so that the mass loss can be neglected. The volume of fluid method is strictly mass conservative but accurate models are harder to implement. The volume of fluid method based on the approach of Hirt \cite{Hirt.1981} can be divided into two categories: One with a finite interface thickness and the other with infinitesimal interface thickness. Models using finite interface thickness are easier to implement and are still used today e.g. in OpenFOAM and Fluent but are often not sufficiently accurate. Models with an infinitesimal thick interface comprise of two steps: Interface reconstruction and interface advection. The most common reconstruction method is the piecewise linear interface reconstruction method (PLIC), where the interface within a cell is represented by a plane cutting the cell into two subcells. The resulting surface lacks $C^{0}$ continuity. Several approaches for the computation of the normal vector of the plane can be found in literature. In the widely used method by Youngs\cite{Youngs.1982}, the normal vector is calculated from the gradient of the volume fraction. This method is easy to implement and works in three dimensions and on unstructured meshes but lacks accuracy. In the efficient least squares volume-of-fluid interface reconstruction algorithm (ELVIRA) \cite{Pilliod.2004} an interface plane is projected into the neighbouring cells. The plane is used to calculate a volume fraction. The error of the calculated volume fraction in the neighbouring cells is minimized by changing the orientation of the plane. This delivers accurate results on all mesh types in three dimensions. However, the method needs a two dimensional minimizer in three dimensions which complicates the implementation and makes the method numerically demanding. Currently, the most accurate approach for the estimation of the normal vector and the curvature is the height function method in which the volume fractions in columns of cells are integrated by direction (x,y,z). The derivative is calculated by finite difference operators. Height functions are second order accurate in curvature and interface normal calculation but only work with hexahedral structured meshes in 3D. In 2D the method has been extended to unstructured meshes\cite{Ito.2014} and in 3D a so-called Embedded Height Function (EHF) approach was invented\cite{Ivey.2015}, where the VOF field is geometrically mapped onto an overlapping hexahedral mesh on which the height function and its derivatives can then be calculated. This method gives second order convergence, albeit at a fairly high price in terms of complexity of the implementation. A simpler approach was investigated by Cummins et al. \cite{Cummins.2005} who reconstructed a signed distance function, denoted RDF, after a first reconstruction step based on the Youngs method \cite{Youngs.1982}. This method works on unstructured meshes in three dimensions, and although it reduces the errors of the Youngs method by roughly one order of magnitude, it does not converge with mesh refinement.

In the present paper we present a further elaboration of Cummins' method where the main novelty lies in the introduction of several iterations over the RDF construction and the interface reconstruction. As demonstrated in the following, this leads to second order convergence with mesh refinement while still retaining the computational expenses at a moderate level.

The paper is organised as follows: In Section~\ref{sec:VOFEqn}, we clarify the role and importance of accurate interface reconstruction in the advection step of interfacial flow simulations. In Section~\ref{sec:recon}, we present the three variants of interface reconstruction procedures under investigation. In Section~\ref{Sec:Results}, we demonstrate the accuracy and convergence properties of these reconstruction schemes and compare with results from literature. In Section~\ref{sec:advectionResults}, we investigate their convergence properties in a number of advection test cases. Finally, in Section~\ref{Sec:Conclusion} we summarize our findings.

\section{The fundamental VOF equation}\label{sec:VOFEqn}

We consider a system of two immiscible, incompressible fluids occupying the space of a computational domain. The state of such a system is described by a velocity field, $\mathbf u(\mathbf x,t)$, a pressure field, $p(\mathbf x,t)$, and finally the instantaneous position of the dynamically evolving sharp fluid interface. The time evolution of the velocity and pressure field are determined by the Navier--Stokes equations together with the incompressibility condition and appropriate boundary conditions on the domain boundaries and on the fluid interface. The fluid interface may be represented by the fluid density field, $\rho(\mathbf x,t)$, taking a constant value in each of the two fluids, and thus jumping from one to the other value at the fluid interface. With this representation, the interface evolution in time is determined by the continuity equation. This is the basis of the Volume--of--Fluid (VOF) method, where, in the spirit of the Finite Volume Method (FVM), the equation is integrated over a control volume in order to account for the change of mass within the volume. With the control volume being the $i$'th cell of a FVM mesh, the volume integrated continuity equation can be written,
\begin{equation}\label{eq:volIntContEqn}
	\frac{d}{dt}\int_{V_i}\rho(\mathbf x,t)dV + \sum_f \int_f\rho(\mathbf x,t)\mathbf u(\mathbf x,t) \cdot d\mathbf S = 0,
\end{equation}
where $V_i$ is the cell volume and the sum is over the polygonal faces comprising the cell boundary. Let us call the reference fluid A and the other fluid B and denote their densities $\rho_A$ and $\rho_B$, respectively. Then we can write the density field in terms of the 3--dimensional Heaviside function, which is $1$ in fluid A and $0$ in fluid B,
\begin{equation}\label{eq:rhoOfH}
	\rho(\mathbf x,t) = \rho_A H(\mathbf x,t) + \rho_B (1 - H(\mathbf x,t)).
\end{equation}
Inserting Eqn.~\eqref{eq:rhoOfH} into Eqn.~\eqref{eq:volIntContEqn} and formally integrating in time, we get after some rearrangement
\begin{equation}\label{eq:alphaEvol}
	\alpha_i(t+\Delta t) = \alpha_i(t) - \frac1{V_i}\sum_f \int_t^{t+\Delta t}\int_f H(\mathbf x,\tau)\mathbf u(\mathbf x,\tau)\cdot d\mathbf S d\tau,
\end{equation}
where we have defined the volume fraction (of fluid A) in cell $i$ as
\begin{equation}\label{eq:alphaDef}
	\alpha_i(t) = \frac{1}{V_i}\int_{V_i}H(\mathbf x,t) dV.
\end{equation}
In the VOF method the cell volume fractions, $\alpha_i$, are used to implicitly represent the fluid interface, and so Eqn.~\eqref{eq:alphaEvol} is the exact evolution equation which is to be discretized and solved approximately as part of our interfacial flow solution algorithm.

If cells are sufficiently small then the interface within a cell can be approximated with a plane cutting the cell into two disjoint subcells filled with fluid A and fluid B, respectively. Clearly the change in $\alpha_i$ during a time interval is closely linked to the motion of the interface inside the cell during the time step. Thus, accurate reconstruction of this interface from the cell VOF data at the beginning of the time step is an essential task in geometric VOF-based interface advection methods. In the following, we describe the interface reconstruction procedures to be investigated in this paper.

\section{Interface reconstruction schemes}\label{sec:recon}

If the interface inside a cell is regarded as approximately planar, the role of the reconstruction step in a geometric VOF algorithm is to provide the interface centre position, $\mathbf x_S$, and unit normal, $\mathbf n_S$, inside each interface cell at the beginning of a time step. Here we will define an interface cell as a cell with $\epsilon < \alpha_i < 1 - \epsilon$, where we typically set $\epsilon = 10^{-8}$. By convention, $\mathbf n_S$ points away from fluid A and into fluid B.

The three procedures for obtaining $\mathbf x_S$ and $\mathbf n_S$ in an interface cell from the available VOF data are described in the following three subsections.

\subsection{\texttt{iso-Alpha}}\label{iso-Alpha}
The first procedure is the $\alpha$-isosurface based method described in \cite{Roenby.2016}, where it forms the reconstruction step in the isoAdvector VOF method. This reconstruction procedure can be described by the following three steps:
\begin{enumerate}
\item\label{ialp:VOFToVertices}If cell $i$ is an interface cell, then for each of its vertices we interpolate the surrounding cell centre $\alpha$-values to the vertex by using the inverse cell--centre--to--vertex distances as interpolation weights.
\item\label{ialp:isovalue}Based on these vertex $\alpha$ values, we calculate an internal $\alpha_0$-isosurface cutting the cell into two subcells. A root finding algorithm is used to find the (cell specific) isovalue, $\alpha_0$, such that $V_A/V_i = \alpha_i$, where $V_A$ is the volume of the subcell filled with fluid A.
\item\label{ialp:xSandnS}Finally, we calculate the face centre, $\mathbf x_S$, and face unit normal, $\hat{\mathbf n}_S$, of the internal isosurface using the same formulae as for ordinary mesh faces.
\end{enumerate}

\subsection{\texttt{iso-RDF-N}}\label{iso-RDF-N}
This method uses the $\mathbf x_S$ and $\hat{\mathbf n}_S$ obtained with \verb|iso-Alpha| as the initial guess for an iterative improvement procedure with $N-1$ subsequent iterations based on RDF-isosurfaces. The method can be described by the following steps:
\begin{enumerate}
\item\label{isoRDF:initxSandnS}Use \verb|iso-Alpha| to obtain an initial guess of $\mathbf x_S$ and $\hat{\mathbf n}_S$ in cell $i$.
\item\label{isoRDF:cellRDF}Using the plane $(\mathbf x - \mathbf x_S)\cdot \hat{\mathbf n}_S = 0$ as the interface approximation, calculate the signed distance function, RDF, in the cell centre of cell $i$ and all its point-neighbours (i.e. cells sharing a vertex with cell $i$).
\item\label{isoRDF:vertexRDF}Interpolate the RDF from the cell centres of cell $i$ and its point neighbours to all the vertices of cell $i$ using the least square interpolation method.
\item\label{isoRDF:RDFisoface}Based on these vertex RDF values, we calculate an internal RDF-isosurface cutting the cell into two subcells. A root finding algorithm is used to find the (cell specific) RDF isovalue such that $V_A/V_i = \alpha_i$, where $V_A$ is the volume of the subcell filled with fluid A.
\item\label{isoRDF:newxSandnS}Calculate the face centre and unit normal of the newly found RDF isosurface and overwrite the old $\mathbf x_S$ and $\hat{\mathbf n}_S$ vectors with these new estimates.
\item\label{isoRDF:RDFiter}Repeat steps \ref{isoRDF:cellRDF}-\ref{isoRDF:newxSandnS} $N-2$ times to improve the calculated $\mathbf x_S$ and $\hat{\mathbf n}_S$.
\end{enumerate}
$N$ is the numbers of times the root finding based interface reconstruction procedure is performed, which is by far the most expensive step in the reconstruction procedure.

\subsection{\texttt{plic-RDF-N}}\label{plic-RDF-N}
Our third and final variant is similar to \verb|iso-RDF-N|, but is based on $\nabla \alpha$ and the PLIC approach for the initial $\mathbf x_S$ and $\hat{\mathbf n}_S$ estimates. Furthermore it uses the gradient of the RDF rather than its isosurface in the iterative improvements of $\hat{\mathbf n}_S$. The method can be described by the following steps:
\begin{enumerate}
\item\label{plicRDF:initnS}If cell $i$ is an interface cell, calculate the cell-centred $\nabla\alpha$ vector and use $\hat{\mathbf n}_S = -\nabla \alpha/|\nabla \alpha|$ as the initial guess of the internal interface unit normal. (In the advection algorithm this initial $\hat{\mathbf n}_S$ estimate will be replaced by an interpolation of $\hat{\mathbf n}_S$ from the previous time step.)
\item\label{plicRDF:initxS}Use a root finding algorithm to find an interface centre, $\mathbf x_S$, inside the cell such that the plane defined by $(\mathbf x - \mathbf x_S)\cdot \hat{\mathbf n}_S = 0$ cuts the cell into two subcells satisfying $V_A/V_i = \alpha_i$, where $V_A$ is the volume of the subcell filled with fluid A.
\item\label{plicRDF:cellRDF}Using the plane $(\mathbf x - \mathbf x_S)\cdot \hat{\mathbf n}_S = 0$ as the interface approximation, reconstruct the signed distance function, RDF, in the centres of cell $i$ and all its point-neighbours.
\item\label{plicRDF:gradRDF}Based on the RDF value in cell $i$ and in all its point neighbours calculate the gradient of the RDF using the least square method.
\item\label{plicRDF:newnS}Assuming the RDF gradient to be constant within cell $i$, normalise it and use the result as our new estiamate of $\hat{\mathbf n}_S$.
\item\label{plicRDF:newxS}Apply step \ref{plicRDF:initxS} to find a new interface centre, $\mathbf x_S$, which is consistent with the new $\hat{\mathbf n}_S$ from step \ref{plicRDF:newnS}.
\item\label{plicRDF:RDFiter}Repeat steps \ref{plicRDF:cellRDF}-\ref{plicRDF:newxS} $N-2$ times to improve the calculated $\mathbf x_S$ and $\hat{\mathbf n}_S$.
\end{enumerate}
Also here $N$ is the number of times the expensive root finding based interface reconstruction procedure is performed. The default value is $N=5$ for the reconstruction test cases and $N=2$ for the advection test cases.

In the following sections we give some further details regarding the implementation of the three proposed reconstruction schemes.

\subsection{Isosurface reconstruction from vertex data}

In all three methods described above there are steps involving reconstruction of an isosurface inside a cell based on values of a function in the cell vertices. In \verb|iso-Alpha| (step \ref{ialp:isovalue}) and \verb|iso-RDF-N| (step \ref{isoRDF:initxSandnS}), this is done with the volume fractions interpolated to the vertices. In \verb|iso-RDF-N| (step \ref{isoRDF:RDFisoface}) and \verb|plic-RDF-N| (step \ref{plicRDF:initxS}), we find planar RDF isosurfaces based on RDF data at the vertices.

Here we describe the procedure to obtain an internal $f_0$-isosurface inside a cell where the values of the function $f$ are known at the cell vertices: Consider the edge between the two vertices, $x_i$ and $x_j$, with function values $f_i$ and $f_j$, respectively. If $\min(f_i,f_j) \leq f_0 < \max(f_i,f_j)$, then the edge is cut by the $f_0$-isosurface at the point
\begin{equation}
		x_{\textrm{cut}} = x_i + \lambda_{ij}(x_j - x_i),\quad \textrm{   where   } \quad \lambda_{ij} = \frac{f_0 - f_i}{f_j - f_i}.
\end{equation}
After finding all such edge cut points for the cell, the cut points can be connected to form the periphery of an internal ``isoface'' as illustrated in Fig.~\ref{fig:slicedCellwithFace}(left). We remark that the face centre, $\mathbf x_S$, and unit normal, $\hat{\mathbf n}_S$, for such an isoface can be calculated the same way as for a mesh face. With $N$ vertex points defining the face, its area vector is defined as the area weighted average of its traingular decomposition based on an average point,
\begin{equation}
	\mathbf n_S = \sum_{k=1}^N \mathbf n_{S,k}, \ \textrm{where} \ \mathbf n_{S,k} = \frac12 (\mathbf x_{k+1} - \mathbf x_k)\times(\bar{\mathbf x} - \mathbf x_k) \ \textrm{and} \
	\bar{\mathbf x} = \frac1N \sum_{k = 1}^N \mathbf x_k.
\end{equation}
Here we use $\mathbf x_{N+1} = \mathbf x_{1}$, and the vertex points are assumed to be ordered such that each $\mathbf n_{S,k}$ points into fluid B. The unit normal vector and face centre point are then calculated as
\begin{equation}
	\hat{\mathbf n}_S = \frac{\mathbf n_S}{|\mathbf n_S|} \quad \textrm{and} \quad
	\mathbf x_S = \sum_{k = 1}^N \frac{|\mathbf n_{S,k}|}{|\mathbf n_S|}\frac{\mathbf x_k + \mathbf x_{k+1} + \bar{\mathbf x}}{3}.
\end{equation}
These definitions are directly adapted from the OpenFOAM code library.
\begin{figure}[tbph]
	\centering
	\includegraphics[width=1.0\textwidth]{slicedCellwithFace}
	\caption{Isosurface inside a hexahedral cell. Left: ``isoface'' constructed from connecting its cut points along the cell edges. Right: Construction of a subface as part of the boundary of the subcell submerged in fluid A.}
	\label{fig:slicedCellwithFace}
\end{figure}

\subsection{The root finding based subcell reconstruction}\label{ssec:rootFinding}
  In all the steps where an isosurface is constructed inside a cell $i$, this must be done multiple times, to find the isovalue yielding a volume $V_A$ of the subcell submerged in fluid A such that $V_A/V_i = \alpha_i$. This subcell's boundary is defined by 1) all the cell faces that are fully submerged in fluid A, 2) the ``submerged'' part of all the faces that are cut by the isosurface (see Fig.~\ref{fig:slicedCellwithFace}(right)), and 3) the internal isosurface itself. The subcell volume, $V_A$, can be calculated based on the face centres, $\mathbf x_f$, and face normals, $\mathbf n_f$, the same way as it is done for a polygonal mesh cell (pyramid decomposition)
\begin{equation}
	V_A = \sum_f \frac13|\mathbf n_f\cdot (\mathbf x_f - \bar{\mathbf x})|, \quad \textrm{ where } \quad \bar{\mathbf x} = \frac1N\sum_f \mathbf x_f,
\end{equation}
and the sums run over all the faces of the subcell. Reconstructing the internal interface and the subfaces can be relatively costly, if done many times to reach $V_A/V_i = \alpha_i$ within a specified tolerance (typically $10^{-8}$). As described in \cite{Roenby.2016} the number of necessary evaluations can be limited by exploiting that $V_A$ is a piecewise cubic polynomial in $f_0$ with polynomial coefficients being constant on the $f_0$ intervals between the vertex values. Thus, once we have done the expensive $V_A$ calculation a few times to find the correct $f_0$ subinterval, we only need two further evaluations (four in all on the subinterval) to obtain the polynomial coefficients of the function $V_A(f_0)$. Then we can use an inexpensive polynomial root finding algorithm to quickly get within the required tolerance.

\subsection{Reconstructing the signed distance function}

A good VOF based interface advection algorithm ensures that the jump from $\alpha = 0$ to $\alpha = 1$ happens over one cell everywhere along the interface. This sharpness of the volume fraction field comprises a significant challenge when derivatives are required e.g. for curvature calculations. As is well known from the Level Set literature, the signed distance function is much better suited for this kind of calculations. This is the motivation for attempting to numerically reconstruct the signed distance function from the VOF data.

\begin{figure}[!tbp]
	\centering
	\begin{minipage}{0.45\linewidth}
		\includegraphics[width=\textwidth]{alphaAndHeight.eps}
		\caption{Representation of a circle with the volume fraction field and the distance function field.}
		\label{fig:CompareAlphaHeight}
	\end{minipage}
	\hfill
	\begin{minipage}{0.45\linewidth}
		\includegraphics[width=\textwidth]{signedDistConceptTri}
		\caption{Reconstruction of the signed distance function.}
		\label{fig:SignedDistConstruct}
	\end{minipage}
\end{figure}

Fig.~\ref{fig:CompareAlphaHeight} illustrates the difference in smoothness of the raw VOF data and an RDF for a circle on a triangular prism mesh. An RDF is similar to a signed distance function, but where the length of the gradient of the signed distance function is equal to 1, this is not necessary exactly fulfilled for an RDF.

When constructing the RDF in a cell we calculate the shortest distance from the cell centre to the plane. However, a cell may in general have several neighbour cells holding an interface as illustrated in Fig.~\ref{fig:SignedDistConstruct}. The question is then which interface to use to calculate the RDF in the cell centre? In the current implementation the RDF in a cell is obtained as a weighted average of the distances to the interfaces in the cell itself and its neighbours. So for cell $j$, the distance to the interface in cell $i$ is calculated as
\begin{equation}
	\tilde F_{ij} =  \hat{\mathbf n}_S \cdot (\mathbf x_j -  \mathbf x_S),
\end{equation}
where $\mathbf x_S$ and $\hat{\mathbf n}_S$ denote interface centre and unit normal in cell $i$. From these distances the RDF in cell $j$ is calculated as
\begin{equation}
	F_j =  \frac{\sum w_{ij} \tilde F_{ij}}{\sum w_{ij}},
\end{equation}
where the sum is over all point neighbours of cell $j$ that are interface cells and the weighting factor is chosen to be
\begin{equation}
	w_{ij} =  \frac{|\hat{\mathbf n}_S \cdot (\mathbf x_j -  \mathbf x_S)|^A}{|\mathbf x_j -  \mathbf x_S|^A},\quad \textrm{with} \quad A = 2.
\end{equation}
Cummins et al. \cite{Cummins.2005} found that their RDF approach worked best with a higher exponent, $A=25$, and their choice is frequently adapted in literature. In our proposed iterative algorithm, we have found the exponent to have no influence on the accuracy for well--resolved interfaces. For under--resolved regions we have found that a higher exponent value can cause the method to diverge. Since $A = 2$ is also less expensive to compute than a higher exponent, this seems like a reasonable choice.

RDF values at the centres of faces belonging to the domain boundary are computed in a similar fashion as described above.

\subsection{Initial estimate of $\hat{\mathbf n}_S$ in \texttt{plic-RDF-N}}\label{ssec:initnS}
In step \ref{plicRDF:initnS} of \verb|plic-RDF-N| (Section \ref{plic-RDF-N}) the most obvious first guess of the interface orientation inside interface cells is in terms of a numerically calculated gradient of the volume fraction data. With this approach the new method is similar to the approach by Cummins et al. \cite{Cummins.2005}. However, utilizing the gradient from the volume fraction data can provide an erroneous initial guess, especially on fine meshes. Hence, this approach will only be applied in the first time step of an interface advection simulation. For subsequent time steps we have found that an orientation estimate based on the previous time step is more accurate and thus requires fewer iterations. In the current implementation we use nearest--neighbour interpolation for this first guess of $\hat{\mathbf n}_S$. If the normal vectors in all surrounding cells differ by more than 30 degrees, the gradient of volume fraction is used as initial guess. On coarse meshes the gradient of the volume fraction does not limit the accuracy.

\subsection{Least square fitting}
In step \ref{plicRDF:gradRDF} of the \verb|plic-RDF-N| algorithm described in Section~\ref{plic-RDF-N} we use the numerically calculated gradient of the RDF for the interface orientation. To ensure that the gradient scheme does not limit the accuracy of our procedure, a least square method is implemented. Sozer et al. \cite{Sozer.2014} demonstrated that cell based least square gradient schemes are superior to the standard Gauss-Green method. OpenFOAM already has a least square gradient scheme implemented, but it is applied to the whole domain making it inefficient for our purpose. The method implemented here only computes the gradient in the cells where it is required. For the gradient calculation, we assume the RDF to be locally trilinear:
\begin{equation}\label{eq:fitFunc}
	F(x,y,z) =  a x + b y + c z + d.
\end{equation}
With $N$ cells in the stencil (comprising of the interface cell and its point neighbours), application of the least squares method to Eqn.~\ref{eq:fitFunc} results in the linear system:
 \begin{equation}\label{eq:fitMatrix}
 \sum_{i=1}^{N} \begin{pmatrix}
 x_i x_i & x_i y_i & x_i z_i & x_i \\
 x_i y_i & y_i y_i & y_i z_i & y_i \\
 z_i x_i  & z_i y_i  &z_i z_i & z_i  \\
 x_i & y_i & z_i & 1
 \end{pmatrix}
 \begin{pmatrix}
 a \\ b \\ c \\ d
 \end{pmatrix}
 =
 \sum_{i=1}^{N} F_i\begin{pmatrix}
 x_i \\ y_i \\ z_i \\ 1
 \end{pmatrix},
 \end{equation}
where $F_i$ are the RDF cell values. The coefficients in the trilinear function are then found by solving the linear system and the gradient is obtained as $\nabla F =  (a,b,c)$. This gradient is assumed to be constant in the cell and so the interface normal estimate at the point $\mathbf x_S$ is simply $\hat{\mathbf n}_S = (a,b,c)/|(a,b,c)|$.

In step \ref{isoRDF:vertexRDF} of the \verb|iso-RDF-N| scheme we also use the least square interpolation method to obtain the volume fraction value in a vertex point based on the surrounding cell centre volume fractions. The calculation of the coefficients $a, b, c$ and $d$ are done with the vertex position as the origin and so the vertex value is simply equal to $d$.

\section{Numerical results for interface reconstruction}\label{Sec:Results}

The performance of the three proposed methods is now investigated by comparison to simple analytical solutions. First we reconstruct a circle on 2D meshes and a sphere on 3D meshes of varying type and resolution, and compare with results from literature. Thereafter, the new reconstruction methods are tested in combination with the isoAdvector advection procedure.

For our convergence tests we use the $L_{\infty}$ error norm which picks out the maximum error in the domain and therefore is the most difficult norm to get to converge. For the interface position, $\mathbf x_S$, and normal, $\hat{\mathbf n}_S$, inside a cell we define the two norms,
\begin{equation}\label{eq:ErrRecon}
	L_{\infty}^x= \max(|\mathbf x_S - \mathbf x_{\textrm{exact}}|) \quad \textrm{and} \quad
	L_{\infty}^n= \max(1 - \hat{\mathbf n}_S \cdot \mathbf n_{\textrm{exact}}),
\end{equation}
where the maximum is taken over all interface cells, $x_{\textrm{exact}}$ is the exact interface position, and $n_{exact}$ is the exact interface unit normal.

\subsection{Reconstruction of a circle on 2D meshes}

\begin{figure}[tb!]
		\includegraphics[width=\textwidth]{Meshes2D}
		\caption{Zoomed view of examples of 2D meshes used for reconstruction test.}
		\label{fig:Meshes2D}
\end{figure}

In the 2D reconstruction test the VOF field is calculated for a circular disc with radius 0.5 m placed at the centre of a square domain $[0, 2]\times [0, 2]$ m$^2$. The $L_\infty^x$ and $L_\infty^n$ error norms are calculated, and the exercise is repeated on increasingly fine meshes. All this is done for three mesh types consisting of square, triangular, and polygonal prisms, respectively, with the two latter mesh types being unstructured. The triangle prism meshes are generated with the open source mesh generator gmsh version 3.06. The polygonal prism meshes are generated by applying the OpenFOAM tool, \verb|polyDualMesh|, to the triangle meshes. Examples of the three types of meshes are shown in Fig.~\ref{fig:Meshes2D}.

\begin{figure}[tb!]
	\subcaptionbox{$L_{\infty}^n$ for triangular prisms
	\label{fig:Recon-Circle-tri-iter}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon-Circle-tri_iter}}
	\subcaptionbox{$L_{\infty}^x$ for triangular prisms\label{fig:Recon-Circle-tri-pos}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon-Circle-tri-pos}}
	\subcaptionbox{$L_{\infty}^n$ for square prisms\label{fig:Recon-Circle-hex}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon-Circle-hex}}
	\subcaptionbox{$L_{\infty}^x$ for square prisms\label{fig:Recon-Circle-hex-pos}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon-Circle-hex-pos}}
	\subcaptionbox{$L_{\infty}^n$ for polygonal prisms\label{fig:Recon-Circle-poly}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon-Circle-poly}}
	\subcaptionbox{$L_{\infty}^x$ for polygonal prisms\label{fig:Recon-Circle-poly-pos}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon-Circle-poly-pos}}
	\caption{$L_\infty$ for circle reconstructed on various mesh types.}
	\label{fig:Recon2D}
\end{figure}

Results are shown in Fig.~\ref{fig:Recon2D} with the $L_\infty^n$ in the left panels (a, c and e) and $L_\infty^x$ in the right panels (b, d and f). The first panel, Fig.~\ref{fig:Recon-Circle-tri-iter}, shows the orientation error on the unstructured triangular prism meshes for all three methods and for different number of iterations in the \verb|iso-RDF-N| and \verb|plic-RDF| methods. The first observation is that \verb|iso-Alpha| does not converge in this norm. The second observation is that both \verb|iso-RDF-N| and \verb|plic-RDF-N| can be brought to second order convergence rate on the tested refinement range if the number of iterations, $N$, is increased to 5. For lower $N$ both methods loose second order convergence at the finest meshes with \verb|iso-RDF-N| exhibiting the largest absolute errors. Based on these observations, we only show the $N=5$ results for the iterative methods in the remaining panels. In Fig.~\ref{fig:Recon-Circle-tri-pos} we see that $L_\infty^x$ for \verb|iso-Alpha|  falls off as $(\Delta x)^2$ on the coarsest triangular prism meshes but that the second order convergence is lost on the two finest of these meshes. Both iterative methods exhibit second order convergence on the full refinement range.

For the orientation error on square prism meshes in Fig.~\ref{fig:Recon-Circle-hex}, \verb|iso-Alpha| performs slightly better than on the triangular prism meshes, now converging albeit at a very slow rate. \verb|plic-RDF-5| and \verb|iso-RDF-5| both exhibit second order convergence on the full resolution range. As for the position error in Fig.~\ref{fig:Recon-Circle-hex-pos} all methods are second order on the full resolution range for square prism meshes.

The errors on unstructured polygonal prism meshes are shown in Figs.~\ref{fig:Recon-Circle-poly} and \ref{fig:Recon-Circle-poly-pos}. All three methods behave in a similar way as for the triangular prism meshes with the $L_\infty^n$ error of \verb|iso-Alpha| not converging and all other errors exhibiting second order convergence except $L_\infty^x$ for \verb|iso-Alpha| on the finest meshes.

In summary, the iterative RDF based methods show second order convergence of both interface orientation and position on all meshes types in 2D.

\begin{figure}[htb!]
	\subcaptionbox{$L_\infty^n$ for tetrahedral cells\label{fig:Recon3D-Circle-tri}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon3D-Circle-tri}}
	\subcaptionbox{$L_\infty^x$ for tetrahedral cells\label{fig:Recon3D-Circle-tri-pos}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon3D-Circle-tri-pos}}
	\subcaptionbox{$L_\infty^n$ for hexahedral cells\label{fig:Recon3D-Circle-hex}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon3D-Circle-hex}}
	\subcaptionbox{$L_\infty^x$ for hexahedral cells\label{fig:Recon3D-Circle-hex-pos}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon3D-Circle-hex-pos}}
	\subcaptionbox{$L_\infty^n$ for polyhedral cells\label{fig:Recon3D-Circle-poly}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon3D-Circle-poly}}
	\subcaptionbox{$L_\infty^x$ for polyhedral cells\label{fig:Recon3D-Circle-poly-pos}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon3D-Circle-poly-pos}}
	\caption{$L_\infty$ for reconstructed sphere on various mesh types.}\label{fig:Recontriandpoly}
\end{figure}

\subsection{Reconstruction of a sphere on 3D meshes}

To test the accuracy in 3D \verb|iso-Alpha|, \verb|plic-RDF-N| and \verb|iso-RDF-N| are applied to a sphere of radius 0.25 m centred in a unit cube domain covered by a range of mesh types and resolutions. The error norms are shown in Fig.~\ref{fig:Recontriandpoly} for unstructured tetrahedral meshes in the two top panels, structured hexahedral meshes in the two middle panels, and unstructured polyhedral meshes in the two bottom panels. The behaviour is very similar to the 2D test above: The orientation error of \verb|iso-Alpha|  converges slowly on the structured mesh but not on the unstructured meshes. The position error of \verb|iso-Alpha| falls off as $(\Delta x)^2$ on all meshes except for the finest unstructured meshes. Both \verb|iso-RDF-5| and \verb|plic-RDF-5| exhibit second order convergence in both interface orientation and position for all three mesh types on the full resolution range.

\begin{figure}[htb!]
	\subcaptionbox{Circle on triangular prisms\label{fig:Recon-Comp-tri}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon-Comp-tri}}
	\subcaptionbox{Sphere on hexahedral cells\label{fig:Recon3D-Comp-hex}}
	[.5\linewidth]{\includegraphics[width=0.5\textwidth]{Recon3D-Comp-hex}}
	\caption{$L_1^n$ reconstruction errors compared with literature.}	\label{fig:Recon3D-Comp}
\end{figure}

\subsection{Comparison with other methods}

We now compare the accuracy of the proposed reconstruction methods with the methods presented in Ito et al \cite{Ito.2014} and Ivey and Moin \cite{Ivey.2015}. In these papers the less strict 1-norm is used for the orientation error:
\begin{equation}\label{eq:ErrReconL1}
	L_1^n= \frac1N \sum_i |1 - \hat{\mathbf n}_{S,i}  \cdot \hat{\mathbf n}_{\textrm{exact},i}|,
\end{equation}
where $\hat{\mathbf n}_{Si}$ and $\hat{\mathbf n}_{\textrm{exact},i}$ are, respectively, the calculated and exact interface unit normal in interface cell $i$, and the sum is over all the $N$ interface cells. Ivey and Moin \cite{Ivey.2015} weighted the error in eq. \eqref{eq:ErrReconL1} with the volume, but on the cubic meshes below in comparison with their results all cell volumes are the same, and so the error definitions are identical for this case.

The 2D case is again a circle of radius 0.5 centred in the domain $[0, 2]\times [0, 2]$, which is covered by an unstructured mesh of triangular prisms. Fig.~\ref{fig:Recon-Comp-tri} depicts the $L_1^n$ error as it varies with the cell size (i.e. the circle radius divided by the root mean square of the average cell area). Also with this error norm \verb|iso-Alpha| shows lack of convergence in the orientation as does The Parker--Youngs method (PY). Both \verb|iso-RDF-5| and \verb|plic-RDF-5| achieve second order convergence as do the height function method and the geometric least square (GLS) method \cite{Shahbazi.2003}. The abolute errors of \verb|iso-RDF-5| and \verb|plic-RDF-5| are in between those of GSL and the height function method.

The 3D case is a sphere of radius 2 in a cubic box of side lenght 8. Fig.~\ref{fig:Recon3D-Comp-hex} shows the $L_1^n$ error variation with the mesh resolution (i.e. the radius of the sphere divided by the cubic root of number of cells). In contrast to PY, the \verb|iso-Alpha| method converges but at a slow rate. Both \verb|iso-RDF-5| and \verb|plic-RDF-5| achieve second order convergence as do the EHF method and LVIRA. The abolute errors of \verb|iso-RDF-5| and \verb|plic-RDF-5| are significantly lower than LVIRA's and similar to those of the height function method.

From a performance point of view the root finding step described in Section~\ref{ssec:rootFinding} is the most costly operation in our procedures due to the subcell reconstruction. \verb|iso-Alpha| only needs one such root finding step and is therefore the fastest of our proposed schemes. The methods \verb|iso-RDF-5| and \verb|plic-RDF-5| are roughly 15-20 times slower since the step is being performed multiple times. Computational performance depends on multiple factors and is as such difficult to compare between different methods and implementations. For structured meshes, the height function method is expected to be comparable in terms of performance to \verb|iso-Alpha| since only one subcell reconstruction step is performed. Ito et al. \cite{Ito.2014} and Ivey and Moin \cite{Ivey.2015} did not report absolute calculation times. However, the methods can be expected to be slower than \verb|plic-RDF-5| and \verb|iso-RDF-5|, since in addition to the subcell reconstruction step, multiple volume intersections between polyhedra are calculated. Ivey and Moin \cite{Ivey.2015} compared the performance of their EHF method to LVIRA and found that LVIRA is 20-100 times slower. Lopez et al. \cite{Lopez.2008} reported that LVIRA is 800 times slower than a method expected to be comparable in computational cost to \verb|iso-Alpha|. It is therefore reasonable to assume that \verb|plic-RDF-5| and \verb|iso-RDF-5| are significantly faster than LVIRA.

\section{Numerical results for interface advection}\label{sec:advectionResults}

We now investigate the effect of replacing \verb|iso-Alpha| with \verb|plic-RDF-N| in the isoAdvector algorithm described in \cite{Roenby.2016}.
As mentioned in Section \ref{ssec:initnS}, we have found that the number of required iterations in \verb|plic-RDF-N|, and thus the computational time, can be significantly reduced by improving the initial guess for $\hat{\mathbf n}_S$. In the following, we therefore use \verb|plic-RDF-N| with $N = 2$ and nearest-neighbour interpolation of $\hat{\mathbf n}_S$ from the previous time step as initial guess instead of the $\nabla\alpha$ values. For the error norm we will use

\begin{equation}\label{eq:E1AdvectError}
	E_{1}(t) = \frac{\sum_{i}V_{i}|\alpha_{i}(t)-\alpha_{i}^{exact}(t)|}{\sum_{i}V_{i}\alpha_{i}^{exact}(t)},
\end{equation}
where the sums are over all cells, and $\alpha_{i}^{exact}(t)$ is the volume fraction in cell $i$ corresponding to the known exact solution.

\subsection{Circle in constant flow}

\begin{figure}[ht]
	\subcaptionbox{Triangular prism mesh\label{fig:Recon-Circle-triIsoVsRDF}}
	{\includegraphics[width=0.5\linewidth]{advectDisc-triIsoVsRDF}}
	\subcaptionbox{Zoom of interface at $T = 4 s$\label{fig:interfaceAdvectCricleTri3}}
	{\includegraphics[width=0.5\linewidth]{interfaceAdvectCircleTri30}}
	\subcaptionbox{Square prism mesh\label{fig:Recon-Circle-hexIsoVsRDF}}
	{\includegraphics[width=0.5\linewidth]{advectDisc-hexIsoVsRDF}}
	\subcaptionbox{Polygonal prism mesh\label{fig:Recon-Circle-polyIsoVsRDF}}
	{\includegraphics[width=0.5\linewidth]{advectDisc-polyIsoVsRDF}}
	\caption{$E_1$ error for circle advected in constant flow.}\label{fig:RecontriIterandHex}
\end{figure}

Our first test case is a circular interface of radius 0.25 m initially centred at $(0.5,0.5)$ and moving with a constant, uniform velocity field without changing its shape. The advection test is performed on triangle, square and polygonal prism meshes. On the triangular and polygonal prism meshes the velocity is set to $(1.0, 0.0)$. On the square prism mesh the velocity is set to (1.0, 0.5) to avoid motion alignment with the mesh directions.

In Fig.~\ref{fig:RecontriIterandHex} we show the $E_{1}$ errors at time $t = 4$ s as functions of the mesh resolution for the different mesh types and performed with three different CFL numbers. Fig.~\ref{fig:Recon-Circle-triIsoVsRDF} shows the results for the triangular prism meshes. It is observed that \verb|iso-Alpha| falls of with a convergence rate between 1 and 2 for the coarsest meshes but that the rate drops to first order at finer meshes. The \verb|iso-Alpha| errors are independent of the time step size but are significantly larger than the \verb|plic-RDF-2| errors on all resolutions. Both the \verb|plic-RDF-2| results with CFL = 0.1 and CFL = 0.2 exhibit second order convergence in the full range of resolution. With CLF = 0.5 \verb|plic-RDF-2| falls off to first order at high resolutions but still has lower absolut error than \verb|iso-Alpha|. Fig.~\ref{fig:interfaceAdvectCricleTri3} shows a close--up of the interface at $t = 4$ s for the simulation with 10 cells per radius. It is evident that the interface segments of \verb|plic-RDF-N| (blue) are significantly improved  compared to \verb|iso-Alpha| (red).

Fig.~\ref{fig:Recon-Circle-hexIsoVsRDF} shows the $E_{1}$ errors as functions of mesh resolution for the square prism meshes. In contrast to the triangle prism results, \verb|iso-Alpha| improves with reduced time step size, especially for high resolutions. At lower resolutions \verb|iso-Alpha|  achieves second order convergence but at higher resolutions the convergence falls off, especially for large time steps. In contrast, \verb|plic-RDF-2| achieves second order convergence on the full resolution range with both CFL = 0.1 and 0.2.

Fig.~\ref{fig:Recon-Circle-polyIsoVsRDF} shows the behaviour on polygonal prism meshes. The behaviour is similar to the trianglar prism results with the \verb|iso-Alpha| giving larger absolute errors and dropping off to first order on finer meshes. Again, \verb|plic-RDF-2| retains second order convergence on the full resolution range with CFL = 0.1 and 0.2. For CFL = 0.5 the convergence rate of \verb|plic-RDF-2| drops at the finest meshes but from an absolute level which is interestingly somewhat smaller than the lower CFL results, and thus ending up with an absolute error very close to the CFL = 0.1 and 0.2 results on the finest mesh.

In summary, \verb|plic-RDF-2| in combination with the isoAdvector advection scheme from \cite{Roenby.2016} consistently lowers the absolute errors and gives second order convergence rate for small enough time step size.

\subsection{Disc in reversed spiral flow}\label{SpirallingDisc}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.9\linewidth]{vortexShearedDisc}
	\cprotect\caption{Spiralling disc test case on triangular prism mesh. Resulting interface shape after one complete reversal at $t = 8$ s is shown for \verb|plic-RDF-2| in blue and for \verb|iso-Alpha| in red toghether with the exact solution in green. The volume fraction field at the time of maximum stretching, $t = 4$ s, is shown in red and blue.}
	\label{fig:vortexShearedDisc}
\end{figure}

\begin{table}[h]
	\centering
	\begin{tabular}{l|c c c | c c c}

		&	& \verb|iso-Alpha| &  & & \verb|plic-RDF-2| & \\
		&	 square &  triangle &  polygon &  square & triangle & polygon\\
		\hline
		\hline
		coarse	&	0.044 	& 0.054 & 0.081 & 0.065 & 0.047 & 0.069 \\
				&	1.9		& 1.4 	&	2.4	&	2.2	&	2.3	&	2.4	\\
		medium	&	0.012 	& 0.02 	&  0.015& 0.014 & 0.0094&  0.013 \\
				&	2.5		& 	1.1	&	2.0	&	2.5	&	2.6	&	2.4 \\
		fine 	&	0.0021 	& 0.0095& 0.0037& 0.0025& 0.0016& 0.0024 \\
	\end{tabular}
	\caption{$E_1$ and convergence rates for a disc in reversed spiral flow with CFL = 0.5.}
	\label{tab:spirallingDisc}
\end{table}

We now move on to a case which has become a standard advection test case in the literature, namely a circular disc of radius 0.15 initialized centred at $(0.5,0.75)$ in the domain $[0,1]\times [0,1]$ advected in the velocity field
\begin{equation}
	\mathbf u(x,y,t) = \cos(2\pi t/16)\left( -\sin^2(\pi x) \sin(2 \pi y), \sin(2 \pi x) \sin^2 (\pi y) \right).
\end{equation}
This flow causes a stretching and spiraling of the disc to a maximum level reached at time $t = 4$ s as depicted in Fig.~\ref{fig:vortexShearedDisc}.  From time $t = 4$ s to 8 s, the flow is running backwards causing the fluid to end up at its initial configuration. The simulation is repeated for square, triangular and polygonal prism meshes of three different resolutions with the linear cell size halved at each refinement.

Table~\ref{tab:spirallingDisc} shows the final errors for both schemes together with the estimated order of convergence \cite{Tryggvason.2011},
\begin{equation}\label{eq:convRate}
	\mathcal O = \frac{\ln[E_1(\Delta x)/E_1(\Delta x/2)]}{\ln 2},
\end{equation}
where $E_1(\Delta x)$ is the error for the simulation with resolution $\Delta x$. Simulations were repeated with CFL = 0.1, 0.2 and 0.5 giving rise to almost identical results so we only show the CFL = 0.5 errors in Table~\ref{tab:spirallingDisc}.

We observe that \verb|plic-RDF-2| increases the absolute errors slightly on the square meshes compared to \verb|iso-Alpha|. On the triangular meshes \verb|plic-RDF-2| significantly reduces the errors compared to  \verb|iso-Alpha|, especially on the finer meshes. Finally, on the polygonal meshes \verb|plic-RDF-2| moderately reduces the absolute errors compared to \verb|iso-Alpha|. As for the convergence rates \verb|plic-RDF-2| is equal or better than \verb|iso-Alpha| for all cases with all rates in the range 2.2-2.6. In particular for the triangular mesh simulations \verb|plic-RDF-2| lowers the errors and increases the convergence rate significantly.

We note that with the resolutions used here and elsewhere in literature, the interface is not well-resolved at the time of maximum deformation of the circle. This is illustrated in Fig.~\ref{fig:vortexShearedDisc}, where pinched off ``droplets'' are clearly visible at the tail of the spiral. In the exact solution and in simulations where the thin tail is well-resolved this does not happen. Inspection of the simulations reveals that the fluid particles in the top part of the circular disc that is most distorted in Fig.~\ref{fig:vortexShearedDisc} are indeed the ones forming the pinching tail at time $t = 4$ s.

\subsection{Ellipsoid in constant flow}\label{ssec:EllipsoidInConstFlow}

Our final advection test case is an ellipsoidal interface with axes $(0.25,0.15,0.2)$ centred at $(0.5,0.5,0.5)$ in the domain $[0, 5] \times [0, 0.6] \times  [0, 0.6]$. The velocity field takes the constant value $(1, 0, 0)$ and the simulation runs until $t = 4$ s, where the exact interface solution has moved without changing its shape to the point $(4.5,0.5,0.5)$. The simulation is performed on three unstructured tetrahedral meshes with the linear cell size halved at every refinement. Simulations were performed with CLF = 0.1, 0.2 and 0.5 giving rise to essentially identical results, so we only show the CLF = 0.5 results here.

\begin{table}[h]
	\centering
	\begin{tabular}{l|c  c | c  c }
	$a / \Delta x$	&	\verb|iso-Alpha| 	&  & \verb|plic-RDF-2| &  \\
		\hline
		5  & 0.767 & - & 0.784 & - \\
		10 & 0.141 & 2.4 & 0.11 & 2.8 \\
		20 & 0.0359 & 2.0 & 0.0256 & 2.1 \\
		40 & 0.012 & 1.6 & 0.0063 & 2.0 \\
	\end{tabular}
	\caption{$E_1$ at time $T = 4$s for ellipsoid advected in constant flow for different mesh resolutions. The principal axes $a=0.25$ is scaled with the cell-to-cell centre distance $\Delta x$ of the mesh.}
	\label{tab:advectEllipsoid4}
\end{table}

The $E_1$ errors are shown in Table~\ref{tab:advectEllipsoid4} for \verb|iso-Alpha| and \verb|plic-RDF-2| together with convergence rates obtained with Eqn.~\ref{eq:convRate}. On the coarsest mesh the \verb|iso-Alpha| and \verb|plic-RDF-2| errors are almost identical. When the resolution is increased \verb|plic-RDF-2| converges faster than \verb|iso-Alpha| ending up with an absolute error which is half of the \verb|iso-Alpha| error on the finest mesh. This is a substantial improvement, but considering the large difference in convergence rate between \verb|iso-Alpha| and \verb|plic-RDF-2| on tetrahedral meshes depicted in Fig.~\ref{fig:Recon3D-Circle-tri} and the large local orientation errors of \verb|iso-Alpha| illustrated in Fig.~\ref{fig:Ellipsoid} it may seem a bit surprising that the difference in convergence rate is not bigger.
\begin{figure}[h]
	\centering
 		\includegraphics[width=0.6\textwidth]{Ellipsoid}
	\cprotect\caption{Ellipsoid reconstructed with \verb|iso-Alpha| (left) and with \verb|plic-RDF-2| (right).}
	\label{fig:Ellipsoid}
\end{figure}
The explanation is to be found in the fact that the $E_1$ error at a given time is a product of both reconstruction and advection error. On tetrahedral meshes the reconstruction errors are random in nature and thus do not accumulate with time. The advection error, on the other hand, is known to be systematic with a small overprediction of the advection speed as shown in \cite{Roenby.2016}. This systematic error decreases with mesh resolution, but it nevertheless leads to a systematic accumulation with simulation time as the exact ellipsoid falls behind the numerical ellipsoid. Thus, when the simulation time becomes long enough, the $E_1$ error will be dominated by the systematic advection error. To illustrate this, we plot the $E_1$ error as a function of time for the two finest meshes in Fig.~\ref{fig:errorOverTime}.
\begin{figure}[h]
	\centering
	\includegraphics[width=0.6\textwidth]{advection-error-time2}
	\cprotect\caption{Time evolution of $E_1$ for \verb|iso-Alpha| and \verb|plic-RDF-2|.}
	\label{fig:errorOverTime}
\end{figure}
On both meshes the error starts of with a significantly lower value when using \verb|plic-RDF-2| compared to \verb|iso-Alpha|. As time goes by the error steadily increases with a rate which is the same for the two methods, and which decreases with mesh refinement. If we calculate the convergence rate at the first time step, $t = 0.1$ s in Fig.~\ref{fig:errorOverTime}, instead of at $t = 4$ s, we get a convergence rate of 1.03 with \verb|iso-Alpha| and a convergence rate of 1.9 with \verb|plic-RDF-2|. This demonstrates the peculiar fact that an interface advection method may obtain second order convergence rate with mesh refinement even with a first order accurate reconstruction scheme. We remark that whereas the noisy interface obtained with \verb|iso-Alpha| may not be devastating for the interface advection step, it most likely will be for the calculation of interface curvature required in surface tension dominated flows.

We now investigate the overall performance of the schemes by using the open source profiling software, Valgrind, on the finest mesh case. We find that the advection step takes 50\% of the time with \verb|iso-Alpha|, 25\% of the time with \verb|plic-RDF-2| and 10\% of the time with \verb|plic-RDF-5|. In other words, the reconstruction plus advection takes twice as long with \verb|plic-RDF-2| and 5 times longer with \verb|plic-RDF-5| compared to \verb|iso-Alpha|. The advection step in Jofre et al. \cite{Jofre.2015} took 81\% of the total simulation time with a reconstruction step which we assume to be comparable in cost to \verb|iso-Alpha|. This means that isoAdvector is roughly twice as fast as Jofre et al. \cite{Jofre.2015}, which may be explained by the absence of polyhedron-polyhedron intersection calculations in isoAdvector. The scheme by Jofre et al \cite{Jofre.2015} is comparable in cost to those of Hernandez et al. \cite{Hernandez.2008}, Maric et al. \cite{Maric.} and Owkes and Desjardins \cite{Owkes.2014}.

\section{Conclusion}\label{Sec:Conclusion}

We have presented a new iterative interface reconstruction procedure based on a reconstructed signed distance function to extract the local interface position and orientation from the raw volume fraction data. The new method is demonstrated to exhibit second order convergence with mesh refinement on both 2D and 3D structured and unstructured meshes. The algorithm has been developed in two variants based on piecewise linear interface construction (PLIC) and RDF isosurface reconstruction. Especially on unstructured meshes and for the local interface orientation both methods have significantly improved convergence properties compared to the reconstruction method presented in \cite{Roenby.2016} which was based on isosurfaces of the volume fraction data. The implementation of the new reconstruction methods is straightforward and the increased computational cost caused by the introduction of iterations is limited to an acceptable level.

The new reconstruction scheme is combined with the isoAdvector advection step \cite{Roenby.2016}, and a number of simple advection test cases are investigated. The overall conclusion is that second order convergence with mesh refinement is achieved on all mesh types with the new reconstruction scheme giving better convergence and lower absolute errors especially on fine unstructured meshes.

The new methods are embedded in an OpenFOAM based code framework, which is released with this paper as open source\cite{Scheufler.2018}. The framework is structured in a way allowing for easy implementation and testing of new reconstruction and advection schemes and with a special emphasis on future extensions to adaptive mesh refinement. It is our hope that the released library will be used, tested and further developed by the CFD community, and eventually result in improved simulation quality in the broad field of applications involving interfacial flows.

\section*{Acknowledgement}

This work was supported by "German Aerospace Center - DLR".

\bibliographystyle{elsarticle-num}


text:
\section{Introduction}
	\setcounter{lemma}{0}
	\setcounter{theorem}{0}
	\setcounter{corollary}{0}
	\setcounter{remark}{0}
	\setcounter{equation}{0}
	\setcounter{conjecture}{0}

For each prime $p$, let $\mathbb{F}_p$ denote the finite field of $p$ elements. For any non-empty subsets $A_1,A_2,\cdots,A_k\subseteq\mathbb{F}_p$. Set
$$A_1+A_2+\cdots+A_k=\{a_1+\cdots+a_k:\ a_j\in A_j\ \text{for}\ j=1,\cdots,k\}.$$
Also, $|A|$ denotes the cardinality of a set $A$.

In 2012, S\'{a}rk\"{o}zy \cite{Sarkozy} investigated the additive decompositions concerning the set $R_p$ of all quadratic residues modulo $p$ and posed the following interesting conjecture in \cite[Conjecture 1.6]{Sarkozy}.

\begin{conjecture}[S\'{a}rk\"{o}zy]
$R_p$ has no decomposition of the form $R_p=A+B$ with $|A|,|B|\ge2$ if $p$ is sufficiently large.
\end{conjecture}
This conjecture seems beyond reach. On the other hand, under the assumption that $A+B=R_p$ with $|A|,|B|\ge2$, the bounds of $|A|$ and $|B|$ have been extensively investigated. For example, S\'{a}rk\"{o}zy \cite[Theorem 2.1]{Sarkozy} showed that if $p$ is large enough and $A+B=R_p$ with $|A|,|B|\ge2$, then
$$\frac{\sqrt{p}}{3\log{p}}<|A|,|B|<\sqrt{p}\log{p}.$$
The factor $\log{p}$ was later removed by Shkredov \cite{Shkredov}. Moreover, in the same paper Shkredov showed that $A+A\neq R_p$ for any $A$ with $|A|\ge2$. In 2021, Chen and Yan \cite[Theorem 1.1]{Chen} improved S\'{a}rk\"{o}zy's result and proved that
$$\frac{7-\sqrt{17}}{6}\sqrt{p}+1\le |A|,|B|\le \frac{7+\sqrt{17}}{4}\sqrt{p}-6.63$$
for each odd prime $p$. As a consequence of this improvement, Chen and Yan \cite[Theorem 1.2]{Chen} showed the following beautiful result:
$$A+B+C\neq R_p$$
for any $A,B,C\subseteq\mathbb{F}_p$ with $|A|,|B|,|C|\ge2$. This refines S\'{a}rk\"{o}zy's result \cite[Theorem 1.2]{Sarkozy}.
Recently, Chen and Xi \cite{Chen and Xi} obtained
$$\frac{1}{4}\sqrt{p}+\frac{1}{8}\le |A|,|B|\le 2\sqrt{p}-1$$
for each odd prime and studied the representation function
$$r_{A,B}(x)=|\{(a,b)\in A\times B:\ a+b=x\}|.$$

Motivated by the above results, in this paper, let $C_p$ be the set of all non-zero cubes in $\mathbb{F}_p$. We study several topics on additive decompositions of $C_p$. Note that if $p\equiv2\pmod3$, then $C_p=\mathbb{F}_p\setminus\{0\}$. Hence we just need consider the case $p\equiv1\pmod3$.

We now state our first result.

\begin{theorem}\label{Thm. A}
	Let $p\equiv1\pmod3$ be a prime with $p\ge9096$. Suppose $A+B=C_p$ with $|A|,|B|\ge2$. Then
	$$\frac{\sqrt{p}}{18}\le |A|,|B|\le 3\sqrt{p}+269.$$
	Moreover, if $p$ is large enough, then we further have
	$$(\frac{1}{9}-o(1))\sqrt{p}\le|A|,|B|.$$
\end{theorem}

As a direct consequence of Theorem \ref{Thm. A}, we obtain the following result concerning the $3$-decompositions of $C_p$.

\begin{theorem}\label{Thm. B}
	Let $p\equiv1\pmod3$ be a prime with $p>184291$. Then $C_p$ does not have a non-trivial $3$-decomposition, i.e.,
	$$A+B+C\neq C_p,$$
	for any $A,B,C\subseteq\mathbb{F}_p$ with $|A|,|B|,|C|\ge2$.
\end{theorem}

On the other hand, in the case $|A|=|B|$, we can obtain the following improvement of Theorem \ref{Thm. A}.
\begin{theorem}\label{Thm. C}
	Let $p\equiv1\pmod3$ be a prime. Suppose $A+B=C_p$ with $|A|=|B|$. Then
	$$\sqrt{\frac{p-1}{3}}\le |A| \le\sqrt{p}.$$
	Moreover, if $A+A=C_p$, then
	$$\frac{-1}{2}+\frac{\sqrt{3}}{6}\sqrt{8p-5}\le|A|\le\sqrt{p}.$$
\end{theorem}

Bachoc, Matolcsi and Ruzsa \cite{Bach}, as well as Shkredov \cite{Shkredov} independently obtained that if $A-A=R_p\cup\{0\}$, then
$$p\ge\begin{cases}
	|A|^2+|A|-1&\mbox{if}\ |A|\ \text{is even},\\
	|A|^2+2|A|-1&\mbox{otherwise}.
\end{cases} $$
Inspired by their results, our next result concerns $A-A=C_p\cup\{0\}$.

\begin{theorem}\label{Thm. D}
	Let $p\equiv1\pmod3$ be a prime. Suppose $A-A=C_p\cup\{0\}$. Then
	$$\sqrt{\frac{p+2}{3}}\le|A|
	\le\begin{cases}
		\sqrt{p}&\mbox{if}\ |A|\equiv 0\pmod 3,\\
		\frac{1+\sqrt{17}}{4}\sqrt{p}+\frac{1+\sqrt{17}}{4}&\mbox{otherwise.}
	\end{cases} $$
\end{theorem}
\begin{remark}
	Note that $\frac{1+\sqrt{17}}{4}=1.28\cdots$.
\end{remark}

The outline of this paper is as follows. In section 2, we will introduce some necessary lemmas for the proof of Theorem \ref{Thm. A}. The proofs of Theorems \ref{Thm. A}--\ref{Thm. B} will be given in section 3. In section 4, we will give some preparations for the proofs of Theorems \ref{Thm. C}--\ref{Thm. D}. The proofs of Theorems \ref{Thm. C}--\ref{Thm. D} will be given in section 5.

\section{Preparations for the Proof of Theorem \ref{Thm. A}}
\setcounter{lemma}{0}
\setcounter{theorem}{0}
\setcounter{corollary}{0}
\setcounter{remark}{0}
\setcounter{equation}{0}
\setcounter{conjecture}{0}

For each prime $p$, let $\mathbb{F}_p^{\times}=\{x\in\mathbb{F}_p:\ x\neq0\}$ and let $\widehat{\mathbb{F}_p^{\times}}$ be the set of all multiplicative characters of $\mathbb{F}_p$. In addition, let $\chi\in\widehat{\mathbb{F}_p^{\times}}$ be a character of order $3$ and define
$$\psi(x)=\chi(x)+\chi(x^2)$$
for any $x\in\mathbb{F}_p$. Then
$$\psi(x)=\begin{cases}
	2&\mbox{if}\ x\in C_p,\\
	0&\mbox{if}\ x=0,\\
	-1&\mbox{otherwise}.
\end{cases} $$
For simplicity, we write $\sum_{x}$ instead of $\sum_{x\in\mathbb{F}_p}$.

We begin with the following result which is known as the Weil Theorem (cf. \cite[Theorem 5.41]{LN}).

\begin{lemma}\label{Lem. Weil's theorem}{\rm (Weil's Theorem)}
	Let $\chi_p\in\widehat{\F_p^{\times}}$ be a character of order $m>1$ and let $f(x)\in\F_p[x]$ be a
	monic polynomial which is not of the form $g(x)^m$ for any $g(x)\in\F_p[x]$. Then for any $a\in\F_p$ we have
	\begin{equation*}
		\left|\sum_{x}\chi_p(af(x))\right|\le (r-1)p^{1/2},
	\end{equation*}
 where $r$ is the number of distinct roots of $f(x)$ in $\overline{\F_p}$.
\end{lemma}

Note that $C_{13}=\{1,5,8,12\}\subseteq\mathbb{F}_{13}$ has a decomposition $C_{13}=A+B,$
where $A=\{1,5\}$ and $B=\{0,7\}$. When $p>13$ we have the following result.

\begin{lemma}\label{Lem. k>2}
	Let $p\equiv1\pmod3$ be a prime greater than $13$. Suppose that $A+B=C_p$ with $|A|\ge|B|=k\ge2$. Then $k\ge3$.
\end{lemma}

\begin{proof}
	We prove this lemma by contradiction. Suppose $k=2$. Note that for any $y\in\mathbb{F}_p$ we have
	$$(A+y)+(B-y)=A+B.$$
	By this we may assume $0\in B$ and hence $A\subseteq C_p$. Set $B=\{0,b\}$. Given an element $x\in C_p$, if $x\in A$, then
	$$x+b\in A+B=C_p.$$
	If $x\not\in A$, since $A+B=C_p$, then
	$$x-b\in A\subseteq C_p.$$
	Hence for any $x\in C_p$ either $x+b\in C_p$ or $x-b\in C_p$. Let
	$$H(x)=(1+\psi(x))\cdot(2-\psi(x+b))\cdot(2-\psi(x-b)).$$
	Then $H(x)=0$ for any $x\in\mathbb{F}_p^{\times}$. Note also that
	$$0\le H(0)=\sum_{x}H(x)=(2-\psi(b))^2\le 9.$$

	On the other hand, $\sum_{x}H(x)$ is equal to
	\begin{align*}
		&\sum_{x}4+4\sum_{x}\psi(x)-2\sum_{x}\psi(x+b)
		-2\sum_{x}\psi(x-b)\\
		-2&\sum_{x}\psi(x)\psi(x+b)-2\sum_{x}\psi(x)\psi(x-b)+\sum_{x}\psi(x-b)\psi(x+b)\\
		+&\sum_{x}\psi(x)\psi(x-b)\psi(x+b).
	\end{align*}
	Via a computation, $\sum_{x}H(x)$ is equal to
	$$4p-4\sum_{x}\psi(x)\psi(x+b)+\sum_{x}\psi(x)\psi(x+2b)+\sum_{x}\psi(x)\psi(x+b)\psi(x+2b).$$
	Recall that $\psi(x)=\chi(x)+\chi(x^2)$. When $p>85$ by Lemma \ref{Lem. Weil's theorem} we obtain
	$$9\ge H(0)=\sum_{x}H(x)\ge 4p-36\sqrt{p}>9,$$
	which is a contradiction. When $p\le 85$ and $p\neq13$, by computations one can verify the desired result.
	This completes the proof.
\end{proof}

\begin{lemma}\label{Lem. another bound}
	Let $p\equiv1\pmod3$ be a prime. Suppose $A+B=C_p$ with $|A|\ge|B|=k\ge2$. Then
	$$\frac{p-1}{3k}\le |A| \le \frac{p+(3^{k-1}\cdot 2k-2^k+1)\sqrt{p}}{3^k}\le \frac{p}{3^k}+\frac{2k}{3}\sqrt{p}.$$
\end{lemma}

\begin{proof}
	As $A+B=C_p$, we have
	$$|A||B|\ge|A+B|=|C_p|=\frac{p-1}{3}$$
	and hence $|A|\ge\frac{p-1}{3k}$.

	On the other hand, set $B=\{b_1,b_2,\cdots,b_k\}$ and let
	$$h(x)=\frac{1}{3^k}\prod_{j=1}^k\left(1+\psi(x+b_j)\right).$$
	Note that $h(x)\ge0$ for any $x\in\mathbb{F}_p$ and $h(x)=1$ for any $x\in A$. By Lemma \ref{Lem. Weil's theorem} we obtain
	\begin{align*}
		|A|\le\sum_{x}h(x)&=\frac{1}{3^k}\sum_{x}\prod_{j=1}^k\left(1+\psi(x+b_j)\right)\\
		   &=\frac{1}{3^k}\sum_{x}\left(1+\sum_{j=1}^{k}\psi(x+b_j)
		   +\sum_{i=2}^{k}\sum_{1\le j_1<\cdots<j_i\le k}\prod_{i=1}^r\psi(x+b_{j_r})\right)\\
		   &=\frac{1}{3^k}\left(p+\sum_{i=2}^k\sum_{1\le j_1<\cdots<j_i\le k}\sum_{x}\prod_{r=1}^i\psi(x+b_{j_r})\right)\\
		   &\le\frac{1}{3^k}\left(p+\sqrt{p}\sum_{i=2}^k\binom{k}{i}2^i(i-1)\right)\\
		   &=\frac{1}{3^k}\left(p+(3^{k-1}\cdot 2k-3^k+1)\sqrt{p}\right)\le\frac{p}{3^k}+\frac{2k}{3}\sqrt{p}.
	\end{align*}

	This completes the proof.
\end{proof}

Note that this lemma implies that
$$\left(\frac{1}{3k}-\frac{1}{3^k}\right)p-\frac{2k}{3}\sqrt{p}-\frac{1}{3k}\le 0.$$
Hence we see that $A$ and $B$ are sufficiently large sets as $p\rightarrow\infty$.

We next introduce some necessary notations. Let
$$L(\mathbb{F}_p)=\{f:\ f\ \text{is a complex function on $\mathbb{F}_p$}\}.$$
For any $f\in L(\mathbb{F}_p)$ and a positive integer $k$, the function $C_{k+1}(f): \mathbb{F}_p^k\rightarrow\mathbb{C}$ is defined by
$$C_{k+1}(f)(x_1,\cdots,x_k)=\sum_{x}f(x)f(x+x_1)\cdots f(x+x_k).$$
In addition, for any $f,g\in L(\mathbb{F}_p)$ we define
$$(f\circ g)(x)=\sum_{y}f(y)g(x+y).$$
Shkredov \cite[Lemma 2.1]{Shkredov} obtained the following result.
\begin{lemma}[Shkredov]\label{Lem. Shkredov on C(f)}
For any $f,g\in L(\mathbb{F}_p)$ and a positive integer $k$,
	\begin{equation*}
		\sum_{x}(f\circ g)^{k+1}(x)
		=\sum_{x_1,\cdots,x_k\in\mathbb{F}_p}C_{k+1}(f)(x_1,\cdots,x_k)\cdot C_{k+1}(g)(x_1,\cdots,x_k).
	\end{equation*}
\end{lemma}

\section{Proofs of Theorems \ref{Thm. A}--\ref{Thm. B}}
\setcounter{lemma}{0}
\setcounter{theorem}{0}
\setcounter{corollary}{0}
\setcounter{remark}{0}
\setcounter{equation}{0}
\setcounter{conjecture}{0}

For any $A\subseteq\mathbb{F}_p$, the characteristic function of $A$ is defined by
$$A(x)=\begin{cases}
	1&\mbox{if}\ x\in A\\
	0&\mbox{otherwise.}
\end{cases}$$

{\bf Proof of Theorem \ref{Thm. A}.} Suppose $A+B=C_p$ with $|A|,|B|\ge2$. Then
\begin{align*}
	\sum_{x\in B}(A\circ\psi)^4(x)
	&=\sum_{x\in B}\sum_{y_1,y_2,y_3,y_4\in\mathbb{F}_p}\prod_{j=1}^4A(y_j)\psi(y_j+x)\\
	&=\sum_{x\in B}\sum_{y_1,y_2,y_3,y_4\in A}\prod_{j=1}^42=16|A|^4|B|.
\end{align*}
Hence
\begin{equation}\label{Eq. A in the proof of Thm. A}
16|A|^4|B|=\sum_{x\in B}(A\circ\psi)^4(x)\le\sum_{x}(A\circ\psi)^4(x).
\end{equation}

On the other hand, by Lemma \ref{Lem. Shkredov on C(f)} we have
$$\sum_{x}(A\circ\psi)^4(x)=\sum_{x_1,x_2,x_3\in\mathbb{F}_p}C_4(A)(x_1,x_2,x_3)\cdot C_4(\psi)(x_1,x_2,x_3).$$
Let
$$E=\{(x,x,0),(x,0,x),(0,x,x):\ x\in\mathbb{F}_p^{\times}\}\cup\{(0,0,0)\}.$$
By Lemma \ref{Lem. Weil's theorem} if $(x_1,x_2,x_3)\not\in E$, then
\begin{equation}\label{Eq. B in the proof of Thm. A}
	|C_4(\psi)(x_1,x_2,x_3)|\le48\sqrt{p}.
\end{equation}
If $(x_1,x_2,x_3)\in E\setminus\{(0,0,0)\}$, then there is an element $t\in\mathbb{F}_p^{\times}$ such that
\begin{align*}
	C_{4}(\psi)(x_1,x_2,x_3)
	&=\sum_{x}\psi^2(x)\psi^2(x+t)\\
	&=\sum_{x}(2\chi_0(x)+\psi(x))(2\chi_0(x+t)+\psi(x+t))\\
	&=4(p-2)+\sum_{x}\psi(x)\psi(x+t).
\end{align*}
By this and Lemma \ref{Lem. another bound} one can verify that
\begin{equation}\label{Eq. B1 in the proof of Thm. A}
	C_4(\psi)(x_1,x_2,x_3)\le\begin{cases}
		4p+4\sqrt{p}&\mbox{if}\ (x_1,x_2,x_3)\in E\setminus\{(0,0,0)\},\\
	6p&\mbox{if}\ (x_1,x_2,x_3)=(0,0,0).	\end{cases}
\end{equation}

By (\ref{Eq. A in the proof of Thm. A})--(\ref{Eq. B1 in the proof of Thm. A}) we obtain
\begin{align*}
	16|A|^4|B|&\le\sum_{x_1,x_2,x_3\in\mathbb{F}_p}C_4(A)(x_1,x_2,x_3)\cdot C_4(\psi)(x_1,x_2,x_3)\\
	          &\le 48\sqrt{p}|A|^4+3(4p+4\sqrt{p})|A|^2+6p|A|,
\end{align*}
i.e.,
\begin{equation}\label{Eq. C in the proof of Thm. A}
	|A|^2|B|\le 3\sqrt{p}|A|^2+\frac{3}{4}(p+\sqrt{p})+\frac{3p}{8|A|}\le 3\sqrt{p}|A|^2+\frac{3}{4}(p+\sqrt{p})+\frac{p}{8}.
\end{equation}
The last inequality follows from Lemma \ref{Lem. k>2}.

By this and $|A||B|\ge|C_p|$ we obtain
\begin{equation}\label{Eq. D in the proof of Thm. A}
	\frac{p-1}{3}|A|\le 3\sqrt{p}|A|^2+\frac{7}{8}p+\frac{3}{4}\sqrt{p}.
\end{equation}
When $p\ge9096$ one can verify that
$$0\le\frac{\frac{p-1}{3}-\sqrt{\Delta}}{6\sqrt{p}}<6,$$
where
$$\Delta=\left(\frac{p-1}{3}\right)^2-3\sqrt{p}\left(\frac{7}{2}p+3\sqrt{p}\right).$$
On the other hand, note that by Lemma \ref{Lem. another bound} we have $|A|,|B|\ge6$ if $p\ge 9096$.
Hence the quadratic inequality (\ref{Eq. D in the proof of Thm. A}) indeed implies that for $p\ge9096$ we have
\begin{equation}\label{Eq. E in the proof of Thm. A}
	|A|\ge\frac{\frac{p-1}{3}+\sqrt{\Delta}}{6\sqrt{p}}\ge\frac{\sqrt{p}}{18}.
\end{equation}
By this and (\ref{Eq. C in the proof of Thm. A}) again we further obtain that if $p\ge9096$, then
\begin{equation}\label{Eq. F in the proof of Thm. A}
	|B|\le 3\sqrt{p}+\frac{3p+3\sqrt{p}}{4|A|^2}+\frac{3p}{8|A|^3}\le 3\sqrt{p}+269.
\end{equation}
With the method as above, we can also obtain that $|B|\ge \frac{\sqrt{p}}{18}$
and $|A|\le 3\sqrt{p}+269$ if $p\ge9096$.

Moreover, if $p$ is large enough, then by (\ref{Eq. E in the proof of Thm. A}) and the above method one can verify that
\begin{equation}\label{Eq. F in the proof of Thm. A}
	|A|,|B|\ge (\frac{1}{9}-o(1))\sqrt{p}.
\end{equation}

In view of the above, we have completed the proof.\qed

To prove our next result, we need the following result which is a special case of \cite[Theorem 1.2]{GMR}.
\begin{lemma}[Gyarmati, Matolcsi and Ruzsa]\label{Lem. GMR}
	Let $p$ be an odd prime and let $A_1,A_2,\cdots, A_k\subseteq\mathbb{F}_p$ be non-empty subsets. Then
	$$\left|\sum_{j=1}^kA_j\right|^{k-1}\le\prod_{j=1}^k\left|\sum_{i\neq j}A_i\right|.$$
\end{lemma}

{\bf Proof of Theorem \ref{Thm. B}.} We prove this theorem by contradiction. Now suppose $A+B+C=C_p$ with $|A|,|B|,|C|\ge2$. Then by Lemma \ref{Lem. GMR} and Theorem \ref{Thm. A} we have
$$\left(\frac{p-1}{3}\right)^2=|A+B+C|^2\le|A+B|\cdot|B+C|\cdot|C+A|\le\left(3\sqrt{p}+379\right)^3.$$
Howeover, when $p>184291$ we have
$$\left(\frac{p-1}{3}\right)^2>\left(3\sqrt{p}+269\right)^3,$$
which deduces a contradiction.

In view of the above, we have completed the proof.\qed

\section{Preparations for the Proofs of Theorem \ref{Thm. C}--\ref{Thm. D}}
\setcounter{lemma}{0}
\setcounter{theorem}{0}
\setcounter{corollary}{0}
\setcounter{remark}{0}
\setcounter{equation}{0}
\setcounter{conjecture}{0}

We first introduce some necessary notations. For any $f,g\in L(\mathbb{F}_p)$, the inner product $\langle f,g\rangle$ of $f,g$ is defined by
$$\langle f,g\rangle=\sum_{x}f(x)\overline{g(x)}.$$
Set
$$||f||_2=(\langle f,f\rangle)^{1/2}=\left(\sum_{x}|f(x)|^2\right)^{1/2},$$
and
$$\langle f\rangle=\langle f,1\rangle=\sum_{x}f(x).$$
Also, the convolution $f*g$ is defined by
$$(f*g)(x)=\sum_{y}f(y)g(x-y).$$
Note that
\begin{equation}\label{Eq. relations between circ and convolution}
(f*g)(x)=(f^c\circ g)(x),
\end{equation}
where $f^c(x)=f(-x)$.

\begin{lemma}\label{Lem. computation of inner product}
	Let $p\equiv1\pmod3$ be a prime and let $f,g\in L(\mathbb{F}_p)$. Then  $\langle f\circ\psi, g\circ\psi\rangle$ is equal to
	\begin{align*}
	 2p\langle f,g\rangle-2\langle f\rangle\cdot\langle \bar{g}\rangle+J(\chi,\chi)\langle\chi^2,\bar{f}\circ g\rangle+J(\chi^2,\chi^2)\langle\chi,\bar{f}\circ g\rangle,
	\end{align*}
where $J(\chi,\chi)$ is the Jacobi sum.
\end{lemma}
\begin{proof}
	By definition
\begin{align*}
		\langle f\circ\psi,  g\circ\psi\rangle
		&=\sum_{x}(f\circ\psi)(x)\cdot\overline{( g\circ\psi)(x)}\\
		&=\sum_{y,z}f(y)\overline{g(z)}\sum_{x}\psi(x+y)\psi(x+z).
\end{align*}
If $y=z$, then
\begin{align*}
	\sum_{x}\psi(x+y)\psi(x+z)=\sum_{x}\psi(x)^2=\frac{4(p-1)}{3}+\frac{2(p-1)}{3}=2(p-1).
\end{align*}
If $y\neq z$, writing $t=y-z\neq0$, then
\begin{align*}
	\sum_{x}\psi(x+y)\psi(x+z)
	&=\sum_{x}\psi(x)\psi(x+t)\\
	&=\chi^2(t)J(\chi,\chi)+\chi(t)\overline{J(\chi,\chi)}+2J(\chi,\chi^2)\\
	&=\chi^2(t)J(\chi,\chi)+\chi(t)\overline{J(\chi,\chi)}-2.
\end{align*}
The last equality follows from $J(\chi,\chi^2)=-1$.
By above results, one can verify that $\langle f\circ\psi, g\circ\psi\rangle$ is equal to
\begin{align*}
	2p\langle f,g\rangle-2\langle f\rangle\cdot\langle \bar{g}\rangle+J(\chi,\chi)\langle\chi^2,\bar{f}\circ g\rangle+J(\chi^2,\chi^2)\langle\chi,\bar{f}\circ g\rangle.
\end{align*}
This completes the proof.
\end{proof}

\section{Proofs of Theorems \ref{Thm. C}--\ref{Thm. D}}
\setcounter{lemma}{0}
\setcounter{theorem}{0}
\setcounter{corollary}{0}
\setcounter{remark}{0}
\setcounter{equation}{0}
\setcounter{conjecture}{0}

{\bf Proof of Theorem \ref{Thm. C}.} Suppose $A+B=C_p$. Set $|A|=a,|B|=b$. Define a function $r\in L(\mathbb{F}_p)$ by
$$(A\circ\psi)(x)=2aB(x)+r(x)$$
for any $x\in\mathbb{F}_p$. For any $x\in B$, we observe that
$$(A\circ\psi)(x)=\sum_{y}A(y)\psi(y+x)=\sum_{y\in A}\psi(x+y)=2a=2aB(x),$$
i.e., $r(x)=0$ for any $x\in B$. Hence
$$||r||_2^2=\sum_{x\not\in B}|r(x)|^2=\sum_{x\not\in B}|(A\circ\psi)(x)|^2=||A\circ\psi||_2^2-\sum_{x\in B}|(A\circ\psi)(x)|^2.$$
By Lemma \ref{Lem. computation of inner product}
\begin{align*}
	||A\circ\psi||_2^2
	&=2p\langle A,A\rangle-2\langle A\rangle^2+2\Re\left(J(\chi,\chi)\langle\chi^2,A\circ A\rangle\right)\\
	&=2pa-2a^2+2\Re\left(J(\chi,\chi)\langle\chi^2,A\circ A\rangle\right),
\end{align*}
where $\Re(z)$ is the real part of a complex number $z$. On the other hand,
$$\left|\Re\left(J(\chi,\chi)\langle\chi^2,A\circ A\rangle\right)\right|\le |J(\chi,\chi)|\cdot|\langle\chi^2,A\circ A\rangle|\le a^2\sqrt{p}.$$
The last inequality follows from $|J(\chi,\chi)|=\sqrt{p}$ and
\begin{align*}
|\langle\chi^2,A\circ A\rangle|
=\left|\sum_{x}\chi^2(x)(A\circ A)(x)\right|
&=\left|\sum_{x}\chi^2(x)\sum_{y}A(y)A(x+y)\right|\\
&=\left|\sum_{x,y\in A}\chi^2(x-y)\right|\le a^2.
\end{align*}
Note also that
\begin{align*}
\sum_{x\in B}|(A\circ\psi)(x)|^2
&=\sum_{x\in B}\sum_{y}A(y)\psi(x+y)\sum_{z}A(z)\psi(x+z)\\
&=\sum_{x\in B}\sum_{y,z\in A}\psi(x+y)\psi(x+z)=4a^2b.
\end{align*}
By the above results, we have
\begin{equation}\label{Eq. a in the proof of Thm. C}
	||r||^2_2=||A\circ\psi||_2^2-\sum_{x\in B}|(A\circ\psi)(x)|^2\le 2ap+2a^2\sqrt{p}-2a^2-4a^2b.
\end{equation}

On the other hand,
\begin{align*}
	\langle r\rangle=\sum_{x}r(x)=\sum_{x}(A\circ\psi)(x)-2a\sum_{x}B(x)=-2ab.
\end{align*}
The last equality follows from
$$\sum_{x}(A\circ\psi)(x)=\sum_{x}\sum_{y}A(y)\psi(x+y)=\sum_{y}A(y)\sum_{x}\psi(x+y)=0.$$
By the Cauchy-Schwarz inequality
$$|\langle r\rangle|^2=\left|\sum_{x\not\in B}r(x)\right|^2\le\sum_{x\not\in B}1^2\sum_{x\not\in B}r(x)^2=(p-b)\cdot||r||_2^2.$$
Combining this with (\ref{Eq. a in the proof of Thm. C}), we obtain
\begin{equation}\label{Eq. b in the proof of Thm. C}
	a+b+2ab\le p+a\sqrt{p}+\left(\frac{1}{p}-\frac{1}{\sqrt{p}}\right)ab\le p+a\sqrt{p}.
\end{equation}

Now suppose $|A|=|B|=a$. By (\ref{Eq. b in the proof of Thm. C}) we obtain
$$2a^2+2a\le p+a\sqrt{p}.$$
This implies that
$$a\le\frac{\sqrt{p}-2+\sqrt{(\sqrt{p}-2)^2+8p}}{4}<\sqrt{p}.$$
Moreover, if $A=B$, then
$$\binom{a}{2}+a\ge |C_p|.$$
This implies
$$a\ge\frac{-1}{2}+\frac{\sqrt{3}}{6}\sqrt{8p-5}.$$

In view of the above, we have completed the proof.\qed

To prove our last result, we need the following result due to Shkredov \cite[Lemma 3.1]{Shkredov}.
\begin{lemma}[Shkredov]\label{Lem. Shkredov trick}
	Let $c$ be an integer and let $f:\mathbb{F}_p\rightarrow\mathbb{Z}$ be a function. Then
	\begin{equation}\label{Eq. a1 of Shkredov}
		||f||_2^2\ge c|\langle f\rangle|-(c-1)\cdot\left|\sum_{x: 0<|f(x)|<c}f(x)\right|.
	\end{equation}
Further,
\begin{equation}\label{Eq. a2 of Shkredov}
	||f||_2^2=c\langle f\rangle+\sum_{k\in\mathbb{Z}}N_k(k^2-ck),
\end{equation}
where $N_k=|\{x\in\mathbb{F}_p:\ f(x)=k\}|$.
\end{lemma}

Now we prove our last result.

{\bf Proof of Theorem \ref{Thm. D}.} Set $|A|=a$. The case $a=1$ is trivial. We now suppose $a\ge2$ and
$A-A=C_p\cup\{0\}$. Note that
for any $y\in\mathbb{F}_p$ we have
$$(A-y)-(A-y)=A-A.$$
By this we may assume $0\in A$ and hence $\emptyset\ne A\setminus\{0\}\subseteq C_p$. Let $s\in L(\mathbb{F}_p)$ defined by
$$(A*\psi)(x)=2(a-1)A(x)+s(x)$$
for any $x\in\mathbb{F}_p$. Then for any $x\in A$ we have
$$(A*\psi)(x)=\sum_{y}A(y)\psi(x-y)=\sum_{y\in A}\psi(x-y)=2(a-1)=2(a-1)A(x),$$
i.e., $s(x)=0$ for any $x\in A$.  By (\ref{Eq. relations between circ and convolution}) we have
$A*\psi=A^c\circ\psi$. Hence by the method appeared in the proof of Theorem \ref{Thm. C} we obtain
\begin{equation}\label{Eq. 2 norm of s}
	||s||_2^2\le 2ap+2a^2\sqrt{p}-2a^2-4(a-1)^2a
\end{equation}
and
\begin{equation}\label{Eq. langle s rangle}
\langle s\rangle=-2a(a-1).
\end{equation}
Moreover, as $||s||_2\ge 0$, we obtain a trivial upper bound of $a$, i.e.,
\begin{equation}\label{Eq. lower bound of p}
	\sqrt{p}\ge \frac{-a+\sqrt{9a^2-12a+8}}{2}= \frac{-a+\sqrt{5a^2+4a^2-12a+8}}{2}\ge \frac{\sqrt{5}-1}{2}a.
\end{equation}

On the other hand, for any $x\not\in A$ we have
\begin{align*}
	s(x)=\sum_{y}A(y)\psi(x-y)
	&=\sum_{y\in A}\psi(x-y)\\
	&=2|C_p\cap(x-A)|-|D_p\cap(x-A)|\\
	&=3|C_p\cap(x-A)|-a,
\end{align*}
where $D_p=\{x\in\mathbb{F}_p:\ \psi(x)=-1\}$. The above result implies that
$$s(x)\equiv -a\pmod3$$
for any $x\not\in A$.

{\bf Case 1.} $a\equiv0\pmod3$.

Applying $c=3$ to (\ref{Eq. a1 of Shkredov}), we obtain $$||s||_2^2\ge 3|\langle s\rangle|.$$
By (\ref{Eq. 2 norm of s}) and (\ref{Eq. langle s rangle}) we have
$$2a^2-\sqrt{p}a-p-1\le 0$$
and hence
$$a\le \frac{\sqrt{p}+\sqrt{9p-8}}{4}<\sqrt{p}.$$

{\bf Case 2.} $a\not\equiv0\pmod3$.

In this case, as $s(x)\equiv-a\pmod3$ for any $x\not\in A$, we have
$$s(x)=0\Leftrightarrow x\in A.$$
Hence $N_0=a$. On the other hand, note that
$$s(x)=-1\Leftrightarrow\sum_{y\in A}\psi(x-y)=-1\Leftrightarrow G(x)=-1,$$
where
$$G(x)=\sum_{y\in A}\left((x-y)^{\frac{p-1}{3}}+(x-y)^{\frac{2(p-1)}{3}}\right).$$
Hence
$$N_{-1}\le \frac{2(p-1)}{3}.$$
Applying $c=-2$ to (\ref{Eq. a2 of Shkredov}) we obtain
\begin{align*}
	||s||_2^2
	&=-2\langle s\rangle+\sum_{k\in\mathbb{Z}}N_k(k^2+2k)\\
	&=4a(a-1)-N_{-1}+\sum_{k\neq -1}N_k(k^2+2k)\\
	&\ge 4a(a-1)-N_{-1}+3\sum_{k\neq -1,0}N_k\\
	&=4a(a-1)-4N_{-1}-3N_0+3p\\
	&\ge 4a(a-1)-\frac{8(p-1)}{3}-3a+3p.
\end{align*}
This gives
$$4a^2-(2+2\sqrt{p})a-2p-3\le \frac{p+8}{3a}.$$
As $a^2\ge 1+|C_p|$, by (\ref{Eq. 2 norm of s}) we obtain
\begin{align*}
4a^2-(2+2\sqrt{p})a-2p-3\le \frac{p+8}{3a}
&\le \frac{p+8}{\sqrt{3}\sqrt{p+2}}\\
&= \frac{\sqrt{p+2}}{\sqrt{3}}+\frac{2\sqrt{3}}{\sqrt{p+2}}\\
&\le \frac{\sqrt{p+2}}{\sqrt{3}}+\frac{2}{\sqrt{3}}.
\end{align*}
By this we finally obtain
\begin{align*}
	a
	&\le \frac{2+2\sqrt{p}+\sqrt{(2+2\sqrt{p})^2+16(2p+3+\frac{\sqrt{p+2}}{\sqrt{3}}+\frac{2}{\sqrt{3}})}}{8}\\
	&\le \frac{2+2\sqrt{p}+\sqrt{(2+2\sqrt{p})^2+16(2+2\sqrt{p})^2}}{8}\\
	&\le \frac{1+\sqrt{17}}{4}\sqrt{p}+\frac{1+\sqrt{17}}{4}.
\end{align*}

In view of the above, we have completed the proof.\qed

\Ack\ We would like to thank Prof. Hao Pan for his helpful suggestions and steadfast encouragement.

redpajama_set_name: RedPajamaArXiv