Actual source code: burgers_spectral.c


  2: static char help[] = "Solves a simple data assimilation problem with one dimensional Burger's equation using TSAdjoint\n\n";

  4: /*

  6:     Not yet tested in parallel

  8: */

 10: /* ------------------------------------------------------------------------

 12:    This program uses the one-dimensional Burger's equation
 13:        u_t = mu*u_xx - u u_x,
 14:    on the domain 0 <= x <= 1, with periodic boundary conditions

 16:    to demonstrate solving a data assimilation problem of finding the initial conditions
 17:    to produce a given solution at a fixed time.

 19:    The operators are discretized with the spectral element method

 21:    See the paper PDE-CONSTRAINED OPTIMIZATION WITH SPECTRAL ELEMENTS USING PETSC AND TAO
 22:    by OANA MARIN, EMIL CONSTANTINESCU, AND BARRY SMITH for details on the exact solution
 23:    used

 25:   ------------------------------------------------------------------------- */

 27: #include <petsctao.h>
 28: #include <petscts.h>
 29: #include <petscdt.h>
 30: #include <petscdraw.h>
 31: #include <petscdmda.h>

 33: /*
 34:    User-defined application context - contains data needed by the
 35:    application-provided call-back routines.
 36: */

 38: typedef struct {
 39:   PetscInt   n;       /* number of nodes */
 40:   PetscReal *nodes;   /* GLL nodes */
 41:   PetscReal *weights; /* GLL weights */
 42: } PetscGLL;

 44: typedef struct {
 45:   PetscInt  N;               /* grid points per elements*/
 46:   PetscInt  E;               /* number of elements */
 47:   PetscReal tol_L2, tol_max; /* error norms */
 48:   PetscInt  steps;           /* number of timesteps */
 49:   PetscReal Tend;            /* endtime */
 50:   PetscReal mu;              /* viscosity */
 51:   PetscReal L;               /* total length of domain */
 52:   PetscReal Le;
 53:   PetscReal Tadj;
 54: } PetscParam;

 56: typedef struct {
 57:   Vec obj;  /* desired end state */
 58:   Vec grid; /* total grid */
 59:   Vec grad;
 60:   Vec ic;
 61:   Vec curr_sol;
 62:   Vec true_solution; /* actual initial conditions for the final solution */
 63: } PetscData;

 65: typedef struct {
 66:   Vec      grid;  /* total grid */
 67:   Vec      mass;  /* mass matrix for total integration */
 68:   Mat      stiff; /* stifness matrix */
 69:   Mat      keptstiff;
 70:   Mat      grad;
 71:   PetscGLL gll;
 72: } PetscSEMOperators;

 74: typedef struct {
 75:   DM                da; /* distributed array data structure */
 76:   PetscSEMOperators SEMop;
 77:   PetscParam        param;
 78:   PetscData         dat;
 79:   TS                ts;
 80:   PetscReal         initial_dt;
 81: } AppCtx;

 83: /*
 84:    User-defined routines
 85: */
 86: extern PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 87: extern PetscErrorCode RHSMatrixLaplaciangllDM(TS, PetscReal, Vec, Mat, Mat, void *);
 88: extern PetscErrorCode RHSMatrixAdvectiongllDM(TS, PetscReal, Vec, Mat, Mat, void *);
 89: extern PetscErrorCode InitialConditions(Vec, AppCtx *);
 90: extern PetscErrorCode TrueSolution(Vec, AppCtx *);
 91: extern PetscErrorCode ComputeObjective(PetscReal, Vec, AppCtx *);
 92: extern PetscErrorCode MonitorError(Tao, void *);
 93: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
 94: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);

 96: int main(int argc, char **argv)
 97: {
 98:   AppCtx       appctx; /* user-defined application context */
 99:   Tao          tao;
100:   Vec          u; /* approximate solution vector */
101:   PetscInt     i, xs, xm, ind, j, lenglob;
102:   PetscReal    x, *wrk_ptr1, *wrk_ptr2;
103:   MatNullSpace nsp;
104:   PetscMPIInt  size;

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:      Initialize program and set problem parameters
108:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

113:   /*initialize parameters */
114:   appctx.param.N     = 10;   /* order of the spectral element */
115:   appctx.param.E     = 10;   /* number of elements */
116:   appctx.param.L     = 4.0;  /* length of the domain */
117:   appctx.param.mu    = 0.01; /* diffusion coefficient */
118:   appctx.initial_dt  = 5e-3;
119:   appctx.param.steps = PETSC_MAX_INT;
120:   appctx.param.Tend  = 4;

122:   PetscOptionsGetInt(NULL, NULL, "-N", &appctx.param.N, NULL);
123:   PetscOptionsGetInt(NULL, NULL, "-E", &appctx.param.E, NULL);
124:   PetscOptionsGetReal(NULL, NULL, "-Tend", &appctx.param.Tend, NULL);
125:   PetscOptionsGetReal(NULL, NULL, "-mu", &appctx.param.mu, NULL);
126:   appctx.param.Le = appctx.param.L / appctx.param.E;

128:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

131:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132:      Create GLL data structures
133:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134:   PetscMalloc2(appctx.param.N, &appctx.SEMop.gll.nodes, appctx.param.N, &appctx.SEMop.gll.weights);
135:   PetscDTGaussLobattoLegendreQuadrature(appctx.param.N, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, appctx.SEMop.gll.nodes, appctx.SEMop.gll.weights);
136:   appctx.SEMop.gll.n = appctx.param.N;
137:   lenglob            = appctx.param.E * (appctx.param.N - 1);

139:   /*
140:      Create distributed array (DMDA) to manage parallel grid and vectors
141:      and to set up the ghost point communication pattern.  There are E*(Nl-1)+1
142:      total grid values spread equally among all the processors, except first and last
143:   */

145:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, lenglob, 1, 1, NULL, &appctx.da);
146:   DMSetFromOptions(appctx.da);
147:   DMSetUp(appctx.da);

