Actual source code: hyperbolic.c
petsc-3.8.4 2018-03-24
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[]="";
113: int main(int argc, char **argv)
114: {
115: PetscErrorCode ierr;
116: Vec x,x0;
117: Tao tao;
118: AppCtx user;
119: IS is_allstate,is_alldesign;
120: PetscInt lo,hi,hi2,lo2,ksp_old;
121: PetscInt ntests = 1;
122: PetscInt i;
123: #if defined(PETSC_USE_LOG)
124: PetscLogStage stages[1];
125: #endif
127: PetscInitialize(&argc, &argv, (char*)0,help);if (ierr) return ierr;
128: user.mx = 32;
129: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);
130: PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);
131: user.nt = 16;
132: PetscOptionsInt("-nt","Number of time steps","",user.nt,&user.nt,NULL);
133: user.ndata = 64;
134: PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL);
135: user.alpha = 10.0;
136: PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);
137: user.T = 1.0/32.0;
138: PetscOptionsReal("-Tfinal","Final time","",user.T,&user.T,NULL);
140: user.m = user.mx*user.mx*user.nt; /* number of constraints */
141: user.n = user.mx*user.mx*3*user.nt; /* number of variables */
142: user.ht = user.T/user.nt; /* Time step */
143: user.gamma = user.T*user.ht / (user.mx*user.mx);
145: VecCreate(PETSC_COMM_WORLD,&user.u);
146: VecCreate(PETSC_COMM_WORLD,&user.y);
147: VecCreate(PETSC_COMM_WORLD,&user.c);
148: VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m);
149: VecSetSizes(user.y,PETSC_DECIDE,user.m);
150: VecSetSizes(user.c,PETSC_DECIDE,user.m);
151: VecSetFromOptions(user.u);
152: VecSetFromOptions(user.y);
153: VecSetFromOptions(user.c);
155: /* Create scatters for reduced spaces.
156: If the state vector y and design vector u are partitioned as
157: [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
158: then the solution vector x is organized as
159: [y_1; u_1; y_2; u_2; ...; y_np; u_np].
160: The index sets user.s_is and user.d_is correspond to the indices of the
161: state and design variables owned by the current processor.
162: */
163: VecCreate(PETSC_COMM_WORLD,&x);
165: VecGetOwnershipRange(user.y,&lo,&hi);
166: VecGetOwnershipRange(user.u,&lo2,&hi2);
168: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);
169: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is);
171: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);
172: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is);
174: VecSetSizes(x,hi-lo+hi2-lo2,user.n);
175: VecSetFromOptions(x);
177: VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter);
178: VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter);
179: ISDestroy(&is_alldesign);
180: ISDestroy(&is_allstate);
182: /* Create TAO solver and set desired solution method */
183: TaoCreate(PETSC_COMM_WORLD,&tao);
184: TaoSetType(tao,TAOLCL);
186: /* Set up initial vectors and matrices */
187: HyperbolicInitialize(&user);
189: Gather(x,user.y,user.state_scatter,user.u,user.design_scatter);
190: VecDuplicate(x,&x0);
191: VecCopy(x,x0);
193: /* Set solution vector with an initial guess */
194: TaoSetInitialVector(tao,x);
195: TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);
196: TaoSetGradientRoutine(tao, FormGradient, (void *)&user);
197: TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);
198: TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, (void *)&user);
199: TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);
200: TaoSetFromOptions(tao);
201: TaoSetStateDesignIS(tao,user.s_is,user.d_is);
203: /* SOLVE THE APPLICATION */
204: PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL);
205: PetscOptionsEnd();
206: PetscLogStageRegister("Trials",&stages[0]);
207: PetscLogStagePush(stages[0]);
208: user.ksp_its_initial = user.ksp_its;
209: ksp_old = user.ksp_its;
210: for (i=0; i<ntests; i++){
211: TaoSolve(tao);
212: PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old);
213: VecCopy(x0,x);
214: TaoSetInitialVector(tao,x);
215: }
216: PetscLogStagePop();
217: PetscBarrier((PetscObject)x);
218: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
219: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);
220: PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests);
221: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its);
222: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations per trial: ");
223: PetscPrintf(PETSC_COMM_WORLD,"%D\n",(user.ksp_its-user.ksp_its_initial)/ntests);
225: TaoDestroy(&tao);
226: VecDestroy(&x);
227: VecDestroy(&x0);
228: HyperbolicDestroy(&user);
229: PetscFinalize();
230: return ierr;
231: }
232: /* ------------------------------------------------------------------- */
233: /*
234: dwork = Qy - d
235: lwork = L*(u-ur).^2
236: f = 1/2 * (dwork.dork + alpha*y.lwork)
237: */
238: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
239: {
241: PetscReal d1=0,d2=0;
242: AppCtx *user = (AppCtx*)ptr;
245: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
246: MatMult(user->Q,user->y,user->dwork);
247: VecAXPY(user->dwork,-1.0,user->d);
248: VecDot(user->dwork,user->dwork,&d1);
250: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
251: VecPointwiseMult(user->uwork,user->uwork,user->uwork);
252: MatMult(user->L,user->uwork,user->lwork);
253: VecDot(user->y,user->lwork,&d2);
254: *f = 0.5 * (d1 + user->alpha*d2);
255: return(0);
256: }
258: /* ------------------------------------------------------------------- */
259: /*
260: state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2
261: design: g_d = alpha*(L'y).*(u-ur)
262: */
263: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
264: {
266: AppCtx *user = (AppCtx*)ptr;
269: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
270: MatMult(user->Q,user->y,user->dwork);
271: VecAXPY(user->dwork,-1.0,user->d);
273: MatMult(user->QT,user->dwork,user->ywork);
275: MatMult(user->LT,user->y,user->uwork);
276: VecWAXPY(user->vwork,-1.0,user->ur,user->u);
277: VecPointwiseMult(user->uwork,user->vwork,user->uwork);
278: VecScale(user->uwork,user->alpha);
280: VecPointwiseMult(user->vwork,user->vwork,user->vwork);
281: MatMult(user->L,user->vwork,user->lwork);
282: VecAXPY(user->ywork,0.5*user->alpha,user->lwork);
284: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
285: return(0);
286: }
288: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
289: {
291: PetscReal d1,d2;
292: AppCtx *user = (AppCtx*)ptr;
295: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
296: MatMult(user->Q,user->y,user->dwork);
297: VecAXPY(user->dwork,-1.0,user->d);
299: MatMult(user->QT,user->dwork,user->ywork);
301: VecDot(user->dwork,user->dwork,&d1);
303: MatMult(user->LT,user->y,user->uwork);
304: VecWAXPY(user->vwork,-1.0,user->ur,user->u);
305: VecPointwiseMult(user->uwork,user->vwork,user->uwork);
306: VecScale(user->uwork,user->alpha);
308: VecPointwiseMult(user->vwork,user->vwork,user->vwork);
309: MatMult(user->L,user->vwork,user->lwork);
310: VecAXPY(user->ywork,0.5*user->alpha,user->lwork);
312: VecDot(user->y,user->lwork,&d2);
314: *f = 0.5 * (d1 + user->alpha*d2);
315: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
316: return(0);
317: }
319: /* ------------------------------------------------------------------- */
320: /* A
321: MatShell object
322: */
323: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
324: {
326: PetscInt i;
327: AppCtx *user = (AppCtx*)ptr;
330: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
331: Scatter_yi(user->u,user->ui,user->ui_scatter,user->nt);
332: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
333: for (i=0; i<user->nt; i++){
334: MatCopy(user->Divxy[0],user->C[i],SAME_NONZERO_PATTERN);
335: MatCopy(user->Divxy[1],user->Cwork[i],SAME_NONZERO_PATTERN);
337: MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
338: MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
339: MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
340: MatScale(user->C[i],user->ht);
341: MatShift(user->C[i],1.0);
342: }
343: return(0);
344: }
346: /* ------------------------------------------------------------------- */
347: /* B */
348: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
349: {
351: AppCtx *user = (AppCtx*)ptr;
354: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
355: return(0);
356: }
358: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
359: {
361: PetscInt i;
362: AppCtx *user;
365: MatShellGetContext(J_shell,(void**)&user);
366: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
367: user->block_index = 0;
368: MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);
370: for (i=1; i<user->nt; i++){
371: user->block_index = i;
372: MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
373: MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
374: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
375: }
376: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
377: return(0);
378: }
380: PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
381: {
383: PetscInt i;
384: AppCtx *user;
387: MatShellGetContext(J_shell,(void**)&user);
388: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
390: for (i=0; i<user->nt-1; i++){
391: user->block_index = i;
392: MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
393: MatMult(user->M,user->yi[i+1],user->ziwork[i+1]);
394: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i+1]);
395: }
397: i = user->nt-1;
398: user->block_index = i;
399: MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
400: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
401: return(0);
402: }
404: PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
405: {
407: PetscInt i;
408: AppCtx *user;
411: MatShellGetContext(J_shell,(void**)&user);
412: i = user->block_index;
413: VecPointwiseMult(user->uxiwork[i],X,user->uxi[i]);
414: VecPointwiseMult(user->uyiwork[i],X,user->uyi[i]);
415: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
416: MatMult(user->Div,user->uiwork[i],Y);
417: VecAYPX(Y,user->ht,X);
418: return(0);
419: }
421: PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y)
422: {
424: PetscInt i;
425: AppCtx *user;
428: MatShellGetContext(J_shell,(void**)&user);
429: i = user->block_index;
430: MatMult(user->Grad,X,user->uiwork[i]);
431: Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
432: VecPointwiseMult(user->uxiwork[i],user->uxi[i],user->uxiwork[i]);
433: VecPointwiseMult(user->uyiwork[i],user->uyi[i],user->uyiwork[i]);
434: VecWAXPY(Y,1.0,user->uxiwork[i],user->uyiwork[i]);
435: VecAYPX(Y,user->ht,X);
436: return(0);
437: }
439: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
440: {
442: PetscInt i;
443: AppCtx *user;
446: MatShellGetContext(J_shell,(void**)&user);
447: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
448: Scatter_uxi_uyi(X,user->uxiwork,user->uxi_scatter,user->uyiwork,user->uyi_scatter,user->nt);
449: for (i=0; i<user->nt; i++){
450: VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
451: VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
452: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
453: MatMult(user->Div,user->uiwork[i],user->ziwork[i]);
454: VecScale(user->ziwork[i],user->ht);
455: }
456: Gather_yi(Y,user->ziwork,user->yi_scatter,user->nt);
457: return(0);
458: }
460: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
461: {
463: PetscInt i;
464: AppCtx *user;
467: MatShellGetContext(J_shell,(void**)&user);
468: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
469: Scatter_yi(X,user->yiwork,user->yi_scatter,user->nt);
470: for (i=0; i<user->nt; i++){
471: MatMult(user->Grad,user->yiwork[i],user->uiwork[i]);
472: Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
473: VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
474: VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
475: Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
476: VecScale(user->uiwork[i],user->ht);
477: }
478: Gather_yi(Y,user->uiwork,user->ui_scatter,user->nt);
479: return(0);
480: }
482: PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
483: {
485: PetscInt i;
486: AppCtx *user;
489: PCShellGetContext(PC_shell,(void**)&user);
490: i = user->block_index;
491: if (user->c_formed) {
492: MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
493: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
494: return(0);
495: }
497: PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
498: {
500: PetscInt i;
501: AppCtx *user;
504: PCShellGetContext(PC_shell,(void**)&user);
506: i = user->block_index;
507: if (user->c_formed) {
508: MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
509: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
510: return(0);
511: }
513: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
514: {
516: AppCtx *user;
517: PetscInt its,i;
520: MatShellGetContext(J_shell,(void**)&user);
522: if (Y == user->ytrue) {
523: /* First solve is done using true solution to set up problem */
524: KSPSetTolerances(user->solver,1e-4,1e-20,PETSC_DEFAULT,PETSC_DEFAULT);
525: } else {
526: KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
527: }
528: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
529: Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
530: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
532: user->block_index = 0;
533: KSPSolve(user->solver,user->yi[0],user->yiwork[0]);
535: KSPGetIterationNumber(user->solver,&its);
536: user->ksp_its = user->ksp_its + its;
537: for (i=1; i<user->nt; i++){
538: MatMult(user->M,user->yiwork[i-1],user->ziwork[i-1]);
539: VecAXPY(user->yi[i],1.0,user->ziwork[i-1]);
540: user->block_index = i;
541: KSPSolve(user->solver,user->yi[i],user->yiwork[i]);
543: KSPGetIterationNumber(user->solver,&its);
544: user->ksp_its = user->ksp_its + its;
545: }
546: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
547: return(0);
548: }
550: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
551: {
553: AppCtx *user;
554: PetscInt its,i;
557: MatShellGetContext(J_shell,(void**)&user);
559: Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
560: Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
561: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
563: i = user->nt - 1;
564: user->block_index = i;
565: KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);
567: KSPGetIterationNumber(user->solver,&its);
568: user->ksp_its = user->ksp_its + its;
570: for (i=user->nt-2; i>=0; i--){
571: MatMult(user->M,user->yiwork[i+1],user->ziwork[i+1]);
572: VecAXPY(user->yi[i],1.0,user->ziwork[i+1]);
573: user->block_index = i;
574: KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);
576: KSPGetIterationNumber(user->solver,&its);
577: user->ksp_its = user->ksp_its + its;
578: }
579: Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
580: return(0);
581: }
583: PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
584: {
586: AppCtx *user;
589: MatShellGetContext(J_shell,(void**)&user);
591: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_shell);
592: MatShellSetOperation(*new_shell,MATOP_MULT,(void(*)(void))StateMatMult);
593: MatShellSetOperation(*new_shell,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
594: MatShellSetOperation(*new_shell,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
595: MatShellSetOperation(*new_shell,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
596: return(0);
597: }
599: PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
600: {
602: AppCtx *user;
605: MatShellGetContext(J_shell,(void**)&user);
606: VecCopy(user->js_diag,X);
607: return(0);
608: }
610: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
611: {
612: /* con = Ay - q, A = [C(u1) 0 0 ... 0;
613: -M C(u2) 0 ... 0;
614: 0 -M C(u3) ... 0;
615: ... ;
616: 0 ... -M C(u_nt)]
617: C(u) = eye + ht*Div*[diag(u1); diag(u2)] */
619: PetscInt i;
620: AppCtx *user = (AppCtx*)ptr;
623: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
624: Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
625: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
627: user->block_index = 0;
628: MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);
630: for (i=1; i<user->nt; i++){
631: user->block_index = i;
632: MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
633: MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
634: VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
635: }
637: Gather_yi(C,user->yiwork,user->yi_scatter,user->nt);
638: VecAXPY(C,-1.0,user->q);
640: return(0);
641: }
644: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
645: {
649: VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
650: VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
651: VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
652: VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
653: return(0);
654: }
656: PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
657: {
659: PetscInt i;
662: for (i=0; i<nt; i++){
663: VecScatterBegin(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
664: VecScatterEnd(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
665: VecScatterBegin(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
666: VecScatterEnd(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
667: }
668: return(0);
669: }
671: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
672: {
676: VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
677: VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
678: VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
679: VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
680: return(0);
681: }
683: PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
684: {
686: PetscInt i;
689: for (i=0; i<nt; i++){
690: VecScatterBegin(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
691: VecScatterEnd(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
692: VecScatterBegin(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
693: VecScatterEnd(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
694: }
695: return(0);
696: }
698: PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
699: {
701: PetscInt i;
704: for (i=0; i<nt; i++){
705: VecScatterBegin(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
706: VecScatterEnd(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
707: }
708: return(0);
709: }
711: PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
712: {
714: PetscInt i;
717: for (i=0; i<nt; i++){
718: VecScatterBegin(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
719: VecScatterEnd(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
720: }
721: return(0);
722: }
724: PetscErrorCode HyperbolicInitialize(AppCtx *user)
725: {
727: PetscInt n,i,j,linear_index,istart,iend,iblock,lo,hi;
728: Vec XX,YY,XXwork,YYwork,yi,uxi,ui,bc;
729: PetscReal h,sum;
730: PetscScalar hinv,neg_hinv,quarter=0.25,one=1.0,half_hinv,neg_half_hinv;
731: PetscScalar vx,vy,zero=0.0;
732: IS is_from_y,is_to_yi,is_from_u,is_to_uxi,is_to_uyi;
735: user->jformed = PETSC_FALSE;
736: user->c_formed = PETSC_FALSE;
738: user->ksp_its = 0;
739: user->ksp_its_initial = 0;
741: n = user->mx * user->mx;
743: h = 1.0/user->mx;
744: hinv = user->mx;
745: neg_hinv = -hinv;
746: half_hinv = hinv / 2.0;
747: neg_half_hinv = neg_hinv / 2.0;
749: /* Generate Grad matrix */
750: MatCreate(PETSC_COMM_WORLD,&user->Grad);
751: MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,2*n,n);
752: MatSetFromOptions(user->Grad);
753: MatMPIAIJSetPreallocation(user->Grad,3,NULL,3,NULL);
754: MatSeqAIJSetPreallocation(user->Grad,3,NULL);
755: MatGetOwnershipRange(user->Grad,&istart,&iend);
757: for (i=istart; i<iend; i++){
758: if (i<n){
759: iblock = i / user->mx;
760: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
761: MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
762: j = iblock*user->mx + ((i+1) % user->mx);
763: MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
764: }
765: if (i>=n){
766: j = (i - user->mx) % n;
767: MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
768: j = (j + 2*user->mx) % n;
769: MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
770: }
771: }
773: MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
774: MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);
776: MatCreate(PETSC_COMM_WORLD,&user->Gradxy[0]);
777: MatSetSizes(user->Gradxy[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
778: MatSetFromOptions(user->Gradxy[0]);
779: MatMPIAIJSetPreallocation(user->Gradxy[0],3,NULL,3,NULL);
780: MatSeqAIJSetPreallocation(user->Gradxy[0],3,NULL);
781: MatGetOwnershipRange(user->Gradxy[0],&istart,&iend);
783: for (i=istart; i<iend; i++){
784: iblock = i / user->mx;
785: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
786: MatSetValues(user->Gradxy[0],1,&i,1,&j,&half_hinv,INSERT_VALUES);
787: j = iblock*user->mx + ((i+1) % user->mx);
788: MatSetValues(user->Gradxy[0],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
789: MatSetValues(user->Gradxy[0],1,&i,1,&i,&zero,INSERT_VALUES);
790: }
791: MatAssemblyBegin(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
792: MatAssemblyEnd(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
794: MatCreate(PETSC_COMM_WORLD,&user->Gradxy[1]);
795: MatSetSizes(user->Gradxy[1],PETSC_DECIDE,PETSC_DECIDE,n,n);
796: MatSetFromOptions(user->Gradxy[1]);
797: MatMPIAIJSetPreallocation(user->Gradxy[1],3,NULL,3,NULL);
798: MatSeqAIJSetPreallocation(user->Gradxy[1],3,NULL);
799: MatGetOwnershipRange(user->Gradxy[1],&istart,&iend);
801: for (i=istart; i<iend; i++){
802: j = (i + n - user->mx) % n;
803: MatSetValues(user->Gradxy[1],1,&i,1,&j,&half_hinv,INSERT_VALUES);
804: j = (j + 2*user->mx) % n;
805: MatSetValues(user->Gradxy[1],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
806: MatSetValues(user->Gradxy[1],1,&i,1,&i,&zero,INSERT_VALUES);
807: }
808: MatAssemblyBegin(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
809: MatAssemblyEnd(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
811: /* Generate Div matrix */
812: MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div);
813: MatTranspose(user->Gradxy[0],MAT_INITIAL_MATRIX,&user->Divxy[0]);
814: MatTranspose(user->Gradxy[1],MAT_INITIAL_MATRIX,&user->Divxy[1]);
816: /* Off-diagonal averaging matrix */
817: MatCreate(PETSC_COMM_WORLD,&user->M);
818: MatSetSizes(user->M,PETSC_DECIDE,PETSC_DECIDE,n,n);
819: MatSetFromOptions(user->M);
820: MatMPIAIJSetPreallocation(user->M,4,NULL,4,NULL);
821: MatSeqAIJSetPreallocation(user->M,4,NULL);
822: MatGetOwnershipRange(user->M,&istart,&iend);
824: for (i=istart; i<iend; i++){
825: /* kron(Id,Av) */
826: iblock = i / user->mx;
827: j = iblock*user->mx + ((i+user->mx-1) % user->mx);
828: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
829: j = iblock*user->mx + ((i+1) % user->mx);
830: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
832: /* kron(Av,Id) */
833: j = (i + user->mx) % n;
834: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
835: j = (i + n - user->mx) % n;
836: MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
837: }
838: MatAssemblyBegin(user->M,MAT_FINAL_ASSEMBLY);
839: MatAssemblyEnd(user->M,MAT_FINAL_ASSEMBLY);
841: /* Generate 2D grid */
842: VecCreate(PETSC_COMM_WORLD,&XX);
843: VecCreate(PETSC_COMM_WORLD,&user->q);
844: VecSetSizes(XX,PETSC_DECIDE,n);
845: VecSetSizes(user->q,PETSC_DECIDE,n*user->nt);
846: VecSetFromOptions(XX);
847: VecSetFromOptions(user->q);
849: VecDuplicate(XX,&YY);
850: VecDuplicate(XX,&XXwork);
851: VecDuplicate(XX,&YYwork);
852: VecDuplicate(XX,&user->d);
853: VecDuplicate(XX,&user->dwork);
855: VecGetOwnershipRange(XX,&istart,&iend);
856: for (linear_index=istart; linear_index<iend; linear_index++){
857: i = linear_index % user->mx;
858: j = (linear_index-i)/user->mx;
859: vx = h*(i+0.5);
860: vy = h*(j+0.5);
861: VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
862: VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
863: }
865: VecAssemblyBegin(XX);
866: VecAssemblyEnd(XX);
867: VecAssemblyBegin(YY);
868: VecAssemblyEnd(YY);
870: /* Compute final density function yT
871: yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
872: yT = yT / (h^2*sum(yT)) */
873: VecCopy(XX,XXwork);
874: VecCopy(YY,YYwork);
876: VecShift(XXwork,-0.25);
877: VecShift(YYwork,-0.25);
879: VecPointwiseMult(XXwork,XXwork,XXwork);
880: VecPointwiseMult(YYwork,YYwork,YYwork);
882: VecCopy(XXwork,user->dwork);
883: VecAXPY(user->dwork,1.0,YYwork);
884: VecScale(user->dwork,-30.0);
885: VecExp(user->dwork);
886: VecCopy(user->dwork,user->d);
888: VecCopy(XX,XXwork);
889: VecCopy(YY,YYwork);
891: VecShift(XXwork,-0.75);
892: VecShift(YYwork,-0.75);
894: VecPointwiseMult(XXwork,XXwork,XXwork);
895: VecPointwiseMult(YYwork,YYwork,YYwork);
897: VecCopy(XXwork,user->dwork);
898: VecAXPY(user->dwork,1.0,YYwork);
899: VecScale(user->dwork,-30.0);
900: VecExp(user->dwork);
902: VecAXPY(user->d,1.0,user->dwork);
903: VecShift(user->d,1.0);
904: VecSum(user->d,&sum);
905: VecScale(user->d,1.0/(h*h*sum));
907: /* Initial conditions of forward problem */
908: VecDuplicate(XX,&bc);
909: VecCopy(XX,XXwork);
910: VecCopy(YY,YYwork);
912: VecShift(XXwork,-0.5);
913: VecShift(YYwork,-0.5);
915: VecPointwiseMult(XXwork,XXwork,XXwork);
916: VecPointwiseMult(YYwork,YYwork,YYwork);
918: VecWAXPY(bc,1.0,XXwork,YYwork);
919: VecScale(bc,-50.0);
920: VecExp(bc);
921: VecShift(bc,1.0);
922: VecSum(bc,&sum);
923: VecScale(bc,1.0/(h*h*sum));
925: /* Create scatter from y to y_1,y_2,...,y_nt */
926: /* TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
927: PetscMalloc1(user->nt*user->mx*user->mx,&user->yi_scatter);
928: VecCreate(PETSC_COMM_WORLD,&yi);
929: VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx);
930: VecSetFromOptions(yi);
931: VecDuplicateVecs(yi,user->nt,&user->yi);
932: VecDuplicateVecs(yi,user->nt,&user->yiwork);
933: VecDuplicateVecs(yi,user->nt,&user->ziwork);
934: for (i=0; i<user->nt; i++){
935: VecGetOwnershipRange(user->yi[i],&lo,&hi);
936: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_yi);
937: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+i*user->mx*user->mx,1,&is_from_y);
938: VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i]);
939: ISDestroy(&is_to_yi);
940: ISDestroy(&is_from_y);
941: }
943: /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
944: /* TODO: reorder for better parallelism */
945: PetscMalloc1(user->nt*user->mx*user->mx,&user->uxi_scatter);
946: PetscMalloc1(user->nt*user->mx*user->mx,&user->uyi_scatter);
947: PetscMalloc1(user->nt*user->mx*user->mx,&user->ux_scatter);
948: PetscMalloc1(user->nt*user->mx*user->mx,&user->uy_scatter);
949: PetscMalloc1(2*user->nt*user->mx*user->mx,&user->ui_scatter);
950: VecCreate(PETSC_COMM_WORLD,&uxi);
951: VecCreate(PETSC_COMM_WORLD,&ui);
952: VecSetSizes(uxi,PETSC_DECIDE,user->mx*user->mx);
953: VecSetSizes(ui,PETSC_DECIDE,2*user->mx*user->mx);
954: VecSetFromOptions(uxi);
955: VecSetFromOptions(ui);
956: VecDuplicateVecs(uxi,user->nt,&user->uxi);
957: VecDuplicateVecs(uxi,user->nt,&user->uyi);
958: VecDuplicateVecs(uxi,user->nt,&user->uxiwork);
959: VecDuplicateVecs(uxi,user->nt,&user->uyiwork);
960: VecDuplicateVecs(ui,user->nt,&user->ui);
961: VecDuplicateVecs(ui,user->nt,&user->uiwork);
962: for (i=0; i<user->nt; i++){
963: VecGetOwnershipRange(user->uxi[i],&lo,&hi);
964: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
965: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
966: VecScatterCreate(user->u,is_from_u,user->uxi[i],is_to_uxi,&user->uxi_scatter[i]);
968: ISDestroy(&is_to_uxi);
969: ISDestroy(&is_from_u);
971: VecGetOwnershipRange(user->uyi[i],&lo,&hi);
972: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
973: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+(2*i+1)*user->mx*user->mx,1,&is_from_u);
974: VecScatterCreate(user->u,is_from_u,user->uyi[i],is_to_uyi,&user->uyi_scatter[i]);
976: ISDestroy(&is_to_uyi);
977: ISDestroy(&is_from_u);
979: VecGetOwnershipRange(user->uxi[i],&lo,&hi);
980: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
981: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_from_u);
982: VecScatterCreate(user->ui[i],is_from_u,user->uxi[i],is_to_uxi,&user->ux_scatter[i]);
984: ISDestroy(&is_to_uxi);
985: ISDestroy(&is_from_u);
987: VecGetOwnershipRange(user->uyi[i],&lo,&hi);
988: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
989: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+user->mx*user->mx,1,&is_from_u);
990: VecScatterCreate(user->ui[i],is_from_u,user->uyi[i],is_to_uyi,&user->uy_scatter[i]);
992: ISDestroy(&is_to_uyi);
993: ISDestroy(&is_from_u);
995: VecGetOwnershipRange(user->ui[i],&lo,&hi);
996: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
997: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
998: VecScatterCreate(user->u,is_from_u,user->ui[i],is_to_uxi,&user->ui_scatter[i]);
1000: ISDestroy(&is_to_uxi);
1001: ISDestroy(&is_from_u);
1002: }
1004: /* RHS of forward problem */
1005: MatMult(user->M,bc,user->yiwork[0]);
1006: for (i=1; i<user->nt; i++){
1007: VecSet(user->yiwork[i],0.0);
1008: }
1009: Gather_yi(user->q,user->yiwork,user->yi_scatter,user->nt);
1011: /* Compute true velocity field utrue */
1012: VecDuplicate(user->u,&user->utrue);
1013: for (i=0; i<user->nt; i++){
1014: VecCopy(YY,user->uxi[i]);
1015: VecScale(user->uxi[i],150.0*i*user->ht);
1016: VecCopy(XX,user->uyi[i]);
1017: VecShift(user->uyi[i],-10.0);
1018: VecScale(user->uyi[i],15.0*i*user->ht);
1019: }
1020: Gather_uxi_uyi(user->utrue,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1022: /* Initial guess and reference model */
1023: VecDuplicate(user->utrue,&user->ur);
1024: for (i=0; i<user->nt; i++){
1025: VecCopy(XX,user->uxi[i]);
1026: VecShift(user->uxi[i],i*user->ht);
1027: VecCopy(YY,user->uyi[i]);
1028: VecShift(user->uyi[i],-i*user->ht);
1029: }
1030: Gather_uxi_uyi(user->ur,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1032: /* Generate regularization matrix L */
1033: MatCreate(PETSC_COMM_WORLD,&user->LT);
1034: MatSetSizes(user->LT,PETSC_DECIDE,PETSC_DECIDE,2*n*user->nt,n*user->nt);
1035: MatSetFromOptions(user->LT);
1036: MatMPIAIJSetPreallocation(user->LT,1,NULL,1,NULL);
1037: MatSeqAIJSetPreallocation(user->LT,1,NULL);
1038: MatGetOwnershipRange(user->LT,&istart,&iend);
1040: for (i=istart; i<iend; i++){
1041: iblock = (i+n) / (2*n);
1042: j = i - iblock*n;
1043: MatSetValues(user->LT,1,&i,1,&j,&user->gamma,INSERT_VALUES);
1044: }
1046: MatAssemblyBegin(user->LT,MAT_FINAL_ASSEMBLY);
1047: MatAssemblyEnd(user->LT,MAT_FINAL_ASSEMBLY);
1049: MatTranspose(user->LT,MAT_INITIAL_MATRIX,&user->L);
1051: /* Build work vectors and matrices */
1052: VecCreate(PETSC_COMM_WORLD,&user->lwork);
1053: VecSetType(user->lwork,VECMPI);
1054: VecSetSizes(user->lwork,PETSC_DECIDE,user->m);
1055: VecSetFromOptions(user->lwork);
1057: MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);
1059: VecDuplicate(user->y,&user->ywork);
1060: VecDuplicate(user->u,&user->uwork);
1061: VecDuplicate(user->u,&user->vwork);
1062: VecDuplicate(user->u,&user->js_diag);
1063: VecDuplicate(user->c,&user->cwork);
1065: /* Create matrix-free shell user->Js for computing A*x */
1066: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->Js);
1067: MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1068: MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
1069: MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
1070: MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
1072: /* Diagonal blocks of user->Js */
1073: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlock);
1074: MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult);
1075: MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockMultTranspose);
1077: /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
1078: D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
1079: This is an SOR preconditioner for user->JsBlock. */
1080: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlockPrec);
1081: MatShellSetOperation(user->JsBlockPrec,MATOP_MULT,(void(*)(void))StateMatBlockPrecMult);
1082: MatShellSetOperation(user->JsBlockPrec,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockPrecMultTranspose);
1084: /* Create a matrix-free shell user->Jd for computing B*x */
1085: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->n-user->m,user,&user->Jd);
1086: MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);
1087: MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);
1089: /* User-defined routines for computing user->Js\x and user->Js^T\x*/
1090: MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsInv);
1091: MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateMatInvMult);
1092: MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatInvTransposeMult);
1094: /* Build matrices for SOR preconditioner */
1095: Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1096: MatShift(user->Divxy[0],0.0); /* Force C[i] and Divxy[0] to share same nonzero pattern */
1097: MatAXPY(user->Divxy[0],0.0,user->Divxy[1],DIFFERENT_NONZERO_PATTERN);
1098: PetscMalloc1(5*n,&user->C);
1099: PetscMalloc1(2*n,&user->Cwork);
1100: for (i=0; i<user->nt; i++){
1101: MatDuplicate(user->Divxy[0],MAT_COPY_VALUES,&user->C[i]);
1102: MatDuplicate(user->Divxy[1],MAT_COPY_VALUES,&user->Cwork[i]);
1104: MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
1105: MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
1106: MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
1107: MatScale(user->C[i],user->ht);
1108: MatShift(user->C[i],1.0);
1109: }
1111: /* Solver options and tolerances */
1112: KSPCreate(PETSC_COMM_WORLD,&user->solver);
1113: KSPSetType(user->solver,KSPGMRES);
1114: KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec);
1115: KSPSetInitialGuessNonzero(user->solver,PETSC_FALSE); /* TODO: why is true slower? */
1116: KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500);
1117: /* KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500); */
1118: KSPGetPC(user->solver,&user->prec);
1119: PCSetType(user->prec,PCSHELL);
1121: PCShellSetApply(user->prec,StateMatBlockPrecMult);
1122: PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMultTranspose);
1123: PCShellSetContext(user->prec,user);
1125: /* Compute true state function yt given ut */
1126: VecCreate(PETSC_COMM_WORLD,&user->ytrue);
1127: VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt);
1128: VecSetFromOptions(user->ytrue);
1129: user->c_formed = PETSC_TRUE;
1130: VecCopy(user->utrue,user->u); /* Set u=utrue temporarily for StateMatInv */
1131: VecSet(user->ytrue,0.0); /* Initial guess */
1132: StateMatInvMult(user->Js,user->q,user->ytrue);
1133: VecCopy(user->ur,user->u); /* Reset u=ur */
1135: /* Initial guess y0 for state given u0 */
1136: StateMatInvMult(user->Js,user->q,user->y);
1138: /* Data discretization */
1139: MatCreate(PETSC_COMM_WORLD,&user->Q);
1140: MatSetSizes(user->Q,PETSC_DECIDE,PETSC_DECIDE,user->mx*user->mx,user->m);
1141: MatSetFromOptions(user->Q);
1142: MatMPIAIJSetPreallocation(user->Q,0,NULL,1,NULL);
1143: MatSeqAIJSetPreallocation(user->Q,1,NULL);
1145: MatGetOwnershipRange(user->Q,&istart,&iend);
1147: for (i=istart; i<iend; i++){
1148: j = i + user->m - user->mx*user->mx;
1149: MatSetValues(user->Q,1,&i,1,&j,&one,INSERT_VALUES);
1150: }
1152: MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1153: MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);
1155: MatTranspose(user->Q,MAT_INITIAL_MATRIX,&user->QT);
1157: VecDestroy(&XX);
1158: VecDestroy(&YY);
1159: VecDestroy(&XXwork);
1160: VecDestroy(&YYwork);
1161: VecDestroy(&yi);
1162: VecDestroy(&uxi);
1163: VecDestroy(&ui);
1164: VecDestroy(&bc);
1166: /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1167: KSPSetFromOptions(user->solver);
1168: return(0);
1169: }
1171: PetscErrorCode HyperbolicDestroy(AppCtx *user)
1172: {
1174: PetscInt i;
1177: MatDestroy(&user->Q);
1178: MatDestroy(&user->QT);
1179: MatDestroy(&user->Div);
1180: MatDestroy(&user->Divwork);
1181: MatDestroy(&user->Grad);
1182: MatDestroy(&user->L);
1183: MatDestroy(&user->LT);
1184: KSPDestroy(&user->solver);
1185: MatDestroy(&user->Js);
1186: MatDestroy(&user->Jd);
1187: MatDestroy(&user->JsBlockPrec);
1188: MatDestroy(&user->JsInv);
1189: MatDestroy(&user->JsBlock);
1190: MatDestroy(&user->Divxy[0]);
1191: MatDestroy(&user->Divxy[1]);
1192: MatDestroy(&user->Gradxy[0]);
1193: MatDestroy(&user->Gradxy[1]);
1194: MatDestroy(&user->M);
1195: for (i=0; i<user->nt; i++){
1196: MatDestroy(&user->C[i]);
1197: MatDestroy(&user->Cwork[i]);
1198: }
1199: PetscFree(user->C);
1200: PetscFree(user->Cwork);
1201: VecDestroy(&user->u);
1202: VecDestroy(&user->uwork);
1203: VecDestroy(&user->vwork);
1204: VecDestroy(&user->utrue);
1205: VecDestroy(&user->y);
1206: VecDestroy(&user->ywork);
1207: VecDestroy(&user->ytrue);
1208: VecDestroyVecs(user->nt,&user->yi);
1209: VecDestroyVecs(user->nt,&user->yiwork);
1210: VecDestroyVecs(user->nt,&user->ziwork);
1211: VecDestroyVecs(user->nt,&user->uxi);
1212: VecDestroyVecs(user->nt,&user->uyi);
1213: VecDestroyVecs(user->nt,&user->uxiwork);
1214: VecDestroyVecs(user->nt,&user->uyiwork);
1215: VecDestroyVecs(user->nt,&user->ui);
1216: VecDestroyVecs(user->nt,&user->uiwork);
1217: VecDestroy(&user->c);
1218: VecDestroy(&user->cwork);
1219: VecDestroy(&user->ur);
1220: VecDestroy(&user->q);
1221: VecDestroy(&user->d);
1222: VecDestroy(&user->dwork);
1223: VecDestroy(&user->lwork);
1224: VecDestroy(&user->js_diag);
1225: ISDestroy(&user->s_is);
1226: ISDestroy(&user->d_is);
1227: VecScatterDestroy(&user->state_scatter);
1228: VecScatterDestroy(&user->design_scatter);
1229: for (i=0; i<user->nt; i++){
1230: VecScatterDestroy(&user->uxi_scatter[i]);
1231: VecScatterDestroy(&user->uyi_scatter[i]);
1232: VecScatterDestroy(&user->ux_scatter[i]);
1233: VecScatterDestroy(&user->uy_scatter[i]);
1234: VecScatterDestroy(&user->ui_scatter[i]);
1235: VecScatterDestroy(&user->yi_scatter[i]);
1236: }
1237: PetscFree(user->uxi_scatter);
1238: PetscFree(user->uyi_scatter);
1239: PetscFree(user->ux_scatter);
1240: PetscFree(user->uy_scatter);
1241: PetscFree(user->ui_scatter);
1242: PetscFree(user->yi_scatter);
1243: return(0);
1244: }
1246: PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1247: {
1249: Vec X;
1250: PetscReal unorm,ynorm;
1251: AppCtx *user = (AppCtx*)ptr;
1254: TaoGetSolutionVector(tao,&X);
1255: Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1256: VecAXPY(user->ywork,-1.0,user->ytrue);
1257: VecAXPY(user->uwork,-1.0,user->utrue);
1258: VecNorm(user->uwork,NORM_2,&unorm);
1259: VecNorm(user->ywork,NORM_2,&ynorm);
1260: PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1261: return(0);
1262: }