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

107:   PetscFunctionBeginUser;
108:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
109:   user.mx = 32;
110:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "hyperbolic example", NULL);
111:   PetscCall(PetscOptionsInt("-mx", "Number of grid points in each direction", "", user.mx, &user.mx, NULL));
112:   user.nt = 16;
113:   PetscCall(PetscOptionsInt("-nt", "Number of time steps", "", user.nt, &user.nt, NULL));
114:   user.ndata = 64;
115:   PetscCall(PetscOptionsInt("-ndata", "Numbers of data points per sample", "", user.ndata, &user.ndata, NULL));
116:   user.alpha = 10.0;
117:   PetscCall(PetscOptionsReal("-alpha", "Regularization parameter", "", user.alpha, &user.alpha, NULL));
118:   user.T = 1.0 / 32.0;
119:   PetscCall(PetscOptionsReal("-Tfinal", "Final time", "", user.T, &user.T, NULL));
120:   PetscCall(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:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user.u));
129:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user.y));
130:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user.c));
131:   PetscCall(VecSetSizes(user.u, PETSC_DECIDE, user.n - user.m));
132:   PetscCall(VecSetSizes(user.y, PETSC_DECIDE, user.m));
133:   PetscCall(VecSetSizes(user.c, PETSC_DECIDE, user.m));
134:   PetscCall(VecSetFromOptions(user.u));
135:   PetscCall(VecSetFromOptions(user.y));
136:   PetscCall(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:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));

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

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

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

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

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

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

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

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

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

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

206:   PetscCall(TaoDestroy(&tao));
207:   PetscCall(VecDestroy(&x));
208:   PetscCall(VecDestroy(&x0));
209:   PetscCall(HyperbolicDestroy(&user));
210:   PetscCall(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:   PetscFunctionBegin;
225:   PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
226:   PetscCall(MatMult(user->Q, user->y, user->dwork));
227:   PetscCall(VecAXPY(user->dwork, -1.0, user->d));
228:   PetscCall(VecDot(user->dwork, user->dwork, &d1));

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

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

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

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

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

259:   PetscCall(VecPointwiseMult(user->vwork, user->vwork, user->vwork));
260:   PetscCall(MatMult(user->L, user->vwork, user->lwork));
261:   PetscCall(VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork));

263:   PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
264:   PetscFunctionReturn(PETSC_SUCCESS);
265: }

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

272:   PetscFunctionBegin;
273:   PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
274:   PetscCall(MatMult(user->Q, user->y, user->dwork));
275:   PetscCall(VecAXPY(user->dwork, -1.0, user->d));

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

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

281:   PetscCall(MatMult(user->LT, user->y, user->uwork));
282:   PetscCall(VecWAXPY(user->vwork, -1.0, user->ur, user->u));
283:   PetscCall(VecPointwiseMult(user->uwork, user->vwork, user->uwork));
284:   PetscCall(VecScale(user->uwork, user->alpha));

286:   PetscCall(VecPointwiseMult(user->vwork, user->vwork, user->vwork));
287:   PetscCall(MatMult(user->L, user->vwork, user->lwork));
288:   PetscCall(VecAXPY(user->ywork, 0.5 * user->alpha, user->lwork));

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

292:   *f = 0.5 * (d1 + user->alpha * d2);
293:   PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

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

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

314:     PetscCall(MatDiagonalScale(user->C[i], NULL, user->uxi[i]));
315:     PetscCall(MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i]));
316:     PetscCall(MatAXPY(user->C[i], 1.0, user->Cwork[i], SUBSET_NONZERO_PATTERN));
317:     PetscCall(MatScale(user->C[i], user->ht));
318:     PetscCall(MatShift(user->C[i], 1.0));
319:   }
320:   PetscFunctionReturn(PETSC_SUCCESS);
321: }

323: /* ------------------------------------------------------------------- */
324: /* B */
325: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
326: {
327:   AppCtx *user = (AppCtx *)ptr;

329:   PetscFunctionBegin;
330:   PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
331:   PetscFunctionReturn(PETSC_SUCCESS);
332: }

334: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
335: {
336:   PetscInt i;
337:   AppCtx  *user;

339:   PetscFunctionBegin;
340:   PetscCall(MatShellGetContext(J_shell, &user));
341:   PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt));
342:   user->block_index = 0;
343:   PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0]));

345:   for (i = 1; i < user->nt; i++) {
346:     user->block_index = i;
347:     PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i]));
348:     PetscCall(MatMult(user->M, user->yi[i - 1], user->ziwork[i - 1]));
349:     PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i - 1]));
350:   }
351:   PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }

355: PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y)
356: {
357:   PetscInt i;
358:   AppCtx  *user;

360:   PetscFunctionBegin;
361:   PetscCall(MatShellGetContext(J_shell, &user));
362:   PetscCall(Scatter_yi(X, user->yi, user->yi_scatter, user->nt));

364:   for (i = 0; i < user->nt - 1; i++) {
365:     user->block_index = i;
366:     PetscCall(MatMultTranspose(user->JsBlock, user->yi[i], user->yiwork[i]));
367:     PetscCall(MatMult(user->M, user->yi[i + 1], user->ziwork[i + 1]));
368:     PetscCall(VecAXPY(user->yiwork[i], -1.0, user->ziwork[i + 1]));
369:   }

371:   i                 = user->nt - 1;
372:   user->block_index = i;
373:   PetscCall(MatMultTranspose(user->JsBlock, user->yi[i], user->yiwork[i]));
374:   PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
375:   PetscFunctionReturn(PETSC_SUCCESS);
376: }

378: PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y)
379: {
380:   PetscInt i;
381:   AppCtx  *user;

383:   PetscFunctionBegin;
384:   PetscCall(MatShellGetContext(J_shell, &user));
385:   i = user->block_index;
386:   PetscCall(VecPointwiseMult(user->uxiwork[i], X, user->uxi[i]));
387:   PetscCall(VecPointwiseMult(user->uyiwork[i], X, user->uyi[i]));
388:   PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
389:   PetscCall(MatMult(user->Div, user->uiwork[i], Y));
390:   PetscCall(VecAYPX(Y, user->ht, X));
391:   PetscFunctionReturn(PETSC_SUCCESS);
392: }

394: PetscErrorCode StateMatBlockMultTranspose(Mat J_shell, Vec X, Vec Y)
395: {
396:   PetscInt i;
397:   AppCtx  *user;

399:   PetscFunctionBegin;
400:   PetscCall(MatShellGetContext(J_shell, &user));
401:   i = user->block_index;
402:   PetscCall(MatMult(user->Grad, X, user->uiwork[i]));
403:   PetscCall(Scatter(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
404:   PetscCall(VecPointwiseMult(user->uxiwork[i], user->uxi[i], user->uxiwork[i]));
405:   PetscCall(VecPointwiseMult(user->uyiwork[i], user->uyi[i], user->uyiwork[i]));
406:   PetscCall(VecWAXPY(Y, 1.0, user->uxiwork[i], user->uyiwork[i]));
407:   PetscCall(VecAYPX(Y, user->ht, X));
408:   PetscFunctionReturn(PETSC_SUCCESS);
409: }

411: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
412: {
413:   PetscInt i;
414:   AppCtx  *user;

416:   PetscFunctionBegin;
417:   PetscCall(MatShellGetContext(J_shell, &user));
418:   PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt));
419:   PetscCall(Scatter_uxi_uyi(X, user->uxiwork, user->uxi_scatter, user->uyiwork, user->uyi_scatter, user->nt));
420:   for (i = 0; i < user->nt; i++) {
421:     PetscCall(VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i]));
422:     PetscCall(VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i]));
423:     PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
424:     PetscCall(MatMult(user->Div, user->uiwork[i], user->ziwork[i]));
425:     PetscCall(VecScale(user->ziwork[i], user->ht));
426:   }
427:   PetscCall(Gather_yi(Y, user->ziwork, user->yi_scatter, user->nt));
428:   PetscFunctionReturn(PETSC_SUCCESS);
429: }

431: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
432: {
433:   PetscInt i;
434:   AppCtx  *user;

436:   PetscFunctionBegin;
437:   PetscCall(MatShellGetContext(J_shell, &user));
438:   PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt));
439:   PetscCall(Scatter_yi(X, user->yiwork, user->yi_scatter, user->nt));
440:   for (i = 0; i < user->nt; i++) {
441:     PetscCall(MatMult(user->Grad, user->yiwork[i], user->uiwork[i]));
442:     PetscCall(Scatter(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
443:     PetscCall(VecPointwiseMult(user->uxiwork[i], user->yi[i], user->uxiwork[i]));
444:     PetscCall(VecPointwiseMult(user->uyiwork[i], user->yi[i], user->uyiwork[i]));
445:     PetscCall(Gather(user->uiwork[i], user->uxiwork[i], user->ux_scatter[i], user->uyiwork[i], user->uy_scatter[i]));
446:     PetscCall(VecScale(user->uiwork[i], user->ht));
447:   }
448:   PetscCall(Gather_yi(Y, user->uiwork, user->ui_scatter, user->nt));
449:   PetscFunctionReturn(PETSC_SUCCESS);
450: }

452: PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y)
453: {
454:   PetscInt i;
455:   AppCtx  *user;

457:   PetscFunctionBegin;
458:   PetscCall(PCShellGetContext(PC_shell, &user));
459:   i = user->block_index;
460:   if (user->c_formed) {
461:     PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y));
462:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed");
463:   PetscFunctionReturn(PETSC_SUCCESS);
464: }