149:   /*
150:      Extract global and local vectors from DMDA; we use these to store the
151:      approximate solution.  Then duplicate these for remaining vectors that
152:      have the same types.
153:   */

155:   DMCreateGlobalVector(appctx.da, &u);
156:   VecDuplicate(u, &appctx.dat.ic);
157:   VecDuplicate(u, &appctx.dat.true_solution);
158:   VecDuplicate(u, &appctx.dat.obj);
159:   VecDuplicate(u, &appctx.SEMop.grid);
160:   VecDuplicate(u, &appctx.SEMop.mass);
161:   VecDuplicate(u, &appctx.dat.curr_sol);

163:   DMDAGetCorners(appctx.da, &xs, NULL, NULL, &xm, NULL, NULL);
164:   DMDAVecGetArray(appctx.da, appctx.SEMop.grid, &wrk_ptr1);
165:   DMDAVecGetArray(appctx.da, appctx.SEMop.mass, &wrk_ptr2);

167:   /* Compute function over the locally owned part of the grid */

169:   xs = xs / (appctx.param.N - 1);
170:   xm = xm / (appctx.param.N - 1);

172:   /*
173:      Build total grid and mass over entire mesh (multi-elemental)
174:   */

176:   for (i = xs; i < xs + xm; i++) {
177:     for (j = 0; j < appctx.param.N - 1; j++) {
178:       x             = (appctx.param.Le / 2.0) * (appctx.SEMop.gll.nodes[j] + 1.0) + appctx.param.Le * i;
179:       ind           = i * (appctx.param.N - 1) + j;
180:       wrk_ptr1[ind] = x;
181:       wrk_ptr2[ind] = .5 * appctx.param.Le * appctx.SEMop.gll.weights[j];
182:       if (j == 0) wrk_ptr2[ind] += .5 * appctx.param.Le * appctx.SEMop.gll.weights[j];
183:     }
184:   }
185:   DMDAVecRestoreArray(appctx.da, appctx.SEMop.grid, &wrk_ptr1);
186:   DMDAVecRestoreArray(appctx.da, appctx.SEMop.mass, &wrk_ptr2);

188:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189:    Create matrix data structure; set matrix evaluation routine.
190:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191:   DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE);
192:   DMCreateMatrix(appctx.da, &appctx.SEMop.stiff);
193:   DMCreateMatrix(appctx.da, &appctx.SEMop.grad);
194:   /*
195:    For linear problems with a time-dependent f(u,t) in the equation
196:    u_t = f(u,t), the user provides the discretized right-hand-side
197:    as a time-dependent matrix.
198:    */
199:   RHSMatrixLaplaciangllDM(appctx.ts, 0.0, u, appctx.SEMop.stiff, appctx.SEMop.stiff, &appctx);
200:   RHSMatrixAdvectiongllDM(appctx.ts, 0.0, u, appctx.SEMop.grad, appctx.SEMop.grad, &appctx);
201:   /*
202:        For linear problems with a time-dependent f(u,t) in the equation
203:        u_t = f(u,t), the user provides the discretized right-hand-side
204:        as a time-dependent matrix.
205:     */

207:   MatDuplicate(appctx.SEMop.stiff, MAT_COPY_VALUES, &appctx.SEMop.keptstiff);

209:   /* attach the null space to the matrix, this probably is not needed but does no harm */
210:   MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nsp);
211:   MatSetNullSpace(appctx.SEMop.stiff, nsp);
212:   MatSetNullSpace(appctx.SEMop.keptstiff, nsp);
213:   MatNullSpaceTest(nsp, appctx.SEMop.stiff, NULL);
214:   MatNullSpaceDestroy(&nsp);
215:   /* attach the null space to the matrix, this probably is not needed but does no harm */
216:   MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nsp);
217:   MatSetNullSpace(appctx.SEMop.grad, nsp);
218:   MatNullSpaceTest(nsp, appctx.SEMop.grad, NULL);
219:   MatNullSpaceDestroy(&nsp);

221:   /* Create the TS solver that solves the ODE and its adjoint; set its options */
222:   TSCreate(PETSC_COMM_WORLD, &appctx.ts);
223:   TSSetProblemType(appctx.ts, TS_NONLINEAR);
224:   TSSetType(appctx.ts, TSRK);
225:   TSSetDM(appctx.ts, appctx.da);
226:   TSSetTime(appctx.ts, 0.0);
227:   TSSetTimeStep(appctx.ts, appctx.initial_dt);
228:   TSSetMaxSteps(appctx.ts, appctx.param.steps);
229:   TSSetMaxTime(appctx.ts, appctx.param.Tend);
230:   TSSetExactFinalTime(appctx.ts, TS_EXACTFINALTIME_MATCHSTEP);
231:   TSSetTolerances(appctx.ts, 1e-7, NULL, 1e-7, NULL);
232:   TSSetFromOptions(appctx.ts);
233:   /* Need to save initial timestep user may have set with -ts_dt so it can be reset for each new TSSolve() */
234:   TSGetTimeStep(appctx.ts, &appctx.initial_dt);
235:   TSSetRHSFunction(appctx.ts, NULL, RHSFunction, &appctx);
236:   TSSetRHSJacobian(appctx.ts, appctx.SEMop.stiff, appctx.SEMop.stiff, RHSJacobian, &appctx);

238:   /* Set Objective and Initial conditions for the problem and compute Objective function (evolution of true_solution to final time */
239:   InitialConditions(appctx.dat.ic, &appctx);
240:   TrueSolution(appctx.dat.true_solution, &appctx);
241:   ComputeObjective(appctx.param.Tend, appctx.dat.obj, &appctx);

