diff --git a/4DSSEConsole.sh b/4DSSEConsole.sh deleted file mode 100755 index a9e2b9f..0000000 --- a/4DSSEConsole.sh +++ /dev/null @@ -1,205 +0,0 @@ -#!/bin/bash - -# --- Functions --- - -# Display the main menu -display_menu() { - clear # Clear the screen for a clean menu - echo "--- Project Install Console ---" - echo "1. Reconfigure (with tests)" - echo "2. Reconfigure (without tests)" - echo "3. Build" - echo "4. Run Tests" - echo "5. Run Test with gdb (Specify test name)" - echo "6. Wipe (Clean build directory)" - echo "7. Generate Documentation" # Added documentation option - echo "8. Exit" - echo "-------------------------------" - echo -n "Enter your choice [1-8]: " -} - -# Reconfigure the project (meson setup) -reconfigure() { - local test_flag="$1" # Use local to avoid global variable pollution - - if [[ "$test_flag" == "true" ]]; then - meson setup build -Dbuild_tests=true --buildtype=debug - else - meson setup build -Dbuild_tests=false --buildtype=release - fi - - if [[ $? -eq 0 ]]; then # Check exit code for success - echo "Reconfiguration successful." - else - echo "Reconfiguration failed." - return 1 # Return a non-zero exit code to signal failure - fi -} - -# Build the project -build_project() { - meson compile -C build - - if [[ $? -eq 0 ]]; then - echo "Build successful." - else - echo "Build failed." - return 1 - fi -} - -# Run tests -run_tests() { - meson test -C build - - if [[ $? -eq 0 ]]; then - echo "Tests passed." # More specific message. - else - echo "Tests failed." - return 1 - fi -} - -# Run a specific test with GDB -run_test_gdb() { - local test_dir="build/tests" # Store the base test directory - local selected_test="" - local test_options=() - local test_index=1 - local choice - - # --- Find test executables --- - # Use find to get a list of potential test executables, *excluding* .p and .cpp.o files - # We sort it for consistent ordering. The '-type f' ensures we only find files. - # The '-executable' ensures the files are actually executable. - find "$test_dir" -type f -executable \ - ! -name "*.p" ! -name "*.cpp.o" ! -name "*.h.gch" | sort | while read -r file; do - # Remove the 'build/tests/' prefix and the file extension - test_name=$(basename "$file") - test_options+=("$test_name") - echo "$test_index. $test_name" - test_index=$((test_index + 1)) - done - - # Check if any tests were found - if [[ ${#test_options[@]} -eq 0 ]]; then - echo "No test executables found in build/tests." - return 1 - fi - - # --- Prompt for test selection --- - echo "-------------------------------" - echo -n "Enter the number of the test to debug (or 0 to cancel): " - read -r choice - - # --- Input Validation --- - if ! [[ "$choice" =~ ^[0-9]+$ ]]; then # Check if input is a number - echo "Invalid input. Please enter a number." - return 1 - fi - - if [[ "$choice" -eq 0 ]]; then - echo "Test selection cancelled." - return 0 - fi - - if [[ "$choice" -lt 1 || "$choice" -gt ${#test_options[@]} ]]; then - echo "Invalid test number. Please enter a number within the valid range." - return 1 - fi - - # --- Get the selected test name --- - selected_test="${test_options[$((choice - 1))]}" # Adjust for 0-based array indexing - - # --- Run GDB --- - # Construct full path. More reliable than changing directories. - local full_test_path="$test_dir/$selected_test" - - if [[ ! -x "$full_test_path" ]]; then - echo "Error: Test executable '$full_test_path' not found or not executable." - return 1 - fi - - gdb --args "$full_test_path" - -} - -# Wipe the build directory -wipe_build() { - echo -n "Are you sure you want to wipe the build directory? (y/N): " - read -r confirm - - if [[ "$confirm" == "y" || "$confirm" == "Y" ]]; then - rm -rf build - echo "Build directory wiped." - else - echo "Wipe cancelled." - fi -} - -generate_docs() { - echo "Generating documentation..." - doxygen - if [[ $? -ne 0 ]]; then # Check if doxygen succeeded - echo "Error: Doxygen failed." - return 1 - fi - - if [[ -d docs/latex ]]; then #Check if directory exist - cd docs/latex || return #cd into docs/latex and exit the function if unsuccesful - make - if [[ $? -ne 0 ]]; then #Check make - echo "Error: make in docs/latex failed." - cd ../.. || exit #go back, or exit the script if it fails - return 1 - fi - cd ../.. || exit # Go back to the project root - else - echo "Warning: docs/latex directory not found. Skipping LaTeX build." - fi - - echo "Documentation generated." -} - -# --- Main Loop --- - -while true; do - display_menu - read -r choice - - case "$choice" in - 1) - reconfigure "true" - ;; - 2) - reconfigure "false" - ;; - 3) - build_project - ;; - 4) - run_tests - ;; - 5) - run_test_gdb - ;; - 6) - wipe_build - ;; - 7) - generate_docs - ;; - 8) - echo "Exiting..." - break # Exit the while loop - ;; - *) - echo "Invalid choice. Please enter a number between 1 and 8." - ;; - esac - - echo -n "Press Enter to continue..." - read -r # Wait for Enter key press -done - -exit 0 #Good practice \ No newline at end of file diff --git a/docs/derivations/laneEmdenVariationalBlockForm.pdf b/docs/derivations/laneEmdenVariationalBlockForm.pdf index 71df7ad..a26eb45 100644 Binary files a/docs/derivations/laneEmdenVariationalBlockForm.pdf and b/docs/derivations/laneEmdenVariationalBlockForm.pdf differ diff --git a/docs/derivations/laneEmdenVariationalBlockForm.tex b/docs/derivations/laneEmdenVariationalBlockForm.tex index 589cbaf..77e0b36 100644 --- a/docs/derivations/laneEmdenVariationalBlockForm.tex +++ b/docs/derivations/laneEmdenVariationalBlockForm.tex @@ -14,7 +14,7 @@ \title{Deriving The Full Lane Emden Weak Form} \author{Emily M. Boudreaux} -\date{March 2025} +\date{April 2025} \begin{document} diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 4b24351..8f290ab 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -34,6 +34,8 @@ #include "resourceManagerTypes.h" #include "operator.h" +#include "debug.h" + #include "quill/LogMacros.h" @@ -149,14 +151,14 @@ void PolySolver::assembleBlockSystem() { // A full derivation of the weak form can be found in the 4DSSE documentation // --- Assemble the MixedBilinear and Bilinear forms (M, D, and Q) --- - auto Mform = std::make_unique(m_feTheta.get(), m_fePhi.get()); - auto Qform = std::make_unique(m_fePhi.get(), m_feTheta.get()); + 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(negOneCoeff)); - Qform->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator(negOneVCoeff)); - Dform->AddDomainIntegrator(new mfem::VectorFEMassIntegrator(oneCoeff)); + Mform->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator()); + Qform->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator()); + Dform->AddDomainIntegrator(new mfem::VectorFEMassIntegrator()); Mform->Assemble(); Mform->Finalize(); @@ -187,47 +189,24 @@ void PolySolver::solve(){ // --- Set the initial guess for the solution --- setInitialGuess(); - // --- Set the essential true dofs for the operator --- - mfem::Array theta_ess_tdof_list, phi_ess_tdof_list; - std::tie(theta_ess_tdof_list, phi_ess_tdof_list) = getEssentialTrueDof(); - m_polytropOperator->SetEssentialTrueDofs(theta_ess_tdof_list, phi_ess_tdof_list); + setupOperator(); - // --- Load configuration parameters --- - double newtonRelTol = m_config.get("Poly:Solver:Newton:RelTol", 1e-7); - double newtonAbsTol = m_config.get("Poly:Solver:Newton:AbsTol", 1e-7); - int newtonMaxIter = m_config.get("Poly:Solver:Newton:MaxIter", 200); - int newtonPrintLevel = m_config.get("Poly:Solver:Newton:PrintLevel", 1); + // 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); + state_vector.GetBlock(0) = *m_theta; + state_vector.GetBlock(1) = *m_phi; - double gmresRelTol = m_config.get("Poly:Solver:GMRES:RelTol", 1e-10); - double gmresAbsTol = m_config.get("Poly:Solver:GMRES:AbsTol", 1e-12); - int gmresMaxIter = m_config.get("Poly:Solver:GMRES:MaxIter", 2000); - int gmresPrintLevel = m_config.get("Poly:Solver:GMRES:PrintLevel", 0); + mfem::Vector zero_rhs(block_offsets.Last()); + zero_rhs = 0.0; - LOG_DEBUG(m_logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel); - LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); - // --- Set up the Newton solver --- - mfem::NewtonSolver newtonSolver; - newtonSolver.SetRelTol(newtonRelTol); - newtonSolver.SetAbsTol(newtonAbsTol); - newtonSolver.SetMaxIter(newtonMaxIter); - newtonSolver.SetPrintLevel(newtonPrintLevel); - newtonSolver.SetOperator(*m_polytropOperator); - mfem::GMRESSolver gmresSolver; - gmresSolver.SetRelTol(gmresRelTol); - gmresSolver.SetAbsTol(gmresAbsTol); - gmresSolver.SetMaxIter(gmresMaxIter); - gmresSolver.SetPrintLevel(gmresPrintLevel); - newtonSolver.SetSolver(gmresSolver); - // newtonSolver.SetAdaptiveLinRtol(); + mfem::NewtonSolver newtonSolver = setupNewtonSolver(); - mfem::Vector B(m_feTheta->GetTrueVSize()); - B = 0.0; - - newtonSolver.Mult(B, *m_theta); + newtonSolver.Mult(zero_rhs, state_vector); // --- Save and view the solution --- - saveAndViewSolution(); + saveAndViewSolution(state_vector); } @@ -295,22 +274,77 @@ void PolySolver::setInitialGuess() { } -void PolySolver::saveAndViewSolution() { +void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) { + mfem::BlockVector x_block(const_cast(state_vector), m_polytropOperator->GetBlockOffsets()); + mfem::Vector& x_theta = x_block.GetBlock(0); + mfem::Vector& x_phi = x_block.GetBlock(1); + bool doView = m_config.get("Poly:Output:View", false); if (doView) { - Probe::glVisView(*m_theta, *m_mesh, "solution"); + Probe::glVisView(x_theta, *m_feTheta, "θ Solution"); + Probe::glVisView(x_phi, *m_fePhi, "ɸ Solution"); } // --- Extract the Solution --- bool write11DSolution = m_config.get("Poly:Output:1D:Save", true); if (write11DSolution) { std::string solutionPath = m_config.get("Poly:Output:1D:Path", "polytropeSolution_1D.csv"); + std::string derivSolPath = "d" + solutionPath; + double rayCoLatitude = m_config.get("Poly:Output:1D:RayCoLatitude", 0.0); double rayLongitude = m_config.get("Poly:Output:1D:RayLongitude", 0.0); int raySamples = m_config.get("Poly:Output:1D:RaySamples", 100); std::vector rayDirection = {rayCoLatitude, rayLongitude}; - Probe::getRaySolution(*m_theta, *m_feTheta, rayDirection, raySamples, solutionPath); + Probe::getRaySolution(x_theta, *m_feTheta, rayDirection, raySamples, solutionPath); + // Probe::getRaySolution(x_phi, *m_fePhi, rayDirection, raySamples, derivSolPath); } +} + +void PolySolver::setupOperator() { + mfem::Array theta_ess_tdof_list, phi_ess_tdof_list; + std::tie(theta_ess_tdof_list, phi_ess_tdof_list) = getEssentialTrueDof(); + m_polytropOperator->SetEssentialTrueDofs(theta_ess_tdof_list, phi_ess_tdof_list); + + // -- Finalize the operator -- + m_polytropOperator->finalize(); + + if (!m_polytropOperator->isFinalized()) { + LOG_ERROR(m_logger, "PolytropeOperator is not finalized. Cannot solve."); + throw std::runtime_error("PolytropeOperator is not finalized. Cannot solve."); + } +} + +mfem::NewtonSolver PolySolver::setupNewtonSolver(){ + // --- Load configuration parameters --- + double newtonRelTol = m_config.get("Poly:Solver:Newton:RelTol", 1e-7); + double newtonAbsTol = m_config.get("Poly:Solver:Newton:AbsTol", 1e-7); + int newtonMaxIter = m_config.get("Poly:Solver:Newton:MaxIter", 200); + int newtonPrintLevel = m_config.get("Poly:Solver:Newton:PrintLevel", 1); + + double gmresRelTol = m_config.get("Poly:Solver:GMRES:RelTol", 1e-10); + double gmresAbsTol = m_config.get("Poly:Solver:GMRES:AbsTol", 1e-12); + int gmresMaxIter = m_config.get("Poly:Solver:GMRES:MaxIter", 2000); + int gmresPrintLevel = m_config.get("Poly:Solver:GMRES:PrintLevel", 0); + + LOG_DEBUG(m_logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel); + LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); + + // --- Set up the Newton solver --- + mfem::NewtonSolver newtonSolver; + newtonSolver.SetRelTol(newtonRelTol); + newtonSolver.SetAbsTol(newtonAbsTol); + newtonSolver.SetMaxIter(newtonMaxIter); + newtonSolver.SetPrintLevel(newtonPrintLevel); + newtonSolver.SetOperator(*m_polytropOperator); + mfem::GMRESSolver gmresSolver; + gmresSolver.SetRelTol(gmresRelTol); + gmresSolver.SetAbsTol(gmresAbsTol); + gmresSolver.SetMaxIter(gmresMaxIter); + gmresSolver.SetPrintLevel(gmresPrintLevel); + newtonSolver.SetSolver(gmresSolver); + // newtonSolver.SetAdaptiveLinRtol(); + + return newtonSolver; } \ No newline at end of file diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 92ceec9..eb89cf1 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -1,6 +1,7 @@ #ifndef POLYSOLVER_H #define POLYSOLVER_H +#include "linalg/solvers.hpp" #include "mfem.hpp" #include #include @@ -53,7 +54,9 @@ private: // Private methods std::pair, mfem::Array> getEssentialTrueDof(); mfem::Array findCenterElement(); void setInitialGuess(); - void saveAndViewSolution(); + void saveAndViewSolution(const mfem::BlockVector& state_vector); + mfem::NewtonSolver setupNewtonSolver(); + void setupOperator(); }; diff --git a/src/poly/utils/private/integrators.cpp b/src/poly/utils/private/integrators.cpp index 8c7e1a5..7590604 100644 --- a/src/poly/utils/private/integrators.cpp +++ b/src/poly/utils/private/integrators.cpp @@ -60,7 +60,8 @@ namespace polyMFEMUtils { double u_safe = std::max(u_val, 0.0); double u_nl = std::pow(u_safe, m_polytropicIndex); - double coeff_val = m_coeff.Eval(Trans, ip); + // double coeff_val = m_coeff.Eval(Trans, ip); + double coeff_val = 1.0; double x2_u_nl = coeff_val * u_nl; for (int i = 0; i < dof; i++){ @@ -94,7 +95,8 @@ namespace polyMFEMUtils { for (int j = 0; j < dof; j++) { u_val += elfun(j) * shape(j); } - double coeff_val = m_coeff.Eval(Trans, ip); + // double coeff_val = m_coeff.Eval(Trans, ip); + double coeff_val = 1.0; // Calculate the Jacobian diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp index 4de7639..7dd3815 100644 --- a/src/poly/utils/private/operator.cpp +++ b/src/poly/utils/private/operator.cpp @@ -3,6 +3,8 @@ #include "linalg/vector.hpp" #include +#include "debug.h" + PolytropeOperator::PolytropeOperator( std::unique_ptr M, @@ -19,8 +21,12 @@ PolytropeOperator::PolytropeOperator( m_Q = std::move(Q); m_D = std::move(D); m_f = std::move(f); +} - +void PolytropeOperator::finalize() { + if (m_isFinalized) { + return; + } m_Mmat = std::make_unique(m_M->SpMat()); m_Qmat = std::make_unique(m_Q->SpMat()); m_Dmat = std::make_unique(m_D->SpMat()); @@ -28,14 +34,14 @@ PolytropeOperator::PolytropeOperator( m_negM_op = std::make_unique(m_Mmat.get(), -1.0); m_negQ_op = std::make_unique(m_Qmat.get(), -1.0); - MFEM_ASSERT(m_Mmat.get() != nullptr, "Matrix m_Mmat is null in PolytropeOperator constructor"); - MFEM_ASSERT(m_Qmat.get() != nullptr, "Matrix m_Qmat is null in PolytropeOperator constructor"); - MFEM_ASSERT(m_Dmat.get() != nullptr, "Matrix m_Dmat is null in PolytropeOperator constructor"); - MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator constructor"); + m_isFinalized = true; } 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); @@ -89,6 +95,9 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { } 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); @@ -111,6 +120,7 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const { void PolytropeOperator::SetEssentialTrueDofs(const mfem::Array &theta_ess_tofs, const mfem::Array &phi_ess_tofs) { + m_isFinalized = false; m_theta_ess_tofs = theta_ess_tofs; m_phi_ess_tofs = phi_ess_tofs; diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h index 00d02e9..b809280 100644 --- a/src/poly/utils/public/operator.h +++ b/src/poly/utils/public/operator.h @@ -20,13 +20,19 @@ public: void SetEssentialTrueDofs(const mfem::Array &theta_ess_tofs, const mfem::Array &phi_ess_tofs); + bool isFinalized() const { return m_isFinalized; } + + void finalize(); + + const mfem::Array& GetBlockOffsets() const { return m_blockOffsets; } + private: 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; + const mfem::Array m_blockOffsets; mfem::Array m_theta_ess_tofs; mfem::Array m_phi_ess_tofs; @@ -39,6 +45,8 @@ private: std::unique_ptr m_negM_op; std::unique_ptr m_negQ_op; mutable std::unique_ptr m_jacobian; + + bool m_isFinalized = false; }; diff --git a/src/probe/public/probe.h b/src/probe/public/probe.h index d9157f6..dfe3427 100644 --- a/src/probe/public/probe.h +++ b/src/probe/public/probe.h @@ -25,7 +25,6 @@ #include #include #include -#include #include #include "mfem.hpp" @@ -65,7 +64,6 @@ namespace Probe { */ void glVisView(mfem::Vector &vec, mfem::FiniteElementSpace &fes, const std::string &windowTitle = "vector", const std::string& keyset=""); - double getMeshRadius(mfem::Mesh& mesh); diff --git a/tests/testsConfig.yaml b/tests/testsConfig.yaml index 0032269..83cbc2a 100644 --- a/tests/testsConfig.yaml +++ b/tests/testsConfig.yaml @@ -3,8 +3,8 @@ Debug: true Probe: GLVis: Visualization: true - Host: "10.28.92.45" - # Host: "localhost" + # Host: "10.28.92.45" + Host: "localhost" Port: 19916 DefaultKeyset: "iimmMcaa" GetRaySolution: @@ -25,6 +25,7 @@ Poly: Constraint: Gamma: 1e2 Output: + View: true 1D: Save: true Path: "output/Poly/1D/poly.csv"