Actual source code: hyperbolic.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #include <petsctao.h>

  3: /*T
  4:    Concepts: TAO^Solving a system of nonlinear equations, nonlinear least squares
  5:    Routines: TaoCreate();
  6:    Routines: TaoSetType();
  7:    Routines: TaoSetInitialVector();
  8:    Routines: TaoSetObjectiveRoutine();
  9:    Routines: TaoSetGradientRoutine();
 10:    Routines: TaoSetConstraintsRoutine();
 11:    Routines: TaoSetJacobianStateRoutine();
 12:    Routines: TaoSetJacobianDesignRoutine();
 13:    Routines: TaoSetStateDesignIS();
 14:    Routines: TaoSetFromOptions();
 15:    Routines: TaoSetHistory(); TaoGetHistory();
 16:    Routines: TaoSolve();
 17:    Routines: TaoGetConvergedReason(); TaoDestroy();
 18:    Processors: 1
 19: T*/

 21: typedef struct {
 22:   PetscInt n; /*  Number of variables */
 23:   PetscInt m; /*  Number of constraints */
 24:   PetscInt mx; /*  grid points in each direction */
 25:   PetscInt nt; /*  Number of time steps */
 26:   PetscInt ndata; /*  Number of data points per sample */
 27:   IS       s_is;
 28:   IS       d_is;
 29:   VecScatter state_scatter;
 30:   VecScatter design_scatter;
 31:   VecScatter *uxi_scatter,*uyi_scatter,*ux_scatter,*uy_scatter,*ui_scatter;
 32:   VecScatter *yi_scatter;

 34:   Mat       Js,Jd,JsBlockPrec,JsInv,JsBlock;
 35:   PetscBool jformed,c_formed;

 37:   PetscReal alpha; /*  Regularization parameter */
 38:   PetscReal gamma;
 39:   PetscReal ht; /*  Time step */
 40:   PetscReal T; /*  Final time */
 41:   Mat Q,QT;
 42:   Mat L,LT;
 43:   Mat Div,Divwork,Divxy[2];
 44:   Mat Grad,Gradxy[2];
 45:   Mat M;
 46:   Mat *C,*Cwork;
 47:   /* Mat Hs,Hd,Hsd; */
 48:   Vec q;
 49:   Vec ur; /*  reference */

 51:   Vec d;
 52:   Vec dwork;

 54:   Vec y; /*  state variables */
 55:   Vec ywork;
 56:   Vec ytrue;
 57:   Vec *yi,*yiwork,*ziwork;
 58:   Vec *uxi,*uyi,*uxiwork,*uyiwork,*ui,*uiwork;

 60:   Vec u; /*  design variables */
 61:   Vec uwork,vwork;
 62:   Vec utrue;

 64:   Vec js_diag;

 66:   Vec c; /*  constraint vector */
 67:   Vec cwork;

 69:   Vec lwork;

 71:   KSP      solver;
 72:   PC       prec;
 73:   PetscInt block_index;

 75:   PetscInt ksp_its;
 76:   PetscInt ksp_its_initial;
 77: } AppCtx;


 80: PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*);
 81: PetscErrorCode FormGradient(Tao, Vec, Vec, void*);
 82: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*);
 83: PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*);
 84: PetscErrorCode FormJacobianDesign(Tao, Vec, Mat,void*);
 85: PetscErrorCode FormConstraints(Tao, Vec, Vec, void*);
 86: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*);
 87: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
 88: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
 89: PetscErrorCode HyperbolicInitialize(AppCtx *user);
 90: PetscErrorCode HyperbolicDestroy(AppCtx *user);
 91: PetscErrorCode HyperbolicMonitor(Tao, void*);

 93: PetscErrorCode StateMatMult(Mat,Vec,Vec);
 94: PetscErrorCode StateMatBlockMult(Mat,Vec,Vec);
 95: PetscErrorCode StateMatBlockMultTranspose(Mat,Vec,Vec);
 96: PetscErrorCode StateMatMultTranspose(Mat,Vec,Vec);
 97: PetscErrorCode StateMatGetDiagonal(Mat,Vec);
 98: PetscErrorCode StateMatDuplicate(Mat,MatDuplicateOption,Mat*);
 99: PetscErrorCode StateMatInvMult(Mat,Vec,Vec);
100: PetscErrorCode StateMatInvTransposeMult(Mat,Vec,Vec);
101: PetscErrorCode StateMatBlockPrecMult(PC,Vec,Vec);

103: PetscErrorCode DesignMatMult(Mat,Vec,Vec);
104: PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);

106: PetscErrorCode Scatter_yi(Vec,Vec*,VecScatter*,PetscInt); /*  y to y1,y2,...,y_nt */
107: PetscErrorCode Gather_yi(Vec,Vec*,VecScatter*,PetscInt);
108: PetscErrorCode Scatter_uxi_uyi(Vec,Vec*,VecScatter*,Vec*,VecScatter*,PetscInt); /*  u to ux_1,uy_1,ux_2,uy_2,...,u */
109: PetscErrorCode Gather_uxi_uyi(Vec,Vec*,VecScatter*,Vec*,VecScatter*,PetscInt);

111: static  char help[]="";

115: int main(int argc, char **argv)
116: {
117:   PetscErrorCode     ierr;
118:   Vec                x,x0;
119:   Tao                tao;
120:   TaoConvergedReason reason;
121:   AppCtx             user;
122:   IS                 is_allstate,is_alldesign;
123:   PetscInt           lo,hi,hi2,lo2,ksp_old;
124:   PetscInt           ntests = 1;
125:   PetscInt           i;
126:   int                stages[1];

128:   PetscInitialize(&argc, &argv, (char*)0,help);

130:   user.mx = 32;
131:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);
132:   PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);
133:   user.nt = 16;
134:   PetscOptionsInt("-nt","Number of time steps","",user.nt,&user.nt,NULL);
135:   user.ndata = 64;
136:   PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL);
137:   user.alpha = 10.0;
138:   PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);
139:   user.T = 1.0/32.0;
140:   PetscOptionsReal("-Tfinal","Final time","",user.T,&user.T,NULL);

