LBM-MHFEM chapter - added problem formulation

\section{Problem formulation}

\label{sec:lbm-mhfem:problem formulation}

\cref{eq:ns,eq:transport}

\inline{equations: copy Navier--Stokes here, add a general advection-diffusion equation (without physical interpretation)}

\inline{define domains $\Omega_1$ and $\Omega_2$ -- use the problem from the appendix in \cite{klinkovsky2022:WT}}

The flow of an incompressible fluid is governed by the Navier--Stokes equations written in the standard form as

\begin{subequations}\label{eq:lbm-mhfem:ns}

    \begin{equation}\label{eq:lbm-mhfem:ns:mass}

        \nabla \cdot \vec{v} = 0,

    \end{equation}

    \begin{equation}\label{eq:lbm-mhfem:ns:momentum}

        %\rho \left( \frac{\partial \vec{v}}{\partial t} + \vec{v}\cdot\nabla\vec{v}\right) + \nabla p = %\mu  \Delta \vec{v} + \vec{g},

        \frac{\partial \vec{v}}{\partial t} + \vec{v}\cdot\nabla\vec{v} = - \frac{1}{\rho} \nabla p + \nu \Delta \vec{v},

    \end{equation}

\end{subequations}

where $\vec v = [v_x, v_y, v_z]$~[\si{\metre\per\second}] is the fluid velocity, $p$~[\si{\pascal}] is the pressure, $\rho$~[\si{\kg\per\cubic\metre}] is the fluid density, and $\nu$~[\si{\metre\squared\per\second}] is the kinematic viscosity.

All these quantities are functions of spatial coordinates $\vec x = [x, y, z] \in \Omega_1\subset\mathbb{R}^3$ and time $t \in (0, t_{\max})$.

The mass and momentum conservation laws \eqref{eq:lbm-mhfem:ns} are coupled with a generic advection--diffusion equation without sources or sinks that may be written in the conservative form

\begin{subequations}\label{eq:ADE}

    \begin{equation}\label{eq:ADE:conservative}

        \frac{\partial \phi}{\partial t} + \nabla \cdot \left( \phi \vec v - D_0 \nabla \phi \right) = 0, %f_\phi.

    \end{equation}

    where $\phi$ is a variable that may be associated to the physical quantity transported by the fluid (e.g., molar or mass concentration for mass transport, or temperature for heat transport) and $D_0$~[\si{\metre\squared\per\second}] is the corresponding diffusion coefficient.

    Combining \cref{eq:ADE:conservative,eq:lbm-mhfem:ns:mass} leads to the non-conservative form

    \begin{equation}\label{eq:ADE:non-conservative}

        \frac{\partial \phi}{\partial t} + \vec v \cdot \nabla \phi - \nabla \cdot \left( D_0 \nabla \phi \right) = 0. %f_\phi,

    \end{equation}

\end{subequations}

\Cref{eq:ADE:conservative,eq:ADE:non-conservative} are assumed to hold in $\Omega_2 \times (0, t_{\max})$, where $\Omega_2 \subset \Omega_1$.

Since \cref{eq:lbm-mhfem:ns,eq:ADE} are coupled only via the velocity field $\vec v$, \cref{eq:lbm-mhfem:ns} can be solved without \cref{eq:ADE} as the values of $\phi$ do not influence the flow field.

In this work, we are interested in modeling constituent transport in a channel flow.

As a benchmark problem, we consider the situation illustrated in \cref{fig:lbm-mhfem:domain} where $\Omega_1 = [0, 1.75] \times [0, 1] \times [0, 1]$ (dimensions are in meters) is a cuboidal channel with inflow, outflow, and solid wall boundaries.

The constituent transport is tracked in the subdomain $\Omega_2 = [0.5, 1.5] \times [0.25, 0.75] \times [0.25, 0.75]$ (in meters).

For simplicity, we assume that $\Omega_2$ is completely immersed in the domain $\Omega_1$ (i.e., none of the domain boundaries coincide: $\partial\Omega_1 \cap \partial\Omega_2 = \emptyset$).

The kinematic viscosity $\nu = \SI{15.52e-6}{\metre\squared\per\second}$ and the diffusion coefficient $D_0 = \SI{25.52e-6}{\metre\squared\per\second}$ are set the same as in \cref{chapter:vapor transport} (see \cref{tab:physical parameters}).

Both \cref{eq:lbm-mhfem:ns,eq:ADE} must be supplemented by suitable initial and boundary conditions.

