diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 6134560..26aa293 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -106,74 +106,65 @@ PolySolver::PolySolver(const double n, const double order) { PolySolver::~PolySolver() = default; void PolySolver::assembleBlockSystem() { + mfem::Array blockOffsets = computeBlockOffsets(); - // Start by defining the block structure of the system - // Block 0: Theta (scalar space, uses m_feTheta) - // Block 1: Phi (vector space, uses m_fePhi) - mfem::Array feSpaces; - feSpaces.Append(m_feTheta.get()); - feSpaces.Append(m_fePhi.get()); + const std::unique_ptr forms = buildIndividualForms(blockOffsets); - // Create the block offsets. These define the start of each block in the combined vector. - // Block offsets will be [0, thetaDofs, thetaDofs + phiDofs] + // --- Build the BlockOperator --- + m_polytropOperator = std::make_unique( + std::move(forms->M), + std::move(forms->Q), + std::move(forms->D), + std::move(forms->f), + blockOffsets); +} + +mfem::Array PolySolver::computeBlockOffsets() const { mfem::Array blockOffsets; blockOffsets.SetSize(3); blockOffsets[0] = 0; - blockOffsets[1] = feSpaces[0]->GetVSize(); - blockOffsets[2] = feSpaces[1]->GetVSize(); + blockOffsets[1] = m_feTheta->GetVSize(); // Get actual number of dofs *before* applying BCs + blockOffsets[2] = m_fePhi->GetVSize(); blockOffsets.PartialSum(); // Cumulative sum to get the offsets + return blockOffsets; +} - // Add integrators to block form one by one - // We add integrators corresponding to each term in the weak form - // The block form of the residual matrix - // ⎡ 0 -M ⎤ ⎡ θ ⎤ + ⎡f(θ)⎤ = ⎡ 0 ⎤ = R(X) - // ⎣ -Q D ⎦ ⎣ Φ ⎦ ⎣ 0 ⎦ ⎣ 0 ⎦ - // This then simplifies to - // ⎡f(θ) - MΘ ⎤ = ⎡ 0 ⎤ = R - // ⎣ -Qɸ Dθ ⎦ ⎣ 0 ⎦ - // Here M, Q, and D are - // M = ∫∇ψᶿ·Nᵠ dV (MixedVectorWeakDivergenceIntegrator) - // D = ∫ψᵠ·Nᵠ dV (VectorFEMassIntegrator) - // Q = ∫ψᵠ·∇Nᶿ dV (MixedVectorGradientIntegrator) - // f(θ) = ∫ψᶿ·θⁿ dV (NonlinearPowerIntegrator) - // Here ψᶿ and ψᵠ are the test functions for the theta and phi spaces, respectively - // Nᵠ and Nᶿ are the basis functions for the theta and phi spaces, respectively - // A full derivation of the weak form can be found in the 4DSSE documentation - +std::unique_ptr PolySolver::buildIndividualForms(const mfem::Array &blockOffsets) { // --- Assemble the MixedBilinear and Bilinear forms (M, D, and Q) --- - auto Mform = std::make_unique(m_fePhi.get(), m_feTheta.get()); - auto Qform = std::make_unique(m_feTheta.get(), m_fePhi.get()); - auto Dform = std::make_unique(m_fePhi.get()); - - // TODO: Check the sign on all of the integrators - Mform->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator()); - Qform->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator()); - Dform->AddDomainIntegrator(new mfem::VectorFEMassIntegrator()); + auto forms = std::make_unique( + std::make_unique(m_fePhi.get(), m_feTheta.get()), + std::make_unique(m_feTheta.get(), m_fePhi.get()), + std::make_unique(m_fePhi.get()), + std::make_unique(m_feTheta.get()) + ); - Mform->Assemble(); - Mform->Finalize(); + // --- -M negation -> M --- + mfem::Vector negOneVec(m_fePhi->GetVDim()); + negOneVec = -1.0; - Qform->Assemble(); - Qform->Finalize(); + m_negationCoeff = std::make_unique(negOneVec); - Dform->Assemble(); - Dform->Finalize(); + // --- Add the integrators to the forms --- + forms->M->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator(*m_negationCoeff)); + forms->Q->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator()); + forms->D->AddDomainIntegrator(new mfem::VectorFEMassIntegrator()); - // --- Assemble the NonlinearForm (f) --- - // Note that the nonlinear form is built here but its essential true dofs (boundary conditions) are - // not set until later (when solve is called). They are applied through the operator rather than directly. - auto fform = std::make_unique(m_feTheta.get()); - fform->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex)); + // --- Assemble and Finalize the forms --- + assembleAndFinalizeForm(forms->M); + assembleAndFinalizeForm(forms->Q); + assembleAndFinalizeForm(forms->D); - // -- Build the BlockOperator -- - m_polytropOperator = std::make_unique( - std::move(Mform), - std::move(Qform), - std::move(Dform), - std::move(fform), - blockOffsets - ); + forms->f->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex)); + return std::move(forms); +} + +void PolySolver::assembleAndFinalizeForm(auto &f) { + // This constructs / ensures the matrix representation for each form + // Assemble => Computes the local element matrices across the domain. Adds these to the global matrix . Adds these to the global matrix. + // Finalize => Builds the sparsity pattern and allows the SparseMatrix representation to be extracted. + f->Assemble(); + f->Finalize(); } void PolySolver::solve() const { @@ -193,6 +184,9 @@ void PolySolver::solve() const { throw std::runtime_error("PolytropeOperator is not finalized. Cannot solve."); } + GetDofCoordinates(*m_feTheta, "dof2posTheta.csv"); + GetDofCoordinates(*m_fePhi, "dof2posPhi.csv"); + // It's safer to get the offsets directly from the operator after finalization const mfem::Array& block_offsets = m_polytropOperator->GetBlockOffsets(); // Assuming a getter exists or accessing member if public/friend mfem::BlockVector state_vector(block_offsets); @@ -213,6 +207,8 @@ void PolySolver::solve() const { // with updating the preconditioner at every newton step as the // changes to the jacobian are automatically propagated through the // solving chain. This is at least true with MFEM 4.8-rc0 + std::string windowTitle = "testWindow"; + std::string keyset = ""; sb.newton.Mult(zero_rhs, state_vector); // --- Save and view the solution --- @@ -279,7 +275,8 @@ void PolySolver::setInitialGuess() const { const double radius = Probe::getMeshRadius(*m_mesh); const double u = 1/radius; - return -std::pow((u*r), 2)+1.0; // The series expansion is a better guess; however, this is cheaper and ensures that the value at the surface is very close to zero in a way that the series expansion does not + return (-1.0/radius) * r + 1; + // return -std::pow((u*r), 2)+1.0; // The series expansion is a better guess; however, this is cheaper and ensures that the value at the surface is very close to zero in a way that the series expansion does not } ); m_theta->ProjectCoefficient(thetaInitGuess); @@ -338,6 +335,43 @@ void PolySolver::LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); } +void PolySolver::GetDofCoordinates(mfem::FiniteElementSpace &fes, const std::string& filename) const { + mfem::Mesh *mesh = fes.GetMesh(); + double r = Probe::getMeshRadius(*mesh); + std::ofstream outputFile(filename, std::ios::out | std::ios::trunc); + outputFile << "dof,R,r,x,y,z" << '\n'; + + const int nElements = mesh->GetNE(); + + mfem::Vector coords; + mfem::IntegrationPoint ipZero; + double p[3] = {0.0, 0.0, 0.0}; + int actual_idx; + ipZero.Set3(p); + for (int i = 0; i < nElements; i++) { + mfem::Array elemDofs; + fes.GetElementDofs(i, elemDofs); + mfem::ElementTransformation* T = mesh->GetElementTransformation(i); + mfem::Vector physCoord(3); + T->Transform(ipZero, physCoord); + for (int dofID = 0; dofID < elemDofs.Size(); dofID++) { + if (elemDofs[dofID] < 0) { + actual_idx = -elemDofs[dofID] - 1; + } else { + actual_idx = elemDofs[dofID]; + } + outputFile << actual_idx; + if (dofID != elemDofs.Size() - 1) { + outputFile << "|"; + } else { + outputFile << ","; + } + } + outputFile << r << "," << physCoord.Norml2() << "," << physCoord[0] << "," << physCoord[1] << "," << physCoord[2] << '\n'; + } + outputFile.close(); +} + solverBundle PolySolver::setupNewtonSolver() const { // --- Load configuration parameters --- double newtonRelTol, newtonAbsTol, gmresRelTol, gmresAbsTol; diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index d496684..738bbe1 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -44,7 +44,14 @@ struct solverBundle { mfem::NewtonSolver newton; // Must be second so that when it is destroyed the solver is still alive preventing a double delete }; -class PolySolver { +struct formBundle { + std::unique_ptr M; + std::unique_ptr Q; + std::unique_ptr D; + std::unique_ptr f; +}; + +class PolySolver final{ public: // Public methods PolySolver(const double n, const double order); ~PolySolver(); @@ -76,9 +83,96 @@ private: // Private Attributes std::unique_ptr m_prec; + std::unique_ptr m_negationCoeff; + private: // Private methods void assembleBlockSystem(); + + /** + * @breif Compute the block offsets for the operator. These are the offsets that define which dofs belong to which variable. + * + * @details + * + * Create the block offsets. These define the start of each block in the combined vector. + * Block offsets will be [0, thetaDofs, thetaDofs + phiDofs]. + * The interpretation of this is that each block tells the operator where in the flattned (1D) vector + * the degrees of freedom or coefficients for that free parameter start and end. I.e. + * we know that in any flattened vector will have a size thetaDofs + phiDofs. The theta dofs will span + * from blockOffsets[0] -> blockOffsets[1] and the phiDofs will span from blockOffsets[1] -> blockOffsets[2]. + * + * This is the same for matrices only in 2D (rows and columns) + * + * The key point here is that this is fundamentally an accounting structure, it is here to keep track of what + * parts of vectors and matrices belong to which variable. + * + * Also note that we use VSize rather than Size. Size referees to the number of true dofs. That is the dofs which + * still are present in the system after eliminating boundary conditions. This is the wrong size to use if we are + * trying to account for the true size of the system. VSize on the other hand refers to the total number of dofs. + * + * @return blockOffsets The offsets for the blocks in the operator + */ + mfem::Array computeBlockOffsets() const; + + /** + * @breif Build the individual forms for the block operator (M, Q, D, and f) + * + * @param blockOffsets The offsets for the blocks in the operator + * @param Mform The mixed bilinear form for the mass matrix + * @param Qform The mixed bilinear form for the gradient matrix + * @param Dform The bilinear form for the divergence matrix + * @param fform The nonlinear form for the source term + * + * @note These forms are build exactly how they are defined in the derivation. This means that Mform -> M not -M and Qform -> Q not -Q. + * + * @details + * Computes the block offsets + * \f$\{0,\;|\theta|,\;|\theta|+|\phi|\}\f$, then builds and finalizes + * the MixedBilinearForms \c Mform, \c Qform and the BilinearForm \c Dform, + * plus the NonlinearForm \c fform. Finally, these are handed off to + * \c PolytropeOperator along with the offsets. + * + * The discretized weak form is + * \f[ + * R(X) + * = \begin{pmatrix} + * f(\theta) - M\,\theta \\[6pt] + * D\,\theta - Q\,\phi + * \end{pmatrix} + * = \mathbf{0}, + * \f] + * with + * \f[ + * M = \int \nabla\psi^\theta \;\cdot\; N^\phi \,dV,\quad + * D = \int \psi^\phi \;\cdot\; N^\phi \,dV, + * \quad + * Q = \int \psi^\phi \;\cdot\; \nabla N^\theta \,dV, + * \quad + * f(\theta) = \int \psi^\theta \;\cdot\; \theta^n \,dV. + * \f] + * + * @note MFEM’s MixedVectorWeakDivergenceIntegrator implements + * \f$ -\nabla\!\cdot \f$, so we supply a –1 coefficient to make + * `Mform` represent the +M from the derivation. The single negation + * in `PolytropeOperator` then restores the final block sign. + * + + * + * @pre \c m_feTheta and \c m_fePhi must be valid, populated FiniteElementSpaces. + * @post \c m_polytropOperator is constructed with assembled forms and offsets. + * + */ + std::unique_ptr buildIndividualForms(const mfem::Array& blockOffsets); + + /** + * @brief Assemble and finalize the form (Must be a form that can be assembled and finalized) + * + * @param f form which is to be assembled and finalized + * + * @pre f is a valid form that can be assembled and finalized (Such as Bilinear or MixedBilinearForm) + * @post f is assembled and finalized + */ + static void assembleAndFinalizeForm(auto &f); SSE::MFEMArrayPairSet getEssentialTrueDof() const; std::pair, mfem::Array> findCenterElement() const; void setInitialGuess() const; @@ -89,4 +183,6 @@ private: // Private methods void LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel, double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const; + void GetDofCoordinates(mfem::FiniteElementSpace &fes, const std::string& filename) const; + }; \ No newline at end of file diff --git a/src/poly/utils/private/integrators.cpp b/src/poly/utils/private/integrators.cpp index a70eac4..dab8c36 100644 --- a/src/poly/utils/private/integrators.cpp +++ b/src/poly/utils/private/integrators.cpp @@ -22,8 +22,11 @@ #include #include "integrators.h" +#include +// static std::ofstream debugOut("gradient.csv", std::ios::trunc); + namespace polyMFEMUtils { NonlinearPowerIntegrator::NonlinearPowerIntegrator(const double n) : m_polytropicIndex(n) {} @@ -40,11 +43,10 @@ namespace polyMFEMUtils { 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(); + const double weight = ip.weight * Trans.Weight(); el.CalcShape(ip, shape); @@ -52,10 +54,10 @@ namespace polyMFEMUtils { 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, m_polytropicIndex); - - double x2_u_nl = u_nl; + 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; @@ -71,18 +73,23 @@ namespace polyMFEMUtils { mfem::DenseMatrix &elmat) { const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); - int dof = el.GetDof(); + 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); - double weight = ip.weight * Trans.Weight(); + const double weight = ip.weight * Trans.Weight(); + el.CalcShape(ip, shape); - + double u_val = 0.0; for (int j = 0; j < dof; j++) { @@ -90,16 +97,37 @@ namespace polyMFEMUtils { } // Calculate the Jacobian - double u_safe = std::max(u_val, 0.0); - double d_u_nl = m_polytropicIndex * std::pow(u_safe, m_polytropicIndex - 1); - double x2_d_u_nl = d_u_nl; + 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 index ecbf6ef..5f16560 100644 --- a/src/poly/utils/private/operator.cpp +++ b/src/poly/utils/private/operator.cpp @@ -89,6 +89,35 @@ void writeDenseMatrixToCSV(const std::string &filename, int precision, const mfe 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, @@ -219,40 +248,6 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { } -// TODO: I was *very* stupid and I accidentally deleted a lot of the code -// which I had written to find and use the preconditioner. This needs -// to be reimplemented. Once that is working you can get back to -// Trying to understand the multimodal hump in the residuals vector. -// There is a jupyter notebook about this. I was thinking that this was -// perhaps related to the non consistent application of boundary conditions. -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 { - invMat = std::make_unique(diag); - } -} void PolytropeOperator::updateInverseNonlinearJacobian(const mfem::Operator &grad) const { if (const auto *sparse_mat = dynamic_cast(&grad); sparse_mat != nullptr) { @@ -325,4 +320,4 @@ void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_td 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/public/operator.h b/src/poly/utils/public/operator.h index bb469a6..d06b978 100644 --- a/src/poly/utils/public/operator.h +++ b/src/poly/utils/public/operator.h @@ -26,7 +26,7 @@ #include "probe.h" -class PolytropeOperator : public mfem::Operator { +class PolytropeOperator final : public mfem::Operator { public: PolytropeOperator( std::unique_ptr M, @@ -84,8 +84,8 @@ private: /* * The schur preconditioner has the form * - * ⎡ḟ(θ)^-1 0 ⎤ - * ⎣ 0 S^-1 ⎦ + * ⎡ḟ(θ)^-1 0 ⎤ + * ⎣ 0 S^-1 ⎦ * * Where S is the Schur compliment of the system *