LBM-MHFEM chapter - updated discussion on the linear vs cubic interpolation...

Note that the spatial discretization of the domain $\Omega_2 \subset \Omega_1$ is generally different than the equidistant lattice on $\Omega_1$; it may be a regular grid with different space steps or even an unstructured mesh.

The interpolation of the velocity field can be requested at any point $\vec x \in \Omega_2$.

The surrounding lattice points $\hat{\vec x} \in \hat\Omega$ can be easily found and the linear or cubic interpolation in $\mathbb R^3$ can be used to obtain the velocity at $\vec x$ from the velocities at $\hat{\vec x}$.

The surrounding lattice points $\hat{\vec x} \in \mathcal L_{\overline{\Omega}_1}$ can be easily found and the linear or cubic interpolation in $\mathbb R^3$ can be used to obtain the velocity at $\vec x$ from the velocities at $\hat{\vec x}$.

Note that linear interpolation can be implemented more efficiently than cubic interpolation as it uses fewer input data points.

Our implementation of the cubic interpolation causes the whole solver to run about three times slower than when using the linear interpolation.

Although cubic interpolation is more accurate, the results presented in \cref{sec:results} indicate that linear interpolation is sufficient for practical use.

Our implementation of the cubic interpolation is not efficient and may cause the whole solver to run multiple times slower compared to the linear interpolation.

The impact of using the linear or cubic interpolation on the accuracy of the numerical solution is investigated in \cref{sec:lbm-mhfem:numerical analysis}.

The finite element space used by MHFEM imposes some requirements on the interpolation of the velocity field.

The finite element space used by MHFEM imposes requirements on the interpolation of the velocity field.

According to \cite{brezzi:1991mixed}, the Raviart--Thomas--Nédélec space of the lowest order $\mathbf{RTN}_0(\mathcal K_h)$ that is used for the finite element--approximation of the velocity field in this work is formed by functions $\vec \omega \in [L^2(\Omega_2)]^3$ such that:

\begin{enumerate}

    \item for any element $K \in \mathcal K_h$, the restriction of $\vec \omega$ to $K$, $\vec \omega|_K$, must belong to the finite element space $\mathbf{RTN}_0(K)$ on the element $K$,

    \item $\vec \omega$ satisfies the balancing condition for normal traces on all interior faces $E \in \mathcal E_h^{\mathrm{int}}$ of the mesh, i.e.,

    \item

        $\vec \omega$ satisfies the balancing condition for normal traces on all interior faces $E \in \mathcal E_h^{\mathrm{int}}$ of the mesh, i.e.,

        $\int_E \vec \omega|_{K_1} \cdot \vec n_{K_1,E} + \int_E \vec \omega|_{K_2} \cdot \vec n_{K_2,E} = 0$

        for all $E \in \mathcal E_h^{\mathrm{int}}$, $E \in \mathcal E_{K_1} \cap \mathcal E_{K_2}$, where $\vec n_{K_\iota,E}$ is the unit normal vector on the face $E$ oriented outward from the element $K_\iota$, $\iota=1,2$.

        for all $E \in \mathcal E_h^{\mathrm{int}}$, $E \in \mathcal E_{K_1} \cap \mathcal E_{K_2}$, where $\vec n_{K_\ell,E}$ is the unit normal vector on the face $E$ oriented outward from the element $K_\ell$, $\ell=1,2$.

\end{enumerate}

An interpolation strategy compatible with these requirements is as follows.

First, velocity is evaluated at the element face centres $\vec x_E$ for all $E \in \mathcal E_h$.

This can be done at any time level $t$ yielding the approximate velocity values $\hat{\vec v}(\vec x_E, t)$ which are then used for the projection into the $\mathbf{RTN}_0(\mathcal K_h)$ space.

The discrete velocity field is assumed to be piecewise constant on the element sides $\mathcal E_h$ and the values $\hat{\vec v}(\vec x_E, t)$ define the components corresponding to the face $E$ in the finite element spaces $\mathbf{RTN}_0(K_j)$ of the elements $K_j$ adjacent to the face $E$.

Finally, it is important to note that numerical schemes for \cref{eq:ADE:non-conservative,eq:ADE:conservative} do not behave equivalently with a general discrete velocity field interpolated to the mesh.

This is because the discrete velocity field computed by LBM may not satisfy \cref{eq:lbm-mhfem:ns:mass} exactly and the field interpolated to the unstructured mesh may be locally non-conservative.

Thus, solving the transport equation in the conservative form \eqref{eq:ADE:conservative} would lead to numerical oscillations and instability of the solver.

