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