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: TaoSetSolution();
  8:    Routines: TaoSetObjective();
  9:    Routines: TaoSetGradient();
 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);
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:   TaoSetSolution(tao,x);
195:   TaoSetObjective(tao, FormFunction, &user);
196:   TaoSetGradient(tao, NULL, 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:     TaoSetSolution(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 0;
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: {
238:   PetscReal      d1=0,d2=0;
239:   AppCtx         *user = (AppCtx*)ptr;

241:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
242:   MatMult(user->Q,user->y,user->dwork);
243:   VecAXPY(user->dwork,-1.0,user->d);
244:   VecDot(user->dwork,user->dwork,&d1);

246:   VecWAXPY(user->uwork,-1.0,user->ur,user->u);
247:   VecPointwiseMult(user->uwork,user->uwork,user->uwork);
248:   MatMult(user->L,user->uwork,user->lwork);
249:   VecDot(user->y,user->lwork,&d2);
250:   *f = 0.5 * (d1 + user->alpha*d2);
251:   return 0;
252: }

254: /* ------------------------------------------------------------------- */
255: /*
256:     state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2
257:     design: g_d = alpha*(L'y).*(u-ur)
258: */
259: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
260: {
261:   AppCtx         *user = (AppCtx*)ptr;

263:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
264:   MatMult(user->Q,user->y,user->dwork);
265:   VecAXPY(user->dwork,-1.0,user->d);

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

269:   MatMult(user->LT,user->y,user->uwork);
270:   VecWAXPY(user->vwork,-1.0,user->ur,user->u);
271:   VecPointwiseMult(user->uwork,user->vwork,user->uwork);
272:   VecScale(user->uwork,user->alpha);

274:   VecPointwiseMult(user->vwork,user->vwork,user->vwork);
275:   MatMult(user->L,user->vwork,user->lwork);
276:   VecAXPY(user->ywork,0.5*user->alpha,user->lwork);

278:   Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
279:   return 0;
280: }

282: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
283: {
284:   PetscReal      d1,d2;
285:   AppCtx         *user = (AppCtx*)ptr;

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

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

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

295:   MatMult(user->LT,user->y,user->uwork);
296:   VecWAXPY(user->vwork,-1.0,user->ur,user->u);
297:   VecPointwiseMult(user->uwork,user->vwork,user->uwork);
298:   VecScale(user->uwork,user->alpha);

300:   VecPointwiseMult(user->vwork,user->vwork,user->vwork);
301:   MatMult(user->L,user->vwork,user->lwork);
302:   VecAXPY(user->ywork,0.5*user->alpha,user->lwork);

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

306:   *f = 0.5 * (d1 + user->alpha*d2);
307:   Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
308:   return 0;
309: }

311: /* ------------------------------------------------------------------- */
312: /* A
313: MatShell object
314: */
315: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
316: {
317:   PetscInt       i;
318:   AppCtx         *user = (AppCtx*)ptr;

320:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
321:   Scatter_yi(user->u,user->ui,user->ui_scatter,user->nt);
322:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
323:   for (i=0; i<user->nt; i++) {
324:     MatCopy(user->Divxy[0],user->C[i],SUBSET_NONZERO_PATTERN);
325:     MatCopy(user->Divxy[1],user->Cwork[i],SAME_NONZERO_PATTERN);

327:     MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
328:     MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
329:     MatAXPY(user->C[i],1.0,user->Cwork[i],SUBSET_NONZERO_PATTERN);
330:     MatScale(user->C[i],user->ht);
331:     MatShift(user->C[i],1.0);
332:   }
333:   return 0;
334: }

336: /* ------------------------------------------------------------------- */
337: /* B */
338: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
339: {
340:   AppCtx         *user = (AppCtx*)ptr;

342:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
343:   return 0;
344: }

346: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
347: {
348:   PetscInt       i;
349:   AppCtx         *user;

351:   MatShellGetContext(J_shell,&user);
352:   Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
353:   user->block_index = 0;
354:   MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);

356:   for (i=1; i<user->nt; i++) {
357:     user->block_index = i;
358:     MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
359:     MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
360:     VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
361:   }
362:   Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
363:   return 0;
364: }

366: PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
367: {
368:   PetscInt       i;
369:   AppCtx         *user;

371:   MatShellGetContext(J_shell,&user);
372:   Scatter_yi(X,user->yi,user->yi_scatter,user->nt);

374:   for (i=0; i<user->nt-1; i++) {
375:     user->block_index = i;
376:     MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
377:     MatMult(user->M,user->yi[i+1],user->ziwork[i+1]);
378:     VecAXPY(user->yiwork[i],-1.0,user->ziwork[i+1]);
379:   }

381:   i = user->nt-1;
382:   user->block_index = i;
383:   MatMultTranspose(user->JsBlock,user->yi[i],user->yiwork[i]);
384:   Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
385:   return 0;
386: }

388: PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
389: {
390:   PetscInt       i;
391:   AppCtx         *user;

393:   MatShellGetContext(J_shell,&user);
394:   i = user->block_index;
395:   VecPointwiseMult(user->uxiwork[i],X,user->uxi[i]);
396:   VecPointwiseMult(user->uyiwork[i],X,user->uyi[i]);
397:   Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
398:   MatMult(user->Div,user->uiwork[i],Y);
399:   VecAYPX(Y,user->ht,X);
400:   return 0;
401: }

403: PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y)
404: {
405:   PetscInt       i;
406:   AppCtx         *user;

408:   MatShellGetContext(J_shell,&user);
409:   i = user->block_index;
410:   MatMult(user->Grad,X,user->uiwork[i]);
411:   Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
412:   VecPointwiseMult(user->uxiwork[i],user->uxi[i],user->uxiwork[i]);
413:   VecPointwiseMult(user->uyiwork[i],user->uyi[i],user->uyiwork[i]);
414:   VecWAXPY(Y,1.0,user->uxiwork[i],user->uyiwork[i]);
415:   VecAYPX(Y,user->ht,X);
416:   return 0;
417: }