243:   TSSetSaveTrajectory(appctx.ts);
244:   TSSetFromOptions(appctx.ts);

246:   /* Create TAO solver and set desired solution method  */
247:   TaoCreate(PETSC_COMM_WORLD, &tao);
248:   TaoSetMonitor(tao, MonitorError, &appctx, NULL);
249:   TaoSetType(tao, TAOBQNLS);
250:   TaoSetSolution(tao, appctx.dat.ic);
251:   /* Set routine for function and gradient evaluation  */
252:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&appctx);
253:   /* Check for any TAO command line options  */
254:   TaoSetTolerances(tao, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
255:   TaoSetFromOptions(tao);
256:   TaoSolve(tao);

258:   TaoDestroy(&tao);
259:   MatDestroy(&appctx.SEMop.stiff);
260:   MatDestroy(&appctx.SEMop.keptstiff);
261:   MatDestroy(&appctx.SEMop.grad);
262:   VecDestroy(&u);
263:   VecDestroy(&appctx.dat.ic);
264:   VecDestroy(&appctx.dat.true_solution);
265:   VecDestroy(&appctx.dat.obj);
266:   VecDestroy(&appctx.SEMop.grid);
267:   VecDestroy(&appctx.SEMop.mass);
268:   VecDestroy(&appctx.dat.curr_sol);
269:   PetscFree2(appctx.SEMop.gll.nodes, appctx.SEMop.gll.weights);
270:   DMDestroy(&appctx.da);
271:   TSDestroy(&appctx.ts);

273:   /*
274:      Always call PetscFinalize() before exiting a program.  This routine
275:        - finalizes the PETSc libraries as well as MPI
276:        - provides summary and diagnostic information if certain runtime
277:          options are chosen (e.g., -log_summary).
278:   */
279:   PetscFinalize();
280:   return 0;
281: }

283: /* --------------------------------------------------------------------- */
284: /*
285:    InitialConditions - Computes the initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve()

287:                        The routine TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function

289:    Input Parameter:
290:    u - uninitialized solution vector (global)
291:    appctx - user-defined application context

293:    Output Parameter:
294:    u - vector with solution at initial time (global)
295: */
296: PetscErrorCode InitialConditions(Vec u, AppCtx *appctx)
297: {
298:   PetscScalar       *s;
299:   const PetscScalar *xg;
300:   PetscInt           i, xs, xn;

302:   DMDAVecGetArray(appctx->da, u, &s);
303:   DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg);
304:   DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL);
305:   for (i = xs; i < xs + xn; i++) s[i] = 2.0 * appctx->param.mu * PETSC_PI * PetscSinScalar(PETSC_PI * xg[i]) / (2.0 + PetscCosScalar(PETSC_PI * xg[i])) + 0.25 * PetscExpReal(-4.0 * PetscPowReal(xg[i] - 2.0, 2.0));
306:   DMDAVecRestoreArray(appctx->da, u, &s);
307:   DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg);
308:   return 0;
309: }