142:   user.m = user.mx*user.mx*user.nt; /*  number of constraints */
143:   user.n = user.mx*user.mx*3*user.nt; /*  number of variables */
144:   user.ht = user.T/user.nt; /*  Time step */
145:   user.gamma = user.T*user.ht / (user.mx*user.mx);

147:   VecCreate(PETSC_COMM_WORLD,&user.u);
148:   VecCreate(PETSC_COMM_WORLD,&user.y);
149:   VecCreate(PETSC_COMM_WORLD,&user.c);
150:   VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m);
151:   VecSetSizes(user.y,PETSC_DECIDE,user.m);
152:   VecSetSizes(user.c,PETSC_DECIDE,user.m);
153:   VecSetFromOptions(user.u);
154:   VecSetFromOptions(user.y);
155:   VecSetFromOptions(user.c);

157:   /* Create scatters for reduced spaces.
158:      If the state vector y and design vector u are partitioned as
159:      [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
160:      then the solution vector x is organized as
161:      [y_1; u_1; y_2; u_2; ...; y_np; u_np].
162:      The index sets user.s_is and user.d_is correspond to the indices of the
163:      state and design variables owned by the current processor.
164:   */
165:   VecCreate(PETSC_COMM_WORLD,&x);

167:   VecGetOwnershipRange(user.y,&lo,&hi);
168:   VecGetOwnershipRange(user.u,&lo2,&hi2);

170:   ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);
171:   ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is);

173:   ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);
174:   ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is);

176:   VecSetSizes(x,hi-lo+hi2-lo2,user.n);
177:   VecSetFromOptions(x);

179:   VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter);
180:   VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter);
181:   ISDestroy(&is_alldesign);
182:   ISDestroy(&is_allstate);

184:   /* Create TAO solver and set desired solution method */
185:   TaoCreate(PETSC_COMM_WORLD,&tao);
186:   TaoSetType(tao,TAOLCL);

188:   /* Set up initial vectors and matrices */
189:   HyperbolicInitialize(&user);

191:   Gather(x,user.y,user.state_scatter,user.u,user.design_scatter);
192:   VecDuplicate(x,&x0);
193:   VecCopy(x,x0);

195:   /* Set solution vector with an initial guess */
196:   TaoSetInitialVector(tao,x);
197:   TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);
198:   TaoSetGradientRoutine(tao, FormGradient, (void *)&user);
199:   TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);

201:   TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, (void *)&user);

203:   TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);

205:   TaoSetFromOptions(tao);
206:   TaoSetStateDesignIS(tao,user.s_is,user.d_is);

208:   /* SOLVE THE APPLICATION */
209:   PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL);
210:   PetscOptionsEnd();
211:   PetscLogStageRegister("Trials",&stages[0]);
212:   PetscLogStagePush(stages[0]);
213:   user.ksp_its_initial = user.ksp_its;
214:   ksp_old = user.ksp_its;
215:   for (i=0; i<ntests; i++){
216:     TaoSolve(tao);
217:     PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old);
218:     VecCopy(x0,x);
219:     TaoSetInitialVector(tao,x);
220:   }
221:   PetscLogStagePop();
222:   PetscBarrier((PetscObject)x);
223:   PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
224:   PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);
225:   PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests);
226:   PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its);
227:   PetscPrintf(PETSC_COMM_WORLD,"KSP iterations per trial: ");
228:   PetscPrintf(PETSC_COMM_WORLD,"%D\n",(user.ksp_its-user.ksp_its_initial)/ntests);

230:   TaoGetConvergedReason(tao,&reason);

232:   if (reason < 0) {
233:     PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge.\n");
234:   } else {
235:     PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);
236:   }

238:   TaoDestroy(&tao);
239:   VecDestroy(&x);
240:   VecDestroy(&x0);
241:   HyperbolicDestroy(&user);
242:   PetscFinalize();
243:   return 0;
244: }
245: /* ------------------------------------------------------------------- */
248: /*
249:    dwork = Qy - d
250:    lwork = L*(u-ur).^2
251:    f = 1/2 * (dwork.dork + alpha*y.lwork)
252: */
253: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
254: {
256:   PetscReal      d1=0,d2=0;
257:   AppCtx         *user = (AppCtx*)ptr;

260:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
261:   MatMult(user->Q,user->y,user->dwork);
262:   VecAXPY(user->dwork,-1.0,user->d);
263:   VecDot(user->dwork,user->dwork,&d1);

265:   VecWAXPY(user->uwork,-1.0,user->ur,user->u);
266:   VecPointwiseMult(user->uwork,user->uwork,user->uwork);
267:   MatMult(user->L,user->uwork,user->lwork);
268:   VecDot(user->y,user->lwork,&d2);
269:   *f = 0.5 * (d1 + user->alpha*d2);
270:   return(0);
271: }

273: /* ------------------------------------------------------------------- */
276: /*
277:     state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2
278:     design: g_d = alpha*(L'y).*(u-ur)
279: */
280: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
281: {
283:   AppCtx         *user = (AppCtx*)ptr;

286:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
287:   MatMult(user->Q,user->y,user->dwork);
288:   VecAXPY(user->dwork,-1.0,user->d);

290:   MatMult(user->QT,user->dwork,user->ywork);

292:   MatMult(user->LT,user->y,user->uwork);
293:   VecWAXPY(user->vwork,-1.0,user->ur,user->u);
294:   VecPointwiseMult(user->uwork,user->vwork,user->uwork);
295:   VecScale(user->uwork,user->alpha);

297:   VecPointwiseMult(user->vwork,user->vwork,user->vwork);
298:   MatMult(user->L,user->vwork,user->lwork);
299:   VecAXPY(user->ywork,0.5*user->alpha,user->lwork);

301:   Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
302:   return(0);
303: }