419: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
420: {
421:   PetscInt       i;
422:   AppCtx         *user;

424:   MatShellGetContext(J_shell,&user);
425:   Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
426:   Scatter_uxi_uyi(X,user->uxiwork,user->uxi_scatter,user->uyiwork,user->uyi_scatter,user->nt);
427:   for (i=0; i<user->nt; i++) {
428:     VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
429:     VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
430:     Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
431:     MatMult(user->Div,user->uiwork[i],user->ziwork[i]);
432:     VecScale(user->ziwork[i],user->ht);
433:   }
434:   Gather_yi(Y,user->ziwork,user->yi_scatter,user->nt);
435:   return 0;
436: }

438: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
439: {
440:   PetscInt       i;
441:   AppCtx         *user;

443:   MatShellGetContext(J_shell,&user);
444:   Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
445:   Scatter_yi(X,user->yiwork,user->yi_scatter,user->nt);
446:   for (i=0; i<user->nt; i++) {
447:     MatMult(user->Grad,user->yiwork[i],user->uiwork[i]);
448:     Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
449:     VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);
450:     VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);
451:     Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);
452:     VecScale(user->uiwork[i],user->ht);
453:   }
454:   Gather_yi(Y,user->uiwork,user->ui_scatter,user->nt);
455:   return 0;
456: }

458: PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
459: {
460:   PetscInt       i;
461:   AppCtx         *user;

463:   PCShellGetContext(PC_shell,&user);
464:   i = user->block_index;
465:   if (user->c_formed) {
466:     MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
467:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
468:   return 0;
469: }

471: PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
472: {
473:   PetscInt       i;
474:   AppCtx         *user;

476:   PCShellGetContext(PC_shell,&user);

478:   i = user->block_index;
479:   if (user->c_formed) {
480:     MatSOR(user->C[i],X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y);
481:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not formed");
482:   return 0;
483: }

485: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
486: {
487:   AppCtx         *user;
488:   PetscInt       its,i;

490:   MatShellGetContext(J_shell,&user);

492:   if (Y == user->ytrue) {
493:     /* First solve is done using true solution to set up problem */
494:     KSPSetTolerances(user->solver,1e-4,1e-20,PETSC_DEFAULT,PETSC_DEFAULT);
495:   } else {
496:     KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
497:   }
498:   Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
499:   Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
500:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

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

505:   KSPGetIterationNumber(user->solver,&its);
506:   user->ksp_its = user->ksp_its + its;
507:   for (i=1; i<user->nt; i++) {
508:     MatMult(user->M,user->yiwork[i-1],user->ziwork[i-1]);
509:     VecAXPY(user->yi[i],1.0,user->ziwork[i-1]);
510:     user->block_index = i;
511:     KSPSolve(user->solver,user->yi[i],user->yiwork[i]);

513:     KSPGetIterationNumber(user->solver,&its);
514:     user->ksp_its = user->ksp_its + its;
515:   }
516:   Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
517:   return 0;
518: }

520: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
521: {
522:   AppCtx         *user;
523:   PetscInt       its,i;

525:   MatShellGetContext(J_shell,&user);

527:   Scatter_yi(X,user->yi,user->yi_scatter,user->nt);
528:   Scatter_yi(Y,user->yiwork,user->yi_scatter,user->nt);
529:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

531:   i = user->nt - 1;
532:   user->block_index = i;
533:   KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);

535:   KSPGetIterationNumber(user->solver,&its);
536:   user->ksp_its = user->ksp_its + its;

538:   for (i=user->nt-2; i>=0; i--) {
539:     MatMult(user->M,user->yiwork[i+1],user->ziwork[i+1]);
540:     VecAXPY(user->yi[i],1.0,user->ziwork[i+1]);
541:     user->block_index = i;
542:     KSPSolveTranspose(user->solver,user->yi[i],user->yiwork[i]);

544:     KSPGetIterationNumber(user->solver,&its);
545:     user->ksp_its = user->ksp_its + its;
546:   }
547:   Gather_yi(Y,user->yiwork,user->yi_scatter,user->nt);
548:   return 0;
549: }

551: PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
552: {
553:   AppCtx         *user;

555:   MatShellGetContext(J_shell,&user);

557:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_shell);
558:   MatShellSetOperation(*new_shell,MATOP_MULT,(void(*)(void))StateMatMult);
559:   MatShellSetOperation(*new_shell,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
560:   MatShellSetOperation(*new_shell,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
561:   MatShellSetOperation(*new_shell,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);
562:   return 0;
563: }

565: PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
566: {
567:   AppCtx         *user;

569:   MatShellGetContext(J_shell,&user);
570:   VecCopy(user->js_diag,X);
571:   return 0;
572: }

574: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
575: {
576:   /* con = Ay - q, A = [C(u1)  0     0     ...   0;
577:                          -M  C(u2)   0     ...   0;
578:                           0   -M   C(u3)   ...   0;
579:                                       ...         ;
580:                           0    ...      -M C(u_nt)]
581:      C(u) = eye + ht*Div*[diag(u1); diag(u2)]       */
582:   PetscInt       i;
583:   AppCtx         *user = (AppCtx*)ptr;

585:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
586:   Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);
587:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

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

592:   for (i=1; i<user->nt; i++) {
593:     user->block_index = i;
594:     MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);
595:     MatMult(user->M,user->yi[i-1],user->ziwork[i-1]);
596:     VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);
597:   }

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

