diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index f0df901..679d8b3 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -54,29 +54,21 @@ namespace laneEmden { // TODO: Come back to this and think of a better way to get the mesh file const std::string SPHERICAL_MESH = std::string(getenv("MESON_SOURCE_ROOT")) + "/src/resources/mesh/core.msh"; -PolySolver::PolySolver(double n, double order) +PolySolver::PolySolver(double n, double order, mfem::Mesh& mesh_) : logger(logManager.getLogger("log")), n(n), order(order), - meshIO(SPHERICAL_MESH, 3.1415), // TODO : Change this from PI (set to PI right now for testing the n = 1 case) - mesh(meshIO.GetMesh()), + mesh(mesh_), feCollection(std::make_unique(order, mesh.SpaceDimension())), feSpace(std::make_unique(&mesh, feCollection.get())), compositeIntegrator(std::make_unique()), nonlinearForm(std::make_unique(feSpace.get())), - C(std::make_unique(feSpace.get())), - u(std::make_unique(feSpace.get())), - diffusionCoeff(std::make_unique([&](){ - mfem::Vector diffusionCoeffVec(mesh.SpaceDimension()); - diffusionCoeffVec = 1.0; - return diffusionCoeffVec; - }())), - nonLinearSourceCoeff(std::make_unique(-1.0)), - gaussianCoeff(std::make_unique(config.get("Poly:Gaussian:Sigma", 0.1))) { - // C_val is the weighted average of the constraint function - C_val = polyMFEMUtils::calculateGaussianIntegral(mesh, *gaussianCoeff); + u(std::make_unique(feSpace.get())) { + + diffusionCoeff = std::make_unique(mesh.SpaceDimension(), polycoeff::diffusionCoeff); + nonlinearSourceCoeff = std::make_unique(polycoeff::nonlinearSourceCoeff); + assembleNonlinearForm(); - assembleConstraintForm(); } @@ -90,18 +82,12 @@ void PolySolver::assembleNonlinearForm() { compositeIntegrator->add_integrator(wrappedDiffusionIntegrator.release()); // Add the \int_{\Omega}v\theta^{n} d\Omega term - auto nonLinearIntegrator = std::make_unique(*nonLinearSourceCoeff, n); - compositeIntegrator->add_integrator(nonLinearIntegrator.release()); + auto nonlinearIntegrator = std::make_unique(*nonlinearSourceCoeff, n); + compositeIntegrator->add_integrator(nonlinearIntegrator.release()); nonlinearForm->AddDomainIntegrator(compositeIntegrator.release()); } -void PolySolver::assembleConstraintForm() { - auto constraintIntegrator = std::make_unique(*gaussianCoeff); - C->AddDomainIntegrator(constraintIntegrator.release()); - C->Assemble(); -} - void PolySolver::solve(){ // --- Set the initial guess for the solution --- mfem::FunctionCoefficient initCoeff ( @@ -112,181 +98,84 @@ void PolySolver::solve(){ } ); u->ProjectCoefficient(initCoeff); - std::string initGuessFilename = "output/Poly/Debug/Newton/1D/initial_guess.csv"; - Probe::getRaySolution(*u, *feSpace->GetMesh(), {0.0, 0.0}, 100, initGuessFilename); - if (config.get("Poly:Solver:ViewInitialGuess", false)) { - Probe::glVisView(*u, mesh, "initial_guess"); + // mfem::DenseMatrix centerPoint(mesh.SpaceDimension(), 7); + mfem::DenseMatrix centerPoint(mesh.SpaceDimension(), 1); + centerPoint(0, 0) = 0.0; + 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); + mfem::Array centerDofs; + mfem::Array tempDofs; + for (int i = 0; i < elementIDs.Size(); i++) { + feSpace->GetElementDofs(elementIDs[i], tempDofs); + centerDofs.Append(tempDofs); } + mfem::Array ess_tdof_list; + mfem::Array ess_brd(mesh.bdr_attributes.Max()); + ess_brd = 1; + feSpace->GetEssentialTrueDofs(ess_brd, ess_tdof_list); + // combine the essential dofs with the center dofs + ess_tdof_list.Append(centerDofs); + nonlinearForm->SetEssentialTrueDofs(ess_tdof_list); + // Set the center elemID to be the Dirichlet boundary - // --- Combine DOFs (u and λ) into a single vector --- - int lambdaDofOffset = feSpace->GetTrueVSize(); // Get the size of θ space - int totalTrueDofs = lambdaDofOffset + 1; - - mfem::Vector U(totalTrueDofs); - U = 0.0; - - mfem::Vector u_view(U.GetData(), lambdaDofOffset); - u->GetTrueDofs(u_view); - - // --- Setup the Augmented Operator --- - polyMFEMUtils::AugmentedOperator aug_op(*nonlinearForm, *C, lambdaDofOffset, C_val); - - // --- Create the RHS of the augmented system --- - mfem::Vector B(totalTrueDofs); - B = 0.0; - B[lambdaDofOffset] = C_val; - - - // --- Custom Newton Solver --- + double alpha = config.get("Poly:Solver:Alpha", 1e2); + 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.SetOperator(*nonlinearForm); mfem::GMRESSolver gmresSolver; - gmresSolver.SetRelTol(config.get("Poly:Solver:GMRES:RelTol", 1e-8)); - gmresSolver.SetAbsTol(config.get("Poly:Solver:GMRES:AbsTol", 1e-10)); - gmresSolver.SetMaxIter(config.get("Poly:Solver:GMRES:MaxIter", 2000)); - gmresSolver.SetPrintLevel(config.get("Poly:Solver:GMRES:PrintLevel", 0)); + gmresSolver.SetRelTol(1e-10); + gmresSolver.SetAbsTol(1e-12); + gmresSolver.SetMaxIter(2000); + gmresSolver.SetPrintLevel(0); + newtonSolver.SetSolver(gmresSolver); + // newtonSolver.SetAdaptiveLinRtol(); - std::cout << "Setting the Block ILU preconditioner size too " << feSpace->GetTypicalFE()->GetDof() << std::endl; - mfem::BlockILU prec(feSpace->GetTypicalFE()->GetDof(), mfem::BlockILU::Reordering::MINIMUM_DISCARDED_FILL); - gmresSolver.SetPreconditioner(prec); + mfem::Vector B(feSpace->GetTrueVSize()); + B = 0.0; - int iteration = 0; - const int maxIter = config.get("Poly:Solver:Newton:MaxIterations", 200); - const double relTol = config.get("Poly:Solver:Newton:RelTol", 1e-8); - const double absTol = config.get("Poly:Solver:Newton:AbsTol", 1e-10); + newtonSolver.Mult(B, *u); - bool writeIntermediate = config.get("Poly:Debug:Newton:1D:WriteIntermediate", false); - double rayLatitude = config.get("Poly:Debug:Newton:1D:lat", 0.0); - double rayLongitude = config.get("Poly:Debug:Newton:1D:lon", 0.0); - int raySamples = config.get("Poly:Debug:Newton:1D:radialPoints", 100); - double rayMin = config.get("Poly:Debug:Newton:1D:radialMin", 0.0); - double rayMax = config.get("Poly:Debug:Newton:1D:radialMax", 3.14); - double rayStep = (rayMax - rayMin) / raySamples; - int stepsPerWrite = config.get("Poly:Debug:Newton:1D:StepsPerWrite", 1); - bool exitAfterWrite = config.get("Poly:Debug:Newton:1D:Exit", false); - std::string outputDirectory = config.get("Poly:Debug:Newton:1D:OutputDir", "output/Poly/Debug/Newton/1D"); - std::pair, std::vector> samples; - std::vector radialPoints; - radialPoints.reserve(raySamples); - for (int i = 0; i < raySamples; i++) { - radialPoints.push_back(rayMin + i * rayStep); - } - std::vector rayDirection = {rayLatitude, rayLongitude}; + Probe::glVisView(*u, mesh, "solution"); - if (writeIntermediate) { - std::filesystem::create_directories(outputDirectory); + // --- Extract the Solution --- + bool write11DSolution = config.get("Poly:Output:1D:Save", true); + if (write11DSolution) { + std::string solutionPath = config.get("Poly:Output:1D:Path", "polytropeSolution_1D.csv"); + double rayCoLatitude = config.get("Poly:Output:1D:RayCoLatitude", 0.0); + double rayLongitude = config.get("Poly:Output:1D:RayLongitude", 0.0); + int raySamples = config.get("Poly:Output:1D:RaySamples", 100); + + std::vector rayDirection = {rayCoLatitude, rayLongitude}; + + Probe::getRaySolution(*u, *feSpace, rayDirection, raySamples, solutionPath); } - std::string keyset = config.get("Poly:Debug:Newton:GLVis:Keyset", ""); - bool view = config.get("Poly:Debug:Newton:GLVis:View", false); - bool doExit = config.get("Poly:Debug:Newton:GLVis:Exit", false); - int stepsPerView = config.get("Poly:Debug:Newton:GLVis:StepsPerView", 1); - - while (iteration < maxIter) { - mfem::Vector F(totalTrueDofs); - F = 0.0; - aug_op.Mult(U, F); // F now holds augOp(U) - F -= B; - double resNorm = F.Norml2(); - std::cout << "Iteration: " << iteration << " Residual Norm: [ " << resNorm << " ] --- "; - if (resNorm < relTol || resNorm < absTol) { - std::cout << "Convergence achieved!" << std::endl; - break; - } - - // --- Retrieve the Jacobian --- - mfem::Operator &gradOp = aug_op.GetGradient(U); - std::cout << "Size of the Jacobian: " << gradOp.Height() << " x " << gradOp.Width() << std::endl; - gmresSolver.SetOperator(gradOp); - mfem::SparseMatrix *J = dynamic_cast(&gradOp); - if (!J) { - MFEM_ABORT("GetGradient did not return a SparseMatrix"); - } - std::cout << "Jacobian: " << J->Height() << " x " << J->Width() << std::endl; - std::cout << "Non-zero entries: " << J->NumNonZeroElems() << std::endl; - - // --- Solve the Newton Step: J * step = -F --- - mfem::Vector minusF(totalTrueDofs); - minusF = F; // MFEM's vector class does not overload the unary minus operator - minusF *= -1.0; - mfem::Vector step(totalTrueDofs); - step = 0.0; - gmresSolver.Mult(minusF, step); - double stepNorm = step.Norml2(); - std::cout << "Step Norm: " << stepNorm << std::endl; - - U += step; - - // Silly, but a way to manually force the central value to = 1 - mfem::Array elementIds; - mfem::Array ips; - mfem::DenseMatrix rayPoints(3, 1); - rayPoints(0, 0) = 0.0; - rayPoints(1, 0) = 0.0; - rayPoints(2, 0) = 0.0; - mesh.FindPoints(rayPoints, elementIds, ips); - mfem::Array dofs; - feSpace->GetElementDofs(elementIds[0], dofs); - for (int dofID : dofs) { - U[dofID] = 1.0; - } - - - - if (view && iteration % stepsPerView == 0) { - std::string s_iteration = std::to_string(iteration); - Probe::glVisView(U, *feSpace, "U at " + s_iteration, keyset); - if (doExit) { - std::raise(SIGINT); - } - } - - if (writeIntermediate && iteration % stepsPerWrite == 0) { - std::string s_iteration = std::to_string(iteration); - std::string filename = outputDirectory + "/U_" + s_iteration + ".csv"; - Probe::getRaySolution(U, *feSpace, rayDirection, raySamples, filename); - if (exitAfterWrite) { - std::raise(SIGINT); - } - } - - bool endOfStepPause = config.get("Poly:Debug:Newton:EndOfStepPause", false); - if (endOfStepPause) { - Probe::pause(); - } - iteration++; - } - - - // // --- Setup the Newton Solver --- - // mfem::NewtonSolver newtonSolver; - // newtonSolver.SetRelTol(1e-8); - // newtonSolver.SetAbsTol(1e-10); - // newtonSolver.SetMaxIter(200); - // newtonSolver.SetPrintLevel(1); - // newtonSolver.SetOperator(aug_op); - - // // --- Setup the GMRES Solver --- - // // --- GMRES is good for indefinite systems --- - // mfem::GMRESSolver gmresSolver; - // gmresSolver.SetRelTol(1e-10); - // gmresSolver.SetAbsTol(1e-12); - // gmresSolver.SetMaxIter(2000); - // gmresSolver.SetPrintLevel(0); - // newtonSolver.SetSolver(gmresSolver); - // // TODO: Change numeric tolerance to grab from the tol module - - - // // --- Solve the augmented system --- - // newtonSolver.Mult(B, U); - - // // --- Extract the Solution --- - // mfem::Vector u_sol_view(U.GetData(), lambdaDofOffset); - - // DEPRECATION_WARNING_OFF // DISABLE DEPRECATION WARNING - // u->SetData(u_sol_view); - // DEPRECATION_WARNING_ON // REENABLE DEPRECATION WARNING - - // double lambda = U[lambdaDofOffset]; - - // std::cout << "λ = " << lambda << std::endl; - // // TODO : Add a way to get the solution out of the solver } \ No newline at end of file diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 443f7c6..0e9a75c 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -25,32 +25,24 @@ private: Probe::LogManager& logManager = Probe::LogManager::getInstance(); quill::Logger* logger; double n, order; - MeshIO meshIO; mfem::Mesh& mesh; std::unique_ptr feCollection; - std::unique_ptr feSpace; std::unique_ptr compositeIntegrator; - std::unique_ptr nonlinearForm; - std::unique_ptr C; // For the constraint equation std::unique_ptr u; - std::unique_ptr diffusionCoeff; - std::unique_ptr nonLinearSourceCoeff; - std::unique_ptr gaussianCoeff; - - double C_val; + std::unique_ptr diffusionCoeff; + std::unique_ptr nonlinearSourceCoeff; void assembleNonlinearForm(); - void assembleConstraintForm(); public: - PolySolver(double n, double order); + PolySolver(double n, double order, mfem::Mesh& mesh_); ~PolySolver(); void solve(); mfem::Mesh& getMesh() { return mesh; } diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp index 0900ae4..40231f7 100644 --- a/src/poly/utils/private/polyMFEMUtils.cpp +++ b/src/poly/utils/private/polyMFEMUtils.cpp @@ -5,6 +5,8 @@ #include #include #include +#include +#include #include "polyMFEMUtils.h" #include "probe.h" @@ -371,4 +373,323 @@ namespace polyMFEMUtils { return gaussianIntegral(one_gf); } + + ZeroSlopeNewtonSolver::ZeroSlopeNewtonSolver(double alpha_, std::vector zeroSlopeCoordinate_) + : alpha(alpha_), zeroSlopeCoordinate(zeroSlopeCoordinate_) {} + + ZeroSlopeNewtonSolver::~ZeroSlopeNewtonSolver() {} + + void ZeroSlopeNewtonSolver::SetOperator(const mfem::Operator &op) { + LOG_INFO(logger, "Setting operator for zero slope constraint..."); + mfem::NewtonSolver::SetOperator(op); // Call the base class method + LOG_INFO(logger, "Setting operator for zero slope constraint...done"); + LOG_INFO(logger, "Building location of zero slope constraint..."); + + mfem::NonlinearForm *nlf = dynamic_cast(const_cast(&op)); + if (!nlf) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::SetOperator: input operator is not a NonlinearForm"); + MFEM_ABORT("ZeroSlopeNewtonSolver::SetOperator: input operator is not a NonlinearForm"); + } + + mfem::FiniteElementSpace *fes = nlf->FESpace(); + + if (!fes) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::SetOperator: input operator does not have a finite element space"); + MFEM_ABORT("ZeroSlopeNewtonSolver::SetOperator: input operator does not have a finite element space"); + } + u_gf = std::make_unique(fes); + + mfem::Mesh *mesh = fes->GetMesh(); + + if (!mesh) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::SetOperator: input operator does not have a mesh"); + MFEM_ABORT("ZeroSlopeNewtonSolver::SetOperator: input operator does not have a mesh"); + } + + if (mesh->SpaceDimension() != static_cast(zeroSlopeCoordinate.size())) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::SetOperator: input operator mesh dimension does not match the zero slope coordinate dimension"); + MFEM_ABORT("ZeroSlopeNewtonSolver::SetOperator: input operator mesh dimension does not match the zero slope coordinate dimension"); + } + + mfem::DenseMatrix zeroSlopeCoordinateMatrix(mesh->SpaceDimension(), 1); + for (int dimID = 0; dimID < mesh->SpaceDimension(); dimID++) { + zeroSlopeCoordinateMatrix(dimID, 0) = zeroSlopeCoordinate[dimID]; + } + + mfem::Array elementsIDs; + mfem::Array ips; + mesh->FindPoints(zeroSlopeCoordinateMatrix, elementsIDs, ips); + + zeroSlopeElemID = elementsIDs[0]; + zeroSlopeIP = ips[0]; + + LOG_INFO(logger, "Getting element dofs for zero slope constraint..."); + fes->GetElementDofs(zeroSlopeElemID, zeroSlopeDofs); + LOG_INFO(logger, "Getting element dofs for zero slope constraint...done"); + LOG_INFO(logger, "Building location of zero slope constraint...done"); + } + + // void ZeroSlopeNewtonSolver::ProcessNewState(const mfem::Vector &x) const { + // LOG_INFO(logger, "Processing new state for zero slope constraint..."); + // if (zeroSlopeElemID < 0) { + // LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: zero slope element ID is not set"); + // MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: zero slope element ID is not set"); + // } + + // mfem::NonlinearForm *nlf = dynamic_cast(const_cast(oper)); + // if (!nlf) { + // LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator is not a NonlinearForm"); + // MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator is not a NonlinearForm"); + // } + + // mfem::FiniteElementSpace *fes = nlf->FESpace(); + // if (!fes) { + // LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a finite element space"); + // MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a finite element space"); + // } + + // mfem::Mesh *mesh = fes->GetMesh(); + // if (!mesh) { + // LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh"); + // MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh"); + // } + + // mfem::ElementTransformation *T = mesh->GetElementTransformation(zeroSlopeElemID); + // if (!T) { + // LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found"); + // MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found"); + // } + + // mfem::Vector grad_u(3); + + // mfem::GridFunction u_gf(fes); + // DEPRECATION_WARNING_OFF + // u_gf.SetData(x.GetData()); + // DEPRECATION_WARNING_ON + + // T->SetIntPoint(&zeroSlopeIP); + // u_gf.GetGradient(*T, grad_u); + + // int dof; + // LOG_DEBUG(logger, "Adjusting the residual to enforce the zero slope constraint by {:0.4E}...", -alpha*grad_u[0]); + // double rNorm = r.Norml2(); + // LOG_INFO(logger, "||r_B|| = {:0.4E}", rNorm); + // for (int i = 0; i < zeroSlopeDofs.Size(); i++) { + // dof = zeroSlopeDofs[i]; + // r[dof] -= alpha * grad_u[0]; + // r[dof] -= alpha * grad_u[1]; + // r[dof] -= alpha * grad_u[2]; + // } + // rNorm = r.Norml2(); + // LOG_INFO(logger, "||r_A|| = {:0.4E}", rNorm); + // // This still is not working; however, I think I am close. I also need to modify the jacobain. + // } + + void ZeroSlopeNewtonSolver::Mult(const mfem::Vector &b, mfem::Vector &x) const { + using namespace mfem; + using namespace std; + + MFEM_VERIFY(oper != NULL, "the Operator is not set (use SetOperator)."); + MFEM_VERIFY(prec != NULL, "the Solver is not set (use SetSolver)."); + + int it; + real_t norm0, norm, norm_goal; + const bool have_b = (b.Size() == Height()); + + if (!iterative_mode) + { + x = 0.0; + } + + ProcessNewState(x); + + oper->Mult(x, r); + if (have_b) + { + r -= b; + } + // ComputeConstrainedResidual(x, r); + + norm0 = norm = initial_norm = Norm(r); + if (print_options.first_and_last && !print_options.iterations) + { + mfem::out << "Zero slope newton iteration " << setw(2) << 0 + << " : ||r|| = " << norm << "...\n"; + } + norm_goal = std::max(rel_tol*norm, abs_tol); + + prec->iterative_mode = false; + + // x_{i+1} = x_i - [DF(x_i)]^{-1} [F(x_i)-b] + for (it = 0; true; it++) + { + MFEM_VERIFY(IsFinite(norm), "norm = " << norm); + if (print_options.iterations) + { + mfem::out << "Zero slope newton iteration " << setw(2) << it + << " : ||r|| = " << norm; + if (it > 0) + { + mfem::out << ", ||r||/||r_0|| = " << norm/norm0; + } + mfem::out << '\n'; + } + Monitor(it, norm, r, x); + + if (norm <= norm_goal) + { + converged = true; + break; + } + + if (it >= max_iter) + { + converged = false; + break; + } + + grad = dynamic_cast(&oper->GetGradient(x)); + if (!grad) + { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::Mult: Operator does not return a SparseMatrix"); + MFEM_ABORT("ZeroSlopeNewtonSolver::Mult: Operator does not return a SparseMatrix"); + } + ComputeConstrainedGradient(x); + prec->SetOperator(*grad); + + if (lin_rtol_type) + { + AdaptiveLinRtolPreSolve(x, it, norm); + } + + prec->Mult(r, c); // c = [DF(x_i)]^{-1} [F(x_i)-b] + + if (lin_rtol_type) + { + AdaptiveLinRtolPostSolve(c, r, it, norm); + } + + const real_t c_scale = ComputeScalingFactor(x, b); + if (c_scale == 0.0) + { + converged = false; + break; + } + add(x, -c_scale, c, x); + + ProcessNewState(x); + + oper->Mult(x, r); + if (have_b) + { + r -= b; + } + // ComputeConstrainedResidual(x, r); + norm = Norm(r); + } + + final_iter = it; + final_norm = norm; + + if (print_options.summary || (!converged && print_options.warnings) || + print_options.first_and_last) + { + mfem::out << "Newton: Number of iterations: " << final_iter << '\n' + << " ||r|| = " << final_norm + << ", ||r||/||r_0|| = " << final_norm/norm0 << '\n'; + } + if (!converged && (print_options.summary || print_options.warnings)) + { + mfem::out << "Newton: No convergence!\n"; + } + } + + void ZeroSlopeNewtonSolver::ComputeConstrainedResidual(const mfem::Vector &x, mfem::Vector &residual) const { + mfem::NonlinearForm *nlf = dynamic_cast(const_cast(oper)); + if (!nlf) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator is not a NonlinearForm"); + MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator is not a NonlinearForm"); + } + + mfem::FiniteElementSpace *fes = nlf->FESpace(); + if (!fes) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a finite element space"); + MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a finite element space"); + } + + mfem::Mesh *mesh = fes->GetMesh(); + if (!mesh) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh"); + MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh"); + } + + mfem::ElementTransformation *T = mesh->GetElementTransformation(zeroSlopeElemID); + if (!T) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found"); + MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found"); + } + + DEPRECATION_WARNING_OFF + u_gf->SetData(x.GetData()); + DEPRECATION_WARNING_ON + + T->SetIntPoint(&zeroSlopeIP); + mfem::Vector grad_u(3); // TODO make this a unique pointer so it can be dimensionally adaptive + u_gf->GetGradient(*T, grad_u); + + for (int i = 0; i < zeroSlopeDofs.Size(); i++) { + int dof = zeroSlopeDofs[i]; + residual[dof] -= alpha * grad_u[0]; + residual[dof] -= alpha * grad_u[1]; + residual[dof] -= alpha * grad_u[2]; + } + + } + + void ZeroSlopeNewtonSolver::ComputeConstrainedGradient(const mfem::Vector &x) const { + mfem::NonlinearForm *nlf = dynamic_cast(const_cast(oper)); + if (!nlf) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator is not a NonlinearForm"); + MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator is not a NonlinearForm"); + } + + mfem::FiniteElementSpace *fes = nlf->FESpace(); + if (!fes) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a finite element space"); + MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a finite element space"); + } + + mfem::Mesh *mesh = fes->GetMesh(); + if (!mesh) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh"); + MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: input operator does not have a mesh"); + } + + mfem::ElementTransformation *T = mesh->GetElementTransformation(zeroSlopeElemID); + if (!T) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found"); + MFEM_ABORT("ZeroSlopeNewtonSolver::ProcessNewState: element transformation is not found"); + } + + const mfem::FiniteElement* fe = fes->GetFE(zeroSlopeElemID); // Get FE *once*. + mfem::DenseMatrix dshape; // For shape function derivatives. + dshape.SetSize(fe->GetDof(), mesh->Dimension()); + T->SetIntPoint(&zeroSlopeIP); + fe->CalcDShape(zeroSlopeIP, dshape); + if (!grad) { + LOG_ERROR(logger, "ZeroSlopeNewtonSolver::ComputeConstrainedGradient: Grad is not set"); + MFEM_ABORT("ZeroSlopeNewtonSolver::ComputeConstrainedGradient: Grad is not set"); + } + // --- Modify Jacobian --- + LOG_INFO(logger, "Adjusting the Jacobian to enforce the zero slope constraint..."); + for (int i = 0; i < zeroSlopeDofs.Size(); i++) { + for (int j = 0; j < zeroSlopeDofs.Size(); j++) { + grad->Add(zeroSlopeDofs[i], zeroSlopeDofs[j], alpha * dshape(j, 0)); + grad->Add(zeroSlopeDofs[i], zeroSlopeDofs[j], alpha * dshape(j, 1)); + grad->Add(zeroSlopeDofs[i], zeroSlopeDofs[j], alpha * dshape(j, 2)); + } + } + LOG_INFO(logger, "Adjusting the Jacobian to enforce the zero slope constraint...done"); + } + } // namespace polyMFEMUtils \ No newline at end of file diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h index d65ddf2..8bb29bc 100644 --- a/src/poly/utils/public/polyMFEMUtils.h +++ b/src/poly/utils/public/polyMFEMUtils.h @@ -3,7 +3,11 @@ #include "mfem.hpp" #include +#include +#include #include "config.h" +#include "probe.h" +#include "quill/LogMacros.h" @@ -213,6 +217,33 @@ namespace polyMFEMUtils { * @return The Gaussian integral. */ double calculateGaussianIntegral(mfem::Mesh &mesh, polyMFEMUtils::GaussianCoefficient &gaussianCoeff); + + class ZeroSlopeNewtonSolver : public mfem::NewtonSolver { + private: + Config& config = Config::getInstance(); + Probe::LogManager& logManager = Probe::LogManager::getInstance(); + quill::Logger* logger = logManager.getLogger("log"); + + double alpha; // The penalty term for the flat slope at zero + std::vector zeroSlopeCoordinate; // The coordinate of the zero slope point + + int zeroSlopeElemID = -1; + mfem::Array zeroSlopeDofs; + mfem::IntegrationPoint zeroSlopeIP; + + std::unique_ptr u_gf; + mutable mfem::SparseMatrix *grad = nullptr; + + void ComputeConstrainedResidual(const mfem::Vector &x, mfem::Vector &r) const; + void ComputeConstrainedGradient(const mfem::Vector &x) const; + + public: + ZeroSlopeNewtonSolver(double alpha_, std::vector zeroSlopeCoordinate_); + ~ZeroSlopeNewtonSolver(); + // virtual void ProcessNewState(const mfem::Vector &x) const; + virtual void SetOperator(const mfem::Operator &op) override; + void Mult(const mfem::Vector &b, mfem::Vector &x) const override; + }; } // namespace polyMFEMUtils #endif // POLYMFEMUTILS_H \ No newline at end of file