LBM-MHFEM chapter - updated computational algorithm section (38a5c354) · Commits · Jakub Klinkovský / Dissertation

content/LBM-MHFEM.tex

+50 −43

Original line number	Diff line number	Diff line
		@@ -77,60 +77,67 @@ The concrete inflow boundary profiles for both $\vec v$ and $\phi$ will be speci
		\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: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 Navier--Stokes equations \eqref{eq:lbm-mhfem:ns} are solved numerically by the lattice Boltzmann method described in \cref{chapter:LBM}.
		The transport equation \eqref{eq:ADE} is incorporated into the mathematical framework of the \emph{NumDwarf} solver described in \cref{chapter:MHFEM}.
		Both \cref{eq:ADE:conservative} and \cref{eq:ADE:non-conservative} can be converted to the form of \cref{eq:NumDwarf} by taking $n=1$, $\vec Z = [\phi]$, $\matrix N = [1]$, $\vec m = [1]$, $\vec D = [D_0]$, $\vec w = \vec 0$, $\vec r = \vec 0$, $\vec f = [0]$, %[f_\phi]$,
		and depending on the equation:
		\begin{itemize}[itemsep=0ex]
		\item $\vec u = \vec 0$, $\vec a = [\vec v]$ for \cref{eq:ADE:conservative},
		\item $\vec u = [\vec v]$, $\vec a = \vec 0$ for \cref{eq:ADE:non-conservative}.
		\end{itemize}

		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.
		The resulting solver couples the computational algorithms of MHFEM (\cref{sec:mhfem:computational algorithm}) and LBM (\cref{sec:lbm:computational algorithm}).
		The initialization of the solver includes setting all physical parameters, discretization of the domain $\Omega_1$ by the regular lattice $\mathcal L_{\overline{\Omega}_1}$ (see \cref{sec:lbm:lattice}), discretization of the domain $\Omega_2$ by the conforming unstructured mesh $\mathcal K_h$ (see \cref{sec:mhfem:numerical scheme}), and decomposition of the lattice and mesh into subdomains for distributed computing (described later in \cref{sec:lbm-mhfem:decomposition}).
		Then the solver allocates all data structures, applies the initial conditions and starts the main time-loop.
		To optimize the efficiency of the solver, we developed an adaptive time-stepping algorithm based on a "CFL-like" condition.
		The time-stepping part of the computational algorithm can be summarized as follows:
		The time-stepping part of the computational algorithm is summarized in \cref{alg:LBM-MHFEM:time-stepping}.

		\begin{enumerate}
		%\item Initialize the solver for given resolution:
		%\begin{enumerate}
		% \item Set all physical parameters, define discretizations of both domains $\Omega_1$ and $\Omega_2$.
		% \item Decompose the lattice and unstructured mesh into $N_p$ subdomains and assign them to the available computational units.
		% \item Set the data for each boundary condition used in the problem.
		% \item Allocate data for the solver and apply the initial conditions.
		%\end{enumerate}
		\item Set the time tracking variables $t_L := 0$ and $t_M := 0$, physical time step $\Delta t$ and final time $t_{\max} = N_t \Delta t$.
		\begin{algorithm}[Time-stepping in the coupled LBM-MHFEM scheme]
		\label{alg:LBM-MHFEM:time-stepping}
		\begin{algsteps}
		\item Set $t_L := \SI{0}{\second}$ and $t_M := \SI{0}{\second}$, physical time step $\delta_t$~[\si{\second}] and final time $t_{\max}$~[\si{\second}].
		\item \label{step:time-loop}
		While $t_L < N_t \Delta t$, repeat these steps:
		\begin{enumerate}
		While $t_L < t_{\max}$:
		\begin{algsteps}
		\item %If $i \pmod{1000} = 0$, %% NOTE: i is undefined...
		After every 1000 iterations, recompute inflow velocity fluctuations that will be used in the next 1000 iterations.
		\item Perform one iteration of LBM. %\todo{Add more details? I.e. compute on boundaries, compute on internal blocks, wait for communication on boundaries, start MPI communication, wait for communication, wait for computation}
		\item Set $t_L := t_L + \Delta t$.
		\todo{Velocity fluctuations were not introduced yet...}
		\item Perform one iteration of LBM on the lattice $\mathcal L_{\overline{\Omega}_1}$ (i.e., perform the steps \ref{algitem:forcing} to \ref{algitem:output} from \cref{alg:LBM} on \cpageref{alg:LBM}).
		\item Set $t_L := t_L + \delta_t$.
		\item \label{step:mesh-comp}
		If $t_M < t_L$, perform these steps:
		\begin{enumerate}
		If $t_M < t_L$:
		\begin{algsteps}
		\item Interpolate the velocity field from the regular lattice to the unstructured mesh.
		See \cref{sec:interpolation} for details regarding this procedure.
		\item Compute $C = \max_E\{\abs{\vec v_E} \Delta t / \abs{E}\}$, where $E \in \mathcal E_h$ goes over all faces of the unstructured mesh, $\abs{\vec v_E}$ is the magnitude of the interpolated velocity on the face $E$ and $\abs{E}$ is the characteristic length of the face $E$.
		\item Set the time step for MHFEM: $\Delta t_M := \Delta t \floor{C_{\max} / C}$ if $C \leq C_{\max}$, else $\Delta t_M := \Delta t / \ceil{C / C_{\max}}$.
		See \cref{sec:lbm-mhfem:interpolation} for details.
		\item Compute $C = \max_E\{\abs{\hat{\vec v}_E} \delta_t / h_E\}$, where $E \in \mathcal E_h$ goes over all faces of the mesh, $\abs{\hat{\vec v}_E}$ is the interpolated velocity magnitude on face $E$ and $h_E$~[\si{\metre}] is the mesh size on face $E$.
		\item Set the time step for MHFEM: $\delta_{t,M} := \delta_t \floor{C_{\max} / C}$ if $C \leq C_{\max}$, else $\delta_{t,M} := \delta_t / \ceil{C / C_{\max}}$.
		\item Set the number of MHFEM iterations: $n_M := 1$ if $C \leq C_{\max}$, else $n_M := \ceil{C/C_{\max}}$.
		\item Perform $n_M$ iterations of MHFEM with the time step $\Delta t_M$.
		The MHFEM algorithm is described in \cite{fucik:2019NumDwarf}.
		\item Set $t_M := t_M + n_M \Delta t_M$.
		\end{enumerate}
		\item If $t_L$ is a multiple of a pre-defined time period, make a snapshot of the current data for visualization/post-processing.
		\end{enumerate}
		\end{enumerate}
		\item Perform $n_M$ iterations of MHFEM using \cref{alg:MHFEM} with the time step $\delta_{t,M}$.
		\item Set $t_M := t_M + n_M \delta_{t,M}$.
		\end{algsteps}
		\item If $t_L$ or $t_M$ reached a pre-defined time period, make a snapshot of the current data for visualization/post-processing.
		\end{algsteps}
		\end{algsteps}
		\end{algorithm}

		Initially, two separate time tracking variables $t_L$ and $t_M$ are created, one for computations on the lattice and the other for computations on the mesh.
		The lattice variable $t_L$ is the main one which is checked in the time loop condition (step~\ref{step:time-loop} above).
		The variable $t_M$ is incremented separately in the step~\ref{step:mesh-comp} that is responsible for the coupling between computations on the lattice and the mesh.
		Here the velocity field is interpolated from the lattice to the mesh and the time-step control factor $C$ is computed.
		The time step $\Delta t_M$ for the mesh computations can be either longer or shorter than $\Delta t$ depending on the value of $C$.
		If $C \leq C_{\max}$, where $C_{\max}$ is a pre-defined constant, the time step $\Delta t_M$ is set as an integral multiple of $\Delta t$ and the time variable $t_M$ is pushed forward to the new time level by just one iteration of MHFEM.
		Here the velocity field is interpolated from the lattice to the mesh (the procedure is described later in \cref{sec:lbm-mhfem:interpolation}) and the time-step control factor $C$ is computed.
		The time step $\delta_{t,M}$ for the mesh computations can be either longer or shorter than $\delta_t$ depending on the value of $C$.
		If $C \leq C_{\max}$, where $C_{\max}$ is a pre-defined constant, the time step $\delta_{t,M}$ is set as an integral multiple of $\delta_t$ and the time variable $t_M$ is pushed forward to the new time level by just one iteration of MHFEM.
		Hence, several following LBM iterations can be performed successively until the lattice time $t_L$ overruns $t_M$ and the condition~\ref{step:mesh-comp} is satisfied.
		On the other hand, if $C > C_{\max}$, the MHFEM solver needs to perform $n_M > 1$ iterations to move to the new time level $t_M + \Delta t$ with a shorter time step $\Delta t_M = \Delta t / n_M$.
		On the other hand, if $C > C_{\max}$, the MHFEM solver needs to perform $n_M > 1$ iterations to move to the new time level $t_M + \delta_t$ with a shorter time step $\delta_{t,M} = \delta_t / n_M$.

		In practice, we found empirically that limiting the time step $\delta_{t,M}$ with $C_{\max} = 1/2$ is necessary to ensure the stability of the coupled solver.
		The overall performance of the solver depends on the concrete values of $\delta_{t,M}$ selected by the adaptive algorithm, which are influenced by the quantities needed to compute the time-step control factor $C$, i.e., local velocity magnitude, lattice time step $\delta_t$ (which is related to the lattice space step $\delta_x$), and the local space step $h_E$ of the unstructured mesh.

		In practice, we found empirically that limiting the time step $\Delta t_M$ with $C_{\max} = \frac{1}{2}$ is necessary to ensure the stability of the coupled solver.
		The overall performance of the solver depends on the concrete values of $\Delta t_M$ selected by the adaptive algorithm, which are influenced by the quantities needed to compute the time-step control factor $C$, i.e., local velocity magnitude, lattice time step (which is related to the lattice space step), and the local space step of the unstructured mesh (expressed by $\abs{E}$ in the algorithm).
		An important choice related to the solver performance is the selection of the algorithm for the solution of large systems of linear equations arising from the MHFEM discretization.
		We achieved the best performance using the BiCGstab method combined with the Jacobi preconditioner, both implemented in TNL \cite{oberhuber:2021tnl}.
		In all simulations presented in this thesis, 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.

		\section{Interpolation of the velocity field}
		\label{sec:interpolation}
		\label{sec:lbm-mhfem:interpolation}

		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.