602:   return 0;
603: }

605: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
606: {
607:   VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
608:   VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);
609:   VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
610:   VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);
611:   return 0;
612: }

614: PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
615: {
616:   PetscInt       i;

618:   for (i=0; i<nt; i++) {
619:     VecScatterBegin(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
620:     VecScatterEnd(scatx[i],u,uxi[i],INSERT_VALUES,SCATTER_FORWARD);
621:     VecScatterBegin(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
622:     VecScatterEnd(scaty[i],u,uyi[i],INSERT_VALUES,SCATTER_FORWARD);
623:   }
624:   return 0;
625: }

627: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
628: {
629:   VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
630:   VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);
631:   VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
632:   VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);
633:   return 0;
634: }

636: PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
637: {
638:   PetscInt       i;

640:   for (i=0; i<nt; i++) {
641:     VecScatterBegin(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
642:     VecScatterEnd(scatx[i],uxi[i],u,INSERT_VALUES,SCATTER_REVERSE);
643:     VecScatterBegin(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
644:     VecScatterEnd(scaty[i],uyi[i],u,INSERT_VALUES,SCATTER_REVERSE);
645:   }
646:   return 0;
647: }

649: PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
650: {
651:   PetscInt       i;

653:   for (i=0; i<nt; i++) {
654:     VecScatterBegin(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
655:     VecScatterEnd(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD);
656:   }
657:   return 0;
658: }

660: PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
661: {
662:   PetscInt       i;

664:   for (i=0; i<nt; i++) {
665:     VecScatterBegin(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
666:     VecScatterEnd(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE);
667:   }
668:   return 0;
669: }

671: PetscErrorCode HyperbolicInitialize(AppCtx *user)
672: {
673:   PetscInt       n,i,j,linear_index,istart,iend,iblock,lo,hi;
674:   Vec            XX,YY,XXwork,YYwork,yi,uxi,ui,bc;
675:   PetscReal      h,sum;
676:   PetscScalar    hinv,neg_hinv,quarter=0.25,one=1.0,half_hinv,neg_half_hinv;
677:   PetscScalar    vx,vy,zero=0.0;
678:   IS             is_from_y,is_to_yi,is_from_u,is_to_uxi,is_to_uyi;

680:   user->jformed = PETSC_FALSE;
681:   user->c_formed = PETSC_FALSE;

683:   user->ksp_its = 0;
684:   user->ksp_its_initial = 0;

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

688:   h = 1.0/user->mx;
689:   hinv = user->mx;
690:   neg_hinv = -hinv;
691:   half_hinv = hinv / 2.0;
692:   neg_half_hinv = neg_hinv / 2.0;

694:   /* Generate Grad matrix */
695:   MatCreate(PETSC_COMM_WORLD,&user->Grad);
696:   MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,2*n,n);
697:   MatSetFromOptions(user->Grad);
698:   MatMPIAIJSetPreallocation(user->Grad,3,NULL,3,NULL);
699:   MatSeqAIJSetPreallocation(user->Grad,3,NULL);
700:   MatGetOwnershipRange(user->Grad,&istart,&iend);

702:   for (i=istart; i<iend; i++) {
703:     if (i<n) {
704:       iblock = i / user->mx;
705:       j = iblock*user->mx + ((i+user->mx-1) % user->mx);
706:       MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
707:       j = iblock*user->mx + ((i+1) % user->mx);
708:       MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
709:     }
710:     if (i>=n) {
711:       j = (i - user->mx) % n;
712:       MatSetValues(user->Grad,1,&i,1,&j,&half_hinv,INSERT_VALUES);
713:       j = (j + 2*user->mx) % n;
714:       MatSetValues(user->Grad,1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
715:     }
716:   }

718:   MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
719:   MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);

721:   MatCreate(PETSC_COMM_WORLD,&user->Gradxy[0]);
722:   MatSetSizes(user->Gradxy[0],PETSC_DECIDE,PETSC_DECIDE,n,n);
723:   MatSetFromOptions(user->Gradxy[0]);
724:   MatMPIAIJSetPreallocation(user->Gradxy[0],3,NULL,3,NULL);
725:   MatSeqAIJSetPreallocation(user->Gradxy[0],3,NULL);
726:   MatGetOwnershipRange(user->Gradxy[0],&istart,&iend);

728:   for (i=istart; i<iend; i++) {
729:     iblock = i / user->mx;
730:     j = iblock*user->mx + ((i+user->mx-1) % user->mx);
731:     MatSetValues(user->Gradxy[0],1,&i,1,&j,&half_hinv,INSERT_VALUES);
732:     j = iblock*user->mx + ((i+1) % user->mx);
733:     MatSetValues(user->Gradxy[0],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
734:     MatSetValues(user->Gradxy[0],1,&i,1,&i,&zero,INSERT_VALUES);
735:   }
736:   MatAssemblyBegin(user->Gradxy[0],MAT_FINAL_ASSEMBLY);
737:   MatAssemblyEnd(user->Gradxy[0],MAT_FINAL_ASSEMBLY);

739:   MatCreate(PETSC_COMM_WORLD,&user->Gradxy[1]);
740:   MatSetSizes(user->Gradxy[1],PETSC_DECIDE,PETSC_DECIDE,n,n);
741:   MatSetFromOptions(user->Gradxy[1]);
742:   MatMPIAIJSetPreallocation(user->Gradxy[1],3,NULL,3,NULL);
743:   MatSeqAIJSetPreallocation(user->Gradxy[1],3,NULL);
744:   MatGetOwnershipRange(user->Gradxy[1],&istart,&iend);

746:   for (i=istart; i<iend; i++) {
747:     j = (i + n - user->mx) % n;
748:     MatSetValues(user->Gradxy[1],1,&i,1,&j,&half_hinv,INSERT_VALUES);
749:     j = (j + 2*user->mx) % n;
750:     MatSetValues(user->Gradxy[1],1,&i,1,&j,&neg_half_hinv,INSERT_VALUES);
751:     MatSetValues(user->Gradxy[1],1,&i,1,&i,&zero,INSERT_VALUES);
752:   }
753:   MatAssemblyBegin(user->Gradxy[1],MAT_FINAL_ASSEMBLY);
754:   MatAssemblyEnd(user->Gradxy[1],MAT_FINAL_ASSEMBLY);

756:   /* Generate Div matrix */
757:   MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div);
758:   MatTranspose(user->Gradxy[0],MAT_INITIAL_MATRIX,&user->Divxy[0]);
759:   MatTranspose(user->Gradxy[1],MAT_INITIAL_MATRIX,&user->Divxy[1]);

761:   /* Off-diagonal averaging matrix */
762:   MatCreate(PETSC_COMM_WORLD,&user->M);
763:   MatSetSizes(user->M,PETSC_DECIDE,PETSC_DECIDE,n,n);
764:   MatSetFromOptions(user->M);
765:   MatMPIAIJSetPreallocation(user->M,4,NULL,4,NULL);
766:   MatSeqAIJSetPreallocation(user->M,4,NULL);
767:   MatGetOwnershipRange(user->M,&istart,&iend);

769:   for (i=istart; i<iend; i++) {
770:     /* kron(Id,Av) */
771:     iblock = i / user->mx;
772:     j = iblock*user->mx + ((i+user->mx-1) % user->mx);
773:     MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
774:     j = iblock*user->mx + ((i+1) % user->mx);
775:     MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);

777:     /* kron(Av,Id) */
778:     j = (i + user->mx) % n;
779:     MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
780:     j = (i + n - user->mx) % n;
781:     MatSetValues(user->M,1,&i,1,&j,&quarter,INSERT_VALUES);
782:   }
783:   MatAssemblyBegin(user->M,MAT_FINAL_ASSEMBLY);
784:   MatAssemblyEnd(user->M,MAT_FINAL_ASSEMBLY);

786:   /* Generate 2D grid */
787:   VecCreate(PETSC_COMM_WORLD,&XX);
788:   VecCreate(PETSC_COMM_WORLD,&user->q);
789:   VecSetSizes(XX,PETSC_DECIDE,n);
790:   VecSetSizes(user->q,PETSC_DECIDE,n*user->nt);
791:   VecSetFromOptions(XX);
792:   VecSetFromOptions(user->q);

794:   VecDuplicate(XX,&YY);
795:   VecDuplicate(XX,&XXwork);
796:   VecDuplicate(XX,&YYwork);
797:   VecDuplicate(XX,&user->d);
798:   VecDuplicate(XX,&user->dwork);

800:   VecGetOwnershipRange(XX,&istart,&iend);
801:   for (linear_index=istart; linear_index<iend; linear_index++) {
802:     i = linear_index % user->mx;
803:     j = (linear_index-i)/user->mx;
804:     vx = h*(i+0.5);
805:     vy = h*(j+0.5);
806:     VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
807:     VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
808:   }

810:   VecAssemblyBegin(XX);
811:   VecAssemblyEnd(XX);
812:   VecAssemblyBegin(YY);
813:   VecAssemblyEnd(YY);

815:   /* Compute final density function yT
816:      yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
817:      yT = yT / (h^2*sum(yT)) */
818:   VecCopy(XX,XXwork);
819:   VecCopy(YY,YYwork);

821:   VecShift(XXwork,-0.25);
822:   VecShift(YYwork,-0.25);

824:   VecPointwiseMult(XXwork,XXwork,XXwork);
825:   VecPointwiseMult(YYwork,YYwork,YYwork);

827:   VecCopy(XXwork,user->dwork);
828:   VecAXPY(user->dwork,1.0,YYwork);
829:   VecScale(user->dwork,-30.0);
830:   VecExp(user->dwork);
831:   VecCopy(user->dwork,user->d);

833:   VecCopy(XX,XXwork);
834:   VecCopy(YY,YYwork);

836:   VecShift(XXwork,-0.75);
837:   VecShift(YYwork,-0.75);

839:   VecPointwiseMult(XXwork,XXwork,XXwork);
840:   VecPointwiseMult(YYwork,YYwork,YYwork);

842:   VecCopy(XXwork,user->dwork);
843:   VecAXPY(user->dwork,1.0,YYwork);
844:   VecScale(user->dwork,-30.0);
845:   VecExp(user->dwork);

847:   VecAXPY(user->d,1.0,user->dwork);
848:   VecShift(user->d,1.0);
849:   VecSum(user->d,&sum);
850:   VecScale(user->d,1.0/(h*h*sum));

852:   /* Initial conditions of forward problem */
853:   VecDuplicate(XX,&bc);
854:   VecCopy(XX,XXwork);
855:   VecCopy(YY,YYwork);

857:   VecShift(XXwork,-0.5);
858:   VecShift(YYwork,-0.5);

860:   VecPointwiseMult(XXwork,XXwork,XXwork);
861:   VecPointwiseMult(YYwork,YYwork,YYwork);

863:   VecWAXPY(bc,1.0,XXwork,YYwork);
864:   VecScale(bc,-50.0);
865:   VecExp(bc);
866:   VecShift(bc,1.0);
867:   VecSum(bc,&sum);
868:   VecScale(bc,1.0/(h*h*sum));

870:   /* Create scatter from y to y_1,y_2,...,y_nt */
871:   /*  TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
872:   PetscMalloc1(user->nt*user->mx*user->mx,&user->yi_scatter);
873:   VecCreate(PETSC_COMM_WORLD,&yi);
874:   VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx);
875:   VecSetFromOptions(yi);
876:   VecDuplicateVecs(yi,user->nt,&user->yi);
877:   VecDuplicateVecs(yi,user->nt,&user->yiwork);
878:   VecDuplicateVecs(yi,user->nt,&user->ziwork);
879:   for (i=0; i<user->nt; i++) {
880:     VecGetOwnershipRange(user->yi[i],&lo,&hi);
881:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_yi);
882:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+i*user->mx*user->mx,1,&is_from_y);
883:     VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i]);
884:     ISDestroy(&is_to_yi);
885:     ISDestroy(&is_from_y);
886:   }

888:   /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
889:   /*  TODO: reorder for better parallelism */
890:   PetscMalloc1(user->nt*user->mx*user->mx,&user->uxi_scatter);
891:   PetscMalloc1(user->nt*user->mx*user->mx,&user->uyi_scatter);
892:   PetscMalloc1(user->nt*user->mx*user->mx,&user->ux_scatter);
893:   PetscMalloc1(user->nt*user->mx*user->mx,&user->uy_scatter);
894:   PetscMalloc1(2*user->nt*user->mx*user->mx,&user->ui_scatter);
895:   VecCreate(PETSC_COMM_WORLD,&uxi);
896:   VecCreate(PETSC_COMM_WORLD,&ui);
897:   VecSetSizes(uxi,PETSC_DECIDE,user->mx*user->mx);
898:   VecSetSizes(ui,PETSC_DECIDE,2*user->mx*user->mx);
899:   VecSetFromOptions(uxi);
900:   VecSetFromOptions(ui);
901:   VecDuplicateVecs(uxi,user->nt,&user->uxi);
902:   VecDuplicateVecs(uxi,user->nt,&user->uyi);
903:   VecDuplicateVecs(uxi,user->nt,&user->uxiwork);
904:   VecDuplicateVecs(uxi,user->nt,&user->uyiwork);
905:   VecDuplicateVecs(ui,user->nt,&user->ui);
906:   VecDuplicateVecs(ui,user->nt,&user->uiwork);
907:   for (i=0; i<user->nt; i++) {
908:     VecGetOwnershipRange(user->uxi[i],&lo,&hi);
909:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
910:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
911:     VecScatterCreate(user->u,is_from_u,user->uxi[i],is_to_uxi,&user->uxi_scatter[i]);

913:     ISDestroy(&is_to_uxi);
914:     ISDestroy(&is_from_u);

916:     VecGetOwnershipRange(user->uyi[i],&lo,&hi);
917:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
918:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+(2*i+1)*user->mx*user->mx,1,&is_from_u);
919:     VecScatterCreate(user->u,is_from_u,user->uyi[i],is_to_uyi,&user->uyi_scatter[i]);

921:     ISDestroy(&is_to_uyi);
922:     ISDestroy(&is_from_u);

924:     VecGetOwnershipRange(user->uxi[i],&lo,&hi);
925:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
926:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_from_u);
927:     VecScatterCreate(user->ui[i],is_from_u,user->uxi[i],is_to_uxi,&user->ux_scatter[i]);

929:     ISDestroy(&is_to_uxi);
930:     ISDestroy(&is_from_u);

932:     VecGetOwnershipRange(user->uyi[i],&lo,&hi);
933:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uyi);
934:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+user->mx*user->mx,1,&is_from_u);
935:     VecScatterCreate(user->ui[i],is_from_u,user->uyi[i],is_to_uyi,&user->uy_scatter[i]);

