diff --git a/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.pdf b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.pdf index 092993d..05ea41c 100644 Binary files a/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.pdf and b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.pdf differ diff --git a/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.tex b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.tex index 6604047..b0ad19f 100644 --- a/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.tex +++ b/docs/derivations/laneEmdenVariationalForm/laneEmdenVariationalBlockForm.tex @@ -84,7 +84,7 @@ Now we exploit the linearity of summation and integration to move the sums out o \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 \\ + 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} @@ -96,18 +96,18 @@ f(\bar{\theta}) \equiv \int_{\Omega}\psi_{k}^{\theta}\left(\theta_{h}\right)^{n} \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_{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{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} \\ + 0 & \mathbf{M} \\ -\mathbf{Q} & \mathbf{D} \end{bmatrix} \begin{bmatrix} @@ -126,7 +126,7 @@ From this form we can easily see that the residual matrix is \begin{align} R &= \begin{bmatrix} - f(\bar{\theta}) - M\bar{\phi} \\ + f(\bar{\theta}) + M\bar{\phi} \\ D\bar{\phi} - Q\bar{\theta} \end{bmatrix} \end{align} @@ -144,11 +144,11 @@ in our Newton-Raphson method. Generally the Jacobian is the matrix of partial de 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(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} \\ + \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} @@ -156,7 +156,7 @@ So then the Jacobian is 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 \\ + \frac{df}{d\theta} & M \\ -Q & D \end{bmatrix} \end{align} @@ -164,9 +164,9 @@ Finally, we know that the matrices $M$, $D$, and $Q$ are constant with respect t \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}} \\ + \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}} \\ + \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}} diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/polytropeOperator.cpp similarity index 92% rename from src/poly/utils/private/operator.cpp rename to src/poly/utils/private/polytropeOperator.cpp index db092a0..e81a032 100644 --- a/src/poly/utils/private/operator.cpp +++ b/src/poly/utils/private/polytropeOperator.cpp @@ -58,38 +58,50 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) { m_Qmat = std::make_unique(m_Q->SpMat()); m_Dmat = std::make_unique(m_D->SpMat()); + saveSparseMatrixBinary(*m_Mmat, "MmatRawO2.bin"); + saveSparseMatrixBinary(*m_Qmat, "QmatRawO2.bin"); + saveSparseMatrixBinary(*m_Dmat, "DmatRawO2.bin"); + // Remove the essential dofs from the constant matrices for (const auto& dof: m_theta_ess_tdofs.first) { + std::cout << "Eliminating dof: " << dof << std::endl; m_Mmat->EliminateRow(dof); m_Qmat->EliminateCol(dof); } + saveSparseMatrixBinary(*m_Mmat, "MmatBCThetaO2.bin"); + saveSparseMatrixBinary(*m_Qmat, "QmatBCThetaO2.bin"); + saveSparseMatrixBinary(*m_Dmat, "DmatBCThetaO2.bin"); + for (const auto& dof: m_phi_ess_tdofs.first) { m_Mmat->EliminateCol(dof); m_Qmat->EliminateRow(dof); m_Dmat->EliminateRowCol(dof); } + saveSparseMatrixBinary(*m_Mmat, "MmatBCO2.bin"); + saveSparseMatrixBinary(*m_Qmat, "QmatBCO2.bin"); + saveSparseMatrixBinary(*m_Dmat, "DmatBCO2.bin"); m_negM_mat = std::make_unique(m_Mmat.get(), -1.0); m_negQ_mat = std::make_unique(m_Qmat.get(), -1.0); - m_schurCompliment = std::make_unique(*m_Qmat, *m_Dmat, *m_Mmat); + // m_schurCompliment = std::make_unique(*m_Qmat, *m_Dmat, *m_Mmat); // Set up the constant parts of the jacobian now // We use the assembled matrices with their boundary conditions already removed for the jacobian m_jacobian = std::make_unique(m_blockOffsets); - m_jacobian->SetBlock(0, 1, m_negM_mat.get()); //<- -M (constant) + m_jacobian->SetBlock(0, 1, m_Mmat.get()); //<- M (constant) m_jacobian->SetBlock(1, 0, m_negQ_mat.get()); //<- -Q (constant) m_jacobian->SetBlock(1, 1, m_Dmat.get()); //<- D (constant) - m_invSchurCompliment = std::make_unique(*m_schurCompliment); + // m_invSchurCompliment = std::make_unique(*m_schurCompliment); m_isFinalized = true; // Build the initial preconditioner based on some initial guess const auto &grad = m_f->GetGradient(initTheta); - updatePreconditioner(grad); + // updatePreconditioner(grad); } @@ -137,7 +149,7 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { mfem::Vector Qtheta_term(phi_size); // Calculate R0 and R1 terms - // R0 = f(θ) - Mɸ + // R0 = f(θ) + Mɸ // R1 = Dɸ - Qθ MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult"); @@ -147,7 +159,7 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { m_D->Mult(x_phi, Dphi_term); m_Q->Mult(x_theta, Qtheta_term); - subtract(f_term, Mphi_term, y_R0); + add(f_term, Mphi_term, y_R0); subtract(Dphi_term, Qtheta_term, y_R1); // -- Apply essential boundary conditions -- @@ -211,7 +223,7 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const { const mfem::Vector& x_theta = x_block.GetBlock(0); auto &grad = m_f->GetGradient(x_theta); - updatePreconditioner(grad); + // updatePreconditioner(grad); m_jacobian->SetBlock(0, 0, &grad); diff --git a/src/poly/utils/private/utilities.cpp b/src/poly/utils/private/utilities.cpp new file mode 100644 index 0000000..227b549 --- /dev/null +++ b/src/poly/utils/private/utilities.cpp @@ -0,0 +1,5 @@ +// +// Created by Emily Boudreaux on 6/4/25. +// + +#include "utilities.h" diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/polytropeOperator.h similarity index 100% rename from src/poly/utils/public/operator.h rename to src/poly/utils/public/polytropeOperator.h