Permalink
Browse files

FDS Source, V&V, Manuals : Rename GLMAT to UGLMAT, GLMAT IBM to GLMAT…

… for consistency with ScaRC.
  • Loading branch information...
marcosvanella committed Jun 27, 2018
1 parent 31e6dc2 commit 1c94956d5605d716cbc3a0525652b4272d04bbbd
@@ -1629,17 +1629,17 @@ \section{Accuracy of the Pressure Solver}
\section{Optional Pressure Solver}
\label{optional_pressure_solver}
The default linear solver in FDS is the CrayFishpak FFT solver ({\ct SOLVER='FFT'}). For non-overlapping, non-stretched meshes at the same refinement level covering a single connected domain (single pressure zone), the solver may be changed to {\ct SOLVER='GLMAT'} on {\ct PRES}. GLMAT enables the solution of the Poisson equation from the global discretization matrix. The solver used in this mode is the Intel\textsuperscript{\textregistered} MKL Sparse Cluster Solver. When using {\ct SOLVER='GLMAT'}, the unknown values of the pseudo-pressure, $H$, live only in gas-phase cells. This allows for the normal components of velocity of a solid surface to be computed exactly (no velocity penetration error). Strictly speaking, in this mode, FDS is no longer using an ``immersed boundary method'' for flow obstructions---the pressure solve is now ``unstructured.'' If the solver is defined as {\ct SOLVER='GLMAT IBM'}, the pressure is solved in both solid and gas cells using an immersed boundary method for flow obstructions, as usual for FDS. This mode will reproduce the exact same pressure solution as a single mesh FFT solve. That is, velocity errors at any mesh boundaries will be eliminated. As GLMAT is based on Intel\textsuperscript{\textregistered} MKL parallel LU decomposition solver, it requires computation and storage of global factorization matrices. Depending on the size of the problem, these operations can be taxing on computing resources. Our tests find it suitable for problems with up to a million cells on a multicore workstation with 64 GB of RAM memory. Note that currently a single discretization matrix is built for all gas cells defined on the model, therefore the solver is not meant to be used in cases of multiple sealed compartments and pressure zones.
The default linear solver in FDS is the CrayFishpak FFT solver ({\ct SOLVER='FFT'}). For non-overlapping, non-stretched meshes at the same refinement level covering a single connected domain (single pressure zone), the solver may be changed to {\ct SOLVER='UGLMAT'} on {\ct PRES}. UGLMAT enables the solution of the Poisson equation from the global discretization matrix. The solver used in this mode is the Intel\textsuperscript{\textregistered} MKL Sparse Cluster Solver. When using {\ct SOLVER='UGLMAT'}, the unknown values of the pseudo-pressure, $H$, live only in gas-phase cells. This allows for the normal components of velocity of a solid surface to be computed exactly (no velocity penetration error). Strictly speaking, in this mode, FDS is no longer using an ``immersed boundary method'' for flow obstructions---the pressure solution is now ``unstructured''. If the solver is defined as {\ct SOLVER='GLMAT'}, the pressure is solved in both solid and gas cells using an immersed boundary method for flow obstructions, as usual for FDS. This mode will reproduce the exact same pressure solution as a single mesh FFT solve. That is, velocity errors at any mesh boundaries will be eliminated. As UGLMAT and GLMAT are based on Intel\textsuperscript{\textregistered} MKL parallel LU decomposition solver, they require computation and storage of global factorization matrices. Depending on the size of the problem, these operations can be taxing on computing resources. Our tests find it suitable for problems with up to a million cells on a multicore workstation with 64 GB of RAM memory. Note that currently a single discretization matrix is built for all gas cells defined on the model, therefore the solver is not meant to be used in cases of multiple sealed compartments and pressure zones.
\subsubsection{Example Case: Pressure\_Solver/duct\_flow}
\label{duct_flow}
To demonstrate how to improve the accuracy of the pressure solver, consider the flow of air through a square duct that crosses several meshes. In the sample input file, {\ct duct\_flow.fds}, air is pushed through a 1 m$^2$ duct at 1~m/s. With no thermal expansion, the volume flow into the duct ought to equal the volume flow out of the duct. Figure~\ref{duct_flow_fig} displays the computed inflow and outflow as a function of time and the number of pressure iterations required. The default pressure solution strategy uses block-wise direct FFT solves combined with a direct forcing immersed boundary method to model flow obstructions. Hence, there are two sources of pressure error: one at a mesh interface and one at a solid boundary. The default pressure iteration scheme tries to drive both sources of error to the specified (or default) velocity tolerance. In the {\ct duct\_flow} case, the {\ct VELOCITY\_TOLERANCE} has been set to 0.001~m/s with {\ct MAX\_PRESSURE\_ITERATIONS} set to 1000 and the grid cell size is 0.2~m. For the default scheme (FFT IBM), the outflow does not match the inflow exactly because of inaccuracies at the solid and mesh boundaries. We have also included the results using the GLMAT solver with an unstructured pressure matrix, which allows the flow solver to enforce the impermeability condition (zero normal component of velocity) at a solid boundary.
To demonstrate how to improve the accuracy of the pressure solver, consider the flow of air through a square duct that crosses several meshes. In the sample input file, {\ct duct\_flow.fds}, air is pushed through a 1 m$^2$ duct at 1~m/s. With no thermal expansion, the volume flow into the duct ought to equal the volume flow out of the duct. Figure~\ref{duct_flow_fig} displays the computed inflow and outflow as a function of time and the number of pressure iterations required. The default pressure solution strategy uses block-wise direct FFT solves combined with a direct forcing immersed boundary method to model flow obstructions. Hence, there are two sources of pressure error: one at a mesh interface and one at a solid boundary. The default pressure iteration scheme tries to drive both sources of error to the specified (or default) velocity tolerance. In the {\ct duct\_flow} case, the {\ct VELOCITY\_TOLERANCE} has been set to 0.001~m/s with {\ct MAX\_PRESSURE\_ITERATIONS} set to 1000 and the grid cell size is 0.2~m. For the default scheme (FFT + IBM), the outflow does not match the inflow exactly because of inaccuracies at the solid and mesh boundaries. We have also included the results using the UGLMAT solver with an unstructured pressure matrix, which allows the flow solver to enforce the impermeability condition (zero normal component of velocity) at a solid boundary.
\begin{lstlisting}
&PRES SOLVER='GLMAT' /
&PRES SOLVER='UGLMAT' /
\end{lstlisting}
As seen in Fig.~\ref{duct_flow_fig}, the volume flow matches perfectly and the number of required pressure iterations is one. Whether this pressure solver strategy is most efficient likely depends on the problem and the required error tolerance. In this particular case, the GLMAT case runs about 25 \% faster than the default FFT case, and it is more accurate---easy choice.
As seen in Fig.~\ref{duct_flow_fig}, the volume flow matches perfectly and the number of required pressure iterations is one. Whether this pressure solver strategy is most efficient likely depends on the problem and the required error tolerance. In this particular case, the UGLMAT case runs about 25 \% faster than the default FFT case, and it is more accurate---easy choice.
\begin{figure}[ht]
\begin{center}
@@ -1654,20 +1654,20 @@ \subsubsection{Example Case: Pressure\_Solver/duct\_flow}
\subsubsection{Example Case: Pressure\_Solver/dancing\_eddies}
\label{dancing_eddies}
In this example, air is pushed through a 30~cm long, two-dimensional channel at 0.5~m/s. A plate obstruction normal to the flow creates a Karman vortex street. The computational domain is divided into 4 meshes. Three simulations are performed: one in which the {\ct VELOCITY\_TOLERANCE} is set to the default value, one in which it is set to a relatively small value, and one in which the {\ct SOLVER='GLMAT'} is used. Figure~\ref{dancing_eddies_figure} shows the downstream pressure histories for these cases compared to a simulation that uses only one mesh. The case with the tighter tolerance produces a better match to the single mesh solution, but at a higher computational cost---about a factor of 2.5 . The {\ct SOLVER='GLMAT'} case takes just a little longer than the default---by a factor of 1.25---but the error is machine precision. These cases used MPI for domain decomposition, but only a single OpenMP thread. Note, however, that the GLMAT solver is threaded, whereas the FFT solver (FDS default) is not.
In this example, air is pushed through a 30~cm long, two-dimensional channel at 0.5~m/s. A plate obstruction normal to the flow creates a Karman vortex street. The computational domain is divided into 4 meshes. Three simulations are performed: one in which the {\ct VELOCITY\_TOLERANCE} is set to the default value, one in which it is set to a relatively small value, and one in which the {\ct SOLVER='UGLMAT'} is used. Figure~\ref{dancing_eddies_figure} shows the downstream pressure histories for these cases compared to a simulation that uses only one mesh. The case with the tighter tolerance produces a better match to the single mesh solution, but at a higher computational cost---about a factor of 2.5 . The {\ct SOLVER='UGLMAT'} case takes just a little longer than the default---by a factor of 1.25---but the error is machine precision. These cases used MPI for domain decomposition, but only a single OpenMP thread. Note, however, that the UGLMAT solver is threaded, whereas currently the FFT solver (FDS default) is not.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=3in]{SCRIPT_FIGURES/dancing_eddies_default}
\includegraphics[width=3in]{SCRIPT_FIGURES/dancing_eddies_tight}
\includegraphics[width=3in]{SCRIPT_FIGURES/dancing_eddies_glmat}
\includegraphics[width=3in]{SCRIPT_FIGURES/dancing_eddies_uglmat}
\includegraphics[width=\textwidth]{SCRIPT_FIGURES/dancing_eddies}
\end{center}
\caption[Results of the {\ct dancing\_eddies} test cases]{(Top) Comparison of pressure traces in the channel for three different settings of {\ct VELOCITY\_TOLERANCE}, the default value (upper-left), a tighter tolerance (upper-right), machine precision from the GLMAT solver with a single iteration (middle). (Bottom) A contour plot of the pressure after 2~s with the default tolerance.}
\caption[Results of the {\ct dancing\_eddies} test cases]{(Top) Comparison of pressure traces in the channel for three different settings of {\ct VELOCITY\_TOLERANCE}, the default value (upper-left), a tighter tolerance (upper-right), machine precision from the UGLMAT solver with a single iteration (middle). (Bottom) A contour plot of the pressure after 2~s with the default tolerance.}
\label{dancing_eddies_figure}
\end{figure}
One strategy for reducing the computational cost associated with a tightened tolerance on the normal velocity at mesh boundaries is to overlap the meshes. The idea is described in detail for a non-structured mesh in Ref.~\cite{Grinberg:JCP2010}. The basic idea is that a slight overlap of meshes can lead to faster convergence of the pressure solver towards the desired level of accuracy. In the {\ct dancing\_eddies} test case, the case called {\ct dancing\_eddies\_tight} is rerun as {\ct dancing\_eddies\_tight\_overlap} where each mesh is extended an additional 10 grid cells in the longitudinal direction. This increases the work by approximately a factor of 85/75, but it reduces the number of iterations. Figure~\ref{dancing_eddies_overlap} shows the reduction in the number of pressure iterations (left) and the CPU time (right). It can be seen that the overlap strategy does not ultimately reduce the CPU cost, whereas the GLMAT solver reduces the cost by a factor of 0.5 by using a single iteration as shown on the left plot.
One strategy for reducing the computational cost associated with a tightened tolerance on the normal velocity at mesh boundaries is to overlap the meshes. The idea is described in detail for a non-structured mesh in Ref.~\cite{Grinberg:JCP2010}. The basic idea is that a slight overlap of meshes can lead to faster convergence of the pressure solver towards the desired level of accuracy. In the {\ct dancing\_eddies} test case, the case called {\ct dancing\_eddies\_tight} is rerun as {\ct dancing\_eddies\_tight\_overlap} where each mesh is extended an additional 10 grid cells in the longitudinal direction. This increases the work by approximately a factor of 85/75, but it reduces the number of iterations. Figure~\ref{dancing_eddies_overlap} shows the reduction in the number of pressure iterations (left) and the CPU time (right). It can be seen that the overlap strategy does not ultimately reduce the CPU cost, whereas the UGLMAT solver reduces the cost by a factor of 0.5 by using a single iteration as shown on the left plot.
\begin{figure}[ht]
\includegraphics[width=3in]{SCRIPT_FIGURES/dancing_eddies_iterations}
@@ -7582,13 +7582,13 @@ SUBROUTINE READ_PRES
IF (SCARC_METHOD == 'NONE') SCARC_METHOD = 'KRYLOV' ! Use Krylov default solver for SCARC
IF (SCARC_PRECON == 'NONE') SCARC_PRECON = 'FFT' ! Use FFT as default preconditioner for SCARC
CASE('GLMAT')
CASE('UGLMAT')
PRES_METHOD = 'GLMAT'
GLMAT_SOLVER = .TRUE.
PRES_ON_WHOLE_DOMAIN = .FALSE.
IF (CHECK_POISSON) GLMAT_VERBOSE=.TRUE.
CASE('GLMAT IBM')
CASE('GLMAT')
PRES_METHOD = 'GLMAT'
GLMAT_SOLVER = .TRUE.
PRES_ON_WHOLE_DOMAIN = .TRUE.
Oops, something went wrong.

0 comments on commit 1c94956

Please sign in to comment.