1 #include <petsctao.h>
2
3 typedef struct {
4 PetscInt n; /* Number of variables */
5 PetscInt m; /* Number of constraints */
6 PetscInt mx; /* grid points in each direction */
7 PetscInt nt; /* Number of time steps */
8 PetscInt ndata; /* Number of data points per sample */
9 IS s_is;
10 IS d_is;
11 VecScatter state_scatter;
12 VecScatter design_scatter;
13 VecScatter *uxi_scatter, *uyi_scatter, *ux_scatter, *uy_scatter, *ui_scatter;
14 VecScatter *yi_scatter;
15
16 Mat Js, Jd, JsBlockPrec, JsInv, JsBlock;
17 PetscBool jformed, c_formed;
18
19 PetscReal alpha; /* Regularization parameter */
20 PetscReal gamma;
21 PetscReal ht; /* Time step */
22 PetscReal T; /* Final time */
23 Mat Q, QT;
24 Mat L, LT;
25 Mat Div, Divwork, Divxy[2];
26 Mat Grad, Gradxy[2];
27 Mat M;
28 Mat *C, *Cwork;
29 /* Mat Hs,Hd,Hsd; */
30 Vec q;
31 Vec ur; /* reference */
32
33 Vec d;
34 Vec dwork;
35
36 Vec y; /* state variables */
37 Vec ywork;
38 Vec ytrue;
39 Vec *yi, *yiwork, *ziwork;
40 Vec *uxi, *uyi, *uxiwork, *uyiwork, *ui, *uiwork;
41
42 Vec u; /* design variables */
43 Vec uwork, vwork;
44 Vec utrue;
45
46 Vec js_diag;
47
48 Vec c; /* constraint vector */
49 Vec cwork;
50
51 Vec lwork;
52
53 KSP solver;
54 PC prec;
55 PetscInt block_index;
56
57 PetscInt ksp_its;
58 PetscInt ksp_its_initial;
59 } AppCtx;
60
61 PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
62 PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
63 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
64 PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void *);
65 PetscErrorCode FormJacobianDesign(Tao, Vec, Mat, void *);
66 PetscErrorCode FormConstraints(Tao, Vec, Vec, void *);
67 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
68 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
69 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
70 PetscErrorCode HyperbolicInitialize(AppCtx *user);
71 PetscErrorCode HyperbolicDestroy(AppCtx *user);
72 PetscErrorCode HyperbolicMonitor(Tao, void *);
73
74 PetscErrorCode StateMatMult(Mat, Vec, Vec);
75 PetscErrorCode StateMatBlockMult(Mat, Vec, Vec);
76 PetscErrorCode StateMatBlockMultTranspose(Mat, Vec, Vec);
77 PetscErrorCode StateMatMultTranspose(Mat, Vec, Vec);
78 PetscErrorCode StateMatGetDiagonal(Mat, Vec);
79 PetscErrorCode StateMatDuplicate(Mat, MatDuplicateOption, Mat *);
80 PetscErrorCode StateMatInvMult(Mat, Vec, Vec);
81 PetscErrorCode StateMatInvTransposeMult(Mat, Vec, Vec);
82 PetscErrorCode StateMatBlockPrecMult(PC, Vec, Vec);
83
84 PetscErrorCode DesignMatMult(Mat, Vec, Vec);
85 PetscErrorCode DesignMatMultTranspose(Mat, Vec, Vec);
86
87 PetscErrorCode Scatter_yi(Vec, Vec *, VecScatter *, PetscInt); /* y to y1,y2,...,y_nt */
88 PetscErrorCode Gather_yi(Vec, Vec *, VecScatter *, PetscInt);
89 PetscErrorCode Scatter_uxi_uyi(Vec, Vec *, VecScatter *, Vec *, VecScatter *, PetscInt); /* u to ux_1,uy_1,ux_2,uy_2,...,u */
90 PetscErrorCode Gather_uxi_uyi(Vec, Vec *, VecScatter *, Vec *, VecScatter *, PetscInt);
91
92 static char help[] = "";
93
main(int argc,char ** argv)94 int main(int argc, char **argv)
95 {
96 Vec x, x0;
97 Tao tao;
98 AppCtx user;
99 IS is_allstate, is_alldesign;
100 PetscInt lo, hi, hi2, lo2, ksp_old;
101 PetscInt ntests = 1;
102 PetscInt i;
103 PetscLogStage stages[1];
104
105 PetscFunctionBeginUser;
106 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
107 user.mx = 32;
108 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "hyperbolic example", NULL);
109 PetscCall(PetscOptionsInt("-mx", "Number of grid points in each direction", "", user.mx, &user.mx, NULL));
110 user.nt = 16;
111 PetscCall(PetscOptionsInt("-nt", "Number of time steps", "", user.nt, &user.nt, NULL));
112 user.ndata = 64;
113 PetscCall(PetscOptionsInt("-ndata", "Numbers of data points per sample", "", user.ndata, &user.ndata, NULL));
114 user.alpha = 10.0;
115 PetscCall(PetscOptionsReal("-alpha", "Regularization parameter", "", user.alpha, &user.alpha, NULL));
116 user.T = 1.0 / 32.0;
117 PetscCall(PetscOptionsReal("-Tfinal", "Final time", "", user.T, &user.T, NULL));
118 PetscCall(PetscOptionsInt("-ntests", "Number of times to repeat TaoSolve", "", ntests, &ntests, NULL));
119 PetscOptionsEnd();
120
121 user.m = user.mx * user.mx * user.nt; /* number of constraints */
122 user.n = user.mx * user.mx * 3 * user.nt; /* number of variables */
123 user.ht = user.T / user.nt; /* Time step */
124 user.gamma = user.T * user.ht / (user.mx * user.mx);
125
126 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.u));
127 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.y));
128 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.c));
129 PetscCall(VecSetSizes(user.u, PETSC_DECIDE, user.n - user.m));
130 PetscCall(VecSetSizes(user.y, PETSC_DECIDE, user.m));
131 PetscCall(VecSetSizes(user.c, PETSC_DECIDE, user.m));
132 PetscCall(VecSetFromOptions(user.u));
133 PetscCall(VecSetFromOptions(user.y));
134 PetscCall(VecSetFromOptions(user.c));
135
136 /* Create scatters for reduced spaces.
137 If the state vector y and design vector u are partitioned as
138 [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
139 then the solution vector x is organized as
140 [y_1; u_1; y_2; u_2; ...; y_np; u_np].
141 The index sets user.s_is and user.d_is correspond to the indices of the
142 state and design variables owned by the current processor.
143 */
144 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
145
146 PetscCall(VecGetOwnershipRange(user.y, &lo, &hi));
147 PetscCall(VecGetOwnershipRange(user.u, &lo2, &hi2));
148
149 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_allstate));
150 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + lo2, 1, &user.s_is));
151
152 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, lo2, 1, &is_alldesign));
153 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, hi + lo2, 1, &user.d_is));
154
155 PetscCall(VecSetSizes(x, hi - lo + hi2 - lo2, user.n));
156 PetscCall(VecSetFromOptions(x));
157
158 PetscCall(VecScatterCreate(x, user.s_is, user.y, is_allstate, &user.state_scatter));
159 PetscCall(VecScatterCreate(x, user.d_is, user.u, is_alldesign, &user.design_scatter));
160 PetscCall(ISDestroy(&is_alldesign));
161 PetscCall(ISDestroy(&is_allstate));
162
163 /* Create TAO solver and set desired solution method */
164 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
165 PetscCall(TaoSetType(tao, TAOLCL));
166
167 /* Set up initial vectors and matrices */
168 PetscCall(HyperbolicInitialize(&user));
169
170 PetscCall(Gather(x, user.y, user.state_scatter, user.u, user.design_scatter));
171 PetscCall(VecDuplicate(x, &x0));
172 PetscCall(VecCopy(x, x0));
173
174 /* Set solution vector with an initial guess */
175 PetscCall(TaoSetSolution(tao, x));
176 PetscCall(TaoSetObjective(tao, FormFunction, &user));
177 PetscCall(TaoSetGradient(tao, NULL, FormGradient, &user));
178 PetscCall(TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user));
179 PetscCall(TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, &user));
180 PetscCall(TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user));
181 PetscCall(TaoSetFromOptions(tao));
182 PetscCall(TaoSetStateDesignIS(tao, user.s_is, user.d_is));
183
184 /* SOLVE THE APPLICATION */
185 PetscCall(PetscLogStageRegister("Trials", &stages[0]));
186 PetscCall(PetscLogStagePush(stages[0]));
187 user.ksp_its_initial = user.ksp_its;
188 ksp_old = user.ksp_its;
189 for (i = 0; i < ntests; i++) {
190 PetscCall(TaoSolve(tao));
191 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP Iterations = %" PetscInt_FMT "\n", user.ksp_its - ksp_old));
192 PetscCall(VecCopy(x0, x));
193 PetscCall(TaoSetSolution(tao, x));
194 }
195 PetscCall(PetscLogStagePop());
196 PetscCall(PetscBarrier((PetscObject)x));
197 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations within initialization: "));
198 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its_initial));
199 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Total KSP iterations over %" PetscInt_FMT " trial(s): ", ntests));
200 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its));
201 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations per trial: "));
202 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", (user.ksp_its - user.ksp_its_initial) / ntests));
203
204 PetscCall(TaoDestroy(&tao));
205 PetscCall(VecDestroy(&x));
206 PetscCall(VecDestroy(&x0));
207 PetscCall(HyperbolicDestroy(&user));
208 PetscCall(PetscFinalize());
209 return 0;
210 }
211 /* ------------------------------------------------------------------- */
212 /*
213 dwork = Qy - d
214 lwork = L*(u-ur).^2
215 f = 1/2 * (dwork.dork + alpha*y.lwork)
216 */
FormFunction(Tao tao,Vec X,PetscReal * f,void * ptr)217 PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr)
218 {
219 PetscReal d1 = 0, d2 = 0;
220 AppCtx *user = (AppCtx *)ptr;
221
222 PetscFunctionBegin;
223 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
224 PetscCall(MatMult(user->Q, user->y, user->dwork));
225 PetscCall(VecAXPY(user->dwork, -1.0, user->d));
226 PetscCall(VecDot(user->dwork, user->dwork, &d1));
227
228 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u));
229 PetscCall(VecPointwiseMult(user->uwork, user->uwork, user->uwork));
230 PetscCall(MatMult(user->L, user->uwork, user->lwork));
231 PetscCall(VecDot(user->y, user->lwork, &d2));
232 *f = 0.5 * (d1 + user->alpha * d2);
233 PetscFunctionReturn(PETSC_SUCCESS);
234 }
235
236 /* ------------------------------------------------------------------- */
237 /*
238 state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2
239 design: g_d = alpha*(L'y).*(u-ur)
240 */
FormGradient(Tao tao,Vec X,Vec G,void * ptr)241 PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr)
242 {
243 AppCtx *user = (AppCtx *)ptr;
244
245 PetscFunctionBegin;
246 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
247 PetscCall(MatMult(user->Q, user->y, user->dwork));
248 PetscCall(VecAXPY(user->dwork, -1.0, user->d));
249
250 PetscCall(MatMult(user->QT, user->dwork, user->ywork));
251
252 PetscCall(MatMult(user->LT, user->y, user->uwork));
253 PetscCall(VecWAXPY(user->vwork, -1.0, user->ur, user->u));
254 PetscCall(VecPointwiseMult(user->uwork, user->vwork, user->uwork));
255 PetscCall(VecScale(user->uwork, user->alpha));
256
257 PetscCall(VecPointwiseMult(user->vwork, user->vwork, user->vwork));
258 PetscCall(MatMult(user->L, user->vwork, user->lwork));
259 PetscCall(VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork));
260
261 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
262 PetscFunctionReturn(PETSC_SUCCESS);
263 }
264
FormFunctionGradient(Tao tao,Vec X,PetscReal * f,Vec G,void * ptr)265 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
266 {
267 PetscReal d1, d2;
268 AppCtx *user = (AppCtx *)ptr;
269
270 PetscFunctionBegin;
271 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
272 PetscCall(MatMult(user->Q, user->y, user->dwork));
273 PetscCall(VecAXPY(user->dwork, -1.0, user->d));
274
275 PetscCall(MatMult(user->QT, user->dwork, user->ywork));
276
277 PetscCall(VecDot(user->dwork, user->dwork, &d1));
278
279 PetscCall(MatMult(user->LT, user->y, user->uwork));
280 PetscCall(VecWAXPY(user->vwork, -1.0, user->ur, user->u));
281 PetscCall(VecPointwiseMult(user->uwork, user->vwork, user->uwork));
282 PetscCall(VecScale(user->uwork, user->alpha));
283
284 PetscCall(VecPointwiseMult(user->vwork, user->vwork, user->vwork));
285 PetscCall(MatMult(user->L, user->vwork, user->lwork));
286 PetscCall(VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork));
287
288 PetscCall(VecDot(user->y, user->lwork, &d2));
289
290 *f = 0.5 * (d1 + user->alpha * d2);
291 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
292 PetscFunctionReturn(PETSC_SUCCESS);
293 }
294
295 /* ------------------------------------------------------------------- */
296 /* A
297 MatShell object
298 */
FormJacobianState(Tao tao,Vec X,Mat J,Mat JPre,Mat JInv,void * ptr)299 PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
300 {
301 PetscInt i;
302 AppCtx *user = (AppCtx *)ptr;
303
304 PetscFunctionBegin;
305 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
306 PetscCall(Scatter_yi(user->u, user->ui, user->ui_scatter, user->nt));
307 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
308 for (i = 0; i < user->nt; i++) {
309 PetscCall(MatCopy(user->Divxy[0], user->C[i], SUBSET_NONZERO_PATTERN));
310 PetscCall(MatCopy(user->Divxy[1], user->Cwork[i], SAME_NONZERO_PATTERN));
311
312 PetscCall(MatDiagonalScale(user->C[i], NULL, user->uxi[i]));
313 PetscCall(MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i]));
314 PetscCall(MatAXPY(user->C[i], 1.0, user->Cwork[i], SUBSET_NONZERO_PATTERN));
315 PetscCall(MatScale(user->C[i], user->ht));
316 PetscCall(MatShift(user->C[i], 1.0));
317 }
318 PetscFunctionReturn(PETSC_SUCCESS);
319 }
320
321 /* ------------------------------------------------------------------- */
322 /* B */
FormJacobianDesign(Tao tao,Vec X,Mat J,void * ptr)323 PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
324 {
325 AppCtx *user = (AppCtx *)ptr;
326
327 PetscFunctionBegin;
328 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
329 PetscFunctionReturn(PETSC_SUCCESS);
330 }
331
StateMatMult(Mat J_shell,Vec X,Vec Y)332 PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
333 {
334 PetscInt i;
335 AppCtx *user;
336
337 PetscFunctionBegin;
338 PetscCall(MatShellGetContext(J_shell, &user));
339 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt));
340 user->block_index = 0;
341 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0]));
342
343 for (i = 1; i < user->nt; i++) {
344 user->block_index = i;
345 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i]));
346 PetscCall(MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1]));
347 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1]));
348 }
349 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
350 PetscFunctionReturn(PETSC_SUCCESS);
351 }
352
StateMatMultTranspose(Mat J_shell,Vec X,Vec Y)353 PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
354 {
355 PetscInt i;
356 AppCtx *user;
357
358 PetscFunctionBegin;
359 PetscCall(MatShellGetContext(J_shell, &user));
360 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt));
361
362 for (i = 0; i < user->nt - 1; i++) {
363 user->block_index = i;
364 PetscCall(MatMultTranspose(user->JsBlock, user->yi[i], user->yiwork[i]));
365 PetscCall(MatMult(user->M, user->yi[i + 1], user->ziwork[i + 1]));
366 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i + 1]));
367 }
368
369 i = user->nt - 1;
370 user->block_index = i;
371 PetscCall(MatMultTranspose(user->JsBlock, user->yi[i], user->yiwork[i]));
372 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
373 PetscFunctionReturn(PETSC_SUCCESS);
374 }
375
StateMatBlockMult(Mat J_shell,Vec X,Vec Y)376 PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
377 {
378 PetscInt i;
379 AppCtx *user;
380
381 PetscFunctionBegin;
382 PetscCall(MatShellGetContext(J_shell, &user));
383 i = user->block_index;
384 PetscCall(VecPointwiseMult(user->uxiwork[i], X, user->uxi[i]));
385 PetscCall(VecPointwiseMult(user->uyiwork[i], X, user->uyi[i]));
386 PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
387 PetscCall(MatMult(user->Div, user->uiwork[i], Y));
388 PetscCall(VecAYPX(Y, user->ht, X));
389 PetscFunctionReturn(PETSC_SUCCESS);
390 }
391
StateMatBlockMultTranspose(Mat J_shell,Vec X,Vec Y)392 PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y)
393 {
394 PetscInt i;
395 AppCtx *user;
396
397 PetscFunctionBegin;
398 PetscCall(MatShellGetContext(J_shell, &user));
399 i = user->block_index;
400 PetscCall(MatMult(user->Grad, X, user->uiwork[i]));
401 PetscCall(Scatter(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
402 PetscCall(VecPointwiseMult(user->uxiwork[i], user->uxi[i], user->uxiwork[i]));
403 PetscCall(VecPointwiseMult(user->uyiwork[i], user->uyi[i], user->uyiwork[i]));
404 PetscCall(VecWAXPY(Y, 1.0, user->uxiwork[i], user->uyiwork[i]));
405 PetscCall(VecAYPX(Y, user->ht, X));
406 PetscFunctionReturn(PETSC_SUCCESS);
407 }
408
DesignMatMult(Mat J_shell,Vec X,Vec Y)409 PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
410 {
411 PetscInt i;
412 AppCtx *user;
413
414 PetscFunctionBegin;
415 PetscCall(MatShellGetContext(J_shell, &user));
416 PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt));
417 PetscCall(Scatter_uxi_uyi(X, user->uxiwork, user->uxi_scatter, user->uyiwork, user->uyi_scatter, user->nt));
418 for (i = 0; i < user->nt; i++) {
419 PetscCall(VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i]));
420 PetscCall(VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i]));
421 PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
422 PetscCall(MatMult(user->Div, user->uiwork[i], user->ziwork[i]));
423 PetscCall(VecScale(user->ziwork[i], user->ht));
424 }
425 PetscCall(Gather_yi(Y, user->ziwork, user->yi_scatter, user->nt));
426 PetscFunctionReturn(PETSC_SUCCESS);
427 }
428
DesignMatMultTranspose(Mat J_shell,Vec X,Vec Y)429 PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
430 {
431 PetscInt i;
432 AppCtx *user;
433
434 PetscFunctionBegin;
435 PetscCall(MatShellGetContext(J_shell, &user));
436 PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt));
437 PetscCall(Scatter_yi(X, user->yiwork, user->yi_scatter, user->nt));
438 for (i = 0; i < user->nt; i++) {
439 PetscCall(MatMult(user->Grad, user->yiwork[i], user->uiwork[i]));
440 PetscCall(Scatter(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
441 PetscCall(VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i]));
442 PetscCall(VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i]));
443 PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
444 PetscCall(VecScale(user->uiwork[i], user->ht));
445 }
446 PetscCall(Gather_yi(Y, user->uiwork, user->ui_scatter, user->nt));
447 PetscFunctionReturn(PETSC_SUCCESS);
448 }
449
StateMatBlockPrecMult(PC PC_shell,Vec X,Vec Y)450 PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
451 {
452 PetscInt i;
453 AppCtx *user;
454
455 PetscFunctionBegin;
456 PetscCall(PCShellGetContext(PC_shell, &user));
457 i = user->block_index;
458 PetscCheck(user->c_formed, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed");
459 PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y));
460 PetscFunctionReturn(PETSC_SUCCESS);
461 }
462
StateMatBlockPrecMultTranspose(PC PC_shell,Vec X,Vec Y)463 PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
464 {
465 PetscInt i;
466 AppCtx *user;
467
468 PetscFunctionBegin;
469 PetscCall(PCShellGetContext(PC_shell, &user));
470
471 i = user->block_index;
472 PetscCheck(user->c_formed, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed");
473 PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y));
474 PetscFunctionReturn(PETSC_SUCCESS);
475 }
476
StateMatInvMult(Mat J_shell,Vec X,Vec Y)477 PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
478 {
479 AppCtx *user;
480 PetscInt its, i;
481
482 PetscFunctionBegin;
483 PetscCall(MatShellGetContext(J_shell, &user));
484
485 if (Y == user->ytrue) {
486 /* First solve is done using true solution to set up problem */
487 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, PETSC_CURRENT, PETSC_CURRENT));
488 } else {
489 PetscCall(KSPSetTolerances(user->solver, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
490 }
491 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt));
492 PetscCall(Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt));
493 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
494
495 user->block_index = 0;
496 PetscCall(KSPSolve(user->solver, user->yi[0], user->yiwork[0]));
497
498 PetscCall(KSPGetIterationNumber(user->solver, &its));
499 user->ksp_its = user->ksp_its + its;
500 for (i = 1; i < user->nt; i++) {
501 PetscCall(MatMult(user->M, user->yiwork[i - 1], user->ziwork[i - 1]));
502 PetscCall(VecAXPY(user->yi[i], 1.0, user->ziwork[i - 1]));
503 user->block_index = i;
504 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i]));
505
506 PetscCall(KSPGetIterationNumber(user->solver, &its));
507 user->ksp_its = user->ksp_its + its;
508 }
509 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
510 PetscFunctionReturn(PETSC_SUCCESS);
511 }
512
StateMatInvTransposeMult(Mat J_shell,Vec X,Vec Y)513 PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
514 {
515 AppCtx *user;
516 PetscInt its, i;
517
518 PetscFunctionBegin;
519 PetscCall(MatShellGetContext(J_shell, &user));
520
521 PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt));
522 PetscCall(Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt));
523 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
524
525 i = user->nt - 1;
526 user->block_index = i;
527 PetscCall(KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i]));
528
529 PetscCall(KSPGetIterationNumber(user->solver, &its));
530 user->ksp_its = user->ksp_its + its;
531
532 for (i = user->nt - 2; i >= 0; i--) {
533 PetscCall(MatMult(user->M, user->yiwork[i + 1], user->ziwork[i + 1]));
534 PetscCall(VecAXPY(user->yi[i], 1.0, user->ziwork[i + 1]));
535 user->block_index = i;
536 PetscCall(KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i]));
537
538 PetscCall(KSPGetIterationNumber(user->solver, &its));
539 user->ksp_its = user->ksp_its + its;
540 }
541 PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
542 PetscFunctionReturn(PETSC_SUCCESS);
543 }
544
StateMatDuplicate(Mat J_shell,MatDuplicateOption opt,Mat * new_shell)545 PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
546 {
547 AppCtx *user;
548
549 PetscFunctionBegin;
550 PetscCall(MatShellGetContext(J_shell, &user));
551
552 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, new_shell));
553 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT, (PetscErrorCodeFn *)StateMatMult));
554 PetscCall(MatShellSetOperation(*new_shell, MATOP_DUPLICATE, (PetscErrorCodeFn *)StateMatDuplicate));
555 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatMultTranspose));
556 PetscCall(MatShellSetOperation(*new_shell, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)StateMatGetDiagonal));
557 PetscFunctionReturn(PETSC_SUCCESS);
558 }
559
StateMatGetDiagonal(Mat J_shell,Vec X)560 PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
561 {
562 AppCtx *user;
563
564 PetscFunctionBegin;
565 PetscCall(MatShellGetContext(J_shell, &user));
566 PetscCall(VecCopy(user->js_diag, X));
567 PetscFunctionReturn(PETSC_SUCCESS);
568 }
569
FormConstraints(Tao tao,Vec X,Vec C,void * ptr)570 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
571 {
572 /* con = Ay - q, A = [C(u1) 0 0 ... 0;
573 -M C(u2) 0 ... 0;
574 0 -M C(u3) ... 0;
575 ... ;
576 0 ... -M C(u_nt)]
577 C(u) = eye + ht*Div*[diag(u1); diag(u2)] */
578 PetscInt i;
579 AppCtx *user = (AppCtx *)ptr;
580
581 PetscFunctionBegin;
582 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
583 PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt));
584 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
585
586 user->block_index = 0;
587 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0]));
588
589 for (i = 1; i < user->nt; i++) {
590 user->block_index = i;
591 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i]));
592 PetscCall(MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1]));
593 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1]));
594 }
595
596 PetscCall(Gather_yi(C, user->yiwork, user->yi_scatter, user->nt));
597 PetscCall(VecAXPY(C, -1.0, user->q));
598 PetscFunctionReturn(PETSC_SUCCESS);
599 }
600
Scatter(Vec x,Vec state,VecScatter s_scat,Vec design,VecScatter d_scat)601 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
602 {
603 PetscFunctionBegin;
604 PetscCall(VecScatterBegin(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD));
605 PetscCall(VecScatterEnd(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD));
606 PetscCall(VecScatterBegin(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD));
607 PetscCall(VecScatterEnd(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD));
608 PetscFunctionReturn(PETSC_SUCCESS);
609 }
610
Scatter_uxi_uyi(Vec u,Vec * uxi,VecScatter * scatx,Vec * uyi,VecScatter * scaty,PetscInt nt)611 PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
612 {
613 PetscInt i;
614
615 PetscFunctionBegin;
616 for (i = 0; i < nt; i++) {
617 PetscCall(VecScatterBegin(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD));
618 PetscCall(VecScatterEnd(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD));
619 PetscCall(VecScatterBegin(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD));
620 PetscCall(VecScatterEnd(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD));
621 }
622 PetscFunctionReturn(PETSC_SUCCESS);
623 }
624
Gather(Vec x,Vec state,VecScatter s_scat,Vec design,VecScatter d_scat)625 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
626 {
627 PetscFunctionBegin;
628 PetscCall(VecScatterBegin(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE));
629 PetscCall(VecScatterEnd(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE));
630 PetscCall(VecScatterBegin(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE));
631 PetscCall(VecScatterEnd(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE));
632 PetscFunctionReturn(PETSC_SUCCESS);
633 }
634
Gather_uxi_uyi(Vec u,Vec * uxi,VecScatter * scatx,Vec * uyi,VecScatter * scaty,PetscInt nt)635 PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
636 {
637 PetscInt i;
638
639 PetscFunctionBegin;
640 for (i = 0; i < nt; i++) {
641 PetscCall(VecScatterBegin(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE));
642 PetscCall(VecScatterEnd(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE));
643 PetscCall(VecScatterBegin(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE));
644 PetscCall(VecScatterEnd(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE));
645 }
646 PetscFunctionReturn(PETSC_SUCCESS);
647 }
648
Scatter_yi(Vec y,Vec * yi,VecScatter * scat,PetscInt nt)649 PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
650 {
651 PetscInt i;
652
653 PetscFunctionBegin;
654 for (i = 0; i < nt; i++) {
655 PetscCall(VecScatterBegin(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD));
656 PetscCall(VecScatterEnd(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD));
657 }
658 PetscFunctionReturn(PETSC_SUCCESS);
659 }
660
Gather_yi(Vec y,Vec * yi,VecScatter * scat,PetscInt nt)661 PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
662 {
663 PetscInt i;
664
665 PetscFunctionBegin;
666 for (i = 0; i < nt; i++) {
667 PetscCall(VecScatterBegin(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE));
668 PetscCall(VecScatterEnd(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE));
669 }
670 PetscFunctionReturn(PETSC_SUCCESS);
671 }
672
HyperbolicInitialize(AppCtx * user)673 PetscErrorCode HyperbolicInitialize(AppCtx *user)
674 {
675 PetscInt n, i, j, linear_index, istart, iend, iblock, lo, hi;
676 Vec XX, YY, XXwork, YYwork, yi, uxi, ui, bc;
677 PetscReal h, sum;
678 PetscScalar hinv, neg_hinv, quarter = 0.25, one = 1.0, half_hinv, neg_half_hinv;
679 PetscScalar vx, vy, zero = 0.0;
680 IS is_from_y, is_to_yi, is_from_u, is_to_uxi, is_to_uyi;
681
682 PetscFunctionBegin;
683 user->jformed = PETSC_FALSE;
684 user->c_formed = PETSC_FALSE;
685
686 user->ksp_its = 0;
687 user->ksp_its_initial = 0;
688
689 n = user->mx * user->mx;
690
691 h = 1.0 / user->mx;
692 hinv = user->mx;
693 neg_hinv = -hinv;
694 half_hinv = hinv / 2.0;
695 neg_half_hinv = neg_hinv / 2.0;
696
697 /* Generate Grad matrix */
698 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Grad));
699 PetscCall(MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, 2 * n, n));
700 PetscCall(MatSetFromOptions(user->Grad));
701 PetscCall(MatMPIAIJSetPreallocation(user->Grad, 3, NULL, 3, NULL));
702 PetscCall(MatSeqAIJSetPreallocation(user->Grad, 3, NULL));
703 PetscCall(MatGetOwnershipRange(user->Grad, &istart, &iend));
704
705 for (i = istart; i < iend; i++) {
706 if (i < n) {
707 iblock = i / user->mx;
708 j = iblock * user->mx + ((i + user->mx - 1) % user->mx);
709 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
710 j = iblock * user->mx + ((i + 1) % user->mx);
711 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
712 }
713 if (i >= n) {
714 j = (i - user->mx) % n;
715 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
716 j = (j + 2 * user->mx) % n;
717 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
718 }
719 }
720
721 PetscCall(MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY));
722 PetscCall(MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY));
723
724 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[0]));
725 PetscCall(MatSetSizes(user->Gradxy[0], PETSC_DECIDE, PETSC_DECIDE, n, n));
726 PetscCall(MatSetFromOptions(user->Gradxy[0]));
727 PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[0], 3, NULL, 3, NULL));
728 PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[0], 3, NULL));
729 PetscCall(MatGetOwnershipRange(user->Gradxy[0], &istart, &iend));
730
731 for (i = istart; i < iend; i++) {
732 iblock = i / user->mx;
733 j = iblock * user->mx + ((i + user->mx - 1) % user->mx);
734 PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
735 j = iblock * user->mx + ((i + 1) % user->mx);
736 PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
737 PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &i, &zero, INSERT_VALUES));
738 }
739 PetscCall(MatAssemblyBegin(user->Gradxy[0], MAT_FINAL_ASSEMBLY));
740 PetscCall(MatAssemblyEnd(user->Gradxy[0], MAT_FINAL_ASSEMBLY));
741
742 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[1]));
743 PetscCall(MatSetSizes(user->Gradxy[1], PETSC_DECIDE, PETSC_DECIDE, n, n));
744 PetscCall(MatSetFromOptions(user->Gradxy[1]));
745 PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[1], 3, NULL, 3, NULL));
746 PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[1], 3, NULL));
747 PetscCall(MatGetOwnershipRange(user->Gradxy[1], &istart, &iend));
748
749 for (i = istart; i < iend; i++) {
750 j = (i + n - user->mx) % n;
751 PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
752 j = (j + 2 * user->mx) % n;
753 PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
754 PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &i, &zero, INSERT_VALUES));
755 }
756 PetscCall(MatAssemblyBegin(user->Gradxy[1], MAT_FINAL_ASSEMBLY));
757 PetscCall(MatAssemblyEnd(user->Gradxy[1], MAT_FINAL_ASSEMBLY));
758
759 /* Generate Div matrix */
760 PetscCall(MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div));
761 PetscCall(MatTranspose(user->Gradxy[0], MAT_INITIAL_MATRIX, &user->Divxy[0]));
762 PetscCall(MatTranspose(user->Gradxy[1], MAT_INITIAL_MATRIX, &user->Divxy[1]));
763
764 /* Off-diagonal averaging matrix */
765 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->M));
766 PetscCall(MatSetSizes(user->M, PETSC_DECIDE, PETSC_DECIDE, n, n));
767 PetscCall(MatSetFromOptions(user->M));
768 PetscCall(MatMPIAIJSetPreallocation(user->M, 4, NULL, 4, NULL));
769 PetscCall(MatSeqAIJSetPreallocation(user->M, 4, NULL));
770 PetscCall(MatGetOwnershipRange(user->M, &istart, &iend));
771
772 for (i = istart; i < iend; i++) {
773 /* kron(Id,Av) */
774 iblock = i / user->mx;
775 j = iblock * user->mx + ((i + user->mx - 1) % user->mx);
776 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
777 j = iblock * user->mx + ((i + 1) % user->mx);
778 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
779
780 /* kron(Av,Id) */
781 j = (i + user->mx) % n;
782 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
783 j = (i + n - user->mx) % n;
784 PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
785 }
786 PetscCall(MatAssemblyBegin(user->M, MAT_FINAL_ASSEMBLY));
787 PetscCall(MatAssemblyEnd(user->M, MAT_FINAL_ASSEMBLY));
788
789 /* Generate 2D grid */
790 PetscCall(VecCreate(PETSC_COMM_WORLD, &XX));
791 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->q));
792 PetscCall(VecSetSizes(XX, PETSC_DECIDE, n));
793 PetscCall(VecSetSizes(user->q, PETSC_DECIDE, n * user->nt));
794 PetscCall(VecSetFromOptions(XX));
795 PetscCall(VecSetFromOptions(user->q));
796
797 PetscCall(VecDuplicate(XX, &YY));
798 PetscCall(VecDuplicate(XX, &XXwork));
799 PetscCall(VecDuplicate(XX, &YYwork));
800 PetscCall(VecDuplicate(XX, &user->d));
801 PetscCall(VecDuplicate(XX, &user->dwork));
802
803 PetscCall(VecGetOwnershipRange(XX, &istart, &iend));
804 for (linear_index = istart; linear_index < iend; linear_index++) {
805 i = linear_index % user->mx;
806 j = (linear_index - i) / user->mx;
807 vx = h * (i + 0.5);
808 vy = h * (j + 0.5);
809 PetscCall(VecSetValues(XX, 1, &linear_index, &vx, INSERT_VALUES));
810 PetscCall(VecSetValues(YY, 1, &linear_index, &vy, INSERT_VALUES));
811 }
812
813 PetscCall(VecAssemblyBegin(XX));
814 PetscCall(VecAssemblyEnd(XX));
815 PetscCall(VecAssemblyBegin(YY));
816 PetscCall(VecAssemblyEnd(YY));
817
818 /* Compute final density function yT
819 yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
820 yT = yT / (h^2*sum(yT)) */
821 PetscCall(VecCopy(XX, XXwork));
822 PetscCall(VecCopy(YY, YYwork));
823
824 PetscCall(VecShift(XXwork, -0.25));
825 PetscCall(VecShift(YYwork, -0.25));
826
827 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
828 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));
829
830 PetscCall(VecCopy(XXwork, user->dwork));
831 PetscCall(VecAXPY(user->dwork, 1.0, YYwork));
832 PetscCall(VecScale(user->dwork, -30.0));
833 PetscCall(VecExp(user->dwork));
834 PetscCall(VecCopy(user->dwork, user->d));
835
836 PetscCall(VecCopy(XX, XXwork));
837 PetscCall(VecCopy(YY, YYwork));
838
839 PetscCall(VecShift(XXwork, -0.75));
840 PetscCall(VecShift(YYwork, -0.75));
841
842 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
843 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));
844
845 PetscCall(VecCopy(XXwork, user->dwork));
846 PetscCall(VecAXPY(user->dwork, 1.0, YYwork));
847 PetscCall(VecScale(user->dwork, -30.0));
848 PetscCall(VecExp(user->dwork));
849
850 PetscCall(VecAXPY(user->d, 1.0, user->dwork));
851 PetscCall(VecShift(user->d, 1.0));
852 PetscCall(VecSum(user->d, &sum));
853 PetscCall(VecScale(user->d, 1.0 / (h * h * sum)));
854
855 /* Initial conditions of forward problem */
856 PetscCall(VecDuplicate(XX, &bc));
857 PetscCall(VecCopy(XX, XXwork));
858 PetscCall(VecCopy(YY, YYwork));
859
860 PetscCall(VecShift(XXwork, -0.5));
861 PetscCall(VecShift(YYwork, -0.5));
862
863 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
864 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));
865
866 PetscCall(VecWAXPY(bc, 1.0, XXwork, YYwork));
867 PetscCall(VecScale(bc, -50.0));
868 PetscCall(VecExp(bc));
869 PetscCall(VecShift(bc, 1.0));
870 PetscCall(VecSum(bc, &sum));
871 PetscCall(VecScale(bc, 1.0 / (h * h * sum)));
872
873 /* Create scatter from y to y_1,y_2,...,y_nt */
874 /* TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
875 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->yi_scatter));
876 PetscCall(VecCreate(PETSC_COMM_WORLD, &yi));
877 PetscCall(VecSetSizes(yi, PETSC_DECIDE, user->mx * user->mx));
878 PetscCall(VecSetFromOptions(yi));
879 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yi));
880 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yiwork));
881 PetscCall(VecDuplicateVecs(yi, user->nt, &user->ziwork));
882 for (i = 0; i < user->nt; i++) {
883 PetscCall(VecGetOwnershipRange(user->yi[i], &lo, &hi));
884 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_yi));
885 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + i * user->mx * user->mx, 1, &is_from_y));
886 PetscCall(VecScatterCreate(user->y, is_from_y, user->yi[i], is_to_yi, &user->yi_scatter[i]));
887 PetscCall(ISDestroy(&is_to_yi));
888 PetscCall(ISDestroy(&is_from_y));
889 }
890
891 /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
892 /* TODO: reorder for better parallelism */
893 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uxi_scatter));
894 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uyi_scatter));
895 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->ux_scatter));
896 PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uy_scatter));
897 PetscCall(PetscMalloc1(2 * user->nt * user->mx * user->mx, &user->ui_scatter));
898 PetscCall(VecCreate(PETSC_COMM_WORLD, &uxi));
899 PetscCall(VecCreate(PETSC_COMM_WORLD, &ui));
900 PetscCall(VecSetSizes(uxi, PETSC_DECIDE, user->mx * user->mx));
901 PetscCall(VecSetSizes(ui, PETSC_DECIDE, 2 * user->mx * user->mx));
902 PetscCall(VecSetFromOptions(uxi));
903 PetscCall(VecSetFromOptions(ui));
904 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uxi));
905 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uyi));
906 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uxiwork));
907 PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uyiwork));
908 PetscCall(VecDuplicateVecs(ui, user->nt, &user->ui));
909 PetscCall(VecDuplicateVecs(ui, user->nt, &user->uiwork));
910 for (i = 0; i < user->nt; i++) {
911 PetscCall(VecGetOwnershipRange(user->uxi[i], &lo, &hi));
912 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi));
913 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u));
914 PetscCall(VecScatterCreate(user->u, is_from_u, user->uxi[i], is_to_uxi, &user->uxi_scatter[i]));
915
916 PetscCall(ISDestroy(&is_to_uxi));
917 PetscCall(ISDestroy(&is_from_u));
918
919 PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi));
920 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi));
921 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + (2 * i + 1) * user->mx * user->mx, 1, &is_from_u));
922 PetscCall(VecScatterCreate(user->u, is_from_u, user->uyi[i], is_to_uyi, &user->uyi_scatter[i]));
923
924 PetscCall(ISDestroy(&is_to_uyi));
925 PetscCall(ISDestroy(&is_from_u));
926
927 PetscCall(VecGetOwnershipRange(user->uxi[i], &lo, &hi));
928 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi));
929 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_from_u));
930 PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uxi[i], is_to_uxi, &user->ux_scatter[i]));
931
932 PetscCall(ISDestroy(&is_to_uxi));
933 PetscCall(ISDestroy(&is_from_u));
934
935 PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi));
936 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi));
937 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + user->mx * user->mx, 1, &is_from_u));
938 PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uyi[i], is_to_uyi, &user->uy_scatter[i]));
939
940 PetscCall(ISDestroy(&is_to_uyi));
941 PetscCall(ISDestroy(&is_from_u));
942
943 PetscCall(VecGetOwnershipRange(user->ui[i], &lo, &hi));
944 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi));
945 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u));
946 PetscCall(VecScatterCreate(user->u, is_from_u, user->ui[i], is_to_uxi, &user->ui_scatter[i]));
947
948 PetscCall(ISDestroy(&is_to_uxi));
949 PetscCall(ISDestroy(&is_from_u));
950 }
951
952 /* RHS of forward problem */
953 PetscCall(MatMult(user->M, bc, user->yiwork[0]));
954 for (i = 1; i < user->nt; i++) PetscCall(VecSet(user->yiwork[i], 0.0));
955 PetscCall(Gather_yi(user->q, user->yiwork, user->yi_scatter, user->nt));
956
957 /* Compute true velocity field utrue */
958 PetscCall(VecDuplicate(user->u, &user->utrue));
959 for (i = 0; i < user->nt; i++) {
960 PetscCall(VecCopy(YY, user->uxi[i]));
961 PetscCall(VecScale(user->uxi[i], 150.0 * i * user->ht));
962 PetscCall(VecCopy(XX, user->uyi[i]));
963 PetscCall(VecShift(user->uyi[i], -10.0));
964 PetscCall(VecScale(user->uyi[i], 15.0 * i * user->ht));
965 }
966 PetscCall(Gather_uxi_uyi(user->utrue, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
967
968 /* Initial guess and reference model */
969 PetscCall(VecDuplicate(user->utrue, &user->ur));
970 for (i = 0; i < user->nt; i++) {
971 PetscCall(VecCopy(XX, user->uxi[i]));
972 PetscCall(VecShift(user->uxi[i], i * user->ht));
973 PetscCall(VecCopy(YY, user->uyi[i]));
974 PetscCall(VecShift(user->uyi[i], -i * user->ht));
975 }
976 PetscCall(Gather_uxi_uyi(user->ur, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
977
978 /* Generate regularization matrix L */
979 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->LT));
980 PetscCall(MatSetSizes(user->LT, PETSC_DECIDE, PETSC_DECIDE, 2 * n * user->nt, n * user->nt));
981 PetscCall(MatSetFromOptions(user->LT));
982 PetscCall(MatMPIAIJSetPreallocation(user->LT, 1, NULL, 1, NULL));
983 PetscCall(MatSeqAIJSetPreallocation(user->LT, 1, NULL));
984 PetscCall(MatGetOwnershipRange(user->LT, &istart, &iend));
985
986 for (i = istart; i < iend; i++) {
987 iblock = (i + n) / (2 * n);
988 j = i - iblock * n;
989 PetscCall(MatSetValues(user->LT, 1, &i, 1, &j, &user->gamma, INSERT_VALUES));
990 }
991
992 PetscCall(MatAssemblyBegin(user->LT, MAT_FINAL_ASSEMBLY));
993 PetscCall(MatAssemblyEnd(user->LT, MAT_FINAL_ASSEMBLY));
994
995 PetscCall(MatTranspose(user->LT, MAT_INITIAL_MATRIX, &user->L));
996
997 /* Build work vectors and matrices */
998 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->lwork));
999 PetscCall(VecSetType(user->lwork, VECMPI));
1000 PetscCall(VecSetSizes(user->lwork, PETSC_DECIDE, user->m));
1001 PetscCall(VecSetFromOptions(user->lwork));
1002
1003 PetscCall(MatDuplicate(user->Div, MAT_SHARE_NONZERO_PATTERN, &user->Divwork));
1004
1005 PetscCall(VecDuplicate(user->y, &user->ywork));
1006 PetscCall(VecDuplicate(user->u, &user->uwork));
1007 PetscCall(VecDuplicate(user->u, &user->vwork));
1008 PetscCall(VecDuplicate(user->u, &user->js_diag));
1009 PetscCall(VecDuplicate(user->c, &user->cwork));
1010
1011 /* Create matrix-free shell user->Js for computing A*x */
1012 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->Js));
1013 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT, (PetscErrorCodeFn *)StateMatMult));
1014 PetscCall(MatShellSetOperation(user->Js, MATOP_DUPLICATE, (PetscErrorCodeFn *)StateMatDuplicate));
1015 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatMultTranspose));
1016 PetscCall(MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)StateMatGetDiagonal));
1017
1018 /* Diagonal blocks of user->Js */
1019 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlock));
1020 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT, (PetscErrorCodeFn *)StateMatBlockMult));
1021 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatBlockMultTranspose));
1022
1023 /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
1024 D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
1025 This is an SOR preconditioner for user->JsBlock. */
1026 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlockPrec));
1027 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT, (PetscErrorCodeFn *)StateMatBlockPrecMult));
1028 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatBlockPrecMultTranspose));
1029
1030 /* Create a matrix-free shell user->Jd for computing B*x */
1031 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->n - user->m, user, &user->Jd));
1032 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT, (PetscErrorCodeFn *)DesignMatMult));
1033 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)DesignMatMultTranspose));
1034
1035 /* User-defined routines for computing user->Js\x and user->Js^T\x*/
1036 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsInv));
1037 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT, (PetscErrorCodeFn *)StateMatInvMult));
1038 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatInvTransposeMult));
1039
1040 /* Build matrices for SOR preconditioner */
1041 PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
1042 PetscCall(PetscMalloc1(5 * n, &user->C));
1043 PetscCall(PetscMalloc1(2 * n, &user->Cwork));
1044 for (i = 0; i < user->nt; i++) {
1045 PetscCall(MatDuplicate(user->Divxy[0], MAT_COPY_VALUES, &user->C[i]));
1046 PetscCall(MatDuplicate(user->Divxy[1], MAT_COPY_VALUES, &user->Cwork[i]));
1047
1048 PetscCall(MatDiagonalScale(user->C[i], NULL, user->uxi[i]));
1049 PetscCall(MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i]));
1050 PetscCall(MatAXPY(user->C[i], 1.0, user->Cwork[i], DIFFERENT_NONZERO_PATTERN));
1051 PetscCall(MatScale(user->C[i], user->ht));
1052 PetscCall(MatShift(user->C[i], 1.0));
1053 }
1054
1055 /* Solver options and tolerances */
1056 PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->solver));
1057 PetscCall(KSPSetType(user->solver, KSPGMRES));
1058 PetscCall(KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec));
1059 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500));
1060 /* PetscCall(KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500)); */
1061 PetscCall(KSPGetPC(user->solver, &user->prec));
1062 PetscCall(PCSetType(user->prec, PCSHELL));
1063
1064 PetscCall(PCShellSetApply(user->prec, StateMatBlockPrecMult));
1065 PetscCall(PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMultTranspose));
1066 PetscCall(PCShellSetContext(user->prec, user));
1067
1068 /* Compute true state function yt given ut */
1069 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ytrue));
1070 PetscCall(VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt));
1071 PetscCall(VecSetFromOptions(user->ytrue));
1072 user->c_formed = PETSC_TRUE;
1073 PetscCall(VecCopy(user->utrue, user->u)); /* Set u=utrue temporarily for StateMatInv */
1074 PetscCall(VecSet(user->ytrue, 0.0)); /* Initial guess */
1075 PetscCall(StateMatInvMult(user->Js, user->q, user->ytrue));
1076 PetscCall(VecCopy(user->ur, user->u)); /* Reset u=ur */
1077
1078 /* Initial guess y0 for state given u0 */
1079 PetscCall(StateMatInvMult(user->Js, user->q, user->y));
1080
1081 /* Data discretization */
1082 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Q));
1083 PetscCall(MatSetSizes(user->Q, PETSC_DECIDE, PETSC_DECIDE, user->mx * user->mx, user->m));
1084 PetscCall(MatSetFromOptions(user->Q));
1085 PetscCall(MatMPIAIJSetPreallocation(user->Q, 0, NULL, 1, NULL));
1086 PetscCall(MatSeqAIJSetPreallocation(user->Q, 1, NULL));
1087
1088 PetscCall(MatGetOwnershipRange(user->Q, &istart, &iend));
1089
1090 for (i = istart; i < iend; i++) {
1091 j = i + user->m - user->mx * user->mx;
1092 PetscCall(MatSetValues(user->Q, 1, &i, 1, &j, &one, INSERT_VALUES));
1093 }
1094
1095 PetscCall(MatAssemblyBegin(user->Q, MAT_FINAL_ASSEMBLY));
1096 PetscCall(MatAssemblyEnd(user->Q, MAT_FINAL_ASSEMBLY));
1097
1098 PetscCall(MatTranspose(user->Q, MAT_INITIAL_MATRIX, &user->QT));
1099
1100 PetscCall(VecDestroy(&XX));
1101 PetscCall(VecDestroy(&YY));
1102 PetscCall(VecDestroy(&XXwork));
1103 PetscCall(VecDestroy(&YYwork));
1104 PetscCall(VecDestroy(&yi));
1105 PetscCall(VecDestroy(&uxi));
1106 PetscCall(VecDestroy(&ui));
1107 PetscCall(VecDestroy(&bc));
1108
1109 /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1110 PetscCall(KSPSetFromOptions(user->solver));
1111 PetscFunctionReturn(PETSC_SUCCESS);
1112 }
1113
HyperbolicDestroy(AppCtx * user)1114 PetscErrorCode HyperbolicDestroy(AppCtx *user)
1115 {
1116 PetscInt i;
1117
1118 PetscFunctionBegin;
1119 PetscCall(MatDestroy(&user->Q));
1120 PetscCall(MatDestroy(&user->QT));
1121 PetscCall(MatDestroy(&user->Div));
1122 PetscCall(MatDestroy(&user->Divwork));
1123 PetscCall(MatDestroy(&user->Grad));
1124 PetscCall(MatDestroy(&user->L));
1125 PetscCall(MatDestroy(&user->LT));
1126 PetscCall(KSPDestroy(&user->solver));
1127 PetscCall(MatDestroy(&user->Js));
1128 PetscCall(MatDestroy(&user->Jd));
1129 PetscCall(MatDestroy(&user->JsBlockPrec));
1130 PetscCall(MatDestroy(&user->JsInv));
1131 PetscCall(MatDestroy(&user->JsBlock));
1132 PetscCall(MatDestroy(&user->Divxy[0]));
1133 PetscCall(MatDestroy(&user->Divxy[1]));
1134 PetscCall(MatDestroy(&user->Gradxy[0]));
1135 PetscCall(MatDestroy(&user->Gradxy[1]));
1136 PetscCall(MatDestroy(&user->M));
1137 for (i = 0; i < user->nt; i++) {
1138 PetscCall(MatDestroy(&user->C[i]));
1139 PetscCall(MatDestroy(&user->Cwork[i]));
1140 }
1141 PetscCall(PetscFree(user->C));
1142 PetscCall(PetscFree(user->Cwork));
1143 PetscCall(VecDestroy(&user->u));
1144 PetscCall(VecDestroy(&user->uwork));
1145 PetscCall(VecDestroy(&user->vwork));
1146 PetscCall(VecDestroy(&user->utrue));
1147 PetscCall(VecDestroy(&user->y));
1148 PetscCall(VecDestroy(&user->ywork));
1149 PetscCall(VecDestroy(&user->ytrue));
1150 PetscCall(VecDestroyVecs(user->nt, &user->yi));
1151 PetscCall(VecDestroyVecs(user->nt, &user->yiwork));
1152 PetscCall(VecDestroyVecs(user->nt, &user->ziwork));
1153 PetscCall(VecDestroyVecs(user->nt, &user->uxi));
1154 PetscCall(VecDestroyVecs(user->nt, &user->uyi));
1155 PetscCall(VecDestroyVecs(user->nt, &user->uxiwork));
1156 PetscCall(VecDestroyVecs(user->nt, &user->uyiwork));
1157 PetscCall(VecDestroyVecs(user->nt, &user->ui));
1158 PetscCall(VecDestroyVecs(user->nt, &user->uiwork));
1159 PetscCall(VecDestroy(&user->c));
1160 PetscCall(VecDestroy(&user->cwork));
1161 PetscCall(VecDestroy(&user->ur));
1162 PetscCall(VecDestroy(&user->q));
1163 PetscCall(VecDestroy(&user->d));
1164 PetscCall(VecDestroy(&user->dwork));
1165 PetscCall(VecDestroy(&user->lwork));
1166 PetscCall(VecDestroy(&user->js_diag));
1167 PetscCall(ISDestroy(&user->s_is));
1168 PetscCall(ISDestroy(&user->d_is));
1169 PetscCall(VecScatterDestroy(&user->state_scatter));
1170 PetscCall(VecScatterDestroy(&user->design_scatter));
1171 for (i = 0; i < user->nt; i++) {
1172 PetscCall(VecScatterDestroy(&user->uxi_scatter[i]));
1173 PetscCall(VecScatterDestroy(&user->uyi_scatter[i]));
1174 PetscCall(VecScatterDestroy(&user->ux_scatter[i]));
1175 PetscCall(VecScatterDestroy(&user->uy_scatter[i]));
1176 PetscCall(VecScatterDestroy(&user->ui_scatter[i]));
1177 PetscCall(VecScatterDestroy(&user->yi_scatter[i]));
1178 }
1179 PetscCall(PetscFree(user->uxi_scatter));
1180 PetscCall(PetscFree(user->uyi_scatter));
1181 PetscCall(PetscFree(user->ux_scatter));
1182 PetscCall(PetscFree(user->uy_scatter));
1183 PetscCall(PetscFree(user->ui_scatter));
1184 PetscCall(PetscFree(user->yi_scatter));
1185 PetscFunctionReturn(PETSC_SUCCESS);
1186 }
1187
HyperbolicMonitor(Tao tao,void * ptr)1188 PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1189 {
1190 Vec X;
1191 PetscReal unorm, ynorm;
1192 AppCtx *user = (AppCtx *)ptr;
1193
1194 PetscFunctionBegin;
1195 PetscCall(TaoGetSolution(tao, &X));
1196 PetscCall(Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
1197 PetscCall(VecAXPY(user->ywork, -1.0, user->ytrue));
1198 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue));
1199 PetscCall(VecNorm(user->uwork, NORM_2, &unorm));
1200 PetscCall(VecNorm(user->ywork, NORM_2, &ynorm));
1201 PetscCall(PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n", (double)unorm, (double)ynorm));
1202 PetscFunctionReturn(PETSC_SUCCESS);
1203 }
1204
1205 /*TEST
1206
1207 build:
1208 requires: !complex
1209
1210 test:
1211 requires: !single
1212 args: -tao_monitor_constraint_norm -tao_max_funcs 10 -tao_type lcl -tao_gatol 1.e-5
1213
1214 test:
1215 suffix: guess_pod
1216 requires: !single
1217 args: -tao_monitor_constraint_norm -tao_max_funcs 10 -tao_type lcl -ksp_guess_type pod -tao_gatol 1.e-5
1218
1219 TEST*/
1220