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