v0.16.3
Loading...
Searching...
No Matches
IncrementalOptimization.cpp
Go to the documentation of this file.
1/**
2 * @file IncrementalOptimization.cpp
3 * @brief Physics-independent incremental optimization orchestration
4 */
5
6#define SINGULARITY
7#include <MoFEM.hpp>
8using namespace MoFEM;
9
10#include <IncrementalOptimization.hpp>
11
12namespace EshelbianPlasticity {
13
14MoFEMErrorCode IncrementalOptimizationContext::rollbackTrial() {
15 MoFEMFunctionBegin;
16 CHKERR VecCopy(committedControl, candidateControl);
17 CHKERR VecGhostUpdateBegin(candidateControl, INSERT_VALUES, SCATTER_FORWARD);
18 CHKERR VecGhostUpdateEnd(candidateControl, INSERT_VALUES, SCATTER_FORWARD);
19 MoFEMFunctionReturn(0);
20}
21
22boost::shared_ptr<IncrementalOptimizationContext>
23createIncrementalOptimizationContext(
24 boost::shared_ptr<IncrementalOptimizationProblem> problem,
25 SmartPetscObj<Vec> state) {
26 if (!problem || !state)
27 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
28 "Incremental optimization requires a problem and state");
29
30 CHK_THROW_MESSAGE(problem->validate(),
31 "Validate incremental-optimization problem");
32 auto context = boost::make_shared<IncrementalOptimizationContext>();
33 context->problem = std::move(problem);
34 CHK_THROW_MESSAGE(context->problem->createControlData(
35 context->committedControl, context->candidateControl),
36 "Create incremental-optimization control data");
37 if (!context->committedControl || !context->candidateControl)
38 CHK_THROW_MESSAGE(
39 MOFEM_DATA_INCONSISTENCY,
40 "Incremental-optimization problem returned incomplete control data");
41#ifndef NDEBUG
42 if (context->committedControl.get() == context->candidateControl.get())
43 CHK_THROW_MESSAGE(
44 MOFEM_DATA_INCONSISTENCY,
45 "Candidate and committed controls must be distinct vectors");
46#endif
47
48 DM control_dm = context->problem->getControlDM();
49 if (!control_dm)
50 CHK_THROW_MESSAGE(
51 MOFEM_DATA_INCONSISTENCY,
52 "Incremental-optimization problem returned no control DM");
53 const Problem *control_problem_ptr = nullptr;
54 CHK_THROW_MESSAGE(DMMoFEMGetProblemPtr(control_dm, &control_problem_ptr),
55 "Get incremental-optimization control problem");
56 if (!control_problem_ptr)
57 CHK_THROW_MESSAGE(
58 MOFEM_DATA_INCONSISTENCY,
59 "Incremental-optimization control DM has no MoFEM problem");
60
61 const MPI_Comm control_comm =
62 PetscObjectComm(reinterpret_cast<PetscObject>(control_dm));
63 context->constraintDM = createDM(control_comm, "DMMOFEM");
64 const std::string constraint_problem_name =
65 control_problem_ptr->getName() + "_CONSTRAINTS";
66 CHK_THROW_MESSAGE(DMMoFEMCreateSubDM(context->constraintDM, control_dm,
67 constraint_problem_name.c_str()),
68 "Create incremental-optimization constraint DM");
69 CHK_THROW_MESSAGE(DMMoFEMSetDestroyProblem(context->constraintDM, PETSC_TRUE),
70 "Set constraint-DM problem ownership");
71 CHK_THROW_MESSAGE(DMMoFEMSetSquareProblem(context->constraintDM, PETSC_FALSE),
72 "Set rectangular constraint-DM layout");
73 CHK_THROW_MESSAGE(
74 context->problem->configureConstraintDM(context->constraintDM),
75 "Configure incremental-optimization constraint DM");
76 CHK_THROW_MESSAGE(DMSetUp(context->constraintDM),
77 "Set up incremental-optimization constraint DM");
78
79 context->inequalityConstraints =
80 createDMVector(context->constraintDM, RowColData::ROW);
81 context->inequalityJacobian = createDMMatrix(context->constraintDM);
82 if (!context->inequalityConstraints || !context->inequalityJacobian)
83 CHK_THROW_MESSAGE(
84 MOFEM_DATA_INCONSISTENCY,
85 "Incremental-optimization problem returned incomplete inequality "
86 "data");
87
88#ifndef NDEBUG
89 PetscInt local_controls = 0;
90 PetscInt global_controls = 0;
91 PetscInt local_constraints = 0;
92 PetscInt global_constraints = 0;
93 PetscInt local_rows = 0;
94 PetscInt local_columns = 0;
95 PetscInt global_rows = 0;
96 PetscInt global_columns = 0;
97 CHK_THROW_MESSAGE(VecGetLocalSize(context->committedControl, &local_controls),
98 "Get local incremental-control size");
99 CHK_THROW_MESSAGE(VecGetSize(context->committedControl, &global_controls),
100 "Get global incremental-control size");
101 CHK_THROW_MESSAGE(
102 VecGetLocalSize(context->inequalityConstraints, &local_constraints),
103 "Get local inequality-constraint size");
104 CHK_THROW_MESSAGE(
105 VecGetSize(context->inequalityConstraints, &global_constraints),
106 "Get global inequality-constraint size");
107 CHK_THROW_MESSAGE(MatGetLocalSize(context->inequalityJacobian, &local_rows,
108 &local_columns),
109 "Get local constraint-Jacobian size");
110 CHK_THROW_MESSAGE(MatGetSize(context->inequalityJacobian, &global_rows,
111 &global_columns),
112 "Get global constraint-Jacobian size");
113 if (local_rows != local_constraints || global_rows != global_constraints ||
114 local_columns != local_controls || global_columns != global_controls)
115 CHK_THROW_MESSAGE(
116 MOFEM_DATA_INCONSISTENCY,
117 "Constraint DM is incompatible with the constraint/control vectors");
118#endif
119
120 context->gradient = vectorDuplicate(context->committedControl);
121 context->objectiveGradient = vectorDuplicate(context->committedControl);
122 context->equilibratedState = vectorDuplicate(state);
123 context->referenceControl = vectorDuplicate(context->committedControl);
124 context->baselineState = vectorDuplicate(state);
125 context->lastValidState = vectorDuplicate(state);
126 context->cachedControl = vectorDuplicate(context->committedControl);
127 context->cachedState = vectorDuplicate(state);
128
129 CHK_THROW_MESSAGE(
130 PetscOptionsGetReal(PETSC_NULLPTR, "",
131 "-incremental_optimization_equilibrium_tolerance",
132 &context->equilibriumTolerance, PETSC_NULLPTR),
133 "Read incremental-optimization equilibrium tolerance");
134 if (!(context->equilibriumTolerance > 0) ||
135 !std::isfinite(context->equilibriumTolerance))
136 CHK_THROW_MESSAGE(
137 MOFEM_INVALID_DATA,
138 "Incremental-optimization equilibrium tolerance must be positive");
139
140 CHK_THROW_MESSAGE(VecCopy(state, context->equilibratedState),
141 "Initialise equilibrated state");
142 CHK_THROW_MESSAGE(context->problem->initialiseReferenceControl(
143 context->committedControl, context->referenceControl),
144 "Initialise incremental-optimization reference control");
145 CHK_THROW_MESSAGE(VecCopy(state, context->baselineState),
146 "Initialise baseline state snapshot");
147 CHK_THROW_MESSAGE(VecCopy(state, context->lastValidState),
148 "Initialise last-valid state snapshot");
149 CHK_THROW_MESSAGE(VecCopy(state, context->cachedState),
150 "Initialise cached state snapshot");
151 context->stateSolve = [problem_ptr = context->problem](
152 Vec trial_state,
153 EquilibriumSolveReport &report) {
154 return problem_ptr->solveState(trial_state, report);
155 };
156 return context;
157}
158
159namespace {
160
161bool isFiniteSuccessfulObjective(
162 const IncrementalObjectiveEvaluation &evaluation) {
163 return evaluation.stateStatus == EquilibriumSolveStatus::Converged &&
164 std::isfinite(evaluation.conservativeValue) &&
165 std::isfinite(evaluation.dissipativeValue) &&
166 std::isfinite(evaluation.equilibriumNorm);
167}
168
169bool isFiniteSuccessfulEvaluation(
170 const IncrementalObjectiveEvaluation &evaluation) {
171 return isFiniteSuccessfulObjective(evaluation) &&
172 evaluation.smoothGradient.get() &&
173 evaluation.objectiveGradient.get() &&
174 std::isfinite(evaluation.smoothGradientNorm) &&
175 std::isfinite(evaluation.objectiveGradientNorm);
176}
177
178MoFEMErrorCode
179clearIncrementalOptimizationForces(IncrementalOptimizationContext &context) {
180 MoFEMFunctionBegin;
181 CHKERR VecZeroEntries(context.gradient);
182 CHKERR VecZeroEntries(context.objectiveGradient);
183 MoFEMFunctionReturn(0);
184}
185
186/** Restore the complete pre-solve transaction on every exit before commit. */
187struct IncrementalOptimizationSolveRollback {
188 boost::shared_ptr<IncrementalOptimizationContext> context;
189 Vec state = nullptr;
190 SmartPetscObj<Vec> stateBefore;
191 bool active = true;
192
193 IncrementalOptimizationSolveRollback(
194 boost::shared_ptr<IncrementalOptimizationContext> context, Vec state,
195 SmartPetscObj<Vec> state_before)
196 : context(std::move(context)), state(state),
197 stateBefore(std::move(state_before)) {}
198
199 MoFEMErrorCode rollback() {
200 MoFEMErrorCode first_error = 0;
201 const auto record_error = [&first_error](const MoFEMErrorCode error) {
202 if (error && !first_error)
203 first_error = error;
204 };
205 record_error(context->rollbackTrial());
206 record_error(VecCopy(stateBefore, context->equilibratedState));
207 record_error(VecCopy(stateBefore, context->baselineState));
208 record_error(VecCopy(stateBefore, context->lastValidState));
209 record_error(context->problem->restoreState(context->equilibratedState));
210 record_error(VecCopy(stateBefore, state));
211 record_error(clearIncrementalOptimizationForces(*context));
212 context->baselineValid = false;
213 context->cacheValid = false;
214 active = false;
215 return first_error;
216 }
217
218 void dismiss() { active = false; }
219
220 ~IncrementalOptimizationSolveRollback() {
221 if (active) {
222 const auto error = rollback();
223 if (error)
224 MOFEM_LOG("EP", Sev::error)
225 << "Incremental-optimization rollback failed with error "
226 << error;
227 }
228 }
229};
230
231struct IncrementalOptimizationTAOCtx : public MoFEM::TaoCtx {
232 boost::shared_ptr<IncrementalOptimizationContext> context;
233 IncrementalObjectiveEvaluation finalEvaluation;
234
235 explicit IncrementalOptimizationTAOCtx(
236 boost::shared_ptr<IncrementalOptimizationContext> context)
237 : MoFEM::TaoCtx(*getInterfacePtr(context->problem->getControlDM()),
238 "INCREMENTAL_OPTIMIZATION"),
239 context(std::move(context)) {
240 getObjectiveHook() =
241 [this](Tao tao, Vec control, PetscReal *objective) {
242 return evaluate(tao, control, objective, nullptr);
243 };
244 getGradientHook() = [this](Tao tao, Vec control, Vec gradient) {
245 MoFEMFunctionBegin;
246 PetscReal objective = 0;
247 CHKERR evaluate(tao, control, &objective, gradient);
248 // A gradient-only callback cannot return the rejected-objective sentinel.
249 if (!isFiniteSuccessfulEvaluation(finalEvaluation))
250 SETERRQ(PetscObjectComm(reinterpret_cast<PetscObject>(tao)),
251 MOFEM_OPERATION_UNSUCCESSFUL,
252 "Cannot evaluate the incremental objective gradient");
253 MoFEMFunctionReturn(0);
254 };
255 getObjectiveAndGradientHook() =
256 [this](Tao tao, Vec control, PetscReal *objective, Vec gradient) {
257 return evaluate(tao, control, objective, gradient);
258 };
259 }
260
261 static PetscErrorCode evaluateInequalityConstraints(Tao, Vec control,
262 Vec constraints,
263 void *data) {
264 MoFEMFunctionBegin;
265 auto *tao_context = static_cast<IncrementalOptimizationTAOCtx *>(data);
266 if (!tao_context)
267 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
268 "Incremental inequality callback has no TAO context");
269 // The control is vector-owned: its scalar slot is Delta kappa, whereas
270 // the corresponding mesh field stores committed kappa_n. Calling the
271 // standard MoFEM TAO wrapper here would scatter Delta kappa to the mesh
272 // before evaluating this vector-native constraint.
273 CHKERR evaluateIncrementalOptimizationInequalityConstraints(
274 tao_context->context, control, constraints);
275 MoFEMFunctionReturn(0);
276 }
277
278 static PetscErrorCode evaluateInequalityJacobian(Tao, Vec control,
279 Mat jacobian,
280 Mat preconditioner,
281 void *data) {
282 MoFEMFunctionBegin;
283 auto *tao_context = static_cast<IncrementalOptimizationTAOCtx *>(data);
284 if (!tao_context)
285 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
286 "Incremental inequality-Jacobian callback has no TAO context");
287 if (jacobian != preconditioner)
288 SETERRQ(PetscObjectComm(reinterpret_cast<PetscObject>(jacobian)),
289 MOFEM_NOT_IMPLEMENTED,
290 "Incremental inequality Jacobian and preconditioner must use "
291 "the same matrix");
292 // As above, bypass the mesh-scattering wrapper and assemble directly from
293 // the TAO control vector.
294 CHKERR assembleIncrementalOptimizationInequalityJacobian(
295 tao_context->context, control, preconditioner);
296 MoFEMFunctionReturn(0);
297 }
298
299 MoFEMErrorCode evaluate(Tao tao, Vec control, PetscReal *objective,
300 Vec gradient) {
301 MoFEMFunctionBegin;
302 finalEvaluation = {};
303 IncrementalObjectiveEvaluation evaluation;
304 CHKERR evaluateIncrementalOptimizationTAOObjectiveAndGradient(
305 context, control, *objective, gradient, evaluation);
306 finalEvaluation = evaluation;
307 if (evaluation.stateStatus != EquilibriumSolveStatus::Converged) {
308 ++context->rejectedTrialCount;
309 MoFEMFunctionReturnHot(0);
310 }
311 if (!std::isfinite(*objective) ||
312 !isFiniteSuccessfulObjective(evaluation) ||
313 (gradient && !isFiniteSuccessfulEvaluation(evaluation))) {
314 ++context->rejectedTrialCount;
315 *objective = std::numeric_limits<PetscReal>::infinity();
316 if (gradient)
317 CHKERR VecZeroEntries(gradient);
318 MoFEMFunctionReturnHot(0);
319 }
320
321 PetscInt iteration = 0;
322 CHKERR TaoGetIterationNumber(tao, &iteration);
323 const double conservative_increment =
324 evaluation.conservativeValue - context->baselineConservativeValue;
325 MOFEM_LOG("EP", Sev::inform)
326 << context->problem->getName() << " TAO evaluation at iteration "
327 << iteration << ": objective " << *objective
328 << " (conservative increment " << conservative_increment
329 << ", dissipation " << evaluation.dissipativeValue << ")"
330 << ", equilibrium residual " << evaluation.equilibriumNorm;
331 if (gradient)
332 MOFEM_LOG("EP", Sev::inform)
333 << context->problem->getName() << " TAO gradient norm "
334 << evaluation.objectiveGradientNorm;
335 MoFEMFunctionReturn(0);
336 }
337};
338
339} // namespace
340
341MoFEMErrorCode evaluateIncrementalObjective(
342 const boost::shared_ptr<IncrementalOptimizationContext> &context,
343 Vec control, IncrementalObjectiveEvaluation &evaluation) {
344 MoFEMFunctionBegin;
345 if (!context || !context->problem || !control || !context->stateSolve)
346 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
347 "Incremental-optimization evaluation has an incomplete context");
348
349 PetscBool is_reference = PETSC_FALSE;
350 CHKERR VecEqual(control, context->referenceControl, &is_reference);
351 if (!context->baselineValid && !is_reference) {
352 IncrementalObjectiveEvaluation baseline_evaluation;
353 CHKERR evaluateIncrementalObjective(
354 context, context->referenceControl, baseline_evaluation);
355 if (baseline_evaluation.stateStatus != EquilibriumSolveStatus::Converged) {
356 evaluation = baseline_evaluation;
357 MoFEMFunctionReturnHot(0);
358 }
359 }
360
361 if (context->cacheValid) {
362 PetscBool cache_match = PETSC_FALSE;
363 CHKERR VecEqual(control, context->cachedControl, &cache_match);
364 if (cache_match) {
365 CHKERR VecCopy(control, context->candidateControl);
366 CHKERR VecGhostUpdateBegin(context->candidateControl, INSERT_VALUES,
367 SCATTER_FORWARD);
368 CHKERR VecGhostUpdateEnd(context->candidateControl, INSERT_VALUES,
369 SCATTER_FORWARD);
370 CHKERR VecCopy(context->cachedState, context->equilibratedState);
371 CHKERR context->problem->restoreState(context->equilibratedState);
372 evaluation = context->cachedEvaluation;
373 evaluation.cacheHit = true;
374 evaluation.stateSolveCount = context->stateSolveCount;
375 MoFEMFunctionReturnHot(0);
376 }
377 }
378 context->cacheValid = false;
379
380 CHKERR VecCopy(control, context->candidateControl);
381 CHKERR VecGhostUpdateBegin(context->candidateControl, INSERT_VALUES,
382 SCATTER_FORWARD);
383 CHKERR VecGhostUpdateEnd(context->candidateControl, INSERT_VALUES,
384 SCATTER_FORWARD);
385 CHKERR VecCopy(context->lastValidState, context->equilibratedState);
386
387 EquilibriumSolveReport solve_report;
388 ++context->stateSolveCount;
389 PetscLogDouble equilibrium_start = 0;
390 PetscLogDouble equilibrium_end = 0;
391 CHKERR PetscTime(&equilibrium_start);
392 const MoFEMErrorCode callback_error =
393 context->stateSolve(context->equilibratedState, solve_report);
394 CHKERR PetscTime(&equilibrium_end);
395 context->equilibriumSeconds += equilibrium_end - equilibrium_start;
396 if (callback_error) {
397 solve_report.status = EquilibriumSolveStatus::Failed;
398 solve_report.errorCode = callback_error;
399 }
400 if (solve_report.status == EquilibriumSolveStatus::Converged &&
401 (!std::isfinite(solve_report.equilibriumNorm) ||
402 solve_report.equilibriumNorm > context->equilibriumTolerance)) {
403 solve_report.status = EquilibriumSolveStatus::Diverged;
404 solve_report.errorCode = MOFEM_OPERATION_UNSUCCESSFUL;
405 }
406 if (solve_report.status != EquilibriumSolveStatus::Converged) {
407 CHKERR VecCopy(context->lastValidState, context->equilibratedState);
408 CHKERR context->problem->restoreState(context->equilibratedState);
409 CHKERR context->rollbackTrial();
410 evaluation = {};
411 evaluation.stateStatus = solve_report.status;
412 evaluation.stateSolveError = solve_report.errorCode;
413 evaluation.equilibriumNorm = solve_report.equilibriumNorm;
414 evaluation.stateSolveCount = context->stateSolveCount;
415 evaluation.smoothGradient = context->gradient;
416 evaluation.objectiveGradient = context->objectiveGradient;
417 CHKERR clearIncrementalOptimizationForces(*context);
418 context->cacheValid = false;
419 MoFEMFunctionReturnHot(0);
420 }
421
422 double conservative_value = 0;
423 CHKERR context->problem->evaluateConservativeValue(conservative_value);
424 double dissipative_value = 0;
425 CHKERR context->problem->evaluateDissipation(context->candidateControl,
426 dissipative_value);
427 if (!std::isfinite(conservative_value) ||
428 !std::isfinite(dissipative_value)) {
429 CHKERR VecCopy(context->lastValidState, context->equilibratedState);
430 CHKERR context->problem->restoreState(context->equilibratedState);
431 CHKERR context->rollbackTrial();
432 evaluation = {};
433 evaluation.stateStatus = EquilibriumSolveStatus::Failed;
434 evaluation.stateSolveError = MOFEM_INVALID_DATA;
435 evaluation.equilibriumNorm = solve_report.equilibriumNorm;
436 evaluation.stateSolveCount = context->stateSolveCount;
437 evaluation.smoothGradient = context->gradient;
438 evaluation.objectiveGradient = context->objectiveGradient;
439 CHKERR clearIncrementalOptimizationForces(*context);
440 context->cacheValid = false;
441 MoFEMFunctionReturnHot(0);
442 }
443
444 evaluation = {};
445 evaluation.conservativeValue = conservative_value;
446 evaluation.dissipativeValue = dissipative_value;
447 evaluation.equilibriumNorm = solve_report.equilibriumNorm;
448 evaluation.stateStatus = solve_report.status;
449 evaluation.stateSolveError = solve_report.errorCode;
450 evaluation.stateSolveCount = context->stateSolveCount;
451
452 CHKERR VecCopy(context->equilibratedState, context->lastValidState);
453 if (is_reference && !context->baselineValid) {
454 context->baselineConservativeValue = conservative_value;
455 CHKERR VecCopy(context->equilibratedState, context->baselineState);
456 context->baselineValid = true;
457 }
458 CHKERR VecCopy(control, context->cachedControl);
459 CHKERR VecCopy(context->equilibratedState, context->cachedState);
460 context->cachedEvaluation = evaluation;
461 context->cachedEvaluation.cacheHit = false;
462 context->cacheValid = true;
463
464 MoFEMFunctionReturn(0);
465}
466
467MoFEMErrorCode evaluateIncrementalObjectiveAndGradient(
468 const boost::shared_ptr<IncrementalOptimizationContext> &context,
469 Vec control, IncrementalObjectiveEvaluation &evaluation) {
470 MoFEMFunctionBegin;
471 CHKERR evaluateIncrementalObjective(context, control, evaluation);
472 if (evaluation.stateStatus != EquilibriumSolveStatus::Converged ||
473 evaluation.smoothGradient.get())
474 MoFEMFunctionReturnHot(0);
475
476 PetscLogDouble gradient_start = 0;
477 PetscLogDouble gradient_end = 0;
478 CHKERR PetscTime(&gradient_start);
479 CHKERR context->problem->assembleConservativeGradient(context->gradient);
480 CHKERR PetscTime(&gradient_end);
481 context->gradientAssemblySeconds += gradient_end - gradient_start;
482 PetscReal gradient_norm = 0;
483 CHKERR VecNorm(context->gradient, NORM_2, &gradient_norm);
484 if (!std::isfinite(gradient_norm)) {
485 CHKERR context->rollbackTrial();
486 CHKERR clearIncrementalOptimizationForces(*context);
487 context->cacheValid = false;
488 evaluation.stateStatus = EquilibriumSolveStatus::Failed;
489 evaluation.stateSolveError = MOFEM_INVALID_DATA;
490 MoFEMFunctionReturnHot(0);
491 }
492 evaluation.smoothGradient = context->gradient;
493 evaluation.smoothGradientNorm = gradient_norm;
494 context->cachedEvaluation = evaluation;
495 MoFEMFunctionReturn(0);
496}
497
498MoFEMErrorCode evaluateIncrementalOptimizationTAOObjectiveAndGradient(
499 const boost::shared_ptr<IncrementalOptimizationContext> &context,
500 Vec control, PetscReal &objective, Vec objective_gradient,
501 IncrementalObjectiveEvaluation &evaluation) {
502 MoFEMFunctionBegin;
503 if (objective_gradient)
504 CHKERR evaluateIncrementalObjectiveAndGradient(context, control, evaluation);
505 else
506 CHKERR evaluateIncrementalObjective(context, control, evaluation);
507 if (evaluation.stateStatus != EquilibriumSolveStatus::Converged) {
508 objective = std::numeric_limits<PetscReal>::infinity();
509 if (objective_gradient)
510 CHKERR VecZeroEntries(objective_gradient);
511 MoFEMFunctionReturnHot(0);
512 }
513
514 objective = evaluation.conservativeValue - context->baselineConservativeValue +
515 evaluation.dissipativeValue;
516 if (!std::isfinite(objective)) {
517 if (objective_gradient)
518 CHKERR VecZeroEntries(objective_gradient);
519 MoFEMFunctionReturnHot(0);
520 }
521 if (!objective_gradient)
522 MoFEMFunctionReturnHot(0);
523
524 // TAO may overwrite its gradient workspace (e.g. with ALMM penalties).
525 // Cache the adjoint gradient only; this local combination is inexpensive.
526 CHKERR context->problem->assembleObjectiveGradient(
527 context->candidateControl, evaluation.smoothGradient,
528 context->objectiveGradient);
529 PetscReal gradient_norm = 0;
530 CHKERR VecNorm(context->objectiveGradient, NORM_2, &gradient_norm);
531 evaluation.objectiveGradient = context->objectiveGradient;
532 evaluation.objectiveGradientNorm = gradient_norm;
533 if (!std::isfinite(gradient_norm)) {
534 objective = std::numeric_limits<PetscReal>::infinity();
535 CHKERR context->rollbackTrial();
536 CHKERR clearIncrementalOptimizationForces(*context);
537 context->cacheValid = false;
538 evaluation.stateStatus = EquilibriumSolveStatus::Failed;
539 evaluation.stateSolveError = MOFEM_INVALID_DATA;
540 }
541 if (objective_gradient != context->objectiveGradient.get())
542 CHKERR VecCopy(context->objectiveGradient, objective_gradient);
543 MoFEMFunctionReturn(0);
544}
545
546MoFEMErrorCode evaluateIncrementalOptimizationInequalityConstraints(
547 const boost::shared_ptr<IncrementalOptimizationContext> &context,
548 Vec control, Vec constraints) {
549 MoFEMFunctionBegin;
550 if (!context || !context->problem || !control || !constraints)
551 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
552 "Cannot evaluate incomplete incremental inequalities");
553 CHKERR VecZeroEntries(constraints);
554 CHKERR context->problem->evaluateInequalityConstraints(control,
555 constraints);
556 MoFEMFunctionReturn(0);
557}
558
559MoFEMErrorCode assembleIncrementalOptimizationInequalityJacobian(
560 const boost::shared_ptr<IncrementalOptimizationContext> &context,
561 Vec control, Mat jacobian) {
562 MoFEMFunctionBegin;
563 if (!context || !context->problem || !control || !jacobian)
564 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
565 "Cannot assemble incomplete incremental inequality Jacobian");
566 IncrementalObjectiveEvaluation evaluation;
567 if (context->cacheValid) {
568 PetscBool cache_match = PETSC_FALSE;
569 CHKERR VecEqual(control, context->cachedControl, &cache_match);
570 if (cache_match)
571 evaluation = context->cachedEvaluation;
572 }
573 CHKERR assembleIncrementalOptimizationInequalityJacobianFromEvaluation(
574 context, control, evaluation, jacobian);
575 MoFEMFunctionReturn(0);
576}
577
578MoFEMErrorCode assembleIncrementalOptimizationInequalityJacobianFromEvaluation(
579 const boost::shared_ptr<IncrementalOptimizationContext> &context,
580 Vec control, const IncrementalObjectiveEvaluation &evaluation,
581 Mat jacobian) {
582 MoFEMFunctionBegin;
583 if (!context || !context->problem || !control || !jacobian)
584 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
585 "Cannot assemble an inequality Jacobian from incomplete "
586 "objective data");
587 // A cached conservative gradient selects the exact-apex subgradient. Away
588 // from the apex the plastic Jacobian depends only on the control.
589 const bool has_finite_objective_data =
590 isFiniteSuccessfulObjective(evaluation) &&
591 evaluation.smoothGradient.get() &&
592 std::isfinite(evaluation.smoothGradientNorm);
593#ifndef NDEBUG
594 if (has_finite_objective_data) {
595 PetscBool matches_candidate = PETSC_FALSE;
596 CHKERR VecEqual(control, context->candidateControl, &matches_candidate);
597 if (!matches_candidate)
598 SETERRQ(PetscObjectComm(reinterpret_cast<PetscObject>(control)),
599 MOFEM_DATA_INCONSISTENCY,
600 "Constraint Jacobian evaluation does not match the current "
601 "equilibrated control");
602 }
603#endif
604 CHKERR MatZeroEntries(jacobian);
605 if (has_finite_objective_data) {
606 CHKERR context->problem->assembleInequalityJacobian(
607 control, evaluation.smoothGradient.get(), jacobian);
608 } else
609 CHKERR context->problem
610 ->assembleInequalityJacobianWithoutSmoothGradient(control, jacobian);
611 MoFEMFunctionReturn(0);
612}
613
614MoFEMErrorCode solveIncrementalOptimizationTAO(
615 const boost::shared_ptr<IncrementalOptimizationContext> &context,
616 Vec state) {
617 MoFEMFunctionBegin;
618 if (!context || !context->problem || !state)
619 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
620 "Cannot solve an incomplete incremental-optimization problem");
621
622 // Duplicate the context state so the first copy also preflights that the
623 // caller's state has the layout needed by rollback and finalization.
624 auto state_before = vectorDuplicate(context->equilibratedState);
625 CHKERR VecCopy(state, state_before);
626 IncrementalOptimizationSolveRollback rollback_guard(context, state,
627 state_before);
628
629 CHKERR VecCopy(context->referenceControl, context->candidateControl);
630 auto tao_context = boost::make_shared<IncrementalOptimizationTAOCtx>(context);
631 DM control_dm = context->problem->getControlDM();
632 const MPI_Comm control_comm =
633 PetscObjectComm(reinterpret_cast<PetscObject>(control_dm));
634 auto tao = createTao(control_comm);
635 auto tao_solution = vectorDuplicate(context->referenceControl);
636 CHKERR VecCopy(context->referenceControl, tao_solution);
637 CHKERR TaoSetOptionsPrefix(tao, "incremental_optimization_");
638 CHKERR TaoRegisterNullSpace();
639 CHKERR TaoSetType(tao, TAOALMM);
640 Tao subsolver = nullptr;
641 CHKERR TaoALMMGetSubsolver(tao, &subsolver);
642 CHKERR TaoSetType(subsolver, TAOBQNKTR);
643 CHKERR TaoSetInitialTrustRegionRadius(subsolver, 1e-3);
644 CHKERR TaoSetTolerances(tao, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
645 CHKERR TaoSetConstraintTolerances(tao, 1e-8, PETSC_DEFAULT);
646 CHKERR ::TaoSetObjective(tao, MoFEM::TaoSetObjective, tao_context.get());
647 CHKERR ::TaoSetGradient(tao, context->objectiveGradient,
648 MoFEM::TaoSetGradient, tao_context.get());
649 CHKERR ::TaoSetObjectiveAndGradient(
650 tao, context->objectiveGradient, MoFEM::TaoSetObjectiveAndGradient,
651 tao_context.get());
652 CHKERR TaoSetInequalityConstraintsRoutine(
653 tao, context->inequalityConstraints,
654 IncrementalOptimizationTAOCtx::evaluateInequalityConstraints,
655 tao_context.get());
656 CHKERR TaoSetJacobianInequalityRoutine(
657 tao, context->inequalityJacobian, context->inequalityJacobian,
658 IncrementalOptimizationTAOCtx::evaluateInequalityJacobian,
659 tao_context.get());
660 CHKERR TaoSetFromOptions(tao);
661 PetscBool is_almm = PETSC_FALSE;
662 PetscBool is_nullspace = PETSC_FALSE;
663 CHKERR PetscObjectTypeCompare(reinterpret_cast<PetscObject>(tao.get()),
664 TAOALMM, &is_almm);
665 CHKERR PetscObjectTypeCompare(reinterpret_cast<PetscObject>(tao.get()),
666 TAONULLSPACE, &is_nullspace);
667 if (!is_almm && !is_nullspace)
668 SETERRQ(control_comm, MOFEM_NOT_IMPLEMENTED,
669 "Incremental constrained optimization currently requires "
670 "-incremental_optimization_tao_type almm or nullspace");
671 if (is_nullspace)
672 CHKERR TaoNullSpaceSetConstraints(
673 tao, context->constraintDM, context->inequalityConstraints,
674 context->inequalityJacobian, context->inequalityJacobian);
675 CHKERR TaoSetSolution(tao, tao_solution);
676 CHKERR TaoSetUp(tao);
677 Vec tao_multipliers = nullptr;
678 if (is_almm) {
679 TaoALMMType almm_type;
680 CHKERR TaoALMMGetType(tao, &almm_type);
681 if (almm_type != TAO_ALMM_PHR)
682 SETERRQ(control_comm, MOFEM_NOT_IMPLEMENTED,
683 "Incremental constrained optimization currently requires "
684 "-incremental_optimization_tao_almm_type phr so that returned "
685 "inequality multipliers use the physical nonnegative sign");
686 CHKERR TaoALMMGetMultipliers(tao, &tao_multipliers);
687 } else {
688 CHKERR TaoNullSpaceGetMultipliers(tao, &tao_multipliers);
689 }
690 CHKERR VecZeroEntries(tao_multipliers);
691 CHKERR context->problem->initialiseInequalityMultipliers(
692 context->constraintDM, tao_multipliers);
693 // TaoSolve_ALMM clears its internal multipliers unless recycle is enabled.
694 // Preserve the physics-specific eta-stationarity initialization above.
695 if (is_almm)
696 CHKERR TaoSetRecycleHistory(tao, PETSC_TRUE);
697 CHKERR TaoSolve(tao);
698 TaoConvergedReason reason;
699 CHKERR TaoGetConvergedReason(tao, &reason);
700 const bool reached_iteration_limit = reason == TAO_DIVERGED_MAXITS;
701 if (reached_iteration_limit) {
702 MOFEM_LOG("EP", Sev::warning)
703 << context->problem->getName()
704 << " incremental optimization reached the TAO iteration limit; "
705 "checking the final iterate before continuing";
706 } else if (reason <= 0) {
707 SETERRQ(control_comm, MOFEM_OPERATION_UNSUCCESSFUL,
708 "Incremental-optimization TAO did not converge: %s",
709 TaoConvergedReasons[reason]);
710 }
711
712 Vec solution = nullptr;
713 CHKERR TaoGetSolution(tao, &solution);
714 PetscReal final_objective = 0;
715 CHKERR tao_context->evaluate(tao, solution, &final_objective,
716 context->objectiveGradient);
717 if (tao_context->finalEvaluation.stateStatus !=
718 EquilibriumSolveStatus::Converged)
719 SETERRQ(control_comm, MOFEM_OPERATION_UNSUCCESSFUL,
720 "Final TAO control does not have an equilibrated state");
721 if (!std::isfinite(final_objective))
722 SETERRQ(control_comm, MOFEM_OPERATION_UNSUCCESSFUL,
723 "Final TAO control has a non-finite objective");
724 if (!isFiniteSuccessfulEvaluation(tao_context->finalEvaluation))
725 SETERRQ(control_comm, MOFEM_OPERATION_UNSUCCESSFUL,
726 "Final TAO control has incomplete or non-finite objective data");
727
728#ifndef NDEBUG
729 PetscBool candidate_matches_solution = PETSC_FALSE;
730 CHKERR VecEqual(context->candidateControl, solution,
731 &candidate_matches_solution);
732 if (!candidate_matches_solution)
733 SETERRQ(control_comm, MOFEM_DATA_INCONSISTENCY,
734 "Final equilibrated control does not match the TAO solution");
735#endif
736
737 PetscReal gradient_absolute_tolerance = 0;
738 CHKERR TaoGetTolerances(tao, &gradient_absolute_tolerance, PETSC_NULLPTR,
739 PETSC_NULLPTR);
740 PetscReal constraint_absolute_tolerance = 0;
741 CHKERR TaoGetConstraintTolerances(tao, &constraint_absolute_tolerance,
742 PETSC_NULLPTR);
743 Vec inequality_multipliers = nullptr;
744 if (is_almm)
745 CHKERR TaoALMMGetMultipliers(tao, &inequality_multipliers);
746 else
747 CHKERR TaoNullSpaceGetMultipliers(tao, &inequality_multipliers);
748 if (!inequality_multipliers)
749 SETERRQ(control_comm, MOFEM_DATA_INCONSISTENCY,
750 "TAO did not provide final inequality multipliers");
751 std::string diagnostics;
752 CHKERR context->problem->validateSolution(
753 context->constraintDM, context->candidateControl, context->gradient,
754 inequality_multipliers, gradient_absolute_tolerance,
755 constraint_absolute_tolerance, diagnostics);
756#ifndef NDEBUG
757 PetscBool validated_control_matches_solution = PETSC_FALSE;
758 CHKERR VecEqual(context->candidateControl, solution,
759 &validated_control_matches_solution);
760 if (!validated_control_matches_solution)
761 SETERRQ(control_comm, MOFEM_DATA_INCONSISTENCY,
762 "Validated control no longer matches the TAO solution");
763#endif
764 auto accepted_solution = vectorDuplicate(context->candidateControl);
765 CHKERR VecCopy(context->candidateControl, accepted_solution);
766 CHKERR VecGhostUpdateBegin(accepted_solution, INSERT_VALUES,
767 SCATTER_FORWARD);
768 CHKERR VecGhostUpdateEnd(accepted_solution, INSERT_VALUES, SCATTER_FORWARD);
769 // Finalize generic state while rollback is still possible. The physics
770 // commit below owns the remaining history/control transition atomically.
771 CHKERR VecCopy(context->equilibratedState, state);
772 CHKERR VecCopy(context->equilibratedState, context->baselineState);
773 CHKERR VecCopy(context->equilibratedState, context->lastValidState);
774 const MoFEMErrorCode commit_error = context->problem->commit(
775 accepted_solution, context->equilibratedState,
776 context->committedControl, context->candidateControl,
777 context->referenceControl);
778 if (commit_error)
779 CHKERR commit_error;
780 context->baselineValid = false;
781 context->cacheValid = false;
782 rollback_guard.dismiss();
783
784 MOFEM_LOG("EP", reached_iteration_limit ? Sev::warning : Sev::inform)
785 << context->problem->getName() << " incremental optimization "
786 << (reached_iteration_limit ? "accepted at the TAO iteration limit: "
787 : "converged: ")
788 << "equilibrium " << tao_context->finalEvaluation.equilibriumNorm
789 << ", TAO gradient "
790 << tao_context->finalEvaluation.objectiveGradientNorm << ", objective "
791 << final_objective << ", state evaluations " << context->stateSolveCount
792 << ", rejected trials " << context->rejectedTrialCount
793 << ", seconds equilibrium/gradient = " << context->equilibriumSeconds
794 << "/" << context->gradientAssemblySeconds
795 << (diagnostics.empty() ? "" : ", ") << diagnostics;
796 MoFEMFunctionReturn(0);
797}
798
799} // namespace EshelbianPlasticity
finalEvaluation
IncrementalObjectiveEvaluation finalEvaluation
Definition IncrementalOptimization.cpp:233
stateBefore
SmartPetscObj< Vec > stateBefore
Definition IncrementalOptimization.cpp:190
context
boost::shared_ptr< IncrementalOptimizationContext > context
Definition IncrementalOptimization.cpp:188
IncrementalOptimization.hpp
Physics-independent incremental optimization API.
ROW
@ ROW
Definition definitions.h:136
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition definitions.h:36
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition DMMoFEM.cpp:422
MoFEM::createDMVector
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1237
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1194
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
EshelbianPlasticity::createIncrementalOptimizationContext
boost::shared_ptr< IncrementalOptimizationContext > createIncrementalOptimizationContext(boost::shared_ptr< IncrementalOptimizationProblem > problem, SmartPetscObj< Vec > state)
Definition IncrementalOptimization.cpp:23
EshelbianPlasticity::solveIncrementalOptimizationTAO
MoFEMErrorCode solveIncrementalOptimizationTAO(const boost::shared_ptr< IncrementalOptimizationContext > &context, Vec state)
Definition IncrementalOptimization.cpp:614
EshelbianPlasticity::evaluateIncrementalObjectiveAndGradient
MoFEMErrorCode evaluateIncrementalObjectiveAndGradient(const boost::shared_ptr< IncrementalOptimizationContext > &context, Vec control, IncrementalObjectiveEvaluation &evaluation)
Definition IncrementalOptimization.cpp:467
EshelbianPlasticity::assembleIncrementalOptimizationInequalityJacobianFromEvaluation
MoFEMErrorCode assembleIncrementalOptimizationInequalityJacobianFromEvaluation(const boost::shared_ptr< IncrementalOptimizationContext > &context, Vec control, const IncrementalObjectiveEvaluation &evaluation, Mat jacobian)
Definition IncrementalOptimization.cpp:578
EshelbianPlasticity::evaluateIncrementalOptimizationTAOObjectiveAndGradient
MoFEMErrorCode evaluateIncrementalOptimizationTAOObjectiveAndGradient(const boost::shared_ptr< IncrementalOptimizationContext > &context, Vec control, PetscReal &objective, Vec objective_gradient, IncrementalObjectiveEvaluation &evaluation)
Definition IncrementalOptimization.cpp:498
EshelbianPlasticity::evaluateIncrementalObjective
MoFEMErrorCode evaluateIncrementalObjective(const boost::shared_ptr< IncrementalOptimizationContext > &context, Vec control, IncrementalObjectiveEvaluation &evaluation)
Definition IncrementalOptimization.cpp:341
EshelbianPlasticity::assembleIncrementalOptimizationInequalityJacobian
MoFEMErrorCode assembleIncrementalOptimizationInequalityJacobian(const boost::shared_ptr< IncrementalOptimizationContext > &context, Vec control, Mat jacobian)
Definition IncrementalOptimization.cpp:559
EshelbianPlasticity::evaluateIncrementalOptimizationInequalityConstraints
MoFEMErrorCode evaluateIncrementalOptimizationInequalityConstraints(const boost::shared_ptr< IncrementalOptimizationContext > &context, Vec control, Vec constraints)
Definition IncrementalOptimization.cpp:546
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition DMMoFEM.cpp:434
MoFEM::TaoSetObjective
PetscErrorCode TaoSetObjective(Tao tao, Vec x, PetscReal *f, void *ctx)
Sets the objective function value for a TAO optimization context.
Definition TaoCtx.cpp:37
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition DeprecatedPetsc.hpp:152
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
MoFEM::getInterfacePtr
auto getInterfacePtr(DM dm)
Get the Interface Ptr object.
Definition DMMoFEM.hpp:1171
MoFEM::TaoSetObjectiveAndGradient
PetscErrorCode TaoSetObjectiveAndGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
Sets the objective function value and gradient for a TAO optimization solver.
Definition TaoCtx.cpp:178
MoFEM::TaoSetGradient
PetscErrorCode TaoSetGradient(Tao tao, Vec x, Vec f, void *ctx)
Sets the gradient vector for a TAO optimization context.
Definition TaoCtx.cpp:92
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition PetscSmartObj.hpp:143
MoFEM::createTao
auto createTao(MPI_Comm comm)
Definition PetscSmartObj.hpp:251
MoFEM::Problem
keeps basic data about problem
Definition ProblemsMultiIndices.hpp:54
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
MoFEM::TaoCtx
Interface for TAO solvers.
Definition TaoCtx.hpp:14
Generated by Doxygen 1.9.8 and hosted at University of Glasgow