fix(poly): working on solving polytrope

src/poly/solver/private/polySolver.cpp

@@ -104,26 +104,6 @@ void PolySolver::solve(){
     centerPoint(1, 0) = 0.0;
     centerPoint(2, 0) = 0.0;
 
-    // double controlPoint = 0.25;
-    // int sign;
-    // for (int i = 1; i < 7; i++) {
-    //     sign = i % 2 == 0 ? -1 : 1; 
-    //     if (i == 1 || i == 2) {
-    //         centerPoint(0, i) = controlPoint * sign;
-    //         centerPoint(1, i) = 0.0;
-    //         centerPoint(2, i) = 0.0;
-    //     } else if (i == 3 || i == 4) {
-    //         centerPoint(0, i) = 0.0;
-    //         centerPoint(1, i) = controlPoint * sign;
-    //         centerPoint(2, i) = 0.0;
-    //     } else {
-    //         centerPoint(0, i) = 0.0;
-    //         centerPoint(1, i) = 0.0;
-    //         centerPoint(2, i) = controlPoint * sign;
-    //     }
-    // }
-
-
     mfem::Array<int> elementIDs;
     mfem::Array<mfem::IntegrationPoint> ips;
     mesh.FindPoints(centerPoint, elementIDs, ips);
@@ -142,19 +122,32 @@ void PolySolver::solve(){
     nonlinearForm->SetEssentialTrueDofs(ess_tdof_list);
     // Set the center elemID to be the Dirichlet boundary
 
-    double alpha = config.get<double>("Poly:Solver: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);
+    int newtonPrintLevel = config.get<int>("Poly:Solver:Newton:PrintLevel", 1);
+
+    double gmresRelTol = config.get<double>("Poly:Solver:GMRES:RelTol", 1e-10);
+    double gmresAbsTol = config.get<double>("Poly:Solver:GMRES:AbsTol", 1e-12);
+    int gmresMaxIter = config.get<int>("Poly:Solver:GMRES:MaxIter", 2000);
+    int gmresPrintLevel = config.get<int>("Poly:Solver:GMRES:PrintLevel", 0);
+
+    LOG_INFO(logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel);
+    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);
-    newtonSolver.SetRelTol(1e-8);
+    newtonSolver.SetRelTol(newtonRelTol);
-    newtonSolver.SetAbsTol(1e-10);
+    newtonSolver.SetAbsTol(newtonAbsTol);
-    newtonSolver.SetMaxIter(200);
+    newtonSolver.SetMaxIter(newtonMaxIter);
-    newtonSolver.SetPrintLevel(1);
+    newtonSolver.SetPrintLevel(newtonPrintLevel);
     newtonSolver.SetOperator(*nonlinearForm);
     mfem::GMRESSolver gmresSolver;
-    gmresSolver.SetRelTol(1e-10);
+    gmresSolver.SetRelTol(gmresRelTol);
-    gmresSolver.SetAbsTol(1e-12);
+    gmresSolver.SetAbsTol(gmresAbsTol);
-    gmresSolver.SetMaxIter(2000);
+    gmresSolver.SetMaxIter(gmresMaxIter);
-    gmresSolver.SetPrintLevel(0);
+    gmresSolver.SetPrintLevel(gmresPrintLevel);
     newtonSolver.SetSolver(gmresSolver);
     // newtonSolver.SetAdaptiveLinRtol();
 

src/poly/utils/public/polyMFEMUtils.h

@@ -228,8 +228,11 @@ namespace polyMFEMUtils {
       std::vector<double> zeroSlopeCoordinate; // The coordinate of the zero slope point
 
       int zeroSlopeElemID = -1;
-      mfem::Array<int> zeroSlopeDofs;
+      double zeroIPReferenceCoord[4] = {0.0, 0.0, 0.0, 1.0};
-      mfem::IntegrationPoint zeroSlopeIP;
+      mfem::IntegrationPoint zeroIP;
+      mfem::Array<int> zeroSlopeConnectedElements;
+      std::vector<mfem::IntegrationPoint> zeroSlopeIPs;
+      std::vector<mfem::Array<int>> zeroSlopeDofs;
 
       std::unique_ptr<mfem::GridFunction> u_gf;
       mutable mfem::SparseMatrix *grad = nullptr;

tests/poly/polyTest.cpp

@@ -23,7 +23,8 @@ TEST_F(polyTest, Solve) {
     LOG_INFO(logger, "Starting polytrope solve test 1...");
     config.loadConfig(CONFIG_FILENAME);
 
-    MeshIO meshIO(SPHERICAL_MESH, std::numbers::pi);
+    double polyRadius = config.get<double>("Tests:Poly:Radius", std::numbers::pi);
+    MeshIO meshIO(SPHERICAL_MESH, polyRadius);
     mfem::Mesh& mesh = meshIO.GetMesh();
     double radius = Probe::getMeshRadius(mesh);
     LOG_INFO(logger, "Mesh radius: {}", radius);

tests/testsConfig.yaml

@@ -1,3 +1,5 @@
+Debug: true
+
 Probe:
   GLVis:
     Visualization: false
@@ -9,9 +11,17 @@ Probe:
 Poly:
   Solver:
     ViewInitialGuess: false
-    Alpha: 1e10
     GMRES:
       MaxIter: 5000
+      RelTol: 1.0e-8
+      AbsTol: 1.0e-10
+      PrintLevel: 0
+    Newton:
+      MaxIter: 200
+      RelTol: 1.0e-8
+      AbsTol: 1.0e-10
+      Alpha: 0.05
+      PrintLevel: 1
     Newton:
   Output:
     1D:
@@ -25,4 +35,5 @@ Poly:
 # ANY OPTIONS NEEDED FOR THE TEST SUITE SHOULD BE PLACED HERE
 Tests:
   Poly:
-    Index: 1.1
+    Index: 1.0
+    # Radius: 3.25