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*/