466: PetscErrorCode StateMatBlockPrecMultTranspose(PC PC_shell, Vec X, Vec Y)
467: {
468:   PetscInt i;
469:   AppCtx  *user;

471:   PetscFunctionBegin;
472:   PetscCall(PCShellGetContext(PC_shell, &user));

474:   i = user->block_index;
475:   if (user->c_formed) {
476:     PetscCall(MatSOR(user->C[i], X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP), 0.0, 1, 1, Y));
477:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not formed");
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
482: {
483:   AppCtx  *user;
484:   PetscInt its, i;

486:   PetscFunctionBegin;
487:   PetscCall(MatShellGetContext(J_shell, &user));

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

499:   user->block_index = 0;
500:   PetscCall(KSPSolve(user->solver, user->yi[0], user->yiwork[0]));

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

510:     PetscCall(KSPGetIterationNumber(user->solver, &its));
511:     user->ksp_its = user->ksp_its + its;
512:   }
513:   PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
514:   PetscFunctionReturn(PETSC_SUCCESS);
515: }

517: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
518: {
519:   AppCtx  *user;
520:   PetscInt its, i;

522:   PetscFunctionBegin;
523:   PetscCall(MatShellGetContext(J_shell, &user));

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

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

533:   PetscCall(KSPGetIterationNumber(user->solver, &its));
534:   user->ksp_its = user->ksp_its + its;

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

542:     PetscCall(KSPGetIterationNumber(user->solver, &its));
543:     user->ksp_its = user->ksp_its + its;
544:   }
545:   PetscCall(Gather_yi(Y, user->yiwork, user->yi_scatter, user->nt));
546:   PetscFunctionReturn(PETSC_SUCCESS);
547: }

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

553:   PetscFunctionBegin;
554:   PetscCall(MatShellGetContext(J_shell, &user));

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

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

568:   PetscFunctionBegin;
569:   PetscCall(MatShellGetContext(J_shell, &user));
570:   PetscCall(VecCopy(user->js_diag, X));
571:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
586:   PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter));
587:   PetscCall(Scatter_yi(user->y, user->yi, user->yi_scatter, user->nt));
588:   PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));

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

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

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

603:   PetscFunctionReturn(PETSC_SUCCESS);
604: }

606: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
607: {
608:   PetscFunctionBegin;
609:   PetscCall(VecScatterBegin(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD));
610:   PetscCall(VecScatterEnd(s_scat, x, state, INSERT_VALUES, SCATTER_FORWARD));
611:   PetscCall(VecScatterBegin(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD));
612:   PetscCall(VecScatterEnd(d_scat, x, design, INSERT_VALUES, SCATTER_FORWARD));
613:   PetscFunctionReturn(PETSC_SUCCESS);
614: }

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

620:   PetscFunctionBegin;
621:   for (i = 0; i < nt; i++) {
622:     PetscCall(VecScatterBegin(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD));
623:     PetscCall(VecScatterEnd(scatx[i], u, uxi[i], INSERT_VALUES, SCATTER_FORWARD));
624:     PetscCall(VecScatterBegin(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD));
625:     PetscCall(VecScatterEnd(scaty[i], u, uyi[i], INSERT_VALUES, SCATTER_FORWARD));
626:   }
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
631: {
632:   PetscFunctionBegin;
633:   PetscCall(VecScatterBegin(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE));
634:   PetscCall(VecScatterEnd(s_scat, state, x, INSERT_VALUES, SCATTER_REVERSE));
635:   PetscCall(VecScatterBegin(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE));
636:   PetscCall(VecScatterEnd(d_scat, design, x, INSERT_VALUES, SCATTER_REVERSE));
637:   PetscFunctionReturn(PETSC_SUCCESS);
638: }

640: PetscErrorCode Gather_uxi_uyi(Vec u, Vec *uxi, VecScatter *scatx, Vec *uyi, VecScatter *scaty, PetscInt nt)
641: {
642:   PetscInt i;

644:   PetscFunctionBegin;
645:   for (i = 0; i < nt; i++) {
646:     PetscCall(VecScatterBegin(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE));
647:     PetscCall(VecScatterEnd(scatx[i], uxi[i], u, INSERT_VALUES, SCATTER_REVERSE));
648:     PetscCall(VecScatterBegin(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE));
649:     PetscCall(VecScatterEnd(scaty[i], uyi[i], u, INSERT_VALUES, SCATTER_REVERSE));
650:   }
651:   PetscFunctionReturn(PETSC_SUCCESS);
652: }