For simplicity, the velocity field is initialized by zero ($\vec v(\vec x, 0) = \vec 0$) and the variable $\phi$ is initialized as $\phi(\vec x, 0) = \phi_0$, where $\phi_0$ is assumed to be constant.

The boundary conditions on $\partial \Omega_1$ are posed as follows:

\begin{itemize}

    \item

        Solid walls (i.e., the top, bottom, front, and back sides of $\Omega_1$) are modeled using the standard no-slip boundary condition.

    \item

        A fixed value of pressure and a zero velocity gradient in the normal direction are prescribed on the outflow side (i.e., right hand side) of domain $\Omega_1$.

    \item

        Fixed values of velocity are prescribed on the inflow side (i.e., left hand side) of $\Omega_1$.

        Details will be described later.

\end{itemize}

The conditions for the variable $\phi$ on the subdomain boundary $\partial \Omega_2$ are:

\begin{itemize}

    \item

        A Dirichlet-type condition is used to prescribe fixed values on the inflow side (i.e., left hand side) of $\Omega_2$.

    \item

        A Neumann-type condition is used to prescribe zero gradient in the normal direction on all remaining sides of $\Omega_2$.

        %($\frac{\partial \phi}{\partial \vec x} \cdot \vec n = 0$)

\end{itemize}

The concrete inflow boundary profiles for both $\vec v$ and $\phi$ will be specified later in \cref{sec:lbm-mhfem:numerical analysis} where the benchmark problem will be studied numerically.

\begin{figure}[!tb]

\begin{figure}[tb]

    \centering

    \includegraphics[width=0.8\textwidth]{figures/lbm-mhfem/ADE2_domain_annotated.png}

    \caption{Schematic configuration of the computational domains $\Omega_1$ and $\Omega_2$ (2D cross-section along the plane $y = 0$).}

    \caption{Schematic configuration of the computational domains $\Omega_1$ and $\Omega_2$ (2D cross-section of a 3D cuboidal channel along the plane $y = 0$).}

    \label{fig:lbm-mhfem:domain}

\end{figure}

\section{Computational algorithm and time adaptivity}

\label{sec:lbm-mhfem:algorithm}

The previous two sections described the numerical methods used for the discretizations of \cref{eq:ns,eq:transport}, respectively.

The previous two sections described the numerical methods used for the discretizations of \cref{eq:lbm-mhfem:ns,eq:ADE}, respectively.

In this section, we describe the main computational algorithm of the coupled scheme focusing on its adaptive time control.

The initialization of the solver includes setting all physical parameters, discretization of the domain $\Omega_1$ by the lattice $\hat \Omega$, discretization of the domain $\Omega_2$ by the unstructured mesh $\mathcal K_h$, and decomposition of the lattice and mesh into subdomains for distributed computing.

\section{Interpolation of the velocity field}

\label{sec:interpolation}

Since \cref{eq:ns,eq:transport} are coupled by the velocity field $\vec v(\vec x, t)$, the numerical approach relies on the interpolation of the approximate velocity field computed by LBM and its projection into the finite element space used by MHFEM.

Since \cref{eq:lbm-mhfem:ns,eq:ADE} are coupled by the velocity field $\vec v(\vec x, t)$, the numerical approach relies on the interpolation of the approximate velocity field computed by LBM and its projection into the finite element space used by MHFEM.

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$.

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:transport:non-conservative,eq:transport: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: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:transport:conservative} would lead to numerical oscillations and instability of the solver.

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:transport:non-conservative} is less sensitive to errors in the velocity field and the post-processing schemes are not necessary in this case.

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.

\section{Implementation remarks}

The result of this decomposition procedure is illustrated in \cref{fig:domain_decomposition}b.

Overall, the decomposition algorithm optimizes the computational cost and memory requirements of each MPI rank at the cost of increased communication due to increased number of lattice subdomains.

\inline{Add the 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}

\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}

\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.

The aim of this section is to study the differences between the conservative and non-conservative formulations of the transport equation.

\Cref{eq:ns} governing the fluid flow is solved in a cuboidal channel $\Omega_1 = [0, 1.75] \times [0, 1] \times [0, 1]$  (dimensions are in meters) with parameters similar to the main problem discussed in the paper (kinematic viscosity $\nu = \SI{15.52e-6}{\metre\squared\per\second}$, mean inflow velocity magnitude $v_{\max} = \SI{1}{\metre\per\second}$).

Note that the channel is free of all obstacles, but we induce turbulent flow using the unsteady (time-varying) inflow boundary condition described in \cref{sec:inflow:fluctuations}.

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}.

The aim of this section is to study the differences between the conservative and non-conservative formulations of the advection--diffusion equation \eqref{eq:ADE}.