307: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
308: {
310:   PetscReal      d1,d2;
311:   AppCtx         *user = (AppCtx*)ptr;

314:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
315:   MatMult(user->Q,user->y,user->dwork);
316:   VecAXPY(user->dwork,-1.0,user->d);

318:   MatMult(user->QT,user->dwork,user->ywork);

320:   VecDot(user->dwork,user->dwork,&d1);

322:   MatMult(user->LT,user->y,user->uwork);
323:   VecWAXPY(user->vwork,-1.0,user->ur,user->u);
324:   VecPointwiseMult(user->uwork,user->vwork,user->uwork);
325:   VecScale(user->uwork,user->alpha);

327:   VecPointwiseMult(user->vwork,user->vwork,user->vwork);
328:   MatMult(user->L,user->vwork,user->lwork);
329:   VecAXPY(user->ywork,0.5*user->alpha,user->lwork);

331:   VecDot(user->y,user->lwork,&d2);

333:   *f = 0.5 * (d1 + user->alpha*d2);
334:   Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
335:   return(0);
336: }

338: /* ------------------------------------------------------------------- */
341: /* A
342: MatShell object
343: */
344: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
345: {
347:   PetscInt       i;
348:   AppCtx         *user = (AppCtx*)ptr;

351:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
352:   Scatter_yi(user->u,user->ui,user->ui_scatter,user->nt);
353:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
354:   for (i=0; i<user->nt; i++){
355:     MatCopy(user->Divxy[0],user->C[i],SAME_NONZERO_PATTERN);
356:     MatCopy(user->Divxy[1],user->Cwork[i],SAME_NONZERO_PATTERN);

358:     MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
359:     MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
360:     MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
361:     MatScale(user->C[i],user->ht);
362:     MatShift(user->C[i],1.0);
363:   }
364:   return(0);
365: }

367: /* ------------------------------------------------------------------- */
370: /* B */
371: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
372: {
374:   AppCtx         *user = (AppCtx*)ptr;

377:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
378:   return(0);
379: }

383: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
384: {
386:   PetscInt       i;
387:   AppCtx         *user;

390:   MatShellGetContext(J_shell,(void**)&user);
391:   Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
392:   user->block_index = 0;
393:   MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);

395:   for (i=1; i<user->nt; i++){
396:     user->block_index = i;
397:     MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
398:     MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
399:     VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
400:   }
401:   Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
402:   return(0);
403: }

407: PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
408: {
410:   PetscInt       i;
411:   AppCtx         *user;

414:   MatShellGetContext(J_shell,(void**)&user);
415:   Scatter_yi(X,user->yi,user->yi_scatter,user->nt);

417:   for (i=0; i<user->nt-1; i++){
418:     user->block_index = i;
419:     MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
420:     MatMult(user->M,user->yi[i+1],user->ziwork[i+1]);
421:     VecAXPY(user->yiwork[i],-1.0,user->ziwork[i+1]);
422:   }

424:   i = user->nt-1;
425:   user->block_index = i;
426:   MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
427:   Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
428:   return(0);
429: }

433: PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
434: {
436:   PetscInt       i;
437:   AppCtx         *user;

440:   MatShellGetContext(J_shell,(void**)&user);
441:   i = user->block_index;
442:   VecPointwiseMult(user->uxiwork[i],X,user->uxi[i]);
443:   VecPointwiseMult(user->uyiwork[i],X,user->uyi[i]);
444:   Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
445:   MatMult(user->Div,user->uiwork[i],Y);
446:   VecAYPX(Y,user->ht,X);
447:   return(0);
448: }

452: PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y)
453: {
455:   PetscInt       i;
456:   AppCtx         *user;

459:   MatShellGetContext(J_shell,(void**)&user);
460:   i = user->block_index;
461:   MatMult(user->Grad,X,user->uiwork[i]);
462:   Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
463:   VecPointwiseMult(user->uxiwork[i],user->uxi[i],user->uxiwork[i]);
464:   VecPointwiseMult(user->uyiwork[i],user->uyi[i],user->uyiwork[i]);
465:   VecWAXPY(Y,1.0,user->uxiwork[i],user->uyiwork[i]);
466:   VecAYPX(Y,user->ht,X);
467:   return(0);
468: }

472: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
473: {
475:   PetscInt       i;
476:   AppCtx         *user;

479:   MatShellGetContext(J_shell,(void**)&user);
480:   Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
481:   Scatter_uxi_uyi(X,user->uxiwork,user->uxi_scatter,user->uyiwork,user->uyi_scatter,user->nt);
482:   for (i=0; i<user->nt; i++){
483:     VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
484:     VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
485:     Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
486:     MatMult(user->Div,user->uiwork[i],user->ziwork[i]);
487:     VecScale(user->ziwork[i],user->ht);
488:   }
489:   Gather_yi(Y,user->ziwork,user->yi_scatter,user->nt);
490:   return(0);
491: }

495: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
496: {
498:   PetscInt       i;
499:   AppCtx         *user;

502:   MatShellGetContext(J_shell,(void**)&user);
503:   Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
504:   Scatter_yi(X,user->yiwork,user->yi_scatter,user->nt);
505:   for (i=0; i<user->nt; i++){
506:     MatMult(user->Grad,user->yiwork[i],user->uiwork[i]);
507:     Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
508:     VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
509:     VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
510:     Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
511:     VecScale(user->uiwork[i],user->ht);
512:   }
513:   Gather_yi(Y,user->uiwork,user->ui_scatter,user->nt);
514:   return(0);
515: }

