11 files changed, 340 insertions, 38 deletions

diff --git a/doc/analytic-k2.png b/doc/analytic-k2.png
new file mode 100644
index 0000000..6389688
--- /dev/null
+++ b/doc/analytic-k2.png
diff --git a/doc/analytic.png b/doc/analytic.png
new file mode 100644
index 0000000..83078d0
--- /dev/null
+++ b/doc/analytic.png
diff --git a/doc/excerpt.cc b/doc/excerpt.cc
new file mode 100644
index 0000000..5a57a0d
--- /dev/null
+++ b/doc/excerpt.cc
@@ -0,0 +1,19 @@
+for (int i = 0; i < ir->GetNPoints(); ++i) {
+        // Get shape functions, coefficients at this integration point.
+        ...
+
+        // Calculate phi_grad.
+        phi_grad = 0.0;
+        for (int b = 0; b < nen; ++b) {
+                for (int j = 0; j < nsd; ++j) {
+                        phi_grad(j) += dshape(b, j) * elfun(b);
+                }
+        }
+
+        // Calculate SUPG term.
+        for (int a = 0; a < nen; ++a) {
+                dshape.GetRow(a, Na_grad);
+                double tau_tmp = (A_ip*Na_grad)*tau_ip;
+                elvec(a) += tau_tmp * (A_ip * phi_grad) * wdetj;
+        }
+}
diff --git a/doc/excerpt2.cc b/doc/excerpt2.cc
new file mode 100644
index 0000000..765e338
--- /dev/null
+++ b/doc/excerpt2.cc
@@ -0,0 +1,14 @@
+for (int i = 0; i < ir->GetNPoints(); ++i) {
+        // Get shape functions, coefficients at this integration point.
+        ...
+
+        // Calculate tangent matrix.
+        for (int a = 0; a < nen; ++a) {
+                dshape.GetRow(a, Na_grad);
+                double tau_tmp = (A_ip*Na_grad)*tau_ip;
+                for (int b = 0; b < nen; ++b) {
+                        dshape.GetRow(b, Nb_grad);
+                        elmat(a, b) += tau_tmp * (A_ip * Nb_grad) * wdetj;
+                }
+        }
+}
diff --git a/doc/line-sol-k2.png b/doc/line-sol-k2.png
new file mode 100644
index 0000000..aba807f
--- /dev/null
+++ b/doc/line-sol-k2.png
diff --git a/doc/line-supg-sol-k2.png b/doc/line-supg-sol-k2.png
new file mode 100644
index 0000000..1acd732
--- /dev/null
+++ b/doc/line-supg-sol-k2.png
diff --git a/doc/paper-excerpt1.cc b/doc/paper-excerpt1.cc
new file mode 100644
index 0000000..d34bce3
--- /dev/null
+++ b/doc/paper-excerpt1.cc
@@ -0,0 +1,24 @@
+for (int i = 0; i < ir->GetNPoints(); ++i) {
+        // Get shape functions, coefficients at this integration point.
+        const IntegrationPoint& ip = ir->IntPoint(i);
+        Tr.SetIntPoint(&ip);
+        el.CalcPhysDShape(Tr, dshape);
+        A_ir.GetColumn(i, A_ip);
+        double tau_ip = tau->Eval(Tr, ip);
+        const double wdetj = ip.weight * Tr.Weight();
+
+        // Calculate phi_grad.
+        phi_grad = 0.0;
+        for (int b = 0; b < nen; ++b) {
+                for (int j = 0; j < nsd; ++j) {
+                        phi_grad(j) += dshape(b, j) * elfun(b);
+                }
+        }
+
+        // Calculate SUPG term.
+        for (int a = 0; a < nen; ++a) {
+                dshape.GetRow(a, Na_grad);
+                double tau_tmp = (A_ip*Na_grad)*tau_ip;
+                elvec(a) += tau_tmp * (A_ip * phi_grad) * wdetj;
+        }
+}
diff --git a/doc/paper-excerpt2.cc b/doc/paper-excerpt2.cc
new file mode 100644
index 0000000..5843c4d
--- /dev/null
+++ b/doc/paper-excerpt2.cc
@@ -0,0 +1,19 @@
+for (int i = 0; i < ir->GetNPoints(); ++i) {
+        // Get shape functions, coefficients at this integration point.
+        const IntegrationPoint& ip = ir->IntPoint(i);
+        Tr.SetIntPoint(&ip);
+        el.CalcPhysDShape(Tr, dshape);
+        A_ir.GetColumn(i, A_ip);
+        double tau_ip = tau->Eval(Tr, ip);
+        double wdetj = ip.weight * Tr.Weight();
+
+        // Calculate tangent matrix.
+        for (int a = 0; a < nen; ++a) {
+                dshape.GetRow(a, Na_grad);
+                double tau_tmp = (A_ip*Na_grad)*tau_ip;
+                for (int b = 0; b < nen; ++b) {
+                        dshape.GetRow(b, Nb_grad);
+                        elmat(a, b) += tau_tmp * (A_ip * Nb_grad) * wdetj;
+                }
+        }
+}
diff --git a/doc/report.pdf b/doc/report.pdf
index c45cdcc..e9709a0 100644
--- a/doc/report.pdf
+++ b/doc/report.pdf
diff --git a/doc/report.tex b/doc/report.tex
index 9c72b59..510637f 100644
--- a/doc/report.tex
+++ b/doc/report.tex
@@ -5,13 +5,15 @@
 \usepackage{amsmath}
 \usepackage{graphicx}
 \usepackage{hyperref}
