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