654: PetscErrorCode Scatter_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
655: {
656:   PetscInt i;

658:   PetscFunctionBegin;
659:   for (i = 0; i < nt; i++) {
660:     PetscCall(VecScatterBegin(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD));
661:     PetscCall(VecScatterEnd(scat[i], y, yi[i], INSERT_VALUES, SCATTER_FORWARD));
662:   }
663:   PetscFunctionReturn(PETSC_SUCCESS);
664: }

666: PetscErrorCode Gather_yi(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
667: {
668:   PetscInt i;

670:   PetscFunctionBegin;
671:   for (i = 0; i < nt; i++) {
672:     PetscCall(VecScatterBegin(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE));
673:     PetscCall(VecScatterEnd(scat[i], yi[i], y, INSERT_VALUES, SCATTER_REVERSE));
674:   }
675:   PetscFunctionReturn(PETSC_SUCCESS);
676: }

678: PetscErrorCode HyperbolicInitialize(AppCtx *user)
679: {
680:   PetscInt    n, i, j, linear_index, istart, iend, iblock, lo, hi;
681:   Vec         XX, YY, XXwork, YYwork, yi, uxi, ui, bc;
682:   PetscReal   h, sum;
683:   PetscScalar hinv, neg_hinv, quarter = 0.25, one = 1.0, half_hinv, neg_half_hinv;
684:   PetscScalar vx, vy, zero = 0.0;
685:   IS          is_from_y, is_to_yi, is_from_u, is_to_uxi, is_to_uyi;

687:   PetscFunctionBegin;
688:   user->jformed  = PETSC_FALSE;
689:   user->c_formed = PETSC_FALSE;

691:   user->ksp_its         = 0;
692:   user->ksp_its_initial = 0;

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

696:   h             = 1.0 / user->mx;
697:   hinv          = user->mx;
698:   neg_hinv      = -hinv;
699:   half_hinv     = hinv / 2.0;
700:   neg_half_hinv = neg_hinv / 2.0;

702:   /* Generate Grad matrix */
703:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Grad));
704:   PetscCall(MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, 2 * n, n));
705:   PetscCall(MatSetFromOptions(user->Grad));
706:   PetscCall(MatMPIAIJSetPreallocation(user->Grad, 3, NULL, 3, NULL));
707:   PetscCall(MatSeqAIJSetPreallocation(user->Grad, 3, NULL));
708:   PetscCall(MatGetOwnershipRange(user->Grad, &istart, &iend));

710:   for (i = istart; i < iend; i++) {
711:     if (i < n) {
712:       iblock = i / user->mx;
713:       j      = iblock * user->mx + ((i + user->mx - 1) % user->mx);
714:       PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
715:       j = iblock * user->mx + ((i + 1) % user->mx);
716:       PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
717:     }
718:     if (i >= n) {
719:       j = (i - user->mx) % n;
720:       PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
721:       j = (j + 2 * user->mx) % n;
722:       PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
723:     }
724:   }

726:   PetscCall(MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY));
727:   PetscCall(MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY));

729:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[0]));
730:   PetscCall(MatSetSizes(user->Gradxy[0], PETSC_DECIDE, PETSC_DECIDE, n, n));
731:   PetscCall(MatSetFromOptions(user->Gradxy[0]));
732:   PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[0], 3, NULL, 3, NULL));
733:   PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[0], 3, NULL));
734:   PetscCall(MatGetOwnershipRange(user->Gradxy[0], &istart, &iend));

736:   for (i = istart; i < iend; i++) {
737:     iblock = i / user->mx;
738:     j      = iblock * user->mx + ((i + user->mx - 1) % user->mx);
739:     PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
740:     j = iblock * user->mx + ((i + 1) % user->mx);
741:     PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
742:     PetscCall(MatSetValues(user->Gradxy[0], 1, &i, 1, &i, &zero, INSERT_VALUES));
743:   }
744:   PetscCall(MatAssemblyBegin(user->Gradxy[0], MAT_FINAL_ASSEMBLY));
745:   PetscCall(MatAssemblyEnd(user->Gradxy[0], MAT_FINAL_ASSEMBLY));

747:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Gradxy[1]));
748:   PetscCall(MatSetSizes(user->Gradxy[1], PETSC_DECIDE, PETSC_DECIDE, n, n));
749:   PetscCall(MatSetFromOptions(user->Gradxy[1]));
750:   PetscCall(MatMPIAIJSetPreallocation(user->Gradxy[1], 3, NULL, 3, NULL));
751:   PetscCall(MatSeqAIJSetPreallocation(user->Gradxy[1], 3, NULL));
752:   PetscCall(MatGetOwnershipRange(user->Gradxy[1], &istart, &iend));