937:     ISDestroy(&is_to_uyi);
938:     ISDestroy(&is_from_u);

940:     VecGetOwnershipRange(user->ui[i],&lo,&hi);
941:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_uxi);
942:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+2*i*user->mx*user->mx,1,&is_from_u);
943:     VecScatterCreate(user->u,is_from_u,user->ui[i],is_to_uxi,&user->ui_scatter[i]);

945:     ISDestroy(&is_to_uxi);
946:     ISDestroy(&is_from_u);
947:   }

949:   /* RHS of forward problem */
950:   MatMult(user->M,bc,user->yiwork[0]);
951:   for (i=1; i<user->nt; i++) {
952:     VecSet(user->yiwork[i],0.0);
953:   }
954:   Gather_yi(user->q,user->yiwork,user->yi_scatter,user->nt);

956:   /* Compute true velocity field utrue */
957:   VecDuplicate(user->u,&user->utrue);
958:   for (i=0; i<user->nt; i++) {
959:     VecCopy(YY,user->uxi[i]);
960:     VecScale(user->uxi[i],150.0*i*user->ht);
961:     VecCopy(XX,user->uyi[i]);
962:     VecShift(user->uyi[i],-10.0);
963:     VecScale(user->uyi[i],15.0*i*user->ht);
964:   }
965:   Gather_uxi_uyi(user->utrue,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

967:   /* Initial guess and reference model */
968:   VecDuplicate(user->utrue,&user->ur);
969:   for (i=0; i<user->nt; i++) {
970:     VecCopy(XX,user->uxi[i]);
971:     VecShift(user->uxi[i],i*user->ht);
972:     VecCopy(YY,user->uyi[i]);
973:     VecShift(user->uyi[i],-i*user->ht);
974:   }
975:   Gather_uxi_uyi(user->ur,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);

977:   /* Generate regularization matrix L */
978:   MatCreate(PETSC_COMM_WORLD,&user->LT);
979:   MatSetSizes(user->LT,PETSC_DECIDE,PETSC_DECIDE,2*n*user->nt,n*user->nt);
980:   MatSetFromOptions(user->LT);
981:   MatMPIAIJSetPreallocation(user->LT,1,NULL,1,NULL);
982:   MatSeqAIJSetPreallocation(user->LT,1,NULL);
983:   MatGetOwnershipRange(user->LT,&istart,&iend);

985:   for (i=istart; i<iend; i++) {
986:     iblock = (i+n) / (2*n);
987:     j = i - iblock*n;
988:     MatSetValues(user->LT,1,&i,1,&j,&user->gamma,INSERT_VALUES);
989:   }

991:   MatAssemblyBegin(user->LT,MAT_FINAL_ASSEMBLY);
992:   MatAssemblyEnd(user->LT,MAT_FINAL_ASSEMBLY);

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

996:   /* Build work vectors and matrices */
997:   VecCreate(PETSC_COMM_WORLD,&user->lwork);
998:   VecSetType(user->lwork,VECMPI);
999:   VecSetSizes(user->lwork,PETSC_DECIDE,user->m);
1000:   VecSetFromOptions(user->lwork);

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

1004:   VecDuplicate(user->y,&user->ywork);
1005:   VecDuplicate(user->u,&user->uwork);
1006:   VecDuplicate(user->u,&user->vwork);
1007:   VecDuplicate(user->u,&user->js_diag);
1008:   VecDuplicate(user->c,&user->cwork);

1010:   /* Create matrix-free shell user->Js for computing A*x */
1011:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->Js);
1012:   MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1013:   MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);
1014:   MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);
1015:   MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);

