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

@@ -19,7 +19,7 @@
 //
 // *********************************************************************** */
 #include "mfem.hpp"
-#include <cmath>
+#include <iostream>
 
 #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;
+  }
+
 }

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

src/poly/solver/private/polySolver.cpp

@@ -1,16 +1,7 @@
 #include "mfem.hpp"
 
 #include <string>
-#include <iostream>
-#include <memory>
-#include <stdexcept>
-#include <csignal>
-#include <filesystem>
-#include <vector>
-#include <array>
-#include <utility>
 
-#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<double>("Poly:Solver:Newton:Alpha", 1e2);
+    // double alpha = config.get<double>("Poly:Solver:Newton:Alpha", 1e2);
     double newtonRelTol = config.get<double>("Poly:Solver:Newton:RelTol", 1e-7);
     double newtonAbsTol = config.get<double>("Poly:Solver:Newton:AbsTol", 1e-7);
     int newtonMaxIter = config.get<int>("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<double> 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);

tests/poly/polyTest.cpp

@@ -1,11 +1,12 @@
 #include <gtest/gtest.h>
 
-#include<numbers>
+#include <iostream>
 
 #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<double>("Tests:Poly:Radius", std::numbers::pi);
+    double polytropicIndex = config.get<double>("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<double>("Tests:Poly:Index", 1);
-    LOG_INFO(logger, "Solving polytrope with n = {:0.2f}", polytropicIndex);
 
     PolySolver polytrope(polytropicIndex, 1, mesh);
     LOG_INFO(logger, "Solving polytrope...");