519: PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
520: {
522:   PetscInt       i;
523:   AppCtx         *user;

526:   PCShellGetContext(PC_shell,(void**)&user);
527:   i = user->block_index;
528:   if (user->c_formed) {
529:     MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
530:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
531:   return(0);
532: }

536: PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
537: {
539:   PetscInt       i;
540:   AppCtx         *user;

543:   PCShellGetContext(PC_shell,(void**)&user);

545:   i = user->block_index;
546:   if (user->c_formed) {
547:     MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
548:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
549:   return(0);
550: }

554: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
555: {
557:   AppCtx         *user;
558:   PetscInt       its,i;

561:   MatShellGetContext(J_shell,(void**)&user);

563:   if (Y == user->ytrue) {
564:     /* First solve is done using true solution to set up problem */
565:     KSPSetTolerances(user->solver,1e-4,1e-20,PETSC_DEFAULT,PETSC_DEFAULT);
566:   } else {
567:     KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
568:   }
569:   Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
570:   Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
571:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

573:   user->block_index = 0;
574:   KSPSolve(user->solver,user->yi[0],user->yiwork[0]);

576:   KSPGetIterationNumber(user->solver,&its);
577:   user->ksp_its = user->ksp_its + its;
578:   for (i=1; i<user->nt; i++){
579:     MatMult(user->M,user->yiwork[i-1],user->ziwork[i-1]);
580:     VecAXPY(user->yi[i],1.0,user->ziwork[i-1]);
581:     user->block_index = i;
582:     KSPSolve(user->solver,user->yi[i],user->yiwork[i]);

584:     KSPGetIterationNumber(user->solver,&its);
585:     user->ksp_its = user->ksp_its + its;
586:   }
587:   Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
588:   return(0);
589: }

593: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
594: {
596:   AppCtx         *user;
597:   PetscInt       its,i;

600:   MatShellGetContext(J_shell,(void**)&user);

602:   Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
603:   Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
604:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

606:   i = user->nt - 1;
607:   user->block_index = i;
608:   KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);

610:   KSPGetIterationNumber(user->solver,&its);
611:   user->ksp_its = user->ksp_its + its;

613:   for (i=user->nt-2; i>=0; i--){
614:     MatMult(user->M,user->yiwork[i+1],user->ziwork[i+1]);
615:     VecAXPY(user->yi[i],1.0,user->ziwork[i+1]);
616:     user->block_index = i;
617:     KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);

619:     KSPGetIterationNumber(user->solver,&its);
620:     user->ksp_its = user->ksp_its + its;
621:   }
622:   Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
623:   return(0);
624: }

