Actual source code: spectraladjointassimilation.c
petsc-3.11.4 2019-09-28
2: static char help[] ="Solves a simple data assimilation problem with one dimensional advection diffusion equation using TSAdjoint\n\n";
4: /*
6: Not yet tested in parallel
8: */
9: /*
10: Concepts: TS^time-dependent linear problems
11: Concepts: TS^heat equation
12: Concepts: TS^diffusion equation
13: Concepts: adjoints
14: Processors: n
15: */
17: /* ------------------------------------------------------------------------
19: This program uses the one-dimensional advection-diffusion equation),
20: u_t = mu*u_xx - a u_x,
21: on the domain 0 <= x <= 1, with periodic boundary conditions
23: to demonstrate solving a data assimilation problem of finding the initial conditions
24: to produce a given solution at a fixed time.
26: The operators are discretized with the spectral element method
28: ------------------------------------------------------------------------- */
30: /*
31: Include "petscts.h" so that we can use TS solvers. Note that this file
32: automatically includes:
33: petscsys.h - base PETSc routines petscvec.h - vectors
34: petscmat.h - matrices
35: petscis.h - index sets petscksp.h - Krylov subspace methods
36: petscviewer.h - viewers petscpc.h - preconditioners
37: petscksp.h - linear solvers petscsnes.h - nonlinear solvers
38: */
40: #include <petsctao.h>
41: #include <petscts.h>
42: #include <petscgll.h>
43: #include <petscdraw.h>
44: #include <petscdmda.h>
46: /*
47: User-defined application context - contains data needed by the
48: application-provided call-back routines.
49: */
51: typedef struct {
52: PetscInt N; /* grid points per elements*/
53: PetscInt E; /* number of elements */
54: PetscReal tol_L2,tol_max; /* error norms */
55: PetscInt steps; /* number of timesteps */
56: PetscReal Tend; /* endtime */
57: PetscReal mu; /* viscosity */
58: PetscReal a; /* advection speed */
59: PetscReal L; /* total length of domain */
60: PetscReal Le;
61: PetscReal Tadj;
62: } PetscParam;
64: typedef struct {
65: Vec reference; /* desired end state */
66: Vec grid; /* total grid */
67: Vec grad;
68: Vec ic;
69: Vec curr_sol;
70: Vec joe;
71: Vec true_solution; /* actual initial conditions for the final solution */
72: } PetscData;
74: typedef struct {
75: Vec grid; /* total grid */
76: Vec mass; /* mass matrix for total integration */
77: Mat stiff; /* stifness matrix */
78: Mat advec;
79: Mat keptstiff;
80: PetscGLL gll;
81: } PetscSEMOperators;
83: typedef struct {
84: DM da; /* distributed array data structure */
85: PetscSEMOperators SEMop;
86: PetscParam param;
87: PetscData dat;
88: TS ts;
89: PetscReal initial_dt;
90: PetscReal *solutioncoefficients;
91: PetscInt ncoeff;
92: } AppCtx;
94: /*
95: User-defined routines
96: */
97: extern PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
98: extern PetscErrorCode RHSLaplacian(TS,PetscReal,Vec,Mat,Mat,void*);
99: extern PetscErrorCode RHSAdvection(TS,PetscReal,Vec,Mat,Mat,void*);
100: extern PetscErrorCode InitialConditions(Vec,AppCtx*);
101: extern PetscErrorCode ComputeReference(TS,PetscReal,Vec,AppCtx*);
102: extern PetscErrorCode MonitorError(Tao,void*);
103: extern PetscErrorCode MonitorDestroy(void**);
104: extern PetscErrorCode ComputeSolutionCoefficients(AppCtx*);
105: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
106: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
108: int main(int argc,char **argv)
109: {
110: AppCtx appctx; /* user-defined application context */
111: Tao tao;
112: Vec u; /* approximate solution vector */
114: PetscInt i, xs, xm, ind, j, lenglob;
115: PetscReal x, *wrk_ptr1, *wrk_ptr2;
116: MatNullSpace nsp;
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Initialize program and set problem parameters
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
125: /*initialize parameters */
126: appctx.param.N = 10; /* order of the spectral element */
127: appctx.param.E = 8; /* number of elements */
128: appctx.param.L = 1.0; /* length of the domain */
129: appctx.param.mu = 0.00001; /* diffusion coefficient */
130: appctx.param.a = 0.0; /* advection speed */
131: appctx.initial_dt = 1e-4;
132: appctx.param.steps = PETSC_MAX_INT;
133: appctx.param.Tend = 0.01;
134: appctx.ncoeff = 2;
136: PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL);
137: PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL);
138: PetscOptionsGetInt(NULL,NULL,"-ncoeff",&appctx.ncoeff,NULL);
139: PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL);
140: PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL);
141: PetscOptionsGetReal(NULL,NULL,"-a",&appctx.param.a,NULL);
142: appctx.param.Le = appctx.param.L/appctx.param.E;
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Create GLL data structures
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: PetscGLLCreate(appctx.param.N,PETSCGLL_VIA_LINEARALGEBRA,&appctx.SEMop.gll);
149: lenglob = appctx.param.E*(appctx.param.N-1);
151: /*
152: Create distributed array (DMDA) to manage parallel grid and vectors
153: and to set up the ghost point communication pattern. There are E*(Nl-1)+1
154: total grid values spread equally among all the processors, except first and last
155: */
157: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da);
158: DMSetFromOptions(appctx.da);
159: DMSetUp(appctx.da);
161: /*
162: Extract global and local vectors from DMDA; we use these to store the
163: approximate solution. Then duplicate these for remaining vectors that
164: have the same types.
165: */
167: DMCreateGlobalVector(appctx.da,&u);
168: VecDuplicate(u,&appctx.dat.ic);
169: VecDuplicate(u,&appctx.dat.true_solution);
170: VecDuplicate(u,&appctx.dat.reference);
171: VecDuplicate(u,&appctx.SEMop.grid);
172: VecDuplicate(u,&appctx.SEMop.mass);
173: VecDuplicate(u,&appctx.dat.curr_sol);
174: VecDuplicate(u,&appctx.dat.joe);
176: DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
177: DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
178: DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
180: /* Compute function over the locally owned part of the grid */
182: xs=xs/(appctx.param.N-1);
183: xm=xm/(appctx.param.N-1);
185: /*
186: Build total grid and mass over entire mesh (multi-elemental)
187: */
189: for (i=xs; i<xs+xm; i++) {
190: for (j=0; j<appctx.param.N-1; j++) {
191: x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i;
192: ind=i*(appctx.param.N-1)+j;
193: wrk_ptr1[ind]=x;
194: wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
195: if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j];
196: }
197: }
198: DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1);
199: DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2);
202: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: Create matrix data structure; set matrix evaluation routine.
204: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
205: DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE);
206: DMCreateMatrix(appctx.da,&appctx.SEMop.stiff);
207: DMCreateMatrix(appctx.da,&appctx.SEMop.advec);
209: /*
210: For linear problems with a time-dependent f(u,t) in the equation
211: u_t = f(u,t), the user provides the discretized right-hand-side
212: as a time-dependent matrix.
213: */
214: RHSLaplacian(appctx.ts,0.0,u,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx);
215: RHSAdvection(appctx.ts,0.0,u,appctx.SEMop.advec,appctx.SEMop.advec,&appctx);
216: MatAXPY(appctx.SEMop.stiff,-1.0,appctx.SEMop.advec,DIFFERENT_NONZERO_PATTERN);
217: MatDuplicate(appctx.SEMop.stiff,MAT_COPY_VALUES,&appctx.SEMop.keptstiff);
219: /* attach the null space to the matrix, this probably is not needed but does no harm */
220: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp);
221: MatSetNullSpace(appctx.SEMop.stiff,nsp);
222: MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL);
223: MatNullSpaceDestroy(&nsp);
225: /* Create the TS solver that solves the ODE and its adjoint; set its options */
226: TSCreate(PETSC_COMM_WORLD,&appctx.ts);
227: TSSetSolutionFunction(appctx.ts,(PetscErrorCode (*)(TS,PetscReal,Vec, void *))ComputeReference,&appctx);
228: TSSetProblemType(appctx.ts,TS_LINEAR);
229: TSSetType(appctx.ts,TSRK);
230: TSSetDM(appctx.ts,appctx.da);
231: TSSetTime(appctx.ts,0.0);
232: TSSetTimeStep(appctx.ts,appctx.initial_dt);
233: TSSetMaxSteps(appctx.ts,appctx.param.steps);
234: TSSetMaxTime(appctx.ts,appctx.param.Tend);
235: TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
236: TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL);
237: TSSetFromOptions(appctx.ts);
238: /* Need to save initial timestep user may have set with -ts_dt so it can be reset for each new TSSolve() */
239: TSGetTimeStep(appctx.ts,&appctx.initial_dt);
240: TSSetRHSFunction(appctx.ts,NULL,TSComputeRHSFunctionLinear,&appctx);
241: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,TSComputeRHSJacobianConstant,&appctx);
242: /* TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
243: TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx); */
244: TSSetSaveTrajectory(appctx.ts);
246: /* Set random initial conditions as initial guess, compute analytic reference solution and analytic (true) initial conditions */
247: ComputeSolutionCoefficients(&appctx);
248: InitialConditions(appctx.dat.ic,&appctx);
249: ComputeReference(appctx.ts,appctx.param.Tend,appctx.dat.reference,&appctx);
250: ComputeReference(appctx.ts,0.0,appctx.dat.true_solution,&appctx);
252: /* Create TAO solver and set desired solution method */
253: TaoCreate(PETSC_COMM_WORLD,&tao);
254: TaoSetMonitor(tao,MonitorError,&appctx,MonitorDestroy);
255: TaoSetType(tao,TAOBQNLS);
256: TaoSetInitialVector(tao,appctx.dat.ic);
257: /* Set routine for function and gradient evaluation */
258: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&appctx);
259: /* Check for any TAO command line options */
260: TaoSetTolerances(tao,1e-8,PETSC_DEFAULT,PETSC_DEFAULT);
261: TaoSetFromOptions(tao);
262: TaoSolve(tao);
264: TaoDestroy(&tao);
265: PetscFree(appctx.solutioncoefficients);
266: MatDestroy(&appctx.SEMop.advec);
267: MatDestroy(&appctx.SEMop.stiff);
268: MatDestroy(&appctx.SEMop.keptstiff);
269: VecDestroy(&u);
270: VecDestroy(&appctx.dat.ic);
271: VecDestroy(&appctx.dat.joe);
272: VecDestroy(&appctx.dat.true_solution);
273: VecDestroy(&appctx.dat.reference);
274: VecDestroy(&appctx.SEMop.grid);
275: VecDestroy(&appctx.SEMop.mass);
276: VecDestroy(&appctx.dat.curr_sol);
277: PetscGLLDestroy(&appctx.SEMop.gll);
278: DMDestroy(&appctx.da);
279: TSDestroy(&appctx.ts);
281: /*
282: Always call PetscFinalize() before exiting a program. This routine
283: - finalizes the PETSc libraries as well as MPI
284: - provides summary and diagnostic information if certain runtime
285: options are chosen (e.g., -log_summary).
286: */
287: PetscFinalize();
288: return ierr;
289: }
291: /*
292: Computes the coefficients for the analytic solution to the PDE
293: */
294: PetscErrorCode ComputeSolutionCoefficients(AppCtx *appctx)
295: {
296: PetscErrorCode ierr;
297: PetscRandom rand;
298: PetscInt i;
301: PetscMalloc1(appctx->ncoeff,&appctx->solutioncoefficients);
302: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
303: PetscRandomSetInterval(rand,.9,1.0);
304: for (i=0; i<appctx->ncoeff; i++) {
305: PetscRandomGetValue(rand,&appctx->solutioncoefficients[i]);
306: }
307: PetscRandomDestroy(&rand);
308: return(0);
309: }
311: /* --------------------------------------------------------------------- */
312: /*
313: InitialConditions - Computes the (random) initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve()
315: Input Parameter:
316: u - uninitialized solution vector (global)
317: appctx - user-defined application context
319: Output Parameter:
320: u - vector with solution at initial time (global)
321: */
322: PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
323: {
324: PetscScalar *s;
325: const PetscScalar *xg;
326: PetscErrorCode ierr;
327: PetscInt i,j,lenglob;
328: PetscReal sum,val;
329: PetscRandom rand;
332: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
333: PetscRandomSetInterval(rand,.9,1.0);
334: DMDAVecGetArray(appctx->da,u,&s);
335: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
336: lenglob = appctx->param.E*(appctx->param.N-1);
337: for (i=0; i<lenglob; i++) {
338: s[i]= 0;
339: for (j=0; j<appctx->ncoeff; j++) {
340: PetscRandomGetValue(rand,&val);
341: s[i] += val*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
342: }
343: }
344: PetscRandomDestroy(&rand);
345: DMDAVecRestoreArray(appctx->da,u,&s);
346: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
347: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
348: VecSum(u,&sum);
349: VecShift(u,-sum/lenglob);
350: return(0);
351: }
354: /*
355: TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function.
357: InitialConditions() computes the initial conditions for the begining of the Tao iterations
359: Input Parameter:
360: u - uninitialized solution vector (global)
361: appctx - user-defined application context
363: Output Parameter:
364: u - vector with solution at initial time (global)
365: */
366: PetscErrorCode TrueSolution(Vec u,AppCtx *appctx)
367: {
368: PetscScalar *s;
369: const PetscScalar *xg;
370: PetscErrorCode ierr;
371: PetscInt i,j,lenglob;
372: PetscReal sum;
375: DMDAVecGetArray(appctx->da,u,&s);
376: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
377: lenglob = appctx->param.E*(appctx->param.N-1);
378: for (i=0; i<lenglob; i++) {
379: s[i]= 0;
380: for (j=0; j<appctx->ncoeff; j++) {
381: s[i] += appctx->solutioncoefficients[j]*PetscSinScalar(2*(j+1)*PETSC_PI*xg[i]);
382: }
383: }
384: DMDAVecRestoreArray(appctx->da,u,&s);
385: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
386: /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
387: VecSum(u,&sum);
388: VecShift(u,-sum/lenglob);
389: return(0);
390: }
391: /* --------------------------------------------------------------------- */
392: /*
393: Sets the desired profile for the final end time
395: Input Parameters:
396: t - final time
397: obj - vector storing the desired profile
398: appctx - user-defined application context
400: */
401: PetscErrorCode ComputeReference(TS ts,PetscReal t,Vec obj,AppCtx *appctx)
402: {
403: PetscScalar *s,tc;
404: const PetscScalar *xg;
405: PetscErrorCode ierr;
406: PetscInt i, j,lenglob;
409: DMDAVecGetArray(appctx->da,obj,&s);
410: DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
411: lenglob = appctx->param.E*(appctx->param.N-1);
412: for (i=0; i<lenglob; i++) {
413: s[i]= 0;
414: for (j=0; j<appctx->ncoeff; j++) {
415: tc = -appctx->param.mu*(j+1)*(j+1)*4.0*PETSC_PI*PETSC_PI*t;
416: s[i] += appctx->solutioncoefficients[j]*PetscSinScalar(2*(j+1)*PETSC_PI*(xg[i] + appctx->param.a*t))*PetscExpReal(tc);
417: }
418: }
419: DMDAVecRestoreArray(appctx->da,obj,&s);
420: DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg);
421: return(0);
422: }
424: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
425: {
427: AppCtx *appctx = (AppCtx*)ctx;
430: MatMult(appctx->SEMop.keptstiff,globalin,globalout);
431: return(0);
432: }
434: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec globalin,Mat A, Mat B,void *ctx)
435: {
437: AppCtx *appctx = (AppCtx*)ctx;
440: MatCopy(appctx->SEMop.keptstiff,A,DIFFERENT_NONZERO_PATTERN);
441: return(0);
442: }
444: /* --------------------------------------------------------------------- */
446: /*
447: RHSLaplacian - matrix for diffusion
449: Input Parameters:
450: ts - the TS context
451: t - current time (ignored)
452: X - current solution (ignored)
453: dummy - optional user-defined context, as set by TSetRHSJacobian()
455: Output Parameters:
456: AA - Jacobian matrix
457: BB - optionally different matrix from which the preconditioner is built
458: str - flag indicating matrix structure
460: Scales by the inverse of the mass matrix (perhaps that should be pulled out)
462: */
463: PetscErrorCode RHSLaplacian(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
464: {
465: PetscReal **temp;
466: PetscReal vv;
467: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
469: PetscInt i,xs,xn,l,j;
470: PetscInt *rowsDM;
473: /*
474: Creates the element stiffness matrix for the given gll
475: */
476: PetscGLLElementLaplacianCreate(&appctx->SEMop.gll,&temp);
478: /* scale by the size of the element */
479: for (i=0; i<appctx->param.N; i++) {
480: vv=-appctx->param.mu*2.0/appctx->param.Le;
481: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
482: }
484: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
485: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
487: if (appctx->param.N-1 < 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Polynomial order must be at least 2");
488: xs = xs/(appctx->param.N-1);
489: xn = xn/(appctx->param.N-1);
491: PetscMalloc1(appctx->param.N,&rowsDM);
492: /*
493: loop over local elements
494: */
495: for (j=xs; j<xs+xn; j++) {
496: for (l=0; l<appctx->param.N; l++) {
497: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
498: }
499: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
500: }
501: PetscFree(rowsDM);
502: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
503: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
504: VecReciprocal(appctx->SEMop.mass);
505: MatDiagonalScale(A,appctx->SEMop.mass,0);
506: VecReciprocal(appctx->SEMop.mass);
508: PetscGLLElementLaplacianDestroy(&appctx->SEMop.gll,&temp);
509: return(0);
510: }
512: /*
513: Almost identical to Laplacian
515: Note that the element matrix is NOT scaled by the size of element like the Laplacian term.
516: */
517: PetscErrorCode RHSAdvection(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx)
518: {
519: PetscReal **temp;
520: PetscReal vv;
521: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
523: PetscInt i,xs,xn,l,j;
524: PetscInt *rowsDM;
527: /*
528: Creates the element stiffness matrix for the given gll
529: */
530: PetscGLLElementAdvectionCreate(&appctx->SEMop.gll,&temp);
532: /* scale by the size of the element */
533: for (i=0; i<appctx->param.N; i++) {
534: vv = -appctx->param.a;
535: for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv;
536: }
538: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
539: DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL);
541: if (appctx->param.N-1 < 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Polynomial order must be at least 2");
542: xs = xs/(appctx->param.N-1);
543: xn = xn/(appctx->param.N-1);
545: PetscMalloc1(appctx->param.N,&rowsDM);
546: /*
547: loop over local elements
548: */
549: for (j=xs; j<xs+xn; j++) {
550: for (l=0; l<appctx->param.N; l++) {
551: rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l;
552: }
553: MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES);
554: }
555: PetscFree(rowsDM);
556: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
557: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
558: VecReciprocal(appctx->SEMop.mass);
559: MatDiagonalScale(A,appctx->SEMop.mass,0);
560: VecReciprocal(appctx->SEMop.mass);
562: PetscGLLElementAdvectionDestroy(&appctx->SEMop.gll,&temp);
563: return(0);
564: }
566: /* ------------------------------------------------------------------ */
567: /*
568: FormFunctionGradient - Evaluates the function and corresponding gradient.
570: Input Parameters:
571: tao - the Tao context
572: ic - the input vector
573: ctx - optional user-defined context, as set when calling TaoSetObjectiveAndGradientRoutine()
575: Output Parameters:
576: f - the newly evaluated function
577: G - the newly evaluated gradient
579: Notes:
581: The forward equation is
582: M u_t = F(U)
583: which is converted to
584: u_t = M^{-1} F(u)
585: in the user code since TS has no direct way of providing a mass matrix. The Jacobian of this is
586: M^{-1} J
587: where J is the Jacobian of F. Now the adjoint equation is
588: M v_t = J^T v
589: but TSAdjoint does not solve this since it can only solve the transposed system for the
590: Jacobian the user provided. Hence TSAdjoint solves
591: w_t = J^T M^{-1} w (where w = M v)
592: since there is no way to indicate the mass matrix as a seperate entitity to TS. Thus one
593: must be careful in initializing the "adjoint equation" and using the result. This is
594: why
595: G = -2 M(u(T) - u_d)
596: below (instead of -2(u(T) - u_d)
599: */
600: PetscErrorCode FormFunctionGradient(Tao tao,Vec ic,PetscReal *f,Vec G,void *ctx)
601: {
602: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
603: PetscErrorCode ierr;
604: Vec temp;
607: TSSetTime(appctx->ts,0.0);
608: TSSetStepNumber(appctx->ts,0);
609: TSResetTrajectory(appctx->ts);
610: TSSetTimeStep(appctx->ts,appctx->initial_dt);
611: VecCopy(ic,appctx->dat.curr_sol);
613: TSSolve(appctx->ts,appctx->dat.curr_sol);
614: VecCopy(appctx->dat.curr_sol,appctx->dat.joe);
616: /* Compute the difference between the current ODE solution and target ODE solution */
617: VecWAXPY(G,-1.0,appctx->dat.curr_sol,appctx->dat.reference);
619: /* Compute the objective/cost function */
620: VecDuplicate(G,&temp);
621: VecPointwiseMult(temp,G,G);
622: VecDot(temp,appctx->SEMop.mass,f);
623: VecDestroy(&temp);
625: /* Compute initial conditions for the adjoint integration. See Notes above */
626: VecScale(G, -2.0);
627: VecPointwiseMult(G,G,appctx->SEMop.mass);
628: TSSetCostGradients(appctx->ts,1,&G,NULL);
630: TSAdjointSolve(appctx->ts);
631: /* VecPointwiseDivide(G,G,appctx->SEMop.mass);*/
632: return(0);
633: }
635: PetscErrorCode MonitorError(Tao tao,void *ctx)
636: {
637: AppCtx *appctx = (AppCtx*)ctx;
638: Vec temp,grad;
639: PetscReal nrm;
641: PetscInt its;
642: PetscReal fct,gnorm;
645: VecDuplicate(appctx->dat.ic,&temp);
646: VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution);
647: VecPointwiseMult(temp,temp,temp);
648: VecDot(temp,appctx->SEMop.mass,&nrm);
649: nrm = PetscSqrtReal(nrm);
650: TaoGetGradientVector(tao,&grad);
651: VecPointwiseMult(temp,temp,temp);
652: VecDot(temp,appctx->SEMop.mass,&gnorm);
653: gnorm = PetscSqrtReal(gnorm);
654: VecDestroy(&temp);
655: TaoGetIterationNumber(tao,&its);
656: TaoGetObjective(tao,&fct);
657: if (!its) {
658: PetscPrintf(PETSC_COMM_WORLD,"%% Iteration Error Objective Gradient-norm\n");
659: PetscPrintf(PETSC_COMM_WORLD,"history = [\n");
660: }
661: PetscPrintf(PETSC_COMM_WORLD,"%3D %g %g %g\n",its,(double)nrm,(double)fct,(double)gnorm);
662: return(0);
663: }
665: PetscErrorCode MonitorDestroy(void **ctx)
666: {
670: PetscPrintf(PETSC_COMM_WORLD,"];\n");
671: return(0);
672: }
675: /*TEST
677: build:
678: requires: !complex
680: test:
681: requires: !single
682: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
684: test:
685: suffix: cn
686: requires: !single
687: args: -ts_type cn -ts_dt .003 -pc_type lu -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none
689: test:
690: suffix: 2
691: requires: !single
692: args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -a .1 -tao_bqnls_mat_lmvm_scale_type none
695: TEST*/