754:   for (i = istart; i < iend; i++) {
755:     j = (i + n - user->mx) % n;
756:     PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &half_hinv, INSERT_VALUES));
757:     j = (j + 2 * user->mx) % n;
758:     PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &j, &neg_half_hinv, INSERT_VALUES));
759:     PetscCall(MatSetValues(user->Gradxy[1], 1, &i, 1, &i, &zero, INSERT_VALUES));
760:   }
761:   PetscCall(MatAssemblyBegin(user->Gradxy[1], MAT_FINAL_ASSEMBLY));
762:   PetscCall(MatAssemblyEnd(user->Gradxy[1], MAT_FINAL_ASSEMBLY));

764:   /* Generate Div matrix */
765:   PetscCall(MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div));
766:   PetscCall(MatTranspose(user->Gradxy[0], MAT_INITIAL_MATRIX, &user->Divxy[0]));
767:   PetscCall(MatTranspose(user->Gradxy[1], MAT_INITIAL_MATRIX, &user->Divxy[1]));

769:   /* Off-diagonal averaging matrix */
770:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->M));
771:   PetscCall(MatSetSizes(user->M, PETSC_DECIDE, PETSC_DECIDE, n, n));
772:   PetscCall(MatSetFromOptions(user->M));
773:   PetscCall(MatMPIAIJSetPreallocation(user->M, 4, NULL, 4, NULL));
774:   PetscCall(MatSeqAIJSetPreallocation(user->M, 4, NULL));
775:   PetscCall(MatGetOwnershipRange(user->M, &istart, &iend));

777:   for (i = istart; i < iend; i++) {
778:     /* kron(Id,Av) */
779:     iblock = i / user->mx;
780:     j      = iblock * user->mx + ((i + user->mx - 1) % user->mx);
781:     PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
782:     j = iblock * user->mx + ((i + 1) % user->mx);
783:     PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));

785:     /* kron(Av,Id) */
786:     j = (i + user->mx) % n;
787:     PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
788:     j = (i + n - user->mx) % n;
789:     PetscCall(MatSetValues(user->M, 1, &i, 1, &j, &quarter, INSERT_VALUES));
790:   }
791:   PetscCall(MatAssemblyBegin(user->M, MAT_FINAL_ASSEMBLY));
792:   PetscCall(MatAssemblyEnd(user->M, MAT_FINAL_ASSEMBLY));

794:   /* Generate 2D grid */
795:   PetscCall(VecCreate(PETSC_COMM_WORLD, &XX));
796:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->q));
797:   PetscCall(VecSetSizes(XX, PETSC_DECIDE, n));
798:   PetscCall(VecSetSizes(user->q, PETSC_DECIDE, n * user->nt));
799:   PetscCall(VecSetFromOptions(XX));
800:   PetscCall(VecSetFromOptions(user->q));

802:   PetscCall(VecDuplicate(XX, &YY));
803:   PetscCall(VecDuplicate(XX, &XXwork));
804:   PetscCall(VecDuplicate(XX, &YYwork));
805:   PetscCall(VecDuplicate(XX, &user->d));
806:   PetscCall(VecDuplicate(XX, &user->dwork));

808:   PetscCall(VecGetOwnershipRange(XX, &istart, &iend));
809:   for (linear_index = istart; linear_index < iend; linear_index++) {
810:     i  = linear_index % user->mx;
811:     j  = (linear_index - i) / user->mx;
812:     vx = h * (i + 0.5);
813:     vy = h * (j + 0.5);
814:     PetscCall(VecSetValues(XX, 1, &linear_index, &vx, INSERT_VALUES));
815:     PetscCall(VecSetValues(YY, 1, &linear_index, &vy, INSERT_VALUES));
816:   }

818:   PetscCall(VecAssemblyBegin(XX));
819:   PetscCall(VecAssemblyEnd(XX));
820:   PetscCall(VecAssemblyBegin(YY));
821:   PetscCall(VecAssemblyEnd(YY));

823:   /* Compute final density function yT
824:      yT = 1.0 + exp(-30*((x-0.25)^2+(y-0.25)^2)) + exp(-30*((x-0.75)^2+(y-0.75)^2))
825:      yT = yT / (h^2*sum(yT)) */
826:   PetscCall(VecCopy(XX, XXwork));
827:   PetscCall(VecCopy(YY, YYwork));

829:   PetscCall(VecShift(XXwork, -0.25));
830:   PetscCall(VecShift(YYwork, -0.25));

832:   PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
833:   PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));

835:   PetscCall(VecCopy(XXwork, user->dwork));
836:   PetscCall(VecAXPY(user->dwork, 1.0, YYwork));
837:   PetscCall(VecScale(user->dwork, -30.0));
838:   PetscCall(VecExp(user->dwork));
839:   PetscCall(VecCopy(user->dwork, user->d));

