Actual source code: hyperbolic.c

  1: #include <petsctao.h>

  3: typedef struct {
  4:   PetscInt    n;     /*  Number of variables */
  5:   PetscInt    m;     /*  Number of constraints */
  6:   PetscInt    mx;    /*  grid points in each direction */
  7:   PetscInt    nt;    /*  Number of time steps */
  8:   PetscInt    ndata; /*  Number of data points per sample */
  9:   IS          s_is;
 10:   IS          d_is;
 11:   VecScatter  state_scatter;
 12:   VecScatter  design_scatter;
 13:   VecScatter *uxi_scatter, *uyi_scatter, *ux_scatter, *uy_scatter, *ui_scatter;
 14:   VecScatter *yi_scatter;

 16:   Mat       Js, Jd, JsBlockPrec, JsInv, JsBlock;
 17:   PetscBool jformed, c_formed;

 19:   PetscReal alpha; /*  Regularization parameter */
 20:   PetscReal gamma;
 21:   PetscReal ht; /*  Time step */
 22:   PetscReal T;  /*  Final time */
 23:   Mat       Q, QT;
 24:   Mat       L, LT;
 25:   Mat       Div, Divwork, Divxy[2];
 26:   Mat       Grad, Gradxy[2];
 27:   Mat       M;
 28:   Mat      *C, *Cwork;
 29:   /* Mat Hs,Hd,Hsd; */
 30:   Vec q;
 31:   Vec ur; /*  reference */

 33:   Vec d;
 34:   Vec dwork;

 36:   Vec  y; /*  state variables */
 37:   Vec  ywork;
 38:   Vec  ytrue;
 39:   Vec *yi, *yiwork, *ziwork;
 40:   Vec *uxi, *uyi, *uxiwork, *uyiwork, *ui, *uiwork;

 42:   Vec u; /*  design variables */
 43:   Vec uwork, vwork;
 44:   Vec utrue;

 46:   Vec js_diag;

 48:   Vec c; /*  constraint vector */
 49:   Vec cwork;

 51:   Vec lwork;

 53:   KSP      solver;
 54:   PC       prec;
 55:   PetscInt block_index;

 57:   PetscInt ksp_its;
 58:   PetscInt ksp_its_initial;
 59: } AppCtx;

 61: PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
 62: PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
 63: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 64: PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void *);
 65: PetscErrorCode FormJacobianDesign(Tao, Vec, Mat, void *);
 66: PetscErrorCode FormConstraints(Tao, Vec, Vec, void *);
 67: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
 68: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
 69: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
 70: PetscErrorCode HyperbolicInitialize(AppCtx *user);
 71: PetscErrorCode HyperbolicDestroy(AppCtx *user);
 72: PetscErrorCode HyperbolicMonitor(Tao, void *);

 74: PetscErrorCode StateMatMult(Mat, Vec, Vec);
 75: PetscErrorCode StateMatBlockMult(Mat, Vec, Vec);
 76: PetscErrorCode StateMatBlockMultTranspose(Mat, Vec, Vec);
 77: PetscErrorCode StateMatMultTranspose(Mat, Vec, Vec);
 78: PetscErrorCode StateMatGetDiagonal(Mat, Vec);
 79: PetscErrorCode StateMatDuplicate(Mat, MatDuplicateOption, Mat *);
 80: PetscErrorCode StateMatInvMult(Mat, Vec, Vec);
 81: PetscErrorCode StateMatInvTransposeMult(Mat, Vec, Vec);
 82: PetscErrorCode StateMatBlockPrecMult(PC, Vec, Vec);

 84: PetscErrorCode DesignMatMult(Mat, Vec, Vec);
 85: PetscErrorCode DesignMatMultTranspose(Mat, Vec, Vec);

 87: PetscErrorCode Scatter_yi(Vec, Vec *, VecScatter *, PetscInt); /*  y to y1,y2,...,y_nt */
 88: PetscErrorCode Gather_yi(Vec, Vec *, VecScatter *, PetscInt);
 89: PetscErrorCode Scatter_uxi_uyi(Vec, Vec *, VecScatter *, Vec *, VecScatter *, PetscInt); /*  u to ux_1,uy_1,ux_2,uy_2,...,u */
 90: PetscErrorCode Gather_uxi_uyi(Vec, Vec *, VecScatter *, Vec *, VecScatter *, PetscInt);

 92: static char help[] = "";

 94: int main(int argc, char **argv)
 95: {
 96:   Vec      x, x0;
 97:   Tao      tao;
 98:   AppCtx   user;
 99:   IS       is_allstate, is_alldesign;
100:   PetscInt lo, hi, hi2, lo2, ksp_old;
101:   PetscInt ntests = 1;
102:   PetscInt i;
103: #if defined(PETSC_USE_LOG)
104:   PetscLogStage stages[1];
105: #endif

108:   PetscInitialize(&argc, &argv, (char *)0, help);
109:   user.mx = 32;
110:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "hyperbolic example", NULL);
111:   PetscOptionsInt("-mx", "Number of grid points in each direction", "", user.mx, &user.mx, NULL);
112:   user.nt = 16;
113:   PetscOptionsInt("-nt", "Number of time steps", "", user.nt, &user.nt, NULL);
114:   user.ndata = 64;
115:   PetscOptionsInt("-ndata", "Numbers of data points per sample", "", user.ndata, &user.ndata, NULL);
116:   user.alpha = 10.0;
117:   PetscOptionsReal("-alpha", "Regularization parameter", "", user.alpha, &user.alpha, NULL);
118:   user.T = 1.0 / 32.0;
119:   PetscOptionsReal("-Tfinal", "Final time", "", user.T, &user.T, NULL);
120:   PetscOptionsInt("-ntests", "Number of times to repeat TaoSolve", "", ntests, &ntests, NULL);
121:   PetscOptionsEnd();