628: PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
629: {
631:   AppCtx         *user;

634:   MatShellGetContext(J_shell,(void**)&user);

636:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_shell);
637:   MatShellSetOperation(*new_shell,MATOP_MULT,(void(*)(void))StateMatMult);
638:   MatShellSetOperation(*new_shell,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
639:   MatShellSetOperation(*new_shell,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
640:   MatShellSetOperation(*new_shell,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
641:   return(0);
642: }

646: PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
647: {
649:   AppCtx         *user;

652:   MatShellGetContext(J_shell,(void**)&user);
653:    VecCopy(user->js_diag,X);
654:   return(0);
655: }

659: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
660: {
661:   /* con = Ay - q, A = [C(u1)  0     0     ...   0;
662:                          -M  C(u2)   0     ...   0;
663:                           0   -M   C(u3)   ...   0;
664:                                       ...         ;
665:                           0    ...      -M C(u_nt)]
666:      C(u) = eye + ht*Div*[diag(u1); diag(u2)]       */
668:   PetscInt       i;
669:   AppCtx         *user = (AppCtx*)ptr;

672:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
673:   Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
674:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

676:   user->block_index = 0;
677:   MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);

679:   for (i=1; i<user->nt; i++){
680:     user->block_index = i;
681:     MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
682:     MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
683:     VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
684:   }

686:   Gather_yi(C,user->yiwork,user->yi_scatter,user->nt);
687:   VecAXPY(C,-1.0,user->q);

689:   return(0);
690: }


695: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
696: {

700:   VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
701:   VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
702:   VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
703:   VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
704:   return(0);
705: }

709: PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
710: {
712:   PetscInt       i;

715:   for (i=0; i<nt; i++){
716:     VecScatterBegin(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
717:     VecScatterEnd(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
718:     VecScatterBegin(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
719:     VecScatterEnd(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
720:   }
721:   return(0);
722: }

726: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
727: {

731:   VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
732:   VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
733:   VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
734:   VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
735:   return(0);
736: }

740: PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
741: {
743:   PetscInt       i;

746:   for (i=0; i<nt; i++){
747:     VecScatterBegin(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
748:     VecScatterEnd(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
749:     VecScatterBegin(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
750:     VecScatterEnd(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
751:   }
752:   return(0);
753: }

757: PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
758: {
760:   PetscInt       i;

763:   for (i=0; i<nt; i++){
764:     VecScatterBegin(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
765:     VecScatterEnd(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
766:   }
767:   return(0);
768: }

772: PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
773: {
775:   PetscInt       i;

778:   for (i=0; i<nt; i++){
779:     VecScatterBegin(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
780:     VecScatterEnd(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
781:   }
782:   return(0);
783: }

787: PetscErrorCode HyperbolicInitialize(AppCtx *user)
788: {
790:   PetscInt       n,i,j,linear_index,istart,iend,iblock,lo,hi;
791:   Vec            XX,YY,XXwork,YYwork,yi,uxi,ui,bc;
792:   PetscReal      h,sum;
793:   PetscScalar    hinv,neg_hinv,quarter=0.25,one=1.0,half_hinv,neg_half_hinv;
794:   PetscScalar    vx,vy,zero=0.0;
795:   IS             is_from_y,is_to_yi,is_from_u,is_to_uxi,is_to_uyi;

798:   user->jformed = PETSC_FALSE;
799:   user->c_formed = PETSC_FALSE;

801:   user->ksp_its = 0;
802:   user->ksp_its_initial = 0;

804:   n = user->mx * user->mx;

806:   h = 1.0/user->mx;
807:   hinv = user->mx;
808:   neg_hinv = -hinv;
809:   half_hinv = hinv / 2.0;
810:   neg_half_hinv = neg_hinv / 2.0;

812:   /* Generate Grad matrix */
813:   MatCreate(PETSC_COMM_WORLD,&user->Grad);
814:   MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,2*n,n);
815:   MatSetFromOptions(user->Grad);
816:   MatMPIAIJSetPreallocation(user->Grad,3,NULL,3,NULL);
817:   MatSeqAIJSetPreallocation(user->Grad,3,NULL);
818:   MatGetOwnershipRange(user->Grad,&istart,&iend);

820:   for (i=istart; i<iend; i++){
821:     if (i<n){
822:       iblock = i / user->mx;
823:       j = iblock*user->mx + ((i+user->mx-1) % user->mx);
824:       MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
825:       j = iblock*user->mx + ((i+1) % user->mx);
826:       MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
827:     }
828:     if (i>=n){
829:       j = (i - user->mx) % n;
830:       MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
831:       j = (j + 2*user->mx) % n;
832:       MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
833:     }
834:   }

836:   MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
837:   MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);

839:   MatCreate(PETSC_COMM_WORLD,&user->Gradxy[0]);
840:   MatSetSizes(user->Gradxy[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
841:   MatSetFromOptions(user->Gradxy[0]);
842:   MatMPIAIJSetPreallocation(user->Gradxy[0],3,NULL,3,NULL);
843:   MatSeqAIJSetPreallocation(user->Gradxy[0],3,NULL);
844:   MatGetOwnershipRange(user->Gradxy[0],&istart,&iend);

846:   for (i=istart; i<iend; i++){
847:     iblock = i / user->mx;
848:     j = iblock*user->mx + ((i+user->mx-1) % user->mx);
849:     MatSetValues(user->Gradxy[0],1,&i,1,&j,&half_hinv,INSERT_VALUES);
850:     j = iblock*user->mx + ((i+1) % user->mx);
851:     MatSetValues(user->Gradxy[0],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
852:     MatSetValues(user->Gradxy[0],1,&i,1,&i,&zero,INSERT_VALUES);
853:   }
854:   MatAssemblyBegin(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
855:   MatAssemblyEnd(user->Gradxy[0],MAT_FINAL_ASSEMBLY);

857:   MatCreate(PETSC_COMM_WORLD,&user->Gradxy[1]);
858:   MatSetSizes(user->Gradxy[1],PETSC_DECIDE,PETSC_DECIDE,n,n);
859:   MatSetFromOptions(user->Gradxy[1]);
860:   MatMPIAIJSetPreallocation(user->Gradxy[1],3,NULL,3,NULL);
861:   MatSeqAIJSetPreallocation(user->Gradxy[1],3,NULL);
862:   MatGetOwnershipRange(user->Gradxy[1],&istart,&iend);

864:   for (i=istart; i<iend; i++){
865:     j = (i + n - user->mx) % n;
866:     MatSetValues(user->Gradxy[1],1,&i,1,&j,&half_hinv,INSERT_VALUES);
867:     j = (j + 2*user->mx) % n;
868:     MatSetValues(user->Gradxy[1],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
869:     MatSetValues(user->Gradxy[1],1,&i,1,&i,&zero,INSERT_VALUES);
870:   }
871:   MatAssemblyBegin(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
872:   MatAssemblyEnd(user->Gradxy[1],MAT_FINAL_ASSEMBLY);

874:   /* Generate Div matrix */
875:   MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div);
876:   MatTranspose(user->Gradxy[0],MAT_INITIAL_MATRIX,&user->Divxy[0]);
877:   MatTranspose(user->Gradxy[1],MAT_INITIAL_MATRIX,&user->Divxy[1]);

879:   /* Off-diagonal averaging matrix */
880:   MatCreate(PETSC_COMM_WORLD,&user->M);
881:   MatSetSizes(user->M,PETSC_DECIDE,PETSC_DECIDE,n,n);
882:   MatSetFromOptions(user->M);
883:   MatMPIAIJSetPreallocation(user->M,4,NULL,4,NULL);
884:   MatSeqAIJSetPreallocation(user->M,4,NULL);
885:   MatGetOwnershipRange(user->M,&istart,&iend);

887:   for (i=istart; i<iend; i++){
888:     /* kron(Id,Av) */
889:     iblock = i / user->mx;
890:     j = iblock*user->mx + ((i+user->mx-1) % user->mx);
891:     MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
892:     j = iblock*user->mx + ((i+1) % user->mx);
893:     MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);

895:     /* kron(Av,Id) */
896:     j = (i + user->mx) % n;
897:     MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
898:     j = (i + n - user->mx) % n;
899:     MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
900:   }
901:   MatAssemblyBegin(user->M,MAT_FINAL_ASSEMBLY);
902:   MatAssemblyEnd(user->M,MAT_FINAL_ASSEMBLY);

904:   /* Generate 2D grid */
905:   VecCreate(PETSC_COMM_WORLD,&XX);
906:   VecCreate(PETSC_COMM_WORLD,&user->q);
907:   VecSetSizes(XX,PETSC_DECIDE,n);
908:   VecSetSizes(user->q,PETSC_DECIDE,n*user->nt);
909:   VecSetFromOptions(XX);
910:   VecSetFromOptions(user->q);

912:   VecDuplicate(XX,&YY);
913:   VecDuplicate(XX,&XXwork);
914:   VecDuplicate(XX,&YYwork);
915:   VecDuplicate(XX,&user->d);
916:   VecDuplicate(XX,&user->dwork);

918:   VecGetOwnershipRange(XX,&istart,&iend);
919:   for (linear_index=istart; linear_index<iend; linear_index++){
920:     i = linear_index % user->mx;
921:     j = (linear_index-i)/user->mx;
922:     vx = h*(i+0.5);
923:     vy = h*(j+0.5);
924:     VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
925:     VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
926:   }

928:   VecAssemblyBegin(XX);
929:   VecAssemblyEnd(XX);
930:   VecAssemblyBegin(YY);
931:   VecAssemblyEnd(YY);

933:   /* Compute final density function yT
934:      yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
935:      yT = yT / (h^2*sum(yT)) */
936:   VecCopy(XX,XXwork);
937:   VecCopy(YY,YYwork);

939:   VecShift(XXwork,-0.25);
940:   VecShift(YYwork,-0.25);

942:   VecPointwiseMult(XXwork,XXwork,XXwork);
943:   VecPointwiseMult(YYwork,YYwork,YYwork);

945:   VecCopy(XXwork,user->dwork);
946:   VecAXPY(user->dwork,1.0,YYwork);
947:   VecScale(user->dwork,-30.0);
948:   VecExp(user->dwork);
949:   VecCopy(user->dwork,user->d);

951:   VecCopy(XX,XXwork);
952:   VecCopy(YY,YYwork);

954:   VecShift(XXwork,-0.75);
955:   VecShift(YYwork,-0.75);

957:   VecPointwiseMult(XXwork,XXwork,XXwork);
958:   VecPointwiseMult(YYwork,YYwork,YYwork);

960:   VecCopy(XXwork,user->dwork);
961:   VecAXPY(user->dwork,1.0,YYwork);
962:   VecScale(user->dwork,-30.0);
963:   VecExp(user->dwork);

965:   VecAXPY(user->d,1.0,user->dwork);
966:   VecShift(user->d,1.0);
967:   VecSum(user->d,&sum);
968:   VecScale(user->d,1.0/(h*h*sum));

970:   /* Initial conditions of forward problem */
971:   VecDuplicate(XX,&bc);
972:   VecCopy(XX,XXwork);
973:   VecCopy(YY,YYwork);

975:   VecShift(XXwork,-0.5);
976:   VecShift(YYwork,-0.5);

978:   VecPointwiseMult(XXwork,XXwork,XXwork);
979:   VecPointwiseMult(YYwork,YYwork,YYwork);

981:   VecWAXPY(bc,1.0,XXwork,YYwork);
982:   VecScale(bc,-50.0);
983:   VecExp(bc);
984:   VecShift(bc,1.0);
985:   VecSum(bc,&sum);
986:   VecScale(bc,1.0/(h*h*sum));

988:   /* Create scatter from y to y_1,y_2,...,y_nt */
989:   /*  TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
990:   PetscMalloc1(user->nt*user->mx*user->mx,&user->yi_scatter);
991:   VecCreate(PETSC_COMM_WORLD,&yi);
992:   VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx);
993:   VecSetFromOptions(yi);
994:   VecDuplicateVecs(yi,user->nt,&user->yi);
995:   VecDuplicateVecs(yi,user->nt,&user->yiwork);
996:   VecDuplicateVecs(yi,user->nt,&user->ziwork);
997:   for (i=0; i<user->nt; i++){
998:     VecGetOwnershipRange(user->yi[i],&lo,&hi);
999:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_yi);
1000:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+i*user->mx*user->mx,1,&is_from_y);
1001:     VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i]);
1002:     ISDestroy(&is_to_yi);
1003:     ISDestroy(&is_from_y);
1004:   }

1006:   /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
1007:   /*  TODO: reorder for better parallelism */
1008:   PetscMalloc1(user->nt*user->mx*user->mx,&user->uxi_scatter);
1009:   PetscMalloc1(user->nt*user->mx*user->mx,&user->uyi_scatter);
1010:   PetscMalloc1(user->nt*user->mx*user->mx,&user->ux_scatter);
1011:   PetscMalloc1(user->nt*user->mx*user->mx,&user->uy_scatter);
1012:   PetscMalloc1(2*user->nt*user->mx*user->mx,&user->ui_scatter);
1013:   VecCreate(PETSC_COMM_WORLD,&uxi);
1014:   VecCreate(PETSC_COMM_WORLD,&ui);
1015:   VecSetSizes(uxi,PETSC_DECIDE,user->mx*user->mx);
1016:   VecSetSizes(ui,PETSC_DECIDE,2*user->mx*user->mx);
1017:   VecSetFromOptions(uxi);
1018:   VecSetFromOptions(ui);
1019:   VecDuplicateVecs(uxi,user->nt,&user->uxi);
1020:   VecDuplicateVecs(uxi,user->nt,&user->uyi);
1021:   VecDuplicateVecs(uxi,user->nt,&user->uxiwork);
1022:   VecDuplicateVecs(uxi,user->nt,&user->uyiwork);
1023:   VecDuplicateVecs(ui,user->nt,&user->ui);
1024:   VecDuplicateVecs(ui,user->nt,&user->uiwork);
1025:   for (i=0; i<user->nt; i++){
1026:     VecGetOwnershipRange(user->uxi[i],&lo,&hi);
1027:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1028:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
1029:     VecScatterCreate(user->u,is_from_u,user->uxi[i],is_to_uxi,&user->uxi_scatter[i]);

1031:     ISDestroy(&is_to_uxi);
1032:     ISDestroy(&is_from_u);

1034:     VecGetOwnershipRange(user->uyi[i],&lo,&hi);
1035:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
1036:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+(2*i+1)*user->mx*user->mx,1,&is_from_u);
1037:     VecScatterCreate(user->u,is_from_u,user->uyi[i],is_to_uyi,&user->uyi_scatter[i]);

1039:     ISDestroy(&is_to_uyi);
1040:     ISDestroy(&is_from_u);

1042:     VecGetOwnershipRange(user->uxi[i],&lo,&hi);
1043:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1044:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_from_u);
1045:     VecScatterCreate(user->ui[i],is_from_u,user->uxi[i],is_to_uxi,&user->ux_scatter[i]);

1047:     ISDestroy(&is_to_uxi);
1048:     ISDestroy(&is_from_u);

1050:     VecGetOwnershipRange(user->uyi[i],&lo,&hi);
1051:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
1052:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+user->mx*user->mx,1,&is_from_u);
1053:     VecScatterCreate(user->ui[i],is_from_u,user->uyi[i],is_to_uyi,&user->uy_scatter[i]);

1055:     ISDestroy(&is_to_uyi);
1056:     ISDestroy(&is_from_u);

1058:     VecGetOwnershipRange(user->ui[i],&lo,&hi);
1059:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
1060:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
1061:     VecScatterCreate(user->u,is_from_u,user->ui[i],is_to_uxi,&user->ui_scatter[i]);

1063:     ISDestroy(&is_to_uxi);
1064:     ISDestroy(&is_from_u);
1065:   }

1067:   /* RHS of forward problem */
1068:   MatMult(user->M,bc,user->yiwork[0]);
1069:   for (i=1; i<user->nt; i++){
1070:     VecSet(user->yiwork[i],0.0);
1071:   }
1072:   Gather_yi(user->q,user->yiwork,user->yi_scatter,user->nt);

1074:   /* Compute true velocity field utrue */
1075:   VecDuplicate(user->u,&user->utrue);
1076:   for (i=0; i<user->nt; i++){
1077:     VecCopy(YY,user->uxi[i]);
1078:     VecScale(user->uxi[i],150.0*i*user->ht);
1079:     VecCopy(XX,user->uyi[i]);
1080:     VecShift(user->uyi[i],-10.0);
1081:     VecScale(user->uyi[i],15.0*i*user->ht);
1082:   }
1083:   Gather_uxi_uyi(user->utrue,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

1085:   /* Initial guess and reference model */
1086:   VecDuplicate(user->utrue,&user->ur);
1087:   for (i=0; i<user->nt; i++){
1088:     VecCopy(XX,user->uxi[i]);
1089:     VecShift(user->uxi[i],i*user->ht);
1090:     VecCopy(YY,user->uyi[i]);
1091:     VecShift(user->uyi[i],-i*user->ht);
1092:   }
1093:   Gather_uxi_uyi(user->ur,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

1095:   /* Generate regularization matrix L */
1096:   MatCreate(PETSC_COMM_WORLD,&user->LT);
1097:   MatSetSizes(user->LT,PETSC_DECIDE,PETSC_DECIDE,2*n*user->nt,n*user->nt);
1098:   MatSetFromOptions(user->LT);
1099:   MatMPIAIJSetPreallocation(user->LT,1,NULL,1,NULL);
1100:   MatSeqAIJSetPreallocation(user->LT,1,NULL);
1101:   MatGetOwnershipRange(user->LT,&istart,&iend);

1103:   for (i=istart; i<iend; i++){
1104:     iblock = (i+n) / (2*n);
1105:     j = i - iblock*n;
1106:     MatSetValues(user->LT,1,&i,1,&j,&user->gamma,INSERT_VALUES);
1107:   }

1109:   MatAssemblyBegin(user->LT,MAT_FINAL_ASSEMBLY);
1110:   MatAssemblyEnd(user->LT,MAT_FINAL_ASSEMBLY);

1112:   MatTranspose(user->LT,MAT_INITIAL_MATRIX,&user->L);

1114:   /* Build work vectors and matrices */
1115:   VecCreate(PETSC_COMM_WORLD,&user->lwork);
1116:   VecSetType(user->lwork,VECMPI);
1117:   VecSetSizes(user->lwork,PETSC_DECIDE,user->m);
1118:   VecSetFromOptions(user->lwork);

1120:   MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);

1122:   VecDuplicate(user->y,&user->ywork);
1123:   VecDuplicate(user->u,&user->uwork);
1124:   VecDuplicate(user->u,&user->vwork);
1125:   VecDuplicate(user->u,&user->js_diag);
1126:   VecDuplicate(user->c,&user->cwork);

1128:   /* Create matrix-free shell user->Js for computing A*x */
1129:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->Js);
1130:   MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1131:   MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
1132:   MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
1133:   MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);

1135:   /* Diagonal blocks of user->Js */
1136:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlock);
1137:   MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult);
1138:   MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockMultTranspose);

1140:   /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
1141:      D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
1142:      This is an SOR preconditioner for user->JsBlock. */
1143:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlockPrec);
1144:   MatShellSetOperation(user->JsBlockPrec,MATOP_MULT,(void(*)(void))StateMatBlockPrecMult);
1145:   MatShellSetOperation(user->JsBlockPrec,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockPrecMultTranspose);

1147:   /* Create a matrix-free shell user->Jd for computing B*x */
1148:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->n-user->m,user,&user->Jd);
1149:   MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);
1150:   MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);

1152:   /* User-defined routines for computing user->Js\x and user->Js^T\x*/
1153:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsInv);
1154:   MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateMatInvMult);
1155:   MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatInvTransposeMult);

1157:   /* Build matrices for SOR preconditioner */
1158:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1159:   MatShift(user->Divxy[0],0.0); /*  Force C[i] and Divxy[0] to share same nonzero pattern */
1160:   MatAXPY(user->Divxy[0],0.0,user->Divxy[1],DIFFERENT_NONZERO_PATTERN);
1161:   PetscMalloc1(5*n,&user->C);
1162:   PetscMalloc1(2*n,&user->Cwork);
1163:   for (i=0; i<user->nt; i++){
1164:     MatDuplicate(user->Divxy[0],MAT_COPY_VALUES,&user->C[i]);
1165:     MatDuplicate(user->Divxy[1],MAT_COPY_VALUES,&user->Cwork[i]);

1167:     MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
1168:     MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
1169:     MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
1170:     MatScale(user->C[i],user->ht);
1171:     MatShift(user->C[i],1.0);
1172:   }

1174:   /* Solver options and tolerances */
1175:   KSPCreate(PETSC_COMM_WORLD,&user->solver);
1176:   KSPSetType(user->solver,KSPGMRES);
1177:   KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec);
1178:   KSPSetInitialGuessNonzero(user->solver,PETSC_FALSE); /*  TODO: why is true slower? */
1179:   KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500);
1180:   /* KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500); */
1181:   KSPGetPC(user->solver,&user->prec);
1182:   PCSetType(user->prec,PCSHELL);

