Actual source code: hyperbolic.c
petsc-3.5.4 2015-05-23
1: #include <petsctao.h>
2: #include "../src/tao/pde_constrained/impls/lcl/lcl.h"
4: /*T
5: Concepts: TAO^Solving a system of nonlinear equations, nonlinear least squares
6: Routines: TaoCreate();
7: Routines: TaoSetType();
8: Routines: TaoSetInitialVector();
9: Routines: TaoSetObjectiveRoutine();
10: Routines: TaoSetGradientRoutine();
11: Routines: TaoSetConstraintsRoutine();
12: Routines: TaoSetJacobianStateRoutine();
13: Routines: TaoSetJacobianDesignRoutine();
14: Routines: TaoSetStateDesignIS();
15: Routines: TaoSetFromOptions();
16: Routines: TaoSetHistory(); TaoGetHistory();
17: Routines: TaoSolve();
18: Routines: TaoGetConvergedReason(); TaoDestroy();
19: Processors: 1
20: T*/
22: typedef struct {
23: PetscInt n; /* Number of variables */
24: PetscInt m; /* Number of constraints */
25: PetscInt mx; /* grid points in each direction */
26: PetscInt nt; /* Number of time steps */
27: PetscInt ndata; /* Number of data points per sample */
28: IS s_is;
29: IS d_is;
30: VecScatter state_scatter;
31: VecScatter design_scatter;
32: VecScatter *uxi_scatter,*uyi_scatter,*ux_scatter,*uy_scatter,*ui_scatter;
33: VecScatter *yi_scatter;
35: Mat Js,Jd,JsBlockPrec,JsInv,JsBlock;
36: PetscBool jformed,c_formed;
38: PetscReal alpha; /* Regularization parameter */
39: PetscReal gamma;
40: PetscReal ht; /* Time step */
41: PetscReal T; /* Final time */
42: Mat Q,QT;
43: Mat L,LT;
44: Mat Div,Divwork,Divxy[2];
45: Mat Grad,Gradxy[2];
46: Mat M;
47: Mat *C,*Cwork;
48: /* Mat Hs,Hd,Hsd; */
49: Vec q;
50: Vec ur; /* reference */
52: Vec d;
53: Vec dwork;
55: Vec y; /* state variables */
56: Vec ywork;
57: Vec ytrue;
58: Vec *yi,*yiwork,*ziwork;
59: Vec *uxi,*uyi,*uxiwork,*uyiwork,*ui,*uiwork;
61: Vec u; /* design variables */
62: Vec uwork,vwork;
63: Vec utrue;
65: Vec js_diag;
67: Vec c; /* constraint vector */
68: Vec cwork;
70: Vec lwork;
72: KSP solver;
73: PC prec;
74: PetscInt block_index;
76: TAO_LCL *lcl;
77: PetscInt ksp_its;
78: PetscInt ksp_its_initial;
79: } AppCtx;
82: PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*);
83: PetscErrorCode FormGradient(Tao, Vec, Vec, void*);
84: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*);
85: PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*);
86: PetscErrorCode FormJacobianDesign(Tao, Vec, Mat,void*);
87: PetscErrorCode FormConstraints(Tao, Vec, Vec, void*);
88: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*);
89: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
90: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
91: PetscErrorCode HyperbolicInitialize(AppCtx *user);
92: PetscErrorCode HyperbolicDestroy(AppCtx *user);
93: PetscErrorCode HyperbolicMonitor(Tao, void*);
95: PetscErrorCode StateMatMult(Mat,Vec,Vec);
96: PetscErrorCode StateMatBlockMult(Mat,Vec,Vec);
97: PetscErrorCode StateMatBlockMultTranspose(Mat,Vec,Vec);
98: PetscErrorCode StateMatMultTranspose(Mat,Vec,Vec);
99: PetscErrorCode StateMatGetDiagonal(Mat,Vec);
100: PetscErrorCode StateMatDuplicate(Mat,MatDuplicateOption,Mat*);
101: PetscErrorCode StateMatInvMult(Mat,Vec,Vec);
102: PetscErrorCode StateMatInvTransposeMult(Mat,Vec,Vec);
103: PetscErrorCode StateMatBlockPrecMult(PC,Vec,Vec);
105: PetscErrorCode DesignMatMult(Mat,Vec,Vec);
106: PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);
108: PetscErrorCode Scatter_yi(Vec,Vec*,VecScatter*,PetscInt); /* y to y1,y2,...,y_nt */
109: PetscErrorCode Gather_yi(Vec,Vec*,VecScatter*,PetscInt);
110: PetscErrorCode Scatter_uxi_uyi(Vec,Vec*,VecScatter*,Vec*,VecScatter*,PetscInt); /* u to ux_1,uy_1,ux_2,uy_2,...,u */
111: PetscErrorCode Gather_uxi_uyi(Vec,Vec*,VecScatter*,Vec*,VecScatter*,PetscInt);
113: static char help[]="";
117: int main(int argc, char **argv)
118: {
119: PetscErrorCode ierr;
120: Vec x,x0;
121: Tao tao;
122: TaoConvergedReason reason;
123: AppCtx user;
124: IS is_allstate,is_alldesign;
125: PetscInt lo,hi,hi2,lo2,ksp_old;
126: PetscBool flag;
127: PetscInt ntests = 1;
128: PetscInt i;
129: int stages[1];
131: PetscInitialize(&argc, &argv, (char*)0,help);
133: user.mx = 32;
134: PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,&flag);
135: user.nt = 16;
136: PetscOptionsInt("-nt","Number of time steps","",user.nt,&user.nt,&flag);
137: user.ndata = 64;
138: PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,&flag);
139: user.alpha = 10.0;
140: PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,&flag);
141: user.T = 1.0/32.0;
142: PetscOptionsReal("-Tfinal","Final time","",user.T,&user.T,&flag);
144: user.m = user.mx*user.mx*user.nt; /* number of constraints */
145: user.n = user.mx*user.mx*3*user.nt; /* number of variables */
146: user.ht = user.T/user.nt; /* Time step */
147: user.gamma = user.T*user.ht / (user.mx*user.mx);
149: VecCreate(PETSC_COMM_WORLD,&user.u);
150: VecCreate(PETSC_COMM_WORLD,&user.y);
151: VecCreate(PETSC_COMM_WORLD,&user.c);
152: VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m);
153: VecSetSizes(user.y,PETSC_DECIDE,user.m);
154: VecSetSizes(user.c,PETSC_DECIDE,user.m);
155: VecSetFromOptions(user.u);
156: VecSetFromOptions(user.y);
157: VecSetFromOptions(user.c);
159: /* Create scatters for reduced spaces.
160: If the state vector y and design vector u are partitioned as
161: [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
162: then the solution vector x is organized as
163: [y_1; u_1; y_2; u_2; ...; y_np; u_np].
164: The index sets user.s_is and user.d_is correspond to the indices of the
165: state and design variables owned by the current processor.
166: */
167: VecCreate(PETSC_COMM_WORLD,&x);
169: VecGetOwnershipRange(user.y,&lo,&hi);
170: VecGetOwnershipRange(user.u,&lo2,&hi2);
172: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);
173: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is);
175: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);
176: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is);
178: VecSetSizes(x,hi-lo+hi2-lo2,user.n);
179: VecSetFromOptions(x);
181: VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter);
182: VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter);
183: ISDestroy(&is_alldesign);
184: ISDestroy(&is_allstate);
186: /* Create TAO solver and set desired solution method */
187: TaoCreate(PETSC_COMM_WORLD,&tao);
188: TaoSetType(tao,TAOLCL);
189: user.lcl = (TAO_LCL*)(tao->data);
191: /* Set up initial vectors and matrices */
192: HyperbolicInitialize(&user);
194: Gather(x,user.y,user.state_scatter,user.u,user.design_scatter);
195: VecDuplicate(x,&x0);
196: VecCopy(x,x0);
198: /* Set solution vector with an initial guess */
199: TaoSetInitialVector(tao,x);
200: TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);
201: TaoSetGradientRoutine(tao, FormGradient, (void *)&user);
202: TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);
204: TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, (void *)&user);
206: TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);
208: TaoSetFromOptions(tao);
209: TaoSetStateDesignIS(tao,user.s_is,user.d_is);
211: /* SOLVE THE APPLICATION */
212: PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,&flag);
213: PetscLogStageRegister("Trials",&stages[0]);
214: PetscLogStagePush(stages[0]);
215: user.ksp_its_initial = user.ksp_its;
216: ksp_old = user.ksp_its;
217: for (i=0; i<ntests; i++){
218: TaoSolve(tao);
219: PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old);
220: VecCopy(x0,x);
221: TaoSetInitialVector(tao,x);
222: }
223: PetscLogStagePop();
224: PetscBarrier((PetscObject)x);
225: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
226: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);
227: PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests);
228: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its);
229: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations per trial: ");
230: PetscPrintf(PETSC_COMM_WORLD,"%D\n",(user.ksp_its-user.ksp_its_initial)/ntests);
232: TaoGetConvergedReason(tao,&reason);
234: if (reason < 0) {
235: PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge.\n");
236: } else {
237: PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);
238: }
240: TaoDestroy(&tao);
241: VecDestroy(&x);
242: VecDestroy(&x0);
243: HyperbolicDestroy(&user);
244: PetscFinalize();
245: return 0;
246: }
247: /* ------------------------------------------------------------------- */
250: /*
251: dwork = Qy - d
252: lwork = L*(u-ur).^2
253: f = 1/2 * (dwork.dork + alpha*y.lwork)
254: */
255: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
256: {
258: PetscReal d1=0,d2=0;
259: AppCtx *user = (AppCtx*)ptr;
262: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
263: MatMult(user->Q,user->y,user->dwork);
264: VecAXPY(user->dwork,-1.0,user->d);
265: VecDot(user->dwork,user->dwork,&d1);
267: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
268: VecPointwiseMult(user->uwork,user->uwork,user->uwork);
269: MatMult(user->L,user->uwork,user->lwork);
270: VecDot(user->y,user->lwork,&d2);
271: *f = 0.5 * (d1 + user->alpha*d2);
272: return(0);
273: }
275: /* ------------------------------------------------------------------- */
278: /*
279: state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2
280: design: g_d = alpha*(L'y).*(u-ur)
281: */
282: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
283: {
285: AppCtx *user = (AppCtx*)ptr;
288: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
289: MatMult(user->Q,user->y,user->dwork);
290: VecAXPY(user->dwork,-1.0,user->d);
292: MatMult(user->QT,user->dwork,user->ywork);
294: MatMult(user->LT,user->y,user->uwork);
295: VecWAXPY(user->vwork,-1.0,user->ur,user->u);
296: VecPointwiseMult(user->uwork,user->vwork,user->uwork);
297: VecScale(user->uwork,user->alpha);
299: VecPointwiseMult(user->vwork,user->vwork,user->vwork);
300: MatMult(user->L,user->vwork,user->lwork);
301: VecAXPY(user->ywork,0.5*user->alpha,user->lwork);
303: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
304: return(0);
305: }
309: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
310: {
312: PetscReal d1,d2;
313: AppCtx *user = (AppCtx*)ptr;
316: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
317: MatMult(user->Q,user->y,user->dwork);
318: VecAXPY(user->dwork,-1.0,user->d);
320: MatMult(user->QT,user->dwork,user->ywork);
322: VecDot(user->dwork,user->dwork,&d1);
324: MatMult(user->LT,user->y,user->uwork);
325: VecWAXPY(user->vwork,-1.0,user->ur,user->u);
326: VecPointwiseMult(user->uwork,user->vwork,user->uwork);
327: VecScale(user->uwork,user->alpha);
329: VecPointwiseMult(user->vwork,user->vwork,user->vwork);
330: MatMult(user->L,user->vwork,user->lwork);
331: VecAXPY(user->ywork,0.5*user->alpha,user->lwork);
333: VecDot(user->y,user->lwork,&d2);
335: *f = 0.5 * (d1 + user->alpha*d2);
336: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
337: return(0);
338: }
340: /* ------------------------------------------------------------------- */
343: /* A
344: MatShell object
345: */
346: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
347: {
349: PetscInt i;
350: AppCtx *user = (AppCtx*)ptr;
353: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
354: Scatter_yi(user->u,user->ui,user->ui_scatter,user->nt);
355: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
356: for (i=0; i<user->nt; i++){
357: MatCopy(user->Divxy[0],user->C[i],SAME_NONZERO_PATTERN);
358: MatCopy(user->Divxy[1],user->Cwork[i],SAME_NONZERO_PATTERN);
360: MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
361: MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
362: MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
363: MatScale(user->C[i],user->ht);
364: MatShift(user->C[i],1.0);
365: }
366: return(0);
367: }
369: /* ------------------------------------------------------------------- */
372: /* B */
373: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
374: {
376: AppCtx *user = (AppCtx*)ptr;
379: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
380: return(0);
381: }
385: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
386: {
388: PetscInt i;
389: AppCtx *user;
392: MatShellGetContext(J_shell,(void**)&user);
393: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
394: user->block_index = 0;
395: MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);
397: for (i=1; i<user->nt; i++){
398: user->block_index = i;
399: MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
400: MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
401: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
402: }
403: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
404: return(0);
405: }
409: PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
410: {
412: PetscInt i;
413: AppCtx *user;
416: MatShellGetContext(J_shell,(void**)&user);
417: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
419: for (i=0; i<user->nt-1; i++){
420: user->block_index = i;
421: MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
422: MatMult(user->M,user->yi[i+1],user->ziwork[i+1]);
423: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i+1]);
424: }
426: i = user->nt-1;
427: user->block_index = i;
428: MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
429: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
430: return(0);
431: }
435: PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
436: {
438: PetscInt i;
439: AppCtx *user;
442: MatShellGetContext(J_shell,(void**)&user);
443: i = user->block_index;
444: VecPointwiseMult(user->uxiwork[i],X,user->uxi[i]);
445: VecPointwiseMult(user->uyiwork[i],X,user->uyi[i]);
446: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
447: MatMult(user->Div,user->uiwork[i],Y);
448: VecAYPX(Y,user->ht,X);
449: return(0);
450: }
454: PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y)
455: {
457: PetscInt i;
458: AppCtx *user;
461: MatShellGetContext(J_shell,(void**)&user);
462: i = user->block_index;
463: MatMult(user->Grad,X,user->uiwork[i]);
464: Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
465: VecPointwiseMult(user->uxiwork[i],user->uxi[i],user->uxiwork[i]);
466: VecPointwiseMult(user->uyiwork[i],user->uyi[i],user->uyiwork[i]);
467: VecWAXPY(Y,1.0,user->uxiwork[i],user->uyiwork[i]);
468: VecAYPX(Y,user->ht,X);
469: return(0);
470: }
474: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
475: {
477: PetscInt i;
478: AppCtx *user;
481: MatShellGetContext(J_shell,(void**)&user);
482: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
483: Scatter_uxi_uyi(X,user->uxiwork,user->uxi_scatter,user->uyiwork,user->uyi_scatter,user->nt);
484: for (i=0; i<user->nt; i++){
485: VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
486: VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
487: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
488: MatMult(user->Div,user->uiwork[i],user->ziwork[i]);
489: VecScale(user->ziwork[i],user->ht);
490: }
491: Gather_yi(Y,user->ziwork,user->yi_scatter,user->nt);
492: return(0);
493: }
497: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
498: {
500: PetscInt i;
501: AppCtx *user;
504: MatShellGetContext(J_shell,(void**)&user);
505: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
506: Scatter_yi(X,user->yiwork,user->yi_scatter,user->nt);
507: for (i=0; i<user->nt; i++){
508: MatMult(user->Grad,user->yiwork[i],user->uiwork[i]);
509: Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
510: VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
511: VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
512: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
513: VecScale(user->uiwork[i],user->ht);
514: }
515: Gather_yi(Y,user->uiwork,user->ui_scatter,user->nt);
516: return(0);
517: }
521: PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
522: {
524: PetscInt i;
525: AppCtx *user;
528: PCShellGetContext(PC_shell,(void**)&user);
529: i = user->block_index;
530: if (user->c_formed) {
531: MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
532: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
533: return(0);
534: }
538: PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
539: {
541: PetscInt i;
542: AppCtx *user;
545: PCShellGetContext(PC_shell,(void**)&user);
547: i = user->block_index;
548: if (user->c_formed) {
549: MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
550: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
551: return(0);
552: }
556: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
557: {
559: AppCtx *user;
560: PetscReal tau;
561: PetscInt its,i;
564: MatShellGetContext(J_shell,(void**)&user);
566: if (Y == user->ytrue) {
567: KSPSetTolerances(user->solver,1e-4,1e-20,PETSC_DEFAULT,PETSC_DEFAULT);
568: } else if (user->lcl) {
569: tau = user->lcl->tau[user->lcl->solve_type];
570: KSPSetTolerances(user->solver,tau,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
571: }
572: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
573: Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
574: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
576: user->block_index = 0;
577: KSPSolve(user->solver,user->yi[0],user->yiwork[0]);
579: KSPGetIterationNumber(user->solver,&its);
580: user->ksp_its = user->ksp_its + its;
581: for (i=1; i<user->nt; i++){
582: MatMult(user->M,user->yiwork[i-1],user->ziwork[i-1]);
583: VecAXPY(user->yi[i],1.0,user->ziwork[i-1]);
584: user->block_index = i;
585: KSPSolve(user->solver,user->yi[i],user->yiwork[i]);
587: KSPGetIterationNumber(user->solver,&its);
588: user->ksp_its = user->ksp_its + its;
589: }
590: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
591: return(0);
592: }
596: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
597: {
599: AppCtx *user;
600: PetscReal tau;
601: PetscInt its,i;
604: MatShellGetContext(J_shell,(void**)&user);
606: if (user->lcl) {
607: tau = user->lcl->tau[user->lcl->solve_type];
608: KSPSetTolerances(user->solver,tau,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
609: }
610: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
611: Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
612: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
614: i = user->nt - 1;
615: user->block_index = i;
616: KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);
618: KSPGetIterationNumber(user->solver,&its);
619: user->ksp_its = user->ksp_its + its;
621: for (i=user->nt-2; i>=0; i--){
622: MatMult(user->M,user->yiwork[i+1],user->ziwork[i+1]);
623: VecAXPY(user->yi[i],1.0,user->ziwork[i+1]);
624: user->block_index = i;
625: KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);
627: KSPGetIterationNumber(user->solver,&its);
628: user->ksp_its = user->ksp_its + its;
629: }
630: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
631: return(0);
632: }
636: PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
637: {
639: AppCtx *user;
642: MatShellGetContext(J_shell,(void**)&user);
644: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_shell);
645: MatShellSetOperation(*new_shell,MATOP_MULT,(void(*)(void))StateMatMult);
646: MatShellSetOperation(*new_shell,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
647: MatShellSetOperation(*new_shell,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
648: MatShellSetOperation(*new_shell,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
649: return(0);
650: }
654: PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
655: {
657: AppCtx *user;
660: MatShellGetContext(J_shell,(void**)&user);
661: VecCopy(user->js_diag,X);
662: return(0);
663: }
667: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
668: {
669: /* con = Ay - q, A = [C(u1) 0 0 ... 0;
670: -M C(u2) 0 ... 0;
671: 0 -M C(u3) ... 0;
672: ... ;
673: 0 ... -M C(u_nt)]
674: C(u) = eye + ht*Div*[diag(u1); diag(u2)] */
676: PetscInt i;
677: AppCtx *user = (AppCtx*)ptr;
680: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
681: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
682: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
684: user->block_index = 0;
685: MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);
687: for (i=1; i<user->nt; i++){
688: user->block_index = i;
689: MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
690: MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
691: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
692: }
694: Gather_yi(C,user->yiwork,user->yi_scatter,user->nt);
695: VecAXPY(C,-1.0,user->q);
697: return(0);
698: }
703: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
704: {
708: VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
709: VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
710: VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
711: VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
712: return(0);
713: }
717: PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
718: {
720: PetscInt i;
723: for (i=0; i<nt; i++){
724: VecScatterBegin(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
725: VecScatterEnd(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
726: VecScatterBegin(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
727: VecScatterEnd(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
728: }
729: return(0);
730: }
734: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
735: {
739: VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
740: VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
741: VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
742: VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
743: return(0);
744: }
748: PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
749: {
751: PetscInt i;
754: for (i=0; i<nt; i++){
755: VecScatterBegin(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
756: VecScatterEnd(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
757: VecScatterBegin(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
758: VecScatterEnd(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
759: }
760: return(0);
761: }
765: PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
766: {
768: PetscInt i;
771: for (i=0; i<nt; i++){
772: VecScatterBegin(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
773: VecScatterEnd(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
774: }
775: return(0);
776: }
780: PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
781: {
783: PetscInt i;
786: for (i=0; i<nt; i++){
787: VecScatterBegin(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
788: VecScatterEnd(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
789: }
790: return(0);
791: }
795: PetscErrorCode HyperbolicInitialize(AppCtx *user)
796: {
798: PetscInt n,i,j,linear_index,istart,iend,iblock,lo,hi;
799: Vec XX,YY,XXwork,YYwork,yi,uxi,ui,bc;
800: PetscReal h,sum;
801: PetscScalar hinv,neg_hinv,quarter=0.25,one=1.0,half_hinv,neg_half_hinv;
802: PetscScalar vx,vy,zero=0.0;
803: IS is_from_y,is_to_yi,is_from_u,is_to_uxi,is_to_uyi;
806: user->jformed = PETSC_FALSE;
807: user->c_formed = PETSC_FALSE;
809: user->ksp_its = 0;
810: user->ksp_its_initial = 0;
812: n = user->mx * user->mx;
814: h = 1.0/user->mx;
815: hinv = user->mx;
816: neg_hinv = -hinv;
817: half_hinv = hinv / 2.0;
818: neg_half_hinv = neg_hinv / 2.0;
820: /* Generate Grad matrix */
821: MatCreate(PETSC_COMM_WORLD,&user->Grad);
822: MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,2*n,n);
823: MatSetFromOptions(user->Grad);
824: MatMPIAIJSetPreallocation(user->Grad,3,NULL,3,NULL);
825: MatSeqAIJSetPreallocation(user->Grad,3,NULL);
826: MatGetOwnershipRange(user->Grad,&istart,&iend);
828: for (i=istart; i<iend; i++){
829: if (i<n){
830: iblock = i / user->mx;
831: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
832: MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
833: j = iblock*user->mx + ((i+1) % user->mx);
834: MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
835: }
836: if (i>=n){
837: j = (i - user->mx) % n;
838: MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
839: j = (j + 2*user->mx) % n;
840: MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
841: }
842: }
844: MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
845: MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);
847: MatCreate(PETSC_COMM_WORLD,&user->Gradxy[0]);
848: MatSetSizes(user->Gradxy[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
849: MatSetFromOptions(user->Gradxy[0]);
850: MatMPIAIJSetPreallocation(user->Gradxy[0],3,NULL,3,NULL);
851: MatSeqAIJSetPreallocation(user->Gradxy[0],3,NULL);
852: MatGetOwnershipRange(user->Gradxy[0],&istart,&iend);
854: for (i=istart; i<iend; i++){
855: iblock = i / user->mx;
856: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
857: MatSetValues(user->Gradxy[0],1,&i,1,&j,&half_hinv,INSERT_VALUES);
858: j = iblock*user->mx + ((i+1) % user->mx);
859: MatSetValues(user->Gradxy[0],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
860: MatSetValues(user->Gradxy[0],1,&i,1,&i,&zero,INSERT_VALUES);
861: }
862: MatAssemblyBegin(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
863: MatAssemblyEnd(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
865: MatCreate(PETSC_COMM_WORLD,&user->Gradxy[1]);
866: MatSetSizes(user->Gradxy[1],PETSC_DECIDE,PETSC_DECIDE,n,n);
867: MatSetFromOptions(user->Gradxy[1]);
868: MatMPIAIJSetPreallocation(user->Gradxy[1],3,NULL,3,NULL);
869: MatSeqAIJSetPreallocation(user->Gradxy[1],3,NULL);
870: MatGetOwnershipRange(user->Gradxy[1],&istart,&iend);
872: for (i=istart; i<iend; i++){
873: j = (i + n - user->mx) % n;
874: MatSetValues(user->Gradxy[1],1,&i,1,&j,&half_hinv,INSERT_VALUES);
875: j = (j + 2*user->mx) % n;
876: MatSetValues(user->Gradxy[1],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
877: MatSetValues(user->Gradxy[1],1,&i,1,&i,&zero,INSERT_VALUES);
878: }
879: MatAssemblyBegin(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
880: MatAssemblyEnd(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
882: /* Generate Div matrix */
883: MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div);
884: MatTranspose(user->Gradxy[0],MAT_INITIAL_MATRIX,&user->Divxy[0]);
885: MatTranspose(user->Gradxy[1],MAT_INITIAL_MATRIX,&user->Divxy[1]);
887: /* Off-diagonal averaging matrix */
888: MatCreate(PETSC_COMM_WORLD,&user->M);
889: MatSetSizes(user->M,PETSC_DECIDE,PETSC_DECIDE,n,n);
890: MatSetFromOptions(user->M);
891: MatMPIAIJSetPreallocation(user->M,4,NULL,4,NULL);
892: MatSeqAIJSetPreallocation(user->M,4,NULL);
893: MatGetOwnershipRange(user->M,&istart,&iend);
895: for (i=istart; i<iend; i++){
896: /* kron(Id,Av) */
897: iblock = i / user->mx;
898: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
899: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
900: j = iblock*user->mx + ((i+1) % user->mx);
901: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
903: /* kron(Av,Id) */
904: j = (i + user->mx) % n;
905: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
906: j = (i + n - user->mx) % n;
907: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
908: }
909: MatAssemblyBegin(user->M,MAT_FINAL_ASSEMBLY);
910: MatAssemblyEnd(user->M,MAT_FINAL_ASSEMBLY);
912: /* Generate 2D grid */
913: VecCreate(PETSC_COMM_WORLD,&XX);
914: VecCreate(PETSC_COMM_WORLD,&user->q);
915: VecSetSizes(XX,PETSC_DECIDE,n);
916: VecSetSizes(user->q,PETSC_DECIDE,n*user->nt);
917: VecSetFromOptions(XX);
918: VecSetFromOptions(user->q);
920: VecDuplicate(XX,&YY);
921: VecDuplicate(XX,&XXwork);
922: VecDuplicate(XX,&YYwork);
923: VecDuplicate(XX,&user->d);
924: VecDuplicate(XX,&user->dwork);
926: VecGetOwnershipRange(XX,&istart,&iend);
927: for (linear_index=istart; linear_index<iend; linear_index++){
928: i = linear_index % user->mx;
929: j = (linear_index-i)/user->mx;
930: vx = h*(i+0.5);
931: vy = h*(j+0.5);
932: VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
933: VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
934: }
936: VecAssemblyBegin(XX);
937: VecAssemblyEnd(XX);
938: VecAssemblyBegin(YY);
939: VecAssemblyEnd(YY);
941: /* Compute final density function yT
942: yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
943: yT = yT / (h^2*sum(yT)) */
944: VecCopy(XX,XXwork);
945: VecCopy(YY,YYwork);
947: VecShift(XXwork,-0.25);
948: VecShift(YYwork,-0.25);
950: VecPointwiseMult(XXwork,XXwork,XXwork);
951: VecPointwiseMult(YYwork,YYwork,YYwork);
953: VecCopy(XXwork,user->dwork);
954: VecAXPY(user->dwork,1.0,YYwork);
955: VecScale(user->dwork,-30.0);
956: VecExp(user->dwork);
957: VecCopy(user->dwork,user->d);
959: VecCopy(XX,XXwork);
960: VecCopy(YY,YYwork);
962: VecShift(XXwork,-0.75);
963: VecShift(YYwork,-0.75);
965: VecPointwiseMult(XXwork,XXwork,XXwork);
966: VecPointwiseMult(YYwork,YYwork,YYwork);
968: VecCopy(XXwork,user->dwork);
969: VecAXPY(user->dwork,1.0,YYwork);
970: VecScale(user->dwork,-30.0);
971: VecExp(user->dwork);
973: VecAXPY(user->d,1.0,user->dwork);
974: VecShift(user->d,1.0);
975: VecSum(user->d,&sum);
976: VecScale(user->d,1.0/(h*h*sum));
978: /* Initial conditions of forward problem */
979: VecDuplicate(XX,&bc);
980: VecCopy(XX,XXwork);
981: VecCopy(YY,YYwork);
983: VecShift(XXwork,-0.5);
984: VecShift(YYwork,-0.5);
986: VecPointwiseMult(XXwork,XXwork,XXwork);
987: VecPointwiseMult(YYwork,YYwork,YYwork);
989: VecWAXPY(bc,1.0,XXwork,YYwork);
990: VecScale(bc,-50.0);
991: VecExp(bc);
992: VecShift(bc,1.0);
993: VecSum(bc,&sum);
994: VecScale(bc,1.0/(h*h*sum));
996: /* Create scatter from y to y_1,y_2,...,y_nt */
997: /* TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
998: PetscMalloc1(user->nt*user->mx*user->mx,&user->yi_scatter);
999: VecCreate(PETSC_COMM_WORLD,&yi);
1000: VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx);
1001: VecSetFromOptions(yi);
1002: VecDuplicateVecs(yi,user->nt,&user->yi);
1003: VecDuplicateVecs(yi,user->nt,&user->yiwork);
1004: VecDuplicateVecs(yi,user->nt,&user->ziwork);
1005: for (i=0; i<user->nt; i++){
1006: VecGetOwnershipRange(user->yi[i],&lo,&hi);
1007: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_yi);
1008: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+i*user->mx*user->mx,1,&is_from_y);
1009: VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i]);
1010: ISDestroy(&is_to_yi);
1011: ISDestroy(&is_from_y);
1012: }
1014: /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
1015: /* TODO: reorder for better parallelism */
1016: PetscMalloc1(user->nt*user->mx*user->mx,&user->uxi_scatter);
1017: PetscMalloc1(user->nt*user->mx*user->mx,&user->uyi_scatter);
1018: PetscMalloc1(user->nt*user->mx*user->mx,&user->ux_scatter);
1019: PetscMalloc1(user->nt*user->mx*user->mx,&user->uy_scatter);
1020: PetscMalloc1(2*user->nt*user->mx*user->mx,&user->ui_scatter);
1021: VecCreate(PETSC_COMM_WORLD,&uxi);
1022: VecCreate(PETSC_COMM_WORLD,&ui);
1023: VecSetSizes(uxi,PETSC_DECIDE,user->mx*user->mx);
1024: VecSetSizes(ui,PETSC_DECIDE,2*user->mx*user->mx);
1025: VecSetFromOptions(uxi);
1026: VecSetFromOptions(ui);
1027: VecDuplicateVecs(uxi,user->nt,&user->uxi);
1028: VecDuplicateVecs(uxi,user->nt,&user->uyi);
1029: VecDuplicateVecs(uxi,user->nt,&user->uxiwork);
1030: VecDuplicateVecs(uxi,user->nt,&user->uyiwork);
1031: VecDuplicateVecs(ui,user->nt,&user->ui);
1032: VecDuplicateVecs(ui,user->nt,&user->uiwork);
1033: for (i=0; i<user->nt; i++){
1034: VecGetOwnershipRange(user->uxi[i],&lo,&hi);
1035: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1036: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
1037: VecScatterCreate(user->u,is_from_u,user->uxi[i],is_to_uxi,&user->uxi_scatter[i]);
1039: ISDestroy(&is_to_uxi);
1040: ISDestroy(&is_from_u);
1042: VecGetOwnershipRange(user->uyi[i],&lo,&hi);
1043: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
1044: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+(2*i+1)*user->mx*user->mx,1,&is_from_u);
1045: VecScatterCreate(user->u,is_from_u,user->uyi[i],is_to_uyi,&user->uyi_scatter[i]);
1047: ISDestroy(&is_to_uyi);
1048: ISDestroy(&is_from_u);
1050: VecGetOwnershipRange(user->uxi[i],&lo,&hi);
1051: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1052: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_from_u);
1053: VecScatterCreate(user->ui[i],is_from_u,user->uxi[i],is_to_uxi,&user->ux_scatter[i]);
1055: ISDestroy(&is_to_uxi);
1056: ISDestroy(&is_from_u);
1058: VecGetOwnershipRange(user->uyi[i],&lo,&hi);
1059: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
1060: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+user->mx*user->mx,1,&is_from_u);
1061: VecScatterCreate(user->ui[i],is_from_u,user->uyi[i],is_to_uyi,&user->uy_scatter[i]);
1063: ISDestroy(&is_to_uyi);
1064: ISDestroy(&is_from_u);
1066: VecGetOwnershipRange(user->ui[i],&lo,&hi);
1067: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1068: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
1069: VecScatterCreate(user->u,is_from_u,user->ui[i],is_to_uxi,&user->ui_scatter[i]);
1071: ISDestroy(&is_to_uxi);
1072: ISDestroy(&is_from_u);
1073: }
1075: /* RHS of forward problem */
1076: MatMult(user->M,bc,user->yiwork[0]);
1077: for (i=1; i<user->nt; i++){
1078: VecSet(user->yiwork[i],0.0);
1079: }
1080: Gather_yi(user->q,user->yiwork,user->yi_scatter,user->nt);
1082: /* Compute true velocity field utrue */
1083: VecDuplicate(user->u,&user->utrue);
1084: for (i=0; i<user->nt; i++){
1085: VecCopy(YY,user->uxi[i]);
1086: VecScale(user->uxi[i],150.0*i*user->ht);
1087: VecCopy(XX,user->uyi[i]);
1088: VecShift(user->uyi[i],-10.0);
1089: VecScale(user->uyi[i],15.0*i*user->ht);
1090: }
1091: Gather_uxi_uyi(user->utrue,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1093: /* Initial guess and reference model */
1094: VecDuplicate(user->utrue,&user->ur);
1095: for (i=0; i<user->nt; i++){
1096: VecCopy(XX,user->uxi[i]);
1097: VecShift(user->uxi[i],i*user->ht);
1098: VecCopy(YY,user->uyi[i]);
1099: VecShift(user->uyi[i],-i*user->ht);
1100: }
1101: Gather_uxi_uyi(user->ur,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1103: /* Generate regularization matrix L */
1104: MatCreate(PETSC_COMM_WORLD,&user->LT);
1105: MatSetSizes(user->LT,PETSC_DECIDE,PETSC_DECIDE,2*n*user->nt,n*user->nt);
1106: MatSetFromOptions(user->LT);
1107: MatMPIAIJSetPreallocation(user->LT,1,NULL,1,NULL);
1108: MatSeqAIJSetPreallocation(user->LT,1,NULL);
1109: MatGetOwnershipRange(user->LT,&istart,&iend);
1111: for (i=istart; i<iend; i++){
1112: iblock = (i+n) / (2*n);
1113: j = i - iblock*n;
1114: MatSetValues(user->LT,1,&i,1,&j,&user->gamma,INSERT_VALUES);
1115: }
1117: MatAssemblyBegin(user->LT,MAT_FINAL_ASSEMBLY);
1118: MatAssemblyEnd(user->LT,MAT_FINAL_ASSEMBLY);
1120: MatTranspose(user->LT,MAT_INITIAL_MATRIX,&user->L);
1122: /* Build work vectors and matrices */
1123: VecCreate(PETSC_COMM_WORLD,&user->lwork);
1124: VecSetType(user->lwork,VECMPI);
1125: VecSetSizes(user->lwork,PETSC_DECIDE,user->m);
1126: VecSetFromOptions(user->lwork);
1128: MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);
1130: VecDuplicate(user->y,&user->ywork);
1131: VecDuplicate(user->u,&user->uwork);
1132: VecDuplicate(user->u,&user->vwork);
1133: VecDuplicate(user->u,&user->js_diag);
1134: VecDuplicate(user->c,&user->cwork);
1136: /* Create matrix-free shell user->Js for computing A*x */
1137: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->Js);
1138: MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1139: MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
1140: MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
1141: MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
1143: /* Diagonal blocks of user->Js */
1144: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlock);
1145: MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult);
1146: MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockMultTranspose);
1148: /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
1149: D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
1150: This is an SOR preconditioner for user->JsBlock. */
1151: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlockPrec);
1152: MatShellSetOperation(user->JsBlockPrec,MATOP_MULT,(void(*)(void))StateMatBlockPrecMult);
1153: MatShellSetOperation(user->JsBlockPrec,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockPrecMultTranspose);
1155: /* Create a matrix-free shell user->Jd for computing B*x */
1156: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->n-user->m,user,&user->Jd);
1157: MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);
1158: MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);
1160: /* User-defined routines for computing user->Js\x and user->Js^T\x*/
1161: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsInv);
1162: MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateMatInvMult);
1163: MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatInvTransposeMult);
1165: /* Build matrices for SOR preconditioner */
1166: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1167: MatShift(user->Divxy[0],0.0); /* Force C[i] and Divxy[0] to share same nonzero pattern */
1168: MatAXPY(user->Divxy[0],0.0,user->Divxy[1],DIFFERENT_NONZERO_PATTERN);
1169: PetscMalloc1(5*n,&user->C);
1170: PetscMalloc1(2*n,&user->Cwork);
1171: for (i=0; i<user->nt; i++){
1172: MatDuplicate(user->Divxy[0],MAT_COPY_VALUES,&user->C[i]);
1173: MatDuplicate(user->Divxy[1],MAT_COPY_VALUES,&user->Cwork[i]);
1175: MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
1176: MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
1177: MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
1178: MatScale(user->C[i],user->ht);
1179: MatShift(user->C[i],1.0);
1180: }
1182: /* Solver options and tolerances */
1183: KSPCreate(PETSC_COMM_WORLD,&user->solver);
1184: KSPSetType(user->solver,KSPGMRES);
1185: KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec);
1186: KSPSetInitialGuessNonzero(user->solver,PETSC_FALSE); /* TODO: why is true slower? */
1187: KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500);
1188: /* KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500); */
1189: KSPGetPC(user->solver,&user->prec);
1190: PCSetType(user->prec,PCSHELL);
1192: PCShellSetApply(user->prec,StateMatBlockPrecMult);
1193: PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMultTranspose);
1194: PCShellSetContext(user->prec,user);
1196: /* Compute true state function yt given ut */
1197: VecCreate(PETSC_COMM_WORLD,&user->ytrue);
1198: VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt);
1199: VecSetFromOptions(user->ytrue);
1200: user->c_formed = PETSC_TRUE;
1201: VecCopy(user->utrue,user->u); /* Set u=utrue temporarily for StateMatInv */
1202: VecSet(user->ytrue,0.0); /* Initial guess */
1203: StateMatInvMult(user->Js,user->q,user->ytrue);
1204: VecCopy(user->ur,user->u); /* Reset u=ur */
1206: /* Initial guess y0 for state given u0 */
1207: user->lcl->solve_type = LCL_FORWARD1;
1208: StateMatInvMult(user->Js,user->q,user->y);
1210: /* Data discretization */
1211: MatCreate(PETSC_COMM_WORLD,&user->Q);
1212: MatSetSizes(user->Q,PETSC_DECIDE,PETSC_DECIDE,user->mx*user->mx,user->m);
1213: MatSetFromOptions(user->Q);
1214: MatMPIAIJSetPreallocation(user->Q,0,NULL,1,NULL);
1215: MatSeqAIJSetPreallocation(user->Q,1,NULL);
1217: MatGetOwnershipRange(user->Q,&istart,&iend);
1219: for (i=istart; i<iend; i++){
1220: j = i + user->m - user->mx*user->mx;
1221: MatSetValues(user->Q,1,&i,1,&j,&one,INSERT_VALUES);
1222: }
1224: MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1225: MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);
1227: MatTranspose(user->Q,MAT_INITIAL_MATRIX,&user->QT);
1229: VecDestroy(&XX);
1230: VecDestroy(&YY);
1231: VecDestroy(&XXwork);
1232: VecDestroy(&YYwork);
1233: VecDestroy(&yi);
1234: VecDestroy(&uxi);
1235: VecDestroy(&ui);
1236: VecDestroy(&bc);
1238: /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1239: KSPSetFromOptions(user->solver);
1240: return(0);
1241: }
1245: PetscErrorCode HyperbolicDestroy(AppCtx *user)
1246: {
1248: PetscInt i;
1251: MatDestroy(&user->Q);
1252: MatDestroy(&user->QT);
1253: MatDestroy(&user->Div);
1254: MatDestroy(&user->Divwork);
1255: MatDestroy(&user->Grad);
1256: MatDestroy(&user->L);
1257: MatDestroy(&user->LT);
1258: KSPDestroy(&user->solver);
1259: MatDestroy(&user->Js);
1260: MatDestroy(&user->Jd);
1261: MatDestroy(&user->JsBlockPrec);
1262: MatDestroy(&user->JsInv);
1263: MatDestroy(&user->JsBlock);
1264: MatDestroy(&user->Divxy[0]);
1265: MatDestroy(&user->Divxy[1]);
1266: MatDestroy(&user->Gradxy[0]);
1267: MatDestroy(&user->Gradxy[1]);
1268: MatDestroy(&user->M);
1269: for (i=0; i<user->nt; i++){
1270: MatDestroy(&user->C[i]);
1271: MatDestroy(&user->Cwork[i]);
1272: }
1273: PetscFree(user->C);
1274: PetscFree(user->Cwork);
1275: VecDestroy(&user->u);
1276: VecDestroy(&user->uwork);
1277: VecDestroy(&user->vwork);
1278: VecDestroy(&user->utrue);
1279: VecDestroy(&user->y);
1280: VecDestroy(&user->ywork);
1281: VecDestroy(&user->ytrue);
1282: VecDestroyVecs(user->nt,&user->yi);
1283: VecDestroyVecs(user->nt,&user->yiwork);
1284: VecDestroyVecs(user->nt,&user->ziwork);
1285: VecDestroyVecs(user->nt,&user->uxi);
1286: VecDestroyVecs(user->nt,&user->uyi);
1287: VecDestroyVecs(user->nt,&user->uxiwork);
1288: VecDestroyVecs(user->nt,&user->uyiwork);
1289: VecDestroyVecs(user->nt,&user->ui);
1290: VecDestroyVecs(user->nt,&user->uiwork);
1291: VecDestroy(&user->c);
1292: VecDestroy(&user->cwork);
1293: VecDestroy(&user->ur);
1294: VecDestroy(&user->q);
1295: VecDestroy(&user->d);
1296: VecDestroy(&user->dwork);
1297: VecDestroy(&user->lwork);
1298: VecDestroy(&user->js_diag);
1299: ISDestroy(&user->s_is);
1300: ISDestroy(&user->d_is);
1301: VecScatterDestroy(&user->state_scatter);
1302: VecScatterDestroy(&user->design_scatter);
1303: for (i=0; i<user->nt; i++){
1304: VecScatterDestroy(&user->uxi_scatter[i]);
1305: VecScatterDestroy(&user->uyi_scatter[i]);
1306: VecScatterDestroy(&user->ux_scatter[i]);
1307: VecScatterDestroy(&user->uy_scatter[i]);
1308: VecScatterDestroy(&user->ui_scatter[i]);
1309: VecScatterDestroy(&user->yi_scatter[i]);
1310: }
1311: PetscFree(user->uxi_scatter);
1312: PetscFree(user->uyi_scatter);
1313: PetscFree(user->ux_scatter);
1314: PetscFree(user->uy_scatter);
1315: PetscFree(user->ui_scatter);
1316: PetscFree(user->yi_scatter);
1317: return(0);
1318: }
1322: PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1323: {
1325: Vec X;
1326: PetscReal unorm,ynorm;
1327: AppCtx *user = (AppCtx*)ptr;
1330: TaoGetSolutionVector(tao,&X);
1331: Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1332: VecAXPY(user->ywork,-1.0,user->ytrue);
1333: VecAXPY(user->uwork,-1.0,user->utrue);
1334: VecNorm(user->uwork,NORM_2,&unorm);
1335: VecNorm(user->ywork,NORM_2,&ynorm);
1336: PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1337: return(0);
1338: }