123:   user.m     = user.mx * user.mx * user.nt;     /*  number of constraints */
124:   user.n     = user.mx * user.mx * 3 * user.nt; /*  number of variables */
125:   user.ht    = user.T / user.nt;                /*  Time step */
126:   user.gamma = user.T * user.ht / (user.mx * user.mx);

128:   VecCreate(PETSC_COMM_WORLD, &user.u);
129:   VecCreate(PETSC_COMM_WORLD, &user.y);
130:   VecCreate(PETSC_COMM_WORLD, &user.c);
131:   VecSetSizes(user.u, PETSC_DECIDE, user.n - user.m);
132:   VecSetSizes(user.y, PETSC_DECIDE, user.m);
133:   VecSetSizes(user.c, PETSC_DECIDE, user.m);
134:   VecSetFromOptions(user.u);
135:   VecSetFromOptions(user.y);
136:   VecSetFromOptions(user.c);

138:   /* Create scatters for reduced spaces.
139:      If the state vector y and design vector u are partitioned as
140:      [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
141:      then the solution vector x is organized as
142:      [y_1; u_1; y_2; u_2; ...; y_np; u_np].
143:      The index sets user.s_is and user.d_is correspond to the indices of the
144:      state and design variables owned by the current processor.
145:   */
146:   VecCreate(PETSC_COMM_WORLD, &x);

148:   VecGetOwnershipRange(user.y, &lo, &hi);
149:   VecGetOwnershipRange(user.u, &lo2, &hi2);

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

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

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

160:   VecScatterCreate(x, user.s_is, user.y, is_allstate, &user.state_scatter);
161:   VecScatterCreate(x, user.d_is, user.u, is_alldesign, &user.design_scatter);
162:   ISDestroy(&is_alldesign);
163:   ISDestroy(&is_allstate);

165:   /* Create TAO solver and set desired solution method */
166:   TaoCreate(PETSC_COMM_WORLD, &tao);
167:   TaoSetType(tao, TAOLCL);

169:   /* Set up initial vectors and matrices */
170:   HyperbolicInitialize(&user);

172:   Gather(x, user.y, user.state_scatter, user.u, user.design_scatter);
173:   VecDuplicate(x, &x0);
174:   VecCopy(x, x0);

176:   /* Set solution vector with an initial guess */
177:   TaoSetSolution(tao, x);
178:   TaoSetObjective(tao, FormFunction, &user);
179:   TaoSetGradient(tao, NULL, FormGradient, &user);
180:   TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user);
181:   TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, &user);
182:   TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user);
183:   TaoSetFromOptions(tao);
184:   TaoSetStateDesignIS(tao, user.s_is, user.d_is);

186:   /* SOLVE THE APPLICATION */
187:   PetscLogStageRegister("Trials", &stages[0]);
188:   PetscLogStagePush(stages[0]);
189:   user.ksp_its_initial = user.ksp_its;
190:   ksp_old              = user.ksp_its;
191:   for (i = 0; i < ntests; i++) {
192:     TaoSolve(tao);
193:     PetscPrintf(PETSC_COMM_WORLD, "KSP Iterations = %" PetscInt_FMT "\n", user.ksp_its - ksp_old);
194:     VecCopy(x0, x);
195:     TaoSetSolution(tao, x);
196:   }
197:   PetscLogStagePop();
198:   PetscBarrier((PetscObject)x);
199:   PetscPrintf(PETSC_COMM_WORLD, "KSP iterations within initialization: ");
200:   PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its_initial);
201:   PetscPrintf(PETSC_COMM_WORLD, "Total KSP iterations over %" PetscInt_FMT " trial(s): ", ntests);
202:   PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its);
203:   PetscPrintf(PETSC_COMM_WORLD, "KSP iterations per trial: ");
204:   PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", (user.ksp_its - user.ksp_its_initial) / ntests);

206:   TaoDestroy(&tao);
207:   VecDestroy(&x);
208:   VecDestroy(&x0);
209:   HyperbolicDestroy(&user);
210:   PetscFinalize();
211:   return 0;
212: }
213: /* ------------------------------------------------------------------- */
214: /*
215:    dwork = Qy - d
216:    lwork = L*(u-ur).^2
217:    f = 1/2 * (dwork.dork + alpha*y.lwork)
218: */
219: PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr)
220: {
221:   PetscReal d1 = 0, d2 = 0;
222:   AppCtx   *user = (AppCtx *)ptr;

224:   Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter);
225:   MatMult(user->Q, user->y, user->dwork);
226:   VecAXPY(user->dwork, -1.0, user->d);
227:   VecDot(user->dwork, user->dwork, &d1);

229:   VecWAXPY(user->uwork, -1.0, user->ur, user->u);
230:   VecPointwiseMult(user->uwork, user->uwork, user->uwork);
231:   MatMult(user->L, user->uwork, user->lwork);
232:   VecDot(user->y, user->lwork, &d2);
233:   *f = 0.5 * (d1 + user->alpha * d2);
234:   return 0;
235: }

237: /* ------------------------------------------------------------------- */
238: /*
239:     state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2
240:     design: g_d = alpha*(L'y).*(u-ur)
241: */
242: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr)
243: {
244:   AppCtx *user = (AppCtx *)ptr;

246:   Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter);
247:   MatMult(user->Q, user->y, user->dwork);
248:   VecAXPY(user->dwork, -1.0, user->d);

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

252:   MatMult(user->LT, user->y, user->uwork);
253:   VecWAXPY(user->vwork, -1.0, user->ur, user->u);
254:   VecPointwiseMult(user->uwork, user->vwork, user->uwork);
255:   VecScale(user->uwork, user->alpha);

