Actual source code: burgers_spectral.c
1: static char help[] = "Solves a simple data assimilation problem with one dimensional Burger's equation using TSAdjoint\n\n";
3: /*
5: Not yet tested in parallel
7: */
9: /* ------------------------------------------------------------------------
11: This program uses the one-dimensional Burger's equation
12: u_t = mu*u_xx - u u_x,
13: on the domain 0 <= x <= 1, with periodic boundary conditions
15: to demonstrate solving a data assimilation problem of finding the initial conditions
16: to produce a given solution at a fixed time.
18: The operators are discretized with the spectral element method
20: See the paper PDE-CONSTRAINED OPTIMIZATION WITH SPECTRAL ELEMENTS USING PETSC AND TAO
21: by OANA MARIN, EMIL CONSTANTINESCU, AND BARRY SMITH for details on the exact solution
22: used
24: ------------------------------------------------------------------------- */
26: #include <petsctao.h>
27: #include <petscts.h>
28: #include <petscdt.h>
29: #include <petscdmda.h>
31: /*
32: User-defined application context - contains data needed by the
33: application-provided call-back routines.
34: */
36: typedef struct {
37: PetscInt n; /* number of nodes */
38: PetscReal *nodes; /* GLL nodes */
39: PetscReal *weights; /* GLL weights */
40: } PetscGLL;
42: typedef struct {
43: PetscInt N; /* grid points per elements*/
44: PetscInt E; /* number of elements */
45: PetscReal tol_L2, tol_max; /* error norms */
46: PetscInt steps; /* number of timesteps */
47: PetscReal Tend; /* endtime */
48: PetscReal mu; /* viscosity */
49: PetscReal L; /* total length of domain */
50: PetscReal Le;
51: PetscReal Tadj;
52: } PetscParam;
54: typedef struct {
55: Vec obj; /* desired end state */
56: Vec grid; /* total grid */
57: Vec grad;
58: Vec ic;
59: Vec curr_sol;
60: Vec true_solution; /* actual initial conditions for the final solution */
61: } PetscData;
63: typedef struct {
64: Vec grid; /* total grid */
65: Vec mass; /* mass matrix for total integration */
66: Mat stiff; /* stiffness matrix */
67: Mat keptstiff;
68: Mat grad;
69: PetscGLL gll;
70: } PetscSEMOperators;
72: typedef struct {
73: DM da; /* distributed array data structure */
74: PetscSEMOperators SEMop;
75: PetscParam param;
76: PetscData dat;
77: TS ts;
78: PetscReal initial_dt;
79: } AppCtx;
81: /*
82: User-defined routines
83: */
84: extern PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
85: extern PetscErrorCode RHSMatrixLaplaciangllDM(TS, PetscReal, Vec, Mat, Mat, void *);
86: extern PetscErrorCode RHSMatrixAdvectiongllDM(TS, PetscReal, Vec, Mat, Mat, void *);
87: extern PetscErrorCode InitialConditions(Vec, AppCtx *);
88: extern PetscErrorCode TrueSolution(Vec, AppCtx *);
89: extern PetscErrorCode ComputeObjective(PetscReal, Vec, AppCtx *);
90: extern PetscErrorCode MonitorError(Tao, void *);
91: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
92: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
94: int main(int argc, char **argv)
95: {
96: AppCtx appctx; /* user-defined application context */
97: Tao tao;
98: Vec u; /* approximate solution vector */
99: PetscInt i, xs, xm, ind, j, lenglob;
100: PetscReal x, *wrk_ptr1, *wrk_ptr2;
101: MatNullSpace nsp;
102: PetscMPIInt size;
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Initialize program and set problem parameters
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: PetscFunctionBeginUser;
108: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
110: /*initialize parameters */
111: appctx.param.N = 10; /* order of the spectral element */
112: appctx.param.E = 10; /* number of elements */
113: appctx.param.L = 4.0; /* length of the domain */
114: appctx.param.mu = 0.01; /* diffusion coefficient */
115: appctx.initial_dt = 5e-3;
116: appctx.param.steps = PETSC_INT_MAX;
117: appctx.param.Tend = 4;
119: PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &appctx.param.N, NULL));
120: PetscCall(PetscOptionsGetInt(NULL, NULL, "-E", &appctx.param.E, NULL));
121: PetscCall(PetscOptionsGetReal(NULL, NULL, "-Tend", &appctx.param.Tend, NULL));
122: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &appctx.param.mu, NULL));
123: appctx.param.Le = appctx.param.L / appctx.param.E;
125: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
126: PetscCheck((appctx.param.E % size) == 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of elements must be divisible by number of processes");
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Create GLL data structures
130: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: PetscCall(PetscMalloc2(appctx.param.N, &appctx.SEMop.gll.nodes, appctx.param.N, &appctx.SEMop.gll.weights));
132: PetscCall(PetscDTGaussLobattoLegendreQuadrature(appctx.param.N, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, appctx.SEMop.gll.nodes, appctx.SEMop.gll.weights));
133: appctx.SEMop.gll.n = appctx.param.N;
134: lenglob = appctx.param.E * (appctx.param.N - 1);
136: /*
137: Create distributed array (DMDA) to manage parallel grid and vectors
138: and to set up the ghost point communication pattern. There are E*(Nl-1)+1
139: total grid values spread equally among all the processors, except first and last
140: */
142: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, lenglob, 1, 1, NULL, &appctx.da));
143: PetscCall(DMSetFromOptions(appctx.da));
144: PetscCall(DMSetUp(appctx.da));
146: /*
147: Extract global and local vectors from DMDA; we use these to store the
148: approximate solution. Then duplicate these for remaining vectors that
149: have the same types.
150: */
152: PetscCall(DMCreateGlobalVector(appctx.da, &u));
153: PetscCall(VecDuplicate(u, &appctx.dat.ic));
154: PetscCall(VecDuplicate(u, &appctx.dat.true_solution));
155: PetscCall(VecDuplicate(u, &appctx.dat.obj));
156: PetscCall(VecDuplicate(u, &appctx.SEMop.grid));
157: PetscCall(VecDuplicate(u, &appctx.SEMop.mass));
158: PetscCall(VecDuplicate(u, &appctx.dat.curr_sol));
160: PetscCall(DMDAGetCorners(appctx.da, &xs, NULL, NULL, &xm, NULL, NULL));
161: PetscCall(DMDAVecGetArray(appctx.da, appctx.SEMop.grid, &wrk_ptr1));
162: PetscCall(DMDAVecGetArray(appctx.da, appctx.SEMop.mass, &wrk_ptr2));
164: /* Compute function over the locally owned part of the grid */
166: xs = xs / (appctx.param.N - 1);
167: xm = xm / (appctx.param.N - 1);
169: /*
170: Build total grid and mass over entire mesh (multi-elemental)
171: */
173: for (i = xs; i < xs + xm; i++) {
174: for (j = 0; j < appctx.param.N - 1; j++) {
175: x = (appctx.param.Le / 2.0) * (appctx.SEMop.gll.nodes[j] + 1.0) + appctx.param.Le * i;
176: ind = i * (appctx.param.N - 1) + j;
177: wrk_ptr1[ind] = x;
178: wrk_ptr2[ind] = .5 * appctx.param.Le * appctx.SEMop.gll.weights[j];
179: if (j == 0) wrk_ptr2[ind] += .5 * appctx.param.Le * appctx.SEMop.gll.weights[j];
180: }
181: }
182: PetscCall(DMDAVecRestoreArray(appctx.da, appctx.SEMop.grid, &wrk_ptr1));
183: PetscCall(DMDAVecRestoreArray(appctx.da, appctx.SEMop.mass, &wrk_ptr2));
185: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: Create matrix data structure; set matrix evaluation routine.
187: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
188: PetscCall(DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE));
189: PetscCall(DMCreateMatrix(appctx.da, &appctx.SEMop.stiff));
190: PetscCall(DMCreateMatrix(appctx.da, &appctx.SEMop.grad));
191: /*
192: For linear problems with a time-dependent f(u,t) in the equation
193: u_t = f(u,t), the user provides the discretized right-hand side
194: as a time-dependent matrix.
195: */
196: PetscCall(RHSMatrixLaplaciangllDM(appctx.ts, 0.0, u, appctx.SEMop.stiff, appctx.SEMop.stiff, &appctx));
197: PetscCall(RHSMatrixAdvectiongllDM(appctx.ts, 0.0, u, appctx.SEMop.grad, appctx.SEMop.grad, &appctx));
198: /*
199: For linear problems with a time-dependent f(u,t) in the equation
200: u_t = f(u,t), the user provides the discretized right-hand side
201: as a time-dependent matrix.
202: */
204: PetscCall(MatDuplicate(appctx.SEMop.stiff, MAT_COPY_VALUES, &appctx.SEMop.keptstiff));
206: /* attach the null space to the matrix, this probably is not needed but does no harm */
207: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nsp));
208: PetscCall(MatSetNullSpace(appctx.SEMop.stiff, nsp));
209: PetscCall(MatSetNullSpace(appctx.SEMop.keptstiff, nsp));
210: PetscCall(MatNullSpaceTest(nsp, appctx.SEMop.stiff, NULL));
211: PetscCall(MatNullSpaceDestroy(&nsp));
212: /* attach the null space to the matrix, this probably is not needed but does no harm */
213: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nsp));
214: PetscCall(MatSetNullSpace(appctx.SEMop.grad, nsp));
215: PetscCall(MatNullSpaceTest(nsp, appctx.SEMop.grad, NULL));
216: PetscCall(MatNullSpaceDestroy(&nsp));
218: /* Create the TS solver that solves the ODE and its adjoint; set its options */
219: PetscCall(TSCreate(PETSC_COMM_WORLD, &appctx.ts));
220: PetscCall(TSSetProblemType(appctx.ts, TS_NONLINEAR));
221: PetscCall(TSSetType(appctx.ts, TSRK));
222: PetscCall(TSSetDM(appctx.ts, appctx.da));
223: PetscCall(TSSetTime(appctx.ts, 0.0));
224: PetscCall(TSSetTimeStep(appctx.ts, appctx.initial_dt));
225: PetscCall(TSSetMaxSteps(appctx.ts, appctx.param.steps));
226: PetscCall(TSSetMaxTime(appctx.ts, appctx.param.Tend));
227: PetscCall(TSSetExactFinalTime(appctx.ts, TS_EXACTFINALTIME_MATCHSTEP));
228: PetscCall(TSSetTolerances(appctx.ts, 1e-7, NULL, 1e-7, NULL));
229: PetscCall(TSSetFromOptions(appctx.ts));
230: /* Need to save initial timestep user may have set with -ts_time_step so it can be reset for each new TSSolve() */
231: PetscCall(TSGetTimeStep(appctx.ts, &appctx.initial_dt));
232: PetscCall(TSSetRHSFunction(appctx.ts, NULL, RHSFunction, &appctx));
233: PetscCall(TSSetRHSJacobian(appctx.ts, appctx.SEMop.stiff, appctx.SEMop.stiff, RHSJacobian, &appctx));
235: /* Set Objective and Initial conditions for the problem and compute Objective function (evolution of true_solution to final time */
236: PetscCall(InitialConditions(appctx.dat.ic, &appctx));
237: PetscCall(TrueSolution(appctx.dat.true_solution, &appctx));
238: PetscCall(ComputeObjective(appctx.param.Tend, appctx.dat.obj, &appctx));
240: PetscCall(TSSetSaveTrajectory(appctx.ts));
241: PetscCall(TSSetFromOptions(appctx.ts));
243: /* Create TAO solver and set desired solution method */
244: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
245: PetscCall(TaoMonitorSet(tao, MonitorError, &appctx, NULL));
246: PetscCall(TaoSetType(tao, TAOBQNLS));
247: PetscCall(TaoSetSolution(tao, appctx.dat.ic));
248: /* Set routine for function and gradient evaluation */
249: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&appctx));
250: /* Check for any TAO command line options */
251: PetscCall(TaoSetTolerances(tao, 1e-8, PETSC_CURRENT, PETSC_CURRENT));
252: PetscCall(TaoSetFromOptions(tao));
253: PetscCall(TaoSolve(tao));
255: PetscCall(TaoDestroy(&tao));
256: PetscCall(MatDestroy(&appctx.SEMop.stiff));
257: PetscCall(MatDestroy(&appctx.SEMop.keptstiff));
258: PetscCall(MatDestroy(&appctx.SEMop.grad));
259: PetscCall(VecDestroy(&u));
260: PetscCall(VecDestroy(&appctx.dat.ic));
261: PetscCall(VecDestroy(&appctx.dat.true_solution));
262: PetscCall(VecDestroy(&appctx.dat.obj));
263: PetscCall(VecDestroy(&appctx.SEMop.grid));
264: PetscCall(VecDestroy(&appctx.SEMop.mass));
265: PetscCall(VecDestroy(&appctx.dat.curr_sol));
266: PetscCall(PetscFree2(appctx.SEMop.gll.nodes, appctx.SEMop.gll.weights));
267: PetscCall(DMDestroy(&appctx.da));
268: PetscCall(TSDestroy(&appctx.ts));
270: /*
271: Always call PetscFinalize() before exiting a program. This routine
272: - finalizes the PETSc libraries as well as MPI
273: - provides summary and diagnostic information if certain runtime
274: options are chosen (e.g., -log_view).
275: */
276: PetscCall(PetscFinalize());
277: return 0;
278: }
280: /* --------------------------------------------------------------------- */
281: /*
282: InitialConditions - Computes the initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve()
284: The routine TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function
286: Input Parameter:
287: u - uninitialized solution vector (global)
288: appctx - user-defined application context
290: Output Parameter:
291: u - vector with solution at initial time (global)
292: */
293: PetscErrorCode InitialConditions(Vec u, AppCtx *appctx)
294: {
295: PetscScalar *s;
296: const PetscScalar *xg;
297: PetscInt i, xs, xn;
299: PetscFunctionBegin;
300: PetscCall(DMDAVecGetArray(appctx->da, u, &s));
301: PetscCall(DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
302: PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));
303: 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));
304: PetscCall(DMDAVecRestoreArray(appctx->da, u, &s));
305: PetscCall(DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*
310: TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function.
312: InitialConditions() computes the initial conditions for the beginning of the Tao iterations
314: Input Parameter:
315: u - uninitialized solution vector (global)
316: appctx - user-defined application context
318: Output Parameter:
319: u - vector with solution at initial time (global)
320: */
321: PetscErrorCode TrueSolution(Vec u, AppCtx *appctx)
322: {
323: PetscScalar *s;
324: const PetscScalar *xg;
325: PetscInt i, xs, xn;
327: PetscFunctionBegin;
328: PetscCall(DMDAVecGetArray(appctx->da, u, &s));
329: PetscCall(DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
330: PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));
331: 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]));
332: PetscCall(DMDAVecRestoreArray(appctx->da, u, &s));
333: PetscCall(DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
336: /* --------------------------------------------------------------------- */
337: /*
338: Sets the desired profile for the final end time
340: Input Parameters:
341: t - final time
342: obj - vector storing the desired profile
343: appctx - user-defined application context
345: */
346: PetscErrorCode ComputeObjective(PetscReal t, Vec obj, AppCtx *appctx)
347: {
348: PetscScalar *s;
349: const PetscScalar *xg;
350: PetscInt i, xs, xn;
352: PetscFunctionBegin;
353: PetscCall(DMDAVecGetArray(appctx->da, obj, &s));
354: PetscCall(DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
355: PetscCall(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: PetscCall(DMDAVecRestoreArray(appctx->da, obj, &s));
360: PetscCall(DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, PetscCtx ctx)
365: {
366: AppCtx *appctx = (AppCtx *)ctx;
368: PetscFunctionBegin;
369: PetscCall(MatMult(appctx->SEMop.grad, globalin, globalout)); /* grad u */
370: PetscCall(VecPointwiseMult(globalout, globalin, globalout)); /* u grad u */
371: PetscCall(VecScale(globalout, -1.0));
372: PetscCall(MatMultAdd(appctx->SEMop.keptstiff, globalin, globalout, globalout));
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: /*
378: K is the discretiziation of the Laplacian
379: G is the discretization of the gradient
381: Computes Jacobian of K u + diag(u) G u which is given by
382: K + diag(u)G + diag(Gu)
383: */
384: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec globalin, Mat A, Mat B, PetscCtx ctx)
385: {
386: AppCtx *appctx = (AppCtx *)ctx;
387: Vec Gglobalin;
389: PetscFunctionBegin;
390: /* A = diag(u) G */
392: PetscCall(MatCopy(appctx->SEMop.grad, A, SAME_NONZERO_PATTERN));
393: PetscCall(MatDiagonalScale(A, globalin, NULL));
395: /* A = A + diag(Gu) */
396: PetscCall(VecDuplicate(globalin, &Gglobalin));
397: PetscCall(MatMult(appctx->SEMop.grad, globalin, Gglobalin));
398: PetscCall(MatDiagonalSet(A, Gglobalin, ADD_VALUES));
399: PetscCall(VecDestroy(&Gglobalin));
401: /* A = K - A */
402: PetscCall(MatScale(A, -1.0));
403: PetscCall(MatAXPY(A, 1.0, appctx->SEMop.keptstiff, SAME_NONZERO_PATTERN));
404: PetscFunctionReturn(PETSC_SUCCESS);
405: }
407: /* --------------------------------------------------------------------- */
409: /*
410: RHSMatrixLaplacian - User-provided routine to compute the right-hand-side
411: matrix for the heat equation.
413: Input Parameters:
414: ts - the TS context
415: t - current time (ignored)
416: X - current solution (ignored)
417: dummy - optional user-defined context, as set by TSetRHSJacobian()
419: Output Parameters:
420: AA - Jacobian matrix
421: BB - optionally different matrix from which the preconditioner is built
423: */
424: PetscErrorCode RHSMatrixLaplaciangllDM(TS ts, PetscReal t, Vec X, Mat A, Mat BB, PetscCtx ctx)
425: {
426: PetscReal **temp;
427: PetscReal vv;
428: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
429: PetscInt i, xs, xn, l, j;
430: PetscInt *rowsDM;
432: PetscFunctionBegin;
433: /*
434: Creates the element stiffness matrix for the given gll
435: */
436: PetscCall(PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
437: /* workaround for clang analyzer warning: Division by zero */
438: PetscCheck(appctx->param.N > 1, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Spectral element order should be > 1");
440: /* scale by the size of the element */
441: for (i = 0; i < appctx->param.N; i++) {
442: vv = -appctx->param.mu * 2.0 / appctx->param.Le;
443: for (j = 0; j < appctx->param.N; j++) temp[i][j] = temp[i][j] * vv;
444: }
446: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
447: PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));
449: xs = xs / (appctx->param.N - 1);
450: xn = xn / (appctx->param.N - 1);
452: PetscCall(PetscMalloc1(appctx->param.N, &rowsDM));
453: /*
454: loop over local elements
455: */
456: for (j = xs; j < xs + xn; j++) {
457: for (l = 0; l < appctx->param.N; l++) rowsDM[l] = 1 + (j - xs) * (appctx->param.N - 1) + l;
458: PetscCall(MatSetValuesLocal(A, appctx->param.N, rowsDM, appctx->param.N, rowsDM, &temp[0][0], ADD_VALUES));
459: }
460: PetscCall(PetscFree(rowsDM));
461: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
462: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
463: PetscCall(VecReciprocal(appctx->SEMop.mass));
464: PetscCall(MatDiagonalScale(A, appctx->SEMop.mass, NULL));
465: PetscCall(VecReciprocal(appctx->SEMop.mass));
467: PetscCall(PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /*
472: RHSMatrixAdvection - User-provided routine to compute the right-hand-side
473: matrix for the Advection equation.
475: Input Parameters:
476: ts - the TS context
477: t - current time
478: global_in - global input vector
479: dummy - optional user-defined context, as set by TSetRHSJacobian()
481: Output Parameters:
482: AA - Jacobian matrix
483: BB - optionally different matrix used to compute the preconditioner
485: */
486: PetscErrorCode RHSMatrixAdvectiongllDM(TS ts, PetscReal t, Vec X, Mat A, Mat BB, PetscCtx ctx)
487: {
488: PetscReal **temp;
489: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
490: PetscInt xs, xn, l, j;
491: PetscInt *rowsDM;
493: PetscFunctionBegin;
494: /*
495: Creates the advection matrix for the given gll
496: */
497: PetscCall(PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
498: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
500: PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));
502: xs = xs / (appctx->param.N - 1);
503: xn = xn / (appctx->param.N - 1);
505: PetscCall(PetscMalloc1(appctx->param.N, &rowsDM));
506: for (j = xs; j < xs + xn; j++) {
507: for (l = 0; l < appctx->param.N; l++) rowsDM[l] = 1 + (j - xs) * (appctx->param.N - 1) + l;
508: PetscCall(MatSetValuesLocal(A, appctx->param.N, rowsDM, appctx->param.N, rowsDM, &temp[0][0], ADD_VALUES));
509: }
510: PetscCall(PetscFree(rowsDM));
511: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
512: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
514: PetscCall(VecReciprocal(appctx->SEMop.mass));
515: PetscCall(MatDiagonalScale(A, appctx->SEMop.mass, NULL));
516: PetscCall(VecReciprocal(appctx->SEMop.mass));
517: PetscCall(PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
520: /* ------------------------------------------------------------------ */
521: /*
522: FormFunctionGradient - Evaluates the function and corresponding gradient.
524: Input Parameters:
525: tao - the Tao context
526: IC - the input vector
527: ctx - optional user-defined context, as set when calling TaoSetObjectiveAndGradient()
529: Output Parameters:
530: f - the newly evaluated function
531: G - the newly evaluated gradient
533: Notes:
535: The forward equation is
536: M u_t = F(U)
537: which is converted to
538: u_t = M^{-1} F(u)
539: in the user code since TS has no direct way of providing a mass matrix. The Jacobian of this is
540: M^{-1} J
541: where J is the Jacobian of F. Now the adjoint equation is
542: M v_t = J^T v
543: but TSAdjoint does not solve this since it can only solve the transposed system for the
544: Jacobian the user provided. Hence TSAdjoint solves
545: w_t = J^T M^{-1} w (where w = M v)
546: since there is no way to indicate the mass matrix as a separate entity to TS. Thus one
547: must be careful in initializing the "adjoint equation" and using the result. This is
548: why
549: G = -2 M(u(T) - u_d)
550: below (instead of -2(u(T) - u_d) and why the result is
551: G = G/appctx->SEMop.mass (that is G = M^{-1}w)
552: below (instead of just the result of the "adjoint solve").
554: */
555: PetscErrorCode FormFunctionGradient(Tao tao, Vec IC, PetscReal *f, Vec G, PetscCtx ctx)
556: {
557: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
558: Vec temp;
559: PetscInt its;
560: PetscReal ff, gnorm, cnorm, xdiff, errex;
561: TaoConvergedReason reason;
563: PetscFunctionBegin;
564: PetscCall(TSSetTime(appctx->ts, 0.0));
565: PetscCall(TSSetStepNumber(appctx->ts, 0));
566: PetscCall(TSSetTimeStep(appctx->ts, appctx->initial_dt));
567: PetscCall(VecCopy(IC, appctx->dat.curr_sol));
569: PetscCall(TSSolve(appctx->ts, appctx->dat.curr_sol));
571: PetscCall(VecWAXPY(G, -1.0, appctx->dat.curr_sol, appctx->dat.obj));
573: /*
574: Compute the L2-norm of the objective function, cost function is f
575: */
576: PetscCall(VecDuplicate(G, &temp));
577: PetscCall(VecPointwiseMult(temp, G, G));
578: PetscCall(VecDot(temp, appctx->SEMop.mass, f));
580: /* local error evaluation */
581: PetscCall(VecWAXPY(temp, -1.0, appctx->dat.ic, appctx->dat.true_solution));
582: PetscCall(VecPointwiseMult(temp, temp, temp));
583: /* for error evaluation */
584: PetscCall(VecDot(temp, appctx->SEMop.mass, &errex));
585: PetscCall(VecDestroy(&temp));
586: errex = PetscSqrtReal(errex);
588: /*
589: Compute initial conditions for the adjoint integration. See Notes above
590: */
592: PetscCall(VecScale(G, -2.0));
593: PetscCall(VecPointwiseMult(G, G, appctx->SEMop.mass));
594: PetscCall(TSSetCostGradients(appctx->ts, 1, &G, NULL));
595: PetscCall(TSAdjointSolve(appctx->ts));
596: PetscCall(VecPointwiseDivide(G, G, appctx->SEMop.mass));
598: PetscCall(TaoGetSolutionStatus(tao, &its, &ff, &gnorm, &cnorm, &xdiff, &reason));
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: PetscErrorCode MonitorError(Tao tao, PetscCtx ctx)
603: {
604: AppCtx *appctx = (AppCtx *)ctx;
605: Vec temp;
606: PetscReal nrm;
608: PetscFunctionBegin;
609: PetscCall(VecDuplicate(appctx->dat.ic, &temp));
610: PetscCall(VecWAXPY(temp, -1.0, appctx->dat.ic, appctx->dat.true_solution));
611: PetscCall(VecPointwiseMult(temp, temp, temp));
612: PetscCall(VecDot(temp, appctx->SEMop.mass, &nrm));
613: PetscCall(VecDestroy(&temp));
614: nrm = PetscSqrtReal(nrm);
615: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error for initial conditions %g\n", (double)nrm));
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*TEST
621: build:
622: requires: !complex
624: test:
625: args: -tao_max_it 5 -tao_gatol 1.e-4
626: requires: !single
628: test:
629: requires: !defined(PETSCTEST_VALGRIND)
630: suffix: 2
631: nsize: 2
632: args: -tao_max_it 5 -tao_gatol 1.e-4
633: requires: !single
635: TEST*/