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: PetscCheck(user->c_formed, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed");
459: PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y));
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
464: {
465: PetscInt i;
466: AppCtx *user;
468: PetscFunctionBegin;
469: PetscCall(PCShellGetContext(PC_shell, &user));
471: i = user->block_index;
472: PetscCheck(user->c_formed, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed");
473: PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y));
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
478: {
479: AppCtx *user;
480: PetscInt its, i;
482: PetscFunctionBegin;
483: PetscCall(MatShellGetContext(J_shell, &user));
485: if (Y == user->ytrue) {
486: /* First solve is done using true solution to set up problem */
487: PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, PETSC_CURRENT, PETSC_CURRENT));
488: } else {
489: PetscCall(KSPSetTolerances(user->solver, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
490: }
491: PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt));
492: PetscCall(Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt));
493: PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
495: user->block_index = 0;
496: PetscCall(KSPSolve(user->solver, user->yi[0], user->yiwork[0]));
498: PetscCall(KSPGetIterationNumber(user->solver, &its));
499: user->ksp_its = user->ksp_its + its;
500: for (i = 1; i < user->nt; i++) {
501: PetscCall(MatMult(user->M, user->yiwork[i - 1], user->ziwork[i - 1]));
502: PetscCall(VecAXPY(user->yi[i], 1.0, user->ziwork[i - 1]));
503: user->block_index = i;
504: PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i]));
506: PetscCall(KSPGetIterationNumber(user->solver, &its));
507: user->ksp_its = user->ksp_its + its;
508: }
509: PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
514: {
515: AppCtx *user;
516: PetscInt its, i;
518: PetscFunctionBegin;
519: PetscCall(MatShellGetContext(J_shell, &user));
521: PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt));
522: PetscCall(Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt));
523: PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
525: i = user->nt - 1;
526: user->block_index = i;
527: PetscCall(KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i]));
529: PetscCall(KSPGetIterationNumber(user->solver, &its));
530: user->ksp_its = user->ksp_its + its;
532: for (i = user->nt - 2; i >= 0; i--) {
533: PetscCall(MatMult(user->M, user->yiwork[i + 1], user->ziwork[i + 1]));
534: PetscCall(VecAXPY(user->yi[i], 1.0, user->ziwork[i + 1]));
535: user->block_index = i;
536: PetscCall(KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i]));
538: PetscCall(KSPGetIterationNumber(user->solver, &its));
539: user->ksp_its = user->ksp_its + its;
540: }
541: PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
542: PetscFunctionReturn(PETSC_SUCCESS);
543: }
545: PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
546: {
547: AppCtx *user;
549: PetscFunctionBegin;
550: PetscCall(MatShellGetContext(J_shell, &user));
552: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, new_shell));
553: PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT, (void (*)(void))StateMatMult));
554: PetscCall(MatShellSetOperation(*new_shell, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate));
555: PetscCall(MatShellSetOperation(*new_shell, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose));
556: PetscCall(MatShellSetOperation(*new_shell, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal));
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
561: {
562: AppCtx *user;
564: PetscFunctionBegin;
565: PetscCall(MatShellGetContext(J_shell, &user));
566: PetscCall(VecCopy(user->js_diag, X));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
571: {
572: /* con = Ay - q, A = [C(u1) 0 0 ... 0;
573: -M C(u2) 0 ... 0;
574: 0 -M C(u3) ... 0;
575: ... ;
576: 0 ... -M C(u_nt)]
577: C(u) = eye + ht*Div*[diag(u1); diag(u2)] */
578: PetscInt i;
579: AppCtx *user = (AppCtx *)ptr;
581: PetscFunctionBegin;
582: PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
583: PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt));
584: PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
586: user->block_index = 0;
587: PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0]));
589: for (i = 1; i < user->nt; i++) {
590: user->block_index = i;
591: PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i]));
592: PetscCall(MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1]));
593: PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1]));
594: }
596: PetscCall(Gather_yi(C, user->yiwork, user->yi_scatter, user->nt));
597: PetscCall(VecAXPY(C, -1.0, user->q));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
602: {
603: PetscFunctionBegin;
604: PetscCall(VecScatterBegin(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD));
605: PetscCall(VecScatterEnd(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD));
606: PetscCall(VecScatterBegin(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD));
607: PetscCall(VecScatterEnd(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD));
608: PetscFunctionReturn(PETSC_SUCCESS);
609: }
611: PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
612: {
613: PetscInt i;
615: PetscFunctionBegin;
616: for (i = 0; i < nt; i++) {
617: PetscCall(VecScatterBegin(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD));
618: PetscCall(VecScatterEnd(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD));
619: PetscCall(VecScatterBegin(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD));
620: PetscCall(VecScatterEnd(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD));
621: }
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
626: {
627: PetscFunctionBegin;
628: PetscCall(VecScatterBegin(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE));
629: PetscCall(VecScatterEnd(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE));
630: PetscCall(VecScatterBegin(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE));
631: PetscCall(VecScatterEnd(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
636: {
637: PetscInt i;
639: PetscFunctionBegin;
640: for (i = 0; i < nt; i++) {
641: PetscCall(VecScatterBegin(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE));
642: PetscCall(VecScatterEnd(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE));
643: PetscCall(VecScatterBegin(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE));
644: PetscCall(VecScatterEnd(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE));
645: }
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
650: {
651: PetscInt i;
653: PetscFunctionBegin;
654: for (i = 0; i < nt; i++) {
655: PetscCall(VecScatterBegin(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD));
656: PetscCall(VecScatterEnd(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD));
657: }
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
662: {
663: PetscInt i;
665: PetscFunctionBegin;
666: for (i = 0; i < nt; i++) {
667: PetscCall(VecScatterBegin(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE));
668: PetscCall(VecScatterEnd(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE));
669: }
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: PetscErrorCode HyperbolicInitialize(AppCtx *user)
674: {
675: PetscInt n, i, j, linear_index, istart, iend, iblock, lo, hi;
676: Vec XX, YY, XXwork, YYwork, yi, uxi, ui, bc;
677: PetscReal h, sum;
678: PetscScalar hinv, neg_hinv, quarter = 0.25, one = 1.0, half_hinv, neg_half_hinv;
679: PetscScalar vx, vy, zero = 0.0;
680: IS is_from_y, is_to_yi, is_from_u, is_to_uxi, is_to_uyi;
682: PetscFunctionBegin;
683: user->jformed = PETSC_FALSE;
684: user->c_formed = PETSC_FALSE;
686: user->ksp_its = 0;
687: user->ksp_its_initial = 0;
689: n = user->mx * user->mx;
691: h = 1.0 / user->mx;
692: hinv = user->mx;
693: neg_hinv = -hinv;
694: half_hinv = hinv / 2.0;
695: neg_half_hinv = neg_hinv / 2.0;
697: /* Generate Grad matrix */
698: PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Grad));
699: PetscCall(MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, 2 * n, n));
700: PetscCall(MatSetFromOptions(user->Grad));
701: PetscCall(MatMPIAIJSetPreallocation(user->Grad, 3, NULL, 3, NULL));
702: PetscCall(MatSeqAIJSetPreallocation(user->Grad, 3, NULL));
703: PetscCall(MatGetOwnershipRange(user->Grad, &istart, &iend));
705: for (i = istart; i < iend; i++) {
706: if (i < n) {
707: iblock = i / user->mx;
708: j = iblock * user->mx + ((i + user->mx - 1) % user->mx);
709: PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
710: j = iblock * user->mx + ((i + 1) % user->mx);
711: PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
712: }
713: if (i >= n) {
714: j = (i - user->mx) % n;
715: PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
716: j = (j + 2 * user->mx) % n;
717: PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
718: }
719: }
721: PetscCall(MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY));
722: PetscCall(MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY));
724: PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[0]));
725: PetscCall(MatSetSizes(user->Gradxy[0], PETSC_DECIDE, PETSC_DECIDE, n, n));
726: PetscCall(MatSetFromOptions(user->Gradxy[0]));
727: PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[0], 3, NULL, 3, NULL));
728: PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[0], 3, NULL));
729: PetscCall(MatGetOwnershipRange(user->Gradxy[0], &istart, &iend));
731: for (i = istart; i < iend; i++) {
732: iblock = i / user->mx;
733: j = iblock * user->mx + ((i + user->mx - 1) % user->mx);
734: PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
735: j = iblock * user->mx + ((i + 1) % user->mx);
736: PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
737: PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &i, &zero, INSERT_VALUES));
738: }
739: PetscCall(MatAssemblyBegin(user->Gradxy[0], MAT_FINAL_ASSEMBLY));
740: PetscCall(MatAssemblyEnd(user->Gradxy[0], MAT_FINAL_ASSEMBLY));
742: PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[1]));
743: PetscCall(MatSetSizes(user->Gradxy[1], PETSC_DECIDE, PETSC_DECIDE, n, n));
744: PetscCall(MatSetFromOptions(user->Gradxy[1]));
745: PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[1], 3, NULL, 3, NULL));
746: PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[1], 3, NULL));
747: PetscCall(MatGetOwnershipRange(user->Gradxy[1], &istart, &iend));
749: for (i = istart; i < iend; i++) {
750: j = (i + n - user->mx) % n;
751: PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
752: j = (j + 2 * user->mx) % n;
753: PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
754: PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &i, &zero, INSERT_VALUES));
755: }
756: PetscCall(MatAssemblyBegin(user->Gradxy[1], MAT_FINAL_ASSEMBLY));
757: PetscCall(MatAssemblyEnd(user->Gradxy[1], MAT_FINAL_ASSEMBLY));
759: /* Generate Div matrix */
760: PetscCall(MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div));
761: PetscCall(MatTranspose(user->Gradxy[0], MAT_INITIAL_MATRIX, &user->Divxy[0]));
762: PetscCall(MatTranspose(user->Gradxy[1], MAT_INITIAL_MATRIX, &user->Divxy[1]));
764: /* Off-diagonal averaging matrix */
765: PetscCall(MatCreate(PETSC_COMM_WORLD, &user->M));
766: PetscCall(MatSetSizes(user->M, PETSC_DECIDE, PETSC_DECIDE, n, n));
767: PetscCall(MatSetFromOptions(user->M));
768: PetscCall(MatMPIAIJSetPreallocation(user->M, 4, NULL, 4, NULL));
769: PetscCall(MatSeqAIJSetPreallocation(user->M, 4, NULL));
770: PetscCall(MatGetOwnershipRange(user->M, &istart, &iend));
772: for (i = istart; i < iend; i++) {
773: /* kron(Id,Av) */
774: iblock = i / user->mx;
775: j = iblock * user->mx + ((i + user->mx - 1) % user->mx);
776: PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
777: j = iblock * user->mx + ((i + 1) % user->mx);
778: PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
780: /* kron(Av,Id) */
781: j = (i + user->mx) % n;
782: PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
783: j = (i + n - user->mx) % n;
784: PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
785: }
786: PetscCall(MatAssemblyBegin(user->M, MAT_FINAL_ASSEMBLY));
787: PetscCall(MatAssemblyEnd(user->M, MAT_FINAL_ASSEMBLY));
789: /* Generate 2D grid */
790: PetscCall(VecCreate(PETSC_COMM_WORLD, &XX));
791: PetscCall(VecCreate(PETSC_COMM_WORLD, &user->q));
792: PetscCall(VecSetSizes(XX, PETSC_DECIDE, n));
793: PetscCall(VecSetSizes(user->q, PETSC_DECIDE, n * user->nt));
794: PetscCall(VecSetFromOptions(XX));
795: PetscCall(VecSetFromOptions(user->q));
797: PetscCall(VecDuplicate(XX, &YY));
798: PetscCall(VecDuplicate(XX, &XXwork));
799: PetscCall(VecDuplicate(XX, &YYwork));
800: PetscCall(VecDuplicate(XX, &user->d));
801: PetscCall(VecDuplicate(XX, &user->dwork));
803: PetscCall(VecGetOwnershipRange(XX, &istart, &iend));
804: for (linear_index = istart; linear_index < iend; linear_index++) {
805: i = linear_index % user->mx;
806: j = (linear_index - i) / user->mx;
807: vx = h * (i + 0.5);
808: vy = h * (j + 0.5);
809: PetscCall(VecSetValues(XX, 1, &linear_index, &vx, INSERT_VALUES));
810: PetscCall(VecSetValues(YY, 1, &linear_index, &vy, INSERT_VALUES));
811: }
813: PetscCall(VecAssemblyBegin(XX));
814: PetscCall(VecAssemblyEnd(XX));
815: PetscCall(VecAssemblyBegin(YY));
816: PetscCall(VecAssemblyEnd(YY));
818: /* Compute final density function yT
819: yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
820: yT = yT / (h^2*sum(yT)) */
821: PetscCall(VecCopy(XX, XXwork));
822: PetscCall(VecCopy(YY, YYwork));
824: PetscCall(VecShift(XXwork, -0.25));
825: PetscCall(VecShift(YYwork, -0.25));
827: PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
828: PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));
830: PetscCall(VecCopy(XXwork, user->dwork));
831: PetscCall(VecAXPY(user->dwork, 1.0, YYwork));
832: PetscCall(VecScale(user->dwork, -30.0));
833: PetscCall(VecExp(user->dwork));
834: PetscCall(VecCopy(user->dwork, user->d));
836: PetscCall(VecCopy(XX, XXwork));
837: PetscCall(VecCopy(YY, YYwork));
839: PetscCall(VecShift(XXwork, -0.75));
840: PetscCall(VecShift(YYwork, -0.75));
842: PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
843: PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));
845: PetscCall(VecCopy(XXwork, user->dwork));
846: PetscCall(VecAXPY(user->dwork, 1.0, YYwork));
847: PetscCall(VecScale(user->dwork, -30.0));
848: PetscCall(VecExp(user->dwork));
850: PetscCall(VecAXPY(user->d, 1.0, user->dwork));
851: PetscCall(VecShift(user->d, 1.0));
852: PetscCall(VecSum(user->d, &sum));
853: PetscCall(VecScale(user->d, 1.0 / (h * h * sum)));
855: /* Initial conditions of forward problem */
856: PetscCall(VecDuplicate(XX, &bc));
857: PetscCall(VecCopy(XX, XXwork));
858: PetscCall(VecCopy(YY, YYwork));
860: PetscCall(VecShift(XXwork, -0.5));
861: PetscCall(VecShift(YYwork, -0.5));
863: PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
864: PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));
866: PetscCall(VecWAXPY(bc, 1.0, XXwork, YYwork));
867: PetscCall(VecScale(bc, -50.0));
868: PetscCall(VecExp(bc));
869: PetscCall(VecShift(bc, 1.0));
870: PetscCall(VecSum(bc, &sum));
871: PetscCall(VecScale(bc, 1.0 / (h * h * sum)));
873: /* Create scatter from y to y_1,y_2,...,y_nt */
874: /* TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
875: PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->yi_scatter));
876: PetscCall(VecCreate(PETSC_COMM_WORLD, &yi));
877: PetscCall(VecSetSizes(yi, PETSC_DECIDE, user->mx * user->mx));
878: PetscCall(VecSetFromOptions(yi));
879: PetscCall(VecDuplicateVecs(yi, user->nt, &user->yi));
880: PetscCall(VecDuplicateVecs(yi, user->nt, &user->yiwork));
881: PetscCall(VecDuplicateVecs(yi, user->nt, &user->ziwork));
882: for (i = 0; i < user->nt; i++) {
883: PetscCall(VecGetOwnershipRange(user->yi[i], &lo, &hi));
884: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_yi));
885: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + i * user->mx * user->mx, 1, &is_from_y));
886: PetscCall(VecScatterCreate(user->y, is_from_y, user->yi[i], is_to_yi, &user->yi_scatter[i]));
887: PetscCall(ISDestroy(&is_to_yi));
888: PetscCall(ISDestroy(&is_from_y));
889: }
891: /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
892: /* TODO: reorder for better parallelism */
893: PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uxi_scatter));
894: PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uyi_scatter));
895: PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->ux_scatter));
896: PetscCall(PetscMalloc1(user->nt * user->mx * user->mx, &user->uy_scatter));
897: PetscCall(PetscMalloc1(2 * user->nt * user->mx * user->mx, &user->ui_scatter));
898: PetscCall(VecCreate(PETSC_COMM_WORLD, &uxi));
899: PetscCall(VecCreate(PETSC_COMM_WORLD, &ui));
900: PetscCall(VecSetSizes(uxi, PETSC_DECIDE, user->mx * user->mx));
901: PetscCall(VecSetSizes(ui, PETSC_DECIDE, 2 * user->mx * user->mx));
902: PetscCall(VecSetFromOptions(uxi));
903: PetscCall(VecSetFromOptions(ui));
904: PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uxi));
905: PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uyi));
906: PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uxiwork));
907: PetscCall(VecDuplicateVecs(uxi, user->nt, &user->uyiwork));
908: PetscCall(VecDuplicateVecs(ui, user->nt, &user->ui));
909: PetscCall(VecDuplicateVecs(ui, user->nt, &user->uiwork));
910: for (i = 0; i < user->nt; i++) {
911: PetscCall(VecGetOwnershipRange(user->uxi[i], &lo, &hi));
912: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi));
913: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u));
914: PetscCall(VecScatterCreate(user->u, is_from_u, user->uxi[i], is_to_uxi, &user->uxi_scatter[i]));
916: PetscCall(ISDestroy(&is_to_uxi));
917: PetscCall(ISDestroy(&is_from_u));
919: PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi));
920: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi));
921: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + (2 * i + 1) * user->mx * user->mx, 1, &is_from_u));
922: PetscCall(VecScatterCreate(user->u, is_from_u, user->uyi[i], is_to_uyi, &user->uyi_scatter[i]));
924: PetscCall(ISDestroy(&is_to_uyi));
925: PetscCall(ISDestroy(&is_from_u));
927: PetscCall(VecGetOwnershipRange(user->uxi[i], &lo, &hi));
928: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi));
929: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_from_u));
930: PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uxi[i], is_to_uxi, &user->ux_scatter[i]));
932: PetscCall(ISDestroy(&is_to_uxi));
933: PetscCall(ISDestroy(&is_from_u));
935: PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi));
936: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi));
937: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + user->mx * user->mx, 1, &is_from_u));
938: PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uyi[i], is_to_uyi, &user->uy_scatter[i]));
940: PetscCall(ISDestroy(&is_to_uyi));
941: PetscCall(ISDestroy(&is_from_u));
943: PetscCall(VecGetOwnershipRange(user->ui[i], &lo, &hi));
944: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi));
945: PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u));
946: PetscCall(VecScatterCreate(user->u, is_from_u, user->ui[i], is_to_uxi, &user->ui_scatter[i]));
948: PetscCall(ISDestroy(&is_to_uxi));
949: PetscCall(ISDestroy(&is_from_u));
950: }
952: /* RHS of forward problem */
953: PetscCall(MatMult(user->M, bc, user->yiwork[0]));
954: for (i = 1; i < user->nt; i++) PetscCall(VecSet(user->yiwork[i], 0.0));
955: PetscCall(Gather_yi(user->q, user->yiwork, user->yi_scatter, user->nt));
957: /* Compute true velocity field utrue */
958: PetscCall(VecDuplicate(user->u, &user->utrue));
959: for (i = 0; i < user->nt; i++) {
960: PetscCall(VecCopy(YY, user->uxi[i]));
961: PetscCall(VecScale(user->uxi[i], 150.0 * i * user->ht));
962: PetscCall(VecCopy(XX, user->uyi[i]));
963: PetscCall(VecShift(user->uyi[i], -10.0));
964: PetscCall(VecScale(user->uyi[i], 15.0 * i * user->ht));
965: }
966: PetscCall(Gather_uxi_uyi(user->utrue, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
968: /* Initial guess and reference model */
969: PetscCall(VecDuplicate(user->utrue, &user->ur));
970: for (i = 0; i < user->nt; i++) {
971: PetscCall(VecCopy(XX, user->uxi[i]));
972: PetscCall(VecShift(user->uxi[i], i * user->ht));
973: PetscCall(VecCopy(YY, user->uyi[i]));
974: PetscCall(VecShift(user->uyi[i], -i * user->ht));
975: }
976: PetscCall(Gather_uxi_uyi(user->ur, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
978: /* Generate regularization matrix L */
979: PetscCall(MatCreate(PETSC_COMM_WORLD, &user->LT));
980: PetscCall(MatSetSizes(user->LT, PETSC_DECIDE, PETSC_DECIDE, 2 * n * user->nt, n * user->nt));
981: PetscCall(MatSetFromOptions(user->LT));
982: PetscCall(MatMPIAIJSetPreallocation(user->LT, 1, NULL, 1, NULL));
983: PetscCall(MatSeqAIJSetPreallocation(user->LT, 1, NULL));
984: PetscCall(MatGetOwnershipRange(user->LT, &istart, &iend));
986: for (i = istart; i < iend; i++) {
987: iblock = (i + n) / (2 * n);
988: j = i - iblock * n;
989: PetscCall(MatSetValues(user->LT, 1, &i, 1, &j, &user->gamma, INSERT_VALUES));
990: }
992: PetscCall(MatAssemblyBegin(user->LT, MAT_FINAL_ASSEMBLY));
993: PetscCall(MatAssemblyEnd(user->LT, MAT_FINAL_ASSEMBLY));
995: PetscCall(MatTranspose(user->LT, MAT_INITIAL_MATRIX, &user->L));
997: /* Build work vectors and matrices */
998: PetscCall(VecCreate(PETSC_COMM_WORLD, &user->lwork));
999: PetscCall(VecSetType(user->lwork, VECMPI));
1000: PetscCall(VecSetSizes(user->lwork, PETSC_DECIDE, user->m));
1001: PetscCall(VecSetFromOptions(user->lwork));
1003: PetscCall(MatDuplicate(user->Div, MAT_SHARE_NONZERO_PATTERN, &user->Divwork));
1005: PetscCall(VecDuplicate(user->y, &user->ywork));
1006: PetscCall(VecDuplicate(user->u, &user->uwork));
1007: PetscCall(VecDuplicate(user->u, &user->vwork));
1008: PetscCall(VecDuplicate(user->u, &user->js_diag));
1009: PetscCall(VecDuplicate(user->c, &user->cwork));
1011: /* Create matrix-free shell user->Js for computing A*x */
1012: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->Js));
1013: PetscCall(MatShellSetOperation(user->Js, MATOP_MULT, (void (*)(void))StateMatMult));
1014: PetscCall(MatShellSetOperation(user->Js, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate));
1015: PetscCall(MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose));
1016: PetscCall(MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal));
1018: /* Diagonal blocks of user->Js */
1019: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlock));
1020: PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT, (void (*)(void))StateMatBlockMult));
1021: PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockMultTranspose));
1023: /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
1024: D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
1025: This is an SOR preconditioner for user->JsBlock. */
1026: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlockPrec));
1027: PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT, (void (*)(void))StateMatBlockPrecMult));
1028: PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockPrecMultTranspose));
1030: /* Create a matrix-free shell user->Jd for computing B*x */
1031: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->n - user->m, user, &user->Jd));
1032: PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT, (void (*)(void))DesignMatMult));
1033: PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT_TRANSPOSE, (void (*)(void))DesignMatMultTranspose));
1035: /* User-defined routines for computing user->Js\x and user->Js^T\x*/
1036: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsInv));
1037: PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT, (void (*)(void))StateMatInvMult));
1038: PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatInvTransposeMult));
1040: /* Build matrices for SOR preconditioner */
1041: PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
1042: PetscCall(PetscMalloc1(5 * n, &user->C));
1043: PetscCall(PetscMalloc1(2 * n, &user->Cwork));
1044: for (i = 0; i < user->nt; i++) {
1045: PetscCall(MatDuplicate(user->Divxy[0], MAT_COPY_VALUES, &user->C[i]));
1046: PetscCall(MatDuplicate(user->Divxy[1], MAT_COPY_VALUES, &user->Cwork[i]));
1048: PetscCall(MatDiagonalScale(user->C[i], NULL, user->uxi[i]));
1049: PetscCall(MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i]));
1050: PetscCall(MatAXPY(user->C[i], 1.0, user->Cwork[i], DIFFERENT_NONZERO_PATTERN));
1051: PetscCall(MatScale(user->C[i], user->ht));
1052: PetscCall(MatShift(user->C[i], 1.0));
1053: }
1055: /* Solver options and tolerances */
1056: PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->solver));
1057: PetscCall(KSPSetType(user->solver, KSPGMRES));
1058: PetscCall(KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec));
1059: PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500));
1060: /* PetscCall(KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500)); */
1061: PetscCall(KSPGetPC(user->solver, &user->prec));
1062: PetscCall(PCSetType(user->prec, PCSHELL));
1064: PetscCall(PCShellSetApply(user->prec, StateMatBlockPrecMult));
1065: PetscCall(PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMultTranspose));
1066: PetscCall(PCShellSetContext(user->prec, user));
1068: /* Compute true state function yt given ut */
1069: PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ytrue));
1070: PetscCall(VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt));
1071: PetscCall(VecSetFromOptions(user->ytrue));
1072: user->c_formed = PETSC_TRUE;
1073: PetscCall(VecCopy(user->utrue, user->u)); /* Set u=utrue temporarily for StateMatInv */
1074: PetscCall(VecSet(user->ytrue, 0.0)); /* Initial guess */
1075: PetscCall(StateMatInvMult(user->Js, user->q, user->ytrue));
1076: PetscCall(VecCopy(user->ur, user->u)); /* Reset u=ur */
1078: /* Initial guess y0 for state given u0 */
1079: PetscCall(StateMatInvMult(user->Js, user->q, user->y));
1081: /* Data discretization */
1082: PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Q));
1083: PetscCall(MatSetSizes(user->Q, PETSC_DECIDE, PETSC_DECIDE, user->mx * user->mx, user->m));
1084: PetscCall(MatSetFromOptions(user->Q));
1085: PetscCall(MatMPIAIJSetPreallocation(user->Q, 0, NULL, 1, NULL));
1086: PetscCall(MatSeqAIJSetPreallocation(user->Q, 1, NULL));
1088: PetscCall(MatGetOwnershipRange(user->Q, &istart, &iend));
1090: for (i = istart; i < iend; i++) {
1091: j = i + user->m - user->mx * user->mx;
1092: PetscCall(MatSetValues(user->Q, 1, &i, 1, &j, &one, INSERT_VALUES));
1093: }
1095: PetscCall(MatAssemblyBegin(user->Q, MAT_FINAL_ASSEMBLY));
1096: PetscCall(MatAssemblyEnd(user->Q, MAT_FINAL_ASSEMBLY));
1098: PetscCall(MatTranspose(user->Q, MAT_INITIAL_MATRIX, &user->QT));
1100: PetscCall(VecDestroy(&XX));
1101: PetscCall(VecDestroy(&YY));
1102: PetscCall(VecDestroy(&XXwork));
1103: PetscCall(VecDestroy(&YYwork));
1104: PetscCall(VecDestroy(&yi));
1105: PetscCall(VecDestroy(&uxi));
1106: PetscCall(VecDestroy(&ui));
1107: PetscCall(VecDestroy(&bc));
1109: /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1110: PetscCall(KSPSetFromOptions(user->solver));
1111: PetscFunctionReturn(PETSC_SUCCESS);
1112: }
1114: PetscErrorCode HyperbolicDestroy(AppCtx *user)
1115: {
1116: PetscInt i;
1118: PetscFunctionBegin;
1119: PetscCall(MatDestroy(&user->Q));
1120: PetscCall(MatDestroy(&user->QT));
1121: PetscCall(MatDestroy(&user->Div));
1122: PetscCall(MatDestroy(&user->Divwork));
1123: PetscCall(MatDestroy(&user->Grad));
1124: PetscCall(MatDestroy(&user->L));
1125: PetscCall(MatDestroy(&user->LT));
1126: PetscCall(KSPDestroy(&user->solver));
1127: PetscCall(MatDestroy(&user->Js));
1128: PetscCall(MatDestroy(&user->Jd));
1129: PetscCall(MatDestroy(&user->JsBlockPrec));
1130: PetscCall(MatDestroy(&user->JsInv));
1131: PetscCall(MatDestroy(&user->JsBlock));
1132: PetscCall(MatDestroy(&user->Divxy[0]));
1133: PetscCall(MatDestroy(&user->Divxy[1]));
1134: PetscCall(MatDestroy(&user->Gradxy[0]));
1135: PetscCall(MatDestroy(&user->Gradxy[1]));
1136: PetscCall(MatDestroy(&user->M));
1137: for (i = 0; i < user->nt; i++) {
1138: PetscCall(MatDestroy(&user->C[i]));
1139: PetscCall(MatDestroy(&user->Cwork[i]));
1140: }
1141: PetscCall(PetscFree(user->C));
1142: PetscCall(PetscFree(user->Cwork));
1143: PetscCall(VecDestroy(&user->u));
1144: PetscCall(VecDestroy(&user->uwork));
1145: PetscCall(VecDestroy(&user->vwork));
1146: PetscCall(VecDestroy(&user->utrue));
1147: PetscCall(VecDestroy(&user->y));
1148: PetscCall(VecDestroy(&user->ywork));
1149: PetscCall(VecDestroy(&user->ytrue));
1150: PetscCall(VecDestroyVecs(user->nt, &user->yi));
1151: PetscCall(VecDestroyVecs(user->nt, &user->yiwork));
1152: PetscCall(VecDestroyVecs(user->nt, &user->ziwork));
1153: PetscCall(VecDestroyVecs(user->nt, &user->uxi));
1154: PetscCall(VecDestroyVecs(user->nt, &user->uyi));
1155: PetscCall(VecDestroyVecs(user->nt, &user->uxiwork));
1156: PetscCall(VecDestroyVecs(user->nt, &user->uyiwork));
1157: PetscCall(VecDestroyVecs(user->nt, &user->ui));
1158: PetscCall(VecDestroyVecs(user->nt, &user->uiwork));
1159: PetscCall(VecDestroy(&user->c));
1160: PetscCall(VecDestroy(&user->cwork));
1161: PetscCall(VecDestroy(&user->ur));
1162: PetscCall(VecDestroy(&user->q));
1163: PetscCall(VecDestroy(&user->d));
1164: PetscCall(VecDestroy(&user->dwork));
1165: PetscCall(VecDestroy(&user->lwork));
1166: PetscCall(VecDestroy(&user->js_diag));
1167: PetscCall(ISDestroy(&user->s_is));
1168: PetscCall(ISDestroy(&user->d_is));
1169: PetscCall(VecScatterDestroy(&user->state_scatter));
1170: PetscCall(VecScatterDestroy(&user->design_scatter));
1171: for (i = 0; i < user->nt; i++) {
1172: PetscCall(VecScatterDestroy(&user->uxi_scatter[i]));
1173: PetscCall(VecScatterDestroy(&user->uyi_scatter[i]));
1174: PetscCall(VecScatterDestroy(&user->ux_scatter[i]));
1175: PetscCall(VecScatterDestroy(&user->uy_scatter[i]));
1176: PetscCall(VecScatterDestroy(&user->ui_scatter[i]));
1177: PetscCall(VecScatterDestroy(&user->yi_scatter[i]));
1178: }
1179: PetscCall(PetscFree(user->uxi_scatter));
1180: PetscCall(PetscFree(user->uyi_scatter));
1181: PetscCall(PetscFree(user->ux_scatter));
1182: PetscCall(PetscFree(user->uy_scatter));
1183: PetscCall(PetscFree(user->ui_scatter));
1184: PetscCall(PetscFree(user->yi_scatter));
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1189: {
1190: Vec X;
1191: PetscReal unorm, ynorm;
1192: AppCtx *user = (AppCtx *)ptr;
1194: PetscFunctionBegin;
1195: PetscCall(TaoGetSolution(tao, &X));
1196: PetscCall(Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
1197: PetscCall(VecAXPY(user->ywork, -1.0, user->ytrue));
1198: PetscCall(VecAXPY(user->uwork, -1.0, user->utrue));
1199: PetscCall(VecNorm(user->uwork, NORM_2, &unorm));
1200: PetscCall(VecNorm(user->ywork, NORM_2, &ynorm));
1201: PetscCall(PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n", (double)unorm, (double)ynorm));
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: /*TEST
1207: build:
1208: requires: !complex
1210: test:
1211: requires: !single
1212: args: -tao_monitor_constraint_norm -tao_max_funcs 10 -tao_type lcl -tao_gatol 1.e-5
1214: test:
1215: suffix: guess_pod
1216: requires: !single
1217: args: -tao_monitor_constraint_norm -tao_max_funcs 10 -tao_type lcl -ksp_guess_type pod -tao_gatol 1.e-5
1219: TEST*/