257:   VecPointwiseMult(user->vwork, user->vwork, user->vwork);
258:   MatMult(user->L, user->vwork, user->lwork);
259:   VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork);

261:   Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter);
262:   return 0;
263: }

265: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
266: {
267:   PetscReal d1, d2;
268:   AppCtx   *user = (AppCtx *)ptr;

270:   Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter);
271:   MatMult(user->Q, user->y, user->dwork);
272:   VecAXPY(user->dwork, -1.0, user->d);

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

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

278:   MatMult(user->LT, user->y, user->uwork);
279:   VecWAXPY(user->vwork, -1.0, user->ur, user->u);
280:   VecPointwiseMult(user->uwork, user->vwork, user->uwork);
281:   VecScale(user->uwork, user->alpha);

283:   VecPointwiseMult(user->vwork, user->vwork, user->vwork);
284:   MatMult(user->L, user->vwork, user->lwork);
285:   VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork);

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

289:   *f = 0.5 * (d1 + user->alpha * d2);
290:   Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter);
291:   return 0;
292: }

294: /* ------------------------------------------------------------------- */
295: /* A
296: MatShell object
297: */
298: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
299: {
300:   PetscInt i;
301:   AppCtx  *user = (AppCtx *)ptr;

303:   Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter);
304:   Scatter_yi(user->u, user->ui, user->ui_scatter, user->nt);
305:   Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt);
306:   for (i = 0; i < user->nt; i++) {
307:     MatCopy(user->Divxy[0], user->C[i], SUBSET_NONZERO_PATTERN);
308:     MatCopy(user->Divxy[1], user->Cwork[i], SAME_NONZERO_PATTERN);

310:     MatDiagonalScale(user->C[i], NULL, user->uxi[i]);
311:     MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i]);
312:     MatAXPY(user->C[i], 1.0, user->Cwork[i], SUBSET_NONZERO_PATTERN);
313:     MatScale(user->C[i], user->ht);
314:     MatShift(user->C[i], 1.0);
315:   }
316:   return 0;
317: }

319: /* ------------------------------------------------------------------- */
320: /* B */
321: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
322: {
323:   AppCtx *user = (AppCtx *)ptr;

325:   Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter);
326:   return 0;
327: }

329: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
330: {
331:   PetscInt i;
332:   AppCtx  *user;

334:   MatShellGetContext(J_shell, &user);
335:   Scatter_yi(X, user->yi, user->yi_scatter, user->nt);
336:   user->block_index = 0;
337:   MatMult(user->JsBlock, user->yi[0], user->yiwork[0]);

339:   for (i = 1; i < user->nt; i++) {
340:     user->block_index = i;
341:     MatMult(user->JsBlock, user->yi[i], user->yiwork[i]);
342:     MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1]);
343:     VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1]);
344:   }
345:   Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt);
346:   return 0;
347: }

349: PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
350: {
351:   PetscInt i;
352:   AppCtx  *user;

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

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

364:   i                 = user->nt - 1;
365:   user->block_index = i;
366:   MatMultTranspose(user->JsBlock, user->yi[i], user->yiwork[i]);
367:   Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt);
368:   return 0;
369: }

371: PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
372: {
373:   PetscInt i;
374:   AppCtx  *user;

376:   MatShellGetContext(J_shell, &user);
377:   i = user->block_index;
378:   VecPointwiseMult(user->uxiwork[i], X, user->uxi[i]);
379:   VecPointwiseMult(user->uyiwork[i], X, user->uyi[i]);
380:   Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]);
381:   MatMult(user->Div, user->uiwork[i], Y);
382:   VecAYPX(Y, user->ht, X);
383:   return 0;
384: }

386: PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y)
387: {
388:   PetscInt i;
389:   AppCtx  *user;

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

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

407:   MatShellGetContext(J_shell, &user);
408:   Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt);
409:   Scatter_uxi_uyi(X, user->uxiwork, user->uxi_scatter, user->uyiwork, user->uyi_scatter, user->nt);
410:   for (i = 0; i < user->nt; i++) {
411:     VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i]);
412:     VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i]);
413:     Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]);
414:     MatMult(user->Div, user->uiwork[i], user->ziwork[i]);
415:     VecScale(user->ziwork[i], user->ht);
416:   }
417:   Gather_yi(Y, user->ziwork, user->yi_scatter, user->nt);
418:   return 0;
419: }

421: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
422: {
423:   PetscInt i;
424:   AppCtx  *user;

426:   MatShellGetContext(J_shell, &user);
427:   Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt);
428:   Scatter_yi(X, user->yiwork, user->yi_scatter, user->nt);
429:   for (i = 0; i < user->nt; i++) {
430:     MatMult(user->Grad, user->yiwork[i], user->uiwork[i]);
431:     Scatter(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]);
432:     VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i]);
433:     VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i]);
434:     Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]);
435:     VecScale(user->uiwork[i], user->ht);
436:   }
437:   Gather_yi(Y, user->uiwork, user->ui_scatter, user->nt);
438:   return 0;
439: }

441: PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
442: {
443:   PetscInt i;
444:   AppCtx  *user;

446:   PCShellGetContext(PC_shell, &user);
447:   i = user->block_index;
448:   if (user->c_formed) {
449:     MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y);
450:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed");
451:   return 0;
452: }

