Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "Moose.h"
11 : #include "MooseError.h"
12 : #include "OptimizeSolve.h"
13 : #include "OptimizationAppTypes.h"
14 : #include "OptimizationReporterBase.h"
15 : #include "Steady.h"
16 :
17 : InputParameters
18 998 : OptimizeSolve::validParams()
19 : {
20 998 : InputParameters params = emptyInputParameters();
21 : MooseEnum tao_solver_enum(
22 : "taontr taobntr taobncg taonls taobnls taobqnktr taontl taobntl taolmvm "
23 1996 : "taoblmvm taonm taobqnls taoowlqn taogpcg taobmrm taoalmm");
24 1996 : params.addRequiredParam<MooseEnum>(
25 : "tao_solver", tao_solver_enum, "Tao solver to use for optimization.");
26 998 : ExecFlagEnum exec_enum = ExecFlagEnum();
27 2994 : exec_enum.addAvailableFlags(EXEC_NONE,
28 : OptimizationAppTypes::EXEC_FORWARD,
29 : OptimizationAppTypes::EXEC_ADJOINT,
30 : OptimizationAppTypes::EXEC_HOMOGENEOUS_FORWARD);
31 : exec_enum = {OptimizationAppTypes::EXEC_FORWARD,
32 : OptimizationAppTypes::EXEC_ADJOINT,
33 4990 : OptimizationAppTypes::EXEC_HOMOGENEOUS_FORWARD};
34 1996 : params.addParam<ExecFlagEnum>(
35 : "solve_on", exec_enum, "List of flags indicating when inner system solve should occur.");
36 2994 : params.addParam<bool>(
37 : "output_optimization_iterations",
38 1996 : false,
39 : "Use the time step as the current iteration for outputting optimization history.");
40 998 : return params;
41 1996 : }
42 :
43 499 : OptimizeSolve::OptimizeSolve(Executioner & ex)
44 : : SolveObject(ex),
45 998 : _my_comm(MPI_COMM_SELF),
46 998 : _solve_on(getParam<ExecFlagEnum>("solve_on")),
47 499 : _verbose(getParam<bool>("verbose")),
48 998 : _output_opt_iters(getParam<bool>("output_optimization_iterations")),
49 998 : _tao_solver_enum(getParam<MooseEnum>("tao_solver").getEnum<TaoSolverEnum>()),
50 1996 : _parameters(std::make_unique<libMesh::PetscVector<Number>>(_my_comm))
51 : {
52 499 : if (libMesh::n_threads() > 1)
53 0 : mooseError("OptimizeSolve does not currently support threaded execution");
54 :
55 499 : if (_output_opt_iters && _problem.isTransient())
56 0 : mooseDocumentedError(
57 : "moose", 27225, "Outputting for transient executioners has not been implemented.");
58 499 : }
59 :
60 : bool
61 499 : OptimizeSolve::solve()
62 : {
63 998 : TIME_SECTION("optimizeSolve", 1, "Optimization Solve");
64 : // Initial solve
65 499 : _inner_solve->solve();
66 :
67 : // Grab objective function
68 998 : if (!_problem.hasUserObject("OptimizationReporter"))
69 0 : mooseError("No OptimizationReporter object found.");
70 499 : _obj_function = &_problem.getUserObject<OptimizationReporterBase>("OptimizationReporter");
71 :
72 : // Initialize solution and matrix
73 499 : _obj_function->setInitialCondition(*_parameters);
74 491 : _ndof = _parameters->size();
75 :
76 : // time step defaults 1, we want to start at 0 for first iteration to be
77 : // consistent with TAO iterations.
78 491 : if (_output_opt_iters)
79 29 : _problem.timeStep() = 0;
80 491 : bool solveInfo = (taoSolve() == 0);
81 485 : return solveInfo;
82 : }
83 :
84 : PetscErrorCode
85 491 : OptimizeSolve::taoSolve()
86 : {
87 : PetscFunctionBegin;
88 : // Initialize tao object
89 491 : LibmeshPetscCallQ(TaoCreate(_my_comm.get(), &_tao));
90 :
91 : #if PETSC_RELEASE_LESS_THAN(3, 21, 0)
92 : LibmeshPetscCallQ(TaoSetMonitor(_tao, monitor, this, nullptr));
93 : #else
94 491 : LibmeshPetscCallQ(TaoMonitorSet(_tao, monitor, this, nullptr));
95 : #endif
96 :
97 491 : switch (_tao_solver_enum)
98 : {
99 0 : case TaoSolverEnum::NEWTON_TRUST_REGION:
100 0 : LibmeshPetscCallQ(TaoSetType(_tao, TAONTR));
101 : break;
102 0 : case TaoSolverEnum::BOUNDED_NEWTON_TRUST_REGION:
103 0 : LibmeshPetscCallQ(TaoSetType(_tao, TAOBNTR));
104 : break;
105 46 : case TaoSolverEnum::BOUNDED_CONJUGATE_GRADIENT:
106 46 : LibmeshPetscCallQ(TaoSetType(_tao, TAOBNCG));
107 : break;
108 45 : case TaoSolverEnum::NEWTON_LINE_SEARCH:
109 45 : LibmeshPetscCallQ(TaoSetType(_tao, TAONLS));
110 : break;
111 0 : case TaoSolverEnum::BOUNDED_NEWTON_LINE_SEARCH:
112 0 : LibmeshPetscCallQ(TaoSetType(_tao, TAOBNLS));
113 : break;
114 101 : case TaoSolverEnum::BOUNDED_QUASI_NEWTON_TRUST_REGION:
115 101 : LibmeshPetscCallQ(TaoSetType(_tao, TAOBQNKTR));
116 : break;
117 0 : case TaoSolverEnum::NEWTON_TRUST_LINE:
118 0 : LibmeshPetscCallQ(TaoSetType(_tao, TAONTL));
119 : break;
120 0 : case TaoSolverEnum::BOUNDED_NEWTON_TRUST_LINE:
121 0 : LibmeshPetscCallQ(TaoSetType(_tao, TAOBNTL));
122 : break;
123 117 : case TaoSolverEnum::QUASI_NEWTON:
124 117 : LibmeshPetscCallQ(TaoSetType(_tao, TAOLMVM));
125 : break;
126 72 : case TaoSolverEnum::BOUNDED_QUASI_NEWTON:
127 72 : LibmeshPetscCallQ(TaoSetType(_tao, TAOBLMVM));
128 : break;
129 :
130 28 : case TaoSolverEnum::NELDER_MEAD:
131 28 : LibmeshPetscCallQ(TaoSetType(_tao, TAONM));
132 : break;
133 :
134 54 : case TaoSolverEnum::BOUNDED_QUASI_NEWTON_LINE_SEARCH:
135 54 : LibmeshPetscCallQ(TaoSetType(_tao, TAOBQNLS));
136 : break;
137 0 : case TaoSolverEnum::ORTHANT_QUASI_NEWTON:
138 0 : LibmeshPetscCallQ(TaoSetType(_tao, TAOOWLQN));
139 : break;
140 0 : case TaoSolverEnum::GRADIENT_PROJECTION_CONJUGATE_GRADIENT:
141 0 : LibmeshPetscCallQ(TaoSetType(_tao, TAOGPCG));
142 : break;
143 0 : case TaoSolverEnum::BUNDLE_RISK_MIN:
144 0 : LibmeshPetscCallQ(TaoSetType(_tao, TAOBMRM));
145 : break;
146 28 : case TaoSolverEnum::AUGMENTED_LAGRANGIAN_MULTIPLIER_METHOD:
147 : #if !PETSC_VERSION_LESS_THAN(3, 15, 0)
148 28 : LibmeshPetscCallQ(TaoSetType(_tao, TAOALMM));
149 : // Need to cancel monitors for ALMM, if not there is a segfault at MOOSE destruction. Setup
150 : // default constraint monitor.
151 : #if PETSC_RELEASE_GREATER_EQUALS(3, 21, 0)
152 28 : LibmeshPetscCallQ(TaoMonitorCancel(_tao));
153 : #else
154 : LibmeshPetscCallQ(TaoCancelMonitors(_tao));
155 : #endif
156 28 : LibmeshPetscCallQ(PetscOptionsSetValue(NULL, "-tao_cmonitor", NULL));
157 : break;
158 : #else
159 : mooseError("ALMM is only compatible with PETSc versions above 3.14. ");
160 : #endif
161 :
162 0 : default:
163 0 : mooseError("Invalid Tao solve type");
164 : }
165 :
166 : // Set objective and gradient functions
167 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
168 491 : LibmeshPetscCallQ(TaoSetObjective(_tao, objectiveFunctionWrapper, this));
169 : #else
170 : LibmeshPetscCallQ(TaoSetObjectiveRoutine(_tao, objectiveFunctionWrapper, this));
171 : #endif
172 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
173 491 : LibmeshPetscCallQ(
174 : TaoSetObjectiveAndGradient(_tao, NULL, objectiveAndGradientFunctionWrapper, this));
175 : #else
176 : LibmeshPetscCallQ(
177 : TaoSetObjectiveAndGradientRoutine(_tao, objectiveAndGradientFunctionWrapper, this));
178 : #endif
179 :
180 : // Set matrix-free version of the Hessian function
181 491 : LibmeshPetscCallQ(MatCreateShell(_my_comm.get(), _ndof, _ndof, _ndof, _ndof, this, &_hessian));
182 : // Link matrix-free Hessian to Tao
183 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
184 491 : LibmeshPetscCallQ(TaoSetHessian(_tao, _hessian, _hessian, hessianFunctionWrapper, this));
185 : #else
186 : LibmeshPetscCallQ(TaoSetHessianRoutine(_tao, _hessian, _hessian, hessianFunctionWrapper, this));
187 : #endif
188 :
189 : // Set initial guess
190 : #if !PETSC_VERSION_LESS_THAN(3, 17, 0)
191 491 : LibmeshPetscCallQ(TaoSetSolution(_tao, _parameters->vec()));
192 : #else
193 : LibmeshPetscCallQ(TaoSetInitialVector(_tao, _parameters->vec()));
194 : #endif
195 :
196 : // Set TAO petsc options
197 491 : LibmeshPetscCallQ(TaoSetFromOptions(_tao));
198 :
199 : // save nonTAO PETSC options to reset before every call to execute()
200 491 : _petsc_options = _problem.getPetscOptions();
201 : // We only use a single system solve at this point
202 491 : _solver_params = _problem.solverParams(0);
203 :
204 : // Set bounds for bounded optimization
205 491 : LibmeshPetscCallQ(TaoSetVariableBoundsRoutine(_tao, variableBoundsWrapper, this));
206 :
207 491 : if (_tao_solver_enum == TaoSolverEnum::AUGMENTED_LAGRANGIAN_MULTIPLIER_METHOD)
208 28 : LibmeshPetscCallQ(taoALCreate());
209 :
210 : // Backup multiapps so transient problems start with the same initial condition
211 491 : _problem.backupMultiApps(OptimizationAppTypes::EXEC_FORWARD);
212 491 : _problem.backupMultiApps(OptimizationAppTypes::EXEC_ADJOINT);
213 491 : _problem.backupMultiApps(OptimizationAppTypes::EXEC_HOMOGENEOUS_FORWARD);
214 :
215 : // Solve optimization
216 491 : LibmeshPetscCallQ(TaoSolve(_tao));
217 :
218 : // Print solve statistics
219 970 : if (getParam<bool>("verbose"))
220 458 : LibmeshPetscCallQ(TaoView(_tao, PETSC_VIEWER_STDOUT_WORLD));
221 :
222 485 : LibmeshPetscCallQ(TaoDestroy(&_tao));
223 :
224 485 : LibmeshPetscCallQ(MatDestroy(&_hessian));
225 :
226 485 : if (_tao_solver_enum == TaoSolverEnum::AUGMENTED_LAGRANGIAN_MULTIPLIER_METHOD)
227 28 : LibmeshPetscCallQ(taoALDestroy());
228 :
229 485 : PetscFunctionReturn(PETSC_SUCCESS);
230 : }
231 :
232 : void
233 2556 : OptimizeSolve::getTaoSolutionStatus(std::vector<int> & tot_iters,
234 : std::vector<double> & gnorm,
235 : std::vector<int> & obj_iters,
236 : std::vector<double> & cnorm,
237 : std::vector<int> & grad_iters,
238 : std::vector<double> & xdiff,
239 : std::vector<int> & hess_iters,
240 : std::vector<double> & f,
241 : std::vector<int> & tot_solves) const
242 : {
243 : const auto num = _total_iterate_vec.size();
244 2556 : tot_iters.resize(num);
245 2556 : obj_iters.resize(num);
246 2556 : grad_iters.resize(num);
247 2556 : hess_iters.resize(num);
248 2556 : tot_solves.resize(num);
249 2556 : f.resize(num);
250 2556 : gnorm.resize(num);
251 2556 : cnorm.resize(num);
252 2556 : xdiff.resize(num);
253 :
254 7659 : for (const auto i : make_range(num))
255 : {
256 5103 : tot_iters[i] = _total_iterate_vec[i];
257 5103 : obj_iters[i] = _obj_iterate_vec[i];
258 5103 : grad_iters[i] = _grad_iterate_vec[i];
259 5103 : hess_iters[i] = _hess_iterate_vec[i];
260 5103 : tot_solves[i] = _function_solve_vec[i];
261 5103 : f[i] = _f_vec[i];
262 5103 : gnorm[i] = _gnorm_vec[i];
263 5103 : cnorm[i] = _cnorm_vec[i];
264 5103 : xdiff[i] = _xdiff_vec[i];
265 : }
266 2556 : }
267 :
268 : void
269 5900 : OptimizeSolve::setTaoSolutionStatus(double f, int its, double gnorm, double cnorm, double xdiff)
270 : {
271 : // set data from TAO
272 5900 : _total_iterate_vec.push_back(its);
273 5900 : _f_vec.push_back(f);
274 5900 : _gnorm_vec.push_back(gnorm);
275 5900 : _cnorm_vec.push_back(cnorm);
276 5900 : _xdiff_vec.push_back(xdiff);
277 : // set data we collect on this optimization iteration and then reset for next iteration
278 5900 : _obj_iterate_vec.push_back(_obj_iterate);
279 5900 : _grad_iterate_vec.push_back(_grad_iterate);
280 5900 : _hess_iterate_vec.push_back(_hess_iterate);
281 : // count total number of FE solves
282 5900 : int solves = _obj_iterate + _grad_iterate + 2 * _hess_iterate;
283 5900 : _function_solve_vec.push_back(solves);
284 5900 : _obj_iterate = 0;
285 5900 : _grad_iterate = 0;
286 5900 : _hess_iterate = 0;
287 :
288 : // Pass down the iteration number if the subapp is of the Steady/SteadyAndAdjoint type.
289 : // This enables exodus per-iteration output.
290 11879 : for (auto & sub_app : _app.getExecutioner()->feProblem().getMultiAppWarehouse().getObjects())
291 : {
292 5979 : if (auto steady = dynamic_cast<Steady *>(sub_app->getExecutioner(0)))
293 4017 : steady->setIterationNumberOutput((unsigned int)its);
294 : }
295 :
296 : // Output the converged iteration outputs
297 5900 : _problem.outputStep(OptimizationAppTypes::EXEC_FORWARD);
298 :
299 : // Increment timestep. In steady problems timestep = time for outputting.
300 : // See Output.C
301 5900 : if (_output_opt_iters)
302 54 : _problem.timeStep() += 1;
303 :
304 : // print verbose per iteration output
305 5900 : if (_verbose)
306 5279 : _console << "TAO SOLVER: iteration=" << its << "\tf=" << f << "\tgnorm=" << gnorm
307 5279 : << "\tcnorm=" << cnorm << "\txdiff=" << xdiff << std::endl;
308 5900 : }
309 :
310 : PetscErrorCode
311 5900 : OptimizeSolve::monitor(Tao tao, void * ctx)
312 : {
313 : TaoConvergedReason reason;
314 : PetscInt its;
315 : PetscReal f, gnorm, cnorm, xdiff;
316 :
317 : PetscFunctionBegin;
318 5900 : LibmeshPetscCallQ(TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason));
319 :
320 : auto * solver = static_cast<OptimizeSolve *>(ctx);
321 5900 : solver->setTaoSolutionStatus((double)f, (int)its, (double)gnorm, (double)cnorm, (double)xdiff);
322 :
323 5900 : PetscFunctionReturn(PETSC_SUCCESS);
324 : }
325 :
326 : PetscErrorCode
327 5740 : OptimizeSolve::objectiveFunctionWrapper(Tao /*tao*/, Vec x, Real * objective, void * ctx)
328 : {
329 : PetscFunctionBegin;
330 : auto * solver = static_cast<OptimizeSolve *>(ctx);
331 :
332 5740 : libMesh::PetscVector<Number> param(x, solver->_my_comm);
333 5740 : solver->_parameters->swap(param);
334 :
335 5740 : (*objective) = solver->objectiveFunction();
336 5740 : solver->_parameters->swap(param);
337 5740 : PetscFunctionReturn(PETSC_SUCCESS);
338 5740 : }
339 :
340 : PetscErrorCode
341 4816 : OptimizeSolve::objectiveAndGradientFunctionWrapper(
342 : Tao /*tao*/, Vec x, Real * objective, Vec gradient, void * ctx)
343 : {
344 : PetscFunctionBegin;
345 : auto * solver = static_cast<OptimizeSolve *>(ctx);
346 :
347 4816 : libMesh::PetscVector<Number> param(x, solver->_my_comm);
348 4816 : solver->_parameters->swap(param);
349 :
350 4816 : (*objective) = solver->objectiveFunction();
351 4812 : libMesh::PetscVector<Number> grad(gradient, solver->_my_comm);
352 4812 : solver->gradientFunction(grad);
353 4810 : solver->_parameters->swap(param);
354 4810 : PetscFunctionReturn(PETSC_SUCCESS);
355 4810 : }
356 :
357 : PetscErrorCode
358 45 : OptimizeSolve::hessianFunctionWrapper(
359 : Tao /*tao*/, Vec /*x*/, Mat /*hessian*/, Mat /*pc*/, void * ctx)
360 : {
361 : PetscFunctionBegin;
362 : // Define Hessian-vector multiplication routine
363 : auto * solver = static_cast<OptimizeSolve *>(ctx);
364 45 : LibmeshPetscCallQ(MatShellSetOperation(
365 : solver->_hessian, MATOP_MULT, (void (*)(void))OptimizeSolve::applyHessianWrapper));
366 45 : PetscFunctionReturn(PETSC_SUCCESS);
367 : }
368 :
369 : PetscErrorCode
370 117 : OptimizeSolve::applyHessianWrapper(Mat H, Vec s, Vec Hs)
371 : {
372 : void * ctx;
373 :
374 : PetscFunctionBegin;
375 117 : LibmeshPetscCallQ(MatShellGetContext(H, &ctx));
376 :
377 117 : auto * solver = static_cast<OptimizeSolve *>(ctx);
378 117 : libMesh::PetscVector<Number> sbar(s, solver->_my_comm);
379 117 : libMesh::PetscVector<Number> Hsbar(Hs, solver->_my_comm);
380 234 : return solver->applyHessian(sbar, Hsbar);
381 117 : }
382 :
383 : PetscErrorCode
384 301 : OptimizeSolve::variableBoundsWrapper(Tao tao, Vec /*xl*/, Vec /*xu*/, void * ctx)
385 : {
386 : PetscFunctionBegin;
387 : auto * solver = static_cast<OptimizeSolve *>(ctx);
388 :
389 301 : LibmeshPetscCallQ(solver->variableBounds(tao));
390 301 : PetscFunctionReturn(PETSC_SUCCESS);
391 : }
392 :
393 : Real
394 10556 : OptimizeSolve::objectiveFunction()
395 : {
396 21112 : TIME_SECTION("objectiveFunction", 2, "Objective forward solve");
397 10556 : _obj_function->updateParameters(*_parameters);
398 :
399 10556 : Moose::PetscSupport::petscSetOptions(_petsc_options, _solver_params);
400 10556 : _problem.execute(OptimizationAppTypes::EXEC_FORWARD);
401 :
402 10556 : _problem.restoreMultiApps(OptimizationAppTypes::EXEC_FORWARD);
403 21108 : if (!_problem.execMultiApps(OptimizationAppTypes::EXEC_FORWARD))
404 : {
405 : // We do this so we can output for failed solves.
406 0 : _problem.outputStep(OptimizationAppTypes::EXEC_FORWARD);
407 0 : mooseError("Forward solve multiapp failed!");
408 : }
409 10552 : if (_solve_on.isValueSet(OptimizationAppTypes::EXEC_FORWARD))
410 7230 : _inner_solve->solve();
411 :
412 10552 : _obj_iterate++;
413 21104 : return _obj_function->computeObjective();
414 : }
415 :
416 : void
417 4812 : OptimizeSolve::gradientFunction(libMesh::PetscVector<Number> & gradient)
418 : {
419 9624 : TIME_SECTION("gradientFunction", 2, "Gradient adjoint solve");
420 4812 : _obj_function->updateParameters(*_parameters);
421 :
422 4812 : Moose::PetscSupport::petscSetOptions(_petsc_options, _solver_params);
423 4812 : _problem.execute(OptimizationAppTypes::EXEC_ADJOINT);
424 4812 : _problem.restoreMultiApps(OptimizationAppTypes::EXEC_ADJOINT);
425 9624 : if (!_problem.execMultiApps(OptimizationAppTypes::EXEC_ADJOINT))
426 0 : mooseError("Adjoint solve multiapp failed!");
427 4812 : if (_solve_on.isValueSet(OptimizationAppTypes::EXEC_ADJOINT))
428 4017 : _inner_solve->solve();
429 :
430 4812 : _grad_iterate++;
431 4812 : _obj_function->computeGradient(gradient);
432 4810 : }
433 :
434 : PetscErrorCode
435 117 : OptimizeSolve::applyHessian(libMesh::PetscVector<Number> & s, libMesh::PetscVector<Number> & Hs)
436 : {
437 : PetscFunctionBegin;
438 234 : TIME_SECTION("applyHessian", 2, "Hessian forward/adjoint solve");
439 : // What happens for material inversion when the Hessian
440 : // is dependent on the parameters? Deal with it later???
441 : // see notes on how this needs to change for Material inversion
442 234 : if (_problem.hasMultiApps() &&
443 234 : !_problem.hasMultiApps(OptimizationAppTypes::EXEC_HOMOGENEOUS_FORWARD))
444 0 : mooseError("Hessian based optimization algorithms require a sub-app with:\n"
445 : " execute_on = HOMOGENEOUS_FORWARD");
446 117 : _obj_function->updateParameters(s);
447 :
448 117 : Moose::PetscSupport::petscSetOptions(_petsc_options, _solver_params);
449 117 : _problem.execute(OptimizationAppTypes::EXEC_HOMOGENEOUS_FORWARD);
450 117 : _problem.restoreMultiApps(OptimizationAppTypes::EXEC_HOMOGENEOUS_FORWARD);
451 234 : if (!_problem.execMultiApps(OptimizationAppTypes::EXEC_HOMOGENEOUS_FORWARD))
452 0 : mooseError("Homogeneous forward solve multiapp failed!");
453 117 : if (_solve_on.isValueSet(OptimizationAppTypes::EXEC_HOMOGENEOUS_FORWARD))
454 117 : _inner_solve->solve();
455 :
456 117 : _obj_function->setMisfitToSimulatedValues();
457 :
458 117 : Moose::PetscSupport::petscSetOptions(_petsc_options, _solver_params);
459 117 : _problem.execute(OptimizationAppTypes::EXEC_ADJOINT);
460 117 : _problem.restoreMultiApps(OptimizationAppTypes::EXEC_ADJOINT);
461 234 : if (!_problem.execMultiApps(OptimizationAppTypes::EXEC_ADJOINT))
462 0 : mooseError("Adjoint solve multiapp failed!");
463 117 : if (_solve_on.isValueSet(OptimizationAppTypes::EXEC_ADJOINT))
464 117 : _inner_solve->solve();
465 :
466 117 : _obj_function->computeGradient(Hs);
467 117 : _hess_iterate++;
468 117 : PetscFunctionReturn(PETSC_SUCCESS);
469 : }
470 :
471 : PetscErrorCode
472 301 : OptimizeSolve::variableBounds(Tao tao)
473 : {
474 : PetscFunctionBegin;
475 301 : unsigned int sz = _obj_function->getNumParams();
476 :
477 301 : libMesh::PetscVector<Number> xl(_my_comm, sz);
478 301 : libMesh::PetscVector<Number> xu(_my_comm, sz);
479 :
480 : // copy values from upper and lower bounds to xl and xu
481 1402 : for (const auto i : make_range(sz))
482 : {
483 1101 : xl.set(i, _obj_function->getLowerBound(i));
484 1101 : xu.set(i, _obj_function->getUpperBound(i));
485 : }
486 : // set upper and lower bounds in tao solver
487 301 : LibmeshPetscCallQ(TaoSetVariableBounds(tao, xl.vec(), xu.vec()));
488 301 : PetscFunctionReturn(PETSC_SUCCESS);
489 301 : }
490 :
491 : PetscErrorCode
492 627 : OptimizeSolve::equalityFunctionWrapper(Tao /*tao*/, Vec /*x*/, Vec ce, void * ctx)
493 : {
494 : PetscFunctionBegin;
495 : // grab the solver
496 : auto * solver = static_cast<OptimizeSolve *>(ctx);
497 627 : libMesh::PetscVector<Number> eq_con(ce, solver->_my_comm);
498 : // use the OptimizationReporterBase class to actually compute equality constraints
499 : OptimizationReporterBase * obj_func = solver->getObjFunction();
500 627 : obj_func->computeEqualityConstraints(eq_con);
501 627 : PetscFunctionReturn(PETSC_SUCCESS);
502 627 : }
503 :
504 : PetscErrorCode
505 627 : OptimizeSolve::equalityGradientFunctionWrapper(
506 : Tao /*tao*/, Vec /*x*/, Mat gradient_e, Mat /*gradient_epre*/, void * ctx)
507 : {
508 : PetscFunctionBegin;
509 : // grab the solver
510 : auto * solver = static_cast<OptimizeSolve *>(ctx);
511 627 : libMesh::PetscMatrix<Number> grad_eq(gradient_e, solver->_my_comm);
512 : // use the OptimizationReporterBase class to actually compute equality
513 : // constraints gradient
514 : OptimizationReporterBase * obj_func = solver->getObjFunction();
515 627 : obj_func->computeEqualityGradient(grad_eq);
516 627 : PetscFunctionReturn(PETSC_SUCCESS);
517 627 : }
518 :
519 : PetscErrorCode
520 504 : OptimizeSolve::inequalityFunctionWrapper(Tao /*tao*/, Vec /*x*/, Vec ci, void * ctx)
521 : {
522 : PetscFunctionBegin;
523 : // grab the solver
524 : auto * solver = static_cast<OptimizeSolve *>(ctx);
525 504 : libMesh::PetscVector<Number> ineq_con(ci, solver->_my_comm);
526 : // use the OptimizationReporterBase class to actually compute equality constraints
527 : OptimizationReporterBase * obj_func = solver->getObjFunction();
528 504 : obj_func->computeInequalityConstraints(ineq_con);
529 504 : PetscFunctionReturn(PETSC_SUCCESS);
530 504 : }
531 :
532 : PetscErrorCode
533 504 : OptimizeSolve::inequalityGradientFunctionWrapper(
534 : Tao /*tao*/, Vec /*x*/, Mat gradient_i, Mat /*gradient_ipre*/, void * ctx)
535 : {
536 : PetscFunctionBegin;
537 : // grab the solver
538 : auto * solver = static_cast<OptimizeSolve *>(ctx);
539 504 : libMesh::PetscMatrix<Number> grad_ineq(gradient_i, solver->_my_comm);
540 : // use the OptimizationReporterBase class to actually compute equality
541 : // constraints gradient
542 : OptimizationReporterBase * obj_func = solver->getObjFunction();
543 504 : obj_func->computeInequalityGradient(grad_ineq);
544 504 : PetscFunctionReturn(PETSC_SUCCESS);
545 504 : }
546 :
547 : PetscErrorCode
548 28 : OptimizeSolve::taoALCreate()
549 : {
550 : PetscFunctionBegin;
551 28 : if (_obj_function->getNumEqCons())
552 : {
553 : // Create equality vector
554 19 : LibmeshPetscCallQ(VecCreate(_my_comm.get(), &_ce));
555 19 : LibmeshPetscCallQ(
556 : VecSetSizes(_ce, _obj_function->getNumEqCons(), _obj_function->getNumEqCons()));
557 19 : LibmeshPetscCallQ(VecSetFromOptions(_ce));
558 19 : LibmeshPetscCallQ(VecSetUp(_ce));
559 :
560 : // Set equality jacobian matrix
561 19 : LibmeshPetscCallQ(MatCreate(_my_comm.get(), &_gradient_e));
562 19 : LibmeshPetscCallQ(MatSetSizes(
563 : _gradient_e, _obj_function->getNumEqCons(), _ndof, _obj_function->getNumEqCons(), _ndof));
564 19 : LibmeshPetscCallQ(MatSetFromOptions(_gradient_e));
565 19 : LibmeshPetscCallQ(MatSetUp(_gradient_e));
566 :
567 : // Set the Equality Constraints
568 19 : LibmeshPetscCallQ(TaoSetEqualityConstraintsRoutine(_tao, _ce, equalityFunctionWrapper, this));
569 :
570 : // Set the Equality Constraints Jacobian
571 19 : LibmeshPetscCallQ(TaoSetJacobianEqualityRoutine(
572 : _tao, _gradient_e, _gradient_e, equalityGradientFunctionWrapper, this));
573 : }
574 :
575 28 : if (_obj_function->getNumInEqCons())
576 : {
577 : // Create inequality vector
578 9 : LibmeshPetscCallQ(VecCreate(_my_comm.get(), &_ci));
579 9 : LibmeshPetscCallQ(
580 : VecSetSizes(_ci, _obj_function->getNumInEqCons(), _obj_function->getNumInEqCons()));
581 9 : LibmeshPetscCallQ(VecSetFromOptions(_ci));
582 9 : LibmeshPetscCallQ(VecSetUp(_ci));
583 :
584 : // Set inequality jacobian matrix
585 9 : LibmeshPetscCallQ(MatCreate(_my_comm.get(), &_gradient_i));
586 9 : LibmeshPetscCallQ(MatSetSizes(_gradient_i,
587 : _obj_function->getNumInEqCons(),
588 : _ndof,
589 : _obj_function->getNumInEqCons(),
590 : _ndof));
591 9 : LibmeshPetscCallQ(MatSetFromOptions(_gradient_i));
592 9 : LibmeshPetscCallQ(MatSetUp(_gradient_i));
593 :
594 : // Set the Inequality constraints
595 9 : LibmeshPetscCallQ(
596 : TaoSetInequalityConstraintsRoutine(_tao, _ci, inequalityFunctionWrapper, this));
597 :
598 : // Set the Inequality constraints Jacobian
599 9 : LibmeshPetscCallQ(TaoSetJacobianInequalityRoutine(
600 : _tao, _gradient_i, _gradient_i, inequalityGradientFunctionWrapper, this));
601 : }
602 28 : PetscFunctionReturn(PETSC_SUCCESS);
603 : }
604 :
605 : PetscErrorCode
606 28 : OptimizeSolve::taoALDestroy()
607 : {
608 : PetscFunctionBegin;
609 28 : if (_obj_function->getNumEqCons())
610 : {
611 19 : LibmeshPetscCallQ(VecDestroy(&_ce));
612 19 : LibmeshPetscCallQ(MatDestroy(&_gradient_e));
613 : }
614 28 : if (_obj_function->getNumInEqCons())
615 : {
616 9 : LibmeshPetscCallQ(VecDestroy(&_ci));
617 9 : LibmeshPetscCallQ(MatDestroy(&_gradient_i));
618 : }
619 :
620 28 : PetscFunctionReturn(PETSC_SUCCESS);
621 : }
|