content/LBM.tex

+1 −0

Original line number	Diff line number	Diff line
		@@ -30,6 +30,7 @@ The detailed derivation of LBM is not within the scope of this thesis.
		Instead, we review the components needed to describe the implementation of the method and refer the reader to \cite{kruger:2017lattice} for further details.

		\subsection{Lattice and velocity set}
		\label{sec:lbm:lattice}

		We consider the discretization of a cuboidal domain $\Omega \subset \mathbb R^D$ using an equidistant orthogonal grid represented by a finite set $\mathcal G_\Omega$ of grid nodes.
		The extended grid $\mathcal G_{\overline \Omega}$ consists of the interior nodes from $\mathcal G_\Omega$ and an additional outer layer of nodes, which is added to facilitate the discretization of the domain boundary.

content/vapor_transport.tex

+0 −12

Original line number	Diff line number	Diff line
		@@ -378,14 +378,6 @@ However, water vapor is still released diffusively from the sides where the no-s

		\inline{Connect to \cref{chapter:MHFEM,chapter:LBM,chapter:LBM-MHFEM}}

		The transport equation \cref{eq:transport} has been incorporated into the mathematical framework of the \emph{NumDwarf} solver described in \cref{chapter:MHFEM}.
		Both \cref{eq:transport:conservative} and \cref{eq:transport:non-conservative} can be converted to the form of \cref{eq:NumDwarf} by taking $n=1$, $\vec Z = [\phi]$, $\matrix N = [1]$, $\vec m = [1]$, $\vec D = [D]$, $\vec w = \vec 0$, $\vec r = \vec 0$, $\vec f = [0]$, %[f_\phi]$,
		and depending on the equation:
		\begin{itemize}
		\item $\vec u = \vec 0$, $\vec a = [\vec v]$ for \cref{eq:transport:conservative},
		\item $\vec u = [\vec v]$, $\vec a = \vec 0$ for \cref{eq:transport:non-conservative}.
		\end{itemize}

		To solve the system of Eqs.~\eqref{eq:NumDwarf}, suitable initial and boundary conditions must be supplied.
		A simple Dirichlet-type boundary condition is used to prescribe fixed values of the unknown function on parts of the domain boundary $\partial \Omega_2$ as described in \cref{sec:boundary conditions}.
		On the remaining parts of $\partial \Omega_2$, which are either impermeable walls where the no-slip condition on velocity is imposed or free-stream boundaries of $\Omega_2$ inside $\Omega_1$, a Neumann-type condition is used to prescribe zero gradient of the unknown function in the normal direction.
		@@ -401,10 +393,6 @@ The boundary condition for velocity on the downstream faces of synthetic plants
		On the inflow, the discrete distribution functions are approximated by the discrete equilibrium functions evaluated from the known macroscopic variables \cite{latt2008,mohamad2009,haussmann2019}.
		On the outflow, the extrapolation outflow boundary condition is used to approximate the discrete distribution functions.

		An important choice related to the solver performance is the selection of the algorithm for the solution of large systems of linear equations.
		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.}