454: PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
455: {
456:   PetscInt i;
457:   AppCtx  *user;

459:   PCShellGetContext(PC_shell, &user);

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

468: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
469: {
470:   AppCtx  *user;
471:   PetscInt its, i;

473:   MatShellGetContext(J_shell, &user);

475:   if (Y == user->ytrue) {
476:     /* First solve is done using true solution to set up problem */
477:     KSPSetTolerances(user->solver, 1e-4, 1e-20, PETSC_DEFAULT, PETSC_DEFAULT);
478:   } else {
479:     KSPSetTolerances(user->solver, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
480:   }
481:   Scatter_yi(X, user->yi, user->yi_scatter, user->nt);
482:   Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt);
483:   Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt);

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

488:   KSPGetIterationNumber(user->solver, &its);
489:   user->ksp_its = user->ksp_its + its;
490:   for (i = 1; i < user->nt; i++) {
491:     MatMult(user->M, user->yiwork[i - 1], user->ziwork[i - 1]);
492:     VecAXPY(user->yi[i], 1.0, user->ziwork[i - 1]);
493:     user->block_index = i;
494:     KSPSolve(user->solver, user->yi[i], user->yiwork[i]);

496:     KSPGetIterationNumber(user->solver, &its);
497:     user->ksp_its = user->ksp_its + its;
498:   }
499:   Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt);
500:   return 0;
501: }

503: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
504: {
505:   AppCtx  *user;
506:   PetscInt its, i;

508:   MatShellGetContext(J_shell, &user);

510:   Scatter_yi(X, user->yi, user->yi_scatter, user->nt);
511:   Scatter_yi(Y, user->yiwork, user->yi_scatter, user->nt);
512:   Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt);

514:   i                 = user->nt - 1;
515:   user->block_index = i;
516:   KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i]);

518:   KSPGetIterationNumber(user->solver, &its);
519:   user->ksp_its = user->ksp_its + its;

521:   for (i = user->nt - 2; i >= 0; i--) {
522:     MatMult(user->M, user->yiwork[i + 1], user->ziwork[i + 1]);
523:     VecAXPY(user->yi[i], 1.0, user->ziwork[i + 1]);
524:     user->block_index = i;
525:     KSPSolveTranspose(user->solver, user->yi[i], user->yiwork[i]);

527:     KSPGetIterationNumber(user->solver, &its);
528:     user->ksp_its = user->ksp_its + its;
529:   }
530:   Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt);
531:   return 0;
532: }

534: PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
535: {
536:   AppCtx *user;

538:   MatShellGetContext(J_shell, &user);

540:   MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, new_shell);
541:   MatShellSetOperation(*new_shell, MATOP_MULT, (void (*)(void))StateMatMult);
542:   MatShellSetOperation(*new_shell, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate);
543:   MatShellSetOperation(*new_shell, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose);
544:   MatShellSetOperation(*new_shell, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal);
545:   return 0;
546: }

548: PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
549: {
550:   AppCtx *user;

552:   MatShellGetContext(J_shell, &user);
553:   VecCopy(user->js_diag, X);
554:   return 0;
555: }

557: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
558: {
559:   /* con = Ay - q, A = [C(u1)  0     0     ...   0;
560:                          -M  C(u2)   0     ...   0;
561:                           0   -M   C(u3)   ...   0;
562:                                       ...         ;
563:                           0    ...      -M C(u_nt)]
564:      C(u) = eye + ht*Div*[diag(u1); diag(u2)]       */
565:   PetscInt i;
566:   AppCtx  *user = (AppCtx *)ptr;

568:   Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter);
569:   Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt);
570:   Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt);

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

575:   for (i = 1; i < user->nt; i++) {
576:     user->block_index = i;
577:     MatMult(user->JsBlock, user->yi[i], user->yiwork[i]);
578:     MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1]);
579:     VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1]);
580:   }

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

585:   return 0;
586: }

588: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
589: {
590:   VecScatterBegin(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD);
591:   VecScatterEnd(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD);
592:   VecScatterBegin(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD);
593:   VecScatterEnd(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD);
594:   return 0;
595: }

597: PetscErrorCode Scatter_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
598: {
599:   PetscInt i;

601:   for (i = 0; i < nt; i++) {
602:     VecScatterBegin(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD);
603:     VecScatterEnd(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD);
604:     VecScatterBegin(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD);
605:     VecScatterEnd(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD);
606:   }
607:   return 0;
608: }

610: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
611: {
612:   VecScatterBegin(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE);
613:   VecScatterEnd(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE);
614:   VecScatterBegin(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE);
615:   VecScatterEnd(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE);
616:   return 0;
617: }

619: PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
620: {
621:   PetscInt i;

623:   for (i = 0; i < nt; i++) {
624:     VecScatterBegin(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE);
625:     VecScatterEnd(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE);
626:     VecScatterBegin(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE);
627:     VecScatterEnd(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE);
628:   }
629:   return 0;
630: }

632: PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
633: {
634:   PetscInt i;

636:   for (i = 0; i < nt; i++) {
637:     VecScatterBegin(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD);
638:     VecScatterEnd(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD);
639:   }
640:   return 0;
641: }

643: PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
644: {
645:   PetscInt i;

647:   for (i = 0; i < nt; i++) {
648:     VecScatterBegin(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE);
649:     VecScatterEnd(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE);
650:   }
651:   return 0;
652: }

654: PetscErrorCode HyperbolicInitialize(AppCtx *user)
655: {
656:   PetscInt    n, i, j, linear_index, istart, iend, iblock, lo, hi;
657:   Vec         XX, YY, XXwork, YYwork, yi, uxi, ui, bc;
658:   PetscReal   h, sum;
659:   PetscScalar hinv, neg_hinv, quarter = 0.25, one = 1.0, half_hinv, neg_half_hinv;
660:   PetscScalar vx, vy, zero = 0.0;
661:   IS          is_from_y, is_to_yi, is_from_u, is_to_uxi, is_to_uyi;

663:   user->jformed  = PETSC_FALSE;
664:   user->c_formed = PETSC_FALSE;

666:   user->ksp_its         = 0;
667:   user->ksp_its_initial = 0;

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

671:   h             = 1.0 / user->mx;
672:   hinv          = user->mx;
673:   neg_hinv      = -hinv;
674:   half_hinv     = hinv / 2.0;
675:   neg_half_hinv = neg_hinv / 2.0;

677:   /* Generate Grad matrix */
678:   MatCreate(PETSC_COMM_WORLD, &user->Grad);
679:   MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, 2 * n, n);
680:   MatSetFromOptions(user->Grad);
681:   MatMPIAIJSetPreallocation(user->Grad, 3, NULL, 3, NULL);
682:   MatSeqAIJSetPreallocation(user->Grad, 3, NULL);
683:   MatGetOwnershipRange(user->Grad, &istart, &iend);

