Actual source code: hyperbolic.c
petsc-3.7.5 2017-01-01
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: 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: AppCtx user;
121: IS is_allstate,is_alldesign;
122: PetscInt lo,hi,hi2,lo2,ksp_old;
123: PetscInt ntests = 1;
124: PetscInt i;
125: #if defined(PETSC_USE_LOG)
126: PetscLogStage stages[1];
127: #endif
128:
129: PetscInitialize(&argc, &argv, (char*)0,help);
131: user.mx = 32;
132: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);
133: PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);
134: user.nt = 16;
135: PetscOptionsInt("-nt","Number of time steps","",user.nt,&user.nt,NULL);
136: user.ndata = 64;
137: PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL);
138: user.alpha = 10.0;
139: PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);
140: user.T = 1.0/32.0;
141: PetscOptionsReal("-Tfinal","Final time","",user.T,&user.T,NULL);
143: user.m = user.mx*user.mx*user.nt; /* number of constraints */
144: user.n = user.mx*user.mx*3*user.nt; /* number of variables */
145: user.ht = user.T/user.nt; /* Time step */
146: user.gamma = user.T*user.ht / (user.mx*user.mx);
148: VecCreate(PETSC_COMM_WORLD,&user.u);
149: VecCreate(PETSC_COMM_WORLD,&user.y);
150: VecCreate(PETSC_COMM_WORLD,&user.c);
151: VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m);
152: VecSetSizes(user.y,PETSC_DECIDE,user.m);
153: VecSetSizes(user.c,PETSC_DECIDE,user.m);
154: VecSetFromOptions(user.u);
155: VecSetFromOptions(user.y);
156: VecSetFromOptions(user.c);
158: /* Create scatters for reduced spaces.
159: If the state vector y and design vector u are partitioned as
160: [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
161: then the solution vector x is organized as
162: [y_1; u_1; y_2; u_2; ...; y_np; u_np].
163: The index sets user.s_is and user.d_is correspond to the indices of the
164: state and design variables owned by the current processor.
165: */
166: VecCreate(PETSC_COMM_WORLD,&x);
168: VecGetOwnershipRange(user.y,&lo,&hi);
169: VecGetOwnershipRange(user.u,&lo2,&hi2);
171: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);
172: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is);
174: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);
175: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is);
177: VecSetSizes(x,hi-lo+hi2-lo2,user.n);
178: VecSetFromOptions(x);
180: VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter);
181: VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter);
182: ISDestroy(&is_alldesign);
183: ISDestroy(&is_allstate);
185: /* Create TAO solver and set desired solution method */
186: TaoCreate(PETSC_COMM_WORLD,&tao);
187: TaoSetType(tao,TAOLCL);
189: /* Set up initial vectors and matrices */
190: HyperbolicInitialize(&user);
192: Gather(x,user.y,user.state_scatter,user.u,user.design_scatter);
193: VecDuplicate(x,&x0);
194: VecCopy(x,x0);
196: /* Set solution vector with an initial guess */
197: TaoSetInitialVector(tao,x);
198: TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);
199: TaoSetGradientRoutine(tao, FormGradient, (void *)&user);
200: TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);
202: TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, (void *)&user);
204: TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);
206: TaoSetFromOptions(tao);
207: TaoSetStateDesignIS(tao,user.s_is,user.d_is);
209: /* SOLVE THE APPLICATION */
210: PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL);
211: PetscOptionsEnd();
212: PetscLogStageRegister("Trials",&stages[0]);
213: PetscLogStagePush(stages[0]);
214: user.ksp_its_initial = user.ksp_its;
215: ksp_old = user.ksp_its;
216: for (i=0; i<ntests; i++){
217: TaoSolve(tao);
218: PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old);
219: VecCopy(x0,x);
220: TaoSetInitialVector(tao,x);
221: }
222: PetscLogStagePop();
223: PetscBarrier((PetscObject)x);
224: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
225: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);
226: PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests);
227: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its);
228: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations per trial: ");
229: PetscPrintf(PETSC_COMM_WORLD,"%D\n",(user.ksp_its-user.ksp_its_initial)/ntests);
231: TaoDestroy(&tao);
232: VecDestroy(&x);
233: VecDestroy(&x0);
234: HyperbolicDestroy(&user);
235: PetscFinalize();
236: return 0;
237: }
238: /* ------------------------------------------------------------------- */
241: /*
242: dwork = Qy - d
243: lwork = L*(u-ur).^2
244: f = 1/2 * (dwork.dork + alpha*y.lwork)
245: */
246: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
247: {
249: PetscReal d1=0,d2=0;
250: AppCtx *user = (AppCtx*)ptr;
253: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
254: MatMult(user->Q,user->y,user->dwork);
255: VecAXPY(user->dwork,-1.0,user->d);
256: VecDot(user->dwork,user->dwork,&d1);
258: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
259: VecPointwiseMult(user->uwork,user->uwork,user->uwork);
260: MatMult(user->L,user->uwork,user->lwork);
261: VecDot(user->y,user->lwork,&d2);
262: *f = 0.5 * (d1 + user->alpha*d2);
263: return(0);
264: }
266: /* ------------------------------------------------------------------- */
269: /*
270: state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2
271: design: g_d = alpha*(L'y).*(u-ur)
272: */
273: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
274: {
276: AppCtx *user = (AppCtx*)ptr;
279: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
280: MatMult(user->Q,user->y,user->dwork);
281: VecAXPY(user->dwork,-1.0,user->d);
283: MatMult(user->QT,user->dwork,user->ywork);
285: MatMult(user->LT,user->y,user->uwork);
286: VecWAXPY(user->vwork,-1.0,user->ur,user->u);
287: VecPointwiseMult(user->uwork,user->vwork,user->uwork);
288: VecScale(user->uwork,user->alpha);
290: VecPointwiseMult(user->vwork,user->vwork,user->vwork);
291: MatMult(user->L,user->vwork,user->lwork);
292: VecAXPY(user->ywork,0.5*user->alpha,user->lwork);
294: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
295: return(0);
296: }
300: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
301: {
303: PetscReal d1,d2;
304: AppCtx *user = (AppCtx*)ptr;
307: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
308: MatMult(user->Q,user->y,user->dwork);
309: VecAXPY(user->dwork,-1.0,user->d);
311: MatMult(user->QT,user->dwork,user->ywork);
313: VecDot(user->dwork,user->dwork,&d1);
315: MatMult(user->LT,user->y,user->uwork);
316: VecWAXPY(user->vwork,-1.0,user->ur,user->u);
317: VecPointwiseMult(user->uwork,user->vwork,user->uwork);
318: VecScale(user->uwork,user->alpha);
320: VecPointwiseMult(user->vwork,user->vwork,user->vwork);
321: MatMult(user->L,user->vwork,user->lwork);
322: VecAXPY(user->ywork,0.5*user->alpha,user->lwork);
324: VecDot(user->y,user->lwork,&d2);
326: *f = 0.5 * (d1 + user->alpha*d2);
327: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
328: return(0);
329: }
331: /* ------------------------------------------------------------------- */
334: /* A
335: MatShell object
336: */
337: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
338: {
340: PetscInt i;
341: AppCtx *user = (AppCtx*)ptr;
344: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
345: Scatter_yi(user->u,user->ui,user->ui_scatter,user->nt);
346: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
347: for (i=0; i<user->nt; i++){
348: MatCopy(user->Divxy[0],user->C[i],SAME_NONZERO_PATTERN);
349: MatCopy(user->Divxy[1],user->Cwork[i],SAME_NONZERO_PATTERN);
351: MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
352: MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
353: MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
354: MatScale(user->C[i],user->ht);
355: MatShift(user->C[i],1.0);
356: }
357: return(0);
358: }
360: /* ------------------------------------------------------------------- */
363: /* B */
364: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
365: {
367: AppCtx *user = (AppCtx*)ptr;
370: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
371: return(0);
372: }
376: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
377: {
379: PetscInt i;
380: AppCtx *user;
383: MatShellGetContext(J_shell,(void**)&user);
384: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
385: user->block_index = 0;
386: MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);
388: for (i=1; i<user->nt; i++){
389: user->block_index = i;
390: MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
391: MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
392: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
393: }
394: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
395: return(0);
396: }
400: PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
401: {
403: PetscInt i;
404: AppCtx *user;
407: MatShellGetContext(J_shell,(void**)&user);
408: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
410: for (i=0; i<user->nt-1; i++){
411: user->block_index = i;
412: MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
413: MatMult(user->M,user->yi[i+1],user->ziwork[i+1]);
414: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i+1]);
415: }
417: i = user->nt-1;
418: user->block_index = i;
419: MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
420: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
421: return(0);
422: }
426: PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
427: {
429: PetscInt i;
430: AppCtx *user;
433: MatShellGetContext(J_shell,(void**)&user);
434: i = user->block_index;
435: VecPointwiseMult(user->uxiwork[i],X,user->uxi[i]);
436: VecPointwiseMult(user->uyiwork[i],X,user->uyi[i]);
437: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
438: MatMult(user->Div,user->uiwork[i],Y);
439: VecAYPX(Y,user->ht,X);
440: return(0);
441: }
445: PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y)
446: {
448: PetscInt i;
449: AppCtx *user;
452: MatShellGetContext(J_shell,(void**)&user);
453: i = user->block_index;
454: MatMult(user->Grad,X,user->uiwork[i]);
455: Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
456: VecPointwiseMult(user->uxiwork[i],user->uxi[i],user->uxiwork[i]);
457: VecPointwiseMult(user->uyiwork[i],user->uyi[i],user->uyiwork[i]);
458: VecWAXPY(Y,1.0,user->uxiwork[i],user->uyiwork[i]);
459: VecAYPX(Y,user->ht,X);
460: return(0);
461: }
465: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
466: {
468: PetscInt i;
469: AppCtx *user;
472: MatShellGetContext(J_shell,(void**)&user);
473: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
474: Scatter_uxi_uyi(X,user->uxiwork,user->uxi_scatter,user->uyiwork,user->uyi_scatter,user->nt);
475: for (i=0; i<user->nt; i++){
476: VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
477: VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
478: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
479: MatMult(user->Div,user->uiwork[i],user->ziwork[i]);
480: VecScale(user->ziwork[i],user->ht);
481: }
482: Gather_yi(Y,user->ziwork,user->yi_scatter,user->nt);
483: return(0);
484: }
488: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
489: {
491: PetscInt i;
492: AppCtx *user;
495: MatShellGetContext(J_shell,(void**)&user);
496: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
497: Scatter_yi(X,user->yiwork,user->yi_scatter,user->nt);
498: for (i=0; i<user->nt; i++){
499: MatMult(user->Grad,user->yiwork[i],user->uiwork[i]);
500: Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
501: VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
502: VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
503: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
504: VecScale(user->uiwork[i],user->ht);
505: }
506: Gather_yi(Y,user->uiwork,user->ui_scatter,user->nt);
507: return(0);
508: }
512: PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
513: {
515: PetscInt i;
516: AppCtx *user;
519: PCShellGetContext(PC_shell,(void**)&user);
520: i = user->block_index;
521: if (user->c_formed) {
522: MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
523: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
524: return(0);
525: }
529: PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
530: {
532: PetscInt i;
533: AppCtx *user;
536: PCShellGetContext(PC_shell,(void**)&user);
538: i = user->block_index;
539: if (user->c_formed) {
540: MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
541: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
542: return(0);
543: }
547: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
548: {
550: AppCtx *user;
551: PetscInt its,i;
554: MatShellGetContext(J_shell,(void**)&user);
556: if (Y == user->ytrue) {
557: /* First solve is done using true solution to set up problem */
558: KSPSetTolerances(user->solver,1e-4,1e-20,PETSC_DEFAULT,PETSC_DEFAULT);
559: } else {
560: KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
561: }
562: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
563: Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
564: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
566: user->block_index = 0;
567: KSPSolve(user->solver,user->yi[0],user->yiwork[0]);
569: KSPGetIterationNumber(user->solver,&its);
570: user->ksp_its = user->ksp_its + its;
571: for (i=1; i<user->nt; i++){
572: MatMult(user->M,user->yiwork[i-1],user->ziwork[i-1]);
573: VecAXPY(user->yi[i],1.0,user->ziwork[i-1]);
574: user->block_index = i;
575: KSPSolve(user->solver,user->yi[i],user->yiwork[i]);
577: KSPGetIterationNumber(user->solver,&its);
578: user->ksp_its = user->ksp_its + its;
579: }
580: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
581: return(0);
582: }
586: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
587: {
589: AppCtx *user;
590: PetscInt its,i;
593: MatShellGetContext(J_shell,(void**)&user);
595: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
596: Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
597: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
599: i = user->nt - 1;
600: user->block_index = i;
601: KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);
603: KSPGetIterationNumber(user->solver,&its);
604: user->ksp_its = user->ksp_its + its;
606: for (i=user->nt-2; i>=0; i--){
607: MatMult(user->M,user->yiwork[i+1],user->ziwork[i+1]);
608: VecAXPY(user->yi[i],1.0,user->ziwork[i+1]);
609: user->block_index = i;
610: KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);
612: KSPGetIterationNumber(user->solver,&its);
613: user->ksp_its = user->ksp_its + its;
614: }
615: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
616: return(0);
617: }
621: PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
622: {
624: AppCtx *user;
627: MatShellGetContext(J_shell,(void**)&user);
629: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_shell);
630: MatShellSetOperation(*new_shell,MATOP_MULT,(void(*)(void))StateMatMult);
631: MatShellSetOperation(*new_shell,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
632: MatShellSetOperation(*new_shell,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
633: MatShellSetOperation(*new_shell,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
634: return(0);
635: }
639: PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
640: {
642: AppCtx *user;
645: MatShellGetContext(J_shell,(void**)&user);
646: VecCopy(user->js_diag,X);
647: return(0);
648: }
652: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
653: {
654: /* con = Ay - q, A = [C(u1) 0 0 ... 0;
655: -M C(u2) 0 ... 0;
656: 0 -M C(u3) ... 0;
657: ... ;
658: 0 ... -M C(u_nt)]
659: C(u) = eye + ht*Div*[diag(u1); diag(u2)] */
661: PetscInt i;
662: AppCtx *user = (AppCtx*)ptr;
665: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
666: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
667: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
669: user->block_index = 0;
670: MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);
672: for (i=1; i<user->nt; i++){
673: user->block_index = i;
674: MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
675: MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
676: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
677: }
679: Gather_yi(C,user->yiwork,user->yi_scatter,user->nt);
680: VecAXPY(C,-1.0,user->q);
682: return(0);
683: }
688: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
689: {
693: VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
694: VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
695: VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
696: VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
697: return(0);
698: }
702: PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
703: {
705: PetscInt i;
708: for (i=0; i<nt; i++){
709: VecScatterBegin(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
710: VecScatterEnd(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
711: VecScatterBegin(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
712: VecScatterEnd(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
713: }
714: return(0);
715: }
719: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
720: {
724: VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
725: VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
726: VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
727: VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
728: return(0);
729: }
733: PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
734: {
736: PetscInt i;
739: for (i=0; i<nt; i++){
740: VecScatterBegin(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
741: VecScatterEnd(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
742: VecScatterBegin(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
743: VecScatterEnd(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
744: }
745: return(0);
746: }
750: PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
751: {
753: PetscInt i;
756: for (i=0; i<nt; i++){
757: VecScatterBegin(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
758: VecScatterEnd(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
759: }
760: return(0);
761: }
765: PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
766: {
768: PetscInt i;
771: for (i=0; i<nt; i++){
772: VecScatterBegin(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
773: VecScatterEnd(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
774: }
775: return(0);
776: }
780: PetscErrorCode HyperbolicInitialize(AppCtx *user)
781: {
783: PetscInt n,i,j,linear_index,istart,iend,iblock,lo,hi;
784: Vec XX,YY,XXwork,YYwork,yi,uxi,ui,bc;
785: PetscReal h,sum;
786: PetscScalar hinv,neg_hinv,quarter=0.25,one=1.0,half_hinv,neg_half_hinv;
787: PetscScalar vx,vy,zero=0.0;
788: IS is_from_y,is_to_yi,is_from_u,is_to_uxi,is_to_uyi;
791: user->jformed = PETSC_FALSE;
792: user->c_formed = PETSC_FALSE;
794: user->ksp_its = 0;
795: user->ksp_its_initial = 0;
797: n = user->mx * user->mx;
799: h = 1.0/user->mx;
800: hinv = user->mx;
801: neg_hinv = -hinv;
802: half_hinv = hinv / 2.0;
803: neg_half_hinv = neg_hinv / 2.0;
805: /* Generate Grad matrix */
806: MatCreate(PETSC_COMM_WORLD,&user->Grad);
807: MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,2*n,n);
808: MatSetFromOptions(user->Grad);
809: MatMPIAIJSetPreallocation(user->Grad,3,NULL,3,NULL);
810: MatSeqAIJSetPreallocation(user->Grad,3,NULL);
811: MatGetOwnershipRange(user->Grad,&istart,&iend);
813: for (i=istart; i<iend; i++){
814: if (i<n){
815: iblock = i / user->mx;
816: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
817: MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
818: j = iblock*user->mx + ((i+1) % user->mx);
819: MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
820: }
821: if (i>=n){
822: j = (i - user->mx) % n;
823: MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
824: j = (j + 2*user->mx) % n;
825: MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
826: }
827: }
829: MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
830: MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);
832: MatCreate(PETSC_COMM_WORLD,&user->Gradxy[0]);
833: MatSetSizes(user->Gradxy[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
834: MatSetFromOptions(user->Gradxy[0]);
835: MatMPIAIJSetPreallocation(user->Gradxy[0],3,NULL,3,NULL);
836: MatSeqAIJSetPreallocation(user->Gradxy[0],3,NULL);
837: MatGetOwnershipRange(user->Gradxy[0],&istart,&iend);
839: for (i=istart; i<iend; i++){
840: iblock = i / user->mx;
841: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
842: MatSetValues(user->Gradxy[0],1,&i,1,&j,&half_hinv,INSERT_VALUES);
843: j = iblock*user->mx + ((i+1) % user->mx);
844: MatSetValues(user->Gradxy[0],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
845: MatSetValues(user->Gradxy[0],1,&i,1,&i,&zero,INSERT_VALUES);
846: }
847: MatAssemblyBegin(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
848: MatAssemblyEnd(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
850: MatCreate(PETSC_COMM_WORLD,&user->Gradxy[1]);
851: MatSetSizes(user->Gradxy[1],PETSC_DECIDE,PETSC_DECIDE,n,n);
852: MatSetFromOptions(user->Gradxy[1]);
853: MatMPIAIJSetPreallocation(user->Gradxy[1],3,NULL,3,NULL);
854: MatSeqAIJSetPreallocation(user->Gradxy[1],3,NULL);
855: MatGetOwnershipRange(user->Gradxy[1],&istart,&iend);
857: for (i=istart; i<iend; i++){
858: j = (i + n - user->mx) % n;
859: MatSetValues(user->Gradxy[1],1,&i,1,&j,&half_hinv,INSERT_VALUES);
860: j = (j + 2*user->mx) % n;
861: MatSetValues(user->Gradxy[1],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
862: MatSetValues(user->Gradxy[1],1,&i,1,&i,&zero,INSERT_VALUES);
863: }
864: MatAssemblyBegin(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
865: MatAssemblyEnd(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
867: /* Generate Div matrix */
868: MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div);
869: MatTranspose(user->Gradxy[0],MAT_INITIAL_MATRIX,&user->Divxy[0]);
870: MatTranspose(user->Gradxy[1],MAT_INITIAL_MATRIX,&user->Divxy[1]);
872: /* Off-diagonal averaging matrix */
873: MatCreate(PETSC_COMM_WORLD,&user->M);
874: MatSetSizes(user->M,PETSC_DECIDE,PETSC_DECIDE,n,n);
875: MatSetFromOptions(user->M);
876: MatMPIAIJSetPreallocation(user->M,4,NULL,4,NULL);
877: MatSeqAIJSetPreallocation(user->M,4,NULL);
878: MatGetOwnershipRange(user->M,&istart,&iend);
880: for (i=istart; i<iend; i++){
881: /* kron(Id,Av) */
882: iblock = i / user->mx;
883: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
884: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
885: j = iblock*user->mx + ((i+1) % user->mx);
886: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
888: /* kron(Av,Id) */
889: j = (i + user->mx) % n;
890: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
891: j = (i + n - user->mx) % n;
892: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
893: }
894: MatAssemblyBegin(user->M,MAT_FINAL_ASSEMBLY);
895: MatAssemblyEnd(user->M,MAT_FINAL_ASSEMBLY);
897: /* Generate 2D grid */
898: VecCreate(PETSC_COMM_WORLD,&XX);
899: VecCreate(PETSC_COMM_WORLD,&user->q);
900: VecSetSizes(XX,PETSC_DECIDE,n);
901: VecSetSizes(user->q,PETSC_DECIDE,n*user->nt);
902: VecSetFromOptions(XX);
903: VecSetFromOptions(user->q);
905: VecDuplicate(XX,&YY);
906: VecDuplicate(XX,&XXwork);
907: VecDuplicate(XX,&YYwork);
908: VecDuplicate(XX,&user->d);
909: VecDuplicate(XX,&user->dwork);
911: VecGetOwnershipRange(XX,&istart,&iend);
912: for (linear_index=istart; linear_index<iend; linear_index++){
913: i = linear_index % user->mx;
914: j = (linear_index-i)/user->mx;
915: vx = h*(i+0.5);
916: vy = h*(j+0.5);
917: VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
918: VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
919: }
921: VecAssemblyBegin(XX);
922: VecAssemblyEnd(XX);
923: VecAssemblyBegin(YY);
924: VecAssemblyEnd(YY);
926: /* Compute final density function yT
927: yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
928: yT = yT / (h^2*sum(yT)) */
929: VecCopy(XX,XXwork);
930: VecCopy(YY,YYwork);
932: VecShift(XXwork,-0.25);
933: VecShift(YYwork,-0.25);
935: VecPointwiseMult(XXwork,XXwork,XXwork);
936: VecPointwiseMult(YYwork,YYwork,YYwork);
938: VecCopy(XXwork,user->dwork);
939: VecAXPY(user->dwork,1.0,YYwork);
940: VecScale(user->dwork,-30.0);
941: VecExp(user->dwork);
942: VecCopy(user->dwork,user->d);
944: VecCopy(XX,XXwork);
945: VecCopy(YY,YYwork);
947: VecShift(XXwork,-0.75);
948: VecShift(YYwork,-0.75);
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);
958: VecAXPY(user->d,1.0,user->dwork);
959: VecShift(user->d,1.0);
960: VecSum(user->d,&sum);
961: VecScale(user->d,1.0/(h*h*sum));
963: /* Initial conditions of forward problem */
964: VecDuplicate(XX,&bc);
965: VecCopy(XX,XXwork);
966: VecCopy(YY,YYwork);
968: VecShift(XXwork,-0.5);
969: VecShift(YYwork,-0.5);
971: VecPointwiseMult(XXwork,XXwork,XXwork);
972: VecPointwiseMult(YYwork,YYwork,YYwork);
974: VecWAXPY(bc,1.0,XXwork,YYwork);
975: VecScale(bc,-50.0);
976: VecExp(bc);
977: VecShift(bc,1.0);
978: VecSum(bc,&sum);
979: VecScale(bc,1.0/(h*h*sum));
981: /* Create scatter from y to y_1,y_2,...,y_nt */
982: /* TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
983: PetscMalloc1(user->nt*user->mx*user->mx,&user->yi_scatter);
984: VecCreate(PETSC_COMM_WORLD,&yi);
985: VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx);
986: VecSetFromOptions(yi);
987: VecDuplicateVecs(yi,user->nt,&user->yi);
988: VecDuplicateVecs(yi,user->nt,&user->yiwork);
989: VecDuplicateVecs(yi,user->nt,&user->ziwork);
990: for (i=0; i<user->nt; i++){
991: VecGetOwnershipRange(user->yi[i],&lo,&hi);
992: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_yi);
993: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+i*user->mx*user->mx,1,&is_from_y);
994: VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i]);
995: ISDestroy(&is_to_yi);
996: ISDestroy(&is_from_y);
997: }
999: /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
1000: /* TODO: reorder for better parallelism */
1001: PetscMalloc1(user->nt*user->mx*user->mx,&user->uxi_scatter);
1002: PetscMalloc1(user->nt*user->mx*user->mx,&user->uyi_scatter);
1003: PetscMalloc1(user->nt*user->mx*user->mx,&user->ux_scatter);
1004: PetscMalloc1(user->nt*user->mx*user->mx,&user->uy_scatter);
1005: PetscMalloc1(2*user->nt*user->mx*user->mx,&user->ui_scatter);
1006: VecCreate(PETSC_COMM_WORLD,&uxi);
1007: VecCreate(PETSC_COMM_WORLD,&ui);
1008: VecSetSizes(uxi,PETSC_DECIDE,user->mx*user->mx);
1009: VecSetSizes(ui,PETSC_DECIDE,2*user->mx*user->mx);
1010: VecSetFromOptions(uxi);
1011: VecSetFromOptions(ui);
1012: VecDuplicateVecs(uxi,user->nt,&user->uxi);
1013: VecDuplicateVecs(uxi,user->nt,&user->uyi);
1014: VecDuplicateVecs(uxi,user->nt,&user->uxiwork);
1015: VecDuplicateVecs(uxi,user->nt,&user->uyiwork);
1016: VecDuplicateVecs(ui,user->nt,&user->ui);
1017: VecDuplicateVecs(ui,user->nt,&user->uiwork);
1018: for (i=0; i<user->nt; i++){
1019: VecGetOwnershipRange(user->uxi[i],&lo,&hi);
1020: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1021: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
1022: VecScatterCreate(user->u,is_from_u,user->uxi[i],is_to_uxi,&user->uxi_scatter[i]);
1024: ISDestroy(&is_to_uxi);
1025: ISDestroy(&is_from_u);
1027: VecGetOwnershipRange(user->uyi[i],&lo,&hi);
1028: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
1029: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+(2*i+1)*user->mx*user->mx,1,&is_from_u);
1030: VecScatterCreate(user->u,is_from_u,user->uyi[i],is_to_uyi,&user->uyi_scatter[i]);
1032: ISDestroy(&is_to_uyi);
1033: ISDestroy(&is_from_u);
1035: VecGetOwnershipRange(user->uxi[i],&lo,&hi);
1036: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1037: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_from_u);
1038: VecScatterCreate(user->ui[i],is_from_u,user->uxi[i],is_to_uxi,&user->ux_scatter[i]);
1040: ISDestroy(&is_to_uxi);
1041: ISDestroy(&is_from_u);
1043: VecGetOwnershipRange(user->uyi[i],&lo,&hi);
1044: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
1045: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+user->mx*user->mx,1,&is_from_u);
1046: VecScatterCreate(user->ui[i],is_from_u,user->uyi[i],is_to_uyi,&user->uy_scatter[i]);
1048: ISDestroy(&is_to_uyi);
1049: ISDestroy(&is_from_u);
1051: VecGetOwnershipRange(user->ui[i],&lo,&hi);
1052: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1053: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
1054: VecScatterCreate(user->u,is_from_u,user->ui[i],is_to_uxi,&user->ui_scatter[i]);
1056: ISDestroy(&is_to_uxi);
1057: ISDestroy(&is_from_u);
1058: }
1060: /* RHS of forward problem */
1061: MatMult(user->M,bc,user->yiwork[0]);
1062: for (i=1; i<user->nt; i++){
1063: VecSet(user->yiwork[i],0.0);
1064: }
1065: Gather_yi(user->q,user->yiwork,user->yi_scatter,user->nt);
1067: /* Compute true velocity field utrue */
1068: VecDuplicate(user->u,&user->utrue);
1069: for (i=0; i<user->nt; i++){
1070: VecCopy(YY,user->uxi[i]);
1071: VecScale(user->uxi[i],150.0*i*user->ht);
1072: VecCopy(XX,user->uyi[i]);
1073: VecShift(user->uyi[i],-10.0);
1074: VecScale(user->uyi[i],15.0*i*user->ht);
1075: }
1076: Gather_uxi_uyi(user->utrue,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1078: /* Initial guess and reference model */
1079: VecDuplicate(user->utrue,&user->ur);
1080: for (i=0; i<user->nt; i++){
1081: VecCopy(XX,user->uxi[i]);
1082: VecShift(user->uxi[i],i*user->ht);
1083: VecCopy(YY,user->uyi[i]);
1084: VecShift(user->uyi[i],-i*user->ht);
1085: }
1086: Gather_uxi_uyi(user->ur,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1088: /* Generate regularization matrix L */
1089: MatCreate(PETSC_COMM_WORLD,&user->LT);
1090: MatSetSizes(user->LT,PETSC_DECIDE,PETSC_DECIDE,2*n*user->nt,n*user->nt);
1091: MatSetFromOptions(user->LT);
1092: MatMPIAIJSetPreallocation(user->LT,1,NULL,1,NULL);
1093: MatSeqAIJSetPreallocation(user->LT,1,NULL);
1094: MatGetOwnershipRange(user->LT,&istart,&iend);
1096: for (i=istart; i<iend; i++){
1097: iblock = (i+n) / (2*n);
1098: j = i - iblock*n;
1099: MatSetValues(user->LT,1,&i,1,&j,&user->gamma,INSERT_VALUES);
1100: }
1102: MatAssemblyBegin(user->LT,MAT_FINAL_ASSEMBLY);
1103: MatAssemblyEnd(user->LT,MAT_FINAL_ASSEMBLY);
1105: MatTranspose(user->LT,MAT_INITIAL_MATRIX,&user->L);
1107: /* Build work vectors and matrices */
1108: VecCreate(PETSC_COMM_WORLD,&user->lwork);
1109: VecSetType(user->lwork,VECMPI);
1110: VecSetSizes(user->lwork,PETSC_DECIDE,user->m);
1111: VecSetFromOptions(user->lwork);
1113: MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);
1115: VecDuplicate(user->y,&user->ywork);
1116: VecDuplicate(user->u,&user->uwork);
1117: VecDuplicate(user->u,&user->vwork);
1118: VecDuplicate(user->u,&user->js_diag);
1119: VecDuplicate(user->c,&user->cwork);
1121: /* Create matrix-free shell user->Js for computing A*x */
1122: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->Js);
1123: MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1124: MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
1125: MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
1126: MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
1128: /* Diagonal blocks of user->Js */
1129: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlock);
1130: MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult);
1131: MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockMultTranspose);
1133: /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
1134: D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
1135: This is an SOR preconditioner for user->JsBlock. */
1136: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlockPrec);
1137: MatShellSetOperation(user->JsBlockPrec,MATOP_MULT,(void(*)(void))StateMatBlockPrecMult);
1138: MatShellSetOperation(user->JsBlockPrec,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockPrecMultTranspose);
1140: /* Create a matrix-free shell user->Jd for computing B*x */
1141: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->n-user->m,user,&user->Jd);
1142: MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);
1143: MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);
1145: /* User-defined routines for computing user->Js\x and user->Js^T\x*/
1146: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsInv);
1147: MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateMatInvMult);
1148: MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatInvTransposeMult);
1150: /* Build matrices for SOR preconditioner */
1151: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1152: MatShift(user->Divxy[0],0.0); /* Force C[i] and Divxy[0] to share same nonzero pattern */
1153: MatAXPY(user->Divxy[0],0.0,user->Divxy[1],DIFFERENT_NONZERO_PATTERN);
1154: PetscMalloc1(5*n,&user->C);
1155: PetscMalloc1(2*n,&user->Cwork);
1156: for (i=0; i<user->nt; i++){
1157: MatDuplicate(user->Divxy[0],MAT_COPY_VALUES,&user->C[i]);
1158: MatDuplicate(user->Divxy[1],MAT_COPY_VALUES,&user->Cwork[i]);
1160: MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
1161: MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
1162: MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
1163: MatScale(user->C[i],user->ht);
1164: MatShift(user->C[i],1.0);
1165: }
1167: /* Solver options and tolerances */
1168: KSPCreate(PETSC_COMM_WORLD,&user->solver);
1169: KSPSetType(user->solver,KSPGMRES);
1170: KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec);
1171: KSPSetInitialGuessNonzero(user->solver,PETSC_FALSE); /* TODO: why is true slower? */
1172: KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500);
1173: /* KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500); */
1174: KSPGetPC(user->solver,&user->prec);
1175: PCSetType(user->prec,PCSHELL);
1177: PCShellSetApply(user->prec,StateMatBlockPrecMult);
1178: PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMultTranspose);
1179: PCShellSetContext(user->prec,user);
1181: /* Compute true state function yt given ut */
1182: VecCreate(PETSC_COMM_WORLD,&user->ytrue);
1183: VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt);
1184: VecSetFromOptions(user->ytrue);
1185: user->c_formed = PETSC_TRUE;
1186: VecCopy(user->utrue,user->u); /* Set u=utrue temporarily for StateMatInv */
1187: VecSet(user->ytrue,0.0); /* Initial guess */
1188: StateMatInvMult(user->Js,user->q,user->ytrue);
1189: VecCopy(user->ur,user->u); /* Reset u=ur */
1191: /* Initial guess y0 for state given u0 */
1192: StateMatInvMult(user->Js,user->q,user->y);
1194: /* Data discretization */
1195: MatCreate(PETSC_COMM_WORLD,&user->Q);
1196: MatSetSizes(user->Q,PETSC_DECIDE,PETSC_DECIDE,user->mx*user->mx,user->m);
1197: MatSetFromOptions(user->Q);
1198: MatMPIAIJSetPreallocation(user->Q,0,NULL,1,NULL);
1199: MatSeqAIJSetPreallocation(user->Q,1,NULL);
1201: MatGetOwnershipRange(user->Q,&istart,&iend);
1203: for (i=istart; i<iend; i++){
1204: j = i + user->m - user->mx*user->mx;
1205: MatSetValues(user->Q,1,&i,1,&j,&one,INSERT_VALUES);
1206: }
1208: MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1209: MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);
1211: MatTranspose(user->Q,MAT_INITIAL_MATRIX,&user->QT);
1213: VecDestroy(&XX);
1214: VecDestroy(&YY);
1215: VecDestroy(&XXwork);
1216: VecDestroy(&YYwork);
1217: VecDestroy(&yi);
1218: VecDestroy(&uxi);
1219: VecDestroy(&ui);
1220: VecDestroy(&bc);
1222: /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1223: KSPSetFromOptions(user->solver);
1224: return(0);
1225: }
1229: PetscErrorCode HyperbolicDestroy(AppCtx *user)
1230: {
1232: PetscInt i;
1235: MatDestroy(&user->Q);
1236: MatDestroy(&user->QT);
1237: MatDestroy(&user->Div);
1238: MatDestroy(&user->Divwork);
1239: MatDestroy(&user->Grad);
1240: MatDestroy(&user->L);
1241: MatDestroy(&user->LT);
1242: KSPDestroy(&user->solver);
1243: MatDestroy(&user->Js);
1244: MatDestroy(&user->Jd);
1245: MatDestroy(&user->JsBlockPrec);
1246: MatDestroy(&user->JsInv);
1247: MatDestroy(&user->JsBlock);
1248: MatDestroy(&user->Divxy[0]);
1249: MatDestroy(&user->Divxy[1]);
1250: MatDestroy(&user->Gradxy[0]);
1251: MatDestroy(&user->Gradxy[1]);
1252: MatDestroy(&user->M);
1253: for (i=0; i<user->nt; i++){
1254: MatDestroy(&user->C[i]);
1255: MatDestroy(&user->Cwork[i]);
1256: }
1257: PetscFree(user->C);
1258: PetscFree(user->Cwork);
1259: VecDestroy(&user->u);
1260: VecDestroy(&user->uwork);
1261: VecDestroy(&user->vwork);
1262: VecDestroy(&user->utrue);
1263: VecDestroy(&user->y);
1264: VecDestroy(&user->ywork);
1265: VecDestroy(&user->ytrue);
1266: VecDestroyVecs(user->nt,&user->yi);
1267: VecDestroyVecs(user->nt,&user->yiwork);
1268: VecDestroyVecs(user->nt,&user->ziwork);
1269: VecDestroyVecs(user->nt,&user->uxi);
1270: VecDestroyVecs(user->nt,&user->uyi);
1271: VecDestroyVecs(user->nt,&user->uxiwork);
1272: VecDestroyVecs(user->nt,&user->uyiwork);
1273: VecDestroyVecs(user->nt,&user->ui);
1274: VecDestroyVecs(user->nt,&user->uiwork);
1275: VecDestroy(&user->c);
1276: VecDestroy(&user->cwork);
1277: VecDestroy(&user->ur);
1278: VecDestroy(&user->q);
1279: VecDestroy(&user->d);
1280: VecDestroy(&user->dwork);
1281: VecDestroy(&user->lwork);
1282: VecDestroy(&user->js_diag);
1283: ISDestroy(&user->s_is);
1284: ISDestroy(&user->d_is);
1285: VecScatterDestroy(&user->state_scatter);
1286: VecScatterDestroy(&user->design_scatter);
1287: for (i=0; i<user->nt; i++){
1288: VecScatterDestroy(&user->uxi_scatter[i]);
1289: VecScatterDestroy(&user->uyi_scatter[i]);
1290: VecScatterDestroy(&user->ux_scatter[i]);
1291: VecScatterDestroy(&user->uy_scatter[i]);
1292: VecScatterDestroy(&user->ui_scatter[i]);
1293: VecScatterDestroy(&user->yi_scatter[i]);
1294: }
1295: PetscFree(user->uxi_scatter);
1296: PetscFree(user->uyi_scatter);
1297: PetscFree(user->ux_scatter);
1298: PetscFree(user->uy_scatter);
1299: PetscFree(user->ui_scatter);
1300: PetscFree(user->yi_scatter);
1301: return(0);
1302: }
1306: PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1307: {
1309: Vec X;
1310: PetscReal unorm,ynorm;
1311: AppCtx *user = (AppCtx*)ptr;
1314: TaoGetSolutionVector(tao,&X);
1315: Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1316: VecAXPY(user->ywork,-1.0,user->ytrue);
1317: VecAXPY(user->uwork,-1.0,user->utrue);
1318: VecNorm(user->uwork,NORM_2,&unorm);
1319: VecNorm(user->ywork,NORM_2,&ynorm);
1320: PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1321: return(0);
1322: }