311: /*
312:    TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function.

314:              InitialConditions() computes the initial conditions for the beginning of the Tao iterations

316:    Input Parameter:
317:    u - uninitialized solution vector (global)
318:    appctx - user-defined application context

320:    Output Parameter:
321:    u - vector with solution at initial time (global)
322: */
323: PetscErrorCode TrueSolution(Vec u, AppCtx *appctx)
324: {
325:   PetscScalar       *s;
326:   const PetscScalar *xg;
327:   PetscInt           i, xs, xn;

329:   DMDAVecGetArray(appctx->da, u, &s);
330:   DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg);
331:   DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL);
332:   for (i = xs; i < xs + xn; i++) s[i] = 2.0 * appctx->param.mu * PETSC_PI * PetscSinScalar(PETSC_PI * xg[i]) / (2.0 + PetscCosScalar(PETSC_PI * xg[i]));
333:   DMDAVecRestoreArray(appctx->da, u, &s);
334:   DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg);
335:   return 0;
336: }
337: /* --------------------------------------------------------------------- */
338: /*
339:    Sets the desired profile for the final end time

341:    Input Parameters:
342:    t - final time
343:    obj - vector storing the desired profile
344:    appctx - user-defined application context

346: */
347: PetscErrorCode ComputeObjective(PetscReal t, Vec obj, AppCtx *appctx)
348: {
349:   PetscScalar       *s;
350:   const PetscScalar *xg;
351:   PetscInt           i, xs, xn;

353:   DMDAVecGetArray(appctx->da, obj, &s);
354:   DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg);
355:   DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL);
356:   for (i = xs; i < xs + xn; i++) {
357:     s[i] = 2.0 * appctx->param.mu * PETSC_PI * PetscSinScalar(PETSC_PI * xg[i]) * PetscExpScalar(-PETSC_PI * PETSC_PI * t * appctx->param.mu) / (2.0 + PetscExpScalar(-PETSC_PI * PETSC_PI * t * appctx->param.mu) * PetscCosScalar(PETSC_PI * xg[i]));
358:   }
359:   DMDAVecRestoreArray(appctx->da, obj, &s);
360:   DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg);
361:   return 0;
362: }

364: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
365: {
366:   AppCtx *appctx = (AppCtx *)ctx;

368:   MatMult(appctx->SEMop.grad, globalin, globalout); /* grad u */
369:   VecPointwiseMult(globalout, globalin, globalout); /* u grad u */
370:   VecScale(globalout, -1.0);
371:   MatMultAdd(appctx->SEMop.keptstiff, globalin, globalout, globalout);
372:   return 0;
373: }

375: /*

377:       K is the discretiziation of the Laplacian
378:       G is the discretization of the gradient

380:       Computes Jacobian of      K u + diag(u) G u   which is given by
381:               K   + diag(u)G + diag(Gu)
382: */
383: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec globalin, Mat A, Mat B, void *ctx)
384: {
385:   AppCtx *appctx = (AppCtx *)ctx;
386:   Vec     Gglobalin;

388:   /*    A = diag(u) G */

390:   MatCopy(appctx->SEMop.grad, A, SAME_NONZERO_PATTERN);
391:   MatDiagonalScale(A, globalin, NULL);

393:   /*    A  = A + diag(Gu) */
394:   VecDuplicate(globalin, &Gglobalin);
395:   MatMult(appctx->SEMop.grad, globalin, Gglobalin);
396:   MatDiagonalSet(A, Gglobalin, ADD_VALUES);
397:   VecDestroy(&Gglobalin);

399:   /*   A  = K - A    */
400:   MatScale(A, -1.0);
401:   MatAXPY(A, 1.0, appctx->SEMop.keptstiff, SAME_NONZERO_PATTERN);
402:   return 0;
403: }

405: /* --------------------------------------------------------------------- */

407: /*
408:    RHSMatrixLaplacian - User-provided routine to compute the right-hand-side
409:    matrix for the heat equation.

411:    Input Parameters:
412:    ts - the TS context
413:    t - current time  (ignored)
414:    X - current solution (ignored)
415:    dummy - optional user-defined context, as set by TSetRHSJacobian()

417:    Output Parameters:
418:    AA - Jacobian matrix
419:    BB - optionally different matrix from which the preconditioner is built
420:    str - flag indicating matrix structure

422: */
423: PetscErrorCode RHSMatrixLaplaciangllDM(TS ts, PetscReal t, Vec X, Mat A, Mat BB, void *ctx)
424: {
425:   PetscReal **temp;
426:   PetscReal   vv;
427:   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
428:   PetscInt    i, xs, xn, l, j;
429:   PetscInt   *rowsDM;

431:   /*
432:    Creates the element stiffness matrix for the given gll
433:    */
434:   PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp);
435:   /* workaround for clang analyzer warning: Division by zero */