685:   for (i = istart; i < iend; i++) {
686:     if (i < n) {
687:       iblock = i / user->mx;
688:       j      = iblock * user->mx + ((i + user->mx - 1) % user->mx);
689:       MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES);
690:       j = iblock * user->mx + ((i + 1) % user->mx);
691:       MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES);
692:     }
693:     if (i >= n) {
694:       j = (i - user->mx) % n;
695:       MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES);
696:       j = (j + 2 * user->mx) % n;
697:       MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES);
698:     }
699:   }

701:   MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY);
702:   MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY);

704:   MatCreate(PETSC_COMM_WORLD, &user->Gradxy[0]);
705:   MatSetSizes(user->Gradxy[0], PETSC_DECIDE, PETSC_DECIDE, n, n);
706:   MatSetFromOptions(user->Gradxy[0]);
707:   MatMPIAIJSetPreallocation(user->Gradxy[0], 3, NULL, 3, NULL);
708:   MatSeqAIJSetPreallocation(user->Gradxy[0], 3, NULL);
709:   MatGetOwnershipRange(user->Gradxy[0], &istart, &iend);

711:   for (i = istart; i < iend; i++) {
712:     iblock = i / user->mx;
713:     j      = iblock * user->mx + ((i + user->mx - 1) % user->mx);
714:     MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &half_hinv, INSERT_VALUES);
715:     j = iblock * user->mx + ((i + 1) % user->mx);
716:     MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES);
717:     MatSetValues(user->Gradxy[0], 1, &i, 1, &i, &zero, INSERT_VALUES);
718:   }
719:   MatAssemblyBegin(user->Gradxy[0], MAT_FINAL_ASSEMBLY);
720:   MatAssemblyEnd(user->Gradxy[0], MAT_FINAL_ASSEMBLY);

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

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

739:   /* Generate Div matrix */
740:   MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div);
741:   MatTranspose(user->Gradxy[0], MAT_INITIAL_MATRIX, &user->Divxy[0]);
742:   MatTranspose(user->Gradxy[1], MAT_INITIAL_MATRIX, &user->Divxy[1]);

744:   /* Off-diagonal averaging matrix */
745:   MatCreate(PETSC_COMM_WORLD, &user->M);
746:   MatSetSizes(user->M, PETSC_DECIDE, PETSC_DECIDE, n, n);
747:   MatSetFromOptions(user->M);
748:   MatMPIAIJSetPreallocation(user->M, 4, NULL, 4, NULL);
749:   MatSeqAIJSetPreallocation(user->M, 4, NULL);
750:   MatGetOwnershipRange(user->M, &istart, &iend);

752:   for (i = istart; i < iend; i++) {
753:     /* kron(Id,Av) */
754:     iblock = i / user->mx;
755:     j      = iblock * user->mx + ((i + user->mx - 1) % user->mx);
756:     MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES);
757:     j = iblock * user->mx + ((i + 1) % user->mx);
758:     MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES);

760:     /* kron(Av,Id) */
761:     j = (i + user->mx) % n;
762:     MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES);
763:     j = (i + n - user->mx) % n;
764:     MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES);
765:   }
766:   MatAssemblyBegin(user->M, MAT_FINAL_ASSEMBLY);
767:   MatAssemblyEnd(user->M, MAT_FINAL_ASSEMBLY);

769:   /* Generate 2D grid */
770:   VecCreate(PETSC_COMM_WORLD, &XX);
771:   VecCreate(PETSC_COMM_WORLD, &user->q);
772:   VecSetSizes(XX, PETSC_DECIDE, n);
773:   VecSetSizes(user->q, PETSC_DECIDE, n * user->nt);
774:   VecSetFromOptions(XX);
775:   VecSetFromOptions(user->q);

777:   VecDuplicate(XX, &YY);
778:   VecDuplicate(XX, &XXwork);
779:   VecDuplicate(XX, &YYwork);
780:   VecDuplicate(XX, &user->d);
781:   VecDuplicate(XX, &user->dwork);

783:   VecGetOwnershipRange(XX, &istart, &iend);
784:   for (linear_index = istart; linear_index < iend; linear_index++) {
785:     i  = linear_index % user->mx;
786:     j  = (linear_index - i) / user->mx;
787:     vx = h * (i + 0.5);
788:     vy = h * (j + 0.5);
789:     VecSetValues(XX, 1, &linear_index, &vx, INSERT_VALUES);
790:     VecSetValues(YY, 1, &linear_index, &vy, INSERT_VALUES);
791:   }

793:   VecAssemblyBegin(XX);
794:   VecAssemblyEnd(XX);
795:   VecAssemblyBegin(YY);
796:   VecAssemblyEnd(YY);

798:   /* Compute final density function yT
799:      yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
800:      yT = yT / (h^2*sum(yT)) */
801:   VecCopy(XX, XXwork);
802:   VecCopy(YY, YYwork);

804:   VecShift(XXwork, -0.25);
805:   VecShift(YYwork, -0.25);

807:   VecPointwiseMult(XXwork, XXwork, XXwork);
808:   VecPointwiseMult(YYwork, YYwork, YYwork);

810:   VecCopy(XXwork, user->dwork);
811:   VecAXPY(user->dwork, 1.0, YYwork);
812:   VecScale(user->dwork, -30.0);
813:   VecExp(user->dwork);
814:   VecCopy(user->dwork, user->d);

816:   VecCopy(XX, XXwork);
817:   VecCopy(YY, YYwork);

819:   VecShift(XXwork, -0.75);
820:   VecShift(YYwork, -0.75);

822:   VecPointwiseMult(XXwork, XXwork, XXwork);
823:   VecPointwiseMult(YYwork, YYwork, YYwork);

825:   VecCopy(XXwork, user->dwork);
826:   VecAXPY(user->dwork, 1.0, YYwork);
827:   VecScale(user->dwork, -30.0);
828:   VecExp(user->dwork);

830:   VecAXPY(user->d, 1.0, user->dwork);
831:   VecShift(user->d, 1.0);
832:   VecSum(user->d, &sum);
833:   VecScale(user->d, 1.0 / (h * h * sum));

835:   /* Initial conditions of forward problem */
836:   VecDuplicate(XX, &bc);
837:   VecCopy(XX, XXwork);
838:   VecCopy(YY, YYwork);

840:   VecShift(XXwork, -0.5);
841:   VecShift(YYwork, -0.5);

843:   VecPointwiseMult(XXwork, XXwork, XXwork);
844:   VecPointwiseMult(YYwork, YYwork, YYwork);

846:   VecWAXPY(bc, 1.0, XXwork, YYwork);
847:   VecScale(bc, -50.0);
848:   VecExp(bc);
849:   VecShift(bc, 1.0);
850:   VecSum(bc, &sum);
851:   VecScale(bc, 1.0 / (h * h * sum));

853:   /* Create scatter from y to y_1,y_2,...,y_nt */
854:   /*  TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
855:   PetscMalloc1(user->nt * user->mx * user->mx, &user->yi_scatter);
856:   VecCreate(PETSC_COMM_WORLD, &yi);
857:   VecSetSizes(yi, PETSC_DECIDE, user->mx * user->mx);
858:   VecSetFromOptions(yi);
859:   VecDuplicateVecs(yi, user->nt, &user->yi);
860:   VecDuplicateVecs(yi, user->nt, &user->yiwork);
861:   VecDuplicateVecs(yi, user->nt, &user->ziwork);
862:   for (i = 0; i < user->nt; i++) {
863:     VecGetOwnershipRange(user->yi[i], &lo, &hi);
864:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_yi);
865:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + i * user->mx * user->mx, 1, &is_from_y);
866:     VecScatterCreate(user->y, is_from_y, user->yi[i], is_to_yi, &user->yi_scatter[i]);
867:     ISDestroy(&is_to_yi);
868:     ISDestroy(&is_from_y);
869:   }

871:   /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
872:   /*  TODO: reorder for better parallelism */
873:   PetscMalloc1(user->nt * user->mx * user->mx, &user->uxi_scatter);
874:   PetscMalloc1(user->nt * user->mx * user->mx, &user->uyi_scatter);
875:   PetscMalloc1(user->nt * user->mx * user->mx, &user->ux_scatter);
876:   PetscMalloc1(user->nt * user->mx * user->mx, &user->uy_scatter);
877:   PetscMalloc1(2 * user->nt * user->mx * user->mx, &user->ui_scatter);
878:   VecCreate(PETSC_COMM_WORLD, &uxi);
879:   VecCreate(PETSC_COMM_WORLD, &ui);
880:   VecSetSizes(uxi, PETSC_DECIDE, user->mx * user->mx);
881:   VecSetSizes(ui, PETSC_DECIDE, 2 * user->mx * user->mx);
882:   VecSetFromOptions(uxi);
883:   VecSetFromOptions(ui);
884:   VecDuplicateVecs(uxi, user->nt, &user->uxi);
885:   VecDuplicateVecs(uxi, user->nt, &user->uyi);
886:   VecDuplicateVecs(uxi, user->nt, &user->uxiwork);
887:   VecDuplicateVecs(uxi, user->nt, &user->uyiwork);
888:   VecDuplicateVecs(ui, user->nt, &user->ui);
889:   VecDuplicateVecs(ui, user->nt, &user->uiwork);
890:   for (i = 0; i < user->nt; i++) {
891:     VecGetOwnershipRange(user->uxi[i], &lo, &hi);
892:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi);
893:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u);
894:     VecScatterCreate(user->u, is_from_u, user->uxi[i], is_to_uxi, &user->uxi_scatter[i]);

896:     ISDestroy(&is_to_uxi);
897:     ISDestroy(&is_from_u);

899:     VecGetOwnershipRange(user->uyi[i], &lo, &hi);
900:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi);
901:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + (2 * i + 1) * user->mx * user->mx, 1, &is_from_u);
902:     VecScatterCreate(user->u, is_from_u, user->uyi[i], is_to_uyi, &user->uyi_scatter[i]);

904:     ISDestroy(&is_to_uyi);
905:     ISDestroy(&is_from_u);

907:     VecGetOwnershipRange(user->uxi[i], &lo, &hi);
908:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi);
909:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_from_u);
910:     VecScatterCreate(user->ui[i], is_from_u, user->uxi[i], is_to_uxi, &user->ux_scatter[i]);

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

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

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

923:     VecGetOwnershipRange(user->ui[i], &lo, &hi);
924:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi);
925:     ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u);
926:     VecScatterCreate(user->u, is_from_u, user->ui[i], is_to_uxi, &user->ui_scatter[i]);

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

932:   /* RHS of forward problem */
933:   MatMult(user->M, bc, user->yiwork[0]);
934:   for (i = 1; i < user->nt; i++) VecSet(user->yiwork[i], 0.0);
935:   Gather_yi(user->q, user->yiwork, user->yi_scatter, user->nt);

937:   /* Compute true velocity field utrue */
938:   VecDuplicate(user->u, &user->utrue);
939:   for (i = 0; i < user->nt; i++) {
940:     VecCopy(YY, user->uxi[i]);
941:     VecScale(user->uxi[i], 150.0 * i * user->ht);
942:     VecCopy(XX, user->uyi[i]);
943:     VecShift(user->uyi[i], -10.0);
944:     VecScale(user->uyi[i], 15.0 * i * user->ht);
945:   }
946:   Gather_uxi_uyi(user->utrue, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt);

948:   /* Initial guess and reference model */
949:   VecDuplicate(user->utrue, &user->ur);
950:   for (i = 0; i < user->nt; i++) {
951:     VecCopy(XX, user->uxi[i]);
952:     VecShift(user->uxi[i], i * user->ht);
953:     VecCopy(YY, user->uyi[i]);
954:     VecShift(user->uyi[i], -i * user->ht);
955:   }
956:   Gather_uxi_uyi(user->ur, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt);

958:   /* Generate regularization matrix L */
959:   MatCreate(PETSC_COMM_WORLD, &user->LT);
960:   MatSetSizes(user->LT, PETSC_DECIDE, PETSC_DECIDE, 2 * n * user->nt, n * user->nt);
961:   MatSetFromOptions(user->LT);
962:   MatMPIAIJSetPreallocation(user->LT, 1, NULL, 1, NULL);
963:   MatSeqAIJSetPreallocation(user->LT, 1, NULL);
964:   MatGetOwnershipRange(user->LT, &istart, &iend);

966:   for (i = istart; i < iend; i++) {
967:     iblock = (i + n) / (2 * n);
968:     j      = i - iblock * n;
969:     MatSetValues(user->LT, 1, &i, 1, &j, &user->gamma, INSERT_VALUES);
970:   }

972:   MatAssemblyBegin(user->LT, MAT_FINAL_ASSEMBLY);
973:   MatAssemblyEnd(user->LT, MAT_FINAL_ASSEMBLY);

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

977:   /* Build work vectors and matrices */
978:   VecCreate(PETSC_COMM_WORLD, &user->lwork);
979:   VecSetType(user->lwork, VECMPI);
980:   VecSetSizes(user->lwork, PETSC_DECIDE, user->m);
981:   VecSetFromOptions(user->lwork);

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

985:   VecDuplicate(user->y, &user->ywork);
986:   VecDuplicate(user->u, &user->uwork);
987:   VecDuplicate(user->u, &user->vwork);
988:   VecDuplicate(user->u, &user->js_diag);
989:   VecDuplicate(user->c, &user->cwork);

991:   /* Create matrix-free shell user->Js for computing A*x */
992:   MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->Js);
993:   MatShellSetOperation(user->Js, MATOP_MULT, (void (*)(void))StateMatMult);
994:   MatShellSetOperation(user->Js, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate);
995:   MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose);
996:   MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal);

998:   /* Diagonal blocks of user->Js */
999:   MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlock);
1000:   MatShellSetOperation(user->JsBlock, MATOP_MULT, (void (*)(void))StateMatBlockMult);
1001:   MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockMultTranspose);

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

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

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

1020:   /* Build matrices for SOR preconditioner */
1021:   Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt);
1022:   PetscMalloc1(5 * n, &user->C);
1023:   PetscMalloc1(2 * n, &user->Cwork);
1024:   for (i = 0; i < user->nt; i++) {
1025:     MatDuplicate(user->Divxy[0], MAT_COPY_VALUES, &user->C[i]);
1026:     MatDuplicate(user->Divxy[1], MAT_COPY_VALUES, &user->Cwork[i]);

1028:     MatDiagonalScale(user->C[i], NULL, user->uxi[i]);
1029:     MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i]);
1030:     MatAXPY(user->C[i], 1.0, user->Cwork[i], DIFFERENT_NONZERO_PATTERN);
1031:     MatScale(user->C[i], user->ht);
1032:     MatShift(user->C[i], 1.0);
1033:   }

1035:   /* Solver options and tolerances */
1036:   KSPCreate(PETSC_COMM_WORLD, &user->solver);
1037:   KSPSetType(user->solver, KSPGMRES);
1038:   KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec);
1039:   KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500);
1040:   /* KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500); */
1041:   KSPGetPC(user->solver, &user->prec);
1042:   PCSetType(user->prec, PCSHELL);

1044:   PCShellSetApply(user->prec, StateMatBlockPrecMult);
1045:   PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMultTranspose);
1046:   PCShellSetContext(user->prec, user);

1048:   /* Compute true state function yt given ut */
1049:   VecCreate(PETSC_COMM_WORLD, &user->ytrue);
1050:   VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt);
1051:   VecSetFromOptions(user->ytrue);
1052:   user->c_formed = PETSC_TRUE;
1053:   VecCopy(user->utrue, user->u); /*  Set u=utrue temporarily for StateMatInv */
1054:   VecSet(user->ytrue, 0.0);      /*  Initial guess */
1055:   StateMatInvMult(user->Js, user->q, user->ytrue);
1056:   VecCopy(user->ur, user->u); /*  Reset u=ur */

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

1061:   /* Data discretization */
1062:   MatCreate(PETSC_COMM_WORLD, &user->Q);
1063:   MatSetSizes(user->Q, PETSC_DECIDE, PETSC_DECIDE, user->mx * user->mx, user->m);
1064:   MatSetFromOptions(user->Q);
1065:   MatMPIAIJSetPreallocation(user->Q, 0, NULL, 1, NULL);
1066:   MatSeqAIJSetPreallocation(user->Q, 1, NULL);

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

