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