1017:   /* Diagonal blocks of user->Js */
1018:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,n,n,user,&user->JsBlock);
1019:   MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult);
1020:   MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockMultTranspose);

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

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

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

1039:   /* Build matrices for SOR preconditioner */
1040:   Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
1041:   PetscMalloc1(5*n,&user->C);
1042:   PetscMalloc1(2*n,&user->Cwork);
1043:   for (i=0; i<user->nt; i++) {
1044:     MatDuplicate(user->Divxy[0],MAT_COPY_VALUES,&user->C[i]);
1045:     MatDuplicate(user->Divxy[1],MAT_COPY_VALUES,&user->Cwork[i]);

1047:     MatDiagonalScale(user->C[i],NULL,user->uxi[i]);
1048:     MatDiagonalScale(user->Cwork[i],NULL,user->uyi[i]);
1049:     MatAXPY(user->C[i],1.0,user->Cwork[i],DIFFERENT_NONZERO_PATTERN);
1050:     MatScale(user->C[i],user->ht);
1051:     MatShift(user->C[i],1.0);
1052:   }

1054:   /* Solver options and tolerances */
1055:   KSPCreate(PETSC_COMM_WORLD,&user->solver);
1056:   KSPSetType(user->solver,KSPGMRES);
1057:   KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec);
1058:   KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500);
1059:   /* KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500); */
1060:   KSPGetPC(user->solver,&user->prec);
1061:   PCSetType(user->prec,PCSHELL);

1063:   PCShellSetApply(user->prec,StateMatBlockPrecMult);
1064:   PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMultTranspose);
1065:   PCShellSetContext(user->prec,user);

1067:   /* Compute true state function yt given ut */
1068:   VecCreate(PETSC_COMM_WORLD,&user->ytrue);
1069:   VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt);
1070:   VecSetFromOptions(user->ytrue);
1071:   user->c_formed = PETSC_TRUE;
1072:   VecCopy(user->utrue,user->u); /*  Set u=utrue temporarily for StateMatInv */
1073:   VecSet(user->ytrue,0.0); /*  Initial guess */
1074:   StateMatInvMult(user->Js,user->q,user->ytrue);
1075:   VecCopy(user->ur,user->u); /*  Reset u=ur */

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

1080:   /* Data discretization */
1081:   MatCreate(PETSC_COMM_WORLD,&user->Q);
1082:   MatSetSizes(user->Q,PETSC_DECIDE,PETSC_DECIDE,user->mx*user->mx,user->m);
1083:   MatSetFromOptions(user->Q);
1084:   MatMPIAIJSetPreallocation(user->Q,0,NULL,1,NULL);
1085:   MatSeqAIJSetPreallocation(user->Q,1,NULL);

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

1089:   for (i=istart; i<iend; i++) {
1090:     j = i + user->m - user->mx*user->mx;
1091:     MatSetValues(user->Q,1,&i,1,&j,&one,INSERT_VALUES);
1092:   }

