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