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