1070:   for (i = istart; i < iend; i++) {
1071:     j = i + user->m - user->mx * user->mx;
1072:     MatSetValues(user->Q, 1, &i, 1, &j, &one, INSERT_VALUES);
1073:   }

1075:   MatAssemblyBegin(user->Q, MAT_FINAL_ASSEMBLY);
1076:   MatAssemblyEnd(user->Q, MAT_FINAL_ASSEMBLY);

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

1080:   VecDestroy(&XX);
1081:   VecDestroy(&YY);
1082:   VecDestroy(&XXwork);
1083:   VecDestroy(&YYwork);
1084:   VecDestroy(&yi);
1085:   VecDestroy(&uxi);
1086:   VecDestroy(&ui);
1087:   VecDestroy(&bc);

1089:   /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1090:   KSPSetFromOptions(user->solver);
1091:   return 0;
1092: }

1094: PetscErrorCode HyperbolicDestroy(AppCtx *user)
1095: {
1096:   PetscInt i;

1098:   MatDestroy(&user->Q);
1099:   MatDestroy(&user->QT);
1100:   MatDestroy(&user->Div);
1101:   MatDestroy(&user->Divwork);
1102:   MatDestroy(&user->Grad);
1103:   MatDestroy(&user->L);
1104:   MatDestroy(&user->LT);
1105:   KSPDestroy(&user->solver);
1106:   MatDestroy(&user->Js);
1107:   MatDestroy(&user->Jd);
1108:   MatDestroy(&user->JsBlockPrec);
1109:   MatDestroy(&user->JsInv);
1110:   MatDestroy(&user->JsBlock);
1111:   MatDestroy(&user->Divxy[0]);
1112:   MatDestroy(&user->Divxy[1]);
1113:   MatDestroy(&user->Gradxy[0]);
1114:   MatDestroy(&user->Gradxy[1]);
1115:   MatDestroy(&user->M);
1116:   for (i = 0; i < user->nt; i++) {
1117:     MatDestroy(&user->C[i]);
1118:     MatDestroy(&user->Cwork[i]);
1119:   }
1120:   PetscFree(user->C);
1121:   PetscFree(user->Cwork);
1122:   VecDestroy(&user->u);
1123:   VecDestroy(&user->uwork);
1124:   VecDestroy(&user->vwork);
1125:   VecDestroy(&user->utrue);
1126:   VecDestroy(&user->y);
1127:   VecDestroy(&user->ywork);
1128:   VecDestroy(&user->ytrue);
1129:   VecDestroyVecs(user->nt, &user->yi);
1130:   VecDestroyVecs(user->nt, &user->yiwork);
1131:   VecDestroyVecs(user->nt, &user->ziwork);
1132:   VecDestroyVecs(user->nt, &user->uxi);
1133:   VecDestroyVecs(user->nt, &user->uyi);
1134:   VecDestroyVecs(user->nt, &user->uxiwork);
1135:   VecDestroyVecs(user->nt, &user->uyiwork);
1136:   VecDestroyVecs(user->nt, &user->ui);
1137:   VecDestroyVecs(user->nt, &user->uiwork);
1138:   VecDestroy(&user->c);
1139:   VecDestroy(&user->cwork);
1140:   VecDestroy(&user->ur);
1141:   VecDestroy(&user->q);
1142:   VecDestroy(&user->d);
1143:   VecDestroy(&user->dwork);
1144:   VecDestroy(&user->lwork);
1145:   VecDestroy(&user->js_diag);
1146:   ISDestroy(&user->s_is);
1147:   ISDestroy(&user->d_is);
1148:   VecScatterDestroy(&user->state_scatter);
1149:   VecScatterDestroy(&user->design_scatter);
1150:   for (i = 0; i < user->nt; i++) {
1151:     VecScatterDestroy(&user->uxi_scatter[i]);
1152:     VecScatterDestroy(&user->uyi_scatter[i]);
1153:     VecScatterDestroy(&user->ux_scatter[i]);
1154:     VecScatterDestroy(&user->uy_scatter[i]);
1155:     VecScatterDestroy(&user->ui_scatter[i]);
1156:     VecScatterDestroy(&user->yi_scatter[i]);
1157:   }
1158:   PetscFree(user->uxi_scatter);
1159:   PetscFree(user->uyi_scatter);
1160:   PetscFree(user->ux_scatter);
1161:   PetscFree(user->uy_scatter);
1162:   PetscFree(user->ui_scatter);
1163:   PetscFree(user->yi_scatter);
1164:   return 0;
1165: }

1167: PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1168: {
1169:   Vec       X;
1170:   PetscReal unorm, ynorm;
1171:   AppCtx   *user = (AppCtx *)ptr;

1173:   TaoGetSolution(tao, &X);
1174:   Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter);
1175:   VecAXPY(user->ywork, -1.0, user->ytrue);
1176:   VecAXPY(user->uwork, -1.0, user->utrue);
1177:   VecNorm(user->uwork, NORM_2, &unorm);
1178:   VecNorm(user->ywork, NORM_2, &ynorm);
1179:   PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n", (double)unorm, (double)ynorm);
1180:   return 0;
1181: }

1183: /*TEST

1185:    build:
1186:       requires: !complex

1188:    test:
1189:       requires: !single
1190:       args: -tao_cmonitor -tao_max_funcs 10 -tao_type lcl -tao_gatol 1.e-5

1192:    test:
1193:       suffix: guess_pod
1194:       requires: !single
1195:       args: -tao_cmonitor -tao_max_funcs 10 -tao_type lcl -ksp_guess_type pod -tao_gatol 1.e-5

1197: TEST*/