v0.16.3
Loading...
Searching...
No Matches
Functions | Variables
tao_ctx_constraints.cpp File Reference

Test TAO constraint pipelines and application evaluation hooks. More...

#include <MoFEM.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "Test TaoCtx constraint callbacks.\n\n"
 

Detailed Description

Test TAO constraint pipelines and application evaluation hooks.

Definition in file tao_ctx_constraints.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 11 of file tao_ctx_constraints.cpp.

11 {
12 MoFEM::Core::Initialize(&argc, &argv, nullptr, help);
13
14 try {
15 moab::Core moab_core;
16 MoFEM::Core mofem_core(moab_core);
17 MoFEM::Interface &m_field = mofem_core;
18
19 CHKERR DMRegister_MoFEM("DMMOFEM");
20
21 auto simple = m_field.getInterface<Simple>();
22 CHKERR simple->getOptions();
23 CHKERR simple->loadFile();
24 CHKERR simple->addDomainField("X", H1, AINSWORTH_LEGENDRE_BASE, 1);
25 CHKERR simple->setFieldOrder("X", 1);
26 CHKERR simple->setUp();
27
28 auto dm = simple->getDM();
29 auto x = createDMVector(dm);
30 auto c = vectorDuplicate(x);
31 auto J = createDMMatrix(dm);
32 auto diagonal = vectorDuplicate(x);
33
34 Tao tao;
35 CHKERR TaoCreate(PETSC_COMM_WORLD, &tao);
36 CHKERR TaoSetType(tao, TAOALMM);
37 CHKERR TaoSetSolution(tao, x);
38
39 TaoCtx tao_ctx(m_field, simple->getProblemName());
40
41 auto check_vector_pipeline = [&](auto context, const double value,
42 auto &pre_process, auto &post_process) {
43 auto pre = boost::make_shared<BasicMethod>();
44 auto post = boost::make_shared<BasicMethod>();
45 auto *pre_ptr = pre.get();
46 auto *post_ptr = post.get();
47
48 pre->preProcessHook = [pre_ptr, context, value, c]() {
49 MoFEMFunctionBegin;
50 if (pre_ptr->tao_ctx != context || pre_ptr->tao_c != c)
51 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
52 "Incorrect TAO constraint context");
53 CHKERR VecSet(pre_ptr->tao_c, value);
54 MoFEMFunctionReturn(0);
55 };
56 post->postProcessHook = [post_ptr, context, c]() {
57 MoFEMFunctionBegin;
58 if (post_ptr->tao_ctx != context || post_ptr->tao_c != c)
59 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
60 "Incorrect TAO constraint context");
61 MoFEMFunctionReturn(0);
62 };
63
64 pre_process.push_back(pre);
65 post_process.push_back(post);
66 };
67
68 auto check_jacobian_pipeline = [&](auto context, const double value,
69 auto &pre_process, auto &post_process) {
70 auto pre = boost::make_shared<BasicMethod>();
71 auto post = boost::make_shared<BasicMethod>();
72 auto *pre_ptr = pre.get();
73 auto *post_ptr = post.get();
74
75 pre->preProcessHook = [pre_ptr, context, value, J]() {
76 MoFEMFunctionBegin;
77 if (pre_ptr->tao_ctx != context || pre_ptr->tao_J != J ||
78 pre_ptr->tao_Jpre != J)
79 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
80 "Incorrect TAO constraint Jacobian context");
81 CHKERR MatShift(pre_ptr->tao_Jpre, value);
82 MoFEMFunctionReturn(0);
83 };
84 post->postProcessHook = [post_ptr, context, J]() {
85 MoFEMFunctionBegin;
86 if (post_ptr->tao_ctx != context || post_ptr->tao_J != J ||
87 post_ptr->tao_Jpre != J)
88 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
89 "Incorrect TAO constraint Jacobian context");
90 MoFEMFunctionReturn(0);
91 };
92
93 pre_process.push_back(pre);
94 post_process.push_back(post);
95 };
96
97 check_vector_pipeline(TaoMethod::CTX_TAO_EQUALITY_CONSTRAINT, 2.0,
98 tao_ctx.getPreProcEqualityConstraint(),
99 tao_ctx.getPostProcEqualityConstraint());
100 check_vector_pipeline(TaoMethod::CTX_TAO_INEQUALITY_CONSTRAINT, 3.0,
101 tao_ctx.getPreProcInequalityConstraint(),
102 tao_ctx.getPostProcInequalityConstraint());
103 check_jacobian_pipeline(TaoMethod::CTX_TAO_EQUALITY_JACOBIAN, 4.0,
104 tao_ctx.getPreProcEqualityJacobian(),
105 tao_ctx.getPostProcEqualityJacobian());
106 check_jacobian_pipeline(TaoMethod::CTX_TAO_INEQUALITY_JACOBIAN, 5.0,
107 tao_ctx.getPreProcInequalityJacobian(),
108 tao_ctx.getPostProcInequalityJacobian());
109
110 CHKERR TaoSetEqualityConstraintsRoutine(
111 tao, c, MoFEM::TaoSetEqualityConstraints, &tao_ctx);
112 CHKERR TaoSetInequalityConstraintsRoutine(
113 tao, c, MoFEM::TaoSetInequalityConstraints, &tao_ctx);
114 CHKERR TaoSetJacobianEqualityRoutine(
115 tao, J, J, MoFEM::TaoSetJacobianEquality, &tao_ctx);
116 CHKERR TaoSetJacobianInequalityRoutine(
117 tao, J, J, MoFEM::TaoSetJacobianInequality, &tao_ctx);
118
119 auto check_vector = [&](const double expected) {
120 MoFEMFunctionBegin;
121 CHKERR VecShift(c, -expected);
122 PetscReal norm;
123 CHKERR VecNorm(c, NORM_INFINITY, &norm);
124 if (norm > 10 * PETSC_MACHINE_EPSILON)
125 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
126 "Incorrect TAO constraint vector");
127 MoFEMFunctionReturn(0);
128 };
129
130 auto check_jacobian = [&](const double expected) {
131 MoFEMFunctionBegin;
132 CHKERR MatGetDiagonal(J, diagonal);
133 CHKERR VecShift(diagonal, -expected);
134 PetscReal norm;
135 CHKERR VecNorm(diagonal, NORM_INFINITY, &norm);
136 if (norm > 10 * PETSC_MACHINE_EPSILON)
137 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
138 "Incorrect TAO constraint Jacobian");
139 MoFEMFunctionReturn(0);
140 };
141
142 CHKERR TaoComputeEqualityConstraints(tao, x, c);
143 CHKERR check_vector(2.0);
144 CHKERR TaoComputeInequalityConstraints(tao, x, c);
145 CHKERR check_vector(3.0);
146 CHKERR TaoComputeJacobianEquality(tao, x, J, J);
147 CHKERR check_jacobian(4.0);
148 CHKERR TaoComputeJacobianInequality(tao, x, J, J);
149 CHKERR check_jacobian(5.0);
150
151 TaoCtx copied_ctx(m_field, simple->getProblemName());
152 CHKERR copied_ctx.copyLoops(tao_ctx);
153 if (copied_ctx.getPreProcEqualityConstraint().empty() ||
154 copied_ctx.getPreProcInequalityConstraint().empty() ||
155 copied_ctx.getPreProcEqualityJacobian().empty() ||
156 copied_ctx.getPreProcInequalityJacobian().empty())
157 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
158 "TAO constraint pipelines were not copied");
159 CHKERR copied_ctx.clearLoops();
160 if (!copied_ctx.getPreProcEqualityConstraint().empty() ||
161 !copied_ctx.getPreProcInequalityConstraint().empty() ||
162 !copied_ctx.getPreProcEqualityJacobian().empty() ||
163 !copied_ctx.getPreProcInequalityJacobian().empty())
164 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
165 "TAO constraint pipelines were not cleared");
166
167 // A reduced objective owns its mesh state and can reuse its gradient.
168 // Neither wrapper may scatter a trial or clear that cached result.
169 CHKERR VecZeroEntries(x);
170 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
171 simple->getProblemName(), COL, x, INSERT_VALUES, SCATTER_REVERSE);
172 CHKERR VecSet(x, 2.0);
173 auto gradient = vectorDuplicate(x);
174 auto expected_gradient = vectorDuplicate(x);
175 CHKERR VecSet(gradient, 7.0);
176 CHKERR VecCopy(gradient, expected_gradient);
177 PetscInt objective_calls = 0, gradient_calls = 0;
178 tao_ctx.getObjectiveHook() = [&](Tao, Vec trial, PetscReal *objective) {
179 MoFEMFunctionBegin;
180 if (trial != x)
181 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
182 "Incorrect TAO objective-hook controls");
183 ++objective_calls;
184 *objective = 3.0;
185 MoFEMFunctionReturn(0);
186 };
187 tao_ctx.getGradientHook() = [&](Tao, Vec trial, Vec g) {
188 MoFEMFunctionBegin;
189 PetscBool equal;
190 CHKERR VecEqual(g, expected_gradient, &equal);
191 if (trial != x || !equal)
192 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
193 "TAO gradient hook did not receive the cached gradient");
194 ++gradient_calls;
195 MoFEMFunctionReturn(0);
196 };
197 CHKERR ::TaoSetObjective(tao, MoFEM::TaoSetObjective, &tao_ctx);
198 CHKERR ::TaoSetGradient(tao, gradient, MoFEM::TaoSetGradient, &tao_ctx);
199 CHKERR ::TaoSetObjectiveAndGradient(
200 tao, gradient, MoFEM::TaoSetObjectiveAndGradient, &tao_ctx);
201 PetscReal objective;
202 CHKERR TaoComputeObjective(tao, x, &objective);
203 if (objective != 3.0 || objective_calls != 1 || gradient_calls != 0)
204 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
205 "Objective-only wrapper did not dispatch exclusively to its hook");
206 CHKERR TaoComputeGradient(tao, x, gradient);
207 if (objective_calls != 1 || gradient_calls != 1)
208 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
209 "Gradient-only wrapper did not dispatch exclusively to its hook");
210 CHKERR TaoComputeObjectiveAndGradient(tao, x, &objective, gradient);
211 PetscBool equal;
212 CHKERR VecEqual(gradient, expected_gradient, &equal);
213 if (objective != 3.0 || objective_calls != 2 || gradient_calls != 2 ||
214 !equal)
215 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
216 "Combined wrapper did not fall back to the separate hooks");
217 auto mesh_controls = vectorDuplicate(x);
218 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
219 simple->getProblemName(), COL, mesh_controls, INSERT_VALUES,
220 SCATTER_FORWARD);
221 PetscReal mesh_norm;
222 CHKERR VecNorm(mesh_controls, NORM_INFINITY, &mesh_norm);
223 if (mesh_norm != 0)
224 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
225 "TAO evaluation hooks implicitly scattered trial controls");
226
227 CHKERR TaoDestroy(&tao);
228 }
229 CATCH_ERRORS;
230
231 CHKERR MoFEM::Core::Finalize();
232 return 0;
233}
diagonal
SmartPetscObj< Vec > diagonal
Definition TaoDiagonalNNLS.cpp:21
gradient
SmartPetscObj< Vec > gradient
Definition TaoDiagonalNNLS.cpp:20
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
COL
@ COL
Definition definitions.h:136
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
H1
@ H1
continuous field
Definition definitions.h:85
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
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::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEM::createDMVector
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1250
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1207
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:40
J
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition MatrixFunction.hpp:64
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:237
MoFEM::TaoSetJacobianInequality
PetscErrorCode TaoSetJacobianInequality(Tao tao, Vec x, Mat J, Mat Jpre, void *ctx)
Assemble the inequality constraint Jacobian.
Definition TaoCtx.cpp:502
MoFEM::TaoSetEqualityConstraints
PetscErrorCode TaoSetEqualityConstraints(Tao tao, Vec x, Vec c, void *ctx)
Evaluate equality constraints for a TAO optimization context.
Definition TaoCtx.cpp:468
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
MoFEM::TaoSetInequalityConstraints
PetscErrorCode TaoSetInequalityConstraints(Tao tao, Vec x, Vec c, void *ctx)
Evaluate inequality constraints for a TAO optimization context.
Definition TaoCtx.cpp:479
MoFEM::TaoSetJacobianEquality
PetscErrorCode TaoSetJacobianEquality(Tao tao, Vec x, Mat J, Mat Jpre, void *ctx)
Assemble the equality constraint Jacobian.
Definition TaoCtx.cpp:490
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:386
MoFEM::TaoSetGradient
PetscErrorCode TaoSetGradient(Tao tao, Vec x, Vec f, void *ctx)
Sets the gradient vector for a TAO optimization context.
Definition TaoCtx.cpp:296
g
constexpr double g
Definition shallow_wave.cpp:63
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:83
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:68
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::Simple
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEM::TaoCtx
Interface for TAO solvers.
Definition TaoCtx.hpp:14
MoFEM::TaoMethod::CTX_TAO_EQUALITY_CONSTRAINT
@ CTX_TAO_EQUALITY_CONSTRAINT
Evaluating equality constraints.
Definition LoopMethods.hpp:347
MoFEM::TaoMethod::CTX_TAO_EQUALITY_JACOBIAN
@ CTX_TAO_EQUALITY_JACOBIAN
Assembling equality Jacobian.
Definition LoopMethods.hpp:349
MoFEM::TaoMethod::CTX_TAO_INEQUALITY_JACOBIAN
@ CTX_TAO_INEQUALITY_JACOBIAN
Assembling inequality Jacobian.
Definition LoopMethods.hpp:350
MoFEM::TaoMethod::CTX_TAO_INEQUALITY_CONSTRAINT
@ CTX_TAO_INEQUALITY_CONSTRAINT
Evaluating inequality constraints.
Definition LoopMethods.hpp:348
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:103
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition VecManager.hpp:23
help
static char help[]
Definition tao_ctx_constraints.cpp:9

Variable Documentation

◆ help

char help[] = "Test TaoCtx constraint callbacks.\n\n"
static

Definition at line 9 of file tao_ctx_constraints.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow