Actual source code: hyperbolic.c
petsc-3.6.4 2016-04-12
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: TaoSetInitialVector();
8: Routines: TaoSetObjectiveRoutine();
9: Routines: TaoSetGradientRoutine();
10: Routines: TaoSetConstraintsRoutine();
11: Routines: TaoSetJacobianStateRoutine();
12: Routines: TaoSetJacobianDesignRoutine();
13: Routines: TaoSetStateDesignIS();
14: Routines: TaoSetFromOptions();
15: Routines: TaoSetHistory(); TaoGetHistory();
16: Routines: TaoSolve();
17: Routines: TaoGetConvergedReason(); TaoDestroy();
18: Processors: 1
19: T*/
21: typedef struct {
22: PetscInt n; /* Number of variables */
23: PetscInt m; /* Number of constraints */
24: PetscInt mx; /* grid points in each direction */
25: PetscInt nt; /* Number of time steps */
26: PetscInt ndata; /* Number of data points per sample */
27: IS s_is;
28: IS d_is;
29: VecScatter state_scatter;
30: VecScatter design_scatter;
31: VecScatter *uxi_scatter,*uyi_scatter,*ux_scatter,*uy_scatter,*ui_scatter;
32: VecScatter *yi_scatter;
34: Mat Js,Jd,JsBlockPrec,JsInv,JsBlock;
35: PetscBool jformed,c_formed;
37: PetscReal alpha; /* Regularization parameter */
38: PetscReal gamma;
39: PetscReal ht; /* Time step */
40: PetscReal T; /* Final time */
41: Mat Q,QT;
42: Mat L,LT;
43: Mat Div,Divwork,Divxy[2];
44: Mat Grad,Gradxy[2];
45: Mat M;
46: Mat *C,*Cwork;
47: /* Mat Hs,Hd,Hsd; */
48: Vec q;
49: Vec ur; /* reference */
51: Vec d;
52: Vec dwork;
54: Vec y; /* state variables */
55: Vec ywork;
56: Vec ytrue;
57: Vec *yi,*yiwork,*ziwork;
58: Vec *uxi,*uyi,*uxiwork,*uyiwork,*ui,*uiwork;
60: Vec u; /* design variables */
61: Vec uwork,vwork;
62: Vec utrue;
64: Vec js_diag;
66: Vec c; /* constraint vector */
67: Vec cwork;
69: Vec lwork;
71: KSP solver;
72: PC prec;
73: PetscInt block_index;
75: PetscInt ksp_its;
76: PetscInt ksp_its_initial;
77: } AppCtx;
80: PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*);
81: PetscErrorCode FormGradient(Tao, Vec, Vec, void*);
82: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*);
83: PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*);
84: PetscErrorCode FormJacobianDesign(Tao, Vec, Mat,void*);
85: PetscErrorCode FormConstraints(Tao, Vec, Vec, void*);
86: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*);
87: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
88: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
89: PetscErrorCode HyperbolicInitialize(AppCtx *user);
90: PetscErrorCode HyperbolicDestroy(AppCtx *user);
91: PetscErrorCode HyperbolicMonitor(Tao, void*);
93: PetscErrorCode StateMatMult(Mat,Vec,Vec);
94: PetscErrorCode StateMatBlockMult(Mat,Vec,Vec);
95: PetscErrorCode StateMatBlockMultTranspose(Mat,Vec,Vec);
96: PetscErrorCode StateMatMultTranspose(Mat,Vec,Vec);
97: PetscErrorCode StateMatGetDiagonal(Mat,Vec);
98: PetscErrorCode StateMatDuplicate(Mat,MatDuplicateOption,Mat*);
99: PetscErrorCode StateMatInvMult(Mat,Vec,Vec);
100: PetscErrorCode StateMatInvTransposeMult(Mat,Vec,Vec);
101: PetscErrorCode StateMatBlockPrecMult(PC,Vec,Vec);
103: PetscErrorCode DesignMatMult(Mat,Vec,Vec);
104: PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);
106: PetscErrorCode Scatter_yi(Vec,Vec*,VecScatter*,PetscInt); /* y to y1,y2,...,y_nt */
107: PetscErrorCode Gather_yi(Vec,Vec*,VecScatter*,PetscInt);
108: PetscErrorCode Scatter_uxi_uyi(Vec,Vec*,VecScatter*,Vec*,VecScatter*,PetscInt); /* u to ux_1,uy_1,ux_2,uy_2,...,u */
109: PetscErrorCode Gather_uxi_uyi(Vec,Vec*,VecScatter*,Vec*,VecScatter*,PetscInt);
111: static char help[]="";
115: int main(int argc, char **argv)
116: {
117: PetscErrorCode ierr;
118: Vec x,x0;
119: Tao tao;
120: TaoConvergedReason reason;
121: AppCtx user;
122: IS is_allstate,is_alldesign;
123: PetscInt lo,hi,hi2,lo2,ksp_old;
124: PetscInt ntests = 1;
125: PetscInt i;
126: int stages[1];
128: PetscInitialize(&argc, &argv, (char*)0,help);
130: user.mx = 32;
131: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);
132: PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);
133: user.nt = 16;
134: PetscOptionsInt("-nt","Number of time steps","",user.nt,&user.nt,NULL);
135: user.ndata = 64;
136: PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL);
137: user.alpha = 10.0;
138: PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);
139: user.T = 1.0/32.0;
140: PetscOptionsReal("-Tfinal","Final time","",user.T,&user.T,NULL);
142: user.m = user.mx*user.mx*user.nt; /* number of constraints */
143: user.n = user.mx*user.mx*3*user.nt; /* number of variables */
144: user.ht = user.T/user.nt; /* Time step */
145: user.gamma = user.T*user.ht / (user.mx*user.mx);
147: VecCreate(PETSC_COMM_WORLD,&user.u);
148: VecCreate(PETSC_COMM_WORLD,&user.y);
149: VecCreate(PETSC_COMM_WORLD,&user.c);
150: VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m);
151: VecSetSizes(user.y,PETSC_DECIDE,user.m);
152: VecSetSizes(user.c,PETSC_DECIDE,user.m);
153: VecSetFromOptions(user.u);
154: VecSetFromOptions(user.y);
155: VecSetFromOptions(user.c);
157: /* Create scatters for reduced spaces.
158: If the state vector y and design vector u are partitioned as
159: [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
160: then the solution vector x is organized as
161: [y_1; u_1; y_2; u_2; ...; y_np; u_np].
162: The index sets user.s_is and user.d_is correspond to the indices of the
163: state and design variables owned by the current processor.
164: */
165: VecCreate(PETSC_COMM_WORLD,&x);
167: VecGetOwnershipRange(user.y,&lo,&hi);
168: VecGetOwnershipRange(user.u,&lo2,&hi2);
170: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);
171: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is);
173: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);
174: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is);
176: VecSetSizes(x,hi-lo+hi2-lo2,user.n);
177: VecSetFromOptions(x);
179: VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter);
180: VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter);
181: ISDestroy(&is_alldesign);
182: ISDestroy(&is_allstate);
184: /* Create TAO solver and set desired solution method */
185: TaoCreate(PETSC_COMM_WORLD,&tao);
186: TaoSetType(tao,TAOLCL);
188: /* Set up initial vectors and matrices */
189: HyperbolicInitialize(&user);
191: Gather(x,user.y,user.state_scatter,user.u,user.design_scatter);
192: VecDuplicate(x,&x0);
193: VecCopy(x,x0);
195: /* Set solution vector with an initial guess */
196: TaoSetInitialVector(tao,x);
197: TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);
198: TaoSetGradientRoutine(tao, FormGradient, (void *)&user);
199: TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);
201: TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, (void *)&user);
203: TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);
205: TaoSetFromOptions(tao);
206: TaoSetStateDesignIS(tao,user.s_is,user.d_is);
208: /* SOLVE THE APPLICATION */
209: PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL);
210: PetscOptionsEnd();
211: PetscLogStageRegister("Trials",&stages[0]);
212: PetscLogStagePush(stages[0]);
213: user.ksp_its_initial = user.ksp_its;
214: ksp_old = user.ksp_its;
215: for (i=0; i<ntests; i++){
216: TaoSolve(tao);
217: PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old);
218: VecCopy(x0,x);
219: TaoSetInitialVector(tao,x);
220: }
221: PetscLogStagePop();
222: PetscBarrier((PetscObject)x);
223: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
224: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);
225: PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests);
226: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its);
227: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations per trial: ");
228: PetscPrintf(PETSC_COMM_WORLD,"%D\n",(user.ksp_its-user.ksp_its_initial)/ntests);
230: TaoGetConvergedReason(tao,&reason);
232: if (reason < 0) {
233: PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge.\n");
234: } else {
235: PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);
236: }
238: TaoDestroy(&tao);
239: VecDestroy(&x);
240: VecDestroy(&x0);
241: HyperbolicDestroy(&user);
242: PetscFinalize();
243: return 0;
244: }
245: /* ------------------------------------------------------------------- */
248: /*
249: dwork = Qy - d
250: lwork = L*(u-ur).^2
251: f = 1/2 * (dwork.dork + alpha*y.lwork)
252: */
253: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
254: {
256: PetscReal d1=0,d2=0;
257: AppCtx *user = (AppCtx*)ptr;
260: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
261: MatMult(user->Q,user->y,user->dwork);
262: VecAXPY(user->dwork,-1.0,user->d);
263: VecDot(user->dwork,user->dwork,&d1);
265: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
266: VecPointwiseMult(user->uwork,user->uwork,user->uwork);
267: MatMult(user->L,user->uwork,user->lwork);
268: VecDot(user->y,user->lwork,&d2);
269: *f = 0.5 * (d1 + user->alpha*d2);
270: return(0);
271: }
273: /* ------------------------------------------------------------------- */
276: /*
277: state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2
278: design: g_d = alpha*(L'y).*(u-ur)
279: */
280: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
281: {
283: AppCtx *user = (AppCtx*)ptr;
286: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
287: MatMult(user->Q,user->y,user->dwork);
288: VecAXPY(user->dwork,-1.0,user->d);
290: MatMult(user->QT,user->dwork,user->ywork);
292: MatMult(user->LT,user->y,user->uwork);
293: VecWAXPY(user->vwork,-1.0,user->ur,user->u);
294: VecPointwiseMult(user->uwork,user->vwork,user->uwork);
295: VecScale(user->uwork,user->alpha);
297: VecPointwiseMult(user->vwork,user->vwork,user->vwork);
298: MatMult(user->L,user->vwork,user->lwork);
299: VecAXPY(user->ywork,0.5*user->alpha,user->lwork);
301: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
302: return(0);
303: }
307: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
308: {
310: PetscReal d1,d2;
311: AppCtx *user = (AppCtx*)ptr;
314: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
315: MatMult(user->Q,user->y,user->dwork);
316: VecAXPY(user->dwork,-1.0,user->d);
318: MatMult(user->QT,user->dwork,user->ywork);
320: VecDot(user->dwork,user->dwork,&d1);
322: MatMult(user->LT,user->y,user->uwork);
323: VecWAXPY(user->vwork,-1.0,user->ur,user->u);
324: VecPointwiseMult(user->uwork,user->vwork,user->uwork);
325: VecScale(user->uwork,user->alpha);
327: VecPointwiseMult(user->vwork,user->vwork,user->vwork);
328: MatMult(user->L,user->vwork,user->lwork);
329: VecAXPY(user->ywork,0.5*user->alpha,user->lwork);
331: VecDot(user->y,user->lwork,&d2);
333: *f = 0.5 * (d1 + user->alpha*d2);
334: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
335: return(0);
336: }
338: /* ------------------------------------------------------------------- */
341: /* A
342: MatShell object
343: */
344: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
345: {
347: PetscInt i;
348: AppCtx *user = (AppCtx*)ptr;
351: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
352: Scatter_yi(user->u,user->ui,user->ui_scatter,user->nt);
353: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
354: for (i=0; i<user->nt; i++){
355: MatCopy(user->Divxy[0],user->C[i],SAME_NONZERO_PATTERN);
356: MatCopy(user->Divxy[1],user->Cwork[i],SAME_NONZERO_PATTERN);
358: MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
359: MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
360: MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
361: MatScale(user->C[i],user->ht);
362: MatShift(user->C[i],1.0);
363: }
364: return(0);
365: }
367: /* ------------------------------------------------------------------- */
370: /* B */
371: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
372: {
374: AppCtx *user = (AppCtx*)ptr;
377: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
378: return(0);
379: }
383: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
384: {
386: PetscInt i;
387: AppCtx *user;
390: MatShellGetContext(J_shell,(void**)&user);
391: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
392: user->block_index = 0;
393: MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);
395: for (i=1; i<user->nt; i++){
396: user->block_index = i;
397: MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
398: MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
399: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
400: }
401: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
402: return(0);
403: }
407: PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
408: {
410: PetscInt i;
411: AppCtx *user;
414: MatShellGetContext(J_shell,(void**)&user);
415: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
417: for (i=0; i<user->nt-1; i++){
418: user->block_index = i;
419: MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
420: MatMult(user->M,user->yi[i+1],user->ziwork[i+1]);
421: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i+1]);
422: }
424: i = user->nt-1;
425: user->block_index = i;
426: MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
427: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
428: return(0);
429: }
433: PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
434: {
436: PetscInt i;
437: AppCtx *user;
440: MatShellGetContext(J_shell,(void**)&user);
441: i = user->block_index;
442: VecPointwiseMult(user->uxiwork[i],X,user->uxi[i]);
443: VecPointwiseMult(user->uyiwork[i],X,user->uyi[i]);
444: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
445: MatMult(user->Div,user->uiwork[i],Y);
446: VecAYPX(Y,user->ht,X);
447: return(0);
448: }
452: PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y)
453: {
455: PetscInt i;
456: AppCtx *user;
459: MatShellGetContext(J_shell,(void**)&user);
460: i = user->block_index;
461: MatMult(user->Grad,X,user->uiwork[i]);
462: Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
463: VecPointwiseMult(user->uxiwork[i],user->uxi[i],user->uxiwork[i]);
464: VecPointwiseMult(user->uyiwork[i],user->uyi[i],user->uyiwork[i]);
465: VecWAXPY(Y,1.0,user->uxiwork[i],user->uyiwork[i]);
466: VecAYPX(Y,user->ht,X);
467: return(0);
468: }
472: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
473: {
475: PetscInt i;
476: AppCtx *user;
479: MatShellGetContext(J_shell,(void**)&user);
480: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
481: Scatter_uxi_uyi(X,user->uxiwork,user->uxi_scatter,user->uyiwork,user->uyi_scatter,user->nt);
482: for (i=0; i<user->nt; i++){
483: VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
484: VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
485: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
486: MatMult(user->Div,user->uiwork[i],user->ziwork[i]);
487: VecScale(user->ziwork[i],user->ht);
488: }
489: Gather_yi(Y,user->ziwork,user->yi_scatter,user->nt);
490: return(0);
491: }
495: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
496: {
498: PetscInt i;
499: AppCtx *user;
502: MatShellGetContext(J_shell,(void**)&user);
503: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
504: Scatter_yi(X,user->yiwork,user->yi_scatter,user->nt);
505: for (i=0; i<user->nt; i++){
506: MatMult(user->Grad,user->yiwork[i],user->uiwork[i]);
507: Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
508: VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
509: VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
510: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
511: VecScale(user->uiwork[i],user->ht);
512: }
513: Gather_yi(Y,user->uiwork,user->ui_scatter,user->nt);
514: return(0);
515: }
519: PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
520: {
522: PetscInt i;
523: AppCtx *user;
526: PCShellGetContext(PC_shell,(void**)&user);
527: i = user->block_index;
528: if (user->c_formed) {
529: MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
530: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
531: return(0);
532: }
536: PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
537: {
539: PetscInt i;
540: AppCtx *user;
543: PCShellGetContext(PC_shell,(void**)&user);
545: i = user->block_index;
546: if (user->c_formed) {
547: MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
548: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
549: return(0);
550: }
554: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
555: {
557: AppCtx *user;
558: PetscInt its,i;
561: MatShellGetContext(J_shell,(void**)&user);
563: if (Y == user->ytrue) {
564: /* First solve is done using true solution to set up problem */
565: KSPSetTolerances(user->solver,1e-4,1e-20,PETSC_DEFAULT,PETSC_DEFAULT);
566: } else {
567: KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
568: }
569: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
570: Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
571: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
573: user->block_index = 0;
574: KSPSolve(user->solver,user->yi[0],user->yiwork[0]);
576: KSPGetIterationNumber(user->solver,&its);
577: user->ksp_its = user->ksp_its + its;
578: for (i=1; i<user->nt; i++){
579: MatMult(user->M,user->yiwork[i-1],user->ziwork[i-1]);
580: VecAXPY(user->yi[i],1.0,user->ziwork[i-1]);
581: user->block_index = i;
582: KSPSolve(user->solver,user->yi[i],user->yiwork[i]);
584: KSPGetIterationNumber(user->solver,&its);
585: user->ksp_its = user->ksp_its + its;
586: }
587: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
588: return(0);
589: }
593: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
594: {
596: AppCtx *user;
597: PetscInt its,i;
600: MatShellGetContext(J_shell,(void**)&user);
602: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
603: Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
604: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
606: i = user->nt - 1;
607: user->block_index = i;
608: KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);
610: KSPGetIterationNumber(user->solver,&its);
611: user->ksp_its = user->ksp_its + its;
613: for (i=user->nt-2; i>=0; i--){
614: MatMult(user->M,user->yiwork[i+1],user->ziwork[i+1]);
615: VecAXPY(user->yi[i],1.0,user->ziwork[i+1]);
616: user->block_index = i;
617: KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);
619: KSPGetIterationNumber(user->solver,&its);
620: user->ksp_its = user->ksp_its + its;
621: }
622: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
623: return(0);
624: }
628: PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
629: {
631: AppCtx *user;
634: MatShellGetContext(J_shell,(void**)&user);
636: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_shell);
637: MatShellSetOperation(*new_shell,MATOP_MULT,(void(*)(void))StateMatMult);
638: MatShellSetOperation(*new_shell,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
639: MatShellSetOperation(*new_shell,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
640: MatShellSetOperation(*new_shell,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
641: return(0);
642: }
646: PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
647: {
649: AppCtx *user;
652: MatShellGetContext(J_shell,(void**)&user);
653: VecCopy(user->js_diag,X);
654: return(0);
655: }
659: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
660: {
661: /* con = Ay - q, A = [C(u1) 0 0 ... 0;
662: -M C(u2) 0 ... 0;
663: 0 -M C(u3) ... 0;
664: ... ;
665: 0 ... -M C(u_nt)]
666: C(u) = eye + ht*Div*[diag(u1); diag(u2)] */
668: PetscInt i;
669: AppCtx *user = (AppCtx*)ptr;
672: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
673: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
674: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
676: user->block_index = 0;
677: MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);
679: for (i=1; i<user->nt; i++){
680: user->block_index = i;
681: MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
682: MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
683: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
684: }
686: Gather_yi(C,user->yiwork,user->yi_scatter,user->nt);
687: VecAXPY(C,-1.0,user->q);
689: return(0);
690: }
695: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
696: {
700: VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
701: VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
702: VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
703: VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
704: return(0);
705: }
709: PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
710: {
712: PetscInt i;
715: for (i=0; i<nt; i++){
716: VecScatterBegin(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
717: VecScatterEnd(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
718: VecScatterBegin(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
719: VecScatterEnd(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
720: }
721: return(0);
722: }
726: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
727: {
731: VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
732: VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
733: VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
734: VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
735: return(0);
736: }
740: PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
741: {
743: PetscInt i;
746: for (i=0; i<nt; i++){
747: VecScatterBegin(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
748: VecScatterEnd(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
749: VecScatterBegin(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
750: VecScatterEnd(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
751: }
752: return(0);
753: }
757: PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
758: {
760: PetscInt i;
763: for (i=0; i<nt; i++){
764: VecScatterBegin(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
765: VecScatterEnd(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
766: }
767: return(0);
768: }
772: PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
773: {
775: PetscInt i;
778: for (i=0; i<nt; i++){
779: VecScatterBegin(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
780: VecScatterEnd(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
781: }
782: return(0);
783: }
787: PetscErrorCode HyperbolicInitialize(AppCtx *user)
788: {
790: PetscInt n,i,j,linear_index,istart,iend,iblock,lo,hi;
791: Vec XX,YY,XXwork,YYwork,yi,uxi,ui,bc;
792: PetscReal h,sum;
793: PetscScalar hinv,neg_hinv,quarter=0.25,one=1.0,half_hinv,neg_half_hinv;
794: PetscScalar vx,vy,zero=0.0;
795: IS is_from_y,is_to_yi,is_from_u,is_to_uxi,is_to_uyi;
798: user->jformed = PETSC_FALSE;
799: user->c_formed = PETSC_FALSE;
801: user->ksp_its = 0;
802: user->ksp_its_initial = 0;
804: n = user->mx * user->mx;
806: h = 1.0/user->mx;
807: hinv = user->mx;
808: neg_hinv = -hinv;
809: half_hinv = hinv / 2.0;
810: neg_half_hinv = neg_hinv / 2.0;
812: /* Generate Grad matrix */
813: MatCreate(PETSC_COMM_WORLD,&user->Grad);
814: MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,2*n,n);
815: MatSetFromOptions(user->Grad);
816: MatMPIAIJSetPreallocation(user->Grad,3,NULL,3,NULL);
817: MatSeqAIJSetPreallocation(user->Grad,3,NULL);
818: MatGetOwnershipRange(user->Grad,&istart,&iend);
820: for (i=istart; i<iend; i++){
821: if (i<n){
822: iblock = i / user->mx;
823: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
824: MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
825: j = iblock*user->mx + ((i+1) % user->mx);
826: MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
827: }
828: if (i>=n){
829: j = (i - user->mx) % n;
830: MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
831: j = (j + 2*user->mx) % n;
832: MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
833: }
834: }
836: MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
837: MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);
839: MatCreate(PETSC_COMM_WORLD,&user->Gradxy[0]);
840: MatSetSizes(user->Gradxy[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
841: MatSetFromOptions(user->Gradxy[0]);
842: MatMPIAIJSetPreallocation(user->Gradxy[0],3,NULL,3,NULL);
843: MatSeqAIJSetPreallocation(user->Gradxy[0],3,NULL);
844: MatGetOwnershipRange(user->Gradxy[0],&istart,&iend);
846: for (i=istart; i<iend; i++){
847: iblock = i / user->mx;
848: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
849: MatSetValues(user->Gradxy[0],1,&i,1,&j,&half_hinv,INSERT_VALUES);
850: j = iblock*user->mx + ((i+1) % user->mx);
851: MatSetValues(user->Gradxy[0],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
852: MatSetValues(user->Gradxy[0],1,&i,1,&i,&zero,INSERT_VALUES);
853: }
854: MatAssemblyBegin(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
855: MatAssemblyEnd(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
857: MatCreate(PETSC_COMM_WORLD,&user->Gradxy[1]);
858: MatSetSizes(user->Gradxy[1],PETSC_DECIDE,PETSC_DECIDE,n,n);
859: MatSetFromOptions(user->Gradxy[1]);
860: MatMPIAIJSetPreallocation(user->Gradxy[1],3,NULL,3,NULL);
861: MatSeqAIJSetPreallocation(user->Gradxy[1],3,NULL);
862: MatGetOwnershipRange(user->Gradxy[1],&istart,&iend);
864: for (i=istart; i<iend; i++){
865: j = (i + n - user->mx) % n;
866: MatSetValues(user->Gradxy[1],1,&i,1,&j,&half_hinv,INSERT_VALUES);
867: j = (j + 2*user->mx) % n;
868: MatSetValues(user->Gradxy[1],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
869: MatSetValues(user->Gradxy[1],1,&i,1,&i,&zero,INSERT_VALUES);
870: }
871: MatAssemblyBegin(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
872: MatAssemblyEnd(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
874: /* Generate Div matrix */
875: MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div);
876: MatTranspose(user->Gradxy[0],MAT_INITIAL_MATRIX,&user->Divxy[0]);
877: MatTranspose(user->Gradxy[1],MAT_INITIAL_MATRIX,&user->Divxy[1]);
879: /* Off-diagonal averaging matrix */
880: MatCreate(PETSC_COMM_WORLD,&user->M);
881: MatSetSizes(user->M,PETSC_DECIDE,PETSC_DECIDE,n,n);
882: MatSetFromOptions(user->M);
883: MatMPIAIJSetPreallocation(user->M,4,NULL,4,NULL);
884: MatSeqAIJSetPreallocation(user->M,4,NULL);
885: MatGetOwnershipRange(user->M,&istart,&iend);
887: for (i=istart; i<iend; i++){
888: /* kron(Id,Av) */
889: iblock = i / user->mx;
890: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
891: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
892: j = iblock*user->mx + ((i+1) % user->mx);
893: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
895: /* kron(Av,Id) */
896: j = (i + user->mx) % n;
897: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
898: j = (i + n - user->mx) % n;
899: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
900: }
901: MatAssemblyBegin(user->M,MAT_FINAL_ASSEMBLY);
902: MatAssemblyEnd(user->M,MAT_FINAL_ASSEMBLY);
904: /* Generate 2D grid */
905: VecCreate(PETSC_COMM_WORLD,&XX);
906: VecCreate(PETSC_COMM_WORLD,&user->q);
907: VecSetSizes(XX,PETSC_DECIDE,n);
908: VecSetSizes(user->q,PETSC_DECIDE,n*user->nt);
909: VecSetFromOptions(XX);
910: VecSetFromOptions(user->q);
912: VecDuplicate(XX,&YY);
913: VecDuplicate(XX,&XXwork);
914: VecDuplicate(XX,&YYwork);
915: VecDuplicate(XX,&user->d);
916: VecDuplicate(XX,&user->dwork);
918: VecGetOwnershipRange(XX,&istart,&iend);
919: for (linear_index=istart; linear_index<iend; linear_index++){
920: i = linear_index % user->mx;
921: j = (linear_index-i)/user->mx;
922: vx = h*(i+0.5);
923: vy = h*(j+0.5);
924: VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
925: VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
926: }
928: VecAssemblyBegin(XX);
929: VecAssemblyEnd(XX);
930: VecAssemblyBegin(YY);
931: VecAssemblyEnd(YY);
933: /* Compute final density function yT
934: yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
935: yT = yT / (h^2*sum(yT)) */
936: VecCopy(XX,XXwork);
937: VecCopy(YY,YYwork);
939: VecShift(XXwork,-0.25);
940: VecShift(YYwork,-0.25);
942: VecPointwiseMult(XXwork,XXwork,XXwork);
943: VecPointwiseMult(YYwork,YYwork,YYwork);
945: VecCopy(XXwork,user->dwork);
946: VecAXPY(user->dwork,1.0,YYwork);
947: VecScale(user->dwork,-30.0);
948: VecExp(user->dwork);
949: VecCopy(user->dwork,user->d);
951: VecCopy(XX,XXwork);
952: VecCopy(YY,YYwork);
954: VecShift(XXwork,-0.75);
955: VecShift(YYwork,-0.75);
957: VecPointwiseMult(XXwork,XXwork,XXwork);
958: VecPointwiseMult(YYwork,YYwork,YYwork);
960: VecCopy(XXwork,user->dwork);
961: VecAXPY(user->dwork,1.0,YYwork);
962: VecScale(user->dwork,-30.0);
963: VecExp(user->dwork);
965: VecAXPY(user->d,1.0,user->dwork);
966: VecShift(user->d,1.0);
967: VecSum(user->d,&sum);
968: VecScale(user->d,1.0/(h*h*sum));
970: /* Initial conditions of forward problem */
971: VecDuplicate(XX,&bc);
972: VecCopy(XX,XXwork);
973: VecCopy(YY,YYwork);
975: VecShift(XXwork,-0.5);
976: VecShift(YYwork,-0.5);
978: VecPointwiseMult(XXwork,XXwork,XXwork);
979: VecPointwiseMult(YYwork,YYwork,YYwork);
981: VecWAXPY(bc,1.0,XXwork,YYwork);
982: VecScale(bc,-50.0);
983: VecExp(bc);
984: VecShift(bc,1.0);
985: VecSum(bc,&sum);
986: VecScale(bc,1.0/(h*h*sum));
988: /* Create scatter from y to y_1,y_2,...,y_nt */
989: /* TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
990: PetscMalloc1(user->nt*user->mx*user->mx,&user->yi_scatter);
991: VecCreate(PETSC_COMM_WORLD,&yi);
992: VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx);
993: VecSetFromOptions(yi);
994: VecDuplicateVecs(yi,user->nt,&user->yi);
995: VecDuplicateVecs(yi,user->nt,&user->yiwork);
996: VecDuplicateVecs(yi,user->nt,&user->ziwork);
997: for (i=0; i<user->nt; i++){
998: VecGetOwnershipRange(user->yi[i],&lo,&hi);
999: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_yi);
1000: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+i*user->mx*user->mx,1,&is_from_y);
1001: VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i]);
1002: ISDestroy(&is_to_yi);
1003: ISDestroy(&is_from_y);
1004: }
1006: /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
1007: /* TODO: reorder for better parallelism */
1008: PetscMalloc1(user->nt*user->mx*user->mx,&user->uxi_scatter);
1009: PetscMalloc1(user->nt*user->mx*user->mx,&user->uyi_scatter);
1010: PetscMalloc1(user->nt*user->mx*user->mx,&user->ux_scatter);
1011: PetscMalloc1(user->nt*user->mx*user->mx,&user->uy_scatter);
1012: PetscMalloc1(2*user->nt*user->mx*user->mx,&user->ui_scatter);
1013: VecCreate(PETSC_COMM_WORLD,&uxi);
1014: VecCreate(PETSC_COMM_WORLD,&ui);
1015: VecSetSizes(uxi,PETSC_DECIDE,user->mx*user->mx);
1016: VecSetSizes(ui,PETSC_DECIDE,2*user->mx*user->mx);
1017: VecSetFromOptions(uxi);
1018: VecSetFromOptions(ui);
1019: VecDuplicateVecs(uxi,user->nt,&user->uxi);
1020: VecDuplicateVecs(uxi,user->nt,&user->uyi);
1021: VecDuplicateVecs(uxi,user->nt,&user->uxiwork);
1022: VecDuplicateVecs(uxi,user->nt,&user->uyiwork);
1023: VecDuplicateVecs(ui,user->nt,&user->ui);
1024: VecDuplicateVecs(ui,user->nt,&user->uiwork);
1025: for (i=0; i<user->nt; i++){
1026: VecGetOwnershipRange(user->uxi[i],&lo,&hi);
1027: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1028: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
1029: VecScatterCreate(user->u,is_from_u,user->uxi[i],is_to_uxi,&user->uxi_scatter[i]);
1031: ISDestroy(&is_to_uxi);
1032: ISDestroy(&is_from_u);
1034: VecGetOwnershipRange(user->uyi[i],&lo,&hi);
1035: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
1036: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+(2*i+1)*user->mx*user->mx,1,&is_from_u);
1037: VecScatterCreate(user->u,is_from_u,user->uyi[i],is_to_uyi,&user->uyi_scatter[i]);
1039: ISDestroy(&is_to_uyi);
1040: ISDestroy(&is_from_u);
1042: VecGetOwnershipRange(user->uxi[i],&lo,&hi);
1043: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1044: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_from_u);
1045: VecScatterCreate(user->ui[i],is_from_u,user->uxi[i],is_to_uxi,&user->ux_scatter[i]);
1047: ISDestroy(&is_to_uxi);
1048: ISDestroy(&is_from_u);
1050: VecGetOwnershipRange(user->uyi[i],&lo,&hi);
1051: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
1052: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+user->mx*user->mx,1,&is_from_u);
1053: VecScatterCreate(user->ui[i],is_from_u,user->uyi[i],is_to_uyi,&user->uy_scatter[i]);
1055: ISDestroy(&is_to_uyi);
1056: ISDestroy(&is_from_u);
1058: VecGetOwnershipRange(user->ui[i],&lo,&hi);
1059: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1060: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
1061: VecScatterCreate(user->u,is_from_u,user->ui[i],is_to_uxi,&user->ui_scatter[i]);
1063: ISDestroy(&is_to_uxi);
1064: ISDestroy(&is_from_u);
1065: }
1067: /* RHS of forward problem */
1068: MatMult(user->M,bc,user->yiwork[0]);
1069: for (i=1; i<user->nt; i++){
1070: VecSet(user->yiwork[i],0.0);
1071: }
1072: Gather_yi(user->q,user->yiwork,user->yi_scatter,user->nt);
1074: /* Compute true velocity field utrue */
1075: VecDuplicate(user->u,&user->utrue);
1076: for (i=0; i<user->nt; i++){
1077: VecCopy(YY,user->uxi[i]);
1078: VecScale(user->uxi[i],150.0*i*user->ht);
1079: VecCopy(XX,user->uyi[i]);
1080: VecShift(user->uyi[i],-10.0);
1081: VecScale(user->uyi[i],15.0*i*user->ht);
1082: }
1083: Gather_uxi_uyi(user->utrue,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1085: /* Initial guess and reference model */
1086: VecDuplicate(user->utrue,&user->ur);
1087: for (i=0; i<user->nt; i++){
1088: VecCopy(XX,user->uxi[i]);
1089: VecShift(user->uxi[i],i*user->ht);
1090: VecCopy(YY,user->uyi[i]);
1091: VecShift(user->uyi[i],-i*user->ht);
1092: }
1093: Gather_uxi_uyi(user->ur,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1095: /* Generate regularization matrix L */
1096: MatCreate(PETSC_COMM_WORLD,&user->LT);
1097: MatSetSizes(user->LT,PETSC_DECIDE,PETSC_DECIDE,2*n*user->nt,n*user->nt);
1098: MatSetFromOptions(user->LT);
1099: MatMPIAIJSetPreallocation(user->LT,1,NULL,1,NULL);
1100: MatSeqAIJSetPreallocation(user->LT,1,NULL);
1101: MatGetOwnershipRange(user->LT,&istart,&iend);
1103: for (i=istart; i<iend; i++){
1104: iblock = (i+n) / (2*n);
1105: j = i - iblock*n;
1106: MatSetValues(user->LT,1,&i,1,&j,&user->gamma,INSERT_VALUES);
1107: }
1109: MatAssemblyBegin(user->LT,MAT_FINAL_ASSEMBLY);
1110: MatAssemblyEnd(user->LT,MAT_FINAL_ASSEMBLY);
1112: MatTranspose(user->LT,MAT_INITIAL_MATRIX,&user->L);
1114: /* Build work vectors and matrices */
1115: VecCreate(PETSC_COMM_WORLD,&user->lwork);
1116: VecSetType(user->lwork,VECMPI);
1117: VecSetSizes(user->lwork,PETSC_DECIDE,user->m);
1118: VecSetFromOptions(user->lwork);
1120: MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);
1122: VecDuplicate(user->y,&user->ywork);
1123: VecDuplicate(user->u,&user->uwork);
1124: VecDuplicate(user->u,&user->vwork);
1125: VecDuplicate(user->u,&user->js_diag);
1126: VecDuplicate(user->c,&user->cwork);
1128: /* Create matrix-free shell user->Js for computing A*x */
1129: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->Js);
1130: MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1131: MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
1132: MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
1133: MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
1135: /* Diagonal blocks of user->Js */
1136: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlock);
1137: MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult);
1138: MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockMultTranspose);
1140: /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
1141: D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
1142: This is an SOR preconditioner for user->JsBlock. */
1143: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlockPrec);
1144: MatShellSetOperation(user->JsBlockPrec,MATOP_MULT,(void(*)(void))StateMatBlockPrecMult);
1145: MatShellSetOperation(user->JsBlockPrec,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockPrecMultTranspose);
1147: /* Create a matrix-free shell user->Jd for computing B*x */
1148: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->n-user->m,user,&user->Jd);
1149: MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);
1150: MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);
1152: /* User-defined routines for computing user->Js\x and user->Js^T\x*/
1153: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsInv);
1154: MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateMatInvMult);
1155: MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatInvTransposeMult);
1157: /* Build matrices for SOR preconditioner */
1158: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1159: MatShift(user->Divxy[0],0.0); /* Force C[i] and Divxy[0] to share same nonzero pattern */
1160: MatAXPY(user->Divxy[0],0.0,user->Divxy[1],DIFFERENT_NONZERO_PATTERN);
1161: PetscMalloc1(5*n,&user->C);
1162: PetscMalloc1(2*n,&user->Cwork);
1163: for (i=0; i<user->nt; i++){
1164: MatDuplicate(user->Divxy[0],MAT_COPY_VALUES,&user->C[i]);
1165: MatDuplicate(user->Divxy[1],MAT_COPY_VALUES,&user->Cwork[i]);
1167: MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
1168: MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
1169: MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
1170: MatScale(user->C[i],user->ht);
1171: MatShift(user->C[i],1.0);
1172: }
1174: /* Solver options and tolerances */
1175: KSPCreate(PETSC_COMM_WORLD,&user->solver);
1176: KSPSetType(user->solver,KSPGMRES);
1177: KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec);
1178: KSPSetInitialGuessNonzero(user->solver,PETSC_FALSE); /* TODO: why is true slower? */
1179: KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500);
1180: /* KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500); */
1181: KSPGetPC(user->solver,&user->prec);
1182: PCSetType(user->prec,PCSHELL);
1184: PCShellSetApply(user->prec,StateMatBlockPrecMult);
1185: PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMultTranspose);
1186: PCShellSetContext(user->prec,user);
1188: /* Compute true state function yt given ut */
1189: VecCreate(PETSC_COMM_WORLD,&user->ytrue);
1190: VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt);
1191: VecSetFromOptions(user->ytrue);
1192: user->c_formed = PETSC_TRUE;
1193: VecCopy(user->utrue,user->u); /* Set u=utrue temporarily for StateMatInv */
1194: VecSet(user->ytrue,0.0); /* Initial guess */
1195: StateMatInvMult(user->Js,user->q,user->ytrue);
1196: VecCopy(user->ur,user->u); /* Reset u=ur */
1198: /* Initial guess y0 for state given u0 */
1199: StateMatInvMult(user->Js,user->q,user->y);
1201: /* Data discretization */
1202: MatCreate(PETSC_COMM_WORLD,&user->Q);
1203: MatSetSizes(user->Q,PETSC_DECIDE,PETSC_DECIDE,user->mx*user->mx,user->m);
1204: MatSetFromOptions(user->Q);
1205: MatMPIAIJSetPreallocation(user->Q,0,NULL,1,NULL);
1206: MatSeqAIJSetPreallocation(user->Q,1,NULL);
1208: MatGetOwnershipRange(user->Q,&istart,&iend);
1210: for (i=istart; i<iend; i++){
1211: j = i + user->m - user->mx*user->mx;
1212: MatSetValues(user->Q,1,&i,1,&j,&one,INSERT_VALUES);
1213: }
1215: MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1216: MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);
1218: MatTranspose(user->Q,MAT_INITIAL_MATRIX,&user->QT);
1220: VecDestroy(&XX);
1221: VecDestroy(&YY);
1222: VecDestroy(&XXwork);
1223: VecDestroy(&YYwork);
1224: VecDestroy(&yi);
1225: VecDestroy(&uxi);
1226: VecDestroy(&ui);
1227: VecDestroy(&bc);
1229: /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1230: KSPSetFromOptions(user->solver);
1231: return(0);
1232: }
1236: PetscErrorCode HyperbolicDestroy(AppCtx *user)
1237: {
1239: PetscInt i;
1242: MatDestroy(&user->Q);
1243: MatDestroy(&user->QT);
1244: MatDestroy(&user->Div);
1245: MatDestroy(&user->Divwork);
1246: MatDestroy(&user->Grad);
1247: MatDestroy(&user->L);
1248: MatDestroy(&user->LT);
1249: KSPDestroy(&user->solver);
1250: MatDestroy(&user->Js);
1251: MatDestroy(&user->Jd);
1252: MatDestroy(&user->JsBlockPrec);
1253: MatDestroy(&user->JsInv);
1254: MatDestroy(&user->JsBlock);
1255: MatDestroy(&user->Divxy[0]);
1256: MatDestroy(&user->Divxy[1]);
1257: MatDestroy(&user->Gradxy[0]);
1258: MatDestroy(&user->Gradxy[1]);
1259: MatDestroy(&user->M);
1260: for (i=0; i<user->nt; i++){
1261: MatDestroy(&user->C[i]);
1262: MatDestroy(&user->Cwork[i]);
1263: }
1264: PetscFree(user->C);
1265: PetscFree(user->Cwork);
1266: VecDestroy(&user->u);
1267: VecDestroy(&user->uwork);
1268: VecDestroy(&user->vwork);
1269: VecDestroy(&user->utrue);
1270: VecDestroy(&user->y);
1271: VecDestroy(&user->ywork);
1272: VecDestroy(&user->ytrue);
1273: VecDestroyVecs(user->nt,&user->yi);
1274: VecDestroyVecs(user->nt,&user->yiwork);
1275: VecDestroyVecs(user->nt,&user->ziwork);
1276: VecDestroyVecs(user->nt,&user->uxi);
1277: VecDestroyVecs(user->nt,&user->uyi);
1278: VecDestroyVecs(user->nt,&user->uxiwork);
1279: VecDestroyVecs(user->nt,&user->uyiwork);
1280: VecDestroyVecs(user->nt,&user->ui);
1281: VecDestroyVecs(user->nt,&user->uiwork);
1282: VecDestroy(&user->c);
1283: VecDestroy(&user->cwork);
1284: VecDestroy(&user->ur);
1285: VecDestroy(&user->q);
1286: VecDestroy(&user->d);
1287: VecDestroy(&user->dwork);
1288: VecDestroy(&user->lwork);
1289: VecDestroy(&user->js_diag);
1290: ISDestroy(&user->s_is);
1291: ISDestroy(&user->d_is);
1292: VecScatterDestroy(&user->state_scatter);
1293: VecScatterDestroy(&user->design_scatter);
1294: for (i=0; i<user->nt; i++){
1295: VecScatterDestroy(&user->uxi_scatter[i]);
1296: VecScatterDestroy(&user->uyi_scatter[i]);
1297: VecScatterDestroy(&user->ux_scatter[i]);
1298: VecScatterDestroy(&user->uy_scatter[i]);
1299: VecScatterDestroy(&user->ui_scatter[i]);
1300: VecScatterDestroy(&user->yi_scatter[i]);
1301: }
1302: PetscFree(user->uxi_scatter);
1303: PetscFree(user->uyi_scatter);
1304: PetscFree(user->ux_scatter);
1305: PetscFree(user->uy_scatter);
1306: PetscFree(user->ui_scatter);
1307: PetscFree(user->yi_scatter);
1308: return(0);
1309: }
1313: PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1314: {
1316: Vec X;
1317: PetscReal unorm,ynorm;
1318: AppCtx *user = (AppCtx*)ptr;
1321: TaoGetSolutionVector(tao,&X);
1322: Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1323: VecAXPY(user->ywork,-1.0,user->ytrue);
1324: VecAXPY(user->uwork,-1.0,user->utrue);
1325: VecNorm(user->uwork,NORM_2,&unorm);
1326: VecNorm(user->ywork,NORM_2,&ynorm);
1327: PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1328: return(0);
1329: }