438:   /* scale by the size of the element */
439:   for (i = 0; i < appctx->param.N; i++) {
440:     vv = -appctx->param.mu * 2.0 / appctx->param.Le;
441:     for (j = 0; j < appctx->param.N; j++) temp[i][j] = temp[i][j] * vv;
442:   }

444:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
445:   DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL);

447:   xs = xs / (appctx->param.N - 1);
448:   xn = xn / (appctx->param.N - 1);

450:   PetscMalloc1(appctx->param.N, &rowsDM);
451:   /*
452:    loop over local elements
453:    */
454:   for (j = xs; j < xs + xn; j++) {
455:     for (l = 0; l < appctx->param.N; l++) rowsDM[l] = 1 + (j - xs) * (appctx->param.N - 1) + l;
456:     MatSetValuesLocal(A, appctx->param.N, rowsDM, appctx->param.N, rowsDM, &temp[0][0], ADD_VALUES);
457:   }
458:   PetscFree(rowsDM);
459:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
460:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
461:   VecReciprocal(appctx->SEMop.mass);
462:   MatDiagonalScale(A, appctx->SEMop.mass, 0);
463:   VecReciprocal(appctx->SEMop.mass);

465:   PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp);
466:   return 0;
467: }

469: /*
470:    RHSMatrixAdvection - User-provided routine to compute the right-hand-side
471:    matrix for the Advection equation.

473:    Input Parameters:
474:    ts - the TS context
475:    t - current time
476:    global_in - global input vector
477:    dummy - optional user-defined context, as set by TSetRHSJacobian()

479:    Output Parameters:
480:    AA - Jacobian matrix
481:    BB - optionally different preconditioning matrix
482:    str - flag indicating matrix structure

484: */
485: PetscErrorCode RHSMatrixAdvectiongllDM(TS ts, PetscReal t, Vec X, Mat A, Mat BB, void *ctx)
486: {
487:   PetscReal **temp;
488:   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
489:   PetscInt    xs, xn, l, j;
490:   PetscInt   *rowsDM;

492:   /*
493:    Creates the advection matrix for the given gll
494:    */
495:   PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp);
496:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

498:   DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL);

500:   xs = xs / (appctx->param.N - 1);
501:   xn = xn / (appctx->param.N - 1);

503:   PetscMalloc1(appctx->param.N, &rowsDM);
504:   for (j = xs; j < xs + xn; j++) {
505:     for (l = 0; l < appctx->param.N; l++) rowsDM[l] = 1 + (j - xs) * (appctx->param.N - 1) + l;
506:     MatSetValuesLocal(A, appctx->param.N, rowsDM, appctx->param.N, rowsDM, &temp[0][0], ADD_VALUES);
507:   }
508:   PetscFree(rowsDM);
509:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
510:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

