PolyFEM
Loading...
Searching...
No Matches
ALSolver.cpp
Go to the documentation of this file.
1#include "ALSolver.hpp"
2
4
5namespace polyfem::solver
6{
8 std::shared_ptr<BCLagrangianForm> lagr_form,
9 std::shared_ptr<BCPenaltyForm> pen_form,
10 const double initial_al_weight,
11 const double scaling,
12 const double max_al_weight,
13 const double eta_tol,
14 const std::function<void(const Eigen::VectorXd &)> &update_barrier_stiffness)
15 : lagr_form(lagr_form),
16 pen_form(pen_form),
17 initial_al_weight(initial_al_weight),
18 scaling(scaling),
19 max_al_weight(max_al_weight),
20 eta_tol(eta_tol),
21 update_barrier_stiffness(update_barrier_stiffness)
22 {
23 }
24
25 void ALSolver::solve_al(std::shared_ptr<NLSolver> nl_solver, NLProblem &nl_problem, Eigen::MatrixXd &sol)
26 {
27 assert(sol.size() == nl_problem.full_size());
28
29 const Eigen::VectorXd initial_sol = sol;
30 Eigen::VectorXd tmp_sol = nl_problem.full_to_reduced(sol);
31 assert(tmp_sol.size() == nl_problem.reduced_size());
32
33 // --------------------------------------------------------------------
34
35 double al_weight = initial_al_weight;
36 int al_steps = 0;
37 const int iters = nl_solver->stop_criteria().iterations;
38
39 const double initial_error = pen_form->compute_error(sol);
40
41 nl_problem.line_search_begin(sol, tmp_sol);
42
43 while (!std::isfinite(nl_problem.value(tmp_sol))
44 || !nl_problem.is_step_valid(sol, tmp_sol)
45 || !nl_problem.is_step_collision_free(sol, tmp_sol))
46 {
47 nl_problem.line_search_end();
48
49 set_al_weight(nl_problem, sol, al_weight);
50 logger().debug("Solving AL Problem with weight {}", al_weight);
51
52 nl_problem.init(sol);
54 tmp_sol = sol;
55
56 try
57 {
58 nl_solver->minimize(nl_problem, tmp_sol);
59 }
60 catch (const std::runtime_error &e)
61 {
62 }
63
64 sol = tmp_sol;
65 set_al_weight(nl_problem, sol, -1);
66
67 const double current_error = pen_form->compute_error(sol);
68 const double eta = 1 - sqrt(current_error / initial_error);
69
70 logger().debug("Current eta = {}", eta);
71
72 if (eta < 0)
73 {
74 logger().debug("Higher error than initial, increase weight and revert to previous solution");
75 sol = initial_sol;
76 }
77
78 tmp_sol = nl_problem.full_to_reduced(sol);
79 nl_problem.line_search_begin(sol, tmp_sol);
80
81 if (eta < eta_tol && al_weight < max_al_weight)
82 al_weight *= scaling;
83 else
84 lagr_form->update_lagrangian(sol, al_weight);
85
86 post_subsolve(al_weight);
87 ++al_steps;
88 }
89 nl_problem.line_search_end();
90 nl_solver->stop_criteria().iterations = iters;
91 }
92
93 void ALSolver::solve_reduced(std::shared_ptr<NLSolver> nl_solver, NLProblem &nl_problem, Eigen::MatrixXd &sol)
94 {
95 assert(sol.size() == nl_problem.full_size());
96
97 Eigen::VectorXd tmp_sol = nl_problem.full_to_reduced(sol);
98 nl_problem.line_search_begin(sol, tmp_sol);
99
100 if (!std::isfinite(nl_problem.value(tmp_sol))
101 || !nl_problem.is_step_valid(sol, tmp_sol)
102 || !nl_problem.is_step_collision_free(sol, tmp_sol))
103 log_and_throw_error("Failed to apply boundary conditions; solve with augmented lagrangian first!");
104
105 // --------------------------------------------------------------------
106 // Perform one final solve with the DBC projected out
107
108 logger().debug("Successfully applied boundary conditions; solving in reduced space");
109
110 nl_problem.init(sol);
112 try
113 {
114 nl_solver->minimize(nl_problem, tmp_sol);
115 }
116 catch (const std::runtime_error &e)
117 {
118 sol = nl_problem.reduced_to_full(tmp_sol);
119 throw e;
120 }
121 sol = nl_problem.reduced_to_full(tmp_sol);
122
123 post_subsolve(0);
124 }
125
126 void ALSolver::set_al_weight(NLProblem &nl_problem, const Eigen::VectorXd &x, const double weight)
127 {
128 if (pen_form == nullptr || lagr_form == nullptr)
129 return;
130 if (weight > 0)
131 {
132 pen_form->enable();
133 lagr_form->enable();
134 pen_form->set_weight(weight);
135 nl_problem.use_full_size();
136 nl_problem.set_apply_DBC(x, false);
137 }
138 else
139 {
140 pen_form->disable();
141 lagr_form->disable();
142 nl_problem.use_reduced_size();
143 nl_problem.set_apply_DBC(x, true);
144 }
145 }
146
147} // namespace polyfem::solver
int x
const double max_al_weight
Definition ALSolver.hpp:43
void set_al_weight(NLProblem &nl_problem, const Eigen::VectorXd &x, const double weight)
Definition ALSolver.cpp:126
std::function< void(const Eigen::VectorXd &)> update_barrier_stiffness
Definition ALSolver.hpp:47
ALSolver(std::shared_ptr< BCLagrangianForm > lagr_form, std::shared_ptr< BCPenaltyForm > pen_form, const double initial_al_weight, const double scaling, const double max_al_weight, const double eta_tol, const std::function< void(const Eigen::VectorXd &)> &update_barrier_stiffness)
Definition ALSolver.cpp:7
std::shared_ptr< BCPenaltyForm > pen_form
Definition ALSolver.hpp:40
std::shared_ptr< BCLagrangianForm > lagr_form
Definition ALSolver.hpp:39
std::function< void(const double)> post_subsolve
Definition ALSolver.hpp:34
const double initial_al_weight
Definition ALSolver.hpp:41
void solve_reduced(std::shared_ptr< NLSolver > nl_solver, NLProblem &nl_problem, Eigen::MatrixXd &sol)
Definition ALSolver.cpp:93
void solve_al(std::shared_ptr< NLSolver > nl_solver, NLProblem &nl_problem, Eigen::MatrixXd &sol)
Definition ALSolver.cpp:25
virtual void line_search_end() override
virtual void init(const TVector &x0) override
void line_search_begin(const TVector &x0, const TVector &x1) override
Definition NLProblem.cpp:80
virtual bool is_step_valid(const TVector &x0, const TVector &x1) override
Definition NLProblem.cpp:90
virtual TVector full_to_reduced(const TVector &full) const
virtual bool is_step_collision_free(const TVector &x0, const TVector &x1) override
Definition NLProblem.cpp:95
virtual TVector reduced_to_full(const TVector &reduced) const
virtual double value(const TVector &x) override
void set_apply_DBC(const TVector &x, const bool val)
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:42
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:71