The problem of compatibility between flow schemes producing a discrete velocity field and transport schemes using the interpolated velocity was extensively researched and several velocity post-processing algorithms were developed to recover the compatibility \cite{larson2004,sun2006,kees2008}.

%\later{This deserves further investigation, but there is not enough time...}

However, such post-processing would incur additional cost to the computational algorithm.

On the other hand, we observed that solving the transport equation in the non-conservative form \eqref{eq:ADE:non-conservative} is less sensitive to errors in the velocity field and the post-processing schemes are not necessary in this case.

First, velocity is evaluated at the element face centers $\vec x_E$ for all $E \in \mathcal E_h$.

This can be done at any time level $t$ yielding the approximate velocity values $\hat{\vec v}_E \equiv \hat{\vec v}(\vec x_E, t)$ which are then used for the projection into the $\mathbf{RTN}_0(\mathcal K_h)$ space.

The discrete velocity field is assumed to be piecewise constant on the element sides $\mathcal E_h$ and the values $\hat{\vec v}_E$ define the components corresponding to the face $E$ in the finite element spaces $\mathbf{RTN}_0(K_j)$ of the elements $K_j$ adjacent to the face $E$.

Finally, it is important to note that the MHFEM schemes for \cref{eq:ADE:non-conservative,eq:ADE:conservative} do not behave equivalently with a general discrete velocity field interpolated to the mesh.

This is because the discrete velocity field computed by LBM may not satisfy \cref{eq:lbm-mhfem:ns:mass} exactly and even if it did, the interpolation scheme combines values from different locations in the flow field on a single element.

Hence, the field interpolated to the unstructured mesh may be locally non-conservative, i.e., the discrete approximation of the velocity divergence $\sum\limits_{E \in \mathcal E_K} \hat{\vec v}_E \cdot \vec n_{K,E}$ on element $K \in \mathcal K_h$ may be non-zero.

The accuracy of the numerical scheme applied to the conservative form of \cref{eq:ADE:conservative} and non-conservative form of \cref{eq:ADE:non-conservative} is investigated in \cref{sec:lbm-mhfem:numerical analysis} on a benchmark problem.

\section{Implementation remarks}

\later[inline]{Future work: problem of mapping MPI ranks to GPUs -- quadratic assignment problem, plus we need to get the weights (communication cost between each pair of GPUs) somehow.}

\section{Numerical analysis of the conservative vs. non-conservative formulation}

\section{Numerical analysis}

\label{sec:lbm-mhfem:numerical analysis}

In this section, we study numerically the convergence of the coupled LBM-MHFEM scheme using an artificial benchmark problem described in \cref{sec:lbm-mhfem:problem formulation}.

Furthermore, it can be noticed in \cref{tab:ADE2:norms} that changing the interpolation or upwind scheme does not have a significant effect on the error when the non-conservative formulation is used.

In the finest resolution RES-A3, using the linear interpolation and explicit upwind is not only advantageous for the performance of the solver, but also leads to a smaller error.

The presented results show that solving the transport equation in the conservative form of \cref{eq:ADE:conservative} with a highly turbulent velocity field may lead to large deviations in the numerical solution, whereas solving the non-conservative form of \cref{eq:ADE:non-conservative} with the same velocity field results in significantly more accurate solution.

An alternative approach to address the problem of non-zero divergence of the discrete velocity field might be to use a post-processing algorithm to recover the discrete divergence-free condition on the given mesh.

The problem of compatibility between flow schemes producing a discrete velocity field and transport schemes using the interpolated velocity was extensively researched and several velocity post-processing algorithms were developed \cite{larson2004,sun2006,kees2008}.

However, such post-processing would incur additional cost to the computational algorithm and the approach is not investigated further in this thesis.

\begin{table}[tb]

    \centering

    \caption{Results of the numerical analysis for different formulations and variants of the MHFEM scheme.}

We have found that the best performance is obtained using the BiCGstab method combined with the Jacobi preconditioner, both implemented in TNL \cite{oberhuber:2021tnl}.

In the simulations presented in this chapter, the BiCGstab method took at most 4 iterations to converge in most of the time steps, so improved performance cannot be expected from stronger preconditioners.

Based on the results from \cref{sec:lbm-mhfem:numerical analysis}, we use the linear interpolation of the velocity field and the explicit upwind scheme in MHFEM.

\todo{Add that the non-conservative formulation of the transport equation is used for the numerical solution.}

\section{Validation results}

\label{sec:results}