841:   PetscCall(VecCopy(XX, XXwork));
842:   PetscCall(VecCopy(YY, YYwork));

844:   PetscCall(VecShift(XXwork, -0.75));
845:   PetscCall(VecShift(YYwork, -0.75));

847:   PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
848:   PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));

850:   PetscCall(VecCopy(XXwork, user->dwork));
851:   PetscCall(VecAXPY(user->dwork, 1.0, YYwork));
852:   PetscCall(VecScale(user->dwork, -30.0));
853:   PetscCall(VecExp(user->dwork));

855:   PetscCall(VecAXPY(user->d, 1.0, user->dwork));
856:   PetscCall(VecShift(user->d, 1.0));
857:   PetscCall(VecSum(user->d, &sum));
858:   PetscCall(VecScale(user->d, 1.0 / (h * h * sum)));

860:   /* Initial conditions of forward problem */
861:   PetscCall(VecDuplicate(XX, &bc));
862:   PetscCall(VecCopy(XX, XXwork));
863:   PetscCall(VecCopy(YY, YYwork));

865:   PetscCall(VecShift(XXwork, -0.5));
866:   PetscCall(VecShift(YYwork, -0.5));

868:   PetscCall(VecPointwiseMult(XXwork, XXwork, XXwork));
869:   PetscCall(VecPointwiseMult(YYwork, YYwork, YYwork));

871:   PetscCall(VecWAXPY(bc, 1.0, XXwork, YYwork));
872:   PetscCall(VecScale(bc, -50.0));
873:   PetscCall(VecExp(bc));
874:   PetscCall(VecShift(bc, 1.0));
875:   PetscCall(VecSum(bc, &sum));
876:   PetscCall(VecScale(bc, 1.0 / (h * h * sum)));

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

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

921:     PetscCall(ISDestroy(&is_to_uxi));
922:     PetscCall(ISDestroy(&is_from_u));

924:     PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi));
925:     PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi));
926:     PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + (2 * i + 1) * user->mx * user->mx, 1, &is_from_u));
927:     PetscCall(VecScatterCreate(user->u, is_from_u, user->uyi[i], is_to_uyi, &user->uyi_scatter[i]));

929:     PetscCall(ISDestroy(&is_to_uyi));
930:     PetscCall(ISDestroy(&is_from_u));

932:     PetscCall(VecGetOwnershipRange(user->uxi[i], &lo, &hi));
933:     PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi));
934:     PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_from_u));
935:     PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uxi[i], is_to_uxi, &user->ux_scatter[i]));

937:     PetscCall(ISDestroy(&is_to_uxi));
938:     PetscCall(ISDestroy(&is_from_u));

940:     PetscCall(VecGetOwnershipRange(user->uyi[i], &lo, &hi));
941:     PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uyi));
942:     PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + user->mx * user->mx, 1, &is_from_u));
943:     PetscCall(VecScatterCreate(user->ui[i], is_from_u, user->uyi[i], is_to_uyi, &user->uy_scatter[i]));

945:     PetscCall(ISDestroy(&is_to_uyi));
946:     PetscCall(ISDestroy(&is_from_u));

948:     PetscCall(VecGetOwnershipRange(user->ui[i], &lo, &hi));
949:     PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo, 1, &is_to_uxi));
950:     PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + 2 * i * user->mx * user->mx, 1, &is_from_u));
951:     PetscCall(VecScatterCreate(user->u, is_from_u, user->ui[i], is_to_uxi, &user->ui_scatter[i]));

953:     PetscCall(ISDestroy(&is_to_uxi));
954:     PetscCall(ISDestroy(&is_from_u));
955:   }

957:   /* RHS of forward problem */
958:   PetscCall(MatMult(user->M, bc, user->yiwork[0]));
959:   for (i = 1; i < user->nt; i++) PetscCall(VecSet(user->yiwork[i], 0.0));
960:   PetscCall(Gather_yi(user->q, user->yiwork, user->yi_scatter, user->nt));

962:   /* Compute true velocity field utrue */
963:   PetscCall(VecDuplicate(user->u, &user->utrue));
964:   for (i = 0; i < user->nt; i++) {
965:     PetscCall(VecCopy(YY, user->uxi[i]));
966:     PetscCall(VecScale(user->uxi[i], 150.0 * i * user->ht));
967:     PetscCall(VecCopy(XX, user->uyi[i]));
968:     PetscCall(VecShift(user->uyi[i], -10.0));
969:     PetscCall(VecScale(user->uyi[i], 15.0 * i * user->ht));
970:   }
971:   PetscCall(Gather_uxi_uyi(user->utrue, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));