The fluid flow is coupled with a transport equation either in the conservative form \cref{eq:transport:conservative} or non-conservative form \cref{eq:transport:non-conservative}, where $\phi$~[-] is exempt from its physical meaning and for the purpose of this section, it is interpreted as the concentration of a generic constituent transported by the fluid.

The diffusion coefficient $D = \SI{25.52e-6}{\metre\squared\per\second}$ is set the same as in the main problem discussed in the paper (see \cref{tab:physical parameters}).

The transport equation is solved in domain $\Omega_2 = [0.5, 1.5] \times [0.25, 0.75] \times [0.25, 0.75]$ (in meters) that is completely immersed in the domain $\Omega_1$ (i.e., none of the domain boundaries coincide: $\partial\Omega_1 \cap \partial\Omega_2 = \emptyset$).

See \cref{fig:lbm-mhfem:domain} for schematic configuration of the domains.

In order to study the differences between the conservative and non-conservative formulations, the initial and boundary conditions outlined in \cref{sec:lbm-mhfem:problem formulation} are specified with the following profiles.

Initially, we set $\phi(\vec x, 0) = \phi_0 = 1$ for all $\vec x \in \Omega_2$.

On the inflow boundary of $\Omega_2$ ($x = \SI{0.5}{\metre}$), we prescribe a fixed value $\phi_{\mathrm{in}} = 1$.

Turbulent flow in $\Omega_1$ is induced using the unsteady (time-varying) inflow boundary condition described in \cref{sec:inflow:fluctuations}\todo{move the section to this chapter to avoid a forward reference} with the mean inflow velocity magnitude $v_{\max} = \SI{1}{\metre\per\second}$.\todo{the symbol $v_{\max}$ is used in the power-law profile, not in the fluctuations -- change it to $\bar{\vec v}_{\mathrm{in}}$?}

Given a divergence-free velocity field due to \cref{eq:lbm-mhfem:ns:mass}, this initial-boundary-value problem has a trivial analytical solution $\phi(\vec x, t) = 1$ for all $\vec x \in \Omega_2$ and $t > 0$.

In order to study the differences between the conservative and non-conservative formulations, the initial and boundary conditions for the transport equation are posed as follows.

Initially, we set $\phi = 1$ uniformly in the whole domain $\Omega_2$.

On the inflow boundary ($x = 0.5$), we prescribe a fixed value $\phi = 1$.

On all remaining parts of $\partial \Omega_2$, we prescribe a zero gradient in the normal direction ($\frac{\partial \phi}{\partial \vec x} \cdot \vec n = 0$).

Given a divergence-free velocity field due to \cref{eq:ns:mass}, this initial-boundary-value problem has a trivial analytical solution $\phi(\vec x, t) = 1$ for all $\vec x \in \Omega_2$ and $t > 0$.

The coupled problem is solved numerically using the LBM-MHFEM scheme as described in \cref{sec:computational approach}.

Several variants of the MHFEM scheme were used, namely explicit or implicit upwind, and linear or cubic interpolation of the velocity field.

Several variants of the MHFEM scheme from \cref{chapter:MHFEM} were used, namely explicit or implicit upwind, and linear or cubic interpolation of the velocity field.

Each variant was computed in three resolutions denoted as RES-A1, RES-A2, and RES-A3, see \cref{tab:ADE2:resolutions}.

Note that single precision was used in the LBM part for fluid flow and double precision was used in the MHFEM part.

To illustrate the turbulent flow field in $\Omega_1$, \cref{fig:ADE2:vx} shows the horizontal velocity ($v_x$) field in the final time $t_{\max} = \SI{10}{\second}$.

\Cref{fig:ADE2:concentrations} shows qualitative differences between the concentration ($\phi$) fields that were computed using different variants of the MHFEM scheme.

\Cref{fig:ADE2:concentrations} shows qualitative differences between the fields of the transported variable $\phi$ that were computed using different variants of the MHFEM scheme.

Since the fields obtained using any variant with the non-conservative formulation were visually indistinguishable from the constant analytical solution on the scale used in \cref{fig:ADE2:concentrations}, only the conservative formulation variants are shown in the figure.

Note that for given resolution, the velocity field is the same in all variants of the MHFEM scheme.

Quantitative comparison is presented in \cref{tab:ADE2:norms} in terms of $L^p$ norms of the differences between the analytical solution $\phi = 1$ and each numerical solution $\phi_h$.