+\usepackage{tabularx}
+\usepackage{listings}
 
 \usepackage[backend=biber]{biblatex}
 \addbibresource{sources.bib}
 
-\title{%
-Term Project Report \\
-MANE 6660: Finite Element Method
+\title{
+        SUPG Stabilization for Advection-Diffusion with MFEM \\
+        \large{MANE 6660: Finite Element Method Final Project Report}
 }
 \author{Aiden Woodruff}
 \date{\today}
@@ -19,67 +21,102 @@ MANE 6660: Finite Element Method
 \begin{document}
 \maketitle
 
-\section*{MFEM SUPG miniapp}
-My project proposal for the FEM term project is a SUPG miniapp implemented
-with MFEM. This is an independent project related to my graduate research
-(which has focused on shock processing).
+\section{Introduction}
 
-\subsection*{Advection-diffusion equation}
-The complete strong form of the advection-diffusion equation is given:
+The advection-diffusion equation, also known as the
+convection-diffusion equation, is a PDE that combines advection and diffusion
+to describe fluid flow:
 
-\begin{gather}
+\begin{gather} \label{eq:pde}
         \frac{\delta\phi}{\delta t} + \nabla\cdot\pmb F(\phi) = s \\
-        \pmb F(\phi) = \pmb F^\text{adv} + \pmb F^\text{diff} =
+        \pmb F(\phi) = \pmb F^\text{adv}(\phi) + \pmb F^\text{diff}(\phi) =
         \pmb a\phi + (-\kappa\nabla\phi)
 \end{gather}
 
-Where \(\phi\) is the solution variable (i.e. pressure or temperature),
-\(\pmb a\) is the advection velocity vector, \(\kappa\) is diffusivity, and
-\(s\) is the source term. I make the following simplifying assumptions:
+The solution variable \(\phi\) may be a number of quantities such as pressure,
+temperature, or species concentration. \(\pmb a\) is the advection velocity,
+\(\kappa\) is diffusivity, and \(s\) is the source term. In this project, a
+simplified form of the equation is solved with the following assumptions:
 
 \begin{itemize}
         \item No time dependent term.
-        \item Constant \(\pmb a\) and \(\kappa\).
-        \item No source term; \(s=0\).
+        \item Constant \(\pmb a\) and \(\kappa\) (w.r.t. \(\phi\)).
+        \item No source term, i.e. \(s=0\).
         \item Only Dirichlet boundary conditions.
         \item Linear finite elements.
 \end{itemize}
 
-Then the simplified strong form residual follows:
+The resulting strong form residual is:
 
-\begin{equation}
+\begin{equation} \label{eq:strong}
         R(\phi) = \nabla\cdot\left(\pmb a\phi - \kappa\nabla\phi\right) = 0
 \end{equation}
 
-Using index notation:
+The strong form residual \eqref{eq:strong} is expressed in index notation as
 
-\begin{equation}
+\begin{equation} \label{eq:strong-idx}
         R(\phi) = a_i \phi,_i - \left(\kappa\phi,_i\right),_i
 \end{equation}
 
-Apply the method of weighted residuals and integration by parts to obtain the
-weak form:
+The solution can be obtained analytically for the 1D case using separation of
+variables:
 
 \begin{equation}
-        \int_\Omega -w,_i a_i \phi + w,_i \kappa \phi,_i d\Omega = 0
+        \phi = \frac{\kappa}{a}e^{\frac{a}{k}x + c_1} + c_2
+\end{equation}
+
+Where \(c_1\) and \(c_2\) are solved based on the boundary conditions
+\((x_L,\phi_L)\) and \((x_R,\phi_R)\):
+
+\begin{equation}
+        c_1 = \ln\left(
+                \frac{a}{\kappa} \cdot
+                \frac{\phi_R - \phi_L}{e^{\frac{a}{\kappa}x_R} - e^{\frac{a}{\kappa}x_L}}
+        \right)
+        \quad
+        c_2 = \phi_L - \frac{\phi_R - \phi_L}{
+                e^{\frac{a}{\kappa}x_R} - e^{\frac{a}{\kappa}x_L}
+        }
 \end{equation}
 
-\subsection*{Stabilization}
+An analytical solution is shown in Figure \ref{fig:analytic}.
+
+\begin{figure}
+        \centering
+        \includegraphics[width=0.45\linewidth]{analytic-k1.png}
+        \caption{
+                Analytic solution with \(a_x = 1\), \(\kappa = 0.1\), and boundary
+                conditions \(\phi(0) = 0.2; \phi(1.6) = 1.4\)
+        }
+        \label{fig:analytic}
+\end{figure}
+
+The finite element method is typically used to solve cases in higher dimensions
+with general domains. The weak form is obtained from \eqref{eq:strong} by
+applying the method of weighted residuals and integration by parts:
+
+\begin{equation}
+        \int_\Omega -w,_i a_i \phi + w,_i \kappa \phi,_i d\Omega = 0
+\end{equation}
 
-Using the standard Galerkin approach in advection dominated cases results in an
-ill-conditioned matrix. See Figure \ref{fig:fem-unstable} for an example of an
-unstable advection-dominated solution.
+Then linear finite elements can be used. However, using the standard Galerkin
+approach in advection dominated cases results in an ill-conditioned matrix and
+numerical instability. See Figure \ref{fig:fem-unstable} for a comparison
+of a stable FEM solution and an unstable advection-dominated FEM solution.
 
 \begin{figure}
         \centering
-        \includegraphics[width=0.45\textwidth]{line-sol-k1e0.png}
-        \includegraphics[width=0.45\textwidth]{line-sol-k1e-3.png}
-        \caption{Stable solution (kappa=1) and unstable solution (kappa=1e-3)}
+        \includegraphics[width=0.45\textwidth]{lecsoln-k1.png}
+        \includegraphics[width=0.45\textwidth]{lecsoln-k2.png}
+        \caption{
+                Stable solution (kappa=0.1) and unstable solution (kappa=0.01) from a
+                reference FEM implementation with exact solution shown
+        }
         \label{fig:fem-unstable}
 \end{figure}
 
 One approach to remedy this instability is to introduce a stabilization term
-(artifical numerical diffusion). A simple formulation is called Streamlined-
+(artificial numerical diffusion). A simple formulation is called Streamlined-
 Upwind Petrov-Galerkin (SUPG). Here, an extra term is added to the weak form:
 \(a({\pmb F}^\text{adv}(w), -\tau R(\phi))\). In the linear finite element
 case, the diffusive term of the strong form residual can be dropped if care is
@@ -96,23 +133,180 @@ form is:
 The stabilizing term depends on the residual so that it goes to zero when the
 solution is achieved. The parameter \(\tau\) is of special importance, as too
 small or too large \(\tau\) can result in an ineffective method. Shakib and
-Hughes \cite{shakib_1991} identify an optimal value based on the cell Peclet
-number. However, this value is computationally inefficient, so an approximation
-is also available:
+Hughes \cite{shakib_1991} identify an optimal value in 1D based on the cell
+P\'{e}clet number.
+
+\begin{equation}
+        \tau_{skb} = \frac{h}{2|a_x|} \cdot \left(\coth(Pe) - \frac{q}{Pe}\right)
+        \quad
+        Pe = \frac{|a_x|h}{2\kappa}
+\end{equation}
+
+However, this value is computationally inefficient, so an approximation
+is also made available:
 
 \begin{equation}
-        \tau_{skb} = \left[
+        \tau_{skb,alg} = \left[
                 \left(\frac{2|a_x|}{h}\right)^2
-                +\left(\frac{4\kappa}{h^2}\right)^2
+                +9\left(\frac{4\kappa}{h^2}\right)^2
         \right]^{-1/2}
 \end{equation}
 
-\begin{figure}
-        \includegraphics[width=\textwidth]{line-supg-sol-k1e-3.png}
+An FEM solution with SUPG stabilization is shown in Figure \ref{fig:fem-supg}.
+
+\begin{figure}[htp]
+        \centering
+        \includegraphics[width=0.45\textwidth]{lecsoln-k2-supg.png}
         \caption{SUPG stabilized solution (kappa=1e-3)}
         \label{fig:fem-supg}
 \end{figure}
 
+In this project I use MFEM, a C++ framework for finite elements
+\cite{mfem-web}, to solve the advection-diffusion equation in an
+advection-dominated case. I compare my results to the exact solution and the
+reference FEM solution.
+
+\section{Methods}
+
+The MFEM library allows specification of the weak form using pre-defined
+\texttt{Integrator}s. Two such classes are readily available for solving
+the advection-diffusion equation: \texttt{ConservativeConvectionIntegrator} and
+\texttt{DiffusionIntegrator}. Implementing the SUPG term requires the addition
+of a custom integrator, \texttt{SUPGIntegrator}. Additionally, because the SUPG
+term includes the solution variable, it is implemented by extending
+\texttt{NonlinearFormIntegrator} instead of \texttt{BilinearFormIntegrator},
+which changes the problem setup for the stabilized case. See Listing
+\ref{inp:ex1} for the code for the nonlinear residual and Listing \ref{inp:ex2}
+for code for the nonlinear jacobian.
+
+\lstinputlisting[
+        basicstyle=\ttfamily\small,tabsize=2,frame=single,
+        caption={Key lines of code for the SUPG nonlinear residual},
+        label={inp:ex1}
+]{paper-excerpt1.cc}
+
+\lstinputlisting[
+        basicstyle=\ttfamily\small,tabsize=2,frame=single,
+        caption={Key lines of code for the SUPG nonlinear jacobian},
+        label={inp:ex2}
+]{paper-excerpt2.cc}
+
+\section{Results}
+
+The custom \texttt{SUPGIntegrator} provided stabilization and accurate results
+for advection-dominated cases, except when the left boundary condition was
+enforced -- note that \ref{tbl:k2} and the right of Figure \ref{fig:mfem-supg}
+start from \(\phi = 0.0\) whereas the left and \ref{tbl:k0} start from \(\phi =
+0.2\). The cause for this is unclear, as the implementation uses physical
+gradients and includes the inverse jacobian determinant.
+
+\begin{figure}[htp]
+        \centering
+        \includegraphics[width=0.45\textwidth]{line-sol-k0.png}
+        \includegraphics[width=0.45\textwidth]{line-supg-sol-k2.png}
+        \caption{
+                MFEM solution for \(a_x=1\), \(\kappa=1\) (left) and \(a_x=1\),
+                \(\kappa=0.01\) (right)
+        }
+        \label{fig:mfem-supg}
+\end{figure}
+
+\begin{table}
+        \centering
+        \begin{tabular}{c|c|c}
+                \small
+                Exact solution & Reference FEM solution & MFEM solution \\
+                \hline
+                0.2 & 0.2 & 0.2 \\
+                \hline
+                0.23192615 & 0.23192615 & 0.2318994842537078 \\
+                \hline
+                0.26721 & 0.26721 & 0.2671570027631416 \\
+                \hline
+                0.30620469 & 0.30620469 & 0.3061264579767903 \\
+                \hline
+                0.34930048 & 0.34930048 & 0.3491981952038976 \\
+                \hline
+                0.3969287 & 0.3969287 & 0.3968037571864139 \\
+                \hline
+                0.44956602 & 0.44956602 & 0.4494211044369776 \\
+                \hline
+                0.50773925 & 0.50773926 & 0.5075777552014797 \\
+                \hline
+                0.57203062 & 0.57203063 & 0.5718563242761415 \\
+                \hline
+                0.64308357 & 0.64308358 & 0.6429009576212408 \\
+                \hline
+                0.72160923 & 0.72160923 & 0.7214237583971538 \\
+                \hline
+                0.8083935 & 0.8083935 & 0.8082118551597843 \\
+                \hline
+                0.90430495 & 0.90430496 & 0.904135215032668 \\
+                \hline
+                1.0103035 & 1.0103035  & 1.010155960623283 \\
+                \hline
+                1.12745001 & 1.12745001 & 1.12733665457002 \\
+                \hline
+                1.25691693 & 1.25691693 & 1.256851800816963 \\
+                \hline
+                1.4 & 1.4 & 1.4
+        \end{tabular}
+        \caption{Comparison of nodal values for \(a_x=1\), \(\kappa=1\)}
+        \label{tbl:k0}
+\end{table}
+
+\begin{table}
+        \centering
+        \begin{tabular}{c|c|c}
+                \small
+                Exact solution & Reference FEM solution & MFEM solution \\
+                \hline
+                7.07789169e-85 & 0.00000000e+00 & 0 \\
+                \hline
+                1.00446783e-65 & 6.72915553e-24 & 1.344027471478775e-21 \\
+                \hline
+                2.21258808e-61 & 2.47544695e-22 & -3.440531560398167e-19 \\
+                \hline
+                4.87354958e-57 & 8.86558356e-21 & -5.957125219217317e-19 \\
+                \hline
+                1.07347073e-52 & 3.17278380e-19 & 5.144414882568559e-17 \\
+                \hline
+                2.36447663e-48 & 1.13544128e-17 & 1.095961397199189e-16 \\
+                \hline
+                5.20810637e-44 & 4.06339111e-16 & -5.302401810502381e-15 \\
+                \hline
+                1.14716177e-39 & 1.45416125e-14 & 4.901903445752921e-15 \\
+                \hline
+                2.52679194e-35 & 5.20399065e-13 & 1.007010439154497e-12 \\
+                \hline
+                5.56562963e-31 & 1.86234634e-11 & 1.920121752194572e-11 \\
+                \hline
+                1.22591151e-26 & 6.66475808e-10 & 6.356196583077191e-10 \\
+                \hline
+                2.70024979e-22 & 2.38510955e-08 & 2.383479327580857e-08 \\
+                \hline
+                5.94769596e-18 & 8.53556496e-07 & 8.549283844700147e-07 \\
+                \hline
+                1.31006722e-13 & 3.05461312e-05 & 3.054579043392884e-05 \\
+                \hline
+                2.88561507e-09 & 1.09315099e-03 & 0.001093115172881016 \\
+                \hline
+                6.35599017e-05 & 3.91204728e-02 & 0.03912052192397422 \\
+                \hline
+                1.40000000e+00 & 1.40000000e+00 & 1.4
+        \end{tabular}
+        \caption{Comparison of nodal values for \(a_x=1\), \(\kappa=0.01\)}
+        \label{tbl:k2}
+\end{table}
+
+\section{Lessons learned}
+
+Using MFEM turned out to be the most difficult part of the project. Initially,
+2D cases were also planned for, but had to be abandoned due to time
+constraints. It was unclear that a \texttt{NonlinearFormIntegrator} was
+required, and since the needs of this project differed from the provided
+examples, it was somewhat unclear what to do.
+
 Code is available at \url{https://git.aidenw.net/rpi/mfem-supg}.
 
 \printbibliography
diff --git a/doc/sources.bib b/doc/sources.bib
index 84d577c..58efee9 100644
--- a/doc/sources.bib
+++ b/doc/sources.bib
@@ -1,3 +1,35 @@
+@article{mfem-2024,
+  title = {High-Performance Finite Elements with {MFEM}},
+  author = {J. Andrej and N. Atallah and J.-P. Bäcker and  J.-S. Camier and
+            D. Copeland and V. Dobrev and Y. Dudouit and T. Duswald and
+            B. Keith and D. Kim and T. Kolev and B. Lazarov and  K. Mittal and
+            W. Pazner and S. Petrides and S. Shiraiwa and M. Stowell and V. Tomov},
+  journal={The International Journal of High Performance Computing Applications},
+  volume={38},
+  number={5},
+  pages={447-467},
+  year={2024},
+  publisher={SAGE Publications Sage UK: London, England}
+}
+
+@article{mfem-2021,
+  title = {{MFEM}: A Modular Finite Element Methods Library},
+  author = {R. Anderson and J. Andrej and A. Barker and J. Bramwell and J.-S. Camier and
+            J. Cerveny and V. Dobrev and Y. Dudouit and A. Fisher and Tz. Kolev and W. Pazner and
+            M. Stowell and V. Tomov and I. Akkerman and J. Dahm and D. Medina and S. Zampini},
+  journal = {Computers \& Mathematics with Applications},
+  doi = {10.1016/j.camwa.2020.06.009},
+  volume = {81},
+  pages = {42-74},
+  year = {2021}
+}
+
+@misc{mfem-web,
+  key = {mfem},
+  title = {{MFEM}: Modular Finite Element Methods {[Software]}},
+  howpublished = {\url{mfem.org}},
+  doi = {10.11578/dc.20171025.1248}
+}
 
 @article{shakib_1991,
         title = {A new finite element formulation for computational fluid dynamics: X. The compressible Euler and Navier-Stokes equations},