1094:   MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1095:   MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);

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

1099:   VecDestroy(&XX);
1100:   VecDestroy(&YY);
1101:   VecDestroy(&XXwork);
1102:   VecDestroy(&YYwork);
1103:   VecDestroy(&yi);
1104:   VecDestroy(&uxi);
1105:   VecDestroy(&ui);
1106:   VecDestroy(&bc);

1108:   /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1109:   KSPSetFromOptions(user->solver);
1110:   return 0;
1111: }

1113: PetscErrorCode HyperbolicDestroy(AppCtx *user)
1114: {
1115:   PetscInt       i;

1117:   MatDestroy(&user->Q);
1118:   MatDestroy(&user->QT);
1119:   MatDestroy(&user->Div);
1120:   MatDestroy(&user->Divwork);
1121:   MatDestroy(&user->Grad);
1122:   MatDestroy(&user->L);
1123:   MatDestroy(&user->LT);
1124:   KSPDestroy(&user->solver);
1125:   MatDestroy(&user->Js);
1126:   MatDestroy(&user->Jd);
1127:   MatDestroy(&user->JsBlockPrec);
1128:   MatDestroy(&user->JsInv);
1129:   MatDestroy(&user->JsBlock);
1130:   MatDestroy(&user->Divxy[0]);
1131:   MatDestroy(&user->Divxy[1]);
1132:   MatDestroy(&user->Gradxy[0]);
1133:   MatDestroy(&user->Gradxy[1]);
1134:   MatDestroy(&user->M);
1135:   for (i=0; i<user->nt; i++) {
1136:     MatDestroy(&user->C[i]);
1137:     MatDestroy(&user->Cwork[i]);
1138:   }
1139:   PetscFree(user->C);
1140:   PetscFree(user->Cwork);
1141:   VecDestroy(&user->u);
1142:   VecDestroy(&user->uwork);
1143:   VecDestroy(&user->vwork);
1144:   VecDestroy(&user->utrue);
1145:   VecDestroy(&user->y);
1146:   VecDestroy(&user->ywork);
1147:   VecDestroy(&user->ytrue);
1148:   VecDestroyVecs(user->nt,&user->yi);
1149:   VecDestroyVecs(user->nt,&user->yiwork);
1150:   VecDestroyVecs(user->nt,&user->ziwork);
1151:   VecDestroyVecs(user->nt,&user->uxi);
1152:   VecDestroyVecs(user->nt,&user->uyi);
1153:   VecDestroyVecs(user->nt,&user->uxiwork);
1154:   VecDestroyVecs(user->nt,&user->uyiwork);
1155:   VecDestroyVecs(user->nt,&user->ui);
1156:   VecDestroyVecs(user->nt,&user->uiwork);
1157:   VecDestroy(&user->c);
1158:   VecDestroy(&user->cwork);
1159:   VecDestroy(&user->ur);
1160:   VecDestroy(&user->q);
1161:   VecDestroy(&user->d);
1162:   VecDestroy(&user->dwork);
1163:   VecDestroy(&user->lwork);
1164:   VecDestroy(&user->js_diag);
1165:   ISDestroy(&user->s_is);
1166:   ISDestroy(&user->d_is);
1167:   VecScatterDestroy(&user->state_scatter);
1168:   VecScatterDestroy(&user->design_scatter);
1169:   for (i=0; i<user->nt; i++) {
1170:     VecScatterDestroy(&user->uxi_scatter[i]);
1171:     VecScatterDestroy(&user->uyi_scatter[i]);
1172:     VecScatterDestroy(&user->ux_scatter[i]);
1173:     VecScatterDestroy(&user->uy_scatter[i]);
1174:     VecScatterDestroy(&user->ui_scatter[i]);
1175:     VecScatterDestroy(&user->yi_scatter[i]);
1176:   }
1177:   PetscFree(user->uxi_scatter);
1178:   PetscFree(user->uyi_scatter);
1179:   PetscFree(user->ux_scatter);
1180:   PetscFree(user->uy_scatter);
1181:   PetscFree(user->ui_scatter);
1182:   PetscFree(user->yi_scatter);
1183:   return 0;
1184: }

1186: PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1187: {
1188:   Vec            X;
1189:   PetscReal      unorm,ynorm;
1190:   AppCtx         *user = (AppCtx*)ptr;

1192:   TaoGetSolution(tao,&X);
1193:   Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1194:   VecAXPY(user->ywork,-1.0,user->ytrue);
1195:   VecAXPY(user->uwork,-1.0,user->utrue);
1196:   VecNorm(user->uwork,NORM_2,&unorm);
1197:   VecNorm(user->ywork,NORM_2,&ynorm);
1198:   PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1199:   return 0;
1200: }

1202: /*TEST

1204:    build:
1205:       requires: !complex

1207:    test:
1208:       requires: !single
1209:       args: -tao_cmonitor -tao_max_funcs 10 -tao_type lcl -tao_gatol 1.e-5

1211:    test:
1212:       suffix: guess_pod
1213:       requires: !single
1214:       args: -tao_cmonitor -tao_max_funcs 10 -tao_type lcl -ksp_guess_type pod -tao_gatol 1.e-5

1216: TEST*/