Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AdjointNLProblem.cpp
Go to the documentation of this file.
2
9#include <polyfem/State.hpp>
11
12#include <list>
13#include <stack>
14
15namespace polyfem::solver
16{
17 namespace
18 {
19
20 Eigen::VectorXd get_updated_mesh_nodes(const VariableToSimulationGroup &variables_to_simulation, const std::shared_ptr<State> &curr_state, const Eigen::VectorXd &x)
21 {
22 Eigen::MatrixXd V;
23 curr_state->get_vertices(V);
24 Eigen::VectorXd X = utils::flatten(V);
25
26 for (auto &p : variables_to_simulation)
27 {
28 for (const auto &state : p->get_states())
29 if (state.get() != curr_state.get())
30 continue;
31 if (p->get_parameter_type() != ParameterType::Shape)
32 continue;
33 auto state_variable = p->get_parametrization().eval(x);
34 auto output_indexing = p->get_output_indexing(x);
35 for (int i = 0; i < output_indexing.size(); ++i)
36 X(output_indexing(i)) = state_variable(i);
37 }
38
39 return X;
40 }
41
42 using namespace std;
43 // Class to represent a graph
44 class Graph
45 {
46 int V; // No. of vertices'
47
48 // adjacency lists
49 vector<list<int>> adj;
50
51 // A function used by topologicalSort
52 void topologicalSortUtil(int v, vector<bool> &visited, stack<int> &Stack);
53
54 public:
55 Graph(int V); // Constructor
56
57 // function to add an edge to graph
58 void addEdge(int v, int w);
59
60 // prints a Topological Sort of the complete graph
61 vector<int> topologicalSort();
62 };
63
64 Graph::Graph(int V)
65 {
66 this->V = V;
67 adj.resize(V);
68 }
69
70 void Graph::addEdge(int v, int w)
71 {
72 adj[v].push_back(w); // Add w to v’s list.
73 }
74
75 // A recursive function used by topologicalSort
76 void Graph::topologicalSortUtil(int v, vector<bool> &visited,
77 stack<int> &Stack)
78 {
79 // Mark the current node as visited.
80 visited[v] = true;
81
82 // Recur for all the vertices adjacent to this vertex
83 list<int>::iterator i;
84 for (i = adj[v].begin(); i != adj[v].end(); ++i)
85 if (!visited[*i])
86 topologicalSortUtil(*i, visited, Stack);
87
88 // Push current vertex to stack which stores result
89 Stack.push(v);
90 }
91
92 // The function to do Topological Sort. It uses recursive
93 // topologicalSortUtil()
94 vector<int> Graph::topologicalSort()
95 {
96 stack<int> Stack;
97
98 // Mark all the vertices as not visited
99 vector<bool> visited(V, false);
100
101 // Call the recursive helper function to store Topological
102 // Sort starting from all vertices one by one
103 for (int i = 0; i < V; i++)
104 if (visited[i] == false)
105 topologicalSortUtil(i, visited, Stack);
106
107 // Print contents of stack
108 vector<int> sorted;
109 while (Stack.empty() == false)
110 {
111 sorted.push_back(Stack.top());
112 Stack.pop();
113 }
114
115 return sorted;
116 }
117 } // namespace
118
119 AdjointNLProblem::AdjointNLProblem(std::shared_ptr<AdjointForm> form, const VariableToSimulationGroup &variables_to_simulation, const std::vector<std::shared_ptr<State>> &all_states, const json &args)
120 : FullNLProblem({form}),
121 form_(form),
122 variables_to_simulation_(variables_to_simulation),
123 all_states_(all_states),
124 save_freq(args["output"]["save_frequency"]),
125 enable_slim(args["solver"]["advanced"]["enable_slim"]),
126 smooth_line_search(args["solver"]["advanced"]["smooth_line_search"]),
127 solve_in_parallel(args["solver"]["advanced"]["solve_in_parallel"])
128 {
129 cur_grad.setZero(0);
130
131 if (enable_slim && args["solver"]["nonlinear"]["advanced"]["apply_gradient_fd"] != "None")
132 adjoint_logger().warn("SLIM may affect the finite difference result!");
133
134 if (enable_slim && smooth_line_search)
135 adjoint_logger().warn("Both in-line-search SLIM and after-line-search SLIM are ON!");
136
137 if (args["output"]["solution"] != "")
138 {
139 solution_ostream.open(args["output"]["solution"].get<std::string>(), std::ofstream::out);
140 if (!solution_ostream.is_open())
141 adjoint_logger().error("Cannot open solution file for writing!");
142 }
143
144 solve_in_order.clear();
145 {
146 Graph G(all_states.size());
147 for (int k = 0; k < all_states.size(); k++)
148 {
149 auto &arg = args["states"][k];
150 if (arg["initial_guess"].get<int>() >= 0)
151 G.addEdge(arg["initial_guess"].get<int>(), k);
152 }
153
154 solve_in_order = G.topologicalSort();
155 }
156
157 active_state_mask.assign(all_states_.size(), false);
158 for (int i = 0; i < all_states_.size(); i++)
159 {
160 for (const auto &v2sim : variables_to_simulation_)
161 {
162 for (const auto &state : v2sim->get_states())
163 {
164 if (all_states_[i].get() == state.get())
165 {
166 active_state_mask[i] = true;
167 break;
168 }
169 }
170 }
171 }
172 }
173
174 AdjointNLProblem::AdjointNLProblem(std::shared_ptr<AdjointForm> form, const std::vector<std::shared_ptr<AdjointForm>> &stopping_conditions, const VariableToSimulationGroup &variables_to_simulation, const std::vector<std::shared_ptr<State>> &all_states, const json &args) : AdjointNLProblem(form, variables_to_simulation, all_states, args)
175 {
176 stopping_conditions_ = stopping_conditions;
177 }
178
179 void AdjointNLProblem::hessian(const Eigen::VectorXd &x, StiffnessMatrix &hessian)
180 {
181 log_and_throw_adjoint_error("Hessian not supported!");
182 }
183
184 double AdjointNLProblem::value(const Eigen::VectorXd &x)
185 {
186 return form_->value(x);
187 }
188
189 void AdjointNLProblem::gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv)
190 {
191 if (cur_grad.size() == x.size())
192 gradv = cur_grad;
193 else
194 {
195 gradv.setZero(x.size());
196
197 {
198 POLYFEM_SCOPED_TIMER("adjoint solve");
199 for (int i = 0; i < all_states_.size(); i++)
200 all_states_[i]->solve_adjoint_cached(form_->compute_reduced_adjoint_rhs(x, *all_states_[i])); // caches inside state
201 }
202
203 {
204 POLYFEM_SCOPED_TIMER("gradient assembly");
205 form_->first_derivative(x, gradv);
206 if (x.size() < 10)
207 {
208 adjoint_logger().trace("x {}", x.transpose());
209 adjoint_logger().trace("gradient {}", gradv.transpose());
210 }
211 }
212
213 cur_grad = gradv;
214 }
215 }
216
217 bool AdjointNLProblem::is_step_valid(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1)
218 {
219 bool need_rebuild_basis = false;
220
221 // update to new parameter and check if the new parameter is valid to solve
222 for (const auto &v : variables_to_simulation_)
223 if (v->get_parameter_type() == ParameterType::Shape || v->get_parameter_type() == ParameterType::PeriodicShape)
224 need_rebuild_basis = true;
225
226 if (need_rebuild_basis && smooth_line_search)
227 {
228 Eigen::MatrixXd X, V0, V1;
229 Eigen::MatrixXi F;
230
231 int state_count = 0;
232 for (auto state_ : all_states_)
233 {
234
236 get_updated_mesh_nodes(variables_to_simulation_, state_, x1),
237 state_->mesh->dimension());
238 state_->get_vertices(V0);
239 state_->get_elements(F);
240
241 Eigen::MatrixXd V_smooth;
242 bool slim_success = polyfem::mesh::apply_slim(V0, F, V1, V_smooth);
243 if (!slim_success)
244 {
245 adjoint_logger().info("SLIM failed, step not valid!");
246 return false;
247 }
248
249 V1 = V_smooth;
250
251 bool flipped = AdjointTools::is_flipped(V1, F);
252 if (flipped)
253 {
254 adjoint_logger().info("Found flipped element in LS, step not valid!");
255 return false;
256 }
257
258 state_count++;
259 }
260 }
261
262 return form_->is_step_valid(x0, x1);
263 }
264
265 bool AdjointNLProblem::is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1)
266 {
267 return form_->is_step_collision_free(x0, x1);
268 }
269
270 double AdjointNLProblem::max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1)
271 {
272 return form_->max_step_size(x0, x1);
273 }
274
275 void AdjointNLProblem::line_search_begin(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1)
276 {
277 form_->line_search_begin(x0, x1);
278 }
279
281 {
282 form_->line_search_end();
283 }
284
285 void AdjointNLProblem::post_step(const polysolve::nonlinear::PostStepData &data)
286 {
287 save_to_file(save_iter++, data.x);
288
289 form_->post_step(data);
290 }
291
292 void AdjointNLProblem::save_to_file(const int iter_num, const Eigen::VectorXd &x0)
293 {
294 int id = 0;
295
296 if (solution_ostream.is_open())
297 {
298 adjoint_logger().debug("Save solution at iteration {} to file...", iter_num);
299 solution_ostream << iter_num << ": " << std::setprecision(16) << x0.transpose() << std::endl;
300 solution_ostream.flush();
301 }
302
303 if (iter_num % save_freq != 0)
304 return;
305 adjoint_logger().info("Saving iteration {}", iter_num);
306 for (const auto &state : all_states_)
307 {
308 bool save_vtu = true;
309 bool save_rest_mesh = true;
310
311 std::string vis_mesh_path = state->resolve_output_path(fmt::format("opt_state_{:d}_iter_{:d}.vtu", id, iter_num));
312 std::string mesh_ext = state->mesh->is_volume() ? ".msh" : ".obj";
313 std::string rest_mesh_path = state->resolve_output_path(fmt::format("opt_state_{:d}_iter_{:d}" + mesh_ext, id, iter_num));
314 id++;
315
316 if (!save_vtu)
317 continue;
318 adjoint_logger().debug("Save final vtu to file {} ...", vis_mesh_path);
319
320 double tend = state->args.value("tend", 1.0);
321 double dt = 1;
322 if (!state->args["time"].is_null())
323 dt = state->args["time"]["dt"];
324
325 Eigen::MatrixXd sol = state->diff_cached.u(-1);
326
327 state->out_geom.save_vtu(
328 vis_mesh_path,
329 *state,
330 sol,
331 Eigen::MatrixXd::Zero(state->n_pressure_bases, 1),
332 tend, dt,
333 io::OutGeometryData::ExportOptions(state->args, state->mesh->is_linear(), state->problem->is_scalar()), state->is_contact_enabled());
334
335 if (!save_rest_mesh)
336 continue;
337 adjoint_logger().debug("Save rest mesh to file {} ...", rest_mesh_path);
338
339 // If shape opt, save rest meshes as well
340 Eigen::MatrixXd V;
341 Eigen::MatrixXi F;
342 state->get_vertices(V);
343 state->get_elements(F);
344 if (state->mesh->is_volume())
345 io::MshWriter::write(rest_mesh_path, V, F, state->mesh->get_body_ids(), true, false);
346 else
347 io::OBJWriter::write(rest_mesh_path, V, F);
348 }
349 }
350
351 void AdjointNLProblem::solution_changed(const Eigen::VectorXd &newX)
352 {
353 bool need_rebuild_basis = false;
354
355 // update to new parameter and check if the new parameter is valid to solve
356 for (const auto &v : variables_to_simulation_)
357 {
358 v->update(newX);
359 if (v->get_parameter_type() == ParameterType::Shape || v->get_parameter_type() == ParameterType::PeriodicShape)
360 need_rebuild_basis = true;
361 }
362
363 if (need_rebuild_basis)
364 {
365 for (const auto &state : all_states_)
366 state->build_basis();
367 }
368
369 // solve PDE
370 solve_pde();
371
372 form_->solution_changed(newX);
373
374 curr_x = newX;
375 }
376
377 bool AdjointNLProblem::after_line_search_custom_operation(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1)
378 {
379 if (!enable_slim)
380 return false;
381
382 for (const auto &v : variables_to_simulation_)
383 v->update(x0);
384
385 std::vector<Eigen::MatrixXd> V_old_list;
386 for (auto state : all_states_)
387 {
388 Eigen::MatrixXd V;
389 state->get_vertices(V);
390 V_old_list.push_back(V);
391 }
392
393 for (const auto &v : variables_to_simulation_)
394 v->update(x1);
395
396 // Apply slim to all states on a frequency
397 int state_num = 0;
398 for (auto state : all_states_)
399 {
400 Eigen::MatrixXd V_new, V_smooth;
401 Eigen::MatrixXi F;
402 state->get_vertices(V_new);
403 state->get_elements(F);
404
405 bool slim_success = polyfem::mesh::apply_slim(V_old_list[state_num++], F, V_new, V_smooth, 50);
406
407 if (!slim_success)
408 log_and_throw_adjoint_error("SLIM cannot succeed");
409
410 for (int i = 0; i < V_smooth.rows(); ++i)
411 state->set_mesh_vertex(i, V_smooth.row(i));
412 }
413 adjoint_logger().debug("SLIM succeeded!");
414
415 return true;
416 }
417
419 {
421 {
422 adjoint_logger().info("Run simulations in parallel...");
423
424 utils::maybe_parallel_for(all_states_.size(), [&](int start, int end, int thread_id) {
425 for (int i = start; i < end; i++)
426 {
427 auto state = all_states_[i];
428 if (active_state_mask[i] || state->diff_cached.size() == 0)
429 {
430 state->assemble_rhs();
431 state->assemble_mass_mat();
432 Eigen::MatrixXd sol, pressure; // solution is also cached in state
433 state->solve_problem(sol, pressure);
434 }
435 }
436 });
437 }
438 else
439 {
440 adjoint_logger().info("Run simulations in serial...");
441
442 Eigen::MatrixXd sol, pressure; // solution is also cached in state
443 for (int i : solve_in_order)
444 {
445 auto state = all_states_[i];
446 if (active_state_mask[i] || state->diff_cached.size() == 0)
447 {
448 state->assemble_rhs();
449 state->assemble_mass_mat();
450
451 state->solve_problem(sol, pressure);
452 }
453 }
454 }
455
456 cur_grad.resize(0);
457 }
458
459 bool AdjointNLProblem::stop(const TVector &x)
460 {
461 if (stopping_conditions_.size() == 0)
462 return false;
463
464 for (auto &obj : stopping_conditions_)
465 {
466 obj->solution_changed(x);
467 if (obj->value(x) > 0)
468 return false;
469 }
470 return true;
471 }
472
473} // namespace polyfem::solver
vector< list< int > > adj
int V
int x
#define POLYFEM_SCOPED_TIMER(...)
Definition Timer.hpp:10
static void write(const std::string &path, const mesh::Mesh &mesh, const bool binary)
saves the mesh
Definition MshWriter.cpp:7
static bool write(const std::string &path, const Eigen::MatrixXd &v, const Eigen::MatrixXi &e, const Eigen::MatrixXi &f)
Definition OBJWriter.cpp:18
bool after_line_search_custom_operation(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
double max_step_size(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
void gradient(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) override
void save_to_file(const int iter_num, const Eigen::VectorXd &x0)
bool stop(const TVector &x) override
bool is_step_valid(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
std::vector< std::shared_ptr< AdjointForm > > stopping_conditions_
void hessian(const Eigen::VectorXd &x, StiffnessMatrix &hessian) override
void post_step(const polysolve::nonlinear::PostStepData &data) override
AdjointNLProblem(std::shared_ptr< AdjointForm > form, const VariableToSimulationGroup &variables_to_simulation, const std::vector< std::shared_ptr< State > > &all_states, const json &args)
void solution_changed(const Eigen::VectorXd &new_x) override
bool is_step_collision_free(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
double value(const Eigen::VectorXd &x) override
std::vector< std::shared_ptr< State > > all_states_
std::shared_ptr< AdjointForm > form_
const VariableToSimulationGroup variables_to_simulation_
void line_search_begin(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1) override
A collection of VariableToSimulation.
bool apply_slim(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXd &V_new, Eigen::MatrixXd &V_smooth, const int max_iters)
bool is_flipped(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F)
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
Eigen::VectorXd flatten(const Eigen::MatrixXd &X)
Flatten rowwises.
void maybe_parallel_for(int size, const std::function< void(int, int, int)> &partial_for)
spdlog::logger & adjoint_logger()
Retrieves the current logger for adjoint.
Definition Logger.cpp:28
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_adjoint_error(const std::string &msg)
Definition Logger.cpp:77
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22