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