973:   /* Initial guess and reference model */
974:   PetscCall(VecDuplicate(user->utrue, &user->ur));
975:   for (i = 0; i < user->nt; i++) {
976:     PetscCall(VecCopy(XX, user->uxi[i]));
977:     PetscCall(VecShift(user->uxi[i], i * user->ht));
978:     PetscCall(VecCopy(YY, user->uyi[i]));
979:     PetscCall(VecShift(user->uyi[i], -i * user->ht));
980:   }
981:   PetscCall(Gather_uxi_uyi(user->ur, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));

983:   /* Generate regularization matrix L */
984:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->LT));
985:   PetscCall(MatSetSizes(user->LT, PETSC_DECIDE, PETSC_DECIDE, 2 * n * user->nt, n * user->nt));
986:   PetscCall(MatSetFromOptions(user->LT));
987:   PetscCall(MatMPIAIJSetPreallocation(user->LT, 1, NULL, 1, NULL));
988:   PetscCall(MatSeqAIJSetPreallocation(user->LT, 1, NULL));
989:   PetscCall(MatGetOwnershipRange(user->LT, &istart, &iend));

991:   for (i = istart; i < iend; i++) {
992:     iblock = (i + n) / (2 * n);
993:     j      = i - iblock * n;
994:     PetscCall(MatSetValues(user->LT, 1, &i, 1, &j, &user->gamma, INSERT_VALUES));
995:   }

997:   PetscCall(MatAssemblyBegin(user->LT, MAT_FINAL_ASSEMBLY));
998:   PetscCall(MatAssemblyEnd(user->LT, MAT_FINAL_ASSEMBLY));

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

1002:   /* Build work vectors and matrices */
1003:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->lwork));
1004:   PetscCall(VecSetType(user->lwork, VECMPI));
1005:   PetscCall(VecSetSizes(user->lwork, PETSC_DECIDE, user->m));
1006:   PetscCall(VecSetFromOptions(user->lwork));

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

1010:   PetscCall(VecDuplicate(user->y, &user->ywork));
1011:   PetscCall(VecDuplicate(user->u, &user->uwork));
1012:   PetscCall(VecDuplicate(user->u, &user->vwork));
1013:   PetscCall(VecDuplicate(user->u, &user->js_diag));
1014:   PetscCall(VecDuplicate(user->c, &user->cwork));

1016:   /* Create matrix-free shell user->Js for computing A*x */
1017:   PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->Js));
1018:   PetscCall(MatShellSetOperation(user->Js, MATOP_MULT, (void (*)(void))StateMatMult));
1019:   PetscCall(MatShellSetOperation(user->Js, MATOP_DUPLICATE, (void (*)(void))StateMatDuplicate));
1020:   PetscCall(MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatMultTranspose));
1021:   PetscCall(MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (void (*)(void))StateMatGetDiagonal));

1023:   /* Diagonal blocks of user->Js */
1024:   PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlock));
1025:   PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT, (void (*)(void))StateMatBlockMult));
1026:   PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockMultTranspose));

1028:   /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
1029:      D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
1030:      This is an SOR preconditioner for user->JsBlock. */
1031:   PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, user, &user->JsBlockPrec));
1032:   PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT, (void (*)(void))StateMatBlockPrecMult));
1033:   PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatBlockPrecMultTranspose));

1035:   /* Create a matrix-free shell user->Jd for computing B*x */
1036:   PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->n - user->m, user, &user->Jd));
1037:   PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT, (void (*)(void))DesignMatMult));
1038:   PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT_TRANSPOSE, (void (*)(void))DesignMatMultTranspose));

1040:   /* User-defined routines for computing user->Js\x and user->Js^T\x*/
1041:   PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsInv));
1042:   PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT, (void (*)(void))StateMatInvMult));
1043:   PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT_TRANSPOSE, (void (*)(void))StateMatInvTransposeMult));

1045:   /* Build matrices for SOR preconditioner */
1046:   PetscCall(Scatter_uxi_uyi(user->u, user->uxi, user->uxi_scatter, user->uyi, user->uyi_scatter, user->nt));
1047:   PetscCall(PetscMalloc1(5 * n, &user->C));
1048:   PetscCall(PetscMalloc1(2 * n, &user->Cwork));
1049:   for (i = 0; i < user->nt; i++) {
1050:     PetscCall(MatDuplicate(user->Divxy[0], MAT_COPY_VALUES, &user->C[i]));
1051:     PetscCall(MatDuplicate(user->Divxy[1], MAT_COPY_VALUES, &user->Cwork[i]));

1053:     PetscCall(MatDiagonalScale(user->C[i], NULL, user->uxi[i]));
1054:     PetscCall(MatDiagonalScale(user->Cwork[i], NULL, user->uyi[i]));
1055:     PetscCall(MatAXPY(user->C[i], 1.0, user->Cwork[i], DIFFERENT_NONZERO_PATTERN));
1056:     PetscCall(MatScale(user->C[i], user->ht));
1057:     PetscCall(MatShift(user->C[i], 1.0));
1058:   }