1184:   PCShellSetApply(user->prec,StateMatBlockPrecMult);
1185:   PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMultTranspose);
1186:   PCShellSetContext(user->prec,user);

1188:   /* Compute true state function yt given ut */
1189:   VecCreate(PETSC_COMM_WORLD,&user->ytrue);
1190:   VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt);
1191:   VecSetFromOptions(user->ytrue);
1192:   user->c_formed = PETSC_TRUE;
1193:   VecCopy(user->utrue,user->u); /*  Set u=utrue temporarily for StateMatInv */
1194:   VecSet(user->ytrue,0.0); /*  Initial guess */
1195:   StateMatInvMult(user->Js,user->q,user->ytrue);
1196:   VecCopy(user->ur,user->u); /*  Reset u=ur */

1198:   /* Initial guess y0 for state given u0 */
1199:   StateMatInvMult(user->Js,user->q,user->y);

1201:   /* Data discretization */
1202:   MatCreate(PETSC_COMM_WORLD,&user->Q);
1203:   MatSetSizes(user->Q,PETSC_DECIDE,PETSC_DECIDE,user->mx*user->mx,user->m);
1204:   MatSetFromOptions(user->Q);
1205:   MatMPIAIJSetPreallocation(user->Q,0,NULL,1,NULL);
1206:   MatSeqAIJSetPreallocation(user->Q,1,NULL);

