From 8dcdf92414420a281a4447126ddf4f6c08982cff Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 18 Mar 2025 10:15:51 -0400 Subject: [PATCH] feat(poly): interpolating polynomial to find polytrope surface Instead of treating the polytrope as a free boundary problem I have defined an interpolating polynominal, accurate to within 0.01 percent over n=[0,5) which is used to set the size of the domain for a given n --- src/poly/coeff/private/polyCoeff.cpp | 12 +++++++++++- src/poly/coeff/public/polyCoeff.h | 13 +++++++++++++ src/poly/solver/private/polySolver.cpp | 15 +++------------ tests/poly/polyTest.cpp | 13 ++++++++----- 4 files changed, 35 insertions(+), 18 deletions(-) diff --git a/src/poly/coeff/private/polyCoeff.cpp b/src/poly/coeff/private/polyCoeff.cpp index aa9a826..44cc0d1 100644 --- a/src/poly/coeff/private/polyCoeff.cpp +++ b/src/poly/coeff/private/polyCoeff.cpp @@ -19,7 +19,7 @@ // // *********************************************************************** */ #include "mfem.hpp" -#include +#include #include "polyCoeff.h" @@ -35,4 +35,14 @@ namespace polycoeff{ v.SetSize(3); for (int i = 0; i < 3; i++) { v(i) = -1; } } + + double x1(const double n) + { + double r = 0; + for (int i = 0; i < x1InterpCoeff::numTerms; i++) { + r += x1InterpCoeff::coeff[i] * std::pow(n, x1InterpCoeff::power[i]); + } + return r; + } + } \ No newline at end of file diff --git a/src/poly/coeff/public/polyCoeff.h b/src/poly/coeff/public/polyCoeff.h index 0812498..4b18066 100644 --- a/src/poly/coeff/public/polyCoeff.h +++ b/src/poly/coeff/public/polyCoeff.h @@ -42,6 +42,19 @@ namespace polycoeff * @param v Output vector to store the computed xi coefficient. */ void diffusionCoeff(const mfem::Vector &x, mfem::Vector &v); + + double x1(const double n); + /** + * @brief Coefficients for the interpolations of the surface location of a polytrope + * @param numTerms Number of terms in the polynomial interpolator + * @param power Array of the powers of the polynomial interpolator + * @param coeff Array of the coefficients of the polynomial interpolator + */ + struct x1InterpCoeff { + constexpr static int numTerms = 51; + constexpr static int power[51] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50}; + constexpr static double coeff[51] = {2.449260049071003e+00, 4.340987403815444e-01, 2.592462347475860e+00, -2.283794512491552e+01, 1.135535208066129e+02, -3.460093673587010e+02, 6.948139716983651e+02, -9.564368645158952e+02, 9.184546798891624e+02, -6.130308987803503e+02, 2.747228735318193e+02, -7.416795903196430e+01, 7.318638538580859e+00, 1.749441306797260e+00, -4.240148582456829e-01, -5.818809544982156e-02, 1.514877217199105e-02, 3.228634707578998e-03, -2.862524323980516e-04, -1.622486968261819e-04, -1.253644717104076e-05, 4.334945141292894e-06, 1.296452565229763e-06, 8.634802209209870e-08, -3.337511676486084e-08, -1.094796628367775e-08, -1.228178760540410e-09, 1.416744125622751e-10, 8.513777351265677e-11, 1.624582561364811e-11, 8.377207519041114e-13, -4.363812865112836e-13, -1.535862757816461e-13, -2.485085045037669e-14, -6.566281276491033e-16, 8.405047965853478e-16, 2.673804441025638e-16, 4.176337890285142e-17, 8.150073570140493e-19, -1.531673805257016e-18, -4.746996933716653e-19, -6.976825828390195e-20, 8.807513368604331e-22, 3.307739310180278e-21, 8.593260940093030e-22, 7.385969061093440e-23, -2.211130577977291e-23, -8.557291048388455e-24, -4.169359901215994e-25, 4.609379358875657e-25, -3.590870335035984e-26}; + }; } // namespace polyCoeff #endif // POLYCOEFF_H \ No newline at end of file diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 258e71d..b568c3a 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -1,16 +1,7 @@ #include "mfem.hpp" #include -#include -#include -#include -#include -#include -#include -#include -#include -#include "meshIO.h" #include "polySolver.h" #include "polyMFEMUtils.h" #include "polyCoeff.h" @@ -19,7 +10,6 @@ #include "quill/LogMacros.h" -#include "warning_control.h" namespace laneEmden { @@ -122,7 +112,7 @@ void PolySolver::solve(){ nonlinearForm->SetEssentialTrueDofs(ess_tdof_list); // Set the center elemID to be the Dirichlet boundary - double alpha = config.get("Poly:Solver:Newton:Alpha", 1e2); + // double alpha = config.get("Poly:Solver:Newton:Alpha", 1e2); double newtonRelTol = config.get("Poly:Solver:Newton:RelTol", 1e-7); double newtonAbsTol = config.get("Poly:Solver:Newton:AbsTol", 1e-7); int newtonMaxIter = config.get("Poly:Solver:Newton:MaxIter", 200); @@ -137,7 +127,8 @@ void PolySolver::solve(){ LOG_INFO(logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); std::vector zeroSlopeCoordinate = {0.0, 0.0, 0.0}; - polyMFEMUtils::ZeroSlopeNewtonSolver newtonSolver(alpha, zeroSlopeCoordinate); + // polyMFEMUtils::ZeroSlopeNewtonSolver newtonSolver(alpha, zeroSlopeCoordinate); + mfem::NewtonSolver newtonSolver; newtonSolver.SetRelTol(newtonRelTol); newtonSolver.SetAbsTol(newtonAbsTol); newtonSolver.SetMaxIter(newtonMaxIter); diff --git a/tests/poly/polyTest.cpp b/tests/poly/polyTest.cpp index a4f5246..9a9639e 100644 --- a/tests/poly/polyTest.cpp +++ b/tests/poly/polyTest.cpp @@ -1,11 +1,12 @@ #include -#include +#include #include "quill/LogMacros.h" #include "mfem.hpp" #include "polySolver.h" +#include "polyCoeff.h" #include "probe.h" #include "config.h" #include "meshIO.h" @@ -23,14 +24,16 @@ TEST_F(polyTest, Solve) { LOG_INFO(logger, "Starting polytrope solve test 1..."); config.loadConfig(CONFIG_FILENAME); - double polyRadius = config.get("Tests:Poly:Radius", std::numbers::pi); + double polytropicIndex = config.get("Tests:Poly:Index", 1); + double polyRadius = polycoeff::x1(polytropicIndex); + std::cout << "Polytropic index: " << polytropicIndex << std::endl; + std::cout << "Polytropic radius: " << polyRadius << std::endl; + LOG_INFO(logger, "Solving polytrope with n = {:0.2f}", polytropicIndex); MeshIO meshIO(SPHERICAL_MESH, polyRadius); mfem::Mesh& mesh = meshIO.GetMesh(); double radius = Probe::getMeshRadius(mesh); - LOG_INFO(logger, "Mesh radius: {}", radius); + LOG_INFO(logger, "Mesh radius: {:0.4f} (target={:0.4f})", radius, polyRadius); - double polytropicIndex = config.get("Tests:Poly:Index", 1); - LOG_INFO(logger, "Solving polytrope with n = {:0.2f}", polytropicIndex); PolySolver polytrope(polytropicIndex, 1, mesh); LOG_INFO(logger, "Solving polytrope...");