1 #include <petsc/private/taoimpl.h>
2
3 typedef struct {
4 PetscInt n; /* Number of variables */
5 PetscInt m; /* Number of constraints per time step */
6 PetscInt mx; /* grid points in each direction */
7 PetscInt nt; /* Number of time steps; as of now, must be divisible by 8 */
8 PetscInt ndata; /* Number of data points per sample */
9 PetscInt ns; /* Number of samples */
10 PetscInt *sample_times; /* Times of samples */
11 IS s_is;
12 IS d_is;
13
14 VecScatter state_scatter;
15 VecScatter design_scatter;
16 VecScatter *yi_scatter;
17 VecScatter *di_scatter;
18
19 Mat Js, Jd, JsBlockPrec, JsInv, JsBlock;
20 PetscBool jformed, dsg_formed;
21
22 PetscReal alpha; /* Regularization parameter */
23 PetscReal beta; /* Weight attributed to ||u||^2 in regularization functional */
24 PetscReal noise; /* Amount of noise to add to data */
25 PetscReal ht; /* Time step */
26
27 Mat Qblock, QblockT;
28 Mat L, LT;
29 Mat Div, Divwork;
30 Mat Grad;
31 Mat Av, Avwork, AvT;
32 Mat DSG;
33 Vec q;
34 Vec ur; /* reference */
35
36 Vec d;
37 Vec dwork;
38 Vec *di;
39
40 Vec y; /* state variables */
41 Vec ywork;
42
43 Vec ytrue;
44 Vec *yi, *yiwork;
45
46 Vec u; /* design variables */
47 Vec uwork;
48
49 Vec utrue;
50 Vec js_diag;
51 Vec c; /* constraint vector */
52 Vec cwork;
53
54 Vec lwork;
55 Vec S;
56 Vec Rwork, Swork, Twork;
57 Vec Av_u;
58
59 KSP solver;
60 PC prec;
61
62 PetscInt ksp_its;
63 PetscInt ksp_its_initial;
64 } AppCtx;
65
66 PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
67 PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
68 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
69 PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void *);
70 PetscErrorCode FormJacobianDesign(Tao, Vec, Mat, void *);
71 PetscErrorCode FormConstraints(Tao, Vec, Vec, void *);
72 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
73 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
74 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
75 PetscErrorCode ParabolicInitialize(AppCtx *user);
76 PetscErrorCode ParabolicDestroy(AppCtx *user);
77 PetscErrorCode ParabolicMonitor(Tao, void *);
78
79 PetscErrorCode StateMatMult(Mat, Vec, Vec);
80 PetscErrorCode StateMatBlockMult(Mat, Vec, Vec);
81 PetscErrorCode StateMatMultTranspose(Mat, Vec, Vec);
82 PetscErrorCode StateMatGetDiagonal(Mat, Vec);
83 PetscErrorCode StateMatDuplicate(Mat, MatDuplicateOption, Mat *);
84 PetscErrorCode StateMatInvMult(Mat, Vec, Vec);
85 PetscErrorCode StateMatInvTransposeMult(Mat, Vec, Vec);
86 PetscErrorCode StateMatBlockPrecMult(PC, Vec, Vec);
87
88 PetscErrorCode DesignMatMult(Mat, Vec, Vec);
89 PetscErrorCode DesignMatMultTranspose(Mat, Vec, Vec);
90
91 PetscErrorCode Gather_i(Vec, Vec *, VecScatter *, PetscInt);
92 PetscErrorCode Scatter_i(Vec, Vec *, VecScatter *, PetscInt);
93
94 static char help[] = "";
95
main(int argc,char ** argv)96 int main(int argc, char **argv)
97 {
98 Vec x, x0;
99 Tao tao;
100 AppCtx user;
101 IS is_allstate, is_alldesign;
102 PetscInt lo, hi, hi2, lo2, ksp_old;
103 PetscInt ntests = 1;
104 PetscInt i;
105 PetscLogStage stages[1];
106
107 PetscFunctionBeginUser;
108 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
109 user.mx = 8;
110 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "parabolic example", NULL);
111 PetscCall(PetscOptionsInt("-mx", "Number of grid points in each direction", "", user.mx, &user.mx, NULL));
112 user.nt = 8;
113 PetscCall(PetscOptionsInt("-nt", "Number of time steps", "", user.nt, &user.nt, NULL));
114 user.ndata = 64;
115 PetscCall(PetscOptionsInt("-ndata", "Numbers of data points per sample", "", user.ndata, &user.ndata, NULL));
116 user.ns = 8;
117 PetscCall(PetscOptionsInt("-ns", "Number of samples", "", user.ns, &user.ns, NULL));
118 user.alpha = 1.0;
119 PetscCall(PetscOptionsReal("-alpha", "Regularization parameter", "", user.alpha, &user.alpha, NULL));
120 user.beta = 0.01;
121 PetscCall(PetscOptionsReal("-beta", "Weight attributed to ||u||^2 in regularization functional", "", user.beta, &user.beta, NULL));
122 user.noise = 0.01;
123 PetscCall(PetscOptionsReal("-noise", "Amount of noise to add to data", "", user.noise, &user.noise, NULL));
124 PetscCall(PetscOptionsInt("-ntests", "Number of times to repeat TaoSolve", "", ntests, &ntests, NULL));
125 PetscOptionsEnd();
126
127 user.m = user.mx * user.mx * user.mx; /* number of constraints per time step */
128 user.n = user.m * (user.nt + 1); /* number of variables */
129 user.ht = (PetscReal)1 / user.nt;
130
131 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.u));
132 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.y));
133 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.c));
134 PetscCall(VecSetSizes(user.u, PETSC_DECIDE, user.n - user.m * user.nt));
135 PetscCall(VecSetSizes(user.y, PETSC_DECIDE, user.m * user.nt));
136 PetscCall(VecSetSizes(user.c, PETSC_DECIDE, user.m * user.nt));
137 PetscCall(VecSetFromOptions(user.u));
138 PetscCall(VecSetFromOptions(user.y));
139 PetscCall(VecSetFromOptions(user.c));
140
141 /* Create scatters for reduced spaces.
142 If the state vector y and design vector u are partitioned as
143 [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
144 then the solution vector x is organized as
145 [y_1; u_1; y_2; u_2; ...; y_np; u_np].
146 The index sets user.s_is and user.d_is correspond to the indices of the
147 state and design variables owned by the current processor.
148 */
149 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
150
151 PetscCall(VecGetOwnershipRange(user.y, &lo, &hi));
152 PetscCall(VecGetOwnershipRange(user.u, &lo2, &hi2));
153
154 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_allstate));
155 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + lo2, 1, &user.s_is));
156
157 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, lo2, 1, &is_alldesign));
158 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, hi + lo2, 1, &user.d_is));
159
160 PetscCall(VecSetSizes(x, hi - lo + hi2 - lo2, user.n));
161 PetscCall(VecSetFromOptions(x));
162
163 PetscCall(VecScatterCreate(x, user.s_is, user.y, is_allstate, &user.state_scatter));
164 PetscCall(VecScatterCreate(x, user.d_is, user.u, is_alldesign, &user.design_scatter));
165 PetscCall(ISDestroy(&is_alldesign));
166 PetscCall(ISDestroy(&is_allstate));
167
168 /* Create TAO solver and set desired solution method */
169 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
170 PetscCall(TaoSetType(tao, TAOLCL));
171
172 /* Set up initial vectors and matrices */
173 PetscCall(ParabolicInitialize(&user));
174
175 PetscCall(Gather(x, user.y, user.state_scatter, user.u, user.design_scatter));
176 PetscCall(VecDuplicate(x, &x0));
177 PetscCall(VecCopy(x, x0));
178
179 /* Set solution vector with an initial guess */
180 PetscCall(TaoSetSolution(tao, x));
181 PetscCall(TaoSetObjective(tao, FormFunction, &user));
182 PetscCall(TaoSetGradient(tao, NULL, FormGradient, &user));
183 PetscCall(TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user));
184
185 PetscCall(TaoSetJacobianStateRoutine(tao, user.Js, user.JsBlockPrec, user.JsInv, FormJacobianState, &user));
186 PetscCall(TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user));
187
188 PetscCall(TaoSetFromOptions(tao));
189 PetscCall(TaoSetStateDesignIS(tao, user.s_is, user.d_is));
190
191 /* SOLVE THE APPLICATION */
192 PetscCall(PetscLogStageRegister("Trials", &stages[0]));
193 PetscCall(PetscLogStagePush(stages[0]));
194 user.ksp_its_initial = user.ksp_its;
195 ksp_old = user.ksp_its;
196 for (i = 0; i < ntests; i++) {
197 PetscCall(TaoSolve(tao));
198 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP Iterations = %" PetscInt_FMT "\n", user.ksp_its - ksp_old));
199 PetscCall(VecCopy(x0, x));
200 PetscCall(TaoSetSolution(tao, x));
201 }
202 PetscCall(PetscLogStagePop());
203 PetscCall(PetscBarrier((PetscObject)x));
204 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations within initialization: "));
205 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its_initial));
206 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Total KSP iterations over %" PetscInt_FMT " trial(s): ", ntests));
207 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its));
208 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP iterations per trial: "));
209 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", (user.ksp_its - user.ksp_its_initial) / ntests));
210
211 PetscCall(TaoDestroy(&tao));
212 PetscCall(VecDestroy(&x));
213 PetscCall(VecDestroy(&x0));
214 PetscCall(ParabolicDestroy(&user));
215
216 PetscCall(PetscFinalize());
217 return 0;
218 }
219 /* ------------------------------------------------------------------- */
220 /*
221 dwork = Qy - d
222 lwork = L*(u-ur)
223 f = 1/2 * (dwork.dork + alpha*lwork.lwork)
224 */
FormFunction(Tao tao,Vec X,PetscReal * f,void * ptr)225 PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr)
226 {
227 PetscReal d1 = 0, d2 = 0;
228 PetscInt i, j;
229 AppCtx *user = (AppCtx *)ptr;
230
231 PetscFunctionBegin;
232 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
233 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt));
234 for (j = 0; j < user->ns; j++) {
235 i = user->sample_times[j];
236 PetscCall(MatMult(user->Qblock, user->yi[i], user->di[j]));
237 }
238 PetscCall(Gather_i(user->dwork, user->di, user->di_scatter, user->ns));
239 PetscCall(VecAXPY(user->dwork, -1.0, user->d));
240 PetscCall(VecDot(user->dwork, user->dwork, &d1));
241
242 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u));
243 PetscCall(MatMult(user->L, user->uwork, user->lwork));
244 PetscCall(VecDot(user->lwork, user->lwork, &d2));
245
246 *f = 0.5 * (d1 + user->alpha * d2);
247 PetscFunctionReturn(PETSC_SUCCESS);
248 }
249
250 /* ------------------------------------------------------------------- */
251 /*
252 state: g_s = Q' *(Qy - d)
253 design: g_d = alpha*L'*L*(u-ur)
254 */
FormGradient(Tao tao,Vec X,Vec G,void * ptr)255 PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr)
256 {
257 PetscInt i, j;
258 AppCtx *user = (AppCtx *)ptr;
259
260 PetscFunctionBegin;
261 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
262 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt));
263 for (j = 0; j < user->ns; j++) {
264 i = user->sample_times[j];
265 PetscCall(MatMult(user->Qblock, user->yi[i], user->di[j]));
266 }
267 PetscCall(Gather_i(user->dwork, user->di, user->di_scatter, user->ns));
268 PetscCall(VecAXPY(user->dwork, -1.0, user->d));
269 PetscCall(Scatter_i(user->dwork, user->di, user->di_scatter, user->ns));
270 PetscCall(VecSet(user->ywork, 0.0));
271 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt));
272 for (j = 0; j < user->ns; j++) {
273 i = user->sample_times[j];
274 PetscCall(MatMult(user->QblockT, user->di[j], user->yiwork[i]));
275 }
276 PetscCall(Gather_i(user->ywork, user->yiwork, user->yi_scatter, user->nt));
277
278 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u));
279 PetscCall(MatMult(user->L, user->uwork, user->lwork));
280 PetscCall(MatMult(user->LT, user->lwork, user->uwork));
281 PetscCall(VecScale(user->uwork, user->alpha));
282 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
283 PetscFunctionReturn(PETSC_SUCCESS);
284 }
285
FormFunctionGradient(Tao tao,Vec X,PetscReal * f,Vec G,void * ptr)286 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
287 {
288 PetscReal d1, d2;
289 PetscInt i, j;
290 AppCtx *user = (AppCtx *)ptr;
291
292 PetscFunctionBegin;
293 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
294 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt));
295 for (j = 0; j < user->ns; j++) {
296 i = user->sample_times[j];
297 PetscCall(MatMult(user->Qblock, user->yi[i], user->di[j]));
298 }
299 PetscCall(Gather_i(user->dwork, user->di, user->di_scatter, user->ns));
300 PetscCall(VecAXPY(user->dwork, -1.0, user->d));
301 PetscCall(VecDot(user->dwork, user->dwork, &d1));
302 PetscCall(Scatter_i(user->dwork, user->di, user->di_scatter, user->ns));
303 PetscCall(VecSet(user->ywork, 0.0));
304 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt));
305 for (j = 0; j < user->ns; j++) {
306 i = user->sample_times[j];
307 PetscCall(MatMult(user->QblockT, user->di[j], user->yiwork[i]));
308 }
309 PetscCall(Gather_i(user->ywork, user->yiwork, user->yi_scatter, user->nt));
310
311 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u));
312 PetscCall(MatMult(user->L, user->uwork, user->lwork));
313 PetscCall(VecDot(user->lwork, user->lwork, &d2));
314 PetscCall(MatMult(user->LT, user->lwork, user->uwork));
315 PetscCall(VecScale(user->uwork, user->alpha));
316 *f = 0.5 * (d1 + user->alpha * d2);
317
318 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
319 PetscFunctionReturn(PETSC_SUCCESS);
320 }
321
322 /* ------------------------------------------------------------------- */
323 /* A
324 MatShell object
325 */
FormJacobianState(Tao tao,Vec X,Mat J,Mat JPre,Mat JInv,void * ptr)326 PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
327 {
328 AppCtx *user = (AppCtx *)ptr;
329
330 PetscFunctionBegin;
331 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
332 PetscCall(VecSet(user->uwork, 0));
333 PetscCall(VecAXPY(user->uwork, -1.0, user->u));
334 PetscCall(VecExp(user->uwork));
335 PetscCall(MatMult(user->Av, user->uwork, user->Av_u));
336 PetscCall(VecCopy(user->Av_u, user->Swork));
337 PetscCall(VecReciprocal(user->Swork));
338 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN));
339 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Swork));
340 if (user->dsg_formed) {
341 PetscCall(MatProductNumeric(user->DSG));
342 } else {
343 PetscCall(MatMatMult(user->Divwork, user->Grad, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &user->DSG));
344 user->dsg_formed = PETSC_TRUE;
345 }
346
347 /* B = speye(nx^3) + ht*DSG; */
348 PetscCall(MatScale(user->DSG, user->ht));
349 PetscCall(MatShift(user->DSG, 1.0));
350 PetscFunctionReturn(PETSC_SUCCESS);
351 }
352
353 /* ------------------------------------------------------------------- */
354 /* B */
FormJacobianDesign(Tao tao,Vec X,Mat J,void * ptr)355 PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
356 {
357 AppCtx *user = (AppCtx *)ptr;
358
359 PetscFunctionBegin;
360 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
361 PetscFunctionReturn(PETSC_SUCCESS);
362 }
363
StateMatMult(Mat J_shell,Vec X,Vec Y)364 PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
365 {
366 PetscInt i;
367 AppCtx *user;
368
369 PetscFunctionBegin;
370 PetscCall(MatShellGetContext(J_shell, &user));
371 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt));
372 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0]));
373 for (i = 1; i < user->nt; i++) {
374 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i]));
375 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i - 1]));
376 }
377 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt));
378 PetscFunctionReturn(PETSC_SUCCESS);
379 }
380
StateMatMultTranspose(Mat J_shell,Vec X,Vec Y)381 PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
382 {
383 PetscInt i;
384 AppCtx *user;
385
386 PetscFunctionBegin;
387 PetscCall(MatShellGetContext(J_shell, &user));
388 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt));
389 for (i = 0; i < user->nt - 1; i++) {
390 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i]));
391 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i + 1]));
392 }
393 i = user->nt - 1;
394 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i]));
395 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt));
396 PetscFunctionReturn(PETSC_SUCCESS);
397 }
398
StateMatBlockMult(Mat J_shell,Vec X,Vec Y)399 PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
400 {
401 AppCtx *user;
402
403 PetscFunctionBegin;
404 PetscCall(MatShellGetContext(J_shell, &user));
405 PetscCall(MatMult(user->Grad, X, user->Swork));
406 PetscCall(VecPointwiseDivide(user->Swork, user->Swork, user->Av_u));
407 PetscCall(MatMult(user->Div, user->Swork, Y));
408 PetscCall(VecAYPX(Y, user->ht, X));
409 PetscFunctionReturn(PETSC_SUCCESS);
410 }
411
DesignMatMult(Mat J_shell,Vec X,Vec Y)412 PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
413 {
414 PetscInt i;
415 AppCtx *user;
416
417 PetscFunctionBegin;
418 PetscCall(MatShellGetContext(J_shell, &user));
419
420 /* sdiag(1./v) */
421 PetscCall(VecSet(user->uwork, 0));
422 PetscCall(VecAXPY(user->uwork, -1.0, user->u));
423 PetscCall(VecExp(user->uwork));
424
425 /* sdiag(1./((Av*(1./v)).^2)) */
426 PetscCall(MatMult(user->Av, user->uwork, user->Swork));
427 PetscCall(VecPointwiseMult(user->Swork, user->Swork, user->Swork));
428 PetscCall(VecReciprocal(user->Swork));
429
430 /* (Av * (sdiag(1./v) * b)) */
431 PetscCall(VecPointwiseMult(user->uwork, user->uwork, X));
432 PetscCall(MatMult(user->Av, user->uwork, user->Twork));
433
434 /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */
435 PetscCall(VecPointwiseMult(user->Swork, user->Twork, user->Swork));
436
437 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt));
438 for (i = 0; i < user->nt; i++) {
439 /* (sdiag(Grad*y(:,i)) */
440 PetscCall(MatMult(user->Grad, user->yi[i], user->Twork));
441
442 /* ht * Div * (sdiag(Grad*y(:,i)) * (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b)))) */
443 PetscCall(VecPointwiseMult(user->Twork, user->Twork, user->Swork));
444 PetscCall(MatMult(user->Div, user->Twork, user->yiwork[i]));
445 PetscCall(VecScale(user->yiwork[i], user->ht));
446 }
447 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt));
448 PetscFunctionReturn(PETSC_SUCCESS);
449 }
450
DesignMatMultTranspose(Mat J_shell,Vec X,Vec Y)451 PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
452 {
453 PetscInt i;
454 AppCtx *user;
455
456 PetscFunctionBegin;
457 PetscCall(MatShellGetContext(J_shell, &user));
458
459 /* sdiag(1./((Av*(1./v)).^2)) */
460 PetscCall(VecSet(user->uwork, 0));
461 PetscCall(VecAXPY(user->uwork, -1.0, user->u));
462 PetscCall(VecExp(user->uwork));
463 PetscCall(MatMult(user->Av, user->uwork, user->Rwork));
464 PetscCall(VecPointwiseMult(user->Rwork, user->Rwork, user->Rwork));
465 PetscCall(VecReciprocal(user->Rwork));
466
467 PetscCall(VecSet(Y, 0.0));
468 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt));
469 PetscCall(Scatter_i(X, user->yiwork, user->yi_scatter, user->nt));
470 for (i = 0; i < user->nt; i++) {
471 /* (Div' * b(:,i)) */
472 PetscCall(MatMult(user->Grad, user->yiwork[i], user->Swork));
473
474 /* sdiag(Grad*y(:,i)) */
475 PetscCall(MatMult(user->Grad, user->yi[i], user->Twork));
476
477 /* (sdiag(Grad*y(:,i)) * (Div' * b(:,i))) */
478 PetscCall(VecPointwiseMult(user->Twork, user->Swork, user->Twork));
479
480 /* (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i)))) */
481 PetscCall(VecPointwiseMult(user->Twork, user->Rwork, user->Twork));
482
483 /* (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */
484 PetscCall(MatMult(user->AvT, user->Twork, user->yiwork[i]));
485
486 /* sdiag(1./v) * (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */
487 PetscCall(VecPointwiseMult(user->yiwork[i], user->uwork, user->yiwork[i]));
488 PetscCall(VecAXPY(Y, user->ht, user->yiwork[i]));
489 }
490 PetscFunctionReturn(PETSC_SUCCESS);
491 }
492
StateMatBlockPrecMult(PC PC_shell,Vec X,Vec Y)493 PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
494 {
495 AppCtx *user;
496
497 PetscFunctionBegin;
498 PetscCall(PCShellGetContext(PC_shell, &user));
499
500 PetscCheck(user->dsg_formed, PETSC_COMM_SELF, PETSC_ERR_SUP, "DSG not formed");
501 PetscCall(MatSOR(user->DSG, X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y));
502 PetscFunctionReturn(PETSC_SUCCESS);
503 }
504
StateMatInvMult(Mat J_shell,Vec X,Vec Y)505 PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
506 {
507 AppCtx *user;
508 PetscInt its, i;
509
510 PetscFunctionBegin;
511 PetscCall(MatShellGetContext(J_shell, &user));
512
513 if (Y == user->ytrue) {
514 /* First solve is done with true solution to set up problem */
515 PetscCall(KSPSetTolerances(user->solver, 1e-8, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
516 } else {
517 PetscCall(KSPSetTolerances(user->solver, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
518 }
519
520 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt));
521 PetscCall(KSPSolve(user->solver, user->yi[0], user->yiwork[0]));
522 PetscCall(KSPGetIterationNumber(user->solver, &its));
523 user->ksp_its = user->ksp_its + its;
524
525 for (i = 1; i < user->nt; i++) {
526 PetscCall(VecAXPY(user->yi[i], 1.0, user->yiwork[i - 1]));
527 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i]));
528 PetscCall(KSPGetIterationNumber(user->solver, &its));
529 user->ksp_its = user->ksp_its + its;
530 }
531 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt));
532 PetscFunctionReturn(PETSC_SUCCESS);
533 }
534
StateMatInvTransposeMult(Mat J_shell,Vec X,Vec Y)535 PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
536 {
537 AppCtx *user;
538 PetscInt its, i;
539
540 PetscFunctionBegin;
541 PetscCall(MatShellGetContext(J_shell, &user));
542
543 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt));
544
545 i = user->nt - 1;
546 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i]));
547
548 PetscCall(KSPGetIterationNumber(user->solver, &its));
549 user->ksp_its = user->ksp_its + its;
550
551 for (i = user->nt - 2; i >= 0; i--) {
552 PetscCall(VecAXPY(user->yi[i], 1.0, user->yiwork[i + 1]));
553 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i]));
554
555 PetscCall(KSPGetIterationNumber(user->solver, &its));
556 user->ksp_its = user->ksp_its + its;
557 }
558
559 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt));
560 PetscFunctionReturn(PETSC_SUCCESS);
561 }
562
StateMatDuplicate(Mat J_shell,MatDuplicateOption opt,Mat * new_shell)563 PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
564 {
565 AppCtx *user;
566
567 PetscFunctionBegin;
568 PetscCall(MatShellGetContext(J_shell, &user));
569
570 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, new_shell));
571 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT, (PetscErrorCodeFn *)StateMatMult));
572 PetscCall(MatShellSetOperation(*new_shell, MATOP_DUPLICATE, (PetscErrorCodeFn *)StateMatDuplicate));
573 PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatMultTranspose));
574 PetscCall(MatShellSetOperation(*new_shell, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)StateMatGetDiagonal));
575 PetscFunctionReturn(PETSC_SUCCESS);
576 }
577
StateMatGetDiagonal(Mat J_shell,Vec X)578 PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
579 {
580 AppCtx *user;
581
582 PetscFunctionBegin;
583 PetscCall(MatShellGetContext(J_shell, &user));
584 PetscCall(VecCopy(user->js_diag, X));
585 PetscFunctionReturn(PETSC_SUCCESS);
586 }
587
FormConstraints(Tao tao,Vec X,Vec C,void * ptr)588 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
589 {
590 /* con = Ay - q, A = [B 0 0 ... 0;
591 -I B 0 ... 0;
592 0 -I B ... 0;
593 ... ;
594 0 ... -I B]
595 B = ht * Div * Sigma * Grad + eye */
596 PetscInt i;
597 AppCtx *user = (AppCtx *)ptr;
598
599 PetscFunctionBegin;
600 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
601 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt));
602 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0]));
603 for (i = 1; i < user->nt; i++) {
604 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i]));
605 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i - 1]));
606 }
607 PetscCall(Gather_i(C, user->yiwork, user->yi_scatter, user->nt));
608 PetscCall(VecAXPY(C, -1.0, user->q));
609 PetscFunctionReturn(PETSC_SUCCESS);
610 }
611
Scatter(Vec x,Vec state,VecScatter s_scat,Vec design,VecScatter d_scat)612 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
613 {
614 PetscFunctionBegin;
615 PetscCall(VecScatterBegin(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD));
616 PetscCall(VecScatterEnd(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD));
617 PetscCall(VecScatterBegin(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD));
618 PetscCall(VecScatterEnd(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD));
619 PetscFunctionReturn(PETSC_SUCCESS);
620 }
621
Scatter_i(Vec y,Vec * yi,VecScatter * scat,PetscInt nt)622 PetscErrorCode Scatter_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
623 {
624 PetscInt i;
625
626 PetscFunctionBegin;
627 for (i = 0; i < nt; i++) {
628 PetscCall(VecScatterBegin(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD));
629 PetscCall(VecScatterEnd(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD));
630 }
631 PetscFunctionReturn(PETSC_SUCCESS);
632 }
633
Gather(Vec x,Vec state,VecScatter s_scat,Vec design,VecScatter d_scat)634 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
635 {
636 PetscFunctionBegin;
637 PetscCall(VecScatterBegin(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE));
638 PetscCall(VecScatterEnd(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE));
639 PetscCall(VecScatterBegin(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE));
640 PetscCall(VecScatterEnd(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE));
641 PetscFunctionReturn(PETSC_SUCCESS);
642 }
643
Gather_i(Vec y,Vec * yi,VecScatter * scat,PetscInt nt)644 PetscErrorCode Gather_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
645 {
646 PetscInt i;
647
648 PetscFunctionBegin;
649 for (i = 0; i < nt; i++) {
650 PetscCall(VecScatterBegin(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE));
651 PetscCall(VecScatterEnd(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE));
652 }
653 PetscFunctionReturn(PETSC_SUCCESS);
654 }
655
ParabolicInitialize(AppCtx * user)656 PetscErrorCode ParabolicInitialize(AppCtx *user)
657 {
658 PetscInt m, n, i, j, k, linear_index, istart, iend, iblock, lo, hi, lo2, hi2;
659 Vec XX, YY, ZZ, XXwork, YYwork, ZZwork, UTwork, yi, di, bc;
660 PetscReal *x, *y, *z;
661 PetscReal h, stime;
662 PetscScalar hinv, neg_hinv, half = 0.5, sqrt_beta;
663 PetscInt im, indx1, indx2, indy1, indy2, indz1, indz2, nx, ny, nz;
664 PetscReal xri, yri, zri, xim, yim, zim, dx1, dx2, dy1, dy2, dz1, dz2, Dx, Dy, Dz;
665 PetscScalar v, vx, vy, vz;
666 IS is_from_y, is_to_yi, is_from_d, is_to_di;
667 /* Data locations */
668 PetscScalar xr[64] = {0.4970, 0.8498, 0.7814, 0.6268, 0.7782, 0.6402, 0.3617, 0.3160, 0.3610, 0.5298, 0.6987, 0.3331, 0.7962, 0.5596, 0.3866, 0.6774, 0.5407, 0.4518, 0.6702, 0.6061, 0.7580, 0.8997,
669 0.5198, 0.8326, 0.2138, 0.9198, 0.3000, 0.2833, 0.8288, 0.7076, 0.1820, 0.0728, 0.8447, 0.2367, 0.3239, 0.6413, 0.3114, 0.4731, 0.1192, 0.9273, 0.5724, 0.4331, 0.5136, 0.3547,
670 0.4413, 0.2602, 0.5698, 0.7278, 0.5261, 0.6230, 0.2454, 0.3948, 0.7479, 0.6582, 0.4660, 0.5594, 0.7574, 0.1143, 0.5900, 0.1065, 0.4260, 0.3294, 0.8276, 0.0756};
671
672 PetscScalar yr[64] = {0.7345, 0.9120, 0.9288, 0.7528, 0.4463, 0.4985, 0.2497, 0.6256, 0.3425, 0.9026, 0.6983, 0.4230, 0.7140, 0.2970, 0.4474, 0.8792, 0.6604, 0.2485, 0.7968, 0.6127, 0.1796, 0.2437,
673 0.5938, 0.6137, 0.3867, 0.5658, 0.4575, 0.1009, 0.0863, 0.3361, 0.0738, 0.3985, 0.6602, 0.1437, 0.0934, 0.5983, 0.5950, 0.0763, 0.0768, 0.2288, 0.5761, 0.1129, 0.3841, 0.6150,
674 0.6904, 0.6686, 0.1361, 0.4601, 0.4491, 0.3716, 0.1969, 0.6537, 0.6743, 0.6991, 0.4811, 0.5480, 0.1684, 0.4569, 0.6889, 0.8437, 0.3015, 0.2854, 0.8199, 0.2658};
675
676 PetscScalar zr[64] = {0.7668, 0.8573, 0.2654, 0.2719, 0.1060, 0.1311, 0.6232, 0.2295, 0.8009, 0.2147, 0.2119, 0.9325, 0.4473, 0.3600, 0.3374, 0.3819, 0.4066, 0.5801, 0.1673, 0.0959, 0.4638, 0.8236,
677 0.8800, 0.2939, 0.2028, 0.8262, 0.2706, 0.6276, 0.9085, 0.6443, 0.8241, 0.0712, 0.1824, 0.7789, 0.4389, 0.8415, 0.7055, 0.6639, 0.3653, 0.2078, 0.1987, 0.2297, 0.4321, 0.8115,
678 0.4915, 0.7764, 0.4657, 0.4627, 0.4569, 0.4232, 0.8514, 0.0674, 0.3227, 0.1055, 0.6690, 0.6313, 0.9226, 0.5461, 0.4126, 0.2364, 0.6096, 0.7042, 0.3914, 0.0711};
679
680 PetscFunctionBegin;
681 PetscCall(PetscMalloc1(user->mx, &x));
682 PetscCall(PetscMalloc1(user->mx, &y));
683 PetscCall(PetscMalloc1(user->mx, &z));
684 user->jformed = PETSC_FALSE;
685 user->dsg_formed = PETSC_FALSE;
686
687 n = user->mx * user->mx * user->mx;
688 m = 3 * user->mx * user->mx * (user->mx - 1);
689 sqrt_beta = PetscSqrtScalar(user->beta);
690
691 user->ksp_its = 0;
692 user->ksp_its_initial = 0;
693
694 stime = (PetscReal)user->nt / user->ns;
695 PetscCall(PetscMalloc1(user->ns, &user->sample_times));
696 for (i = 0; i < user->ns; i++) user->sample_times[i] = (PetscInt)(stime * ((PetscReal)i + 1.0) - 0.5);
697
698 PetscCall(VecCreate(PETSC_COMM_WORLD, &XX));
699 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->q));
700 PetscCall(VecSetSizes(XX, PETSC_DECIDE, n));
701 PetscCall(VecSetSizes(user->q, PETSC_DECIDE, n * user->nt));
702 PetscCall(VecSetFromOptions(XX));
703 PetscCall(VecSetFromOptions(user->q));
704
705 PetscCall(VecDuplicate(XX, &YY));
706 PetscCall(VecDuplicate(XX, &ZZ));
707 PetscCall(VecDuplicate(XX, &XXwork));
708 PetscCall(VecDuplicate(XX, &YYwork));
709 PetscCall(VecDuplicate(XX, &ZZwork));
710 PetscCall(VecDuplicate(XX, &UTwork));
711 PetscCall(VecDuplicate(XX, &user->utrue));
712 PetscCall(VecDuplicate(XX, &bc));
713
714 /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */
715 h = 1.0 / user->mx;
716 hinv = user->mx;
717 neg_hinv = -hinv;
718
719 PetscCall(VecGetOwnershipRange(XX, &istart, &iend));
720 for (linear_index = istart; linear_index < iend; linear_index++) {
721 i = linear_index % user->mx;
722 j = ((linear_index - i) / user->mx) % user->mx;
723 k = ((linear_index - i) / user->mx - j) / user->mx;
724 vx = h * (i + 0.5);
725 vy = h * (j + 0.5);
726 vz = h * (k + 0.5);
727 PetscCall(VecSetValues(XX, 1, &linear_index, &vx, INSERT_VALUES));
728 PetscCall(VecSetValues(YY, 1, &linear_index, &vy, INSERT_VALUES));
729 PetscCall(VecSetValues(ZZ, 1, &linear_index, &vz, INSERT_VALUES));
730 if ((vx < 0.6) && (vx > 0.4) && (vy < 0.6) && (vy > 0.4) && (vy < 0.6) && (vz < 0.6) && (vz > 0.4)) {
731 v = 1000.0;
732 PetscCall(VecSetValues(bc, 1, &linear_index, &v, INSERT_VALUES));
733 }
734 }
735
736 PetscCall(VecAssemblyBegin(XX));
737 PetscCall(VecAssemblyEnd(XX));
738 PetscCall(VecAssemblyBegin(YY));
739 PetscCall(VecAssemblyEnd(YY));
740 PetscCall(VecAssemblyBegin(ZZ));
741 PetscCall(VecAssemblyEnd(ZZ));
742 PetscCall(VecAssemblyBegin(bc));
743 PetscCall(VecAssemblyEnd(bc));
744
745 /* Compute true parameter function
746 ut = 0.5 + exp(-10*((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)) */
747 PetscCall(VecCopy(XX, XXwork));
748 PetscCall(VecCopy(YY, YYwork));
749 PetscCall(VecCopy(ZZ, ZZwork));
750
751 PetscCall(VecShift(XXwork, -0.5));
752 PetscCall(VecShift(YYwork, -0.5));
753 PetscCall(VecShift(ZZwork, -0.5));
754
755 PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
756 PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));
757 PetscCall(VecPointwiseMult(ZZwork, ZZwork, ZZwork));
758
759 PetscCall(VecCopy(XXwork, user->utrue));
760 PetscCall(VecAXPY(user->utrue, 1.0, YYwork));
761 PetscCall(VecAXPY(user->utrue, 1.0, ZZwork));
762 PetscCall(VecScale(user->utrue, -10.0));
763 PetscCall(VecExp(user->utrue));
764 PetscCall(VecShift(user->utrue, 0.5));
765
766 PetscCall(VecDestroy(&XX));
767 PetscCall(VecDestroy(&YY));
768 PetscCall(VecDestroy(&ZZ));
769 PetscCall(VecDestroy(&XXwork));
770 PetscCall(VecDestroy(&YYwork));
771 PetscCall(VecDestroy(&ZZwork));
772 PetscCall(VecDestroy(&UTwork));
773
774 /* Initial guess and reference model */
775 PetscCall(VecDuplicate(user->utrue, &user->ur));
776 PetscCall(VecSet(user->ur, 0.0));
777
778 /* Generate Grad matrix */
779 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Grad));
780 PetscCall(MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, m, n));
781 PetscCall(MatSetFromOptions(user->Grad));
782 PetscCall(MatMPIAIJSetPreallocation(user->Grad, 2, NULL, 2, NULL));
783 PetscCall(MatSeqAIJSetPreallocation(user->Grad, 2, NULL));
784 PetscCall(MatGetOwnershipRange(user->Grad, &istart, &iend));
785
786 for (i = istart; i < iend; i++) {
787 if (i < m / 3) {
788 iblock = i / (user->mx - 1);
789 j = iblock * user->mx + (i % (user->mx - 1));
790 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES));
791 j = j + 1;
792 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES));
793 }
794 if (i >= m / 3 && i < 2 * m / 3) {
795 iblock = (i - m / 3) / (user->mx * (user->mx - 1));
796 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1)));
797 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES));
798 j = j + user->mx;
799 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES));
800 }
801 if (i >= 2 * m / 3) {
802 j = i - 2 * m / 3;
803 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES));
804 j = j + user->mx * user->mx;
805 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES));
806 }
807 }
808
809 PetscCall(MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY));
810 PetscCall(MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY));
811
812 /* Generate arithmetic averaging matrix Av */
813 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Av));
814 PetscCall(MatSetSizes(user->Av, PETSC_DECIDE, PETSC_DECIDE, m, n));
815 PetscCall(MatSetFromOptions(user->Av));
816 PetscCall(MatMPIAIJSetPreallocation(user->Av, 2, NULL, 2, NULL));
817 PetscCall(MatSeqAIJSetPreallocation(user->Av, 2, NULL));
818 PetscCall(MatGetOwnershipRange(user->Av, &istart, &iend));
819
820 for (i = istart; i < iend; i++) {
821 if (i < m / 3) {
822 iblock = i / (user->mx - 1);
823 j = iblock * user->mx + (i % (user->mx - 1));
824 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES));
825 j = j + 1;
826 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES));
827 }
828 if (i >= m / 3 && i < 2 * m / 3) {
829 iblock = (i - m / 3) / (user->mx * (user->mx - 1));
830 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1)));
831 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES));
832 j = j + user->mx;
833 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES));
834 }
835 if (i >= 2 * m / 3) {
836 j = i - 2 * m / 3;
837 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES));
838 j = j + user->mx * user->mx;
839 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES));
840 }
841 }
842
843 PetscCall(MatAssemblyBegin(user->Av, MAT_FINAL_ASSEMBLY));
844 PetscCall(MatAssemblyEnd(user->Av, MAT_FINAL_ASSEMBLY));
845
846 /* Generate transpose of averaging matrix Av */
847 PetscCall(MatTranspose(user->Av, MAT_INITIAL_MATRIX, &user->AvT));
848
849 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->L));
850 PetscCall(MatSetSizes(user->L, PETSC_DECIDE, PETSC_DECIDE, m + n, n));
851 PetscCall(MatSetFromOptions(user->L));
852 PetscCall(MatMPIAIJSetPreallocation(user->L, 2, NULL, 2, NULL));
853 PetscCall(MatSeqAIJSetPreallocation(user->L, 2, NULL));
854 PetscCall(MatGetOwnershipRange(user->L, &istart, &iend));
855
856 for (i = istart; i < iend; i++) {
857 if (i < m / 3) {
858 iblock = i / (user->mx - 1);
859 j = iblock * user->mx + (i % (user->mx - 1));
860 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES));
861 j = j + 1;
862 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES));
863 }
864 if (i >= m / 3 && i < 2 * m / 3) {
865 iblock = (i - m / 3) / (user->mx * (user->mx - 1));
866 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1)));
867 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES));
868 j = j + user->mx;
869 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES));
870 }
871 if (i >= 2 * m / 3 && i < m) {
872 j = i - 2 * m / 3;
873 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES));
874 j = j + user->mx * user->mx;
875 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES));
876 }
877 if (i >= m) {
878 j = i - m;
879 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &sqrt_beta, INSERT_VALUES));
880 }
881 }
882
883 PetscCall(MatAssemblyBegin(user->L, MAT_FINAL_ASSEMBLY));
884 PetscCall(MatAssemblyEnd(user->L, MAT_FINAL_ASSEMBLY));
885
886 PetscCall(MatScale(user->L, PetscPowScalar(h, 1.5)));
887
888 /* Generate Div matrix */
889 PetscCall(MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div));
890
891 /* Build work vectors and matrices */
892 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->S));
893 PetscCall(VecSetSizes(user->S, PETSC_DECIDE, user->mx * user->mx * (user->mx - 1) * 3));
894 PetscCall(VecSetFromOptions(user->S));
895
896 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->lwork));
897 PetscCall(VecSetSizes(user->lwork, PETSC_DECIDE, m + user->mx * user->mx * user->mx));
898 PetscCall(VecSetFromOptions(user->lwork));
899
900 PetscCall(MatDuplicate(user->Div, MAT_SHARE_NONZERO_PATTERN, &user->Divwork));
901 PetscCall(MatDuplicate(user->Av, MAT_SHARE_NONZERO_PATTERN, &user->Avwork));
902
903 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->d));
904 PetscCall(VecSetSizes(user->d, PETSC_DECIDE, user->ndata * user->nt));
905 PetscCall(VecSetFromOptions(user->d));
906
907 PetscCall(VecDuplicate(user->S, &user->Swork));
908 PetscCall(VecDuplicate(user->S, &user->Av_u));
909 PetscCall(VecDuplicate(user->S, &user->Twork));
910 PetscCall(VecDuplicate(user->S, &user->Rwork));
911 PetscCall(VecDuplicate(user->y, &user->ywork));
912 PetscCall(VecDuplicate(user->u, &user->uwork));
913 PetscCall(VecDuplicate(user->u, &user->js_diag));
914 PetscCall(VecDuplicate(user->c, &user->cwork));
915 PetscCall(VecDuplicate(user->d, &user->dwork));
916
917 /* Create matrix-free shell user->Js for computing A*x */
918 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m * user->nt, user, &user->Js));
919 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT, (PetscErrorCodeFn *)StateMatMult));
920 PetscCall(MatShellSetOperation(user->Js, MATOP_DUPLICATE, (PetscErrorCodeFn *)StateMatDuplicate));
921 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatMultTranspose));
922 PetscCall(MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)StateMatGetDiagonal));
923
924 /* Diagonal blocks of user->Js */
925 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsBlock));
926 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT, (PetscErrorCodeFn *)StateMatBlockMult));
927 /* Blocks are symmetric */
928 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatBlockMult));
929
930 /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
931 D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
932 This is an SSOR preconditioner for user->JsBlock. */
933 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsBlockPrec));
934 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT, (PetscErrorCodeFn *)StateMatBlockPrecMult));
935 /* JsBlockPrec is symmetric */
936 PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatBlockPrecMult));
937 PetscCall(MatSetOption(user->JsBlockPrec, MAT_SYMMETRIC, PETSC_TRUE));
938 PetscCall(MatSetOption(user->JsBlockPrec, MAT_SYMMETRY_ETERNAL, PETSC_TRUE));
939
940 /* Create a matrix-free shell user->Jd for computing B*x */
941 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m, user, &user->Jd));
942 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT, (PetscErrorCodeFn *)DesignMatMult));
943 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)DesignMatMultTranspose));
944
945 /* User-defined routines for computing user->Js\x and user->Js^T\x*/
946 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m * user->nt, user, &user->JsInv));
947 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT, (PetscErrorCodeFn *)StateMatInvMult));
948 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatInvTransposeMult));
949
950 /* Solver options and tolerances */
951 PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->solver));
952 PetscCall(KSPSetType(user->solver, KSPCG));
953 PetscCall(KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec));
954 PetscCall(KSPSetInitialGuessNonzero(user->solver, PETSC_FALSE));
955 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500));
956 PetscCall(KSPSetFromOptions(user->solver));
957 PetscCall(KSPGetPC(user->solver, &user->prec));
958 PetscCall(PCSetType(user->prec, PCSHELL));
959
960 PetscCall(PCShellSetApply(user->prec, StateMatBlockPrecMult));
961 PetscCall(PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMult));
962 PetscCall(PCShellSetContext(user->prec, user));
963
964 /* Create scatter from y to y_1,y_2,...,y_nt */
965 PetscCall(PetscMalloc1(user->nt * user->m, &user->yi_scatter));
966 PetscCall(VecCreate(PETSC_COMM_WORLD, &yi));
967 PetscCall(VecSetSizes(yi, PETSC_DECIDE, user->mx * user->mx * user->mx));
968 PetscCall(VecSetFromOptions(yi));
969 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yi));
970 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yiwork));
971
972 PetscCall(VecGetOwnershipRange(user->y, &lo2, &hi2));
973 istart = 0;
974 for (i = 0; i < user->nt; i++) {
975 PetscCall(VecGetOwnershipRange(user->yi[i], &lo, &hi));
976 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_yi));
977 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo2 + istart, 1, &is_from_y));
978 PetscCall(VecScatterCreate(user->y, is_from_y, user->yi[i], is_to_yi, &user->yi_scatter[i]));
979 istart = istart + hi - lo;
980 PetscCall(ISDestroy(&is_to_yi));
981 PetscCall(ISDestroy(&is_from_y));
982 }
983 PetscCall(VecDestroy(&yi));
984
985 /* Create scatter from d to d_1,d_2,...,d_ns */
986 PetscCall(PetscMalloc1(user->ns * user->ndata, &user->di_scatter));
987 PetscCall(VecCreate(PETSC_COMM_WORLD, &di));
988 PetscCall(VecSetSizes(di, PETSC_DECIDE, user->ndata));
989 PetscCall(VecSetFromOptions(di));
990 PetscCall(VecDuplicateVecs(di, user->ns, &user->di));
991 PetscCall(VecGetOwnershipRange(user->d, &lo2, &hi2));
992 istart = 0;
993 for (i = 0; i < user->ns; i++) {
994 PetscCall(VecGetOwnershipRange(user->di[i], &lo, &hi));
995 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_di));
996 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo2 + istart, 1, &is_from_d));
997 PetscCall(VecScatterCreate(user->d, is_from_d, user->di[i], is_to_di, &user->di_scatter[i]));
998 istart = istart + hi - lo;
999 PetscCall(ISDestroy(&is_to_di));
1000 PetscCall(ISDestroy(&is_from_d));
1001 }
1002 PetscCall(VecDestroy(&di));
1003
1004 /* Assemble RHS of forward problem */
1005 PetscCall(VecCopy(bc, user->yiwork[0]));
1006 for (i = 1; i < user->nt; i++) PetscCall(VecSet(user->yiwork[i], 0.0));
1007 PetscCall(Gather_i(user->q, user->yiwork, user->yi_scatter, user->nt));
1008 PetscCall(VecDestroy(&bc));
1009
1010 /* Compute true state function yt given ut */
1011 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ytrue));
1012 PetscCall(VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt));
1013 PetscCall(VecSetFromOptions(user->ytrue));
1014
1015 /* First compute Av_u = Av*exp(-u) */
1016 PetscCall(VecSet(user->uwork, 0));
1017 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue)); /* Note: user->utrue */
1018 PetscCall(VecExp(user->uwork));
1019 PetscCall(MatMult(user->Av, user->uwork, user->Av_u));
1020
1021 /* Symbolic DSG = Div * Grad */
1022 PetscCall(MatProductCreate(user->Div, user->Grad, NULL, &user->DSG));
1023 PetscCall(MatProductSetType(user->DSG, MATPRODUCT_AB));
1024 PetscCall(MatProductSetAlgorithm(user->DSG, "default"));
1025 PetscCall(MatProductSetFill(user->DSG, PETSC_DETERMINE));
1026 PetscCall(MatProductSetFromOptions(user->DSG));
1027 PetscCall(MatProductSymbolic(user->DSG));
1028
1029 user->dsg_formed = PETSC_TRUE;
1030
1031 /* Next form DSG = Div*Grad */
1032 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN));
1033 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Av_u));
1034 if (user->dsg_formed) {
1035 PetscCall(MatProductNumeric(user->DSG));
1036 } else {
1037 PetscCall(MatMatMult(user->Div, user->Grad, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &user->DSG));
1038 user->dsg_formed = PETSC_TRUE;
1039 }
1040 /* B = speye(nx^3) + ht*DSG; */
1041 PetscCall(MatScale(user->DSG, user->ht));
1042 PetscCall(MatShift(user->DSG, 1.0));
1043
1044 /* Now solve for ytrue */
1045 PetscCall(StateMatInvMult(user->Js, user->q, user->ytrue));
1046
1047 /* Initial guess y0 for state given u0 */
1048
1049 /* First compute Av_u = Av*exp(-u) */
1050 PetscCall(VecSet(user->uwork, 0));
1051 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); /* Note: user->u */
1052 PetscCall(VecExp(user->uwork));
1053 PetscCall(MatMult(user->Av, user->uwork, user->Av_u));
1054
1055 /* Next form DSG = Div*S*Grad */
1056 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN));
1057 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Av_u));
1058 if (user->dsg_formed) {
1059 PetscCall(MatProductNumeric(user->DSG));
1060 } else {
1061 PetscCall(MatMatMult(user->Div, user->Grad, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &user->DSG));
1062
1063 user->dsg_formed = PETSC_TRUE;
1064 }
1065 /* B = speye(nx^3) + ht*DSG; */
1066 PetscCall(MatScale(user->DSG, user->ht));
1067 PetscCall(MatShift(user->DSG, 1.0));
1068
1069 /* Now solve for y */
1070 PetscCall(StateMatInvMult(user->Js, user->q, user->y));
1071
1072 /* Construct projection matrix Q, a block diagonal matrix consisting of nt copies of Qblock along the diagonal */
1073 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Qblock));
1074 PetscCall(MatSetSizes(user->Qblock, PETSC_DECIDE, PETSC_DECIDE, user->ndata, n));
1075 PetscCall(MatSetFromOptions(user->Qblock));
1076 PetscCall(MatMPIAIJSetPreallocation(user->Qblock, 8, NULL, 8, NULL));
1077 PetscCall(MatSeqAIJSetPreallocation(user->Qblock, 8, NULL));
1078
1079 for (i = 0; i < user->mx; i++) {
1080 x[i] = h * (i + 0.5);
1081 y[i] = h * (i + 0.5);
1082 z[i] = h * (i + 0.5);
1083 }
1084
1085 PetscCall(MatGetOwnershipRange(user->Qblock, &istart, &iend));
1086 nx = user->mx;
1087 ny = user->mx;
1088 nz = user->mx;
1089 for (i = istart; i < iend; i++) {
1090 xri = xr[i];
1091 im = 0;
1092 xim = x[im];
1093 while (xri > xim && im < nx) {
1094 im = im + 1;
1095 xim = x[im];
1096 }
1097 indx1 = im - 1;
1098 indx2 = im;
1099 dx1 = xri - x[indx1];
1100 dx2 = x[indx2] - xri;
1101
1102 yri = yr[i];
1103 im = 0;
1104 yim = y[im];
1105 while (yri > yim && im < ny) {
1106 im = im + 1;
1107 yim = y[im];
1108 }
1109 indy1 = im - 1;
1110 indy2 = im;
1111 dy1 = yri - y[indy1];
1112 dy2 = y[indy2] - yri;
1113
1114 zri = zr[i];
1115 im = 0;
1116 zim = z[im];
1117 while (zri > zim && im < nz) {
1118 im = im + 1;
1119 zim = z[im];
1120 }
1121 indz1 = im - 1;
1122 indz2 = im;
1123 dz1 = zri - z[indz1];
1124 dz2 = z[indz2] - zri;
1125
1126 Dx = x[indx2] - x[indx1];
1127 Dy = y[indy2] - y[indy1];
1128 Dz = z[indz2] - z[indz1];
1129
1130 j = indx1 + indy1 * nx + indz1 * nx * ny;
1131 v = (1 - dx1 / Dx) * (1 - dy1 / Dy) * (1 - dz1 / Dz);
1132 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES));
1133
1134 j = indx1 + indy1 * nx + indz2 * nx * ny;
1135 v = (1 - dx1 / Dx) * (1 - dy1 / Dy) * (1 - dz2 / Dz);
1136 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES));
1137
1138 j = indx1 + indy2 * nx + indz1 * nx * ny;
1139 v = (1 - dx1 / Dx) * (1 - dy2 / Dy) * (1 - dz1 / Dz);
1140 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES));
1141
1142 j = indx1 + indy2 * nx + indz2 * nx * ny;
1143 v = (1 - dx1 / Dx) * (1 - dy2 / Dy) * (1 - dz2 / Dz);
1144 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES));
1145
1146 j = indx2 + indy1 * nx + indz1 * nx * ny;
1147 v = (1 - dx2 / Dx) * (1 - dy1 / Dy) * (1 - dz1 / Dz);
1148 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES));
1149
1150 j = indx2 + indy1 * nx + indz2 * nx * ny;
1151 v = (1 - dx2 / Dx) * (1 - dy1 / Dy) * (1 - dz2 / Dz);
1152 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES));
1153
1154 j = indx2 + indy2 * nx + indz1 * nx * ny;
1155 v = (1 - dx2 / Dx) * (1 - dy2 / Dy) * (1 - dz1 / Dz);
1156 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES));
1157
1158 j = indx2 + indy2 * nx + indz2 * nx * ny;
1159 v = (1 - dx2 / Dx) * (1 - dy2 / Dy) * (1 - dz2 / Dz);
1160 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES));
1161 }
1162 PetscCall(MatAssemblyBegin(user->Qblock, MAT_FINAL_ASSEMBLY));
1163 PetscCall(MatAssemblyEnd(user->Qblock, MAT_FINAL_ASSEMBLY));
1164
1165 PetscCall(MatTranspose(user->Qblock, MAT_INITIAL_MATRIX, &user->QblockT));
1166 PetscCall(MatTranspose(user->L, MAT_INITIAL_MATRIX, &user->LT));
1167
1168 /* Add noise to the measurement data */
1169 PetscCall(VecSet(user->ywork, 1.0));
1170 PetscCall(VecAYPX(user->ywork, user->noise, user->ytrue));
1171 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt));
1172 for (j = 0; j < user->ns; j++) {
1173 i = user->sample_times[j];
1174 PetscCall(MatMult(user->Qblock, user->yiwork[i], user->di[j]));
1175 }
1176 PetscCall(Gather_i(user->d, user->di, user->di_scatter, user->ns));
1177
1178 /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1179 PetscCall(KSPSetFromOptions(user->solver));
1180 PetscCall(PetscFree(x));
1181 PetscCall(PetscFree(y));
1182 PetscCall(PetscFree(z));
1183 PetscFunctionReturn(PETSC_SUCCESS);
1184 }
1185
ParabolicDestroy(AppCtx * user)1186 PetscErrorCode ParabolicDestroy(AppCtx *user)
1187 {
1188 PetscInt i;
1189
1190 PetscFunctionBegin;
1191 PetscCall(MatDestroy(&user->Qblock));
1192 PetscCall(MatDestroy(&user->QblockT));
1193 PetscCall(MatDestroy(&user->Div));
1194 PetscCall(MatDestroy(&user->Divwork));
1195 PetscCall(MatDestroy(&user->Grad));
1196 PetscCall(MatDestroy(&user->Av));
1197 PetscCall(MatDestroy(&user->Avwork));
1198 PetscCall(MatDestroy(&user->AvT));
1199 PetscCall(MatDestroy(&user->DSG));
1200 PetscCall(MatDestroy(&user->L));
1201 PetscCall(MatDestroy(&user->LT));
1202 PetscCall(KSPDestroy(&user->solver));
1203 PetscCall(MatDestroy(&user->Js));
1204 PetscCall(MatDestroy(&user->Jd));
1205 PetscCall(MatDestroy(&user->JsInv));
1206 PetscCall(MatDestroy(&user->JsBlock));
1207 PetscCall(MatDestroy(&user->JsBlockPrec));
1208 PetscCall(VecDestroy(&user->u));
1209 PetscCall(VecDestroy(&user->uwork));
1210 PetscCall(VecDestroy(&user->utrue));
1211 PetscCall(VecDestroy(&user->y));
1212 PetscCall(VecDestroy(&user->ywork));
1213 PetscCall(VecDestroy(&user->ytrue));
1214 PetscCall(VecDestroyVecs(user->nt, &user->yi));
1215 PetscCall(VecDestroyVecs(user->nt, &user->yiwork));
1216 PetscCall(VecDestroyVecs(user->ns, &user->di));
1217 PetscCall(PetscFree(user->yi));
1218 PetscCall(PetscFree(user->yiwork));
1219 PetscCall(PetscFree(user->di));
1220 PetscCall(VecDestroy(&user->c));
1221 PetscCall(VecDestroy(&user->cwork));
1222 PetscCall(VecDestroy(&user->ur));
1223 PetscCall(VecDestroy(&user->q));
1224 PetscCall(VecDestroy(&user->d));
1225 PetscCall(VecDestroy(&user->dwork));
1226 PetscCall(VecDestroy(&user->lwork));
1227 PetscCall(VecDestroy(&user->S));
1228 PetscCall(VecDestroy(&user->Swork));
1229 PetscCall(VecDestroy(&user->Av_u));
1230 PetscCall(VecDestroy(&user->Twork));
1231 PetscCall(VecDestroy(&user->Rwork));
1232 PetscCall(VecDestroy(&user->js_diag));
1233 PetscCall(ISDestroy(&user->s_is));
1234 PetscCall(ISDestroy(&user->d_is));
1235 PetscCall(VecScatterDestroy(&user->state_scatter));
1236 PetscCall(VecScatterDestroy(&user->design_scatter));
1237 for (i = 0; i < user->nt; i++) PetscCall(VecScatterDestroy(&user->yi_scatter[i]));
1238 for (i = 0; i < user->ns; i++) PetscCall(VecScatterDestroy(&user->di_scatter[i]));
1239 PetscCall(PetscFree(user->yi_scatter));
1240 PetscCall(PetscFree(user->di_scatter));
1241 PetscCall(PetscFree(user->sample_times));
1242 PetscFunctionReturn(PETSC_SUCCESS);
1243 }
1244
ParabolicMonitor(Tao tao,void * ptr)1245 PetscErrorCode ParabolicMonitor(Tao tao, void *ptr)
1246 {
1247 Vec X;
1248 PetscReal unorm, ynorm;
1249 AppCtx *user = (AppCtx *)ptr;
1250
1251 PetscFunctionBegin;
1252 PetscCall(TaoGetSolution(tao, &X));
1253 PetscCall(Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
1254 PetscCall(VecAXPY(user->ywork, -1.0, user->ytrue));
1255 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue));
1256 PetscCall(VecNorm(user->uwork, NORM_2, &unorm));
1257 PetscCall(VecNorm(user->ywork, NORM_2, &ynorm));
1258 PetscCall(PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n", (double)unorm, (double)ynorm));
1259 PetscFunctionReturn(PETSC_SUCCESS);
1260 }
1261
1262 /*TEST
1263
1264 build:
1265 requires: !complex
1266
1267 test:
1268 args: -tao_monitor_constraint_norm -tao_type lcl -ns 1 -tao_gatol 1.e-4 -ksp_max_it 30
1269 requires: !single
1270
1271 TEST*/
1272