1208:   MatGetOwnershipRange(user->Q,&istart,&iend);

1210:   for (i=istart; i<iend; i++){
1211:     j = i + user->m - user->mx*user->mx;
1212:     MatSetValues(user->Q,1,&i,1,&j,&one,INSERT_VALUES);
1213:   }

1215:   MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1216:   MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);

1218:   MatTranspose(user->Q,MAT_INITIAL_MATRIX,&user->QT);

1220:   VecDestroy(&XX);
1221:   VecDestroy(&YY);
1222:   VecDestroy(&XXwork);
1223:   VecDestroy(&YYwork);
1224:   VecDestroy(&yi);
1225:   VecDestroy(&uxi);
1226:   VecDestroy(&ui);
1227:   VecDestroy(&bc);

1229:   /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1230:   KSPSetFromOptions(user->solver);
1231:   return(0);
1232: }

1236: PetscErrorCode HyperbolicDestroy(AppCtx *user)
1237: {
1239:   PetscInt       i;

1242:   MatDestroy(&user->Q);
1243:   MatDestroy(&user->QT);
1244:   MatDestroy(&user->Div);
1245:   MatDestroy(&user->Divwork);
1246:   MatDestroy(&user->Grad);
1247:   MatDestroy(&user->L);
1248:   MatDestroy(&user->LT);
1249:   KSPDestroy(&user->solver);
1250:   MatDestroy(&user->Js);
1251:   MatDestroy(&user->Jd);
1252:   MatDestroy(&user->JsBlockPrec);
1253:   MatDestroy(&user->JsInv);
1254:   MatDestroy(&user->JsBlock);
1255:   MatDestroy(&user->Divxy[0]);
1256:   MatDestroy(&user->Divxy[1]);
1257:   MatDestroy(&user->Gradxy[0]);
1258:   MatDestroy(&user->Gradxy[1]);
1259:   MatDestroy(&user->M);
1260:   for (i=0; i<user->nt; i++){
1261:     MatDestroy(&user->C[i]);
1262:     MatDestroy(&user->Cwork[i]);
1263:   }
1264:   PetscFree(user->C);
1265:   PetscFree(user->Cwork);
1266:   VecDestroy(&user->u);
1267:   VecDestroy(&user->uwork);
1268:   VecDestroy(&user->vwork);
1269:   VecDestroy(&user->utrue);
1270:   VecDestroy(&user->y);
1271:   VecDestroy(&user->ywork);
1272:   VecDestroy(&user->ytrue);
1273:   VecDestroyVecs(user->nt,&user->yi);
1274:   VecDestroyVecs(user->nt,&user->yiwork);
1275:   VecDestroyVecs(user->nt,&user->ziwork);
1276:   VecDestroyVecs(user->nt,&user->uxi);
1277:   VecDestroyVecs(user->nt,&user->uyi);
1278:   VecDestroyVecs(user->nt,&user->uxiwork);
1279:   VecDestroyVecs(user->nt,&user->uyiwork);
1280:   VecDestroyVecs(user->nt,&user->ui);
1281:   VecDestroyVecs(user->nt,&user->uiwork);
1282:   VecDestroy(&user->c);
1283:   VecDestroy(&user->cwork);
1284:   VecDestroy(&user->ur);
1285:   VecDestroy(&user->q);
1286:   VecDestroy(&user->d);
1287:   VecDestroy(&user->dwork);
1288:   VecDestroy(&user->lwork);
1289:   VecDestroy(&user->js_diag);
1290:   ISDestroy(&user->s_is);
1291:   ISDestroy(&user->d_is);
1292:   VecScatterDestroy(&user->state_scatter);
1293:   VecScatterDestroy(&user->design_scatter);
1294:   for (i=0; i<user->nt; i++){
1295:     VecScatterDestroy(&user->uxi_scatter[i]);
1296:     VecScatterDestroy(&user->uyi_scatter[i]);
1297:     VecScatterDestroy(&user->ux_scatter[i]);
1298:     VecScatterDestroy(&user->uy_scatter[i]);
1299:     VecScatterDestroy(&user->ui_scatter[i]);
1300:     VecScatterDestroy(&user->yi_scatter[i]);
1301:   }
1302:   PetscFree(user->uxi_scatter);
1303:   PetscFree(user->uyi_scatter);
1304:   PetscFree(user->ux_scatter);
1305:   PetscFree(user->uy_scatter);
1306:   PetscFree(user->ui_scatter);
1307:   PetscFree(user->yi_scatter);
1308:   return(0);
1309: }

1313: PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1314: {
1316:   Vec            X;
1317:   PetscReal      unorm,ynorm;
1318:   AppCtx         *user = (AppCtx*)ptr;

1321:   TaoGetSolutionVector(tao,&X);
1322:   Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1323:   VecAXPY(user->ywork,-1.0,user->ytrue);
1324:   VecAXPY(user->uwork,-1.0,user->utrue);
1325:   VecNorm(user->uwork,NORM_2,&unorm);
1326:   VecNorm(user->ywork,NORM_2,&ynorm);
1327:   PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1328:   return(0);
1329: }