From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001
From: stainer_t <thomas.stainer@oecd-nea.org>
Date: Mon, 8 Sep 2025 13:48:49 +0200
Subject: Initial commit from Polytechnique Montreal

---
 doc/IGE335/Section3.04_sn.tex | 220 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 220 insertions(+)
 create mode 100644 doc/IGE335/Section3.04_sn.tex

(limited to 'doc/IGE335/Section3.04_sn.tex')

diff --git a/doc/IGE335/Section3.04_sn.tex b/doc/IGE335/Section3.04_sn.tex
new file mode 100644
index 0000000..6e77e6d
--- /dev/null
+++ b/doc/IGE335/Section3.04_sn.tex
@@ -0,0 +1,220 @@
+\subsubsection{The {\tt SNT:} tracking module}\label{sect:SNData}
+
+The {\tt SNT:} module provides an implementation of the discrete ordinates method. The {\tt SNT:} module can process one-, two- and three-dimensional regular
+geometries of type \moc{CAR1D}, \moc{TUBE}, \moc{SPHERE}, \moc{CAR2D}, \moc{TUBEZ}, \moc{CAR3D}, \moc{HEX2D} and \moc{HEX3D}.
+
+\vskip 0.2cm
+
+The calling specification for this module is:
+
+\begin{DataStructure}{Structure \dstr{SNT:}}
+\dusa{TRKNAM}
+\moc{:=} \moc{SNT:} $[$ \dusa{TRKNAM} $]$ 
+\dusa{GEONAM} \moc{::}  \dstr{desctrack} \dstr{descsn}
+\end{DataStructure}
+
+\noindent  where
+\begin{ListeDeDescription}{mmmmmmm}
+
+\item[\dusa{TRKNAM}] {\tt character*12} name of the \dds{tracking} data
+structure that will contain region volume and surface area vectors in
+addition to region identification pointers and other tracking information.
+If \dusa{TRKNAM} also appears on the RHS, the previous tracking 
+parameters will be applied by default on the current geometry.
+
+\item[\dusa{GEONAM}] {\tt character*12} name of the \dds{geometry} data
+structure.
+
+\item[\dstr{desctrack}] structure describing the general tracking data (see
+\Sect{TRKData})
+
+\item[\dstr{descsn}] structure describing the transport tracking data
+specific to \moc{SNT:}.
+
+\end{ListeDeDescription}
+
+\vskip 0.2cm
+
+The \moc{SNT:} specific tracking data in \dstr{descsn} is defined as
+
+\begin{DataStructure}{Structure \dstr{descsn}}
+$[~\{$ \moc{ONEG} $|$ \moc{ALLG} $\}~]~[$ \moc{KBA} \dusa{m} $]$ \\
+$[$ \moc{SCHM} \dusa{ischm} $]~[$ \moc{DIAM} \dusa{mm} $]$ \\
+\moc{SN} \dusa{n} $~[$ \moc{SCAT} \dusa{iscat} $]~$ \\
+$[~\{$ \moc{LIVO} \dusa{icl1} \dusa{icl2} $|$ \moc{NLIVO} $|$ \moc{GMRES} \dusa{nstart} $\}~]$ \\
+$[~\{$ \moc{DSA} \dusa{ndsa} \dusa{mdsa} \dusa{sdsa} $|$ \moc{NDSA} $\}~]$ \\
+$[$ \moc{NSHT} $]$ \\
+$[$ \moc{FOUR} \dusa{nfou} $]$ \\
+$[$ \moc{MAXI} \dusa{maxi} $]~[$ \moc{EPSI} \dusa{epsi} $]$ \\
+$[$ \moc{QUAD} \dusa{iquad} $]$ \\
+$[~\{$ \moc{BTE} $|$ \moc{BFPG} $|$ \moc{BFPL}$ \}~]$ \\
+$[$ \moc{ESCHM} \dusa{eschm} $]~[$ \moc{EDIAM} \dusa{emm} $]$ \\
+$[$ \moc{GQ} \dusa{igq} $]$ \\
+$[~[$ \moc{QUAB} \dusa{iquab} $]~[~\{$ \moc{SAPO} $|$ \moc{HEBE} $|$ \moc{SLSI} $[$ \dusa{frtm} $]~\}~]~]$ \\
+{\tt ;}
+\end{DataStructure}
+
+\noindent where
+
+\begin{ListeDeDescription}{mmmmmmm}
+
+\item[\dstr{desctrack}] structure describing the general tracking data (see
+\Sect{TRKData})
+
+\item[\moc{ONEG}] keyword to specify that the multigroup flux is computed as a sequence of one-group solutions using Gauss-Seidel iterations. This is the default option.
+
+\item[\moc{ALLG}] keyword to specify that the multigroup flux is computed in parallel for a set of energy groups, using the Jacobi method.
+
+\item[\moc{KBA}] keyword to specify that nested loops over both macrocells (in each wavefront) and energy (dependent on whether \moc{ONEG} or \moc{ALLG} is selected) are used for multithreading in 2D and 3D geometries, somewhat similar to the Koch-Baker-Alcouffe (KBA) algorithm.\cite{kba,domino}
+
+\item[\dusa{m}] for Cartesian geometries, subdivide the domain into $m\times m$ (2D) or $m \times m \times m$ (3D) macrocells for the sweep algorithm. This is irrespective of the meshing of the domain, as long as $m$ is smaller than the mesh length along each axis. For 2D hexagonal geometry, the subdivision of the domain is limited to the hexagons being macrocells; hence, any value of $m>1$ triggers this method but results in exactly the same behaviour. For 3D hexagonal geometry, the $z$-axis is divided into $m$ macrocells. 
+
+%Update to tabular format?
+\item[\moc{SCHM}] keyword to specify the spatial discretisation scheme. 
+
+\item[\dusa{ischm}] index to specify the spatial discretisation scheme. \dusa{ischm} $=1$ is used for High-Order Diamond Differencing (HODD) (default value). \dusa{ischm} $=2$ is used for the Discontinuous Galerkin finite element method (DG) currently available only in 1D slab, and 2D/3D Cartesian/hexagonal geometries. \dusa{ischm} $=3$ is used for the Adaptive Weighted Difference method (AWD), only available for Cartesian geometries.
+
+\item[\moc{DIAM}] keyword to fix the spatial approximation order.
+
+\item[\dusa{mm}] order of the Legendre polynomial expansion used in the spatial discretisation method. For HODD, \dusa{mm} $=0$ is the default, while for DG, it is \dusa{mm} $=1$.
+For Cartesian geometries, any order \dusa{mm} $\geq0$ is available. For 2D hexagonal geometries, linear and parabolic orders are available. Classical diamond difference (\dusa{mm} $=0$ with \dusa{ischm} $=1$) are available for 1D tube and 1D sphere geometries. Adaptive schemes (\dusa{ischm} $=3$) are only available with constant order.
+\begin{displaymath}
+\dusa{mm} = \left\{
+\begin{array}{rl}
+ 0 & \textrm{Constant (classical diamond scheme (HODD) or step scheme (DG))} \\
+ 1 & \textrm{Linear} \\
+ 2 & \textrm{Parabolic} \\
+ >3  & \textrm{Higher-orders}
+\end{array} \right.
+\end{displaymath}
+
+\item[\moc{SN}] keyword to fix the angular approximation order of the flux.
+
+\item[\dusa{n}] order of the $S_N$ approximation (even number).
+
+\item[\moc{SCAT}] keyword to limit the anisotropy of scattering sources.
+
+\item[\dusa{iscat}] number of terms in the scattering sources. \dusa{iscat} $=1$ is used for
+isotropic scattering in the laboratory system. \dusa{iscat} $=2$ is used for
+linearly anisotropic scattering in the laboratory system. The default value is set to $n$.
+
+\item[\moc{LIVO}] keyword to enable Livolant acceleration of the scattering iterations (default value).
+\item[\dusa{icl1},~\dusa{icl2}] Numbers of respectively free and accelerated iterations in the Livolant method.
+\item[\moc{NLIVO}] keyword to disable acceleration method and to perform free scattering iterations
+
+\item[\moc{GMRES}] keyword to set the GMRES(m) acceleration of the scattering iterations. The default value,
+equivalent to \dusa{nstart}=0, corresponds to a one-parameter Livolant acceleration.\cite{gmres}
+
+\item[\dusa{nstart}] restarts the GMRES method every \dusa{nstart} iterations.
+
+\item[\moc{DSA}] keyword to enable diffusion synthetic acceleration using BIVAC or TRIVAC.
+\item[\dusa{ndsa}] apply the synthetic acceleration every \dusa{ndsa} number of inner flux iterations. Depending on the test case, if the DSA is enabled too soon or enabled at every inner iteration, instabilities and convergence failure can occur. A value of $0$ can be set to start the DSA immediately and have the acceleration applied every inner iteration thereafter. The default is \dusa{nsdsa}~$=1000$, indicating the DSA will not be applied. Benchmarks suggests that the optimal values are $3$ and $5$ for Cartesian and hexagonal geometries respectively.
+\item[\dusa{mdsa}] order of the Raviart-Thomas spatial approximation used in the DSA resolution. Sometimes, using the same order as the transport calculation does not provide any benefit to the solution, and ends up being a drain on computational resources. Hence, there is the option of using a different order than the transport approximation.
+\begin{displaymath}
+\dusa{mdsa} = \left\{
+\begin{array}{rl}
+ 0 & \textrm{Constant} \\
+ 1 & \textrm{Linear} \\
+ 2 & \textrm{Parabolic} \\
+\end{array} \right.
+\end{displaymath}
+\item[\dusa{sdsa}] choose the solver to use for the synthetic acceleration. Note that TRIVAC generally works better and is faster with hexagonal geometries for the matrix assemblies. Also, for 3D geometries, TRIVAC \emph{has} to be chosen.
+\begin{displaymath}
+\dusa{sdsa} = \left\{
+\begin{array}{rl}
+ 1 & \textrm{BIVAC} \\
+ 2 & \textrm{TRIVAC} \\
+\end{array} \right.
+\end{displaymath}
+
+\item[\moc{NDSA}] keyword to disable diffusion synthetic acceleration (default).
+
+\item[\moc{NSHT}] keyword to disable the shooting method for 1D cases -- can be useful for debugging purposes.
+
+\item[\moc{FOUR}] keyword to pass the number of frequencies to be investigated in Fourier analysis (only works in 1D Cartesian geometry).
+\item[\dusa{nfou}] number of frequencies to be investigated in 1D Fourier analysis along the range $[0, \frac{2\pi}{L})$ where $L$ is the length of the slab.
+
+\item[\moc{MAXI}] keyword to set the maximum number of inner iterations (or GMRES iterations if actived).
+\item[\dusa{maxi}] maximum number of inner iterations. Default value: $100$.
+
+\item[\moc{EPSI}] set the convergence criterion on inner iterations (or GMRES iterations if actived).
+\item[\dusa{epsi}] convergence criterion on inner iterations. The default value is $1\times 10^{-5}$.
+\item[\moc{QUAD}] keyword to set the type of angular quadrature.
+
+\item[\dusa{iquad}] type of quadrature. For 1D quadrature (over $\mu$)
+\begin{displaymath}
+    \dusa{iquad} = \left\{
+    \begin{array}{rl}
+        1 & \textrm{Gauss-Lobatto} \\
+        2 & \textrm{Gauss-Legendre (default option in 1D)} \\
+    \end{array} \right.
+\end{displaymath}
+and for 2D/3D quadrature (over the unit sphere):
+\begin{displaymath}
+    \dusa{iquad} = \left\{
+    \begin{array}{rl}
+        1 & \textrm{Lathrop-Carlson level-symmetric quadrature} \\
+        2 & \textrm{$\mu_1$--optimi\-zed level-symmetric quadrature (default option in 2D and in 3D)} \\
+        3 & \textrm{Snow-code level-symmetric quadrature (obsolete)} \\
+        4 & \textrm{Legendre-Chebyshev quadrature (variable number of base points per axial level)} \\
+        5 & \textrm{symmetric Legendre-Chebyshev quadrature} \\
+        6 & \textrm{quadruple range (QR) quadrature\cite{quadrupole}} \\
+        10 & \textrm{Legendre-Chebyshev quadrature (equal number of base points per axial level)} \\
+    \end{array} \right.
+\end{displaymath}
+
+\item[\moc{BTE}] solution of the Boltzmann transport equation (default option).
+
+\item[\moc{BFPG}] solution of the Boltzmann Fokker-Planck equation with Galerkin energy propagation factors.
+
+\item[\moc{BFPL}] solution of the Boltzmann Fokker-Planck equation with Przybylski and Ligou energy propagation factors.\cite{ligou}
+
+\item[\moc{ESCHM}] keyword to specify the energy discretisation scheme to use for the continuous slowing-down term of the Boltzmann Fokker-Planck equation. 
+
+\item[\dusa{eschm}] index to specify the energy discretisation scheme. \dusa{ischm} $=1$ is used for High-Order Diamond Differencing (HODD) (default value). \dusa{ischm} $=2$ is
+used for the Discontinuous Galerkin finite element method (DG). \dusa{ischm} $=3$ is used for the Adaptive Weighted Difference method (AWD). All of these schemes are available
+only for Cartesian geometries.
+
+\item[\moc{EDIAM}] keyword to fix the energy approximation order.
+
+\item[\dusa{emm}] order of the Legendre polynomial expansion used in the energy discretisation method. For HODD, \dusa{emm} $=0$ is the default, while for DG,
+it is \dusa{emm} $=1$. For Cartesian geometries, any order \dusa{emm} $\geq0$ is available. Adaptive schemes (\dusa{ischm} $=3$) are only available with
+constant order.
+\begin{displaymath}
+\dusa{emm} = \left\{
+\begin{array}{rl}
+ 0 & \textrm{Constant (classical diamond scheme (HODD) or step scheme (DG))} \\
+ 1 & \textrm{Linear} \\
+ 2 & \textrm{Parabolic} \\
+ >3  & \textrm{Higher-orders}
+\end{array} \right.
+\end{displaymath}
+
+\item[\moc{GQ}] keyword to specify the angular interpolation method.
+
+\item[\dusa{igq}] type of angular interpolation. =0: classical SN method, with all spherical harmonics up to the given Legendre truncation order; =1: Galerkin
+quadrature method, where the moment-to-discrete matrix ($\shadowM$) is given by the inverse of the discrete-to-moment ($\shadowD$) matrix; =2: Galerkin quadrature
+method, where the discrete-to-moment ($\shadowD$) matrix is given by the inverse of the moment-to-discrete ($\shadowM$) matrix \cite{morel1989,shands2024}.
+
+\item[\moc{QUAB}] keyword to specify the number of basis point for the
+numerical integration of each micro-structure in cases involving double
+heterogeneity (Bihet).
+
+\item[\dusa{iquab}] the number of basis point for the numerical integration of
+the collision probabilities in the micro-volumes using the  Gauss-Jacobi
+formula. The values permitted are: 1 to 20, 24, 28, 32 or  64. The default value
+is \dusa{iquab}=5. If \dusa{iquab} is negative, its absolute value will be used in the She-Liu-Shi approach to determine the
+split level in the tracking used to compute the probability collisions.
+
+\item[\moc{SAPO}] use the Sanchez-Pomraning double-heterogeneity model.\cite{sapo}
+
+\item[\moc{HEBE}] use the Hebert double-heterogeneity model (default option).\cite{BIHET}
+
+\item[\moc{SLSI}] use the She-Liu-Shi double-heterogeneity model without shadow effect.\cite{She2017}
+
+\item[\dusa{frtm}] the minimum microstructure volume fraction used to compute the size of the equivalent cylinder in She-Liu-Shi approach. The default value is
+\dusa{frtm} $=0.05$.
+
+\end{ListeDeDescription}
+
+\eject
-- 
cgit v1.2.3

