From b2ddbc41e7c730cae426df4c85c7e4fe2c95e65a Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 26 Mar 2025 12:34:30 -0400 Subject: [PATCH] fix(poly): increased default gamma and changed initial guess The default gamma value has been upped to 1e4 which is enough to strongly constrain the solution to have zero slope at the core region. Further, the initial guess has been changed from a series expansion of theta to a simple quadratic that is one at origin and zero at the polytrope radius. This is faster to evaluate and seems to work just as well. --- src/poly/solver/private/polySolver.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 2ff1633..3a73522 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -104,7 +104,7 @@ void PolySolver::assembleNonlinearForm() { compositeIntegrator->add_integrator(nonlinearIntegrator.release()); // Add the contraint term \gamma(\nabla \theta(0)\cdot\nabla v(0))^{2} - double gamma = config.get("Poly:Solver:Constraint:Gamma", 1e2); + double gamma = config.get("Poly:Solver:Constraint:Gamma", 1e4); auto constraintIntegrator = std::make_unique( new polyMFEMUtils::ConstraintIntegrator(gamma, &mesh) ); @@ -118,12 +118,12 @@ void PolySolver::solve(){ mfem::FunctionCoefficient initCoeff ( [this](const mfem::Vector &x) { double r = x.Norml2(); - double theta = laneEmden::thetaSerieseExpansion(r, n, 10); - return theta; - // double radius = Probe::getMeshRadius(mesh); - // double u = 1/radius; + // double theta = laneEmden::thetaSerieseExpansion(r, n, 10); + // return theta; + double radius = Probe::getMeshRadius(mesh); + double u = 1/radius; - // return -std::pow((u*r), 2)+1.0; + return -std::pow((u*r), 2)+1.0; } ); u->ProjectCoefficient(initCoeff);