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