diff --git a/docs/derivations/laneEmdenVariationalBlockForm.aux b/docs/derivations/laneEmdenVariationalBlockForm.aux new file mode 100644 index 0000000..5dcb826 --- /dev/null +++ b/docs/derivations/laneEmdenVariationalBlockForm.aux @@ -0,0 +1,11 @@ +\relax +\@writefile{toc}{\contentsline {section}{\numberline {1}Continuous Variational Form}{1}{}\protected@file@percent } +\@writefile{toc}{\contentsline {section}{\numberline {2}Discritized Variational Form}{1}{}\protected@file@percent } +\@writefile{toc}{\contentsline {subsection}{\numberline {2.1}The Jacobian}{3}{}\protected@file@percent } +\@writefile{toc}{\contentsline {subsection}{\numberline {2.2}Preconditioner}{4}{}\protected@file@percent } +\@writefile{toc}{\contentsline {subsection}{\numberline {2.3}A Few Quick Notes}{5}{}\protected@file@percent } +\@writefile{toc}{\contentsline {section}{\numberline {3}Representation in FEM}{5}{}\protected@file@percent } +\@writefile{toc}{\contentsline {subsection}{\numberline {3.1}MFEM Integrators}{5}{}\protected@file@percent } +\@writefile{lot}{\contentsline {table}{\numberline {1}{\ignorespaces Selection of MFEM Mixed Bilinear Form Integrators}}{6}{}\protected@file@percent } +\newlabel{tab:mfem_mixed_integrators}{{1}{6}{}{table.1}{}} +\gdef \@abspage@last{6} diff --git a/docs/derivations/laneEmdenVariationalBlockForm.pdf b/docs/derivations/laneEmdenVariationalBlockForm.pdf new file mode 100644 index 0000000..092993d Binary files /dev/null and b/docs/derivations/laneEmdenVariationalBlockForm.pdf differ diff --git a/docs/derivations/laneEmdenVariationalBlockForm.tex b/docs/derivations/laneEmdenVariationalBlockForm.tex new file mode 100644 index 0000000..6604047 --- /dev/null +++ b/docs/derivations/laneEmdenVariationalBlockForm.tex @@ -0,0 +1,291 @@ +\documentclass{article}[twocolumn] +\usepackage{graphicx} % Required for inserting images +\usepackage{booktabs} +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{geometry} + \geometry{ + letterpaper, + left=20mm, + top=20mm, + right=20mm, + bottom=30mm + } + +\title{Deriving The Full Lane Emden Weak Form} +\author{Emily M. Boudreaux} +\date{April 2025} + +\begin{document} + +\maketitle + +\section{Continuous Variational Form} +We start with the strong form of the Lane-Emden equation in three dimensions. +\begin{align} + \Delta \theta + \theta^{n} = 0 +\end{align} +We put this into weak form by multiplying by some scalar test function $v^{\theta}$ which lives in the Sobolev space $H^{1}(\Omega)$ +\begin{align} + \int_{\Omega}v^{\theta}\Delta\theta dV + \int_{\Omega}v^{\theta}\theta^{n}dV = 0 +\end{align} +Applying Green's first identity +\begin{align} +\oint_{\partial \Omega}v^{\theta}\left(\nabla\theta\cdot\hat{n}\right)dA - \int_{\Omega}\nabla v^{\theta}\cdot\nabla\theta dV + \int_{\Omega}v^{\theta}\theta^{n} dV = 0 +\end{align} +We let the surface integral go to zero as the value of $\theta$, and therefore $v_{\theta}$, at the surface of the domain is physically constrained to equal 0. +\begin{align} + - \int_{\Omega}\nabla v^{\theta}\cdot\nabla\theta dV + \int_{\Omega}v^{\theta}\theta^{n} dV = 0 +\end{align} +Now we define a new variable $\phi\equiv\nabla\theta$ so that we can eventually apply essential boundary conditions to both $\theta$ and $\nabla\theta$. We must also then find the variational form of this expression. For that we multiply by some vector test function, $\vec{v}^{\phi}$ which will live in some vector space (In MFEM we will eventually use a Raviart-Thomas space, denoted $RT^{0}(\Omega)$) +\begin{align} + \int_{\Omega}\vec{v}^{\phi}\cdot\vec{\phi}dV - \int_{\Omega}\vec{v}^{\phi}\cdot\nabla\theta = 0 +\end{align} +So then the final, continuous variational system of equations which we have is +\begin{align} + -\int_{\Omega}\nabla v^{\theta}\cdot\phi dV + \int_{\Omega}v^{\theta}\theta^{n}dV &= 0 \\ + \int_{\Omega}\vec{v}^{\phi}\cdot\vec{\phi}dV - \int_{\Omega}\vec{v}^{\phi}\cdot\nabla\theta &= 0 +\end{align} +\section{Discritized Variational Form} +In order to work with this in FEM we need to discritize this. First, Let $\theta_{h}$ be some discrete approximation of $\theta$ which lives on $v^{\theta}_{h}\subset v^{\theta}$ such that +\begin{align} + \theta_{h} = \sum_{i=1}^{N^{\theta}_{dof}}\theta_{i}N_{i}^{\theta} +\end{align} +Where $\{N_{i}^{\theta}\}_{i=1}^{N^{\theta}_{dof}}$ is a set of basis functions which span $v_{h}^{\theta}$ and $\theta_{i}$ are scalar degrees of freedom associated to each basis function. Similarly we can discritize $\phi$ by first letting $\phi_{h}$ be a discrete approximation of $\phi$ which live on $v_{h}^{\phi}\subset v^{\phi}$ +\begin{align} +\phi_{h} = \sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}N_{j}^{\phi} +\end{align} +where $\{N_{j}^{\phi}\}_{j=1}^{N_{dof}^{\phi}}$ is a set of vector basis functions which span $v_{h}^{\phi}$. Let us further let the column vectors $\bar{\theta}$ and $\bar{\phi}$ be +\begin{align} + \bar\theta &\equiv \left[\theta_{1},\ ...\ ,\theta_{N^{\theta}_{dof}}\right]^{T}\\ + \bar\phi &\equiv \left[\phi_{1},\ ...\ ,\phi_{N^{\phi}_{dof}}\right]^{T} +\end{align} + +In order to discritize the weak form we need to adopt a method of selecting test functions for $\theta$ and $\phi$. In the Galerkin method the test functions are chosen from the same finite dimensional subspaces which the approximate solutions are defined on. This is typically done by selecting each basis function to be a test function. This means then that we approximate +\begin{align} +v^{\theta}_{h} &= N_{k}^{\theta} \ \ \ \forall \ \ \ k=1,\ ...\ ,N^{\theta}_{dof} \\ +v^{\phi}_{h} &= N_{\ell}^{\phi} \ \ \ \forall \ \ \ \ell=1,\ ...\ ,N^{\phi}_{dof} +\end{align} +We can now substitute these discritzed expressions for $\theta$, $\phi$, $v^{\theta}$, and $v^{\phi}$ back into the weak form... +\begin{align} + -\int_{\Omega}\nabla N_{k}^{\theta}\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV + \int_{\Omega}N_{k}^{\theta}\left(\sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}N_{i}^{\theta}\right)^{n}dV &= 0 \\ + \int_{\Omega}\vec{N}_{\ell}^{\phi}\cdot\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV - \int_{\Omega}\vec{N}_{\ell}^{\phi}\cdot\nabla\left(\sum_{i=1}^{N^{\theta}_{dof}}\theta_{i}N_{i}^{\theta}\right)dV &= 0 +\end{align} + +I want to pause here and make a note of a point of symbolics which might become confusing latter. We are going to be substituting the basis function, $N^{a}_{b}$, into various places in these equations. However, depending on if we substitute them in for the test functions, $v^{a}$, or the trial functions, $\theta$ and $\phi$, the semantic meaning of those basis functions changes. Any basis function set, $N^{a}_{b}$, used to represented a test function will eventually represent the \textbf{range} of the operator; whereas, any basis function set used to represent a trial function will eventually represent the \textbf{domain} of the operator. This becomes confusing since we use the same symbolics for them. Therefore, for the rest of this derivation I will use $N^{a}_{b}$ to represent the trial function basis set and $\psi^{a}_{b}$ to represent the test function basis set. Using this new symbology we can rewrite the previous two equations as the equivalent forms +\begin{align} + -\int_{\Omega}\nabla \psi_{k}^{\theta}\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV + \int_{\Omega}\psi_{k}^{\theta}\left(\sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}N_{i}^{\theta}\right)^{n}dV &= 0 \\ + \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV - \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\nabla\left(\sum_{i=1}^{N^{\theta}_{dof}}\theta_{i}N_{i}^{\theta}\right)dV &= 0 +\end{align} +Now we exploit the linearity of summation and integration to move the sums out of the integrals +\begin{align} +-\sum_{j=1}^{N_{dof}^{\phi}}\phi_{j}\int_{\Omega}\nabla \psi_{k}^{\theta}\cdot\vec{N}_{j}^{\phi}dV + \int_{\Omega}\psi_{k}^{\theta}\left(\theta_{h}\right)^{n} &= 0 \\ +\sum_{j=1}^{N_{dof}^{\phi}}\phi_{j}\int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot \vec{N}_{j}^{\phi}dV - \sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}\int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\nabla N_{i}^{\theta} dV &=0 +\end{align} +We will now define $M_{kj}$, $D_{\ell j}$, and $Q_{\ell i}$ such that +\begin{align} + M_{kj} &\equiv \int_{\Omega}\nabla \psi_{k}^{\theta}\cdot \vec{N}_{j}^{\phi}dV \\ + D_{\ell j} &\equiv \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\vec{N}_{j}^{\phi}dV \\ + Q_{\ell i} &\equiv \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\nabla N_{i}^{\theta} dV +\end{align} +Further we will define $\mathbf{M},\ \mathbf{D},\ \text{and} \ \mathbf{Q}$ such that they are the matrices associated to $M_{kj}$, $D_{\ell j}$, and $Q_{\ell j}$. + +Note that we do not define a matrix for the non-linear term. This is because we need to treat it as a separate term in computational FEM software, so it is useful for us to split it out now. Instead, let us define $f(\theta)$ to handle the non linear term such that +\begin{align} +f(\bar{\theta}) \equiv \int_{\Omega}\psi_{k}^{\theta}\left(\theta_{h}\right)^{n}dV +\end{align} +We can write the variational form of our system of equations as +\begin{align} +-\sum_{j=1}^{N_{dof}^{\phi}}\phi_{j}M_{kj} + f(\bar{\theta)} &= 0 \\ +\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}D_{\ell j} - \sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}Q_{\ell i} &= 0 +\end{align} +Or using the notation we defined +\begin{align} +-\mathbf{M}\bar{\phi} + f(\bar{\theta}) &= 0 \\ +\mathbf{D}\bar{\phi} - \mathbf{Q}\bar{\theta} &= 0 +\end{align} +We can then set this up as a matrix operation +\begin{align} + \begin{bmatrix} + 0 & -\mathbf{M} \\ + -\mathbf{Q} & \mathbf{D} + \end{bmatrix} + \begin{bmatrix} + \bar{\theta} \\ + \bar{\phi} + \end{bmatrix} + \begin{bmatrix} + f(\bar{\theta}) \\ + 0 + \end{bmatrix} = \begin{bmatrix} + 0 \\ + 0 + \end{bmatrix} +\end{align} + +From this form we can easily see that the residual matrix is + +\begin{align} + R &= \begin{bmatrix} + f(\bar{\theta}) - M\bar{\phi} \\ + D\bar{\phi} - Q\bar{\theta} + \end{bmatrix} +\end{align} + +\subsection{The Jacobian} +We need to define the Jacobian of this system of equations so that we can use it +in our Newton-Raphson method. Generally the Jacobian is the matrix of partial derivitives wrt. the state vector. We will let our state vector, $X$, be +\begin{align} + X = \begin{bmatrix} + \bar{\theta} \\ + \bar{\phi} + \end{bmatrix} +\end{align} + +So then the Jacobian is +\begin{align} + J &= \begin{bmatrix} + \frac{\partial}{\partial \theta}\left(f(\theta) - M\phi\right) & \frac{\partial}{\partial \phi}\left(f(\theta) - M\phi\right) \\ + \frac{\partial}{\partial \theta}\left(D\phi - Q\theta\right) & \frac{\partial}{\partial \phi}\left(D\phi - Q\theta\right) + \end{bmatrix} \\ + J &= \begin{bmatrix} + \frac{df}{d\theta} - \phi\frac{\partial M}{\partial \theta} & -M-\phi\frac{\partial M}{\partial \phi} \\ + -Q - \theta\frac{\partial Q}{\partial \theta} & D + \phi\frac{\partial D}{\partial \phi} - \theta\frac{\partial Q}{\partial \phi} + \end{bmatrix} +\end{align} + +Finally, we know that the matrices $M$, $D$, and $Q$ are constant with respect to $\theta$ and $\phi$. Therefore, we can drop the partial derivatives with respect to $\theta$ and $\phi$ from the Jacobian. This gives us +\begin{align} + \mathbf{J} &= \begin{bmatrix} + \frac{df}{d\theta} & -M \\ + -Q & D + \end{bmatrix} +\end{align} + +\noindent In a fully assembled, distritized form this will look like + +\begin{align} + \mathbf{J} = \begin{bmatrix} \frac{df}{d\theta}_{00} & \dots & \frac{df}{d\theta}_{0n_{\theta}} & -M_{00} & \dots & -M_{0n_{\phi}} \\ + \vdots & \ddots & & \vdots & \ddots & \\ + \frac{df}{d\theta}_{n_{\theta}0} & & \frac{df}{d\theta}_{n_{\theta}n_{\theta}} & -M_{n_{\theta}0} & & -M_{n_{\theta}n_{\phi}} \\ + -Q_{00} & \dots & -Q_{0n_{\theta}} & D_{00} & \dots & D_{0n_{\phi}} \\ + \vdots & \ddots & & \vdots & \ddots & \\ + -Q_{n_{\phi}0} & & -Q_{n_{\phi}n_{\theta}} & D_{n_{\phi}0} & & D_{n_{\phi}n_{\phi}} + \end{bmatrix} +\end{align} + +\noindent Where $N_{dof}^{\theta} = n_{\theta}$ is the number of degrees of freedom on $\theta$, and $N_{dof}^{\phi} = n_{\phi}$ is the number of degrees of freedom on $\phi$. Note how the Jacobian is a matrix of size $\left(N_{dof}^{\theta} + N_{dof}^{\phi} \times N_{dof}^{\theta} + N_{dof}^{\phi}\right)$ + +\subsection{Preconditioner} +Due to the eventual size of these matrices we would like to be able to solve +each step in this using a memory efficiet approach. Krylov solvers, such as +GMRES, allow for matrix free iterative solutions (as long the concept of +multiplication is defined). However, for these systems to be well formed for +such solvers it is useful for us to use a preconditioner. However, this is a +somewhat strongly coupled system where we cannot simply use the inverse +diagonals of the matrix as a preconditioner. Instead, to encode the coupling we +will use Schur's Compliment. Each Newton iteration we solve the equation +\begin{align} + \mathbf{J}\Delta \vec{x} = \vec{b} +\end{align} +If we expand this out +\begin{align} + \begin{bmatrix} \mathbf{\dot{f}} & -\mathbf{M} \\ + -\mathbf{Q} & \mathbf{D} + \end{bmatrix}\begin{bmatrix} \theta \\ + \phi + \end{bmatrix} = \begin{bmatrix} b_{0} \\ + b_{1}\end{bmatrix} +\end{align} +We can pull out the first equation from this system +\begin{align} + \mathbf{\dot{f}}\theta - \mathbf{M}\phi &= b_{0} \\ + \theta &= \mathbf{\dot{f}}^{-1}b_{0} + \mathbf{\dot{f}}^{-1}\mathbf{M}\phi +\end{align} +Then if we pull out the second equation from the system +\begin{align} + -\mathbf{Q}\theta + \mathbf{D}\phi &= b_{1} \\ + -\mathbf{Q}\left(\mathbf{\dot{f}}^{-1}b_{0} + \mathbf{\dot{f}}^{-1}\mathbf{M}\phi\right) + \mathbf{D}\phi &= b_{1} +\end{align} +rearanging terms a bit +\begin{align} + -\mathbf{Q}\mathbf{\dot{f}}^{-1}b_{0}-\mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}\phi + \mathbf{D}\phi &= b_{1} \\ + \left(\mathbf{D} - \mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}\right)\phi &= b_{1} + \mathbf{Q}\mathbf{\dot{f}}^{-1}b_{0} +\end{align} +The term $\mathbf{D}-\mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}$ is Schur's Compliment for this system, and we will represent this by the symbol $\mathbf{\tilde{S}}$. We can use Schur's Compilment to precondition our equation if we let the preconditioner be of the form +\begin{align} + \mathbf{P} = \begin{bmatrix} \mathbf{\dot{f}}^{-1} & 0 \\ + 0 & \mathbf{\tilde{S}}^{-1} + \end{bmatrix} +\end{align} +So then the preconditioned equation which can be more easily solved by some Krylov solver (such as GMRES) is +\begin{align} + \mathbf{P}\mathbf{J}\Delta \vec{x} = \mathbf{P}\vec{b} +\end{align} + +It is easy to see here that for this system to be solvable / well defined both +$\mathbf{\tilde{S}}$ and $\mathbf{\dot{f}}$ need to be invertable +matrices. They are both easily shown to be square (with $\mathbf{\tilde{S}}$ +having a size $\left(N_{dof}^{\phi}\times N_{dof}^{\phi}\right)$ and +$\mathbf{\dot{f}}$ having a size $\left(N_{dof}^{\theta}\times +N_{dof}^{\theta}\right)$). + + +\subsection{A Few Quick Notes} +A few notes on the dimensions of $\mathbf{M}$, $\mathbf{Q}$, $\mathbf{D}$, and $f(\bar{\theta})$. +\begin{itemize} + \item $\mathbf{M}$ is a matrix of size $\left(N_{dof}^{\theta}\ \times\ N_{dof}^{\phi}\right)$. + \item $\mathbf{Q}$ is a matrix of size $\left(N_{dof}^{\phi}\ \times\ N_{dof}^{\theta}\right)$. + \item $\mathbf{D}$ is a matrix of size $\left(N_{dof}^{\phi}\ \times\ N_{dof}^{\phi}\right)$. + \item $f(\bar{\theta})$ is a vector of size $N_{dof}^{\theta}$. +\end{itemize} + +\section{Representation in FEM} +We will make use of the MFEM library\footnote{https://mfem.org/} to encode this system of equations. This document is not intended to be a comprehensive guide to using MFEM; rather, here we will provide an explanation for how to translate $\mathbf{M}$, $\mathbf{D}$, and $\mathbf{Q}$ into pre-existing MFEM integrators. The non linear term must be written as a custom integrator and an explanation of this process is outside the scope of this document. + +\subsection{MFEM Integrators} +MFEM provides an extensive set of integrators. Of interest here are the BilinearFormIntegrators and MixedBilinearFormIntegrators. We will explain how to use these by following the process of deciding how $\mathbf{M}$ should be represented. Recall that +\begin{align} + \mathbf{M} &= \left[M_{kj}\right] \\ + M_{kj} &= \int_{\Omega}\nabla\psi_{k}^{\theta}\cdot\vec{N}_{j}^{\phi}dV +\end{align} +Also recall that $\psi$ denotes the test space while $N$ denotes the trial space. MFEM provides a robust set of integrators. Because $\mathbf{M}$ is composed of terms from the $\theta$ and $\psi$ spaces it is what is called a Mixed form. Therefore, we will look at the mixed form integrators provided by MFEM. +\begin{table}[htbp] +\scriptsize +\centering +\caption{Selection of MFEM Mixed Bilinear Form Integrators} +\label{tab:mfem_mixed_integrators} +% Adjust column specifications (l, c, r, p{}) as needed for alignment and width +\begin{tabular}{@{}lllllll@{}} +\toprule +Class Name & Domain & Range & Coef. & Operator & Continuous Op. & Dimension \\ +\midrule +MixedDotProductIntegrator & ND, RT & H1, L2 & V & $(\vec{\lambda} \cdot \vec{u}, v)$ & $\vec{\lambda} \cdot \vec{u}$ & 2D, 3D \\ +MixedScalarCrossProductIntegrator & ND, RT & H1, L2 & V & $(\vec{\lambda} \times \vec{u}, v)$ & $\vec{\lambda} \times \vec{u}$ & 2D \\ +MixedVectorWeakDivergenceIntegrator & ND, RT & H1 & S, D, M & $(-\lambda \vec{u}, \nabla v)$ & $\nabla \cdot (\lambda \vec{u})$ & 2D, 3D \\ +MixedWeakDivCrossIntegrator & ND, RT & H1 & V & $(-\vec{\lambda} \times \vec{u}, \nabla v)$ & $\nabla \cdot (\vec{\lambda} \times \vec{u})$ & 3D \\ +MixedVectorMassIntegrator & ND, RT & ND, RT & S, D, M & $(\lambda \vec{u}, \vec{v})$ & $\lambda \vec{u}$ & 2D, 3D \\ +MixedCrossProductIntegrator & ND, RT & ND, RT & V & $(\vec{\lambda} \times \vec{u}, \vec{v})$ & $\vec{\lambda} \times \vec{u}$ & 3D \\ +MixedVectorWeakCurlIntegrator & ND, RT & ND & S, D, M & $(\lambda \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\lambda \vec{u})$ & 3D \\ +MixedWeakCurlCrossIntegrator & ND, RT & ND & V & $(\vec{\lambda} \times \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\vec{\lambda} \times \vec{u})$ & 3D \\ +MixedScalarWeakCurlCrossIntegrator & ND, RT & ND & V & $(\vec{\lambda} \times \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\vec{\lambda} \times \vec{u})$ & 2D \\ % Note: Same Op as previous? Verify source. +MixedWeakGradDotIntegrator & ND, RT & RT & V & $(-\vec{\lambda} \cdot \vec{u}, \nabla \cdot \vec{v})$ & $\nabla(\vec{\lambda} \cdot \vec{u})$ & 2D, 3D \\ +MixedScalarCurlIntegrator & ND & H1, L2 & S & $(\lambda \nabla \times \vec{u}, v)$ & $\lambda \nabla \times \vec{u}$ & 2D \\ +MixedCrossCurlGradIntegrator & ND & H1 & V & $(\vec{\lambda} \times (\nabla \times \vec{u}), \nabla v)$ & $-\nabla \cdot (\vec{\lambda} \times (\nabla \times \vec{u}))$ & 3D \\ +MixedVectorCurlIntegrator & ND & ND, RT & S, D, M & $(\lambda \nabla \times \vec{u}, \vec{v})$ & $\lambda \nabla \times \vec{u}$ & 3D \\ +MixedCrossCurlIntegrator & ND & ND, RT & V & $(\vec{\lambda} \times (\nabla \times \vec{u}), \vec{v})$ & $\vec{\lambda} \times (\nabla \times \vec{u})$ & 3D \\ +MixedScalarCrossCurlIntegrator & ND & ND, RT & V & $(\vec{\lambda} \times \hat{z} (\nabla \times \vec{u}), \vec{v})$ & $\vec{\lambda} \times \hat{z} (\nabla \times \vec{u})$ & 2D \\ +MixedCurlCurlIntegrator & ND & ND & S, D, M & $(\lambda \nabla \times \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\lambda \nabla \times \vec{u})$ & 3D \\ +MixedCrossCurlCurlIntegrator & ND & ND & V & $(\vec{\lambda} \times (\nabla \times \vec{u}), \nabla \times \vec{v})$ & $\nabla \times (\vec{\lambda} \times (\nabla \times \vec{u}))$ & 3D \\ +MixedScalarDivergenceIntegrator & RT & H1, L2 & S & $(\lambda \nabla \cdot \vec{u}, v)$ & $\lambda \nabla \cdot \vec{u}$ & 2D, 3D \\ +MixedDivGradIntegrator & RT & H1 & V & $(\vec{\lambda} (\nabla \cdot \vec{u}), \nabla v)$ & $-\nabla \cdot (\vec{\lambda} (\nabla \cdot \vec{u}))$ & 2D, 3D \\ +MixedVectorDivergenceIntegrator & RT & ND, RT & V & $(\vec{\lambda} (\nabla \cdot \vec{u}), \vec{v})$ & $\vec{\lambda} (\nabla \cdot \vec{u})$ & 2D, 3D \\ +\bottomrule +\end{tabular} +\end{table} + +There is a lot of information in this table so we will break it down. First we need to identify what spaces the domain and range of our operator exist in. The range of the operator is that which contains the test function while the domain is the space containing the trial function. For $\mathbf{M}$ the test function is in the $\theta$ space, or H1, while the trial function is in the $\phi$ space, or RT. + +Next we look at the Operator column. These define the operation within the integral where $(a,b)$ is the inner product of $a$ and $b$. More specifically, the MFEM documentation provides that $u$ is the trial function and $v$ is the test function. So then we are looking for an integrator which has the operator of the inner product of the trial function and the gradient of the test function while also satisfying the domain and range constraints. Upon investigation of this table we can see that the \texttt{MixedVectorWeakDivegenceIntegrator} has a range defined on H1 and domain defined on RT, just like we need. Further, its operator is given as $(-\lambda \vec{u},\nabla v)$. This is the same form as we have. $\left(-\lambda\vec{u} \rightarrow 1\times\vec{N}_{j}^{\phi}, \nabla v\rightarrow \nabla \psi_{k}^{\theta}\right)$. Therefore, \texttt{MixedVectorWeakDivergenceIntegrator} will compute the matrix $\textbf{M}$ over our domain without any modifications (as long as we are careful about the sign on the coefficient $\lambda$). + +The other integrators map to \texttt{VectorFEMassIntegrator} and \texttt{MixedVectorGradientIntegrator} for $\mathbf{D}$ and $\mathbf{Q}$ respectively. From here one would assemble these, along with the non-linear term into a block form; however, that is beyond the scope of this document. +\end{document} diff --git a/docs/derivations/laneEmdenVariationalForm/.ipynb_checkpoints/JacobianEigenValuesPostPreconditioning-checkpoint.pdf b/docs/derivations/laneEmdenVariationalForm/.ipynb_checkpoints/JacobianEigenValuesPostPreconditioning-checkpoint.pdf new file mode 100644 index 0000000..b4cce6e Binary files /dev/null and b/docs/derivations/laneEmdenVariationalForm/.ipynb_checkpoints/JacobianEigenValuesPostPreconditioning-checkpoint.pdf differ diff --git a/docs/derivations/laneEmdenVariationalForm/JacobianEigenValuesPostPreconditioning.pdf b/docs/derivations/laneEmdenVariationalForm/JacobianEigenValuesPostPreconditioning.pdf new file mode 100644 index 0000000..b4cce6e Binary files /dev/null and b/docs/derivations/laneEmdenVariationalForm/JacobianEigenValuesPostPreconditioning.pdf differ diff --git a/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.aux b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.aux new file mode 100644 index 0000000..5dcb826 --- /dev/null +++ b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.aux @@ -0,0 +1,11 @@ +\relax +\@writefile{toc}{\contentsline {section}{\numberline {1}Continuous Variational Form}{1}{}\protected@file@percent } +\@writefile{toc}{\contentsline {section}{\numberline {2}Discritized Variational Form}{1}{}\protected@file@percent } +\@writefile{toc}{\contentsline {subsection}{\numberline {2.1}The Jacobian}{3}{}\protected@file@percent } +\@writefile{toc}{\contentsline {subsection}{\numberline {2.2}Preconditioner}{4}{}\protected@file@percent } +\@writefile{toc}{\contentsline {subsection}{\numberline {2.3}A Few Quick Notes}{5}{}\protected@file@percent } +\@writefile{toc}{\contentsline {section}{\numberline {3}Representation in FEM}{5}{}\protected@file@percent } +\@writefile{toc}{\contentsline {subsection}{\numberline {3.1}MFEM Integrators}{5}{}\protected@file@percent } +\@writefile{lot}{\contentsline {table}{\numberline {1}{\ignorespaces Selection of MFEM Mixed Bilinear Form Integrators}}{6}{}\protected@file@percent } +\newlabel{tab:mfem_mixed_integrators}{{1}{6}{}{table.1}{}} +\gdef \@abspage@last{6} diff --git a/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.pdf b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.pdf new file mode 100644 index 0000000..092993d Binary files /dev/null and b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.pdf differ diff --git a/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.tex b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.tex new file mode 100644 index 0000000..6604047 --- /dev/null +++ b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.tex @@ -0,0 +1,291 @@ +\documentclass{article}[twocolumn] +\usepackage{graphicx} % Required for inserting images +\usepackage{booktabs} +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{geometry} + \geometry{ + letterpaper, + left=20mm, + top=20mm, + right=20mm, + bottom=30mm + } + +\title{Deriving The Full Lane Emden Weak Form} +\author{Emily M. Boudreaux} +\date{April 2025} + +\begin{document} + +\maketitle + +\section{Continuous Variational Form} +We start with the strong form of the Lane-Emden equation in three dimensions. +\begin{align} + \Delta \theta + \theta^{n} = 0 +\end{align} +We put this into weak form by multiplying by some scalar test function $v^{\theta}$ which lives in the Sobolev space $H^{1}(\Omega)$ +\begin{align} + \int_{\Omega}v^{\theta}\Delta\theta dV + \int_{\Omega}v^{\theta}\theta^{n}dV = 0 +\end{align} +Applying Green's first identity +\begin{align} +\oint_{\partial \Omega}v^{\theta}\left(\nabla\theta\cdot\hat{n}\right)dA - \int_{\Omega}\nabla v^{\theta}\cdot\nabla\theta dV + \int_{\Omega}v^{\theta}\theta^{n} dV = 0 +\end{align} +We let the surface integral go to zero as the value of $\theta$, and therefore $v_{\theta}$, at the surface of the domain is physically constrained to equal 0. +\begin{align} + - \int_{\Omega}\nabla v^{\theta}\cdot\nabla\theta dV + \int_{\Omega}v^{\theta}\theta^{n} dV = 0 +\end{align} +Now we define a new variable $\phi\equiv\nabla\theta$ so that we can eventually apply essential boundary conditions to both $\theta$ and $\nabla\theta$. We must also then find the variational form of this expression. For that we multiply by some vector test function, $\vec{v}^{\phi}$ which will live in some vector space (In MFEM we will eventually use a Raviart-Thomas space, denoted $RT^{0}(\Omega)$) +\begin{align} + \int_{\Omega}\vec{v}^{\phi}\cdot\vec{\phi}dV - \int_{\Omega}\vec{v}^{\phi}\cdot\nabla\theta = 0 +\end{align} +So then the final, continuous variational system of equations which we have is +\begin{align} + -\int_{\Omega}\nabla v^{\theta}\cdot\phi dV + \int_{\Omega}v^{\theta}\theta^{n}dV &= 0 \\ + \int_{\Omega}\vec{v}^{\phi}\cdot\vec{\phi}dV - \int_{\Omega}\vec{v}^{\phi}\cdot\nabla\theta &= 0 +\end{align} +\section{Discritized Variational Form} +In order to work with this in FEM we need to discritize this. First, Let $\theta_{h}$ be some discrete approximation of $\theta$ which lives on $v^{\theta}_{h}\subset v^{\theta}$ such that +\begin{align} + \theta_{h} = \sum_{i=1}^{N^{\theta}_{dof}}\theta_{i}N_{i}^{\theta} +\end{align} +Where $\{N_{i}^{\theta}\}_{i=1}^{N^{\theta}_{dof}}$ is a set of basis functions which span $v_{h}^{\theta}$ and $\theta_{i}$ are scalar degrees of freedom associated to each basis function. Similarly we can discritize $\phi$ by first letting $\phi_{h}$ be a discrete approximation of $\phi$ which live on $v_{h}^{\phi}\subset v^{\phi}$ +\begin{align} +\phi_{h} = \sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}N_{j}^{\phi} +\end{align} +where $\{N_{j}^{\phi}\}_{j=1}^{N_{dof}^{\phi}}$ is a set of vector basis functions which span $v_{h}^{\phi}$. Let us further let the column vectors $\bar{\theta}$ and $\bar{\phi}$ be +\begin{align} + \bar\theta &\equiv \left[\theta_{1},\ ...\ ,\theta_{N^{\theta}_{dof}}\right]^{T}\\ + \bar\phi &\equiv \left[\phi_{1},\ ...\ ,\phi_{N^{\phi}_{dof}}\right]^{T} +\end{align} + +In order to discritize the weak form we need to adopt a method of selecting test functions for $\theta$ and $\phi$. In the Galerkin method the test functions are chosen from the same finite dimensional subspaces which the approximate solutions are defined on. This is typically done by selecting each basis function to be a test function. This means then that we approximate +\begin{align} +v^{\theta}_{h} &= N_{k}^{\theta} \ \ \ \forall \ \ \ k=1,\ ...\ ,N^{\theta}_{dof} \\ +v^{\phi}_{h} &= N_{\ell}^{\phi} \ \ \ \forall \ \ \ \ell=1,\ ...\ ,N^{\phi}_{dof} +\end{align} +We can now substitute these discritzed expressions for $\theta$, $\phi$, $v^{\theta}$, and $v^{\phi}$ back into the weak form... +\begin{align} + -\int_{\Omega}\nabla N_{k}^{\theta}\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV + \int_{\Omega}N_{k}^{\theta}\left(\sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}N_{i}^{\theta}\right)^{n}dV &= 0 \\ + \int_{\Omega}\vec{N}_{\ell}^{\phi}\cdot\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV - \int_{\Omega}\vec{N}_{\ell}^{\phi}\cdot\nabla\left(\sum_{i=1}^{N^{\theta}_{dof}}\theta_{i}N_{i}^{\theta}\right)dV &= 0 +\end{align} + +I want to pause here and make a note of a point of symbolics which might become confusing latter. We are going to be substituting the basis function, $N^{a}_{b}$, into various places in these equations. However, depending on if we substitute them in for the test functions, $v^{a}$, or the trial functions, $\theta$ and $\phi$, the semantic meaning of those basis functions changes. Any basis function set, $N^{a}_{b}$, used to represented a test function will eventually represent the \textbf{range} of the operator; whereas, any basis function set used to represent a trial function will eventually represent the \textbf{domain} of the operator. This becomes confusing since we use the same symbolics for them. Therefore, for the rest of this derivation I will use $N^{a}_{b}$ to represent the trial function basis set and $\psi^{a}_{b}$ to represent the test function basis set. Using this new symbology we can rewrite the previous two equations as the equivalent forms +\begin{align} + -\int_{\Omega}\nabla \psi_{k}^{\theta}\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV + \int_{\Omega}\psi_{k}^{\theta}\left(\sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}N_{i}^{\theta}\right)^{n}dV &= 0 \\ + \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\left(\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}\vec{N_{j}^{\phi}}\right)dV - \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\nabla\left(\sum_{i=1}^{N^{\theta}_{dof}}\theta_{i}N_{i}^{\theta}\right)dV &= 0 +\end{align} +Now we exploit the linearity of summation and integration to move the sums out of the integrals +\begin{align} +-\sum_{j=1}^{N_{dof}^{\phi}}\phi_{j}\int_{\Omega}\nabla \psi_{k}^{\theta}\cdot\vec{N}_{j}^{\phi}dV + \int_{\Omega}\psi_{k}^{\theta}\left(\theta_{h}\right)^{n} &= 0 \\ +\sum_{j=1}^{N_{dof}^{\phi}}\phi_{j}\int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot \vec{N}_{j}^{\phi}dV - \sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}\int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\nabla N_{i}^{\theta} dV &=0 +\end{align} +We will now define $M_{kj}$, $D_{\ell j}$, and $Q_{\ell i}$ such that +\begin{align} + M_{kj} &\equiv \int_{\Omega}\nabla \psi_{k}^{\theta}\cdot \vec{N}_{j}^{\phi}dV \\ + D_{\ell j} &\equiv \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\vec{N}_{j}^{\phi}dV \\ + Q_{\ell i} &\equiv \int_{\Omega}\vec{\psi}_{\ell}^{\phi}\cdot\nabla N_{i}^{\theta} dV +\end{align} +Further we will define $\mathbf{M},\ \mathbf{D},\ \text{and} \ \mathbf{Q}$ such that they are the matrices associated to $M_{kj}$, $D_{\ell j}$, and $Q_{\ell j}$. + +Note that we do not define a matrix for the non-linear term. This is because we need to treat it as a separate term in computational FEM software, so it is useful for us to split it out now. Instead, let us define $f(\theta)$ to handle the non linear term such that +\begin{align} +f(\bar{\theta}) \equiv \int_{\Omega}\psi_{k}^{\theta}\left(\theta_{h}\right)^{n}dV +\end{align} +We can write the variational form of our system of equations as +\begin{align} +-\sum_{j=1}^{N_{dof}^{\phi}}\phi_{j}M_{kj} + f(\bar{\theta)} &= 0 \\ +\sum_{j=1}^{N^{\phi}_{dof}}\phi_{j}D_{\ell j} - \sum_{i=1}^{N_{dof}^{\theta}}\theta_{i}Q_{\ell i} &= 0 +\end{align} +Or using the notation we defined +\begin{align} +-\mathbf{M}\bar{\phi} + f(\bar{\theta}) &= 0 \\ +\mathbf{D}\bar{\phi} - \mathbf{Q}\bar{\theta} &= 0 +\end{align} +We can then set this up as a matrix operation +\begin{align} + \begin{bmatrix} + 0 & -\mathbf{M} \\ + -\mathbf{Q} & \mathbf{D} + \end{bmatrix} + \begin{bmatrix} + \bar{\theta} \\ + \bar{\phi} + \end{bmatrix} + \begin{bmatrix} + f(\bar{\theta}) \\ + 0 + \end{bmatrix} = \begin{bmatrix} + 0 \\ + 0 + \end{bmatrix} +\end{align} + +From this form we can easily see that the residual matrix is + +\begin{align} + R &= \begin{bmatrix} + f(\bar{\theta}) - M\bar{\phi} \\ + D\bar{\phi} - Q\bar{\theta} + \end{bmatrix} +\end{align} + +\subsection{The Jacobian} +We need to define the Jacobian of this system of equations so that we can use it +in our Newton-Raphson method. Generally the Jacobian is the matrix of partial derivitives wrt. the state vector. We will let our state vector, $X$, be +\begin{align} + X = \begin{bmatrix} + \bar{\theta} \\ + \bar{\phi} + \end{bmatrix} +\end{align} + +So then the Jacobian is +\begin{align} + J &= \begin{bmatrix} + \frac{\partial}{\partial \theta}\left(f(\theta) - M\phi\right) & \frac{\partial}{\partial \phi}\left(f(\theta) - M\phi\right) \\ + \frac{\partial}{\partial \theta}\left(D\phi - Q\theta\right) & \frac{\partial}{\partial \phi}\left(D\phi - Q\theta\right) + \end{bmatrix} \\ + J &= \begin{bmatrix} + \frac{df}{d\theta} - \phi\frac{\partial M}{\partial \theta} & -M-\phi\frac{\partial M}{\partial \phi} \\ + -Q - \theta\frac{\partial Q}{\partial \theta} & D + \phi\frac{\partial D}{\partial \phi} - \theta\frac{\partial Q}{\partial \phi} + \end{bmatrix} +\end{align} + +Finally, we know that the matrices $M$, $D$, and $Q$ are constant with respect to $\theta$ and $\phi$. Therefore, we can drop the partial derivatives with respect to $\theta$ and $\phi$ from the Jacobian. This gives us +\begin{align} + \mathbf{J} &= \begin{bmatrix} + \frac{df}{d\theta} & -M \\ + -Q & D + \end{bmatrix} +\end{align} + +\noindent In a fully assembled, distritized form this will look like + +\begin{align} + \mathbf{J} = \begin{bmatrix} \frac{df}{d\theta}_{00} & \dots & \frac{df}{d\theta}_{0n_{\theta}} & -M_{00} & \dots & -M_{0n_{\phi}} \\ + \vdots & \ddots & & \vdots & \ddots & \\ + \frac{df}{d\theta}_{n_{\theta}0} & & \frac{df}{d\theta}_{n_{\theta}n_{\theta}} & -M_{n_{\theta}0} & & -M_{n_{\theta}n_{\phi}} \\ + -Q_{00} & \dots & -Q_{0n_{\theta}} & D_{00} & \dots & D_{0n_{\phi}} \\ + \vdots & \ddots & & \vdots & \ddots & \\ + -Q_{n_{\phi}0} & & -Q_{n_{\phi}n_{\theta}} & D_{n_{\phi}0} & & D_{n_{\phi}n_{\phi}} + \end{bmatrix} +\end{align} + +\noindent Where $N_{dof}^{\theta} = n_{\theta}$ is the number of degrees of freedom on $\theta$, and $N_{dof}^{\phi} = n_{\phi}$ is the number of degrees of freedom on $\phi$. Note how the Jacobian is a matrix of size $\left(N_{dof}^{\theta} + N_{dof}^{\phi} \times N_{dof}^{\theta} + N_{dof}^{\phi}\right)$ + +\subsection{Preconditioner} +Due to the eventual size of these matrices we would like to be able to solve +each step in this using a memory efficiet approach. Krylov solvers, such as +GMRES, allow for matrix free iterative solutions (as long the concept of +multiplication is defined). However, for these systems to be well formed for +such solvers it is useful for us to use a preconditioner. However, this is a +somewhat strongly coupled system where we cannot simply use the inverse +diagonals of the matrix as a preconditioner. Instead, to encode the coupling we +will use Schur's Compliment. Each Newton iteration we solve the equation +\begin{align} + \mathbf{J}\Delta \vec{x} = \vec{b} +\end{align} +If we expand this out +\begin{align} + \begin{bmatrix} \mathbf{\dot{f}} & -\mathbf{M} \\ + -\mathbf{Q} & \mathbf{D} + \end{bmatrix}\begin{bmatrix} \theta \\ + \phi + \end{bmatrix} = \begin{bmatrix} b_{0} \\ + b_{1}\end{bmatrix} +\end{align} +We can pull out the first equation from this system +\begin{align} + \mathbf{\dot{f}}\theta - \mathbf{M}\phi &= b_{0} \\ + \theta &= \mathbf{\dot{f}}^{-1}b_{0} + \mathbf{\dot{f}}^{-1}\mathbf{M}\phi +\end{align} +Then if we pull out the second equation from the system +\begin{align} + -\mathbf{Q}\theta + \mathbf{D}\phi &= b_{1} \\ + -\mathbf{Q}\left(\mathbf{\dot{f}}^{-1}b_{0} + \mathbf{\dot{f}}^{-1}\mathbf{M}\phi\right) + \mathbf{D}\phi &= b_{1} +\end{align} +rearanging terms a bit +\begin{align} + -\mathbf{Q}\mathbf{\dot{f}}^{-1}b_{0}-\mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}\phi + \mathbf{D}\phi &= b_{1} \\ + \left(\mathbf{D} - \mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}\right)\phi &= b_{1} + \mathbf{Q}\mathbf{\dot{f}}^{-1}b_{0} +\end{align} +The term $\mathbf{D}-\mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}$ is Schur's Compliment for this system, and we will represent this by the symbol $\mathbf{\tilde{S}}$. We can use Schur's Compilment to precondition our equation if we let the preconditioner be of the form +\begin{align} + \mathbf{P} = \begin{bmatrix} \mathbf{\dot{f}}^{-1} & 0 \\ + 0 & \mathbf{\tilde{S}}^{-1} + \end{bmatrix} +\end{align} +So then the preconditioned equation which can be more easily solved by some Krylov solver (such as GMRES) is +\begin{align} + \mathbf{P}\mathbf{J}\Delta \vec{x} = \mathbf{P}\vec{b} +\end{align} + +It is easy to see here that for this system to be solvable / well defined both +$\mathbf{\tilde{S}}$ and $\mathbf{\dot{f}}$ need to be invertable +matrices. They are both easily shown to be square (with $\mathbf{\tilde{S}}$ +having a size $\left(N_{dof}^{\phi}\times N_{dof}^{\phi}\right)$ and +$\mathbf{\dot{f}}$ having a size $\left(N_{dof}^{\theta}\times +N_{dof}^{\theta}\right)$). + + +\subsection{A Few Quick Notes} +A few notes on the dimensions of $\mathbf{M}$, $\mathbf{Q}$, $\mathbf{D}$, and $f(\bar{\theta})$. +\begin{itemize} + \item $\mathbf{M}$ is a matrix of size $\left(N_{dof}^{\theta}\ \times\ N_{dof}^{\phi}\right)$. + \item $\mathbf{Q}$ is a matrix of size $\left(N_{dof}^{\phi}\ \times\ N_{dof}^{\theta}\right)$. + \item $\mathbf{D}$ is a matrix of size $\left(N_{dof}^{\phi}\ \times\ N_{dof}^{\phi}\right)$. + \item $f(\bar{\theta})$ is a vector of size $N_{dof}^{\theta}$. +\end{itemize} + +\section{Representation in FEM} +We will make use of the MFEM library\footnote{https://mfem.org/} to encode this system of equations. This document is not intended to be a comprehensive guide to using MFEM; rather, here we will provide an explanation for how to translate $\mathbf{M}$, $\mathbf{D}$, and $\mathbf{Q}$ into pre-existing MFEM integrators. The non linear term must be written as a custom integrator and an explanation of this process is outside the scope of this document. + +\subsection{MFEM Integrators} +MFEM provides an extensive set of integrators. Of interest here are the BilinearFormIntegrators and MixedBilinearFormIntegrators. We will explain how to use these by following the process of deciding how $\mathbf{M}$ should be represented. Recall that +\begin{align} + \mathbf{M} &= \left[M_{kj}\right] \\ + M_{kj} &= \int_{\Omega}\nabla\psi_{k}^{\theta}\cdot\vec{N}_{j}^{\phi}dV +\end{align} +Also recall that $\psi$ denotes the test space while $N$ denotes the trial space. MFEM provides a robust set of integrators. Because $\mathbf{M}$ is composed of terms from the $\theta$ and $\psi$ spaces it is what is called a Mixed form. Therefore, we will look at the mixed form integrators provided by MFEM. +\begin{table}[htbp] +\scriptsize +\centering +\caption{Selection of MFEM Mixed Bilinear Form Integrators} +\label{tab:mfem_mixed_integrators} +% Adjust column specifications (l, c, r, p{}) as needed for alignment and width +\begin{tabular}{@{}lllllll@{}} +\toprule +Class Name & Domain & Range & Coef. & Operator & Continuous Op. & Dimension \\ +\midrule +MixedDotProductIntegrator & ND, RT & H1, L2 & V & $(\vec{\lambda} \cdot \vec{u}, v)$ & $\vec{\lambda} \cdot \vec{u}$ & 2D, 3D \\ +MixedScalarCrossProductIntegrator & ND, RT & H1, L2 & V & $(\vec{\lambda} \times \vec{u}, v)$ & $\vec{\lambda} \times \vec{u}$ & 2D \\ +MixedVectorWeakDivergenceIntegrator & ND, RT & H1 & S, D, M & $(-\lambda \vec{u}, \nabla v)$ & $\nabla \cdot (\lambda \vec{u})$ & 2D, 3D \\ +MixedWeakDivCrossIntegrator & ND, RT & H1 & V & $(-\vec{\lambda} \times \vec{u}, \nabla v)$ & $\nabla \cdot (\vec{\lambda} \times \vec{u})$ & 3D \\ +MixedVectorMassIntegrator & ND, RT & ND, RT & S, D, M & $(\lambda \vec{u}, \vec{v})$ & $\lambda \vec{u}$ & 2D, 3D \\ +MixedCrossProductIntegrator & ND, RT & ND, RT & V & $(\vec{\lambda} \times \vec{u}, \vec{v})$ & $\vec{\lambda} \times \vec{u}$ & 3D \\ +MixedVectorWeakCurlIntegrator & ND, RT & ND & S, D, M & $(\lambda \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\lambda \vec{u})$ & 3D \\ +MixedWeakCurlCrossIntegrator & ND, RT & ND & V & $(\vec{\lambda} \times \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\vec{\lambda} \times \vec{u})$ & 3D \\ +MixedScalarWeakCurlCrossIntegrator & ND, RT & ND & V & $(\vec{\lambda} \times \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\vec{\lambda} \times \vec{u})$ & 2D \\ % Note: Same Op as previous? Verify source. +MixedWeakGradDotIntegrator & ND, RT & RT & V & $(-\vec{\lambda} \cdot \vec{u}, \nabla \cdot \vec{v})$ & $\nabla(\vec{\lambda} \cdot \vec{u})$ & 2D, 3D \\ +MixedScalarCurlIntegrator & ND & H1, L2 & S & $(\lambda \nabla \times \vec{u}, v)$ & $\lambda \nabla \times \vec{u}$ & 2D \\ +MixedCrossCurlGradIntegrator & ND & H1 & V & $(\vec{\lambda} \times (\nabla \times \vec{u}), \nabla v)$ & $-\nabla \cdot (\vec{\lambda} \times (\nabla \times \vec{u}))$ & 3D \\ +MixedVectorCurlIntegrator & ND & ND, RT & S, D, M & $(\lambda \nabla \times \vec{u}, \vec{v})$ & $\lambda \nabla \times \vec{u}$ & 3D \\ +MixedCrossCurlIntegrator & ND & ND, RT & V & $(\vec{\lambda} \times (\nabla \times \vec{u}), \vec{v})$ & $\vec{\lambda} \times (\nabla \times \vec{u})$ & 3D \\ +MixedScalarCrossCurlIntegrator & ND & ND, RT & V & $(\vec{\lambda} \times \hat{z} (\nabla \times \vec{u}), \vec{v})$ & $\vec{\lambda} \times \hat{z} (\nabla \times \vec{u})$ & 2D \\ +MixedCurlCurlIntegrator & ND & ND & S, D, M & $(\lambda \nabla \times \vec{u}, \nabla \times \vec{v})$ & $\nabla \times (\lambda \nabla \times \vec{u})$ & 3D \\ +MixedCrossCurlCurlIntegrator & ND & ND & V & $(\vec{\lambda} \times (\nabla \times \vec{u}), \nabla \times \vec{v})$ & $\nabla \times (\vec{\lambda} \times (\nabla \times \vec{u}))$ & 3D \\ +MixedScalarDivergenceIntegrator & RT & H1, L2 & S & $(\lambda \nabla \cdot \vec{u}, v)$ & $\lambda \nabla \cdot \vec{u}$ & 2D, 3D \\ +MixedDivGradIntegrator & RT & H1 & V & $(\vec{\lambda} (\nabla \cdot \vec{u}), \nabla v)$ & $-\nabla \cdot (\vec{\lambda} (\nabla \cdot \vec{u}))$ & 2D, 3D \\ +MixedVectorDivergenceIntegrator & RT & ND, RT & V & $(\vec{\lambda} (\nabla \cdot \vec{u}), \vec{v})$ & $\vec{\lambda} (\nabla \cdot \vec{u})$ & 2D, 3D \\ +\bottomrule +\end{tabular} +\end{table} + +There is a lot of information in this table so we will break it down. First we need to identify what spaces the domain and range of our operator exist in. The range of the operator is that which contains the test function while the domain is the space containing the trial function. For $\mathbf{M}$ the test function is in the $\theta$ space, or H1, while the trial function is in the $\phi$ space, or RT. + +Next we look at the Operator column. These define the operation within the integral where $(a,b)$ is the inner product of $a$ and $b$. More specifically, the MFEM documentation provides that $u$ is the trial function and $v$ is the test function. So then we are looking for an integrator which has the operator of the inner product of the trial function and the gradient of the test function while also satisfying the domain and range constraints. Upon investigation of this table we can see that the \texttt{MixedVectorWeakDivegenceIntegrator} has a range defined on H1 and domain defined on RT, just like we need. Further, its operator is given as $(-\lambda \vec{u},\nabla v)$. This is the same form as we have. $\left(-\lambda\vec{u} \rightarrow 1\times\vec{N}_{j}^{\phi}, \nabla v\rightarrow \nabla \psi_{k}^{\theta}\right)$. Therefore, \texttt{MixedVectorWeakDivergenceIntegrator} will compute the matrix $\textbf{M}$ over our domain without any modifications (as long as we are careful about the sign on the coefficient $\lambda$). + +The other integrators map to \texttt{VectorFEMassIntegrator} and \texttt{MixedVectorGradientIntegrator} for $\mathbf{D}$ and $\mathbf{Q}$ respectively. From here one would assemble these, along with the non-linear term into a block form; however, that is beyond the scope of this document. +\end{document} diff --git a/docs/derivations/laneEmdenVariationalForm/makefile b/docs/derivations/laneEmdenVariationalForm/makefile new file mode 100644 index 0000000..cf1ffac --- /dev/null +++ b/docs/derivations/laneEmdenVariationalForm/makefile @@ -0,0 +1,8 @@ +TC=pdflatex +JOBNAME=laneEmdenVariationalBlockForm +SRC=laneEmdenVariationalBlockForm.tex + +default: all + +all: laneEmdenVariationalBlockForm.tex + $(TC) -jobname=$(JOBNAME) $(SRC) \ No newline at end of file diff --git a/docs/derivations/makefile b/docs/derivations/makefile new file mode 100644 index 0000000..cf1ffac --- /dev/null +++ b/docs/derivations/makefile @@ -0,0 +1,8 @@ +TC=pdflatex +JOBNAME=laneEmdenVariationalBlockForm +SRC=laneEmdenVariationalBlockForm.tex + +default: all + +all: laneEmdenVariationalBlockForm.tex + $(TC) -jobname=$(JOBNAME) $(SRC) \ No newline at end of file diff --git a/src/meson.build b/src/meson.build index c4db238..2946320 100644 --- a/src/meson.build +++ b/src/meson.build @@ -3,6 +3,7 @@ # as there are dependencies which exist between them. # Utility Libraries +subdir('types') subdir('misc') subdir('config') subdir('probe') @@ -21,3 +22,4 @@ subdir('resource') # Physics Libraries subdir('network') +subdir('poly') diff --git a/src/poly/coeff/meson.build b/src/poly/coeff/meson.build index 112a1ff..e69de29 100644 --- a/src/poly/coeff/meson.build +++ b/src/poly/coeff/meson.build @@ -1,23 +0,0 @@ -polyCoeff_sources = files( - 'private/coeff.cpp' -) - -polyCoeff_headers = files( - 'public/coeff.h' -) - -libPolyCoeff = static_library('polyCoeff', - polyCoeff_sources, - include_directories : include_directories('.'), - cpp_args: ['-fvisibility=default'], - dependencies: [mfem_dep], - install: true -) - - -polyCoeff_dep = declare_dependency( - include_directories : include_directories('.'), - link_with : libPolyCoeff, - sources : polyCoeff_sources, - dependencies : [mfem_dep] -) \ No newline at end of file diff --git a/src/poly/coeff/private/polyCoeff.cpp b/src/poly/coeff/private/polyCoeff.cpp deleted file mode 100644 index 36a5ab9..0000000 --- a/src/poly/coeff/private/polyCoeff.cpp +++ /dev/null @@ -1,60 +0,0 @@ -/* *********************************************************************** -// -// Copyright (C) 2025 -- The 4D-STAR Collaboration -// File Author: Emily Boudreaux -// Last Modified: February 14, 2025 -// -// 4DSSE is free software; you can use it and/or modify -// it under the terms and restrictions the GNU General Library Public -// License version 3 (GPLv3) as published by the Free Software Foundation. -// -// 4DSSE is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -// See the GNU Library General Public License for more details. -// -// You should have received a copy of the GNU Library General Public License -// along with this software; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// *********************************************************************** */ -#include "mfem.hpp" -#include - -#include "coeff.h" - -/** - * @brief Computes the xi coefficient function. - * - * @param x Input vector. - * @return double The computed xi coefficient. - */ -double xi_coeff_func(const mfem::Vector &x) -{ - return std::pow(x(0), 2); -} - -/** - * @brief Computes the vector xi coefficient function. - * - * @param x Input vector. - * @param v Output vector to store the computed xi coefficient. - */ -void vec_xi_coeff_func(const mfem::Vector &x, mfem::Vector &v) -{ - v.SetSize(1); - v[0] = -std::pow(x(0), 2); -} - -/** - * @brief Computes the initial guess for theta. - * - * @param x Input vector. - * @param root Root value used in the computation. - * @return double The initial guess for theta. - */ -double theta_initial_guess(const mfem::Vector &x, double root) -{ - double xi = x[0]; - return 1 - std::pow(xi / root, 2); -} \ No newline at end of file diff --git a/src/poly/utils/meson.build b/src/poly/utils/meson.build index fc072eb..9b73678 100644 --- a/src/poly/utils/meson.build +++ b/src/poly/utils/meson.build @@ -1,24 +1,48 @@ +# *********************************************************************** +# +# Copyright (C) 2025 -- The 4D-STAR Collaboration +# File Author: Emily Boudreaux +# Last Modified: March 19, 2025 +# +# 4DSSE is free software; you can use it and/or modify +# it under the terms and restrictions the GNU General Library Public +# License version 3 (GPLv3) as published by the Free Software Foundation. +# +# 4DSSE is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +# See the GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General Public License +# along with this software; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# *********************************************************************** # polyutils_sources = files( - 'private/polyIO.cpp', - 'private/polyMFEMUtils.cpp' + 'private/integrators.cpp', + 'private/operator.cpp' ) -polyutils_headers = files( - 'public/polyIO.h', - 'public/polyMFEMUtils.h' -) +dependencies = [ + mfem_dep, + macros_dep, + probe_dep, + quill_dep, + config_dep, + types_dep, +] libpolyutils = static_library('polyutils', polyutils_sources, - include_directories : include_directories('.'), + include_directories : include_directories('./public'), cpp_args: ['-fvisibility=default'], - dependencies: [mfem_dep], + dependencies: dependencies, install: true ) -libpolyutils_dep = declare_dependency( - include_directories : include_directories('.'), +polyutils_dep = declare_dependency( + include_directories : include_directories('./public'), link_with : libpolyutils, sources : polyutils_sources, - dependencies : [mfem_dep] + dependencies : dependencies ) diff --git a/src/poly/utils/private/integrators.cpp b/src/poly/utils/private/integrators.cpp new file mode 100644 index 0000000..e433d62 --- /dev/null +++ b/src/poly/utils/private/integrators.cpp @@ -0,0 +1,135 @@ +/* *********************************************************************** +// +// Copyright (C) 2025 -- The 4D-STAR Collaboration +// File Author: Emily Boudreaux +// Last Modified: April 21, 2025 +// +// 4DSSE is free software; you can use it and/or modify +// it under the terms and restrictions the GNU General Library Public +// License version 3 (GPLv3) as published by the Free Software Foundation. +// +// 4DSSE is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU Library General Public License for more details. +// +// You should have received a copy of the GNU Library General Public License +// along with this software; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// *********************************************************************** */ +#include "mfem.hpp" +#include + +#include "integrators.h" +#include + + +// static std::ofstream debugOut("gradient.csv", std::ios::trunc); + +namespace polyMFEMUtils { + NonlinearPowerIntegrator::NonlinearPowerIntegrator(const double n) : + m_polytropicIndex(n) {} + + void NonlinearPowerIntegrator::AssembleElementVector( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::Vector &elvect) { + + const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); + int dof = el.GetDof(); + elvect.SetSize(dof); + elvect = 0.0; + + mfem::Vector shape(dof); + for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { + mfem::IntegrationPoint ip = ir->IntPoint(iqp); + Trans.SetIntPoint(&ip); + const double weight = ip.weight * Trans.Weight(); + + el.CalcShape(ip, shape); + + double u_val = 0.0; + for (int j = 0; j < dof; j++) { + u_val += elfun(j) * shape(j); + } + const double u_safe = std::max(u_val, 0.0); + const double u_nl = std::pow(u_safe, m_polytropicIndex); + + const double x2_u_nl = u_nl; + + for (int i = 0; i < dof; i++){ + elvect(i) += shape(i) * x2_u_nl * weight; + } + } + } + + + void NonlinearPowerIntegrator::AssembleElementGrad ( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::DenseMatrix &elmat) { + + const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); + const int dof = el.GetDof(); + elmat.SetSize(dof); + elmat = 0.0; + mfem::Vector shape(dof); + mfem::DenseMatrix dshape(dof, 3); + mfem::DenseMatrix invJ(3, 3); + mfem::Vector gradPhys(3); + mfem::Vector physCoord(3); + + for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { + const mfem::IntegrationPoint &ip = ir->IntPoint(iqp); + Trans.SetIntPoint(&ip); + const double weight = ip.weight * Trans.Weight(); + + + el.CalcShape(ip, shape); + + double u_val = 0.0; + + for (int j = 0; j < dof; j++) { + u_val += elfun(j) * shape(j); + } + + // Calculate the Jacobian + + // TODO: Check for when theta ~ 0? + const double u_safe = std::max(u_val, 0.0); + const double d_u_nl = m_polytropicIndex * std::pow(u_safe, m_polytropicIndex - 1); + const double x2_d_u_nl = d_u_nl; + + for (int i = 0; i < dof; i++) { + for (int j = 0; j < dof; j++) { + elmat(i, j) += shape(i) * x2_d_u_nl * shape(j) * weight; + } + } + + // // --- Debug Code to write out gradient --- + // Trans.Transform(ip,physCoord); + // el.CalcDShape(ip, dshape); + // + // mfem::CalcInverse(Trans.Jacobian(), invJ); + // + // mfem::DenseMatrix dshapePhys; + // dshapePhys.SetSize(dof, physCoord.Size()); + // mfem::Mult(dshape, invJ, dshapePhys); + // + // gradPhys = 0.0; + // for (int j = 0; j < dof; j++) { + // for (int d = 0; d < gradPhys.Size(); d++) { + // gradPhys(d) += elfun(j)*dshapePhys(j, d); + // } + // } + // + // debugOut + // << physCoord(0) << ", " << physCoord(1) << ", " << physCoord(2) + // << ", " << gradPhys(0) << ", " << gradPhys(1) << ", " << gradPhys(2) << '\n'; + } + // debugOut.flush(); + } +} // namespace polyMFEMUtils \ No newline at end of file diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp new file mode 100644 index 0000000..ea363d5 --- /dev/null +++ b/src/poly/utils/private/operator.cpp @@ -0,0 +1,336 @@ +/* *********************************************************************** +// +// Copyright (C) 2025 -- The 4D-STAR Collaboration +// File Author: Emily Boudreaux +// Last Modified: April 21, 2025 +// +// 4DSSE is free software; you can use it and/or modify +// it under the terms and restrictions the GNU General Library Public +// License version 3 (GPLv3) as published by the Free Software Foundation. +// +// 4DSSE is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU Library General Public License for more details. +// +// You should have received a copy of the GNU Library General Public License +// along with this software; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// *********************************************************************** */ +#include "operator.h" +#include "4DSTARTypes.h" +#include "mfem.hpp" +#include + +static int s_PolyOperatorMultCalledCount = 0; + +void writeDenseMatrixToCSV(const std::string &filename, int precision, const mfem::DenseMatrix *mat) { + if (!mat) { + throw std::runtime_error("The operator is not a SparseMatrix."); + } + + std::ofstream outfile(filename); + if (!outfile.is_open()) { + throw std::runtime_error("Failed to open file: " + filename); + } + + + int height = mat->Height(); + int width = mat->Width(); + + // Set precision for floating-point output + outfile << std::fixed << std::setprecision(precision); + + for (int i = 0; i < width; i++) { + outfile << i; + if (i < width - 1) { + outfile << ","; + } + else { + outfile << "\n"; + } + } + + // Iterate through rows + for (int i = 0; i < height; ++i) { + for (int j = 0; j < width; ++j) { + outfile << mat->Elem(i, j); + if (j < width - 1) { + outfile << ","; + } + } + outfile << std::endl; + } + + outfile.close(); +} + +/** + * @brief Writes the dense representation of an MFEM Operator (if it's a SparseMatrix) to a CSV file. + * + * @param op The MFEM Operator to write. + * @param filename The name of the output CSV file. + * @param precision Number of decimal places for floating-point values. + */ + void writeOperatorToCSV(const mfem::Operator &op, + const std::string &filename, + const int precision = 6) // Add precision argument +{ + // Attempt to cast the Operator to a SparseMatrix + const auto *sparse_mat = dynamic_cast(&op); + if (!sparse_mat) { + throw std::runtime_error("The operator is not a SparseMatrix."); + } + const mfem::DenseMatrix *mat = sparse_mat->ToDenseMatrix(); + writeDenseMatrixToCSV(filename, precision, mat); +} + +void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr& invMat, const std::string& name="matrix") { + // PERF: This likely can be made much more efficient and will probably be called in tight loops, a good + // PERF: place for some easy optimization might be here. + + // Confirm that mat is a square matrix + MFEM_ASSERT(mat.Height() == mat.Width(), "Matrix " + name + " is not square, cannot invert."); + + mfem::Vector diag; + mat.GetDiag(diag); + + // Invert the diagonal + for (int i = 0; i < diag.Size(); i++) { + MFEM_ASSERT(diag(i) != 0, "Diagonal element (" + std::to_string(i) +") in " + name + " is zero, cannot invert."); + diag(i) = 1.0 / diag(i); + } + + // If the matrix is already inverted, just set the diagonal to avoid reallocation + if (invMat != nullptr) { + MFEM_ASSERT(invMat->Height() == invMat->Width(), "invMat (result matrix) is not square, cannot invert " + name + " into it."); + MFEM_ASSERT(invMat->Height() == mat.Height(), "Incompatible matrix sizes for inversion of " + name + ", expected " + std::to_string(mat.Height()) + " but got " + std::to_string(invMat->Height())); + for (int i = 0; i < diag.Size(); i++) { + MFEM_ASSERT(diag(i) != 0, "Diagonal element (" + std::to_string(i) +") in " + name + " is zero, resulting matrix would be singular."); + invMat->Elem(i, i) = diag(i); + } + } else { // The matrix has not been allocated yet so that needs to be done. Sparse Matrix has a constructor that can build from the diagonals + invMat = std::make_unique(diag); + } +} + +PolytropeOperator::PolytropeOperator( + + std::unique_ptr M, + std::unique_ptr Q, + std::unique_ptr D, + std::unique_ptr f, + const mfem::Array &blockOffsets) : + + mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector + m_blockOffsets(blockOffsets), + m_jacobian(nullptr) { + + m_M = std::move(M); + m_Q = std::move(Q); + m_D = std::move(D); + m_f = std::move(f); +} + +void PolytropeOperator::finalize(const mfem::Vector &initTheta) { + if (m_isFinalized) { + return; // do nothing if already finalized + } + + m_negM_op = std::make_unique(m_M.get(), -1.0); + m_negQ_op = std::make_unique(m_Q.get(), -1.0); + + // Set up the constant parts of the jacobian now + m_jacobian = std::make_unique(m_blockOffsets); + m_jacobian->SetBlock(0, 1, m_negM_op.get()); //<- -M (constant) + m_jacobian->SetBlock(1, 0, m_negQ_op.get()); //<- -Q (constant) + m_jacobian->SetBlock(1, 1, m_D.get()); //<- D (constant) + + // Build the initial preconditioner based on some initial guess + const auto &grad = m_f->GetGradient(initTheta); + updatePreconditioner(grad); + + m_isFinalized = true; +} + +const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const { + if (m_jacobian == nullptr) { + MFEM_ABORT("Jacobian has not been initialized before GetJacobianOperator() call."); + } + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator not finalized prior to call to GetJacobianOperator()."); + } + + return *m_jacobian; +} + +mfem::BlockDiagonalPreconditioner& PolytropeOperator::GetPreconditioner() const { + if (m_schurPreconditioner == nullptr) { + MFEM_ABORT("Schur preconditioner has not been initialized before GetPreconditioner() call."); + } + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator not finalized prior to call to GetPreconditioner()."); + } + return *m_schurPreconditioner; +} + +void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator::Mult called before finalize"); + } + // -- Create BlockVector views for input x and output y + mfem::BlockVector x_block(const_cast(x), m_blockOffsets); + mfem::BlockVector y_block(y, m_blockOffsets); + + // -- Get Vector views for individual blocks + const mfem::Vector &x_theta = x_block.GetBlock(0); + const mfem::Vector &x_phi = x_block.GetBlock(1); + mfem::Vector &y_R0 = y_block.GetBlock(0); // Residual Block 0 (theta) + mfem::Vector &y_R1 = y_block.GetBlock(1); // Residual Block 1 (phi) + + int theta_size = m_blockOffsets[1] - m_blockOffsets[0]; + int phi_size = m_blockOffsets[2] - m_blockOffsets[1]; + + mfem::Vector f_term(theta_size); + mfem::Vector Mphi_term(theta_size); + mfem::Vector Dphi_term(phi_size); + mfem::Vector Qtheta_term(phi_size); + + // Calculate R0 and R1 terms + // R0 = f(θ) - Mɸ + // R1 = Dɸ - Qθ + + MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult"); + + m_f->Mult(x_theta, f_term); // fixme: this may be the wrong way to assemble m_f? + m_M->Mult(x_phi, Mphi_term); + m_D->Mult(x_phi, Dphi_term); + m_Q->Mult(x_theta, Qtheta_term); + + subtract(f_term, Mphi_term, y_R0); + subtract(Dphi_term, Qtheta_term, y_R1); + + // -- Apply essential boundary conditions -- + for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) { + if (int idx = m_theta_ess_tdofs.first[i]; idx >= 0 && idx < y_R0.Size()) { + const double &targetValue = m_theta_ess_tdofs.second[i]; + // y_block.GetBlock(0)[idx] = targetValue - x_theta(idx); // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising + y_block.GetBlock(0)[idx] = 0; // Zero out the essential theta dofs in the bi-linear form + } + } + + // TODO: Are the residuals for φ being calculated correctly? + + std::ofstream outfileθ("PolyOperatorMultCallTheta_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc); + outfileθ << "dof,R_θ\n"; + for (int i = 0; i < y_R0.Size(); i++) { + outfileθ << i << "," << y_R0(i) << "\n"; + } + outfileθ.close(); + std::ofstream outfileφ("PolyOperatorMultCallPhi_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc); + outfileφ << "dof,R_ɸ\n"; + for (int i = 0; i < y_R1.Size(); i++) { + outfileφ << i << "," << y_R1(i) << "\n"; + } + outfileφ.close(); + ++s_PolyOperatorMultCalledCount; + +} + + +void PolytropeOperator::updateInverseNonlinearJacobian(const mfem::Operator &grad) const { + if (const auto *sparse_mat = dynamic_cast(&grad); sparse_mat != nullptr) { + approxJacobiInvert(*sparse_mat, m_invNonlinearJacobian, "Nonlinear Jacobian"); + } else { + MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear jacobian where nonlinear jacobian is not dynamically castable to a sparse matrix"); + } +} + +void PolytropeOperator::updateInverseSchurCompliment() const { + // TODO Add a flag in to make sure this tracks in parallel (i.e. every time the non linear jacobian inverse is updated set the flag to true and then check if the flag is true here and if so do work (if not throw error). then at the end of this function set it to false. + if (m_invNonlinearJacobian == nullptr) { + MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseNonlinearJacobian"); + } + mfem::SparseMatrix* schurCompliment(&m_D->SpMat()); // This is now a copy of D + + // Calculate S = D - Q df^{-1} M + mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(m_Q->SpMat(), *m_invNonlinearJacobian); // Q * df^{-1} + const mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, m_M->SpMat()); // Q * df^{-1} * M + + // PERF: Could potentially add some caching in here to only update the preconditioner when some condition has been met + schurCompliment->Add(-1, *temp_QxdfInvxM); // D - Q * df^{-1} * M + approxJacobiInvert(*schurCompliment, m_invSchurCompliment, "Schur Compliment"); + + if (m_schurPreconditioner == nullptr) { + m_schurPreconditioner = std::make_unique(m_blockOffsets); + } + + // ⎡ḟ(θ)^-1 0⎤ + // ⎣0 S^-1⎦ + m_schurPreconditioner->SetDiagonalBlock(0, m_invNonlinearJacobian.get()); + m_schurPreconditioner->SetDiagonalBlock(1, m_invSchurCompliment.get()); +} + +void PolytropeOperator::updatePreconditioner(const mfem::Operator &grad) const { + updateInverseNonlinearJacobian(grad); + updateInverseSchurCompliment(); +} + + +mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const { + if (!m_isFinalized) { + MFEM_ABORT("PolytropeOperator::GetGradient called before finalize"); + } + // --- Get the gradient of f --- + mfem::BlockVector x_block(const_cast(x), m_blockOffsets); + const mfem::Vector& x_theta = x_block.GetBlock(0); + + auto &grad = m_f->GetGradient(x_theta); + updatePreconditioner(grad); + + m_jacobian->SetBlock(0, 0, &grad); + // The other blocks are set up in finalize since they are constant. Only J00 depends on the current state of theta + + return *m_jacobian; +} + +void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs) { + m_isFinalized = false; + m_theta_ess_tdofs = theta_ess_tdofs; + m_phi_ess_tdofs = phi_ess_tdofs; + + if (m_f) { + m_f->SetEssentialTrueDofs(theta_ess_tdofs.first); + } else { + MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs"); + } + + if (m_M) { + m_M->EliminateTestEssentialBC(theta_ess_tdofs.first); + m_M->EliminateTrialEssentialBC(phi_ess_tdofs.first); + } else { + MFEM_ABORT("m_M is null in PolytropeOperator::SetEssentialTrueDofs"); + } + + if (m_Q) { + m_Q->EliminateTrialEssentialBC(theta_ess_tdofs.first); + m_Q->EliminateTestEssentialBC(phi_ess_tdofs.first); + } else { + MFEM_ABORT("m_Q is null in PolytropeOperator::SetEssentialTrueDofs"); + } + + if (m_D) { + m_D->EliminateEssentialBC(phi_ess_tdofs.first); + } else { + MFEM_ABORT("m_D is null in PolytropeOperator::SetEssentialTrueDofs"); + } +} + +void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set) { + SetEssentialTrueDofs(ess_tdof_pair_set.first, ess_tdof_pair_set.second); +} + +SSE::MFEMArrayPairSet PolytropeOperator::GetEssentialTrueDofs() const { + return std::make_pair(m_theta_ess_tdofs, m_phi_ess_tdofs); +} \ No newline at end of file diff --git a/src/poly/utils/private/polyIO.cpp b/src/poly/utils/private/polyIO.cpp deleted file mode 100644 index 007cda8..0000000 --- a/src/poly/utils/private/polyIO.cpp +++ /dev/null @@ -1,43 +0,0 @@ -/* *********************************************************************** -// -// Copyright (C) 2025 -- The 4D-STAR Collaboration -// File Author: Emily Boudreaux -// Last Modified: February 14, 2025 -// -// 4DSSE is free software; you can use it and/or modify -// it under the terms and restrictions the GNU General Library Public -// License version 3 (GPLv3) as published by the Free Software Foundation. -// -// 4DSSE is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -// See the GNU Library General Public License for more details. -// -// You should have received a copy of the GNU Library General Public License -// along with this software; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// *********************************************************************** */ -#include "mfem.hpp" -#include -#include - -#include "polyIO.h" - -void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename) { - std::ofstream file(filename); - if (!file.is_open()) { - std::cerr << "Error: Could not open " << filename << " for writing." << std::endl; - return; - } - - file << "xi,u\n"; // CSV header - - for (int i = 0; i < u.Size(); i++) { - double xi = mesh.GetVertex(i)[0]; // Get spatial coordinate - file << xi << "," << u[i] << "\n"; // Write to CSV - } - - file.close(); - std::cout << "Solution written to " << filename << std::endl; -} \ No newline at end of file diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp deleted file mode 100644 index 65636e2..0000000 --- a/src/poly/utils/private/polyMFEMUtils.cpp +++ /dev/null @@ -1,195 +0,0 @@ -/* *********************************************************************** -// -// Copyright (C) 2025 -- The 4D-STAR Collaboration -// File Author: Emily Boudreaux -// Last Modified: February 12, 2025 -// -// 4DSSE is free software; you can use it and/or modify -// it under the terms and restrictions the GNU General Library Public -// License version 3 (GPLv3) as published by the Free Software Foundation. -// -// 4DSSE is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -// See the GNU Library General Public License for more details. -// -// You should have received a copy of the GNU Library General Public License -// along with this software; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// *********************************************************************** */ -#include "mfem.hpp" -#include -#include -#include - -#include "polyMFEMUtils.h" - -NonlinearPowerIntegrator::NonlinearPowerIntegrator( - mfem::FunctionCoefficient &coeff, - double n) : coeff_(coeff), polytropicIndex(n) { - -} - -void NonlinearPowerIntegrator::AssembleElementVector( - const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::Vector &elvect) { - - const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); - int dof = el.GetDof(); - elvect.SetSize(dof); - elvect = 0.0; - - mfem::Vector shape(dof); - - for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { - mfem::IntegrationPoint ip = ir->IntPoint(iqp); - Trans.SetIntPoint(&ip); - double weight = ip.weight * Trans.Weight(); - - el.CalcShape(ip, shape); - - double u_val = 0.0; - for (int j = 0; j < dof; j++) { - u_val += elfun(j) * shape(j); - } - double u_safe = std::max(u_val, 0.0); - double u_nl = std::pow(u_safe, polytropicIndex); - - double coeff_val = coeff_.Eval(Trans, ip); - double x2_u_nl = coeff_val * u_nl; - - for (int i = 0; i < dof; i++){ - elvect(i) += shape(i) * x2_u_nl * weight; - } - } -} - -void NonlinearPowerIntegrator::AssembleElementGrad ( - const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::DenseMatrix &elmat) { - - const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); - int dof = el.GetDof(); - elmat.SetSize(dof); - elmat = 0.0; - mfem::Vector shape(dof); - - for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { - mfem::IntegrationPoint ip = ir->IntPoint(iqp); - Trans.SetIntPoint(&ip); - double weight = ip.weight * Trans.Weight(); - - el.CalcShape(ip, shape); - - double u_val = 0.0; - - for (int j = 0; j < dof; j++) { - u_val += elfun(j) * shape(j); - } - double coeff_val = coeff_.Eval(Trans, ip); - - - // Calculate the Jacobian - double u_safe = std::max(u_val, 0.0); - double d_u_nl = coeff_val * polytropicIndex * std::pow(u_safe, polytropicIndex - 1); - double x2_d_u_nl = d_u_nl; - - for (int i = 0; i < dof; i++) { - for (int j = 0; j < dof; j++) { - elmat(i, j) += shape(i) * x2_d_u_nl * shape(j) * weight; - } - } - - } -} - -BilinearIntegratorWrapper::BilinearIntegratorWrapper( - mfem::BilinearFormIntegrator *integratorInput - ) : integrator(integratorInput) { } - -BilinearIntegratorWrapper::~BilinearIntegratorWrapper() { - delete integrator; -} - -void BilinearIntegratorWrapper::AssembleElementVector( - const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::Vector &elvect) { - int dof = el.GetDof(); - mfem::DenseMatrix elMat(dof); - integrator->AssembleElementMatrix(el, Trans, elMat); - elvect.SetSize(dof); - elvect = 0.0; - for (int i = 0; i < dof; i++) - { - double sum = 0.0; - for (int j = 0; j < dof; j++) - { - sum += elMat(i, j) * elfun(j); - } - elvect(i) = sum; - } -} - -void BilinearIntegratorWrapper::AssembleElementGrad(const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::DenseMatrix &elmat) { - int dof = el.GetDof(); - elmat.SetSize(dof, dof); - elmat = 0.0; - integrator->AssembleElementMatrix(el, Trans, elmat); -} - -CompositeNonlinearIntegrator::CompositeNonlinearIntegrator() { } - - -CompositeNonlinearIntegrator::~CompositeNonlinearIntegrator() { - for (size_t i = 0; i < integrators.size(); i++) { - delete integrators[i]; - } -} - -void CompositeNonlinearIntegrator::add_integrator(mfem::NonlinearFormIntegrator *integrator) { - integrators.push_back(integrator); -} - -void CompositeNonlinearIntegrator::AssembleElementVector( - const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::Vector &elvect) { - int dof = el.GetDof(); - elvect.SetSize(dof); - elvect = 0.0; - mfem::Vector temp(dof); - - for (size_t i = 0; i < integrators.size(); i++) { - temp= 0.0; - integrators[i]->AssembleElementVector(el, Trans, elfun, temp); - elvect.Add(1.0, temp); - } -} - -void CompositeNonlinearIntegrator::AssembleElementGrad( - const mfem::FiniteElement &el, - mfem::ElementTransformation &Trans, - const mfem::Vector &elfun, - mfem::DenseMatrix &elmat) { - int dof = el.GetDof(); - elmat.SetSize(dof, dof); - elmat = 0.0; - mfem::DenseMatrix temp(dof); - temp.SetSize(dof, dof); - for (size_t i = 0; i < integrators.size(); i++) { - temp = 0.0; - integrators[i] -> AssembleElementGrad(el, Trans, elfun, temp); - elmat.Add(1.0, temp); - } -} \ No newline at end of file diff --git a/src/poly/utils/public/integrators.h b/src/poly/utils/public/integrators.h new file mode 100644 index 0000000..8a625d6 --- /dev/null +++ b/src/poly/utils/public/integrators.h @@ -0,0 +1,77 @@ +/* *********************************************************************** +// +// Copyright (C) 2025 -- The 4D-STAR Collaboration +// File Author: Emily Boudreaux +// Last Modified: April 21, 2025 +// +// 4DSSE is free software; you can use it and/or modify +// it under the terms and restrictions the GNU General Library Public +// License version 3 (GPLv3) as published by the Free Software Foundation. +// +// 4DSSE is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU Library General Public License for more details. +// +// You should have received a copy of the GNU Library General Public License +// along with this software; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// *********************************************************************** */ +#pragma once + +#include "mfem.hpp" +#include +#include "config.h" +#include "probe.h" + + +/** + * @file integrators.h + * @brief A collection of utilities for working with MFEM and solving the lane-emden equation. + */ + + +/** + * @namespace polyMFEMUtils + * @brief A namespace for utilities for working with MFEM and solving the lane-emden equation. + */ +namespace polyMFEMUtils { + /** + * @brief A class for nonlinear power integrator. + */ + class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator { + public: + /** + * @brief Constructor for NonlinearPowerIntegrator. + * + * @param coeff The function coefficient. + * @param n The polytropic index. + */ + NonlinearPowerIntegrator(double n); + + /** + * @brief Assembles the element vector. + * + * @param el The finite element. + * @param Trans The element transformation. + * @param elfun The element function. + * @param elvect The element vector to be assembled. + */ + virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; + /** + * @brief Assembles the element gradient. + * + * @param el The finite element. + * @param Trans The element transformation. + * @param elfun The element function. + * @param elmat The element matrix to be assembled. + */ + virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; + private: + Config& m_config = Config::getInstance(); + Probe::LogManager& m_logManager = Probe::LogManager::getInstance(); + quill::Logger* m_logger = m_logManager.getLogger("log"); + double m_polytropicIndex; + }; +} // namespace polyMFEMUtils \ No newline at end of file diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h new file mode 100644 index 0000000..1c84db0 --- /dev/null +++ b/src/poly/utils/public/operator.h @@ -0,0 +1,98 @@ +/* *********************************************************************** +// +// Copyright (C) 2025 -- The 4D-STAR Collaboration +// File Author: Emily Boudreaux +// Last Modified: April 21, 2025 +// +// 4DSSE is free software; you can use it and/or modify +// it under the terms and restrictions the GNU General Library Public +// License version 3 (GPLv3) as published by the Free Software Foundation. +// +// 4DSSE is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU Library General Public License for more details. +// +// You should have received a copy of the GNU Library General Public License +// along with this software; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// *********************************************************************** */ +#pragma once + +#include "mfem.hpp" +#include "4DSTARTypes.h" +#include + +#include "probe.h" + +class PolytropeOperator final : public mfem::Operator { +public: + PolytropeOperator( + std::unique_ptr M, + std::unique_ptr Q, + std::unique_ptr D, + std::unique_ptr f, + const mfem::Array &blockOffsets); + ~PolytropeOperator() override = default; + + void Mult(const mfem::Vector &x, mfem::Vector &y) const override; + + + mfem::Operator& GetGradient(const mfem::Vector &x) const override; + + void SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs); + void SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set); + + SSE::MFEMArrayPairSet GetEssentialTrueDofs() const; + + bool isFinalized() const { return m_isFinalized; } + + void finalize(const mfem::Vector &initTheta); + + const mfem::Array& GetBlockOffsets() const { return m_blockOffsets; } + + const mfem::BlockOperator &GetJacobianOperator() const; + + mfem::BlockDiagonalPreconditioner &GetPreconditioner() const; + + +private: + Probe::LogManager& m_logManager = Probe::LogManager::getInstance(); + quill::Logger* m_logger = m_logManager.getLogger("log"); + std::unique_ptr m_M; + std::unique_ptr m_Q; + std::unique_ptr m_D; + std::unique_ptr m_f; + + const mfem::Array m_blockOffsets; + + SSE::MFEMArrayPair m_theta_ess_tdofs; + SSE::MFEMArrayPair m_phi_ess_tdofs; + + std::unique_ptr m_negM_op; + std::unique_ptr m_negQ_op; + mutable std::unique_ptr m_jacobian; + + mutable std::unique_ptr m_invSchurCompliment; + mutable std::unique_ptr m_invNonlinearJacobian; + + /* + * The schur preconditioner has the form + * + * ⎡ḟ(θ)^-1 0 ⎤ + * ⎣ 0 S^-1 ⎦ + * + * Where S is the Schur compliment of the system + * + */ + + mutable std::unique_ptr m_schurPreconditioner; + + bool m_isFinalized = false; + +private: + void updateInverseNonlinearJacobian(const mfem::Operator &grad) const; + void updateInverseSchurCompliment() const; + void updatePreconditioner(const mfem::Operator &grad) const; +}; diff --git a/src/poly/utils/public/polyIO.h b/src/poly/utils/public/polyIO.h deleted file mode 100644 index acf2237..0000000 --- a/src/poly/utils/public/polyIO.h +++ /dev/null @@ -1,36 +0,0 @@ -/* *********************************************************************** -// -// Copyright (C) 2025 -- The 4D-STAR Collaboration -// File Author: Emily Boudreaux -// Last Modified: February 14, 2025 -// -// 4DSSE is free software; you can use it and/or modify -// it under the terms and restrictions the GNU General Library Public -// License version 3 (GPLv3) as published by the Free Software Foundation. -// -// 4DSSE is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -// See the GNU Library General Public License for more details. -// -// You should have received a copy of the GNU Library General Public License -// along with this software; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// *********************************************************************** */ -#ifndef POLY_IO_H -#define POLY_IO_H - -#include "mfem.hpp" -#include - -/** - * @brief Writes the solution to a CSV file. - * - * @param u The GridFunction containing the solution. - * @param mesh The mesh associated with the solution. - * @param filename The name of the CSV file to write to. - */ -void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename); - -#endif // POLY_IO_H \ No newline at end of file diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h deleted file mode 100644 index d48a78d..0000000 --- a/src/poly/utils/public/polyMFEMUtils.h +++ /dev/null @@ -1,148 +0,0 @@ -/* *********************************************************************** -// -// Copyright (C) 2025 -- The 4D-STAR Collaboration -// File Author: Emily Boudreaux -// Last Modified: February 14, 2025 -// -// 4DSSE is free software; you can use it and/or modify -// it under the terms and restrictions the GNU General Library Public -// License version 3 (GPLv3) as published by the Free Software Foundation. -// -// 4DSSE is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. -// See the GNU Library General Public License for more details. -// -// You should have received a copy of the GNU Library General Public License -// along with this software; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// *********************************************************************** */ -#include "mfem.hpp" -#include - - -void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename); - -/** - * @brief A class for nonlinear power integrator. - */ -class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator { -private: - mfem::FunctionCoefficient coeff_; - double polytropicIndex; -public: - /** - * @brief Constructor for NonlinearPowerIntegrator. - * - * @param coeff The function coefficient. - * @param n The polytropic index. - */ - NonlinearPowerIntegrator(mfem::FunctionCoefficient &coeff, double n); - - /** - * @brief Assembles the element vector. - * - * @param el The finite element. - * @param Trans The element transformation. - * @param elfun The element function. - * @param elvect The element vector to be assembled. - */ - virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; - - /** - * @brief Assembles the element gradient. - * - * @param el The finite element. - * @param Trans The element transformation. - * @param elfun The element function. - * @param elmat The element matrix to be assembled. - */ - virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; -}; - -/** - * @brief A wrapper class for bilinear integrator. - */ -class BilinearIntegratorWrapper : public mfem::NonlinearFormIntegrator -{ -private: - mfem::BilinearFormIntegrator *integrator; -public: - /** - * @brief Constructor for BilinearIntegratorWrapper. - * - * @param integratorInput The bilinear form integrator input. - */ - BilinearIntegratorWrapper(mfem::BilinearFormIntegrator *integratorInput); - - /** - * @brief Destructor for BilinearIntegratorWrapper. - */ - virtual ~BilinearIntegratorWrapper(); - - /** - * @brief Assembles the element vector. - * - * @param el The finite element. - * @param Trans The element transformation. - * @param elfun The element function. - * @param elvect The element vector to be assembled. - */ - virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; - - /** - * @brief Assembles the element gradient. - * - * @param el The finite element. - * @param Trans The element transformation. - * @param elfun The element function. - * @param elmat The element matrix to be assembled. - */ - virtual void AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; -}; - -/** - * @brief A class for composite nonlinear integrator. - */ -class CompositeNonlinearIntegrator: public mfem::NonlinearFormIntegrator { - private: - std::vector integrators; - public: - /** - * @brief Constructor for CompositeNonlinearIntegrator. - */ - CompositeNonlinearIntegrator(); - - /** - * @brief Destructor for CompositeNonlinearIntegrator. - */ - virtual ~CompositeNonlinearIntegrator(); - - /** - * @brief Adds an integrator to the composite integrator. - * - * @param integrator The nonlinear form integrator to add. - */ - void add_integrator(mfem::NonlinearFormIntegrator *integrator); - - /** - * @brief Assembles the element vector. - * - * @param el The finite element. - * @param Trans The element transformation. - * @param elfun The element function. - * @param elvect The element vector to be assembled. - */ - virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; - - /** - * @brief Assembles the element gradient. - * - * @param el The finite element. - * @param Trans The element transformation. - * @param elfun The element function. - * @param elmat The element matrix to be assembled. - */ - virtual void AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; -}; \ No newline at end of file diff --git a/src/types/meson.build b/src/types/meson.build new file mode 100644 index 0000000..5b71ae8 --- /dev/null +++ b/src/types/meson.build @@ -0,0 +1 @@ +types_dep = declare_dependency(include_directories: include_directories('public')) \ No newline at end of file diff --git a/src/poly/coeff/public/polyCoeff.h b/src/types/public/4DSTARTypes.h similarity index 73% rename from src/poly/coeff/public/polyCoeff.h rename to src/types/public/4DSTARTypes.h index 240b6af..687d0a1 100644 --- a/src/poly/coeff/public/polyCoeff.h +++ b/src/types/public/4DSTARTypes.h @@ -2,7 +2,7 @@ // // Copyright (C) 2025 -- The 4D-STAR Collaboration // File Author: Emily Boudreaux -// Last Modified: February 12, 2025 +// Last Modified: April 09, 2025 // // 4DSSE is free software; you can use it and/or modify // it under the terms and restrictions the GNU General Library Public @@ -18,11 +18,16 @@ // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // // *********************************************************************** */ +#ifndef _4DSTAR_TYPES_H +#define _4DSTAR_TYPES_H + +#include #include "mfem.hpp" -#include -double xi_coeff_func(const mfem::Vector &x); +// TODO : Need a better namespace name for these types +namespace SSE { + typedef std::pair, mfem::Array> MFEMArrayPair; + typedef std::pair MFEMArrayPairSet; +} -void vec_xi_coeff_func(const mfem::Vector &x, mfem::Vector &v); - -double theta_initial_guess(const mfem::Vector &x, double root); \ No newline at end of file +#endif // _4DSTAR_TYPES_H \ No newline at end of file