Both qualitative and quantitative results in \cref{fig:ADE2:concentrations,tab:ADE2:norms} indicate that for the conservative formulation, changing linear interpolation to cubic, as well as changing the explicit upwind discretization to implicit upwind, leads to smoother and more accurate results.

Furthermore, all these variants converge to the analytical solution as the lattice and grid are refined.

However, even the most accurate numerical solution obtained using the conservative formulation exhibits an error that is larger by orders of magnitude compared to the non-conservative formulation, even in the coarsest resolution.

The only difference between the discretizations of the non-conservative and conservative formulations is in \cref{eq:mhfem:advection terms discrete} where the former contains a term corresponding to the discrete divergence of velocity.

The results indicate that this extra term can be understood as a compensation for the non-zero divergence of the discrete velocity field interpolated on the mesh.

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.

\begin{table}[!tb]

\begin{table}[tb]

    \centering

    \caption{Characteristics of lattice and grid resolutions used for the numerical analysis.}

    \label{tab:ADE2:resolutions}

    \begin{tabular}{lrrr}

        \toprule

        & RES-A1 & RES-A2 & RES-A3 \\

        No. of lattice sites & approx. $3.5 \cdot 10^6$ & approx. $29 \cdot 10^6$ & approx. $233 \cdot 10^6$ \\

        No. of grid cells & approx. $0.5 \cdot 10^6$ & approx. $4 \cdot 10^6$ & approx. $33 \cdot 10^6$ \\

        Base time step $\Delta t$ & \SI{1.39e-3}{\second} & \SI{3.38e-4}{\second} & \SI{8.32e-5}{\second} \\

        \multicolumn{1}{p{13em}}{Average no. of LBM iters per MHFEM step ($\floor{C_{\max} / C}$)} & 2 & 4 & 9 \\

        \multicolumn{1}{p{12em}}{Average no. of LBM iters per MHFEM step ($\floor{C_{\max} / C}$)} & 2 & 4 & 9 \\

        \bottomrule

    \end{tabular}

    \caption{Characteristics of lattice and grid resolutions used for the numerical analysis.}

    \label{tab:ADE2:resolutions}

\end{table}

\begin{table}[!tb]

Both qualitative and quantitative results in \cref{fig:ADE2:concentrations,tab:ADE2:norms} indicate that for the conservative formulation, changing linear interpolation to cubic, as well as changing the explicit upwind discretization to implicit upwind, leads to smoother and more accurate results.

Furthermore, all these variants converge to the analytical solution as the lattice and grid are refined.

However, even the most accurate numerical solution obtained using the conservative formulation exhibits an error that is larger by orders of magnitude compared to the non-conservative formulation, even in the coarsest resolution.

The only difference between the discretizations of the non-conservative and conservative formulations is in \cref{eq:mhfem:advection terms discrete}\todo{fix reference when the term is added to \cref{sec:mhfem:numerical scheme}} where the former contains a term corresponding to the discrete divergence of velocity.

The results indicate that this extra term can be understood as a compensation for the non-zero divergence of the discrete velocity field interpolated on the mesh.

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.

\begin{table}[tb]

    \centering

    \input{data/lbm-mhfem/norms_concentration.tex}

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

    \label{tab:ADE2:norms}

    \input{data/lbm-mhfem/norms_concentration.tex}

\end{table}

\begin{figure}[!tb]

\begin{figure}[tb]

    \centering

    \includegraphics[width=0.8\textwidth]{data/lbm-mhfem/vx/lbmres=8_ictype=turb_interp=linear_ML=enabled_advection=explicit-upwind_transport=conservative.png}

    \\\vspace{1ex}

    \label{fig:ADE2:vx}

\end{figure}

\begin{figure}[!tb]

\begin{figure}[tb]

    \centering

    \begin{subfigure}[b]{0.495\textwidth}

        \centering

    \\\vspace{1ex}

    \includegraphics[width=\textwidth]{data/lbm-mhfem/cbar_concentration.png}

    \caption{

        Simulated concentration field ($\phi$) along the plane $y = 0$ in $\Omega_2$ in the benchmark problem using the \emph{conservative formulation} of the transport equation \eqref{eq:transport:conservative}.

        Simulated fields of the transported variable $\phi$ along the plane $y = 0$ in $\Omega_2$ in the benchmark problem using the \emph{conservative formulation} of the transport equation \eqref{eq:ADE:conservative}.

        Several configurations of the numerical scheme are compared: linear and cubic interpolation of the velocity from LBM to MHFEM, and discretization of the advection term in the MHFEM scheme based on explicit and implicit upwind.

        Only the first two resolutions RES-A1 and RES-A2 are shown here.

    }