1060:   /* Solver options and tolerances */
1061:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->solver));
1062:   PetscCall(KSPSetType(user->solver, KSPGMRES));
1063:   PetscCall(KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec));
1064:   PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500));
1065:   /* PetscCall(KSPSetTolerances(user->solver,1e-8,1e-16,1e3,500)); */
1066:   PetscCall(KSPGetPC(user->solver, &user->prec));
1067:   PetscCall(PCSetType(user->prec, PCSHELL));

1069:   PetscCall(PCShellSetApply(user->prec, StateMatBlockPrecMult));
1070:   PetscCall(PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMultTranspose));
1071:   PetscCall(PCShellSetContext(user->prec, user));

1073:   /* Compute true state function yt given ut */
1074:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ytrue));
1075:   PetscCall(VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt));
1076:   PetscCall(VecSetFromOptions(user->ytrue));
1077:   user->c_formed = PETSC_TRUE;
1078:   PetscCall(VecCopy(user->utrue, user->u)); /*  Set u=utrue temporarily for StateMatInv */
1079:   PetscCall(VecSet(user->ytrue, 0.0));      /*  Initial guess */
1080:   PetscCall(StateMatInvMult(user->Js, user->q, user->ytrue));
1081:   PetscCall(VecCopy(user->ur, user->u)); /*  Reset u=ur */

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

1086:   /* Data discretization */
1087:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Q));
1088:   PetscCall(MatSetSizes(user->Q, PETSC_DECIDE, PETSC_DECIDE, user->mx * user->mx, user->m));
1089:   PetscCall(MatSetFromOptions(user->Q));
1090:   PetscCall(MatMPIAIJSetPreallocation(user->Q, 0, NULL, 1, NULL));
1091:   PetscCall(MatSeqAIJSetPreallocation(user->Q, 1, NULL));

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

1095:   for (i = istart; i < iend; i++) {
1096:     j = i + user->m - user->mx * user->mx;
1097:     PetscCall(MatSetValues(user->Q, 1, &i, 1, &j, &one, INSERT_VALUES));
1098:   }

1100:   PetscCall(MatAssemblyBegin(user->Q, MAT_FINAL_ASSEMBLY));
1101:   PetscCall(MatAssemblyEnd(user->Q, MAT_FINAL_ASSEMBLY));

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

1105:   PetscCall(VecDestroy(&XX));
1106:   PetscCall(VecDestroy(&YY));
1107:   PetscCall(VecDestroy(&XXwork));
1108:   PetscCall(VecDestroy(&YYwork));
1109:   PetscCall(VecDestroy(&yi));
1110:   PetscCall(VecDestroy(&uxi));
1111:   PetscCall(VecDestroy(&ui));
1112:   PetscCall(VecDestroy(&bc));

1114:   /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1115:   PetscCall(KSPSetFromOptions(user->solver));
1116:   PetscFunctionReturn(PETSC_SUCCESS);
1117: }

1119: PetscErrorCode HyperbolicDestroy(AppCtx *user)
1120: {
1121:   PetscInt i;

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

1193: PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
1194: {
1195:   Vec       X;
1196:   PetscReal unorm, ynorm;
1197:   AppCtx   *user = (AppCtx *)ptr;

1199:   PetscFunctionBegin;
1200:   PetscCall(TaoGetSolution(tao, &X));
1201:   PetscCall(Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter));
1202:   PetscCall(VecAXPY(user->ywork, -1.0, user->ytrue));
1203:   PetscCall(VecAXPY(user->uwork, -1.0, user->utrue));
1204:   PetscCall(VecNorm(user->uwork, NORM_2, &unorm));
1205:   PetscCall(VecNorm(user->ywork, NORM_2, &ynorm));
1206:   PetscCall(PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n", (double)unorm, (double)ynorm));
1207:   PetscFunctionReturn(PETSC_SUCCESS);
1208: }

1210: /*TEST

1212:    build:
1213:       requires: !complex

1215:    test:
1216:       requires: !single
1217:       args: -tao_cmonitor -tao_max_funcs 10 -tao_type lcl -tao_gatol 1.e-5

1219:    test:
1220:       suffix: guess_pod
1221:       requires: !single
1222:       args: -tao_cmonitor -tao_max_funcs 10 -tao_type lcl -ksp_guess_type pod -tao_gatol 1.e-5

1224: TEST*/