512:   VecReciprocal(appctx->SEMop.mass);
513:   MatDiagonalScale(A, appctx->SEMop.mass, 0);
514:   VecReciprocal(appctx->SEMop.mass);
515:   PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp);
516:   return 0;
517: }
518: /* ------------------------------------------------------------------ */
519: /*
520:    FormFunctionGradient - Evaluates the function and corresponding gradient.

522:    Input Parameters:
523:    tao - the Tao context
524:    IC   - the input vector
525:    ctx - optional user-defined context, as set when calling TaoSetObjectiveAndGradient()

527:    Output Parameters:
528:    f   - the newly evaluated function
529:    G   - the newly evaluated gradient

531:    Notes:

533:           The forward equation is
534:               M u_t = F(U)
535:           which is converted to
536:                 u_t = M^{-1} F(u)
537:           in the user code since TS has no direct way of providing a mass matrix. The Jacobian of this is
538:                  M^{-1} J
539:           where J is the Jacobian of F. Now the adjoint equation is
540:                 M v_t = J^T v
541:           but TSAdjoint does not solve this since it can only solve the transposed system for the
542:           Jacobian the user provided. Hence TSAdjoint solves
543:                  w_t = J^T M^{-1} w  (where w = M v)
544:           since there is no way to indicate the mass matrix as a separate entity to TS. Thus one
545:           must be careful in initializing the "adjoint equation" and using the result. This is
546:           why
547:               G = -2 M(u(T) - u_d)
548:           below (instead of -2(u(T) - u_d) and why the result is
549:               G = G/appctx->SEMop.mass (that is G = M^{-1}w)
550:           below (instead of just the result of the "adjoint solve").

552: */
553: PetscErrorCode FormFunctionGradient(Tao tao, Vec IC, PetscReal *f, Vec G, void *ctx)
554: {
555:   AppCtx            *appctx = (AppCtx *)ctx; /* user-defined application context */
556:   Vec                temp;
557:   PetscInt           its;
558:   PetscReal          ff, gnorm, cnorm, xdiff, errex;
559:   TaoConvergedReason reason;

561:   TSSetTime(appctx->ts, 0.0);
562:   TSSetStepNumber(appctx->ts, 0);
563:   TSSetTimeStep(appctx->ts, appctx->initial_dt);
564:   VecCopy(IC, appctx->dat.curr_sol);

566:   TSSolve(appctx->ts, appctx->dat.curr_sol);

568:   VecWAXPY(G, -1.0, appctx->dat.curr_sol, appctx->dat.obj);

570:   /*
571:      Compute the L2-norm of the objective function, cost function is f
572:   */
573:   VecDuplicate(G, &temp);
574:   VecPointwiseMult(temp, G, G);
575:   VecDot(temp, appctx->SEMop.mass, f);

577:   /* local error evaluation   */
578:   VecWAXPY(temp, -1.0, appctx->dat.ic, appctx->dat.true_solution);
579:   VecPointwiseMult(temp, temp, temp);
580:   /* for error evaluation */
581:   VecDot(temp, appctx->SEMop.mass, &errex);
582:   VecDestroy(&temp);
583:   errex = PetscSqrtReal(errex);

585:   /*
586:      Compute initial conditions for the adjoint integration. See Notes above
587:   */

589:   VecScale(G, -2.0);
590:   VecPointwiseMult(G, G, appctx->SEMop.mass);
591:   TSSetCostGradients(appctx->ts, 1, &G, NULL);
592:   TSAdjointSolve(appctx->ts);
593:   VecPointwiseDivide(G, G, appctx->SEMop.mass);

595:   TaoGetSolutionStatus(tao, &its, &ff, &gnorm, &cnorm, &xdiff, &reason);
596:   return 0;
597: }

599: PetscErrorCode MonitorError(Tao tao, void *ctx)
600: {
601:   AppCtx   *appctx = (AppCtx *)ctx;
602:   Vec       temp;
603:   PetscReal nrm;

605:   VecDuplicate(appctx->dat.ic, &temp);
606:   VecWAXPY(temp, -1.0, appctx->dat.ic, appctx->dat.true_solution);
607:   VecPointwiseMult(temp, temp, temp);
608:   VecDot(temp, appctx->SEMop.mass, &nrm);
609:   VecDestroy(&temp);
610:   nrm = PetscSqrtReal(nrm);
611:   PetscPrintf(PETSC_COMM_WORLD, "Error for initial conditions %g\n", (double)nrm);
612:   return 0;
613: }

615: /*TEST

617:     build:
618:       requires: !complex

620:     test:
621:       args: -tao_max_it 5 -tao_gatol 1.e-4
622:       requires: !single

624:     test:
625:       suffix: 2
626:       nsize: 2
627:       args: -tao_max_it 5 -tao_gatol 1.e-4
628:       requires: !single

630: TEST*/