diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 679d8b3..258e71d 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/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 elementIDs; mfem::Array 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("Poly:Solver: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); + int newtonPrintLevel = config.get("Poly:Solver:Newton:PrintLevel", 1); + + double gmresRelTol = config.get("Poly:Solver:GMRES:RelTol", 1e-10); + double gmresAbsTol = config.get("Poly:Solver:GMRES:AbsTol", 1e-12); + int gmresMaxIter = config.get("Poly:Solver:GMRES:MaxIter", 2000); + int gmresPrintLevel = config.get("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 zeroSlopeCoordinate = {0.0, 0.0, 0.0}; polyMFEMUtils::ZeroSlopeNewtonSolver newtonSolver(alpha, zeroSlopeCoordinate); - newtonSolver.SetRelTol(1e-8); - newtonSolver.SetAbsTol(1e-10); - newtonSolver.SetMaxIter(200); - newtonSolver.SetPrintLevel(1); + newtonSolver.SetRelTol(newtonRelTol); + newtonSolver.SetAbsTol(newtonAbsTol); + newtonSolver.SetMaxIter(newtonMaxIter); + newtonSolver.SetPrintLevel(newtonPrintLevel); newtonSolver.SetOperator(*nonlinearForm); mfem::GMRESSolver gmresSolver; - gmresSolver.SetRelTol(1e-10); - gmresSolver.SetAbsTol(1e-12); - gmresSolver.SetMaxIter(2000); - gmresSolver.SetPrintLevel(0); + gmresSolver.SetRelTol(gmresRelTol); + gmresSolver.SetAbsTol(gmresAbsTol); + gmresSolver.SetMaxIter(gmresMaxIter); + gmresSolver.SetPrintLevel(gmresPrintLevel); newtonSolver.SetSolver(gmresSolver); // newtonSolver.SetAdaptiveLinRtol(); diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h index 8bb29bc..fa76e55 100644 --- a/src/poly/utils/public/polyMFEMUtils.h +++ b/src/poly/utils/public/polyMFEMUtils.h @@ -228,8 +228,11 @@ namespace polyMFEMUtils { std::vector zeroSlopeCoordinate; // The coordinate of the zero slope point int zeroSlopeElemID = -1; - mfem::Array zeroSlopeDofs; - mfem::IntegrationPoint zeroSlopeIP; + double zeroIPReferenceCoord[4] = {0.0, 0.0, 0.0, 1.0}; + mfem::IntegrationPoint zeroIP; + mfem::Array zeroSlopeConnectedElements; + std::vector zeroSlopeIPs; + std::vector> zeroSlopeDofs; std::unique_ptr u_gf; mutable mfem::SparseMatrix *grad = nullptr; diff --git a/tests/poly/polyTest.cpp b/tests/poly/polyTest.cpp index 92cf25f..a4f5246 100644 --- a/tests/poly/polyTest.cpp +++ b/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("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); diff --git a/tests/testsConfig.yaml b/tests/testsConfig.yaml index e9f6b19..7685ff5 100644 --- a/tests/testsConfig.yaml +++ b/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 \ No newline at end of file + Index: 1.0 + # Radius: 3.25 \ No newline at end of file